**Charles S. Peskin and David M. McQueen
Courant Institute of Mathematical Sciences
New York University
251 Mercer Street
New York, NY 10012 USA
peskin@cims.nyu.edu
mcqueen@cims.nyu.edu
Biological Fluid Dynamics
ed. by C. P. Ellington and T. J. Pedley
The Company of Biologists Limited
Cambridge, 1995**

** **

At this Symposium on Biological Fluid Dynamics, it is appropriate to ask whether there is any common theme that unites the diverse problems that arise in the study of living systems interacting with fluids. The answer that immediately comes to mind is this: biological fluid dynamics invariably involves the interaction of elastic flexible tissue with viscous incompressible fluid. (In many cases the tissue is not only elastic, it is also active, i.e., capable of doing work on the fluid.) This paper describes the immersed boundary method, which is a general framework for the computer simulation of biofluid dynamic systems. This method has already been applied to blood flow in the heart (including the computer-assisted design of prosthetic cardiac valves), platelet aggregation during blood clotting, aquatic animal locomotion, wave propagation along the basilar membrane of the inner ear, and flow in collapsible tubes. In the immersed boundary method, the elastic (and possibly active) biological tissue is treated as a part of the fluid in which additional forces (derived from the tissue stresses) are applied. Because the tissue is represented in terms of its force field, the method remains straightforward even when the geometry of the biological tissue is complicated, dynamic, and not known in advance.

A GENERAL METHOD FOR THE COMPUTER SIMULATION OF BIOLOGICAL SYSTEMS INTERACTING WITH FLUIDS

Introduction

This paper is intended as a tutorial on the immersed boundary method, which has been applied to a wide range of problems, mostly in biofluid dynamics, including blood flow in the heart (Peskin, 1972, 1977, 1980, 1992; Peskin and McQueen, 1980, 1989, 1992a, 1992b; McQueen and Peskin, 1983, 1985, 1989, 1990, 1991; McQueen, Peskin, and Yellin, 1982; McCracken and Peskin, 1980; Meisner, McQueen, et al., 1985; Peskin and Printz, 1992), platelet aggregation during blood clotting (Fogelson, 1984, 1985; Fauci and Fogelson, 1993), flow of suspensions (Fogelson and Peskin, 1988; Sulsky and Brackbill, 1991), fluid dynamics of the inner ear (Beyer, 1992), aquatic animal locomotion (Fauci and Peskin, 1988; Fauci, 1990; Fauci and Fogelson, 1993), and flow in collapsible tubes (Rosar, 1994). Generally speaking, the method is applicable to any problem in which a fluid interacts with an elastic material. The elastic material may be an active one, as in the case of cardiac muscle or a swimming fish.

The philosophy of the immersed boundary method is to treat the elastic material as a part of the fluid in which additional forces (arising from the elastic stresses) are applied. The fluid equations are solved on a regular cubic lattice, the structure of which is not modified in any way by the presence of the immersed elastic bodies, the geometry of which may be quite complicated. The elastic material is tracked in Lagrangian fashion, by following a collection of representative material points. The spatial configuration of these points is used to compute elastic forces, which are applied to the nearby lattice points of the fluid. The fluid velocity is updated under the influence of these forces, and the new velocity is then interpolated at the elastic material points, which are moved at the interpolated velocity to complete the time step.

The various applications of the immersed boundary method have been listed above (with references), but the following remarks may give the reader a more detailed picture of what has already been done. The method has been used to design a prosthetic cardiac valve, for which a patent has been issued. It has been used to create a detailed model of the platelet aggregation process during blood clotting, a model that includes the platelets themselves, the blood plasma in which they are immersed, elastic interplatelet links as well as links between a platelet and the injured vessel wall, secretion of adenosine diphosphate (ADP) that occurs upon platelet activation, the convection and diffusion of ADP in the flowing blood, and the activation of platelets that experience an above-threshold level of ADP as well as platelets that encounter the injured region of the vessel wall. Because the immersed boundary method is used, the fluid dynamics automatically adapt to the changing geometry of the flow region as the platelet aggregate grows. The immersed boundary method has also been used to study the swimming of various creatures, from undulating eels to algae that do the breast stroke. Interactions of swimming spermatoza with each other and with elastic channel walls have also been studied. Because the immersed boundary method is used, these studies automatically include the mechanical feedback influence of the hydrodynamics on the swimming motions of the creature in question. The deformations of the creature are not rigidly prescribed but are influenced by hydrodynamic forces. In sedimentation problems, the immersed boundary method makes it possible to study the effects of particle elasticity, and also to study large numbers of interacting particles, since the computational effort grows only linearly with the number of immersed particles. Wave propagation along the basilar membrane of the inner ear and flow in collapsible tubes are two problems for which a variety of idealized models have been introduced and studied. In such problems, the immersed boundary method opens up the prospect of raw simulation, almost akin to experiment, in which something approaching the full complexity of the system is allowed to operate in the model.

