---
name: grn-dynamics
description: Gene Regulatory Network dynamics simulator — Boolean attractor search, ODE dynamics, GRN inference, Waddington landscape, and perturbation analysis. Use when asked to simulate GRN dynamics, find attractors, analyze cell fate decisions, or knock out/overexpress genes in silico.
tags: [gene-regulatory-network, boolean-networks, ode-dynamics, attractor-analysis, systems-biology, waddington-landscape, grn-inference, cell-fate]
---

# GRNDynamics: Gene Regulatory Network Dynamics Simulator

## Trigger Conditions
Use when the user wants to:
- Simulate gene regulatory network (GRN) dynamics from interaction data
- Find attractors (stable states) of a regulatory network: fixed points and limit cycles
- Model cell fate decisions, differentiation, and bistability
- Reconstruct a GRN from gene expression data using mutual information or correlation
- Analyze network topology: feedback loops, motifs, master regulators
- Simulate what happens when you knock out or overexpress a gene
- Compute Waddington epigenetic landscape for a GRN

Example triggers:
- "Simulate the dynamics of this gene regulatory network"
- "Find all attractors in this Boolean network"
- "Which genes are master regulators in this network?"
- "Model the EMT transition using ODE dynamics"
- "Infer a GRN from this expression matrix and simulate it"
- "What happens if I knock out NANOG in this pluripotency network?"

---

## Overview

**GRNDynamics** is a complete gene regulatory network simulation engine built from scratch in pure NumPy/SciPy. It combines three complementary modeling frameworks under one pipeline:

1. **Boolean Network Dynamics** (discrete, exact) — exhaustive attractor enumeration for N ≤ 22
2. **ODE Continuous Dynamics** (continuous, mechanistic) — Hill-function-based with adaptive RK4
3. **Network Inference** from expression data — ARACNE (MI) and GENIE3 (random forest)

Plus: topology analysis, perturbation screens, Waddington pseudo-potential landscape, and interactive Plotly visualization.

**GitHub:** https://github.com/junior1p/GRNDynamics
**clawRxiv:** https://clawrxiv.io/abs/2604.01530
**Dependencies:** numpy ≥ 1.24, scipy ≥ 1.10, pandas ≥ 1.5, plotly ≥ 5.15, networkx ≥ 3.1, scikit-learn ≥ 1.3
**Python:** 3.9+, CPU only

---

## Step-by-Step Instructions

### Step 0: Environment Setup

```bash
pip install numpy scipy pandas plotly networkx matplotlib scikit-learn \
    --break-system-packages -q
```

Verify:
```python
python3 -c "import numpy, scipy, pandas, plotly, networkx, sklearn; print('OK')"
```

### Step 1: Network Representation

```python
import numpy as np
import pandas as pd
import networkx as nx
from dataclasses import dataclass, field

@dataclass
class GRN:
    genes: list
    W: np.ndarray                          # signed weight matrix
    boolean_rules: dict = field(default_factory=dict)
    alpha: np.ndarray = None
    gamma: np.ndarray = None
    hill_n: float = 2.0
    hill_k: float = 0.5

    def __post_init__(self):
        N = len(self.genes)
        if self.alpha is None:
            self.alpha = np.ones(N)
        if self.gamma is None:
            self.gamma = np.ones(N)

    @property
    def N(self):
        return len(self.genes)

    def to_networkx(self) -> nx.DiGraph:
        G = nx.DiGraph()
        for i, target in enumerate(self.genes):
            for j, regulator in enumerate(self.genes):
                w = self.W[i, j]
                if abs(w) > 1e-10:
                    G.add_edge(regulator, target, weight=abs(w),
                               sign="+" if w > 0 else "-")
        return G


def grn_from_dict(genes: list, edges: list) -> GRN:
    """Build GRN from (regulator, target, weight) triples. weight>0=activation, <0=repression."""
    N = len(genes)
    gene_idx = {g: i for i, g in enumerate(genes)}
    W = np.zeros((N, N))
    for reg, tgt, w in edges:
        W[gene_idx[tgt], gene_idx[reg]] = w
    return GRN(genes=genes, W=W)


def grn_from_csv(path: str) -> GRN:
    df = pd.read_csv(path)
    genes = sorted(set(df["regulator"]) | set(df["target"]))
    edges = list(zip(df["regulator"], df["target"], df["weight"]))
    return grn_from_dict(genes, edges)
```

