Home > RAVEN > getEssentialRxns.m

getEssentialRxns

PURPOSE ^

getEssentialRxns

SYNOPSIS ^

function [essentialRxns, essentialRxnsIndexes]=getEssentialRxns(model,ignoreRxns)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005