← Minimal Surfaces Next: Packing & Tessellation →

Module 6: Morphogenesis & Turing Patterns

In 1952 Alan Turing asked: how can a uniform ball of cells spontaneously organise into an animal? His answer — a reaction-diffusion instability driving a uniform state into spatial patterns of stripes and spots — is now the paradigmatic mechanism of morphogenesis. This module derives Turing's instability condition, derives morphogen gradients (Wolpert's French-flag model, the Bicoid gradient in Drosophila), and shows how differential adhesion produces cell sorting. Zebrafish stripes, leopard spots, and seashell pigmentation all emerge from these same mathematical rules.

1. Turing's 1952 Insight

In “The Chemical Basis of Morphogenesis” (Phil. Trans. R. Soc. B, 1952), Alan Turing considered a system of two interacting chemicals — an activator \(u\) and an inhibitor \(v\) — diffusing on a tissue:

\[ \partial_t u = D_u \nabla^2 u + f(u,v), \qquad \partial_t v = D_v \nabla^2 v + g(u,v) \]

Naively, diffusion smooths inhomogeneities. Turing showed, astonishingly, that diffusion can destabilise a homogeneous steady state when \(D_v \gg D_u\): the inhibitor spreads faster than the activator, producing local “peaks of activation surrounded by halos of inhibition” that tile the tissue into spots or stripes.

1.1 Linear Stability Analysis

Let \((u_0, v_0)\) be a spatially uniform fixed point: \(f(u_0, v_0)=g(u_0, v_0)=0\). Perturb \(u = u_0 + \delta u\), \(v = v_0 + \delta v\) and linearise:

\[ \partial_t \begin{pmatrix}\delta u\\\delta v\end{pmatrix} = \left(\underbrace{\begin{pmatrix}f_u & f_v\\g_u & g_v\end{pmatrix}}_{J} + \begin{pmatrix}D_u\nabla^2 & 0\\0 & D_v\nabla^2\end{pmatrix}\right)\begin{pmatrix}\delta u\\\delta v\end{pmatrix} \]

Expand in Fourier modes \(\propto e^{\lambda t + i\mathbf{k}\cdot\mathbf{x}}\):

\[ \lambda = \tfrac{1}{2}(\mathrm{tr}\,M_k) \pm \tfrac{1}{2}\sqrt{(\mathrm{tr}\,M_k)^2 - 4\det M_k} \]

where \(M_k = J - k^2 \mathrm{diag}(D_u, D_v)\).

1.2 Instability Conditions

For Turing instability we need:

Stable without diffusion: \(\mathrm{tr}\,J < 0\) and \(\det J > 0\).
Unstable with diffusion for some \(k\neq 0\): \(\det M_k < 0\) for that \(k\).

The critical wavenumber emerges from \(d(\det M_k)/dk^2 = 0\):

\[ k_c^2 = \frac{D_v f_u + D_u g_v}{2 D_u D_v} \]

The instability condition becomes:

\[ (D_v f_u + D_u g_v)^2 > 4 D_u D_v \det J \]

This requires \(D_v/D_u\) to exceed a critical ratio — typically 10 to 100. The characteristic pattern wavelength is \(\lambda \sim 2\pi/k_c \sim \sqrt{D_v/|g_v|}\).

2. Biological Turing Patterns

2.1 Zebrafish Stripes

Zebrafish stripes (Danio rerio) arise from interactions between melanophores (black, \(u\)) and xanthophores (yellow, \(v\)). Melanophores inhibit nearby xanthophores; xanthophores at distance promote melanophores. This long-range activation / short-range inhibition matches Turing's model precisely (Nakamasu et al. 2009 used laser ablation to confirm the interaction topology in vivo).

2.2 Leopard Spots and Mammalian Coats

Murray (1981) showed that the typical coat patterns — spots on leopards, stripes on zebras, reticulation on giraffes — arise naturally from a single Turing system with mild variations in tissue geometry during embryogenesis. Large embryos (e.g. elephants) are outside the unstable regime — hence plain coats.

2.3 Seashell Patterns (Meinhardt)

Hans Meinhardt's The Algorithmic Beauty of Sea Shells (2003) models mollusc pigmentation as 1-D reaction-diffusion on the growing shell edge, with time as the second axis. This reproduces oblique stripes, triangles, and checkerboard patterns observed in real species.

2.4 Cat Fur Patterns

