← Allometric Scaling Next: Morphogenesis →

Module 5: Minimal Surfaces & Soap Films

A soap film stretched across a wire is the physical instantiation of a deep mathematical problem: find the surface of minimum area bounded by a given curve. Plateau's 19th-century experiments revealed universal geometric rules obeyed by all soap-film junctions — rules that also govern cell membranes, radiolarian skeletons, beetle cuticle photonic structures, and the biconcave red blood cell. This module derives the minimal-surface equation (\(H=0\)) from an area functional, presents Plateau's laws, and shows how the Helfrich bending energy explains the RBC shape and how Schwarz P/D/G triply periodic minimal surfaces generate iridescent biological nanostructures.

1. Minimal Surfaces

A minimal surface is a smooth surface that locally minimises area. Given a domain \(D\subset\mathbb{R}^2\) and a parametrisation\(\mathbf{r}(u,v)\), the area functional is:

\[ A[\mathbf r] = \iint_D \sqrt{EG - F^2}\,du\,dv \]

\(E, F, G\) are the first fundamental form coefficients.

1.1 Euler-Lagrange Derivation

For a graph \(z = u(x,y)\), the area functional becomes:

\[ A[u] = \iint \sqrt{1 + u_x^2 + u_y^2}\,dx\,dy \]

The Euler-Lagrange equation gives the minimal surface equation:

\[ (1+u_y^2) u_{xx} - 2 u_x u_y u_{xy} + (1+u_x^2) u_{yy} = 0 \]

Equivalent to \(H = 0\): the mean curvature vanishes everywhere.

1.2 The Mean Curvature Vanishes

The mean curvature \(H = (\kappa_1 + \kappa_2)/2\) is the average of the two principal curvatures. For a minimal surface, \(\kappa_1 = -\kappa_2\): every point is a saddle. This is because the Laplace-Young pressure jump\(\Delta P = 2\gamma H\) vanishes for a soap film of surface tension \(\gamma\)spanning two points at equal pressure.

1.3 Classical Examples

Plane

Trivial; H = 0 everywhere.

Catenoid

Only surface of revolution that is minimal; r(z)=a cosh(z/a).

Helicoid

Ruled minimal surface; related to catenoid by Bonnet transformation.

Enneper surface

Self-intersecting; first discovered.

Scherk surface

Periodic; z=ln(cos y/cos x).

Schwarz P/D/G

Triply periodic; biology's favorite.

2. Plateau's Laws

Joseph Plateau (1873) catalogued the rules that soap films obey at junctions. Jean Taylor (1976) proved them rigorously as consequences of area minimisation:

Smooth surfaces: each film piece is a smooth surface of constant mean curvature (vanishing for open films; non-zero in bubble clusters).
Edge rule: films meet three at a time, along edges, at angles of exactly \(120°\).
Vertex rule: edges meet four at a time at vertices, at the tetrahedral angle \(\arccos(-1/3) = 109.47°\).
Pressure rule: at every interface,\(\Delta P = 2\gamma H\) (Young-Laplace).

2.1 Derivation of the 120° rule

At an edge where three films meet, the surface tensions \(\mathbf T_1, \mathbf T_2, \mathbf T_3\)(each equal in magnitude to \(\gamma\)) must balance. Setting\(\mathbf T_1 + \mathbf T_2 + \mathbf T_3 = 0\) for three co-planar unit vectors of equal magnitude yields:

\[ \cos\theta = -\tfrac{1}{2} \;\Rightarrow\; \theta = 120° \]

2.2 Derivation of the Tetrahedral Vertex Angle

Four edges meet at a vertex. The four edge-tension vectors must sum to zero. By symmetry (Platonic tetrahedron) the angle between any two edges is \(\cos^{-1}(-1/3) \approx 109.47°\).

2.3 Biological Relevance

