Automatic Mapping of Atrial Fiber Orientations for Patient-Specific Modeling of Cardiac Electromechanics using Image-Registration

Knowledge of appropriate local fiber architecture is necessary to simulate patient-specific electromechanics in the human heart. However, it is not yet possible to reliably measure in-vivo fiber directions, especially in human atria. Thus, we present a method which defines the fiber architecture in arbitrarily shaped atria using image registration and reorientation methods based on atlas atria with fibers predefined from detailed histological observations. Thereby, it is possible to generate detailed fiber families in every new patient-specific geometry in an automated, time-efficient process. We demonstrate the good performance of the image registration and fiber definition on ten differently shaped human atria. Additionally, we show that characteristics of the electrophysiological activation pattern which appear in the atlas atria also appear in the patients' atria. We arrive at analogous conclusions for coupled electro-mechano-hemodynamical computations.

the left and right atrial wall. 7 Because of individual fiber bundles running in different directions, it is not straightforward to use rule-based approaches, as it can be done in the ventricles. 8,9 Besides rule-based methods, another promising approach to define the fiber direction in ventricles is to use diffusion tensor MRI (DTMRI), which is capable to measure noninvasively the fiber architecture of the left ventricular myocardium. 10 This technology is often used ex vivo, [11][12][13] since in vivo measurements are challenging 14,15 and only few slices can be acquired. Furthermore, it requires sophisticated reconstructions of the fibers for the whole ventricles. 16,17 However, until now, it is not possible to obtain in vivo fiber directions in the atria since their thickness is smaller than current DTMRI voxel size. Precisely, the atrial wall is around 2 mm thick, 18 while in comparison, the left ventricle is around 8 mm thick. 19 Only recently, ex vivo fiber orientation in eight different atria could be analyzed with submillimeter DTMRI, 20 which could be used as additional information for fiber definitions in the future.
Until now, only few methods have been proposed to create fiber directions in patient-specific atria. [21][22][23] The approach in Krueger et al 23 is a first step towards realistic atrial models. The importance of fiber orientations for patient-specific simulations is demonstrated with different geometrical models, to which the semiautomatic method is applied and additional electrophysiological simulations are performed. The method uses voxel-based atrial images and manually defined seed points to specify different fiber bundles using marching level sets methods. It is shown that this method works very robust, and the results correspond well with reported literature. However, the semiautomatic approach is strongly depending on user input of the 22 seed points, which have to be set in an accurate way. Variation of the seed points can lead to different results. Strong shape variations, which are common in human atria, are difficult to handle with this algorithm, since it depends on shortest paths between seed and auxiliary points and subdivisions at fixed relations. For example, to incorporate absent or additional pulmonary vein orifices, an adaption of the algorithm is necessary. Additionally, these models are limited to the information about anatomical fiber orientations known at the time of the development of the methods. New information can be only incorporated into the model by changing the methodology, which in turn would lead to more user interaction by defining additional seed points.
To the authors' best knowledge, until now, image registration techniques are not yet used for fiber estimation in the atria. We propose a method to map atlas atria to a patient's atria using image registration techniques. Using the computed deformation map, the fiber directions are reoriented and then transferred to the patient's atria. The initial user input of seed point is reduced in comparison with rule-based models and is only needed for a first general alignment of the geometries. The influence of user variations is therefore greatly decreased. Additionally, the benefit of using registration methods is that the accuracy of the fiber orientation can be easily improved by adapting only the atlas model. More complex data and details have to be included only in the atlas model without changing the registration algorithm. Additionally, geometry variations as the number of pulmonary orifices, which is a challenge for the rule-based models, can be handled by the usage of different atlas atria. Another benefit is the possibility of using ex vivo measured DTMRI atria with fiber data as atlas atria, thus, improving the accuracy of the model. In Vadakkumpadan et al, 24 a method has been proposed to register an atlas ventricle to the patient's ventricle using large deformation diffeomorphic metric mapping. The goal of this paper is to propose an image registration and reorientation method for atrial geometries and fibers and to demonstrate its performance solving an electromechanical problem of patient-specific atria using our fiber definition.
The reminder of the paper is organized as follows. In Section 2, we characterize the method to define the fibers of the atlas atria. Then we describe the process of image registration, fiber reorientation, and fiber lifting to generate the fibers in the patient's atria. Additionally, we shortly describe the electro-mechano-hemodynamical model, with which we simulate the functions of the atria with defined and mapped fibers. In Section 3, we describe the results of the registration and mapping procedure in comparison with the defined fibers. Additionally, we compare, analyze, and discuss the results of the electromechanical simulation with mapped fibers in all atria.