### Step 2: Built-in Demo Networks

```python
def demo_pluripotency_network() -> GRN:
    """OCT4-SOX2-NANOG pluripotency + CDX2/GATA6/FGF4. Bistable."""
    return grn_from_dict(
        genes=["OCT4", "SOX2", "NANOG", "CDX2", "GATA6", "FGF4"],
        edges=[
            ("OCT4",  "SOX2",  +2.0), ("SOX2",  "OCT4",  +2.0),
            ("OCT4",  "NANOG", +1.5), ("SOX2",  "NANOG", +1.5),
            ("NANOG", "OCT4",  +1.0), ("NANOG", "SOX2",  +1.0),
            ("NANOG", "NANOG", +0.8), ("OCT4",  "OCT4",  +0.5),
            ("CDX2",  "CDX2",  +1.5),
            ("NANOG", "CDX2",  -2.5), ("CDX2",  "NANOG", -2.5),
            ("OCT4",  "CDX2",  -1.0),
            ("OCT4",  "GATA6", +0.5), ("GATA6", "NANOG", -1.0),
            ("NANOG", "GATA6", -1.0), ("GATA6", "GATA6", +1.2),
            ("OCT4",  "FGF4",  +1.0), ("SOX2",  "FGF4",  +1.0),
            ("FGF4",  "GATA6", +1.5), ("FGF4",  "NANOG", -0.5),
        ]
    )

def demo_EMT_network() -> GRN:
    """EMT circuit: SNAI1/ZEB1/miR200. Tristable (epithelial/mesenchymal/hybrid)."""
    return grn_from_dict(
        genes=["SNAI1", "ZEB1", "ZEB2", "miR200", "miR34", "TWIST1", "ECAD", "VIM"],
        edges=[
            ("SNAI1", "miR34",  -2.5), ("miR34",  "SNAI1", -2.5),
            ("SNAI1", "ECAD",   -2.0), ("SNAI1",  "VIM",   +2.0),
            ("SNAI1", "ZEB1",   +1.0), ("SNAI1",  "SNAI1", +0.5),
            ("ZEB1",  "miR200", -2.5), ("miR200", "ZEB1",  -2.5),
            ("ZEB2",  "miR200", -2.0), ("miR200", "ZEB2",  -2.0),
            ("ZEB1",  "ECAD",   -2.0), ("ZEB1",   "VIM",   +1.5),
            ("ZEB1",  "SNAI1",  +0.8), ("ZEB1",   "ZEB1",  +1.0),
            ("miR34", "miR200", +0.5), ("miR200", "miR34", +0.3),
            ("TWIST1","SNAI1",  +1.2), ("TWIST1", "ZEB1",  +1.0),
            ("TWIST1","ECAD",   -1.5), ("TWIST1", "VIM",   +1.5),
        ]
    )

def demo_cell_cycle_network() -> GRN:
    """Cell cycle G1/S switch: Rb/E2F/CyclinE. Bistable."""
    return grn_from_dict(
        genes=["CycD", "CDK4", "Rb", "E2F", "CycE", "CDK2", "p21", "p27", "Myc"],
        edges=[
            ("CycD",  "CDK4",  +2.0), ("CDK4",  "Rb",    -2.0),
            ("Rb",    "E2F",   -2.5), ("E2F",   "CycE",  +2.0),
            ("CycE",  "CDK2",  +2.0), ("CDK2",  "Rb",    -1.5),
            ("E2F",   "E2F",   +1.0),
            ("p21",   "CDK2",  -2.0), ("p27",   "CDK4",  -1.5),
            ("p27",   "CDK2",  -1.5),
            ("Myc",   "CycD",  +1.5), ("Myc",   "CDK4",  +1.0),
            ("Myc",   "p27",   -1.0), ("E2F",   "Myc",   +0.8),
            ("Rb",    "Myc",   -1.0), ("CycD",  "p27",   -0.5),
        ]
    )
```

