doubleProductionEnvelope plots maximum growth rate as a function of the output of two specified products [x1,x2,y] = doubleProductionEnvelope(model,deletions,prod1,prod2,biomassRxn,geneDelFlag,nPts) INPUTS model COBRA model structure deletions The reactions or genes to knockout of the model prod1 One of the two products to investigate prod2 The other product to investigate OPTIONAL INPUTS biomassRxn The biomass objective function rxn name (Default = 'biomass_SC4_bal') geneDelFlag Perform gene and not reaction deletions (Default = false) nPts Number of points to plot for each product (Default = 20) OUTPUTS x1 The range of rates plotted for prod1 x2 The range of rates plotted for prod2 y The plotted growth rates at each (x1,x2) Jeff Orth 9/12/07
0001 function [x1,x2,y] = doubleProductionEnvelope(model,deletions,prod1,prod2,biomassRxn,geneDelFlag,nPts) 0002 %doubleProductionEnvelope plots maximum growth rate as a function of the 0003 %output of two specified products 0004 % 0005 % [x1,x2,y] = doubleProductionEnvelope(model,deletions,prod1,prod2,biomassRxn,geneDelFlag,nPts) 0006 % 0007 %INPUTS 0008 % model COBRA model structure 0009 % deletions The reactions or genes to knockout of the model 0010 % prod1 One of the two products to investigate 0011 % prod2 The other product to investigate 0012 % 0013 %OPTIONAL INPUTS 0014 % biomassRxn The biomass objective function rxn name 0015 % (Default = 'biomass_SC4_bal') 0016 % geneDelFlag Perform gene and not reaction deletions 0017 % (Default = false) 0018 % nPts Number of points to plot for each product 0019 % (Default = 20) 0020 % 0021 %OUTPUTS 0022 % x1 The range of rates plotted for prod1 0023 % x2 The range of rates plotted for prod2 0024 % y The plotted growth rates at each (x1,x2) 0025 % 0026 % Jeff Orth 9/12/07 0027 0028 if (nargin < 5) 0029 biomassRxn = 'biomass_SC4_bal'; 0030 end 0031 if (nargin < 6) 0032 geneDelFlag = false; 0033 end 0034 if (nargin < 7) 0035 nPts = 20; 0036 end 0037 0038 % Create model with deletions 0039 if (length(deletions) > 0) 0040 if (geneDelFlag) 0041 model = deleteModelGenes(model,deletions); 0042 else 0043 model = changeRxnBounds(model,deletions,zeros(size(deletions)),'b'); 0044 end 0045 end 0046 0047 % find range for prod1 0048 modelFixed1 = changeRxnBounds(model,biomassRxn,0,'b'); 0049 modelFixed1 = changeObjective(modelFixed1,prod1); 0050 fbasol1 = optimizeCbModel(modelFixed1,'max'); 0051 max1 = fbasol1.f; 0052 0053 % find range for prod2 0054 modelFixed2 = changeRxnBounds(model,biomassRxn,0,'b'); 0055 modelFixed2 = changeObjective(modelFixed2,prod2); 0056 fbasol2 = optimizeCbModel(modelFixed2,'max'); 0057 max2 = fbasol2.f; 0058 0059 % vary both ranges, find max growth rate 0060 x1 = linspace(0,max1,nPts); 0061 x2 = linspace(0,max2,nPts); 0062 y = zeros(nPts); 0063 0064 for i = 1:nPts 0065 for j = 1:nPts 0066 prod1val = x1(i); 0067 prod2val = x2(j); 0068 modelY = changeRxnBounds(model,prod1,prod1val,'b'); 0069 modelY = changeRxnBounds(modelY,prod2,prod2val,'b'); 0070 fbasol = optimizeCbModel(modelY,'max'); 0071 y(j,i) = fbasol.f; %no, this isn't a mistake (probably) 0072 end 0073 end 0074 0075 % plot 0076 surf(x1,x2,y); 0077 alpha(.4); 0078 xlabel([strrep(prod1,'_','\_'),' (mmol/gDW h)']); 0079 ylabel([strrep(prod2,'_','\_'),' (mmol/gDW h)']); 0080 zlabel('Growth Rate (1/h)'); 0081 0082