Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Ann Biomed Eng. Author manuscript; available in PMC 2012 December 11.
Published in final edited form as:
PMCID: PMC3518842

A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models


Electrical waves traveling throughout the myocardium elicit muscle contractions responsible for pumping blood throughout the body. The shape and direction of these waves depend on the spatial arrangement of ventricular myocytes, termed fiber orientation. In computational studies simulating electrical wave propagation or mechanical contraction in the heart, accurately representing fiber orientation is critical so that model predictions corroborate with experimental data. Typically, fiber orientation is assigned to heart models based on Diffusion Tensor Imaging (DTI) data, yet few alternative methodologies exist if DTI data is noisy or absent. Here we present a novel Laplace–Dirichlet Rule-Based (LDRB) algorithm to perform this task with speed, precision, and high usability. We demonstrate the application of the LDRB algorithm in an image-based computational model of the canine ventricles. Simulations of electrical activation in this model are compared to those in the same geometrical model but with DTI-derived fiber orientation. The results demonstrate that activation patterns from simulations with LDRB and DTI-derived fiber orientations are nearly indistinguishable, with relative differences ≤6%, absolute mean differences in activation times ≤3.15 ms, and positive correlations ≥0.99. These results convincingly show that the LDRB algorithm is a robust alternative to DTI for assigning fiber orientation to computational heart models.

Keywords: Diffusion tensor, Magnetic resonance, Imaging, Finite element, Mesh, Laplace–Dirichlet


Electrical waves traveling throughout the myocardium trigger the muscle contractions necessary to pump blood throughout the body. The shape and direction of these waves are functions of the spatial arrangement of myocytes within the myocardium,28 which is termed fiber orientation since it macroscopically resembles a network of fibers within a laminar sheet architecture.32,22 Furthermore, ventricular mechanics also depends on fiber orientation.7 Thus, in computational studies that simulate electrical wave propagation or mechanical contraction in the myocardium,33 fiber orientation should be represented accurately so that model predictions corroborate with experimental data.18,27,35,36

Common approaches to perform this task include image-based and rule-based methods. The most utilized image-based method is Diffusion Tensor Imaging (DTI) since it provides non-destructive three-dimensional information on myocardial fiber orientation. Specifically, DTI provides three orthotropic axes for every voxel in a data set pertaining to an imaged heart,19,30 with the axes being in the longitudinal (along the fibers), transverse (along laminar sheets and perpendicular to the fibers) and sheet normal (orthogonal to the longitudinal and transverse directions) directions. These axes can be directly interpolated into a computational heart model constructed from geometrical data acquired on the same heart with magnetic resonance imaging (MRI).15,34 However, in the presence of noise due to partial volume effects1 on cardiac surfaces and in anatomical structures with a large surface-area-to-volume ratio,29 fiber orientation axes are difficult to interpolate accurately. Furthermore, DTI is performed ex vivo. Therefore, in models constructed from geometrical data acquired in vivo using computed tomography or MRI, which are not able to image fiber orientation well, fiber orientation needs to be approximated or mapped from an alternate source.

For these reasons, rule-based algorithms6,25 serve as an alternative method to DTI for assigning myocardial fiber orientation to computational heart models. These algorithms are termed “rule-based” since they generate mathematical descriptions for fiber orientation with rules formulated on the basis of observations from histology10,11,21,22,32 and DTI.30,31 For all rule-based algorithms, the transmural and apicobasal directions throughout the entire myocardium of a model need to be parametrized in order to systematically assign orthotropic fiber orientation.

Traditionally, rule-based algorithms parameterize the transmural direction based on the minimal distance between the endocardial and epicardial surfaces, and define the apicobasal direction parallel to the long axis of the heart.4 This approach has shown some success in ventricular models with basic geometries.6,25 However, minimal distance parameterizations do not guarantee the absence of singularities in the minimal distance function throughout the entire myocardium,3 particularly in the septum, and in endocardial structures such as papillary muscles. These singularities can yield unpredictable results for the fiber orientation, and the problem can persist even when weighted distance mapping is performed to minimize their occurrence. For example, when using minimal distance to determine transmural depth from the endocardium to epicardium, the normalized distance is zero on the left ventricle (LV) and right ventricle (RV) surfaces, 1.0 on the epicardial surface, and 0.5 in the middle of the septum and at its junction with the LV and RV walls. The mean value theorem from calculus says that with this definition, there must be a point interior to the septum and at its junction with the LV and RV walls where the gradient is zero. Weighted averages will only exacerbate this problem, and redefining the minimal distance function in the septum to eliminate these singularities, and in a manner that is also continuous with the rest of the myocardium, is not trivial. Furthermore, minimal distance algorithms rely on a fixed apicobasal axis throughout the myocardium to assign fiber orientation, which can produce inaccuracies in the fiber orientation towards the apex due to portions of the LV and RV walls not being oriented along the long axis of the heart.

Thus, we have developed a novel Laplace–Dirichlet Rule-Based (LDRB) algorithm for assigning myocardial fiber orientation to computational heart models with speed, accuracy, and high usability. The first major innovation in this new algorithm is the use of a Laplace–Dirichlet method3,5 to simultaneously define the transmural and apicobasal directions for the entire myocardium. By using the continuity properties of the harmonic solutions to Laplace’s equation with Dirichlet boundary conditions, instead of minimal distance parameterizations, our new algorithm ensures that the transmural and apicobasal directions change smoothly and continuously throughout the entire myocardium, even in models with complex geometries. The second major innovation in the algorithm is the use of bi-directional spherical linear interpolation (bislerp) to interpolate fiber orientations within the myocardium. bislerp guarantees that fiber orientation will change smoothly and continuously throughout the entire myocardium, particularly in the septum and at its junctions with the LV and RV.

