sbml_to_net constant-folds synthesized rate parameters, so sensitivities from the .net are silently zero

[wshlavacek]

sbml_to_net writes a .net in which a synthesized rate-law parameter has been
constant-folded to a number, dropping its symbolic reference to the model
parameter it was built from. The parameter is still declared in the .net, so
nothing downstream notices — but the ODE right-hand side no longer mentions it,
and a forward-sensitivity solve with respect to it returns identically zero,
finite, with no warning.

The same model loaded with Model.from_sbml keeps the reference and gives the
correct answer. So the two loaders disagree about the same model.

Reproducer

One species, one parameter, one reaction. The trigger is a kinetic law with no
explicit compartment factor, which makes bngsim synthesize a volume-normalised
_rateLaw_* parameter.

decay.xml:

<?xml version="1.0" encoding="UTF-8"?>
<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
  <model id="decay">
    <listOfCompartments>
      <compartment id="c" size="1" constant="true"/>
    </listOfCompartments>
    <listOfSpecies>
      <species id="S" compartment="c" initialConcentration="10"
               hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
    </listOfSpecies>
    <listOfParameters>
      <parameter id="k" value="0.3" constant="true"/>
    </listOfParameters>
    <listOfReactions>
      <reaction id="v1" reversible="false">
        <listOfReactants>
          <speciesReference species="S" stoichiometry="1" constant="true"/>
        </listOfReactants>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply><times/><ci>k</ci><ci>S</ci></apply>
          </math>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>

import numpy as np, bngsim
bngsim.sbml_to_net("decay.xml", "decay.net")

for label, m in (("from_sbml", bngsim.Model.from_sbml("decay.xml")),
                 ("from_net ", bngsim.Model.from_net("decay.net"))):
    d = m._core.codegen_data()
    rl = [(p["name"], p["is_const"], p["expression"])
          for p in d["parameters"] if p["name"].startswith("_rateLaw")]
    sim = bngsim.Simulator(m, method="ode", sensitivity_params=["k"],
                           sensitivity_method="staggered")
    s = np.asarray(sim.run(t_span=(0, 5), n_points=20, rtol=1e-10, atol=1e-12).sensitivities)
    print(label, rl, "max|dS/dk| =", np.max(np.abs(s)))

Actual

from_sbml [('_rateLaw_v1', False, 'k / (c / _V0_c)')]  max|dS/dk| = 12.258476
from_net  [('_rateLaw_v1', True,  '0.3')]              max|dS/dk| = 0.0

Expected

Both should agree. For S(t) = S0 e^{-kt}, dS/dk = -S0 t e^{-kt}, whose
magnitude peaks at t = 1/k = 3.33 at 10 * (1/0.3) * e^{-1} = 12.2626.
from_sbml's 12.258476 is that value to the resolution of a 20-point grid.
from_net's 0.0 is wrong, and is returned as a finite array with no error.

Why the zero is the serious part

k is still a declared parameter of the network, so
sensitivity_params=["k"] is accepted, the solve runs, the tensor has the right
shape and np.isfinite(...).all() is True. Nothing in the result says the
model no longer depends on the parameter. Anything built on the .net path —
forward sensitivities, Fisher information, gradient-based fitting, or any
parameter-dependency analysis — gets a confidently wrong answer.

Scope

.net files written by BNG2.pl are unaffected: this is specific to what
sbml_to_net writes. Only kinetic laws that need a synthesized _rateLaw_*
parameter are affected; a law that already carries its compartment factor
(k*S*c) uses k directly and survives the round trip.

Spot-checked against four PEtab benchmark problems, counting _rateLaw_*
parameters that come back is_const=True from the .net:

problem	folded / total
Crauste_CellSystems2017	12 / 12
Raia_CancerResearch2011	2 / 2
Sneyd_PNAS2002	1 / 1
Borghans_BiophysChem1997	1 / 1
Brannmark_JBC2010	0 / 0 (synthesizes none)

For Crauste_CellSystems2017 every one of the 12 fitted rate constants is
folded away, and _rateLaw_v1_v_0 comes back as 0.739907308603256 — exactly
mu_N's value.

Found downstream in bngsim-paper: a parameter-classification pass that asks
"which fitted parameters reach a reaction rate or an initial condition" reads
the .net and now answers none for several SBML problems, where the same pass
on the same inputs answered correctly before. That tool will move to
Model.from_sbml, but the silent-zero behaviour seemed worth reporting on its
own.

Environment

bngsim 0.15.1, build 458bf7d72354
checkout v0.15.1-4-g49dc939
macOS-26.6-arm64-arm-64bit, Python 3.12

Not a prescription, but the two candidate directions

1. Keep the symbolic expression in the .net for a synthesized rate parameter,
   as from_sbml already does, so the two loaders agree.
2. If folding is deliberate, make the loss loud: refuse sensitivity_params
   naming a parameter no rate law or initial condition references, rather than
   returning zeros for it.

(2) is worth having regardless of (1) — a sensitivity request for a parameter
the RHS does not depend on is almost always a mistake, and zero is the least
useful way to say so.