Skip to main content
Physics LibreTexts

Mathematical Continuum Descriptions of Membranes

  • Page ID
  • Membranes play an important role in cellular functions. They act as barriers to the outside environment, but also must allow the exchange of nutrients and waste. They must also be physically flexible to allow for cellular growth and movement, such as during cell division and migrations [1]. Membranes are composed of phospholipid bilayers. The phospholipids are amphipathic, containing a hydrophilic head and a hydrophobic tail. The membrane is composed of two adjacent layers of lipids, with the tails facing each other, and the hydrophilic heads facing the outside aqueous environment (Figure 1A). These membranes have a high aspect ratio with a low thickness of about 5 nm [1]. This is since the membrane is only two lipid molecules thick. However, the hydrophobic interior of the bilayer is a strong barrier against the transport of hydrophilic molecules across the membrane. Diffusion across the membrane is very slow for hydrophilic molecules and ions (transport requires transporter proteins) (Figure 1B). In addition, lipids can flip between leaflets, but this is very slow, as it requires the hydrophilic head pass through the hydrophobic region (Figure 1B). Lateral diffusion can be quite quick for lipids (~ 1 µm/sand even transmembrane proteins (0.1 µm/s2) (Figure 1B). Due to the described properties of the cellular membranes, they are commonly modeled as a 2-dimensional planar, elastic sheet in a 3-dimensional space (Figure 2). Biological structures like vesicles and even entire cells can be described as a closed form made by these sheets. These descriptions are known as continuum descriptions of the membrane. This is in contrast to molecular- and particle-based descriptions of the membrane, which focus on smaller time scales and length scales of membranes. 



    Figure 1: Structure of Biological Membranes (Adapted from [1]). A) The cell membrane is composed of a phospholipid bilayer. B) Lipids and proteins can undergo lateral diffusion, while hydrophilic molecules and ions require transporters. Lipids can also flip between leaflets.


    To understand the shape and dynamics of membranes, we will need to describe the geometry of the membrane surface itself. In addition, we will describe how the energy changes for different membrane configurations. In particular, we will look at how the energy changes when a membrane is stretched, bent, sheared and changed in the thickness (Figure 2).


    Figure 2: Biological membranes can be modeled as a two-dimensional sheet that stretches, bends, and changes in thickness. Adapted from [1].

    Area Changes:

    The area deformation can easily be parameterized by a single variable \(\Delta a\) or a function \(\Delta a(x, y) \) if the area deformation is dependent on position. We can treat the lipids in the membrane as a system of connected elastic springs, with increasing the area equivalent to stretching the springs (Figure 3).

    As such, the energy of stretching a piece of membrane is:

    \[E_{a}=\frac{1}{2} \frac{K_{a} \Delta a^{2}}{a_{0}}\]

    where \(K_a\) is the area modulus, which is intrinsic to the membrane, and \(a_0\) is some reference area from which change in area is calculated (i.e. \(\Delta a=a-a_{0}\)). This equation is equivalent to a Hookean spring. Of course, when looking at a whole system, this term has to be integrated over the total area:

    \[E_{a}=\frac{K_{a}}{2} \int (\frac{ \Delta a}{a_{0}})^2 da \]

    If the area change is constant over the surface of the membrane, then the integral reduces to Equation (1).

    The area modulus has values of about 55-70 kBT or 230-290 mN/m [1].


    Figure 3: The area change of a membrane can be modeled by treating the membrane as a network of springs. Adapted from [1].

    Thickness Changes

    Membrane compression can simply be described by a function \(w(x,y)\) that outputs the thickness at point \((x,y)\). Again, we can model the thickness change as a system of springs (Figure 4). If the half-width of the membrane (essentially the lipid length in the membrane at equilibrium) is defined as \(w_0\), and the thickness is \(2w_0\) (Figure 4), then we get a similar equation for the energy penalty of changing the thickness:

    \[E_{\text { thickness }}=\frac{K_{t}}{2} \int\left(\frac{w(x, y)-w_{0}}{w_{0}}\right)^{2} d a\]


    where \(K_t\) is the stiffness for thickness variation, which is intrinsic to the membrane. It has units of energy/area and is approximately 60 kBT/nm[1].


    Figure 4: The thickness change of a membrane can be modeled by treating the membrane as a network of springs. Adapted from [1].

    Shear Changes:

    Pure phospholipid membranes are essentially two-dimensional fluids. Due to the fluid nature of phospholipid bilayer, the membrane cannot support a shear strain and does not have a shear modulus. Therefore, as seen in Figure 2, static shear change does not deform a fluid membrane and does not cost any energy. However, most biological systems are not pure fluid membranes. For example, red blood cells have an underlying cytoskeletal network, which can support shear. If the shear elasticity is an important effect that needs to be modeled, it must be accounted for to model the equilibrium shapes and the dynamics of such cells. An example of this is described later.

    An Interlude on Membrane Geometry

    To understand how the energy costs of bending are determined, it is important to describe the geometry of the membrane.

    If we take a reference two-dimensional plane, we can define a function \(h(x,y)\) that characterizes the height of the membrane above the reference plane at point \((x,y)\). In most situations, this can accurately represent the shape of the membrane. This representation is known as the Monge parameterization (Figure 5). From the Monge parameterization, we can now calculate the curvature of a membrane. This is important as the free energy cost of bending the membrane is dependent on the curvature.


    Figure 5: Monge representation. Adapted from [1].

    In a single dimension, if we have a smooth function \(y=h(x)\), then at a point \(x\) we can define the curvature \(\kappa (x) \) to be the inverse of the radius of a circle that closely approximates the curve at point \(x\) (Figure 6A):


    It might be expected that a straight line has a curvature of 0. Indeed, since the circle that approximates the straight line is that of an infinite radius, the curvature is \( \kappa=\frac{1}{\infty}=0 \) (Figure 7).



    Figure 6: A) The curvature is defined as the inverse of the radius of the circle that best approximates the curve at a given point. B) For a surface, the curvature is defined by two values, which are the curvatures of the curve when the surface is cut by orthogonal planes.  Adapted from [1].

    Image result for circle of infinite radius

    Figure 7: A line has an infinite radius of curvature and therefore a curvature of zero.


    To generalize this to the two-dimensional case, we could cut through the two-dimensional surface with a plane, obtaining a curve, from which we can calculate the curvature. However, depending on the orientation of the plane used to cut the surface, the value of the curvature for that orientation is different. It can be shown that the curvatures in two directions describe a single point on a surface. Specifically, there is one particular choice of orthogonal planes that the curvatures take extreme values. These values \(c_1 , c_2 \) are known as the principal curvatures (Figure 6B). In fact, one can calculate the curvature in any other direction using the formula:

    \[c_{X}=c_{1} \cos ^{2} \alpha_{X}+c_{2} \sin ^{2} \alpha_{X}\]

    where \(c_{x}\) is the curvature in the direction X, and \(\alpha_{X}\) is the angle between the direction X and the direction belonging to curvature \(c_{1}\). These statements are known as the Euler Curvature Formula.

    As an example of principal curvatures, the principal curvatures of a cylinder would be \(1/R\) around the circumference of the cylinder (as the radius of curvature is \( R \), but \(0\) along the length of the cylinder.

    Other curvature quantities calculated from the principal curvatures are also commonly defined:

    \[ \begin{array}{c}{K=c_{1}+c_{2}} \\ {K_{G}=c_{1} c_{2}} \\ {H=\frac{K}{2}=\frac{c_{1}+c_{2}}{2}}\end{array} \]

    where \( K \) is the is the extrinsic curvature, \(G\) is the Gaussian curvature, and \(H\) is the mean curvature.

    With the Monge parameterization, the curvatures can be expressed in terms of the height function. It can be shown that the extrinsic curvature is:

    \[K=\nabla \cdot\left(\frac{\nabla h(\boldsymbol{r})}{\sqrt{1+(\nabla h(\boldsymbol{r}))^{2}}}\right)\]

    For a small gradient, the extrinsic curvature can be approximated as \( K \approx \Delta h(\boldsymbol{r}) \).

    Why is the Monge parameterization and expressing the curvature using the parameterization important? These are important as the bending energy is dependent on the curvature, and to determine the equilibrium shapes of a membrane system, the height function that minimizes the total free energy of the system will be determined to represent the final membrane shape.

    It is important to note that Monge parameterization is not the only way of describing a surface. It is not the most general parameterization or the most convenient parameterization in many cases. It will not work if there are folds or overlaps of the surface.

    This is only an introduction to membrane/surface geometry, which itself is only a small part of the field of differential geometry. For more information on differential geometry in relation to membranes, the reader is referred to [3]. For a more general and advanced introduction to differential geometry, the reader is referred to [4].

    Bending Changes – Helfrich theory:

    The bending energy of a membrane is described by the well-known Helfrich theory. The Helfrich bending energy is defined in terms of the extrinsic and Gaussian curvature:

    \[ E_{\text {bend}}=\int\left(\frac{1}{2} \kappa_{b}\left(c_{1}+c_{2}-c_{0}\right)^{2}+\overline{\kappa}_{b} c_{1} c_{2}\right) d A=\int\left(\frac{1}{2} \kappa_{b}\left(K-c_{0}\right)+\overline{\kappa}_{b} G\right) d A \]

    where \(\kappa_b\) is the bending rigidity and \(\overline{\kappa}_{b}\) is the Gaussian bending rigidity. This equation also contains the spontaneous curvature \(c_0\). It is nonzero if there is some curvature in the membrane at an equilibrium state, which occurs if there is an asymmetric bilayer (different lipid composition for each monolayer). Typical values of \(\kappa_b\) is 10-20 kBT [1].

    It is important to note that that Gaussian bending term can be ignored in most cases due to the Gauss-Bonnet theorem. The Gauss-Bonnet theorem says that if the topology of the system remains unchanged, the integral over the Gaussian curvature will be a constant. The topology (genus) can be determined by the number of holes/handles the surface has. For example, a sphere has no holes and therefore has a genus of 0, but a donut has one hole and a genus of 1. Therefore, unless the genus of the system changes, the Gaussian term vanishes. This is the case for many biological membranes. However, for complicated morphologies of the endoplasmic reticulum and mitochondria, Gaussian curvature could have a significant contribution as membrane holes, tubes, and tori are created and destroyed.

    We can simplify the above equation for Monge parameterization, with no spontaneous curvature. It can be shown for the Monge parameterization that the area element is \( d A=d x d y \sqrt{1+(\nabla h)^{2}} \). For the small gradient approximation, this is equivalent to \( d A=\left(1+\frac{1}{2}(\nabla h)^{2}\right) d x d y \) (based on a simple Taylor expansion). We can also penalize area stretching with a term \( \frac{1}{2} \sum(\nabla h)^{2} \) where \(\Sigma \) is membrane tension. Therefore, a Helfrich energy equation for a Monge parameterization, under the small gradient approximation, and includes bending and tension is:

    \[ E_{b e n d}=\frac{1}{2} \int d x d y\left[\kappa_{b}(\Delta h)^{2}+\Sigma(\nabla h)^{2}\right] \]

    Vesicles and Free Energy Cost

    Vesicles are spherical lipid bilayer structures. Vesicles are used for many cellular processes for transporting large amounts of cargo within cells or even between cells. An example of the latter is the transport of neurotransmitters between neurons at the synaptic cleft.

    We can calculate the free energy cost of the vesicles. The principal curvatures are the same: \( 1/R \) Therefore, the free energy cost is:

    \[ E_{v e s i c l e}=\frac{\kappa_{b}}{2} \int\left(\frac{2}{R}\right)^{2} d A=\frac{\kappa_{b}}{2}\left(\frac{2}{R}\right)^{2} 4 \pi R^{2}=8 \pi \kappa_{b} \]

    Interestingly, this means there is a fixed energy cost to create a vesicle regardless of the size of the vesicle. If \(\kappa_{b} \approx 10 k_{B} T \) then \( E_{\text {vesicle}} \approx 250 \mathrm{k}_{B} T \), which means there must be an active mechanism for the formation of vesicles in cells. Indeed vesicle-coating proteins like clathrin decrease this energy barrier such that vesicle formation is favorable.


    Membrane tethers:

    When a pulling force is applied to a single point on a vesicle, it will generate an extruded membrane tubule known as membrane tethers. Membrane tethers have many biological implications. For example, neutrophils can use membrane tethers to attach to the endothelial wall when rolling. Membrane tethers also have an important role in the endomembrane system. In addition, the energetics and dynamics of membrane tethers are similar to that of other biological processes, such as endocytosis, and actin projections during cell motility.

    Experiments use optical tweezers or molecular motors to apply a force at a single point on the membrane. We will use a simplified model (Figure 8) of the process where we have a sphere of radius \(R\) attached to a cylinder of radius \(r\) and length \(L\). At the end of the cylinder is a hemispherical cap.


    Figure 8: A simplified model of membrane tethers. Adapted from [1].

    We can write the free energy of the vesicle+tether system by calculating the contributions of bending, area stretching, and compression. The bending energy of the system is:

    \[ \begin{array}{c}{E_{\text {bending}}=8 \pi \kappa_{b}+\pi \kappa_{b} \frac{L}{r}+4 \pi \kappa_{b}} \\ {=12 \pi \kappa_{b}+\pi \kappa_{b} \frac{L}{r}}\end{array} \]

    The first term corresponds to the spherical vesicle (as demonstrated previously). The second term corresponds to the cylinder, which is calculated by squaring the sum of the principal curvatures \( \left(0+\frac{1}{r}\right)^{2} \)and multiplying by the area. The last term corresponds to the hemispherical cap.

    Area stretching also contributes to the free energy of the system:

    \[ E_{\text {area}}=\frac{K_{a}}{2} \frac{\left(a-a_{0}\right)^{2}}{a_{0}} \]

    where \( a_0 \) is the reference area and \( a=4 \pi R^{2}+2 \pi r L \) is the area of the system. The first term corresponds to a sphere of radius \( R\) and the second term corresponds to a cylinder of raidus \( r \). Note that the area of the hemispherical cap \( 4 \pi r^2 \) is ignored as it is assumed \( r \ll R, L \).

    Note that this system has no thickness change of the membrane and therefore does not contribute to the free energy of the system.

    In the case of a closed bilayer like a vesicle, the pressure difference across the membrane \( \Delta p \) also contributes to the energy of the system:

    \[ E_{p V}=-\Delta p \int d V=-\Delta p\left(\frac{4}{3} \pi R^{3}+r^{2} \pi L\right) \]

    Again, the contribution of the hemispherical cap is ignored. The minus sign is present because excess internal pressure increases the volume.

    Finally, the force pulling the membrane must also contribute to the free energy of the system:

    \[ E_{l o a d}=-f L \]

    where \( f \) is the force applied. The minus sign is present because increased force leads to a high tether length.

    Summing equations 11, 12, 13, and 14 together, we get:

    \[ E_{t o t}=12 \pi \kappa_{b}+\pi \kappa_{b} \frac{L}{r}+\frac{K_{a}}{2} \frac{\left(a-a_{0}\right)^{2}}{a_{0}}-\Delta p\left(\frac{4}{3} \pi R^{3}+r^{2} \pi L\right)-f L \]

    At equilibrium, the total energy is minimized. Therefore, the equilibrium shape can be determined by minimizing the free energy.

    In this case, we can take the derivatives of equation 15 with respect to each of the variables to minimize the energy.

    Note that

    \[ \frac{\partial E_{\text {area}}}{\partial R}=K_{a} \frac{a-a_{0}}{a_{0}} \frac{\partial a}{\partial R}=\tau 8 \pi R \]

    where \( \tau \) is the membrane surface tension, which can be related to area stretching by the formula \( \tau = \frac{K_{a}\left(a-a_{0}\right)}{a_{0}} \). Similarly, it can be shown that \( \partial E_{\text {area}} / \partial r=\tau 2 \pi L \) and \( \partial E_{\text {area}} / \partial L=\tau 2 \pi r \). Therefore:

    \[ \begin{array}{l}{\frac{\partial E_{t o t}}{\partial R}=0 \rightarrow 8 \pi \tau R-4 \pi \Delta p R^{2}=0 \rightarrow \Delta p=\frac{2 \tau}{R}} \\ {\frac{\partial E_{t o t}}{\partial r}=0 \rightarrow-\pi \kappa_{b} \frac{L}{r^{2}}+2 \pi \tau L-2 \pi \Delta p r L=0} \\ {\frac{\partial E_{t o t}}{\partial L}=0 \rightarrow \pi \kappa_{b} \frac{1}{r}+2 \pi \tau r-\pi \Delta p r^{2}-f=0}\end{array} \]

    The first minimization leads to an equation known as the Laplace-Young relation, which is a well-known equation that relates the pressure difference \( \Delta p \) across a fluid interface (here the membrane) and the surface tension \( \tau \). A more general form takes into account both radii of curvature:

    \[ \Delta p=\tau\left(\frac{1}{R_{1}}+\frac{1}{R_{2}}\right) \]

    We can plug in the Laplace-Young relation into the other minimized equations:

    \[ \begin{array}{l}{-\pi \kappa_{b} \frac{L}{r^{2}}+2 \pi \tau L-\frac{4 \pi \tau r L}{R}=0} \\ {\pi \kappa_{b} \frac{1}{r}+2 \pi \tau r-f-\frac{2 \pi \tau r^{2}}{R}=0}\end{array} \]

    Since \( r \ll R \) (Figure 8) the last terms can be ignored:

    \[ 2 \pi \tau L=\pi \kappa_{b} \frac{L}{r^{2}} \]

    \[ r=\sqrt{\frac{\kappa_{b}}{2 \tau}} \]

    \[ \pi \kappa_{b} \frac{1}{r}+2 \pi \tau r-f=0 \rightarrow \pi \kappa_{b} \sqrt{\frac{2 \tau}{\kappa_{b}}}+2 \pi \tau \sqrt{\frac{\kappa_{b}}{2 \tau}}=f \]

    \[ f=2 \pi \sqrt{2 K_{b} \tau} \]

    We now get a relation between the radius of tether, bending energy, and membrane tension. In addition, we get a relation between the force applied, bending energy and tension. Similar relations can actually be used to determining the bending rigidity \( \kappa_{b} \) using a micropipette aspiration and optical tweezers used to make membrane tethers [5].

    A general framework of determining equilibrium shapes (variational calculus):

    The previous section looked at how relationships between meaningful physical quantities can be derived if the shape of the system is known. This was done by calculating the free energy cost for the membrane shape and minimizing it with respect to the variables of the system. However, a different approach can also be taken. What if the physical quantities of the system are known, but the equilibrium shapes need to be determined? This can be done by minimizing an energy functional of the system, which is the basis of variational calculus.

    We can define a functional as a function that assigns a scalar value for a given function. In this case, we can think of the free energy of the system \( E_{tot} \) as a functional that is dependent on the height function \( h(x,y) \). Variational calculus provides tools to find functions that minimize the functionals. This is exactly what we need to determine the equilibrium shapes! Variational calculus can be used to find the shape function \( h(x,y) \) that minimizes the energy \( E_{tot} \)

    In variational calculus, in a one-dimensional case, the minimization of the functional:

    \[ J[y]=\int L\left(x, y, y^{\prime}\right) d x \]

    is solved by the solution of the differential equation:

    \[ \frac{\partial L}{\partial f}-\frac{d}{d x} \frac{\partial L}{\partial f^{\prime}}=0 \]

    where \( y=f(x) \) minimizes the functional. This equation is known as the Euler-Lagrange equation. Indeed, the equation for \( E_{tot} \) has this form.

    As a side note, many equations in physics have this form and the Euler-Lagrange equation is used. In fact, the dynamics of a physical system can be represented as a functional known as the Lagrangian, and solving the Euler-Lagrange equations for this Lagrangian lead to the equations of motion. This approach to classical mechanics is known as Lagrangian mechanics.

    Going back to membranes, if we use Equation (9), we can get the following Euler-Lagrange equation:

    \[ \kappa \Delta \Delta h-\Sigma \Delta h=0 \]

    This partial differential equation can be solved for the height function that describes the shape of the membrane at equilibrium. This is called the shape equation.

    Minimization of red blood cells:

    While solving the shape equations are non-trivial and won’t be shown here, some results in the case of red blood cells will instead be shown.

    It is important to note that there are often two additional energy terms added to the Helfrich bending energy to obtain the total free energy for red blood cells.

    When an initially flat membrane is bent, the outer monolayer is stretched, and the inner monolayer is compressed and the number of lipids per unit area is higher for the inner surface and lower for the outer surface. The energy cost of this is added to the model as the following model:

    \[ E_{A D E}=\frac{\kappa_{A D E}}{2} \frac{\pi}{A D^{2}}\left(\Delta A-\Delta A_{0}\right)^{2} \]

    where \(\kappa_{ADE} \) is a separate bending energy, which a material parameter, \(A\) is the membrane area, and \(D\) is the membrane thickness. This is known as the area-difference-elasticity (ADE) model.

    In addition, an elastic term is added of the form:

    \[ E_{e}=\frac{K_{A}}{2} \oint \alpha^{2} d A+\mu \oint \beta d A \]

    where \( \alpha^{2} \) and \(\beta^{2} \) are the local area and strain invariants, \( K_{A} \) and \( \mu \) are elastic moduli for stretching and strain, respectively.  

    When the total free energy \( E_{b e n d}+E_{A D E}+E_{e} \) is minimized, the equilibrium shapes of the red blood cells can be determined. The below figure shows a comparison of the theoretically determined shapes and those that have been found experimentally.


    Figure 9: A comparison of experimentally observed red blood cell shapes (left side) and theoretically determined equilibrium shapes (right side). Adapted from [7]


    1. Phillips, R.; Kondev, J.; Theriot, J.; Garcia, H. G. Physical Biology of the Cell, 2013, Garland Science, Taylor & Francis Group, LLC: Abingdon, UK.
    2. Boal, D. Mechanics of the Cell. Cambridge University Press: Cambridge, UK. 2002
    3. Deserno, Markus. Fluid Lipid Membranes – A Primer. Max-Planck-Institut fu ̈r Polymerforschung. Mainz, Germany. 2006
    4. Deserno, Markus. Notes on Differential Geometry. Max-Planck-Institut fu ̈r Polymerforschung. Mainz, Germany. 2004
    5. Bo, Lin, and Richard E. Waugh. “Determination of Bilayer Membrane Bending Stiffness by Tether Formation from Giant, Thin-Walled Vesicles.” Biophysical Journal, vol. 55, doi:10.1016/S0006-3495(89)82844-9. Accessed 13 June 2019.
    6. Li, Xuejin, et al. Continuum- and Particle-Based Modeling of Shapes and Dynamics of Red Blood Cells in Health and Disease. Soft Matter, vol. 9, no. 1, 2013, pp. 28–37., doi:10.1039/c2sm26891d.

    7. W., G. Lim H., et al. “Stomatocyte-Discocyte-Echinocyte Sequence of the Human Red Blood Cell: Evidence for the Bilayer- Couple Hypothesis from Membrane Mechanics.” Proceedings of the National Academy of Sciences, vol. 99, no. 26, 2002, pp. 16766–16769., doi:10.1073/pnas.202617299.