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 (opt, default model.rxns)
   params        parameter structure as used by getMILPParams (opt)
   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) (opt, 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 (opt, default model.rxns)
0011 %   params        parameter structure as used by getMILPParams (opt)
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) (opt, 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