### Step 3: Boolean Attractor Search

```python
from itertools import product
from collections import defaultdict

def build_boolean_rules_from_W(grn: GRN, threshold: float = 0.5) -> dict:
    """Auto-generate Boolean update rules from weight matrix W.
    Gene i is ON if any activator is ON AND no strong repressor is ON."""
    rules = {}
    for i, gene in enumerate(grn.genes):
        activators = [grn.genes[j] for j in range(grn.N) if grn.W[i, j] > threshold]
        repressors  = [grn.genes[j] for j in range(grn.N) if grn.W[i, j] < -threshold]
        has_self_act = grn.W[i, i] > threshold

        def make_rule(acts, reps, self_act, idx):
            def rule(state: dict) -> bool:
                if any(state[r] for r in reps):
                    return False
                if acts:
                    return any(state[a] for a in acts) or (self_act and state[grn.genes[idx]])
                if self_act:
                    return state[grn.genes[idx]]
                return False
            return rule

        rules[gene] = make_rule(activators, repressors, has_self_act, i)
    return rules


def synchronous_update(state: dict, rules: dict) -> dict:
    return {gene: rules[gene](state) for gene in state}


def state_to_int(state: dict, genes: list) -> int:
    result = 0
    for i, g in enumerate(genes):
        if state[g]:
            result |= (1 << i)
    return result


def int_to_state(n: int, genes: list) -> dict:
    return {g: bool(n & (1 << i)) for i, g in enumerate(genes)}


def find_all_attractors_boolean(grn: GRN, rules: dict = None, max_steps: int = 200) -> dict:
    """Exhaustive Boolean attractor search. O(2^N). Feasible for N ≤ 25."""
    if rules is None:
        rules = build_boolean_rules_from_W(grn)

    N = grn.N
    total_states = 2 ** N
    if total_states > 2 ** 25:
        raise ValueError(f"Network too large: {N} genes = {total_states:,} states. Use random sampling.")

    genes = grn.genes
    print(f"  Boolean exhaustive: {N} genes, {total_states:,} states")

    successor = np.zeros(total_states, dtype=np.int64)
    for s in range(total_states):
        state_dict = int_to_state(s, genes)
        successor[s] = state_to_int(synchronous_update(state_dict, rules), genes)

    visited = np.full(total_states, -1, dtype=np.int64)
    attractor_id = np.full(total_states, -1, dtype=np.int64)
    attractors = []
    basins = defaultdict(set)

    for start in range(total_states):
        if visited[start] >= 0:
            continue
        path = []
        seen = {}
        s = start
        while s not in seen:
            if visited[s] >= 0:
                aid = attractor_id[s]
                for ps in path:
                    visited[ps] = aid
                    attractor_id[ps] = aid
                    basins[aid].add(ps)
                break
            seen[s] = len(path)
            path.append(s)
            s = successor[s]
        else:
            cycle_start_idx = seen[s]
            cycle = path[cycle_start_idx:]
            pre_cycle = path[:cycle_start_idx]
            aid = len(attractors)
            cycle_states = [int_to_state(c, genes) for c in cycle]
            if len(cycle) == 1:
                attractors.append({"type": "fixed_point", "state": cycle_states[0], "id": aid, "int": cycle[0]})
            else:
                attractors.append({"type": "limit_cycle", "length": len(cycle), "states": cycle_states, "id": aid})
            for c in cycle:
                visited[c] = aid
                attractor_id[c] = aid
                basins[aid].add(c)
            for ps in pre_cycle:
                visited[ps] = aid
                attractor_id[ps] = aid
                basins[aid].add(ps)

    fixed_points = [a for a in attractors if a["type"] == "fixed_point"]
    limit_cycles = [a for a in attractors if a["type"] == "limit_cycle"]

    print(f"  Found {len(fixed_points)} fixed points, {len(limit_cycles)} limit cycles")
    for a in fixed_points:
        on_genes = [g for g, v in a["state"].items() if v]
        basin_size = len(basins[a["id"]])
        print(f"    FP{a['id']}: ON={on_genes}  basin={basin_size}/{total_states} ({100*basin_size/total_states:.1f}%)")

    return {"attractors": attractors, "fixed_points": fixed_points,
            "limit_cycles": limit_cycles, "basins": dict(basins), "n_states": total_states}
```

