The UFZ services GitLab and Mattermost will be unavailable on Monday, October 25 from 06:00 AM to 08:00 AM due to maintenance work.

Commit 2b1a22da authored by Florian Centler's avatar Florian Centler
Browse files

Initial microbialSim version

parent 391b28c3
......@@ -631,7 +631,7 @@ to attach them to the start of each source file to most effectively
state the exclusion of warranty; and each file should have at least
the "copyright" line and a pointer to where the full notice is found.
microbialSim
microbialSimDev
Copyright (C) 2019 UMBSysBio
This program is free software: you can redistribute it and/or modify
......@@ -652,7 +652,7 @@ Also add information on how to contact you by electronic and paper mail.
If the program does terminal interaction, make it output a short
notice like this when it starts in an interactive mode:
microbialSim Copyright (C) 2019 UMBSysBio
microbialSimDev Copyright (C) 2019 UMBSysBio
This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
......
# microbialSim
# µbialSim
This is µbialSim, a dynamic Flux-Balance-Analysis simulator for complex microbial communities.
This is µbialSim (pronounced *microbialsim*), a dynamic Flux-Balance-Analysis-based simulator for complex microbial communities. Batch and chemostat operation can be simulated. Simulation results detail the temporal evolution of compound concentrations and biomass concentrations for each microbial species in the bioreactor, and internal flux distributions. For low diversity microbiomes Matlab's ODE solvers can be used, while for more complex communities an augmented forward Euler method is implemented in µbialSim.
Developed by Florian Centler and Denny Popp at UFZ Leipzig, Germany.
\ No newline at end of file
## Getting Started
The latest version of µbialSim can be obtained via git:
```
git clone https://git.ufz.de/UMBSysBio/microbialSim.git
```
### Prerequisites
µbialSim is implemented in Matlab and requires the COBRA Toolbox and/or CellNetAnalyzer for FBA computations. The examples make use of the former.
Information on the COBRA Toolbox is available here:
https://opencobra.github.io/cobratoolbox/stable/index.html
Information on CellNetAnalyzer is available here:
https://www2.mpi-magdeburg.mpg.de/projects/cna/cna.html
### Installing
For using the COBRA Toolbox, its location must be specified in the main simulator file `microbialSimMain.m` in lines 87ff, note that multiple locations can be specified, easing operation in different environments:
```
placesToLookForCobraToolbox = {
'C:\Users\centlerf\cobratoolbox', ...
'Y:\Home\centlerf\Projects\MetabolicModeling\FBATools\cobratoolbox', ...
'/data/cobratoolbox' ...
'..\cobratoolbox'
};
```
## Running the Examples
Three examples are included with µbialSim and can be run as follows.
### Example 1
The first example in which batch-culture growth of a single hydrogenotrophic species (*Methanococcus maripaludis*) is simulated can be run with the command:
```
microbialSimMain(1)
```
The FBA model is taken from *M. A. Richards, T. J. Lie, J. Zhang, S. W. Ragsdale, J. A. Leigh, and N. D. Price, “Exploring hydrogenotrophic methanogenesis: A genome scale metabolic reconstruction of Methanococcus maripaludis,” J. Bacteriol., vol. 198, no. 24, pp. 3379–3390, 2016.*
### Example 2
Batch-culture growth of a binary syntrophic, methanogenic community (*Syntrophobacter fumaroxidans* and *Methanospirillium hungatei*) transforming propionate to methane is started by:
```
microbialSimMain(2)
```
The FBA models are taken from *J. J. Hamilton, M. Calixto Contreras, and J. L. Reed, “Thermodynamics and H2 Transfer in a Methanogenic, Syntrophic Community,” PLOS Comput. Biol., vol. 11, no. 7, p. e1004364, 2015.*
### 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:
```
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.
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.
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
* [Matlab](https://www.mathworks.com/products/matlab.html) - Implementation language and numerical ODE solver
* [COBRA Toolbox](https://opencobra.github.io/cobratoolbox/stable/index.html) - For loading SBML models and running FBA computations
* [CellNetAnalyzer](https://www2.mpi-magdeburg.mpg.de/projects/cna/cna.html) - For loading SBML models and running FBA computations
## Contributing
Please read [CONTRIBUTING.md](https://gist.github.com/PurpleBooth/b24679402957c63ec426) for details on our code of conduct, and the process for submitting pull requests to us.
## Versioning
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).
## Authors
* **Florian Centler**, UFZ - Helmholtz Centre for Environmental Research, Leipzig, Germany
* **Denny Popp**, UFZ - Helmholtz Centre for Environmental Research, Leipzig, Germany
Contact: florian.centler@ufz.de
## License
This project is licensed under the GNU General Public License v3.0 - see the [LICENSE](LICENSE) file for details
## Acknowledgments
* This code was made possible by funding from the German Federal Ministry of Education and Research to Florian Centler (e:Bio project McBiogas, FKZ 031A317)
......@@ -11,7 +11,7 @@ initFile = solverPars.readInitialStateFrom;
if isempty(initFile)
startTime = 0;
stronglyConsumedFlag=ones(length(reactor.metabolites),1)*0; % this flags metabolites that in one time step were more consumed than produced leading to time step size reduction to avoid negative concentrations
stronglyConsumedFlag=ones(length(reactor.compounds),1)*0; % this flags compounds that in one time step were more consumed than produced leading to time step size reduction to avoid negative concentrations
else
myInit = load(initFile);
startTime = myInit.recordResult.time;
......@@ -23,7 +23,7 @@ else
end
timeSteps = floor((solverPars.tend - startTime) / solverPars.timeStepSize);
numberOfMetabolites = length(reactor.metabolites);
numberOfCompounds = length(reactor.compounds);
if length(reactor.biomassInit) ~= length(reactor.biomassInflow) || length(metabolicModel) ~= length(reactor.biomassInit)
error('Number of biomasses, inflow, and FBA models mismatches!');
......@@ -50,11 +50,11 @@ end
if isempty(initFile)
result.time(1) = 0;
result.metabolites(1,:) = reactor.metabolitesInit;
result.compounds(1,:) = reactor.compoundsInit;
result.biomass(1,:) = reactor.biomassInit;
else
result.time(1) = myInit.recordResult.time;
result.metabolites(1,:) = myInit.recordResult.metabolites;
result.compounds(1,:) = myInit.recordResult.compounds;
result.biomass(1,:) = myInit.recordResult.biomass;
end
......@@ -70,19 +70,19 @@ if type == 1
options = odeset('OutputFcn', @(t,y,flag) integrator_output(t,y,flag,metabolicModel,solverPars), 'RelTol', solverPars.relTol, 'AbsTol', solverPars.absTol);
%ODE Simulation
if solverPars.nonNegative == 1
options = odeset(options, 'NonNegative', 1:1:(length(metabolicModel)+numberOfMetabolites));
options = odeset(options, 'NonNegative', 1:1:(length(metabolicModel)+numberOfCompounds));
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.metabolitesInit], options);
[fba.t,fba.y] = ode45(@(t,y) getDyWrapper(t,y,metabolicModel,reactor, solverPars), [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.metabolitesInit], options);
[fba.t,fba.y] = ode15s(@(t,y) getDyWrapper(t,y,metabolicModel,reactor, solverPars), [0:solverPars.timeStepSize:solverPars.tend], [reactor.biomassInit reactor.compoundsInit], options);
otherwise
error('Unknown solver.')
end
result = ODEresult;
result.biomass = fba.y(:,1:length(metabolicModel));
result.metabolites = fba.y(:,(length(metabolicModel)+1:end));
result.compounds = fba.y(:,(length(metabolicModel)+1:end));
result.time(2:(timeSteps+1)) = linspace(solverPars.timeStepSize, solverPars.tend, timeSteps)';
else
n = 1;
......@@ -93,45 +93,45 @@ else
while result.time(n) < solverPars.tend
% get FBA predictions for all models
result.metabolites(n+1, :) = result.metabolites(n, :);
result.compounds(n+1, :) = result.compounds(n, :);
result.biomass(n+1, :) = result.biomass(n, :);
if solveFullFBA == 1
clearvars dbio dmetab flux success
clearvars dbio dcompound flux success
if solverPars.parallel == 1
ticBytes(gcp);
parfor i = 1:length(metabolicModel)
currentModel = metabolicModel(i);
[dbio(i), dmetabRAW{i}, flux{i}, success(i)] = getDy(currentModel, reactor, result.metabolites(n, :), result.biomass(n, i), reactor.biomassInflow(i), solverPars);
[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)
else
for i = 1:length(metabolicModel)
[dbio(i), dmetabRAW{i}, flux{i}, success(i)] = getDy(metabolicModel(i), reactor, result.metabolites(n, :), result.biomass(n, i), reactor.biomassInflow(i), solverPars);
[dbio(i), dcompoundRAW{i}, flux{i}, success(i)] = getDy(metabolicModel(i), reactor, result.compounds(n, :), result.biomass(n, i), reactor.biomassInflow(i), solverPars);
end
end
for i = 1:length(metabolicModel)
metabolicModel(i).previousFlux = metabolicModel(i).lastFlux;
metabolicModel(i).lastFlux = flux{i};
dmetab(i,:) = mapExchangeToReactorMetabolites(dmetabRAW{i}, metabolicModel(i).reactorMetaboliteIDs, length(reactor.metabolites));
dcompound(i,:) = mapExchangeToReactorCompounds(dcompoundRAW{i}, metabolicModel(i).reactorCompoundIDs, length(reactor.compounds));
end
if length(metabolicModel) > 1
myDmetab = sum(dmetab);
myDcompound = sum(dcompound);
else
myDmetab = dmetab;
myDcompound = dcompound;
end
dbioTotal = dbio + considerFlow(reactor.biomassInflow, result.biomass(n, :), reactor);
dmetabFlow = considerFlow(reactor.metabolitesInflow, result.metabolites(n, :), reactor);
dmetabTotal = myDmetab + dmetabFlow;
dcompoundFlow = considerFlow(reactor.compoundsInflow, result.compounds(n, :), reactor);
dcompoundTotal = myDcompound + dcompoundFlow;
end
%Update reactor with changes due to biotic activity and flow
result.biomass(n+1, :) = result.biomass(n+1, :) + dbioTotal * myTimeStepSize;
result.metabolites(n+1, :) = result.metabolites(n+1, :) + dmetabTotal * myTimeStepSize;
result.compounds(n+1, :) = result.compounds(n+1, :) + dcompoundTotal * myTimeStepSize;
% check if time step size is ok or needs to be reduced due to negative
% concentrations, or as exchanged, strongly consumed metabolites
% concentrations, or as exchanged, strongly consumed compounds
% deviate too much
if solveFullFBA == 1
......@@ -150,31 +150,31 @@ else
end
end
% do negative metabolite concentrations occur?
for i = 1:numberOfMetabolites
if result.metabolites(n+1, i) <= 0 && result.metabolites(n, i) ~= 0
if result.metabolites(n+1, i) < 0 && result.metabolites(n, i) <= 0
% do negative compound concentrations occur?
for i = 1:numberOfCompounds
if result.compounds(n+1, i) <= 0 && result.compounds(n, i) ~= 0
if result.compounds(n+1, i) < 0 && result.compounds(n, i) <= 0
error('In TimeStepAdjustment: previous concenctration was zero or negative! How was uptake possible? Something is wrong. Aborting')
end
deltaConcentration = result.metabolites(n+1, i) - result.metabolites(n, i);
if result.metabolites(n+1, i) < 0 && deltaConcentration == 0
deltaConcentration = result.compounds(n+1, i) - result.compounds(n, i);
if result.compounds(n+1, i) < 0 && deltaConcentration == 0
error('In TimeStepAdjustment: deltaConcentration is zero. I am confused. Aborting')
end
requiredTimeStepSize = getTimeToSteadyState(metabolicModel, i, result.metabolites(n, i), dmetab(:,i), dmetabFlow(i), deltaConcentration/solverPars.timeStepSize, result.biomass(n,:), reactor, solverPars);
requiredTimeStepSize = getTimeToSteadyState(metabolicModel, i, result.compounds(n, i), dcompound(:,i), dcompoundFlow(i), deltaConcentration/solverPars.timeStepSize, result.biomass(n,:), reactor, solverPars);
if requiredTimeStepSize < myTimeStepSize
myTimeStepSize = requiredTimeStepSize;
end
stronglyConsumedFlag(i) = 1;
else
% avoid fluctuating concentrations for metabolites produced and
% avoid fluctuating concentrations for compounds produced and
% strongly consumed
if stronglyConsumedFlag(i) == 1 && result.metabolites(n, i) ~= 0
if abs(result.metabolites(n+1,i) - result.metabolites(n, i)) > result.metabolites(n,i)/100 * solverPars.maxDeviation;
requiredTimeStepSize = result.metabolites(n, i) / 100.0 * solverPars.maxDeviation / abs(result.metabolites(n+1,i)-result.metabolites(n,i)) * solverPars.timeStepSize;
if stronglyConsumedFlag(i) == 1 && result.compounds(n, i) ~= 0
if abs(result.compounds(n+1,i) - result.compounds(n, i)) > result.compounds(n,i)/100 * solverPars.maxDeviation;
requiredTimeStepSize = result.compounds(n, i) / 100.0 * solverPars.maxDeviation / abs(result.compounds(n+1,i)-result.compounds(n,i)) * solverPars.timeStepSize;
if requiredTimeStepSize < myTimeStepSize
myTimeStepSize = requiredTimeStepSize;
end
% Uncommenting the next two lines speeds up computation, but leads to non-smooth metabolite trajectories. Runtime and accuracy will be in the middle between
% Uncommenting the next two lines speeds up computation, but leads to non-smooth compound trajectories. Runtime and accuracy will be in the middle between
% setting solverPars.maxDeviation to "5.0" (or other finite values) (slow computation, high accurracy) and "inf" (fast computation, low accurracy).
% else
% stronglyConsumedFlag(i) = 0;
......@@ -182,7 +182,7 @@ else
end
end
if stronglyConsumedFlag(i) == 1
if min(dmetab(:,i)) >= 0 || (dmetabFlow(i) <= 0 && max(dmetab(:,i)) <= 0)
if min(dcompound(:,i)) >= 0 || (dcompoundFlow(i) <= 0 && max(dcompound(:,i)) <= 0)
stronglyConsumedFlag(i) = 0;
end
end
......@@ -196,10 +196,10 @@ else
end
end
% for small biomass and metabolite concentrations: set to zero
for i = 1:numberOfMetabolites
if abs(result.metabolites(n+1, i)) < solverPars.myAccuracy
result.metabolites(n+1, i) = 0.0;
% for small biomass and compound concentrations: set to zero
for i = 1:numberOfCompounds
if abs(result.compounds(n+1, i)) < solverPars.myAccuracy
result.compounds(n+1, i) = 0.0;
end
end
for i = 1:length(metabolicModel)
......@@ -219,6 +219,8 @@ else
else
result.FBA(i).fluxes(n, :) = flux{i};
end
%store current specific growth rate
result.mu(n, i) = result.FBA(i).fluxes(n, metabolicModel(i).biomassReac);
end
result.time(n+1) = result.time(n) + myTimeStepSize;
......@@ -248,48 +250,49 @@ for i = 1:length(metabolicModel)
result.FBA(i).ngamReac = metabolicModel(i).ngamReac;
result.FBA(i).exchangeReactions = metabolicModel(i).exchangeReactions;
result.FBA(i).coupledReactions = metabolicModel(i).coupledReactions;
result.FBA(i).reactorMetaboliteIDs = metabolicModel(i).reactorMetaboliteIDs;
result.FBA(i).reactorCompoundIDs = metabolicModel(i).reactorCompoundIDs;
result.FBA(i).modelName = metabolicModel(i).modelName;
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));
end
end
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));
% copy data to coupledReactions table
result.FBA(i).coupledReactions.MassBalance = 0 * ones(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),:));
result.FBA(i).exchangeReactions.finalFluxLimit = fluxLimits(result.FBA(i).exchangeReactions.ReacID)';
result.FBA(i).coupledReactions.finalFluxLimit = fluxLimits(result.FBA(i).coupledReactions.ReacID)';
end
% copy data to coupledReactions table
result.FBA(i).coupledReactions.MassBalance = 0 * ones(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.metabolites(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
result.timePars = solverPars;
result.metaboliteNames = reactor.metabolites;
result.compoundNames = reactor.compounds;
result.modelNames = {};
for i = 1:length(metabolicModel)
if ischar(metabolicModel(i).modelName)
......
function [ dbio, dmetab, flux, success ] = getDy(metabolicModel, reactor, metabolites, biomass, inflowBiomass, timePars)
function [ dbio, dcompound, flux, success ] = getDy(metabolicModel, reactor, compounds, biomass, inflowBiomass, timePars)
%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, dmetab, flux, success ] = getDy(metabolicModel, reactor, metabo
end
success = 1;
dbio = 0.0;
dmetab = 0.0 * ones(1, length(metabolicModel.coupledReactions.ReacID));
dcompound = 0.0 * ones(1, length(metabolicModel.coupledReactions.ReacID));
return
end
% SET UPTAKE LIMITS
for i = 1:length(metabolicModel.coupledReactions.ReacID)
limit = getLimitFlux(metabolicModel, i, metabolites(metabolicModel.reactorMetaboliteIDs(i)));
limit = getLimitFlux(metabolicModel, i, compounds(metabolicModel.reactorCompoundIDs(i)));
if metabolicModel.coupledReactions{i, 'SecretionSense'} > 0
switch metabolicModel.FBAsolver
case 1
......@@ -432,9 +432,9 @@ function [ dbio, dmetab, flux, success ] = getDy(metabolicModel, reactor, metabo
dbio = flux(metabolicModel.biomassReac) * biomass; % g/l/hour
for i = 1:length(metabolicModel.coupledReactions.ReacID)
dmetab(i) = flux(metabolicModel.coupledReactions.ReacID(i)) * metabolicModel.coupledReactions.SecretionSense(i) * biomass;
% if dmetab(i) > 0 % don't feed back production to reactor to test what happens without metabolite exchange!
% dmetab(i) = 0;
dcompound(i) = flux(metabolicModel.coupledReactions.ReacID(i)) * metabolicModel.coupledReactions.SecretionSense(i) * biomass;
% if dcompound(i) > 0 % don't feed back production to reactor to test what happens without compound exchange!
% dcompound(i) = 0;
% end
end
end
......
......@@ -8,7 +8,7 @@ global ODEresult
numberOfSpecies = length(metabolicModel);
biomass = y(1:numberOfSpecies);
metabolites = y((numberOfSpecies+1):end);
compounds = y((numberOfSpecies+1):end);
for i = 1:numberOfSpecies
metabolicModel(i).lastFlux = ODEresult.metabolicModel(i).lastFlux;
......@@ -16,31 +16,31 @@ end
if timePars.parallel == 1
parfor i = 1:numberOfSpecies
[ dbio(i), dmetabRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), reactor, transpose(metabolites), biomass(i), reactor.biomassInflow(i), timePars);
[ dbio(i), dcompoundRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), reactor, transpose(compounds), biomass(i), reactor.biomassInflow(i), timePars);
end
else
for i = 1:numberOfSpecies
[ dbio(i), dmetabRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), reactor, transpose(metabolites), biomass(i), reactor.biomassInflow(i), timePars);
[ dbio(i), dcompoundRAW{i}, flux{i}, success(i) ] = getDy(metabolicModel(i), reactor, transpose(compounds), biomass(i), reactor.biomassInflow(i), timePars);
end
end
for i = 1:numberOfSpecies
dmetab(i,:) = mapExchangeToReactorMetabolites(dmetabRAW{i}, metabolicModel(i).reactorMetaboliteIDs, length(reactor.metabolites));
dcompound(i,:) = mapExchangeToReactorCompounds(dcompoundRAW{i}, metabolicModel(i).reactorCompoundIDs, length(reactor.compounds));
% consider flow (biomasses)
dbio(i) = dbio(i) + considerFlow(reactor.biomassInflow(i), biomass(i), reactor);
end
ODEresult.ODElastFlux = flux;
% consider flow (metabolites)
% consider flow (compounds)
if numberOfSpecies > 1
myDmetab = sum(dmetab);
myDcompound = sum(dcompound);
else
myDmetab = dmetab;
myDcompound = dcompound;
end
dmetab = myDmetab + considerFlow(reactor.metabolitesInflow, metabolites', reactor);
dcompound = myDcompound + considerFlow(reactor.compoundsInflow, compounds', reactor);
dydt = transpose([dbio dmetab]);
dydt = transpose([dbio dcompound]);
end
function [ result ] = getFluxLimits( flux, metabolicModel, metabolites )
function [ result ] = getFluxLimits( flux, metabolicModel, compounds )
%GETFLUXLIMITS 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
......@@ -29,7 +29,7 @@ end
% get current uptake limits
for i = 1:length(metabolicModel.coupledReactions.ReacID)
limit = getLimitFlux(metabolicModel, i, metabolites(metabolicModel.reactorMetaboliteIDs(i)));
limit = getLimitFlux(metabolicModel, i, compounds(metabolicModel.reactorCompoundIDs(i)));
if metabolicModel.coupledReactions{i, 'SecretionSense'} > 0
lb(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
else
......
function [ deltaTime ] = getTimeToSteadyState(metabolicModel, metaboliteID, metabConc, deltaMetab, deltaFlow, slope, biomass, reactor, solverPars)
function [ deltaTime ] = getTimeToSteadyState(metabolicModel, compoundID, metabConc, deltaCompound, deltaFlow, slope, biomass, reactor, solverPars)
%GETTIMETOSTEADYSTATE Summary of this function goes here
% Detailed explanation goes here
......@@ -9,11 +9,11 @@ productionRate = 0;
consumptionRate = 0;
for i = 1:length(metabolicModel)
if deltaMetab(i) > 0
productionRate = productionRate + deltaMetab(i);
else if deltaMetab(i) < 0
if deltaCompound(i) > 0
productionRate = productionRate + deltaCompound(i);
else if deltaCompound(i) < 0
consumers(end + 1) = i;
consumptionRate = consumptionRate + deltaMetab(i);
consumptionRate = consumptionRate + deltaCompound(i);
end
end
end
......@@ -23,7 +23,7 @@ end
reacIndexStr = cell(1,1);
for i = 1:length(consumers)
tmp = find(~(metabolicModel(consumers(i)).reactorMetaboliteIDs - metaboliteID)); % find index at which position value metaboliteID is located (finding the index of the coupled reaction which refers to the reactor metabolite)
tmp = find(~(metabolicModel(consumers(i)).reactorCompoundIDs - compoundID)); % find index at which position value compoundID is located (finding the index of the coupled reaction which refers to the reactor compound)
reacIndexStr{i} = int2str(metabolicModel(consumers(i)).coupledReactions{tmp, 'ReacID'});
end
......@@ -44,7 +44,7 @@ else % consumption & producing processes present
consRate = consRate * biomass(consumers(i));
consumptionRate = consumptionRate + consRate;
end
consumptionRate = consumptionRate + considerFlow(reactor.metabolitesInflow(metaboliteID), targetConcentration, reactor);
consumptionRate = consumptionRate + considerFlow(reactor.compoundsInflow(compoundID), targetConcentration, reactor);
if productionRate + consumptionRate < 0
topBound = targetConcentration;
else
......@@ -60,4 +60,4 @@ else % consumption & producing processes present
end
% compute deltaTime at which steady state concentration is reached
deltaTime = (targetConcentration - metabConc) / slope;
end
\ No newline at end of file
end
......@@ -42,6 +42,7 @@ function [ status ] = integrator_output(t, y, flag, metabolicModel, timePars)
else
ODEresult.FBA(i).fluxes(n, :) = flux{i};
end
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);
......
function [ dcompound ] = mapExchangeToReactorCompounds( dcompoundRAW, reactorMapping, numOfReactorCompounds )
%MAPEXCHANGETOREACTORMETABOLITES Summary of this function goes here
% Detailed explanation goes here
% dcompound(i,:) = mapExchangeToReactorCompounds(dcompoundRAWE, metabolicModel(i).reactorCompoundIDs);
dcompound = zeros(1,numOfReactorCompounds);
for i = 1:length(reactorMapping)
dcompound(reactorMapping(i)) = dcompoundRAW(i);
end
end
function [ dmetab ] = mapExchangeToReactorMetabolites( dmetabRAW, reactorMapping, numOfReactorMetabolites )
%MAPEXCHANGETOREACTORMETABOLITES Summary of this function goes here
% Detailed explanation goes here
%dmetab(i,:) = mapExchangeToReactorMetabolites(dmetabRAW, metabolicModel(i).reactorMetabolitesIDs);
dmetab = zeros(1,numOfReactorMetabolites);
for i = 1:length(reactorMapping)
dmetab(reactorMapping(i)) = dmetabRAW(i);
end
end