Commit 3aa0879a authored by Florian Centler's avatar Florian Centler

Updating to version v1.1.1

parent 1dc84383
......@@ -222,7 +222,7 @@ else
if abs(result.biomass(n+1, i)) < solverPars.myBioAccuracy
requiredTimeStepSize = deltaTimeToZero;
else
requiredTimeStepSize = deltaTimeToZero / 2.0;
requiredTimeStepSize = deltaTimeToZero / solverPars.biomassReductionFactor;
end
if requiredTimeStepSize < myTimeStepSize
myTimeStepSize = requiredTimeStepSize;
......@@ -306,7 +306,7 @@ else
result.time(n+1) = result.time(n) + myTimeStepSize;
if solverPars.recording == 1
recordFilename = strcat(solverPars.trajectoryFile, '_', int2str(n),'.mat');
writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, 0);
writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, 0, solverPars, workerModels);
end
n = n + 1;
solveFullFBA = 1;
......@@ -315,12 +315,12 @@ else
end
if solverPars.recording == 1 % record final step, there are no fluxes then!
recordFilename = strcat(solverPars.trajectoryFile, '_', int2str(n),'.mat');
writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, 1);
writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, 1, solverPars, workerModels);
end
% record restart file
recordFilename = strcat(solverPars.trajectoryFile, '_restartInit.mat');
writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, 1);
writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, 1, solverPars, workerModels);
end
result.time = transpose(result.time);
......
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);
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/L, 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
......@@ -15,6 +15,7 @@ function [ dbio, dcompound, flux, success, resultLimitingFluxes ] = getDy(metabo
success = 1;
dbio = 0.0;
dcompound = zeros(1, length(metabolicModel.coupledReactions.ReacID));
resultLimitingFluxes = '';
return
end
......@@ -454,7 +455,7 @@ if solverPars.recordLimitingFluxes == 1
resultLimitingFluxes(:, 2) = [flux(fixedFluxes); flux(limitedFluxes)];
resultLimitingFluxes(:, 3) = [ones(length(fixedFluxes), 1); zeros(length(limitedFluxes), 1)];
else
resultLimitingFluxes = 0;
resultLimitingFluxes = '';
end
end
......@@ -4,7 +4,7 @@ function [exchangeFluxTable, limitedFluxesIDs] = getExchangeFluxes(metabolicMode
% Detailed explanation goes here
if nargin > 4
error("Too many arguments in getExchangeFluxes()!");
error('Too many arguments in getExchangeFluxes()!');
end
exchangeFluxTable = metabolicModel.exchangeReactions;
......
......@@ -5,7 +5,7 @@ function [ result, lb, ub ] = getLimitedFluxes( flux, metabolicModel, solverPars
% reactions are marked as 100.
if nargin > 4
error("Too many arguments in getLimitedFluxes()!");
error('Too many arguments in getLimitedFluxes()!');
end
switch metabolicModel.FBAsolver
......
......@@ -3,9 +3,9 @@ function [ returnValue ] = initCNA
% This function calls the init function of CellNetAnalyzer
directoriesToSearch=char( ...
'W:\MetabolicModeling\CellNetAnalyzer\CellNetAnalyzer_for_MATLAB75_or_higher', ...
'C:\Users\centlerf\Nextcloud2\Home\Projects\MetabolicModeling\FBATools\CellNetAnalyzer\CellNetAnalyzer_for_MATLAB75_or_higher', ...
'U:\DBFZ\Florian_Kopplung\FBACoupling\CNA2012.1_for_MATLAB75_or_higher\CellNetAnalyzer_for_MATLAB75_or_higher' ...
'W:/MetabolicModeling/CellNetAnalyzer/CellNetAnalyzer_for_MATLAB75_or_higher', ...
'C:/Users/centlerf/Nextcloud2/Home/Projects/MetabolicModeling/FBATools/CellNetAnalyzer/CellNetAnalyzer_for_MATLAB75_or_higher', ...
'U:/DBFZ/Florian_Kopplung/FBACoupling/CNA2012.1_for_MATLAB75_or_higher/CellNetAnalyzer_for_MATLAB75_or_higher' ...
);
for i = 1:length(directoriesToSearch)
......
......@@ -48,14 +48,23 @@ function [ status ] = integrator_output(t, y, flag, metabolicModel, solverPars)
fprintf('\nTime %f (%2.1f%% done)\n\n', t(1), t(1) / solverPars.tend * 100);
for k = 2:size(t, 2)
n = ODEresult.counter;
ODEresult.counter = ODEresult.counter + 1;
% to have values at each time step ...
for i = 1:length(metabolicModel)
ODEresult.FBA(i).fluxes(n, :) = ODEresult.FBA(i).fluxes(n-1, :);
ODEresult.mu(n, i) = ODEresult.mu(n-1, i);
end
% ... but more correct would be this, but requiring the filtering of
% negative values during plotting.
% for i = 1:length(metabolicModel)
% ODEresult.FBA(i).fluxes(n, :) = -ones(1, size(ODEresult.FBA(i).fluxes, 2));
% ODEresult.mu(n, i) = -1;
% end
fprintf('\n(Time %f (%2.1f%% done))\n\n', t(k), t(k) / solverPars.tend * 100);
end
......
function [ trajectory ] = microbialSimMain(scenarioID)
%microbialSimMain v1.1.0
%microbialSimMain v1.1.1
% 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 three batch cultures:
......@@ -49,6 +49,9 @@ 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.biomassReductionFactor = 2.0; % If negative biomass concentrations
% occur, reduce time step size such that at new time step, biomass
% becomes old biomass divided by this factor.
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
......@@ -88,10 +91,12 @@ switch solverPars.FBAsolver
currDir = pwd;
% ADD PATH TO YOUR COBRA INSTALLATION HERE (REPLACE OR APPEND AT END OF LIST)
placesToLookForCobraToolbox = {
'C:\Users\centlerf\cobratoolbox', ...
'Y:\Home\centlerf\Projects\MetabolicModeling\FBATools\cobratoolbox', ...
'W:/cobratoolbox', ...
'C:/Users/centlerf/cobratoolbox', ...
'Y:/Home/centlerf/Projects/MetabolicModeling/FBATools/cobratoolbox', ...
'/data/cobratoolbox' ...
'..\cobratoolbox'
'../cobratoolbox' ...
'/Users/flori/cobratoolbox'
};
for i = 1:length(placesToLookForCobraToolbox)
......@@ -124,7 +129,7 @@ switch scenarioID
solverPars.timeStepSize = 0.002;
solverPars.solverType = 0; % 0 dFBA or 1 for ODE
model = prepareFBAmodel_iMR539('models\MODELRICHARDS 2016.xml', solverPars.FBAsolver);
model = prepareFBAmodel_iMR539('models/MODELRICHARDS 2016.xml', solverPars.FBAsolver);
model = convertToOnlyIrrevReac(model);
model = parametrizeFBAmodel_iMR539(model);
models(1) = model;
......@@ -141,10 +146,10 @@ switch scenarioID
solverPars.timeStepSize = 0.002;
solverPars.solverType = 0; % 0 dFBA or 1 ODE
iSfu648Model = prepareFBAmodel_iSfu648Model('models\S5_Dataset_SDH_Coupling.xml', solverPars.FBAsolver);
iSfu648Model = prepareFBAmodel_iSfu648Model('models/S5_Dataset_SDH_Coupling.xml', solverPars.FBAsolver);
iSfu648Model = convertToOnlyIrrevReac(iSfu648Model);
iSfu648Model = parametrizeFBAmodel_iSfu648(iSfu648Model);
iMhu428Model = prepareFBAmodel_iMhu428Model('models\S3_Dataset.xml', solverPars.FBAsolver);
iMhu428Model = prepareFBAmodel_iMhu428Model('models/S3_Dataset.xml', solverPars.FBAsolver);
iMhu428Model = convertToOnlyIrrevReac(iMhu428Model);
iMhu428Model = parametrizeFBAmodel_iMhu428(iMhu428Model);
models(1) = iSfu648Model;
......@@ -162,9 +167,9 @@ switch scenarioID
%speciesToConsider = [1, 2, 3];
speciesToConsider = [46, 125, 164, 178, 245, 270, 362, 482];
[models, externalCompounds] = prepareFBAmodel_773Agora(speciesToConsider, 'models\AGORA-1.01-Western-Diet\', solverPars.FBAsolver);
[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)
%load('.\loadedModels_26-Mar-2019 11_34_00.mat');
%load('./loadedModels_26-Mar-2019 11_34_00.mat');
solverPars.readInitialStateFrom = '';
......@@ -191,7 +196,7 @@ if solverPars.parallel == 1
global CBT_LP_SOLVER
solver = CBT_LP_SOLVER;
parfor i = 1:length(models)
changeCobraSolver(solver, 'LP');
changeCobraSolver(solver, 'LP', 1, -1);
end
% each worker gets his share of FBA models
......
......@@ -36,7 +36,7 @@ function [ f ] = plotExchangeFluxes(trajectory)
xlabel('Time (hours)')
ylabel('Flux (mmol/gDW/h)')
title({string(trajectory.modelNames(i)), 'Coupled (non-zero) fluxes'})
title({char(trajectory.modelNames(i)), 'Coupled (non-zero) fluxes'}, 'Interpreter', 'none')
legend(myLegend, 'Interpreter', 'none'); % no super/subscripting
end
subplot(ysize, xsize, panelCounter)
......@@ -59,7 +59,7 @@ function [ f ] = plotExchangeFluxes(trajectory)
plot(time, trajectory.FBA(i).fluxes(:,finalIds))
xlabel('Time (hours)')
ylabel('Flux (mmol/gDW/h)')
title({string(trajectory.modelNames(i)), 'Uncoupled (non-zero) exchange fluxes'})
title({char(trajectory.modelNames(i)), 'Uncoupled (non-zero) exchange fluxes'}, 'Interpreter', 'none')
legend(myLegend, 'Interpreter', 'none'); % no super/subscripting
end
end
......
function [f] = plotTopBiomassCompounds(trajectory, maxNumberOfBiomass, maxNumberOfCompounds)
%PLOTTOPBIOMASSCOMPOUNDS Summary of this function goes here
% Detailed explanation goes here
f = figure
subplot(1,2,1)
%xlim([0 1])
numberOfData = min(maxNumberOfBiomass, size(trajectory.biomass, 2));
colors = colorcube(numberOfData);
[dummy, idx] = sort(trajectory.biomass(end,:), 'descend');
for i=1:numberOfData
hold all
h(i)=plot(trajectory.time,trajectory.biomass(:,idx(i)),'Color',colors(i,:));
end
l = legend(trajectory.modelNames(idx(1:numberOfData)), 'Interpreter', 'none', 'Location', 'southoutside');
title('Biomass')
xlabel('Time (hours)')
ylabel('Biomass concentration (gDW/L)')
subplot(1,2,2)
%xlim([0 1])
numberOfData = min(maxNumberOfCompounds, size(trajectory.compounds, 2));
colors = colorcube(numberOfData);
[dummy, idx] = sort(trajectory.compounds(end,:), 'descend');
for i=1:numberOfData
hold all
h(i)=plot(trajectory.time,trajectory.compounds(:,idx(i)),'Color',colors(i,:));
end
l = legend(trajectory.compoundNames(idx(1:numberOfData)), 'Interpreter', 'none', 'Location', 'southoutside');
title('Compounds')
xlabel('Time (hours)')
ylabel('Concentration (mmol/L)')
end
......@@ -42,7 +42,7 @@ function [ f ] = plotTrajectoryCmp( varargin )
end
xlabel('Time (hours)')
ylabel('Concentration (mmol/L)')
title(varargin{1,1}.compoundNames(i))
title(varargin{1,1}.compoundNames(i), 'Interpreter', 'none')
hold off
end
......
This diff is collapsed.
......@@ -16,8 +16,8 @@ for i = 1:size(exchangeTable, 1)
end
end
stringsToReplace = ["[", "]"];
replacementStrings = ["_", "_"];
stringsToReplace = ['[', ']'];
replacementStrings = ['_', '_'];
cmpNames = exchangeTable.Properties.RowNames;
spcsNames = exchangeTable.Properties.VariableNames;
% leading numbers not allowed as variable name in table
......
function [success] = writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, isLastRecord)
function [success] = writeStateToFile(recordFilename, n, result, metabolicModel, stronglyConsumedFlag, isLastRecord, solverPars, workerModels)
%WRITESTATETOFILE Summary of this function goes here
% Detailed explanation goes here
recordResult.time = result.time(n);
recordResult.compounds = result.compounds(n,:);
recordResult.biomass = result.biomass(n,:);
recordResult.stronglyConsumedFlag = stronglyConsumedFlag;
if solverPars.parallel == 1
for i = 1:numel(workerModels)
tmpWorker = workerModels{i};
for j = 1:numel(tmpWorker)
metabolicModel(i+numel(workerModels)*(j-1)).lastFlux = tmpWorker(j).lastFlux;
metabolicModel(i+numel(workerModels)*(j-1)).previousFlux = tmpWorker(j).previousFlux;
end
end
end
for m = 1:length(metabolicModel)
if exist('result.FBA')
recordResult.FBA(m).fluxes = result.FBA(m).fluxes(n,:);
......@@ -19,4 +29,3 @@ function [success] = writeStateToFile(recordFilename, n, result, metabolicModel,
clear recordResult
success = 1;
end
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment