Analytic sensitivity RHS for Functional rate laws (df/dp), so forward sensitivity leaves the CVODES difference-quotient path

[wshlavacek]

Summary

generate_sens_rhs_c / generate_sens_from_model (python/bngsim/_codegen.py) return None on the first non-Elementary reaction, so any model with a Functional rate law gets no compiled sensitivity RHS and CVODES falls back to its internal difference quotient. CVodeSensInit1 takes one sens-RHS callback for all columns, so this is all-or-nothing per model.

Every if()-gated rate law is Functional, so the switch-time feature from #48 runs permanently on that fallback — but the problem is generic, not a #48 regression: beta alone is worse than all three switch times on the Lin2021 exemplar.

Across the BNGL-Models corpus:

27 / 65 .net files have >=1 Functional reaction  ->  whole model on the DQ path

including mek_isoform_erk_feedback_miller2026, where 1 Functional law out of 689 reactions puts all 689 on the difference-quotient path.

Why this is worth doing: a spike measured the mechanism

The DQ fallback causes a step-size collapse, not just expensive steps

Lin2021, plain solve vs a one-parameter (beta) sensitivity solve, day 100, 51 output points:

	rtol	wall	n_steps	n_rhs	n_jac	err-test fails
plain	1e-8	0.72 ms	425	577	130	45
sens beta (DQ)	1e-8	5.64 ms	1 555	3 398	290	40
plain	1e-10	1.27 ms	871	1 065	189	58
sens beta (DQ)	1e-10	12 008 ms	5 559 838	12 852 513	1 761 155	1 472

The 9474x wall-clock ratio tracks a 6383x step-count ratio. Error-test failures stay low, i.e. the solver is succeeding at every minuscule step — the signature of an error estimate it can always meet and never grow.

Mechanism. CVODES' internal difference quotient perturbs the parameter by sqrt(rtol)*|p|, so the source term it returns is only ~sqrt(rtol) relatively accurate — 1e-5 at rtol 1e-10. Error control is asked for 1e-10 from a source term carrying 1e-5 of error, so it shrinks h chasing an accuracy the RHS cannot deliver. At rtol 1e-8 the DQ delivers ~1e-4 and the solver tolerates it. The cliff should therefore appear exactly when the requested tolerance passes below the DQ's own accuracy, which is what the table shows.

Controlled A/B: an exact source term removes the cliff

An all-Elementary synthetic SEIR (exponential growth from a single seed in a population of 2e7 over 100 days — the stiffness property blamed for the cliff) does get an analytic sens RHS. Forcing generate_sens_from_model / generate_sens_rhs_c to return None reproduces the Functional situation exactly (no bngsim_codegen_sens_rhs symbol => CVODES internal DQ) on identical dynamics:

	rtol	wall	n_steps
plain	1e-8	0.07 ms	254
analytic sens RHS	1e-8	0.15 ms	272
forced DQ sens RHS	1e-8	0.15 ms	272
plain	1e-10	0.13 ms	492
analytic sens RHS	1e-10	0.29 ms	649
forced DQ sens RHS	1e-10	26.03 ms	76 935

DQ / analytic = 90.8x at rtol 1e-10, 1.01x at rtol 1e-8. Both gradients agree to rel 1.4e-8 — the DQ reaches the right answer, it just pays ~91x for it. With the analytic source term there is no cliff at all: 649 steps against the plain solve's 492.

And it flips the comparison against finite differences

Removing the cliff is necessary but not sufficient — forward sensitivity also has to beat central FD's 2N plain solves. Each analytic sens column is one linear variational solve (no Newton iteration), so it should cost ~1 plain-solve equivalent rather than 2. Measured on the same synthetic model:

rtol	N	analytic sens	per param	central FD (2N)	sens/FD
1e-8	1	0.123 ms	0.057 ms	0.131 ms	0.94x
1e-8	2	0.163 ms	0.049 ms	0.263 ms	0.62x
1e-8	3	0.215 ms	0.050 ms	0.394 ms	0.55x
1e-10	3	0.431 ms	0.106 ms	0.676 ms	0.64x

An analytic sens column costs 0.75–0.95 plain-solve equivalents against central FD's 2, so forward sensitivity wins from N=2 and the margin widens with N.

Lin2021 on the DQ path costs ~4 plain-solve equivalents per column today (2.6 ms/param against a 0.66 ms plain solve), which is why a tuned central FD is currently 2.3–5.1x cheaper per gradient there. Landing Lin2021 in the synthetic model's 0.75–1 band would be a ~4x per-gradient improvement at rtol 1e-8 and would turn a 2.3x loss to FD into a ~1.6–2x win.

