The following information in the environment file must be specified:
• You can either use the y variables directly, or assign y values to more relevant names such as X for biomass and S for substrates.
• Feed rates and feed compositions.
• Mass transfer expressions.
• You can map your LP fluxes to the v vector or use them directly for dy.
• The dy vector carries the dynamics information. In this case, the last elements of the dy vector contain the DAE information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dy = DRHS(t, y, INFO)

% 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 +

% Assign values from states
Vol = y(1);
X(1) = y(2);
X(2) = y(3);
for i=1:13
S(i) = y(3+i);
end

% Feed rates
Fin = 1;
Fout = 1;

% Biomass Feed concentrations
Xfeed(1) = 0;
Xfeed(2) = 0;

% Mass transfer coefficients
KhO2 = 0.0013;
KhCO2 = 0.035;

% Mass transfer expressions
MT(1) = 0;
MT(2) = 0.6*(KhO2*0.21*1000 - S(2));
MT(3) = 0.58*(KhCO2*0.00035*1000 - S(10));
MT(4) = 0;
MT(5) = 0;
MT(6) = 0;

% Substrate feed concentrations
Sfeed(1) = 15.01;
Sfeed(2) = KhO2*0.21*1000;
Sfeed(3) = KhCO2*0.00035*1000;
Sfeed(4) = 0;
Sfeed(5) = 40;
Sfeed(6) = 0;

% The elements of the flux matrix have the sign given to them by the
% coefficients in the Cost vector in main.
% 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:
% flux(k,i) = 3*v_144 + v_832 - v_931
% The penalty is an array containing the feasibility LP problem optimal
% objective function value for each model.

%% Update bounds and solve for fluxes
INFO.t = t;
[lbx,ubx]=RHS(t,y,INFO);
[ flux,penalty ] = CPLEXObjsSolve( INFO, lbx, ubx);
%%

% Yeast fluxes
for i=1:4
v(1,i) = flux(1,i+1);
end
v(1,5) = 0;
v(1,6) = 0;

% Algae fluxes
v(2,1) = 0;
v(2,2) = flux(2,3);
v(2,3) = flux(2,4);
v(2,4) = 0;
v(2,5) = flux(2,2);

%% Dynamics
dy = zeros(6,1); % a column vector
dy(1) = Fin-Fout; % Volume
dy(2) = flux(1,1)*y(2) + (Xfeed(1)*Fin - y(2)*Fout)/y(1); % Biomass yeast
dy(3) = flux(2,1)*y(3) + (Xfeed(2)*Fin - y(3)*Fout)/y(1); % Biomass algae
for i = 1:5
dy(i+3) = v(1,i)*X(1) + v(2,i)*X(2) + MT(i) + (Sfeed(i)*Fin - S(i)*Fout)/y(1) ;
end

if (S(2)/1000 > KhO2 && dy(3+2)>0)
dy(3+2) = 0;
end

if (S(3)/1000 > KhCO2 && dy(3+3)>0)
dy(3+3) = 0;
end

dy(9) = 0; % Leave total nitrogen constant
dy(10) = penalty(1) + penalty(2); %Penalty function

% Algebraic Equations
Ka = 10^-9.4003;
K1c = 10^-6.3819;
K2c = 10^-10.3767;
Nt = S(6);
Ct = S(3);

x = y(11:16);
F = [-x(1) + Nt*x(6)/(x(6) + Ka) ;
-x(4) + x(1)/(1 + 2*K2c/x(6));
-x(3) + x(6)*x(4)/K1c;
-x(5) + x(4)*K2c/x(6);
-Nt + x(1) + x(2)
-Ct + x(3) + x(4) + x(5)];
dy(11:16) = F;

end