getMinNrFluxes

 getMinNrFluxes
   Returns the minimal set of fluxes that satisfy the model using
   mixed integer linear programming.

    model         a model structure
   toMinimize    either a cell array of reaction IDs, a logical vector
                 with the same number of elements as reactions in the model,
                 of a vector of indexes for the reactions that should be
                 minimized (optional, default model.rxns)
   params        *obsolete option*
   scores        vector of weights for the reactions. Negative scores
                 should not have flux. Positive scores are not possible in this
                 implementation, and they are changed to max(scores(scores<0)).
                 Must have the same dimension as toMinimize (find(toMinimize)
                 if it is a logical vector) (optional, default -1 for all reactions)

   x             the corresponding fluxes for the full model
   I             the indexes of the reactions in toMinimize that were used
                 in the solution
   exitFlag      1: optimal solution found
                -1: no feasible solution found
                -2: optimization time out

   NOTE: Uses 1000 mmol/gDW/h as an arbitary large flux. Could possibly
   cause problems if the fluxes in the model are larger than that.

 Usage: [x,I,exitFlag]=getMinNrFluxes(model, toMinimize, params, scores)

0001 function [x,I,exitFlag]=getMinNrFluxes(model, toMinimize, params,scores)
0002 % getMinNrFluxes
0003 %   Returns the minimal set of fluxes that satisfy the model using
0004 %   mixed integer linear programming.
0005 %
0006 %    model         a model structure
0007 %   toMinimize    either a cell array of reaction IDs, a logical vector
0008 %                 with the same number of elements as reactions in the model,
0009 %                 of a vector of indexes for the reactions that should be
0010 %                 minimized (optional, default model.rxns)
0011 %   params        *obsolete option*
0012 %   scores        vector of weights for the reactions. Negative scores
0013 %                 should not have flux. Positive scores are not possible in this
0014 %                 implementation, and they are changed to max(scores(scores<0)).
0015 %                 Must have the same dimension as toMinimize (find(toMinimize)
0016 %                 if it is a logical vector) (optional, default -1 for all reactions)
0017 %
0018 %   x             the corresponding fluxes for the full model
0019 %   I             the indexes of the reactions in toMinimize that were used
0020 %                 in the solution
0021 %   exitFlag      1: optimal solution found
0022 %                -1: no feasible solution found
0023 %                -2: optimization time out
0024 %
0025 %   NOTE: Uses 1000 mmol/gDW/h as an arbitary large flux. Could possibly
0026 %   cause problems if the fluxes in the model are larger than that.
0027 %
0028 % Usage: [x,I,exitFlag]=getMinNrFluxes(model, toMinimize, params, scores)
0029 
0030 exitFlag=1;
0031 
0032 if nargin<2
0033     toMinimize=model.rxns;
0034 elseif ~islogical(toMinimize) && ~isnumeric(toMinimize)
0035     toMinimize=convertCharArray(toMinimize);
0036 else
0037     toMinimize=model.rxns(toMinimize);
0038 end
0039 
0040 %For passing parameters to the solver
0041 if nargin<3
0042     params=struct();
0043 end
0044 
0045 if nargin<4
0046     %It says that the default is -1, but that is to fit with other code
0047     scores=ones(numel(toMinimize),1)*1;
0048 else
0049     if numel(scores)~=numel(toMinimize)
0050         EM='The number of scores must be the same as the number of reactions to minimize';
0051         dispEM(EM);
0052     end
0053     
0054     %Change positive scores to have a small negative weight. This is a
0055     %temporary solution.
0056     scores(scores>=0)=max(scores(scores<0));
0057     
0058     %It says that the default is -1, but that is to fit with other code
0059     scores=scores*-1;
0060 end
0061 
0062 %Check if the model is in irreversible format
0063 if any(model.rev)
0064     %Convert the model to irreversible format
0065     irrevModel=convertToIrrev(model);
0066     
0067     %Find the indexes for the reactions in toMinimize
0068     [indexes, I]=ismember(strrep(irrevModel.rxns,'_REV',''),toMinimize);
0069 else
0070     irrevModel=model;
0071     
0072     %Find the indexes for the reactions in toMinimize
0073     [indexes, I]=ismember(irrevModel.rxns,toMinimize);
0074 end
0075 
0076 indexes=find(indexes);
0077 %Adjust scores to fit with reversible
0078 scores=scores(I(indexes));
0079 
0080 %Add binary constraints in the following manner: -  Add one unique
0081 %"metabolite" for each integer reaction as a substrate.
0082 %   These metabolites can have net production
0083 %-  Add reactions for the production of each of those metabolites. The
0084 %   amount produced in one reaction unit must be larger than the largest
0085 %   possible flux in the model (but not too large to avoid bad scaling)
0086 
0087 %Calculate a solution to the problem without any constraints. This is to
0088 %get an estimate about the magnitude of fluxes in the model and to get a
0089 %feasible start solution.
0090 sol=solveLP(irrevModel,1);
0091 
0092 %Return an empty solution if the non-constrained problem couldn't be solved
0093 if isempty(sol.x)
0094     x=[];
0095     I=[];
0096     exitFlag=-1;
0097     return;
0098 end
0099 
0100 %Take the maximal times 5 to have a safe margin. If it's smaller than 1000,
0101 %then use 1000 instead.
0102 maxFlux=max(max(sol.x)*5,1000);
0103 
0104 intArray=speye(numel(irrevModel.rxns))*-1;
0105 intArray=intArray(indexes,:);
0106 prob.a=[irrevModel.S;intArray];
0107 a=[sparse(numel(irrevModel.mets),numel(indexes));speye(numel(indexes))*maxFlux];
0108 prob.a=[prob.a a];
0109 prob.ints.sub=numel(irrevModel.rxns)+1:numel(irrevModel.rxns)+numel(indexes);
0110 
0111 prob.c=[zeros(numel(irrevModel.rxns),1);scores(:);zeros(size(prob.a,1),1)]; %Minimize the number of fluxes
0112 prob.A=[prob.a -speye(size(prob.a,1))];
0113 prob.blc=[irrevModel.b(:,1);zeros(numel(indexes),1)];
0114 if size(irrevModel.b,2)==2
0115     prob.buc=[irrevModel.b(:,2);inf(numel(indexes),1)];
0116 else
0117     prob.buc=[irrevModel.b(:,1);inf(numel(indexes),1)];
0118 end
0119 prob.blx=[irrevModel.lb;zeros(numel(indexes),1)];
0120 prob.bux=[irrevModel.ub;ones(numel(indexes),1)];
0121 prob.lb = [prob.blx; prob.blc];
0122 prob.ub = [prob.bux; prob.buc];
0123 prob.osense=1;
0124 prob.csense=repmat('E', 1, size(prob.a,1),1);
0125 prob.b=zeros(size(prob.a,1), 1);
0126 
0127 %Use the output from the linear solution as starting point. Only the values
0128 %for the integer variables will be used, but all are supplied.
0129 prob.sol.int.xx=zeros(numel(prob.c),1);
0130 prob.sol.int.xx(prob.ints.sub(sol.x(indexes)>10^-12))=1;
0131 prob.x0=[];
0132 prob.vartype=repmat('C', size(prob.A,2), 1);
0133 prob.vartype(prob.ints.sub) = 'I'; % with .lb = 0 and .ub = 1, they are binary
0134 % integers (glpk in octave only allows 'continuous' or '', not 'binary')
0135 prob=rmfield(prob,{'blx','bux','blc','buc'});
0136 
0137 % Optimize the problem
0138 res = optimizeProb(prob,params,false);
0139 isFeasible=checkSolution(res);
0140 
0141 if ~isFeasible
0142     x=[];
0143     I=[];
0144     exitFlag=-1;
0145     return;
0146 end
0147 
0148 xx=res.full(1:numel(irrevModel.rxns));
0149 I=res.full(numel(xx)+1:end);
0150 
0151 %Check if Mosek aborted because it reached the time limit
0152 %TODO: modify for cobra/gurobi
0153 % if strcmp('MSK_RES_TRM_MAX_TIME',res.rcode)
0154 %     exitFlag=-2;
0155 % end
0156 
0157 %Map back to original model from irrevModel
0158 x=xx(1:numel(model.rxns));
0159 if numel(irrevModel.rxns)>numel(model.rxns)
0160     x(model.rev~=0)=x(model.rev~=0)-xx(numel(model.rxns)+1:end);
0161 end
0162 
0163 I=ismember(toMinimize,strrep(irrevModel.rxns(indexes(I>10^-12)),'_REV',''));
0164 end