METHODS
Several steps are necessary for the estimation of the fiber architecture in patient's atria. Firstly, an atlas atria with a physiologically detailed fiber architecture has to be created. This is done once at the beginning, and the atlas is then used for fiber estimation for all atrial geometries. For each patient, the following procedure is performed; see Figure 1. From medical imaging data, a geometry is created, which is registered to the atlas. Then the deformed atlas with reoriented fibers is used to generate the fibers of the patient. Finally, the fibers are realigned so that they are tangential to the surface and at last a harmonic lifting is performed.
At the end of this chapter, we will also briefly describe how computations using an electro-mechano-hemodynamical model are performed with the results of the fiber generation as geometrical input.

Geometry creation
To construct the geometry of a patient's atria, we use segmentation, design modification, and meshing tools. First, a surface representation is generated using the software Mimics (Materialise, Leuven, Belgium). For this, the lumen of both atria is segmented manually so that the endocardial surfaces are obtained. To generate the epicardial surface, we extrude the endocardial surface by 2 mm, which corresponds to an average thickness of the atrial wall. 18 Additionally, we add an interatrial muscular bridge between the right and the left atrium, the Bachmann bundle, to allow a physiological propagation of the electrical signal. Finally, we create a 3D volume mesh with tetrahedral elements with a maximal element size of 0.9 mm using Gmsh. 25 This leads to about two to three elements through the atrial wall. Note that this does not pose accuracy issues in the computations, since we use higher order spatial discretizations.

Fiber definition for atlas atria
For the atlas atria, we define the fiber orientation using reported detailed knowledge of atrial fiber structure. 7 The geometry of the atlas are the atria of a 71-old individual without known cardiac pathological findings. The model was obtained from a cardiac CT image with a resolution of 0.4 × 0.4 × 0.8 mm. The atlas geometry is then created identically as any other patients' atria and as has just been described in Section 2.1.
To define the fibers in the atlas atria, we divide the epicardial and endocardial wall of the left and right atrium into regions with a constant fiber direction (see colored regions in Figure 2). In doing so, we manually set the main fiber bundles crista terminalis, pectinate muscles, circumferential vestibule fibers, and the Bachmann bundle. Additionally, also other prominent bundles and fiber alignment in the right and left atrium are defined. In the right atrium, endocardial and epicardial fiber directions coincide, ie, the fiber bundle direction is constant throughout the whole thickness of the wall. In contrary, different fiber bundles throughout the thickness of the wall run in the left atrium at the posterior wall. To model this, we assign different fiber directions on the epicardial and endocardial surface. After defining a fiber direction on each node on the surface, we interpolate the fibers into the volume using harmonic lifting techniques. 13 The atlas atria with fiber directions and the different regions can be seen in Figure 2.

Atlas geometry registration
The deformation of the atlas atria to the shape of the patients' atria is done in two major steps. First, an affine transformation is calculated and second the actual elastic registration process is computed. The registration process is performed in MATLAB (Release 2017a, The MathWorks, Inc, Natick, Massachusetts)

