C. elegans Connectomics + Transcriptomics

Multi-modal analysis asking how much of local wiring in a small nervous system is determined by topology vs. gene expression.

Mon Sep 01 2025 00:00:00 GMT+0000 (Coordinated Universal Time)

Thesis

Local connectomic topology in C. elegans is largely self-determining, while gene expression acts as a modulatory refinement layer rather than a primary architect of wiring.

Interactive: the hermaphrodite connectome

Click any neuron class to highlight its immediate neighborhood — outgoing chemical synapses and gap junctions. Colors encode the 3-class behavioral role (sensory / inter / motor) that the expression-only classifier predicts at 92.7% accuracy in this project.

What I built

A multi-modal pipeline over the hermaphrodite C. elegans connectome (N=91 somatic neurons with both connectivity and single-cell expression data):

  • GNCA (Graph Neural Cellular Automata) predicting outgoing synaptic strength — r = 0.987 on held-out edges. Controls: shuffled-expression r = 0.934 (topology carries the signal), graph-only r = 0.861, expression-only r = 0.27.
  • Neuron-role classifier — 92.7% accuracy predicting 3-class behavioral role (sensory / inter / motor) from expression alone.
  • CCA between expression space and motif space: canonical correlations ρ = [0.875, 0.836, 0.809].
  • Honest null — a PhaseD hybrid gene-graph LSTM showed gene_causal_by_shuffle: False: dynamics prediction doesn’t improve with gene expression on this dataset.
Genetic collapse curve: model correlation (red) vs. cumulative top-driver gene knockouts, alongside average probability assigned to hub neurons (navy). Performance plateaus well above chance even after 200 genes removed — evidence of genetic redundancy supporting the topology-dominant thesis.

Methodology notes

  • Shuffle / ablation / null controls throughout — the thesis rests on the gap between topology-only and expression-augmented models.
  • Triangulation across CCA (linear), GNCA (nonlinear message passing), and classifier probes.
  • Cross-species Caenorhabditis and parasitic-nematode genomes loaded for a planned comparative follow-up.

Status

In progress. Draft manuscript structured around the topology-dominant / gene-modulatory framing. Target venues: eLife, PLOS Computational Biology, or Network Neuroscience for the main paper; NeurIPS GRL / ICLR LMRL workshops for the GNCA methods paper.

Repo is private while the work is unfinished.

Digital worm — embodied simulator in progress

A long-term effort to build a browser-native, expression-aware C. elegans body simulator that plugs into the connectome-constrained brain model the analyses above feed into. The body is a MuJoCo MJCF (/data/wormbody.xml) — 21 bodies, 20 joints, 19 muscle actuators — and three controllers are running on it below:

  • Kinematic (live) — browser-side wave-equation model with interactive sliders. Fast, smooth, not physics-constrained.
  • Physics (MuJoCo CPG) — the MJCF body driven by a tuned central-pattern-generator controller with anisotropic resistive-force-theory drag (c⊥/c∥ = 2.0). The anisotropy is what converts lateral undulation into forward thrust. Crawls at ~145 µm/s, in the biological range.
  • Imitation (RL-trained) — MLP policy trained with PPO against a biologically parameterised reference (Fang-Yen 2010 / Boyle-Berri-Cohen 2012 crawling parameters) using MSE-of-segment-positions as reward. The action is a residual on top of the CPG baseline, so training only has to learn the physics-vs-kinematics correction, not the whole gait. 120k env steps of training, final rollout crawls at ~117 µm/s with peak per-segment curvature of ~22°.

Next phases: (3) swap the controller for a Brian2 LIF brain built from the Cook 2019 connectome + the GNCA-predicted synaptic weights (the C. elegans analog of Shiu et al. 2024’s fly brain model); (4) connectome-and-expression-constrained perturbation experiments against known knockout phenotypes from the C. elegans behavioral genetics literature.

Brain-driven: closed-loop demo (Phase 3)

Below is the Phase 3 closed loop. A 300-neuron Brian2 LIF brain (Cook 2019 hermaphrodite wiring, Loer & Rand 2022 neurotransmitter identity) drives the MuJoCo body through event classifiers fit from Atanas et al. 2023’s paired whole-brain calcium + behavior recordings (DANDI 000776).

Architecture mirrors the Shiu et al. 2024 fly-brain pattern: sensory input → LIF brain → motor-neuron readout → event classifier bank → behavioural state machine → state-specific CPGs. The brain-body sync cadence follows the Eon Systems 2026 embodied-fly integration.

