Plotting Kinetic Model Monod, Haldane Equation

Hello! I am trying to plot graphs (Substrate vs time, Biomass vs time, and Product vs time) from these kinetic equations below:

Where the parameter kinetic are shown in the following table

And the maintenance coefficient of cells (ms) gave a value less than 0.005 g g−1 h−1, it was neglected in the model.

I want to plot this paramaeter to the model kinetic to get plot data like this

where the data shown in the table below

I have tried code like this, but it is still error. Any little suggestion or solutions or comments are very welcomed! Thank you so much!

function CobaLunaFlores_ode23

Rentang = [0:10:80];
C0 = [28.90 0.62 0.07586];
Miumax=0.145;
Ks=1.6607;
Ki=50.3151;
Kp=8798.15;
gamma=3.6658;
Yxs=0.3941;
alfa=371.45;
betha=1.5502;
ms=0.005;

[t,C]=ode23(@LunaFlores,Rentang,C0)

plot(t,C)

end

function dCdt=LunaFlores(t,C)


dCdt(1,:)=((Miumax*C(1))/(Ks+C(1)+(C(1)^2/Ki))*(1-P(1)/Kp)^gamma)*C(2);
dCdt(2,:)= -C(2)*((((Miumax*C(1))/(Ks+C(1)+(C(1)^2/Ki))*(1-P(1)/Kp)^gamma)/Yxs)+ms);
dCdt(3,:)= C(2)*(alfa*((Miumax*C(1))/(Ks+C(1)+(C(1)^2/Ki))*(1-P(1)/Kp)^gamma)+betha);
end

Answer 1

[Star Strider]

I got it to run correctly, however I made no other changes (and there need to be other changes), leaving those for you.

It is necessary to pass all the parameters to ‘LunaFlores’ as extra parameters (see the documentation section on Passing Extra Parameters for details). Also, ‘P(1)’ should be ‘C(3)’. I made those corrections and corrected the resulting function call in the ode23 call, however it is still not workign correctly.

Those problems are fairly obvious, since the code does not match the posted symbolic code, and likely needs to. Code ‘C(1)’ as ‘X’, ‘C(2)’ as ‘S’, and ‘C(3) as ‘P’ throughout the code. With those changes and other corrections required for the code to match the posted symbolic differential equations, it should work.

CobaLunaFlores_ode23

function CobaLunaFlores_ode23

Rentang = linspace(0, 80, 50);
C0 = [28.90 0.62 0.07586];
Miumax=0.145;
Ks=1.6607;
Ki=50.3151;
Kp=8798.15;
gamma=3.6658;
Yxs=0.3941;
alfa=371.45;
betha=1.5502;
ms=0.005;

[t,C]=ode23(@(t,C)LunaFlores(t,C,Miumax,gamma,ms,betha,Kp,Ki,alfa,Ks,Yxs),Rentang,C0);

figure
plot(t,C)
% set(gca, 'YScale','log')
legend('C_1','C_2','C_3', 'Location','best')

end

function dCdt=LunaFlores(t,C,Miumax,gamma,ms,betha,Kp,Ki,alfa,Ks,Yxs)


dCdt(1,:)=((Miumax*C(1))/(Ks+C(1)+(C(1)^2/Ki))*(1-C(3)/Kp)^gamma)*C(2);
dCdt(2,:)= -C(2)*((((Miumax*C(1))/(Ks+C(1)+(C(1)^2/Ki))*(1-C(3)/Kp)^gamma)/Yxs)+ms);
dCdt(3,:)= C(2)*(alfa*((Miumax*C(1))/(Ks+C(1)+(C(1)^2/Ki))*(1-C(3)/Kp)^gamma)+betha);
end

I will provide further help if you require it, however the changes needed to get this running correctly are obvious and should be easy to fix.

.

[Nadya Shalsabila salman]

hello @Star Strider, thanks for the answer i really appreciate it, yes I admit that forgot to change the P(1) to C(3) etc.

I want to ask further about how to make 2 y axis (at the left and right), because the C(1) and C(2) has the same unit in g/L (left y-axis), whereas C(3) has different unit in mikro g/L (right y-axis). Could you help me for the code to make 2 y axis? Thank you so much.

[Star Strider]

My pleasure!

Change the plot call to:

figure

yyaxis left

plot(t,C(:,[1 2]))

ylabel('Concentration (g/L)')

yyaxis right

plot(t, C(:,3))

ylabel('Concentration (\mu g/L')

grid

legend('C_1','C_2','C_3', 'Location','best')

That should do what you want.