### Step 4: ODE Dynamics

```python
from scipy.integrate import solve_ivp
from scipy.linalg import eigvals

def hill_activation(x: float, k: float = 0.5, n: float = 2.0) -> float:
    xn = max(x, 0.0) ** n
    kn = k ** n
    return xn / (kn + xn)

def hill_repression(x: float, k: float = 0.5, n: float = 2.0) -> float:
    return 1.0 - hill_activation(x, k, n)

def grn_odes(t: float, x: np.ndarray, grn: GRN) -> np.ndarray:
    N = grn.N
    dxdt = np.zeros(N)
    x_pos = np.maximum(x, 0.0)
    for i in range(N):
        reg_input = 1.0
        has_any = False
        for j in range(N):
            w = grn.W[i, j]
            if abs(w) < 1e-10:
                continue
            has_any = True
            k_eff = grn.hill_k / abs(w)
            if w > 0:
                reg_input *= hill_activation(x_pos[j], k=k_eff, n=grn.hill_n)
            else:
                reg_input *= hill_repression(x_pos[j], k=k_eff, n=grn.hill_n)
        if not has_any:
            reg_input = 0.0
        dxdt[i] = grn.alpha[i] * reg_input - grn.gamma[i] * x_pos[i]
    return dxdt

def integrate_ode_trajectory(grn: GRN, x0: np.ndarray, t_end: float = 50.0, n_timepoints: int = 500) -> tuple:
    t_eval = np.linspace(0, t_end, n_timepoints)
    sol = solve_ivp(fun=lambda t, x: grn_odes(t, x, grn), t_span=(0, t_end), y0=x0,
                     method="RK45", t_eval=t_eval, rtol=1e-6, atol=1e-9, max_step=t_end/100)
    return sol.t, sol.y.T

def find_ode_attractors(grn: GRN, n_random_starts: int = 200, t_end: float = 100.0,
                         cluster_tol: float = 0.05, rng_seed: int = 42) -> list:
    rng = np.random.default_rng(rng_seed)
    N = grn.N
    endpoints = []
    print(f"  ODE attractor search: {n_random_starts} trajectories")
    for trial in range(n_random_starts):
        x0 = rng.uniform(0, 2, N)
        try:
            t, X = integrate_ode_trajectory(grn, x0, t_end=t_end, n_timepoints=100)
            endpoints.append(X[-1])
        except Exception:
            pass

    endpoints = np.array(endpoints)
    clusters = []
    assigned = np.zeros(len(endpoints), dtype=bool)
    for i, ep in enumerate(endpoints):
        if assigned[i]:
            continue
        dists = np.linalg.norm(endpoints - ep, axis=1)
        close = dists < cluster_tol * N
        if close.sum() > 0:
            center = endpoints[close].mean(axis=0)
            clusters.append({"center": center, "members": np.where(close)[0], "count": int(close.sum())})
            assigned[close] = True

    def numerical_jacobian(x_eq: np.ndarray, eps: float = 1e-5) -> np.ndarray:
        N = len(x_eq)
        J = np.zeros((N, N))
        f0 = grn_odes(0, x_eq, grn)
        for j in range(N):
            x_pert = x_eq.copy()
            x_pert[j] += eps
            f_pert = grn_odes(0, x_pert, grn)
            J[:, j] = (f_pert - f0) / eps
        return J

    attractors = []
    for c in sorted(clusters, key=lambda x: -x["count"]):
        x_eq = np.maximum(c["center"], 0)
        J = numerical_jacobian(x_eq)
        eigs = eigvals(J)
        max_real_eig = float(np.max(eigs.real))
        stability = "stable" if max_real_eig < -1e-3 else ("unstable" if max_real_eig > 1e-3 else "saddle")
        on_genes = [grn.genes[i] for i in range(N) if x_eq[i] > 0.5]
        attractors.append({"state": x_eq, "eigenvalues": eigs, "max_real_eigenvalue": max_real_eig,
                           "stability": stability, "on_genes": on_genes, "basin_fraction": c["count"]/len(endpoints)})

    stable = [a for a in attractors if a["stability"] == "stable"]
    print(f"  {len(attractors)} attractors ({len(stable)} stable)")
    for a in attractors:
        print(f"    [{a['stability'].upper()}] ON={a['on_genes']}  basin~{100*a['basin_fraction']:.1f}%")
    return attractors
```

