checkProduction
Checks which metabolites that can be produced from a model using the
specified constraints.
model a model structure
checkNeededForProduction for each of the metabolites that could not
be produced, include an artificial
production reaction and calculate which new
metabolites that could be produced as en
effect of this (optional, default false)
excretionFromCompartments cell array with compartment ids from which
metabolites can be excreted (optional, default
model.comps)
printDetails print details to the screen (optional, default
true)
notProduced cell array with metabolites that could not
be produced
notProducedNames cell array with names and compartments for
metabolites that could not be produced
neededForProductionMat matrix where n x m is true if metabolite n
allows for production of metabolite m
minToConnect structure with the minimal number of
metabolites that need to be connected in
order to be able to produce all other
metabolites and which metabolites each of
them connects
model updated model structure with excretion
reactions added
The function is intended to be used to identify which metabolites must
be connected in order to have a fully connected network. It does so by
first identifying which metabolites could have a net production in the
network. Then it calculates which other metabolites must be able to
have net production in order to have production of all metabolites in
the network. So, if a network contains the equations A[external]->B,
C->D, and D->E it will identify that production of C will connect
the metabolites D and E.
Usage: [notProduced, notProducedNames,neededForProductionMat,minToConnect,model]=...
checkProduction(model,checkNeededForProduction,...
excretionFromCompartments,printDetails)0001 function [notProduced, notProducedNames, neededForProductionMat,minToConnect,model]=checkProduction(model,checkNeededForProduction,excretionFromCompartments,printDetails) 0002 % checkProduction 0003 % Checks which metabolites that can be produced from a model using the 0004 % specified constraints. 0005 % 0006 % model a model structure 0007 % checkNeededForProduction for each of the metabolites that could not 0008 % be produced, include an artificial 0009 % production reaction and calculate which new 0010 % metabolites that could be produced as en 0011 % effect of this (optional, default false) 0012 % excretionFromCompartments cell array with compartment ids from which 0013 % metabolites can be excreted (optional, default 0014 % model.comps) 0015 % printDetails print details to the screen (optional, default 0016 % true) 0017 % 0018 % notProduced cell array with metabolites that could not 0019 % be produced 0020 % notProducedNames cell array with names and compartments for 0021 % metabolites that could not be produced 0022 % neededForProductionMat matrix where n x m is true if metabolite n 0023 % allows for production of metabolite m 0024 % minToConnect structure with the minimal number of 0025 % metabolites that need to be connected in 0026 % order to be able to produce all other 0027 % metabolites and which metabolites each of 0028 % them connects 0029 % model updated model structure with excretion 0030 % reactions added 0031 % 0032 % The function is intended to be used to identify which metabolites must 0033 % be connected in order to have a fully connected network. It does so by 0034 % first identifying which metabolites could have a net production in the 0035 % network. Then it calculates which other metabolites must be able to 0036 % have net production in order to have production of all metabolites in 0037 % the network. So, if a network contains the equations A[external]->B, 0038 % C->D, and D->E it will identify that production of C will connect 0039 % the metabolites D and E. 0040 % 0041 % Usage: [notProduced, notProducedNames,neededForProductionMat,minToConnect,model]=... 0042 % checkProduction(model,checkNeededForProduction,... 0043 % excretionFromCompartments,printDetails) 0044 0045 if nargin<2 0046 checkNeededForProduction=false; 0047 end 0048 0049 if nargin<3 0050 excretionFromCompartments=model.comps; 0051 else 0052 excretionFromCompartments=convertCharArray(excretionFromCompartments); 0053 end 0054 0055 if nargin<4 0056 printDetails=true; 0057 end 0058 0059 %Add an exchange reaction for each metabolite in the allowed compartments 0060 %and see if it can carry a flux 0061 allowedMetIds=ismember(model.comps(model.metComps),excretionFromCompartments); 0062 allowedMetIndexes=find(allowedMetIds); 0063 [model, addedRxns]=addExchangeRxns(model,'out',allowedMetIndexes); 0064 0065 canProduce=haveFlux(model,10^-5,addedRxns); 0066 0067 notProduced=find(~canProduce); 0068 minToConnect={}; 0069 if checkNeededForProduction==true 0070 %For each of the metabolites that couldn't be produced allow uptake and 0071 %check which of the other metabolites that couldn't be produced that 0072 %can be produced 0073 neededForProductionMat=false(numel(notProduced)); 0074 for i=1:numel(notProduced) 0075 %Add uptake for this metabolite 0076 if i>1 0077 %Reset last iteration 0078 model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i-1))=model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i-1))*-1; 0079 end 0080 %Change the production reaction to an uptake reaction 0081 model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i))=model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i))*-1; 0082 0083 %Test which of the metabolites that couldn't be produced that can 0084 %be produced now 0085 neededForProductionMat(i,:)=haveFlux(model,10^-5,addedRxns(notProduced)); 0086 end 0087 %Calculate the smallest number of metabolites that must be connected to 0088 %make everything connected and return their names 0089 0090 %The algorithm is relatively straight forward. It finds the metabolite 0091 %that connects the most unconnected metabolites (iteratively), adds it 0092 %and removes the now connected metabolites until all are connected. 0093 %This is not guaranteed to find the global minimum 0094 neededForProdTemp=neededForProductionMat; 0095 while 1==1 0096 %Get all metabolites a metabolite makes connected 0097 totalConnected=false(size(neededForProdTemp)); 0098 for i=1:numel(notProduced) 0099 totalConnected(i,:)=neededForProdTemp(i,:); 0100 0101 lastIter=0; 0102 while 1==1 0103 [~, a]=find(neededForProdTemp(totalConnected(i,:),:)); 0104 totalConnected(i,a)=true; 0105 if numel(a)==lastIter 0106 break; %No more connections were possible 0107 else 0108 lastIter=numel(a); 0109 end 0110 end 0111 end 0112 [connections, mostConnected]=max(sum(totalConnected,2)); 0113 0114 if connections>0 0115 %Add the most connected metabolite to the list and remove all 0116 %metabolites that it's connected to 0117 metID=allowedMetIndexes(notProduced(mostConnected)); 0118 entry=[model.metNames{metID},'[',model.comps{model.metComps(metID)},'] (connects ' num2str(connections) ' metabolites)']; 0119 minToConnect=[minToConnect;entry]; 0120 neededForProdTemp(totalConnected(mostConnected,:),:)=false; 0121 neededForProdTemp(:,totalConnected(mostConnected,:))=false; 0122 else 0123 break; 0124 end 0125 end 0126 else 0127 neededForProductionMat=[]; 0128 end 0129 0130 notProducedNames=strcat(model.metNames(allowedMetIndexes(notProduced)),'[',model.comps(model.metComps(allowedMetIndexes(notProduced))),']'); 0131 0132 if printDetails==true 0133 fprintf('The following metabolites could not be produced:\n'); 0134 [notProducedNamesTemp,perm]=sort(notProducedNames); 0135 0136 if checkNeededForProduction==true 0137 neededForProdTemp=neededForProductionMat(:,perm); 0138 neededForProdTemp=neededForProdTemp(perm,:); 0139 fprintf('\tIf the production of a metabolite is dependent on some other metabolites then those are printed under the name\n\n'); 0140 end 0141 for i=1:numel(notProducedNamesTemp) 0142 fprintf([notProducedNamesTemp{i} '\n']); 0143 neededForProdTemp(i,i)=false; %Not neat to do this here. Prevent printing itself 0144 if checkNeededForProduction==true 0145 enablesProduction=find(neededForProdTemp(:,i)); 0146 if any(enablesProduction) 0147 for j=1:numel(enablesProduction) 0148 fprintf(['\t' notProducedNamesTemp{enablesProduction(j)} '\n']); 0149 end 0150 end 0151 end 0152 end 0153 end 0154 end