dynamicFBA

PURPOSE ^

dynamicFBA Perform dynamic FBA simulation using the static optimization

SYNOPSIS ^

function [concentrationMatrix,excRxnNames,timeVec,biomassVec] =dynamicFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns)

DESCRIPTION ^

dynamicFBA Perform dynamic FBA simulation using the static optimization
approach

 [concentrationMatrix,excRxnNames,timeVec,biomassVec]
 dynamicFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns)

INPUTS
 model                 COBRA model structure
 substrateRxns         List of exchange reaction names for substrates
                       initially in the media that may change (e.g. not
                       h2o or co2)
 initConcentrations    Initial concentrations of substrates (in the same
                       structure as substrateRxns)
 initBiomass           Initial biomass (must be non zero)
 timeStep              Time step size
 nSteps                Maximum number of time steps

OPTIONAL INPUTS
 plotRxns              Reactions to be plotted (Default =
                       {'EX_glc(e)','EX_ac(e)','EX_for(e)'})
 exclUptakeRxns        List of uptake reactions whose substrate
                       concentrations do not change (Default =
                       {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)'})
 
OUTPUTS
 concentrationMatrix   Matrix of extracellular metabolite concentrations
 excRxnNames           Names of exchange reactions for the EC metabolites
 timeVec               Vector of time points
 biomassVec            Vector of biomass values

 If no initial concentration is given for a substrate that has an open
 uptake in the model (i.e. model.lb < 0) the concentration is assumed to
 be high enough to not be limiting. If the uptake rate for a nutrient is
 calculated to exceed the maximum uptake rate for that nutrient specified
 in the model and the max uptake rate specified is > 0, the maximum uptake 
 rate specified in the model is used instead of the calculated uptake
 rate.

 NOTE: The dynamic FBA method implemented in this function is essentially 
 the same as the method described in
 [Varma, A., and B. O. Palsson. Appl. Environ. Microbiol. 60:3724 (1994)].
 This function does not implement the dynamic FBA using dynamic optimization approach
 described in [Mahadevan, R. et al. Biophys J, 83:1331-1340 (2003)].

 Markus Herrgard 8/22/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [concentrationMatrix,excRxnNames,timeVec,biomassVec] = ...
0002     dynamicFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns)
0003 %dynamicFBA Perform dynamic FBA simulation using the static optimization
0004 %approach
0005 %
0006 % [concentrationMatrix,excRxnNames,timeVec,biomassVec]
0007 % dynamicFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns)
0008 %
0009 %INPUTS
0010 % model                 COBRA model structure
0011 % substrateRxns         List of exchange reaction names for substrates
0012 %                       initially in the media that may change (e.g. not
0013 %                       h2o or co2)
0014 % initConcentrations    Initial concentrations of substrates (in the same
0015 %                       structure as substrateRxns)
0016 % initBiomass           Initial biomass (must be non zero)
0017 % timeStep              Time step size
0018 % nSteps                Maximum number of time steps
0019 %
0020 %OPTIONAL INPUTS
0021 % plotRxns              Reactions to be plotted (Default =
0022 %                       {'EX_glc(e)','EX_ac(e)','EX_for(e)'})
0023 % exclUptakeRxns        List of uptake reactions whose substrate
0024 %                       concentrations do not change (Default =
0025 %                       {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)'})
0026 %
0027 %OUTPUTS
0028 % concentrationMatrix   Matrix of extracellular metabolite concentrations
0029 % excRxnNames           Names of exchange reactions for the EC metabolites
0030 % timeVec               Vector of time points
0031 % biomassVec            Vector of biomass values
0032 %
0033 % If no initial concentration is given for a substrate that has an open
0034 % uptake in the model (i.e. model.lb < 0) the concentration is assumed to
0035 % be high enough to not be limiting. If the uptake rate for a nutrient is
0036 % calculated to exceed the maximum uptake rate for that nutrient specified
0037 % in the model and the max uptake rate specified is > 0, the maximum uptake
0038 % rate specified in the model is used instead of the calculated uptake
0039 % rate.
0040 %
0041 % NOTE: The dynamic FBA method implemented in this function is essentially
0042 % the same as the method described in
0043 % [Varma, A., and B. O. Palsson. Appl. Environ. Microbiol. 60:3724 (1994)].
0044 % This function does not implement the dynamic FBA using dynamic optimization approach
0045 % described in [Mahadevan, R. et al. Biophys J, 83:1331-1340 (2003)].
0046 %
0047 % Markus Herrgard 8/22/06
0048 
0049 if (nargin < 7)
0050     plotRxns = {'EX_glc(e)','EX_ac(e)','EX_for(e)'};
0051 end
0052 
0053 % Uptake reactions whose substrate concentrations do not change
0054 if (nargin < 8)
0055     exclUptakeRxns = {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)'};
0056 end
0057 
0058 % Find exchange rxns
0059 excInd = findExcRxns(model,false);
0060 excInd = excInd & ~ismember(model.rxns,exclUptakeRxns);
0061 excRxnNames = model.rxns(excInd);
0062 length(excRxnNames)
0063 % Figure out if substrate reactions are correct
0064 missingInd = find(~ismember(substrateRxns,excRxnNames));
0065 if (~isempty(missingInd))
0066     for i = 1:length(missingInd)
0067         fprintf('%s\n',substrateRxns{missingInd(i)});
0068     end
0069     error('Invalid substrate uptake reaction!');
0070 end
0071 
0072 % Initialize concentrations
0073 substrateMatchInd = ismember(excRxnNames,substrateRxns);
0074 concentrations = zeros(length(excRxnNames),1);
0075 concentrations(substrateMatchInd) = initConcentrations;
0076 
0077 % Deal with reactions for which there are no initial concentrations
0078 originalBound = -model.lb(excInd);
0079 noInitConcentration = (concentrations == 0 & originalBound > 0);
0080 concentrations(noInitConcentration) = 1000;
0081 
0082 biomass = initBiomass;
0083 
0084 % Initialize bounds
0085 uptakeBound =  concentrations/(biomass*timeStep);
0086 
0087 % Make sure bounds are not higher than what are specified in the model
0088 aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0);
0089 uptakeBound(aboveOriginal) = originalBound(aboveOriginal);
0090 model.lb(excInd) = -uptakeBound;
0091 
0092 concentrationMatrix = sparse(concentrations);
0093 biomassVec = biomass;
0094 timeVec(1) = 0;
0095 
0096 fprintf('Step number\tBiomass\n');
0097 h = waitbar(0,'Dynamic FBA analysis in progress ...');
0098 for stepNo = 1:nSteps
0099     % Run FBA
0100     sol = optimizeCbModel(model,'max','one');
0101     mu = sol.f;
0102     if (sol.stat ~= 1 || mu == 0)
0103         fprintf('No feasible solution - nutrients exhausted\n');
0104         break;
0105     end
0106     uptakeFlux = sol.x(excInd);
0107     biomass = biomass*exp(mu*timeStep);
0108     %biomass = biomass*(1+mu*timeStep);
0109     biomassVec(end+1) = biomass;
0110     
0111     % Update concentrations
0112     concentrations = concentrations - uptakeFlux/mu*biomass*(1-exp(mu*timeStep));
0113     %concentrations = concentrations + uptakeFlux*biomass*timeStep;
0114     concentrations(concentrations <= 0) = 0;
0115     concentrationMatrix(:,end+1) = sparse(concentrations);
0116     
0117     % Update bounds for uptake reactions
0118     uptakeBound =  concentrations/(biomass*timeStep);
0119     % This is to avoid any numerical issues
0120     uptakeBound(uptakeBound > 1000) = 1000;
0121     % Figure out if the computed bounds were above the original bounds
0122     aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0);
0123     % Revert to original bounds if the rate was too high
0124     uptakeBound(aboveOriginal) = originalBound(aboveOriginal);
0125     uptakeBound(abs(uptakeBound) < 1e-9) = 0;
0126 
0127     model.lb(excInd) = -uptakeBound;  
0128     
0129     fprintf('%d\t%f\n',stepNo,biomass);
0130     waitbar(stepNo/nSteps,h);
0131     timeVec(stepNo+1) = stepNo*timeStep;
0132 end
0133 if ( regexp( version, 'R20') )
0134         close(h);
0135 end
0136 
0137 selNonZero = any(concentrationMatrix>0,2);
0138 concentrationMatrix = concentrationMatrix(selNonZero,:);
0139 excRxnNames = excRxnNames(selNonZero);
0140 selPlot = ismember(excRxnNames,plotRxns);
0141 
0142 % Plot concentrations as a function of time
0143 clf
0144 subplot(1,2,1);
0145 plot(timeVec,biomassVec);
0146 axis tight
0147 title('Biomass');
0148 subplot(1,2,2);
0149 plot(timeVec,concentrationMatrix(selPlot,:));
0150 axis tight
0151 legend(strrep(excRxnNames(selPlot),'EX_',''));

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