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

1. Waddington, C.H. (1957). The Strategy of the Genes.
2. Nobis, G. et al. (2025). "Fractional Diffusion Bridge Models." NeurIPS.
3. Huang, S. et al. (2005). "Cell fates as high-dimensional attractor states."
4. Ferrell, J.E. (2012). "Bistability, bifurcations, and Waddington's landscape."
5. 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.

SDF Interleaving

This skill connects to Software Design for Flexibility (Hanson & Sussman, 2021):

Primary Chapter: 8. Degeneracy

Concepts: redundancy, fallback, multiple strategies, robustness

GF(3) Balanced Triad

waddington-landscape (+) + SDF.Ch8 (−) + [balancer] (○) = 0

Skill Trit: 1 (PLUS - generation)

Secondary Chapters

• Ch7: Propagators
• Ch1: Flexibility through Abstraction
• Ch6: Layering
• Ch4: Pattern Matching

Connection Pattern

Degeneracy provides fallbacks. This skill offers redundant strategies.