← The living map prototypes / a0-regulatory / README

A0: regulatory closure demonstrator

A small mathematical precursor to Track A and the M3 identifiability gate in the CAKE roadmap. Three conductance-based neurons drive a single reduced body variable, which feeds back into their sensory inputs. Each neuron regulates two positive conductances from its own filtered activity.

Scope: synthetic, dimensionless equations with chosen parameters. This is a test of solver assumptions and experimental identifiability, not a model of a worm, an AFD positive control, an H* verdict, or completion of M3. There are no experimental recordings or fitted biological channel kinetics here.

Run

cd prototypes/a0-regulatory
uv sync --locked
uv run pytest -q
uv run cake-a0 --out results

The CLI exits nonzero if any numerical demonstration gate fails. It writes report.md, a detailed report.json, overview.svg, overview.png, and sampled states in trajectories.npz. Runtime dependencies and test dependencies are pinned by uv.lock. NumPy/SciPy run on the CPU; no GPU or data download is needed.

A checked run is stored in evidence/reference/report.md. Its JSON records parameters, schedules, noise seeds, solver tolerances, package versions, and SHA-256 hashes of the implementation, tests and dependency lock. The saved artifacts are numerical evidence, not a promise of bitwise identity across machines or library versions.

Equations and units

All variables and times are dimensionless. Voltage is normalized to three reversal potentials: leak at 0, channel g at 1, channel h at 2. Both regulated channels are inward at the chosen operating voltages. The model has 13 states: [v(3), s(3), log(g)(3), log(h)(3), b]. Log coordinates enforce positive conductances without clipping.

For neuron i:

tau_v dv_i/dt = -v_i + g_i(1-v_i) + h_i(2-v_i)
                 + sum_j W_ij v_j + I_i + q_i b
tau_s ds_i/dt = v_i - s_i
e_i           = r_i - s_i
d log(g_i)/dt = k e_i
d log(h_i)/dt = k (rho_i e_i + eta e_i^2)

tau_b db/dt   = -b + tanh(sum_i m_i v_i)

Here r is the activity target, k is the shared regulation rate, rho is a neuron-specific ratio of regulatory gains, and eta controls history dependence. W is a signed, phenomenological recurrent current coupling, not a receptor or release model. The first-order body variable represents an abstract posture or position signal; it is not validated mechanics. Its motor-to-sensory loop is active during development, perturbation, fitting and evaluation.

Quantity Default
r (0.25, 0.30, 0.35)
rho (0.7, 1.0, 1.3)
k 0.025
tau_v, tau_s, tau_b 0.08, 0.8, 2.0
Initial g, h 0.08, 0.05 per neuron
Baseline sensory drive I 0.03 per neuron

The full coupling arrays are in Circuit and every run's JSON. The local regulatory rule is constructed for this experiment; the squared-error term is a counterexample permitted by a broad local-rule family, not an asserted biological homeostatic mechanism.

The history rule uses deterministic sensors. Persistent sensor fluctuations can keep driving its squared-error term; stochastic stationarity is outside this experiment. Observation noise is added only to the synthetic fitting dataset.

What determines the equilibrium?

At equilibrium, v = s = r and b = tanh(m dot r). Each neuron satisfies only

(1-r_i) g_i + (2-r_i) h_i = r_i - (W r)_i - I_i - q_i b

This is one constraint on two conductances: three neurons leave three free directions. The full 13-state Jacobian has three zero eigenvalues. The remaining eigenvalues in the tested configuration are stable, so these are a continuum of transversely stable equilibria, not separated attractor basins or hysteresis. equilibrium() rejects targets for which no positive-conductance solution exists.

Case 1: invariant rule, eta = 0.

c_i = log(h_i) - rho_i log(g_i)
dc_i/dt = 0

Initial conditions fix c. Substituting h = exp(c) g^rho into the equilibrium constraint leaves a strictly monotone scalar root per neuron. The solver uses Brent's method in log coordinates and returns the unique equilibrium on that specified branch. It does not choose a least-norm point or silently use a pseudoinverse. Direct development from the same initial c must reach that root. Changing the starting c changes the adult conductances despite equal targets.

Case 2: local history rule, eta = 1.

dc_i/dt = k eta (r_i-s_i)^2
c_i(T) = c_i(0) + integral_0^T k eta (r_i-s_i(t))^2 dt

The branch now depends on experienced activity. Two rearing protocols start from the same state and use the same rule. One stays at baseline drive; the other reduces drive by 0.25 for the first 200 time units. Both then spend 6000 time units in the same adult environment. Their final voltages agree while their final conductances differ. A control with eta = 0 loses this rearing-history difference.

