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)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 if nargin<2 0044 ignoreMets=[]; 0045 elseif ~islogical(ignoreMets) && ~isnumeric(ignoreMets) 0046 ignoreMets=convertCharArray(ignoreMets); 0047 end 0048 if nargin<3 0049 isNames=false; 0050 end 0051 if nargin<4 0052 minNrFluxes=false; 0053 end 0054 if nargin<5 0055 params.relGap=0.8; 0056 end 0057 if nargin<6 0058 ignoreIntBounds=false; 0059 end 0060 0061 if isNames==true && ~isempty(ignoreMets) 0062 %Check that metsToRemove is a cell array 0063 if iscellstr(ignoreMets)==false 0064 EM='Must supply a cell array of strings if isNames=true'; 0065 dispEM(EM); 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 [~, 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=1:numel(model.rev); 0114 padding=num2cell(padding)'; 0115 padding=cellfun(@num2str,padding,'uniformoutput',false); 0116 model.rxns=padding; 0117 model.rxnNames=padding; 0118 model.eccodes=padding; 0119 model.rxnMiriams=padding; 0120 model.grRules=padding; 0121 if isfield(model,'genes') 0122 model.rxnGeneMat=sparse(numel(model.rev),numel(model.genes)); 0123 end 0124 model.subSystems=padding; 0125 model.rxnFrom=padding; 0126 model.rxnComps=ones(numel(model.rev),1); 0127 model.rxnNotes=padding; 0128 model.rxnReferences=padding; 0129 model.rxnConfidenceScores=NaN(numel(model.rev),1); 0130 0131 sol=solveLP(model,1); 0132 if any(sol.x) 0133 %It could be that several metabolites were consumed in order to get the 0134 %best solution. The setdiff is to avoid including the last fake 0135 %metabolite 0136 I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1)); 0137 0138 if any(I) %This should always be true 0139 %Change the coefficients so that only the first is consumed. This 0140 %is not always possible, but it is tested for since it it results 0141 %in more easily interpretable results 0142 0143 oldS=model.S; 0144 foundSingle=false; 0145 %Test if any of the metabolites could be consumed on their own 0146 for i=1:numel(I) 0147 model.S=oldS; 0148 J=nRxns+1:numel(model.lb)-1; 0149 J(I(i))=[]; 0150 model.S(:,J)=0; 0151 sol=solveLP(model); 0152 if any(sol.x) 0153 foundSingle=true; 0154 break; 0155 end 0156 end 0157 %This means that there was no metabolite which could be consumed on 0158 %its own. Then let all the consumeable metabolites be consumed. 0159 if foundSingle==false 0160 model.S=oldS; 0161 end 0162 if minNrFluxes==true 0163 %Has to add names for the rxns to prevent an error in 0164 %minNrFluxes 0165 model.rxns=cell(numel(model.lb),1); 0166 model.rxns(:)={'TEMP'}; 0167 model.mets=cell(size(model.b,1),1); 0168 model.mets(:)={'TEMP'}; 0169 sol=solveLP(model,3,params); 0170 else 0171 sol=solveLP(model,1); 0172 end 0173 solution=sol.x(1:nRxns); 0174 metabolite=find(sol.x(nRxns+1:end-1)>0.1); 0175 end 0176 end 0177 end