Affine transformation
For the affine transformation, 13 landmarks on the atlas and the patients' geometry are defined around the orifices of pulmonary veins and around the valvular orifices to the ventricles and the Bachmann bridge ( Figure 3A). First, we compute a rigid motion that aligns optimally in a least square sense of the two sets of landmark points using singular value decomposition. Note that the landmarks are only used for the rigid transformation, which is needed to generally align the two geometries. Furthermore, three landmarks are sufficient for a unique rigid transformation. Second, we perform an iterative closest point affine registration for the sets of all mesh points of the geometries. Figure 4A shows the atria of patient 1 and the atlas atria after calculating and applying the rigid and the affine transformation to the atlas atria.

Registration
For the registration, we use an algorithm that we already successfully used for material modeling and parameter identification of biological materials. 26 To match the atlas geometry to the patients geometry, we first create a 3D binary voxel representation of the atria with a voxel size of 0.6 mm (see Figure 3B). The walls in the binary image are slightly thicker than the original walls in the mesh. Nevertheless, this does not yield a problem, since the interpolation is performed with a weighted nearest neighbor principle. This is done for both the atlas atria and the patient's atria. We denote the image of the patient's atria by I p ∶ Ω → {0, 1} and the image of the atlas atria by I a ∶ Ω → {0, 1}, where Ω ⊂ R 3 is the 3D image cube. We want to find a transformation such that the deformed atlas I a is similar to the patient's atria I p . We use the standard distance measure sum of squared differences (SSD), which is given by with (x) = x + u(x) and u(x) ∶ R 3 → R 3 a spatially varying displacement field. To overcome the inherent ill-posedness of the image registration problem, 27 an elastic potential as regularization  is added. 28 Therefore, we formulate the registration problem as the following: Find a displacement field u * such that with the regularization parameter > 0. We use a regularization parameter of = 0.1 in this work as suggested in Bel-Brunon et al. 26 To solve the optimization problem, we use a Gauss-Newton method. 28 Additionally, we use a multiresolution approach. 29 From the binary images, several levels of coarse grids are consecutively created using convolution with a Gaussian smoothing function. The smoothed image is hereby resampled to an image at a coarser scale with doubled voxelsize. We start the minimization at the coarsest level. The result of the minimization at level n is then linearly interpolated to the grid at level n − 1, and this result is used as a starting point at the new level. This minimization-interpolation steps are performed at each level until at the finest level, the original binary image, the final transformation u * is determined.
In our case, we use n = 4 levels of resolution.

Interpolation and atlas fiber reorientation
The last step in the fiber estimation process is to transfer the fibers of the atlas atria to the patient's atria. To do this, we first calculate the transformation of each mesh node from the optimal transformation u * of the voxels. For each mesh node, the nearest voxels are determined. Using the transformation defined in these voxels, the transformation of the mesh node is computed using the inverse distance weighting average of the voxel centers. Additionally, the deformation gradient is calculated at each node using the affine transformation matrix and the Jacobian of the displacement field u. To reorient the fiber directions, two strategies are possible, the finite strain strategy and the strategy of principal directions. 30 The finite strain strategy uses the rotational component of the deformation gradient to reorient the fibers. Whereas the strategy of preservation of principal directions takes also into account the deformation component of the affine transformation. Thus, the whole deformation gradient is used for the reorientation. In this paper, we use the strategy of principal direction. To map the fiber orientation of the deformed atlas atria to the patient's atria, we use again an inverse distance weighting of the mesh nodes of the geometries. We map only the fiber orientation to the surface of the patient's atria. Then we project the fiber orientation to the tangential plane of the surface nodes and perform afterwards a harmonic lift step to compute the fiber orientation in the interior of the atrial wall. 13 To improve readability, for now on, we will use the keyword mapped for the fibers, which are generated on the patients' atria through registration and interpolation techniques.
In all our computations, we used n = 4 levels of resolution in the registration algorithm.

