consumeSomething Tries to consume any metabolite using as few reactions as possible. The intended use is when you want to make sure that you model cannot consume anything without producing something. It is intended to be used with no active exchange reactions. model a model structure ignoreMets either a cell array of metabolite IDs, a logical vector with the same number of elements as metabolites in the model, of a vector of indexes for metabolites to exclude from this analysis (opt, default []) isNames true if the supplied mets represent metabolite names (as opposed to IDs). This is a way to delete metabolites in several compartments at once without knowing the exact IDs. This only works if ignoreMets is a cell array (opt, default false) minNrFluxes solves the MILP problem of minimizing the number of fluxes instead of the sum. Slower, but can be used if the sum gives too many fluxes (opt, default false) params parameter structure as used by getMILPParams (opt) solution flux vector for the solution metabolite the index of the metabolite(s) which was consumed. If possible only one metabolite is reported, but there are situations where metabolites can only be consumed in pairs (or more) NOTE: This works by forcing at least 1 unit of "any metabolites" to be consumed and then minimize for the sum of fluxes. If more than one metabolite is consumed, it picks one of them to be consumed and then minimizes for the sum of fluxes. Usage: [solution metabolite]=consumeSomething(model,ignoreMets,isNames,... minNrFluxes,params) Rasmus Agren, 2013-02-06
0001 function [solution metabolite]=consumeSomething(model,ignoreMets,isNames,minNrFluxes,params) 0002 % consumeSomething 0003 % Tries to consume any metabolite using as few reactions as possible. 0004 % The intended use is when you want to make sure that you model cannot 0005 % consume anything without producing something. It is intended to be used 0006 % with no active exchange reactions. 0007 % 0008 % model a model structure 0009 % ignoreMets either a cell array of metabolite IDs, a logical vector 0010 % with the same number of elements as metabolites in the model, 0011 % of a vector of indexes for metabolites to exclude from 0012 % this analysis (opt, default []) 0013 % isNames true if the supplied mets represent metabolite names 0014 % (as opposed to IDs). This is a way to delete 0015 % metabolites in several compartments at once without 0016 % knowing the exact IDs. This only works if ignoreMets 0017 % is a cell array (opt, default false) 0018 % minNrFluxes solves the MILP problem of minimizing the number of 0019 % fluxes instead of the sum. Slower, but can be 0020 % used if the sum gives too many fluxes (opt, default 0021 % false) 0022 % params parameter structure as used by getMILPParams (opt) 0023 % 0024 % solution flux vector for the solution 0025 % metabolite the index of the metabolite(s) which was consumed. If 0026 % possible only one metabolite is reported, but there are 0027 % situations where metabolites can only be consumed in 0028 % pairs (or more) 0029 % 0030 % NOTE: This works by forcing at least 1 unit of "any metabolites" to be 0031 % consumed and then minimize for the sum of fluxes. If more than one 0032 % metabolite is consumed, it picks one of them to be consumed and then 0033 % minimizes for the sum of fluxes. 0034 % 0035 % Usage: [solution metabolite]=consumeSomething(model,ignoreMets,isNames,... 0036 % minNrFluxes,params) 0037 % 0038 % Rasmus Agren, 2013-02-06 0039 % 0040 0041 if nargin<2 0042 ignoreMets=[]; 0043 end 0044 0045 if nargin<3 0046 isNames=false; 0047 end 0048 0049 if nargin<4 0050 minNrFluxes=false; 0051 end 0052 0053 if nargin<5 0054 params.relGap=0.8; 0055 end 0056 0057 if isNames==true && ~isempty(ignoreMets) 0058 %Check that metsToRemove is a cell array 0059 if iscellstr(ignoreMets)==false 0060 throw(MException('','Must supply a cell array of strings if isNames=true')); 0061 end 0062 end 0063 0064 if isNames==false 0065 indexesToIgnore=getIndexes(model,ignoreMets,'mets'); 0066 else 0067 indexesToIgnore=[]; 0068 for i=1:numel(ignoreMets) 0069 indexesToIgnore=[indexesToIgnore;find(strcmp(ignoreMets(i),model.metNames))]; 0070 end 0071 end 0072 0073 solution=[]; 0074 metabolite=[]; 0075 0076 nRxns=numel(model.rxns); 0077 nMets=numel(model.mets); 0078 0079 %Add uptake reactions for all metabolites 0080 model.S=[model.S speye(nMets)]; 0081 0082 %Add so that they all consume a fake metabolite 0083 model.S=[model.S;[sparse(1,nRxns) ones(1,nMets)*-1]]; 0084 0085 %Change so that the ignoreMets have a coefficient 0 in their reactions. 0086 %Does not remove the actual reaction to make mapping easier later 0087 model.S(:,indexesToIgnore+nRxns)=0; 0088 0089 %Add an uptake reaction for this last fake metabolite 0090 model.S(size(model.S,1),size(model.S,2)+1)=1; 0091 model.b=[model.b;zeros(1,size(model.b,2))]; 0092 model.lb=[model.lb;zeros(nMets,1);1]; 0093 model.ub=[model.ub;inf(nMets+1,1)]; 0094 model.rev=[model.rev;zeros(nMets+1,1)]; 0095 model.c=zeros(size(model.S,2),1); 0096 0097 sol=solveLP(model,1); 0098 if any(sol.x) 0099 %It could be that several metabolites were consumed in order to get the 0100 %best solution. 0101 %The setdiff is to avoid including the last fake metabolite 0102 I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1)); 0103 0104 if any(I) %This should always be true 0105 %Change the coefficients so that only the first is 0106 %consumed. This is not always possible, but it is tested for since it it 0107 %results in more easily interpretable results 0108 0109 oldS=model.S; 0110 foundSingle=false; 0111 %Test if any of the metabolites could be consumed on their own 0112 for i=1:numel(I) 0113 model.S=oldS; 0114 J=nRxns+1:numel(model.lb)-1; 0115 J(I(i))=[]; 0116 model.S(:,J)=0; 0117 sol=solveLP(model); 0118 if any(sol.x) 0119 foundSingle=true; 0120 break; 0121 end 0122 end 0123 %This means that there was no metabolite which could be consumed on 0124 %its own. Then let all the consumeable metabolites be consumed. 0125 if foundSingle==false 0126 model.S=oldS; 0127 end 0128 if minNrFluxes==true 0129 %Has to add names for the rxns to prevent an error in 0130 %minNrFluxes 0131 model.rxns=cell(numel(model.lb),1); 0132 model.rxns(:)={'TEMP'}; 0133 model.mets=cell(size(model.b,1),1); 0134 model.mets(:)={'TEMP'}; 0135 sol=solveLP(model,3,params); 0136 else 0137 sol=solveLP(model,1); 0138 end 0139 solution=sol.x(1:nRxns); 0140 metabolite=find(sol.x(nRxns+1:end-1)>0.1); 0141 end 0142 end 0143 end