getEssentialRxns Calculate the essential reactions for a model to be solvable model a model structure ignoreRxns cell array of reaction IDs which should not be checked (opt, default {}) essentialRxns cell array with the IDs of the essential reactions essentialRxnsIndexes vector with the indexes of the essential reactions Essential reactions are those which, when constrained to 0, result in an infeasible problem. Usage: [essentialRxns, essentialRxnsIndexes]=getEssentialRxns(model,ignoreRxns) Rasmus Agren, 2013-05-15
0001 function [essentialRxns, essentialRxnsIndexes]=getEssentialRxns(model,ignoreRxns) 0002 % getEssentialRxns 0003 % Calculate the essential reactions for a model to be solvable 0004 % 0005 % model a model structure 0006 % ignoreRxns cell array of reaction IDs which should not be 0007 % checked (opt, default {}) 0008 % 0009 % essentialRxns cell array with the IDs of the essential reactions 0010 % essentialRxnsIndexes vector with the indexes of the essential reactions 0011 % 0012 % Essential reactions are those which, when constrained to 0, result in an 0013 % infeasible problem. 0014 % 0015 % Usage: [essentialRxns, essentialRxnsIndexes]=getEssentialRxns(model,ignoreRxns) 0016 % 0017 % Rasmus Agren, 2013-05-15 0018 % 0019 0020 if nargin<2 0021 ignoreRxns={}; 0022 end 0023 0024 %Too make sure that it doesn't try to optimize for something 0025 model.c=zeros(numel(model.rxns),1); 0026 0027 %First check that the problem is solvable 0028 [sol hsSolOut]=solveLP(model,1); 0029 0030 if sol.stat~=1 || isempty(sol.x) 0031 throw(MException('','No feasible solution to the full model')); 0032 end 0033 0034 %Check which reactions have flux. Only those can be essential. This 0035 %is not the smallest list of reactions, but it's a fast way 0036 rxnsToCheck=setdiff(model.rxns(abs(sol.x)>10^-8),ignoreRxns); 0037 nToCheck=numel(rxnsToCheck); 0038 minimize=true; 0039 while 1 0040 if minimize==true 0041 sol=solveLP(setParam(model,'obj',rxnsToCheck,-1)); 0042 else 0043 sol=solveLP(setParam(model,'obj',rxnsToCheck,1)); 0044 end 0045 rxnsToCheck=intersect(rxnsToCheck,model.rxns(abs(sol.x)>10^-8)); 0046 if numel(rxnsToCheck)>=nToCheck 0047 if minimize==true 0048 minimize=false; 0049 else 0050 break; 0051 end 0052 else 0053 nToCheck=numel(rxnsToCheck); 0054 end 0055 end 0056 0057 essentialRxns={}; 0058 for i=1:numel(rxnsToCheck) 0059 sol=solveLP(setParam(model,'eq',rxnsToCheck(i),0),0,[],hsSolOut); 0060 if sol.stat~=1 || isempty(sol.x) 0061 essentialRxns=[essentialRxns;rxnsToCheck(i)]; 0062 end 0063 end 0064 0065 [crap essentialRxnsIndexes]=ismember(essentialRxns,model.rxns);