### Step 5: GRN Inference (from expression data)

```python
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler

def _mutual_information_discrete(x: np.ndarray, y: np.ndarray, n_bins: int) -> float:
    joint = np.zeros((n_bins, n_bins))
    for xi, yi in zip(x, y):
        joint[min(xi, n_bins-1), min(yi, n_bins-1)] += 1
    joint /= (joint.sum() + 1e-300)
    px = joint.sum(axis=1)
    py = joint.sum(axis=0)
    mi = 0.0
    for i in range(n_bins):
        for j in range(n_bins):
            if joint[i, j] > 1e-300 and px[i] > 1e-300 and py[j] > 1e-300:
                mi += joint[i, j] * np.log2(joint[i, j] / (px[i] * py[j]))
    return max(mi, 0.0)

def infer_grn_mutual_information(expr_matrix: np.ndarray, gene_names: list,
                                  n_bins: int = 10, dpi_threshold: float = 0.1,
                                  mi_threshold_quantile: float = 0.85) -> GRN:
    """ARACNE-style GRN inference via mutual information + DPI."""
    n_samples, N = expr_matrix.shape
    print(f"  ARACNE: {n_samples} samples × {N} genes")
    expr_disc = np.zeros_like(expr_matrix, dtype=int)
    for g in range(N):
        percentiles = np.percentile(expr_matrix[:, g], np.linspace(0, 100, n_bins + 1))
        expr_disc[:, g] = np.digitize(expr_matrix[:, g], percentiles[1:-1])

    MI = np.zeros((N, N))
    for i in range(N):
        for j in range(i + 1, N):
            mi = _mutual_information_discrete(expr_disc[:, i], expr_disc[:, j], n_bins)
            MI[i, j] = MI[j, i] = mi

    W_mi = MI.copy()
    for i in range(N):
        for j in range(i + 1, N):
            for k in range(N):
                if k == i or k == j:
                    continue
                if MI[i, j] < min(MI[i, k], MI[j, k]) - dpi_threshold:
                    W_mi[i, j] = W_mi[j, i] = 0.0
                    break

    nonzero_mi = W_mi[W_mi > 0]
    if len(nonzero_mi) > 0:
        threshold = np.quantile(nonzero_mi, mi_threshold_quantile)
        W_mi[W_mi < threshold] = 0.0

    corr = np.corrcoef(expr_matrix.T)
    W_signed = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            if abs(W_mi[i, j]) > 1e-10:
                sign = 1.0 if corr[i, j] > 0 else -1.0
                W_signed[i, j] = sign * W_mi[i, j]
    print(f"  ARACNE: {int((W_signed != 0).sum())} edges")
    return GRN(genes=gene_names, W=W_signed)

def infer_grn_genie3(expr_matrix: np.ndarray, gene_names: list,
                      n_estimators: int = 100, weight_threshold_quantile: float = 0.80) -> GRN:
    """GENIE3: Random Forest feature importance for each gene predicted from all others."""
    n_samples, N = expr_matrix.shape
    print(f"  GENIE3: {N} genes")
    W = np.zeros((N, N))
    scaler = StandardScaler()
    expr_scaled = scaler.fit_transform(expr_matrix)

    for i, target_gene in enumerate(gene_names):
        y = expr_scaled[:, i]
        X_reg = np.delete(expr_scaled, i, axis=1)
        regulator_idx = [j for j in range(N) if j != i]
        rf = RandomForestRegressor(n_estimators=n_estimators, max_features="sqrt", random_state=42, n_jobs=1)
        rf.fit(X_reg, y)
        for k, j in enumerate(regulator_idx):
            W[i, j] = rf.feature_importances_[k]

    threshold = np.quantile(W[W > 0], weight_threshold_quantile) if (W > 0).any() else 0
    corr = np.corrcoef(expr_scaled.T)
    W_signed = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            if W[i, j] >= threshold:
                sign = 1.0 if corr[i, j] > 0 else -1.0
                W_signed[i, j] = sign * W[i, j]
    print(f"  GENIE3: {int((W_signed != 0).sum())} edges")
    return GRN(genes=gene_names, W=W_signed)
```