Here we demonstrate the utility of the LDRB algorithm in an anatomically accurate computational model of the canine ventricles. To validate this approach, simulations of electrical activation with the model of the canine ventricles containing LDRB fiber orientation are compared to simulation results with the same geometrical model but containing DTI-derived fiber orientation.


The following methods describe the approach for implementing and testing the LDRB algorithm. In the first part of the methods, the implementation of the LDRB algorithm in a generic computational model of the two-chambered mammalian ventricles is presented. The second part of the methods describes the approach for applying the LDRB algorithm to a MRI-based model of the canine ventricles, and for comparing the calculated electrical activation maps in the canine ventricular model with LDRB versus that with DTI-derived fiber orientation.

LDRB Algorithm

Pseudocode to the full LDRB algorithm for assigning fiber orientation to a model of the mammalian ventricles, consisting of a finite element mesh denoted Ω, is listed as Algorithm 1 in the online supplement. The inputs and functions of the algorithm are described in detail below.

LDRB Rules

The LDRB algorithm is based on the following properties of myocardial fiber orientation derived from histological and DTI data:

  • R1: The longitudinal fiber direction in the ventricular walls is parallel to the endocardial and epicardial surfaces.30
  • R2: The longitudinal fiber direction rotates clockwise throughout the ventricular walls from the endocardium (+α) to the epicardium (−α), where α is the helical angle with respect to the counterclockwise circumferential direction in the heart when looking from the base towards the apex.32
  • R3: The longitudinal fiber direction in papillary muscles and trabeculae is parallel to the long axis of these structures.11,21
  • R4: The transverse fiber direction is perpendicular to the longitudinal fiber direction and is defined by the angle β, where β is the angle with respect to the outward transmural axis of the heart.10 Please note, our definition of β is not equivalent to α3 in Streeter et al.32 For simplicity, β is assumed to vary transmurally from −β on the endocardium to +β on the epicardium.
  • R5: The sheet normal is orthonormal to the longitudinal and transverse fiber directions.22
  • R6: Fiber orientation in the septum is continuous with the ventricular walls.31

LDRB Inputs

To assign fiber orientation throughout the myocardium according to the rules R1–R6 above, the LDRB algorithm takes four functions as inputs, representing the desired α and β angles within the septum (s) and the ventricular walls (w). The angles α and β are in degrees and d is the transmural depth normalized from 0 to 1.





If non-linear interpolation functions are desired, these functions need to satisfy αw(0) = αs(0) = −αs(1) modulo 180° and βw(0) = βs(0) = −βs(1) modulo 180°. These conditions ensure that the coordinate systems used for assigning fiber orientation on the endocardium of Ω will be appropriately oriented, regardless of whether orientations are parametrized in the LV or RV wall.

The LDRB algorithm also requires the definition of the following surfaces in order to assign Dirichlet boundary conditions:

  • [partial differential]Ωapex: the apex of the ventricles
  • [partial differential]Ωbase: the base of the ventricles
  • [partial differential]Ωepi: the epicardial surface of the ventricles
  • [partial differential]Ωlv: the endocardial surface of the LV
  • [partial differential]Ωrv: the endocardial surface of the RV

For most studies, the surface [partial differential]Ωbase can be extracted by taking a cutting plane at the apicobasal junction, and the surface [partial differential]Ωapex can be extracted by finding the point lying closest to the ventricular apex of Ω. The surfaces [partial differential]Ωepi, [partial differential]Ωlv, and [partial differential]Ωrv can then be extracted by choosing an arbitrary point on the epicardium, LV endocardium, and RV endocardium, then iteratively expanding from these three points along the surface of Ω until intersection with [partial differential]Ωbase.

LDRB Functions

The LDRB algorithm uses four distinct functions in sequential order to compute the fiber orientation. Pseudocode of these are listed as Functions 1–4 in the online supplement, and each of them is explained below.

Function 1: Laplace–Dirichlet Parameterization of the Myocardium

Laplace(Ω, A, B): Returns the scalar values of the solution to Laplace’s equation on a mesh Ω with Dirichlet boundary conditions defined by surfaces A (maximum) and B (minimum).

Using Laplace, the LDRB algorithm solves one Laplace equation for a “potential energy” function ψab between [partial differential]Ωapex and [partial differential]Ωbase, and 3 more functions (ϕlv, ϕrv, ϕepi) between [partial differential]Ωepi, [partial differential]Ωlv, and [partial differential]Ωrv. The system of equations for ψ or ϕ, which are solved for by constructing the weak form of Laplace’s equation and applying the Galerkin weighted residual method outlined in Chapter 2 of Reddy,26 as well as equations 1–10 in Bayer et al.,3 are listed in Function 1 of the online supplement. We found that a matrix solver which employs the conjugate gradient method preconditioned with block Jacobi containing blocks ILU(0) and a stop condition relative error of 10−7 works well for computing ψ and ϕ.