The spots of a cheetah and the stripes of a tiger emerge from the same mechanism. See the/feline-biophysicscourse, Module 6, for detailed modelling of cat-specific coat pigmentation.

3. Morphogen Gradients: The French Flag Model

Wolpert (1969) proposed an alternative to Turing's self-organisation: localised positional information. A morphogen source at one end of a tissue produces a spatial gradient; cells read their local [M] and adopt one of several fates based on thresholds. Wolpert named this the French flag model — three zones (blue, white, red).

3.1 Diffusion-Degradation Gradient

At steady state, morphogen \(c(x)\) satisfies:

\[ D \frac{d^2 c}{dx^2} - k c = 0 \]

with \(c(0) = c_0\) (source) and \(c(\infty) = 0\). Solution:

\[ c(x) = c_0\, e^{-x/\lambda}, \quad \lambda = \sqrt{D/k} \]

The decay length \(\lambda\) sets the scale of positional information.

3.2 Bicoid in Drosophila

The Bicoid (Bcd) protein gradient in the Drosophila embryo is the prototype morphogen gradient. Measurements (Driever & Nusslein-Volhard 1988; Gregor et al. 2007) show:

Embryo length: \(L \approx 500\,\mu\text{m}\).
Decay length: \(\lambda \approx 100\,\mu\text{m} \approx L/5\).
Diffusion coefficient: \(D \approx 0.3\,\mu\text{m}^2/\text{s}\).
Half-life: \(t_{1/2} \approx 30\,\text{min}\).

3.3 Scale Invariance Problem

A naive diffusion gradient has decay length \(\lambda = \sqrt{D/k}\) set by molecular parameters — it does not scale with embryo size \(L\). But real embryos vary in size within species; pattern boundaries must scale. Mechanisms that restore scaling:

Expansion-repression: a secondary molecule modulates\(k\) so \(\lambda \propto L\) (Ben-Zvi & Barkai 2010).
Pre-steady-state readout: cells sample the gradient before it equilibrates; time-dependent patterns scale.

4. Differential Adhesion & Mechanical Self-Organisation

Morphogenesis is not only chemical. Malcolm Steinberg (1963) observed that when embryonic cells of different types are dissociated and mixed, they spontaneously re-sort into concentric layers — like oil and water separating. He proposed the differential adhesion hypothesis: cells are like liquids with effective surface tensions set by their adhesion energies.

4.1 Steinberg's Inequality

Let \(J_{AA}, J_{BB}, J_{AB}\) be the bond energies between pairs of cell types. The interfacial energy is:

\[ \gamma_{AB} = J_{AB} - \tfrac{1}{2}(J_{AA} + J_{BB}) \]

If \(\gamma_{AB} > 0\), types segregate. If further \(J_{BB} > J_{AA}\), the \(B\) cells form a core and \(A\) cells form the outer layer — reproducing the germ-layer organisation of early embryos.

4.2 Cadherin Expression

Cadherins — homophilic adhesion molecules on cell surfaces — provide the molecular substrate: different cadherin types have different binding affinities, implementing the \(J\) matrix. Foty & Steinberg (2005) quantitatively confirmed that cadherin expression levels predict sorted configurations.

5. Scale-Invariant Patterns

For Turing patterns, the fundamental wavelength is \(\lambda \sim 2\pi/k_c\), set by diffusion constants. On a finite domain of length \(L\), the number of stripes is therefore \(N \sim L/\lambda\). Larger animals have more stripes:

Cheetah (small, short fur): many spots.
Jaguar (mid-size): rosettes (doublet spots).
Tiger (large): horizontal stripes.
Giraffe: reticulated patches.
Elephant (outside unstable regime): plain.

Murray (1988) showed that animal tail patterns are particularly diagnostic: the tail is 1-dimensional-like, and narrower than the body. Stripes on a broad body sometimes become spots on a narrow tail — consistent with a 1D-vs-2D instability geometry.

6. Turing Pattern Formation Stages

7. Zebrafish Stripes & Drosophila Segmentation

8. Simulation 1 — Gray-Scott Reaction-Diffusion

Pearson (1993) classified the Gray-Scott model's behaviour in parameter space. By tuning feed rate\(F\) and kill rate \(k\), one obtains spots, stripes, holes, spirals, or chaos — all from two interacting chemicals. Click “Run” to see six parameter regimes.

Python

