Commit 50d389c0 authored by Florian Centler's avatar Florian Centler

This is microbialSim, version 1.0

parent 308c8b66
function [ dflow ] = considerFlow(inflowConc, current, reactor)
%CONSIDERFLOW Compute the change of reactor concentrations due to
%in-/outflow
% The current flow rate is computed based on set flow rate and reactor
% volume, and the rate of change of reactor concentrations due to flow
% conditions is then computed and returned. The function is called for
% biomass concentrations and compound concentrations separately.
flowrate = reactor.flowRate / reactor.volume;
dflow = flowrate * (inflowConc - current);
end
function [ model ] = convertToOnlyIrrevReac( model )
%CONVERTTOONLYIRREVREAC Modifying a model such that it only contains
%irreversible reactions
% For reversible reactions, these reactions are converted to only
% proceeding in the forward direction. The backward reaction is appended
% to the list of reactions. Exchange reactions are not modified.
% Collecting the ids of originally reversible reactions
model.reversibleReacIDs = [];
switch model.FBAsolver
case 1
numberOfReactions = length(model.CNAmodel.objFunc);
case 2
numberOfReactions = length(model.COBRAmodel.c);
otherwise
error('In convertToOnlyIrrecReac: unknown FBA solver')
end
for i = 1:numberOfReactions
% don't modify exchange reactions
if (~ismember(i, model.exchangeReactions.ReacID))
% check whether reversible
switch model.FBAsolver
case 1
isReversible = sign(model.CNAmodel.reacMin(i)) * sign(model.CNAmodel.reacMax(i)) < 0;
case 2
isReversible = sign(model.COBRAmodel.lb(i)) * sign(model.COBRAmodel.ub(i)) < 0;
otherwise
error('In convertToOnlyIrrecReac: unknown FBA solver')
end
if isReversible
model.reversibleReacIDs(length(model.reversibleReacIDs) + 1) = i;
% copy reaction in stoichiometric matrix and bounds and
% optimization coefficient
switch model.FBAsolver
case 1
model.CNAmodel.stoichMat(:,size(model.CNAmodel.stoichMat, 2) + 1) = model.CNAmodel.stoichMat(:,i);
model.CNAmodel.reacMin(end + 1) = model.CNAmodel.reacMin(i);
model.CNAmodel.reacMax(end + 1) = model.CNAmodel.reacMax(i);
model.CNAmodel.objFunc(end + 1) = model.CNAmodel.objFunc(i);
model.CNAconstraints(end + 1) = model.CNAconstraints(i);
model.CNAmodel.numr = model.CNAmodel.numr + 1;
% fix directions & optimization coefficient
model.CNAmodel.reacMin(i) = 0;
model.CNAmodel.reacMax(length(model.CNAmodel.reacMin)) = 0;
if (model.CNAmodel.objFunc(length(model.CNAmodel.objFunc)) ~=0)
error('Case not considered yet!')
end
if ~isnan(model.CNAconstraints(i))
if model.CNAconstraints(i) >= 0
model.CNAconstraints(length(model.CNAconstraints)) = 0;
else
model.CNAconstraints(i) = 0;
end
end
newName = strcat(model.CNAmodel.reacID(i,:),'-reverse');
newName = strcat(newName, repmat('_',1,length(model.CNAmodel.reacID(1,:))-length(newName)));
model.CNAmodel.reacID(end + 1,:) = newName;
case 2
model.COBRAmodel.S(:,size(model.COBRAmodel.S, 2) + 1) = model.COBRAmodel.S(:,i);
model.COBRAmodel.lb(end + 1) = model.COBRAmodel.lb(i);
model.COBRAmodel.ub(end + 1) = model.COBRAmodel.ub(i);
model.COBRAmodel.c(end + 1) = model.COBRAmodel.c(i);
% fix directions & optimization coefficient
model.COBRAmodel.lb(i) = 0;
model.COBRAmodel.ub(length(model.COBRAmodel.lb)) = 0;
if (model.COBRAmodel.c(length(model.COBRAmodel.c)) ~=0)
error('Case not considered yet!')
end
model.COBRAmodel.rxns(end + 1) = strcat(model.COBRAmodel.rxns(i),'-reverse');
otherwise
error('In convertToOnlyIrrecReac: unknown FBA solver')
end
end
end
end
end
This diff is collapsed.
This diff is collapsed.
function [ dydt ] = getDyWrapper(t, y, metabolicModel, reactor, timePars)
%GETDYWRAPPER This is a wrapper function for getDy which is used when the
% Matlab ODE solver are used
% Detailed explanation goes here
global ODEresult
numberOfSpecies = length(metabolicModel);
biomass = y(1:numberOfSpecies);
metabolites = y((numberOfSpecies+1):end);
for i = 1:numberOfSpecies
metabolicModel(i).lastFlux = ODEresult.metabolicModel(i).lastFlux;
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);
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);
end
end
for i = 1:numberOfSpecies
dmetab(i,:) = mapExchangeToReactorMetabolites(dmetabRAW{i}, metabolicModel(i).reactorMetaboliteIDs, length(reactor.metabolites));
% consider flow (biomasses)
dbio(i) = dbio(i) + considerFlow(reactor.biomassInflow(i), biomass(i), reactor);
end
ODEresult.ODElastFlux = flux;
% consider flow (metabolites)
if numberOfSpecies > 1
myDmetab = sum(dmetab);
else
myDmetab = dmetab;
end
dmetab = myDmetab + considerFlow(reactor.metabolitesInflow, metabolites', reactor);
dydt = transpose([dbio dmetab]);
end
function [ result ] = getFluxLimits( flux, metabolicModel, metabolites )
%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
% reactions are marked as 100.
switch metabolicModel.FBAsolver
case 1
lb = metabolicModel.CNAmodel.reacMin;
ub = metabolicModel.CNAmodel.reacMax;
case 2
lb = metabolicModel.COBRAmodel.lb;
ub = metabolicModel.COBRAmodel.ub;
otherwise
error('getFluxLimits: Unkown model type, aborting')
end
% map back to reversible reactions
numberOfReactions = length(lb);
if isfield(metabolicModel, 'reversibleReacIDs')
numberOfReactions = numberOfReactions - length(metabolicModel.reversibleReacIDs);
for i = 1:length(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, metabolites(metabolicModel.reactorMetaboliteIDs(i)));
if metabolicModel.coupledReactions{i, 'SecretionSense'} > 0
lb(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
else
ub(metabolicModel.coupledReactions{i, 'ReacID'}) = limit;
end
end
% for fixed CNA reactions
if metabolicModel.FBAsolver == 1
for i = 1:numberOfReactions
if ~isnan(metabolicModel.CNAconstraints(i))
lb = metabolicModel.CNAconstraints(i);
ub = metablicModel.CNAconstraints(i);
end
end
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)
result(i) = -1;
else
if flux(i) == ub(i)
result(i) = 1;
else
result(i) = 0;
end
end
if result(i) == -1 && lb(i) == ub(i)
result(i) = 100; % indicating a fixed reaction
end
else
result(i) = 0;
end
end
end
function [ fluxes ] = getFluxes( fluxVector, indices )
%GETFLUXES Get flux values for fluxes with index from indices.
% Detailed explanation goes here
fluxes = fluxVector(indices);
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 [ deltaTime ] = getTimeToSteadyState(metabolicModel, metaboliteID, metabConc, deltaMetab, deltaFlow, slope, biomass, reactor, solverPars)
%GETTIMETOSTEADYSTATE Summary of this function goes here
% Detailed explanation goes here
% get production rate & record consuming species
consumers = [];
productionRate = 0;
consumptionRate = 0;
for i = 1:length(metabolicModel)
if deltaMetab(i) > 0
productionRate = productionRate + deltaMetab(i);
else if deltaMetab(i) < 0
consumers(end + 1) = i;
consumptionRate = consumptionRate + deltaMetab(i);
end
end
end
% for each consuming species, get index of coupled reaction
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)
reacIndexStr{i} = int2str(metabolicModel(consumers(i)).coupledReactions{tmp, 'ReacID'});
end
% get steady state concentration
if productionRate + max(deltaFlow, 0) == 0 % only consuming processes present
targetConcentration = 0;
else % consumption & producing processes present
topBound = metabConc;
lowBound = 0;
consumptionRate = - productionRate - 1000; % just to ensure to enter the while loop
oldConsumptionRate = Inf;
while abs(productionRate + consumptionRate) > solverPars.SteadyStateAccuracy
targetConcentration = (topBound + lowBound) / 2;
consumptionRate = 0;
for i = 1:length(consumers)
consRate = getLimitFlux(metabolicModel(consumers(i)), reacIndexStr(i), targetConcentration);
consRate = consRate * biomass(consumers(i));
consumptionRate = consumptionRate + consRate;
end
consumptionRate = consumptionRate + considerFlow(reactor.metabolitesInflow(metaboliteID), targetConcentration, reactor);
if productionRate + consumptionRate < 0
topBound = targetConcentration;
else
lowBound = targetConcentration;
end
if oldConsumptionRate == consumptionRate % rate does not depend on concentration, e.g. if Km = 0!
targetConcentration = 0;
break
else
oldConsumptionRate = consumptionRate;
end
end
end
% compute deltaTime at which steady state concentration is reached
deltaTime = (targetConcentration - metabConc) / slope;
end
\ No newline at end of file
function [ returnValue ] = initCNA
%INITCNA Initializing CellNetAnalyzer
% 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' ...
);
for i = 1:length(directoriesToSearch)
if (exist(directoriesToSearch(i,:), 'dir'))
cwd = pwd();
cd(directoriesToSearch(i,:));
startcna(1);
cd(cwd);
returnValue = 0;
return
end
end
% failure to initialize CNA
returnValue = 1;
end
function [ status ] = integrator_output(t, y, flag, metabolicModel, timePars)
%INTEGRATOR_OUTPUT Collecting results when using Matlab ODE solver
% Detailed explanation goes here
global ODEresult
switch flag
case 'init'
t = t(1); % first time this function is called, t is [0 tmax](?)
disp('Starting.');
status = 0;
return
case 'done'
t = timePars.tend; % last time this function is called, t is empty
disp('Done.');
status = 0;
return
end
for i = 1:length(metabolicModel)
numberOfReversibleReactions = length(metabolicModel(i).reversibleReacIDs);
switch metabolicModel(i).FBAsolver
case 1
numberOfReactions(i) = length(metabolicModel(i).CNAmodel.reacID) - numberOfReversibleReactions;
case 2
numberOfReactions(i) = length(metabolicModel(i).COBRAmodel.rxns) - numberOfReversibleReactions;
end
end
flux = ODEresult.ODElastFlux;
n = ODEresult.counter;
ODEresult.counter = ODEresult.counter + 1;
for i = 1:length(metabolicModel)
ODEresult.metabolicModel(i).lastFlux = flux{i};
if isfield(metabolicModel(i), 'reversibleReacIDs')
ODEresult.FBA(i).fluxes(n, :) = flux{i}(1:numberOfReactions(i));
% add reverse reactions
for j = 1:length(metabolicModel(i).reversibleReacIDs)
ODEresult.FBA(i).fluxes(n, metabolicModel(i).reversibleReacIDs(j)) = ODEresult.FBA(i).fluxes(n, metabolicModel(i).reversibleReacIDs(j)) + flux{i}(numberOfReactions(i)+j);
end
else
ODEresult.FBA(i).fluxes(n, :) = flux{i};
end
end
fprintf('\nTime %f (%2.1f%% done)\n\n', t(1), t(1) / timePars.tend * 100);
for k = 2:size(t, 2)
n = ODEresult.counter;
ODEresult.counter = ODEresult.counter + 1;
for i = 1:length(metabolicModel)
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);
end
status = 0;
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
function [ trajectory ] = microbialSimMain(scenarioID)
%microbialSimMain v1.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:
% 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)
%% 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)
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*").
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.
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.
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,
% 1 for direct approach (using Matlab's ODE solver)
% Only used if augmented forward Euler method is used:
solverPars.myAccuracy = 1e-15; % Metabolite concentrations below this value
% are evaluated to be zero.
solverPars.myBioAccuracy = 1e-15; % Biomass concentrations below this
% vaue are evaluated to be zero.
solverPars.SteadyStateAccuracy = 1e-5; % Allow for this deviation between
% production and consumption rate when computing steady state during
% time step size reduction
solverPars.maxDeviation = 5.0; % For metabolites of high demand, reduce
% time step size so that metabolite 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.recording = 0; % Set to 1 to save reactor state at each iteration
% to separate files.
% Only used if COBRA is used:
solverPars.minRelaxValue = 1e-11; % If minimization of deviation to
% previous flux fails, relax restrictions, starting with this value
% (allow small zero-boundary violations by this amount) ...
solverPars.maxRelaxValue = 1; % ... and stop at this value.
% Only used for direct approach:
solverPars.solver = 2; % Set to 1 for ode45 solver, to 2 for ode15s solver.
solverPars.nonNegative = 0; % Set to 1 to activate NonNegative option for
% solver.
solverPars.relTol = 1e-9; % Set relative tolerance
solverPars.absTol = 1e-9; % Set absolute tolerance
% Filenames
fileSuffix = strrep(char(datetime('now')),':', '_');
modelFile = strcat('loadedModels_', fileSuffix, '.mat');
solverPars.trajectoryFile = strcat('simulatedTrajectory_', fileSuffix, '.mat');
for name = {modelFile, solverPars.trajectoryFile}
if exist(char(name), 'file') == 2
error(strcat('Filename ', char(name), ' already exists. Aborting.'))
end
end
%% FBA TOOL BOX SETUP
% initialize FBA solver
switch solverPars.FBAsolver
case 1
returnValue = initCNA;
case 2
if 2 ~= exist('optimizeCbModel')
currDir = pwd;
% ADD PATH TO YOUR COBRA INSTALLATION HERE (REPLACE OR APPEND AT END OF LIST)
placesToLookForCobraToolbox = {
'C:\Users\centlerf\Nextcloud2\Home\Projects\MetabolicModeling\FBATools\cobratoolbox', ...
'Y:\Home\centlerf\Projects\MetabolicModeling\FBATools\cobratoolbox', ...
'/data/cobratoolbox' ...
'..\cobratoolbox'
};
for i = 1:length(placesToLookForCobraToolbox)
if 7 == exist(char(placesToLookForCobraToolbox(i)), 'dir') % directory exists
break
end
if i == length(placesToLookForCobraToolbox)
error('COBRA Toolbox not found. Aborting.')
end
end
cd(char(placesToLookForCobraToolbox(i)))
initCobraToolbox
cd(currDir);
else
changeCobraSolver('glpk', 'LP');
end
otherwise
error('FBA solver type unkown!')
end
%% SCENARIO SPECIFIC SETUP
clearvars models
switch scenarioID
case 1 % M. maripaluids monoculture
solverPars.tend = 1;
solverPars.timeStepSize = 0.02;
solverPars.solverType = 0; % 0 dFBA or 1 for ODE
model = prepareFBAmodel_iMR539('models\MODELRICHARDS 2016.xml', solverPars.FBAsolver);
model = convertToOnlyIrrevReac(model);
model = parametrizeiMR539(model);
models(1) = model;
reactor = reactorDefinition_iMR539;
case 2 % Reed Consortium
solverPars.tend = 1;
solverPars.timeStepSize = 0.02;
solverPars.solverType = 0; % 0 dFBA or 1 ODE
iSfu648Model = prepareFBAmodel_iSfu648Model('models\S5_Dataset_SDH_Coupling.xml', solverPars.FBAsolver);
iSfu648Model = convertToOnlyIrrevReac(iSfu648Model);
iSfu648Model = parametrizeiSfu648Model(iSfu648Model);
iMhu428Model = prepareFBAmodel_iMhu428Model('models\S3_Dataset.xml', solverPars.FBAsolver);
iMhu428Model = convertToOnlyIrrevReac(iMhu428Model);
iMhu428Model = parametrizeiMhu428Model(iMhu428Model);
models(1) = iSfu648Model;
models(2) = iMhu428Model;
reactor = reactorDefinitionReed;
case 3 % 773 gut microbiome
% consider all species
speciesToConsider = 1:773;
% consider only selected species
%speciesToConsider = [1, 2, 3];
%speciesToConsider = [46, 125, 164, 178, 245, 270, 362, 482];
[models, externalMetabolites] = prepare773Model(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');
solverPars.tend = 0.3;
solverPars.timeStepSize = 0.002;
solverPars.readInitialStateFrom = '';
solverPars.parallel = 0;
solverPars.recording = 0;
% set KS to nonzero
ksvalue = 0.01;
for i = 1:length(models)
for j = 1:length(models(i).coupledReactions.ks)
models(i).coupledReactions.ks(j) = ksvalue;
end
end
reactor = reactorDefinition773(length(models), externalMetabolites);
otherwise
error('Unknown model scenario. Aborting')
end
if solverPars.parallel == 1
global CBT_LP_SOLVER
solver = CBT_LP_SOLVER;
parfor i = 1:length(models)
changeCobraSolver(solver, 'LP');
end
end
if solverPars.saveLoadedModelsToFile == 1
save(modelFile, 'models', 'externalMetabolites', '-v7.3'); % for > 2GB
end
tic
trajectory = dFBASimulator(models, reactor, solverPars.solverType, solverPars);
elapsedSeconds = toc;
elapsedHours = elapsedSeconds / 60 / 60;
save(solverPars.trajectoryFile, 'trajectory', 'reactor', 'elapsedHours', '-v7.3');
%% plotting
% Quick-n-dirty-Plot
figure
subplot(2,1,1)
plot(trajectory.time, trajectory.biomass)
title('Biomass')
xlabel('Time (hours)')
ylabel('Concentration (gDW/L)')
subplot(2,1,2)
plot(trajectory.time, trajectory.metabolites)
title('Compounds')
xlabel('Time (hours)')
ylabel('Concentration (mmol/L)')
plotTrajectoryCmp(trajectory, 'Simulation')
plotExchangeFluxes(trajectory)
end
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
function [ model ] = parametrizeiSfu648Model(model)
%PARAMETRIZE Setting the parameters (NGAM and uptake kinetics)for the model
% Detailed explanation goes here
%% NGAM
NGAMValue = 5.1176; % mmolATP/gDW/h
switch model.FBAsolver
case 1
model.CNAconstraints(model.ngamReac) = NGAMValue;
case 2
model.COBRAmodel.lb(model.ngamReac) = NGAMValue;
model.COBRAmodel.ub(model.ngamReac) = NGAMValue;
otherwise
display('parametrizeModel(): unkown type selected')
end
%% UPTAKE KINETICS
%CO2, H2, Acetate, Methane
model.coupledReactions.vmax = [50 50 50 0]';
model.coupledReactions.ks = [0.1 0.1 0.1 0]';
function [ model ] = parametrizeiMhu428(model)
%PARAMETRIZE Setting the parameters (NGAM and uptake kinetics)for the model
% Detailed explanation goes here
%% NGAM
NGAMValue = 0.6/24.0; % mmolATP/gDW/h
switch model.FBAsolver
case 1
model.CNAconstraints(model.ngamReac) = NGAMValue;
case 2
model.COBRAmodel.lb(model.ngamReac) = NGAMValue;
model.COBRAmodel.ub(model.ngamReac) = NGAMValue;
otherwise
display('parametrizeModel(): unkown type selected')
end
%% UPTAKE KINETICS
% acetate, co2, h2, Formate, methane
model.coupledReactions.vmax = [0 75.7 102.5*2/2*2 955*0 0]';
model.coupledReactions.ks = [1 1 1 1 1]';
function [ model ] = parametrizeiSfu648Model(model)
%PARAMETRIZE Setting the parameters (NGAM and uptake kinetics)for the model
% Detailed explanation goes here
%% NGAM
NGAMValue = 3.36/24.0; % mmolATP/gDW/h
switch model.FBAsolver
case 1
model.CNAconstraints(model.ngamReac) = NGAMValue;
case 2
model.COBRAmodel.lb(model.ngamReac) = NGAMValue;
model.COBRAmodel.ub(model.ngamReac) = NGAMValue;
otherwise
display('parametrizeModel(): unkown type selected')
end
%% UPTAKE KINETICS
%prop, ac, fum, CO2, H2, for
model.coupledReactions.vmax = [37.7 0 0 0 0 0 0]';
model.coupledReactions.ks = [1 1 1 1 1 1 1]';