dynamicRFBA

PURPOSE ^

dynamicRFBA - perform dynamic rFBA simulation using the static optimization

SYNOPSIS ^

function [concentrationMatrix,excRxnNames,timeVec,biomassVec,drGenes,constrainedRxns,states] =dynamicRFBA(model,substrateRxns,initConcentrations,initBiomass,timeStep,nSteps,plotRxns,exclUptakeRxns)

DESCRIPTION ^

 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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