Michaelis-Menten: the compiled RHS has no sFree clamp, so it contradicts both the interpreted RHS and its own emitted Jacobian

[wshlavacek]

Summary

compute_rxn_rate's Michaelis–Menten branch clamps the free substrate:

double sFree = mm_tqssa(Km, E, S).sFree;
if (sFree < 0.0)
    sFree = 0.0;

The compiled RHS does not. _mm_rate_lines in python/bngsim/_codegen.py
emits the root and goes straight to rate = kcat * sFree * E / (Km + sFree) with
no clamp, so both generate_rhs_c and generate_rhs_from_model evaluate the rate
at a negative sFree.

Meanwhile the derivative emitters in the same file — _mm_jacobian_groups and
_mm_v_lines — both guard on sFree > 0.0 && Dmm > 0.0 and return 0, which is
the clamp seen from the derivative side. So does the C++
mm_tqssa_derivatives.

Two consequences, and the second is the one that bites:

1. The interpreted and compiled backends disagree on the same model at the
   same state.
2. The emitted Jacobian contradicts the emitted RHS beside it in the same
   artifact: the Jacobian says the rate is flat (zero derivative) at points
   where the RHS it is differentiating returns a varying nonzero value. CVODE's
   Newton iteration is handed a Jacobian asserting no dependence for a region
   where its own RHS does vary.

Pre-existing; unrelated to #89, which changed how the root is computed, not
whether it is clamped. Found while fixing #89, filed separately because the fix
needs a decision rather than an edit — see below.

Reproduction

tests/data/mm_tqssa.net (kcat=1, Km=50), enzyme E=10, sweeping the
substrate through zero. Interpreted value computed from the shipped formula;
compiled value from generate_rhs_from_model compiled and called through ctypes;
Jacobian from NetworkModel::_dense_analytical_jacobian. On main at d43607e:

emitted RHS clamps sFree? False
emitted Jac guards sFree? True

         S          sFree    interpreted   compiled dP/dt  native dJac dP/dE
     0.001   8.333356e-04   1.666644e-04     1.666644e-04       1.388877e-05
    -1e-09  -8.333333e-10   0.000000e+00    -1.666667e-10       0.000000e+00
    -1e-06  -8.333333e-07   0.000000e+00    -1.666667e-07       0.000000e+00
      -1.0  -8.309926e-01   0.000000e+00    -1.690074e-01       0.000000e+00
      -5.0  -4.105458e+00   0.000000e+00    -8.945417e-01       0.000000e+00

The last column is the contradiction: the analytical Jacobian reports
d(dP/dt)/dE = 0 at every row where the compiled RHS reports a rate that is
plainly a function of E.

sFree < 0 requires S < 0 (for Km, E >= 0, D >= |delta| gives
sFree >= 0), so this is reachable exactly when the integrator probes a negative
concentration. That is not exotic for MM: it is the substrate-exhaustion endgame,
where the state sits against zero and CVODE routinely oversteps it. S = -1e-9
is a realistic probe and already diverges.

Why this is not obviously "add the clamp to codegen"

Both directions are defensible and they are not equivalent, which is why this is
an issue rather than a patch:

(a) Clamp everywhere. Add if (sFree < 0.0) sFree = 0.0; to
_mm_sfree_c_lines, matching the interpreted path, and everything already agrees
with the existing derivative guards. But the clamp puts a kink in the RHS at
sFree = 0 — the rate is flat on one side and varying on the other — and a
non-smooth RHS is exactly what Newton and the error controller handle badly. It
also means the model's dynamics depend on a numerical guard rather than on the
rate law.

(b) Drop the clamp everywhere. Remove it from compute_rxn_rate, and drop
the sFree > 0.0 guard from mm_tqssa_derivatives, _mm_jacobian_groups and
_mm_v_lines, so the RHS and its derivatives are the same smooth function
everywhere the discriminant is real. This is self-consistent and smooth, and it
changes interpreted-path behaviour on negative-concentration probes — the
trajectory is the same in exact arithmetic, but the transient values a solver
sees are not. The S = 0 boundary case has to stay correct: at S = 0 the rate
really is 0 for every kcat/Km, which is what the guards currently also
encode, and test_both_partials_vanish_where_the_rate_does pins.

Whichever way it goes, the fix must move the RHS emitters and the derivative
guards together — the current state is the one combination that is internally
inconsistent.

Scope

• src/model.cpp — compute_rxn_rate's MM branch (has the clamp) and the SSA
  propensity emitter (has its own if (sf < 0) sf = 0).
• python/bngsim/_codegen.py — _mm_sfree_c_lines / _mm_rate_lines (no
  clamp), against _mm_jacobian_groups and _mm_v_lines (guarded).
• include/bngsim/mm_jacobian.hpp — mm_tqssa_derivatives' sFree <= 0 guard.
• python/bngsim/_jax_rhs.py — no clamp either; whichever answer wins applies
  here too.
• A _CODEGEN_VERSION bump if the emitted C changes.

A regression test wants a state with S < 0 asserted through both backends, plus
a check that the emitted Jacobian and the emitted RHS agree about whether the rate
varies there — the property that is violated today.