The history experiment also passes the terminal c from direct integration to the endpoint solver. Agreement is an oracle check of the endpoint representation. It is not an accelerated way of finding c(T). An endpoint solver using the initial c gives the wrong answer in this case.

Gradients

The branch-constrained root is differentiable with c held fixed. Its analytic derivatives include body feedback and are checked against both root finite differences and independent finite differences of the directly integrated developmental endpoint.

recovery_sensitivities() integrates the full tangent equations dS/dt = (df/dy) S + df/dp, with p = (r0, r1, r2, log(k)). Initial tangents come from the adult equilibrium; the instantaneous depletion is then applied to log(g). LSODA integrates the tangent system; Radau provides independent direct trajectories and finite-difference checks. The analytic state Jacobian is also checked away from equilibrium for both regulatory rules.

For eta > 0, holding terminal c fixed omits its derivative with respect to the program and history. A total derivative needs the additional chain-rule term through c(T), or differentiation of the developmental trajectory. The transient sensitivity API therefore rejects eta > 0; that total history derivative is not implemented. There is no valid ordinary inverse of the unconstrained equilibrium Jacobian in either case.

Identifiability experiment

The unknowns are only three activity targets and one shared regulation rate. Wiring, body parameters, rho, the adult branch, intervention strength and observation model are known. This intentionally separates mathematical identifiability from the much larger biological inference problem.

Endpoint and transient ranks use the same number of voltage observations, noise scale, and parameter coordinates. Rank is a local property; three noise seeds are a small recovery check, not a coverage study or proof of global identifiability. Different unknown nuisance parameters or sparse/noisy calcium observations may remove the apparent identifiability.

Verification and numerical gates

The tests cover the conserved quantity, equilibrium multiplicity and transverse stability, developmental agreement, analytic and transient gradients, rearing history with an eta = 0 control, depletion-induced branch changes, body feedback, independent Radau/BDF integration, tolerance refinement, noisy recovery and prediction of an unseen intervention. Invalid and unreachable inputs fail explicitly. The experiment also evaluates 18 named numerical gates and records every result, including failures, in JSON.

These gates check this constructed model. They were developed with the demonstrator and are not preregistered biological falsification criteria. The 6000-unit burn-in is checked against residuals and the known branch root; merely reaching a simulation time limit is not evidence of convergence.

The integrators and fitting algorithm use existing SciPy implementations: solve_ivp, brentq, and least_squares.

Consequences for the CAKE spec

  1. Sensor equality alone is an incomplete developmental contract. Specify initial states, the actual rule, developmental inputs and the reachable branch.
  2. A warm start can change the selected endpoint when branch information is not fixed. Branch identity or sufficient history must be part of a burn-in cache key.
  3. Implicit differentiation requires a nonsingular selected branch and an account of how that branch changes with inferred parameters. Comparing endpoints alone does not validate those derivatives.
  4. Equilibrium fitting cannot identify every compensation timescale. Recovery measurements must enter the parameter-identifiability gate.
  5. The next scientific gate is a measured single-cell regulatory example with observation and intervention models, followed by the small closed-loop circuit. A0 establishes none of the organism-level E1–E5 or F1–F7 acceptance criteria.

Integral-control checks (cake-a0-slice)

src/cake_a0/slice.py implements the general theory in spec L7 and the endpoint solver in spec 06 §5 for any number of regulated quantities n and sensors k: conserved quantities, the slice equation solved by pseudo-transient continuation, loop-gain stability, the projector P = I − G (JG)⁻¹ J, the predicted covariance, and a Lie-bracket involutivity test. src/cake_a0/slice_experiment.py applies it to two synthetic models:

uv run cake-a0-slice --out results/slice

The command exits nonzero if any of its 24 gates fails. A checked run is in evidence/slice/report.md. What it shows:

Claim (spec L7 / 06 §5) Result
Slice endpoints equal directly developed endpoints A0 circuit 2×10⁻¹⁰; two-sensor neuron < 10⁻⁶
Endpoint spread over initial states follows δz* = P δz₀ Relative RMS 2.9% (A0, nonlinearity at sd 0.05), 0.07% (two-sensor)
Cov(z*) = P Σ₀ Pᵀ (+ jitter term) Within the 95% sampling band for 300 samples, with and without set-point jitter
Sensor-controlled combinations carry only set-point jitter J Cov Jᵀ matches Σ_r (5%); without jitter, their variance is 5×10⁻⁴ of the initial
A root of the sensor equations can be unstable Root residual 6×10⁻¹⁷, loop-gain eigenvalues 0.74 and −0.73; development moves away from it, and continuation does not report it
Constant gain: no lasting effect of a transient manipulation Adult difference 9×10⁻¹⁵; conserved quantities drift < 10⁻⁸
Non-integrable gain: lasting effect at identical set points Adult difference 0.024 in log conductance for a strong transient excitation, 1×10⁻⁴ for mild silencing
Gauge G → G C Endpoints identical (6×10⁻¹¹); recovery half-time 0.75 → 0.20
Implicit gradients (set points and initial state) Agree with finite differences of developed endpoints to 1.3×10⁻⁶
det[J; Nᵀ] = ±det(JG)/√det(GᵀG) Holds to rounding over 500 random draws; zero when JG is singular

