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

12/05/2018 ∙ by Julia M. Hoermann, et al. ∙ 0

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.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 3

page 5

page 6

page 7

page 8

page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 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 e.g. studies of the left ventricle show that different fiber architecture lead to different results in the mechanical contraction due to 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]. Due to individual fiber bundles running in different directions it is not straight-forward 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 non-invasively 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 thick [18], while in comparison the left ventricle is around 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 [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 semi-automatic 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 semi-automatic 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 to 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

[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 to the defined fibers. Additionally, we compare, analyze and discuss the results of the electromechanical simulation with mapped fibers in all atria.

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

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

Figure 1: Processing pipeline for the registration and fiber estimation for the patients’ atria

2.2 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 . 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, i.e. 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.

(a) Posterior view
(b) Anterior view
Figure 2: Fiber direction of atlas atria. The colors represent the regions with different fiber directions.

2.3 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, United States)

2.3.1 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 (a)a

). First, we compute a rigid motion that aligns optimally in a least square sense 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

(a)a shows the atria of Patient 1 and the atlas atria after calculating and applying the rigid and the affine transformation to the atlas atria.

2.3.2 Registration

For the registration we use an algorithm which 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 (see Figure (b)b). 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 and the image of the atlas atria by , where is the 3D image cube. We want to find a transformation such that the deformed atlas is similar to the patient’s atria . We use the standard distance measure sum of squared differences (SSD), which is given by

(1)

with 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: find a displacement field such that

(2)

with the regularization parameter . We use a regularization parameter of in this work as suggested in [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 is then linearly interpolated to the grid at level 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 is determined. In our case we use levels of resolution.

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

(a)
(b)
Figure 3: (a) The landmarks defined on the atlas atria. (b) The binary image of the atlas.
(a)
(b)
(c)
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.

2.4 Electro-mechanico-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 [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 [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) methods [34] with a maximal order of five and a stabilization parameter of . For time discretization we use a semi-implicit discretization [35, 36], i.e. linear terms implicit and the non-linear parts explicit in time with at time step of .

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 modelled as nearly incompressible Mooney-Rivlin material, where a volumetric penalization technique for the incompressibility was used. Rayleigh-type damping was also included. The active part is described by an active stress component [37, 38]. The maximal active tension is defined as . The computation of the mechanical model is performed using tetrahedral quadratic continuous finite elements for space discretization and the generalized- method for time discretization.

For the hemodynamical model we use a three-element Windkessel model [39], which is coupled monolithically with the elastic myocardial wall [40]. 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.

3 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 ten different patients’ atria. We enumerate the atria of the patients by Patient 1 to Patient 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 ten differently shaped atria are generated and the electrophysiological activation pattern and the contractions are computed.

3.1 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 (Figures (a)a and (b)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 to four pulmonary veins of Patient 1, the mapped fibers in this area run in circumferential direction around the orifices in both cases (i.e. in case of the atlas registered to Patient 1 as well as in case of the Patient 1 registered to the atlas atria).

Figure (c)c 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, i.e. it runs a few millimeters further to the left in the mapped case (see Figures (a)a and (b)b), causing big differences in the error calculation (see Figures (c)c, 6). The error is calculated as . Some differences are also visible at the boarders of the circumferential fibers around the veins.

(a)
(b)
(c)
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 .
(a)
(b)
Figure 6: Atlas fibers. Difference in fiber direction between defined and mapped fiber orientation.

3.1.2 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 due to the aforementioned differences in the fiber orientation. Characteristics in the activation pattern which 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 Figures (a)a and (e)e).

The transfer of the activation characteristics from the defined to the registered atria are 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 and for Patient 1 . 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 due to 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).

(a) Patient 1 mapped fibers
(b) Patient 1 defined fibers
(c) Patient 1 activation time difference
(d) Atlas mapped fibers
(e) Atlas defined fibers
(f) Atlas activation time difference
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 . (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.

3.1.3 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 9-10. 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 8 shows the position of the slices. Slice 1, visualized in Figure (a)a, cuts the atria in the plane of the standard 4-chamber-view and Slice 2 (see Figure (b)b) 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 due to the difference in electrical propagation speed.

To see the influence of the different activation we plotted the volume of actively contracting tissue, i.e. 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 (a)a shows the right atrium and Figure (b)b 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 (a)a). The slight shift of the fibers of the crista terminalis on the posterior side of the right atrium (see Paragraph 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 10). Especially in Figure (d)d 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 9 and 10). 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.