By taking the gradient of the Laplace solutions, [nabla]ψab is used to define the apicobasal direction, and [nabla]ϕlv, [nabla]ϕrv, and [nabla]ϕepi are used to find the transmural direction. Notice that in Algorithm 1 and Function 1 in the online supplement, the boundary conditions chosen for ϕlv, ϕrv, and ϕepi allow the following conservation laws to be written:



Within the septum, ϕepi will be approximately zero, which implies that within the septum [nabla]ϕlv≈ −[nabla]ϕrv. Thus, in the septum either [nabla]ϕlv or [nabla]ϕrv can be used to obtain the transmural direction. The same argument applies for the pairs (ϕlv, ϕepi) and (ϕrv, ϕepi) within the LV and RV walls, respectively.

Function 2: Constructing the Fiber Orientation Coordinate System

axis(u, v): Given two vectors u (apicobasal) and v (transmural), returns a 3 × 3 orthogonal matrix representing the coordinate system for assigning fiber orientation within the myocardium.

The function axis takes as inputs, for a given point in Ω, a vector from [nabla]ψab along with a vector from either [nabla]ϕlv, [nabla]ϕrv or [nabla]ϕepi, then yields a right-handed axis system (Q = [ê0, ê1, ê2]) where ê0 is oriented in the circumferential direction, ê1 in the apicobasal direction, and ê2 in the transmural direction. The Q generated by axis is consistent with the coordinate system used in histology studies to quantify myocardial fiber orientation.2,14,22,32 Lastly, due to the manner in which the Laplace gradients are solved, [nabla]ψ and [nabla]ϕ are predominately perpendicular and aligned in the apicobasal and transmural directions, respectively. However, when [nabla]ψ and [nabla]ϕ are not perpendicular, as in endocardial structures (trabeculations or papillary muscles), the direction ê1 is preserved over ê2 in order to satisfy R3.

Function 3: Defining the Orthotropic Fiber Orientation

orient(Q, α, β): Takes a coordinate system Q and the fiber orientation angles α and β from Eqs. (1)(4) at a given point in Ω, then returns an orthonormal coordinate system (F  S  T), where F is the longitudinal direction, S is the sheet normal, and T is the transverse direction, as outlined in the Introduction. Refer to Fig. 1a for further detail.

The coordinate system used for assigning fiber orientation to computational models of the ventricles. (a) The longitudinal fiber (F) and transverse (T) directions with respect to the circumferential (ê0), apicobasal (ê1) and transmural ...

For a given point in Ω, the direction F is obtained by rotating Q around ê2 by the angle α at that point (Fig. 1b):

Rê2α=(ê0ê1ê2)(cos αsin α0sin αcos α0001)(ê0Tê1Tê2T)


The same logic is used to obtain the S and T directions. Given a point in Ω, S and T are generated by rotating the axis system (F êf ê2) from Eq. (7) around F by the angle β at that point (Fig. 1c):

RFβ=(Fêfê2)(1000cos βsin β0sin βcos β)(FêfTê2T)

(F S T)=RFβ(Fêfê2)

The combination of these equations into one operation can be found in Function 3 of the online supplement.

Function 4: Interpolating Fiber Orientations

bislerp(QA, QB, t): Linearly interpolates two orthogonal matrices QA and QB to produce a new orthogonal matrix determined by the interpolation factor t. When t is 0, bislerp returns QA, and when t is 1, bislerp returns QB.

The function bislerp is used by the LDRB algorithm to interpolate fiber orientation axes continuously throughout the myocardium. To satisfy R6, the LDRB algorithm uses bislerp to: (i) perform a weighted average between the axis systems Qlv and Qrv by using ϕlv and ϕrv to determine t, in order to obtain continuous fiber orientation between the RV and LV (Qendo), and then (ii) perform a second weighted average between the axis systems Qendo and Qepi by using ϕepi to determine t, in order to obtain continuous fiber orientation between Qendo and the rest of the ventricles. Note that in Algorithm 1 of the online supplement, when one of the energy functions goes to zero, the remaining two fiber axis systems become identical and the weighted averaging has no effect. Thus, weighted averaging occurs primarily at the junctions of the LV, RV, and septum to continuously blend fiber orientation according to R6. A more detailed description on the inner workings of bislerp can be found in the online supplement under Function 4.

Application of the LDRB Algorithm to a MRI-Based Model of the Canine Ventricles

Development of the Model of the Canine Ventricles

The LDRB algorithm was tested in a MRI-based computational model of the structurally normal canine ventricles. The geometry and DTI-derived fiber orientation of the model were constructed using the methods of Vadakkumpadan et al.34 applied to MRI and DTI data collected by Helm et al.16 To generate the same model but with LDRB fiber orientation, the input functions in Eqs. (1)(4) were optimized so that the mean angle between the LDRB and DTI-derived fiber orientations was minimal. The optimal values for the input parameters αendo, αepi, βendo and βepi of the LDRB algorithm were determined by varying each parameter in the range of ±90° in intervals of 5°, then choosing the parameter combination that produced the smallest mean angle (θmean) between the LDRB and DTI vectors (F  S  T) calculated over the total number of elements in Ω(Nelem). Since fiber orientation is bi-directional, θmean was calculated as



