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