Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2837787

Formats

Article sections

Authors

Related links

Biomech Model Mechanobiol. Author manuscript; available in PMC 2011 April 1.

Published in final edited form as:

Published online 2009 September 2. doi: 10.1007/s10237-009-0170-5

PMCID: PMC2837787

NIHMSID: NIHMS159333

J.P. Carson and A.P Kuprat

Biological Monitoring & Modeling, Pacific Northwest National Laboratory, Richland, WA. Email: vog.lnp@nosrac.semaj

X. Jiao and V. Dyedov

Department of Applied Mathematics & Statistics, Stony Brook University, Stony Brook, NY. Email: ude.bsynus.sma@oaij, Email: ude.bsynus.sma@rimidalv

J.M. Guccione and M.B. Ratcliffe

UCSF Dept. of Surgery & San Francisco VA Medical Center, San Francisco, CA. Email: vog.av@effilctar.kram, Email: vog.av@enoiccug.suiluj

Biological Monitoring & Modeling, Pacific Northwest National Laboratory, Richland, WA. Email: vog.lnp@nietsnie.leinad

The publisher's final edited version of this article is available at Biomech Model Mechanobiol

See other articles in PMC that cite the published article.

Spatial discretization of complex imaging-derived fluid-solid geometries, such as the cardiac environment, is a critical but often overlooked challenge in biomechanical computations. This is particularly true in problems with Lagrangian interfaces, where the fluid and solid phases share a common interface geometrically. For simplicity and better accuracy, it is also highly desirable for the two phases to have a matching surface mesh at the interface between them. We outline a method for solving this problem, and illustrate the approach with a 3D fluid-solid mesh of the mouse heart. An MRI perfusion-fixed dataset of a mouse heart with 50μm isotropic resolution was semi-automatically segmented using a customized multimaterial connected-threshold approach that divided the volume into non-overlapping regions of blood, tissue and background. Subsequently, a multimaterial marching cubes algorithm was applied to the segmented data to produce two detailed, compatible isosurfaces, one for blood and one for tissue. Both isosurfaces were simultaneously smoothed with a multimaterial smoothing algorithm that exactly conserves the volume for each phase. Using these two isosurfaces, we developed and applied novel automated meshing algorithms to generate anisotropic hybrid meshes on arbitrary biological geometries with the number of layers and the desired element anisotropy for each phase as the only input parameters. Since our meshes adapt to the local feature sizes and include boundary layer prisms, they are more efficient and accurate than non-adaptive, isotropic meshes, and the fluid-structure interaction computations will tend to have relative error equilibrated over the whole mesh.

The heart and its immediate environs are intimately linked with the movement of blood. From a mechanical point of view, for example, the resistance to blood flow through the cardiac valves feeds back into the ventricular afterload and contractility, and disturbances in the flow can lead to thrombus formation and coronary obstruction. On a more biological scale, shear stress and flow conditions at the interface are intimately linked to cell life/death and indeed gene expression. For example, it is well known that wall shear stresses regulate coronary endothelium and that inappropriate wall shear stresses can lead to build up of atherosclerotic plaque, e.g. [5,6,9,10,13]. In a similar vein, but less well known, wall shear stresses have been shown to regulate the phenotypes of valve endothelium and by extension valve interstitial cells on the ventricular and aortic sides of the aortic valve [2–4]. In fact, the connection between heart function and the movement and properties of blood is so intimate that the normal morphogenesis of the heart is dependent on proper blood dynamics [3,14,29,35], to the point where with interrupted flow the embryonic heart may develop three chambers rather than four [14].

Because the interface is so important from a biological point of view, it is highly desirable to treat the interface in a Lagrangian framework, where fluid shear stresses can be resolved more accurately. However, this poses two fundamental challenges, the first of which is that the discretization of each phase must be consistent at the interface despite each phase being driven by different physical considerations. For example, the boundary layer is a consideration on the fluid side of the interface but not on the solid side. In addition, the solid side may require only a few layers as a minimum accommodation for bending or lamination, while the fluid phase requires a dedicated grid for resolving the boundary layer, and a compatible fluid interior grid. Due to the complexity of the geometry, all algorithms must be fully automatic with few meaningful user-selectable parameters. The second fundamental challenge is that the grid in both fluid and solid phases may evolve separately. A canonical example is a closing — and opening — heart valve around which the fluid grid must change topology over the course of a simulation, while the grid requirements of the valve tissue are largely unchanged. Here we focus on the first of these challenges. We have addressed, in part, the second of these challenges elsewhere [12]; though work in that area is still ongoing. In short, the strategy is to add an arbitrary Lagrangian-Eulerian (ALE) term to the the Particle Finite Element Method (PFEM) equations of motion (cite ), while preserving a Lagrangian interface by efficiently updating the mesh to avoid overly distorted or squeezed elements. Thus, rather than stopping the simulation to perform a remesh, the mesh is managed dynamically throughout the course of the simulation by invoking the usual operations of nodal smoothing and augementing them with nodal point insertion and deletion as well as edge bisection to allow for topological changes in the mesh.

