In the simulation file the following information must be specified.
• The number of microbial species (nmodel).
• The FBA model names and default bounds. These bounds are equivalent to infinity in the model.
• The reactions in the FBA models that involve exchange fluxes of interest (exID).
• The objectives for each model (C). After maximizing growth, each of the exchange fluxes used in the DRHS file must have an objective.
• The initial conditions (Y0).
• The time of simulation (tspan).
• DFBAlab options:
o Linear program feasibility and optimality tolerance (tol).
o Penalty function tolerance (tolPh1).
o Method: either the primal or the dual LP can be solved (Method).
o Solver: simplex, dual simplex, barrier (solver).
o Do you want the first solve to be displayed? (disp).
o Do you want to use larger phase 1 models? (optionslacks).
o Do you want your simulation to stop if there is an error in the LP section? (error).
o What finite number should be considered as infinity? (Infinity).
• MATLAB integrator options:
o Integrating tolerances
o Mass matrices: useful to solve differential-algebraic equation (DAE) systems.
o Nonnegative variables.
• Plotting section: MATLAB plotting options.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DFBAlab: Dynamic Flux Balance Analysis laboratory %
% Process Systems Engineering Laboratory, Cambridge, MA, USA %
% July 2014 %
% Written by Jose A. Gomez and Kai Höffner %
% %
% This code can only be used for academic purposes. When using this code %
% please cite: %
% %
% Gomez, J.A., Höffner, K. and Barton, P. I. %
% DFBAlab: A fast and reliable MATLAB code for Dynamic Flux Balance %
% Analysis. Submitted. %
% %
% COPYRIGHT (C) 2014 MASSACHUSETTS INSTITUTE OF TECHNOLOGY %
% %
% Read the LICENSE.txt file for more details. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example is based on the High-rate algal pond (HRAP) proposed in
% Buhr, H.O., Miller, S.B. A dynamic model of the high-rate algal-bacterial
% wastewater treatment pond. Water Res. 17, 29-37 (1983), and in Yang, A.
% Modeling and evaluation of CO2 supply and utilization in algal ponds.
% Industrial and Engineering Chemistry Research 50, 11181-11192 (2011). It
% includes a pH balance.
% The algae model used is iRC1080 from Chang, R.L. et al. Metabolic network
% reconstruction of Chlamydnomnas offers insight into light-driven algal
% metabolism. Molecular Systems Biology 7(518) (2011).
% The yeast model is iND750 from Duarte, N.C., Herrgard, M.J., Palsson,
% B.O. Reconstruction and validation of Saccharomyces cerevisiae iND750, a
% fully compartmentalized genome-scale metabolic model. Genome Research
% 14(7), 1298-1309 (2004)

clear all
global CPLEXobjs
nmodel = 2; % Number of models

% Load models. These should be .mat files generated by the COBRA toolbox.
% When generating these files using the COBRA toolbox, a big number is used
% as infinity. This number should be fed to the DB vector (Default bound).

load iND750.mat
model{1} = iND750;
DB(1) = 1000;
load iRC1080.mat
model{2} = iRC1080;
DB(2) = 1000;

%% exID array
% You can either search the reaction names by name or provide them directly
% in the exID array.
% RxnNames = {'EX_glc(e)', 'EX_ac(e)', 'biomass'};
% for i = 1:length(RxnNames)
% [a,exID(i)] = ismember(RxnNames(i),model.rxns);
% end

exID{1}=[428,458,407,420];
exID{2}=[26,28,24,25,1,2,3,4,5,6,7,8,9,10,11,12,13,62,63,64,65,27,81,47,13];

% Lower bounds and upper bounds for these reactions should be provided in
% the RHS code.

% This codes solves the LPs in standard form. Bounds on exchange fluxes in
% the exID array can be modified directly on the first 2*n rows where n is
% the number of exchange fluxes. Order will be lower bound, upper bound,
% lower bound, upper bound in the same order as exID.
%
% NOTE: All bounds on fluxes in the exID arrays are relaxed to -Inf and
% + Inf. These bounds need to be updated if needed in the RHS file.

%% Cost vectors
% Usually the first cost vector will be biomass maximization, but it can
% be any other objective. The CPLEX objects will minimize by default.
% Report only nonzero elements.
% The structure should be:
% C{model} = struct
% Element C{k}(i) of C, is the cost structure i for model k.
% C{k}(i).sense = +1 for minimize, or -1 for maximize.
% C{k}(i).rxns = array containing the reactions in this objective.
% C{k}(i).wts = array containing coefficients for reactions reported in
% rxns. Both arrays should have the same length.
% Example, if:
% C{k}(i).rxns = [144, 832, 931];
% C{k}(i).wts = [3, 1, -1];
% Then the cost vector for this LP will be:
% Cost{k}(i) = 3*v_144 + v_832 - v_931 (fluxes for model k).
% This cost vector will be either maximized or minimized depending on the
% value of C{k}(i).sense.