The next section of this paper describes a mathematical formulation of the equations of motion of a viscous incompressible fluid containing an immersed system of elastic fibers. This formulation forms the foundation of the immersed boundary method, which is then described in the subsequent section. At the core of the immersed boundary method is a smoothed approximation to the Dirac delta-function, the construction of which is given next. The paper concludes with an example application which is qualitatively discussed (with computer-generated figures): blood flow in a three-dimensional computer model of the mammalian heart.

Mathematical Formulation

Consider a viscous incompressible fluid containing an immersed system (a continuum) of elastic fibers. The fibers are pure force-generators. By themselves, they have neither mass nor volume, but together with the fluid in which they are immersed, they form an incompressible visco- elastic material. This material is highly anisotropic, since the fiber stress points always in the fiber direction.

Let the unknown motion of the fibers be described in Lagrangian form:

in which (*q*,*r*,*s*) defines a material point and (*q*,*r*) defines a fiber.
The unit tangent to the fibers is given by

and the fiber tension *T* is given by a generalized Hooke's law:

Here is the force transmitted by the bundle of fibers
, and determines the
local fiber strain, from which this force is computed. The explicit
dependence of on *q*,*r*,*s*,*t* allows for different stress-strain
relations at different locations and times. In particular, the explicit
time dependence makes possible the representation of active materials,
such as muscle.

Let *B* be some arbitrary region of the parameter
plane. This defines a bundle of fibers. Consider the segments of these
fibers defined by . Since the fibers are massless,
the total force acting on these segments must be zero. This total force
includes both the force of the fluid on the fibers and also the force
transmitted across the surfaces and by the fibers
themselves. Thus

Using Newton's law of equal and opposite forces, and also the fundamental theorem of calculus, we may rewrite this equation as follows:

Since *B*, , and are arbitrary, this shows that the density
of the fiber force with respect to the measure is given by

Expanding the derivative, we see that the fiber force density has a tangential component and a component in the direction of the principal normal to the fibers:

where is the fiber curvature and is the principal normal. There is no component of the fiber force in the binormal direction .

We may now formulate equations of motion for the fiber-fluid system.
Let the fluid have density and viscosity , and let , be the velocity and
pressure of the fluid. Note that these are expressed in Eulerian form;
they are functions of the fixed position and the time *t*.
The equations of motion are as follows:

Equations 8-9 are the familiar Navier-Stokes equations
of a viscous incompressible fluid (in Eulerian form). The applied force
density **F** is here used to express the effect of the fibers acting
in the fluid, see below. Equations 12-14 determine
the fiber forces. (These equations have been derived and discussed
above and are repeated here in order to display the entire system of
fiber-fluid equations in one place.) Note that these fiber equations
are in Lagrangian form: , , and are
functions of the Lagrangian parameters *q*, *r*, *s* (and of the time
*t*).

Equations 10-11 are the interaction equations which couple the fibers and the fluid. In these equations denotes the three-dimensional Dirac delta function. Equation 10 expresses the fiber force density as a sum (integral) of point forces applied at the locations of the fiber points. To see more clearly the meaning of this equation, we integrate both sides over an arbitrary region which may or may not contain some of the fibers:

which simply says that the total force of the fibers on the region
at time *t* can be found by integrating over all of the fiber
segments (if any) that happen to lie within the region at time
*t*.