(a) Slice 1
(b) Slice 2
Figure 8: Illustrating the slices through the atria which are used to show the contraction. (a) Slice through the atria in the plane of the 4-chamber-view. (b) Slice through the atria in a plane parallel to the heart skeleton.
(a) Mapped fibers
(b) Defined fibers
(c)
(c) View 1
(d) View 2
Figure 9: Deformation difference between mapped and defined fibers in the atria for Patient 1 at the time of maximal contraction of the left atrium (time = ). (a) and (b) Displacement of the atria for mapped and defined fibers. (c) and (d) Slices through the right and left atrium in the plane shown in Figure (a)a and Figure (b)b, respectively. The black contour shows the atrium in the relaxed state, the green contour the contracted atria with defined fibers and the red contour line the contracted atria with mapped fibers.
(a) Mapped fibers
(b) Defined fibers
(c)
(c) View 1
(d) View 2
Figure 10: 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 = ). (a) and (b) Displacement of the atria for mapped and defined fibers. (c) and (d) Slices through the right and left atrium in the plane shown in Figure (a)a and Figure (b)b, respectively. The black contour shows the atrium in the relaxed state, the green contour the contracted atria with defined fibers and the red contour line the contracted atria with mapped fibers.
(a) Volume Atlas
(b) Pressure Atlas
(c) Volume Patient 1
(d) Pressure Patient 1
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).
(a) Right Atrium
(b) Left Atrium
Figure 12: Volume of active contracting tissue over time for the atlas atria with defined and mapped fibers
(a) Volume Left Atrium
(b) Pressure Left Atrium
(c) Volume Right Atrium
(d) Pressure Right Atrium
Figure 13: Volume (V) and pressure (p) curves plotted over time of the right (R) and left (L) atrium for Patient 2-10.
Activation Time [s] Start Activation [s] Max/Min Volume [ml] EF [%] Max Pressure [mmHg] Time Min. Volume [s]
Left
Atlas map. /
Atlas def. /
Patient 1 map. /
Patient 1 def. /
Right
Atlas map. /
Atlas def. /
Patient 1 map. /
Patient 1 def. /
Table 1: Summary of electrophysiology and contraction outputs for defined and mapped fiber orientations in atlas and Patient 1 geometries.

3.2 Performance Study

In Figure 14 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, i.e. every atria has the same fibers bundles which run in the same directions, for example the crista terminalis, the pectinate muscles, the circumferential vestibule fibers etc. Unusual shapes as in Patient 9 are handled well (see Figure (h)h). 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 (d)d). It is important to remark that no unphysiological fiber orientations were detected.

