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) ignoreIntBounds true if internal bounds (including reversibility) should be ignored. Exchange reactions are not affected. This can be used to find unbalanced solutions which are not possible using the default constraints (opt, default false) 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,ignoreIntBounds) Rasmus Agren, 2013-05-66
0001 function [solution metabolite]=consumeSomething(model,ignoreMets,isNames,minNrFluxes,params,ignoreIntBounds) 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 % ignoreIntBounds true if internal bounds (including reversibility) 0024 % should be ignored. Exchange reactions are not affected. 0025 % This can be used to find unbalanced solutions which are 0026 % not possible using the default constraints (opt, 0027 % default false) 0028 % 0029 % solution flux vector for the solution 0030 % metabolite the index of the metabolite(s) which was consumed. If 0031 % possible only one metabolite is reported, but there are 0032 % situations where metabolites can only be consumed in 0033 % pairs (or more) 0034 % 0035 % NOTE: This works by forcing at least 1 unit of "any metabolites" to be 0036 % consumed and then minimize for the sum of fluxes. If more than one 0037 % metabolite is consumed, it picks one of them to be consumed and then 0038 % minimizes for the sum of fluxes. 0039 % 0040 % Usage: [solution metabolite]=consumeSomething(model,ignoreMets,isNames,... 0041 % minNrFluxes,params,ignoreIntBounds) 0042 % 0043 % Rasmus Agren, 2013-05-66 0044 % 0045 0046 if nargin<2 0047 ignoreMets=[]; 0048 end 0049 if nargin<3 0050 isNames=false; 0051 end 0052 if nargin<4 0053 minNrFluxes=false; 0054 end 0055 if nargin<5 0056 params.relGap=0.8; 0057 end 0058 if nargin<6 0059 ignoreIntBounds=false; 0060 end 0061 0062 if isNames==true && ~isempty(ignoreMets) 0063 %Check that metsToRemove is a cell array 0064 if iscellstr(ignoreMets)==false 0065 throw(MException('','Must supply a cell array of strings if isNames=true')); 0066 end 0067 end 0068 0069 if isNames==false 0070 indexesToIgnore=getIndexes(model,ignoreMets,'mets'); 0071 else 0072 indexesToIgnore=[]; 0073 for i=1:numel(ignoreMets) 0074 indexesToIgnore=[indexesToIgnore;find(strcmp(ignoreMets(i),model.metNames))]; 0075 end 0076 end 0077 0078 %Change all internal reactions to be unbounded in both directions 0079 if ignoreIntBounds==true 0080 [crap I]=getExchangeRxns(model); 0081 nonExc=true(numel(model.rxns),1); 0082 nonExc(I)=false; 0083 model=setParam(model,'lb',nonExc,-1000); 0084 model=setParam(model,'ub',nonExc,1000); 0085 model=setParam(model,'rev',nonExc,1); 0086 end 0087 0088 solution=[]; 0089 metabolite=[]; 0090 0091 nRxns=numel(model.rxns); 0092 nMets=numel(model.mets); 0093 0094 %Add uptake reactions for all metabolites 0095 model.S=[model.S speye(nMets)]; 0096 0097 %Add so that they all consume a fake metabolite 0098 model.S=[model.S;[sparse(1,nRxns) ones(1,nMets)*-1]]; 0099 0100 %Change so that the ignoreMets have a coefficient 0 in their reactions. 0101 %Does not remove the actual reaction to make mapping easier later 0102 model.S(:,indexesToIgnore+nRxns)=0; 0103 0104 %Add an uptake reaction for this last fake metabolite 0105 model.S(size(model.S,1),size(model.S,2)+1)=1; 0106 model.b=[model.b;zeros(1,size(model.b,2))]; 0107 model.lb=[model.lb;zeros(nMets,1);1]; 0108 model.ub=[model.ub;inf(nMets+1,1)]; 0109 model.rev=[model.rev;zeros(nMets+1,1)]; 0110 model.c=zeros(size(model.S,2),1); 0111 0112 sol=solveLP(model,1); 0113 if any(sol.x) 0114 %It could be that several metabolites were consumed in order to get the 0115 %best solution. 0116 %The setdiff is to avoid including the last fake metabolite 0117 I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1)); 0118 0119 if any(I) %This should always be true 0120 %Change the coefficients so that only the first is 0121 %consumed. This is not always possible, but it is tested for since it it 0122 %results in more easily interpretable results 0123 0124 oldS=model.S; 0125 foundSingle=false; 0126 %Test if any of the metabolites could be consumed on their own 0127 for i=1:numel(I) 0128 model.S=oldS; 0129 J=nRxns+1:numel(model.lb)-1; 0130 J(I(i))=[]; 0131 model.S(:,J)=0; 0132 sol=solveLP(model); 0133 if any(sol.x) 0134 foundSingle=true; 0135 break; 0136 end 0137 end 0138 %This means that there was no metabolite which could be consumed on 0139 %its own. Then let all the consumeable metabolites be consumed. 0140 if foundSingle==false 0141 model.S=oldS; 0142 end 0143 if minNrFluxes==true 0144 %Has to add names for the rxns to prevent an error in 0145 %minNrFluxes 0146 model.rxns=cell(numel(model.lb),1); 0147 model.rxns(:)={'TEMP'}; 0148 model.mets=cell(size(model.b,1),1); 0149 model.mets(:)={'TEMP'}; 0150 sol=solveLP(model,3,params); 0151 else 0152 sol=solveLP(model,1); 0153 end 0154 solution=sol.x(1:nRxns); 0155 metabolite=find(sol.x(nRxns+1:end-1)>0.1); 0156 end 0157 end 0158 end