steady_state(): no way to exclude write-only accumulator species from the convergence test

[wshlavacek]

Summary

Simulator.steady_state() has no way to restrict the convergence test to a
subset of species. When a model carries a write-only accumulator species — a
"degraded"/"produced" pool that some reaction produces and no reaction consumes —
that species has a constant non-zero derivative forever, so ||f(y)||_2/n can
never reach tol and the solve reports failure even though every other species
has settled.

AMICI has this knob (Model.set_steadystate_mask). BNGsim currently does not, so
there is no way to express "solve f(y) = 0 on this subspace", and no symmetric
way to compare the two engines on such a model.

Concrete case

beta_catenin_destruction_complex_barua2013 (BNGL-Models), generated with the
model's own bound max_stoich=>{APC=>1,AXIN=>1,bCat=>1} → 409 species, 2737
reactions. Exactly four species are pure sinks — produced by some reaction,
consumed by none:

bCat(ARM34,ARM59,s33s37~U,s45~U,ss~d)
bCat(ARM34,ARM59,s33s37~U,s45~P,ss~d)
bCat(ARM34,ARM59,s33s37~P,s45~P,ss~d)
bCat(ARM34,ARM59,s33s37~P,s45~U,ss~d)

Integrating decades past the settling time leaves the residual pinned while the
state grows linearly:

| max_time | converged | residual | max|y| |
|---|---|---|---|
| 2.5e6 | False | 7.494e-3 | 7.484e6 |
| 2.5e7 | False | 7.494e-3 | 7.484e7 |
| 2.5e8 | False | 7.494e-3 | 7.484e8 |

The residual is unchanged across two decades — a constant non-zero derivative,
not slow convergence. The other 405 species are settled to 1e-10 by the first
horizon.

import bngsim
m = bngsim.Model.from_net("B05_Barua_2013.net")
s = bngsim.Simulator(m, method="ode")
r = s.steady_state(method="integration", max_time=2.5e8)
r.converged        # False
r.residual         # 0.007494031774924183

Proposal

A steady_state_mask on the model (or a mask=/species= argument to
steady_state()) selecting which species enter the convergence norm. Everything
still integrates; only the test is restricted:

sim.steady_state(method="newton", mask=~model.is_pure_sink())

Two things that would make it usable without the caller hand-listing species:

1. Detect pure sinks from the network. A species that never appears as a
   reactant is structurally identifiable at load time — no user annotation
   needed. Exposing that as e.g. Model.pure_sink_species() would let the mask
   default to something sensible, and would be independently useful for
   diagnostics.
2. Say so in the failure. Today a model like this returns
   converged=False with no indication that the cause is structural rather than
   numerical. If a pure sink with non-zero flux is present, the result (or the
   raised error) could name it — that alone would have saved a fair amount of
   time here.

Why it matters beyond one model

Accumulator species are a common BNGL idiom for counting cumulative flux
(degradation, secretion, transcription events). Any such model is currently
outside steady_state()'s reach, even when its biologically meaningful subsystem
has a perfectly good steady state. The workaround — rewriting the accumulator's
producing reactions to have no product — changes the model file to work around a
solver limitation, and is only provably safe because the species is a pure
sink, which is exactly the property the library could check itself.

Context

Found while measuring BNGsim's KINSOL solve against AMICI's damped Newton over
dose scans of six published models for the BNGsim paper (Section 5.7,
Table 15). It is the one model of the six whose steady-state question is
currently ill-posed for BNGsim; the other five are measured and agree.