Automatic Mapping of Atrial Fiber Orientations for Patient-Specific Modeling of Cardiac Electromechanics using Image-Registration
Julia M. Hoermann, Martin R. Pfaller, Linda Avena, Crist\'obal, Bertoglio, Wolfgang A. Wall
TL;DR
This paper introduces an automated image registration method to map detailed atrial fiber orientations onto patient-specific heart models, enabling personalized electromechanical simulations without invasive measurements.
Contribution
The study presents a novel, automated approach for defining atrial fiber architecture in patient-specific models using image registration based on an atlas with histologically detailed fibers.
Findings
Accurately mapped fiber orientations in ten patient-specific atria.
Electrophysiological activation patterns in patient models matched atlas characteristics.
Method enabled efficient, personalized electromechanical simulations.
Abstract
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…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37
Figure 38
Figure 39
Figure 40| 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. | / |
| 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 | / | / | ||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsElasticity and Material Modeling · Cardiac electrophysiology and arrhythmias · Cardiovascular Function and Risk Factors
Automatic Mapping of Atrial Fiber Orientations for Patient-Specific Modeling of Cardiac Electromechanics using Image-Registration
Julia M. Hoermann
Martin R. Pfaller
Linda Avena
Cristóbal Bertoglio*†*
and Wolfgang A. Wall*†* Julia M. Hoermann, Martin R. Pfaller and Wolfgang A. Wall are with the Institute of Computational Mechanics, Technical University of Munich, Boltzmannstraße 15, 85748 Garching bei München, Germany. e-mail: [email protected]óbal Bertoglio is with the Bernoulli Institute, University of Groningen, Nijenborgh 9, 9747 AG Groningen, Netherlands.Linda Avena is with the Department of Electrophysiology, German Heart Center Munich, Technical University of Munich, Lazarettstr. 36, 80636 München, Germany.† Joint last authors
Abstract
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 to analogous conclusions for coupled electro-mechano-hemodynamical computations.
{paperlayoutkeywords}
atrial fiber orientation; patient-specific modeling; cardiac electromechanics.
\paperlayoutpeerreviewmaketitle
1 Introduction
\paperlayoutPARstart
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.
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.
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 3a). 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 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.
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 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 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
[TABLE]
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
[TABLE]
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.
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 0.1\text{,}\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{m}\mathrm{s} in fiber direction and a diffusion coefficient of $\sigma_{2}=$0.01\text{\,}\mathrm{m}\mathrm{m}^{2}\mathrm{/}\mathrm{m}\mathrm{s} 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\text{,}\mathrm{mm}\text{/}\mathrm{ms}$$. 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 5a and 5b). 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 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, i.e. it runs a few millimeters further to the left in the mapped case (see Figures 5a and 5b), causing big differences in the error calculation (see Figures 5c, 6). The error is calculated as . Some differences are also visible at the boarders of the circumferential fibers around the veins.
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 7a and 7e).
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).
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 8a, cuts the atria in the plane of the standard 4-chamber-view and Slice 2 (see Figure 8b) 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 12a shows the right atrium and Figure 12b 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 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 10d 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.
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 14h). 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 14d). 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 15a, 15d, and 15i). In other atria, the activation through the Bachmann bundle is less important, since other interatrial connections are nearby (see e.g. Figure 15e, 15g, and 15h). 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).
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.
\appendices
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.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[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[2] J. Zhao, T. D. Butters, H. Zhang, A. J. Pullan, I. J. Le Grice, 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[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[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[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[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[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[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.