haveFlux

 haveFlux
   Checks which reactions can carry a (positive or negative) flux. Is used
   as a faster version of getAllowedBounds if it is only interesting
   whether the reactions can carry a flux or not

 Input:
   model       a model structure
   cutOff      the flux value that a reaction has to carry to be
               identified as positive (optional, default 10^-8)
   rxns        either a cell array of IDs, a logical vector with the
               same number of elements as metabolites in the model, or a
               vector of indexes (optional, default model.rxns)

 Output:
   I           logical array with true if the corresponding reaction can
               carry a flux

 If a model has +/- Inf bounds then those are replaced with an arbitary
 large value of +/- 10000 prior to solving

 Usage: I = haveFlux(model, cutOff, rxns)

0001 function I=haveFlux(model,cutOff,rxns)
0002 % haveFlux
0003 %   Checks which reactions can carry a (positive or negative) flux. Is used
0004 %   as a faster version of getAllowedBounds if it is only interesting
0005 %   whether the reactions can carry a flux or not
0006 %
0007 % Input:
0008 %   model       a model structure
0009 %   cutOff      the flux value that a reaction has to carry to be
0010 %               identified as positive (optional, default 10^-8)
0011 %   rxns        either a cell array of IDs, a logical vector with the
0012 %               same number of elements as metabolites in the model, or a
0013 %               vector of indexes (optional, default model.rxns)
0014 %
0015 % Output:
0016 %   I           logical array with true if the corresponding reaction can
0017 %               carry a flux
0018 %
0019 % If a model has +/- Inf bounds then those are replaced with an arbitary
0020 % large value of +/- 10000 prior to solving
0021 %
0022 % Usage: I = haveFlux(model, cutOff, rxns)
0023 
0024 if nargin<2
0025     cutOff=10^-6;
0026 end
0027 if isempty(cutOff)
0028     cutOff=10^-6;
0029 end
0030 if nargin<3
0031     rxns=model.rxns;
0032 elseif ~islogical(rxns) && ~isnumeric(rxns)
0033     rxns=convertCharArray(rxns);
0034 end
0035 
0036 %This is since we're maximizing for the sum of fluxes, which isn't possible
0037 %when there are infinite bounds
0038 model.lb(model.lb==-inf)=-10000;
0039 model.ub(model.ub==inf)=10000;
0040 
0041 %Get the reaction IDs. A bit of an awkward way, but fine.
0042 indexes=getIndexes(model,rxns,'rxns');
0043 rxns=model.rxns(indexes);
0044 
0045 %First remove all dead-end reactions
0046 smallModel=simplifyModel(model,false,false,true,true);
0047 
0048 %Get the indexes of the reactions in the connected model
0049 indexes=getIndexes(smallModel,intersect(smallModel.rxns,rxns),'rxns');
0050 J=false(numel(indexes),1);
0051 mixIndexes=indexes(randperm(numel(indexes)));
0052 
0053 %Maximize for all fluxes first in order to get fewer rxns to test
0054 smallModel.c=ones(numel(smallModel.c),1);
0055 sol=solveLP(smallModel);
0056 if ~isempty(sol.x)
0057     J(abs(sol.x(mixIndexes))>cutOff)=true;
0058 end
0059 
0060 %Loop through and maximize then minimize each rxn if it does not already
0061 %have a flux
0062 Z=zeros(numel(smallModel.c),1);
0063 hsSolOut=[];
0064 for i=[1 -1]
0065     for j=1:numel(J)
0066         if J(j)==false
0067             %Only minimize for reversible reactions
0068             if i==1 || smallModel.rev(mixIndexes(j))~=0
0069                 smallModel.c=Z;
0070                 smallModel.c(mixIndexes(j))=i;
0071                 [sol, hsSolOut]=solveLP(smallModel,0,[],hsSolOut);
0072                 if any(sol.x)
0073                     J(abs(sol.x(mixIndexes))>cutOff)=true;
0074                 end
0075             end
0076         end
0077     end
0078 end
0079 
0080 I=ismember(rxns,smallModel.rxns(mixIndexes(J)));
0081 end