New Visualization Techniques
Vol.34 No.1 February 2000
A Three-Dimensional Computer Model of the Human Heart for Studying Cardiac Fluid Dynamics
David M. McQueen and Charles S. Peskin
In all areas of computational fluid dynamics (CFD), proper treatment of the boundary conditions is essential to computing fluid behavior correctly. In many engineering problems, CFD is simplified by a priori knowledge of the motion of the boundary. The well-known parabolic velocity profile in fully-developed flow of an incompressible Newtonian fluid in a pipe of circular cross-section is easily computed because the boundary (the pipe wall) is known to be in a fixed location. Even in more complex settings, such as flow around a ship’s propeller, the motion of the boundary (the propeller) can be specified in advance.
By contrast, in most biological fluid dynamics problems the boundaries are not rigid and their motions are the result of forces imposed on them by the motion of the surrounding fluid. The motion of the fluid, of course, cannot be known without knowledge of the boundary motion. The motion of the boundary and the motion of the fluid form a coupled system; both motions must be computed simultaneously, which makes biological CFD difficult.
A particular problem of interest is the flow of blood in the chambers of the human heart. The heart is an organ whose muscular contractions pump blood around the body. Simplifying somewhat, the heart consists of two main pumping chambers that contract simultaneously. One chamber, the left ventricle, accepts oxygen-enriched blood from the lungs and pumps it to the body. The other chamber, the right ventricle, accepts oxygen-depleted blood from the body and pumps it to the lungs. The inlet and outlet of each ventricle are guarded by valves whose opening and closing guarantee one-directional flow around the circulatory system. There are a total of four valves. The valves generally consist of two or three leaflets - membranes made of very flexible but inextensible material. Familiar examples of materials with this property would be paper or fabric which can be easily bent or twisted but which are not easily stretched. One edge of each valve leaflet is securely attached to the wall of the heart, but the other edge is free of attachment and can move with the flow. Structures analogous to a valve leaflet are a shirt pocket, with one edge (three sides of a rectangular patch pocket) securely stitched to the shirt and one edge free of attachment, or a flag, one edge attached to the flag pole, the other edge free to wave in the wind. When flow is passing through the valve in the forward direction, the valve’s leaflets are positioned out of the way, permitting flow. When flow attempts to pass in the reverse direction, the leaflets come together, their free edges pressing against the free edges of their neighbors to occlude the flow passage.
The motion of the leaflets is not caused by muscles in the valve. The outflow valves are entirely passive structures with no muscular tissue whatsoever. Even in the case of the inflow valves, whose free edges are connected to the heart muscle by a sparse network of tendons, the opening and closing motions result from an interaction with the surrounding fluid. The forward motion of the fluid pushes the leaflets aside out of the main flow stream, but the inextensibilty of the leaflet material prevents free motion of the fluid near the leaflet, affecting the entire flow field. Reverse motion of the fluid causes the leaflets to move back into the flow passage where contact between neighboring leaflets and the inextensibilty of the leaflet material halts the flow. The highly interactive nature of the fluid and leaflet motions makes this an especially interesting and challenging CFD problem.
Commercially available software packages intended for engineering CFD are not equipped to handle this type of dynamic interaction between boundary and fluid. We have developed a numerical method (the "Immersed Boundary Method") which simultaneously computes the motion of a fluid and the motion of an elastic boundary immersed in, and interacting with, that fluid. In the Immersed Boundary Method, the fluid is represented by Eulerian velocities and pressures that are stored on a regular three-dimensional computational lattice. The scale of the heart chambers is such that blood can be treated as a Newtonian fluid. Fluid dynamics is computed by numerical solution of the Navier-Stokes equations, including a body force. The boundary is represented by elastic structures that are free to move continuously in the space sampled by the computational lattice. The essence of the method is to replace the elastic boundary by the forces that result from its deformations. These forces are applied to the lattice in the neighborhood of the elastic boundary with the aid of a numerical approximation to the Dirac delta function. The fluid moves under the action of this body force. The numerical delta function is then used again, to interpolate the newly computed lattice velocities to the locations of the boundary, and then the boundary is moved at the interpolated velocity to a new location (the no-slip condition). The process of computing forces, then fluid motion and then new boundary location is repeated cyclically in a time-stepping procedure with a suitably chosen time step. The only requirements for the method are the physical properties of the fluid, the (possibly time-dependent) elastic properties of the boundary, and the initial geometry of the boundary. A complete description of the Immersed Boundary Method can be found in [1, 2].
Construction of the Model Heart
We are interested in studying blood flow in the heart and have devised a numerical method that uses the (possibly time-dependent) forces arising from an immersed boundary. It seems natural then to construct a computer model of the heart based on the anatomy and mechanical properties of heart muscle fibers. The properties of heart muscle fibers have been extensively documented by physiologists and can be easily incorporated into a computer model.
A remarkable study of the fiber architecture of the heart is . In this work, the hearts of dogs and hogs were bathed in a substance that dissolves the connective tissue between the muscle fibers. Small bundles of the fibers were then pulled away by tweezers in order to determine their paths in three-dimensional space. The paths were reported in a series of black and white drawings and an accompanying verbal description. Figure 1 is indicative of the anatomical drawings in this study. The point is worth emphasizing: two-dimensional drawings like Figure 1 are the data that form the basis for an anatomically correct computer model of the heart. Fortunately  also states verbally the character of the fiber trajectories in different muscle layers, and that, in general, irrespective of the other characteristics of any particular layer, all muscle fibers in the heart originate at some valve ring and terminate at some valve ring. The four valve rings are roughly co-planar, the plane being referred to as the "base of the heart."
We have selected three of the layer types described in  and used them as the basis of a model of the heart. These layers are briefly described below.
We have found that geodesic paths drawn on the surfaces of cones taken in pairs have spatial trajectories similar to those of the fibers in these layers. Other researchers have suggested that muscle fibers follow geodesic paths on surfaces within the heart wall, albeit a different choice of surfaces than cones . The following observation about geodesics on cones is useful: take any cone, oriented with its vertex downward, and then take a cut perpendicular to the axis of the cone. Starting anywhere on the cut, draw any downward-bound geodesic other than a ray of the cone. The geodesic will spiral downward around the cone, approaching the vertex. At some point along its path, however, the geodesic will reach a point of closest approach to the vertex and will then spiral upward away from the vertex, returning to and passing through the plane of the cut. The only geodesics which pass through the vertex are the rays of the cone. (This is a consequence of the intrinsic geometry of cones. A cone may be sliced along a ray and then unrolled onto a plane. All geodesics are straight lines in this plane. Rays are straight lines which pass through the vertex point. All other straight lines have a non-zero minimum distance from the vertex.) If the cut represents a plane near the base of the heart and the vertex represents the apex of the heart, such geodesics have the qualitative character, irrespective of layer, of the fibers of the heart described above. The particular pairs of cones which we use to model each of the three types of layers of the heart described above are shown in Figure 2.
The Outer\Inner Layer is modeled using two cones of differing included angle, one inside the other, vertices coincident. Fibers are modeled by rays of the cones. Rays of the larger cone travel down the exterior, converge at the vertex and then diverge as they travel up the interior. The larger cone represents the epicardium and the smaller cone represents the LV endocardium. See Figure 2a.
The Right-Inner\Left-Outer Layer is modeled using two cones side-by-side with a common ray and coincident vertices. A geodesic which is normal to the common ray is at its point of closest approach to the vertex. Starting at that point, such geodesics will spiral away from the vertex on either cone. In other words, a path may be found which (looking down toward the vertex) follows a geodesic spiraling clockwise downward from the cut on one cone, becomes horizontal at the intersection with the common ray and then follows a geodesic spiraling counter-clockwise upward on the other cone. Each cone models one of the two ventricles. In practice, the cone modeling the right ventricle has a cross-section like a crescent moon, the concave surface coincident with part of the convex surface of the neighboring cone to model the interventricular septum. Even in this changed geometry the geodesic path spirals away from the vertex. See Figure 2b.
Each Left-Ventricular-Cylinder Layer is modeled using two truncated cones, one inside the other with coincident axes (but not coincident vertices) and with a common curve in the plane of truncation. The plane of truncation is chosen normal to the common axis of the cones. A geodesic of either cone which begins tangent to the common curve spirals away from the vertex. Such a geodesic spiraling clockwise downward from the upper cut on one cone can make a seamless transition to the other cone at the common curve and then spiral clockwise upward. There is a V-shaped annular space between any such pair of cones into which another such pair, a scale model of the first pair, may fit. Several such pairs of cones fit together in a nest, like a stack of cookie cutters, discretely modeling the continuous interior of the left ventricular wall. See Figure 2c.
The geodesics on all of these pairs of cone surfaces are computed from their respective minima up to the height of a plane which represents the equatorial plane of the heart. This is the widest part of the heart, slightly below the base of the heart. As stated earlier, all cardiac muscle fibers begin and end on valve rings. In order to terminate the model fibers, the computed geodesics must somehow be extended to the valve rings, which do not lie on the conical surfaces. We have developed a procedure for constructing a surface which smoothly interpolates between an ellipse on a plane at some height Z0 and one or more nonintersecting ellipses on a plane at some other height Z1. The form of the surface is f = 0 where f is some function of the three spatial coordinates. Figure 3 illustrates the procedure for one ellipse at Z0 and two ellipses at Z1; the generalization to more than two ellipses at Z1 is obvious. The restriction to nonintersecting ellipses insures that fk is negative inside (and positive outside) each ellipse. The ellipse is used only because that is the type of conic section which arises in our construction. (The crescent-shaped cross-section of the right ventricle is composed of two elliptical arcs.) If the lower ellipse is taken to be the section of any individual cone in the equatorial plane, and the upper ellipses represent the valve rings appropriate to the layer of which that cone is part, then the interpolating surface models the portion of that layer between the equator and the base of the heart. It is possible to numerically compute geodesic paths on such a surface, given the initial direction. The initial direction is obtained by projecting onto the interpolating surface the direction of the geodesic on the cone as the geodesic passes through the equatorial plane. The geodesic on the interpolating surface is then computed numerically, terminating at any ring in the base. Occasionally, the computed path does not proceed to a valve ring, but turns back to the equator. In that event, the initial direction is perturbed and the path recomputed. A search is performed to find the smallest perturbation which leads to a valve ring.
Figure 4a-4d shows the various layers of the heart model constructed by this procedure. The construction described above refers only to the ventricles. A similar approach is used to construct model atria and portions of the nearby great blood vessels, but the fiber architecture in those parts is somewhat ad hoc. Likewise, the construction of the inflow valves is not based on detailed fiber anatomy, but rather attempts only to capture the gross anatomy of the valves. The outflow valves, on the other hand, are constructed by careful consideration of mechanical equilibrium and numerical solution of a resulting partial differential equation. The details may be found in . Figure 4e shows the four valves of the heart model. Figure 4f shows the entire model. It would be desirable to have the left side of the model heart eject blood into a model systemic circulation which returned blood to the right side of the heart, and the right side eject into a model pulmonary circulation which returned to the left side. The computational demands of that approach would be substantial. Instead, we deal with outflow and inflow to the model heart by closing off the ends of the nearby arteries and veins, and placing sources and sinks near those blind ends. The sources and sinks model the circulatory systems as constant pressure reservoirs connected to the heart by fixed-resistance pathways. In practice, the computed fiber paths are discretized. The entire model is composed of about 4,000 fiber paths requiring about 600,000 discrete points. (The numbers and spacings of fibers, points and layers are parameters which depend on the spatial resolution of the fluid dynamical computation.)
Cardiac Fluid Dynamics and Computer Graphics
Computation of blood flow in the model heart requires specifying the anatomy (that is, tailoring the construction described above to a particular shape), specifying the physiology (the muscular properties of the model fibers, including any time-dependencies) and then the repeated application of the Immersed Boundary Method. Currently, this involves a finite difference solution on a 128x128x128 grid repeated approximately 57,000 times to model a single beat of the heart. Such a computation takes about 250 CPU hours on the Cray T90, a large shared-memory multiprocessor, at the San Diego Supercomputer Center. The result is a prediction of the flow patterns of blood throughout the heart and a simultaneous prediction of the motion of the heart and its valves. Some selected results from one such computation are shown in Figure 5.
It seems almost superfluous to remark that computer graphics is essential to every facet of this work from model design to analysis of results. It would be unimaginable to attempt construction of so complex a three-dimensional object without computer graphics, particularly interactive computer graphics. How else would one assess whether the designed model resembles the object being modeled, particularly when the target design is only available pictorially? It is equally unimaginable to attempt to make sense of the CFD results without interactive computer graphics. The data sets from three-dimensional models are often gigantic, and ours are no exception. The interesting phenomena are often localized, but not necessarily in locations the researcher knows in advance. Without powerful interactive graphics which permit the researcher to see both the large picture and the details, from any desired point of view, simply locating the interesting phenomena, let alone quantifying them, becomes a major challenge.
Results are visualized by a custom display program which uses the GL libraries on SGI workstations. This program allows us to animate the results of a computation, turn on or off any fibers of the heart model, interactively change our point of view and magnify any region of interest. The flow in the heart is represented by the computed motions of fluid markers. Optionally, streaklines showing the recent positions of the markers can be turned on, but this only permits visualization of the flow field where there happen to be markers. One of the great problems we face is in characterizing the overall velocity field.
The computed velocity is a vector field on a grid with more than 2,000,000 storage locations, each storing three velocity components. When all these vectors are displayed, the density of lines obscures the image of the model heart. Even when only a small region of the vector field is displayed, it is often difficult to assess the component of a vector in the direction normal to the screen of the display. The velocity field could be visualized by interactively introducing fluid markers in any arbitrary location and transporting them in the vector field. This is a computationally demanding task which needs the full vector field for any time of interest. To visualize by this approach requires a very fast processor, substantial RAM and a significant investment in disk storage. This is the kind of computation which is more properly done on a supercomputer, as it makes demands on the workstation which are likely to reduce interactive graphics performance. Visualization of the velocity field is an area we are still actively researching.
Although the muscle fibers of the model heart are critical to its function, it can be difficult to make sense out of the wire-frame presentation, which is most natural. Rendering the layers of the model as surfaces is sometimes useful, but there are some technical difficulties resulting from the complicated paths of the geodesics, especially on the interpolating surfaces. Our technique is: on each layer, find the intersections of the model fibers with a set of closely spaced planes parallel to the base of the heart. The difficulty is that these points of intersection are not ordered in any particular way. To find the natural ordering on each plane, we solve the Traveling Salesman Problem. (For the heart model, the largest such problem involves about 200 points. Several hundred such problems need to be solved to surface-render the entire heart. We use the simulated-annealing approach of .) Between the contours on each pair of neighboring planes, the surface to be rendered is a closed ribbon. We choose the triangulation of this ribbon which minimizes the aggregate edge-length. The edges between contours generally follow the fibers, and we apply the triangulation of the initial data to all subsequent data sets. Figure 6 shows a surface rendering of the outer layer of the heart model.
We are currently adjusting the mechanical properties of the model heart to achieve physiologically normal pressures and flows. Our medium-range goals are to investigate certain disease states and to design some prosthetic devices, particularly prosthetic valves. Our long-range goal is to be able to predict flow patterns in the heart in something approaching real time, and to be able to customize the anatomy and physiology of the model to match arbitrary data. It would then be possible to use this approach as an aid to medical diagnosis and treatment. This goal requires significant advances in computer power, algorithm design and interactive graphics.
David M. McQueen is Research Scientist at the Courant Institute. For more than 20 years, he has been applying the immersed boundary method to the computational study of blood flow in the heart.
Charles S. Peskin is Professor of Mathematics at the Courant Institute. He is the creator of the immersed boundary method and has won numerous awards for his innovative work in mathematical biology, including the MacArthur Foundation Prize Fellowship and the Sidney Fernbach Award of the IEEE Computer Society.
The copyright of articles and images printed remains with the author unless otherwise indicated.