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 (optional, 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 (optional, 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,...
0018 %               targetRxn, biomassRxn, nPts)
0019 
0020 if nargin < 4
0021     nPts = 20;
0022 end
0023 
0024 % Run FBA to get upper bound for biomass
0025 model   = setParam(model,'obj',biomassRxn,1);
0026 [solMax,hsSol]  = solveLP(model);
0027 solMax  = solMax.x(logical(model.c));
0028 model.c = -model.c;
0029 solMin  = solveLP(model,0,[],hsSol);
0030 solMin  = solMin.x(logical(model.c));
0031 
0032 % Create biomass range vector
0033 biomassValues = linspace(solMin,solMax,nPts);
0034 targetUpperBound = nan(1,numel(biomassValues));
0035 targetLowerBound = nan(1,numel(biomassValues));
0036 
0037 PB = ProgressBar2(length(biomassValues),'Creating production envelope','cli');
0038 % Max/min for target production
0039 model = setParam(model,'obj',targetRxn,1);
0040 for i = 1:length(biomassValues)
0041     model1 = setParam(model,'eq',biomassRxn,biomassValues(i));
0042     sol = solveLP(model1,0,[],hsSol);
0043     if (sol.stat > 0)
0044         targetUpperBound(i) = sol.x(logical(model.c));
0045     else
0046         targetUpperBound(i) = NaN;
0047     end
0048     model1.c = -model1.c;
0049     sol = solveLP(model1,0,[],hsSol);
0050     if (sol.stat > 0)
0051         targetLowerBound(i) = sol.x(logical(model1.c));
0052     else
0053         targetLowerBound(i) = NaN;
0054     end
0055     count(PB);
0056 end
0057 
0058 % Plot results
0059 plot([biomassValues fliplr(biomassValues)],[targetUpperBound fliplr(targetLowerBound)],'blue','LineWidth',2);
0060 axis tight;
0061 ylabel([strrep(targetRxn,'_','-') ' (mmol/gDW h)']);
0062 xlabel('Growth rate (1/h)');
0063 targetValues = [targetLowerBound; targetUpperBound];
0064 end