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-03-14
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-03-14 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 throw(MException('','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 if res.rcode==1001 0137 throw(MException('','The Mosek licence has expired')); 0138 end 0139 if res.rcode==1010 0140 throw(MException('','The Mosek licence used only supports small problems (up to 300 variables). Have you requested the correct licence?')); 0141 end 0142 0143 xx=res.sol.int.xx(1:numel(irrevModel.rxns)); 0144 I=res.sol.int.xx(numel(xx)+1:end); 0145 0146 %Check that no bounds are violated 0147 if any(xx<irrevModel.lb-10^-7) || any(xx>irrevModel.ub+10^-7) 0148 x=[]; 0149 I=[]; 0150 exitFlag=-1; 0151 return; 0152 end 0153 0154 %Check if Mosek aborted because it reached the time limit 0155 if strcmp('MSK_RES_TRM_MAX_TIME',res.rcodestr) 0156 exitFlag=-2; 0157 end 0158 0159 %Map back to original model from irrevModel 0160 x=xx(1:numel(model.rxns)); 0161 if numel(irrevModel.rxns)>numel(model.rxns) 0162 x(model.rev~=0)=x(model.rev~=0)-xx(numel(model.rxns)+1:end); 0163 end 0164 I=ismember(toMinimize,strrep(irrevModel.rxns(indexes(I>10^-7)),'_REV','')); 0165 end