What the spike does not establish: the magnitude on Lin2021 (synthetic cliff 91x vs Lin2021's 9474x — same mechanism, different severity), and whether the per-column constant transfers from a 4-species/3-reaction model to a 26-species/61-reaction one. The ratio argument is structural; the constant is not measured on a real Functional model.

Scope

The maths is simpler than the analytical Jacobian, not harder

For a Functional reaction, rate_r = func_r(obs, p, t) * prod(R) (.net, apply_species_factor=true) or rate_r = func_r(...) (SBML, per-species). prod(R) is species-only, so

df_i/dp  =  sum_r  stat_r * netstoich_ir * (dfunc_r/dp) * prod(R)_r

There is no product-rule term and no chain rule through observable groups. The dfunc/dobs_k -> dobs_k/dx_j machinery in FunctionalJacobianData::ObservableTerm exists only to turn observable-derivatives into species columns; for df/dp the parameter is already the differentiation variable.

Most of the machinery already exists

piece	where	reuse
_inline_functions	_jacobian.py	as-is
_inline_derived_param_refs	_codegen.py	as-is — flattens sigma = t0+t_delta to primaries before diff
_exprtk_to_sympy / _build_local_dict	_jacobian.py	as-is
_is_emittable, sympy_to_c, resolve_symbol	_jacobian.py / _codegen.py	as-is — generate_jacobian_from_model's c_ref already maps both params (p[i]) and observables (obs[j])
_compute_derived_param_jacobian	_codegen.py	as-is — dp_d/dprimary; _emit_sens_rhs_body's derived_terms already threads it into the switch (iP)
_net_affected	generate_jacobian_from_model	as-is — (row i, coeff = stat*netstoich) is exactly the df/dp scatter
_emit_observable_lines / _emit_function_lines	_codegen.py	as-is — the obs[]/func[] recomputation the derivatives read
bngsim_output_sens_dfdp	generate_output_sens_from_model (#198)	closest analogue: already emits switch (iP) { dfdp_out[m] += ... } over per-function dfunc_m/dp_k with the #15 derived chain folded in, with an (iP, t, y, p, obs, func, out) signature. It is bngsim_dfdp one level up.

Measured on Lin2021: differentiating every Functional rate law w.r.t. every parameter it references takes 23 ms, yields 19 non-zero (law, param) partials, and emits valid C with zero emit failures.

The correctness hazard this introduces

sympy.diff of an if()/Piecewise w.r.t. a condition-only parameter returns a clean 0 (no Dirac delta, so _is_emittable never rejects it). Within a smooth segment that is correct, and #48's jump s+ = s- + (f- - f+)*dt*/dp supplies the whole contribution — so the switch jump stays the mechanism and needs no codegen change.

But that same 0 is a new way to be silently wrong. compute_switch_time_sens only recognises time / counter-clock thresholds. For a state-dependent condition (if(X >= thresh, ...) with thresh fitted) an analytic 0 with no jump is plainly wrong, whereas today's DQ probe at least perturbs the crossing and returns something. Switching to analytic df/dp must therefore be gated: decline the analytic sens RHS (per model — it is all-or-nothing) whenever a requested parameter occurs in a condition that is not a recognised clock threshold.

This is where the correctness risk concentrates and where the test effort goes.

Related: #198's _UNSUPPORTED_EXPR_CONSTRUCTS bans if(), comparisons and logicals outright ("fail loudly"). df/dp has to lift those three bans — every if()-gated law is Functional, which is the entire population — while keeping abs/min/max/floor/ceil/round banned. Note _is_emittable will not catch Min/Max, since both are in _SYMPY_FUNC_TO_EXPRTK; the construct gate has to do it.

What still has to bail

• Michaelis–Menten reactions — out of scope for "Functional" but also on the DQ path today (the gate is != "elementary"). Closed-form and easy; ~1 day if bundled.
• Table functions (tfun_*) — deliberately value-only, never differentiated.
• rateOf — needs the live dx/dt buffer, unavailable in a sens-RHS call.
• abs/min/max/floor/ceil/round — kinks with no jump machinery to compensate.
• Hidden state — any free symbol that is not an observable, parameter or time (e.g. an SBML rate-rule target). Already caught by the existing free-symbol check.
• Cycles / inline depth > 64 — existing _inline_functions guard.
• State-dependent switch conditions — see the hazard above.
• Variable-volume parameters (#171) — the Jacobian has volume_terms for d/dV_live; df/dV needs the analogue. Decline a fitted volume parameter initially.
• Genome-scale — the fix(codegen): bound the sensitivity df/dp derivation with a #95-style budget (#90) #95/#187 derivation budget applies, and an N-parameter df/dp multiplies the SymPy work by N. Needs the same budget guard.

Estimate: ~2 working weeks (7–15 days)

piece	days	notes
df/dp emitter (per-species + per-observable)	2–3	assembly of the table above; bngsim_dfdp gains obs[]/func[], and _analyze_output_sens's pruning must add reaction-bound _rateLawN as roots
Functional bngsim_jac_vec	2–3	mechanical transform of generate_jacobian_from_model's balanced {...} groups into Jv_out[i] += coeff*val*v[j]; touches chunking/sharding and the dense/sparse split
switch-condition gating + #48 pinning interaction	2	mostly test design; the correctness risk lives here
MM df/dp + jac_vec	1	optional
validation: FD cross-checks on the 27 slow-path files, golden + #212/#48 regression	2–3
integration / CI / perf measurement	1

Low end 7 days if the Jacobian groups transform cleanly and no C++ layout change is needed; high end 15 if the sens-RHS user-data change ripples through the sharded-compile path and the SBML per-species path needs its own gating.

Two design points worth settling early:

1. bngsim_dfdp's current signature is (iP, t, y, p, dfdp_out) — it has no obs[]/func[], which the Functional derivatives read. Either extend the signature and recompute them the way generate_jacobian_from_model does, or pass them through CodegenSensUserData. Whichever is chosen, check whether it forces a _CODEGEN_VERSION bump and whether it ripples into compile_rhs's sharded path.
2. bngsim_jac_vec is Elementary-only. Emitting a Functional J*v mirrors the existing Jacobian groups; the alternative (reuse bngsim_codegen_jac/_jac_sparse and matvec) is cheaper to write but needs an n*n or CSC buffer live in the sens RHS.

Related

• Gradient-based fitting of switch-time parameters: forward sensitivity drops the discontinuity-crossing term at if(t>=σ) steps #48 — switch-time forward sensitivity (merged); its write-up flagged this as the orthogonal follow-up.
• Forward sensitivity w.r.t. a switch time encoded as an SBML event (event-time sensitivity dt*/dp != 0) #49 — the same switch encoded as an SBML event. Broader correctness reach (24 corpus models currently refused outright) but it would extend the feature over the same DQ path; the spike above is the argument for doing this issue first.
• fix: accept ExprTk's symbolic logical operators in the symbolic differentiator #53 — accepted ExprTk's symbolic logical operators in the symbolic differentiator, which gives Lin2021 (and 2 other corpus files) an analytical Jacobian. Prerequisite for df/dp coverage of the Lambda/P laws.
• PSA: expose per-reaction scaling diagnostics and a cost/precision decision signal #15, Steady-state FD ∂f/∂p uses an absolute step floor, so a small rate constant is probed at 1500% (15.9x error on before_bunching) #76, fix(codegen): bound the sensitivity df/dp derivation with a #95-style budget (#90) #95, #151, #162, #171, #187, #198 — the machinery being reused.

[wshlavacek]

Staged into four children (2026-07-26)

Scoped against the in-repo 585-model corpus after #56/#51/#52/#54 merged. The
evidence in this issue stands — the children re-stage the work, they do not
replace the argument here. Ordered so that every behavior change is small, late,
and independently revertable:

• #55a — sens RHS plumbing: give bngsim_dfdp observable/function context (no behavior change) #65 — sens RHS plumbing: give bngsim_dfdp observable/function context (no behavior change)
• #55b — analytic df/dp for Functional rate laws (emitted, gated off) #66 — analytic df/dp for Functional rate laws (emitted, gated off; solver-free FD oracle)
• #55c — enable the analytic sens RHS for condition-free Functional models (140/585) #67 — enable the analytic sens RHS for condition-free Functional models (140/585 — first behavior change)
• #55d — extend the switch-condition guard so condition-bearing Functional laws can use the analytic sens RHS #68 — extend the switch-condition guard so condition-bearing Functional laws can use it (28 models, highest risk, done alone)

Three measurements from the scoping pass that re-shape this issue's own
component breakdown:

1. 185/585 models (32%) get no analytic sens RHS, blocked by 1.1% of
   reactions (1736 functional / 11 mm out of ~161k). The all-or-nothing
   CVodeSensInit1 callback is the amplifier. Same phenomenon as this issue's
   "27/65" from the external corpus, measured on what CI and the benchmarks
   actually run.

2. Of those 185, 140 (76%) are CLEAN — Functional laws with no condition and
   no banned construct. Only 28 carry the switch-time correctness hazard; 17 are
   permanently out of scope. Hence staging by risk class rather than by
   component: #55c — enable the analytic sens RHS for condition-free Functional models (140/585) #67 unblocks 140 models with none of the risk, and #55d — extend the switch-condition guard so condition-bearing Functional laws can use the analytic sens RHS #68 handles the
   28 alone.

3. The Functional analytical Jacobian already exists and works —
   generate_jacobian_from_model emitted one for 20/20 sampled CLEAN models. So
   J*v is a matvec over an existing correct matrix, and this issue's 2–3 day
   "Functional bngsim_jac_vec" line leaves the critical path and becomes a pure
   optimization. df/dp is the only genuinely new derivative work.

Michaelis–Menten is deferred: 11 reactions across 2 models, one also
Functional-blocked, so it unblocks exactly one additional model. Done
anyway — that measurement is about the corpus, not the feature. MM(kcat, Km)
is a first-class BNGL rate law and a modeller who writes one should keep the
analytic gradient.

The highest-risk item is not the math — it is duplicating a guard. #52's merged
event-time guard and #68's needed check answer the same question, and two
spellings will drift. See #68 for why that is not hypothetical.

---

All four stages merged, plus Michaelis–Menten (2026-07-28)

stage	merged	corpus with an analytic sens RHS
#65	2c981fa	— (plumbing, byte-identical output)
#66	4ab9e72	— (emitted, gated off)
#67	7ebfdfe (#84)	398 → 535
#68	1dbbba1 (#86)	535 → 543
MM	eeecc4b (#88)	543 → 544

This issue's performance claim held. It predicted 0.75–0.95 plain-solve
equivalents per analytic column against central FD's 2; #67 measured 0.59–1.04
at k=16 against the difference quotient's 9.0–36.7, with DQ needing a median
295× more steps and, on 2 of 6 benchmarked models, never converging at all —
exhausting mxstep and returning a gradient of exactly zero.

Two of the three scoping measurements above held exactly. The third did not:
§3's "J*v is a matvec over an existing correct matrix" was true about the
math and wrong about the cost — 133 of 139 CLEAN models take the dense
Jacobian branch, so a matvec would have put an n×n buffer (max 267 KiB) inside a
CVODES callback and forced a CodegenSensUserData widening. #67 emitted a fused
bngsim_jac_vec instead, reusing the Jacobian's reconstruction through an
add(col, row, value_c, prefix) scatter callable: zero buffer, O(nnz), no ABI
change, one derivation for both consumers. The original 2–3 day estimate for
that line was right.

#68 shipped a stricter rule than it specified. "Every requested sensitivity
parameter's occurrence in a condition must be a clock threshold" would admit
if(Xtot > 0, k*A, 0) — no parameter in the condition — but that crossing still
moves with every parameter through the trajectory, which is #52's neuron
hazard in rate-law clothing, and #52 refuses it for events. Two guards answering
one question opposite ways is exactly the drift #68 was filed against. Measured
first: CVODES' DQ returns exactly 0 for a state-dependent crossing too, so
admitting them was never a regression — just an extension into a class we know
is incomplete. 8 models instead of 14.

What is left

• The sensitivity df/dp path has no #95 derivation budget, so a genome-scale Functional model can hang the build instead of declining #90 — no fix(codegen): bound the sensitivity df/dp derivation with a #95-style budget (#90) #95 derivation budget on the df/dp path. The last unaddressed
  item from this issue's own "What still has to bail" list: "an N-parameter
  df/dp multiplies the SymPy work by N. Needs the same budget guard."
  _codegen.py has zero deadline references. Never tripped on the corpus
  (2.6 s for all 185 models, 0.31 s worst), but the failure mode is a build that
  hangs rather than declines. This issue should stay open on it.

• tQSSA Michaelis-Menten: sFree loses ~2 digits per decade of |delta|/sqrt(4*Km*S), in the RHS as well as the Jacobian #89 — adjacent, and worth knowing before trusting an MM gradient. The tQSSA
  sFree = ½(δ + D) cancels catastrophically once |δ| ≫ √(4·Km·S), ~2 digits
  per decade. That is in the RHS, not just the derivative, so in that regime
  this issue's premise — that the analytic path beats the DQ — is a comparison
  between two wrong answers. No corpus MM model is anywhere near it (both sit at
  ratio ~1), so nothing measured here changes; but "MM has an analytic
  sensitivity RHS" is a weaker claim than it sounds until tQSSA Michaelis-Menten: sFree loses ~2 digits per decade of |delta|/sqrt(4*Km*S), in the RHS as well as the Jacobian #89 lands.

Permanently out of scope, unchanged: tfun, rateOf,
abs/min/max/floor/ceil/round. Already guarded: hidden state, cycles
and inline depth. Declined with a reason: variable-volume parameters (#171),
state-dependent switch conditions (#68), and any Functional law whose derivative
is not representable in C.

Next, per this issue's own argument for doing it first: #49 — the same switch
encoded as an SBML event.