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