Protonation invariant: detect and repair rows that are not protonations; four root causes fixed; the 26.1 run rolled out

[freiburgermsu]

Follow-up to #296, which merged into dev on 2026-09-23 (f7527bce); the repaired rows depend on its a7200948. It was held on the fork as freiburgermsu#7 while #296 was pending and is proposed here unchanged, as its own PR. It merges cleanly onto dev (git merge-tree, verified) and is independent of the manuscript PRs #297–#299.

What was found

The critique named 37 compounds carrying chemically impossible formulas — hydrogens gained with no charge change — and recommended keeping the run and fixing them separately. Investigating it found the class is larger than 37, has four independent root causes, three of them not in Marvin at all, and can be detected and closed by one rule rather than by a list.

Root cause 1 — Marvin is protonating fragments that aren't molecules

InChI disconnects metal–ligand bonds by design. For Fe₄S₄ the InChI is 4Fe.4S; Marvin 26.1 protonates four free sulfides to H₂S and the row comes back H8Fe4S4, charge unchanged. It also protonates bare atoms in both representations: elemental S/Se/P/O → H₂S/H₂Se/PH₃/H₂O. None of these is a protonation state.

A protonation moves protons, so between a row and its source ΔH must equal Δcharge with heavy atoms conserved. That single, representation-agnostic invariant reproduces the critique's count exactly and shows what it missed:

	impossible rows	compounds	shipped formulas that would change impossibly
23.4 bundle	42	22	(already shipping)
26.1 bundle	139	96	69 — 26 H-gained + 11 charge-only = the 37, plus 32 mixed

"Prefer the SMILE row" fixes most but not all, and for the wrong reason on the Fe–S cubanes: MetaCyc's own SMILES writes [SH] on the bridging sulfides, so its smiles.tsv says H4Fe4S4 where its inchi.tsv says Fe4S4 — a source-level disagreement (39 compounds), not a Marvin one.

Root cause 2 — parse_structure drops charge when its parsers fail the string

Checking each inchi.tsv row against what its InChI declares found 58 rows disagreeing with their own /q and /p layers:

• RDKit rejects hypervalent halogen oxides (chlorate, InChI=1S/ClO3/c2-1(3)4/q-1); the OpenBabel fallback warns "Charge(s): Do not match" and returns 0. Nine anions stored as neutral radicals.
• On multi-component strings with both /q and /p (the Mg porphyrins, .../q-1;+2/p-1), RDKit keeps the pre-/p hydrogen count — 49 rows stored one H and one charge unit too many. A neutral chlorophyll shipped as a cation.

This predates 26.1 and propagated into every bundle refreshed through that function (202 bundle InChI rows).

Root cause 3 — the picker never expected two bundles

