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-04-21
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-04-21 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 no solution could be found, it could possible have to do with numerical 0031 %issues from minimizing for the sum of fluxes. Try a normal solution as 0032 %well. This should give the same output, just take a little longer. 0033 if sol.stat~=1 || isempty(sol.x) 0034 [sol hsSolOut]=solveLP(model); 0035 end 0036 if sol.stat~=1 || isempty(sol.x) 0037 throw(MException('','No feasible solution to the full model')); 0038 end 0039 0040 %Check which reactions have flux. Only those can be essential. This 0041 %is not the smallest list of reactions, but it's a fast way 0042 rxnsToCheck=setdiff(model.rxns(abs(sol.x)>10^-8),ignoreRxns); 0043 nToCheck=numel(rxnsToCheck); 0044 minimize=true; 0045 while 1 0046 if minimize==true 0047 sol=solveLP(setParam(model,'obj',rxnsToCheck,-1)); 0048 else 0049 sol=solveLP(setParam(model,'obj',rxnsToCheck,1)); 0050 end 0051 rxnsToCheck=intersect(rxnsToCheck,model.rxns(abs(sol.x)>10^-8)); 0052 if numel(rxnsToCheck)>=nToCheck 0053 if minimize==true 0054 minimize=false; 0055 else 0056 break; 0057 end 0058 else 0059 nToCheck=numel(rxnsToCheck); 0060 end 0061 end 0062 0063 essentialRxns={}; 0064 for i=1:numel(rxnsToCheck) 0065 sol=solveLP(setParam(model,'eq',rxnsToCheck(i),0),0,[],hsSolOut); 0066 if sol.stat~=1 || isempty(sol.x) 0067 essentialRxns=[essentialRxns;rxnsToCheck(i)]; 0068 end 0069 end 0070 0071 [crap essentialRxnsIndexes]=ismember(essentialRxns,model.rxns);