Gapfilling failed in context manager (infeasible)

[wshao1]

Is there an existing issue for this?

• [X] I have searched the existing issues

Problem description

What I tried to achieve:

• I have been trying to apply gapfill to context-specific human GEMs, but have been encountering an infeasibility error.

How I went about it:

• To showcase this error, I tried removing an essential reaction from Human1, and then applying gapfill using the original Human1 GEM as the universal model. I expect this to simply return the removed reaction.
• I also tried lowering the model.solver.configuration.tolerances.feasibility threshold to 1e-9, as previous error threads showed this helped. But this yields the same error.

Why the current behavior is a problem:

• Gapfill fails, and lowering the feasibility threshold does not resolve the issue.

Code sample

Code run:

# Load Human1 GEM
human1 = read_sbml_model('Human-GEM.xml')

# Try gapfilling an essential reaction (MAR06669)
with human1 as test_model:
    rxn = test_model.reactions.get_by_id('MAR06669')
    test_model.remove_reactions([rxn])
    test_model.solver.configuration.tolerances.feasibility = 1e-9
    solution = gapfill(test_model, human1, demand_reactions=False)
    for reaction in solution[0]:
        print(reaction.id)

Traceback:

---------------------------------------------------------------------------
Infeasible                                Traceback (most recent call last)
[~\AppData\Local\Temp\ipykernel_26524\1945771835.py] in <module>
      5     print(test_model.slim_optimize())
      6     test_model.solver.configuration.tolerances.feasibility = 1e-9
----> 7     solution = gapfill(test_model, human1, demand_reactions=False)
      8     for reaction in solution[0]:
      9         print(reaction.id)

[~\anaconda3\lib\site-packages\cobra\flux_analysis\gapfilling.py] in gapfill(model, universal, lower_bound, penalties, demand_reactions, exchange_reactions, iterations)
    405         exchange_reactions=exchange_reactions,
    406     )
--> 407     return gapfiller.fill(iterations=iterations)

[~\anaconda3\lib\site-packages\cobra\flux_analysis\gapfilling.py] in fill(self, iterations)
    294         used_reactions = []
    295         for _ in range(iterations):
--> 296             self.model.slim_optimize(
    297                 error_value=None, message="gap filling optimization failed"
    298             )

[~\anaconda3\lib\site-packages\cobra\core\model.py] in slim_optimize(self, error_value, message)
   1198             return error_value
   1199         else:
-> 1200             assert_optimal(self, message)
   1201 
   1202     def optimize(

[~\anaconda3\lib\site-packages\cobra\util\solver.py] in assert_optimal(model, message)
    588     if status != OPTIMAL:
    589         exception_cls = OPTLANG_TO_EXCEPTIONS_DICT.get(status, OptimizationError)
--> 590         raise exception_cls(f"{message} ({status}).")
    591 
    592 

Infeasible: gap filling optimization failed (infeasible).

Environment

### Package Information | Package | Version | |:--------|--------:| | cobra | 0.29.0 | ### Dependency Information | Package | Version | |:--------------------|------------:| | appdirs | 1.4.4 | | black | 22.6.0 | | bumpversion | **missing** | | depinfo | 2.2.0 | | diskcache | 5.4.0 | | future | 0.18.2 | | httpx | 0.25.2 | | importlib-resources | 5.12.0 | | isort | 5.9.3 | | numpy | 1.21.6 | | optlang | 1.8.1 | | pandas | 1.3.5 | | pydantic | 1.10.7 | | python-libsbml | 5.19.7 | | rich | 13.3.2 | | ruamel.yaml | 0.17.21 | | scipy | 1.7.3 | | swiglpk | 5.0.10 | | tox | **missing** | ### Build Tools Information | Package | Version | |:-----------|--------:| | conda | 23.1.0 | | pip | 24.0 | | setuptools | 69.0.2 | | wheel | 0.37.1 | ### Platform Information | | | |:--------|---------:| | Windows | 10-AMD64 | | CPython | 3.9.13 |

Anything else?

No response

[cdiener]

Hi, thanks for providing a reproducible example. In this case not an issue with the gapfiller (at least for GLPK) but rather an edge case due to the use of then context manager.

In this line

# ...
with human1 as test_model:
# ...

human1 and test_model point to the same object/model. So when you remove the reaction from test_model you also remove it from human1. The with block only reverts all changes after the exiting but it does not take a copy of the model. Substituting this with a copy fixes it, at least for me:

In [70]: # Load Human1 GEM
    ...: human1 = read_sbml_model('Human-GEM.xml')
    ...: human1.solver = "glpk"
    ...: 
    ...: # Try gapfilling an essential reaction (MAR06669)
    ...: with human1.copy() as test_model:  # <-- note the copy here
    ...:     rxn = test_model.reactions.get_by_id('MAR06669')
    ...:     test_model.remove_reactions([rxn])
    ...:     solution = gapfill(test_model, human1, demand_reactions=False)
    ...:     for reaction in solution[0]:
    ...:         print(reaction.id)
    ...: 
MAR06669

I did notice an issue with CPLEX due to a low integer primal and the lack of an option to set this via gapfill. I will send a PR for this part.

[wshao1]

Ah thank you for catching that! I am still having issues with trying to use gapfill however. My goal is to identify alternative reactions/pathways to an essential reaction in a context-specific model. My reasoning was this: if a reaction is essential in the context-specific model but not in the universal model, then there should be alternative reactions/pathways. To identify these, I removed the essential reaction from both the context-specific and universal models, then tried to gap-fill. But I still get the infeasibility error. Is there an error in this reasoning?

Code:

# Load Human1 GEM and context-specific model
human1 = read_sbml_model('Human-GEM.xml')
context_specific_model= read_sbml_model('context_specific_model.xml')

# find alternative reactions to compensate for MAR01820
with context_specific_model as pruned_model:
    with human1 as univ_model:
        pruned_model.remove_reactions([pruned_model.reactions.MAR01820])
        univ_model.remove_reactions([univ_model.reactions.MAR01820])
        test_model.solver.configuration.tolerances.feasibility = 1e-9
        solution = gapfill(pruned_model, univ_model, demand_reactions=False)
        for reaction in solution[0]:
            print(reaction.id)

Context-specific model: context_specific_model.zip

[wshao1]

Also, the original Human1 check still breaks when I try to gapfill some other essential reactions (e.g. MAR06657, MAR00578, MAR06702).