The results of the electrophysiological simulations are shown in Figure 15, where the activation patterns of all patients’ atria are shown. It is visible, that the characteristics of the electrophysiological activation pattern which appear in the atlas atria also appear in the patients’ atria. One example is the increased propagation speed due to 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 due to a more or less distinct triangular shape of the activation pattern in this region in each patient’s atria (see Figure 15). 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 Figures (a)a, (d)d, and (i)i). In other atria, the activation through the Bachmann bundle is less important, since other interatrial connections are nearby (see e.g. Figure (e)e, (g)g, and (h)h). Both scenarios are physiological [41]. The region which 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 16 shows the contour lines of the atria of Patient 2-10 at end-systole. All atria perform a physiologically reasonable contraction. The volume and pressure curves until maximal contraction are plotted in Figure 13. 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 due to sustained atrial arrhythmia and atrial volume is high in some of the patients’ atria (see Figure 13 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 due to a high initial volume (see Table 2).

Left Right
Max/Min Volume [ml] EF [%] Max/Min Volume [ml] EF [%]
Patient 2 / /
Patient 3 / /
Patient 4 / /
Patient 5 / /
Patient 6 / /
Patient 7 / /
Patient 8 / /
Patient 9 / /
Patient 10 / /
Table 2: Maximal and minimal volume and ejection fraction for Patient 2-10
(a) Patient 2
(b) Patient 3
(c) Patient 4
(d) Patient 5
(e) Patient 6
(f) Patient 7
(g) Patient 8
(h) Patient 9
(i) Patient 10
Figure 14: Mapped fiber orientation for Patient 2-10.
(a) Patient 2
(b) Patient 3
(c) Patient 4
(d) Patient 5
(e) Patient 6
(f) Patient 7
(g) Patient 8
(h) Patient 9
(i) Patient 10
(j)
Figure 15: Activation times for the patients’ atria Patient 2 to Patient 10 with mapped fiber orientation.
(a) Patient 2
(b) Patient 3
(c) Patient 4
(d) Patient 5
(e) Patient 6
(f) Patient 7
(g) Patient 8
(h) Patient 9
(i) Patient 10
Figure 16: Displacements at maximal contraction of the left ventricle for Patient 2-10 in the 2-chamber-view. The black contour shows the relaxed atria and red contour shows the contracted atria.

4 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 ten 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, i.e. 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 is 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, i.e. build a new atlas atria with the correct wall thickness and segment the wall of patients’ atria more detailed.

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

Acknowledgment

The authors would like to thank the Centers for Radiology of Klinikum rechts der Isar, German Heart Center, Munich, and Kings Collage London, for image data used in this work. Thanks goes also to Martina Weigl and Jonas Niemeijer for their help in segmenting the geometries.

References

  • [1] L. Clerc, “Directional differences of impulse spread in trabecular muscle from mammalian heart.” The Journal of Physiology, vol. 255, no. 2, pp. 335–346, Feb. 1976.
  • [2] J. Zhao, T. D. Butters, H. Zhang, A. J. Pullan, I. J. LeGrice, G. B. Sands, and B. H. Smaill, “An Image-Based Model of Atrial Muscular Architecture,” Circulation: Arrhythmia and Electrophysiology, vol. 5, no. 2, pp. 361–370, Apr. 2012.
  • [3] T. Eriksson, A. Prassl, G. Plank, and G. Holzapfel, “Influence of myocardial fiber/sheet orientations on left ventricular mechanical contraction,” Mathematics and Mechanics of Solids, vol. 18, no. 6, pp. 592–606, Aug. 2013.
  • [4] A. Nikou, R. C. Gorman, and J. F. Wenk, “Sensitivity of left ventricular mechanics to myofiber architecture: A finite element study,” Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, vol. 230, no. 6, pp. 594–598, Jun. 2016.
  • [5] D. D. Streeter, H. M. Spotnitz, D. P. Patel, J. Ross, and E. H. Sonnenblick, “Fiber orientation in the canine left ventricle during diastole and systole,” Circulation research, vol. 24, no. 3, pp. 339–347, 1969.
  • [6] D. B. Ennis, T. C. Nguyen, J. C. Riboh, L. Wigström, K. B. Harrington, G. T. Daughters, N. B. Ingels, and D. C. Miller, “Myofiber angle distributions in the ovine left ventricle do not conform to computationally optimized predictions,” Journal of Biomechanics, vol. 41, no. 15, pp. 3219–3224, Nov. 2008.
  • [7] S. Y. Ho, R. H. Anderson, and D. Sánchez-Quintana, “Atrial structure and fibres: Morphologic bases of atrial conduction,” Cardiovascular research, vol. 54, no. 2, pp. 325–336, May 2002.
  • [8] J. D. Bayer, R. C. Blake, G. Plank, and N. A. Trayanova, “A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models,” Annals of Biomedical Engineering, vol. 40, no. 10, pp. 2243–2254, Oct. 2012.
  • [9] J. Wong and E. Kuhl, “Generating fiber orientation maps in human heart models using Poisson interpolation,” Computer methods in biomechanics and biomedical engineering, vol. 17, no. 11, pp. 1217–1226, Aug. 2014.
  • [10] E. W. Hsu, A. L. Muzikant, S. A. Matulevicius, R. C. Penland, and C. S. Henriquez, “Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation,” American Journal of Physiology - Heart and Circulatory Physiology, vol. 274, no. 5, pp. H1627–H1634, May 1998.
  • [11] P. A. Helm, H.-J. Tseng, L. Younes, E. R. McVeigh, and R. L. Winslow, “Ex vivo 3D diffusion tensor imaging and quantification of cardiac laminar structure,” Magnetic Resonance in Medicine, vol. 54, no. 4, pp. 850–859, Oct. 2005.
  • [12] J. M. Peyrat, M. Sermesant, X. Pennec, H. Delingette, C. Xu, E. R. McVeigh, and N. Ayache, “A Computational Framework for the Statistical Analysis of Cardiac Diffusion Tensors: Application to a Small Database of Canine Hearts,” IEEE Transactions on Medical Imaging, vol. 26, no. 11, pp. 1500–1514, Nov. 2007.
  • [13]

    A. Nagler, C. Bertoglio, M. Gee, and W. Wall, “Personalization of cardiac fiber orientations from image data using the unscented kalman filter,” in

    Functional Imaging and Modeling of the Heart.   Springer, 2013, pp. 132–140.
  • [14] D. E. Sosnovik, R. Wang, G. Dai, T. G. Reese, and V. J. Wedeen, “Diffusion MR tractography of the heart,” Journal of Cardiovascular Magnetic Resonance, vol. 11, no. 1, p. 47, Nov. 2009.
  • [15] C. T. Stoeck, C. von Deuster, M. Genet, D. Atkinson, and S. Kozerke, “Second-order motion-compensated spin echo diffusion tensor imaging of the human heart,” Magnetic Resonance in Medicine, vol. 75, no. 4, pp. 1669–1676, Apr. 2016.
  • [16] N. Toussaint, C. T. Stoeck, T. Schaeffter, S. Kozerke, M. Sermesant, and P. G. Batchelor, “In vivo human cardiac fibre architecture estimation using shape-based diffusion tensor processing,” Medical Image Analysis, vol. 17, no. 8, pp. 1243–1255, Dec. 2013.
  • [17] A. Nagler, C. Bertoglio, C. T. Stoeck, S. Kozerke, and W. A. Wall, “Maximum likelihood estimation of cardiac fiber bundle orientation from arbitrarily spaced diffusion weighted images,” Medical Image Analysis, vol. 39, pp. 56–77, Jul. 2017.
  • [18] R. Beinart, S. Abbara, A. Blum, M. Ferencik, K. Heist, J. Ruskin, and M. Mansour, “Left Atrial Wall Thickness Variability Measured by CT Scans in Patients Undergoing Pulmonary Vein Isolation,” Journal of Cardiovascular Electrophysiology, vol. 22, no. 11, pp. 1232–1236, Nov. 2011.
  • [19] N. Kawel, E. B. Turkbey, J. J. Carr, J. Eng, A. S. Gomes, W. G. Hundley, C. Johnson, S. C. Masri, M. R. Prince, R. J. van der Geest, J. A. C. Lima, and D. A. Bluemke, “Normal Left Ventricular Myocardial Thickness for Middle-Aged and Older Subjects With Steady-State Free Precession Cardiac Magnetic Resonance,” Circulation: Cardiovascular Imaging, vol. 5, no. 4, pp. 500–508, Jul. 2012.
  • [20] F. Pashakhanloo, D. A. Herzka, H. Ashikaga, S. Mori, N. Gai, D. A. Bluemke, N. A. Trayanova, and E. R. McVeigh, “Myofiber architecture of the human atria as revealed by submillimeter diffusion tensor imaging,” Circulation: Arrhythmia and Electrophysiology, vol. 9, no. 4, p. e004133, 2016.
  • [21] F. B. Sachse, R. Frech, C. D. Werner, and O. Dossel, “A model based approach to assignment of myocardial fibre orientation,” in Computers in Cardiology, 1999.   IEEE, 1999, pp. 145–148.
  • [22] B. D. F. Hermosillo, “Semi-automatic enhancement of atrial models to include atrial architecture and patient specific data: For biophysical simulations,” in 2008 Computers in Cardiology, Sep. 2008, pp. 633–636.
  • [23] M. W. Krueger, V. Schmidt, C. Tobón, F. M. Weber, C. Lorenz, D. U. J. Keller, H. Barschdorf, M. Burdumy, P. Neher, G. Plank, K. Rhode, G. Seemann, D. Sanchez-Quintana, J. Saiz, R. Razavi, and O. Dössel, “Modeling Atrial Fiber Orientation in Patient-Specific Geometries: A Semi-automatic Rule-Based Approach,” in Functional Imaging and Modeling of the Heart, ser. Lecture Notes in Computer Science.   Springer, Berlin, Heidelberg, May 2011, pp. 223–232.
  • [24] F. Vadakkumpadan, H. Arevalo, C. Ceritoglu, M. Miller, and N. Trayanova, “Image-Based Estimation of Ventricular Fiber Orientations for Personalized Modeling of Cardiac Electrophysiology,” IEEE Transactions on Medical Imaging, vol. 31, no. 5, pp. 1051–1060, May 2012.
  • [25] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, Sep. 2009.
  • [26] A. Bel-Brunon, S. Kehl, C. Martin, S. Uhlig, and W. A. Wall, “Numerical identification method for the non-linear viscoelastic compressible behavior of soft tissue using uniaxial tensile tests and image registration – Application to rat lung parenchyma,” Journal of the Mechanical Behavior of Biomedical Materials, vol. 29, pp. 360–374, Jan. 2014.
  • [27] B. Fischer and J. Modersitzki, “Ill-posed medicine—an introduction to image registration,” Inverse Problems, vol. 24, no. 3, p. 034008, Jun. 2008.
  • [28] E. Haber and J. Modersitzki, “A Multilevel Method for Image Registration,” SIAM Journal on Scientific Computing, vol. 27, no. 5, pp. 1594–1607, Jan. 2006.
  • [29] H. Lester and S. R. Arridge, “A survey of hierarchical non-linear medical image registration,” Pattern Recognition, vol. 32, no. 1, pp. 129–149, Jan. 1999.
  • [30] D. C. Alexander, C. Pierpaoli, P. J. Basser, and J. C. Gee, “Spatial transformations of diffusion tensor magnetic resonance images,” IEEE Transactions on Medical Imaging, vol. 20, no. 11, pp. 1131–1139, Nov. 2001.
  • [31] J. M. Hörmann, C. Bertoglio, A. Nagler, M. R. Pfaller, F. Bourier, M. Hadamitzky, I. Deisenhofer, and W. A. Wall, “Multiphysics Modeling of the Atrial Systole under Standard Ablation Strategies,” Cardiovascular Engineering and Technology, vol. 8, no. 2, pp. 205–218, Jun. 2017.
  • [32] A. Bueno-Orovio, E. M. Cherry, and F. H. Fenton, “Minimal model for human ventricular action potentials in tissue,” Journal of Theoretical Biology, vol. 253, no. 3, pp. 544–560, Aug. 2008.
  • [33] C. Lenk, F. M. Weber, M. Bauer, M. Einax, P. Maass, and G. Seeman, “Initiation of atrial fibrillation by interaction of pacemakers with geometrical constraints,” Journal of Theoretical Biology, vol. 366, pp. 13–23, Feb. 2015.
  • [34] J. M. Hoermann, C. Bertoglio, M. Kronbichler, M. R. Pfaller, R. Chabiniok, and W. A. Wall, “An adaptive hybridizable discontinuous Galerkin approach for cardiac electrophysiology,” International Journal for Numerical Methods in Biomedical Engineering, Feb. 2018.
  • [35] J. P. Whiteley, “An Efficient Numerical Technique for the Solution of the Monodomain and Bidomain Equations,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 11, pp. 2139–2147, Nov. 2006.
  • [36] M. A. Fernández and N. Zemzemi, “Decoupled time-marching schemes in computational cardiac electrophysiology and ECG numerical simulation,” Mathematical Biosciences, vol. 226, no. 1, pp. 58–75, Jul. 2010.
  • [37] J. Bestel, F. Clément, and M. Sorine, “A Biomechanical Model of Muscle Contraction,” in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2001, ser. Lecture Notes in Computer Science, W. J. Niessen and M. A. Viergever, Eds.   Springer Berlin Heidelberg, Oct. 2001, no. 2208, pp. 1159–1161.
  • [38] D. Chapelle, P. Le Tallec, P. Moireau, and M. Sorine, “Energy-preserving muscle tissue model: Formulation and compatible discretizations,” International Journal for Multiscale Computational Engineering, vol. 10, no. 2, 2012.
  • [39] Y. Shi, P. Lawford, and R. Hose, “Review of Zero-D and 1-D Models of Blood Flow in the Cardiovascular System,” BioMedical Engineering OnLine, vol. 10, p. 33, Apr. 2011.
  • [40] M. Hirschvogel, M. Bassilious, L. Jagschies, S. Wildhirt, and M. Gee, “A monolithic 3d-0d coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics,” International Journal for Numerical Methods in Biomedical Engineering, vol. 33, no. 8, p. e2842, 2016.
  • [41] V. Markides, R. J. Schilling, S. Y. Ho, A. W. C. Chow, D. W. Davies, and N. S. Peters, “Characterization of Left Atrial Activation in the Intact Human Heart,” Circulation, vol. 107, no. 5, pp. 733–739, Feb. 2003.