The dashboard runs five pre-rendered scenarios (spontaneous / head-touch / osmotic-shock / food / chemotaxis) with a synchronised 20-segment MuJoCo body, 300-neuron 3D brain view, 9-modulator concentration strip, FSM behavioural timeline with event-classifier carets, and a live arena for the chemotaxis scenario. Interactive hooks: click a neuron to lock its 1-hop connectivity + firing-rate history; hover a modulator row or circuit badge to highlight its neurons; scrub the timeline or use the per-scenario “jump to moment” buttons to seek to biologically-relevant windows (reversal onset, FLP-11 surge, dwelling). Keyboard: space/arrows/1–5/?/R. Press ? inside the dashboard for the full guide.

Four scenarios are pre-rendered so the closed-loop simulation runs offline (a live Brian2 + MuJoCo loop in the browser would need MuJoCo-wasm + Brian2Cuda-wasm, not there yet):

  • Spontaneous — no stimulus; baseline behaviour.
  • Head touch — ALM/AVM mechanoreceptor drive at t=5s (Chalfie 1981).
  • Osmotic shock — ASH polymodal avoidance drive at t=5s (Hart 1995). Should produce the most reversal of any scenario.
  • Food — ASI/ASJ/ADF feeding-state tonic drive from t=2s. Should produce the least reversal.

Honest calibration notes are in the component’s “Sources & attribution” fold. v1.5 added a per-neuron affine distribution calibration that matches Brian2 synthetic-calcium moments to the Atanas ΔF/F training distribution. v3 adds an actual slow-neuromodulation layer — 9 peptidergic + monoaminergic concentrations (FLP-11, FLP-1, FLP-2, NLP-12, PDF-1, 5-HT, dopamine, tyramine, octopamine) with releaser + receptor tables extracted directly from CeNGEN single-cell expression (Taylor et al. 2021). This gives the model a causal substrate for sleep-like quiescence (RIS→FLP-11→broad inhibition via NPR-1/NPR-22/DMSR-1, Turek 2016) rather than just a statistical correlate via NSM.

Perturbation validation and reproducibility audit (Phases 3d-3 through 3d-6). An initial v3.0 perturbation suite reported two apparent phenotype reproductions (RIS → quiescence collapse ΔQUI=−0.53; AVA → reversal abolished ΔREV=−0.15), and subsequent biological corrections (5HT pharyngeal exclusion in 3d-5; CeNGEN-derived per-edge glutamate receptors in 3d-4) appeared to regress both results. Rather than iterate on recalibration, the next phase (3d-6) ran a proper reproducibility audit: 3 seeds × 3 configs × 6 ablations × 2 conditions (108 runs, ~97 min compute) with brian2.seed() locked alongside np.random.seed(). The audit overturned the initial story in both directions.

What survives statistical scrutiny (v3.0, default sign mode, n=3 seeds × 20 s): AVA ablation under touch reliably reduces reversal — ΔREV = −0.57 ± 0.37 across 3 seeds, all three correctly negative. This was originally framed as a genuine Chalfie 1985 reproduction. The RIS result, previously reported as ΔQUI=−0.53 flagship, averages to −0.24 ± 0.33 across seeds — directionally consistent with Turek 2016 (2/3 seeds negative) but not statistically robust at 20 s runs. Subsequent horizontal-rebase work (May 2026) reframes both findings — see “Horizontal rebase: honest reframing of phenotype claims” below.

What the biological corrections actually did: the audit confirms they systematically degrade phenotype reproduction. The AVA effect weakens to ΔREV=−0.29±0.52 in v3.1 and disappears entirely in v3.2 (−0.03±0.13). This is statistically confirmed regression across multiple seeds, not a single-run artifact. Adding biologically-accurate modulation shifts the network out of the regime where calibrated phenotypes emerge.

v3.3 scope from the audit: all phenotype claims will require n≥5 seeds × 60 s runs. Modulation strengths need recalibration to recover AVA reproduction under the corrected-biology configurations. NSM, RIM, PDE, AVB ablation claims are dropped — they were all noise at 20 s runs. Full audit in artifacts/ensemble_report.md.

Remaining work: classifier joint/correlation structure mismatch between Brian2 and Atanas; per-edge glutamate receptor mapping available as a toggle pending the v3.3 recal; 5HT target weights need circuit-role reweighting (not just pharyngeal filtering) to capture dwelling; longer simulation runs and more seeds for all future phenotype tests.

