`loopless_solution` does not respect objective

[cfrioux]

Problem description

Hi all, I am running loopless FBA and I noticed that when I run several times the optimisation I get different results. This does not occur systematically though: some model objective functions seem to give constant results.

I observed this behaviour using glpk and cplex solvers. The differences between the results seem worse in cplex - at least in the tests I performed.

I would have expected the results of loopless FBA to be stable - or more stable at least - since for some tests I can get a positive value right after a 0 value.

Am I doing something wrong?

Code Sample

Create a minimal, complete, verifiable example.

from cobra.io import load_model
from cobra.flux_analysis.loopless import add_loopless, loopless_solution

salmonella = load_model('salmonella')

salmonella.objective = "GLCptspp"

for i in range(8):
    print(f"Loopless FBA {i}: {loopless_solution(salmonella).fluxes['GLCptspp']}")

Gives the following output with cplex 22.1.0.0:

Loopless FBA 0: 5.0
Loopless FBA 1: 3.611666666666687
Loopless FBA 2: 0.0
Loopless FBA 3: 3.333333333333346
Loopless FBA 4: 5.0
Loopless FBA 5: 5.0
Loopless FBA 6: 1.714285714285893
Loopless FBA 7: 0.0

And this output with the default glpk solver:

Loopless FBA 0: 5.000000000000094
Loopless FBA 1: 3.10171428571428
Loopless FBA 2: 2.7662337662337544
Loopless FBA 3: 3.0495930047098097
Loopless FBA 4: 3.043293240513111
Loopless FBA 5: 3.4537261236503722
Loopless FBA 6: 4.999999999999877
Loopless FBA 7: 3.0432932405131314

Context

``` System Information ================== OS Darwin OS-release 20.6.0 Python 3.7.7 Package Versions ================ appdirs 1.4.4 black ; extra == 'development' not installed bumpversion ; extra == 'development' not installed cobra 0.25.0 depinfo 1.7.0 diskcache 5.2.1 future 0.18.2 httpx 0.18.2 importlib-resources 5.2.0 isort ; extra == 'development' not installed numpy 1.21.1 optlang 1.5.2 pandas 1.3.1 pip 22.2.2 pydantic 1.8.2 python-libsbml 5.19.0 rich 10.6.0 ruamel.yaml 0.17.10 scipy 1.7.1 setuptools 57.4.0 swiglpk 5.0.3 tox 3.24.3 wheel 0.36.2 ```

[cdiener]

What do you mean by unstable? Uniqueness in FBA is only guaranteed for the objective values, not for the variables/fluxes. There may be many possible flux solutions giving you the same objective value. You can run loopless FVA to investigate the possible spread in that particular reaction flux in the absence of cycles.

[cfrioux]

Yes I agree about the uniqueness but the reaction flux I am looking at is the objective value.

[cdiener]

So sorry, should have read your code more carefully. Looks like there is a bug in loopless_solution where it does not constrain the objective when not passing in reference fluxes.

One thing to note is that the "fast" method will not remove fluxes from the objective itself.

[cfrioux]

The way I looked as the results with loopless.fluxes['GLCptspp'] could be misleading at first sight. It is because loopless.objective_value gives values of 0, which is already documented in https://github.com/opencobra/cobrapy/issues/935.

Thanks for the explanation. So if I understand correctly, the "fast" method constrains the objective to have the nominal flux, and the other reactions fluxes will be modified by removing loops, while the exact Schellenberger implementation is likely to have a reduced objective value should no loopless solution permit the nominal flux?

[cdiener]

Yes exactly. The CycleFreeFlux used in loopless_solution algorithm essentially locks in the objective and all exchanges to its previous flux value, constrains the direction of internal fluxes to its previous direction, and then minimizes the absolute sum of internal fluxes.

Yep, exactly, the Schellenberger implementation adds in constraints that make loops infeasible. That will also apply to the objective. But it is really slow (can take hours for larger models).

I will try to send in a PR for this issue and #935.