0001 function [biomassValues,targetValues] = multiProductionEnvelope(model,deletions,biomassRxn,geneDelFlag,nPts,plotAllFlag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 if (nargin < 2)
0026 deletions = {};
0027 end
0028 if (nargin < 3)
0029
0030 biomassRxn = model.rxns(model.c==1);
0031 end
0032 if (nargin < 4)
0033
0034 geneDelFlag = false;
0035 end
0036 if (nargin < 5)
0037 nPts = 20;
0038 end
0039 if (nargin < 6)
0040 plotAllFlag = false;
0041 end
0042
0043
0044
0045 if (length(deletions) > 0)
0046 if (geneDelFlag)
0047 model = deleteModelGenes(model,deletions);
0048 else
0049 model = changeRxnBounds(model,deletions,zeros(size(deletions)),'b');
0050 end
0051 end
0052
0053
0054 excRxns = model.rxns(findExcRxns(model,false,false));
0055 CRxns = findCarbonRxns(model,1);
0056 CExcRxns = intersect(excRxns,CRxns);
0057 substrateIDs = find(model.lb(findRxnIDs(model,CExcRxns))<0);
0058
0059 for i = 1:length(substrateIDs)
0060 j = substrateIDs(i);
0061 if j == 1
0062 CExcRxns = CExcRxns(2:length(CExcRxns));
0063 elseif j == length(CExcRxns)
0064 CExcRxns = CExcRxns(1:length(CExcRxns)-1);
0065 else
0066 CExcRxns = cat(1,CExcRxns(1:j-i),CExcRxns(j-i+2:length(CExcRxns)));
0067 end
0068 end
0069
0070
0071
0072 model = changeObjective(model,biomassRxn,1);
0073 solMax = optimizeCbModel(model,'max');
0074 solMin = optimizeCbModel(model,'min');
0075
0076
0077 biomassValues = linspace(solMin.f,solMax.f,nPts);
0078
0079 plottedRxns = [];
0080 targetUpperBound = [];
0081 targetLowerBound = [];
0082
0083 for i = 1:length(CExcRxns)
0084 model = changeObjective(model,CExcRxns(i),1);
0085 model2 = changeRxnBounds(model,biomassRxn,max(biomassValues),'b');
0086 fbasol2 = optimizeCbModel(model2,'max');
0087 maxRate = fbasol2.f;
0088 if (plotAllFlag)||(maxRate > 0.5)
0089 plottedRxns = [plottedRxns,i];
0090 for j = 1:length(biomassValues)
0091 model = changeRxnBounds(model,biomassRxn,biomassValues(j),'b');
0092 sol = optimizeCbModel(model,'max');
0093 if (sol.stat > 0)
0094 targetUpperBound(i,j) = sol.f;
0095 else
0096 targetUpperBound(i,j) = NaN;
0097 end
0098 sol = optimizeCbModel(model,'min');
0099 if (sol.stat > 0)
0100 targetLowerBound(i,j) = sol.f;
0101 else
0102 targetLowerBound(i,j) = NaN;
0103 end
0104 end
0105 end
0106 end
0107
0108
0109 colors = {'b','g','r','c','m','y','k'};
0110 for i = 1:length(plottedRxns)
0111 plot([biomassValues fliplr(biomassValues)],[targetUpperBound(plottedRxns(i),:) fliplr(targetLowerBound(plottedRxns(i),:))],colors{mod(i-1,length(colors))+1},'LineWidth',2)
0112 axis tight;
0113 hold on;
0114 end
0115 hold off;
0116 legend(CExcRxns(plottedRxns));
0117 legend off;
0118 ylabel('Production Rate (mmol/gDW h)');
0119 xlabel('Growth Rate (1/h)');
0120 plottools, plotbrowser('on'), figurepalette('hide'), propertyeditor('off');
0121
0122 biomassValues = biomassValues';
0123 targetValues = [targetLowerBound' targetUpperBound'];
0124
0125