Limits. The history effect under a non-integrable gain is real but its size depends on how far the manipulation pushes the error trajectory. With mild silencing it is about 10⁻⁴. The size of the effect in neurons is unknown; F7b measures it. Sensors here are steady-state functions of the conductances, not time averages over stochastic activity. A separate scratch test with noisy sensors found no growth of within-type variance with age under a non-integrable gain, so the spec does not assume it.

Proposals T4 and T5 (cake-a0-fluctuation, cake-a0-oneshot)

Prototypes of two proposals in spec/13-proposals.md. Both use the two-sensor neuron above (n = 3 conductances, k = 2 sensors, constant gain); src/cake_a0/population.py evaluates many such neurons at once, optionally coupled through their mean voltage, with analytic Jacobians.

uv run cake-a0-fluctuation --out results/fluctuation   # ~1 min CPU
uv run cake-a0-oneshot --out results/oneshot           # seconds

Each command exits nonzero if any gate fails. Checked runs: evidence/fluctuation/report.md (15 gates), evidence/oneshot/report.md (15 gates).

T4: regulation rates from spontaneous fluctuations. Each neuron receives an Ornstein–Uhlenbeck input current (correlation time τ_mix = 1), and regulation integrates the instantaneous sensor error instead of its average. 256 replicas per condition, with three rate scales c (ε = τ_mix × fastest rate from 0.034 to 0.135).

Claim Result
Stationary mean does not depend on the rate Averaged equilibria identical (< 10⁻⁹)
Stationary spread does depend on it, so G → G C is not a symmetry of the stationary distribution Variance on the slice 3.7× larger for 4× the rate; matches the Lyapunov prediction of the linearized system within 1%
Loop-gain eigenvalues from unperturbed fluctuations (lagged covariances at 5 and 15 τ_mix) Within 6% of the true values at every rate, and within 3 SE
Onsager regression: fluctuations and recovery after a 10% depletion give the same rates Agree within 10%
With observation noise as large as the fluctuations Within 3 SE of the true values; the slowest mode at the lowest rate has SE ≈ 25%
Non-scalar C Same mean, covariances differ by 55 SE
Negative control: a hidden expression cascade, observed only through its output Rate estimate changes by 77% across lag pairs, against 6% with direct observation, so the missing state is detected

T5: regulation as a constrained prior, fitted all at once. 200 neurons of one type, coupled through their mean voltage (strength 0.4), develop from random initial states with set-point jitter. The data are noisy expression (σ = 0.02) and voltage (σ = 0.005) per neuron. The unknowns are θ = (s*, Nᵀμ) plus three adult coordinates per neuron. The components of μ in range(G) affect no endpoint and are not fitted.

Claim Result
The nested fit (slice solve inside every evaluation, 06 §5.3) and the one-shot fit minimize the same objective Estimates agree to 6×10⁻¹², adult states to 2×10⁻¹⁰
One-shot needs far fewer network evaluations 7 against 122, with nested given warm starts; nested also needs 201 fixed-field single-cell evaluations
Truth recovered Within 1.7 SE for every component of θ
The certified optimum is a developmental outcome Loop gain stable per neuron and for the coupled network (min eigenvalue 0.315); development from points 2.4 away on each slice returns to it within 6×10⁻¹²
Regulation-violation map: per-neuron misfit (χ²₄ under H*) flags idiosyncratic neurons AUC 0.65, 0.86, 0.99 for violator spreads 0.03, 0.06, 0.1; no false positives at the 0.999 threshold
Change-of-variables density with |det[J; Nᵀ]| Normalized (1.000 ± 0.004) and reproduces forward-sampled means; without the determinant the integral is 29

Limits. Both are synthetic. In T4 the fluctuations come from an imposed input current, not from simulated network activity, and the observed quantity is the regulated state itself; real expression reporters add maturation and readout kinetics, which the cascade control shows the method must model. The T4 rates are estimated with the slice basis known. In T5 the gain G, initial spread, jitter and noise levels are known, the slice equilibrium is unique, and coupling is mean-field. The evaluation counts compare algorithms on this problem; they are not a measured speedup for organism-scale inference. With small violations (spread comparable to the expression noise) the violation map has little power; its sensitivity is set by the data, not by the fit.