### Step 6: Topology Analysis

```python
def analyze_network_topology(grn: GRN) -> dict:
    N = grn.N
    W = grn.W
    genes = grn.genes
    out_degree = (np.abs(W) > 0.1).sum(axis=0)
    in_degree  = (np.abs(W) > 0.1).sum(axis=1)
    master_reg_score = out_degree - in_degree
    master_regulators = [(genes[i], int(out_degree[i]), int(in_degree[i]))
                         for i in np.argsort(master_reg_score)[::-1][:5]]
    autoact = [genes[i] for i in range(N) if W[i, i] > 0.1]
    autorep = [genes[i] for i in range(N) if W[i, i] < -0.1]

    pos_feedback, neg_feedback = [], []
    for i in range(N):
        for j in range(i + 1, N):
            if abs(W[i, j]) > 0.1 and abs(W[j, i]) > 0.1:
                pair = (genes[i], genes[j]) if i < j else (genes[j], genes[i])
                if np.sign(W[i, j]) * np.sign(W[j, i]) > 0:
                    pos_feedback.append(pair)
                else:
                    neg_feedback.append(pair)

    results = {
        "n_genes": N, "n_edges": int((np.abs(W) > 0.1).sum()),
        "out_degree": {genes[i]: int(out_degree[i]) for i in range(N)},
        "in_degree":  {genes[i]: int(in_degree[i])  for i in range(N)},
        "master_regulators": master_regulators,
        "autoactivation": autoact, "autorepression": autorep,
        "positive_feedback_loops": pos_feedback, "negative_feedback_loops": neg_feedback,
        "density": (np.abs(W) > 0.1).sum() / max(N * (N - 1), 1),
    }
    print(f"  Topology: {N} genes, {results['n_edges']} edges, "
          f"{len(pos_feedback)} pos FB, {len(neg_feedback)} neg FB, "
          f"master={master_regulators[0][0] if master_regulators else 'N/A'}")
    return results
```

### Step 7: Waddington Pseudo-Potential Landscape

```python
from sklearn.decomposition import PCA
from scipy.stats import gaussian_kde

def compute_waddington_potential(grn: GRN, resolution: int = 40) -> tuple:
    """Compute 2D quasi-potential landscape via ODE trajectory density + PCA projection."""
    rng = np.random.default_rng(0)
    N = grn.N
    endpoints = []
    for _ in range(resolution * resolution // 2):
        x0 = rng.uniform(0, 2, N)
        try:
            t, X = integrate_ode_trajectory(grn, x0, t_end=60.0, n_timepoints=60)
            endpoints.append(X[-1])
        except:
            pass
    if len(endpoints) < 10:
        return None, None, None
    endpoints = np.array(endpoints)
    pca = PCA(n_components=2)
    coords_2d = pca.fit_transform(endpoints)
    kde = gaussian_kde(coords_2d.T, bw_method=0.3)
    x_min, x_max = coords_2d[:, 0].min() - 0.5, coords_2d[:, 0].max() + 0.5
    y_min, y_max = coords_2d[:, 1].min() - 0.5, coords_2d[:, 1].max() + 0.5
    xi = np.linspace(x_min, x_max, resolution)
    yi = np.linspace(y_min, y_max, resolution)
    Xi, Yi = np.meshgrid(xi, yi)
    density = kde(np.vstack([Xi.ravel(), Yi.ravel()])).reshape(resolution, resolution)
    potential = -np.log(density + 1e-10)
    potential -= potential.min()
    return Xi, Yi, potential
```