Electro-mechano-hemodynamical computations
The details of the electrophysiological, mechanical, and hemodynamical models are described in previous work. 31 For the electrophysiological part, we calculate the monodomain equations with the minimal cell model from Bueno-Orovio et al. 32 The parameter set of the cell model is adapted to reproduce atrial activation and a physiological action potential curve for the atrial cell according to Lenk et al. 33 The diffusivity is assumed transverse isotropic with a diffusion coefficient of 1 = 0.1mm 2 ∕ms in fiber direction and a diffusion coefficient of 2 = 0.01 mm 2 ∕ms perpendicular to it. To receive physiological propagation, we use high-order hybridizable discontinuous Galerkin (HDG) methods 34 with a maximal order of five and a stabilization parameter of = 1 mm/ ms. For time discretization, we use a semi-implicit discretization, 35,36 ie, linear terms implicit and the nonlinear parts explicit in time with at time step of 0.1 ms. The mechanical model is coupled to the electrical model via the electrical signal, which triggers the contraction. The model is based on nonlinear elastodynamic equations with a passive and an active part.
The passive material is modeled as nearly incompressible Mooney-Rivlin material, with constants C 1 = 20 kPa and C 2 = 40 Pa where a volumetric penalization technique for the incompressibility was used, 37 ie, (J + 1∕J − 2) 2 with J the Jacobian of the deformation gradient and = 10 4 kPa. This lead to a volume change of 0.3 to 0.75 % from diastole to systole, depending on the atrial model.
Rayleigh-type damping was also included. The Rayleigh damping parameters are 5.0 for the stiffness and 0.0 for the mass damping. The active part is described by an active stress component. 38,39 The maximal active tension is defined as 100 kPa.
The computation of the mechanical model is performed using tetrahedral quadratic continuous finite elements for space discretization and the generalized-method for time discretization, with parameters f = 0.5, m = 0.5, = 0.25, and = 0.5 according to Pfaller et al. 40 For the hemodynamical model, we use a three-element Windkessel model, 41 which is coupled monolithically with the elastic myocardial wall. 42 We used Windkessel models at the orifices to the ventricles. For all atrial models, we used for the sake of simplicity the same set of parameters, which are adjusted for the patient 1 using ventricular ejection fraction captured using cine MRI images.

RESULTS
To demonstrate the performance of our method, we investigate on one hand the difference between mapped and defined fibers, and on the other hand, we demonstrate the functionality of our method on 10 different patients' atria. We enumerate the atria of the patients by patients 1 to 10.
In the following part A, we show a comparison between mapped fibers and defined fibers. For that, we manually define for patient 1 the fiber orientation, performing the same steps as for the atlas atria described in Section 2.2. Then we first map the atlas fibers to patient 1, and second, we map the fibers of patient 1 to the atlas atria. Afterwards, we compare the fiber orientations, the electrophysiological activation, and the mechanical contraction between the two pairs of mapped fibers and manually defined fibers.
In the second part, the fiber orientation for 10 differently shaped atria are generated and the electrophysiological activation pattern and the contractions are computed.

Comparison of defined and mapped fibers 3.1.1 Fibers
The results of the fiber estimation in patient 1 and the atlas are shown in Figures 5 and 6, respectively. The mapped fibers are overall arranged in quite a similar way as the defined fibers ( Figure 5A,B). Between the superior vena cava and the inferior vena cava, the fibers in the crista terminalis are aligned longitudinally, while the pectinate muscles lay perpendicular to them. Around the vestibulum and the orifices of the veins, the fibers run in circumferential direction. Although the atlas atria has three pulmonary veins compared with four pulmonary veins of patient 1, the mapped fibers in this area run in circumferential direction around the orifices in both cases (ie, in case of the atlas registered to patient 1 as well as in case of the patient 1 registered to the atlas atria). Figure 5C shows that the error between mapped and manually defined fibers is small in general (blue color), where larger differences are concentrated at the boundaries between different fiber bundles. For example, the fiber bundle of the crista terminalis is shifted, ie, it runs a few millimeters further to the left in the mapped case (see Figure 5A,B), causing big differences in the error calculation (see Figures 5c and 6). The error is calculated as err = 1 − | T mapped defined |. Some differences are also visible at the boarders of the circumferential fibers around the veins.