script.py60 lines

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

# Gray-Scott reaction-diffusion system (Pearson 1993)
# dU/dt = D_u * Lap U - U V^2 + F (1 - U)
# dV/dt = D_v * Lap V + U V^2 - (F + k) V

def laplacian(Z):
    return (np.roll(Z,  1, axis=0) + np.roll(Z, -1, axis=0)
          + np.roll(Z,  1, axis=1) + np.roll(Z, -1, axis=1)
          - 4*Z)

def run_gs(F, k, Du=0.16, Dv=0.08, N=200, T=6000):
    U = np.ones((N, N))
    V = np.zeros((N, N))
    # Seed a small blob of V
    r = 10
    cx, cy = N//2, N//2
    U[cx-r:cx+r, cy-r:cy+r] = 0.5
    V[cx-r:cx+r, cy-r:cy+r] = 0.25
    U += 0.02*np.random.randn(N, N)
    V += 0.02*np.random.randn(N, N)
    U = np.clip(U, 0, 1); V = np.clip(V, 0, 1)
    dt = 1.0
    for t in range(T):
        Lu = laplacian(U); Lv = laplacian(V)
        UV2 = U*V*V
        U += dt*(Du*Lu - UV2 + F*(1-U))
        V += dt*(Dv*Lv + UV2 - (F+k)*V)
    return V

# Pearson's classic parameter settings
patterns = {
    'spots (kappa)':  (0.0545, 0.062),
    'stripes (mu)':   (0.046, 0.065),
    'spirals (eta)':  (0.03, 0.0625),
    'holes (alpha)':  (0.010, 0.047),
    'waves (theta)':  (0.03, 0.057),
    'chaotic (lambda)': (0.039, 0.065),
}

fig, axes = plt.subplots(2, 3, figsize=(16, 10))
fig.patch.set_facecolor('#030712')
fig.suptitle("Gray-Scott Reaction-Diffusion: Pearson's Parameter Space",
             fontsize=15, color='#fde68a', fontweight='bold')

for ax, (name, (F, k)) in zip(axes.flatten(), patterns.items()):
    V = run_gs(F, k, N=150, T=4000)
    ax.imshow(V, cmap='inferno', vmin=0, vmax=0.5)
    ax.set_title(f'{name}: F={F:.4f}, k={k:.4f}',
                 color='#fde68a', fontsize=10)
    ax.axis('off')

plt.tight_layout()
plt.savefig('output.png', dpi=140, bbox_inches='tight', facecolor='#030712')
print("Gray-Scott generates spots, stripes, spirals, holes -- all from 2 chemicals")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Simulation 2 — French Flag & Bicoid

The French flag model: an exponential morphogen gradient with two thresholds. We also compare Bicoid gradients at different embryo lengths (poor scale invariance) and illustrate how cooperative binding (Hill functions) sharpens cell-fate boundaries.

Python

script.py88 lines

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

# French flag model (Wolpert 1969): morphogen M diffuses from source, forms gradient
# Cells read local [M] and adopt fate A, B, or C based on two thresholds
# Bicoid gradient in Drosophila embryo: L ~ 500 micron

L = 500.0   # embryo length (microns)
D = 1.0     # diffusion constant (micron^2/s)
k = 1/300.0 # degradation (1/s)
J = 1.0     # flux (arbitrary)

# Steady state solution: c(x) = c0 * exp(-x / lambda), lambda = sqrt(D/k)
lam = np.sqrt(D/k)
print(f"Decay length lambda = {lam:.2f} microns")

x = np.linspace(0, L, 400)
c = np.exp(-x/lam)

# Two thresholds
theta_A = 0.5
theta_B = 0.1

# Also show the effect of noise + averaging
np.random.seed(42)
c_noisy = c * (1 + 0.1*np.random.randn(len(x)))

fig, axes = plt.subplots(1, 3, figsize=(18, 6))
fig.patch.set_facecolor('#030712')
fig.suptitle("The French Flag Model (Wolpert 1969) + Bicoid Gradient",
             fontsize=15, color='#fde68a', fontweight='bold')

