PEtab export reads SBML initial values from the wrong place when an initial assignment sets them

[wshlavacek]

SBML lets a model set a parameter's initial value in two places: the value
attribute on the parameter element, and an initialAssignment. When both are
present the initial assignment wins, and the value attribute is often left as a
placeholder. The same is true for a species, where the placeholder is the
initialAmount or initialConcentration attribute.

PyBNF's SBML reader, pybnf/petab/_sbml.py, only ever reads the attributes. It
has no branch for listOfInitialAssignments. So when a model uses an initial
assignment, the reader reports a number the model does not actually start from,
and says nothing about it.

That number is used. pybnf/petab/export.py:1421 (_numeric_nominal) hands it to
pybnf/petab/conditions.py:183-195, which turns a relative condition, such as
multiply this parameter by two, into an absolute number written to
conditions.tsv. The exported problem then carries a wrong number with no
diagnostic.

Repro. Write this file as m.xml:

<?xml version="1.0" encoding="UTF-8"?>
<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
  <model id="m">
    <listOfCompartments><compartment id="c" size="1" constant="true"/></listOfCompartments>
    <listOfParameters>
      <parameter id="koff" value="1" constant="true"/>
      <parameter id="Kd" value="2" constant="true"/>
      <parameter id="kon" value="0" constant="true"/>
    </listOfParameters>
    <listOfInitialAssignments>
      <initialAssignment symbol="kon">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply><divide/><ci>koff</ci><ci>Kd</ci></apply>
        </math>
      </initialAssignment>
    </listOfInitialAssignments>
  </model>
</sbml>

Then run:

from pathlib import Path
from pybnf.petab.export import _read_model, _numeric_nominal
from pybnf.petab.conditions import mutation_target_value
v = _read_model(Path('m.xml').read_text(), Path('m.xml'), 'sbml')
print(_numeric_nominal(v, 'kon'))
print(mutation_target_value('*', 2, nominal=_numeric_nominal(v, 'kon')))

It prints 0.0 and then 0. The model's effective kon is koff/Kd, which is 0.5,
so doubling it should export 1.0.

BioNetGen models are not affected, and only by accident. A BioNetGen expression
is read as a string, so float() on it raises and the relative condition is
refused at conditions.py:183-188 rather than computed from a wrong number.

How much this bites is not yet known, which is why this is an issue rather than
a change folded into other work. Two SBML models in the repository contain
initial assignments: tests/petab_fixtures/boehm_v2, where STAT5A and STAT5B
carry initialConcentration="1" while the initial assignments set them from
207.6*ratio, and tests/full_tests/T2-ade-abcd, which assigns species A. Neither
is currently the target of a relative condition, so no committed artifact is
wrong today. Someone needs to work out how far this reaches, including the
species path, before choosing between the two plausible fixes: read
listOfInitialAssignments and evaluate an assignment when it refers to nothing
else, or refuse the relative condition with a clear error, which is what the
BioNetGen path already does in effect.

For context, the petab library is making the same call for BioNetGen models in
PEtab-dev/libpetab-python#517: a value the model file alone cannot settle is not
reported at all, because a PEtab parameter table may override or estimate what
it depends on. Our SBML reader should end up following the same rule.

[wshlavacek]

Scoped. Summary first: nothing in this repository is affected today, the exposure is in models users bring, and the useful part of the fix is refusing a value rather than computing one.

Who actually uses the values the reader collects

There are two live consumers and one dead one.

pybnf/config.py:2691 puts the reader's values into constants, which the measurement layer binds when an observable or noise formula names a parameter (pybnf/measurement/base.py:117). A wrong value here goes straight into the score. This applies to .ant models too, since an antimony model is converted to SBML and read by the same scanner.

pybnf/petab/export.py:1427 (_numeric_nominal) passes the value to pybnf/petab/conditions.py:183-195, which turns a relative condition into an absolute number in conditions.tsv. This is the path the original report described.

The species_initial field has no consumer in pybnf/ at all. It is only asserted in tests/test_petab_sbml_scanner.py:135. So the species case, which is the more common one in real models, is latent rather than live. Compartment sizes are folded in with parameters and behave the same way.

Two symptoms, depending on whether the file writes a placeholder attribute

When the model writes both an attribute and an initial assignment, the reader hands out the attribute silently. That is the wrong-number case.

When the model writes no attribute, the symbol is still in the namespace but has no value, so a formula naming it fails at materialize time with "references X, which is neither a simulation-output column nor a fit/model parameter". That message is wrong, since X is a model parameter, and it sits in a branch commented "Validation at compile time should make this unreachable". Antimony produces exactly this shape. Checked: k_derived = k_base + 1 compiles to <parameter id="k_derived" constant="true"/> with no value attribute plus an initial assignment.

Exposure in this repository: none

20 SBML files, 3 initial assignments in total, all of them on species: STAT5A and STAT5B in tests/petab_fixtures/boehm_v2, and A in tests/full_tests/T2-ade-abcd. Since nothing consumes species_initial, none of them reach a live path. 5 antimony models, none with a derived parameter. So there is no wrong committed artifact and no failing case in the suite.

Exposure outside it: common

I sampled 25 models from the public PEtab benchmark collection. 16 contain at least one initial assignment. 31 of those are on a parameter or a compartment, and for 28 of them the reader still hands out a stale attribute value. 54 more are on species.

Two concrete cases. Bertozzi_PNAS2020 defines beta_N as (R0_ * gamma_) / N_ while the attribute says 0.0, so a formula over beta_N would be scored with zero today. Laske_PLOSComputBiol2019 has 27 aliases of the form ModelValue_82 = D_rib, where the attribute currently matches the aliased parameter, so those are right until someone estimates D_rib and then quietly stop being right.

Which fix to pick

I measured the shape of those 31 assignments: 30 refer to other entities and exactly 1 is arithmetic over numbers alone. So a fix that evaluates self-contained assignments recovers almost nothing. The value is in the refusal.

Proposed shape. Read listOfInitialAssignments. If the right hand side has no identifiers, evaluate it and use that as the value. Otherwise the entity has no value the model file can settle, so leave it out of parameter_values. The exporter then refuses the relative condition through the path that already exists at conditions.py:183-188, and the measurement layer no longer binds a wrong constant. The MathML serializer needed for this is already in the file, as _assignment_rule_formula and _serialize_mathml, written for assignment rules in #465.

Two loose ends for whoever picks this up. First, once the value is dropped, a formula naming the symbol hits the misleading error above, so the namespace side has to be decided at the same time: either exclude the symbol and inline the right hand side the way #465 does for an assignment rule, or keep it and raise an error that names the real cause. Inlining is only sound when the right hand side is over parameters and compartments, because a species initial value is not that species value at time t. Second, the refusal message at conditions.py:184-188 is worded for BNGL expressions and would need an SBML wording.

Both paths are testable without a simulator, which is what the repro in the first comment does.

For what it is worth this lands on the same rule the petab maintainers chose for BioNetGen models in PEtab-dev/libpetab-python#517: do not report a value the model file alone cannot settle.