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 effects^{1} 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 algorithms^{6,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 histology^{10,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 method^{3,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.