waddington-landscape

# 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

```mathematica
(* 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)

```mathematica
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)

```mathematica
(* 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

```python
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

```python
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.