Select reversibility heuristics per thermodynamic data source (eQuilibrator gets its own cutoff)

[Cooper-Taylor]

Summary

The reaction-direction cascade had exactly one rule set — DEFAULT_HEURISTICS, derived from the Group-Contribution era — applied to every thermodynamic source. This splits it into a registry keyed by source:

set	cascade	basis
GC (default)	ATPS → ABCT → MdeltaG bounds → mMdeltaG band → lowE → =	Jankowski 2008
EQ	ATPS → ABCT → undecomposable → uncorrected transport → ln Γ ± 1σ → =	Noor 2012 index, Beber 2022 uncertainty
EQ2	same, index as a bare point estimate	Flamholz 2012

GC is the default for any source without a set of its own, so both dGPredictor variants are unchanged. The GC cascade itself is untouched.

Why

Three defects in the eQuilibrator path. (Points 1 and 2 were re-characterised after review — see the "Corrected after review" note at the bottom. The code did not change; the explanation of why it was needed was wrong.)

1. The EQ run stopped reading eQuilibrator energies in 2023. The design was correct. From c263e233 (2019) through 17c9739b (2023-09-11), Update_Reaction_eQuilibrator_Energies.py overwrote the canonical energy and stamped the note that gates the reversibility step:

reactions_dict[rxn]['deltag']    = float(eq_reactions[rxn]['dg'])
reactions_dict[rxn]['deltagerr'] = float(eq_reactions[rxn]['dge'])
if('EQU' not in notes_list): notes_list.append('EQU')

So reading deltag gated on EQU was reading eQuilibrator energies, exactly as intended.

3e50646b (2023-09-13) removed the write half and left the read half in place — the updater moved to thermodynamics['eQuilibrator'] only, dropping both the deltag overwrite and the EQU stamp. Later GC rebuilds (Convention A, ad34d6ab) and the PR #265 promotion then rewrote deltag from other sources. Today only 1,797 of 25,028 reactions with an eQuilibrator record have deltag == thermodynamics['eQuilibrator'][0]; restricted to the 17,094 still carrying EQU, only 1,586 do. EQU is now a fossil of the pre-2023 pipeline. Pointing the EQ run at the eQuilibrator sublist restores the original intent rather than inventing a new one.

2. For most σ-flagged reactions eQuilibrator returns no energy at all — we stored a transform term as if it were one. GibbsEnergyPredictor.standard_dg short-circuits when a component falls outside both the reactant- and group-contribution spans:

mu, sigma_fin, sigma_inf, residual = self.get_reaction_prediction(reaction)
if residual:
    return Q_(0, "kJ/mol").plus_minus(self.preprocess.RMSE_inf)

Literally zero — the computed mean is discarded, not projected. standard_dg_prime, which our retrieval script calls, then adds the Legendre/pH transform on top, so the nonzero ΔG′° we store is only the transform. rxn00017's −20.46 kcal/mol is the pH-7 transform of an energy eQuilibrator declined to estimate.

Of the 4,607 σ-flagged reactions: 2,895 are this zeroed case (σ exactly RMSE_inf; 61 still exactly 0.000000, having no net proton change for the transform to act on — e.g. rxn02677, rxn04059), and 1,712 carry RMSE_inf·‖σ_inf‖ and do have a real mean but with an undetermined covariance direction. Neither supports a directional call; the GC bounds rule cannot fire at that width, so all of them fell through to a permissive =.

3. Transport energies are compartment-collapsed. Retrieve_eQuilibrator_Reactions_Energies.py keys its MetaNetX formula on compound id and discards compartment, so any species on both sides nets out — 1,102 transport reactions carry a ΔG′° for a different reaction (rxn12518 at +524 kcal/mol). A further 76 non-transport reactions collapse via stereo-neutral InChIKey matching (rxn00816: D-glucose and galactose merge). Beber 2022 separately notes the transformed framework needs a −N_H·RT·ln(10^ΔpH) − Q·FΔΦ term we never apply.

Separately: the reversibility index eQuilibrator has been computing for us — column 4 of MetaNetX_Reaction_Energies.tbl — was parsed and discarded on every run, and make_ln_reversibility_index_heuristic was never invoked. It now decides 8,944 reactions, and the bare no-evidence default fallback drops from 7,120 to zero.

Results

Per-method operators vs current dev:

source	records	changed
Group contribution	56,002	0
eQuilibrator	25,028	7,103
dGPredictor	27,715	0
dGPredictor-ModelSEED	31,924	0

Only eQuilibrator moves. The three columns scored with the untouched GC rule set are byte-identical, which isolates the change to the new cutoff rather than the refactor.

Canonical reversibility: > 12,670→9,584 (−3,086) · < 2,323→1,603 (−720) · = 15,059→14,002 (−1,057) · ? 25,960→30,823 (+4,863)

The database gets less confident, and that is the point. The +4,863 unknowns track the σ-flagged records almost exactly — reactions that had been getting a permissive = on the strength of a number that, for 2,895 of them, is eQuilibrator's literal zero plus a pH transform.

Verification

• New Scripts/Tests/test_eq_heuristics.py — the ln Γ implementation reproduces eQuilibrator's own published ln_reversibility_index on 17,771 reactions; every residual is a MetaNetX-collapsed reaction, asserted rather than assumed.
• Scripts/Tests/test_reaction_direction.py now separates invariant sources (GC, dGPredictor — must match exactly) from intentionally re-scored ones (eQuilibrator), with --strict to require equality everywhere. Passes.
• Rerun_Thermodynamics.sh is idempotent: a second run reproduces every file byte-for-byte.
• Spot checks: rxn00001 → > (ln Γ −9.18) · rxn00017 → ? · rxn12518 → ? · rxn08173 ATP synthase → =.

Full write-up: Papers/NAR_Update_2026/data/equilibrator_reversibility_heuristics_2026-08-18.md.

Scope notes for review

• Canonical deltag is untouched and mostly remains the Group-Contribution value. Reconciling it is Promote_Reaction_Thermodynamics_to_Canonical.py, deliberately not part of Rerun_Thermodynamics.sh and not run here.
• The ~3,086 lost > calls are the main review surface: some are the index genuinely disagreeing with the GC bounds rule, others are reactions whose only evidence was the zeroed value. Separable from the generated report on request.
• Follow-ups worth considering: is_using_group_contribution() would replace the σ-threshold with eQuilibrator's own verdict; standard_dg_prime_multi() returns the covariance the current per-reaction loop cannot; and the MetaNetX collapse in the retrieval script is fixable at source.

Corrected after review

Commit a86f6adf revises the framing of defects 1 and 2 following review. The original text claimed the EQ run "never" read eQuilibrator energies (it did, correctly, until 2023) and that the σ marker was "read as an error bar" (understated — the energy itself is a discarded mean replaced by zero). No code changed.

🤖 Generated with Claude Code