test(immune): five population-dynamics models with closed forms, and an observable count-predicate defect pinned

[akutuva21]

Immune population-dynamics models with closed forms, and an observable count-predicate defect pinned.

No C++ in this diff. git diff 6889fba -- cpp/ is empty. The one defect I found is already fixed by @correctness in PR #59 (commit 1ddef49), on the identical line, and I withdrew my own fix rather than duplicate it — details below.

What this adds

Five small immune models (models/immune/, all <200 species) plus a runnable harness. Every model is answered by a number derived independently of the model text, or by an exact conservation identity. Nothing here is a plausibility check.

model	the claim it pins	measured
clonal_expansion	logistic closed form K/(1+((K-E0)/E0)e^(-rt)) for dE/dt = r*E*(K-E)/K	max rel deviation 1.681e-09 over 9 steps
	clone approaches declared K=1000, monotone, never overshoots	E(8)=931.738460
	seeded SSA concentrates near the ODE mean	932 vs 931.7385 (0.03%)
cytotoxic_killing	k(E)=Vmax*E/(Kh+E) at E==Kh is exactly Vmax/2, so the decay is T0*exp(-Vmax*t/2) with no free parameter	max rel deviation 1.163e-10 over 11 steps
	saturation signature, swept over E ∈ {20, 200, 2000}	see test_kill_flux_saturates_in_the_effector_count
	validity condition: E is unconsumed, else the closed form does not apply	E == 200 at every step, asserted
exhaustion_switch	stimulus is exactly linear S(t)=S0+a*t	max deviation 0.0
	reaches THETA=400 at exactly t* = (THETA-S0)/a	S==THETA at t=5.0, analytic t*=5.0
	switch flux is kmem*Eff*Hill(1,THETA,n,S)	max rel dev 8.296e-03 (= my trapezoid differencing error at dt=0.05, not a rate error)
	switch-time tolerance is derived, not chosen: the Hill turns on over THETA*(10^(1/n)-1)/a	band ±2.22 d at n=8; measured switch t=5.025
antigen_binding	free+bound antigen conserved at every sampled step, not just t_end	ODE max drift 5.0e-11 over 601 steps; SSA drift exactly 0.0 over 61 steps (counts are integers)
	the split kon/koff rates, pinned by the physical root of x²-(Ag0+R0+Kd)x+Ag0*R0=0	bound=184.112765606 vs analytic 184.112765606, rel 5.897e-14
	.net carries a distinct reverse reaction consuming/reproducing the same species	asserted by token, not by substring

The binding case is the one that matters most for the checklist's §4 trap: conservation alone cannot catch a mis-split rate, because a rule that dropped or doubled koff still conserves the pool. Only the equilibrium root distinguishes them, which is why both are asserted.

Defect found: count predicate on a Molecules observable

OdeIntegrator::compileGroups() applied the predicate to the number of pattern→species embeddings, which is 1 for every single-node pattern.

Smallest reproducer — models/immune/threshold_observable.bngl:

  Molecules Eff  Eff()
  Molecules Gt50 Eff()>50        with Eff() seeded at 100

  BNG2 2.9.3 -> Gt50 = 1.000000000000e+02   (predicate inert; equals Plain)
  BNG3       -> Gt50 = 0.0                 (1 > 50 compared against 1)

Same in SSA with seed 42. This is the "runs and is wrong" shape: no error, a plausible column of zeros.

Oracle, and why the defect is bounded to Molecules. BNG2 2.9.3 at /Users/akutuva/Documents/BioNetGen/bionetgen/bionetgen/bng2/BNG2.pl. Perl2/Observable.pm inspects $patt->Quantifier only inside the Type eq "Species" branch; the Molecules branch never evaluates it. I verified Species-typed observables already agree with BNG2 term for term (SEff/SGt0/SGt1/SCh all to 1e-9), which is what confines the defect to Molecules.

Why the fix is not in this PR

@correctness independently found the same defect and landed 1ddef49 in PR #59, on the same line, gating the filter on the observable type. I had a fail-closed refusal written and ready. I withdrew it:

• Their landing is BNG2-compatible term for term.
• Mine would refuse. I confirmed afterwards that "inert" is BNG2's defined behaviour, not an unexercised branch — Perl2/Observable.pm simply never inspects the quantifier on that path. An unexercised branch is a weaker foundation for refusal than for parity, and the oracle is what settles it.
• perfEngine cleared compileGroups() to me by name, so the collision was avoidable but still real.

My tests therefore assert the parity outcome, and I verified them against a binary that actually carries their fix rather than only against my own branch:

binary WITHOUT the observable fix (6889fba):  3 failed, 19 passed
binary WITH  correctness 1ddef49 applied:     22 passed in 0.88 s
models/immune/verify_immune.py, fixed binary: all checks passed

The second number is the one that matters. A duplicate fix that passes only against its own branch is precisely what breaks when two lanes touch one file, so the models were run through the merged semantics, not each half separately.

Reported, not fixed, not in this diff

BNG3 lowercases a Hill call on emission, and the emitted text does not resolve on the way back in:

  source   Effector() -> Memory()  kmem*Hill(1,THETA,n,Stimulus())
  emitted  _rateLaw1() kmem*hill(1,THETA,n,Stimulus())       <- lowercase
  re-read  "action execution failed: ... hill: unresolved expression symbol: hill"

The trajectory is correct through the normal path (flux tracked to 8.3e-3 against the closed form), so this is a portability defect in emitted text, not a wrong number — which is why it survived. It is exactly docs/DEVELOPMENT_CHECKLIST.md §4: a claim about emitted text that is not a claim about behaviour, except here the behaviour is the opposite of what would be claimed.

Sent to @sciPkPd (owns compile()'s Sat/Hill block) and @swarmCompiler (owns Expression.*) with a two-file reproducer. I have not located the lowercasing site and am not claiming which file emits it. Expression.cpp:333 accepts both spellings at evaluation time, which points at the emission side, but that is a lead and I did not confirm it. I also have no measurement on whether Sat/MM share the problem and claim nothing either way.

What I could not assess here

• PNGL corpus parity. I did not run the validation corpus; no claim about it.
• Sat/MM round-trip. Not exercised by my models, so untested rather than clean.
• NFsim. No model here uses it. My base binary is behind HEAD by exactly one commit (75b22a7, cpp/nfsim/NFinput/NFinput_fromCompiled.cpp only) — git log --since='2026-09-29 15:08:00' --format=%h -- cpp/ returns that alone — which this lane never reaches.
• Multi-compartment immune models. All five are single-compartment by design, to keep the closed forms exact.

Verification

BNG_CPP=build/cpp/bng_cpp python3 -m pytest tests/python/test_immune_models.py -q
  -> 22 passed in 0.88 s
BNG_CPP=build/cpp/bng_cpp python3 models/immune/verify_immune.py
  -> all checks passed
ruff: clean.  black --check --target-version py39: clean.

These tests drive the bng_cpp CLI and never import bionetgen, so neither conftest-guard variant (#53/#56) nor the importorskip hazard affects them; no UNDER TEST pkg/ext line applies because no extension is loaded.

Metric class, declared per Main's rule: every number above is a residual against a closed form or an exact conservation identity. None is a duration, none is load-sensitive, so no benchmark slot is claimed or needed. No timing is reported.

No git stash was used at any point. All writes used absolute paths. git status --porcelain is empty; zero unmerged files; the three pre-existing stash entries were not touched.

Merge note: this branch adds no lines to tests/cpp/CMakeLists.txt, so it does not add to the eight-way contention in OdeIntegrator.cpp. Rebase onto current main is expected but trivial.

🤖 Generated with Claude Code