runRobustnessAnalysis

 runRobustnessAnalysis
   Performs robustness analysis for a reaction of interest and an objective
   of interest. Modified from the COBRA robustnessAnalysis function.

 Input:
   model           a model structure
   controlRxn      reaction of interest whose value is to be controlled
   nPoints         number of points to show on plot (optional, default 20)
   objRxn          reaction identifier of objective to be maximized (optional,
                   default it uses the objective defined in the model)
   plotRedCost     logical whether reduced cost should also be plotted
                   (optional, default false)

 Output:
   controlFlux     flux values of the reaction of interest, ranging from
                   its minimum to its maximum value
   objFlux         optimal values of objective reaction at each control
                   reaction flux value

 Modified from COBRA Toolbox robustnessAnalysis.m

 Usage: runRobustnessAnalysis(model, controlRxn, nPoints, objRxn)

0001 function [controlFlux, objFlux] = runRobustnessAnalysis(model, controlRxn, nPoints, objRxn, plotRedCost)
0002 % runRobustnessAnalysis
0003 %   Performs robustness analysis for a reaction of interest and an objective
0004 %   of interest. Modified from the COBRA robustnessAnalysis function.
0005 %
0006 % Input:
0007 %   model           a model structure
0008 %   controlRxn      reaction of interest whose value is to be controlled
0009 %   nPoints         number of points to show on plot (optional, default 20)
0010 %   objRxn          reaction identifier of objective to be maximized (optional,
0011 %                   default it uses the objective defined in the model)
0012 %   plotRedCost     logical whether reduced cost should also be plotted
0013 %                   (optional, default false)
0014 %
0015 % Output:
0016 %   controlFlux     flux values of the reaction of interest, ranging from
0017 %                   its minimum to its maximum value
0018 %   objFlux         optimal values of objective reaction at each control
0019 %                   reaction flux value
0020 %
0021 % Modified from COBRA Toolbox robustnessAnalysis.m
0022 %
0023 % Usage: runRobustnessAnalysis(model, controlRxn, nPoints, objRxn)
0024 
0025 if nargin < 3
0026     nPoints = 20;
0027 end
0028 if nargin < 4
0029     baseModel = model;
0030 else
0031     baseModel = setParam(model,'obj',objRxn,1);
0032 end
0033 if nargin < 5
0034     plotRedCost = false;
0035 end
0036 
0037 if any(ismember(model.rxns,controlRxn))
0038     controlRxnIdx = getIndexes(model,controlRxn,'rxns');
0039     tmpModel = setParam(model,'obj',controlRxnIdx,1);
0040     solMax = solveLP(tmpModel);
0041     solMax = solMax.x(logical(tmpModel.c));
0042     tmpModel.c = -tmpModel.c;
0043     solMin = solveLP(tmpModel);
0044     solMin = solMin.x(logical(tmpModel.c));
0045 else
0046     error('Control reaction does not exist!');
0047 end
0048 
0049 objFlux = zeros(nPoints,1);
0050 redCost = zeros(nPoints,1);
0051 controlFlux = linspace(solMin,solMax,nPoints)';
0052 
0053 
0054 PB = ProgressBar2(length(controlFlux),'Running robustness analysis','cli');
0055 for i=1:length(controlFlux)
0056     modelControlled = setParam(baseModel,'eq',controlRxnIdx,controlFlux(i));
0057     solControlled = solveLP(modelControlled);
0058     objFlux(i) = solControlled.x(logical(modelControlled.c));
0059     redCost(i) = solControlled.rCost(controlRxnIdx);
0060     count(PB);
0061 end
0062 
0063 if plotRedCost
0064     yyaxis right
0065     plot(controlFlux,redCost)
0066     ylabel([strrep(controlRxn,'_','-') ' reduced cost']);
0067     yyaxis left
0068 end
0069 plot(controlFlux,objFlux)
0070 xlabel(strrep(controlRxn,'_','-'));
0071 ylabel('Objective');
0072 end