Electrophysiology
In Figure 7, the results of the electrophysiological simulations are shown for the atlas atria and patient 1 atria. Here, the activation times obtained from mapped fibers and manually defined fibers are compared. The overall activation pattern is very similar for both fiber architectures: The signal travels fast along the crista terminalis; at the same spots, the activation travels from the right atrium to the left atrium; and the tip of the left appendage is activated at last. Moreover, in both cases, it is clearly visible that the activation of the left atrium occurs through the Bachmann bundle. Differences in the activation can be seen at the borders of the fiber bundles (compare with Figures 2 and 5), as expected because of the aforementioned differences in the fiber orientation. Characteristics in the activation pattern that appear in the atria with defined fibers also appear in the atria with fibers mapped from it. This can be seen for patient 1, where the activation sequence is analogous to the atlas atria with defined fibers (compare Figure 7a and 7e).  The transfer of the activation characteristics from the defined to the registered atria is also visible in the activation time. The atlas atria with defined fibers and the patient 1 atria with mapped fibers take both longer to activate than the patient 1 atria with defined fibers and the atlas with mapped fibers (see Table 1). The maximal difference in activation time is for the atlas atria around 18 % and for patient 1 16 %. In the atlas atria, the region with the biggest difference is the tip of the right appendage, while for patient 1, it is around the left appendage (see Figure 7). These differences are acceptable because of the fact that they appear in the appendages of the atria, which are less important for the ejection fraction (see next paragraph and Table 1).

Mechanical contraction
The mechanical contraction is coupled with the electrophysiology via the transmembrane potential. The results of the contraction between mapped and defined fibers for the atlas atria and Patient 1 are shown in Figures 8 to 9. At the top of each figure, the displacements are shown at the time of maximal contraction of the right and left atrium, respectively. The contour plots in the bottom of each figure show the contour of the atria at differently positioned slices through the atria. Figure 10 shows the position of the slices. Slice 1, visualized in Figure 10A, cuts the atria in the plane of the standard four-chamber view, and slice 2 (see Figure 10B) cuts the atria in a plane parallel to the heart skeleton. We analyzed also the volume curves of the left and right atrium over time until maximal contraction of each chamber (see Figure 11). Additionally, in Figure 11, also the pressures in the right and left atrium are plotted. Only small differences in the volume curves for the atlas atria with mapped and defined fibers are visible. The contraction of the atlas with defined fibers is slightly faster than the contraction of the atlas with mapped fibers because of the difference in electrical propagation speed.  Figure 10A,B, respectively. The black contour shows the atrium in the relaxed state, the green contour shows the contracted atria with defined fibers, and the red contour line shows the contracted atria with mapped fibers  Figure 10A,B, respectively. The black contour shows the atrium in the relaxed state, the green contour shows the contracted atria with defined fibers, and the red contour line shows the contracted atria with mapped fibers To see the influence of the different activation, we plotted the volume of actively contracting tissue, ie, the tissue where the action potential is above a threshold, over time. In the plot, the active tissue is compared for the atlas atria with defined fibers to the atlas atria with mapped fibers. Figure 12A shows the right atrium and Figure 12B the left atrium. It is visible in the plot that for the right atrium, slightly more tissue is active in the atrium with defined fibers than with mapped fibers at the same time. This is consistent with the faster decrease in the volume curve (see Figure 11A). The slight shift of the fibers of the crista terminalis on the posterior side of the right atrium (see Section 3.1.1) in the mapped case leads to fibers in the thickened wall of the crest not running along the crest. This is the reason for the folding near the terminal crest during contraction of the right atrium (see Figure 9). Especially in Figure 9D, the folding of the terminal crest is visible. However, both mapped and defined fibers fold in the regions of fiber orientation change. Differences in the displacement of the left atrial wall of the atlas atria during contraction exist near characteristic shapes, which only exist in individual atria, for example, the pronounced buckle in the inferior wall of the right appendage and the distinctive shape of the tip of the right appendage. Since these characteristics are different and only present in individual atria, the registration process, especially the fiber interpolation, cannot map the correct fiber directions in this regions. However, they are also not known from anatomical studies. In patient 1, we have a very smooth shaped atria. The difference in displacements during contraction is smaller between mapped and defined fibers (see Figures 8 and 9). Also the volume and pressure curves, the ejection fraction, and the time of maximal contraction are similar. Only the displacement of the pulmonary veins, the caval veins, and the left appendage differ slightly. For the atlas atria and the atria of patient 1, characteristic values of the activation and contraction are listed in Table 1.

