|Home | About | Journals | Submit | Contact Us | Français|
Technological limitations pose a major challenge to acquisition of myocardial fiber orientations for patient-specific modeling of cardiac (dys)function and assessment of therapy. The objective of this project was to develop a methodology to estimate cardiac fiber orientations from in vivo images of patient heart geometries. An accurate representation of ventricular geometry and fiber orientations was reconstructed, respectively, from high-resolution ex vivo structural magnetic resonance (MR) and diffusion tensor (DT) MR images of a normal human heart, referred to as the atlas. Ventricular geometry of a patient heart was extracted, via semiautomatic segmentation, from an in vivo computed tomography (CT) image. Using image transformation algorithms, the atlas ventricular geometry was deformed to match that of the patient. Finally, the deformation field was applied to the atlas fiber orientations to obtain an estimate of patient fiber orientations. The accuracy of the fiber estimates was assessed using six normal and three failing canine hearts. The mean absolute difference between inclination angles of acquired and estimated fiber orientations was 15.4°. Computational simulations of ventricular activation maps and pseudo-ECGs in sinus rhythm and ventricular tachycardia indicated that there are no significant differences between estimated and acquired fiber orientations at a clinically observable level.
The iterative interaction between experimentation and simulation has long played a central role in the advancement of biological sciences. Among computational models of various biological systems, the heart is the most highly advanced example of a virtual organ, capable of integrating data at multiple scales, from genes all the way through the whole organ . State-of-the-art whole-heart models of electrophysiology and electromechanics are currently being used to study a wide range of phenomena, such as normal ventricular propagation, arrhythmia, defibrillation, electromechanical coupling, and cardiac resynchronization .
Recent years have witnessed revolutionary advances in imaging, including ex vivo structural and diffusion tensor (DT) magnetic resonance imaging (MRI), that facilitate acquisition of the intact structure of explanted hearts with hitherto unavailable resolution. Leveraging on these advances, a new generation of whole-heart models with unprecedented detail have emerged –. Such ex vivo models are being used in basic research to uncover specific mechanisms of heart dysfunction in diseases such as myocardial infarction and heart failure.
Translation of this ex vivo approach to the in vivo setting of the clinic can yield far reaching consequences. Ultimately, doctors will be able to use simulation tools to predict outcomes of various electrophysiological interventions, and select the optimal cardiac therapy specific to a patient. To this end however, patient-specific ventricular fiber orientation data are necessary, as these orientations determine directions of electrical propagation and strain distributions in the heart. While ventricular geometry and structural remodeling data can now be routinely acquired in vivo with imaging techniques such as late gadolinium enhanced cardiac MR , physicians do not currently have access to a method to acquire myocardial orientations data in vivo.
Myocardial fiber orientation data necessary for computational studies of cardiac electrophysiology and mechanics have been predominantly acquired using histological techniques; examples include the canine  and swine  ventricular models developed at Auckland University, and the rabbit ventricular model developed at the University of California, San Diego . While histological sectioning has traditionally been the gold standard for fiber data acquisition, this technique is exceedingly time consuming, prone to tissue deformation errors introduced due to cutting, restricted to 2-D measurements, and not applicable in vivo. Recently, DTMRI has emerged as a powerful alternative to histology, and DTMRI cardiac fiber data have been acquired in many species , , . DTMRI-acquired fiber orientations have been shown to match histological data , and have begun to be used in whole-heart computational electrophysiology and electromechanics , –, , . Unlike histology, DTMRI offers faster data acquisition, 3-D data, and nondestructive imaging. However, there are several drawbacks to the DTMRI technology, most importantly those that limit its application in the in vivo setting, where the bulk motion of the heart interferes with diffusion measurements.
Attempts to translate the DTMRI technology to the in vivo setting were first made by Edelman et al. , but it was later found that their diffusion encoding interacted with myocardial strain, and therefore, a method to correct for strain effects was proposed , . But this correction was achieved at the cost of increased acquisition noise and complexity. Although the strain correction method in combination with interpolation has recently been used to perform 3-D tractography of ventricular fiber orientations , the data prior to interpolation was very sparse, scan time was about 40 min, and breathholding, a difficult task for patients which typically causes inter-slice misalignments, was necessary. Tseng et al.  and Dou et al.  have proposed techniques that do not require strain correction, but the former is applicable to only those phases of the cardiac cycle called sweet spots, where the strain effect is negligible, whereas the signal-to-noise ratio of the latter is limited by the gradient intensity available in clinical scanners. Recently, a spin-echo technique to measure cardiac diffusion in vivo was proposed, and the utility of this method in measuring 2-D in-plane diffusion along short axis slices demonstrated . Toussaint et al.  have used this method in combination with curvilinear interpolation to perform 3-D tractography of ventricular fiber orientations, but their acquisition was limited to the equatorial part of LV, scan time was about 1.5 h, and breath-hold correction, a very challenging postprocessing step, was necessary. All in all, despite the above mentioned efforts over almost two decades, in vivo translation of cardiac DTMRI remains challenging due to technological inadequacies, e.g., in gradient strength, whole-heart imaging, radio-frequency, and navigators .
Previously, there have been a number of efforts to estimate fiber orientations in the ventricles. Among these, rule-based methods have gained some popularity; they assign fiber orientations to the ventricular walls such that the inclination angles vary smoothly in the transmural direction, as a mathematical function of distance from the endocardium , , , . Though there have been two preliminary studies comparing fiber orientations derived using rule-based methods with those obtained from DTMRI , , how well the idealized rules of myocyte orientation quantitatively correlate with true anatomy in normal and diseased hearts remains largely unknown. A different category of estimation methods utilizes geometric correspondences between two hearts. Panfilov’s group assigned canine heart fiber orientations to the ventricular surfaces of a human heart by computing a landmark-based surface mapping between the human heart and a canine heart . The intramural fiber orientations of the human heart, however, were derived from the acquired surface directions and a rule-based scheme. Sundar et al.  presented a technique to estimate fiber orientations via an elastic image transformation of an atlas heart. They tested the methodology using normal and failing canine hearts, and showed that estimated fiber orientations correlate with DTMRI data better than fiber orientations synthesized using a rule-based method. But the elastic mapping employed by them is not guaranteed to be diffeomorphic, and therefore may not preserve the topology of anatomical structures. More importantly, Sundar et al. did not measure the effect of error in their estimated fiber orientations on the results of simulations of cardiac electrophysiology, or any other functions of heart.
We seek to address the need for a new methodology to acquire myocardial fiber orientations for patient-specific simulations of cardiac function. Importantly, we hypothesize that ventricular fiber orientations of a normal or a failing heart can be estimated, given the geometry of the heart and an atlas, where the atlas is a heart whose geometry and fiber orientations are already available. Should our hypothesis be proved, it would pave the way for the use of estimated fiber orientations in patient-specific models of cardiac function.
In this study, we use state of the art techniques to develop the methodology for estimation of cardiac fiber orientations in vivo, and test the hypothesis in normal and failing ventricles. We quantify the estimation error in each heart, and then assess the effect of this error on simulations of cardiac electrophysiology in the specific heart, by computationally simulating local electrical activation maps as well as pseudo-electrocardiograms (pseudo-ECGs).
The central idea of our fiber estimation methodology is to exploit similarities in fiber orientations, relative to geometry, between different hearts in order to approximate the fiber structure of a (target) heart for which only the geometry information is available. Fig. 1 illustrates the processing pipeline of our methodology. In the following subsections, we describe the various components of the pipeline by demonstrating how the estimation was performed for an example patient who was scanned using in vivo computed tomography (CT). The numbers inside some of the blocks in Fig. 1 refer to the corresponding subsections below.
Geometry and fiber orientations of the atlas were extracted from raw data acquired in a previous study  using high resolution (0.4297 × 0.4297 × 1 mm3) ex vivo structural MRI and DTMRI, respectively, of a normal human heart in diastole. From the acquired structural image, only the ventricular myocardium of the atlas was extracted by fitting, for each short-axis slice, closed splines through a set of landmark points placed semi-automatically along the epicardial and endocardial boundaries in the slice [Fig. 2(A)]. The placement of landmark points was performed manually for a number of slices that were evenly distributed, with an inter-slice spacing of about 10 mm in the image. The landmark points for the remaining slices were obtained automatically by linearly interpolating the manually identified points. A generic mathematical formulation of this interpolation technique has been reported elsewhere , where the technique was used for identifying the atrio-ventricular boundary in an ex-vivo structural MR heart image. Fig. 2 (B) and (C) shows the 3-D geometry and fiber orientations of the reconstructed atlas ventricles; the fibers forming a left-handed helix on the epicardium are clearly visible.
Ventricular geometry of the patient heart in diastole was reconstructed from an in vivo CT image using a segmentation method that is similar to the one used for the atlas. The patient image was re-sampled prior to reconstruction such that the in-plane resolution was 0.4297 × 0.4297 mm3 as in the atlas. Similarly, the number of slices for which landmarks were manually picked, and the interval of out-of-plane interpolation was adjusted so that the segmented patient heart image had a slice thickness of 1 mm. Fig. 3(A) and (B) illustrates the patient ventricular geometry reconstruction.
Following the reconstruction of the patient ventricles, the atlas ventricular image was deformed using the tools of computational anatomy  to match the patient geometry image. The deformation was performed in two steps. The first was an affine transformation based on a set of thirteen landmarks points: the left ventricular (LV) apex, the two right ventricular (RV) insertion points at the base, the two RV insertion points midway between base and apex, and four sets of two points that evenly divide RV and LV epicardial contours at base, and midway between base and apex. Fig. 3 (C) and (D) illustrates the affine transformation of the atlas to match the patient ventricles. In the second step of atlas deformation, the affine-transformed atlas ventricles were further deformed to match the patient geometry, using large deformation diffeomorphic metric mapping (LDDMM).
LDDMM is an image registration algorithm that computes diffeomorphic (invertible, smooth, and with a smooth inverse) transformations between images . Given an atlas image Ia : Ω → R and a patient image Ip : Ω → R, LDDMM computes a flow of diffeomorphisms to transform Ia to match Ip, where Ω R3 is the 3-D cube in which the image data are defined, t [0,1], and υ is a smooth, compactly supported, time-dependent velocity vector field such that
The initial diffeomorphism ϕ0 is the identity transformation. The final diffeomorphism = ϕ1 is calculated by integrating the optimal vector field given by
The norm ‖ · ‖V is defined as ‖f‖V = ‖Lf‖L2 to enforce smoothness on the vector fields υ V, where L = (−α2 + γ)2 I3×3 is a differential operator of the Cauchy–Navier type. The parameters α and γ control the elasticity of the transformation, and ‖ · ‖ L2 is the L2 norm of square integrable functions defined on Ω. The optimal solution is computed by means of a gradient descent search, with the gradient
where and |Df| is the determinant of the Jacobian matrix. The diffeomorphic property of LDDMM guarantees that the atlas does not “fold over” itself during deformation, thereby preserving the integrity of anatomical structures. The above equations also ensure that the computed transformations ϕt define a geodesic path in the space of diffeomorphisms, thereby providing a metric quantification of the difference between the atlas and patient. While geodicity is not essential for fiber orientations estimation, from the standpoint of our long-term goals in patient-specific cardiac image analysis, it is a desirable property. For detailed mathematical descriptions of the LDDMM algorithm and implementation, the reader is referred to previous publications , . Fig. 3(E) illustrates the LDDMM transformation of the atlas to match the patient geometry; the match is remarkable. Note that the registration was performed on 3-D binary images, and not surfaces.
In the final phase of the processing pipeline, the transformation matrix of the affine matching and the deformation field of the LDDMM transformation were applied in sequence to morph the DTMRI image of the atlas to obtain an estimate of the patient ventricular fiber orientations. The morphing of the atlas DTMRI image consisted of spatial repositioning of image voxels in accordance with the spatial transformation of the geometry images, and reorientation of the DTs. For the reorientation of the DTs, we considered two methods, namely the preservation of principal directions (PPD) and finite strain (FS) . The FS method models transformations as rigid rotations, and reorients DTs accordingly. On the other hand, the PPD method models the deformation field at an image voxel as an affine transformation, and reorients the DT at that voxel such that the primary eigenvector of the new DT is , where e1 is the primary eigenvector of the original DT and is the affine transformation . The secondary eigenvector of the new DT is orthogonal to e1 and lies in the plane spanned by and Fe2, where e2 is the secondary eigenvector of the original DT. Thus, the PPD method incorporates, in addition to rotations, deformations such as shearing and nonuniform scaling. This method preserves not only the principal DT direction, but also the plane spanned by the two largest eigenvectors. Therefore, in this study, we employed the PDD method. We also performed some quantitative comparisons of the two methods to ensure that the PPD method performs better.
The estimate of the patient fiber orientations was obtained from the morphed atlas DTMRI image by computing the primary eigenvector of the DTs . Fig. 3(F) illustrates a streamlined visualization of the estimated fiber orientations in the patient ventricles. The fiber orientations vary smoothly across the myocardial wall, form a left-handed helix on the epicardium, and appear as a right-handed helix on the endocardium .
Tests of efficacy of the proposed methodology were performed on canine hearts due to the unavailability of human hearts. Ventricles segmented from a total of six normal and three failing canine hearts, the fiber orientations of all of which were acquired in diastole with ex vivo DTMRI at a resolution of 312.5 × 312.5 × 800 µm3 were used. The datasets have been employed in previous studies , , , , where the reader can find details of the acquisition. In brief, the experimental model of cardiac dyssynchrony and heart failure was generated in canines via radio-frequency ablation of the left bundle-branch block followed by three weeks of tachypacing (210 min−1). The animals were then euthanized and hearts explanted. Images of the explanted hearts were acquired and the ventricles segmented similarly to the human atlas heart, as described under Section III-A1. In the following, ventricles segmented from normal canine hearts are referred to as hearts 1–6, and those segmented from failing canine hearts as hearts 7–9.
To determine if the fiber estimation is accurate in the structurally normal ventricles, five different estimates of ventricular fiber orientations of heart 1 were obtained by using each of hearts 2–6 as an atlas. To test the accuracy of ourmethodology in the structurally remodeled ventricles, fiber orientations for each of the failing ventricles were estimated using heart 1 as the atlas. The purpose of using multiple atlases and a single target in the normal case was twofold. Firstly, it allowed us to confirm that the performance of the proposed methodology was independent of the choice of the atlas. Secondly, it made averaging of ventricular distributions of error across different estimates easier.
Estimation error was computed by comparing inclination angles of estimated and acquired fiber orientations. The effects of this error on outcomes of computational electrophysiological simulations, specifically activation times and ECGs, were examined. Since our estimation methodology was developed for the purpose of modeling, the performance analysis was focused on simulation outcomes.
Fiber orientations are typically characterized using inclination angles  following the tradition of histology, where angular measurements are performed on tissue sections that are cut parallel to the epicardial surface. Inclination angle of a fiber direction is the angle between the circumferencial direction and the projection of the fiber direction on the epicardial tangent plane. Since the angle between the fiber direction and epicardial tangent plane is generally small , , the information loss in describing a fiber direction entirely using its inclination angle is insignificant. We used the technique presented in Scollan et al.  to calculate inclination angles directly from image data. Briefly, each fiber direction was projected onto the tangential plane at the nearest epicardial point, and the angle between the projection and the circumferential direction measured. The tangential plane at an epicardial point was computed as the plane perpendicular to the short-axis imaging plane, and tangential to the short-axis epicardial contour passing through that point. The circumferential direction was calculated as the vector cross product between the normal to the short-axis imaging plane that points toward the base of the heart, and the normal to the tangential plane that points away from the heart. At each voxel that belongs to the myocardium of each target heart, we computed the estimation error as |θe − θa|, where θe and θa are the inclination angles of estimated and acquired fiber orientations, respectively. Note that inclination angles have values between −9° and +90°, and therefore, the estimation error ranges between 0° and 180°. To demonstrate that the inclination angles characterize the fiber orientations well, we also computed the acute angle between estimated and acquired fiber directions in 3-D by means of the vector dot product.
We quantified the effect of estimation error on the outcome of computational cardiac electrophysiology simulations by comparing ventricular activation maps and pseudo-ECGs calculated using estimated fiber orientations with those using acquired fiber orientations.
From heart 1, six models, one with the DTMRI-acquired fiber orientations of heart 1 (referred to as model 1), and five with the five estimated fiber orientations datasets (models 2–6), were constructed. For each of the three failing heart geometries, two ventricular models, one with the DTMRI-acquired fiber orientations and the other with the estimated fiber orientations were also constructed. The heart failure models with DTMRI-acquired fibers were denoted as models 7–9, and those with estimated fibers as models 10–12. Computational meshes were generated for each model. Our methodology for mesh generation, given segmented geometry and fiber orientations, involved discretization of the computational domain into finite elements using the software Tarantula,1 and assignment of fiber orientations to the elements using interpolation. Details of the mesh generation procedure have been reported in a previous publication .
Mathematical description of cardiac tissue was based on the monodomain representation, which is a reaction-diffusion system composed of a partial differential equation coupled to system of ordinary differential equations . The governing equations are
where σi is the intracellular conductivity tensor; Vm is the transmembrane potential; Cm is the membrane specific capacitance; and Iion is the density of the transmembrane current, which in turn depends on Vm and a set of state variables μ describing the dynamics of ionic fluxes across the membrane. For Cm, we have used a value of 1 µF/cm2. For σi in normal canine heart models, we used longitudinal and transverse conductivity values of 0.34 S/m and 0.06 S/m, respectively, as described by Roberts et al. . Iion was represented by the Greenstein et al.  ionic model of the canine ventricular myocyte. This model consisted of canine specific representations of membrane ion channels and intracellular calcium cycling. For the failing canine heart, we used the modifications in canine myocyte membrane dynamics as described by Winslow et al. . Electrical conductivities in canine heart failure ventricular models were also decreased by 30%, consistent with experimental data on electrical propagation in failing hearts . For more details on our cardiac tissue representation including the various parameters, the reader is referred to the relevant publications –.
Sinus rhythm was simulated by replicating activation originating from the Purkinje network, The ventricles were activated at six locations on the endocardium: one on the RV free wall, three on the LV septum and two on the LV free wall, as described in . Appropriate timings of the stimuli were chosen such that the resultant 3-D electrical propagation matched with experimental data , . Since failing hearts are arrhythmogenic, in addition to sinus rhythm, reentrant ventricular tachycardia (VT) was induced in the six failing models to assess the utility of the fiber estimation methodology in models of cardiac arrhythmias. An S1–S2 pacing protocol  was used to induce VT. The timing between S1 and S2 was chosen such that VT activity was sustained for 2 s after S2 delivery. If VT was not induced for any S1–S2 timing, the conductivities were further decreased by up to 70% so that VT was induced, as done in experimental protocols that utilize pharmacological intervention to decrease conduction velocity and thus increase the likelihood of inducing VT .
To assess the effects of estimation error on clinically observable electrophysiological data, pseudo-ECGs were calculated by taking the difference of extracellular potentials between two points in an isotropic bath surrounding the hearts. The two points were placed near the base of the heart separated by 18 cm, and the line connecting themwas perpendicular to the base-apex axis. The extracellular potentials at these points were approximated using the integral equation introduced by Gima et al. . To quantify difference between two ECG waveforms, the so-called mean absolute deviation (MAD) metric was computed as
where X and Y are waveforms of length n . MAD varies between 0% and 100%, corresponding to identical and completely different waveforms, respectively. The MAD metric was suitable because it has been utilized in clinical studies, to compare ECGs of reentrant activity and paced propagation for localization of the organizing centers of reentrant circuits . A MAD score of less than 12% implies that the two underlying propagation patterns are clinically equivalent. All simulations were performed using the software package CARP (CardioSolv, LLC); the numerical methods are described in previous publications , .
Fig. 4(A)–(C) displays streamlined visualizations of estimated as well as DTMRI-derived fiber orientations in normal and failing hearts. Qualitative examination shows that estimated fiber orientations align well with DTMRI-derived ones. Fig. 4(D) illustrates, overlaid on the geometry of heart 1, the distribution of error in normal hearts’ inclination angles, averaged across all five estimates. Fig. 4(E) shows the mean distribution of error in failing hearts’ inclination angles, overlaid on the geometry of heart 1. Fig. 4(F) and (G) presents sections of tissue from the distributions in Fig. 4(D) and (E), respectively. In Fig. 4(E), the errors calculated in the ventricles of each of the three failing hearts were mapped onto the geometry of heart 1 based on the point-to-point correspondences obtained from affine and LDDMM transformations; the mapped errors were then averaged. Fig. 4(F) and (G) highlights the transmural variation of error. The histograms of errors in Fig. 4(H) suggest that most myocardial voxels have small error values. About 80% and 75% of the voxels have errors less than 20 in normal and failing ventricles, respectively. It was found that the mean error, averaged across all estimated datasets, and all image voxels that belonged to the myocardium, were 14.4° and 16.9° in normal and failing ventricles, respectively. The mean error in the entire myocardium, in normal and failing cases combined, was 15.4°. The mean 3-D acute angle between estimated and acquired fiber directions were 17.5° and 18.8° in normal and failing ventricles, respectively. The 3-D angles are comparable to the estimation errors, demonstrating that the information loss in describing fiber orientations by means of inclination angles is insignificant. The standard deviation of error across the five different estimates of the fiber orientations of heart 1 was only 1.9° indicating that the variation in estimation quality from one atlas to another was small.
To compare the performance of the PDD and FS reorientation strategies in our methodology, the processing pipeline was modified by replacing the PPD method with the FS technique in the affine transformation step. The modified pipeline was then applied to estimate fiber orientations in failing ventricles, and errors were recalculated. The new mean error in the failing ventricles was 18°, which is higher than the error that corresponds to the PPD method. This result indicates that the PPD method outperforms the FS method in predicting fiber orientations in failing ventricles.
In the following, we describe our findings regarding the effect of estimation error on the results of electrophysiological simulations. Figs. 5 and and66 present the simulated activation maps of one beat of sinus rhythm activation in normal and failing ventricular models, respectively. In both normal and failing ventricles, the activation maps in models with acquired fiber orientations match well experimental activation maps , with the earliest epicardial activation at the anterior base, anterior apex, and LV free wall. Models with estimated fiber orientations produce activation maps very similar to those of models with acquired orientations; the earliest epicardial activations occur at the same sites, and the directions of propagation match as well. Among the normal ventricles, the total activation time of model 1 was 151 ms, while for models 2–6, the mean total activation time was 154 ms. The overall mean difference in total activation times between the acquired and estimated fiber orientation cases in the normal ventricular models was 5.7 ms, which is a small fraction (3.7%) of the total activation time. Fig. 5(B) indicates that regions with largest differences in local activation times in the normal ventricles are near the base, especially near the RV outflow tract, where the wall is very thin. We believe that this wall thinness caused imperfect registration and large estimation errors, which in turn led to large differences in local activation times. In models 5 and 6 in Fig. 5(D), the large estimation errors have caused significant conduction slowdown near the RV outflow tract. These errors, however, remained local, and did not have a significant bearing on the global differences in activation times. Fig. 5(C) demonstrates that pseudo-ECGs obtained for sinus rhythm simulations with models 1 and 3 have identical morphologies. The MAD score between these two waveforms was 4.14%. On average, the MAD score between sinus rhythm pseudo-ECGs with each of models 2–6 and model 1 was 10.9%.
In simulations of sinus rhythm with failing ventricular models, total activation times for models 7–9 were 183, 144, and 166 ms, respectively, while for models 10–12, they were 187, 139, and 162 ms, respectively. The total activation times in the heart failure models are longer than those in normal ventricular models due to the lower tissue conductivity values. The failing ventricles were also on average 29% larger than normal ventricles resulting in longer paths for wavefront propagation and thus longer total activation times. The mean difference in total activation times between heart failure models with acquired and estimated fiber orientations was only 5.2 ms (3.1%), while the mean MAD score was 4.68%. These results indicate that the outcomes of simulation of ventricular activation in sinus rhythm in normal and failing canine ventricular models with fiber orientations estimated with the present methodology closely match those with acquired orientations. In particular, presence of heart failure did not diminish the accuracy of the estimation. Incidentally, we did not average the three ventricular distributions of difference in activation maps displayed in column 3 of Fig. 6 because these distributions are not similar.
Fig. 7 shows simulated activation maps, in apical views of the ventricles, during one cycle of induced VT in the heart failure models, and corresponding pseudo-ECGs. Simulations with acquired and estimated fiber orientation both exhibit similar figure-of-eight reentrant patterns. The VT cycle lengths for models 7–12 were 392, 393, 436, 391, 408, and 432 ms, respectively. Accordingly, the percentage difference in VT cycle length between models with acquired and estimated fiber orientations, averaged across the three failing hearts, was 1.7%. In hearts 7 and 9, the ECG morphologies corresponding to estimated and acquired fiber orientations were in excellent agreement. In heart 8, the longer VT cycle length in the model with estimated fibers was manifested as a minor phase shift in the ECG, but the overall morphology of the two ECG recordings was still the same. The mean MAD score was 9.3%. These results indicate that canine heart failure models with estimated fiber orientations can closely replicate outcomes of VT simulations performed using acquired fiber orientations.
The goal of this project was to develop a methodology for estimating ventricular fiber orientations from in vivo images of ventricular geometry, and to assess the accuracy of this estimation and its utility in simulations of cardiac electrophysiology. We hypothesized that the fiber orientations can be estimated in normal ventricles and those with structural remodeling, such as heart failure. The results described in the previous section support our hypothesis. They show that inclination angles of predicted fiber orientations are comparable to those acquired by ex vivo DTMRI, the state-of-the-art technique , . Furthermore, activation maps and pseudo-ECGs generated by electrophysiological simulations using estimated fiber orientations closely match those using acquired fiber orientations. This research demonstrates quantitatively that, in the absence of DTMRI, myocardial fiber orientations of normal and failing ventricles can be estimated from in vivo images of their geometries for use in simulations of cardiac electrophysiology. The proposed methodology is demonstrated with in vivo CT data, but it is equally applicable to in vivo MR images of ventricular geometry, addressing the lack of ability to directly acquire patient fiber orientations. It is thus an important step towards the development of personalized models of ventricular electrophysiology for clinical applications. The methodology can also be used to estimate the fiber orientations in ex vivo hearts with high resolution. This is especially helpful when acquiring submillimeter resolution DTMRI images is difficult or prohibitively expensive, due to the very long acquisition times, in the order of days .
Our results are in agreement with research conducted by Helm et al. , who showed that heart failure does not significantly alter fiber orientations relative to geometry. Furthermore, our results demonstrate that the largest errors in estimation occur near the epicardial and endocardial surfaces, and the RV insertion points. The decreased accuracy near the surfaces is due to the partial volume effect in DTMRI  which introduces error and imperfect LDDMM registration, while the larger errors near RV insertion points are due to the complex branching of fiber bundles and the difference in fiber orientations between LV and RV.
For the reorientation of DTs, we chose the PPD strategy, as it accounts for shearing and nonuniform scaling, and therefore is more generic than the FS strategy . On the other hand, Peryrat et al.  have preferred the FS method because it is consistent with the Log-Euclidian framework that they used to compute statistics on DTs, and is not affected by errors in eigenvector extraction. They also argued that FS method is better at preserving transmural variation of fiber structure, and their quantitative experiments comparing all three eigenvectors of DTs indicated that the FS strategy registers images better than PPD. However, FS outperformed PPD only by a very small margin in their work, and their experiments were limited to normal canine hearts. Form the standpoint of patient-specific modeling of cardiac function, it was important for us to adopt a reorientation strategy that performs well in the presence of structural remodeling such as in heart failure, which involves shape changes that PPD may better account for. Furthermore, only the primary eigenvector was relevant in this study. Indeed, when we replaced PPD with FS as part of our pipeline, the fiber estimation error for failing hearts increased. A more comprehensive comparison of various reorientation strategies may be helpful in the future, but such an analysis is outside the scope of this research.
Our electrophysiological simulations suggested that activation maps were not very sensitive to changes in fiber orientations. More importantly, we demonstrated that effects of fiber estimation errors on gross electrophysiology were insignificant at a clinically observable level by means of the MAD score of the pseudo-ECGs. The MAD scores in our results were less than 12%, implying that the propagation patterns generated using estimated and acquired fiber orientations were clinically equivalent . Note that the similarity of propagation patterns will translate to low differences in mechanical activation patterns as well, as experiments show that local electrical and mechanical activation times during sinus rhythm are highly correlated , . In summary, our research will facilitate ventricular simulation studies of any species in health and disease when it is not feasible to acquire fiber orientations using DTMRI. In particular, the proposed methodology paves the way for patient-specific modeling of ventricular whole-heart electrophysiology (and possibly) electromechanics based only on in vivo clinical imaging data. Simulations with such models may ultimately aid physicians to arrive at highly personalized decisions for therapeutic interventions as well as prophylaxis. For example, it has recently been shown that such models can be used to predict the organizing centers of reentrant activity in infarcted hearts, which could aid in guiding ablation therapy for infract-related VT . Also, electromechanical simulations using models of failing hearts have shown promise in noninvasively determining optimal cardiac resynchronization therapy .
The current study has a number of limitations. Firstly, human heart image data were not available to us, and therefore the proposed estimation methodology was validated with images of canine hearts. We expect that the methodology will accurately estimate fiber orientations in human hearts as well because, just as in canine hearts, fiber orientations relative to geometry have been shown to be similar between different human hearts . Also, we tested our methodology in normal and failing hearts only. It would be important to test it under conditions such as myocardial infarction and hypertrophy, where fiber disorganizations are known to occur , . Furthermore, the present study is a proof-of-concept for estimating ventricular fiber orientations for use in personalized simulations of cardiac electrophysiology. Patient-specific acquisition of other electrophysilogically relevant data such as action potential heterogeneity will also be necessary for the development of personalized models of cardiac electrophysiology; our methodology is only one step towards achieving this goal. Finally, the accuracy of the fiber orientation estimation in simulations of cardiac mechanics will also need to be assessed.
The authors would like to thank Dr. R. Winslow, Dr. E. McVeigh, and Dr. P Helm at Johns Hopkins University for providing the ex vivo datasets online.
This work was supported in part by the National Institutes of Health under Grant R01-HL082729 and Grant R01-HL091036, and in part by the National Science Foundation under Grant CBET-0933029.
Fijoy Vadakkumpadan, Institute for Computational Medicine and the Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD 21205 USA.
Hermenegild Arevalo, Institute for Computational Medicine and the Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD 21205 USA.
Can Ceritoglu, Center for Imaging Science, Johns Hopkins University, Baltimore, MD 21218 USA.
Michael Miller, Center for Imaging Science, Institute for Computational Medicine, Baltimore, MD 21218 USA, and also with the Department of Biomedical Engineering, Johns Hopkins University, MD 21205 USA.
Natalia Trayanova, Institute for Computational Medicine and the Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD 21205 USA.