multiProductionEnvelopeInorg

PURPOSE ^

multiProductionEnvelopeInorg calculates the byproduct secretion envelopes

SYNOPSIS ^

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

DESCRIPTION ^

multiProductionEnvelopeInorg calculates the byproduct secretion envelopes 
for every product, including inorganic compounds

 [biomassValues,targetValues] = multiProductionEnvelopeInorg(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 5/1/08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [biomassValues,targetValues] = multiProductionEnvelopeInorg(model,deletions,biomassRxn,geneDelFlag,nPts,plotAllFlag)
0002 %multiProductionEnvelopeInorg calculates the byproduct secretion envelopes
0003 %for every product, including inorganic compounds
0004 %
0005 % [biomassValues,targetValues] = multiProductionEnvelopeInorg(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 5/1/08
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 exchange reactions
0054 [selExc,selUpt] = findExcRxns(model,false,false);
0055 excRxns = [];
0056 for i = 1:length(model.rxns)
0057    if selExc(i) && ~selUpt(i)
0058        excRxns = [excRxns,model.rxns(i)];
0059    end
0060 end
0061 
0062 % Run FBA to get upper bound for biomass
0063 model = changeObjective(model,biomassRxn,1);
0064 solMax = optimizeCbModel(model,'max');
0065 solMin = optimizeCbModel(model,'min');
0066 
0067 % Create biomass range vector
0068 biomassValues = linspace(solMin.f,solMax.f,nPts);
0069 
0070 plottedRxns = [];
0071 targetUpperBound = [];
0072 targetLowerBound = [];
0073 % Max/min for target production
0074 for i = 1:length(excRxns)
0075     model = changeObjective(model,excRxns(i),1);
0076     model2 = changeRxnBounds(model,biomassRxn,max(biomassValues),'b');
0077     fbasol2 = optimizeCbModel(model2,'max');
0078     maxRate = fbasol2.f; %find max production at max growth rate
0079     if (plotAllFlag)||(maxRate > 0.5) %only plot growth coupled solutions
0080         plottedRxns = [plottedRxns,i];
0081         for j = 1:length(biomassValues)
0082             model = changeRxnBounds(model,biomassRxn,biomassValues(j),'b');
0083             sol = optimizeCbModel(model,'max');
0084             if (sol.stat > 0)
0085                 targetUpperBound(i,j) = sol.f;
0086             else
0087                 targetUpperBound(i,j) = NaN;
0088             end
0089             sol = optimizeCbModel(model,'min');    
0090             if (sol.stat > 0)
0091                 targetLowerBound(i,j) = sol.f;
0092             else
0093                 targetLowerBound(i,j) = NaN;
0094             end
0095         end
0096     end
0097 end
0098 
0099 % Plot results
0100 colors = {'b','g','r','c','m','y','k'};
0101 for i = 1:length(plottedRxns)
0102     plot([biomassValues fliplr(biomassValues)],[targetUpperBound(plottedRxns(i),:) fliplr(targetLowerBound(plottedRxns(i),:))],colors{mod(i-1,length(colors))+1},'LineWidth',2)
0103     axis tight;
0104     hold on;
0105 end
0106 hold off;
0107 legend(excRxns(plottedRxns));
0108 legend off;
0109 ylabel('Production Rate (mmol/gDW h)');
0110 xlabel('Growth Rate (1/h)');
0111 plottools, plotbrowser('on'), figurepalette('hide'), propertyeditor('off'); 
0112 
0113 biomassValues = biomassValues';
0114 targetValues = [targetLowerBound' targetUpperBound'];
0115 
0116

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