haveFlux

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

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

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

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