### Step 8: Full Pipeline

```python
import json, os

def run_grndynamics(grn: GRN = None, demo: str = "pluripotency",
                    out_dir: str = "grn_output",
                    run_boolean: bool = True, run_ode: bool = True,
                    n_ode_starts: int = 150) -> dict:
    """Run the full GRNDynamics pipeline."""
    os.makedirs(out_dir, exist_ok=True)
    if grn is None:
        demos = {"pluripotency": demo_pluripotency_network,
                 "emt": demo_EMT_network,
                 "cell_cycle": demo_cell_cycle_network}
        grn = demos[demo]()
    print(f"\n{'='*50}\n  GRNDynamics — {grn.N} genes\n{'='*50}")

    results = {}
    print("\n[1/3] Topology..."); topology = analyze_network_topology(grn); results["topology"] = topology

    bool_result = {}
    if run_boolean:
        print(f"\n[2/3] Boolean ({2**grn.N:,} states)...")
        rules = build_boolean_rules_from_W(grn)
        bool_result = find_all_attractors_boolean(grn, rules)
        results["boolean"] = {"n_fixed_points": len(bool_result.get("fixed_points", []))}

    ode_attractors = []
    if run_ode:
        print(f"\n[3/3] ODE ({n_ode_starts} trajectories)...")
        ode_attractors = find_ode_attractors(grn, n_random_starts=n_ode_starts)
        results["ode"] = {"n_stable": len([a for a in ode_attractors if a["stability"] == "stable"])}

    # Visualization
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots
    Xi, Yi, Zi = compute_waddington_potential(grn) if ode_attractors else (None, None, None)

    fig = make_subplots(rows=2, cols=2,
        subplot_titles=["Network", "Boolean Attractor Basins", "ODE Trajectories", "Waddington Landscape"])
    G = grn.to_networkx()
    pos = nx.spring_layout(G, seed=42, k=2.0)
    for edge in G.edges(data=True):
        x0, y0 = pos[edge[0]]; x1, y1 = pos[edge[1]]
        fig.add_trace(go.Scatter(x=[x0, x1, None], y=[y0, y1, None], mode="lines",
                                 line=dict(color="#e74c3c" if edge[2]["sign"] == "+" else "#3498db"),
                                 showlegend=False), row=1, col=1)
    fig.add_trace(go.Scatter(x=[pos[g][0] for g in grn.genes], y=[pos[g][1] for g in grn.genes],
                             mode="markers+text", marker=dict(size=14, color="steelblue"),
                             text=grn.genes, textposition="top center"), row=1, col=1)
    if bool_result.get("fixed_points"):
        fps = bool_result["fixed_points"]
        basins = bool_result.get("basins", {})
        ns = bool_result.get("n_states", 1)
        fig.add_trace(go.Bar(x=[f"FP{k}" for k in range(len(fps))],
                            y=[100*len(basins.get(fps[k]["id"], []))/ns for k in range(len(fps))],
                            marker_color=["#e74c3c","#3498db","#2ecc71","#f39c12","#9b59b6"]),
                         row=1, col=2)
    rng = np.random.default_rng(1)
    for trial in range(4):
        x0 = rng.uniform(0, 2, grn.N)
        try:
            t, X = integrate_ode_trajectory(grn, x0, t_end=40.0)
            for gi in range(min(3, grn.N)):
                fig.add_trace(go.Scatter(x=t, y=X[:, gi], mode="lines", showlegend=(trial==0),
                                        name=grn.genes[gi], line=dict(width=1.5)), row=2, col=1)
        except: pass
    if Xi is not None:
        fig.add_trace(go.Surface(x=Xi, y=Yi, z=Zi, colorscale="RdBu_r", showscale=False), row=2, col=2)
    fig.update_layout(height=800, title=f"GRNDynamics: {grn.N} genes — {demo}", template="plotly_white")
    fig.write_html(f"{out_dir}/grn_dynamics.html")
    print(f"\n  Saved: {out_dir}/grn_dynamics.html")

    with open(f"{out_dir}/summary.json", "w") as f:
        json.dump(results, f, indent=2, default=str)
    return results
```