# Panel 1: Gradient + thresholds
ax = axes[0]
ax.set_facecolor('#1a1200')
ax.plot(x, c, color='#fbbf24', linewidth=2.5, label='morphogen [M](x)')
ax.axhline(theta_A, color='#f87171', linestyle='--', linewidth=1.5, label=f'theta_A = {theta_A}')
ax.axhline(theta_B, color='#86efac', linestyle='--', linewidth=1.5, label=f'theta_B = {theta_B}')
# Fill regions
x_A = x[c >= theta_A]
x_B = x[(c < theta_A) & (c >= theta_B)]
x_C = x[c < theta_B]
ax.axvspan(0, x_A[-1], alpha=0.2, color='#22d3ee', label='fate A')
ax.axvspan(x_A[-1], x_B[-1] if len(x_B) else L, alpha=0.2, color='#ffffff', label='fate B')
ax.axvspan(x_B[-1] if len(x_B) else x_A[-1], L, alpha=0.2, color='#f87171', label='fate C')
ax.set_xlabel('position along embryo (um)', color='white')
ax.set_ylabel('[M] normalized', color='white')
ax.set_title('French flag', color='#fde68a')
ax.legend(facecolor='#1a1200', edgecolor='#fbbf24', labelcolor='white', fontsize=8)
ax.grid(True, alpha=0.2, color='#fbbf24'); ax.tick_params(colors='white')
for sp in ax.spines.values(): sp.set_color('#fbbf24')

# Panel 2: Bicoid at different embryo sizes (scale invariance?)
ax = axes[1]
ax.set_facecolor('#1a1200')
for L_i, col in zip([300, 500, 700, 900], ['#f87171', '#fbbf24', '#86efac', '#22d3ee']):
    xi = np.linspace(0, L_i, 200)
    ci = np.exp(-xi/lam)
    ax.plot(xi/L_i, ci, color=col, linewidth=2, label=f'L={L_i} um')
ax.set_xlabel('relative position x/L', color='white')
ax.set_ylabel('[Bicoid]', color='white')
ax.set_title('Bicoid gradient does NOT auto-scale with L', color='#fde68a')
ax.legend(facecolor='#1a1200', edgecolor='#fbbf24', labelcolor='white', fontsize=8)
ax.grid(True, alpha=0.2, color='#fbbf24'); ax.tick_params(colors='white')
for sp in ax.spines.values(): sp.set_color('#fbbf24')

# Panel 3: Hunchback target gene sharp response via cooperative binding (Hill function, n=5)
ax = axes[2]
ax.set_facecolor('#1a1200')
nuclei = np.linspace(0, 1, 100)
bcd = np.exp(-nuclei*L/lam)
for n in [1, 3, 5, 10]:
    K = 0.2
    hb = bcd**n / (bcd**n + K**n)
    ax.plot(nuclei, hb, linewidth=2, label=f'Hill n={n}')
ax.set_xlabel('x/L', color='white')
ax.set_ylabel('Hunchback activation', color='white')
ax.set_title('Cooperative binding sharpens boundary', color='#fde68a')
ax.legend(facecolor='#1a1200', edgecolor='#fbbf24', labelcolor='white', fontsize=8)
ax.grid(True, alpha=0.2, color='#fbbf24'); ax.tick_params(colors='white')
for sp in ax.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

10. Simulation 3 — Differential Adhesion Sorting

A Metropolis Monte Carlo implementation of Steinberg's differential adhesion hypothesis. Two cell types with distinct binding energies sort into concentric layers, recapitulating germ-layer organisation.

Python

script.py70 lines

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

# Differential adhesion Monte Carlo (Steinberg 1963)
# Two cell types A and B with adhesion energies J_AA, J_BB, J_AB
# Stable sorted state when J_AB > (J_AA + J_BB)/2 (interfacial energy minimised)

np.random.seed(1)
N = 50
grid = np.random.choice([1, 2], size=(N, N))  # 1=A, 2=B

# Energy matrix (higher = more favorable bond)
J = {(1,1): 3.0, (2,2): 4.0, (1,2): 1.0, (2,1): 1.0,
     (0,0): 0.0, (0,1): 0.0, (0,2): 0.0, (1,0): 0.0, (2,0): 0.0}

def neighbor_sum(g, i, j):
    N = g.shape[0]
    s = 0
    for di, dj in [(1,0), (-1,0), (0,1), (0,-1)]:
        ni, nj = (i+di) % N, (j+dj) % N
        s += J[(g[i,j], g[ni,nj])]
    return s