Also, I changed the ‘LunaFlores’ code to make it a bit more readable, corrected parentheses pairing in ‘mu’, and changed the initial conditions vector so that they are now correct.

With those changes, it works —

CobaLunaFlores_ode23

function CobaLunaFlores_ode23

Rentang = linspace(0, 80, 50);
C0 = [28.90 0.62 0.07586];
C0 = [C0(3) C0(1:2)];
Miumax=0.145;
Ks=1.6607;
Ki=50.3151;
Kp=8798.15;
gamma=3.6658;
Yxs=0.3941;
alfa=371.45;
betha=1.5502;
ms=0.005;

[t,C]=ode23(@(t,C)LunaFlores(t,C,Miumax,gamma,ms,betha,Kp,Ki,alfa,Ks,Yxs),Rentang,C0);

figure
yyaxis left
plot(t,C(:,[1 2]))
ylabel('Concentration (g/L)')
yyaxis right
plot(t, C(:,3))
ylabel('Concentration (\mu g/L')
grid
legend('C_1','C_2','C_3', 'Location','W')

end

function dCdt=LunaFlores(t,C,Miumax,gamma,ms,betha,Kp,Ki,alfa,Ks,Yxs)
mu = ((Miumax*C(2))/(Ks+C(2)+(C(2)^2/Ki)))*((1-C(3)/Kp)^gamma);
dCdt(1,:) = mu*C(1);
dCdt(2,:) = -C(1)*((mu/Yxs)+ms);
dCdt(3,:) = C(1)*((alfa*mu)+betha);
end

.

[Nadya Shalsabila salman]

Wow! Thank's a lot for your big help to my problem @Star Strider. That's very kind of you. I wish you always health and happiness :)

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Star Strider]

The yyaxis left lower limit of -5 appears to result from ‘C(2)’ becoming slightly negative for some reason. It does not do that when plotted alone, and neither do the data when plotted alone. I fixed that by simply defining the yyaxis left values for ylim to begin at 0. However I was uncomfortable with that, so instead defined it as:

ylim([min(C(:,2)) max(ylim)])

so that range is now [-0.2027 30] and the plot appears to be complete.

There was a problem with the legend that I also fixed, so that it plots the concentrations calculated by the differential equations, and does not include the data. (It is obvious which data points are modeled by each differential equation, so identifying them is likely not necessary anyway.)

It also works out that the numbers of ticks on both y-axes are the same. (That is just luck!)

CobaLunaFlores_ode23

function CobaLunaFlores_ode23
Rentang = linspace(0, 80, 50);
C0 = [28.90 0.617 0.75]; %ga sesuai webplot
C0 = [C0(3) C0(1:2)];
Miumax=0.145;
Ks=1.6607;
Ki=50.3151;
Kp=8798.15;
gamma=3.6658;
Yxs=0.3941;
alfa=371.45;
betha=1.5502;
ms=0.005;
[t,C]=ode23(@(t,C)LunaFlores(t,C,Miumax,gamma,ms,betha,Kp,Ki,alfa,Ks,Yxs),Rentang,C0);
tt=[0
    12
    24
    36
    48
    60
    72];
X = [0.6178489703 28.90160183 75.86206897
    0.823798627 28.83295195 404.5977011
    2.402745995 21.62471396 682.7586207
    5.972540046 13.04347826 1858.62069
    8.306636156 6.178489703 3401.149425
    9.885583524 2.402745995 4488.505747
    9.542334096 2.059496568 5044.827586];

figure
yyaxis left
plot(tt,X(:,[1 2]),'p')
hold on
hlp=plot(t,C(:,[1 2]));
hold off
ylabel('Concentration (g/L)')
ylim([min(C(:,2)) max(ylim)])
yyaxis right
plot(tt,X(:,3),'p')
hold on
hlpp=plot(t,C(:,3));
hold off
ylabel('Concentration (\mug/L)')
grid
legend([hlp([1 2]); hlpp],'C_1','C_2','C_3', 'Location','W')
end
function dCdt=LunaFlores(t,C,Miumax,gamma,ms,betha,Kp,Ki,alfa,Ks,Yxs)
mu = ((Miumax*C(2))/(Ks+C(2)+(C(2)^2/Ki)))*((1-C(3)/Kp)^gamma);
dCdt(1,:) = mu*C(1);
dCdt(2,:) = -C(1)*((mu/Yxs)+ms);
dCdt(3,:) = C(1)*((alfa*mu)+betha);
end

.

[Nadya Shalsabila salman]

Thank you very much @Star Strider I really aprreciate it!

[Star Strider]

As always, my pleasure!

.