[jplfaria]

Why this matters in our workflow

We re-gap-fill models and add the results with MSGapfill.integrate_gapfill_solution. That method calls
MSModelUtil.test_solution(solution, keep_changes=True) to remove gap-filled reactions that are "not needed". If
knocking out a reaction makes the LP infeasible, the reaction is treated as not needed and deleted, although the model
cannot work without it. This happens when the model has a forced flux, such as a non-zero ATP maintenance lower bound
(as in e_coli_core), or when the knock-out LP fails for any other reason.

What happens

msmodelutl.py L541-L606 at 062714b:

rxnobj.upper_bound = 0
objective = tempmodel.slim_optimize()
if objective < solution["minobjective"]:
    ...                      # needed: restore the bound
else:
    removed_rxns.append(rxnobj)
    unneeded.append(...)     # "not needed"

slim_optimize() returns NaN for any non-optimal status (its default error_value). Then NaN < minobjective is
False, so the reaction counts as not needed. With keep_changes=True, which is what
MSGapfill.integrate_gapfill_solution passes at L369,
the reaction is removed from the model and from the solution. The model then cannot grow on the medium it was
gap-filled on. The "<" branch at L583 has the same problem.

How to reproduce

Use the package's test fixtures. e_coli_core without GLCpts_c0 is gap-filled on glucose with
template_core_bigg.json. Its ATPM_c0 has a lower bound of 8.39, so without glucose uptake the problem is infeasible.

model = get_model(["GLCpts_c0"])          # fixtures from tests/core/test_msgapfill.py
gf = MSGapfill(model, [template])
sol = gf.run_gapfilling(media_glucose_aerobic, "BIOMASS_Ecoli_core_w_GAM_c0", minimum_obj=0.1)
# sol["new"] == {"GLCpts_c0": ">"}
gf.integrate_gapfill_solution(sol, [], link_gaps_to_objective=False)
"GLCpts_c0" in model.reactions, model.slim_optimize()

	knock-out objective of GLCpts_c0	GLCpts_c0 after integrate	growth on the gap-fill medium
dev	nan	removed	nan (infeasible)
this branch	nan	kept	0.8739

The fix

+import math
 ...
-                    if objective < solution["minobjective"]:
+                    if math.isnan(objective) or objective < solution["minobjective"]:

The change is the same in both directions. An infeasible knock-out means the model cannot meet its constraints without
the reaction, so the reaction is needed. I used an explicit math.isnan instead of slim_optimize(error_value=0),
because with error_value=0 an infeasible knock-out would still count as not needed whenever minobjective is 0,
which MSGapfill(minimum_obj=0) allows.

Callers. test_solution is called only from MSGapfill.integrate_gapfill_solution (msgapfill.py:369), which
is called from MSATPCorrection.apply_growth_media_gapfilling (msatpcorrection.py:423), MSGapfill.gapfill
(msgapfill.py:523) and MSCommunity (mscommunity.py:651). For each of them, the only change is that reactions
whose knock-out is infeasible stay in the model. A feasible knock-out behaves exactly as before.

Not changed here. With keep_changes=True, a "reversed" reaction that is found not needed is removed from the
model entirely (tempmodel.remove_reactions(removed_rxns)), although it was in the model before the gap-fill.
Restoring its original bound may be the intended behaviour. That is a separate question, so this PR leaves it alone.

Tests

tests/core/test_msgapfill.py::test_integrate_gapfill_solution_keeps_reaction_with_infeasible_knockout is the
reproduction above. It checks that GLCpts_c0 stays in the model and in the solution, and that the model grows above
the gap-fill minimum.

Against dev's code:

E       AssertionError: assert 'GLCpts_c0' in [<Reaction PFK_c0 at 0x...>, <Reaction PFL_c0 at 0x...>, ...]
1 failed

With the fix, it passes. Full suite (Python 3.10, cobra 0.32.1, CPLEX 22.1.1): tests/core goes from 56 to 57
passed, and tests/ from 58 to 59. black 22.8.0 and the other pre-commit hooks pass.

Real-model check

We took a genome-scale E. coli K-12 model built with MSBuilder.build() (RAST annotations, GramNeg template v7.1).
We removed its only glucose uptake reaction, rxn05226_c0 (D-glucose transport via PEP:Pyr PTS), re-gap-filled on
minimal glucose with MSGapfill(model, [template], default_target="bio1"), and integrated with
integrate_gapfill_solution:

ATPM_c0 lower bound	gap-fill adds	dev after integrate	this branch after integrate
0	rxn05226_c0	kept; growth 0.5451	kept; growth 0.5451
3	rxn05226_c0	removed; infeasible on glucose	kept; growth 0.5319

Two full MSBuilder.build() runs (E. coli K-12 and a Pantoea strain) gave identical models on dev and on this
branch. In those builds test_solution is never reached, because the ATP correction integrated no core gap-fill.

Related

• Issue: MSModelUtil.test_solution treats an infeasible knock-out (NaN) as 'not needed' and deletes required gap-filled reactions jplfaria/ModelSEEDpy#6

🤖 Generated with Claude Code