Commit 5d79ee9b authored by Florian Centler's avatar Florian Centler
Browse files

This is version v1.1.0

parent d04e0e7b
......@@ -61,13 +61,13 @@ The FBA models are taken from *J. J. Hamilton, M. Calixto Contreras, and J. L. R
### Example 3
Simulating batch growth of a 773 species human gut microbiome first requires the unpacking of the file `AGORA-1.01-Western-Diet.zip` containing the AGORA model collection. Note that for running the simulation with all species, 64GB of RAM are necessary (loading models as a Matlab data structure after the initial loading as SBML files cuts memory demand in half). The simulation of batch growth can then be started by:
Simulating batch growth of an eight species human gut microbiome (SIHUMIx consortium). Models are taken from the AGORA model collection containing 773 human gut microbiome species. Running this simulation first requires the unpacking of the file `AGORA-1.01-Western-Diet.zip` containing the AGORA models. The simulation can then be started by:
```
microbialSimMain(3)
```
Simulation time can considerably be reduced by setting the parameter `maxDeviation` to `inf` in `microbialSimMain.m` at the expense of numerical accuracy. Note that also arbitrary subsets of the model collection can be selected for the simulation.
Any other subset of AGORA species, and even all 773 species together can be simulated. Note that for running the simulation with all 773 species, 64GB of RAM are necessary (loading models as a Matlab data structure after the initial loading as SBML files cuts memory demand in half). For complex microbiomes, consider deactivating the option to record limiting fluxes `recordLimitingFluxes` for faster simulation times. Speed can considerably be further improved by setting the parameter `maxDeviation` to `inf` in `microbialSimMain.m` at the expense of numerical accuracy.
The FBA models stem from *S. Magnúsdóttir et al., “Generation of genome-scale metabolic reconstructions for 773 members of the human gut microbiota,” Nat. Publ. Gr., vol. 35, no. 1, 2016.*
......@@ -90,6 +90,11 @@ Please read [CONTRIBUTING.md](https://gist.github.com/PurpleBooth/b24679402957c6
We use [SemVer](http://semver.org/) for versioning. For the versions available, see the [tags on this repository](https://git.ufz.de/UMBSysBio/microbialSim/tags).
## Version History
* v1.0.0 - July 2019 - initial version
* v1.1.0 - December 2019 - improved parallel performance, new function to estimate maximal uptake flux required to achieve given specific growth rate, new functions to collect, filter, and write compound exchange within the community at a given simulation time point, new feature to record growth-limiting compounds for each species during the simulation, bugfixes
## Authors
* **Florian Centler**, UFZ - Helmholtz Centre for Environmental Research, Leipzig, Germany
......
......@@ -69,6 +69,7 @@ for i = 1:numberOfReactions
error('Case not considered yet!')
end
model.COBRAmodel.rxns(end + 1) = strcat(model.COBRAmodel.rxns(i),'-reverse');
model.COBRAmodel.rxnNames(end + 1) = strcat(model.COBRAmodel.rxnNames(i),'-reverse');
otherwise
error('In convertToOnlyIrrecReac: unknown FBA solver')
......
function [ result ] = dFBASimulator( metabolicModel, reactor, type, solverPars)
function [ result ] = dFBASimulator( metabolicModel, reactor, type, solverPars, workerModels)
%DFBASIMULATOR Running a dynamic FBA simulation with bialSim
% Running the dFBA Simulation. For type == 0, the step-wise iteration is
% used, type == 1 uses the Matlab ODE solver.
......@@ -16,10 +16,20 @@ else
myInit = load(initFile);
startTime = myInit.recordResult.time;
stronglyConsumedFlag = myInit.recordResult.stronglyConsumedFlag;
for m = 1:length(metabolicModel)
metabolicModel(m).lastFlux = myInit.recordResult.metabolicModel(m).lastFlux;
if solverPars.parallel == 1
for i = 1:length(myInit.recordResult.metabolicModel)
myInit_recordResult_metabolicModel_lastFlux{i} = myInit.recordResult.metabolicModel(i).lastFlux;
end
spmd
for i = 1:length(workerModels)
workerModels(i).lastFlux = myInit_recordResult_metabolicModel_lastFlux{labindex+(i-1)*numlabs};
end
end
else
for m = 1:length(metabolicModel)
metabolicModel(m).lastFlux = myInit.recordResult.metabolicModel(m).lastFlux;
end
end
end
timeSteps = floor((solverPars.tend - startTime) / solverPars.timeStepSize);
......@@ -41,13 +51,20 @@ for i = 1:length(metabolicModel)
% Flux reporting will be done based on number of reactions in
% original model.
end
numberOfCoupledReactions(i) = length(metabolicModel(i).coupledReactions.ReacID);
result.FBA(i).fluxes = zeros(timeSteps, numberOfReactions(i));
if isempty(initFile)
metabolicModel(i).lastFlux=[];
end
end
if solverPars.parallel == true
spmd
for i = 1:length(workerModels)
workerModels(i).lastFlux=[];
end
end
end
if isempty(initFile)
result.time(1) = 0;
result.compounds(1,:) = reactor.compoundsInit;
......@@ -58,6 +75,10 @@ else
result.biomass(1,:) = myInit.recordResult.biomass;
end
if (solverPars.recordLimitingFluxes == 1)
result.limitingFluxes = cell2table(cell(0,9), 'VariableNames', {'ModelID', 'ModelName', 'Time', 'FluxID', 'FluxName', 'FluxValue', 'IsFixed', 'Compound', 'ExchangeCoupledFlux'});
end
% DO SIMULATION
if type == 1
......@@ -74,9 +95,9 @@ if type == 1
end
switch solverPars.solver
case 1
[fba.t,fba.y] = ode45(@(t,y) getDyWrapper(t,y,metabolicModel,reactor, solverPars), [0:solverPars.timeStepSize:solverPars.tend], [reactor.biomassInit reactor.compoundsInit], options);
[fba.t,fba.y] = ode45(@(t,y) getDyWrapper(y,metabolicModel,reactor, solverPars, workerModels), [0:solverPars.timeStepSize:solverPars.tend], [reactor.biomassInit reactor.compoundsInit], options);
case 2
[fba.t,fba.y] = ode15s(@(t,y) getDyWrapper(t,y,metabolicModel,reactor, solverPars), [0:solverPars.timeStepSize:solverPars.tend], [reactor.biomassInit reactor.compoundsInit], options);
[fba.t,fba.y] = ode15s(@(t,y) getDyWrapper(y,metabolicModel,reactor, solverPars, workerModels), [0:solverPars.timeStepSize:solverPars.tend], [reactor.biomassInit reactor.compoundsInit], options);
otherwise
error('Unknown solver.')
end
......@@ -98,23 +119,45 @@ else
if solveFullFBA == 1
clearvars dbio dcompound flux success
if solverPars.parallel == 1
ticBytes(gcp);
parfor i = 1:length(metabolicModel)
currentModel = metabolicModel(i);
[dbio(i), dcompoundRAW{i}, flux{i}, success(i)] = getDy(currentModel, reactor, result.compounds(n, :), result.biomass(n, i), reactor.biomassInflow(i), solverPars);
end
tocBytes(gcp)
%ticBytes(gcp);
myResult.compounds = result.compounds(n, :);
myResult.biomass = result.biomass(n, :);
spmd
for i = 1:length(workerModels)
[mydbio(i), mydcompoundRAW{i}, myflux{i}, mysuccess(i), mylimitedFluxes{i}] = getDy(workerModels(i), myResult.compounds, myResult.biomass(labindex + numlabs * (i-1)), solverPars);
end
end
%tocBytes(gcp)
for i = 1:length(mydbio)
dbio(i:length(mydbio):length(metabolicModel)) = mydbio{i};
dcompoundRAW(i:length(mydbio):length(metabolicModel)) = mydcompoundRAW{i};
flux(i:length(mydbio):length(metabolicModel)) = myflux{i};
success(i:length(mydbio):length(metabolicModel)) = mysuccess{i};
limitedFluxes(i:length(mydbio):length(metabolicModel)) = mylimitedFluxes{i};
end
else
for i = 1:length(metabolicModel)
[dbio(i), dcompoundRAW{i}, flux{i}, success(i)] = getDy(metabolicModel(i), reactor, result.compounds(n, :), result.biomass(n, i), reactor.biomassInflow(i), solverPars);
[dbio(i), dcompoundRAW{i}, flux{i}, success(i), limitedFluxes{i}] = getDy(metabolicModel(i), result.compounds(n, :), result.biomass(n, i), solverPars);
end
end
for i = 1:length(metabolicModel)
metabolicModel(i).previousFlux = metabolicModel(i).lastFlux;
metabolicModel(i).lastFlux = flux{i};
dcompound(i,:) = mapExchangeToReactorCompounds(dcompoundRAW{i}, metabolicModel(i).reactorCompoundIDs, length(reactor.compounds));
end
if solverPars.parallel == 1
spmd
for i = 1:length(workerModels)
workerModels(i).previousFlux = workerModels(i).lastFlux;
workerModels(i).lastFlux = flux{labindex+(i-1)*numlabs};
end
end
else
for i = 1:length(metabolicModel)
metabolicModel(i).previousFlux = metabolicModel(i).lastFlux;
metabolicModel(i).lastFlux = flux{i};
end
end
if length(metabolicModel) > 1
myDcompound = sum(dcompound);
else
......@@ -124,6 +167,43 @@ else
dbioTotal = dbio + considerFlow(reactor.biomassInflow, result.biomass(n, :), reactor);
dcompoundFlow = considerFlow(reactor.compoundsInflow, result.compounds(n, :), reactor);
dcompoundTotal = myDcompound + dcompoundFlow;
% record limiting fluxes
if(solverPars.recordLimitingFluxes == 1)
for i = 1:length(metabolicModel)
for j = 1:size(limitedFluxes{i}, 1)
%'ModelID', 'ModelName', 'Time', 'FluxID', 'FluxName',
%'FluxValue', 'IsFixed', 'Compound', 'Exchange/Coupled Flux'
newEntry = {};
newEntry{1,1} = i;
newEntry{1,2} = metabolicModel(i).modelName;
newEntry{1,3} = result.time(n);
newEntry{1,4} = limitedFluxes{i}(j,1);
switch metabolicModel(i).FBAsolver
case 1
newEntry{1,5} = metabolicModel(i).CNAmodel.reacID(limitedFluxes{i}(j, 1));
case 2
newEntry{1,5} = metabolicModel(i).COBRAmodel.rxnNames(limitedFluxes{i}(j, 1));
end
newEntry{1,6} = limitedFluxes{i}(j,2);
newEntry{1,7} = limitedFluxes{i}(j,3);
if ismember(newEntry{1,4}, metabolicModel(i).exchangeReactions.ReacID)
newEntry{1,9} = 1;
else
newEntry{1,9} = 0;
end
myIndex = find(newEntry{1,4} == metabolicModel(i).coupledReactions.ReacID);
if length(myIndex) == 1
newEntry{1,8} = reactor.compounds(metabolicModel(i).reactorCompoundIDs(myIndex));
newEntry{1,9} = 2;
else
newEntry{1,8} = '-';
end
result.limitingFluxes = [result.limitingFluxes;newEntry];
end
end
end
end
%Update reactor with changes due to biotic activity and flow
......@@ -256,42 +336,86 @@ end
% do mass balance for all exchange fluxes
if solverPars.doMassBalance == 1
for i = 1:length(metabolicModel) % for all models
result.FBA(i).exchangeReactions.MassBalance = 0 * ones(height(result.FBA(i).exchangeReactions), 1); % init with zeros
for j = 1:height(result.FBA(i).exchangeReactions) % for all exchange reactions
for k = 1:size(result.FBA(i).fluxes, 1) % for all time steps / recorded rates
result.FBA(i).exchangeReactions{j, 'MassBalance'} = result.FBA(i).exchangeReactions{j, 'MassBalance'} + result.FBA(i).fluxes(k, result.FBA(i).exchangeReactions{j, 'ReacID'}) * result.biomass(k, i) * (result.time(k+1)-result.time(k));
if solverPars.parallel == true
modelBiomass = result.biomass;
modelTime = result.time;
% ticBytes(gcp);
modelFluxes = Composite();
modelBiomass = Composite();
for i = 1:numel(modelFluxes)
modelFluxes{i} = cat(2, result.FBA(i:numel(modelFluxes):length(metabolicModel)).fluxes);
modelBiomass{i} = result.biomass(:,i:numel(modelFluxes):length(metabolicModel));
end
spmd
fluxStartIndex = 1;
for i = 1:length(workerModels)
fluxEndIndex = fluxStartIndex - 1 + length(workerModels(i).lastFlux) - length(workerModels(i).reversibleReacIDs);
for j = 1:height(workerModels(i).exchangeReactions)
myMassBalance(i, j) = 0;
for k = 1:size(modelFluxes, 1)
myMassBalance(i, j) = myMassBalance(i, j) + modelFluxes(k, fluxStartIndex - 1 + workerModels(i).exchangeReactions{j, 'ReacID'}) * modelBiomass(k, i) * (modelTime(k+1)-modelTime(k));
end
end
abs_sum(i) = sum(sum(abs(modelFluxes(:,fluxStartIndex:fluxEndIndex))));
dev_sum(i) = 0;
for j = 2:size(modelFluxes, 1)
dev_sum(i) = dev_sum(i) + sum(abs(modelFluxes(j,fluxStartIndex:fluxEndIndex) - modelFluxes(j-1,fluxStartIndex:fluxEndIndex)));
end
fluxStartIndex = fluxEndIndex + 1;
end
end
% tocBytes(gcp);
for i = 1:length(myMassBalance)
tmp = abs_sum{i};
tmp2 = dev_sum{i};
k = 1;
for j = i:length(myMassBalance):length(metabolicModel)
result.FBA(j).flux_sum = tmp(k);
result.FBA(j).flux_dev_sum = tmp2(k);
k = k + 1;
end
myMB = myMassBalance{i};
myIdx = 1;
for j = i:length(myMassBalance):length(metabolicModel)
result.FBA(j).exchangeReactions = addvars(result.FBA(j).exchangeReactions, transpose(myMB(myIdx, 1:height(result.FBA(j).exchangeReactions))), 'NewVariableNames', 'MassBalance');
myIdx = myIdx + 1;
end
end
else
for i = 1:length(metabolicModel) % for all models
result.FBA(i).exchangeReactions.MassBalance = zeros(height(result.FBA(i).exchangeReactions), 1); % init with zeros
for j = 1:height(result.FBA(i).exchangeReactions) % for all exchange reactions
for k = 1:size(result.FBA(i).fluxes, 1) % for all time steps / recorded rates
result.FBA(i).exchangeReactions{j, 'MassBalance'} = result.FBA(i).exchangeReactions{j, 'MassBalance'} + result.FBA(i).fluxes(k, result.FBA(i).exchangeReactions{j, 'ReacID'}) * result.biomass(k, i) * (result.time(k+1)-result.time(k));
end
end
% compute absolute fluxes & flux difference
result.FBA(i).flux_sum = sum(sum(abs(result.FBA(i).fluxes)));
result.FBA(i).flux_dev_sum = 0;
for j = 2:size(result.FBA(i).fluxes,1)
result.FBA(i).flux_dev_sum = result.FBA(i).flux_dev_sum + sum(abs(result.FBA(i).fluxes(j,:) - result.FBA(i).fluxes(j-1,:)));
end
end
end
for i = 1:length(metabolicModel) % for all models
% copy data to coupledReactions table
result.FBA(i).coupledReactions.MassBalance = 0 * ones(height(result.FBA(i).coupledReactions), 1); % init with zeros
result.FBA(i).coupledReactions.MassBalance = zeros(height(result.FBA(i).coupledReactions), 1); % init with zeros
for j = 1:height(result.FBA(i).coupledReactions)
result.FBA(i).coupledReactions{j, 'MassBalance'} = result.FBA(i).exchangeReactions{int2str(result.FBA(i).coupledReactions{j, 'ReacID'}), 'MassBalance'};
end
% compute absolute fluxes & flux difference
abs_sum = 0;
dev_sum = 0;
for j = 2:size(result.FBA(i).fluxes,1)
abs_sum = abs_sum + sum(sum(abs(result.FBA(i).fluxes)));
dev_sum = dev_sum + sum(abs(result.FBA(i).fluxes(j,:) - result.FBA(i).fluxes(j-1,:)));
end
result.FBA(i).flux_sum = abs_sum;
result.FBA(i).flux_dev_sum = dev_sum;
% Add final flux info to reaction tables
result.FBA(i).exchangeReactions.finalFlux = getFluxes(result.FBA(i).fluxes(size(result.FBA(i).fluxes, 1),:), result.FBA(i).exchangeReactions.ReacID)';
result.FBA(i).coupledReactions.finalFlux = getFluxes(result.FBA(i).fluxes(size(result.FBA(i).fluxes, 1),:), result.FBA(i).coupledReactions.ReacID)';
fluxLimits = getFluxLimits(result.FBA(i).fluxes(size(result.FBA(i).fluxes, 1),:), metabolicModel(i), result.compounds(size(result.FBA(i).fluxes, 1),:));
fluxLimits = getLimitedFluxes(result.FBA(i).fluxes(size(result.FBA(i).fluxes, 1),:), metabolicModel(i), solverPars, result.compounds(size(result.FBA(i).fluxes, 1),:));
result.FBA(i).exchangeReactions.finalFluxLimit = fluxLimits(result.FBA(i).exchangeReactions.ReacID)';
result.FBA(i).coupledReactions.finalFluxLimit = fluxLimits(result.FBA(i).coupledReactions.ReacID)';
end
end
result.timePars = solverPars;
result.solverPars = solverPars;
result.compoundNames = reactor.compounds;
result.modelNames = {};
for i = 1:length(metabolicModel)
......@@ -301,5 +425,10 @@ for i = 1:length(metabolicModel)
result.modelNames{i} = metabolicModel(i).modelName;
end
result.modelNames = [result.modelNames{:}];
if solverPars.parallel == 1
p = gcp('nocreate');
result.numberOfWorkers = p.NumWorkers;
end
end
function [vmax, flux, exchangeFluxes] = estimateVmax(metabolicModel, reactor, solverPars, coupledReactionID, targetMu)
%ESTIMATEVMAX Given a growth rate, this function computes which uptake flux
%of the growth limiting substrate gives rise to this growth rate, assuming
%that all other uptake fluxes are set to reasonable / problem-specific
%values. If
% Detailed explanation goes here
orig.vmax = metabolicModel.coupledReactions.vmax;
orig.ks = metabolicModel.coupledReactions.ks;
numberOfReactions = length(metabolicModel.coupledReactions.vmax);
%% report on growth rate under conditions specified in model
metabolicModel.coupledReactions.ks = zeros(numberOfReactions, 1);
[ dbio, dcompound, flux, success ] = getDy(metabolicModel, ones(1,length(reactor.compounds)), -1, solverPars);
if success == true
fprintf("Growth rate under original model uptake conditions = %d 1/h.\n", flux(metabolicModel.biomassReac))
else
fprintf("Original conditions do not allow for growth.\n");
end
%% check if selected compound is required for growth; if not: adjust all
% flux boundaries
growthLimiting = true;
metabolicModel.coupledReactions.vmax(coupledReactionID) = 0;
[ dbio, dcompound, flux, success ] = getDy(metabolicModel, ones(1,length(reactor.compounds)), -1, solverPars);
if success == true
growthLimiting = false;
fprintf("Selected compound not required for growth, adjusting all flux limits together.\n");
fprintf("Check if uncoupled (and coupled) exchange reactions still contain\ngrowth-relevant fluxes (such as c-sources).\n")
else
fprintf("Selected compound is required for growth.\n");
end
if growthLimiting == true
stepSize = 10;
subStepSize = 1;
else
stepSize = 0.1;
subStepSize = 0.01;
switch metabolicModel.FBAsolver
case 1
orig.lb = metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID);
orig.ub = metabolicModel.CNAmodel.reacMax(metabolicModel.exchangeReactions.ReacID);
case 2
orig.lb = metabolicModel.COBRAmodel.lb(metabolicModel.exchangeReactions.ReacID);
orig.ub = metabolicModel.COBRAmodel.ub(metabolicModel.exchangeReactions.ReacID);
otherwise
error('In estimateVmax: unknown FBA solver')
end
end
%% increase vmax until growth rate is larger then target growth rate
currentMu = 0;
if growthLimiting == true
vmax = 0;
else
vmax = stepSize/1000.0 - stepSize; % factor applied to all uptake fluxes
end
while currentMu < targetMu && vmax < 1000.0
vmax = vmax + stepSize;
if growthLimiting == true
metabolicModel.coupledReactions.vmax(coupledReactionID) = vmax;
else
for i = 1:length(metabolicModel.exchangeReactions.ReacID)
if metabolicModel.exchangeReactions.SecretionSense(i) == 1
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID(i)) = orig.lb(i) * vmax;
case 2
metabolicModel.COBRAmodel.lb(metabolicModel.exchangeReactions.ReacID(i)) = orig.lb(i) * vmax;
otherwise
error('In estimateVmax: unknown FBA solver')
end
else
if metabolicModel.exchangeReactions.SecretionSense(i) == -1
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMax(metabolicModel.exchangeReactions.ReacID(i)) = orig.ub(i) * vmax;
case 2
metabolicModel.COBRAmodel.ub(metabolicModel.exchangeReactions.ReacID(i)) = orig.ub(i) * vmax;
otherwise
error('In estimateVmax: unknown FBA solver')
end
end
end
end
metabolicModel.coupledReactions.vmax = orig.vmax * vmax;
end
[ dbio, dcompound, flux, success ] = getDy(metabolicModel, ones(1,length(reactor.compounds)), -1, solverPars);
if success == true
currentMu = flux(metabolicModel.biomassReac);
else
currentMu = 0;
end
end
if currentMu < targetMu
fprintf('Target (%d) not reachable, just reached %d 1/h\n', targetMu, currentMu);
vmax = -1;
return
end
%% establish lower bound with growth
lowBound = vmax - stepSize - subStepSize;
solvable = false;
while solvable == false
lowBound = lowBound + subStepSize;
if growthLimiting == true
metabolicModel.coupledReactions.vmax(coupledReactionID) = lowBound;
else
for i = 1:length(metabolicModel.exchangeReactions.ReacID)
if metabolicModel.exchangeReactions.SecretionSense(i) == 1
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID(i)) = orig.lb(i) * lowBound;
case 2
metabolicModel.COBRAmodel.lb(metabolicModel.exchangeReactions.ReacID(i)) = orig.lb(i) * lowBound;
otherwise
error('In estimateVmax: unknown FBA solver')
end
else
if metabolicModel.exchangeReactions.SecretionSense(i) == -1
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMax(metabolicModel.exchangeReactions.ReacID(i)) = orig.ub(i) * lowBound;
case 2
metabolicModel.COBRAmodel.ub(metabolicModel.exchangeReactions.ReacID(i)) = orig.ub(i) * lowBound;
otherwise
error('In estimateVmax: unknown FBA solver')
end
end
end
end
metabolicModel.coupledReactions.vmax = orig.vmax * lowBound;
end
[ dbio, dcompound, flux, success ] = getDy(metabolicModel, ones(1,length(reactor.compounds)), -1, solverPars);
if success == true
solvable = true;
end
end
if flux(metabolicModel.biomassReac) > targetMu
fprintf('Failed to find proper lower bound. Chose smaller subStepSize\n');
vmax = -1;
return
end
%% now approach target mu
highBound = vmax;
while abs(currentMu - targetMu) > solverPars.minimalGrowth
vmax = (lowBound + highBound) / 2.0;
if growthLimiting == true
metabolicModel.coupledReactions.vmax(coupledReactionID) = vmax;
else
for i = 1:length(metabolicModel.exchangeReactions.ReacID)
if metabolicModel.exchangeReactions.SecretionSense(i) == 1
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID(i)) = orig.lb(i) * vmax;
case 2
metabolicModel.COBRAmodel.lb(metabolicModel.exchangeReactions.ReacID(i)) = orig.lb(i) * vmax;
otherwise
error('In estimateVmax: unknown FBA solver')
end
else
if metabolicModel.exchangeReactions.SecretionSense(i) == -1
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID(i)) = orig.ub(i) * vmax;
case 2
metabolicModel.COBRAmodel.ub(metabolicModel.exchangeReactions.ReacID(i)) = orig.ub(i) * vmax;
otherwise
error('In estimateVmax: unknown FBA solver')
end
end
end
end
metabolicModel.coupledReactions.vmax = orig.vmax * vmax;
end
[ dbio, dcompound, flux, success ] = getDy(metabolicModel, ones(1,length(reactor.compounds)), -1, solverPars);
if success == true
currentMu = flux(metabolicModel.biomassReac);
if currentMu < targetMu
lowBound = vmax;
else
highBound = vmax;
end
else
fprintf('Fail in approaching target !\n');
vmax = -1;
return
end
end
[exchangeFluxes, limitedFluxesIDs] = getExchangeFluxes(metabolicModel, flux, solverPars, ones(1,length(reactor.compounds)));
if growthLimiting == true
fprintf('For getting = %d 1/h, choose vmax = %d mmol/gDW/h \nfor the %i th coupled reaction.\n', targetMu, vmax, coupledReactionID);
else
fprintf('Applying a factor of %d to all uptake limits results in = %d 1/h.\n', vmax, targetMu);
fprintf('vmax values for coupled reactions:\n');
for i = 1:numberOfReactions
fprintf('%i: %d\n', i, metabolicModel.coupledReactions.vmax(i));
end
end
fprintf('\nFluxes at limit:\n');
for i = 1:length(limitedFluxesIDs)
switch metabolicModel.FBAsolver
case 1
fprintf("ID %i (%s): %d ", limitedFluxesIDs(i), metabolicModel.CNAmodel.reacID(limitedFluxesIDs(i),:), flux(limitedFluxesIDs(i)));
case 2
fprintf("ID %i (%s): %d ", limitedFluxesIDs(i), metabolicModel.COBRAmodel.rxnNames{limitedFluxesIDs(i)}, flux(limitedFluxesIDs(i)));
otherwise
error('In estimateVmax: unknown FBA solver')
end
if ismember(limitedFluxesIDs(i), metabolicModel.exchangeReactions.ReacID)
fprintf("(Exchange reaction, ");
if ismember(limitedFluxesIDs(i), metabolicModel.coupledReactions.ReacID)
fprintf("coupled)");
else
fprintf("uncoupled)");
end
end
fprintf('\n');
end
fprintf('\n');
%likely not really necessary; just for peace of mind
metabolicModel.coupledReactions.vmax = orig.vmax;
metabolicModel.coupledReactions.ks = orig.ks;
if growthLimiting == false
switch metabolicModel.FBAsolver
case 1
metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID) = orig.lb;
metabolicModel.CNAmodel.reacMax(metabolicModel.exchangeReactions.ReacID) = orig.ub;
case 2
metabolicModel.COBRAmodel.lb(metabolicModel.exchangeReactions.ReacID) = orig.lb;
metabolicModel.COBRAmodel.ub(metabolicModel.exchangeReactions.ReacID) = orig.ub;
otherwise
error('In estimateVmax: unknown FBA solver')
end
end
end
function [filteredExchangeTable] = filterCmpndExchangeTable(exchangeTable, mode, threshold)
%FILTERCMPNDEXCHANGEMATRIX Filter table of metabolite exchange according to
%specific criteria
% mode == 0: remove compounds with zero entries only
% mode == 1: true exchange
% mode == 2: only uptake/competition
% mode == 3: only production
% threshold indicates a value, below which fluxes are neglected
rowsToKeep = [];
tmp = table2array(exchangeTable);
tmp(abs(tmp) < threshold) = 0;
tmp2 = array2table(tmp, 'VariableNames', exchangeTable.Properties.VariableNames, 'RowNames', exchangeTable.Properties.RowNames);
exchangeTable = tmp2;