Plateau's 120° rule determines the structure of honeybee combs (Module 7), epithelial cell sheets, and foam-based lightweight biological materials. The tetrahedral vertex is found in the 4-cell stage of embryonic development and in skeletal radiolarian junctions.

3. Radiolaria & Haeckel's Forms

Radiolaria are single-celled marine protists with elaborate silica skeletons. Ernst Haeckel'sKunstformen der Natur (1904) illustrated thousands of species. D'Arcy Thompson (1917) and later Thompson's followers recognised that these skeletons are close approximations to minimal surfaces spanning polyhedral frames — soap-film mathematics at the level of a single cell.

3.1 The Aulonia hexagona Problem

Haeckel described species (e.g. Aulonia hexagona) whose skeletons appeared to be pure hexagonal tilings of the sphere. By Euler's formula, this is impossible:

\[ V - E + F = 2 \]

Euler's polyhedron formula rules out pure hexagonal sphere tilings — at least 12 pentagons are required (as in a football).

Close inspection of Aulonia skeletons confirms this: they contain exactly 12 pentagonal “defects” among mostly hexagonal cells. Biology obeys the theorem.

4. Triply Periodic Minimal Surfaces

Schwarz (1890), Riemann and Schoen (1970) discovered families of minimal surfaces that are periodic in all three spatial directions. They divide space into two interpenetrating labyrinths of equal volume. The simplest explicit approximations use nodal surfaces of trigonometric functions:

\( \text{Schwarz P: }\; \cos x + \cos y + \cos z = 0 \)

\( \text{Schwarz D: }\; \sin x\sin y\sin z + \sin x\cos y\cos z + \cos x\sin y\cos z + \cos x\cos y\sin z = 0 \)

\( \text{Schoen G: }\; \sin x\cos y + \sin y\cos z + \sin z\cos x = 0 \)

4.1 In Biology

Butterfly wings

Papilio palinurus: gyroid (G) surfaces of chitin produce iridescent green; periodicity ~300 nm.

Beetle cuticle

Scarab beetles (Lamprocyphus augustus): photonic gyroid from self-assembled block copolymers.

Mitochondrial cristae

Cubic membrane phases observed in amoeba starved mitochondria (Landh 1995).

Block copolymers

Synthetic analogues of biological TPMS self-assemble in lab at 30-50 nm scale.

Chloroplast thylakoids

Stacked lamellar sheets interconvertible with cubic (Schwarz-D) phases.

Bone trabecula

Approximate TPMS geometry optimises stiffness/mass (gyroid lattice).

5. Helfrich Bending Energy

Lipid bilayer membranes are nearly incompressible. Their shape is dominated by bending. Helfrich (1973) postulated the membrane energy functional:

\[ E = \frac{\kappa}{2}\oint (2H - c_0)^2\,dA \;+\; \bar\kappa\oint K\,dA \;+\; \sigma\oint dA \;+\; \Delta P\int dV \]

\(\kappa \sim 10\,k_B T\) (bilayer bending modulus);\(c_0\) = spontaneous curvature;\(\bar\kappa\) = Gaussian bending modulus;\(\sigma\) = tension.

5.1 Gauss-Bonnet

By the Gauss-Bonnet theorem, the Gaussian curvature integral\(\oint K dA = 4\pi(1-g)\) depends only on the topological genus \(g\). For closed non-fissioning vesicles, this term is a topological constant and can be dropped.

5.2 Shape Equation

Varying \(E\) at fixed area and volume leads to the shape equation (Ou-Yang & Helfrich 1989):

\[ \Delta P - 2\sigma H + \kappa\bigl[\nabla^2(2H) + (2H - c_0)(2H^2 + c_0 H - 2K)\bigr] = 0 \]

5.3 Red Blood Cell Biconcave Shape

Canham (1970) and Evans & Fung (1972) found experimentally that the RBC has a reduced volume\(v = V/[(4\pi/3)(A/4\pi)^{3/2}] \approx 0.59\) (highly oblate). Minimising Helfrich energy at this volume-area constraint reproduces the observed biconcave disc:

