Applications: Diffusion in Biology

Learning Objectives

By the end of this section, you will be able to:

1. Understand Fick's laws of diffusion and their connection to the heat equation
2. Explain how morphogen gradients create spatial patterns in developing embryos
3. Derive the characteristic length scale for diffusion with decay
4. Analyze Turing's reaction-diffusion theory of biological pattern formation
5. Apply diffusion models to drug delivery and pharmacokinetics
6. Calculate oxygen diffusion profiles using the Krogh cylinder model
7. Recognize diffusion processes in neural signaling and synaptic transmission

The Big Picture: Life Runs on Diffusion

"Biology is the study of complicated things that give the appearance of having been designed for a purpose." — Richard Dawkins

Behind this apparent design lies a simple physical process: diffusion. From the spreading of oxygen in your lungs to the formation of stripes on a zebra, the heat equation (or its biological equivalent, Fick's law) governs life at the molecular scale.

Fick's Laws of Diffusion

In 1855, Adolf Fick, a German physiologist, published his laws of diffusion while studying how gases move through membranes. These laws are the biological equivalent of Fourier's law for heat.

Fick's First Law: Flux is Proportional to Gradient

This states that molecules flow from regions of high concentration to low concentration, at a rate proportional to how steep the concentration gradient is. The negative sign ensures flow is "downhill" in concentration.

Fick's Second Law: The Diffusion Equation

Combining Fick's First Law with conservation of mass yields Fick's Second Law:

Morphogen Gradients: How Embryos Know Left from Right

One of the great mysteries of biology is how a single fertilized egg develops into a complex organism with distinct organs in the right places. The answer lies in morphogen gradients — concentration fields of signaling molecules that provide positional information to cells.

The Diffusion-Decay Model

A morphogen is produced at a source (say, one end of the embryo), diffuses through tissue, and is degraded along the way. The steady-state concentration satisfies:

where $k$ is the decay rate. This has the beautiful exponential solution:

The French Flag Model

Lewis Wolpert's famous "French Flag" model explains how a smooth gradient creates sharp boundaries. Cells respond to the local concentration:

• High concentration (C > T₁): Activate Gene A (blue zone)
• Medium concentration (T₂ < C < T₁): Activate Gene B (white zone)
• Low concentration (C < T₂): Activate Gene C (red zone)

1import numpy as np
2import matplotlib.pyplot as plt
3from scipy.integrate import odeint
4
5def morphogen_gradient(y, t, D, k, source_rate):
6    """
7    1D morphogen gradient with diffusion and decay.
8
9    The steady-state solution is:
10    C(x) = C_0 * exp(-x/λ)  where λ = sqrt(D/k)
11
12    This exponential decay creates the "French Flag" model
13    of developmental biology.
14    """
15    n = len(y)
16    dx = 1.0 / (n - 1)
17    dydt = np.zeros(n)
18
19    # Diffusion term: D * d²C/dx²
20    for i in range(1, n-1):
21        dydt[i] = D * (y[i+1] - 2*y[i] + y[i-1]) / dx**2
22
23    # Decay term: -k * C
24    dydt -= k * y
25
26    # Source at x=0
27    dydt[0] = source_rate - k * y[0]
28
29    return dydt
30
31# Parameters
32D = 0.1       # Diffusion coefficient
33k = 0.05      # Decay rate
34lambda_char = np.sqrt(D/k)  # Characteristic length
35
36# Spatial grid
37n_points = 100
38x = np.linspace(0, 2, n_points)
39
40# Initial condition: zero everywhere
41y0 = np.zeros(n_points)
42
43# Time points
44t = np.linspace(0, 50, 200)
45
46# Solve ODE system
47source_rate = 1.0
48solution = odeint(morphogen_gradient, y0, t, args=(D, k, source_rate))
49
50# Plot evolution
51plt.figure(figsize=(12, 5))
52
53plt.subplot(1, 2, 1)
54times = [0, 5, 10, 20, 50]
55colors = plt.cm.viridis(np.linspace(0, 1, len(times)))
56for i, time_idx in enumerate([0, 25, 50, 100, -1]):
57    plt.plot(x, solution[time_idx], color=colors[i],
58             linewidth=2, label=f't = {t[time_idx]:.0f}')
59
60# Analytical steady state
61C_steady = np.exp(-x / lambda_char)
62plt.plot(x, C_steady, 'k--', linewidth=2, label='Steady state (analytical)')
63
64plt.xlabel('Position (x)')
65plt.ylabel('Morphogen Concentration')
66plt.title('Morphogen Gradient Formation')
67plt.legend()
68plt.grid(True, alpha=0.3)
69
70# Show threshold regions (French Flag model)
71plt.subplot(1, 2, 2)
72thresholds = [0.6, 0.3, 0.1]
73colors = ['#3b82f6', '#ffffff', '#ef4444']  # Blue, White, Red
74labels = ['High (Gene A)', 'Medium (Gene B)', 'Low (Gene C)']
75
76for i in range(len(thresholds)):
77    lower = thresholds[i] if i < len(thresholds) else 0
78    upper = 1.0 if i == 0 else thresholds[i-1]
79    mask = (C_steady >= lower) & (C_steady < upper)
80    plt.fill_between(x, 0, 1, where=mask, color=colors[i],
81                     alpha=0.7, label=labels[i])
82
83plt.plot(x, C_steady, 'k-', linewidth=2, label='Concentration')
84plt.xlabel('Position (x)')
85plt.ylabel('Concentration / Threshold')
86plt.title('"French Flag" Model: Gene Activation Zones')
87plt.legend(loc='upper right')
88plt.ylim(0, 1.1)
89plt.grid(True, alpha=0.3)
90
91plt.tight_layout()
92plt.show()
93
94print(f"Characteristic length λ = √(D/k) = {lambda_char:.3f}")
95print(f"At x = λ, concentration drops to 1/e ≈ 37% of source")

Turing Patterns: Spots, Stripes, and the Mathematics of Life

In 1952, Alan Turing — better known for his work on computation — published a revolutionary paper titled "The Chemical Basis of Morphogenesis." He proposed that biological patterns emerge from reaction-diffusion systems.

Turing's Key Insight

Consider two chemicals: an activator (A) that promotes its own production, and an inhibitor (B) that suppresses the activator. Turing showed that if the inhibitor diffuses faster than the activator, spatial patterns spontaneously emerge!

Why Patterns Form

1. Local activation: A small fluctuation in activator concentration triggers positive feedback — more A makes more A.
2. Lateral inhibition: The activator also produces inhibitor B, which diffuses outward and suppresses A in neighboring regions.
3. Scale separation: Because B diffuses faster, it creates a "halo" of inhibition around each activation peak.
4. Pattern selection: The balance of activation and inhibition selects a characteristic wavelength — this determines whether you get spots, stripes, or labyrinths.

1import numpy as np
2import matplotlib.pyplot as plt
3from matplotlib.animation import FuncAnimation
4
5def gray_scott_step(u, v, Du, Dv, f, k, dx, dt):
6    """
7    Gray-Scott reaction-diffusion model step.
8
9    Reactions:
10    U + 2V → 3V  (autocatalytic conversion)
11    V → P        (decay of activator)
12    Feed → U     (replenishment of substrate)
13
14    Equations:
15    ∂u/∂t = Du∇²u - uv² + f(1-u)
16    ∂v/∂t = Dv∇²v + uv² - (f+k)v
17    """
18    # Laplacian using periodic boundary conditions
19    laplacian_u = (
20        np.roll(u, 1, axis=0) + np.roll(u, -1, axis=0) +
21        np.roll(u, 1, axis=1) + np.roll(u, -1, axis=1) - 4*u
22    ) / dx**2
23
24    laplacian_v = (
25        np.roll(v, 1, axis=0) + np.roll(v, -1, axis=0) +
26        np.roll(v, 1, axis=1) + np.roll(v, -1, axis=1) - 4*v
27    ) / dx**2
28
29    # Reaction terms
30    uvv = u * v**2
31
32    # Update equations
33    du = Du * laplacian_u - uvv + f * (1 - u)
34    dv = Dv * laplacian_v + uvv - (f + k) * v
35
36    return u + dt * du, v + dt * dv
37
38# Parameters for different patterns
39patterns = {
40    'spots': {'f': 0.035, 'k': 0.065},
41    'stripes': {'f': 0.055, 'k': 0.062},
42    'labyrinth': {'f': 0.029, 'k': 0.057}
43}
44
45# Choose pattern
46pattern = 'spots'
47params = patterns[pattern]
48
49# Grid setup
50n = 128
51dx = 1.0
52dt = 1.0
53
54# Diffusion coefficients (key: Dv < Du for pattern formation!)
55Du = 0.16
56Dv = 0.08  # Inhibitor diffuses SLOWER than activator
57
58# Initial conditions: uniform u, seed v
59u = np.ones((n, n))
60v = np.zeros((n, n))
61
62# Random initial seeds
63for _ in range(20):
64    cx, cy = np.random.randint(0, n, 2)
65    radius = np.random.randint(3, 8)
66    y, x = np.ogrid[-cx:n-cx, -cy:n-cy]
67    mask = x**2 + y**2 <= radius**2
68    u[mask] = 0.5
69    v[mask] = 0.25
70
71# Run simulation
72steps = 5000
73for step in range(steps):
74    u, v = gray_scott_step(u, v, Du, Dv, params['f'], params['k'], dx, dt)
75
76# Plot result
77plt.figure(figsize=(10, 5))
78
79plt.subplot(1, 2, 1)
80plt.imshow(v, cmap='RdBu_r')
81plt.colorbar(label='Activator (v)')
82plt.title(f'Turing Pattern: {pattern.capitalize()}')
83
84plt.subplot(1, 2, 2)
85plt.imshow(u, cmap='RdBu')
86plt.colorbar(label='Substrate (u)')
87plt.title('Substrate Distribution')
88
89plt.tight_layout()
90plt.show()
91
92print(f"Pattern type: {pattern}")
93print(f"Diffusion ratio Du/Dv = {Du/Dv:.2f}")
94print(f"Key insight: Pattern forms because Dv < Du")

Neural Signaling: Diffusion at the Speed of Thought

Every thought, movement, and sensation depends on electrical signals jumping across tiny gaps called synapses. The gap is only about 20 nm wide — but why so small?

The Synaptic Cleft Problem

When an action potential arrives at a synapse, vesicles release neurotransmitters (like glutamate or dopamine) into the synaptic cleft. These molecules must:

1. Diffuse across the gap to the post-synaptic membrane
2. Bind to receptors before being cleared by reuptake
3. Do this in microseconds to enable rapid signaling

Why Timing Matters

The characteristic diffusion time scales as $t \sim L^2/D$. For the synaptic cleft:

• Gap width: L ≈ 20 nm = 2 × 10⁻⁸ m
• Diffusion coefficient: D ≈ 10⁻⁵ cm²/s
• Diffusion time: t ≈ L²/D ≈ 4 μs

Drug Delivery: Designing Medicines That Reach Their Target

Pharmaceutical effectiveness depends critically on pharmacokinetics— how drugs move through the body. Diffusion governs this process, and the heat equation helps design better drug delivery systems.

The Therapeutic Window

Every drug has a "therapeutic window" — concentrations between the minimum effective concentration (MEC) and maximum tolerated concentration (MTC). The goal is to keep drug levels in this window as long as possible.

Delivery Strategies

Oxygen Transport: The Krogh Cylinder Model

Every cell in your body needs oxygen, delivered by a network of capillaries. But oxygen can only diffuse so far before being consumed — this sets fundamental limits on tissue architecture.

Krogh's Model (1919)

August Krogh modeled tissue as a cylinder of radius $R$centered on a capillary of radius $R_0$. Oxygen diffuses radially outward while being consumed by metabolism. The steady-state equation is:

where $M$ is the metabolic rate (oxygen consumption per unit volume).

The Critical Radius

The critical radius is where oxygen concentration drops to zero — cells beyond this become hypoxic. It depends on:

• Capillary PO₂: Higher arterial oxygen extends the reach
• Metabolic rate: More active tissues consume oxygen faster, shrinking the critical radius
• Diffusion coefficient: Varies with tissue type and temperature

Python Implementation

The Python examples above demonstrate how to simulate biological diffusion systems. Key points:

1. Morphogen gradients reach steady state when diffusion balances decay; the characteristic length λ = √(D/k) sets the scale.
2. Turing patterns require specific parameter ratios (activator/inhibitor diffusion) and emerge from tiny initial perturbations.
3. Numerical methods use the same finite difference schemes we developed for the heat equation — the Laplacian is approximated identically.

Test Your Understanding

Summary

The heat equation, disguised as Fick's law, is the mathematical backbone of countless biological processes. From embryonic development to neural signaling, from drug delivery to tissue oxygenation, diffusion governs life at the molecular scale.

Key Equations

Key Takeaways

1. Fick's laws are the heat equation with concentration instead of temperature; all our mathematical tools transfer directly.
2. Morphogen gradients provide positional information in embryos; cells read concentration to know their location.
3. Turing patterns arise when an inhibitor diffuses faster than an activator, creating spots, stripes, and spirals.
4. Diffusion time scales as L²/D, which is why the synaptic cleft is only 20 nm wide and why capillaries must be densely spaced.
5. Drug delivery design uses diffusion models to maintain therapeutic concentrations over time and space.
6. The Krogh cylinder model predicts oxygen profiles around capillaries and explains why tissues have a maximum thickness.
7. Numerical methods from heat equation analysis apply directly to biological simulations.

Looking Ahead: In the next chapter, we'll explore the wave equation — another fundamental PDE that describes vibrations, sound, and electromagnetic waves. You'll see how hyperbolic PDEs differ from the parabolic heat equation in their behavior and solutions.