pha_cont

function f=pha_cont(t,y)

%Pure culture PHA accumulation model written for CHEE4001
%Damien Batstone 27/8/2019
%Mass basis describes growth of relevant biomass and accumulation
%In fixed volume two-stage chemostat system

%local variable assignment
%reactor 1
Ss_1 = y(1) ;   %substrate
Snh_1 = y(2);   %ammonia
Xb_1 = y(3);    %bacteria
Xphb_1 = y(4);  %PHB

%reactor 2
Ss_2 = y(5) ;   %substrate
Snh_2 = y(6);   %ammonia
Xb_2 = y(7);    %bacteria
Xphb_2 = y(8);  %PHB

%Inputs
qin_1 = 78161;          %L/d fresh feed to R1
qin_2 =  16332;         %L/d fresh feed to R2
Ss_in = 206;            %g/L glucose feed concentration
Snh_in_1 = 3;         %g/L ammonia feed concentration to R1

%Parameters
%physical
V=[78161,94493] ;        %L
%V=[10000,15000] ;        %L

%Biochemical original units
mumax_gs = 0.4321*24;      %d-1	Maximum biomass growth rate
mumax_ps = 0.2628*24;      %d-1	Maximum PHA production rate
mumax_gp = 0.1566*24;      %d-1    Maximum biomass growth rate
KS_s_g	= 1.2;       %kg/m3	Sat coeff for growth
KI_s_g = 16.728;    %kg/m3   Sat inhib coeff for growth
KS_s_p	= 4.13;      %kg/m3	Sat coeff for PHA production
KS_nh_g	= 0.254;	%kg/m3	Ammonia sat coeff for growth
KI_nh_g = 1.691;    %kg/m3    Ammonia inhib coeff for growth
%KI_nh_p	= 80;	    %kg/m3	Ammonia inhib coeff for PHA uptake
KI_nh_p	= 0.2576;	    %kg/m3	Ammonia inhib coeff for PHA uptake
Xmax = 68;       %g/L    Maximum concentration of Biomass
PHAmax	= 0.75;      %max fraction of PHA content
theta = 5.8;        %Inhibitory effect index

%local constitutive eqs
frac_PHA_1 = Xphb_1/(Xb_1+Xphb_1);
frac_PHA_2 = Xphb_2/(Xb_2+Xphb_2);

%stoichiometry from Excel
%Masters 1 values (fed batch)
%v = [-2.083333333         -0.806293019	-0.150442478	1	0	1.108652901
%-1.886792453	-0.338160012	0          0	1	0.72078397];

%2S-CSTR values
v = [-1.779878826	-0.482608211	-0.150442478	1	0	0.66358629
-2.777777778	-1.288544358	0	0	1	2.027562446];


%rates
% reactor 1
r1_1 = mumax_gs*Xb_1*(Ss_1/(KS_s_g+Ss_1))*(Snh_1/(KS_nh_g+Snh_1))
%r1_1 = mumax_gs * (1-((Xb_1/Xmax)^theta)) * (Ss_1/(KS_s_g+Ss_1+((Ss_1^2)/KI_s_g))) * (Snh_1/(KS_nh_g+Snh_1+((Snh_1^2)/KI_nh_g)))
r2_1 = mumax_ps * Xb_1*Ss_1/(KS_s_p+Ss_1)*KI_nh_p/(KI_nh_p+Snh_1)*PHAmax/(PHAmax+frac_PHA_1)

% reactor 2
r1_2 = mumax_gs*Xb_2*Ss_2/(KS_s_g+Ss_2)*Snh_2/(KS_nh_g+Snh_2)
%r1_2 = mumax_gs * (1-((Xb_2/Xmax)^theta)) * (Ss_2/(KS_s_g+Ss_1+((Ss_1^2)/KI_s_g))) * (Snh_2/(KS_nh_g+Snh_2+((Snh_2^2)/KI_nh_g)))
r2_2 = mumax_ps*Xb_2*Ss_2/(KS_s_p+Ss_2)*KI_nh_p/(KI_nh_p+Snh_2)*PHAmax/(PHAmax+frac_PHA_2)

% SS  SNH     XB	XPHB  Add SO2 SCO2 if you need them.

%derivs
f(1) = qin_1/V(1)*(Ss_in-Ss_1) + v(1,1)*r1_1 + v(2,1)*r2_1;
%skip O2
%f(9) = KLa_O*(O_in-O)*(mumax_gs
f(2) = qin_1/V(1)*(Snh_in_1-Snh_1) + v(1,3)*r1_1 + v(2,3)*r2_1;
f(3) = qin_1/V(1)*(0-Xb_1) + v(1,4)*r1_1 + v(2,4)*r2_1;
f(4) = qin_1/V(1)*(0-Xphb_1) + v(1,5)*r1_1 + v(2,5)*r2_1;
%skip co2
f(5) = qin_1/V(2)*(Ss_1-Ss_2)+qin_2/V(2)*(Ss_in-Ss_2)+v(1,1)*r1_2+v(2,1)*r2_2;
%skip O2
f(6) = qin_1/V(2)*(Snh_1-Snh_2)+qin_2/V(2)*(0-Snh_2)+v(1,3)*r1_2+v(2,3)*r2_2;
f(7) = qin_1/V(2)*(Xb_1-Xb_2)+qin_2/V(2)*(0-Xb_2)+v(1,4)*r1_2+v(2,4)*r2_2;
f(8) = qin_1/V(2)*(Xphb_1-Xphb_2)+qin_2/V(2)*(0-Xphb_2)+v(1,5)*r1_2+v(2,5)*r2_2;
%skip co2

f=f';
end