Structural identifiability as a pre-fit gate: decide symbolically whether a parameter can EVER be identified from these observables, instead of inferring it numerically from a flat profile

[wshlavacek]

ProfileLikelihoodAlgorithm._classify reports structural as a verdict — a profile explored with Δχ² < 1e-3 is called structurally unidentifiable. That verdict is inferred numerically, a posteriori, and expensively: it costs a full polish plus a two-directional profile walk with a reduced-dimension re-optimization at every grid point, and it can only ever be evidence, never proof. A flat profile over the explored range is consistent with structural unidentifiability, but also with a very shallow-but-finite curvature the walk did not reach.

The question it is approximating has an exact, data-free, a priori answer.

The gap

Structural identifiability asks: given only the model

  ẋ = f(x, p),    y = h(x, p),    x(0) = x₀(p)

can p be determined from perfect, noise-free, continuous observation of y? It needs no data, no fit, and no simulation — it is a property of the model-and-observable pair alone, decidable symbolically in seconds for models of PyBNF's size.

PyBNF today has no such check. A user can spend hours of CPU on a multistart fit, then a profile-likelihood run, to learn that two rate constants only ever appear as a product — a fact derivable before the first solve.

Why PyBNF is well positioned

• sympy is already a dependency, used across petab/formula.py, analytical_model.py, and measurement/base.py.
• The observables are already symbolic. compile_petab_formula holds each observable as a sympy expression before lambdifying (ADR-0035/0036). h(x, p) needs no new parsing.
• The RHS is mechanically recoverable. A .net expansion is stoichiometry plus rate laws; bngsim_model/net_model.py and bngsim_model/expressions.py already handle both. ẋ = S · v(x, p) is assembly, not inference.
• The verdict has a home. profile_likelihood_summary.txt already has a classification column. A structural pre-check writes into the same vocabulary and disambiguates it: structural becomes a proven claim rather than an inference, and a parameter that is structurally identifiable but flat in the profile is correctly labelled practically unidentifiable — which is a completely different instruction to the user (get better/more data vs. change the model or the observables).

Method options

1. Observability-rank / Lie-derivative condition (Hermann–Krener). Augment the state with ṗ = 0, build successive Lie derivatives of each observable along f, and test the rank of the Jacobian of [h, L_f h, L_f² h, …] w.r.t. (x, p). Implementable natively in sympy, no new dependency. Rank is evaluated at a random rational point to dodge symbolic-rank blowup — which makes it probabilistic (correct with probability 1, not certainty), the standard and accepted tradeoff.
2. Differential algebra (DAISY / input–output equations, Gröbner bases). Gives global identifiability, not just local. Much more expensive and prone to blowing up on models with more than a handful of states.
3. Shell out to StructuralIdentifiability.jl or SIAN. Best-in-class, but imports a Julia/Maple dependency into a Python tool — a hard sell for a default-on pre-check.

Recommendation: (1), native sympy. It answers the question that matters here — local structural identifiability, which is what a profile-likelihood verdict is about anyway — with no new dependency, and it degrades gracefully (a rank test that gets too large can time out and report "undetermined" rather than block the fit).

Worth noting that (1) also distinguishes the two failure modes usefully: the null space of the rank-deficient Jacobian names the offending combination (e.g. "only the product k1·k2 is identifiable"), which is far more actionable than a flat curve.

Scope pins

• ODE models only. A network-free / NFsim model (bngsim_model/nf_model.py) has no finite ODE system to differentiate, so there is nothing to test. Refuse by name, in the style of noise_profiling's refusals, rather than silently skipping.
• Network expansion must terminate. A model whose .net expansion is infinite or impractically large is out for the same reason.
• Initial conditions are part of the question. Identifiability depends on x₀, including which initials are themselves fitted (the seed/condition-target machinery, ADR-0095/0096/0101). A parameter identifiable from a generic x₀ can be unidentifiable from the specific x₀ the experiment uses — the check must use the fit's initial conditions, not a generic point, or it will give confidently wrong answers.
• Observables must be the fit's observables, including transforms and normalization — a scale-normalized column destroys the amplitude information that would otherwise identify a parameter, and the check must see that.

Proposed surface

Cheap enough to run automatically as a pre-fit gate with a warning (not a hard failure — a user may knowingly fit an unidentifiable model), plus a standalone entry point for asking the question directly. job_type = check (#563's model checker) is a plausible host, since it is already the "interrogate this model without fitting it" job.

Evaluation

1. Reproduce known textbook results — the classic two-compartment and Michaelis–Menten unidentifiable cases must come out unidentifiable, with the correct offending combination named.
2. Run it across the Grein subset-I corpus and compare its verdicts against what profile_likelihood reports numerically. Disagreements are the interesting output either way: a parameter called structural by the profile but structurally identifiable by the rank test is a case where more data would actually help, and the user is currently being told the opposite.
3. Measure runtime as a function of state count to find where the rank test stops being seconds-cheap, and set the timeout there.