Equation 11 is the no-slip condition of a viscous fluid, which here plays the role of an equation of motion for the fibers. It states that the fibers move at the local fluid velocity. The reason for writing this out in terms of the Dirac delta function (the second form of Equation 11) will become clear when we consider the numerical scheme that is used to solve the fiber-fluid equations.

It is important to note that the formulation given above still makes
sense when the fibers are confined to a surface (a heart valve leaflet,
for example) instead of filling up a region of space. In that case, we
merely drop one of the Lagrangian parameters (say *q*), and
Equations 10-11 become

The delta function appearing in these equations is still three-dimensional, but Equation 16 now involves only a double integral. It therefore defines as a delta-function layer with support along the immersed elastic boundary only. Although is now infinite on this internal boundary, its integral over any (finite) volume is finite. (Equation 17, by contrast, involves a triple integral which gives directly a finite result without any further integration.)

The Immersed Boundary Method

We now consider the numerical solution of Equations 8-14. Let time proceed in steps of duration , and let superscripts denote the time step index. Thus . At the beginning of the time step and are known. We shall describe the computation of and .

The fluid equations are discretized on a cubic lattice of mesh width *h*.
The following difference operators are applied to functions defined on
this lattice:

where and where , , and are unit vectors in the coordinate directions.

Similarly, the fiber equations are discretized on a rectangular lattice in the space of the Lagrangian parameters (but note that this is some curvilinear lattice in physical space, and that the physical locations of the fiber lattice points are not, typically, those of the fluid lattice.) Let , , be the three fiber-lattice spacings. For any function , let the central difference operator be defined by

Each time step proceeds as follows. First, compute the fiber force from the fiber configuration . This is done using the following equations:

Suppose is defined at integer points . Then and are defined at half-integer points . Finally is defined once again at the integer locations.

The next step is to evaluate the fiber force density as follows:

Here denotes the sum over the values of *q*,*r*,*s* given by
,
where *i*, *j*, and *k* are integers. The function is a
smoothed approximation to the three-dimensional Dirac delta function.
It will be constructed in the next section.

With known, we are ready to update the fluid velocity by integrating the Navier-Stokes equations. This is done by solving the following linear system for , :

In Equation 27, the notation allows for a choice between the forward and backward difference operators. This choice is to be made by upwind differencing. Thus

Because of the upwind differencing, the scheme given by Equations 27-28 is stable, provided that

at all of the points of the fluid lattice.

Note that Equations 27-28 are linear constant-coefficient difference equations in the unknowns . They can therefore be solved by Fourier methods, and this can be made very efficient through the use of the Fast Fourier Transform (FFT) algorithm. It is a consequence of our use of a regular cubic lattice in the discretization of the fluid equations that we are able to exploit the FFT in this way. If the fluid grid had been made to conform to the geometry of the immersed elastic bodies, Fourier methods would not be applicable.

Once has been determined, the next (and last) step is to move the fibers. This is done through a discretization of Equation 11:

where denotes the sum over the cubic lattice of the fluid computation. This is the set of points of the form , where , , and are integers. As in Equation 26, is a smoothed approximation to the Dirac delta function. The construction of is the subject of the next section. Since and have been updated, the time step is complete.

Construction of

Let

where , , and are the components of , and where has the following properties:

**i)**- is continuuous
**ii)**- for
**iii)**- , all
*r* **iv)**- , all
*r* **v)**- , all
*r*

The motivation for imposing these requirements is as follows. Continuity of is needed to avoid sudden jumps in the fiber velocity or force as the fiber points move through the computational lattice of the fluid. Bounded support of is needed to keep down the computational cost of the algorithm: it would be very expensive to have each fiber point interact directly with all of the lattice points of the fluid. (The specific bound 2 which appears in property (ii) defines the smallest support that is consistent with the other criteria.)

Property (iii) has the immediate consequence that
for all *r*. When is used for
interpolation, this ensures that the interpolation of a constant is
exactly correct. When is used for applying the fiber force
to the fluid, this same property assures that the correct total force
is applied by each fiber point (conservation of momentum).

The stronger conditions of property (iii) concerning the separate sums over even and odd points are helpful because they prevent the appearance of a point-to-point oscillation that would otherwise be generated by the application of localized forces. This is related to the use of the central difference operator in the discretization of the Navier-Stokes equations.

