runProductionEnvelope

 runProductionEnvelope
   Calculates the byproduct secretion envelope

 Input:
   model            a model structure
   targetRxn        identifier of target metabolite production reaction
   biomassRxn       identifier of biomass reaction
   nPts             number of points in the plot (opt, default 20)

 Output:
   biomassValues    Biomass values for plotting
   targetValues     Target upper and lower bounds for plotting

 Modified from COBRA Toolbox productionEnvelope.m

 Usage: [biomassValues, targetValues] = runProductionEnvelope(model, targetRxn, biomassRxn, nPts)

0001 function [biomassValues, targetValues] = runProductionEnvelope(model, targetRxn, biomassRxn, nPts)
0002 % runProductionEnvelope
0003 %   Calculates the byproduct secretion envelope
0004 %
0005 % Input:
0006 %   model            a model structure
0007 %   targetRxn        identifier of target metabolite production reaction
0008 %   biomassRxn       identifier of biomass reaction
0009 %   nPts             number of points in the plot (opt, default 20)
0010 %
0011 % Output:
0012 %   biomassValues    Biomass values for plotting
0013 %   targetValues     Target upper and lower bounds for plotting
0014 %
0015 % Modified from COBRA Toolbox productionEnvelope.m
0016 %
0017 % Usage: [biomassValues, targetValues] = runProductionEnvelope(model, targetRxn, biomassRxn, nPts)
0018 
0019 if nargin < 4
0020     nPts = 20;
0021 end
0022 
0023 % Run FBA to get upper bound for biomass
0024 model   = setParam(model,'obj',biomassRxn,1);
0025 [solMax,hsSol]  = solveLP(model);
0026 solMax  = solMax.x(logical(model.c));
0027 model.c = -model.c;
0028 solMin  = solveLP(model,0,[],hsSol);
0029 solMin  = solMin.x(logical(model.c));
0030 
0031 % Create biomass range vector
0032 biomassValues = linspace(solMin,solMax,nPts);
0033 targetUpperBound = nan(1,numel(biomassValues));
0034 targetLowerBound = nan(1,numel(biomassValues));
0035 
0036 PB = ProgressBar2(length(biomassValues),'Creating production envelope','cli');
0037 % Max/min for target production
0038 model = setParam(model,'obj',targetRxn,1);
0039 for i = 1:length(biomassValues)
0040     model1 = setParam(model,'eq',biomassRxn,biomassValues(i));
0041     sol = solveLP(model1,0,[],hsSol);
0042     if (sol.stat > 0)
0043         targetUpperBound(i) = sol.x(logical(model.c));
0044     else
0045         targetUpperBound(i) = NaN;
0046     end
0047     model1.c = -model1.c;
0048     sol = solveLP(model1,0,[],hsSol);
0049     if (sol.stat > 0)
0050         targetLowerBound(i) = sol.x(logical(model1.c));
0051     else
0052         targetLowerBound(i) = NaN;
0053     end
0054     count(PB);
0055 end
0056 
0057 % Plot results
0058 plot([biomassValues fliplr(biomassValues)],[targetUpperBound fliplr(targetLowerBound)],'blue','LineWidth',2);
0059 axis tight;
0060 ylabel([strrep(targetRxn,'_','-') ' (mmol/gDW h)']);
0061 xlabel('Growth rate (1/h)');
0062 targetValues = [targetLowerBound; targetUpperBound];
0063 end