Found by the new key check on the first regeneration, before anything was committed: 5,567 InChIKey rows in the regenerated pick file were not the key of their own InChI row, and 2,006 regenerated compound records carried a formula/charge from one bundle and a SMILES from the other (95 in the shipped records). List_ModelSEED_Structures.py resolves each structure type independently, and BiochemPy.loadStructures globbed every protonations/*.tsv into the Charged stage — so with 23.4 and 26.1 both present, at different protonation states for 6,893 compounds, it took the InChI from one vintage and the key or SMILES from the other. The shipped pick file was last regenerated 2026-07-04, before the 26.1 bundle existed. Any regeneration after #289 would have shipped this, with or without the rest of this branch.

Root cause 4 — the reaction refresh left old charges in every stoichiometry

Rebuild_Stoichiometry.py refreshes the formula embedded in each reaction's stoichiometry entries but never the embedded charge, and balanceReaction() reads that embedded pair. Running the documented Refresh_Reactions.sh after the records changed, every recharged compound left its old charge behind in every reaction it appears in: the rebalance saw a phantom hydrogen imbalance, Adjust_Reaction_Protons.py "fixed" it by adding a proton, and 5,935 reactions that were balanced against their records went OK → CI:1, with proton coefficients rewritten on 17,457. Latent until a refresh recharged thousands of compounds at once — which is exactly what keeping this run does. Found by recomputing the imbalance from the records rather than trusting the status column.

And Rebalance_Reactions.py rejected its own documented invocation — its argument guard admitted verbose but not the save that Refresh_Reactions.sh passes — so the rebalance step of the sequence has silently not been running at all.

Two more, found along the way

• Print_Structure_Formula_Charge.py had never refreshed inchi.tsv/smiles.tsv since the layout migration — it looks for a structure column those files don't have, so every source row was silently skipped. Found because a refresh that should have moved 58 rows moved none.
• Its OpenBabel path strips the R group the run script puts on wildcard SMILES — 8 rows of the 26.1 MetaCyc bundle lose their R on any refresh.

What changed (commits 1–7: the fix)

1. parse_structure derives a standard InChI's formula and charge from its own layers (inchi_layers) — deterministic, parser-version-independent. Over 44,297 source InChIs it agrees with the parsed result on 44,239 and differs on exactly the 58 above. Source-file refresh fixed; R-stripping fixed; --types scopes a refresh. This commit carries the InChI-only refresh (58 source rows, 202 bundle rows). SMILE rows are deliberately not refreshed here — a full refresh also fills 109 previously-empty 23.4 rows the current parsers can read, and is its own change.
2. Validate_Protonations.py — five checks (row, bundle, source, layers, key), strict vs INFO kinds, --fail-on-violation for CI, an allowlist file, --tsv with ModelSEED ids. Repair_Protonation_Rows.py — the rule: a row failing the invariant is replaced by its unprotonated source row, its InChIKey row re-hashed, every replacement recorded in _reports/<bundle>_passthrough_<source>.tsv. Dry-run by default.
3. 26.1 bundle repaired: 130 rows + 96 InChIKey re-hashes. Validates with zero strict failures.
4. 23.4 bundle repaired: 38 rows — at that point it was still consumed (loadStructures globbed both); commit 7 stops consuming it, and it stays validated. Separable: drop this commit and 23.4 stays as it was. Plus the regression test, Scripts/Tests/test_protonation_invariant.py.
5. The gate in Run_Marvin_Protonations.py: the same compare() applied to each row as it's written (passthrough_invariant in the stats; InChIKey hashed from whichever InChI is actually written). Smoke-tested on KEGG --limit 100: elemental sulfur passes through in both rows, output validates clean. A bundle repaired post hoc and one gated at generation are now identical by construction.
6. Docs: README, per-bundle invariant_gate notes in sources.yaml, a section in the run report.
7. One consumed protonation bundle per source. sources.yaml already marks pKa tables consumed_by_production and Add_Compound_pKa_Sources.py is manifest-driven on it; protonation bundles now carry the same flag (26.1 true, 23.4 false — kept for provenance and still validated) and loadStructures honours it. Absent flag or manifest reads as consumed, so nothing else changes behaviour.

Why "pass the source through" is safe: the rule only touches rows the invariant rejects. The 297 disconnected-InChI compounds Marvin legitimately protonates — the /p-2 chlorophylls and cobalamins that the "protonate rather than pass through" decision in #296 exists to protect — all pass, and are untouched by construction.

One defect in my own first cut, caught by the gate before it shipped: the first repair replaced InChI rows but left their InChIKey rows hashed from the discarded protonated string, so elemental sulfur regenerated keyed as hydrosulfide. The key check exists because of it, the repair now re-hashes, and the commits were rebuilt.

What changed (commits 8–15: the rollout — "keep the run")

Each of the three documented steps run unchanged, in order, with an acceptance check between:

step	what moves	acceptance
List_ModelSEED_Structures.py	picked InChI formula/charge for 6,893 compounds	0 impossible against any source row of the compound; 0 InChIKey rows not the key of their InChI row
Update_Compound_Structures_Formulas_Charge.py	formula/charge on 7,544 records (7,434 a protonation vs shipped, 110 corrections of values that were already wrong); SMILES/InChIKey only on 3,738	0 impossible against source
commit 10: fix Rebalance_Reactions.py's guard and Rebuild_Stoichiometry.py's missing charge refresh	—	(the two pipeline bugs above; prerequisite for the next row)
Refresh_Reactions.sh → Refresh_Aliases.sh → Reprint_Biochemistry.py	stoichiometry formulas/charges on 33,005 reactions (57,280 embedded charges refreshed); proton coefficients on 17,457 and water on 23, each resolving a paired H/charge imbalance; status on 8,490 — 5,288 of them charge-imbalanced → OK (26.1's states are more mutually consistent than 23.4's); 63 reactions newly flagged obsolete as duplicates, 18 restored, links/aliases propagated on 81 — flagged per the README, not removed	every embedded formula/charge equals its record; Reprint is a fixed point after one pass (it re-renders 7,412 equation/definition strings); Check_Database_Consistency.py passes on 56,012 reactions and 45,708 compounds

The three compounds the critique named: cpd24682 stays Fe4S4/0 and cpd10112 stays C18H15ClSn/0 (both would have become H8Fe4S4 and C18H18ClSn); cpd33420 goes from Mo7O24/+12 to Mo7O24/−6. Elemental Se and P stop shipping as HSe⁻ and PH₃.

One number to read carefully

Each record's formula/charge checked against its own SMILES field, shipped → regenerated: disagreeing by more than a protonation 96 → 27, heavy-atom disagreement 161 → 151, off by exactly a protonation 95 → 224. Of the 224, 161 are new and every one is a compound whose InChI is more fragmented than its SMILES: Marvin protonated ligands InChI had detached from their metal, the protonation is charge-consistent so it passes the invariant, and the pick order (InChI@Charged first) feeds that formula into the record while the SMILES field carries the connected molecule — ferricyanide records as C6H3FeN6/0 with a SMILES of C6FeN6/−3. 160 of the 161 were consistent in the shipped records only because the 23.4 InChI row was unprotonated. This is #296's "protonate the InChI form" decision surfacing in the records; #296's own regeneration ships the same 161. Netting them out, the pre-existing class falls from 95 to 63.

The invariant cannot see this class by design. Adopted as commits 12–15 (the pick rule): when a compound's InChI is the more fragmented of its two representations, its formula and charge come from its least-fragmented SMILE structure — the test Run_Marvin_Protonations.py already applies to choose its input — applied in List_ModelSEED_Structures.py, where the compound-level values are decided. One guard, found on the first run: only switch when the two sources agree up to a protonation; where they don't (MetaCyc's [SH] on the Fe₄S₄ bridging sulfides), the SMILE row is no more trustworthy and the compound keeps the InChI-derived formula. 338 compounds take the route (_reports/Formula_From_SMILE_Row.txt; tagged in Pick_Reasons.txt), 226 records change, 0 are inconsistent with any source row, and the formula-vs-own-SMILES class goes 224 → 10 (below the shipped baseline of 95). The reaction refresh then touches 479 reactions and leaves the balanced count at 33,496.

What this costs, stated plainly

Keeping the run recharges ~7,500 compounds. The reaction refresh in commit 9 changes H⁺ stoichiometry and balance status on every reaction touched by them. The thermodynamic layer, the evidence grades, and every count in the manuscript PRs #297–#300 were measured before this and are stale against it; they need their own regeneration after merge (re-run Papers/NAR_Update_2026/analysis/population_basis_table.py first — it will say exactly what moved). Not attempted here.

Commits 8, 9, 11, 14 and 15 are tens of thousands of lines of regenerated data; commits 1–7, 10 and 13 are the ~20-file fix. They can be split into two PRs on request.

Left for curation, listed by the validator

• 39 compounds whose inchi.tsv and smiles.tsv disagree beyond a protonation (source check) — the [SH] Fe–S cubanes and similar. The invariant cannot pick a side; a curator can.
• 108 SMILE rows with no source formula at all (unchecked) and 67 wildcard rows whose H count is convention-dependent — reported, never failed.
• The 27 rows that were impossible in both bundles (chlorate as a neutral radical, Se as H₂Se) are now fixed in both, but they had been the shipped values for years; downstream users of the 23.4-era formulas may want to know.

Verification

• Scripts/Tests/test_protonation_invariant.py passes: invariant cases, inchi_layers cases, both bundles with zero strict row failures, every InChIKey row the key of its InChI row, inchi.tsv agreeing with its own layers.
• Every repaired file was rebuilt from the base commit in order and compared byte-for-byte against the working copy before committing.
• Acceptance projection on the regenerated pick file and compound records: zero formula changes impossible against any source row of the compound.
• The reaction refresh was run three times before it was right: once revealing that the rebalance step never runs, once revealing the embedded-charge bug (the OK → CI:1 mass), once clean. Only the clean run is committed.

Relationship to the manuscript PRs

#297, #298 and #299 are text-only and based on dev; nothing in them depends on this PR. Once this merges, a handful of their counts move — most by the 19 reactions the refresh here retires (48,403 → 48,384 live), and M13's mass-and-charge-balanced share from 58% to 69% because the refresh now persists the proton-rebalance step and reads refreshed embedded charges. That re-measurement is prepared as one commit per branch on the fork (nar2026-transport-and-template-scope-after-invariant-gate, nar2026-llm-ensemble-role-after-invariant-gate, nar2026-figures-and-growth-basis-after-invariant-gate) and is applied after this lands, not before.

🤖 Generated with Claude Code

[samseaver]

Looks good per a local review, thanks @freiburgermsu !