When combined with property (iii), property (iv) guarantees that the interpolation of linear functions will be exact and that the correct total torque (with respect to any chosen origin) will be applied to the fluid by each fiber point (conservation of angular momentum).

Property (v) is the most interesting one. It is applied because it guarantees the following inequality:

which holds for all and which follows from (v) by an application of the Schwarz inequality.

Expressions like those on the left-hand side of this inequality arise if the function is used first to apply the fiber forces to the fluid lattice and then to interpolate the resulting force field back to the fiber points themselves. In this two-step process, one would like the interaction of a fiber point with itself to be constant, independent of where that point sits with respect to the fluid lattice, and this is precisely what property (v) asserts. Moreover, one would like the interaction of one fiber point with another be be bounded by the interaction of a fiber point with itself, which is precisely the content of the above inequality.

We leave it as an exercise for the reader to show that the five properties listed above uniquely determine the function , as follows:

(Note that the second line defines on
in terms of the values of on . One could,
of course, write out a formula for this, but it is more efficient to
compute as indicated above, since the typical situation is that
is simultaneously needed at four points of the form
*r*-2, *r*-1, *r*, and *r*+1, where .)

Even though we have only postulated the continuity of itself and said nothing about its derivatives, it is a remarkable fact that the continuity of the first derivative of turns out to be implicit in the five properties proposed above. This can be checked by differentiating the above formulae for to obtain , , and .

For readers who are familliar with our previous work, it should be noted that the function used here is different from the trigonometric function used previously:

One can check that satisfies four of the five properties that were used to define , the one exception being the first-moment or torque condition, property (iv). Although the expressions for and may look very different, their numerical values are startlingly close: plot their graphs with a thick pen and you will have trouble distinguishing them. Thus, the change from to is not likely to have any practical effect on the computed results. We have gone ahead and made the change nonetheless, since it turns out that is slightly cheaper to compute than , and since it is satisfying to have uniquely determined by a reasonable list of axioms.

Example: A Computer Model of the Heart and its Valves

We have used the methodology described above to develop a three-dimensional computer model of the mammalian heart, including the heart valves and the nearby great vessels. The model is constructed as a collection of fibers immersed in a (periodic) box of fluid. These fibers are arranged in such a manner as to simulate the true fiber architecture of the heart, as described by Thomas (1957) and by Streeter et al. (1969,1978). The general mode of construction that we use is to define surfaces within the heart wall on which the fibers are to run, and then to lay out the fibers as geodesics (curves of shortest length) on these surfaces. In the special case of the arterial valves, we have introduced a mathematical theory (Peskin and McQueen, 1994) that can be used to determine the fiber architecture of the valve leaflets from first principles, and we have used that theory to generate the arterial valves of the model heart. It is an open problem how to do the same for the atrioventricular valves, and indeed for the heart as a whole. (There is such a theory available for the left ventricle (Peskin, 1989), but it relies on axial symmetry and we have not used it here.)

The model that we have constructed is anatomically complete. It includes left and right atria; left and right ventricles; aortic, pulmonic, mitral, and tricuspid valves; ascending aorta and main pulmonary artery; superior and inferior vena cavae; and four pulmonary veins. The model arteries and veins are equipped with sources and sinks which apply realistic pressure loads to the model heart. An external source/sink allows for the changes in volume of the heart that occur during the cardiac cycle. The fibers that comprise the model vessels and valve leaflets have fixed elastic parameters, but those of the atrial and ventricular walls have time-dependent stiffnesses and rest lengths so that they can contract and relax at appropriate times in the cardiac cycle. All fibers have nonlinear stress-strain relations with increasing stiffness at greater loads.

All physical parameters of the model with the exception of the blood viscosity have been set appropriately for the human heart. The viscosity has been increased by a factor of 25, which means that the Reynolds number has been reduced by the same factor. Smaller viscosity requires a finer computational mesh, as least near the boundary. The present form of the immersed boundary method uses a uniform mesh, which is therefore expensive to refine, but future research may make it possible to perform local mesh refinement or indeed to use a grid-free method that automatically concentrates the computational effort in regions of high vorticity. Even with uniform meshes, future progress in computer hardware will gradually make it feasible to use finer and finer meshes and hence to reduce the viscosity of the simulated blood.