% In SBML files, usually production fluxes are positive and uptake fluxes
% are negative. Keep in mind that maximizing a negative flux implies
% minimizing its absolute value.
% Different models can have different number of objectives.

min = 1;
max = -1;
% Yeast
% Maximize growth
C{1}(1).sense = max;
C{1}(1).rxns = [1266];
C{1}(1).wts = [1];
% Glucose
C{1}(2).sense = max;
C{1}(2).rxns = [428];
C{1}(2).wts = [1];
% O2
C{1}(3).sense = max;
C{1}(3).rxns = [458];
C{1}(3).wts = [1];
% CO2
C{1}(4).sense = max;
C{1}(4).rxns = [407];
C{1}(4).wts = [1];
% Ethanol
C{1}(5).sense = max;
C{1}(5).rxns = [420];
C{1}(5).wts = [1];

% Algae
% Maximize growth
C{2}(1).sense = max;
C{2}(1).rxns = [63];
C{2}(1).wts = [1];
% Acetate
C{2}(2).sense = min;
C{2}(2).rxns = [28];
C{2}(2).wts = [1];
% O2
C{2}(3).sense = max;
C{2}(3).rxns = [24,81];
C{2}(3).wts = [1,1];
% CO2
C{2}(4).sense = min;
C{2}(4).rxns = [25];
C{2}(4).wts = [1];

% Initial conditions
% Y1 = Volume (L)
% Y2 = Biomass Yeast (gDW/L)
% Y3 = Biomass Algae (gDW/L)
% Y4 = Glucose (mmol/L)
% Y5 = O2 (mmol/L)
% Y6 = Total Carbon (mmol/L)
% Y7 = Ethanol (mmol/L)
% Y8 = Acetate (mmol/L)
% Y9 = Total Nitrogen (mmol/L)
% Y10 = Penalty
% Y11 = NH4+
% Y12 = NH3
% Y13 = CO2
% Y14 = HCO3-
% Y15 = CO3 -2
% Y16 = H +

Y0 = [140 1.1048 1.8774 0.0140, 0.00065156 1.2211 8.2068 0.0237 0.1643 0 0.1643 2.4476E-5 1.0568 0.1643 2.5842E-6 2.6701E-6]';

% Time of simulation
tspan = [0,24];

% CPLEX Objects construction parameters
tol = 1E-9; % Feasibility, optimality and convergence tolerance for Cplex (tol>=1E-9).
% It is recommended it is at least 2 orders of magnitude
% tighter than the integrator tolerance.
% If problems with infeasibility messages, tighten this
% tolerance.
tolPh1 = tol; % Tolerance to determine if a solution to phaseI equals zero.
% It is recommended to be the same as CPLEX.
Method = 1; % Select 0 to work on dual problem (DEFAULT)
% Select 1 to work on primal problem
solver = 2; % Primal Simplex = 1 (RECOMMENDED FOR Method = 0 ONLY)
% Automatic = 0 (CHOOSE IF PROBLEMS ARE BEING ENCOUNTERED -
% SLOWER)
% Dual Simplex = 2 (RECOMMENDED FOR Method = 1 ONLY)
% Barrier Method = 4
disp = 1; % Set to one if you want the solution of the first LP to be
% displayed on MATLAB screen.
optionslacks = zeros(nmodel,1); % Zero is the default value. It will be faster
% because it's using a reduced amounts of slacks. If it
% doesn't work, set this option to 1.
error = 0; % Set to 1 if you want the simulation to be stopped if the LP is not optimal
% during integration (RECOMMENDED).
Infinity = 9E25; % When solving the primal LP, CPLEX does not like infinity in the rhs or
% lhs vectors. The following number should be a very big number (but
% finite) that will play the role of infinity.

% You can modify the integration tolerances here.
% If some of the flows become negative after running the simulation once
% you can add the 'Nonnegative' option.
% The mass matrix M can help you solve DAEs.
M = [eye(10) zeros(10,6); zeros(6,16)];
options = odeset('AbsTol',1E-6,'RelTol',1E-6,'Mass',M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[CPLEXobjs,INFO,pair] = CplexSetup(model,nmodel,DB,C,exID,Y0,tol,tolPh1,solver,disp,Infinity,Method,optionslacks,error);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;
% Look at MATLAB documentation if you want to change solver.
% ode15s is more or less accurate for stiff problems.
[T,Y] = ode15s(@DRHS,tspan,Y0,options,INFO);
display(toc);
% %
% % % %% Plotting
figure(1)
plot(T,Y(:,2:3))
figure(2)
plot(T,Y(:,4))
figure(3)
plot(T,Y(:,5),T,Y(:,13))
figure(4)
plot(T,Y(:,8))
X = -log10(Y(:,16));
figure(5)
plot(T,X);