def step(g, T=0.3):
    N = g.shape[0]
    for _ in range(N*N):
        i, j = np.random.randint(0, N, 2)
        i2, j2 = np.random.randint(0, N, 2)
        if g[i,j] == g[i2,j2]: continue
        E_before = neighbor_sum(g, i, j) + neighbor_sum(g, i2, j2)
        g[i,j], g[i2,j2] = g[i2,j2], g[i,j]
        E_after = neighbor_sum(g, i, j) + neighbor_sum(g, i2, j2)
        dE = -(E_after - E_before)  # we maximised bonds, so energy lowers if E_after higher
        if dE > 0 and np.random.rand() > np.exp(-dE/T):
            g[i,j], g[i2,j2] = g[i2,j2], g[i,j]
    return g

# Record snapshots
snapshots = [grid.copy()]
g = grid.copy()
steps_per_snap = [0, 10, 50, 200, 1000, 3000]
for target in steps_per_snap[1:]:
    while snapshots[-1].any():  # dummy loop with inner counter
        break
    # Actually run
    total = target - len(snapshots)*200  # coarse
    for _ in range(max(target - sum(steps_per_snap[:len(snapshots)]), 100)):
        step(g, T=0.5)
    snapshots.append(g.copy())
    if len(snapshots) >= 6: break

fig, axes = plt.subplots(2, 3, figsize=(15, 10))
fig.patch.set_facecolor('#030712')
fig.suptitle("Differential Adhesion Hypothesis (Steinberg 1963)",
             fontsize=15, color='#fde68a', fontweight='bold')

from matplotlib.colors import ListedColormap
cmap = ListedColormap(['#1a1200', '#fbbf24', '#22d3ee'])
for ax, snap, t in zip(axes.flatten(), snapshots, [0, 200, 400, 800, 2000, 5000]):
    ax.imshow(snap, cmap=cmap, vmin=0, vmax=2, interpolation='nearest')
    ax.set_title(f'step {t}', color='#fde68a')
    ax.axis('off')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#030712')
print("Cells sort spontaneously: B (cyan) wraps around A (amber) (J_BB > J_AA)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Module Summary

Turing instability

Diffusion can destabilize a uniform state when D_v/D_u is large.

Instability condition

(D_v f_u + D_u g_v)^2 > 4 D_u D_v det J.

Pattern wavelength

lambda ~ 2pi/k_c; wavelength set by diffusion constants.

Zebrafish stripes

Nakamasu 2009: short-range activation, long-range inhibition confirmed.

Wolpert French flag

Positional information via morphogen gradient + thresholds.

Bicoid in Drosophila

lambda approx L/5; gradient + cooperative binding -> sharp boundaries.

Steinberg adhesion

gamma_AB = J_AB - (J_AA+J_BB)/2; sorting -> germ layers.

Scale-invariant patterns

Larger embryos -> more stripes; elephants fall out of unstable regime.

References

Turing, A. M. (1952). The chemical basis of morphogenesis. Phil. Trans. R. Soc. B, 237, 37–72.
Wolpert, L. (1969). Positional information and the spatial pattern of cellular differentiation. J. Theor. Biol., 25, 1–47.
Meinhardt, H. (2003). The Algorithmic Beauty of Sea Shells. Springer.
Murray, J. D. (1981). A pre-pattern formation mechanism for animal coat markings. J. Theor. Biol., 88, 161–199.
Nakamasu, A., Takahashi, G., Kanbe, A., & Kondo, S. (2009). Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. PNAS, 106, 8429–8434.
Driever, W. & Nusslein-Volhard, C. (1988). A gradient of bicoid protein in Drosophila embryos. Cell, 54, 83–93.
Gregor, T. et al. (2007). Probing the limits to positional information. Cell, 130, 153–164.
Ben-Zvi, D. & Barkai, N. (2010). Scaling of morphogen gradients by an expansion-repression integral feedback control. PNAS, 107, 6924–6929.
Steinberg, M. S. (1963). Reconstruction of tissues by dissociated cells. Science, 141, 401–408.
Foty, R. A. & Steinberg, M. S. (2005). The differential adhesion hypothesis: a direct evaluation. Dev. Biol., 278, 255–263.
Pearson, J. E. (1993). Complex patterns in a simple system. Science, 261, 189–192.
Kondo, S. & Miura, T. (2010). Reaction-diffusion model as a framework for understanding biological pattern formation. Science, 329, 1616–1620.

← Module 5: Minimal Surfaces Next: Module 7 — Packing & Tessellation →