Figures 1-3 show the heart model in action. In all three of these figures, part (a) shows the model during ventricular diastole and part (b) during ventricular systole. The figures all show the heart in the same orientation, seen from the front, with the left atrial appendage in the foreground. Figure 1 shows the external view of the model heart. In Figure 2, the front of the heart has been cut away to reveal the blood flow inside. The sectioning plane passes through the mitral, aortic, and pulmonic valves, and also cuts close to the apex of the left ventricle. Blood flow is depicted in terms of fluid markers which leave trails behind them so that their motion during the immediate past can be appreciated. Figure 3 shows a somewhat transparent view of the entire heart model, with the blood flow seen through the fibers of the heart walls. Only a small subset of the model fibers are shown; otherwise the walls would appear solid and the blood flow could not be seen.

The model heart beats, but it is not yet a healthy heart. The most obvious pathology is mitral regurgitation, and aortic regurgitation is also present. Further work on the valves is needed, and other changes in the anatomy and physiology of the model heart may be required as well. It is not an easy task to imitate nature's design, but we are doing our best. Once the model has been perfected, it will be used as a test chamber for the computer-assisted design of prosthetic cardiac valves, in the simulation of disease processes affecting the mechanical function of the heart or its valves, and as a tool for achieving improved understanding of the normal physiology of the heart.

Acknowledgments

The work described in this paper is supported by the National Science Foundation (USA) under research grant BIR-9302545. Computation is performed at the Pittsburgh Supercomputing Center under a grant (MCA93S004P) of Cray C-90 computer time from the MetaCenter Allocation Committee (USA).

REFERENCES

Beyer, R.P. (1992). A computational model of the cochlea using the immersed boundary method. J. Comput. Phys. 98:145-162.

Fauci, L.J. (1990). Interaction of oscillating filaments - a computational study. J. Comput. Phys. 86:294-313.

Fauci, L.J. and Fogelson, A.L. (1993). Truncated Newton methods and the modeling of complex immersed elastic structures. Comm. Pure Appld. Math. 46:787-818.

Fauci, L.J. and Peskin, C.S. (1988). A computational model of aquatic animal locomotion. J. Comput. Phys. 77:85-108.

Fogelson, A.L. (1984). A mathematical model and numerical method for studying platelet adhesion and aggregation during blood clotting. J. Comput. Phys. 56:111-134.

Fogelson, A.L. (1985). Mathematical and computational aspects of blood clotting. In Proceedings of the 11th IMACS World Congress on System Simulation and Scientific Computation, Vol. 3 (ed. B. Wahlstrom), pp. 5-8, North Holland, New York.

Fogelson, A.L. and Peskin, C.S. (1988). A fast numerical method for solving the three-dimensional Stokes' equations in the presence of suspended particles. J. Comput. Phys. 79:50-69.

McCracken, M.F. and Peskin, C.S. (1980). A vortex method for blood flow through heart valves. J. Comput. Phys. 35:183-205.

McQueen, D.M. and Peskin, C.S. (1983). Computer-assisted design of pivoting-disc prosthetic mitral valves. J. Thorac. Cardiovasc. Surg. 86:126-135.

McQueen, D.M. and Peskin, C.S. (1985). Computer-assisted design of butterfly bileaflet valves for the mitral position. Scand. J. Thor. Cardiovasc. Surg. 19:139-148.

McQueen, D.M. and Peskin, C.S. (1989). A three-dimensional computational method for blood flow in the heart: (II) contractile fibers. J. Comput. Phys. 82:289-297.

McQueen, D.M. and Peskin, C.S. (1990). A heart valve prosthesis. European Patent Publication Number EP 0 211 576 B1.

McQueen, D.M. and Peskin, C.S. (1991). Curved butterfly bileaflet prosthetic cardiac valve. U.S. Patent Number 5026391.

McQueen, D.M., Peskin, C.S., and Yellin, E.L. (1982). Fluid dynamics of the mitral valve: physiological aspects of a mathematical model. Am. J. Physiol. 242:H1095-H1110.