### Step 9: Quick Start

```python
# Demo run
grn = demo_pluripotency_network()   # or demo_EMT_network(), demo_cell_cycle_network()
results = run_grndynamics(grn, demo="pluripotency", out_dir="grn_output")

# Custom GRN
grn = grn_from_dict(
    genes=["MYC", "BCL2", "BAX"],
    edges=[("MYC", "BCL2", +1.5), ("BCL2", "BAX", -1.0), ("BAX", "MYC", -0.5)]
)
results = run_grndynamics(grn, out_dir="custom_output")

# From CSV
# CSV needs columns: regulator, target, weight
# grn = grn_from_csv("my_network.csv")
# results = run_grndynamics(grn)
```

---

## Output Files

| File | Description |
|------|-------------|
| `grn_dynamics.html` | Interactive 4-panel Plotly visualization |
| `summary.json` | Machine-readable results (topology + attractors) |

---

## Scientific Background

**Waddington epigenetic landscape**: Cells differentiate by rolling from a pluripotent hilltop into differentiated valleys (attractors). Each attractor = a stable gene expression pattern = a cell type. GRNDynamics quantifies this.

**Bistability**: Positive feedback loops (A→B, B→A) or double-negative feedback (A⊣B, B⊣A) create two stable states. OCT4-NANOG vs CDX2 mutual repression = canonical example.

**Limit cycles**: Negative feedback loops (odd number of repressions in a cycle) can produce oscillations — basis of circadian rhythms, cell cycle progression, somitogenesis.

**Master regulators**: Genes with high out-degree and low in-degree — they control many targets but are themselves weakly controlled. NANOG is the top master regulator in the pluripotency network.

---

## Limitations & Pitfalls

1. **Boolean enumeration scales as 2^N** — feasible up to N=25. For larger networks use random sampling (not yet in main pipeline, add `n_trajectories=5000` to `random_boolean_attractor_search`).
2. **ODE parameters (α, γ, k, n) need tuning** — default values work for demo networks but may need adjustment for custom networks.
3. **Network inference methods are correlation-based** — they cannot distinguish direct from indirect regulation without perturbation/intervention data.
4. **Waddington landscape is a 2D projection** — the true landscape is N-dimensional; the projection may obscure some structure.
5. **Asynchronous update** — current implementation uses synchronous update for exhaustive search; asynchronous gives different results for some networks.

---

## References

1. Kauffman, S. (1969). Metabolic stability and epigenesis in randomly constructed genetic nets. *Journal of Theoretical Biology*.
2. Chickarmane, V. et al. (2006). Transcriptional dynamics of the embryonic stem cell switch. *PLoS Computational Biology*.
3. Lu, M. et al. (2013). MicroRNA-based regulation of EMT and cancer stem cells. *PNAS*.
4. Huang, S. et al. (2007). Bifurcation dynamics in lineage-commitment. *Developmental Cell*.
5. Huynh-Thu, V.A. et al. (2010). Inferring regulatory networks from expression data. *PLoS ONE*.
6. Tyson, J.J. & Novak, B. (2015). Bistability, oscillations, and excitability in the vertebrate cell cycle. *Current Biology*.
7. Boyer, L.A. et al. (2005). Core transcriptional regulatory circuitry in embryonic stem cells. *Cell*.