Tier 1 biological-accuracy upgrades (Phases T1a–T1e, shipped). After the audit, I prioritised five architectural corrections aimed at closing the biggest first-principles gaps in the v3 model:

  • T1a — Graded non-spiking dynamics (Kunert-Graf 2014). Replaces LIF with the worm-specific σ(V) = 1/(1+exp(-(V-V_half)/k)) graded-release formulation. Most C. elegans neurons don’t fire vertebrate-style sodium action potentials (Goodman, Hall & Avery 1998) — the LIF spiking model was always biologically wrong at first principles. New GradedBrain class with current-based graded synapses.
  • T1b — L-type Ca channels for plateau potentials. Voltage-gated Ca conductance (egl-19-like) added to 14 command/plateau neurons: AVA, AVE, AVD, AVB, PVC, AIY, RIS, DVA. This makes reversal-bout sustainment an emergent property of neural dynamics (Wicks 1996 plateau mechanism) rather than something imposed by FSM state-hold minimums.
  • T1c — Volume-transmission modulator geometry. 3D soma coordinates extracted from OpenWorm morphology data; per-modulator concentration now a distance-weighted field from each releaser with diffusion length scales λ of 400–700 µm for peptides and 150–250 µm for monoamines. RIS’s FLP-11 now preferentially inhibits nearby head neurons rather than affecting the whole connectome uniformly.
  • T1d — Real closed-loop proprioception. Fixed the v1 O(n²) recompile bug by using persistent Brian2 PoissonGroup with mutable rates. Body curvature now drives PDE/PDA/DVA stretch receptors every 50 ms sync; the loop is genuinely closed for the first time.
  • T1e — 2D chemical-gradient environment. Environment + ChemoGradient classes track worm head position from MuJoCo CoM and sample a 2D Gaussian attractant field. ASE + AWA fire as ON-cells (tonic + positive dC/dt), AWC as OFF-cells (Chalasani 2007). Enables Pierce-Shimomura 1999 chemotaxis-index validation — the canonical worm navigation test. No other connectome-constrained simulator has this for any organism.

Full Tier 1 stack available via ClosedLoopEnv(brain_class="graded", environment=…). Shipped demo scenarios remain on the v3 LIF brain.

v3.3 ensemble audit (graded stack vs LIF baseline). Ran the 6-ablation suite through the Tier 1 graded brain at 3 seeds × 20 s per run. All 36 runs produced identical control vs ablated distributions (FWD 0.03, OMG 0.05, QUI 0.92 — FSM locked in QUIESCENT). The graded brain produces σ(V) patterns that fall below the classifier bank’s detection threshold, because the bank was trained on LIF-derived synthetic calcium. The biological correctness of graded dynamics doesn’t translate to phenotype reproduction until the classifier is retrained on graded-derived σ signals. This is v3.4 work (~2 weeks: generate synthetic-Atanas from graded brain under scripted sensory schedules, retrain 8-event classifier bank, retune FSM thresholds, re-run ensemble). Full audit in artifacts/v33_audit_report.md.

T0 sign-convention resolution (April 25 2026). The April 21 framing — that the v3 LIF brain failed to propagate touch through ALM → AIB → AVA, that synaptic weight calibration was the prerequisite fix — turned out to be wrong on both the cascade description and the fix category. AIB has zero chemical edges to AVD in this connectome (Cook 2019 + Loer & Rand 2022); ALM/AVM also have zero direct chemical edges to AIB. The actual operative cascade is ALM/AVM → PVC → AVD/AVE → AVA, and the cascade failure was a sign-assignment problem, not a weight-calibration problem. Default per-presynaptic-neuron sign convention treated glutamate edges to iGluR-dominant cells as inhibitory; per-edge mode (CeNGEN-derived postsynaptic-receptor signs) flips ~518 chemical edges (14% of total) and resolves cascade firing. Under per-edge mode, AVD/AVA ΔHz +60 on touch (n=10, σ < 1.5 Hz). The historical audit findings above (v3.0/3.1/3.2/3.3) stand as observations under default mode; the per-edge mode opens new questions (PVC/AVB over-activation, FSM/classifier recalibration, RIS silencing under per-edge) that are themselves load-bearing. Per-edge is currently opt-in pending resolution.

Wave 2 cellular layer (April 26-28 2026). A separate trajectory built a NEURON→Brian2 channel-translation pipeline for biophysically-validated cellular dynamics beyond what LIF can represent. Four cells production-grade (matching Nicoletti 2024 published phenotypes to 5-decimal-place agreement on full voltage-clamp + current-clamp validation panels): AVAL (4 channels), AIY (7 channels), RIM (7 channels), AVAR (5 channels). 14 NMODL channel translations, F1-F18 NMODL gotcha catalog, cython codegen baseline (22.71× pure-Brian2 speedup). Wave 2 caught a Mellem-2008-as-AVA-anchor citation misattribution that had carried through earlier docs (Mellem 2008 actually characterizes RMD plateau dynamics — direct quote: “we never observed action potentials in AVA”; AVA experimental data traces to Liu/Chen/Wang 2020 Nat Commun via Nicoletti’s references). Stage IV cross-validation independently confirms the §5 T0 resolution claim under conservative drive (touch protocol: ALM/AVM activate ~60 Hz peak; AVA ΔHz +7.5; the cascade fires under per-edge mode regardless of stim parameters within reasonable range). Wave 2’s value proposition is mechanistic biological resolution beyond LIF’s spike-count abstraction, not a fix for a broken cascade — AVAL and AVAR are biologically distinguishable under Wave 2 cellular detail (different rest, different K-channel set including UNC-103, different plateau amplitudes at same drive) where they are interchangeable in LIF.

