function [dxdt,mmax,Pmax,bP,m,Kd,Qin] = EthanolModel(t,x,u)
%Definition of model parameters
%Kinetic parameters
a1 = 0.05;    % Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50;    % Growth-associated parameter for ethanol production, [-]
AP1 = 6.0;    % Activation energy parameter for ethanol production, [oC]
AP2 = 20.3;   % Activation energy parameter for ethanol production, [oC]
b1 = 0.035;   % Parameter in the exponential expression of the maximum specific growth rate expression, [oC-1]
b2 = 0.15;    % Parameter in the exponential expression of the growth inhibitory ethanol concentration expression, [oC-1]
b3 = 0.40;    % Parameter in the exponential expression of the specific death rate expression,[oC-1]
c1 = 0.38;    % Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29;    % Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00;    % Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0;    % Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0;    % Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0;    % Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90;   % Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90;   % Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025;  % Basal specific cellular biomass death rate, [h-1]
KdT = 30.00;  % Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5;      % Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005; % Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05; % Maximum specific oxygen consumption rate, [h-1]

%Metabolic parameters
YPS = 0.51; % Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97; % Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53; % Theoretical yield of biomass on glucose, [gX gS-1]

%Physicochemical and thermodynamic parameters
Chbr = 4.18;     % Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18;      % Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000; % Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20;       % Reference temperature, [oC]
KH = 200;        % Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792;       % Oxygen compressibility factor, [-]
R = 0.082;       % Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 ;     % Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000;     % Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 ;     % Density of the fermentation broth, [g L-1]
rhoc = 1000 ;    % Density of the cooling agent, [g L-1]
MO = 32.0;       % Molecular weight of oxygen (O2), [g mol-1]

%Bioreactor design data
AT = 1.0;       % Bioreactor heat transfer area, [m2]
V = 1800;       % Bioreactor working volume, [L]
Vcj = 50;       % Cooling jacket volume, [L]
Ogasin = 0.305; % Oxygen concentration in airflow inlet, [g L-1]

%Definition of model variables
%State variables
Xt = x(1);      % Total cellular biomass, [g L-1]
Xv = x(2);      % Viable cellular biomass, [g L-1]
S = x(3);       % Substrate/Glucose concentration, [g L-1]
P = x(4);       % Product/Ethanol concentration, [g L-1]
Oliq = x(5);    % Dissolved oxygen concentration, [g L-1]
Ogas = x(6);    % Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x(7);       % Temperature in the bioreactor, [oC]
Tc = x(8);      % Temperature of the cooling agent in the jacket, [oC]
Vl = x(9);      % Culture volume in the bioreactor, [L]
Sf_cum = x(10); % Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x(11);   % Batch time, [h]

%Initial Conditions
Xt0 = u(10);     % Initial total cellular biomass, [g L-1]
Xv0 = u(11);     % Initial viable cellular biomass, [g L-1]
S0 = u(12);      % Initial substrate/Glucose concentration, [g L-1]
P0 = u(13);      % Initial product/Ethanol concentration, [g L-1]
Oliq0 = u(14);   % Initial Dissolved oxygen concentration, [g L-1]
Ogas0 = u(15);   % Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T0 = u(16);      % Initial Temperature in the bioreactor, [oC]
Tc0 = u(17);     % Initial Temperature of the cooling agent in the jacket, [oC]
Vl0 = u(18);     % Initial Culture volume in the bioreactor, [L]
Sf_cum0 = u(19); % Initial Cumulative substrate/glucose fed to the bioreactor, [g]
Time0 = u(20);   % Initial batch time, [h]

%Control variables
Qin0 = u(1);        % Volumetric inflow rate, [l/h-1]
Qin= Qin0 + 15*heaviside(t-5) + 5*heaviside(t-10) - 6*heaviside(t-20) - 14*heaviside(t-35);
Xtin = u(2);        % Total biomass concentration in the bioreactor feed, [g L-1]
Xvin = u(3);        % Viable biomass concentration in the bioreactor feed, [g L-1]
Qe = u(4);          % Volumetric outflow rate, [l/h-1]
Sin = u(5);         % Substrate/Glucose concentration in bioreactor feed, [g L-1]
Fc = u(6);          % Cooling agent inlet volumetric flowrate, [L h-1]
Fair = u(7);        % Airflow inlet volumetric flowrate, [L h-1]
Tin = u(8);         % Temperature of bioreactor feed, [oC]
Tcin = u(9);        % Temperature of cooling agent inlet, [oC]

% Definition of model equations
% Kinetic rates
% -----------------------------
% Specific growth rate, [h-1] 
mmax = ((a1*(T-k1))*(1-exp(b1*(T-k2))))^2;
Pmax = Pmaxb + PmaxT/(1-exp(-b2*(T-k3)));
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+exp(-(100-S)/1)); % Specific growth rate, [h-1] 
if m1 >= 0
m = m1;
else
m=0.0;
end
% Non-growth-associated ethanol specific production rate, [h-1]
if S > 0
bP = c1 * exp(-AP1/T) - c2 * exp(-AP2/T) ; % Non-growth-associated ethanol specific production rate, [h-1]
else
bP = 0.0;
end
qP = aP*m + bP;
% Ethanol consumption specific rate
qS = m/YXS + qP/YPS;
% Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq);
% Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+exp(-b3*(T-k4)));
% Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH;
% Oxygen mass transfer coefficient
kla = kla0*1.2^(T-20);
% Volume of the gas phase in the bioreactor
Vg = V - Vl;

%Material balances
%-----------------
% Volume of liquid culture
dVl = Qin - Qe;
% Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt);
% Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv);
% Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv;
% Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv;
% Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv;
% Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg;

% Energy balances
%---------------
% Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr;
% Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc;

% Yields & Productivity
%---------------------
% Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin;
dTime = 1;

% Definition of state derivatives vector
% State derivatives
dxdt = [dXt; dXv; dS; dP; dOliq; dOgas; dT; dTc; dVl; dSf_cum; dTime];
return