Lipid bilayers constitute one of the critical parts of all biological membranes, including cell membranes. A nice description of lipid bilayers and their function in biological membranes can be found here. They can be exceptionally complex and contain hundreds of different constitutents, so simpler model lipid bilayers are often produced in the laboratory and studied experimentally. They form closed spheroidal structures, called liposomes, with a thickness of a few nm, and characteristic linear dimensions up to several microns. Larger such structures are usually called Giant Unilamellar Vesicles, or GUVs. Why should we care about these structures as mechanicians? For a number of reasons, the elastic properties of lipid membranes are thought to play a crucial role in governing their potential configurations. Recent experimental studies of the role of membrane curvature on domain formation in biomembranes, for example, provide testament to this notion. Images of their work are reproduced (with permission) here.

Although one could certainly debate the applicability of continuum models to structures such as these, in fact a great deal of insight has been provided by the same. Further, while molecular dynamics and other discrete methods certainly have a lot to offer toward understanding these structures, even coarse-grained methods are impractical for handling the billions of lipids present in GUVs.

Not that things are much easier for continuum models! One of the main issues for simulations based on these concerns the relatively high order of even the simplest evolution equations. For example, the classical bilayer mechanics theory developed by Helfrich (see also the work of Canham and Evans) involves a bending energy that is quadratic in the principle curvatures of the membrane surface. This gives rise to Euler-Lagrange equations that are fourth-order. Proposed models for surface composition, akin to Cahn-Hilliard equations, are similarly fourth-order in the composition field. Given the complex geometry and possible configurations attainable with giant unilamellar vesicles, classical high-order finite difference and spectral methods are essentially impractical.

Classical finite-element methods are also problematic for fourth-order problems. The parent space of functions in this case is the second-order Sobolev space H2(S), which, roughly speaking, requires the use of smooth C1 shape functions on the membrane surface. With good reason, recent work by Feng and Klug (2) has moved away from classical shell finite elements and toward C1-conforming subdivision surfaces. Other recent work incorporating smooth basis functions includes that of Ayton et al. (1), who relied on a smoothed-particle hydrodynamics (SPH)-based method.

Although these approaches have their advantages, neither is capable of reasonably simulating the topology changes (e.g. sphere to toroid) required to predict equilibrium configurations in GUVs (though Ayton may claim otherwise). The only method with such capability appears to be that of Qiang Du and coworkers (3), employing a phase-field regularization of (essentially) the Canham-Helfrich-Evans biomembrane theory. This treatment is appealing, but it is not a particularly efficient approach, as it requires the use of a uniform Cartesian grid and embeds the interface within a higher-dimensional domain. Further, while their most recent work has allowed for spatial variations in elastic moduli, they have yet to couple their approach with composition dynamics.

The three papers included here for discussion are:

(1) G. S. Ayton, J. L. McWhirter, P. McMurthy, & G. A. Voth, 2005. Coupling Field Theory with Continuum Mechanics: A Simulation of Domain Formation in Giant Unilamellar Vesicles, Biophysical Journal, vol. 88, 3855--3869.

(2) F. Feng & W.S. Klug, 2006. Finite element modeling of lipid bilayer membranes, Journal of Computational Physics, vol. 220, 394-408.

(3) Q. Du, C. Liu, & X. Wang, 2006. Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions, Journal of Computational Physics, vol. 212, 757--777.

I recommend beginning with the Klug paper, as it provides some excellent background and a nice introduction to the bending aspects. The Ayton paper touches on the coupled chemo-mechanical aspects and the formation of domains in GUVs. Finally, the Du paper is a very different take on this problem, and one that I believe to be of interest. I welcome discussion on these papers, and will focus my responses on the computational aspects.