checkRxn

 checkRxn
   Checks which reactants in a reaction that can be synthesized and which
   products that can be consumed. This is primarily for debugging
   reactions which cannot have flux

   model       a model structure
   rxn         the id of one reaction to check
   cutoff      minimal flux for successful production/consumption (optional,
               default 10^-7)
   revDir      true if the reaction should be reversed (optional, default
               false)
   printReport print a report (optional, default true)

   report
       reactants   array with reactant indexes
       canMake     boolean array, true if the corresponding reactant can be
                   synthesized
       products    array with product indexes
       canConsume  boolean array, true if the corresponding reactant can
                   be consumed

 Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport)

0001 function report=checkRxn(model,rxn,cutoff,revDir,printReport)
0002 % checkRxn
0003 %   Checks which reactants in a reaction that can be synthesized and which
0004 %   products that can be consumed. This is primarily for debugging
0005 %   reactions which cannot have flux
0006 %
0007 %   model       a model structure
0008 %   rxn         the id of one reaction to check
0009 %   cutoff      minimal flux for successful production/consumption (optional,
0010 %               default 10^-7)
0011 %   revDir      true if the reaction should be reversed (optional, default
0012 %               false)
0013 %   printReport print a report (optional, default true)
0014 %
0015 %   report
0016 %       reactants   array with reactant indexes
0017 %       canMake     boolean array, true if the corresponding reactant can be
0018 %                   synthesized
0019 %       products    array with product indexes
0020 %       canConsume  boolean array, true if the corresponding reactant can
0021 %                   be consumed
0022 %
0023 % Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport)
0024 
0025 rxn=char(rxn);
0026 if nargin<3
0027     cutoff=10^-7;
0028 end
0029 if nargin<4
0030     revDir=false;
0031 end
0032 if isempty(cutoff)
0033     cutoff=10^-7;
0034 end
0035 if nargin<5
0036     printReport=true;
0037 end
0038 
0039 [I, rxnID]=ismember(rxn,model.rxns);
0040 
0041 if ~I
0042     EM='Reaction ID not found';
0043     dispEM(EM);
0044 end
0045 
0046 if revDir==false
0047     report.reactants=find(model.S(:,rxnID)<0);
0048     report.products=find(model.S(:,rxnID)>0);
0049 else
0050     report.reactants=find(model.S(:,rxnID)>0);
0051     report.products=find(model.S(:,rxnID)<0);
0052 end
0053 report.canMake=false(numel(report.reactants),1);
0054 report.canConsume=false(numel(report.products),1);
0055 
0056 %Remove this field as it would give an annoying note otherwise
0057 if isfield(model,'rxnComps')
0058     model=rmfield(model,'rxnComps');
0059 end
0060 
0061 %There are several ways to do this. Here I choose to add the reactions one
0062 %by one and checking their bounds. This might not be optimal
0063 for i=1:numel(report.reactants)
0064     [tempModel, testRxn]=addExchangeRxns(model,'out',report.reactants(i));
0065     tempModel=setParam(tempModel,'obj',testRxn,1);
0066     sol=solveLP(tempModel);
0067     if sol.f>cutoff
0068         report.canMake(i)=true;
0069     else
0070         if printReport==true
0071             fprintf(['Failed to make ' model.metNames{report.reactants(i)} '[' model.comps{model.metComps(report.reactants(i))} ']\n']);
0072         end
0073     end
0074 end
0075 
0076 for i=1:numel(report.products)
0077     [tempModel, testRxn]=addExchangeRxns(model,'in',report.products(i));
0078     tempModel=setParam(tempModel,'obj',testRxn,1);
0079     sol=solveLP(tempModel);
0080     if sol.f>cutoff
0081         report.canConsume(i)=true;
0082     else
0083         if printReport==true
0084             fprintf(['Failed to consume ' model.metNames{report.products(i)} '[' model.comps{model.metComps(report.products(i))} ']\n']);
0085         end
0086     end
0087 end
0088 end