getEssentialRxns

 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)

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 if nargin<2
0018     ignoreRxns={};
0019 else
0020     ignoreRxns=convertCharArray(ignoreRxns);
0021 end
0022 
0023 %Too make sure that it doesn't try to optimize for something
0024 model.c=zeros(numel(model.rxns),1);
0025 
0026 %First check that the problem is solvable
0027 [sol, hsSolOut]=solveLP(model,1);
0028 
0029 if sol.stat==-1 || isempty(sol.x)
0030     EM='No feasible solution to the full model';
0031     dispEM(EM);
0032 end
0033 
0034 %Check which reactions have flux. Only those can be essential. This is not
0035 %the smallest list of reactions, but it's a fast way
0036 rxnsToCheck=setdiff(model.rxns(abs(sol.x)>10^-12),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 [~, essentialRxnsIndexes]=ismember(essentialRxns,model.rxns);
0066 end