opencobra/cobratoolbox — pull request #15 by bbarker (2013-11-20 19:12:06+00:00)
Title: Improving convertIrrevFluxDistribution.m and convertToIrreversible.m

I've revised based on your suggestions - thanks again for those, it has
been helpful. Please see the previous pull request for the original intent
of the features:
https://github.com/opencobra/cobratoolbox/pull/7

First the good news: testOptKnock works identically with the submitted
code. Now some bad news:

In the E. coli core model, there are several reactions (see below)
whose bounds disagree with the reversibility flag; I realize
that at the model annotation level, this is currently OK - and (unfortunately)
it doesn't seem to be standardized. Correcting for this in
convertToIrreversible.m causes the SMILEY algorithm (run from testGrowthExpMatch) to find no solution.

As a non-expert SMILEY user, I really do not know what to make of this,
but I suppose it is possible that the SMILEY code or algorithm may lack some robustness to this form of model annotation.

>> load('./ecoli_core_model.mat');
>> lbRxns = find(model.rev & (model.lb >= 0)) % These are all exchange reactions.

lbRxns =

20
21
22
24
25
26
27
29
30
33
34
38
39
>> ubRxns = find(model.rev & (model.ub <= 0))

ubRxns =

Empty matrix: 0-by-1

A related problem will occur in isSameCobraModel.m (which is used in a test) if
the models have discrepancies between lb/ub and rev.

A possible solution (that doesn't involve looking too much at SMILEY) is to have convertToIrreversible ignore exchange reactions by calling findExcRxns.

Another possible solution is to take the opposite approach and assign all
reversible reactions to have negative lb and positive ub. But I feel we are
getting in to annotation problems here. Also I haven't yet tested to see how growthExpMatch would work in the test case in this scenario.

Yet another solution is to leave the .rev annotation alone, but it is already
changed in convertToIrreversible.m in the case where the reaction is operating strictly in the negative direction (but NOT strictly in the positive direction, which seems a bit strange). The pull request I've submitted removes these inconsistencies but causes the above problems.

Truthfully, these rev versus lb/ub issues have caused me a little
consternation over the last year and I really wonder if it would be better that .rev didn't exist at all. Take as an example from programming the OOP or abstract
interface approach with get/set methods. You really don't want different interfaces
to change the behavior of the same variable (or in this case, different variables
basically describing the same object); it increases the code complexity
and likelihood of introducing errors.

In the short term, maybe the right thing to do is to never use it in an algorithm implementation and deprecate its use. In the long term, maybe remove the flag altogether from new models? Either way I'd really like to hear what other people think.

Going a little further down the road of discussing annotation (and further afield from the current problems stated above): it seems that .rev and .lb/.ub are condition dependent anyway; in microbes it might change depending on the media, in animals it might change depending on the
tissue, etc. I would think that a better approach would be to have multiple sets
of lbs/ubs published for each model - the "base" lb/ub should only contain directionality constraints for non-biological reactions (e.g. biomass, IsA, sources/sinks not normally present biologically).