\( z(\rho) = \pm \tfrac{1}{2} D_0\sqrt{1-\rho^2/R_0^2}(C_0 + C_2\rho^2/R_0^2 + C_4\rho^4/R_0^4) \)

D0 = 2.84 μm, R0 = 3.91 μm, C0 = 0.21, C2 = 2.03, C4 = -1.12 (Evans & Fung 1972).

This shape optimises: (a) surface-to-volume ratio for O\(_2\) exchange, (b) elastic deformability through capillaries of diameter \(\sim 3\,\mu\text{m}\) (smaller than the RBC diameter of 8 μm). It is also a stationary point of the Helfrich functional with\(c_0 \approx 0\) and reduced volume 0.59.

6. Classical Minimal Surfaces

7. Red Blood Cell & Plateau Junction

8. Simulation 1 — Plateau's Problem

Numerically solve for the catenoid spanning two parallel rings separated by\(2h\) with rim radius \(R\). There are two solutions (stable + unstable) for \(h/R < 0.528\); above this, the soap film jumps to the “Goldschmidt discontinuous” minimiser: two separate discs.

Python

script.py104 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# Numerical minimisation of surface area for a catenoid spanning two rings
# Analytic solution: z in [-h, h], rim radius R -> catenary profile
# r(z) = a*cosh(z/a) with a determined by a*cosh(h/a) = R

R_rim = 1.0
h = 0.6

# Solve transcendental a*cosh(h/a) = R
from scipy.optimize import brentq
def eq(a):
    return a*np.cosh(h/a) - R_rim

# Two possible roots (stable + unstable) for h < h_crit approx 0.663 R
a_small = brentq(eq, 0.01, 0.55)
a_big = brentq(eq, 0.56, 2.0)

# Surface areas
def area_catenoid(a, h):
    # A = 2*pi * integral r * sqrt(1 + r'^2) dz
    # r = a cosh(z/a), r' = sinh(z/a), sqrt(1+r'^2) = cosh(z/a)
    # A = 2*pi * a * integral cosh^2(z/a) dz from -h to h
    return 2*np.pi*a * (h + a*np.sinh(h/a)*np.cosh(h/a))

A_small = area_catenoid(a_small, h)
A_big = area_catenoid(a_big, h)
# Competing "Goldschmidt" solution: two separate discs
A_disk = 2*np.pi*R_rim**2
print(f"a_small = {a_small:.4f}, Area = {A_small:.4f}")
print(f"a_big   = {a_big:.4f}, Area = {A_big:.4f} (unstable)")
print(f"Two-disc (Goldschmidt) = {A_disk:.4f}")

# Plot profiles
fig = plt.figure(figsize=(18, 7))
fig.patch.set_facecolor('#030712')
fig.suptitle("Plateau's Problem: Catenoid vs Two-Disc Minimum",
             fontsize=15, color='#fde68a', fontweight='bold')

z = np.linspace(-h, h, 100)
r_small = a_small * np.cosh(z/a_small)
r_big = a_big * np.cosh(z/a_big)

ax1 = fig.add_subplot(1, 3, 1)
ax1.set_facecolor('#1a1200')
ax1.plot(r_small, z, color='#fbbf24', linewidth=2.5, label=f'stable, A={A_small:.3f}')
ax1.plot(-r_small, z, color='#fbbf24', linewidth=2.5)
ax1.plot(r_big, z, '--', color='#f87171', linewidth=2, label=f'unstable, A={A_big:.3f}')
ax1.plot(-r_big, z, '--', color='#f87171', linewidth=2)
ax1.axhline(h, color='#22d3ee', linewidth=1)
ax1.axhline(-h, color='#22d3ee', linewidth=1)
ax1.set_xlabel('r', color='white'); ax1.set_ylabel('z', color='white')
ax1.set_title('Catenoid profiles', color='#fde68a')
ax1.legend(facecolor='#1a1200', edgecolor='#fbbf24', labelcolor='white', fontsize=9)
ax1.set_aspect('equal')
ax1.tick_params(colors='white')
for sp in ax1.spines.values(): sp.set_color('#fbbf24')