Performance study
In Figure 13, all patient's atria with mapped fibers are shown. Although the geometry of the atria has a different shape in every patient, the registration is able to deform the atlas shape to the correct patient's shape. Additionally, the fiber (A) (B)

FIGURE 10
Illustrating the slices through the atria which are used to show the contraction. A, Slice through the atria in the plane of the four-chamber view. B, Slice through the atria in a plane parallel to the heart skeleton architecture of the atria is well defined, ie, every atria has the same fibers bundles which run in the same directions, for example, the crista terminalis, the pectinate muscles, and the circumferential vestibule fibers. Unusual shapes as in patient 9 are handled well (see Figure 13H). Although the geometry has a bump on the left atrium, the fiber direction around it is smooth and also the pump is equipped with fibers. Patient 5, which has an additional right pulmonary vein orifice, has physiological fiber orientation in the region of the pulmonary orifices (see Figure 13D). It is important to remark that no unphysiological fiber orientations were detected.
The results of the electrophysiological simulations are shown in Figure 14, where the activation patterns of all patients' atria are shown. It is visible that the characteristics of the electrophysiological activation pattern that appear in the atlas atria also appear in the patients' atria. One example is the increased propagation speed because of circumferential vestibule fibers in the posterior wall of the left atrium. When the signal reaches the fibers around the vestibule on the posterior left atrial wall, the activation in circumferential direction around the vestibule increases, which is recognizable because of a more or less distinct triangular shape of the activation pattern in this region in each patient's atria (see Figure 14). The signal travels from the right atrium to the left atrium through the Bachmann bundle if the septum is small and the distance between Bachmann bundle and other interatrial connections is big (see Figure 14A,D,I). In other atria, the activation through the Bachmann bundle is less important, since other interatrial connections are nearby (see, eg, Figure 14E,G,H). Both scenarios are physiological. 43 The region that is activated last is the left atrial appendage. However, the overall activation time depends of course on the size and shape of the atria.  Figure 15 shows the contour lines of the atria of patients 2 to 10 at end systole. All atria perform a physiologically reasonable contraction. The volume and pressure curves until maximal contraction are plotted in Figure 16. Despite the differences in volume of the atria, all atria show physiological volume and pressure curves over time. The atrial geometries have been obtained from patients receiving treatment for atrial arrhythmia; thus, some of the geometries are dilated because of sustained atrial arrhythmia, and atrial volume is high in some of the patients' atria (see Figure 16 and Table 2). Hence, these cases are a good test for our approach. However, our method is able to register the atria, and to transfer the fiber orientations although the volume of atlas and patient's atria varies significantly. Additionally, also with such dilated geometry, mechanical contraction can be computed without problems. The ejected volume increases slightly with increasing atrial volume, but the ejection fraction decreases in dilated atria because of a high initial volume (see Table 2).

