Development of an anatomically detailed MRI-derived rabbit ventricular model and assessment of its impact on simulations of electrophysiological function

Martin J. Bishop, Gernot Plank, Rebecca A. B. Burton, Jürgen E. Schneider, David J. Gavaghan, Vicente Grau, Peter Kohl

Abstract

Recent advances in magnetic resonance (MR) imaging technology have unveiled a wealth of information regarding cardiac histoanatomical complexity. However, methods to faithfully translate this level of fine-scale structural detail into computational whole ventricular models are still in their infancy, and, thus, the relevance of this additional complexity for simulations of cardiac function has yet to be elucidated. Here, we describe the development of a highly detailed finite-element computational model (resolution: ∼125 μm) of rabbit ventricles constructed from high-resolution MR data (raw data resolution: 43 × 43 × 36 μm), including the processes of segmentation (using a combination of level-set approaches), identification of relevant anatomical features, mesh generation, and myocyte orientation representation (using a rule-based approach). Full access is provided to the completed model and MR data. Simulation results were compared with those from a simplified model built from the same images but excluding finer anatomical features (vessels/endocardial structures). Initial simulations showed that the presence of trabeculations can provide shortcut paths for excitation, causing regional differences in activation after pacing between models. Endocardial structures gave rise to small-scale virtual electrodes upon the application of external field stimulation, which appeared to protect parts of the endocardium in the complex model from strong polarizations, whereas intramural virtual electrodes caused by blood vessels and extracellular cleft spaces appeared to reduce polarization of the epicardium. Postshock, these differences resulted in the genesis of new excitation wavefronts that were not observed in more simplified models. Furthermore, global differences in the stimulus recovery rates of apex/base regions were observed, causing differences in the ensuing arrhythmogenic episodes. In conclusion, structurally simplified models are well suited for a large range of cardiac modeling applications. However, important differences are seen when behavior at microscales is relevant, particularly when examining the effects of external electrical stimulation on tissue electrophysiology and arrhythmia induction. This highlights the utility of histoanatomically detailed models for investigations of cardiac function, in particular for future patient-specific modeling.

  • cardiac arrhythmias
  • magnetic resonance imaging

computational modeling is playing increasingly important roles in cardiovascular research. The direct combination of computational simulations with experimental measurements has provided important insights into the function of the heart in normal and pathological conditions (5, 22, 23, 55). One of the advantages of computational models is that they can be exploited to obtain information regarding multiple functional parameters throughout the three-dimensional (3-D) volume of the tissue, with resolution limited only by the level of detail represented in the model and the computational resources available to solve it. In contrast, experimental measurements are often restricted to recording a limited range of parameters, and observations are often affected by the very technique used to measure them (5, 45). Current state-of-the-art computational cardiac models incorporate gross cardiac anatomy as well as regional variations in prevailing myocardial cell orientation (usually, if incorrectly, referred to as “fiber orientation”). Such models have been developed for ventricles from the rabbit (54, 59), mouse (17), dog (32), pig (50), and human (39), to name but a few.

Although these models have proved sufficiently realistic to answer important questions regarding the relationship between global cardiac structure and function, they share a number of limitations. First, they follow a “one heart fits all” approach, lacking representation of interindividual variations in structure and anatomy that have important implications when experimental measurements from a particular preparation are compared with computational simulation results from a “generic” model (43). Second, these models are highly simplified in their histoanatomical complexity. By and large, they lack blood vessels, papillary muscles, valves, trabeculations, etc. The inclusion of such detail is essential both in providing an accurate like-for-like comparison between simulations and experiments as well as for providing an improved understanding of the roles played by these anatomical structures in cardiac function.

The requirement to consider anatomical complexity has been highlighted in a series of recent experimental investigations that have shown that fine-scale anatomical structures such as papillary muscles (24, 58) and blood vessels (41, 58) may act as important stabilizers of arrhythmias. These studies suggest that such structures can anchor phase singularities, promoting reentrant activity, and thus play an important role during the initiation and maintenance of heart rhythm disturbances. Furthermore, the role of cardiac structure at macro- and microscales has also been shown to be a key determinant of the cardiac response to external electrical stimulation. High-resolution optical mapping studies, for example, have demonstrated the existence of microscopic intramural virtual electrodes after the application of strong defibrillation shocks. These micro virtual electrodes have been structurally linked to the specific distribution of intramural collagenous planes within the tissue (18), which may facilitate successful defibrillation (19). At the macroscopic scale, the distribution of shock-end virtual electrode polarizations and corresponding vulnerability to strong electrical shocks has been shown to strongly depend on global cardiac ventricular geometry (43). Understanding the specific mechanisms by which anatomical obstacles stabilize arrhythmias, and how macro- or microscopic structures interact with electric field stimuli, could help in the development of improved therapies against arrhythmias, such as monophasic “unpinning” shocks (42), resonance drift (31), or low-voltage defibrillation (3).

To be able to probe the specific roles of fine-scale anatomical structures in cardiac function, it is necessary to include a matching level of structural detail in computational cardiac models. Recently, efforts have been focussed on developing techniques to construct 3-D computational cardiac models directly from magnetic resonance (MR) data, benefitting from the nondestructive nature of that technique (21). However, detail and anatomical complexity of the resulting models have been limited by the inherent resolution limits of MR datasets, with voxel dimensions usually exceeding 100-μm edge length. In the last few years, the advent of stronger magnets and refined scanning protocols has significantly increased the resolution of anatomical MR scans, such that smaller mammalian hearts now can have MR voxel dimensions of 20 μm or less (7, 36).

First attempts have been made to use these high-resolution MR datasets for the creation of computational models, both of the left ventricular (LV) free wall (7, 40, 57) and more recently of the entire ventricles (36, 56). When developing techniques to generate high-resolution structure-function computational models from such MR data, it is important that the methods are automated and as efficient as possible. The development of a computational MR-to-model pipeline provides the potential to generate models on a case by case basis, potentially leading to the generation of individualized models, a precondition for building probabilistic high-resolution MR cardiac atlases (7, 36). The efficient generation of individual models for separate preparations will facilitate the comparison of structure-function interactions in different hearts, thus assessing the link between interindividual heterogeneity in fine-scale anatomical structures and arrhythmia susceptibility. Publication of the computational methods used in the development of the MR-to-model pipeline, in full reproducible detail, is therefore an essential component in the wide-scale application of these methods to other datasets.

A vital feature yet to be addressed when developing this “next generation” of cardiac models is the assessment of the impact that additional anatomical complexity may have on simulations of cardiac function relative to more simplified ventricular models (54, 59). Such an assessment will put into context the conclusions drawn from the previous simulation studies and help to assess the cost-to-benefit ratio of the added overhead involved in work using more detailed models.

The main goal of this study was to present the methodology involved in the construction of a computational model of the rabbit ventricles directly from high-resolution MR data, in full detail, to investigate the effects of fine-scale anatomical heterogeneity (such as endocardial structures, blood vessels, etc) upon electrical activation sequences during different stimulation protocols. In addition, we provide the research community with full access to the new anatomical model and MR dataset. The model generation pipeline first involved the construction of a geometrically faithful representation of the ventricular geometry, obtained by segmentation of the MR dataset followed by finite-element mesh generation of the segmented voxel image stack. Functional electrical properties (in the form of anistropic conduction) were then incorporated into the geometrical model, which was then used for the simulation of electrical activity. Comparison with a more simplified model (derived from the same MR dataset but with endocardial structures and intramyocardial vessels and cleft spaces removed) allowed inquiry into key instances in which fine-scale anatomical complexity affects 1) activation sequences and electrical wavefront morphology under simple single-site pacing protocols; 2) the magnitude and distribution of electrical polarization levels, both intramurally and on epicardial/endocardial surface/structures, in response to field stimulation and how this impacts arrhythmia susceptibility; and 3) the complexity and duration of episodes of arrhythmogenesis. We focused on the resulting findings of key relevance from a model development point of view, highlighting circumstances in which the additional anatomical complexity allowed us to uncover mechanisms that may have important electrophysiological implications as well as those cases in which more simplified models were sufficient. We identified a number of applications in which fine-scale anatomical complexity is important, thus identifying instances in which it is beneficial to use detailed models to investigate such phenomena in future studies. By providing the community with full access to all the datasets acquired and developed in this study, we hope to facilitate more widespread investigation into the use of anatomically complex computational models in all aspects of cardiac electrophysiological research.

MATERIALS AND METHODS

MR Data Acquisition and Sample Preparation

All experiment investigations in this study had the necessary institutional ethical approval and were conducted in accordance with the UK Home Office guidance on the Operation of Animals (Scientific Procedures) Act of 1986.

The sample preparation of the heart before scanning has been described in detail previously (7). Briefly, female rabbit (≈1.2 kg) hearts were isolated and fixed during cardioplegic arrest (using elevated potassium levels, 20 mM) via coronary perfusion with 50 ml of fast-acting Karnovsky's fixative (2% formaldehyde and 2.5% glutaraldehyde mix) containing 2 mM gadodiamide (contrast agent). For MR imaging (MRI), hearts were stabilized in a nuclear MR tube using low-melting 1% agar containing 2 mM paramagnetic gadodiamide contrast agent.