# 3D catenoid
ax2 = fig.add_subplot(1, 3, 2, projection='3d')
ax2.set_facecolor('#1a1200')
theta = np.linspace(0, 2*np.pi, 60)
Z, TH = np.meshgrid(z, theta)
Rm = a_small * np.cosh(Z/a_small)
X = Rm*np.cos(TH); Y = Rm*np.sin(TH)
ax2.plot_surface(X, Y, Z, cmap='inferno', alpha=0.8, edgecolor='none')
ax2.set_title('Catenoid surface (stable)', color='#fde68a')
ax2.set_xlabel('x', color='white'); ax2.set_ylabel('y', color='white'); ax2.set_zlabel('z', color='white')
ax2.tick_params(colors='white')

# Area versus h/R phase diagram
ax3 = fig.add_subplot(1, 3, 3)
ax3.set_facecolor('#1a1200')
h_vals = np.linspace(0.05, 0.7, 60)
A_cat_s = []; A_cat_b = []; A_disc = []
for hv in h_vals:
    try:
        a_s = brentq(lambda a: a*np.cosh(hv/a) - 1.0, 0.01, 0.55)
        a_b = brentq(lambda a: a*np.cosh(hv/a) - 1.0, 0.56, 2.0)
        A_cat_s.append(area_catenoid(a_s, hv))
        A_cat_b.append(area_catenoid(a_b, hv))
        A_disc.append(2*np.pi*1.0**2)
    except Exception:
        A_cat_s.append(np.nan); A_cat_b.append(np.nan); A_disc.append(2*np.pi)
ax3.plot(h_vals, A_cat_s, color='#fbbf24', linewidth=2, label='catenoid (stable)')
ax3.plot(h_vals, A_cat_b, '--', color='#f87171', linewidth=2, label='catenoid (unstable)')
ax3.plot(h_vals, A_disc, color='#22d3ee', linewidth=2, label='two discs')
ax3.axvline(0.528, color='#86efac', linestyle=':', linewidth=2,
            label='jump to Goldschmidt')
ax3.set_xlabel('h/R', color='white'); ax3.set_ylabel('area', color='white')
ax3.set_title('Phase diagram of soap-film minima', color='#fde68a')
ax3.legend(facecolor='#1a1200', edgecolor='#fbbf24', labelcolor='white', fontsize=8)
ax3.grid(True, alpha=0.2, color='#fbbf24')
ax3.tick_params(colors='white')
for sp in ax3.spines.values(): sp.set_color('#fbbf24')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#030712')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Simulation 2 — Helfrich Shape Families

A sketch of the phase diagram of vesicle shapes: sphere, prolate, oblate, and biconcave. At reduced volume \(v \approx 0.59\), the biconcave RBC shape minimises the Helfrich bending energy — matching human erythrocytes.

Python

script.py75 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Helfrich bending energy for axisymmetric vesicles (Seifert 1997 review)
# E = (kappa/2) * int (2H - c0)^2 dA + Delta P V - sigma A
# Solve Euler-Lagrange equations for axisymmetric shapes.
# We parametrize by arclength s, tangent psi(s), radius r(s), height z(s).

kappa = 1.0    # bending modulus (units: k_B T)
c0 = 0.0       # spontaneous curvature
sigma_t = 0.0  # tension
dp = 0.0       # pressure difference

# Use stripped shape-equation in terms of psi(s)
# Shape equation for axisymmetric vesicle (Deuling-Helfrich 1976):
# d/ds (sin(psi)/r) * ... -> instead, use direct ODE from Jenkins 1977:
# We'll illustrate the family of solutions: prolate, oblate, biconcave (RBC)
# via the reduced volume v = V/((4pi/3)(A/(4pi))^(3/2))