Membrane kinetics in this model of the canine ventricles were described by the Greenstein–Winslow myocyte model.12 Orthotropic tissue conductivities of 0.5 (S/m) along F, 0.3 (S/m) along T, and 0.16 (S/m) along S were assigned to produce conduction velocities within the range of 20–80 cm/s as observed in experiments.8 Monodomain simulations were performed with the model of the canine ventricles using the Cardiac Arrhythmia Research Package37 (CardioSolv LLC) running on 16 compute nodes, each with four Dual Core AMD Opteron processors (Model 2222) and 8GB of memory. All simulations were executed with a 10 µs time step.

Simulation Protocol and Activation Maps

Since the pattern of electrical activation in the myocardium, which underlies cardiac pump function, is dependent on fiber orientation, ventricular activation maps were generated and used to analyze the correspondence between the simulation results with LDRB and DTI-derived fiber orientations. Two sets of simulations were performed. The first set was performed by pacing the canine ventricles with LDRB and DTI-derived fiber orientations at the LV epicardium midway between the apex and base to elicit transmural electrical wave propagation, and the second set was performed by pacing at the LV endocardial apex to elicit apicobasal propagation. At each location, the myocardium was paced at a cycle length of 600 ms for 10 beats with stimuli of 2 ms duration and twice diastolic threshold amplitude. Activation maps were generated for the 10th beat by computing, at each location within the myocardium, the moment in time when the action potential upstroke velocity reached maximum. Differences in the activation maps between simulations with LDRB and DTI-derived fiber orientation were quantified using the methods of Han et al.,13 in which a relative difference (RD), root mean square difference (RMSD), and correlation coefficient (CC) were computed for each pacing protocol.


