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
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_',''));