Meisner, J.S., McQueen, D.M., Ishida, Y., Vetter, H.O., Bortolotti, U., Strom, J.A., Frater, R.W.M., Peskin, C.S., and Yellin, E.L. (1985). Effects of timing of atrial systole on LV filling and mitral valve closure: computer and dog studies. Am. J. Physiol. 249:H604-H619.

Peskin, C.S. (1972). Flow Patterns Around Heart Valves: A Digital Computer Method for Solving the Equations of Motion. Ph.D. thesis, Physiology, Albert Einstein College of Medicine. University Microfilms #72-30, 378. 211pp.

Peskin, C.S. (1977). Numerical analysis of blood flow in the heart. J. Comput. Phys. 25:220-252.

Peskin, C.S. (1980). Fluid dynamics of the heart and the ear. In Computing Methods in Applied Sciences and Engineering (ed. R. Glowinski and J.L. Lions), pp. 587-613, North Holland, Amsterdam.

Peskin, C.S. (1989). Fiber architecture of the left ventricular wall: an asymptotic analysis. Comm. Pure Appld. Math. 42:79-113.

Peskin, C.S. (1992). Two examples of mathematics and computing in the biological sciences: blood flow in the heart and molecular dynamics. In American Mathematical Society Centennial Publications vol 2, pp. 395-416. American Mathematical Society, Providence, Rhode Island.

Peskin, C.S. and McQueen, D.M. (1980). Modeling prosthetic heart valves for numerical analysis of blood flow in the heart. J. Comput. Phys. 37:113-132.

Peskin, C.S. and McQueen, D.M. (1989). A three-dimensional computational method for blood flow in the heart: (I) immersed elastic fibers in a viscous incompressible fluid. J. Comput. Phys. 81:372-405.

Peskin, C.S. and McQueen, D.M. (1992). Computational biofluid dynamics. Contemp. Math. 141:161-186.

Peskin, C.S. and McQueen, D.M. (1992). Cardiac Fluid Dynamics. Crit. Rev. Biomed. Eng. 20:451-459.

Peskin, C.S. and McQueen, D.M. (1994). Mechanical equilibrium determines the fractal fiber architecture of aortic heart valve leaflets. Am. J. Physiol. 266:H319-H328.

Peskin, C.S. and Printz, B.F. (1992). Improved volume conservation in the computation of flows with immersed elastic boundaries. J. Comput. Phys. 105:33-46.

Rosar, M.E. (1994). A Three-Dimensional Computer Model for Fluid Flow Through a Collapsible Tube. Ph.D. Thesis, New York University.

Streeter, D.D. Jr., Powers, W.E., Ross, M.A., and Torrent-Guasp, F. (1978). Three-dimensional fiber orientation in the mammalian left ventricular wall. In Cardiovascular System Dynamics (ed. J. Baan, A. Noordergraaf and J. Raines), pp. 73-84, MIT Press, Cambridge, MA.

Streeter, D.D. Jr., Spotnitz, H.M., Patel, D.P., Ross, J. Jr., and Sonnenblick, E.H. (1969). Fiber orientation in the canine left ventricle during diastole and systole. Circ. Res. 24:339-347.

Sulsky, D. and Brackbill, J.U. (1991). A numerical method for suspension flow. J. Comput. Phys. 96:339-368.

Thomas, C.E. (1957). The muscular architecture of the ventricles of hog and dog hearts. Am. J. Anatomy 101:17-57.

Legends

**Figure 1:**- External view of the model heart (a) during ventricular diastole and (b) during ventricular systole. The heart is shown from the front, so the left side of the heart is on the right side of the figure, with the left atrial appendage (auricle) in the foreground. Structures seen at the top of the figure are (from left to right in the figure) main pulmonary artery, ascending aorta, and the left atrium with three of the four pulmonary veins visible. The right atruim and the vena cavae are in back and cannot be seen. Thin black lines indicate some of the fiber trajectories on the surface of the model heart.
**Figure 2:**- Cut-away view of the model heart in the same orientation and at the same times in the cardiac cycle as in Figure 1. Blood flow is depicted in terms of streak-lines. Mitral regurgitation is clearly visible during ventricular systole (b).
**Figure 3:**- Transparent view of the model heart showing selected
fibers and the blood flow within the chambers. Frames shown correspond
to those of Figures 1-2.

Tue Dec 16 11:42:05 EST 1997