Home > RAVEN > getMinNrFluxes.m

getMinNrFluxes

PURPOSE ^

getMinNrFluxes

SYNOPSIS ^

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

DESCRIPTION ^

 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]=getMinNrFluxes(model, toMinimize, params, scores)

   Rasmus Agren, 2013-08-01

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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]=getMinNrFluxes(model, toMinimize, params, scores)
0029 %
0030 %   Rasmus Agren, 2013-08-01
0031 %
0032 
0033 exitFlag=1;
0034 
0035 if nargin<2
0036     toMinimize=model.rxns;
0037 else
0038     if ~iscell(toMinimize)
0039        toMinimize=model.rxns(toMinimize); 
0040     end
0041 end
0042 
0043 %For passing parameters to the solver
0044 if nargin<3
0045     params=[];
0046 end
0047 echo=0;
0048 if isfield(params,'printReport')
0049     if params.printReport==true
0050         echo=3;
0051     end
0052 end
0053 
0054 if nargin<4
0055     %It says that the default is -1, but that is to fit with other code
0056     scores=ones(numel(toMinimize),1)*1;
0057 else
0058     if numel(scores)~=numel(toMinimize)
0059         dispEM('The number of scores must be the same as the number of reactions to minimize');
0060     end
0061     
0062     %Change positive scores to have a small negative weight. This is a
0063     %temporary solution.
0064     scores(scores>=0)=max(scores(scores<0));
0065 
0066     %It says that the default is -1, but that is to fit with other code
0067     scores=scores*-1;
0068 end
0069 
0070 %Check if the model is in irreversible format
0071 if any(model.rev)
0072     %Convert the model to irreversible format
0073     irrevModel=convertToIrrev(model);
0074     
0075     %Find the indexes for the reactions in toMinimize
0076     [indexes I]=ismember(strrep(irrevModel.rxns,'_REV',''),toMinimize);
0077 else
0078     irrevModel=model;
0079     
0080     %Find the indexes for the reactions in toMinimize
0081     [indexes I]=ismember(irrevModel.rxns,toMinimize);
0082 end
0083 
0084 indexes=find(indexes);
0085 %Adjust scores to fit with reversible
0086 scores=scores(I(indexes));
0087 
0088 %Add binary constraints in the following manner:
0089 %-  Add one unique "metabolite" for each integer reaction as a substrate.
0090 %   These metabolites can have net production
0091 %-  Add reactions for the production of each of those metabolites. The
0092 %   amount produced in one reaction unit must be larger than the largest
0093 %   possible flux in the model (but not too large to avoid bad scaling)
0094 
0095 %Calculate a solution to the problem without any constraints. This is to
0096 %get an estimate about the magnitude of fluxes in the model and to get a
0097 %feasible start solution.
0098 sol=solveLP(irrevModel,1);
0099 
0100 %Return an empty solution if the non-constrained problem couldn't be solved
0101 if isempty(sol.x)
0102     x=[];
0103     I=[];
0104     exitFlag=-1;
0105     return;
0106 end
0107 
0108 %Take the maximal times 5 to have a safe margin. If it's smaller than 1000,
0109 %then use 1000 instead.
0110 maxFlux=max(max(sol.x)*5,1000);
0111 
0112 prob.c=[zeros(numel(irrevModel.rxns),1);scores(:)]; %Minimize the number of fluxes
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 
0122 intArray=speye(numel(irrevModel.rxns))*-1;
0123 intArray=intArray(indexes,:);
0124 prob.a=[irrevModel.S;intArray];
0125 a=[sparse(numel(irrevModel.mets),numel(indexes));speye(numel(indexes))*maxFlux];
0126 prob.a=[prob.a a];
0127 prob.ints.sub=numel(irrevModel.rxns)+1:numel(irrevModel.rxns)+numel(indexes);
0128 
0129 %Use the output from the linear solution as starting point. Only the values
0130 %for the integer variables will be used, but all are supplied.
0131 prob.sol.int.xx=zeros(numel(prob.c),1);
0132 prob.sol.int.xx(prob.ints.sub(sol.x(indexes)>10^-7))=1;
0133 
0134 % Optimize the problem
0135 [crap,res] = mosekopt(['minimize echo(' num2str(echo) ')'], prob,getMILPParams(params));
0136 isFeasible=checkSolution(res);
0137 
0138 if ~isFeasible
0139     x=[];
0140     I=[];
0141     exitFlag=-1;
0142     return;
0143 end
0144 
0145 xx=res.sol.int.xx(1:numel(irrevModel.rxns));
0146 I=res.sol.int.xx(numel(xx)+1:end);
0147 
0148 %Check if Mosek aborted because it reached the time limit
0149 if strcmp('MSK_RES_TRM_MAX_TIME',res.rcodestr)
0150     exitFlag=-2;
0151 end
0152 
0153 %Map back to original model from irrevModel
0154 x=xx(1:numel(model.rxns));
0155 if numel(irrevModel.rxns)>numel(model.rxns)
0156     x(model.rev~=0)=x(model.rev~=0)-xx(numel(model.rxns)+1:end);
0157 end
0158 
0159 I=ismember(toMinimize,strrep(irrevModel.rxns(indexes(I>10^-7)),'_REV',''));
0160 end

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005