dynamicRFBA - perform dynamic rFBA simulation using the static optimization approach [concentrationMatrix,excRxnNames,timeVec,biomassVec,drGenes,constrainedRxns,states] = ... dynamicRFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns) model a regulatory COBRA model substrateRxns list of exchange reaction names for substrates initially in the media that may change (i.e. not h2o or co2) initConcentrations initial concentrations of substrates (in the same structure as substrateRxns) initBiomass initial biomass timeStep time step size nSteps maximum number of time steps plotRxns reactions to be plotted exclUptakeRxns list of uptake reactions whose substrate concentrations do not change (opt, default {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)'}) 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 drGenes vector of downregulated genes constrainedRxns vector of downregulated reactions states vector of regulatory network states 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)]. Jeff Orth 9/15/08 (modified dynamicFBA by Markus Herrgard 8/22/06)
0001 function [concentrationMatrix,excRxnNames,timeVec,biomassVec,drGenes,constrainedRxns,states] = ... 0002 dynamicRFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns) 0003 % dynamicRFBA - perform dynamic rFBA simulation using the static optimization 0004 % approach 0005 % 0006 % [concentrationMatrix,excRxnNames,timeVec,biomassVec,drGenes,constrainedRxns,states] = ... 0007 % dynamicRFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns) 0008 % 0009 % model a regulatory COBRA model 0010 % substrateRxns list of exchange reaction names for substrates 0011 % initially in the media that may change (i.e. not 0012 % h2o or co2) 0013 % initConcentrations initial concentrations of substrates (in the same 0014 % structure as substrateRxns) 0015 % initBiomass initial biomass 0016 % timeStep time step size 0017 % nSteps maximum number of time steps 0018 % plotRxns reactions to be plotted 0019 % exclUptakeRxns list of uptake reactions whose substrate 0020 % concentrations do not change (opt, default 0021 % {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)'}) 0022 % 0023 % concentrationMatrix matrix of extracellular metabolite concentrations 0024 % excRxnNames names of exchange reactions for the EC metabolites 0025 % timeVec vector of time points 0026 % biomassVec vector of biomass values 0027 % drGenes vector of downregulated genes 0028 % constrainedRxns vector of downregulated reactions 0029 % states vector of regulatory network states 0030 % 0031 % If no initial concentration is given for a substrate that has an open 0032 % uptake in the model (i.e. model.lb < 0) the concentration is assumed to 0033 % be high enough to not be limiting. If the uptake rate for a nutrient is 0034 % calculated to exceed the maximum uptake rate for that nutrient specified 0035 % in the model and the max uptake rate specified is > 0, the maximum uptake 0036 % rate specified in the model is used instead of the calculated uptake rate. 0037 % 0038 % NOTE: The dynamic FBA method implemented in this function is essentially 0039 % the same as the method described in 0040 % [Varma, A., and B. O. Palsson. Appl. Environ. Microbiol. 60:3724 (1994)]. 0041 % This function does not implement the dynamic FBA using dynamic optimization approach 0042 % described in [Mahadevan, R. et al. Biophys J, 83:1331-1340 (2003)]. 0043 % 0044 % Jeff Orth 9/15/08 (modified dynamicFBA by Markus Herrgard 8/22/06) 0045 0046 if (nargin < 7) 0047 plotRxns = {'EX_glc(e)','EX_ac(e)','EX_for(e)'}; 0048 end 0049 0050 % Uptake reactions whose substrate concentrations do not change 0051 if (nargin < 8) 0052 exclUptakeRxns = {'EX_co2(e)','EX_o2(e)','EX_h2o(e)','EX_h(e)','EX_nh4(e)','EX_pi(e)'}; 0053 end 0054 0055 % Find exchange rxns 0056 excInd = findExcRxns(model,false); 0057 excInd = excInd & ~ismember(model.rxns,exclUptakeRxns); 0058 excRxnNames = model.rxns(excInd); 0059 0060 % Figure out if substrate reactions are correct 0061 missingInd = find(~ismember(substrateRxns,excRxnNames)); 0062 if (~isempty(missingInd)) 0063 for i = 1:length(missingInd) 0064 fprintf('%s\n',substrateRxns{missingInd(i)}); 0065 end 0066 error('Invalid substrate uptake reaction!'); 0067 end 0068 0069 % Initialize concentrations 0070 substrateMatchInd = ismember(excRxnNames,substrateRxns); 0071 concentrations = zeros(length(excRxnNames),1); 0072 concentrations(substrateMatchInd) = initConcentrations; 0073 0074 % Deal with reactions for which there are no initial concentrations 0075 originalBound = -model.lb(excInd); 0076 noInitConcentration = (concentrations == 0 & originalBound > 0); 0077 concentrations(noInitConcentration) = 1000; 0078 0079 biomass = initBiomass; 0080 0081 % Initialize bounds 0082 uptakeBound = concentrations/(biomass*timeStep); 0083 0084 % Make sure bounds are not higher than what are specified in the model 0085 aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0); 0086 uptakeBound(aboveOriginal) = originalBound(aboveOriginal); 0087 model.lb(excInd) = -uptakeBound; 0088 0089 % Initialize the regulatory network with optimizeRegModel 0090 [FBAsols,DRgenes,constrainedRxns,cycleStart,modStates] = optimizeRegModel(model); 0091 %find a growth solution 0092 for i = 1:length(FBAsols) 0093 sol = FBAsols{i}; 0094 if sol.f ~= 0 0095 break 0096 end 0097 end 0098 genes = DRgenes{i}; 0099 rxns = constrainedRxns{i}; 0100 iniState = modStates{i+cycleStart}(1:length(model.regulatoryGenes)); 0101 inputs1state = modStates{i+cycleStart}((length(model.regulatoryGenes)+1):(length(model.regulatoryGenes)+length(model.regulatoryInputs1))); 0102 inputs2state = modStates{i+cycleStart}((length(model.regulatoryGenes)+length(model.regulatoryInputs1)+1):(length(model.regulatoryGenes)+length(model.regulatoryInputs1)+length(model.regulatoryInputs2))); 0103 modelDR = deleteModelGenes(model,genes); % set rxns to 0 0104 0105 concentrationMatrix = sparse(concentrations); 0106 biomassVec = biomass; 0107 timeVec(1) = 0; 0108 % regulatory states 0109 drGenes{1} = genes; 0110 constrainedRxns{1} = rxns; 0111 states{1} = iniState; 0112 0113 noGrowthCount = 0; 0114 %fprintf('Step number\tBiomass\n'); 0115 h = waitbar(0,'Dynamic regulatory FBA analysis in progress ...'); 0116 for stepNo = 1:nSteps 0117 % Run FBA 0118 sol = optimizeCbModel(modelDR,'max',true); 0119 mu = sol.f; 0120 0121 if (sol.stat ~= 1 | mu == 0) % end if no growth for 10 steps 0122 noGrowthCount = noGrowthCount+1; 0123 biomass = biomassVec(end); % no growth 0124 if noGrowthCount >= 20 0125 fprintf('No feasible solution - nutrients exhausted\n'); 0126 break; 0127 end 0128 else 0129 uptakeFlux = sol.x(excInd); 0130 biomass = biomass*exp(mu*timeStep); 0131 % Update concentrations only if growth occurs 0132 concentrations = concentrations - uptakeFlux/mu*biomass*(1-exp(mu*timeStep)); 0133 %concentrations = concentrations + uptakeFlux*biomass*timeStep; 0134 concentrations(concentrations <= 0) = 0; 0135 end 0136 0137 biomassVec(end+1) = biomass; 0138 concentrationMatrix(:,end+1) = sparse(concentrations); 0139 0140 % Update bounds for uptake reactions 0141 uptakeBound = concentrations/(biomass*timeStep); 0142 % This is to avoid any numerical issues 0143 uptakeBound(uptakeBound > 1000) = 1000; 0144 % Figure out if the computed bounds were above the original bounds 0145 aboveOriginal = (uptakeBound > originalBound) & (originalBound > 0); 0146 % Revert to original bounds if the rate was too high 0147 uptakeBound(aboveOriginal) = originalBound(aboveOriginal); 0148 uptakeBound(abs(uptakeBound) < 1e-9) = 0; 0149 0150 model.lb(excInd) = -uptakeBound; 0151 0152 % get current regulatory state and downregulate reactions 0153 [finalState,finalInputs1States,finalInputs2States] = solveBooleanRegModel(model,states{end},inputs1state,inputs2state); 0154 KOgenes = {}; 0155 for i = 1:length(model.regulatoryGenes) 0156 if finalState(i) == false 0157 KOgenes{end+1,1} = model.regulatoryGenes{i}; 0158 end 0159 end 0160 genes = intersect(model.genes,KOgenes); % remove genes not associated with rxns 0161 [modelDR,he,rxns] = deleteModelGenes(model,genes); % set rxns to 0 0162 inputs1state = finalInputs1States; 0163 inputs2state = finalInputs2States; %reset inputs 1 and 2 states to current conditions 0164 0165 drGenes{end+1} = genes; 0166 constrainedRxns{end+1} = rxns; 0167 states{end+1} = finalState; 0168 0169 fprintf('%d\t%f\n',stepNo,biomass); 0170 waitbar(stepNo/nSteps,h); 0171 timeVec(stepNo+1) = stepNo*timeStep; 0172 end 0173 close(h); 0174 0175 selNonZero = any(concentrationMatrix>0,2); 0176 concentrationMatrix = concentrationMatrix(selNonZero,:); 0177 excRxnNames = excRxnNames(selNonZero); 0178 selPlot = ismember(excRxnNames,plotRxns); 0179 0180 % Plot concentrations as a function of time 0181 clf 0182 subplot(1,2,1); 0183 plot(timeVec,biomassVec,'LineWidth',2); 0184 axis tight 0185 title('Biomass'); 0186 subplot(1,2,2); 0187 plot(timeVec,concentrationMatrix(selPlot,:),'LineWidth',2); 0188 axis tight 0189 legend(strrep(excRxnNames(selPlot),'EX_','')); 0190 0191 0192 0193