runDynamicFBA
Performs dynamic FBA simulation using the static optimization approach
Input:
model a model structure
substrateRxns cell array with exchange reaction identifiers for
substrates that are initially in the media, whose
concentration may change (e.g. not h2o or co2)
initConcentrations numeric initial concentrations of substrates
(matching substrateRxns)
initBiomass numeric initial biomass (must be non-zero)
timeStep numeric time step size
nSteps numeric maximum number of time steps
plotRxns cell array with exchange reaction identifiers for
substrates whose concentration should be plotted
exclUptakeRxns cell array with exchange reaction identifiers for
substrates whose concentration does not change
(e.g. co2, o2, h2o, h)
Output:
concentrationMatrix numeric matrix with extracellular metabolite
concentrations
excRxnNames cell array with exchange reaction identifiers that
match the metabolites included in the
concentrationMatrix
timeVec numeric vector of time points
biomassVec numeric vector with biomass concentrations
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.
Modified from COBRA Toolbox dynamicFBA.m
Usage: [concentrationMatrix, excRxnNames, timeVec, biomassVec] = runDynamicFBA(model, substrateRxns, initConcentrations, initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns)0001 function [concentrationMatrix, excRxnNames, timeVec, biomassVec] = runDynamicFBA(model, substrateRxns, initConcentrations, initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns) 0002 % runDynamicFBA 0003 % Performs dynamic FBA simulation using the static optimization approach 0004 % 0005 % Input: 0006 % model a model structure 0007 % substrateRxns cell array with exchange reaction identifiers for 0008 % substrates that are initially in the media, whose 0009 % concentration may change (e.g. not h2o or co2) 0010 % initConcentrations numeric initial concentrations of substrates 0011 % (matching substrateRxns) 0012 % initBiomass numeric initial biomass (must be non-zero) 0013 % timeStep numeric time step size 0014 % nSteps numeric maximum number of time steps 0015 % plotRxns cell array with exchange reaction identifiers for 0016 % substrates whose concentration should be plotted 0017 % exclUptakeRxns cell array with exchange reaction identifiers for 0018 % substrates whose concentration does not change 0019 % (e.g. co2, o2, h2o, h) 0020 % 0021 % Output: 0022 % concentrationMatrix numeric matrix with extracellular metabolite 0023 % concentrations 0024 % excRxnNames cell array with exchange reaction identifiers that 0025 % match the metabolites included in the 0026 % concentrationMatrix 0027 % timeVec numeric vector of time points 0028 % biomassVec numeric vector with biomass concentrations 0029 % 0030 % If no initial concentration is given for a substrate that has an open 0031 % uptake in the model (i.e. `model.lb < 0`) the concentration is assumed to 0032 % be high enough to not be limiting. If the uptake rate for a nutrient is 0033 % calculated to exceed the maximum uptake rate for that nutrient specified 0034 % in the model and the max uptake rate specified is > 0, the maximum uptake 0035 % rate specified in the model is used instead of the calculated uptake 0036 % rate. 0037 % 0038 % Modified from COBRA Toolbox dynamicFBA.m 0039 % 0040 % Usage: [concentrationMatrix, excRxnNames, timeVec, biomassVec] = runDynamicFBA(model, substrateRxns, initConcentrations, initBiomass, timeStep, nSteps, plotRxns, exclUptakeRxns) 0041 0042 % Find exchange rxns 0043 excRxnNames = getExchangeRxns(model); 0044 excRxnNames(ismember(excRxnNames,exclUptakeRxns))=[]; 0045 excInd = getIndexes(model,excRxnNames,'rxns'); 0046 % Figure out if substrate reactions are correct 0047 missingInd = find(~ismember(substrateRxns,excRxnNames)); 0048 if (~isempty(missingInd)) 0049 for i = 1:length(missingInd) 0050 fprintf('%s\n',substrateRxns{missingInd(i)}); 0051 end 0052 error('Invalid substrate uptake reaction!'); 0053 end 0054 0055 % Initialize concentrations 0056 [~, substrateMatchInd] = ismember(substrateRxns,excRxnNames); 0057 concentrations = zeros(length(excRxnNames),1); 0058 concentrations(substrateMatchInd) = initConcentrations; 0059 0060 % Deal with reactions for which there are no initial concentrations 0061 originalBound = -model.lb(excInd); 0062 noInitConcentration = (concentrations == 0 & originalBound > 0); 0063 concentrations(noInitConcentration) = 1000; 0064 0065 biomass = initBiomass; 0066 0067 % Initialize bounds 0068 uptakeBound = concentrations/(biomass*timeStep); 0069 0070 % Make sure bounds are not higher than what are specified in the model 0071 aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0); 0072 uptakeBound(aboveOriginal) = originalBound(aboveOriginal); 0073 model.lb(excInd) = -uptakeBound; 0074 0075 concentrationMatrix = concentrations; 0076 biomassVec = biomass; 0077 timeVec(1) = 0; 0078 [~,hsSol]=solveLP(model,1); 0079 PB = ProgressBar2(nSteps,'Running dynamic FBA analysis','cli'); 0080 for stepNo = 1:nSteps 0081 % Run FBA 0082 [sol,hsSol] = solveLP(model,1,[],hsSol); 0083 mu = -sol.f; 0084 if (sol.stat ~= 1 || mu == 0) 0085 fprintf('\nNo feasible solution - nutrients exhausted. Biomass:\t %f\n', biomass); 0086 break; 0087 end 0088 uptakeFlux = sol.x(excInd); 0089 biomass = biomass*exp(mu*timeStep); 0090 biomassVec(end+1) = biomass; 0091 0092 % Update concentrations 0093 concentrations = concentrations - uptakeFlux/mu*biomass*(1-exp(mu*timeStep)); 0094 concentrations(concentrations <= 0) = 0; 0095 concentrationMatrix(:,end+1) = concentrations; 0096 0097 % Update bounds for uptake reactions 0098 uptakeBound = concentrations/(biomass*timeStep); 0099 % This is to avoid any numerical issues 0100 uptakeBound(uptakeBound > 1000) = 1000; 0101 % Figure out if the computed bounds were above the original bounds 0102 aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0); 0103 % Revert to original bounds if the rate was too high 0104 uptakeBound(aboveOriginal) = originalBound(aboveOriginal); 0105 uptakeBound(abs(uptakeBound) < 1e-9) = 0; 0106 0107 model.lb(excInd) = -uptakeBound; 0108 timeVec(stepNo+1) = stepNo*timeStep; 0109 0110 count(PB); 0111 end 0112 0113 selNonZero = any(concentrationMatrix>0,2); 0114 concentrationMatrix = concentrationMatrix(selNonZero,:); 0115 excRxnNames = excRxnNames(selNonZero); 0116 selPlot = ismember(excRxnNames,plotRxns); 0117 0118 % Plot concentrations as a function of time 0119 clf 0120 subplot(1,2,1); 0121 plot(timeVec,biomassVec); 0122 axis tight 0123 title('Biomass'); 0124 xlabel('Time (h)'); 0125 ylabel('Concentration (g/L)'); 0126 subplot(1,2,2); 0127 plot(timeVec,concentrationMatrix(selPlot,:)); 0128 axis tight 0129 title('Substrates and/or products'); 0130 xlabel('Time (h)'); 0131 ylabel('Concentration (mmol/L)'); 0132 legend(strrep(excRxnNames(selPlot),'EX_','')); 0133 end