All vectors in (F  S  T) were normalized to have unitary magnitude, then visualized as streamlines using the software SpiderView ( Streamlines were colorized according to transmural depth, which was determined using the Laplace solutions for ϕepi, ϕlv and ϕrv. Please note, when streamlining fibers on the epicardial and endocardial surfaces in DTI data sets, streamlines in noisy regions will either appear missing, or be similar to the underlying fiber directions not containing partial volume effects. Activation maps from the simulation protocols described above were visualized using the software Meshalyzer (


LDRB Fiber Orientation in the Model of the Canine Ventricles

The conformal boundary-fitted mesh of the canine ventricles contained 823,077 vertices and 4,588,291 tetrahedra with a mean edge length of 600 µm. Within this mesh, the solutions for Laplace’s equation with Dirichlet boundary conditions applied to the ventricle surfaces (Fig. 2a) were both smooth and harmonic (Fig. 2b). Using the gradients of the Laplace solutions, a coordinate system consistent with histology studies (Fig. 1a) was automatically generated to accurately perform the α and β rotations (Figs. 1b and 1c) necessary to generate fiber orientation according to the rules R1–R6. The optimized LDRB parameters for these rotations in the canine ventricles were αendo = 40°, αepi = 250°, βendo = 265° and βepi = 25°, which are within the range reported by Helm et al.,15 and were used in the model of the canine ventricles with LDRB fiber orientation presented in Figs. 3, ,4,4, ,5,5, and and6.6. The computation time needed to calculate LDRB fiber orientation for any combination of α and β did not exceed 5 min on a single CPU (a maximum of 249 s to compute the Laplace functions, and 48 s to compute the orient, axis, and bislerp functions).

Solutions to the Laplace–Dirichlet scalar fields in the canine ventricles. (a) The Dirichlet boundary conditions (0 or 1) applied to the surfaces of the ventricles. Please note, the boundary condition [partial differential]Ωapex is a single point located ...
The LDRB and DTI-derived longitudinal fiber direction (F) in the model of the canine ventricles. (a) The streamlined LDRB and DTI-derived longitudinal fiber directions defined by the angle α. (b) Streamlines peeled away to visualize the internal ...
The LDRB and DTI-derived transverse direction (T) in the model of the canine ventricles. (a) The streamlined LDRB and DTI-derived transverse directions defined by the angle β. (b) Streamlines peeled away to visualize the internal transverse fiber ...
The LDRB and DTI-derived sheet normal direction (S) in the model of the canine ventricles. (a) The streamlined LDRB and DTI-derived sheet normal directions defined by the cross-product of the fiber directions T and F. (b) Streamlines peeled away to visualize ...
Simulation results from the model of the canine ventricles. (a) Activation maps obtained with the model of the canine ventricles with LDRB and DTI-derived fiber orientations during pacing at the LV epicardium and LV apex. Isochrone lines have a 10 ms ...

Panels (a) and (b) in Figs. 3, ,4,4, and and55 show the LDRB fiber orientation obtained for the canine ventricles next to the DTI-derived fiber orientation (refer to Section 2 of the online supplement for full-size images of the streamlined fiber orientations), while Figs. 3c, ,4c,4c, and and5c5c present the angles between the LDRB and DTI vector fields for the fiber orientation directions F, S, and T. As seen in Figs. 3a and 3b, the DTI-derived fiber orientation direction F displayed a profound transmural rotation parallel to the ventricular surfaces, which is consistent with rules R1 and R2, and was captured well by the LDRB algorithm. The largest differences between the LDRB and DTI-derived F (Fig. 3c) were localized near the endocardial and epicardial surfaces due to partial volume effects. As shown in Figs. 4a and 4b, the DTI-derived T displayed a predominant pattern of −β on the endocardium transitioning to +β on the epicardium, which is consistent with R4 and which the LDRB algorithm was able to capture. However, in contrast to θmean(F), the distribution of differences between the LDRB and DTI-derived T was more homogeneous (Fig. 4c), but this was expected since the DTI-derived T is susceptible to noise throughout the entire myocardium because its eigenvalues are similar in magnitude to S.24 Therefore, S in Fig. 5 displayed results similar to T since S is just the cross-product of T with F (R5). Lastly, the use of quaternions to interpolate fiber orientation proved successful for preserving R6 in all three fiber directions.

Simulations of Electrical Activation with the Model of the Canine Ventricles

Shown in Fig. 6a are the calculated electrical activation maps for pacing at the LV epicardium and apex in the model of the canine ventricles with LDRB and DTI-derived fiber orientations. Despite the differences in S and T between the LDRB and DTI-derived fiber orientations (Figs. 4c and and5c),5c), the LDRB model activation maps were nearly indistinguishable from the DTI model activation maps for both LV epicardial and apical pacing. The magnitude of the absolute difference in activation times between simulations in the model with the LDRB and DTI-derived fiber orientations did not exceed 11 ms (Fig. 6b). The activation time statistics at the top of Table 1 show that the LDRB and DTI activation maps were in excellent agreement, with small RD ≤ 6% and RMSD ≤ 3.15 ms, along with strong positive CC ≥ 0.99.

Activation time statistics for the model of the canine ventricles.


The goal of this study was to develop a novel rule-based algorithm for incorporating myocardial fiber orientation into computational heart models with speed, accuracy, and high usability. To achieve this task, we developed an algorithm termed LDRB, then demonstrated the algorithm’s capabilities by assigning fiber orientation to a realistic image-based computational model of the canine ventricles. The results presented in this study show that the LDRB algorithm is (i) fast, even for large meshes, (ii) accurate, by producing activation maps strikingly similar to those obtained from the model with DTI-derived fiber orientation, and (iii) user friendly, by providing a robust framework to easily navigate throughout the model’s geometry in order to automatically assign continuous fiber orientation. These findings provide strong evidence that the LDRB algorithm is a robust and inexpensive alternative to the use of DTI-derived fiber orientation in computational models of the heart.

LDRB Computation Time

In modeling pipelines designed for clinical use,34 the recent focus has been on fully automating and speeding up the model construction. We have shown that the LDRB algorithm is fast for large meshes. Specifically, the computation time needed to execute the LDRB algorithm in the image-based model of the canine ventricles, on a single CPU, was slightly less than 5 min. For comparison, when applying the minimal-distance rule-based (MDRB) algorithm by Potse et al.25 to generate fiber orientation (though only for F since the algorithm does not include rotations for T and S) in our model of the canine ventricles, the computation time was just shy of 13 min on the same CPU (See Section 4 in the online supplement for details). Thus, run times for the LDRB algorithm are just as fast or faster than existing rule-based methods, making it ideal for rapidly generating fiber orientation in geometrical models constructed from both in vivo and ex vivo MRI or CT scans.

LDRB Accuracy

If activation patterns are significantly altered by the underlying fiber orientation, both the dynamics of electrical activity and the sequence of contraction would be altered. Thus, if fiber orientation is unrealistic for a given model, it is unlikely the electrical and mechanical behaviors in simulations with this model would be physiologic. Therefore, we conducted tests to make sure that activation patterns generated by the LDRB algorithm compare well to those in models with the gold standard, DTI-derived fiber orientation. The LDRB algorithm compared extremely well, in that only small differences existed in the activation maps resulting from LV epicardium and LV apex pacing (small RD = 4–6% and RMSD = 1.48–3.15 ms, high CC ≥ 0.99).

To ensure these values are dependent on changes in activation patterns due to fiber orientation, and not the mesh resolution (which can potentially alter wavefront shape at larger discretizations9), the mesh of the canine ventricles was re-discretized to have a smaller mean edge length of 300 µm (6,344,174 vertices and 36,706,328 tetrahedra). The orthotropic tissue conductivities for this mesh were set to 0.36 (S/m) along F, 0.18 (S/m) along T, and 0.08 (S/m) along S to produce the conduction velocities observed in experiments.8 Then with this model, the same simulation protocol outlined in the methods was performed. As shown at the bottom of Table 1, activation time statistics hardly changed from the values determined at the 600 µm discretization (|ΔRD| ≤ 2%, |ΔRMSD| ≤ 0.91 ms, |ΔCC| ≤ 0.01). These results indicate that wavefront distortion due to the mesh resolution did not significantly factor into the activation time statistics for our study.

Lastly, to understand how activation time statistics change when the LDRB inputs α and β in Eqs. (1)(4) are varied, we performed additional simulations with the model of the canine ventricles (at the 600 µm discretization since computation times were more than 10× faster than with the 300 µm discretization) using a wide range of values for α and β. These results are included in Section 3 of the online supplement and show that activation patterns calculated from simulations with the model of the canine ventricles including unrealistic fiber orientation can produce RD as high as 54%, RMSD as large as 26 ms, and CC as low as 0.69.

LDRB Usability

The rule-based algorithms used currently for mammalian ventricles6,25 are based on minimal distance parameterizations and rely on nearest neighbor averaging of the minimal distance parameterizations to prevent discontinuities in fiber orientation, specifically in the middle of the septum, at the junctions of the LV and RV walls with the septum, and at the junctions of papillary muscles and trabeculae with the ventricular walls (See S-Figures 13 and 18 in the online supplement). Repairing anomalous fiber orientations, which could potentially alter electrical activation or mechanical deformation patterns, is not trivial and may require time-consuming manual intervention. In contrast, the LDRB’s use of bislerp guarantees that fiber orientation changes smoothly and continuously at these locations while avoiding the need to develop complicated conditions in order to tune the algorithm to individual models. This renders the LDRB algorithm virtually automatic after the definition of the input parameters (which can be directly taken from histology and DTI studies), thereby making it very user-friendly, and with results that are consistent for any model.

Alternatively, the use of a large deformation diffeomorphic metric mapping (LDDMM) approach34 to map ventricular fiber orientation into computational meshes of the ventricles could also potentially solve the problems encountered with minimal-distance based methods. However, execution times are very long (at the order of hours), and the LDDMM algorithm does not work well when the resolution of the target heart model is higher than that of the atlas template. For example, if the target model is of high resolution with endocardial structures that cannot be imaged with DTI, LDDMM will have trouble mapping fibers into and around these regions, where the LDRB algorithm would not.

Additional LDRB Capabilities

First, the LDRB algorithm is not limited by the model’s resolution (the mean edge length). However, the full capability of the LDRB algorithm to effortlessly incorporate fiber orientation into complex endocardial structures could not be exploited with the model of the canine ventricles employed in this study. This is due to the fact that available DTI data sets of entire mammalian ventricles are too low in resolution to capture trabeculations and papillary muscles; achieving the sub-millimeter resolution necessary to accurately capture these features with DTI is inherently expensive due to the long scan times.15 In Section 5 of the online supplement, we present the application of the LDRB algorithm to a high-resolution (25 µm isotropic resolution) MRI data set of rabbit LV. The LDRB algorithm was able to successfully define a continuous mathematical description for fiber orientation from the myocardium into the trabeculations and the papillary muscle according to R5 in the methods, while the MDRB algorithm we implemented could not (See S-Figures 18–22 in the online supplement). To the author’s knowledge, no other existing rule-based algorithm is capable of performing this task with automaticity.

Secondly, the angle β does not necessarily change linearly within the ventricular walls, and may vary from apex to base. Depending on the location within the myocardium (RV, LV, apex, base), the transmural function for β may be uniform,2 bimodal,14 or sigmoidal.22 Since the results in Table 1 revealed small differences between simulations with LDRB fiber orientation, generated using a simple linear transmural formulation for β throughout the entire myocardium, to the simulations with the DTI-derived fiber orientation, we did not explore the impact of such nonlinearities in β. However, with slight modification to the LDRB algorithm inputs in Eqs. (1)(4), the algorithm can accommodate differences in fiber orientation between the LV and RV by checking for ϕlv > ϕrv, as well as between the apex and base by expanding Eqs. (1)(4) to be functions of ψab. This approach would be especially useful in capturing changes in fiber orientation resulting from disease such as heart failure, where the sheet angle gradually steepens from apex to base in the septum of the failing heart.16

Thirdly, the LDRB algorithm can be easily modified to test various hypotheses regarding fiber orientation and its role in electrical activation and mechanical contraction. As presented above, the LDRB algorithm assumes the transverse angle, i.e., the imbrication angle of F around the ê1 axis, is zero. However, this may not always be the case for all hearts and could potentially affect transmural activation sequences.20,23 To incorporate the transverse angle into the orient function of the LDRB algorithm, only one additional rotation (as a function of ψ and ϕ) around the ê1 axis would need to be included between the α and β rotations. Furthermore, it is feasible that fiber orientation in the septum may only be continuous with the LV wall in order to optimize LV work,17 which would contradict R6. In this case, one needs to modify Eqs. (1)(4) to ensure that functions for the fiber orientation in the septum are continuous with only the fiber orientation in the LV wall.

Lastly, the LDRB algorithm is very versatile. In addition to computational modeling studies requiring accurate fiber orientation, it is becoming increasingly important to model spatial heterogeneity in ion channel density and conductivity throughout the myocardium.38 Assigning such properties throughout models with complex cardiac geometries poses the same challenges as encountered when assigning fiber orientation. From the results presented in this study, it is evident that the LDRB algorithm provides the accurate parameterization of the myocardium that is necessary to easily perform this otherwise daunting task.


The first limitation of the LDRB algorithm is that it is not applicable to computational models with surfaces that are not well defined. Well defined boundaries make it possible to apply Dirichlet boundary conditions and solve Laplace’s equation. However, this is typically only a problem when constructing models from imaging data sets acquired at low-resolution where surfaces bounding the thin RV wall may contain holes or be extremely jagged. Secondly, the LDRB algorithm produces artificially smooth fiber orientation and lacks natural variability in fiber orientation found in histology10,11,21,22,32 and DTI30,31 studies. It is feasible to add this variation to α and β throughout the myocardium using a random number generator, but in its absence we found the differences between activation maps with LDRB and DTI-derived fiber orientation to be subtle, so the effects of incorporating such variation are likely negligible for the pacing protocols typically used in simulation studies. Lastly, the LDRB fiber orientations need to be validated in simulations of mechanical contraction.


The LDRB algorithm is inexpensive, efficient, and easy to implement. It is capable of rapidly producing fiber orientation based on a-priori knowledge, and the fiber orientation obtained when using this method compares well both qualitatively and quantitatively with DTI-derived fiber orientation. The LDRB algorithm is ideal for cardiac electrophysiology studies, both in vivo and ex vivo, that utilize computational models with realistic descriptions of myocardial fiber orientation.

Supplementary Material



The authors would like to thank Dr. Edward Vigmond at the University of Bordeaux for his software Meshalyzer. This work was supported by grants AHA 10PRE3650037 to Jason Bayer, NIH R01 HL082729 and HL103428, and NSF CBET-0933029 to Natalia Trayanova, and FWF F3210-N18 and NIH R01 HL10119601 to Gernot Plank.



The online version of this article (doi:10.1007/s10439-012-0593-5) contains supplementary material, which is available to authorized users.


1. Alexander AL, Hasan KM, Lazar M, Tsuruda JS, Parker DL. Analysis of partial volume effects in diffusion-tensor MRI. Magn. Reson. Med. 2001;45(5):770–780. [PubMed]
2. Ashikaga H, Criscione JC, Omens JH, Covell JW, Ingels NB. Transmural left ventricular mechanics underlying torsional recoil during relaxation. Am. J. Physiol. Heart Circ. Physiol. 2004;286(2):H640–H647. [PMC free article] [PubMed]
3. Bayer JD, Beaumont J, Krol A. Laplace–Dirichlet energy field specification for deformable models. An FEM approach to active contour fitting. Ann. Biomed. Eng. 2005;33(9):1175–1186. [PubMed]
4. Beyar R, Sideman S. A computer study of the left ventricular performance based on fiber structure, sarcomere dynamics, and transmural electrical propagation velocity. Circ. Res. 1984;55(3):358–375. [PubMed]
5. Bishop MJ, Boyle PM, Plank G, Welsh DG, Vigmond EJ. Modeling the role of the coronary vasculature during external field stimulation. IEEE Trans. Biomed. Eng. 2010;57(10):2335–2345. [PMC free article] [PubMed]
6. Bishop MJ, Plank G, Burton RA, Schneider JE, Gavaghan DJ, Grau V, Kohl P. Development of an anatomically detailed MRI-derived rabbit ventricular model and assessment of its impact on simulations of electrophysiological function. Am. J. Physiol. Heart. Circ. Physiol. 2010;298(2):H699–H718. [PMC free article] [PubMed]
7. Bovendeerd PH, Arts T, Huyghe JM, van Campen DH, Reneman RS. Dependence of local left ventricular wall mechanics on myocardial fiber orientation: a model study. J. Biomech. 1992;25(10):1129–1140. [PubMed]
8. Caldwell BJ, Trew ML, Sands GB, Hooks DA, LeGrice IJ, Smaill BH. Three distinct directions of intramural activation reveal nonuniform side-to-side electrical coupling of ventricular myocytes. Circ. Arrhythm. Electrophysiol. 2009;2(4):433–440. [PubMed]
9. Cherry EM, Greenside HS, Henriquez CS. Efficient simulation of three-dimensional anisotropic cardiac tissue using an adaptive mesh refinement method. Choas. 2003;13(3):853–865. [PubMed]
10. Costa KD, Takayama Y, McCulloch AD, W J. Laminar fiber architecture and three-dimensional systolic mechanics in canine ventricular myocardium. Am. J. Physiol. 1999;276(2 Pt 2):H595–H607. [PubMed]
11. Fernandez-Teran MA, Hurle JM. Myocardial fiber architecture of the human heart ventricles. Anat. Rec. 1982;204(2):137–147. [PubMed]
12. Greenstein JL, Winslow RL. An integrative model of the cardiac ventricular myocyte incorporating local control of Ca2+ release. Biophys. J. 2002;83(6):2918–2945. [PubMed]
13. Han C, Pogwizd SM, Killingsworth CR, He B. Noninvasive reconstruction of the three-dimensional ventricular activation sequence during pacing and ventricular tachycardia in the canine heart. Am. J. Physiol. Heart Circ. Physiol. 2012;302(1):H244–H252. [PubMed]
14. Harrington KB, Rodriguez F, Cheng A, Langer F, Ashikaga H, Daughters GT, Criscione JC, Ingels NB, Miller DC. Direct measurement of transmural laminar architecture in the anterolateral wall of the ovine left ventricle: new implications for wall thickening mechanics. Am. J. Physiol. Heart Circ. Physiol. 2005;288(3):H1324–H1330. [PMC free article] [PubMed]
15. Helm P, Beg MF, Miller MI, Winslow RL. Measuring and mapping cardiac fiber and laminar architecture using diffusion tensor MR imaging. Ann. N. Y. Acad. Sci. 2005;1047:296–307. [PubMed]
16. Helm PA, Younes L, Beg MF, Ennis DB, Leclercq C, Faris OP, McVeigh E, Kass D, Miller MI, Winslow RL. Evidence of structural remodeling in the dyssynchronous failing heart. Circ. Res. 2006;98(1):125–132. [PubMed]
17. Hristov N, Liakopoulos OJ, Buckberg GD, Trummer G. Septal structure and function relationships parallel the left ventricular free wall ascending and descending segments of the helical heart. Eur. J. Cardiothorac. Surg. 2006;29S:S115–S125. [PubMed]
18. Hooks DA, Trew ML, Caldwell BJ, Sands GB, LeGrice IJ, Smaill BH. Laminar arrangement of ventricular myocytes influences electrical behavior of the heart. Circ. Res. 2007;101(10):e103–e112. [PubMed]
19. Hsu EW, Muzikant AL, Matulevicius SA, Penland RC, Henriquez CS. Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation. Am. J. Physiol. 1998;274(5 Pt 2):H1627–H1634. [PubMed]
20. Keller DUJ, Weiss DL, Dossel O, Seemann G. Influence of IKs heterogeneities on the genesis of the T-wave: a computational evaluation. IEEE Trans. Biomed. Eng. 2012;59(2):311–322. [PubMed]
21. Kim YH, Xie F, Yashima M, Wu TJ, Valderrbano M, Lee MH, Ohara T, Voroshilovsky O, Doshi RN, Fishbein MC, Qu Z, Garfinkel A, Weiss JN, Karagueuzian HS, Chen PS. Role of papillary muscle in the generation and maintenance of reentry during ventricular tachycardia and fibrillation in isolated swine right ventricle. Circulation. 1999;100(13):1450–1459. [PubMed]
22. LeGrice IJ, Smaill BH, Chai LZ, Edgar SG, Gavin JB, Hunter PJ. Laminar structure of the heart: Ventricular myocyte arrangement and connective tissue architecture in the dog. Am. J. Physiol. 1995;269(2 Pt 2):H571–H582. [PubMed]
23. Lombaert H, Peyrat J-M, Croisille P, Rapacchi S, Fanton L, Clarysse P, Delingette H, Ayache N. Statistical analysis of the human cardiac fiber architecture from DT-MRI. Proceedings of the 6th International Conference on Functional Imaging and Modeling of the Heart (FIMH’11); 2011. pp. 171–179.
24. Pierpaoli C, Basser PJ. Toward a quantitative assessment of diffusion anisotropy. Magn. Reson. Med. 1996;36(6):893–906. [PubMed]
25. Potse M, Dub B, Richer J, Vinet A, Gulrajani RM. A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart. IEEE Trans. Biomed. Eng. 2006;53(12 Pt 1):2425–2435. [PubMed]
26. Reddy JN. An Introduction to the Finite Element Method. 3rd ed. New York: McGraw-Hill; 2006. p. 766.
27. Rijcken J, Bovendeerd PH, Schoofs AJ, van Campen DH, Arts T. Optimization of cardiac fiber orientation for homogeneous fiber strain during ejection. Ann. Biomed. Eng. 1999;27(3):289–297. [PubMed]
28. Roberts DE, Hersh LT, Scher AM. Influence of cardiac fiber orientation on wavefront voltage, conduction velocity, and tissue resistivity in the dog. Circ. Res. 1979;44(5):701–712. [PubMed]
29. Rohmer D, Sitek A, Gullberg GT. Reconstruction and visualization of fiber and laminar structure in the normal human heart from ex vivo diffusion tensor magnetic resonance imaging (DTMRI) data. Invest. Radiol. 2007;42(11):777–789. [PubMed]
30. Scollan DF, Holmes A, Winslow R, Forder J. Histological validation of myocardial microstructure obtained from diffusion tensor magnetic resonance imaging. Am. J. Physiol. 1998;275(6 Pt 2):H2308–H2318. [PubMed]
31. Scollan DF, Holmes A, Zhang J, Winslow RL. Reconstruction of cardiac ventricular geometry and fiber orientation using magnetic resonance imaging. Ann. Biomed. Eng. 2000;28(8):934–944. [PMC free article] [PubMed]
32. Streeter DD, Spotnitz HM, Patel DP, Ross J, Sonnenblick EH. Fiber orientation in the canine left ventricle during diastole and systole. Circ. Res. 1969;24(3):339–347. [PubMed]
33. Trayanova NA. Whole-heart modeling: applications to cardiac electrophysiology and electro mechanics. Circ. Res. 2011;108(1):113–128. [PMC free article] [PubMed]
34. Vadakkumpadan F, Arevalo H, Prassl AJ, Chen J, Kickinger F, Kohl P, Plank G, Trayanova N. Image-based models of cardiac structure in health and disease. Wiley Interdiscip. Rev. Syst. Biol. Med. 2010;2(4):489–506. [PMC free article] [PubMed]
35. Vendelin M, Bovendeerd PH, Engelbrecht J, Arts T. Optimizing ventricular fibers: uniform strain or stress, but not ATP consumption, leads to high efficiency. Am. J. Physiol. Heart Circ. Physiol. 2002;283(3):H1072–H1081. [PubMed]
36. Vetter FJ, Simons SB, Mironov S, Hyatt CJ, Pertsov AM. Epicardial fiber organization in swine right ventricle and its impact on propagation. Circ. Res. 2005;96(2):244–251. [PubMed]
37. Vigmond EJ, Weber dos Santos R, Prassl AJ, Deo M, Plank G. Solvers for the cardiac bidomain equations. Prog. Biophys. Mol. Biol. 2008;96(1–3):3–18. [PMC free article] [PubMed]
38. Weiss DL, Seemann G, Keller DUJ, Farina D, Sachse FB, Dossel O. Modeling of heterogeneous electrophysiology in the human heart with respect to ECG genesis. In Proceedings of Computers in Cardiology. 2007:49–52.