Using box constraints (+-1000) versus using infinity constraints (+-inf) for unknown bounds can lead to different results

[sergarcia]

For example, consider model_M with a large scalar M for unknown bounds:

model_INF = model_M; 
model_INF.lb(model_M.lb ==-M) = -inf;
model_INF.ub(model_M.ub ==M) = inf;

The outcome of findBlockedReaction for model_M is different than that of model_INF (in my case the later having more blocked reactions). The reason behind this is that when the LP solver maximizes/minimizes those reactions with unknown bounds, which also happen to be unbounded, it returns a status of unbounded which the fluxVariablity function treats as 0. In fact, the following code solves the problem (i.e. the blocked reactions are the same in both models), but the issue can lead erroneous calculations if someone is unaware of this:

model_INF = model_M;
[minFlux,maxFlux] = fluxVariability(model_M,0);
model_INF.lb(model_M.lb ==-M & minFlux ~=-M) = -inf;
model_INF.ub(model_M.ub ==M & maxFlux ~=M) = inf;

I hereby confirm that I have:

• [y] Tried to solve the issue on my own
• [y] Retried to run my code with the latest version of The COBRA Toolbox
• [y] Checked that a similar issue has not already been opened

[rmtfleming]

Thanks for contributing this Sergio.