# Heuristic parametric shapes
def shape(kind, n=200):
    s = np.linspace(0, 1, n)
    if kind == 'sphere':
        theta = s*np.pi
        r = np.sin(theta); z = 1 - np.cos(theta)
    elif kind == 'prolate':
        theta = s*np.pi
        r = 0.55*np.sin(theta); z = 1 - np.cos(theta)
    elif kind == 'oblate':
        theta = s*np.pi
        r = 1.2*np.sin(theta); z = 0.6*(1 - np.cos(theta))
    elif kind == 'biconcave':
        # RBC-like formula (Evans & Fung 1972)
        rh = np.linspace(0, 1, n)
        # z = +/- 0.5*D0*sqrt(1-rho^2/R^2)*(C0 + C2*rho^2 + C4*rho^4)
        D0 = 0.8; C0 = 0.21; C2 = 2.03; C4 = -1.12
        rho = rh
        with np.errstate(invalid='ignore'):
            z = 0.5*D0*np.sqrt(np.clip(1 - rho**2, 0, 1))*(C0 + C2*rho**2 + C4*rho**4)
        r = rh; z = z
    return r, z

fig = plt.figure(figsize=(20, 6))
fig.patch.set_facecolor('#030712')
fig.suptitle("Helfrich Membrane Shapes: The RBC Problem",
             fontsize=15, color='#fde68a', fontweight='bold')

shapes = ['sphere', 'prolate', 'oblate', 'biconcave']
for i, k in enumerate(shapes):
    ax = fig.add_subplot(1, 4, i+1)
    ax.set_facecolor('#1a1200')
    r, z = shape(k)
    # Mirror to make full cross-section
    ax.plot(r, z, color='#fbbf24', linewidth=2.5)
    ax.plot(-r, z, color='#fbbf24', linewidth=2.5)
    if k == 'biconcave':
        ax.plot(r, -z, color='#fbbf24', linewidth=2.5)
        ax.plot(-r, -z, color='#fbbf24', linewidth=2.5)
    ax.set_title(f'{k}', color='#fde68a')
    ax.set_aspect('equal')
    ax.tick_params(colors='white')
    for sp in ax.spines.values(): sp.set_color('#fbbf24')
    # Curvatures & energy computation (schematic)
    H_mean = 1.0/0.5  # not exact, schematic
    ax.text(0.02, 0.98, f'v = {0.95 if k=="sphere" else 0.8 if k=="prolate" else 0.85 if k=="oblate" else 0.59}',
            transform=ax.transAxes, color='#fde68a', fontsize=9, va='top')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#030712')
