...
 
Commits (3)
......@@ -61,19 +61,19 @@ 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.*
### Simulation Output
Two files are generated at the end of the simulation with a date and time stamp in the filename indicating the start of the simulation. Both files hold Matlab data structures. The file `*_restartInit.mat` records the final state of the simulator and can be used as the initial conditions to continue the simulation in a subsequent run of µbialSim. The other file holds the simulated trajectory in the Matlab structure trajectory. The fields `time`, `compounds`, `biomass`, and `mu` hold the time, compound concentrations, biomass concentrations, and specific growth rates for each integration step. The field `FBA` stores data for each FBA model, including the temporal dynamics of all metabolic fluxes, and the mass balance for all exchange reactions.
Two files are generated at the end of the simulation with a date and time stamp in the filename indicating the start of the simulation. Both files hold Matlab data structures. The file `*_restartInit.mat` records the final state of the simulator and can be used as the initial conditions to continue the simulation in a subsequent run of µbialSim. The other file holds the simulated trajectory in the Matlab structure trajectory. The fields `time`, `compounds`, `biomass`, and `mu` hold the time, compound concentrations, biomass concentrations, and specific growth rates for each integration step. The field `FBA` stores data for each FBA model, including the temporal dynamics of all metabolic fluxes, and the mass balance for all exchange reactions. The field `limitingFluxes` stores over time those reactions, for which thecurrent flux is identical to its current flux boundary (upper or lower). These are typically the growth-limiting fluxes.
Additionally, once the simulation is finished, the trajectory is automatically visualized in three Matlab figures detailing the evolution of compound and biomass concentrations as well as exchange fluxes for all species.
## Built With
......@@ -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 - January 2020 - 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')
......
This diff is collapsed.
This diff is collapsed.
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;
%% iterate across table rows (== compounds)
for i = 1:size(exchangeTable, 1)
keepIt = false;
switch mode
case 0
if sum(abs(exchangeTable{i,:})) ~= 0.0
keepIt = true;
end
case 1
if min(exchangeTable{i,:}) < 0 && max(exchangeTable{i,:}) > 0
keepIt = true;
end
case 2
if min(exchangeTable{i,:}) < 0 && max(exchangeTable{i,:}) <= 0
keepIt = true;
end
case 3
if min(exchangeTable{i,:}) >= 0 && max(exchangeTable{i,:}) > 0
keepIt = true;
end
otherwise
error('Unknown mode.');
end
if keepIt
rowsToKeep(end+1) = i;
end
end
%% set up table
filteredExchangeTable = exchangeTable(rowsToKeep, :);
%% don't include species with zero values only
for i = size(filteredExchangeTable, 2):-1:1
if sum(abs(filteredExchangeTable{:,i})) == 0
filteredExchangeTable(:,i) = [];
end
end
%writetable(filteredExchangeTable, 'filename.txt', 'WriteRowNames', true, 'Delimiter', '\t')
function [exchangeTable, reportTime] = getCmpndExchangeTable(trajectory, time)
%GETCMPEXCHANGEMATRIX Summary of this function goes here
% Detailed explanation goes here
%% get index to desired report time
if trajectory.time(end-1) < time %"-1" as for the last index, there is no corresponding flux!
exchangeTable = 0;
reportTime = -1;
return
end
i = 1;
while trajectory.time(i) < time
i = i + 1;
end
if i > 1 && (time - trajectory.time(i-1)) < (trajectory.time(i) - time)
i = i - 1;
end
index = i;
reportTime = trajectory.time(index);
% truncate too long model names (Matlab restriction)
for i = 1:length(trajectory.modelNames)
if length(trajectory.modelNames{i}) > 63
trajectory.modelNames{i} = trajectory.modelNames{i}(1:63);
trajectory.modelNames{i}(63) = 'x';
end
end
%% set up exchange matrix, entries have unit mM/h, negative values indicate uptake
exchangeMatrix = zeros(length(trajectory.compoundNames), length(trajectory.modelNames));
for i = 1:length(trajectory.modelNames)
exchangeMatrix(:, i) = mapExchangeToReactorCompounds(getFluxes(trajectory.FBA(i).fluxes(index, :), trajectory.FBA(i).coupledReactions.ReacID) * trajectory.biomass(index, i) .* transpose(trajectory.FBA(i).coupledReactions.SecretionSense), trajectory.FBA(i).reactorCompoundIDs, length(trajectory.compoundNames));
end
%% prepare exchange table
exchangeTable = array2table(exchangeMatrix, 'VariableNames', trajectory.modelNames, 'RowNames', trajectory.compoundNames);
end
\ No newline at end of file
function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, compounds, biomass, inflowBiomass, timePars)
function [ dbio, dcompound, flux, success, resultLimitingFluxes ] = getDy(metabolicModel, compounds, biomass, solverPars)
%GETDY Compute the current rate of change for biomass and reactor compounds
%based on the current reactor state
% Detailed info
......@@ -14,13 +14,13 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
end
success = 1;
dbio = 0.0;
dcompound = 0.0 * ones(1, length(metabolicModel.coupledReactions.ReacID));
dcompound = zeros(1, length(metabolicModel.coupledReactions.ReacID));
return
end
% SET UPTAKE LIMITS
for i = 1:length(metabolicModel.coupledReactions.ReacID)
limit = getLimitFlux(metabolicModel, i, compounds(metabolicModel.reactorCompoundIDs(i)));
limit = getMaxUptakeFlux(metabolicModel, i, compounds(metabolicModel.reactorCompoundIDs(i)));
if metabolicModel.coupledReactions{i, 'SecretionSense'} > 0
switch metabolicModel.FBAsolver
case 1
......@@ -57,7 +57,14 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
error('In getDy: unknown FBA solver')
end
if success == false || flux(metabolicModel.biomassReac) < timePars.minimalGrowth
% not a simulation, but just called by estimateVmax.m
if biomass == -1
dbio = 0;
dcompound = 0;
return
end
if success == false || flux(metabolicModel.biomassReac) < solverPars.minimalGrowth
disp('No growth conditions.');
switch metabolicModel.FBAsolver
case 1
......@@ -70,7 +77,7 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
% perform parsimonious FBA, fix growth rate and minimize internal fluxes,
% however not allowing reactions to switch direction
originalSum=(sum(abs(flux)));
if timePars.dopFBA
if solverPars.dopFBA
% record current constraints etc.
switch metabolicModel.FBAsolver
case 1
......@@ -168,7 +175,7 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
% Minimize deviation to previous flux
if (timePars.doMin2PrevFlux == 1 && sum(abs(metabolicModel.lastFlux)) ~= 0) % previous flux available, try to minimize deviation
if (solverPars.doMin2PrevFlux == 1 && sum(abs(metabolicModel.lastFlux)) ~= 0) % previous flux available, try to minimize deviation
originalDev = sum(abs(flux-metabolicModel.lastFlux));
% retain current state
switch metabolicModel.FBAsolver
......@@ -238,13 +245,13 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
case 1
if isnan(metabolicModel.CNAconstraints(i))
%metabolicModel.CNAconstraints(i) = flux(i); % fixing does not work unfortunately
metabolicModel.CNAmodel.reacMin(i) = flux(i) - timePars.fluxTolerance;
metabolicModel.CNAmodel.reacMax(i) = flux(i) + timePars.fluxTolerance;
metabolicModel.CNAmodel.reacMin(i) = flux(i) - solverPars.fluxTolerance;
metabolicModel.CNAmodel.reacMax(i) = flux(i) + solverPars.fluxTolerance;
metabolicModel.CNAmodel.objFunc(i) = 0;
end
case 2
metabolicModel.COBRAmodel.lb(i) = flux(i) - timePars.fluxTolerance;
metabolicModel.COBRAmodel.ub(i) = flux(i) + timePars.fluxTolerance;
metabolicModel.COBRAmodel.lb(i) = flux(i) - solverPars.fluxTolerance;
metabolicModel.COBRAmodel.ub(i) = flux(i) + solverPars.fluxTolerance;
metabolicModel.COBRAmodel.c(i) = 0;
end
end
......@@ -345,7 +352,7 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
disp('Minimization to previous flux did not work out. Relaxing restrictions ...')
end
tryValue = timePars.minRelaxValue;
tryValue = solverPars.minRelaxValue;
switch metabolicModel.FBAsolver
case 1
myUB = metabolicModel.CNAmodel.reacMax;
......@@ -354,7 +361,7 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
myUB = metabolicModel.COBRAmodel.ub;
myLB = metabolicModel.COBRAmodel.lb;
end
while (~mysuccess && tryValue < timePars.maxRelaxValue && metabolicModel.FBAsolver == 2)
while (~mysuccess && tryValue < solverPars.maxRelaxValue && metabolicModel.FBAsolver == 2)
fprintf('Trying %d\n', tryValue)
for i = 1:numberOfReactions
if myLB(i) ~= myUB(i)
......@@ -437,5 +444,17 @@ function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, com
% dcompound(i) = 0;
% end
end
% recording limiting fluxes
if solverPars.recordLimitingFluxes == 1
limitingFluxes = getLimitedFluxes(flux, metabolicModel, solverPars)';
fixedFluxes = find(limitingFluxes == 100);
limitedFluxes = find(abs(limitingFluxes) == 1);
resultLimitingFluxes = [fixedFluxes; limitedFluxes];
resultLimitingFluxes(:, 2) = [flux(fixedFluxes); flux(limitedFluxes)];
resultLimitingFluxes(:, 3) = [ones(length(fixedFluxes), 1); zeros(length(limitedFluxes), 1)];
else
resultLimitingFluxes = 0;
end
end
function [ dydt ] = getDyWrapper(t, y, metabolicModel, reactor, timePars)
function [ dydt ] = getDyWrapper(y, metabolicModel, reactor, solverPars, workerModels)
%GETDYWRAPPER This is a wrapper function for getDy which is used when the
% Matlab ODE solver are used
% Detailed explanation goes here
......@@ -10,18 +10,39 @@ numberOfSpecies = length(metabolicModel);
biomass = y(1:numberOfSpecies);
compounds = y((numberOfSpecies+1):end);
for i = 1:numberOfSpecies
metabolicModel(i).lastFlux = ODEresult.metabolicModel(i).lastFlux;
if solverPars.parallel == 1
for i = 1:length(ODEresult.metabolicModel)
myODEresult_metabolicModel_lastFlux{i} = ODEresult.metabolicModel(i).lastFlux;
end
spmd
for i = 1:length(workerModels)
workerModels(i).lastFlux = myODEresult_metabolicModel_lastFlux{labindex+(i-1)*numlabs};
end
end
else
for i = 1:numberOfSpecies
metabolicModel(i).lastFlux = ODEresult.metabolicModel(i).lastFlux;
end
end
if timePars.parallel == 1
parfor i = 1:numberOfSpecies
[ dbio(i), dcompoundRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), reactor, transpose(compounds), biomass(i), reactor.biomassInflow(i), timePars);
end
if solverPars.parallel == 1
%ticBytes(gcp);
spmd
for i = 1:length(workerModels)
[mydbio(i), mydcompoundRAW{i}, myflux{i}, mysuccess(i)] = getDy(workerModels(i), transpose(compounds), 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};
end
else
for i = 1:numberOfSpecies
[ dbio(i), dcompoundRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), reactor, transpose(compounds), biomass(i), reactor.biomassInflow(i), timePars);
end
for i = 1:numberOfSpecies
[ dbio(i), dcompoundRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), transpose(compounds), biomass(i), solverPars);
end
end
for i = 1:numberOfSpecies
......
function [exchangeFluxTable, limitedFluxesIDs] = getExchangeFluxes(metabolicModel, flux, solverPars, varargin)
%UNTITLED Report on the exchange fluxes given a flux distribution, mark
%coupled reactions, mark reactions at limit
% Detailed explanation goes here
if nargin > 4
error("Too many arguments in getExchangeFluxes()!");
end
exchangeFluxTable = metabolicModel.exchangeReactions;
% mapping back reversible reactions
switch metabolicModel.FBAsolver
case 1
numberOfReactions = length(metabolicModel.CNAmodel.objFunc);
case 2
numberOfReactions = length(metabolicModel.COBRAmodel.c);
otherwise
error('In getExchangeFluxes: unknown FBA solver')
end
if isfield(metabolicModel, 'reversibleReacIDs')
numberOfReactions = numberOfReactions - length(metabolicModel.reversibleReacIDs);
for i = 1:length(metabolicModel.reversibleReacIDs)
flux(metabolicModel.reversibleReacIDs(i)) = flux(metabolicModel.reversibleReacIDs(i)) + flux(numberOfReactions + i);
end
end
fluxes = flux(1:numberOfReactions);
FluxValues = getFluxes(fluxes, metabolicModel.exchangeReactions.ReacID);
if nargin == 4
compounds = varargin{1};
[limitedFluxes, evalLB, evalUB] = getLimitedFluxes(fluxes, metabolicModel, solverPars, compounds);
else
[limitedFluxes, evalLB, evalUB] = getLimitedFluxes(fluxes, metabolicModel, solverPars);
end
limitedFluxesIDs = find(limitedFluxes);
ProductionConsumption = FluxValues .* metabolicModel.exchangeReactions.SecretionSense;
AtFluxLimit = limitedFluxes(metabolicModel.exchangeReactions.ReacID)';
coupledReac = zeros(1, length(fluxes))';
coupledReac(metabolicModel.coupledReactions.ReacID) = 1;
IsCoupled = coupledReac(metabolicModel.exchangeReactions.ReacID);
exchangeFluxTable = addvars(exchangeFluxTable, FluxValues);
exchangeFluxTable = addvars(exchangeFluxTable, ProductionConsumption);
exchangeFluxTable = addvars(exchangeFluxTable, AtFluxLimit);
exchangeFluxTable = addvars(exchangeFluxTable, IsCoupled);
exchangeFluxTable = addvars(exchangeFluxTable, evalLB(metabolicModel.exchangeReactions.ReacID), 'NewVariableNames', 'evalLB');
exchangeFluxTable = addvars(exchangeFluxTable, evalUB(metabolicModel.exchangeReactions.ReacID), 'NewVariableNames', 'evalUB');
switch metabolicModel.FBAsolver
case 1
lb = metabolicModel.CNAmodel.reacMin(metabolicModel.exchangeReactions.ReacID);
ub = metabolicModel.CNAmodel.reacMax(metabolicModel.exchangeReactions.ReacID);
case 2
lb = metabolicModel.COBRAmodel.lb(metabolicModel.exchangeReactions.ReacID);
ub = metabolicModel.COBRAmodel.ub(metabolicModel.exchangeReactions.ReacID);
otherwise
error('In getExchangeFluxes: unknown FBA solver')
end
exchangeFluxTable = addvars(exchangeFluxTable, lb);
exchangeFluxTable = addvars(exchangeFluxTable, ub);
end
function [ limitFlux ] = getLimitFlux(metabolicModel, i, metabConcentration)
%GETLIMITFLUX Compute current uptake limits based on Monod kinetics.
% Detailed explanation goes here
if metabolicModel.coupledReactions{i, 'inhibited'} == 1
limitFlux = 0.0;
else
limitFlux = metabolicModel.coupledReactions{i, 'vmax'};
limitFlux = limitFlux * -1 * metabolicModel.coupledReactions{i, 'SecretionSense'};
limitFlux = limitFlux * metabConcentration;
if limitFlux ~= 0 % avoid division by zero if conc is zero and ks is zero
limitFlux = limitFlux / (metabolicModel.coupledReactions{i, 'ks'} + metabConcentration);
end
end
end
function [ result ] = getFluxLimits( flux, metabolicModel, compounds )
%GETFLUXLIMITS Identify reactions, for which the current flux is identical
function [ result, lb, ub ] = getLimitedFluxes( flux, metabolicModel, solverPars, varargin )
%GETLIMITEDFLUXES Identify reactions, for which the current flux is identical
%to a set flux boundary
% Fluxes at lower limit are indicated by -1, at upper limit by +1, fixed
% reactions are marked as 100.
if nargin > 4
error("Too many arguments in getLimitedFluxes()!");
end
switch metabolicModel.FBAsolver
case 1
......@@ -13,7 +16,7 @@ switch metabolicModel.FBAsolver
lb = metabolicModel.COBRAmodel.lb;
ub = metabolicModel.COBRAmodel.ub;
otherwise
error('getFluxLimits: Unkown model type, aborting')
error('getLimitedFluxes: Unkown model type, aborting')
end
% map back to reversible reactions
......@@ -26,17 +29,19 @@ if isfield(metabolicModel, 'reversibleReacIDs')
lb(metabolicModel.reversibleReacIDs(i)) = lb(numberOfReactions + i);
end
end
% get current uptake limits
for i = 1:length(metabolicModel.coupledReactions.ReacID)
limit = getLimitFlux(metabolicModel, i, compounds(metabolicModel.reactorCompoundIDs(i)));
if metabolicModel.coupledReactions{i, 'SecretionSense'} > 0
lb(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
else
ub(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
% get current uptake limits
if nargin == 4
compounds = varargin{1};
for i = 1:length(metabolicModel.coupledReactions.ReacID)
limit = getMaxUptakeFlux(metabolicModel, i, compounds(metabolicModel.reactorCompoundIDs(i)));
if metabolicModel.coupledReactions{i, 'SecretionSense'} > 0
lb(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
else
ub(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
end
end
end
end
% for fixed CNA reactions
if metabolicModel.FBAsolver == 1
......@@ -51,11 +56,11 @@ end
% find non-zero reactions at the upper or lower limit
result = [];
for i = 1:numberOfReactions
if flux(i) ~= 0
if flux(i) == lb(i)
if abs(flux(i)) > solverPars.fluxTolerance % only check non-zero fluxes
if abs(flux(i) - lb(i)) < solverPars.fluxTolerance
result(i) = -1;
else
if flux(i) == ub(i)
if abs(flux(i) - ub(i)) < solverPars.fluxTolerance
result(i) = 1;
else
result(i) = 0;
......
function [ maxFlux ] = getMaxUptakeFlux(metabolicModel, i, metabConcentration)
%GETMAXUPTAKEFLUX Compute current uptake limits based on Monod kinetics.
% Detailed explanation goes here
if metabolicModel.coupledReactions{i, 'inhibited'} == 1
maxFlux = 0.0;
else
maxFlux = metabolicModel.coupledReactions{i, 'vmax'};
maxFlux = maxFlux * -1 * metabolicModel.coupledReactions{i, 'SecretionSense'};
maxFlux = maxFlux * metabConcentration;
if maxFlux ~= 0 % avoid division by zero if conc is zero and ks is zero
maxFlux = maxFlux / (metabolicModel.coupledReactions{i, 'ks'} + metabConcentration);
end
end
end
......@@ -40,7 +40,7 @@ else % consumption & producing processes present
targetConcentration = (topBound + lowBound) / 2;
consumptionRate = 0;
for i = 1:length(consumers)
consRate = getLimitFlux(metabolicModel(consumers(i)), reacIndexStr(i), targetConcentration);
consRate = getMaxUptakeFlux(metabolicModel(consumers(i)), reacIndexStr(i), targetConcentration);
consRate = consRate * biomass(consumers(i));
consumptionRate = consumptionRate + consRate;
end
......
function [ status ] = integrator_output(t, y, flag, metabolicModel, timePars)
function [ status ] = integrator_output(t, y, flag, metabolicModel, solverPars)
%INTEGRATOR_OUTPUT Collecting results when using Matlab ODE solver
% Detailed explanation goes here
......@@ -11,7 +11,7 @@ function [ status ] = integrator_output(t, y, flag, metabolicModel, timePars)
status = 0;
return
case 'done'
t = timePars.tend; % last time this function is called, t is empty
t = solverPars.tend; % last time this function is called, t is empty
disp('Done.');
status = 0;
return
......@@ -45,7 +45,7 @@ function [ status ] = integrator_output(t, y, flag, metabolicModel, timePars)
ODEresult.mu(n, i) = ODEresult.FBA(i).fluxes(n, metabolicModel(i).biomassReac);
end
fprintf('\nTime %f (%2.1f%% done)\n\n', t(1), t(1) / timePars.tend * 100);
fprintf('\nTime %f (%2.1f%% done)\n\n', t(1), t(1) / solverPars.tend * 100);
for k = 2:size(t, 2)
......@@ -56,7 +56,7 @@ function [ status ] = integrator_output(t, y, flag, metabolicModel, timePars)
ODEresult.FBA(i).fluxes(n, :) = ODEresult.FBA(i).fluxes(n-1, :);
end
fprintf('\n(Time %f (%2.1f%% done))\n\n', t(k), t(k) / timePars.tend * 100);
fprintf('\n(Time %f (%2.1f%% done))\n\n', t(k), t(k) / solverPars.tend * 100);
end
status = 0;
......
function [ trajectory ] = microbialSimMain(scenarioID)
%microbialSimMain v1.0.0
%microbialSimMain v1.1.0
% This is the main routine of the dFBA simulator �bialSim. Different
% simulation scenarios can be set up by configuring "scenarioID"s. Out of
% the box, �bialSim can simulate two microbiomes:
% the box, �bialSim can simulate three batch cultures:
% scenarioID == 1: a Methanococcus maripaludis single culture (Richards et al., 2016)
% scenarioID == 2: a syntrophic methanogenesis biculture (Hamilton et al., 2015)
% scenarioID == 3: a 773 species human gut microbiome (or selection
% thereof) (Magnusdottir et al., 2017)
% scenarioID == 3: 8-species human gut consortium (SIHUMIx) with models taken from the
% collection of 773 human gut microbiome species (Magnusdottir et al., 2017)
%% PARAMETER
solverPars.FBAsolver = 2; % Select FBA solver, set to 1 for CellNetAnalyzer
% and to 2 for COBRA Toolbox
sovlerPars.tend = 1; % Simulation end time (in hours)
sovlerPars.timeStepSize=0.002; % Default simulation step size (in hours)
solverPars.tend = 1; % Simulation end time (in hours)
solverPars.timeStepSize=0.002; % Default simulation step size (in hours)
solverPars.saveLoadedModelsToFile = 0; % Set to 1 to save loaded SBML
% models in Matlab format for speeding up subsequent simulation runs.
% Models are stored in "loadedModels<date><time>.mat".
solverPars.readInitialStateFrom = ''; % Provide trajectory file name to
% continue a previous simulation run ("simulatedTrajectory_*.mat_restartInit.mat").
solverPars.parallel = 0; % Set to 1 to compute FBA models in parallel using
% Matlab's parfor command. Note that model data need to be tansferred
% to workers at each individual iteration again and again, limiting
% efficiency.
% Matlab's spmd environment
solverPars.maxWorkers = 12; % maximum number of workers to recruit for
% parallel pool
solverPars.dopFBA = 1; % Minimize sum of total flux in each FBA computation
% (pFBA). Set to 0 to deactivate.
solverPars.doMin2PrevFlux = 1; % Minimize deviation to flux distribution of
% previous time step for each model. Set to 0 to deactivate.
solverPars.fluxTolerance = 1e-6; % During minimization of deviation to
% previous flux, allow for this deviation for unaltered fluxes.
% previous flux, allow for this deviation for unaltered fluxes. Also
% used as the tolerance when reporting limiting fluxes.
solverPars.minimalGrowth = 1e-6; % If predicted growth rate is below this
% threshold, assume that no growth is possible
solverPars.solverType = 0; % Set to 0 for augmented forward Euler method,
......@@ -48,6 +49,8 @@ solverPars.maxDeviation = 5.0; % For compounds of high demand, reduce
% time step size so that compound change is not bigger than this value
% (specified in %). This helps to avoid oscillatory behavior. To disable
% this feature and speed up simulation, set to "inf"
solverPars.recordLimitingFluxes = 1; % record during the simulation which
% fluxes are at their boundary
solverPars.recording = 0; % Set to 1 to save reactor state at each iteration
% to separate files.
......@@ -101,7 +104,7 @@ switch solverPars.FBAsolver
end
cd(char(placesToLookForCobraToolbox(i)))
initCobraToolbox
initCobraToolbox(false); % no updating
cd(currDir);
else
changeCobraSolver('glpk', 'LP');
......@@ -113,6 +116,7 @@ end
%% SCENARIO SPECIFIC SETUP
clearvars models
workerModels = 0;
switch scenarioID
case 1 % M. maripaluids monoculture
......@@ -126,6 +130,12 @@ switch scenarioID
models(1) = model;
reactor = reactorDefinition_iMR539;
%Examplary usage of function estimateVmax:
%[myVmax, myFlux, myExchangeFluxes] = estimateVmax(models(1), reactor, solverPars, 2, 0.0875);
%trajectory = myExchangeFluxes;
%return
case 2 % syntrophic propionate degrading binary consortium
solverPars.tend = 1.0;
solverPars.timeStepSize = 0.002;
......@@ -139,17 +149,18 @@ switch scenarioID
iMhu428Model = parametrizeFBAmodel_iMhu428(iMhu428Model);
models(1) = iSfu648Model;
models(2) = iMhu428Model;
reactor = reactorDefinition_synProp;
case 3 % 773 gut microbiome
solverPars.tend = 1.0;
solverPars.timeStepSize = 0.002;
% consider all species
speciesToConsider = 1:773;
%speciesToConsider = 1:773;
% consider only selected species
%speciesToConsider = [1, 2, 3];
%speciesToConsider = [46, 125, 164, 178, 245, 270, 362, 482];
speciesToConsider = [46, 125, 164, 178, 245, 270, 362, 482];
[models, externalCompounds] = prepareFBAmodel_773Agora(speciesToConsider, 'models\AGORA-1.01-Western-Diet\', solverPars.FBAsolver);
solverPars.saveLoadedModelsToFile = 1; % look for "loadedModels_*.mat" in current directory, can be used in subsequent runs for speed-up (see next line)
......@@ -175,11 +186,21 @@ switch scenarioID
end
if solverPars.parallel == 1
myPool = parpool(min(length(models), solverPars.maxWorkers));
global CBT_LP_SOLVER
solver = CBT_LP_SOLVER;
parfor i = 1:length(models)
changeCobraSolver(solver, 'LP');
end
% each worker gets his share of FBA models
% (memory requirements could be halfed if models were directly created
% within workers)
workerModels = Composite();
for i = 1:numel(workerModels)
workerModels{i} = models(i:numel(workerModels):length(models));
end
end
if solverPars.saveLoadedModelsToFile == 1
......@@ -192,7 +213,11 @@ if solverPars.saveLoadedModelsToFile == 1
end
tic
trajectory = dFBASimulator(models, reactor, solverPars.solverType, solverPars);
if solverPars.parallel == 1
trajectory = dFBASimulator(models, reactor, solverPars.solverType, solverPars, workerModels);
else
trajectory = dFBASimulator(models, reactor, solverPars.solverType, solverPars, 0);
end
elapsedSeconds = toc;
elapsedHours = elapsedSeconds / 60 / 60;
......@@ -203,17 +228,23 @@ save(solverPars.trajectoryFile, 'trajectory', 'reactor', 'elapsedHours', '-v7.3'
figure
subplot(2,1,1)
plot(trajectory.time, trajectory.biomass)
legend(trajectory.modelNames)
legend(trajectory.modelNames, 'Interpreter', 'none')
title('Biomass')
xlabel('Time (hours)')
ylabel('Concentration (gDW/L)')
subplot(2,1,2)
plot(trajectory.time, trajectory.compounds)
legend(trajectory.compoundNames)
legend(trajectory.compoundNames, 'Interpreter', 'none')
title('Compounds')
xlabel('Time (hours)')
ylabel('Concentration (mmol/L)')
plotTrajectoryCmp(trajectory, 'Simulation')
plotExchangeFluxes(trajectory)
if solverPars.parallel == 1
%delete(gcp('nocreate'))
delete(myPool);
end
end
......@@ -31,7 +31,7 @@ function [ f ] = plotTrajectoryCmp( varargin )
xlabel('Time (hours)')
ylabel('Concentration offset (gDW/L)')
title('Biomass')
legend([myLegend{:}])
legend([myLegend{:}], 'Interpreter', 'none')
hold off
for i = 1:(numberOfPlots - 1)
......
......@@ -29,7 +29,8 @@ switch FBAsolver
case 1
ReactionIDName = CNAmodel.reacID(IDs, :);
case 2
ReactionIDName = COBRAmodel.rxns(IDs, :);
%ReactionIDName = COBRAmodel.rxns(IDs, :);
ReactionIDName = COBRAmodel.rxnNames(IDs, :);
end
%prohibit acetate uptake
......@@ -50,7 +51,8 @@ switch FBAsolver
case 1
ReacNames = CNAmodel.reacID(IDs, :);
case 2
ReacNames = COBRAmodel.rxns(IDs, :);
%ReacNames = COBRAmodel.rxns(IDs, :);
ReacNames = COBRAmodel.rxnNames(IDs, :);
end
coupledReactions = table(IDs', cellstr(ReacNames), exchangeReactions.SecretionSense(strtrim(cellstr(num2str(IDs'))')), vmax', ks', inhibitedCoupledReactions' ,'VariableNames',{'ReacID' 'ReacName' 'SecretionSense' 'vmax' 'ks' 'inhibited'}, 'RowNames', strtrim(cellstr(num2str(IDs'))'));
......
......@@ -11,7 +11,7 @@ function [ reactor ] = reactorDefinition_773Agora( numberOfSpecies, externalComp
reactor.compounds = externalCompounds';
reactor.compoundsInit = ones(1, numberOfCompounds);
reactor.compoundsInit = 0.01 * ones(1, numberOfCompounds);
reactor.biomassInit = 0.1 * ones(1, numberOfSpecies);
......
function [exchangeMatrix] = writeCmpndExchangeTableToFile(exchangeTable, filename)
%WRITECMPNDEXCHANGETABLETOFILE Summary of this function goes here
% Detailed explanation goes here
dimension = sum(size(exchangeTable));
exchangeMatrix = zeros(dimension, dimension);
for i = 1:size(exchangeTable, 1)
for j = 1:size(exchangeTable, 2)
if exchangeTable{i, j} < 0
exchangeMatrix(i, j) = -exchangeTable{i,j};
else
exchangeMatrix(j + size(exchangeTable, 1), i + size(exchangeTable, 2)) = exchangeTable{i,j};
end
end
end
stringsToReplace = ["[", "]"];
replacementStrings = ["_", "_"];
cmpNames = exchangeTable.Properties.RowNames;
spcsNames = exchangeTable.Properties.VariableNames;
% leading numbers not allowed as variable name in table
cmpNames = regexprep(cmpNames, '^([0-9])', 'C_$1');
for i = 1:length(stringsToReplace)
cmpNames = strrep(cmpNames, stringsToReplace(i), replacementStrings(i));
spcsNames = strrep(spcsNames, stringsToReplace(i), replacementStrings(i));
end
spcsNames = transpose(spcsNames);
newRowNames = cat(1, cmpNames, spcsNames);
newColNames = cat(1, spcsNames, cmpNames);
%% positive integer values required for circos plots
%exchangeMatrix = exchangeMatrix / min(exchangeMatrix(exchangeMatrix > 0));
fullExchangeTable = array2table(exchangeMatrix, 'VariableNames', cellstr(newColNames), 'RowNames', cellstr(newRowNames));
writetable(fullExchangeTable, strcat(filename, '.txt'), 'WriteRowNames', true, 'Delimiter', '\t')
end