steady_state(method="newton") can return the saddle on a bistable model and report converged

[wshlavacek]

Summary

On a bistable model, steady_state(method="newton") can return the saddle
between the two branches — an unstable equilibrium no trajectory can rest on —
and report converged=True with a small residual. method="integration" on the
same model and dose returns the correct branch, as does integrating first and
then polishing.

The seed-stability guard does not catch it, and I think the reason is structural
rather than a tuning problem: seed stability asks whether refining the seed
changes the root, which is not the same question as whether the root is
dynamically stable. Near a separatrix the trajectory slows down, so two
successively tighter bursts return nearly the same seed, both polish to the same
saddle, and they agree — the guard is satisfied precisely where it is most
needed.

Reproducer

Gardner toggle (genetic_switch_gardner2000) at alpha_2 = 53.526315789473685.
Self-contained .net:

# Created by BioNetGen 2.9.3
begin parameters
    1 alpha_1     156.25  # Constant
    2 alpha_2     53.526315789473685  # Constant
    3 beta        2.5  # Constant
    4 gamma       1.0  # Constant
    5 k_deg       1.0  # Constant
end parameters
begin functions
    1 syn_R1() alpha_1/(1+(Obs_Tot_R2^beta))
    2 syn_R2() alpha_2/(1+(Obs_Tot_R1^gamma))
end functions
begin species
    1 R1() 0
    2 R2() 0
end species
begin reactions
    1 0 1 syn_R1 #_R1
    2 0 2 syn_R2 #_R2
    3 1 0 k_deg #_R3
    4 2 0 k_deg #_R4
end reactions
begin groups
    1 Obs_Tot_R1           1
    2 Obs_Tot_R2           2
end groups

import bngsim
m = bngsim.Model.from_net("toggle.net")
s = bngsim.Simulator(m, method="ode")

r = s.steady_state(method="newton")
r.converged, r.method_used, r.residual, r.concentrations
# True, 'newton', 2.82e-10, array([28.24521577,  1.83025888])   <- the saddle

m.reset()
s.steady_state(method="integration").concentrations
# array([7.59620790e-03, 5.31227841e+01])                        <- correct

m.reset()                      # integrate to the model's own horizon first
s.run(t_span=(0.0, 20.0), n_points=2)
s.steady_state(method="newton").concentrations
# array([7.59620790e-03, 5.31227841e+01])                        <- correct

It is the saddle, not the other branch

Perturbing the returned state by one part in 1e6 sends it to different
attractors depending on sign:

returned  [28.24521577,  1.83025888]
+1e-6  -> [ 7.5962e-03, 53.1228]     (the branch integration reaches)
-1e-6  -> [144.3734,      0.3682]    (the other stable branch)

The ODEs are

dR1/dt = alpha_1/(1 + R2^beta) - k_deg*R1
dR2/dt = alpha_2/(1 + R1^gamma) - k_deg*R2

so the Jacobian is available in closed form. Its spectrum at the three fixed
points:

fixed point	eigenvalues	max Re
[28.2452, 1.8303] ← returned	+0.40643, -2.40643	+0.40643
[0.0076, 53.1228]	-0.86272, -1.13728	−0.86272
[144.3734, 0.3682]	-0.56558, -1.43442	−0.56558

The returned root has a positive eigenvalue. It is unstable, by a wide margin,
and not a borderline case.

Why the burst does not save it

The trajectory from this model's initial condition [0, 0] heads at the saddle
before committing to a branch:

t	R1	R2	rel. distance to saddle
1	26.31	1.79	6.8e-2
2	27.38	1.86	3.1e-2
5	25.77	1.95	8.8e-2
10	13.16	3.08	5.3e-1
15	0.41	27.58	9.9e-1
30	0.0076	53.12	settled

It comes within 3% of the saddle by t = 2 and stays within 10% for about 4.5
time units. A burst terminating anywhere in that window hands KINSOL a seed a few
percent from the saddle, and the polish resolves to the nearest root. Because the
state is barely moving there, a tighter burst lands in nearly the same place —
so the two polishes agree and seed stability is satisfied.

Suggested fix

The polish already has an analytical Jacobian in hand. Checking its spectrum at
the accepted root costs one eigenvalue solve on a matrix that is already
assembled, and would have rejected this root outright (+0.406 against a
tolerance of 0). Concretely: after the polish converges, reject and continue
integrating if max Re(eig(J)) > 0 (with a small tolerance for the marginal
case, and a cheaper sufficient test — e.g. a few steps of power iteration on
J + sigma*I — if a full eigendecomposition is unattractive on large networks).

A cheaper alternative that does not need the spectrum: perturb the accepted root
and integrate briefly. If the state leaves, it was not a steady state. That is
what I ended up doing outside the library to classify these.

Either way the useful invariant seems to be that seed stability and dynamical
stability are different properties, and only the second is what a caller
asking for a steady state actually wants.

Scope

One dose out of 20 on this model, and BNGsim was correct at the other 19 —
method="integration" and burst-then-polish are both right everywhere. So this
is narrow, but it is silent: converged=True with a 1e-10 residual on a state
the system cannot occupy.

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).
Related but distinct from #74, which is about accumulator species rather than
stability.