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-08-01
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-08-01 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 dispEM('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 %Add padding to the reaction annotation to prevent an error in solveLP 0113 padding=cell(numel(model.rev),1); 0114 padding(:)={''}; 0115 model.rxns=padding; 0116 model.rxnNames=padding; 0117 model.eccodes=padding; 0118 model.rxnMiriams=padding; 0119 model.grRules=padding; 0120 if isfield(model,'genes') 0121 model.rxnGeneMat=sparse(numel(model.rev),numel(model.genes)); 0122 end 0123 model.subSystems=padding; 0124 model.rxnFrom=padding; 0125 model.rxnComps=ones(numel(model.rev),1); 0126 0127 sol=solveLP(model,1); 0128 if any(sol.x) 0129 %It could be that several metabolites were consumed in order to get the 0130 %best solution. 0131 %The setdiff is to avoid including the last fake metabolite 0132 I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1)); 0133 0134 if any(I) %This should always be true 0135 %Change the coefficients so that only the first is 0136 %consumed. This is not always possible, but it is tested for since it it 0137 %results in more easily interpretable results 0138 0139 oldS=model.S; 0140 foundSingle=false; 0141 %Test if any of the metabolites could be consumed on their own 0142 for i=1:numel(I) 0143 model.S=oldS; 0144 J=nRxns+1:numel(model.lb)-1; 0145 J(I(i))=[]; 0146 model.S(:,J)=0; 0147 sol=solveLP(model); 0148 if any(sol.x) 0149 foundSingle=true; 0150 break; 0151 end 0152 end 0153 %This means that there was no metabolite which could be consumed on 0154 %its own. Then let all the consumeable metabolites be consumed. 0155 if foundSingle==false 0156 model.S=oldS; 0157 end 0158 if minNrFluxes==true 0159 %Has to add names for the rxns to prevent an error in 0160 %minNrFluxes 0161 model.rxns=cell(numel(model.lb),1); 0162 model.rxns(:)={'TEMP'}; 0163 model.mets=cell(size(model.b,1),1); 0164 model.mets(:)={'TEMP'}; 0165 sol=solveLP(model,3,params); 0166 else 0167 sol=solveLP(model,1); 0168 end 0169 solution=sol.x(1:nRxns); 0170 metabolite=find(sol.x(nRxns+1:end-1)>0.1); 0171 end 0172 end 0173 end