multiProductionEnvelope

PURPOSE ^

multiProductionEnvelope Calculates the byproduct secretion envelopes for

SYNOPSIS ^

function [biomassValues,targetValues] = multiProductionEnvelope(model,deletions,biomassRxn,geneDelFlag,nPts,plotAllFlag)

DESCRIPTION ^

multiProductionEnvelope Calculates the byproduct secretion envelopes for
every product (excreted metabolites with 1 or more Carbons)

 [biomassValues,targetValues] = multiProductionEnvelope(model,deletions,biomassRxn,geneDelFlag,nPts,plotAllFlag)

INPUT
 model         COBRA model structure

OPTIONAL INPUT
 deletions     List of reaction or gene deletions (empty if wild type)
               (Default = {})
 biomassRxn    Biomass rxn name (Default = whatever is defined in model)
 geneDelFlag   Perform gene and not reaction deletions (Default = false)
 nPts          Number of points in the plot (Default = 20)
 plotAllFlag   plot all envelopes, even ones that are not growth coupled
               (Default = false)

OUTPUT
 biomassValues Biomass values for plotting
 targetValues  Target upper and lower bounds for plotting

 Jeff Orth 8/16/07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [biomassValues,targetValues] = multiProductionEnvelope(model,deletions,biomassRxn,geneDelFlag,nPts,plotAllFlag)
0002 %multiProductionEnvelope Calculates the byproduct secretion envelopes for
0003 %every product (excreted metabolites with 1 or more Carbons)
0004 %
0005 % [biomassValues,targetValues] = multiProductionEnvelope(model,deletions,biomassRxn,geneDelFlag,nPts,plotAllFlag)
0006 %
0007 %INPUT
0008 % model         COBRA model structure
0009 %
0010 %OPTIONAL INPUT
0011 % deletions     List of reaction or gene deletions (empty if wild type)
0012 %               (Default = {})
0013 % biomassRxn    Biomass rxn name (Default = whatever is defined in model)
0014 % geneDelFlag   Perform gene and not reaction deletions (Default = false)
0015 % nPts          Number of points in the plot (Default = 20)
0016 % plotAllFlag   plot all envelopes, even ones that are not growth coupled
0017 %               (Default = false)
0018 %
0019 %OUTPUT
0020 % biomassValues Biomass values for plotting
0021 % targetValues  Target upper and lower bounds for plotting
0022 %
0023 % Jeff Orth 8/16/07
0024 
0025 if (nargin < 2)
0026     deletions = {};
0027 end
0028 if (nargin < 3)
0029     % Biomass flux
0030     biomassRxn = model.rxns(model.c==1);
0031 end
0032 if (nargin < 4)
0033     % Gene or rxn deletions
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 % Create model with deletions
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 %get all C exchange reactions
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 %remove the substrate reactions
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 % Run FBA to get upper bound for biomass
0072 model = changeObjective(model,biomassRxn,1);
0073 solMax = optimizeCbModel(model,'max');
0074 solMin = optimizeCbModel(model,'min');
0075 
0076 % Create biomass range vector
0077 biomassValues = linspace(solMin.f,solMax.f,nPts);
0078 
0079 plottedRxns = [];
0080 targetUpperBound = [];
0081 targetLowerBound = [];
0082 % Max/min for target production
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; %find max production at max growth rate
0088     if (plotAllFlag)||(maxRate > 0.5) %only plot growth coupled solutions
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 % Plot results
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

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003