This is my code:
function simulate_sscf()
% Set initial conditions and reaction time
initial_conditions = [30; 0; 0; 7.5; 3; 0]; % initial concentrations of cellulose, cellobiose, glucose, xylose, biomass and ethanol
tspan = [0 36]; % duration of simulation
[t, concentrations] = ode45(@rate_equation, tspan, initial_conditions);
% Plot the Results
plot(t, concentrations)
xlabel('Time (hr)')
ylabel('Concentration (g/L)')
legend('Cellulose','Cellobiose','Glucose','Xylose','Biomass', 'Ethanol');
end
function dxdt = rate_equation(~,x)
global k1r E1B Rs S0 K1IG K1IG2 G XY K1IX k2r E2B K2IG2 K2IG K2IX k3r E2F K3M K3IG K3IX K1ad E1TZ E1max K2ad E2TZ E2max
C = x(1); % Cellulose concentration
G2 = x(2); % Cellobiose concentration
SG = x(3); % Glucose concentration
SX = x(4); % Xylose concentration
X = x(5); % Biomass concentration
P = x(6); % Product concentration
% Define constants for Saccharification
% Substrate reactivity is an assumption based on Figure 2 of ref.paper
K1ad=0.4;
K2ad=0.1;
E1max=0.06;
E2max=0.01;
Ea=-5540;
k1r=22.3;
k2r=7.18;
k3r=285.5;
K1IG2=0.015;
K1IG=0.1;
K1IX=0.1;
K1IE=0.15;
K2IG2=132;
K2IG=0.04;
K2IX=0.2;
K3M=24.3;
K3IG=3.9;
K3IX=201;
XY=0;
Rs=1.987; % where C0 is the initial concentration of cellulose
R=1.987;
E1F=1/2*(-1/K1ad + E1TZ-E1max*C + sqrt((-1/K1ad + E1TZ - E1max*C)^2 + 4*E1TZ/K1ad));
E1B=(E1max*K1ad*E1F*C)/(1+K1ad*E1F);
E2F=1/2*(-1/K2ad + E2TZ - E2max*C + sqrt((-1/K2ad + E2TZ - E2max*C)^2 + 4*E2TZ/K2ad));
E2B=(E2max*K2ad*E2F*C)/(1+K2ad*E2F);
% Define constants for glucose fermentation
mu_mG = 0.662; % Maximum specific growth rate in glucose fermentation
KsG = 0.565; % Monod constant, for growth on glucose
KiG = 283.7; % Inhibition constant, for growth on glucose
betaG = 1.29; % constants in product inhibition model in glucose fermentation
PmG = 95.4; % Ethanol concentration above which cells don't grow in glucose
vmG = 2.005; % Maximum specific rate of glucose formation
Ks2G = 1.342; % Monod constant, for product formation from glucose
Ki2G = 4890; % Inhibition constant, for product formation from glucose
gammaG = 1.42; % Maximum specific rate of glucose formation
Pm2G = 103; % Ethanol concentration above which cells do not produce ethanol in glucose fermentation
YEtG_G = 0.470; % product yield constant (g-ethanol/g-glucose)
YXG_G = 0.115; % cell yield constant from glucose (g-cells/g-substrate)
mG = 0.097; % maintenance coefficient in glucose fermentation
% Define constants for xylose fermentation
mu_mX = 0.190; % Maximum specific growth rate in xylose fermentation
KsX = 3.4; % Monod constant, for growth on xylose
KiX = 18.1; % Inhibition constant, for growth on xylose
betaX = 1.036; % constant in product inhibition model in xylose fermentation
PmX = 59.04; % Ethanol concentration above which cells don't grow in xylose
vmX = 0.250; % Maximum specific rate of xylose formation
Ks2X = 3.4; % Monod constant, for product formation from glucose
Ki2X = 81.3; % Inhibition constant, for product formation from glucose
gammaX = 0.608; % Maximum specific rate of glucose formation
Pm2X = 60.2; % Ethanol concentration above which cells do not produce ethanol in glucose fermentation
YEtX_X = 0.400; % product yield constant (g-ethanol/g-xylose)
YXX_X = 0.162; % cell yield constant from xylose (g-cells/g-substrate)
mX = 0.067; % maintenance coefficient in xylose fermentation
n = 0.5; % weighting factor for glucose consumption
% Define rates of reaction
% cellulose to cellobiose
r(1)=(k1r*E1B*Rs*C)/(1+(G2/K1IG2)+(SG/K1IG)+(SX/K1IX)+(P/K1IE));
% cellulose to glucose
r(2)=(k2r*(E1B+E2B)*Rs*C)/(1+(G2/K2IG2)+(SG/K2IG)+(SX/K2IX));
% cellobiose to glucose
r(3)=(k3r*E2F*G2)/K3M*(1+(SG/K3IG)+(SX/K3IX))+G2;
% Biomass Production from Glucose
mu_G = mu_mG * (SG / (KsG + SG + (SG^2/KiG))) * (1 - (P/PmG)^betaG);
% Biomass Production from Xylose
mu_X = mu_mX * (SX / (KsX + SX + (SX^2/KiX))) * (1 - (P/PmX)^betaX);
% Total Biomass Production from Both
mu_T = n*mu_G + (1-n)*mu_X;
% Ethanol Formation from Glucose
vG = vmG * (SG / (Ks2G + SG + (SG^2/Ki2G))) * (1 - (P/Pm2G)^gammaG);
% Ethanol Formation from Xylose
vX = vmX * (SX / (Ks2X + SX + (SX^2/Ki2X))) * (1 - (P/Pm2X)^gammaX);
% Total Ethanol Formation from both
vT = vG + vX;
% Define rate equations
dCdt = -r(1)-r(2); % Rate of change of cellulose
dG2dt = 1.056*r(1)-r(3); % Rate of change of cellobiose
dXdt = mu_T * X; % Rate of change of biomass
dPdt = vT * X; % Rate of change of product
dSGdt = 1.111*r(2)+1.053*r(3)-(1/YEtG_G) * vG * X; % Rate of change of glucose in SSCF
dSXdt = -(1/YEtX_X) * vX * X; % Rate of change of xylose in SSCF
dxdt = [dCdt; dG2dt; dSGdt; dSXdt; dXdt; dPdt];
end
Any help will be greatly appreciated!