In this manuscript, we describe the underpinnings of our methods, namely multimaterial segmentation, multimaterial isosurface extraction, generation of the gradient-limited feature size, construction of layered tissue tetrahedra, generation of a prismatic boundary layer, and generation of a boundary constrained Delaunay fluid interior mesh. We propose new methods for generating quality, hybrid, Lagrangian fluid-solid grids for complex biomedical geometries from imaging data. The Lagrangian fluid-solid grid has its own set of challenges insofar as a fluid-solid grid is multi-material, consisting of — at least — fluid and solid phases. In addition, for simplicity and for reducing numerical error over the course of a simulation, we require that the fluid and solid boundaries match exactly along their common interface. These challenges are present at every stage, including image segmentation, surface extraction, cleanup and grid generation. Our mesh-generation methods are by design automatic, meaning that they require no user intervention. They are also efficient in the sense that, in most cases, a computable grid may be obtained in few hours, given a proper segmentation. We focus on a perfusion fixed dataset of a mouse heart, though the methods are applicable to other biomedical geometries. Some of these ideas have been developed elsewhere, and are thus described in more general terms. Novel ideas presented herein are described in detail with pseudo code.

We describe our methods starting from a brief description of the acquisition of MRI data. Following that we describe our approach to multimaterial segmentation of the heart tissue and blood. We then describe our approach to multimaterial isosurface extraction. Finally, we highlight our newest contributions to fluid-solid grid generation.

An adult male C57BL6 mouse commercially obtained (Jackson Laboratory, Bar Harbor, ME) was actively stained and imaged as previously described in publication [21]. A 1:20 mixture of gadopentetate dimeglumine (Magnevist; Berlex Laboratories, Wayne, NJ) and10% buffered formalin (Fisher Scientific, Fairlawn, NJ) were systemically perfused into the mouse through the left carotid artery and the right jugular vein. A 3D MR image (512 × 512 × 2048) was acquired at 50μ*m* isotropic resolution with a modified MR apparatus (Signa; GE Medical Systems, Milwaukee, WI). Image reconstruction (Figure 1) was performed on a dedicated 64-bit reconstruction engine (Silicon Graphics, Mountain View, CA). A 3D image subset (249 × 216 × 240) encompassing the heart and nearby tissue was cropped from the whole volume as the focus of this work.

Voxel assignment to different materials for the 3D MRI dataset was performed using a semi-automated segmentation procedure consisting of four steps. Software for performing segmentation routines was written in Python (www.python.org) with the graphical user interface (GUI) implemented using Visual Python (www.vpython.org). The overall segmentation process is articulated in four steps.

Every voxel was first automatically assigned to one of four material types: the perfused blood material (*M _{b}*), the heart tissue material (

Histogram of the pixel intensity across all images in the volume. Each voxel was assigned to four different material types based on empirically determined thresholds.

This step begins by automatically reassigning fully encapsulated disputed material *M _{d}* voxel groups to the surrounding material when the enclosing material is all of the same material type (blood

Next, disputed material *M _{d}* voxels representing coronary vasculature were automatically reassigned to blood material

To complete the reassignment of *M _{d}* voxels to

The application of Step 2 is represented for a sample image and illustrated in Figure 3. The procedure is detailed in Algorithm 2.

Semi-automatic fuzzy determination of disputed voxels. The tissue is shown in dark grey; blood is shown in white, and disputed voxels are shown in light grey. The process of resolving the fuzzy set is shown from left to right.

The non-cardiovascular tissues adjacent to the heart and vessel walls often appear at a similar intensity in the MRI images, thus making them distinguishable only by incorporating anatomical knowledge. A semi-automated masking process was implemented to facilitate the isolation of heart and blood from the remainder of the image (Algorithm 3). First, a user-generated mask indicating the cardiovascular region for the first image in the volume was created. Then, for each subsequent image, a new mask was generated from its preceding image by applying a warp function to the preceding mask. This warp function represents a transformation that deforms an image to match its subsequent image in the volume. The warp function was calculated using multiresolution elastic image registration [23] on two adjacent images in the volume, and then the function was applied to the mask of the first image to generate the mask of the second image. After generating a new mask, a minimal amount of manual augmentation to the mask was performed to address new features emerging from one image to the next.

Several topological criteria were iteratively enforced to ensure the segmentation would generate an anatomically accurate multimaterial mesh (Algorithm 4). First, blood material was forbidden from contacting any material other than heart material. This was accomplished by replacing all *M _{o}* voxels 26-adjacency connected to

Multimaterial imaging-based geometries, such as that represented by the mouse heart shown in Figure 4, represent clear challenges for generating any kind of quality computational grid, let alone one for Lagrangian FSI, where we require exact node matching at the interface. Our approach to this challenge is based on the gradient-limited feature size (GLFS) [25]. The GLFS serves three functions in our grid generation algorithms: 1) a field for tangential adaptation of the surface grid [25], 2) a metric for creating layered tetrahedra [25], and 3) a speed function for construction of a prismatic boundary layer [11] by application of the generalized Huygens' principle [17]. Details of these methods are beyond the present scope and were presented elsewhere [25,11]. Herein, we review some of the key points of those references as background to the novel refinements and enhancements that we will describe shortly. These novel contributions can be summarized as:

