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 that 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 10 differently shaped human atria. Additionally, we show that characteristics of the electrophysiological activation pattern that appear in the atlas atria also appear in the patients' atria. We arrive to analogous conclusions for coupled electro-mechano-hemodynamical computations.

INTRODUCTION

Patient‐specific mathematical and computational models can contribute to understand the (patho‐)physiological function of the heart. These models require not only an accurate geometrical representation of the heart, usually obtained from computed tomography (CT) or from magnetic resonance imaging (MRI), but also a description of the fiber directions. The term fiber is referred to myofiber bundles, which are similarly oriented myocytes running along a certain direction denoted as fiber direction. For a correct representation of the electrophysiological behavior, knowledge about the fiber direction is necessary, since the electrical signal travels faster in fiber direction than perpendicular to it,1 and this anisotropy influences the electrical activation.2 But a correct fiber representation is obviously also important from the mechanical point of view, as, eg, studies of the left ventricle show that different fiber architecture lead to different results in the mechanical contraction because of active and passive anisotropy.3, 4

The fiber architecture of the atria differs from the one of the ventricles. While in the ventricles, the fibers are aligned in an almost regular way,5, 6 the fibers in the atria are aligned in individual bundles, which run in different directions through 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 challenging14, 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 al23 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.

Figure 1

Processing pipeline for the registration and fiber estimation for the patients' atria

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.

Figure 2

Fiber direction of atlas atria. The colors represent the regions with different fiber directions

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.

Figure 3

A, The landmarks defined on the atlas atria. B, The binary image of the atlas

Figure 4

Atlas atria (blue) and patient's atria (yellow) at different steps of the registration process: A, meshed geometry after affine transformation; B, binary images after registration; and C, the meshed geometry after the registration and interpolation to the tetrahedral nodes

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 :Ω→{0,1} and the image of the atlas atria by I :Ω→{0,1}, where is the 3D image cube. We want to find a transformation φ such that the deformed atlas I is similar to the patient's atria I . We use the standard distance measure sum of squared differences (SSD), which is given by with φ(x)  =  x + u(x) and 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 in fiber direction and a diffusion coefficient of perpendicular to it. To receive physiological propagation, we use high‐order hybridizable discontinuous Galerkin (HDG) methods34 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 κ  =  104 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 α   =  0.5, α   =  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

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 5

Patient 1 fibers. A, Defined fiber orientations in patient 1. B, Mapped fiber orientation in patient 1. C, Difference in fiber direction between defined and mapped fiber orientation. The difference is calculated as

Figure 6

Atlas fibers. Difference in fiber direction between defined and mapped fiber orientation

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 . 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).

Figure 7

Results of electrophysiology simulations. Activation times of patient 1 atria and atlas atria with mapped fibers and defined fibers and activation time differences. The difference is calculated as diff  =  |(activation time)mapped − (activation time)defined|. A, Activation time with mapped fibers for patient 1. B, Activation time with defined fibers for patient 1. C, Difference between the activation times with mapped fibers to defined fibers of Patient 1. D, Activation time with mapped fibers for the atlas. E, Activation time with defined fibers for the atlas. F, Difference between the activation times with mapped fibers to defined fibers of the atlas atria

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).

Table 1

Summary of electrophysiology and contraction outputs for defined and mapped fiber orientations in atlas and patient 1 geometries

		Activation time, s	Start activation, s	Max/min volume, mL	EF, %	Max pressure, mmHg	Time min volume, s
Left
	Atlas map	0.338	0.098	90.237613 / 65.7015717816899	27.5732948824884	17.8143418467805	0.349
	Atlas def	0.376	0.09	90.237613 / 62.6711304704444	30.8908730610759	19.026517701898	0.352
	Patient 1 map	0.317	0.076	43.08893 / 26.3452440603676	39.5473506623582	14.6974119280259	0.313
	Patient 1 def	0.288	0.062	43.08893 / 26.8412097214347	38.2789084919675	14.4990183282687	0.303
Right
	Atlas map	0.262	0	92.1740 / 64.0666398285335	31.0299927373123	6.4061280870252	0.277
	Atlas def	0.277	0	92.1740 / 62.2938025487618	33.1849331976425	6.49477659640426	0.274
	Patient 1 map	0.212	0	52.2500 / 33.708274894174	36.0112309123391	5.92783762174759	0.242
	Patient 1 def	0.212	0	52.2500 / 33.3374851816822	36.7043375625503	5.94637927986155	0.236

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 8

Deformation difference between mapped and defined fibers in the atria for patient 1 at the time of maximal contraction of the left atrium (time = 0.31 s). A,B, Displacement of the atria for mapped and defined fibers. C,D, Slices through the right and left atrium in the plane shown in 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 9

Deformation difference between mapped and defined fibers in the atria for the atlas atria at the time of maximal contraction of the left atrium (time = 0.35 s). A,B, Displacement of the atria for mapped and defined fibers. C,D, Slices through the right and left atrium in the plane shown in 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 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

Figure 11

Volume (V) and pressure (p) curves plotted over time of the right (R) and left (L) atrium of patient 1 and the atlas geometry for mapped fibers (map) and defined fibers (def)

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.

Figure 12

Volume of active contracting tissue over time for the atlas atria with defined and mapped fibers

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 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.

Figure 13

Mapped fiber orientation for patients 2 to 10

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 14

Activation times for the patients' atria patients 2 to 10 with mapped fiber orientation

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).

Figure 15

Displacements at maximal contraction of the left ventricle for patients 2 to 10 in the two‐chamber view. The black contour shows the relaxed atria, and red contour shows the contracted atria

Figure 16

Volume (V) and pressure (p) curves plotted over time of the right (R) and left (L) atrium for patients 2 to 10

Table 2

Maximal and minimal volume and ejection fraction for patients 2 to 10

	Left		Right
	Max/min volume, mL	EF, %	Max/min volume, ml	EF, %
Patient 2	114.556698310895 / 83.536649356947	27.4392876676423	126.637118441487 / 96.658209475075	24.1814963949822
Patient 3	128.04763883753 / 95.701506216093	25.5692891827454	161.743118686833 / 119.40086204681	26.5429447826017
Patient 4	157.009621214239 / 121.39546503040	22.9586536344844	107.960005404864 / 76.763407051905	29.3700939172109
Patient 5	190.970827833812 / 153.11325215410	20.1183536960575	118.767502059093 / 85.445548383752	28.5108321344938
Patient 6	145.902620178197 / 114.04613013089	22.1017887345541	146.579127704434 / 107.08356339605	27.4397235027787
Patient 7	136.207083329857 / 103.49447908437	24.2328496493578	150.236727659456 / 112.82118863064	25.5681164928298
Patient 8	127.973774311196 / 96.511726565111	24.9034604194291	176.428094330608 / 136.61722359473	23.1800871639543
Patient 9	96.8237672271759 / 69.480117115129	28.7722180439427	147.920611025809 / 110.33602711163	25.8570572105612
Patient 10	158.223428094065 / 126.53332184517	20.3736621858871	136.887713049734 / 100.18882307724	27.3501772965607

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 differences in the shape of the atria, as, for example, the number of pulmonary veins, are handled well. The fibers can be reoriented according to the deformation calculated in the registration process and applied to the patient's atria.

The activation sequence in the patient's atria is similar to the sequence in the atlas atria; thus, it is very important to have an exact fiber definition in the atlas atria.

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.