MRI experiments were conducted on an 11.7-T (500 MHz) MR system consisting of a vertical magnet (bore size: 123 mm, Magnex Scientific, Oxon, UK), a Bruker Avance console (Bruker Medical, Ettlingen, Germany), and a shielded gradient system (548 mT/m, rise time: 160 μs, Magnex Scientific). Quadrature-driven birdcage coils with inner diameters of 28 and 40 mm (Rapid Biomedical, Wurzburg, Germany) were used to transmit/receive the MR signals. MR scans were performed using an established fast gradient echo technique (echo time/repetition time = 7.5/30 ms) allowing the acquisition of gap-free 3-D images (46).

Scans were acquired from the fixed heart with an in-plane resolution of 43 × 43 μm and an out-of-plane resolution of 36 μm. Reconstruction of the MR data was performed as described in Ref. 46. Two-dimensional (2-D) image planes were collated and used to generate a 3-D MRI-based dataset, which subsequently could be sectioned at any desired angle for visualization and analysis. Experimentally, 608 × 608 × 1,408 voxels were acquired in the x-, y-, and z-directions, respectively, which were zero filled to create the full MR data stack containing 1,024 × 1,024 × 2,048 voxels with corresponding edge lengths of 26.4 × 26.4 × 24.4 μm (zero filling is a commonly used MR technique to visually improve an image without in fact adding additional structural information.) The full dataset can be accessed and downloaded from the online Supplemental Material.1 Figure 1 shows cross sections through the MR dataset obtained from the fixed heart along the transverse (A), frontal (B), and sagittal (C) planes (orientations all relative to the organ's long axis, i.e., ignoring that the heart is positioned in the body at an oblique angle). Note that the transverse plane is normal to the global apex-base direction, the frontal plane is normal to the posterior-anterior direction, and the sagittal plane is normal to the global right-left direction.

Fig. 1.

Three-dimensional (3-D) magnetic resonance (MR) dataset, visualizing the data in the transverse (A), frontal (B), and sagittal (C) planes through the voxel stack. For clarity, the region of interest containing only the heart is shown (902 × 832 × 1,368 voxels). D: enlarged regions of the transmural plane image corresponding to the dashed square in A comparing the full-resolution MR image (left) with the downsampled image (right).

Segmentation of the Dataset

Segmentation is the critical first stage in extracting geometrical information from the MR dataset to facilitate the construction of a high-resolution anatomically detailed computational cardiac model. The process of segmenting a medical image may be thought of as labeling voxels in the image with a specific tag to define their association with different regions, objects, or boundaries. Initially, those voxels in the MR dataset that represented “tissue” were differentiated from those that represented “nontissue.” Discrimination between different tissue types, due to the similarity in intensity values in the MR images, must be done at a later stage, either exploiting the specific geometries of different structures (as described below) or using histological slices previously coregistered with the MR dataset (30). Thus, the initial goal is to create a binary (0, 1) mask allowing discrimination of the heart's tissue volume from the surrounding background volume and cavities as well as the extracellular cleft space and vessel lumens.

Dataset preprocessing.

Before segmentation, the MR data were first downsampled by a factor of 2 in each direction to produce a more manageable dataset containing 512 × 512 × 1,024 voxels. In addition to reducing the memory footprint of the data, downsampling also introduced a uniform smoothing that reduced “salt and pepper” noise artifacts, improving the performance of the segmentation algorithms. Figure 1D shows how downsampling removes noise without having an overly pronounced effect on the level of structural detail present in the image. Finally, the dataset was further reduced in size to 451 × 416 × 714 voxels by selecting only the heart as the region of interest, converted from 16-bit to 8-bit, and then rescaled to occupy the full range of grayscale intensities.

Segmentation pipeline.

Due to the relatively large size of the rabbit heart compared with the MR tube (28-mm diameter), parts of the MR images were overly affected by the bias field. Such an effect is caused by spatial variations in the coil sensitivity and largely affects tissue close to the edges of the MR tube, away from the “MR sweet spot” in the center of the tube. Consequently, a significant overlap was present in the intensity of tissue and background voxels in these regions, thus making simple segmentation approaches such as intensity thresholding impractical. Initial testing with bias field correction algorithms (52) did not completely address the issue and later proved unnecessary. Successful segmentation of the dataset was conducted using a combination of level-set segmentation filters, combined in a sequential processing pipeline to refine the segmentation at each stage.

Active contour/active surface methods are among the most powerful techniques for medical image segmentation. In particular, level-set methods use an implicit mathematical framework that solves many of the instabilities linked to the explicit formulation and allow easy generalization to N dimensions (48). The surface is embedded as the zero level set of a higher dimensional function [Ψ(X, t)], which can be evolved under the control of a differential equation. The zero level set Γ(X, t) = [(X, t) = 0] corresponds to the surface at any time during the evolution. A generic level-set equation to compute the update of the solution Ψ is given by the following: Math(1) where A, P, and Z are the advection, propagation, and expansion terms, respectively, and κ is the mean curvature. The scalar constants α, β, and γ weight the relative influence of each of the terms on the movement of the surface. Equation 1 represents a generic structure that can be adapted to the particular characteristics of the problem by selecting the appropriate weights and expressions for A, P, and Z. These can include different edge detection methods, intensity-based terms, surface shape constraints, etc. Equation 1 uses an iterative algorithm to reach the local minimum of a related cost function, and thus the procedure is sensitive to local minima of the cost function, where the algorithm may stop before reaching the desired global minimum. To avoid the effect of local minima, one can use an initialization close to the desired result. Otherwise, it is possible to choose a level-set function with few local minima; however, this usually means also a less accurate placement of the segmented surface. We used a combination of both approaches by building a pipeline of level-set algorithms, starting with the one most robust to poor initialization and finishing with the most accurate. In this way we obtained an algorithm that was both robust and accurate.

The full segmentation pipeline consisted of the following steps. First, an initial approximate segmentation was performed using a level-set threshold filter, following prior automated selection of seed points throughout the volume of the dataset (Fig. 2B). Seeds were automatically chosen based on explicit local intensity values using a very conservative threshold to avoid placing seeds within tissue voxels. The propagation term of the evolving contour (P in Eq. 1) in the threshold level-set method is dependent on user-defined lower and upper threshold limits. By selecting relatively high values for these limits, an approximate segmentation could be achieved without the fronts “leaking” into regions of tissue. However, as shown by the highlighted regions in the top left and bottom right corners of Fig. 2B, this also has the effect of leaving many areas of the background outside the volume of the heart incorrectly defined (due to the above-mentioned effects of the spatial variations in coil sensitivity) and of missing smaller structures within the myocardium.

Fig. 2.

Results of the automated sequential segmentation pipeline shown in the transverse (top) and frontal (bottom) slices. A: unsegmented MR dataset. B: output from the first stage in the segmentation pipeline, the threshold level-set filter, which acts as a good approximate initial segmentation. C: output from the geodesic active contour filter. D: final result of the segmentation pipeline following the Laplacian level-set filter.

To remove such artifacts, the second stage of the sequential segmentation processing pipeline was to use the output from the threshold level-set filter (Fig. 2B) as an initial level set for a geodesic active contour level-set filter. The geodesic active contour level-set algorithm attracts the contour to object boundaries (high values of the intensity gradient) and so is less affected by the intensity mismatches introduced by the bias field. The output from the geodesic active contour level-set filter is shown in Fig. 2C. The final stage of the segmentation pipeline used the output of the geodesic filter as an initial level set for a Laplacian level-set filter (using the second derivative of the intensity to detect edges), which is useful for fine tuning an initial segmentation but less so as an explicit region growing algorithm, as discussed above. Figure 2D shows the output of the Laplacian filter, the final result of the automated segmentation pipeline. Further details of the level-set filters used, along with the specific parameter values, can be found in the appendix (Segmentation Details).

Segmentation postprocessing.

The binary mask was further improved by a series of morphological operations. Isolated segmented regions were removed using a connected component algorithm to ensure that neighboring tissue voxels were directly connected to one another. This is important as all cardiomyocytes in the heart are electrically connected. Thin-walled blood vessels on the outer surface of the heart volume that were not sufficiently defined during the segmentation process were manually corrected. Finally, despite the care taken in the sample preparation, a small number of air bubbles were visible in the MR dataset and were also manually removed. Atrial tissue was then identified and removed from the segmented dataset, leaving only the ventricles, from which a detailed computational model was produced, in line with the goal of our study.

To achieve this, a series of 20–30 points was manually selected in every fifth frontal slice through the dataset, along the line believed to separate the ventricles from the atria, as shown in Fig. 3A. In most regions, this delineation was made by identifying a band of connective tissue along the atrioventricular border that appeared darker than the surrounding myocardium in the MR data. A smooth surface was then fitted to the points. This 2-D surface, dissecting the image volume, defined a secondary binary mask that separated the ventricles from the atria and other structures, such as great vessels. Figure 3B shows a 3-D representation of the binary mask defined by the 2-D surface within the 3-D object space. Figure 3C shows the result of convoluting the segmented voxel stack with the secondary binary mask of Fig. 3B to generate the final dataset, which contained only the ventricles.

Fig. 3.

Removal of the atrial tissue from the segmented voxel stack. A: frontal slice through the MR dataset showing an example of the manually placed points along the line believed to separate the atria and ventricles (pink dots). B: 3-D representation of the binary mask generated by the two-dimensional surface separating the ventricles and atria. C: frontal slice through the final segmented image stack with the atria removed.

The final stage of postprocessing involved the tagging of different structures in the heart based on known differences in electrical, mechanical, and/or functional properties, which are important for subsequent simulations of heart function. Manual approximate binary masks were created separately for the heart valves, papillary muscles, and larger trabecular structures (in places where they were sufficiently clearly separated from the endocardial surface) and superimposed on the segmented image stack to numerically tag each separate region with an identifier different to the normal myocardium. These numerical tags are preserved in the meshing stage (described in Finite-Element Mesh Generation below) and thus can be used to assign different functional properties to regions before simulation. Finally, the fragments of the free-running Purkinje system that were sufficiently well defined to be segmented (using the pipeline described in Segmentation Pipeline) were also manually tagged using the method described above. However, these structures were removed from the segmented image stack before mesh generation as incorporation of a functional anatomical representation of the Purkinje system requires further work to characterize the junctions between free-running and tissue-embedded components.

Computational tools used.

The automatic algorithms used within the sequential segmentation pipeline above were implemented using the Insight Toolkit library (ITK, www.itk.org), which offers a range of techniques for 3-D medical image segmentation, including several level-set algorithms. The ITK library was used as it provides a comprehensive suite of tested algorithms in an open source environment, which minimized implementation of new code and facilitated reproducibility. Further details of the filters and algorithms used can be found in the appendix (Segmentation Details). Where manual interaction was required (landmark selection, manual labeling, air bubble removal, and epicardial blood vessel patchup), Seg3D software was used (part of the SciRun toolkit, Scientific Computing and Imaging Institute, University of Utah). This was also used for the image visualization shown in Figs. 1 and 2. To ensure a high level of accuracy where manual interaction was required with the dataset, a Cintiq 12WX (Wacom, www.wacom.com/cintiq/12wx.cfm) interactive graphics tablet was used as a user interface. Matlab (The MathWorks) was used for a variety of tasks, including image downsampling and connected component analysis. In addition, an open source Matlab function (gridfit.m) was used to fit a smooth surface to the points defining the plane separating the atria and ventricles (for more information, see http://www.mathworks.com/matlabcentral/fileexchange/8998).

Finite-Element Mesh Generation

The numerical solution of the different systems of equations that model the heart's electrical and mechanical behavior requires the cardiac geometry to be discretized into a set of grid or mesh points. The majority of cardiac electrophysiology simulators (such as the one used in this study and described in Electrophysiological Simulations) use finite-element meshes as an input (23). Therefore, the final segmented and postprocessed dataset was used to generate a tetrahedral finite-element mesh, as described below.

Meshing software: Tarantula.

Tarantula (www.meshing.at) is an Octree-based mesh generator that builds unstructured, boundary fitted, locally refined, conformal finite-element meshes and has been custom designed for cardiac modeling. The underlying algorithm is similar to a recently published image-based unstructured mesh generation technique (40). Briefly, the method uses a modified dual mesh of an Octree applied directly to segmented 3-D image stacks. The algorithm operates fully automatically (no requirement for user interaction) and generates accurate volume-preserving representations of arbitrarily complex geometries with smooth surfaces. Directly meshing the segmented image stack removes the need to generate tessellated representations of the cardiac surface, as required by Delaunay-based mesh generators (38). Avoiding this additional step is a distinct advantage in the development of high-throughput cardiac model generation pipelines (36). A full description of the parameters and their numerical values used is provided in the appendix (Mesh Generation Details). In addition to meshing the myocardial volume, a mesh of the volume conductor surrounding the ventricular mesh was also produced. This extracellular mesh consists not only of the space within the ventricular cavities and the vessels, etc., but also defines an extended bath region outside of the heart, which can be used to model an in vitro situation where the heart is surrounded by a conductive bath, such as during optical mapping experiments (13). Spatial adaptivity within the surrounding volume mesh reduces the computational load by reducing the spatial resolution of the mesh with the distance to the cardiac surfaces (40).

Results of the mesh generation.

The total mesh produced by Tarantula from the final segmented voxel stack consisted of 6,985,108 node points, making up 41,489,283 tetrahedral elements. The myocardial mesh (representing the “intracellular” or “cytosolic domain”) consisted of 4,306,770 node points, making up 24,199,055 tetrahedral elements with 1,061,480 bounding triangular faces. The myocardial mesh had a mean tetrahedral edge length of 125.7 μm. It took 39 min to generate the mesh using two Intel Xeon processors operating at 2.66 GHz (36). Note that Tarantula was used to produce the mesh using a node discretization of ≈3 voxels (see the appendix, Mesh Generation Details). Using a finer discretization (i.e., 1 voxel) would result in a mesh with a vastly increased number of degrees of freedom (i.e., more nodes and elements), which would, in turn, prove problematic for computational tractability, both at the meshing stage and for the simulation of electrical activity. However, this means that fine structures (only 1–2 voxels thick) would become degraded. To prevent this, a “halo” of voxels was added to all outer surface structures (2 voxels in the case of the epicardium, to prevent the lumina of thin-walled epicardial blood vessels from becoming connected to the exterior).

Figure 4 shows the final ventricular mesh. Here, and throughout, Meshalzyer software (carp.meduni-graz.at/) was used for finite-element visualization. Figure 4A shows a standard anterior view of the outer (epicardial) surface of the mesh (left) along with cuts along both frontal (middle) and transverse (right) clipping planes to expose transmural and endocardial structures. Notable features include the large blood vessels close to the epicardial surface and within the midmyocardium as well as the papillary muscles and valves. Figure 4B shows a zoomed-in region of the interior of the LV (in a posterior to frontal view), where the triangular faces of the intracellular finite-element mesh are shown in blue, highlighting the level of detail with which the complex endocardial surfaces are defined. Finally, segmentation tags that were used to label the different regions in the heart were directly transferred from the regular image stack to the unstructured grid. These tags can be subsequently used to define electrophysiological properties on a per region basis. Figure 4C shows tagged structures such as the papillary muscles (green) and valves or cordae tendinae (blue) in a cut along a frontal clipping plane (posterior view) to enable the visualization of important endocardial structures. For the computation of electrical activity within the ventricles, an additional high-resolution mesh was constructed in which the valves were removed from the tagged image stack before meshing, as the valves are considered electrically silent. This mesh is available for download via the online Supplemental Material.

Fig. 4.

Tetrahedral finite-element rabbit ventricular mesh. A: visualization of the final ventricular finite-element mesh from a standard anterior view (left) along with cuts along the frontal (middle) and transveral (right) planes to expose endocardial structures. B: highlighted region from an exposed clipping plane in a posterior view demonstrating the level of detail in the finite-element mesh on the endocardial surfaces. C: cut along a frontal clipping plane in a posterior view showing the tagged structures of the papillary muscles (green) and valves or cordae tendinae (blue). Note that the valves do not retain their in vivo shape due to the preparation of the heart. D: simplified rabbit ventricular finite-element model shown from a standard anterior view (left) along with cuts along frontal (middle) and transverse (right) clipping planes. The helix angle α and vectors z, u, v, and af (defining the global apex-base, transmural, circumferential, and fiber directions, respectively) were used in the calculation of fiber orientation explained in Incorporation of Fiber Orientation Information.

Simplified Ventricular Model

A major goal of this study was to directly assess the impact of structural detail in the high-resolution finite-element rabbit ventricular model on electrical propagation dynamics throughout the volume of the ventricles. To achieve this goal, a simplified model was constructed for comparison. The simplified model was developed from the same MR dataset (and thus was geometrically similar in terms of overall size, cavity dimensions, etc.) using the pipeline described in Segmentation of the Dataset. However, in the segmentation pipeline, seed points were only placed within the two cavities and outside the volume of the heart; thus, only the epicardium and endocardium were segmented, but not higher-level details such as the blood vessels or extracellular cleft spaces. In addition, the manually created binary masks used to delineate the papillary muscles, larger trabeculae, and valves (see Segmentation postprocessing) were used to remove such structures from the binary segmented voxel image stack before meshing (atrial tissue was removed using the same binary mask as shown in Fig. 3B). The result was a simplified model with a level of anatomical complexity similar to that of previously established rabbit ventricular models (54, 59).

Tarantula was used with similar parameters as described in the appendix (Mesh Generation Details) to ensure that the simplified and anatomically complex meshes had similar nodal discretizations. The total simplified mesh produced by Tarantula consisted of 6,449,098 node points, making up 38,299,065 tetrahedral elements. The intracellular mesh consisted of 4,233,856 node points, making up 24,348,259 tetrahedral elements with 642,282 bounding triangular faces. The mean element edge length of the intracellular mesh was 128.1 μm. Figure 4D shows the simplified mesh (left), including views where the clipping planes have been used to expose ventricular cavities (middle and right; compare with Fig. 4A).

Incorporation of Fiber Orientation Information

The preferential arrangement of cardiac cells within the heart (“fiber orientation”) is known to play crucial roles in electrical and mechanical activity (23). Ideally, quantitative measured data should be used to account for the structural orthotropy of the cardiac tissue architecture. Methods such as diffusion tensor MRI (DT-MRI) are well suited for this purpose, and efforts to include DT-MRI into the processing pipeline as an additional imaging modality are currently being undertaken (36). However, for processing the data used in this study this method was not yet available. In such cases, fiber information can be incorporated into the model using a rule-based approach based on a priori assumptions. Here, we present a variation of such a mathematical model for ventricular fiber orientation, as originally presented by Potse et al. (39). The model is based on the experimental data collected by Streeter et al. (51), who found that cardiomyocytes generally tend to be aligned tangential to the surface and with an angle relative to the equatorial plane, which depends on the depth within the myocardial wall. The following method was used to incorporate fiber architecture into both the anatomically complex model and simplified model.

Definition of surfaces.

The algorithm developed by Potse et al. (39) uses a normalized transmural distance across the myocardial wall to define both the plane in which the myofibers lie and also the angle of inclination of the fiber within the plane (termed the “helix” angle, as shown in Fig. 4D). Therefore, a vital postprocessing step on the finite-element mesh (generated as described above in Finite-Element Mesh Generation) was to discriminate between the epicardial, LV endocardial, and right ventricular (RV) endocardial surfaces, such that the distance of each point within the tissue to each respective nearest surface could be determined to allow the construction of a transmural distance map.

Attempts to discriminate between the different surfaces within the list based on their discrete connectivities were unsuccessful. After segmentation, small connections remained between the exterior surfaces and coronary vasculature surfaces and were also formed by any slight leaks or holes introduced in the segmentation or meshing stages. Instead, discrimination between the different surfaces within the final ventricular mesh was obtained by performing a secondary segmentation of the MR data where only information regarding the two cavities and the outer surface was obtained and relating this back to the surface triangles in the original mesh. The specific details of this method can be found in the appendix (Surface Discrimination Algorithm).

Rule-based algorithm for incorporating fiber architecture.

After successful discrimination between the epicardial and LV and RV endocardial surfaces, the first stage in the fiber architecture algorithm of Potse et al. (39) involved constructing a distance map, which was used to define the transmural direction at every point within the tissue. To do so, the minimum distance of each node point within the mesh to both the epicardium (depi) and endocardium (dendo) was computed and used to define a normalized thickness parameter (e) as follows: Math(2) In the septum, depi was taken to be the minimum distance to the endocardial surface of the RV (dRVendo) and dendo taken to be the minimum distance to the endocardial surface of the LV (dLVendo), replicating the apparent structural continuity of the LV and septum seen in previous studies (27, 47). Note that the points labeled as residing within the septum were determined as those for which depi > dLVendo and depi > dRVendo and those within the LV and RV free walls were determined as those for which dRVendo > dLVendo and dLVendo > dRVendo, respectively. To avoid sudden discontinuities in e at boundaries between regions, the value of e at each node was averaged with all of its immediate neighbors (i.e., those neighbours sharing a common element).

The gradient of e within each tetrahedral element was used to define the transmural direction (u). The vector product of u with a unit vector in the z-direction (the global apex-base direction) then defined the circumferential direction (v) or, in other words, the orientation of cardiac fibers in the midmyocardium with zero helix angle (see Fig. 4D for vector definitions). The direction w was then defined as being perpendicular to v and u, i.e., w = u × v. To account for the transmural variation in the fiber vector direction through the myocardial wall, the fiber direction at the center of each tetrahedral element was further determined by a helix angle (α) (39) as follows:Math(3) where eav was taken to be the average of e of the nodal values making up the tetrahedra. A rotation of v and w of α radians about the u-axis finally defined the real fiber vector direction (af) and the normal to the myolaminar (an) at that point. R was set equal to π/3 for the LV and septum and π/4 for the RV in accordance with Ref. 39. Figure 5A shows the orientation of fibers through sections of the LV wall, RV wall, and septum. The fine arrows in Fig. 5A show the explicit direction of the fibers, whereas the color bar indicates the out-of-plane rotation of the fiber vector.

Fig. 5.

A: vector representation of fiber directions (defined via the algorithm described in Rule-based algorithm for incorporating fiber architecture) within a section of the left ventricular (LV) wall, right ventricular (RV) wall, and septum. The color bar represents the out-of-plane rotation of the fiber vectors. B: Visualization of fiber vectors within the papillary muscle following assignment via the method described in Anatomically based model of fiber architecture within the papillary muscles. Fiber vectors are shown here as white arrows for clarity. Visualization was with Meshalzyer software.

Anatomically Based Model of Fiber Architecture within the Papillary Muscles

Although the parametric description of fiber architecture by Streeter et al. (51) [described in the implementation of the algorithm (39) above] provides a detailed representation of the transmural rotation of cardiac fibers throughout the myocardial wall, it does not address fiber architecture within the papillary muscle and large trabecular structures. It is known that the histoarchitecture of the papillary muscles greatly differs from that seen in the bulk myocardium, with cells tending to align parallel to the muscle axis (i.e., in the primary direction of contraction). Therefore, for tetrahedral elements tagged as papillary muscle, a fiber orientation parallel to the axis of the muscle was assigned. To determine the overall global axial direction at any point within the papillary muscle, the structure tensor was used.

Structure tensor.

The structure tensor is a robust method for determining localized structure based on neighboring gradient fields (15). Given a function f, which at the point x, y, z has local gradient ▽f = (fx, fy,fz), the structure tensor S is defined as follows: Math(4) where Kσ is a Gaussian function with standard deviation σ and T represents the transpose of the vector ∇f. The use of the Gaussian effectively performs a spatial average of neighboring gradient vectors for every entry in the total matrix. The eigenvalues (λ1, λ2, λ3) and corresponding eigenvectors (e1,e2,e3) of S indicate the local distribution of the gradient vectors. Depending on the relative sizes of the eigenvalues, overall global structures can be characterized (15). For the specific case that λ12 > λ3 ≈ 0, the overall global structure is characterized by a line-like structure (or cylinder), with the eigenvector e3 defining the cylinder axis. In the case of papillary muscles, which are approximately cylindrical in structure, this axis is used to define the papillary muscle fiber direction.

Implementation.

Before using the structure tensor method described above, it was necessary to define a gradient field associated with the structure of the papillary muscle to define the axial papillary muscle direction. A radial vector, i.e., a vector that would exit the structure normal to the surface, was defined for each tetrahedral finite element within the papillary muscles. The most robust way to compute such a radial vector (rr) was to define the vector as the line connecting the centroid of the tetrahedral element (position vector rc) with the geometrically closest surface node point (position vector rn) taken from the list of all LV and RV endocardial points calculated as described in Definition of surfaces. Thus, rr was simply given as follows: rr = rcrn.

As our finite-element ventricular mesh was unstructured, a Gaussian smoothing of the components of the matrix [▿f × (▿f)T] was performed on an unstructured domain. To do this, the size of the Gaussian kernel associated with each element was taken to include all neighboring tetrahedral elements that were directly connected to that element (i.e., directly shared node points). The value of the Gaussian function at each neighboring tetrahedra was then computed relative to the geometrical distance between the centroids of this and the central tetrahedra. The value of σ was optimized to 0.3 mm, based on visual inspection of the resulting fiber vectors. Although this method can potentially lead to artifacts due to the irregular distance between element centers, in our case this was not a significant concern as the good quality of the mesh prevented this from affecting the result. In addition, the use of a directly connected method ensured that the smoothing did not include neighboring elements associated with close by, but unconnected, papillary structures.

After the smoothing, structure tensor S was formed for each element and eigenvalues and eigenvectors were computed. The eigenvector with the corresponding smallest eigenvalue then defined the required axial fiber direction. Figure 5B shows a cut through the surface of a papillary muscle with the fiber directions calculated as described above, shown as fine white arrows for clarity. The same method for incorporating fiber vectors aligned with the axial direction of a cylindrical structure was also used to assign fiber orientation to the larger trabeculations present in the dataset, which were also tagged as described in Segmentation postprocessing.

Fiber vector smoothing.

To avoid sudden changes in fiber orientation in the insertion point of the papillary muscles with the myocardial wall as well as at the junctions between different wall regions (LV, RV, and septum), all fiber vectors were additionally smoothed by constructing a Gaussian-smoothed structure tensor, as described above. In this case, the gradient function f in Eq. 4 was defined by the unsmoothed fiber vector at each element, and the eigenvector corresponding to the largest eigenvalue of the resulting structure tensor was used to define the smoothed fiber direction within each element.

Electrophysiological Simulations

Governing equations of electrical activation.

Electrical activation throughout the ventricular models was simulated using the following bidomain equations (37):Math(5) Math(6) Math(7) Math(8)where φi and φe are the intracellular and extracellular potentials, respectively, σ̂i and σ̂e are the intracellular and extracellular conductivity tensors, respectively, β is the membrane surface-to-volume ratio, Im is the transmembrane current density, Istim is the current density of the transmembrane stimulus used to initiate an action potential, Ie is the current density of the extracellular stimulus, Cm is the membrane capacitance per unit area, Vm is the transmembrane potential, and Iion is the density of the total current flowing through the membrane ionic channels, pumps, and exchangers, which depends on the transmembrane potential and on a set of state variables (η). At the tissue boundaries, electrical isolation is assumed, which is accounted for by imposing no-flux boundary conditions on φe and φi.

In those cases where the cardiac tissue is surrounded by a conductive medium, such as a perfusing bath in which the heart is submerged, an additional Poisson equation has to be solved:Math(9)where σb is the isotropic conductivity of the conductive medium and Ie is the stimulation current injected into the conductive medium. In this case, no-flux boundary conditions are assumed at the boundaries of the conductive medium, whereas continuity of the normal component of the extracellular current and the continuity of φe are enforced at the tissue-bath interface. The no-flux boundary conditions for φi remain the same.

Numerical solution techniques.

The Cardiac Arrhythmia Research Package (CARP) (61) was used to discretize and numerically solve the bidomain equations. The specifics of the spatiotemporal discretization strategy and the numerical solution scheme have been described in detail elsewhere (34, 36, 61). Membrane dynamics were modeled using the recent Mahajan-Shiferaw model of the rabbit ventricular cell (29), which is equipped with an electroporation current and a hypothetical K+ current that activates at larger positive polarizations beyond +160 mV to reproduce the asymmetry of the membrane response to strong shocks delivered during the plateau phase of the action potential (1). Simulations were performed using resources at the Oxford Supercomputing Center (http://www.oerc.ox.ac.uk/resources/osc) as well as the Distributed European Infrastructure for Supercomputing Applications Extreme Computing Initiative project muHeart. For computational performance and benchmarking, see Ref. 36.

Representing electrical tissue anisotropy.

The information regarding locally prevailing cell alignment (the fiber vector af, see Rule-based algorithm for incorporating fiber architecture) was used to define the intracellular conductivity tensor σ̂i using the formulation in Ref. 12. Intracellular and extracellular conductivities were taken to be 0.174 and 0.625 S/m parallel to the myofiber and 0.019 and 0.236 S/m transverse to the fiber, respectively (11). Conductivity in the bath surrounding the ventricles was taken to be 1.0 S/m.

Stimulation Protocols

Pacing.

The stimulation protocol used to elicit wavefront propagation in both the complex and simplified rabbit ventricular models included stimulation at a site in the RV free wall close to the apex to elicit wavefront propagation in approximately an apicobasal direction (termed apical pacing). A transmembrane current pulse of 50 ×10−4 μA/cm3 was delivered over 1-ms duration. The instant at which Vm crossed a level of −20 mV at each node was recorded and used to construct an activation map plot.

Induction of arrhythmia.

An additional protocol was used to induce an arrhythmia in the ventricles. First, the ventricles were paced at the apex at the same site as in the apical pacing protocol. Then, at a given coupling interval (CI) of 195 ms, when the refractory isoline appeared roughly halfway between the apex and base, an external electrical stimulus of a shock strength of 2 V/cm and 5-ms duration was applied via two plate electrodes in an anterior-posterior configuration. The cathodal stimulus was applied via the anterior plate electrode by current injection. The posterior plate electrode was grounded and served as the anode. Such a protocol was used due its well-documented propensity for arrhythmia induction (43).

RESULTS

Activation Sequences and Wavefront Morphology

Figure 6 shows a comparison of simulated Vm distributions after apical stimulation in both the complex (top half) and simplified (bottom half) ventricular models at different times after the initial stimulus. Figure 6 shows the full ventricular epicardial surface (top rows) and the ventricles cut along the transverse (middle rows) and frontal planes (bottom rows). We will refer to these as the epicardial, transverse, and frontal views, respectively, in the explanation below. The posterior view was chosen to allow the visualization of important endocardial structures, in particular, the major papillary muscles.

Fig. 6.

Propagation of activation wavefronts after an apical stimulus. Shown are snapshots of membrane potential (Vm) distributions within the anatomically complex model (top half) and simplified model (bottom half) at different instances in time after stimulation close to the apex. Each image shows the Vm distribution in an epicardial view (top row) and where clipping planes have been used in transverse (middle row) and frontal (bottom row) plane cuts.

Although the global patterns of activation are similar in the two models, there were a number of important differences in activation sequences due to differences in the underlying additional anatomical detail of the models. Immediately after the stimulation in the complex model, propagation spread not only up the LV and RV free walls but also along the small trabeculations close to the apex of the RV where the stimulation occurred (13 ms, frontal view). These trabeculations linked the RV to the septum and thus caused the septum to depolarize in advance of the global propagating wavefront (32 ms, frontal view).

As a result, activation propagated significantly faster along the septum in the complex model compared with the simplified model (compare 32–58 ms, frontal views). For example, the activation time of point 1 (shown in the 13-ms frontal view of Fig. 6) in the complex model occurred at 36 ms after the stimulus compared with 46 ms for the corresponding point in the simplified model. In contrast, the propagation speed within the LV/RV myocardial walls of the two models were similar: the activation times of points 2 and 3 were, respectively, 64 and 42 ms for the complex model and 60 and 40 ms for the simplified model. In addition to the similarity in speeds, wavefront morphology across the LV/RV walls in both models exhibited a similar characteristic “V” shape (47 and 58 ms in the frontal views, for example). This is due to the inclusion of a realistic representation of fiber architecture in both models, resulting in the activation advancing more quickly close to the epicardial/endocardial surfaces, as fibers in these regions have a greater helix angle and thus are more closely aligned with the global direction of propagation (the apex-base direction, in this case). The corresponding effect was not present in the septum of the complex model due to the premature depolarization of the center of the septum, as described above, causing asymmetrical propagation in this region.

As the wavefront reached the base of the ventricles, the difference in the activation times between the septa of the two models was still apparent: the activation times of point 1′ were 77 and 85 ms for the complex and simplified models, respectively. However, the close similarity of the activation times within the LV and RV walls seen at points 2 and 3 was less apparent, as the faster depolarization of the septum in the complex model subsequently spread out from the center of the septum to the LV/RV free walls (evident in the 58-ms transverse image), resulting in the slightly earlier activation of these regions in the complex model: activation times of points 2′ and 3′ are 100 ms and 81 ms in the complex model compared with 104 and 85 ms in the simplified model (Fig. 6).

A further consequence of the faster propagation within the septum was focal depolarizion of the epicardial surface, appearing in advance of the main propagating wavefront (47-ms epicardial view). This focal activation spread out from the septum and caused the wavefront on the epicardium to be significantly less planar (increased curvature) throughout the remainder of the activation sequence, relative to the more uniform and planar wavefront in the simplified model, particularly evident in the epicardial 58- to 72-ms images.

Response of the Ventricles to Electrical Field Stimulation

Figure 7 shows a comparison of Vm distributions in the complex (top half) and simplified (bottom half) ventricular models just before (195-ms images) and after (200- to 215-ms images) the application of an external electrical stimulus (CI = 195 ms, SS = 2 V/cm). The views in Fig. 7 represent the clipping planes in the transverse (top row) and frontal (bottom row) planes.

Fig. 7.

Distribution of electrical polarization after an external stimulus. Shown are snapshots of Vm distributions within the anatomically complex model (top half) and simplified model (bottom half) at different instances in time after the application of an external electrical stimulus with a coupling interval (CI) = 195 ms and SS = 2 V/cm. Each image shows the Vm distribution in a posterior view where clipping planes have been used in transverse (top row) and frontal (bottom row) plane cuts. Note that white represents polarization levels >40 mV.

Although the distribution of potential was similar between both models just before the shock (195-ms images), immediately after the shock (200 ms) significant differences between polarization states existed. Most noticeably, when we compared the 200-ms images, regions of the anterior endocardial surface of the complex model were appreciably less positively polarized than those of the simplified model. This is because at shock end the endocardial structures in the complex model (papillary muscles/trabeculations) have formed individual virtual electrodes, where part of the structure has become positively charged and the other half negatively. Where these structures then attach to the endocardial surface they introduce a heterogeneous distribution of Vm in the subendocardial layers of the myocardial wall (examples are shown by the enlarged dashed rectangular boxes in Fig. 7). The overall effect is that the endocardial surface of the complex model has significantly more heterogeneity in the polarization levels across the surface compared with the simple model, and, in addition, parts of the endocardial surface (where the endocardial structures attach to the wall) are less strongly polarized in the complex model relative to the corresponding regions in the simplified model. The complex endocardial structures therefore appear to “protect” the endocardial surface from becoming strongly polarized (note that a similar effect was also seen on the posterior endocardial surface). Furthermore, at the end of the stimulus (200-ms images), the existence of small-scale virtual electrodes were found around the blood vessels and extracellular cleft spaces in the complex model. These microvirtual electrodes introduced heterogeneities in the Vm distribution within the bulk of the myocardial wall, causing a further disparity between the complex and simplified models.

Due to saturation of the color scales, an effect that was not immediately evident from the results shown in Fig. 7 (200-ms images), the posterior epicardial surface of the complex model was significantly less positively polarized than that of the simplified model at shock end. For example, considering the polarization states of epicardial points lying within a square 1.2 × 1.2-cm region on the surface [roughly corresponding to the field of view in an optical mapping experiment (5); shown as a white square on the top left image in Fig. 6], we found that just 9% (88%) were polarized > 150 (>100 mV) in the complex model compared with 99% (100%) in the simplified model for respective polarization levels. Such an effect was also present on the anterior epicardium, which experienced a less negatively polarized state at shock end in the complex model: 23% (81%) of points within a corresponding region on the anterior epicardial surface had polarization levels less than −200 (less than −150 mV) compared with 79% (95%) in the simplified model.

Consequently, 5 ms after the end of the stimulus (205-ms images), the posterior epicardial surface of the complex model recovered to a more negative state than the simplified model, which was now evident in the transverse slice images: 0% of points within the square field of view had polarization levels > 50 mV compared with 43% in the simplified model. We also noted that the endocardial surface of the complex model recovered to a lower potential than the simplified model, as it also reached a less positively polarized state at shock end (as described above).

Intramurally, despite the previously noted existence of microvirtual electrodes in the complex model at shock end surrounding intramural structures, postshock the distribution of Vm appeared similar between the two models, as these transient effects dispersed. However, important to note in the 205-ms image of the complex model (frontal slice) is the spreading of activity along a trabeculation that links the depolarized RV wall with the septum, which itself is largely at resting state. By the 215-ms image, this activity has caused the septum to depolarize, introducing an additional wavefront of activity that is not present in the simplified model. At this time, the 215-ms image also highlights the existence of other new shock-induced activation wavefronts within the ventricles, which have different locations between the two models due to differences in anatomical geometry.

Arrhythmia Induction and Dynamics

Postshock evolution of Vm within the ventricles progressed into sustained episodes of reentrant activity in both models. Figure 8 shows the distribution of Vm from a posterior epicardial view for both the complex (top half) and simplified (bottom half) models at different points in time.

Fig. 8.

Progression of electrical activation dynamics during stimulus-induced arrhythmogenesis. Shown is the Vm distribution on the posterior epicardial surface of the ventricles of the anatomically complex model (top half) and simplified model (bottom half) at different time instances after arrhythmia induction for the externally applied stimulus shown in Fig. 7 (CI = 195 ms, SS = 2 V/cm).

The differences in the development of new activation wavefronts formed after the stimulus (highlighted in the 215-ms image in Fig. 7) between the two models led to the formation of different initial reentrant patterns of activity at specific locations in the model. However, despite these localized differences, the global effect of the stimulus was such that the initial reentrant activity on the posterior epicardial surface in both models was dominated by a large propagating wavefront moving from the LV to the RV, although the morphology of the wavefront was different in each case, as were trajectory and speed.

For example, the 260- and 275-ms images of Fig. 8 demonstrate how the posterior epicardial surface of the simplified model was more strongly positively polarized by the shock than that in the more complex model. Here, we found that large areas of the epicardium (particularly at the apex, the base, and toward the RV) were still refractory in the simplified model, whereas they were almost fully recovered in the complex model. Consequently, in the complex model, activity spread rapidly, unhindered, across the posterior surface of the ventricles. In the simplified model, however, activation spread more slowly, experiencing conduction block in regions close to the apex or base that were still refractory. On the second cycle of the resulting reentry (from ∼400 ms onward), the heterogeneous distribution of refractory tissue present in the simplified model caused the wavefronts to fractionate and break up (clearly evident in the 590-ms images). The arrhythmia subsequently progressed into a more disorganised state with the existence of multiple wavefronts, which broke frequently. These patterns were also clear in the movies provided in the Supplemental Material.

In the complex model, the initial wavefront propagating across the posterior epicardium after the shock (260- and 275-ms images) was not hindered by refractory tissue to the same extent as in the simplified model. As a result, on the second reentrant cycle, there was no such dispersion of refractory tissue, and the reentrant wavefront continued to propagate. Consequently, a stable spiral wave reentrant arrhythmia pattern resulted, which was maintained for seven reentrant cycles in the present example. During this time, the phase singularities did not appear to be anchored to any significant anatomical structure. However, after ∼1,200 ms, the rapid activation rate of the tissue eventually led to wavebreak (clearly evident in the 1,430- and 1,535-ms images), and the arrhythmia also progressed into a more disorganized state with multiple wavebreaks, which appeared to be of a similarly disorganized state as that witnessed in the simplified model (although they were brought about by at least partially different mechanisms).

Finally, the episodes of arrhymias in the two models were found to terminate (i.e., extinction of the final remaining wavefront in the ventricular tissue) at 1,265 ms for the simplified model compared with 1,610 ms for the complex model.

DISCUSSION

In this study, we present a full computational pipeline for the generation of a highly detailed, anatomically complex rabbit ventricular model, derived directly from high-resolution MR data, and demonstrate functional implications of the additional level of anatomical heterogeneity on electrical behavior of the ventricular model during simple pacing protocols and after the application of external electrical field stimulation. In comparing the results from the anatomically complex model with a more simplified representation of the same dataset, we provide a first glimpse into differences made by fine structural detail on the electrical activation dynamics as well as elucidate instances where global activation patterns are similar.

Relevance of Methodological Model Development

The differences seen in the activation sequences and wavefront morphologies during pacing and in the shock-induced polarization states and arrhythmia induction after the application of external stimuli support the introduction of higher levels of anatomical complexity in certain aspects of computational studies of cardiac electrophysiology. The validation and comparison of computational simulation results with experimental data is an essential component of cardiac electrophysiological research and development. This is of particular relevance for the use of models desgined to represent cardiac anatomical structure on a case by case basis for individual experimental preparations and will assist in the like-for-like comparison between model and experiment. In addition, these allow quantitative prediction of cardiac electrical activation at very high spatial (∼100 μm) and temporal (∼50 μs) resolutions and allow the identification of underlying mechanisms in a manner that cannot be achieved in experimental investigations or using less detailed computational models.

The MR-to-model pipeline described in this study provides the potential for individualized computational model construction, representing experimentally obtained functional (e.g., optical mapping) and structural (MR) data on individual preparations. This may help to elucidate the role of individual anatomical heterogeneity in episodes of anatomical microreentry or shock-induced arrhythmogenesis in like-for-like wet and dry observations. This may require even higher-resolution representation of cell and tissue type differences throughout the model, for example, using histological identification of tissue components and morphology. Furthermore, the availability of detailed models opens the door for the quantitative analysis of coupled multiphysics problems, such as electromechanical coupling and the resulting fluid-structure interactions. The prescence of fine-scale structures such as the papillary muscles and myocardial valves, which play important roles in electromechanical and fluidic function, will be essential to represent fully coupled whole heart behavior in the future.

Methodological Improvements to the Computational Pipeline

One of the major goals of this project was to establish a generic computational pipeline for the generation of individual cardiac models from MR data (36). Although manual interaction was still required at certain stages, most of the pipeline is now fully automated, providing the potential for its direct application to high-throughput data analysis. Despite the fact that the specific values of parameters used [for example, in the segmentation filters of the appendix (Segmentation Details)] will inevitably need adjustment, the general steps in the pipeline are sufficiently robust to be used on other datasets, even for different species [for example, on rat MR data (6)].

The next stage of the computational model generation pipeline is the incorporation of whole heart serial histology data into the model (7, 36) obtained for the same preparation as the MR data. Whole rabbit ventricle histology data are currently available at a resolution of 1 × 1 × 10 μm and provide information on different cell-type distributions as well as cardiomyocyte alignment (fiber orientation). Current work focuses on the full registration of 3-D histological tissue volumes with the original MR data to remove distortion artifacts introduced during histological processing (30). Such registered histological information will provide important information regarding the distribution of nonexcitable connective tissue within the ventricles as well as a detailed representation of 3-D cardiac fiber orientation at a subcellular resolution. This will support more accurate representation of the specific fiber architecture around fine-scale anatomical structures such as blood vessels.

This distribution of fibers around coronary vessels is known to be complex, and it has been suggested to play a key role in arrhythmias (16). Rule-based methods for assigning fiber directions, such as those used in this study, do not account for the presence of vessels. Specific rule-based strategies to represent the negotiation of fibers around blood vessels are currently being developed (16), and methods for incorporating these into realistic geometries are under way.

The rule-based methods used here play important roles in assigning fiber vectors within computational ventricular models, where DT-MRI data are not available (i.e., the vast majority of patient studies). Rule-based methods have been shown to agree well with DT-MRI measurements in wedge preparations (33) as well as in whole ventricular models (6) and have been shown to offer advantages in areas where DT-MRI measurements are typically noisy, for example, close to epicardial/endocardial surfaces (6). However, a major drawback of rule-based methods when applied to high-resolution ventricular models is their limited representation of anatomical structures that are not part of the myocardial wall/septum (e.g., papillary muscles and large trabeculations). The novel method used in this study to represent fiber architecture within these tube-like structures represents a step forward in the use of anatomical MR data to construct ventricular models in the absence of DT-MRI information (for example, for data obtained in vivo). However, for the assessment of inter-individual variability, access to DT-MRI data will remain highly desirable for further progression of this work, and coregistration mechanisms for multiple MR-based and histological datasets are being developed (36).

Assessment of the Importance of Including Fine-Scale Anatomical Structure in Ventricular Models

Following simple single-site pacing protocols, our results show that the global appearance of activation sequences was similar between the complex and simplified models. In particular, both the initial activation wavefront speeds and wavefront morphologies were similar within the LV and RV walls after apical stimuli. This reconfirms the utility of the generation of ventricular models with moderate anatomical detail for the application to electrophysiological research, for example, in current studies of cardiac resynchronization therapy, where only approximate global activation sequences through the cavity walls are generally required.

However, despite the overall similarities, a number of significant differences between activation dynamics within the two models were identified, which highlights the limits of applying conclusions from generalized simple architecture models to individual cases. As shown in our results, of primary importance is the presence of large trabeculation structures that link the ventricular walls to the septum (particularly found close to the apex). These trabeculae provide a “shortcut” pathway along which activation rapidly propagates, causing depolarization of the tissue in advance of the global propagating activation wavefront. These give rise to individual focal sources of excitation on the epicardium, leading to a less planar wavefront. During sinus rhythm, these trabeculations might very well have a physiological function, as the shortcut pathways may assist the natural activation sequence. The details of this role will be elucidated in the context of incorporation of an anatomically realistic representation of the Purkinje system into the rabbit ventricular model (work in progress) to simulate sinus rhythm.

Conversely, the additional trabecular pathways may be of importance for the initiation and/or sustenance of reentrant activity, for example, in situations where unidirectional conduction block within the bulk of the myocardial wall would otherwise terminate the reentry. In addition, the existence of focal sources of excitation, breaking through to the epicardial surface following conduction along trabecular structures (witnessed during apical pacing), may also have important proarrhythmic consequences. Disparities in the conduction velocities of waves propagating out from focal sources, as has been shown in optical imaging experiments (8, 25), could introduce a heterogeneous distribution of tissue refractoriness around the site of the focus, which in turn could lead to destabilization of wavefronts and wavebreak in subsequent reentrant cycles.

Here, significant differences between the complex and simplified models can be found in the response of the tissue to the shock. First, the epicardial (both posterior and anterior) surfaces were appreciably less strongly polarized by the shock in the complex model than in the simplified one. Second, the fine-scale anatomical detail present in the complex model led to a more localized heterogeneous distribution of polarization levels, both on the endocardial surface (surrounding the junctions of the endocardial structures with the myocardial wall) and, to a lesser extent, intramurally (surrounding blood vessels). Such heterogeneous potential distribution has been witnessed in experimental recordings from intact hearts using optrodes (26) as well as in optical mapping measurements from transmural wedge preparations after intermediate-strength shocks, which are less affected by the well-documented measurement artifacts during such recordings (35, 49). However, such effects had not been reproduced in the wealth of prior computational studies examining the effect of electrical shocks on whole ventricular models, due to the inherent lack of structural detail in those models. Our findings thus illustrate the utility of more complex models in certain scenarios.

It is widely known that the preshock state of the tissue is a key determinant of the response of the ventricles to the shock after pacing (43). In our simulations, although the differences in activation times and wavefront morphologies seen after apical pacing in Activation Sequences and Wavefront Morphology (due chiefly to wavefront propagation along trabeculae) did lead to small differences in the preshock state of certain areas of the ventricles, these differences were largely confined to the septum, most noticeably close to the base (195-ms images in Fig. 7). Importantly, the LV/RV free walls had very similar preshock states (195-ms images in Fig. 7), due to the similar activation times of these regions noted in Activation Sequences and Wavefront Morphology. The differences in the response of the ventricles to the shock, introduced through the inclusion of additional anatomical complexity in our model, were witnessed throughout the volume of the ventricles and were not confined only to regions with (minute) differences in preshock states. Therefore, the way in which fine-scale anatomical structures interact with external electrical stimulation, highlighted in this study, represents a novel finding that may have important implications for defibrillation, which requires further investigation.

Specifically, in this experiment, it appears that endocardial structures may “protect” the bulk myocardium from becoming as strongly polarized as seen in simplified models. Furthermore, our simulations suggest that the intramurally heterogeneous distributions of Vm witnessed here may provide a plausible explanation for the weaker polarization of the epicardial surfaces seen in the complex model relative to the simplified model. This phenomenon may also help to explain the previous disparity between the magnitude of shock end epicardial virtual electrode polarization seen experimentally and those predicted by simulations using simplified ventricular models (5, 43). Although these differences can be accounted for, in part at least, by the distortion effects of photon scattering in optical recordings (5) and the inclusion of electroporation currents (1), simulations with simplified models still overestimate experimentally recorded measurements of the shock-induced change in Vm. Thus, the results from the present study suggest that the use of more realistic models in this context could help to bridge the gap between simulation and experiment.

The highlighted differences in polarization states at shock end lead to discrepancies in activity after external electrical shock application. These include the presence of new excitation wavefronts in different locations between the models as well as global differences in the recovery rates of the apex and base regions from the stimulus. Consequently, these differences were amplified in the overall postshock response of the two models, with different ensuing episodes of arrhythmogeneis observed in each model. The specific arrhythmias induced in these simulations were neither sufficiently long nor complex enough for the analysis of the particular roles played by anatomical structures in their proliferation or maintenance. For example, phase singularities were not seen to be anchored to structures such as blood vessels or papillary muscle junctions, as has been suggested experimentally (24, 41, 58). However, the major goal of this study was to compare the behavior of complex and simple models after the application of identical shock-induced arrhythmia protocols. This was achieved, highlighting significant differences in the response of the models to external stimuli and subsequent arrhythmia induction. This further emphasizes the need to include similar levels of complexity in future studies in this area to understand the mechanisms and consequences of electrical field interaction with myocardial tissue, as well as the role played by heterogeneities in histoanatomical complexity between individuals. It is hoped that identifying this modeling requirement will encourage further research, which may be facilitated by open access to the model used in this study.

Study Limitations

Although downsampling of the original MR data was necessary to reduce the memory footprint of the dataset as well as to introduce a uniform smoothing, which improved the performance of the segmentation algorithms, it also inevitably introduced small distortions to geometrical structures. However, the high resolution of the images was due in part to zero filling of the raw images in the MR acquisition protocol, in fact upsampling the originally collected data, to visually improve an image without adding structural information. As a result, the downsampling, as demonstrated in materials and methods (Fig. 1D), did not cause significant loss of important structural information. However, the reduced definition of some of the finer structures (such as thin blood vessel walls on the epicardial surface) lead to problems in accurately delineating them via automated segmentation. This was addressed in the manual postprocessing stage of the segmentation process, as described in Segmentation postprocessing. Recent attempts have been made to adapt the segmentation algorithms to segment small regions of the full-resolution dataset. However, transferring this resolution over to the meshing stage involved the generation of finite-element meshes with a number of degrees of freedom that exceed current limits for practical application to computational modeling of whole ventricular structure and function. Although some recent studies have generated meshes of ventricular wedges or free wall preparations (40, 56) for use in simulations at slightly higher resolutions (up to ∼60 μm), this study provides a high-resolution whole ventricular computational model that can be used relatively easily in normally available cardiac simulation environments.

Finally, we wish to highlight additional features that have not been included in the model at this stage. First, isotropic conductivity transverse to the myocyte axis has been assumed. Although a widely held assumption (11), a recent study (20) has suggested that there is a discrepency between conductivity in the sheet and sheet-normal directions in the heart. However, due to lack of information on the specific sheet structure in the ventricles for this preparation, and the absence of corresponding rule-based methods for incorporation of such information into computational models, sheet architecture was not represented in the model. In future studies, it is hoped to derive sheet structure information directly from either histology slices previously coregistered with MR data (30) or from DT-MRI data (47) to represent the tissue as electrically orthotropic.

The Purkinje system is not represented in the model. To date, there is only one study (60) representing a Purkinje model formulation that can be used in conjunction with a ventricular bidomain model, but this was used in conjunction with a more simplified model geometry. Although it is possible to extract and represent computationally the free-running Purkinje system in our model, implementation of the junctions between Purkinje fibers and the myocardium and delineation of the endocardial-bound Purkinje system will require further histological and functional data. Although the presence of a Purkinje system changes activation sequences, which will alter some of the effects seen during single-site pacing in this study, the absence of fast electrical impulse propagation through the Purkinje system allows one to dissect the effects of trabeculae-mediated breakthroughs during pacing, which represents a novel finding. Furthermore, a direct comparison between the two models in the presence of a Purkinje system is difficult, as many Purkinje-ventricular junctions attach to finer endocardial structures, such as the papillary muscle, which are not represented in the simplified geometry.

Regarding the electrophysiological detail incorporated into the model, our aim was to match the current state of the art, (e.g. Refs. 2, 10, and 53, all of which were performed with anatomically simplified ventricular models). In line with these studies, our model also does not use a discontinuous approach to account for the discrete effects due to gap junctions, nor do we specify regional heterogeneity in cell membrane dynamics. We used the most recent Mahajan-Shiferaw model of the rabbit ventricular cell. This was tailored to reproduce rabbit cellular dynamics at fast heart rates and was augmented with two additional currents to accurately represent the effects of strong shocks on tissue. At this point the cell model has not been specialized to account for regional heterogeneity, but the framework is in place to do so when appropriate cell model modifications become available. Although the inclusion of higher levels of electrophysiological detail may affect the results found in this study, we believe that the fundamental findings presented here will hold true. For example, the trabeculae-mediated epicardial breakthrough is a robust finding that depends on the activation front of the wave of depolarization and so will be little affected by changes in the repolarization gradients introduced by spatial heterogeneity in action potential shape and duration during single pacing. Although heterogeneity in repolarization would alter the recovery isoline before the delivery of the external stimulus in Response of the Ventricles to Electrical Field Stimulation, the large differences in the effect of the shock on the complex model relative to the simplified model were seen in all regions of the ventricles, thus appearing to be independent of the exact timing of recovery of the tissue. Finally, electrophysiological heterogeneity is largely attenuated at cycle lengths as short as those used in this study (<300 ms) (28). Overall, we believe that despite these limitations, our study shows that incorporation of fine-scale anatomical complexity, such as blood vessels and endocardial structures, can provide insight into novel electrophysiological mechanisms that cannot be elucidated with previously used simplified models.

Conclusions

In this study, we have presented the full methodological pipeline involved in the construction of an anatomically complex computational model of the rabbit ventricles and made both the high-resolution MR data from which the model was constructed and the final finite-element computational model itself available to the research community. We have demonstrated that simple and complex models share a large extent of similarities in the functional modeling of cardiac electrical behavior. This is an important illustration of the utility of simple models for general basic cardiac research applications. However, significant differences exist between standard and high-resolution models in activation sequences and particularly in the response of tissue to external stimuli and the ensuing patterns, complexity, and duration of the induced arrhythmias. This indicates that for a range of applications, sample- or patient-specific modeling will require higher structural resolution than that currently established as a standard.

The previously established standard of computational ventricular models, particularly those developed from the University of California-San Diego rabbit (59) and Auckland canine (32) and pig (50) datasets, have given rise to advances in our knowledge of cardiac electrophysiology over the past decade or so, in particular thanks to their open access policy. Here, we provide an additional high-resolution MR-derived ventricular model of the rabbit to the community, describing in full detail the methods and techniques to allow other researchers to repeat the steps used for the model generation pipeline, and invite investigation of other anatomically influenced differences in micro- and macroscale behavior to be assessed in more detail over a broader range of pathophysiological scenarios.

APPENDIX

Segmentation Details

This section contains specific details of the level-set segmentation methods outlined in materials and methods (Segmentation pipeline), including parameter values where appropriate. For specific details of the nature of individual image processing filters, please refer to the ITK User Guide or the extensive online documentation (www.itk.org). It should be noted that all parameter values given below were optimized based on extensive testing and visual inspection of the resulting segmentations.

Threshold level-set filter.

As described briefly in Segmentation pipeline, the first stage of the segmentation pipeline involved the use of the threshold level-set filter from the ITK library. The threshold level-set filter takes as an input an initial level-set surface generated by the fast marching method (48). The ITK fast marching image filter requires the specification of initial seed points to compute a distance map. The initial selection of seed points throughout the image volume was chosen based on local intensity values. Specifically, this involved moving a 3 × 3 × 3-voxel kernel throughout the image volume and comparing both the mean intensity of all 27 voxels within the kernel as well as the intensity of the central voxel to user-defined limits (160 and 165, respectively). If both the mean and central voxel intensities were above these limits, a seed point was placed at the center of the kernel. The kernel was moved in steps of 2 voxels in each direction to avoid placing seed points within neighboring voxels. A user-provided initial distance (in this case, 3) was subtracted from the distance map to obtain a level set in which the zero set represented the initial contour. This initial level set was given as input to the threshold level-set filter along with the feature image itself (MR dataset).

The other user-defined limits in the threshold level-set filter were the upper (U) and lower (L) threshold limits, which govern the propagation of the evolving contour. The propagation term P in Eq. 1, at a given point in space with intensity level g, has a functional form rendering P positive for L < g < U and negative otherwise. The threshold limits used here were L = 180 and U = 230. Scaling parameters were then used to balance the relative influence of the propagation (inflation) and curvature (surface smoothing) terms from Eq. 1. Respective values of 1.0 and 5.0 were used for the parameters β and γ in Eq. 1. The advection term was not used in this filter. The threshold level-set filter was run for maximum of 2,000 iterations with an root mean square tolerance of 0.02.

Geodesic active contour level-set filter.

The output of the threshold level-set filter, described above, acted as an initial level set for the geodesic active contour level-set filter. This filter uses an advection term to attract the level set to object boundaries. In addition to the initial level set from the threshold filter, the geodesic filter takes as input the raw MR data stack. However, before its input, the raw MR image stack was firstly smoothed with a curvature anisotropic diffusion filter with a time step of 0.0625, number of iterations = 5, and conductance of 9.0. An edge potential image was then produced by passing the smoothed image through a gradient magnitude recursive Gaussian filter (using a σ of 0.05), followed by a sigmoid image filter with parameters of α = −1.0 and β = 2.0. The edge potential image along with the initial level set was then used directly in the geodesic active contours level-set filter. The relative weights of the propagation, curvature, and advection terms of the geodesic active contours level-set filter were set to −10, 1.0, and 1.0, respectively. The filter was run for a maximum of 2,000 iterations with a root mean square error of 0.001.

Laplacian level-set filter.

Finally, the output of the geodesic active contours level-set filter acted as an initial level-set input for the Laplacian level-set filter. In this filter, the propagation (or speed) term in Eq. 1 was constructed by applying a Laplacian image filter onto the raw MR image stack. The Laplacian filter calculates the second derivatives of the image and thus amplifies noise. To minimize this effect, we applied an initial smoothing anisotropic diffusion filter to the MR dataset. In this case, the gradient anisotropic diffusion filter was used with 10 iterations, a time step of 0.0625, and a conductance of 2.0. In addition to the smoothed raw MR image, the second input to the Laplacian level-set filter was an initial level set, which in this case was the direct result from the geodesic level-set filter above. The Laplacian level-set filter was then run with the curvature scaling term (in Eq. 1) set to 10.0 and the propagation scaling term set to 1.0 (the advection term was not used in this filter). Finally, the filter was run with a maximum root mean square error of 0.002 and the maximum number of iterations of 250.

Mesh Generation Details

Here, we provide details of the parameters taken as input by Tarantula in the generation of the finite-element mesh (Results of the mesh generation). Parameter values along with descriptions are shown in Table 1. For further details, please see www.meshing.at.

View this table:
Table 1.

Meshing parameters used in Tarantula

Surface Discrimination Algorithm

Identification of different surfaces within the final ventricular finite-element mesh was achieved by performing a secondary segmentation of the MR data using the same segmentation pipeline described in Segmentation of the Dataset. However, in this case, seed points were manually placed only within the two ventricular cavities and outside the volume of the heart (i.e., excluding the blood vessels and extracellular cleft spaces). In addition, after segmentation, the atria were left attached to the ventricles, and a “cap” was put on top of the aorta to ensure that the cavities were fully isolated from the outside space and from each other. For this particular application, the segmented voxel stack was meshed using the meshing software Tetgen (www.tetgen.berlios.de/) after prior computation of a STL surface representation using a marching cubes algorithm, using the methods described in Ref. 38. Tetgen is capable of producing tetrahedral finite-element meshes, where the elements corresponding to separate enclosed regions are automatically tagged with different numerical tags, a feature that is unfortunately not currently available in Tarantula. However, due to the reduced level of detail present in the segmented voxel stack, the results produced by Tetgen in this case were sufficiently robust.

Meshing of the resulting segmented dataset produced elements with three separate numerical tags representing the myocardial tissue volume, the LV cavity, and the RV cavity, as shown in Fig. 9A. Note that in this case the ventricular cavities were directly connected to the atrial cavities as well. Examination of the numerical tags of each pair of elements bordering each bounding triangle face allowed discrimination between the three individual surfaces to be made. For example, a triangle that is part of the LV (or RV) endocardium also forms part of two tetrahedral elements: one in the myocardium (green) and one in the LV cavity (red) [or RV cavity (blue)]. A triangle that is part of the epicardium only forms part of one tetrahedral element. A schematic diagram demonstrating this identification process is shown in Fig. 9B.

Fig. 9.

Identification of discrete surfaces. A: whole cardiac mesh [including atrial tissue] produced from Tetgen showing the tagged tetrahedral elements: myocardium (green), LV/right atrial (RA) cavities (red), and RV/RA cavities (blue). B: schematic diagram showing the identification of different surfaces based on the tags of bordering elements. The example shown is of an LV endocardial face triangle. C: visualization of the tagged surface nodes in the final ventricular finite-element mesh, representing the epicardium (red), LV endocardium (white), and RV endocardium (yellow).

A list of all node points (and their spatial coordinates) was then formed for the epicardial, LV endocardial, and RV endocardial surfaces in the secondary segmentation. The spatial coordinates of these surface node points were then mapped onto the final ventricular mesh (shown in Fig. 4) using a nearest neighbor mapping algorithm to tag nodes in the final ventricular mesh as being one of the three surfaces. To avoid mapping nodes associated with the atrial tissue, all nodes with x,y,z positions residing above the plane used to remove the atria in Segmentation postprocessing were removed (Fig. 3B). Figure 9C shows the surface nodes with the epicardium highlighted in red, the LV endocardium in white, and the RV endocardium in yellow.

GRANTS

This work was supported by Biotechnology and Biological Sciences Research Council Grant ∼BB/F004834/1 for support of the Oxford 3D Heart Project. G. Plank was the receipient of Austrian Science Fund (Fonds zur Förderungder Wissenschaftlichen Forschung) Grant F3210-N18 and a Distributed European Infrastructure for Supercomputing Applications Extreme Computing Initiative grant (muHeart).

DISCLOSURES

No conflicts of interest are declared by the author(s).

ACKNOWLEDGMENTS

The authors thank Dr. Christian Bollensdorf, Dr. Blanca Rodriguez, Michal Plotkowiak, and Tahir Mansoori as well as Dr. Susana Garcia for the access to Wacom hardware (funded by the Technical Computing Initiative of Microsoft Corporation) and Dr. Edward Vigmond for the use of Meshalyzer software. M. J. Bishop is a Sir Henry Wellcome Postdoctoral Fellow, V. Ggrau is a Research Councils UK Academic Fellow, and P. Kohl is a British Heart Foundation Senior Research Fellow.

Footnotes

  • 1Supplemental Material for this article is available online at the American Journal of Physiology-Heart and Circulatory Physiology website.

REFERENCES

  1. 1.
  2. 2.
  3. 3.
  4. 5.
  5. 6.
  6. 7.
  7. 8.
  8. 9.
  9. 10.
  10. 11.
  11. 12.
  12. 13.
  13. 14.
  14. 15.
  15. 16.
  16. 17.
  17. 18.
  18. 19.
  19. 20.
  20. 21.
  21. 22.
  22. 23.
  23. 24.
  24. 25.
  25. 26.
  26. 27.
  27. 28.
  28. 29.
  29. 30.
  30. 31.
  31. 32.
  32. 33.
  33. 34.
  34. 35.
  35. 36.
  36. 37.
  37. 38.
  38. 39.
  39. 40.
  40. 41.
  41. 42.
  42. 43.
  43. 44.
  44. 45.
  45. 46.
  46. 47.
  47. 48.
  48. 49.
  49. 50.
  50. 51.
  51. 52.
  52. 53.
  53. 54.
  54. 55.
  55. 56.
  56. 57.
  57. 58.
  58. 59.
  59. 60.
  60. 61.
View Abstract