Behavior of Gene.functional and Reaction.functional is unclear

[cdiener]

Problem description

Setting Gene.functional does currently not trigger side effects. This leads to reactions with rxn.functional == False that still carry flux. The use case for that should be described better and documented or the attribute should be made private.

Code Sample

In [31]: mod = create_test_model("textbook")

In [32]: g = list(mod.reactions.ENO.genes)[0]

In [33]: g.functional = False

In [34]: mod.reactions.ENO.functional
Out[34]: False

In [35]: mod.optimize()
Out[35]: <Solution 0.874 at 0x122beb710>

In [36]: mod.reactions.ENO.flux
Out[36]: 14.716139568742832

Actual Output

ENO carries flux.

Expected Output

ENO should not carry flux since the reactions is not functional.

Dependency Information

Please paste the output of python -c "from cobra import show_versions; show_versions()" in the details block below.

System Information

OS Darwin OS-release 17.7.0 Python 3.6.7

Package Versions

cobra 0.15.1 depinfo 1.4.0 future 0.17.1 numpy 1.16.4 optlang 1.4.3 pandas 0.24.2 pip 19.1.1 python-libsbml-experimental 5.17.2 ruamel.yaml 0.15.98 setuptools 41.0.1 six 1.12.0 swiglpk 1.5.1 tabulate 0.8.3 wheel 0.33.4

[akaviaLab]

It looks like functional just sets the flag, but doesn't propagate it upwards. gene.knock_out does propagate it upwards We can copy the code in functional from

    @functional.setter
    @resettable
    def functional(self, value: bool) -> None:
        if not isinstance(value, bool):
            raise ValueError("expected boolean")
        self._functional = value

To call knock_out

    @functional.setter
    @resettable
    def functional(self, value: bool) -> None:
        if not isinstance(value, bool):
            raise ValueError("expected boolean")
        self._functional = value
      if not value:
            self.knock_out

[akaviaLab]

Okay, what I wrote won't work, because knock_out calls functional and thus we end up in a recursion loop.

Actually, this is much more complicated, because when we make the gene functional again, we'd have to reset the bound to whatever they were. Since reversible is calculated from the bounds, we wouldn't know if the bounds of the reaction that is now functional should be (0,1000), (-1000, 1000), or perhaps a different value. It is kept in context, since bounds is resettable, but I'm not sure how to undo the changes. Also, that would only work in context, so if we're not using with model I think it wouldn't work.

Any suggestions?

[akaviaLab]

The other alternative is to modify behavior a bit - if reaction is non-functional, bounds stay the same. However, when optimizing, the reaction is temporarily set to zero (optimize looks at the functional flag), and then the model is optimized. When you reset the functionality of the gene, the functional flag on reaction is reset, so it goes back to normal. Maybe change bounds so it returns (0, 0) when reaction is non-functional without changing the actual values?

I'm not sure how it would interact with optlang and the optimization - where would we have to change it.