- a modification of the functional relationship between the GLFS and the first principal curvature
- a further adjustment of the GLFS to resolve thin tube-like structures of the heart
- the method by which points are cast along normal rays to create quality layered tetrahedra
- an explicit time-dependent computation of the GLFS during interior surface evolution for the generation of prismatic boundary layer elements

Multimaterial segmentation of the mouse heart. Initial example MR slice is shown in (A) with the cropped region highlighted by the box and shown segmented into blood (blue), tissue (red), and other (white) in (B). Various renderings of the segmentation **...**

Following segmentation and topology validation, surfaces were extracted from the segmented mask with a variant of the marching cubes algorithm [34]. Let *S* be the oriented closed triangulated surface or network of surfaces, obtained in Section 2.2.4. For any point **x** *S*, we define the “raw feature size” or “local diameter” *F*[**x**] as the length of the line segment formed by first shooting a ray from **x** in the direction of [**x**], the inward normal at **x**, and then truncating the ray at its first intersection with *S*. That is

$$F\left[\mathbf{x}\right]\equiv \mathit{min}\left\{\lambda >0\mid \mathbf{x}+\lambda \widehat{\mathbf{n}}\left[\mathbf{x}\right]\in S\right\}$$

(1)

Since *S* is closed, the ray proceeding from **x** in the surface normal direction [**x**] will intersect *S* at least once, and hence *F*[**x**] is well-defined. Similarly, we also perform an “outward” interrogation of the geometry to compute another raw feature size field *F _{out}* [

In addition, *F _{out}* [

Once *F*[**x**] has been modified to reflect this additional geometric information, it is gradient limited to render it relatively insensitive to abrupt changes in *S* [25]. First, we initialize *f* [**x**] to *F*[**x**]. Given a bound *G* on the surface gradient of *f* [**x**], the algorithm places the directed edges that violate the gradient limit into a maxpriority queue, ranked by the key

$$f\left[{x}_{1}\right]-(f\left[{x}_{2}\right]+G\parallel {x}_{1}-{x}_{2}\parallel )$$

(2)

which measures how much the gradient violates the gradient limit for a directed edge $\overrightarrow{{x}_{1}{x}_{2}}$ on *S*. Let $\overrightarrow{{x}_{i}{x}_{j}}$ be the directed edge with the largest measure in the queue. We relax *f* [**x**] to satisfy the gradient limit, recompute the measures for the incident edges of *x _{i}*, and update the priority queue accordingly. The process continues until the queue is empty. For computational efficiency, ray-triangle intersections are queried within an axis-aligned bounding box (AABB) tree [22,32] that contains at its leaf nodes the bounding box for each triangle. This algorithm has a complexity of (

GLFS computed on tissue (A) and interior blood (B) surfaces. Note the different dynamic ranges. The coronary vasculature apparent in (B) is implicit in (A) where surface rays detect the involuted boundary of the arteries. Computation of *f* [*x*] for these **...**

Having defined *f* [**x**] and *f _{out}* [

Imaging based isosurfaces, such as those generated in Section 2.2.4, are highly irregular, noisy, and often over-resolved in some areas and conversely under-resolved in some others. To generate a volume mesh based on those isosurfaces would be wasteful of computational resources at best. At worse, error in the computed fields would be concentrated at the smaller scale, for example, in the coronary vasculature or around the valves. To address this, we modify *S* to produce a high-quality surface mesh *S*′ by performing the operations of smoothing, refinement and de-refinement while limiting perturbations to a small fraction of a voxel. To assure scale-invariance, the operations are registered to the GLFS [25]. Specifically, we perform the operations of (1) Rivara edge bisection, (2) node merging, (3) edge swaps, (4) in order to create surface triangles whose edges are close to *c _{t}f* [

Instead, we focus on a necessary further modification of the GLFS in the presence of tube-like structures such as the chordae tendinea and the coronary arteries. In these areas, modulating the GLFS with the first principle curvature, as described in Algorithm 5 and Equation 2, is not sufficient to guarantee a minimum number of elements circumferentially. This is particularly true in instances where these features are close to the resolution of the voxel.

For these cases, we introduce a novel curvature correction to *f* [**x**], such that the surface edge-length is sufficiently small to accommodate an angular spacing of NRING nodes around the circumference of the tube. This implies a surface edge length of

$$\frac{2\pi \times \mathit{tube}\phantom{\rule{thinmathspace}{0ex}}\mathit{radius}}{\mathit{NRING}}$$

(3)

This new edge-length may be less than the previous 'ray-shoot' edge length of ${\scriptstyle \frac{K}{M}}f\left[\mathbf{x}\right]$, where *K* is desired volume-element anisotropy (defined as the ratio of edge-length to volume element depth), and *M* is the desired number of layers. This is especially true where the tube diameter is a single voxel wide, or equal to *L _{min}*. Thus, we dictate that the surface edge length should be as small as ${\scriptstyle \frac{2\pi}{{\kappa}_{1}\mathrm{NRING}}}$. This is done by reducing

$$f\left[\mathbf{x}\right]\leftarrow \mathit{min}\left({f}_{1}\left[\mathbf{x}\right],\frac{K}{M}\frac{2\pi}{{\kappa}_{1}\mathit{NRING}}\right)$$

(4)

where *f*_{1}[**x**] has been previously computed. We again perform a gradient-limiting operation on the resulting GLFS field (Equation 2) and then employ surface mesh manipulation operations to establish an edge length of ${c}_{t}f\left[x\right]\equiv \frac{K}{M}f\left[x\right]$.

*Remark 1* In this section we have specified a functional value for *c _{t}*, and the reader may wonder why this was not defined earlier. The reason is simply that

As observed before, cardiac tissue is laminated. Thus, we desire the computational grid of the cardiac tissue to be similarly laminated. We achieve this by creating layers of tetrahedral elements on the interior of *S*′. Points are cast along “seeding rays” from each point **x*** _{i}* on the surface

If *M* is the target number of layers across the cross-section of the geometry, then points ${\mathbf{x}}_{i}^{m},0\le m\le {\scriptstyle \frac{M}{2}}$, are equally or ratio spaced between ${\mathbf{x}}_{i}^{0}\equiv {\mathbf{x}}_{i}$ on the surface to ${\mathbf{x}}_{i}+{\scriptstyle \frac{1}{2}}f\left[{\mathbf{x}}_{i}\right]\widehat{n}\left[{\mathbf{x}}_{i}\right]$. Due to gradient-limiting, we may have *f* [**x*** _{i}*] <

Having generated a quality layered tetrahedral grid for the heart tissue, we wish to create several layers of prisms at the Lagrangian interface between tissue and blood. There are several reasons for this. First, prismatic boundary layers are desirable for viscous flows, which tend to have strong gradients at the boundary due to wall shear stress, an effect that – as we have already mentioned – has been shown to be extremely important in cardiac/cardiovascular biomechanics. In addition, prismatic layers have been shown to be effective in reducing errors in turbulence, particle deposition, and mass-transfer problems [33,27,7,24]. In biological fluid-structure interaction problems with a Lagrangian interface, prismatic boundary layers are also desirable near the interfaces, such as closing cardiac valves, because the direction of impaction is locally normal to the boundary, so it would be aligned with axis of the prisms. This will not only provide more accurate results (since the larger gradients will be aligned with the mesh) but it is a desirable mesh property for any mesh movement algorithm in the framework of Arbitrary Lagrangian-Eulerian (ALE) formulations.

Our method advances a surface layer by solving the Lagrangian evolution equation,

$$\frac{\partial \mathbf{x}}{\partial t}=\phi (\mathbf{x},t)\widehat{\mathbf{n}}$$

(5)

where *t* denotes time, denotes the normal velocity, and denotes unit surface normal. Generating a layer of prisms reduces to propagating the surface for some time interval by discretizing Equation (5). Traditional Lagrangian methods are prone to the problem of “swallowtails” in strongly concave regions with large curvatures. To avoid such issues, we apply the *face offsetting method* in [17], which is based on a geometric construction called the *generalized Huygens' principle*. Unlike traditional Lagrangian methods, which propagate vertices along some vertex normals, the face-offsetting method first propagates faces and then reconstructs the vertices using numerical techniques of least-squares approximation and eigenvalue analysis. The method also adapts the time step to prevent mesh folding. To achieve mesh quality of the prisms, we also developed a novel variational smoothing procedure that simultaneously improves the shapes of the base triangles and the orthogonality of the side edges, as described in [11].

We have described our general approach of prismatic boundary layer meshes in [11]. In our initial implementation, we used the GLFS of the surface as the speed function in (5), and propagated the surface by a desired percentage of the GLFS. For more complex geometries like the mouse heart, however, it is important for the speed function to adapt to the changing curvatures as the surface evolves. In the present paper, we enhanced our implementation to recompute the GLFS of the evolved surface periodically. This enhancement allows areas with large curvatures to slow down gradually, while allowing the remainder of the surface to propagate as far as possible. This enhancement significantly improved the robustness of our algorithm, and allowed the boundary-layer mesh to cover about 30% of the volume for the blood region of the mouse heart, nearly twice as much as our original method. The reevaluation of the GLFS introduces only moderate overhead to our method. For the mouse heart mesh, each reevaluation requires about 3.5 seconds on a moderate computer. After generating the prismatic boundary layer meshes, we tetrahedralize the interior with a boundary constrained Delaunay method implemented in the code TetGen [31].

Figure 4 illustrates the highly detailed segmentation of tissue and blood, with the segmented blood volumes clearly separated into right (deoxygenated) and left (oxygenated) volumes. Figure 5 depicts the gradient limited feature sized field on exterior and interior surfaces of the reconstructed mouse heart. Note the detection of the coronary vasculature on the exterior surface (Figure 5A) and its revealed structure on the interior heart (Figure 5B). Figure 6 shows the generated layered tetrahedral mesh, which has 15,571,963 elements. Five prism layers, with 6,882,265 prisms in total, are shown in Figures 7 and and8.8. The Delaunay blood mesh, with 5,735,350 tetrahedra, is shown in Figure 8.

Layered tetrahedral mesh of the cardiac tissue. (A) Apex-to-base cut-away of the entire mouse heart. (B) Left ventricular chords. Note the smooth transition of scales from chords to the myocardium. (C) Detail of the left atrium. (D) Tricuspid valve and **...**

Five layers of boundary prisms, propagated by a dynamically updated GLFS field. Note the adaption to scale and the preservation of orthogonality.

Prisms with enclosed Delaunay tetrahedra. Panel A shows a number of coronary arteries traveling through the myocardium (empty space in picture). The arrow indicates the termination of a coronary vessel. Since the capillaries are not resolved, correct **...**

Mesh quality statistics for both tetrahedra and prisms are detailed in Figure 9. For tetrahedra, we report the *aspect ratio* which is proportional to the ratio of the inscribed radius to the length of the longest edge. For prisms, we report the socalled scaled aspect ratio (see [11]), whose definition is somewhat more involved. In effect, the scaled aspect ratio combines the measures of triangle shapes and edge orthogonality. Both quality metrics vary between 0 and 1, where 1 is optimal. The mesh quality for all three regions is excellent. Boundary prisms are layered by design, and layered tetrahedra resolve tissue lamination at all scales, from the ventricular septum to the coronary vasculature.

The role of computational continuum physics in cardiac research has expanded beyond an engineering analysis of devices or even idealized organs, to become a platform for discovery and integration. This thrust has been paralleled by a major research focus on biomedical image processing (www.itk.org, www.slicer.org, mipav.cit.nih.gov, rsbweb.nih.gov/ij). Not surprisingly, biomedical geometries are increasingly derived from imaging modalities such as magnetic resonance or computed tomography. Not only can a connection to imaging data provide critical geometric information — atherosclerosis for example is arguably a disease with a strong geometric component — but also critical physiologic data, such as diffusion anisotropy, genetic expression or even metabolism. Biomedical geometries, however, are complex and hierarchical. Despite this, little attention in the biomedical literature has been paid to the enormous challenges posed by geometric computation of biomedical geometries.

Efficient unstructured mesh generation of these geometries in a way that is tuned to the physics of the problem is challenging. Finite computational resources dictate that the grid must be adapted to the geometry in order to minimize both computational cost and discretization error. At the same time, many biophysics problems require that the grid be organized both by scale and by intrinsic properties. The wall of the heart, for example, is a laminated muscle consisting of three separate layers, each with a separate family of fiber and myocyte orientations. Thus, a grid of the heart muscle must be similarly layered. That same layering persists at the smaller scale of the heart geometry, where, for example the musculature transitions into the connective tissue of the cardiac valves; or at an even smaller scale where the heart muscle involutes to become a network of coronary arteries. These transitions of scale are of course mirrored in the blood. The lumen of a coronary artery, for example, may be several orders of magnitude smaller than the span across a ventricle, and thus there is a need to equilibrate error over a range of meaningful scales. These challenges of course are not limited to the heart.

However, the heart does pose some unique challenges. First, the effects of blood flow on cardiac mechanics and physiology are profound. From a biological viewpoint the interface is the zone of greatest activity and communication; from a medical point of view it is the site of the greatest vulnerability to disease. Thus, biomedical fluid-structure problems must aim to resolve continuum fields accurately at the interface. To address these challenges, we have presented the numerical framework to accomplish a Lagrangian treatment of tissue-blood interfaces [12]. Work in that area continues. Such an approach is greatly simplified if exact nodal compatibility exists between each region of the grid. Secondly, unlike arteries or veins in isolation, most of the heart structures cannot be approximated by a plane-stress membrane assumption. This implies the explicit discretization with volume elements such as tetrahedra, hexahedra and prisms. Quality discretization of complex geometries with hexahedra remains forbidding.

We have presented a method for computing high-quality scale-invariant grids from such imaging data, and have exemplified the approach with a demanding, perfusion fixed micro MRI dataset of the *in situ* mouse heart and its structure. The mesh quality achieved for this geometry is excellent by any measure. This quality owes to three factors: a careful and topologically correct segmentation to produce initial surfaces, the adoption of the gradient-limited feature size field, and variational optimization of marched prisms based on the face-offsetting method.

The concept of scale is directly embedded into all aspects of our approach by efficient computation of the gradient-limited feature size field. Were edge lengths based solely on some geometry-independent values, the number of elements required would increase by two or three orders of magnitude to capture the finest structures. In contrast, adoption of the GLFS field enables the efficient placement of mesh only where it is needed. Further savings are achieved by decoupling the mesh densities in the surface tangential and normal directions, though tangential density is locally isotropic. Adoption of the GLFS also enables adequate circumferential discretization of tubular regions of the geometry (Equation 4), and a finite boundary layer to resolve the stronger gradients that are found at the wall.

To explicate our method, we have presented detailed pseudo-code for novel aspects of our methods. Further details may be found in earlier references [11, 25]. Our methods are applicable to a wide variety of multimaterial data sets and are particularly well-suited to segmentation of perfused tissues. The objective here is to enable others to reproduce our results on related data sets.

Despite the success of our approach, there remain at least four areas of potentially substantial improvement. First, the grid density tangential to the surface is isotropic. On structures such as the cordae or coronary arteries, where fluid flow is largely axial, it would be considerably more efficient to introduce anisotropy into the grid, such that the greater density were arranged in the direction with the highest curvature. Such an approach has been successfully applied to vascular trees in the past [30,16]. A promising approach that is compatible with our current methods is described in [18]. A second area of improvement would be to relax the requirement of constrained Delaunay tetrahedralization for the blood core. Discretization of this zone can be considered preliminary since adaption in our fluid solver is driven by an error metric, and thus discretization will evolve. Nevertheless, the initial quality of this interior grid is important. Although the current approach will typically generate high-quality hybrid meshes for accurate and efficient numerical simulations, the use of constrained Delaunay tetrahedralization can limit the quality of the tetrahedral core for some complex biological geometries. Thirdly, it must be noted that while the overall method does not depend on “good” surfaces, it does depend on “correct” surfaces. In most cases, surface adaptation and volume conserving smoothing improve the input surfaces sufficiently. Application of the latter is critical for mesh improvement, and has the advantage over conventional smoothing approaches that the volume conservation constraint limits drift. However, our smoothing algorithm does not guarantee against self-intersection, and therefore needs to be applied with care. With complex geometries, a robust volume-conserving smoothing with theoretical guarantees against self-intersection would be an important asset. Finally, prisms have inherent advantages with respect to tetrahedra in terms of orthogonality and their avoidance of numerical pitfalls, such as locking^{1}. Additionally, prisms can potentially reduce the element count in a given zone by a factor of 2–3. Thus, the face-offsetting method has additional potential as a discretization tool, if select nodes can be merged, collapsing prisms into a small number of pyramids and tetrahedra. Addressing these four areas of improvement is part of our future work.

Though computational geometry has received little attention in the biomedical literature, there are some notable exceptions. A number of references to discretization of biomedical geometries can be found in the numerical methods literature. We cite only a few here that have taken a comprehensive view of biomedical geometries that spans imaging and geometric analysis. For example, a method described by Labelle for “isosurface stuffing” [26], fills an isosurface with tetrahedral mesh with excellent quality statistics. Though this method has not been applied to biomedical geometries, that we are aware of, it is certainly applicable. While this approach may be appropriate for some applications, it does not organize the grid by layer or scale and is generally unsuitable for fluid-flow simulations. A somewhat related approach that produces multiple material lattice-like tetrahedra on biological domains has been presented in [37]. This approach is similarly suited to a range of applications that fall outside cardiac fluid-structure interactions because the hexahedra are not aligned with the surface and there is no recognition of scale. A promising alternative is the isogeometric approach pioneered by Hughes et al. [15], and adapted for fluid-structure analysis of treed biomedical geometries by Zhang et al. [36]. This method has the dual advantages of discretizing volumes with quality hexahedra and being fully integrated with the numerical solver. However, it is unsuited to domains that cannot be reduced to a one-dimensional network of curves. Thus, it is well suited to arterial trees, but unsuited to the entire cardiac geometry.

In contrast, our method uses a gradient-limited feature size (GLFS) to guide the adaptation of the surface mesh, the generation of layered tetrahedra, the generation of the prismatic mesh at the boundary layer and the generation of tetrahedral mesh in the interior. This achieves scale-invariance such that even the smallest resolved coronary vessel is guaranteed to have a preset number of layers across its diameter. Furthermore, guided by the GLFS, we are able to laminate the cardiac tissue, seamlessly transitioning from myocaridum, to cardiac valve to coronary lumen. In addition, this approach enables the automatic adaption of the boundary layer by scale. Our method is also automatic in the sense that the user must only specify a few meaningful parameters: number of layers, degree of anisotropy, growth ratio for the prisms, and minimum and maximum resolved feature sizes.

Research funded by the National Heart Lung and Blood Institute Awards 1RO1HL073598-01A1 and 1R01HL084431-01A1, by the National Science Foundation under award number DMS-0809285, and by the U.S Department of Energy under LDRD DE-AC05-76RL01830.

^{1}Modern linear tetraheda have overcome historic issues with volumetric and shear locking in the large-deformation, nearly incompressible limit, e.g. see [8]

1. Amenta N, Bern M. Surface reconstruction by Voronoi filtering. Discr. and Comp. Geom. 1999;22:481–504.

2. Butcher JT, Markwald RR. Valvulogenesis: the moving target. Philos Trans R Soc Lond B Biol Sci. 2007;362(1484):1489–503. [PMC free article] [PubMed]

3. Butcher JT, Nerem RM. Valvular endothelial cells and the mechanoregulation of valvular pathology. Philos Trans R Soc Lond B Biol Sci. 2007;362(1484):1445–57. [PMC free article] [PubMed]

4. Butcher JT, Tressel S, Johnson T, Turner D, Sorescu G, Jo H, Nerem RM. Transcriptional profiles of valvular and vascular endothelial cells reveal phenotypic differences: influence of shear stress. Arterioscler Thromb Vasc Biol. 2006;26(1):69–77. [PubMed]

5. Chen LC, Nadziejko C. Effects of subchronic exposures to concentrated ambient particles (caps) in mice. v. caps exacerbate aortic plaque development in hyperlipidemic mice. Inhal Toxicol. 2005;17(4–5):217–24. [PubMed]

6. Chien S, Li S, Shiu YT, Li YS. Molecular basis of mechanical modulation of endothelial cell migration. Front Biosci. 2005;10:1985–2000. [PubMed]

7. Chiou MC. Particle deposition from natural convection boundary layer flow onto an isothermal vertical cylinder. Acta Mechanica. 1998;129:1619–6937.

8. Cisloiu R, Lovell M, Wang J. Astabilized mixed formulation for finite strain deformation for low-order tetrahedral solid elements. Finite Elements in Analysis and Design. 2008;44:472–482.

9. Cummins PM, Cotter EJ, Cahill PA. Hemodynamic regulation of metallopeptidases within the vasculature. Protein Pept Lett. 2004;11(5):433–42. [PubMed]

10. Dardik A, Yamashita A, Aziz F, Asada H, Sumpio BE. Shear stress-stimulated endothelial cells induce smooth muscle cell chemotaxis via platelet-derived growth factor-bb and interleukin-1alpha. J Vasc Surg. 2005;41(2):321–31. [PubMed]

11. Dyedov V, Einstein DR, Jiao X, Kuprat AP, Carson JP, del Pin F. Variational generation of prismatic boundary-layer meshes. Int. J. Numer. Meth. Engrg. 2009 to appear. [PMC free article] [PubMed]

12. Einstein DR, Pin FD, Kuprat AP, Jiao X, Carson JP, Kunzelman KS, Cochran RP, Guccione J, Ratcliffe M. Fluid-structure interactions of the mitral valve and left heart: Comprehensive strategies, past, present and future. Communications in Numerical Methods in Engineering. 2009 to appear. [PMC free article] [PubMed]

13. Ganguli A, Persson L, Palmer IR, Evans I, Yang L, Smallwood R, Black R, Qwarnstrom EE. Distinct nf-kappab regulation by shear stress through ras-dependent ikappa-balpha oscillations: real-time analysis of flow-mediated activation in live cells. Circ Res. 2005;96(6):626–34. [PubMed]

14. Hove JR, Koster RW, Forouhar AS, Acevedo-Bolton G, Fraser SE, Gharib M. Intracardiac fluid forces are an essential epigenetic factor for embryonic cardiogenesis. Nature. 2003;421(6919):172–7. [PubMed]

15. Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement. Computer Methods in Applied Mechanics and Engineering. 2005;194:4135–4195.

16. Jansen KE, Shephard MS, Beall MW. On anisotropic mesh generation and quality control in complex flow problems. 10th International Meshing Roundtable. 2001

17. Jiao X. Face offsetting: A unified approach for explicit moving interfaces. J. Comput. Phys. 2007;220:612–625.

18. Jiao X, Colombi A, Ni X, Hart J. Anisotropic mesh adaptation for evolving triangulated surfaces. 15th International Meshing Roundtable. 2006

19. Jiao X, Einstein DR, Dyedov V, Carson JP. Automatic identification and truncation of boundary outlets in complex imaging-derived biomedical geometries. Medical & Biological Engineering & Computing. 2009 Submitted. [PMC free article] [PubMed]

20. Jiao X, Zha H. Consistent computation of first- and second-order differential quantities for surface meshes. ACM Solid and Physical Modeling Symposium. 2008

21. Johnson GA, Cofer GP, Gewalt SL, Hedlund LW. Morphologic phenotyping with MR microscopy: the visible mouse. Radiology. 2002;222(3):789–93. [PubMed]

22. Khamayseh A, Hansen G. Use of the spatial kd-tree in computational physics applications. Commun. Comput. Phys. 2007;2:545–576.

23. Kostelec P, Weaver J, D. H. Multiresolution elastic image registration. Medical Physics. 1998;25(9):1593–1604. [PubMed]

24. Kroger C, Drossinos Y. Particle deposition in a turbulent boundary layer over a large particle size spectrum. J. Aerosol Sci. 1997;28:631–632.

25. Kuprat AP, Einstein DR. An anisotropic scale-invariant unstructured mesh generator suitable for volumetric imaging data. J. Comput. Phys. 2009;228:619–640. [PMC free article] [PubMed]

26. Labelle F, Shewchuk JR. Isosurface stuffing: Fast tetrahedral meshes with good dihedral angles. ACM Transactions on Graphics (ACM SIGGRAPH Conference Proceedings) 2007;26(3):1–57.

27. Longest WP. Comparison of blood particle deposition models for non-parallel flow domains. J. Biomech. 2003;36(3):421–430. [PubMed]

28. Lorensen WE, Cline HE. Marching cubes: A high resolution 3d surface construction algorithm. Computer Graphics. 1987;21:163–169.

29. Mironov V, Visconti RP, Markwald RR. On the role of shear stress in cardiogenesis. Endothelium. 2005;12(5–6):259–61. [PubMed]

30. Remacle JF, Li X, Shephard MS, Flaherty JE. Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods. Int. J. Num. Meth. Eng. 2005;62:899–923.

31. Si H. Adaptive tetrahedral mesh generation by constrained Delaunay refinement. Int. J. Numer. Meth. Engrg. 2008 http://dx.doi.org/10.1002/nme.2318.

32. Smits B. Efficiency issues for ray tracing. SIGGRAPH '05: ACM SIGGRAPH 2005 Courses. 2005:6. http://doi.acm.org/10.1145/1198555.1198745.

33. Taylor CA, Hughes TJR, Zarins CK. Finite element modeling of blood flow in arteries. Comput. Methods Appl. Mech. Engrg. 1998;158:155–196.

34. Treece GM, Prager RW, Gee AH. Regularised marching tetrahedra: improved iso-surface extraction. Computers & Graphics. 1999;23(4):583–598.

35. Yashiro K, Shiratori H, Hamada H. Haemodynamics determined by a genetic programme govern asymmetric development of the aortic arch. Nature. 2007;450(7167):285–8. [PubMed]

36. Zhang Y, Bazilevs Y, Goswami S, Bajaj C, Hughes TJR. Patient-specific vascular nurbs modeling for isogeometric analysis of blood flow. Computer Methods in Applied Mechanics and Engineering. 2007;196(29–30):2943–2959. [PMC free article] [PubMed]

37. Zhang Y, Hughes TJR, Bajaj CL. An automatic 3d mesh generation method for domains with multiple materials. Computer Methods in Applied Mechanics and Engineering. 2009 to appear. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |