SIAM News Blog

Beating in Fluid: Hearts and Cilia by the Immersed Boundary Method

By Debbie Sniderman

Heart valves are thin flexible membranes immersed in fluid that profoundly affect blood flow. Each side of the heart uses valves to aid in the movement of blood from the atrium to the ventricle, and from the ventricle to the aorta or the pulmonary artery. Each ventricle has an inflow and an outflow  valve, which allow the alternate contraction and relaxation of that ventricle to generate unidirectional flow. The valves open and close in coordination with the pumping of the heart in its two main states, the diastolic or relaxed phase and the systolic or contracted phase (see Figure 1). Heart valves let flow through when open, but the valve leaflets also shear the flow and create vortices that help the valves close efficiently. When closed, the valves prevent backflow.

Figure 1. Diagram showing the physiology of the left side of the heart. The right heart is qualitatively the same but quantitatively different. Upper left. Heart in its two main states, diastole (relaxed) and systole (contracted). LA=left atrium, LV=left ventricle, Ao=aorta. Bottom left. Pressure/volume diagram of the cardiac cycle. Each corner of the rectangular path denotes the opening or closing of a valve. Inflow and outflow pressures are idealized as constant. ∆V denotes the stroke volume. Upper right. Pressures in the Ao, LV, and LA. The ventricle switches between atrial pressure when relaxed and aortic pressure when contracted. Bottom right. Flows through the Ao and mitral (Mi) valves, with the heart sounds indicated. Image credit: Charles Peskin.

Heart valve function presents a fluid-structure interaction problem, and the immersed boundary (IB) method offers a solution. It unifies the fields of elasticity and fluid dynamics, enabling a wide variety of applications in biological and engineering mechanics. Charles Peskin (New York University) first introduced the IB method in his Ph.D. thesis at the Albert Einstein College of Medicine, and further developed it with his students and colleagues at New York University’s Courant Institute of Mathematical Sciences.

IB Method

The IB method [8] treats immersed elastic bodies or boundaries as a part of the fluid in which they are immersed. Structures are represented as elastic, and material points of the structure are tracked so that their spatial configurations can be used to compute elastic forces applied to the fluid. Fluid velocity and pressure are computed on fixed Cartesian grids. 

Image of an aortic valve leaflet stained for collagen, one of the three leaflets that regulates outflow from the left ventricle. It displays a branching and braided fiber architecture. Image Credit: A. Sauren.
The moving structure grid and the fixed fluid grid do not conform to each other, so the IB method’s key problems are the following: determining how to apply the forces generated by the structure to the fluid and deciding how to evaluate the fluid velocity at the structure points. Use of a smoothed approximation to the Dirac delta function, constructed so that important physical quantities such as momentum and angular momentum are preserved during fluid-structure interaction, solves both of these problems. A recently-developed delta function with a constant second moment and three continuous derivatives enables a new version of the IB method, in which the interpolated velocity field is exactly divergence-free.

The principal advantage of the IB method over other methods for fluid-structure interaction is that there is no need to adapt the fluid mesh to that of the structure, or vice versa. There is also no need to deal separately with the many ways that the structure mesh may intersect with the fluid mesh.

The IB method is modular; it separates the task of modeling the fluid from that of modeling the immersed structure. Any structural model that can generate elastic forces from the spatial configuration of the structure mesh’s nodes can be used. Peskin used networks of points connected by springs for structural models in his early work. More recently, finite element structural models are being immersed in fluid by the IB method. Since there is fluid everywhere in an IB computation, even inside an immersed structure, the structural model does not need to enforce incompressibility, a task that is handled by the fluid solver [1].

Immersed boundary software known as IBAMR, which was created by Boyce Griffith (University of North Carolina at Chapel Hill) and is publically available, incorporates both distributed parallelization and adaptive mesh refinement of Cartesian grids [2].

Beating Hearts

Figure 3. Snapshot of a visualization of regions of high velocity flow on the front view of an IB heart model stimulated by electrical activity. The red color indicates regions in which the magnitude of the velocity is above an arbitrarily chosen threshold. A is in diastole with blood entering the ventricles through the mitral and tricuspid valves. B is in systole, with blood being ejected from the ventricles through the aortic and pulmonic valves. In both frames the thin-walled right ventricle is on the left side of the figure. Image credit: Boyce Griffith, David McQueen, and Charles Peskin.
Peskin’s work on the heart began with the goal of improving the design of artificial valves; for recent work of this kind, see [3]. The next step, carried out jointly with David McQueen (New York University), was to build a model heart around the valves [5, 6].

Heart muscle is made of fibers that take geodesic paths on nested toroidal surfaces within the heart walls. Peskin introduced equations for this fiber architecture [7] and found asymptotic solutions in the special case of axial symmetry (applicable to the left ventricle). But solving these equations in the non-axial symmetric case of the whole heart remains an open problem. 

Numerical solutions by McQueen and Peskin [9] of similar equations for the collagen fiber architecture of the aortic valve reveal a fractal structure with dimension 2.2 [10]. Alexander Kaiser (New York University) recently created a fiber-based model of the mitral valve with intricate chordae tendineae. In the IB model heart, all mathematical fibers are immersed or embedded in a viscous and incompressible fluid, creating a fiber-reinforced fluid that has mass, volume, and incompressibility. The stress/strain relationship of the muscle fibers is both nonlinear and time-dependent, with much higher stiffness during systole, or contraction of the heart, than diastole, the relaxation phase. This time dependence of material properties drives the model heart through the cardiac cycle (see Figures 2 and 3).

