Waddington Landscape Skill
"The cell is like a ball rolling down a landscape of valleys. Once it enters a valley, it is canalized toward a particular fate." — Conrad Hal Waddington (1957)
Overview
Waddington's epigenetic landscape is the foundational metaphor for developmental biology:
| Concept | Landscape Metaphor | Mathematical Structure | |---------|-------------------|----------------------| | Cell | Ball/marble | State point θ(t) | | Differentiation | Rolling downhill | Gradient descent | | Cell fate | Valley bottom | Attractor basin | | Fate decision | Ridge/bifurcation | Critical point | | Landscape shape | Epigenetic regulation | Potential V(θ) |
The Mathematics
Gradient Flow on Potential Landscape
dθ/dt = -∇V(θ) + √(2T) dW(t)
Where:
θ = cell state (gene expression profile)
V(θ) = epigenetic potential (landscape height)
T = temperature (stochastic fluctuations)
dW = Brownian motion (noise)
This is Langevin dynamics on the Waddington landscape!
Connection to FDBM (NeurIPS 2025)
The Fractional Diffusion Bridge Models paper (Nobis et al., 2025) provides the modern framework:
Standard Brownian (H=0.5):
Memoryless diffusion
No long-range correlations
Fractional Brownian (H≠0.5):
H > 0.5: Superdiffusive, smooth paths, PERSISTENT
H < 0.5: Subdiffusive, rough paths, ANTI-PERSISTENT
For cell differentiation: H > 0.5 (cells remember their history)
The Hurst index H encodes the epigenetic memory of the cell!
Schrödinger Bridge Formulation
Transport cells from distribution Π₀ (pluripotent) to Π₁ (differentiated):
P^SB = argmin { D_KL(P || Q) ; P₀ = Π₀, P₁ = Π₁ }
P∈Paths
Where Q is the reference process (fBM with Hurst index H).
This is OPTIMAL TRANSPORT on the Waddington landscape!
Modelica Implementation
Basic Landscape with Pitchfork Bifurcation
(* Waddington landscape: pluripotent → differentiated *)
CreateSystemModel["Waddington.PitchforkFate", {
(* Pitchfork bifurcation: V(x,r) = rx²/2 + x⁴/4 *)
(* r > 0: single valley (pluripotent) *)
(* r < 0: two valleys (differentiated fates) *)
x'[t] == -(r[t] * x[t] + x[t]^3) - gamma * x'[t],
r'[t] == -beta, (* Development drives bifurcation *)
V[t] == r[t] * x[t]^2 / 2 + x[t]^4 / 4
}, t, <|
"ParameterValues" -> {beta -> 0.1, gamma -> 0.2},
"InitialValues" -> {x -> 0.001, r -> 1.0}
|>]
(* Simulate cell fate determination *)
sim = SystemModelSimulate["Waddington.PitchforkFate", 30];
SystemModelPlot[sim, {"x", "r", "V"}]
Three-Fate Landscape (Stem Cell → Neuron/Muscle/Blood)
CreateSystemModel["Waddington.ThreeFate", {
(* Potential with 3 minima at 120° angles *)
(* V = -a(x³ - 3xy²) + b(x² + y²)² + c(x² + y²) *)
x'[t] == -(3 a[t] * (x[t]^2 - y[t]^2) +
4 b * x[t] * (x[t]^2 + y[t]^2) + 2 c[t] * x[t]),
y'[t] == -(-6 a[t] * x[t] * y[t] +
4 b * y[t] * (x[t]^2 + y[t]^2) + 2 c[t] * y[t]),
(* Development changes landscape *)
a'[t] == rateA, (* Increasing asymmetry *)
c'[t] == -rateC, (* Decreasing central stability *)
diff[t] == Sqrt[x[t]^2 + y[t]^2] (* Differentiation progress *)
}, t, <|
"ParameterValues" -> {b -> 1.0, rateA -> 0.05, rateC -> 0.1},
"InitialValues" -> {x -> 0.1, y -> 0.05, a -> 0.0, c -> 1.0}
|>]
Fractional Dynamics (MA-fBM)
(* Markov Approximation of Fractional Brownian Motion *)
(* Following Daems et al. (2026) / Harms & Stefanovits *)
CreateSystemModel["Waddington.FractionalFate", {
(* Augmented state: Z = (X, Y₁, ..., Yₖ) *)
(* X = cell state, Yₖ = OU processes for memory *)
(* Main dynamics *)
x'[t] == Sum[omega[k] * y[k][t], {k, 1, K}] - gamma * x[t],
(* K Ornstein-Uhlenbeck processes (memory) *)
y[1]'[t] == -gamma1 * y[1][t], (* + dB *)
y[2]'[t] == -gamma2 * y[2][t],
y[3]'[t] == -gamma3 * y[3][t],
y[4]'[t] == -gamma4 * y[4][t],
y[5]'[t] == -gamma5 * y[5][t],
(* Hurst index H encoded in omega weights *)
(* H > 0.5: superdiffusive (cell memory) *)
(* H < 0.5: subdiffusive (exploratory) *)
}, t, <|
"ParameterValues" -> {
K -> 5,
gamma -> 0.1,
(* Geometric spacing: γₖ = r^(k-n) *)
gamma1 -> 0.25, gamma2 -> 0.5, gamma3 -> 1.0,
gamma4 -> 2.0, gamma5 -> 4.0,
(* L² optimal weights for H = 0.7 (superdiffusive) *)
omega1 -> 0.15, omega2 -> 0.25, omega3 -> 0.30,
omega4 -> 0.20, omega5 -> 0.10
}
|>]
Python Implementation
Waddington Landscape Class
import numpy as np
from dataclasses import dataclass
from typing import Callable, Tuple
@dataclass
class WaddingtonLandscape:
"""
Waddington's epigenetic landscape as potential function.
Connects to:
- Langevin dynamics (gradient descent + noise)
- FDBM (fractional diffusion bridges)
- Schrödinger bridges (optimal transport)
"""
n_fates: int = 2
bifurcation_param: float = 1.0
temperature: float = 0.01
hurst_index: float = 0.5 # H > 0.5 for cell memory
def potential(self, x: np.ndarray, t: float) -> float:
"""
Time-dependent potential V(x, t).
As t increases (development), bifurcations emerge.
"""
r = self.bifurcation_param * (1 - 2*t) # r: + → -
if self.n_fates == 2:
# Pitchfork: V = rx²/2 + x⁴/4
return r * x[0]**2 / 2 + x[0]**4 / 4
else:
# Three-fate: V = -a(x³-3xy²) + b(x²+y²)² + c(x²+y²)
x0, x1 = x[0], x[1]
a = t * 0.5 # Increasing asymmetry
b = 1.0
c = 1.0 - t # Decreasing central stability
return (-a * (x0**3 - 3*x0*x1**2) +
b * (x0**2 + x1**2)**2 +
c * (x0**2 + x1**2))
def gradient(self, x: np.ndarray, t: float) -> np.ndarray:
"""∇V(x, t) - drives cell toward fate."""
eps = 1e-6
grad = np.zeros_like(x)
for i in range(len(x)):
x_plus = x.copy(); x_plus[i] += eps
x_minus = x.copy(); x_minus[i] -= eps
grad[i] = (self.potential(x_plus, t) -
self.potential(x_minus, t)) / (2 * eps)
return grad
def simulate_langevin(
self,
x0: np.ndarray,
n_steps: int = 1000,
dt: float = 0.01
) -> np.ndarray:
"""
Simulate cell trajectory via Langevin dynamics.
dx = -∇V dt + √(2T) dW
"""
trajectory = [x0.copy()]
x = x0.copy()
for step in range(n_steps):
t = step * dt
grad = self.gradient(x, t)
noise = np.random.randn(*x.shape)
x = x - grad * dt + np.sqrt(2 * self.temperature * dt) * noise
trajectory.append(x.copy())
return np.array(trajectory)
def find_cell_fates(landscape: WaddingtonLandscape, t: float = 1.0):
"""Find attractor basins (cell fates) at time t."""
from scipy.optimize import minimize
fates = []
# Start from multiple initial conditions
for x0 in [[-1, 0], [1, 0], [0, 1], [0, -1], [0, 0]]:
result = minimize(
lambda x: landscape.potential(np.array(x), t),
x0, method='L-BFGS-B'
)
if result.success:
fates.append(result.x)
# Remove duplicates
unique_fates = []
for fate in fates:
is_new = all(np.linalg.norm(fate - f) > 0.1 for f in unique_fates)
if is_new:
unique_fates.append(fate)
return unique_fates
FDBM Integration
class FractionalWaddington(WaddingtonLandscape):
"""
Waddington landscape with Fractional Brownian Motion.
Based on: Nobis et al. "Fractional Diffusion Bridge Models" (NeurIPS 2025)
"""
def __init__(self, hurst_index: float = 0.7, n_ou: int = 5, **kwargs):
super().__init__(**kwargs)
self.hurst_index = hurst_index
self.n_ou = n_ou # Number of OU processes for MA-fBM
# Geometric spacing of mean-reversion rates
r = 2.0 # Spacing ratio
n = (n_ou + 1) / 2
self.gammas = [r**(k - n) for k in range(1, n_ou + 1)]
# L² optimal weights (approximate for given H)
self.omegas = self._compute_optimal_weights()
def _compute_optimal_weights(self) -> np.ndarray:
"""
Compute L² optimal approximation coefficients.
See Daems et al. Proposition 3.
"""
H = self.hurst_index
K = self.n_ou
# Simplified: weights favor different scales based on H
# H > 0.5: emphasize slow processes (long memory)
# H < 0.5: emphasize fast processes (rough paths)
weights = np.zeros(K)
for k in range(K):
# Higher H → more weight on slower (smaller γ) processes
weights[k] = self.gammas[k] ** (H - 0.5)
return weights / np.sum(weights)
def simulate_fdbm(
self,
x0: np.ndarray,
x1: np.ndarray, # Target fate
n_steps: int = 1000,
dt: float = 0.01
) -> np.ndarray:
"""
Simulate fractional diffusion bridge from x0 to x1.
This is the Schrödinger bridge on Waddington landscape!
"""
d = len(x0)
K = self.n_ou
# Augmented state: Z = (X, Y₁, ..., Yₖ)
Z = np.zeros(d * (K + 1))
Z[:d] = x0 # X = cell state
# Y_k initialized to 0
trajectory = [x0.copy()]
for step in range(n_steps):
t = step * dt
s = 1.0 - t # Time to terminal
X = Z[:d]
Y = Z[d:].reshape(K, d)
# Drift toward x1 (bridge condition)
mu_1_t = X + sum(self.omegas[k] * Y[k] *
(np.exp(self.gammas[k] * s) - 1) / self.gammas[k]
for k in range(K))
sigma_sq = s # Approximate
# Update drift
bridge_drift = (x1 - mu_1_t) / (sigma_sq + 1e-6)
# Landscape gradient
grad = self.gradient(X, t)
# Combined dynamics
dB = np.random.randn(d) * np.sqrt(dt)
# Update X
dX = -grad * dt + bridge_drift * dt + np.sqrt(2 * self.temperature) * dB
# Update Y (OU processes driven by same dB)
dY = np.zeros((K, d))
for k in range(K):
dY[k] = -self.gammas[k] * Y[k] * dt + dB
Z[:d] += dX
Z[d:] += dY.flatten()
trajectory.append(Z[:d].copy())
return np.array(trajectory)
Connection to Other Skills
The KOH-Opalescence-Waddington Chain
┌─────────────────────────────────────────────────────────────────────┐
│ CRITICAL PHENOMENA ACROSS SCALES │
├─────────────────────────────────────────────────────────────────────┤
│ │
│ KOLMOGOROV-ONSAGER-HURST CRITICAL OPALESCENCE WADDINGTON │
│ (Turbulence/Scaling) (Phase Transitions) (Development) │
│ ───────────────────── ────────────────── ──────────── │
│ H = 1/3 ξ → ∞ at T_c Cell fate │
│ E(k) ~ k^(-5/3) Scale-free fluct. Bifurcations │
│ Anomalous dissipation Light scattering Canalization │
│ │
│ ────────────────── UNIFYING PRINCIPLE ────────────────────────── │
│ │
│ CRITICAL SLOWING DOWN: Near bifurcation/transition, │
│ correlation length diverges, fluctuations at all scales, │
│ small perturbations determine fate (sensitivity to initial cond) │
│ │
│ HURST INDEX H: Memory parameter controlling path roughness │
│ and long-range dependence in ALL these systems │
│ │
└─────────────────────────────────────────────────────────────────────┘
GF(3) Triads
Trit: +1 (PLUS/Generator)
Waddington landscape GENERATES cell fates through gradient flow.
It creates new states from potential function dynamics.
GF(3) Triads:
kolmogorov-onsager-hurst (-1) ⊗ critical-opalescence (0) ⊗ waddington-landscape (+1) = 0 ✓
langevin-dynamics (-1) ⊗ fokker-planck-analyzer (0) ⊗ waddington-landscape (+1) = 0 ✓
lyapunov-function (-1) ⊗ bifurcation (0) ⊗ waddington-landscape (+1) = 0 ✓
Biological Connections
Chreods (Waddington's Term)
Chreod = "necessary path" (Greek χρή + ὁδός)
A chreod is a canalized developmental pathway:
- Genetic perturbations don't derail development
- The valley walls RESIST deviation
- Robustness emerges from landscape topology
Reprogramming (iPSC)
Yamanaka factors (Oct4, Sox2, Klf4, c-Myc):
- Reverse the landscape!
- Push cells UPHILL from differentiated → pluripotent
- Requires energy input (against gradient)
In landscape terms:
Normal: dx/dt = -∇V (downhill to fate)
Reprogram: dx/dt = +∇V + driving (uphill forcing)
Critical Transitions in Development
At bifurcation points:
- Correlation length ξ → ∞ (critical opalescence!)
- Variance increases (early warning signal)
- Recovery time slows (critical slowing down)
- Small noise → fate determination
References
- Waddington, C.H. (1957). The Strategy of the Genes.
- Nobis, G. et al. (2025). "Fractional Diffusion Bridge Models." NeurIPS.
- Huang, S. et al. (2005). "Cell fates as high-dimensional attractor states."
- Ferrell, J.E. (2012). "Bistability, bifurcations, and Waddington's landscape."
- Mojtahedi, M. et al. (2016). "Cell fate decision as high-dimensional critical state transition."
Invocation
/waddington-landscape
Model developmental processes as gradient flows on epigenetic potential landscapes, with connections to Schrödinger bridges and fractional diffusion.