FSEOF

 FSEOF: implements the algorithm of Flux Scanning based on Enforced Objective Flux.

   model           a model structure
   biomassRxn      string with reaction ID of the biomass formation or
                   growth reaction
   targetRxn       string with reaction ID of target reaction
   iterations      double indicating number of iterations (opt, default 10)
   coefficient     double indicating ratio of optimal target reaction
                   flux, must be less than 1 (opt, default 0.9)
   outputFile      string with output filename (opt, default prints to
                   command window)

   targets         structure with target identifying information for each reaction
       logical     logical array indicating FSEOF identified target reaction
       slope       double array with FSEOF calculated slopes for each reaction

0001 function targets=FSEOF(model,biomassRxn,targetRxn,iterations,coefficient,outputFile)
0002 % FSEOF: implements the algorithm of Flux Scanning based on Enforced Objective Flux.
0003 %
0004 %   model           a model structure
0005 %   biomassRxn      string with reaction ID of the biomass formation or
0006 %                   growth reaction
0007 %   targetRxn       string with reaction ID of target reaction
0008 %   iterations      double indicating number of iterations (opt, default 10)
0009 %   coefficient     double indicating ratio of optimal target reaction
0010 %                   flux, must be less than 1 (opt, default 0.9)
0011 %   outputFile      string with output filename (opt, default prints to
0012 %                   command window)
0013 %
0014 %   targets         structure with target identifying information for each reaction
0015 %       logical     logical array indicating FSEOF identified target reaction
0016 %       slope       double array with FSEOF calculated slopes for each reaction
0017 
0018 %OUTPUTS
0019 %   This function writes an tab-delimited file or prints to command window. If an output
0020 %   has been specified (targets), it will also generate a structure indicating for
0021 %   each reaction whether it is identified by FSEOF as a target and the slope of the
0022 %   reaction when switching from biomass formation to product formation.
0023 %
0024 %   Usage: targets=FSEOF(model,biomassRxn,targetRxn,iterations,coefficient,outputFile)
0025 
0026 biomassRxn=char(biomassRxn);
0027 targetRxn=char(targetRxn);
0028 
0029 if nargin<4
0030     iterations=10;
0031     coefficient=0.9;
0032 end
0033 
0034 if nargin <5
0035     coefficient=0.9;
0036 end
0037 
0038 if nargin == 6
0039     output=1;
0040 else
0041     output=0;
0042 end
0043 
0044 %Find out the maximum theoretical yield of target reaction
0045 model=setParam(model,'obj',targetRxn,1);
0046 sol=solveLP(model,1);
0047 targetMax=abs(sol.f*coefficient);   % 90 percent of the theoretical yield
0048 
0049 model=setParam(model,'obj',biomassRxn,1);
0050 
0051 fseof.results=zeros(length(model.rxns),iterations);
0052 fseof.target=zeros(length(model.rxns),1);
0053 rxnDirection=zeros(length(model.rxns),1);
0054 
0055 %Enforce objective flux iteratively
0056 for i=1:iterations
0057     n=i*targetMax/iterations;
0058     model=setParam(model,'lb',targetRxn,n);
0059     
0060     sol=solveLP(model,1);
0061     
0062     fseof.results(:,i)=sol.x;
0063     
0064     %Loop through all fluxes and identify the ones that increased upon the
0065     %enforced objective flux
0066     for j=1:length(fseof.results)
0067         if fseof.results(j,1) > 0   %Check the positive fluxes
0068             
0069             if i == 1   %The initial round
0070                 rxnDirection(j,1)=1;
0071                 fseof.target(j,1)=1;
0072             else
0073                 
0074                 if (fseof.results(j,i) > fseof.results(j,i-1)) & fseof.target(j,1)
0075                     fseof.target(j,1)=1;
0076                 else
0077                     fseof.target(j,1)=0;
0078                 end
0079             end
0080             
0081         elseif fseof.results(j,1) < 0 %Check the negative fluxes
0082             
0083             if i == 1   %The initial round
0084                 rxnDirection(j,1)=-1;
0085                 fseof.target(j,1)=1;
0086             else
0087                 if (fseof.results(j,i) < fseof.results(j,i-1)) & fseof.target(j,1)
0088                     fseof.target(j,1)=1;
0089                 else
0090                     fseof.target(j,1)=0;
0091                 end
0092             end
0093             
0094         end
0095         
0096     end
0097 end
0098 
0099 %Generating output
0100 formatSpec='%s\t%s\t%s\t%s\t%s\t%s\t%s\n';
0101 if output == 1    %Output to a file
0102     outputFile=char(outputFile);
0103     fid=fopen(outputFile,'w');
0104     fprintf(fid,formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystems','Direction','Gr Rule');
0105 else              %Output to screen
0106     fprintf(formatSpec,'Slope','rowID','Enzyme ID','Enzyme Name','Subsystems','Direction','Gr Rule');
0107 end
0108 
0109 for num=1:length(fseof.target)
0110     if fseof.target(num,1) == 1
0111         A0=num2str(abs(fseof.results(num,iterations)-fseof.results(num,1))/abs(targetMax-targetMax/iterations)); %Slope calculation
0112         A1=num2str(num);                                  %row ID
0113         A2=char(model.rxns(num));                         %enzyme ID
0114         A3=char(model.rxnNames(num));                     %enzyme Name
0115         if isfield(model,'subSystems') && ~isempty(model.subSystems{num});
0116             A4=char(strjoin(model.subSystems{num,1},';'));                   %Subsystems
0117         else
0118             A4='';
0119         end
0120         A5=num2str(model.rev(num)*rxnDirection(num,1));   %reaction Dirction
0121         A6=char(model.grRules(num));                      %Gr Rule
0122         if output == 1    %Output to a file
0123             fprintf(fid,formatSpec,A0,A1,A2,A3,A4,A5,A6);
0124         else              %Output screen
0125             fprintf(formatSpec,A0,A1,A2,A3,A4,A5,A6);
0126         end
0127     end
0128 end
0129 
0130 if output == 1    %Output to a file
0131     fclose(fid);
0132 end
0133 
0134 if nargout == 1
0135     targets.logical=logical(fseof.target);
0136     targets.slope=abs(fseof.results(:,iterations)-fseof.results(:,1))/abs(targetMax-targetMax/iterations); %Slope calculation
0137 end
0138 end