Need help plotting modified kinetic equation

Hi, nice to meet you. I am trying to plot graphs (Substrate vs time, Biomass vs time, and Product vs time) from these modified kinetic equations below:

• EDIT (the dX/dt equation should be like this:)

dX/dt = miumax * S/(Ks*S*(S^2/Kxi)) * I/(KXI+I) * (1-X/Xm) *(1-P/Pm)

-dS/dt = 1/Yxs*(dX/dt) + 1/Yps*(dP/dt) + miusx*X

dP/dt = (Ypx*(dX/dt) + miupx*X) * S/(Kps+S+(S^2/Kpi)) * P/(Kp+P) * I/(KPI+I)

where miumax, Ks, Kxi, KXI, Xm, Pm, Yxs, Yps, miusx, Ypx, miupx, Kps, Kpi, Kp, and KPI are kinetic parameters; S, X, P are substrate concentration, biomass concentration, and product concentration respectively.

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

• EDIT

I attached a new MATLAB file in the comment section named "trycurvefit.m", which is a lot more understandable than this code below

function equation1
t = [1:1:20]; %time: 20 days
C0 = [0.1143; 0.9659; 0.004406]; %intial concentration of biomass, susbtrate, and product respectively
[t C] = ode15s(@fungsi,t,C0);
plot(t,C)
end

function dCdt = fungsi(t,C)
Yxs = 2.6872; Yps = 60.3616; miusx = 1.8216; Ypx = 20.7502; miupx = 6.0679; Kps = 0.0017; Kpi = 63.9593; Kp = 19.1078; KPI = 72.5525; %kinetic parameter values
dCdt(1,:) = (Yxs*(-dCdt(2,:) - (1/Yps)*dCdt(3,:) - miusx*C(1))); %dX/dt, i got this changing the dS/dt equation
dCdt(2,:) = (-1/Yxs)*dCdt(1,:)-(1/Yps)*dCdt(3,:)-miusx*C(1);%dS/dt
dCdt(3,:) = (Ypx*dCdt(1,:) + miupx*C(1))*(C(2)/(Kps+C(2)+C(2)^2/Kpi))*(C(3)/(Kp+C(3)))*(45/(KPI+45)); %dP/dt
end 

%i also attach the substrate (S), biomass (X), and product (P) data, so the
%plot sholud be close to these datas
% S = [0.9659 0.6314 0.3819 0.2725 0.2177 0.1934 0.1812 0.1630 0.1326 0.1204 0.1022 0.0778 0.0657 0.0535 0.0413 0.0353 0.0292 0.0170 0.0048 0.0047]
% X = [0.1143 0.1508 0.1873 0.2603 0.5462 0.8138 1.0754 1.2032 1.3126 1.3856 1.6350 1.8600 2.0243 2.1277 2.3224 2.4258 2.5048 2.5353 2.5900 2.6386]
% P = [0.4406 0.7244 0.8622 1.0002 1.1380 1.4217 1.8510 2.5716 3.5832 6.0514 7.5001 8.8027 11.7077 15.1958 18.5373 21.0054 30.3182 36.5723 41.8075 53.5960]

[Star Strider]

I have no idea what result you want.

If your intent is to estimate the parameters, there are several options, one of which is Parameter Estimation for a System of Differential Equations

Start there. If you have problems, post back with a detailed description of them.

[Regina Maria]

@Star Strider thank you so much for answering my question! The link you gave me is definitely a big help and that is actually how I want the result to be. I have tried the code with my data and equation; however, there is 1 curve fit line which is not quite fitting to the data (C2(t)). Here is the plot result of the code:

It will be my pleasure if you could help me out about the C2(t) curve. I have attached the matlab file here. Thank you!!

[James Tursa]

"... there is no dX/dt equation ..."

But you need an independent dX/dt equation. Having three state variables and only two differential equations isn't enough. Consider these lines from your code:

dCdt(1,:) = (Yxs*(-dCdt(2,:) - (1/Yps)*dCdt(3,:) - miusx*C(1))); %dX/dt, i got this changing the dS/dt equation
dCdt(2,:) = (-1/Yxs)*dCdt(1,:)-(1/Yps)*dCdt(3,:)-miusx*C(1);%dS/dt

You can't use the same equation twice as two different differential equations. It looks like you simply used the dS/dt equation to get an expression for dX/dt and then called that your third differential equation. This is not correct. It is not an independent equation for dX/dt. I.e., you can't have dCdt(1,:) depend on dCdt(2,:) and dCdt(3,:) before dCdt(2,:) and dCdt(3,:) are even calculated. This is going to give garbage results.

You need to find that third independent equation for dX/dt for your problem.

[Star Strider]

I am using the genetic algorithm (ga) function to fit the differential equations, since it generally does better than gradient descent optimisers for these sorts of problems.

I note that your code did not constrain the parameters. I constrain kinetics parameters to be positive, since kinetic parameters generally are positive. (If this is not true for your system, please post which are allowed to be negative.)

I am fitting the initial conditions as parameters as well in my atttempt to fit the data. This is turning oout the be a difficult problem, specifically with respect to . Please check to be certain is coded correctly.

I will post the code for it when I get it to work correctly with your system of differential equations.

EDIT — (17 Jun 2021 at 3:08)

I have been working on this most of the day and have not been able to get to fit. (Part of the problem is fitting 15 parameters with 20 data.) Please check to be certain it is coded correctly.

Stopping here.

EDIT — (17 Jun 2021 at 12:34)

Even with the change in I cannot get this to work, either with lsqcurvefit or ga, and I seriously doubt that any optimisation routine could correctly estimate the parameters. It is likely that the model needs to be simplified or changed in some way so that it correctly describes the data. Just now it does not appear to do so.

If you succeed in correctly estimating the parameters, please post back here with them and the plotted results.

Stopping here.

.

[Regina Maria]

@Star Strider Hello once again thank so much you for putting a lot of effort to my problem, I really appreciate that. It is true that all of the kinetic parameters are positive. So sorry to inform you late, but I notice something was wrong with my C2 (dcdt (2)) formula. The whole formula should be multiplied by -1, so it can became like this:

dcdt(2) = -1.*((1/theta(7)).*dcdt(1) + (1/theta(8)).*dcdt(3) + theta(9).*c(1));

or like this

dcdt(2) = -((1/theta(7)).*dcdt(1) - (1/theta(8)).*dcdt(3) theta(9).*c(1));

However, I tried to run the code with this new dcdt(2) code but I just cannot get the result since MATLAB it only shows "Busy" without any result (it takes a really long time to process, I think it is not normal). Do you have any idea or suggestion tho this matter?

[Regina Maria]

@James Tursa Hello, thank you for stopping by and give some advices! Yes, after posting this problem, I realised that my dX/dt equation should not be dependent to dS/dt like that. Therefore, I have revised the dX/dt equation and attached the code file in the previous comment I've posted. Thank you for the info!