print("Biconcave RBC shape minimises bending energy at reduced volume v approx 0.59")
print("(Deuling & Helfrich 1976; Canham 1970)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

10. Simulation 3 — Schwarz P / D / G

Marching-cubes extraction of the triply periodic minimal surfaces Schwarz P, Schwarz D and Schoen G (the gyroid). The gyroid appears in beetle exoskeletons and butterfly wing scales as a photonic crystal (Saranathan et al. 2010).

Python

script.py57 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# Schwarz P triply periodic minimal surface (approximation): cos(x)+cos(y)+cos(z) = 0
# Construct by marching cubes on scalar field

from skimage import measure

N = 50
L = 2*np.pi
x = np.linspace(-L, L, N)
y = np.linspace(-L, L, N)
z = np.linspace(-L, L, N)
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

# Schwarz P
F_P = np.cos(X) + np.cos(Y) + np.cos(Z)
# Schwarz D
F_D = (np.sin(X)*np.sin(Y)*np.sin(Z)
       + np.sin(X)*np.cos(Y)*np.cos(Z)
       + np.cos(X)*np.sin(Y)*np.cos(Z)
       + np.cos(X)*np.cos(Y)*np.sin(Z))
# Schoen G (Gyroid)
F_G = np.sin(X)*np.cos(Y) + np.sin(Y)*np.cos(Z) + np.sin(Z)*np.cos(X)

def plot_surface(ax, F, title, color):
    verts, faces, _, _ = measure.marching_cubes(F, level=0.0,
                         spacing=(x[1]-x[0], y[1]-y[0], z[1]-z[0]))
    ax.plot_trisurf(verts[:,0]-L, verts[:,1]-L, faces, verts[:,2]-L,
                     color=color, alpha=0.85, edgecolor='none',
                     linewidth=0)
    ax.set_title(title, color='#fde68a')
    ax.set_facecolor('#1a1200')
    ax.tick_params(colors='white')
    ax.set_box_aspect([1,1,1])

fig = plt.figure(figsize=(20, 7))
fig.patch.set_facecolor('#030712')
fig.suptitle("Triply Periodic Minimal Surfaces (Schoen 1970)",
             fontsize=15, color='#fde68a', fontweight='bold')

ax1 = fig.add_subplot(1, 3, 1, projection='3d')
plot_surface(ax1, F_P, 'Schwarz P', '#fbbf24')
ax2 = fig.add_subplot(1, 3, 2, projection='3d')
plot_surface(ax2, F_D, 'Schwarz D', '#a78bfa')
ax3 = fig.add_subplot(1, 3, 3, projection='3d')
plot_surface(ax3, F_G, 'Schoen G (gyroid)', '#86efac')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#030712')
print("Schwarz P, D and G all found in butterfly wings, beetle cuticle")
print("(Saranathan et al. 2010)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Module Summary

Minimal surface equation

H = 0 everywhere; Euler-Lagrange of area functional.

Plateau's laws

120° edges, 109.47° vertices; universal soap-film rules.

Catenoid / Helicoid

Canonical pair; related by Bonnet transformation.

Radiolaria

Aulonia hexagona obeys Euler's polyhedron formula.

Schwarz P / D / gyroid

TPMS in butterfly wings, beetle cuticle, mitochondria.

Helfrich energy

E = (kappa/2) int (2H-c0)^2 dA; kappa ~ 10 k_B T.

RBC biconcave shape

Reduced volume 0.59 minimises Helfrich energy.

Gauss-Bonnet

Int K dA = 4pi(1-g) is topological; drops for fixed topology.

References

Plateau, J. (1873). Statique expérimentale et théorique des liquides soumis aux seules forces moléculaires.
Taylor, J. E. (1976). The structure of singularities in soap-bubble-like and soap-film-like minimal surfaces. Ann. Math., 103, 489–539.
Canham, P. B. (1970). The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell. J. Theor. Biol., 26, 61–81.
Helfrich, W. (1973). Elastic properties of lipid bilayers: theory and possible experiments. Z. Naturforsch., 28c, 693–703.
Evans, E. & Fung, Y.-C. (1972). Improved measurements of the erythrocyte geometry. Microvasc. Res., 4, 335–347.
Ou-Yang, Z.-C. & Helfrich, W. (1989). Bending energy of vesicle membranes. Phys. Rev. A, 39, 5280.
Schoen, A. H. (1970). Infinite periodic minimal surfaces without self-intersections. NASA Tech. Note D-5541.
Saranathan, V. et al. (2010). Structure, function, and self-assembly of single-network gyroid (I4132) photonic crystals in butterfly wing scales. PNAS, 107, 11676.
Seifert, U. (1997). Configurations of fluid membranes and vesicles. Adv. Phys., 46, 13–137.
Thompson, D'A. W. (1917). On Growth and Form. Cambridge University Press.
Haeckel, E. (1904). Kunstformen der Natur.
Landh, T. (1995). From entangled membranes to eclectic morphologies: cubic membranes as subcellular space organizers. FEBS Lett., 369, 13–17.

← Module 4: Allometric Scaling Next: Module 6 — Morphogenesis →