Wave2HybridBrain integration — production-grade under M2-pure (May 2026). Initial Phase δ WB3 cross-coupling work in April surfaced a load-bearing readout artifact (σ > 0.5 rising-threshold pseudo-spike emission silently zeros under saturated σ regime). The D7-followup implemented σ-magnitude continuous readout matching graded_brain.py output_rates() precedent. Under that corrected readout, the Wave2HybridBrain initially appeared to NOT activate AVA on touch — leading to a temporary “demote Wave 2 hybrid to research substrate” framing. The horizontal rebase’s W2 investigation thread refuted that interpretation: under M2-pure (per-edge + sign_exceptions={}, the locked production sign mode), Wave2HybridBrain fires the cascade identically to pure LIFBrain M2-pure (AVDL Δ +60.5 Hz; AVAL σ Δ +10% peri-touch with clean Δ_post ≈ 0 separation from drift). The initial CP2 finding was a sign-mode question (DOCUMENTED_SIGN_EXCEPTIONS in non-M2-pure modes collapses the cascade upstream of Wave 2 cells), not a cellular-substitution question. Wave 2 + M2-pure adds biophysical resolution beyond LIF without breaking cascade dynamics.

Horizontal rebase: honest reframing of phenotype claims (May 2026)

A five-phase horizontal rebase was conducted to settle the validation methodology questions the T0 resolution had left open. Repository deliverables: docs/state_of_claims_2026-05-02.md (state-of-claims catalog), docs/brain_v3.5_locked.md (Phase 1 brain lock), docs/phase2_preflight.md (architecture sign-off).

Phase 1 ran a sign-mode decision gauntlet across 4 candidate modes. Only M2-pure (per-edge CeNGEN signs with sign_exceptions={}) fired the +60 Hz touch cascade documented in §5. M1 (default, current production) and M2-current (per-edge + 7 DOCUMENTED_SIGN_EXCEPTIONS) both showed AVDL/AVAL DROPPING on touch by 2-4 Hz. Brain locked at M2-pure.

Phase 2 rebuilt the readout layer from scratch: A2-balanced 21-cell readout (legacy 18 + AVAL + AVAR + AVDL — AVB and PVC pairs absent from Atanas data), CV-trained classifier with leave-one-worm-out 10-fold validation, M2-pure-specific Brian2→Atanas calibration, recalibrated FSM thresholds.

Phase 2.5 ran the validation gauntlet under the recalibrated stack at default tier (n=10 × 60s):

  • AVA → dREV = +0.229 ± 0.137 (2/10 negative) under M2-pure. Wrong direction. The Chalfie 1985 phenotype is NOT reproduced under correct cascade dynamics.
  • AVA → dPIR = −0.005 ± 0.005 (1/10 negative). The dPIR channel-shift hypothesis from T0 §5 (-0.117, 9/10 negative under legacy stack) does not survive recalibration — that finding was a stack artifact, not a real channel-shifted phenotype.
  • RIS → dQUI = −0.007 ± 0.026 (4/10 negative). Turek 2016 RIS quiescence reproduction is null at full power.
  • NEW finding (Direct, awaiting literature precedent): AVA → dFWD = −0.302 ± 0.102 (7/10 negative, Cohen’s d ≈ 0.93). Forward-locomotion suppression rather than reversal abolition — biologically interpretable via AVA-AVB gap-junction coupling (Wang/Liu/Chen 2020) but not the textbook Chalfie phenotype. Headline-result decision deferred pending literature precedent verification.

Honest current state: under correct cascade dynamics + recalibrated readout, this simulator does not reproduce Chalfie 1985 (AVA → reversal abolition) or Turek 2016 (RIS → quiescence) as originally claimed. It does produce a robust AVA → forward-suppression signal that may correspond to a different published phenotype or may be a methodological artifact awaiting external scrutiny. The previously-cited reproductions were real as statistical measurements but mechanistically misleading: they ran via Mode 3 tonic-shift on broken sign convention (default mode dREV) and via a now-disproved channel-shift hypothesis (per-edge dPIR under legacy classifier).

Statistical robustness is not the same as mechanistic validity. The rebase is an exercise in catching that distinction within one’s own work.