Figure 2. Mitral valve model. Three perpendicular views of a mitral valve model simulated in a test chamber by the IB method. Colors indicate the vertical component of fluid velocity. Credit: Alexander Kaiser.

Model heart experiencing the cardiac cycle. Credit: Boyce Griffith, David McQueen, and Charles Peskin.

Beyond Mechanics and Fluids: Electrical Integration

The IB method’s flexibility also allows for the integration of the electrical component of the heart [4] into the mechanical component/fluid modeling. The heart has an autonomous system that generates electrical waves to control beating. Mechanical-to-electrical feedback operates through stretch-activated channels and passive changes in electrical resistance as cardiac tissue deforms during contraction and relaxation.

In joint work with Griffith, Peskin uses the formalism of the IB method to simulate the electrical and mechanical activity of the heart within the same framework and in the same software. Their work is based on the bidomain model, which tracks both intracellular and extracellular voltage and current. In the IB formulation of the bidomain equations, the extracellular domain is analogous to the fluid and the intracellular domain is analogous to the structure.

In the electrical equations of the extracellular and intracellular spaces, the currents leaving the cells appear as a source/sink term. This can be written mathematically as a distributed delta function, which spreads the current onto the extracellular domain identically to how forces are applied to fluid in the IB method. The transmembrane current depends on both intracellular and extracellular voltage, and the evaluation of the extracellular voltage at a given intracellular point is analogous to the evaluation of fluid velocity at a structure point in the IB method. Thus, the entire electrical problem looks like an IB problem.

In the IB formulation, the extracellular space extends beyond the myocardium and into the blood and surrounding tissues, both of which are electrically conducting media. The model heart thus acquires an electrocardiogram.

Beating Cilia

Another example of beating in fluid on a microscopic scale occurs in motile cilia, which are the active hair-like appendages of biological cells. The beat of a single cilium is driven by hundreds of dynein molecular motors that coordinate their activity to produce wavelike bending of the cilium. 

In a model proposed recently by Peskin and Jihun Han (New York University), the coordination emerges spontaneously as a result of a simple dynamical law governing the tension of each dynein motor, together with a geometrical constraint based on the microstructure of the cilium. The IB method immerses the model cilia in fluid, facilitating the study of their fluid-mediated interaction (see Figure 4). These studies reveal a striking tendency towards synchronization of nearby cilia regardless of initial conditions. 

Figure 4. Snapshot of two rows of modeled cilia in a periodic box of fluid. Red particles are fluid markers leaving trails to show recent trajectories. Despite in-phase initial conditions, the two rows temporarily go into anti-phase and later recover synchrony before computations end. Long-term behavior is unknown. The more vertical cilia are moving left (power stroke) and the others are moving right (recovery stroke). Credit: Jihun Han.

The mathematical problems surrounding hearts and cilia are so different in scale, yet both involve spontaneous oscillation and fluid-structure interaction. The immersed boundary method is applicable to both.

This article is based on an invited lecture by Charles Peskin at the SIAM Conference on the Life Sciences, which was held in Boston this July.

[1] Devendran, D., & Peskin, C.S. (2012). An immersed boundary energy-based method for incompressible viscoelasticity. Journal of Computational Physics, 231(14), 4613-4642.

[2] Griffith, B.E., Hornung, R.D., McQueen, D.M., & Peskin, C.S. (2007). An adaptive, formally second order accurate version of the immersed boundary method. Journal of Computational Physics, 223(1), 10-49.

[3] Griffith, B.E., Luo, X., McQueen, D.M., & Peskin, C.S. (2009). Simulating the fluid dynamics of natural and prosthetic heart valves using the immersed boundary method. International Journal of Applied Mechanics, 1(1), 137-177.

[4] Griffith, B.E., & Peskin, C.S. (2013). Electrophysiology. Communications on Pure and Applied Mathematics, 66(12), 1837-1913.

[5] McQueen, D.M., & Peskin, C.S. (2002). Heart simulation by an immersed boundary method with formal second-order accuracy and reduced numerical viscosity. In H. Aref and J.W. Phillips (Eds.), Mechanics for a New Millennium, Proceedings of the 20th International Congress on Theoretical and Applied Mechanics (pp. 429-444). Kluwer Academic Publishers.

[6] McQueen, D.M., O’Donnell T., Griffith, B.E., & Peskin, C.S. (2015). Constructing a patient-specific model heart from CT data. In N. Paragios, J. Duncan, and N. Ayache (Eds.), The Handbook of Biomedical Imaging: Methodologies and Clinical Research (pp. 183-198). New York, NY: Springer.

[7] Peskin, C.S. (1989). Fiber architecture of the left ventricular wall: an asymptotic analysis. Communications on Pure and Applied Mathematics, 42, 79-113.

[8] Peskin, C.S. (2002). The immersed boundary method. Acta Numerica, 11, 479-517.

[9] Peskin, C.S., & McQueen, D.M. (1994). Mechanical equilibrium determines the fractal fiber architecture of the aortic heart valve leaflets. American Journal of Physiology, 266, H319-H328.

[10] Stern, J.V., & Peskin, C.S. (1994). Fractal dimension of an aortic heart valve leaflet. Fractals – Complex Geometry Patterns and Scaling in Nature and Society, 2(3), 461-464. 

Debbie Sniderman is an applied physics engineer, materials scientist, and CEO of VI Ventures LLC, an engineering consulting company. She can be reached at

blog comments powered by Disqus