DISCUSSION
The proposed approach is able to register differently shaped atria to each other and to create physiological reasonable fiber architecture for patient-specific geometries. We show this using 10 different human atria with various shapes. Geometric  Using the proposed pipeline to define the fiber orientation on the atria, an exact atlas atria is needed. As the shape of the atria in patients varies very much, an atlas atria should be used which itself is similar to the patient's atria. Although, the method was able to handle a varying number of pulmonary veins, for a better accuracy of the fiber definition at the pulmonary orifices, it is suggested to use an atlas with the common number of pulmonary orifices. As for atrial arrhythmia simulations the pulmonary veins play an important role, it could be more convenient to use different atlases with the appropriate number of pulmonary orifices.
Small differences in fiber orientation lead to small differences in the electrical activation pattern. However, the mechanical contraction showed to be more sensitive to the fiber orientation. Individual shapes of the atria, for example, additional bulges, cause the registration to deform the atlas atria more to match the individual patient atria, which in turn could lead to strongly varying fiber orientation in the bulges. However, since it is not possible to visualize the fiber orientation in vivo in the atria, it is not measurable which is the correct fiber definition in these bulges, which appear only in single individuals. Our registration algorithm is able to define fiber orientation everywhere in the atria, ie, also in bulges. However, to improve the performance of the registration and fiber reorientation, different atlas atria could be used for different shaped atria, for example, from ex vivo DTMRI. 20 With a similarity measurement between the atlas atria and the patients' atria, for instance, the size of the deformation map of the registration, one could find the most similar atlas to be used along with the proposed approach, which would then lead to the most physiological results regarding the fiber orientation and, thus, the activation and contraction.
The segmentation and the generation of the atrial geometry are restricted by the medical image resolution, which is not enough to capture details of the atrial wall thickness. In the future, a more accurate reconstruction of the atrial wall can be performed during the segmentation process if improved medical imaging is available. In that case, the proposed method should probably be adapted in order to deal with that level of detail, ie, build a new atlas atria with the correct wall thickness and segment the wall of patients' atria more detailed.
Concerning the computational time, the whole registration and fiber definition take at the moment between a few hours to around a day, depending on the size of the atria since the voxel size was the same for all geometries. However, our computer code is so far neither optimized nor parallelized. This time also needs to be put in reference to available alternatives, namely, using manual assignment of the fibers, which can take up to a few days of manual work.
Finally, we want to remark that we only reproduced the kinematics at peak atrial systole for an effective assessment of the impact of the errors in the reconstructed fibers. Therefore, we chose the Mooney-Rivlin material for the sake of computational efficiency because of the large amount of electro-fluid-mechanical simulations that needed to be run, since in our experience, exponential materials require considerably more nonlinear iterations (about twice, in cardiac mechanics testcases). In any case, there is only little evidence on the passive behavior of the human atria, 44 showing moreover great variability among the same specimen and among different subjects.

CONCLUSION
A local fiber definition is important for patient-specific electrophysiological and mechanical modeling of the atria. We presented a method to automatically define the fiber orientation on arbitrarily shaped atria. To do so, we use a single pair of atlas atria with a detailed predefined fiber orientation. Using image registration techniques and reorientation methods, we are able to map the fibers to different atria, even if the shape of the atria differ significantly. The method needs as input only the segmented geometry together with a few landmarks and does not require additional user interaction. Our method is thus very convenient, especially if a study with many individual atria is carried out. We compared the result of the fiber mapping with manually defined fibers. The same fiber bundles were visible in defined and mapped atria and the performed electrophysiology and contraction simulation showed similar results for defined and mapped fibers. Using our registration method for fiber mapping to patients' atria, the activation pattern of the patients showed the same characteristics as the atlas. Thus, with a detailed atlas, we have the same activation properties also in patient's atria. We demonstrated the good performance of our method with 10 different human atria. Different shapes of the atria are handled very well during the registration process. The electrophysiological, contraction, and hemodynamical simulation with atria with mapped fibers showed physiologically reasonable results in terms of activation pattern, spatial contraction, volume and pressure curves, and ejection fraction.