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

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.

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.

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

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).

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.

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.

References

[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.