Plotting solution of ODE System

Wondering why I get such a different graph for the first variable of the system I've attached, together with the values of the initial conditions and of the parameter of the model, and with what I get when I plot the solution and what the authors of the paper got.

I've first implemented this function

    function dydt = syst2( t, y )
dydt = zeros (5, 1);
a = 1.2763;
b = 1.2423;
c = 1.2763; 
d = 1.1529;
e = 1.355;
A = 5.622;
B = 93.023;
C = 0.001;
D = 0.0018; 
E = 0.000101;
F = 0.01;
G = 13.2;
H = 1440;
I = 1.8;
mi = 0.0667;
Ks = 2.15;
beta = 0.0045;
CO2295 = 92;
Yxs = 0.04;
Yco2e = 0.775;
dydt(2) = (1/Yxs)*(-y(1)*(mi * y(2)/(y(2)+Ks+B^b)-E*y(1)))-F*y(1);
dydt(1) = ((exp(-(y(3)-CO2295)))/(exp(y(3)-CO2295)+ exp(-(y(3)-CO2295))))*A*mi*(y(2)/(y(2)+Ks*B^a))*y(1)*(1-(y(1)/(A*mi*(y(2)/(beta*(y(2)+Ks*B^a))))))+(1-((exp(-(y(3)-CO2295)))/(exp(y(3)-CO2295)+ exp(-(y(3)-CO2295)))))*(C*y(1)*dydt(2)-D*y(1));
dydt(5) = H*mi*(y(2)^2 /((y(2)+Ks*B^d)*(y(2)+Ks*B^e)))*y(1)+I*y(1);
dydt(3) = G*mi*(y(2)/(y(2)+Ks*B^c))*y(1)+ dydt(5);
dydt(4) = (1/Yco2e)* dydt(3);
  end

and then this script:

    tspan = [0 500];
yzero= [ 0.215 208.5 0 0 1.6081];
[T, Y] = ode45 (@syst2, tspan, yzero);
plot (T,Y(:,1),'m');
legend ('X');

What do I do wrong?

%

Answer 1

[Star Strider]

It’s difficult for me to follow your code, because I’m not certain what the variables are. (It’s always a good idea to use comments to document your code so you know what everything represents, and how you adapted the differential equations in the paper.)

You have five differential equations, one of which, dydt(5), is likely the second term in the dCO2/dt equation, while the information you posted has four. I haven’t coded biochemical differential equations in decades, so I’m not certain how I would code dCO2/dt. How did these or other authors code dCO2/dt? How did they deal with the derivative term?

Also, check the paper to see the solver the authors used. Different solvers give different results. The MATLAB solvers are likely the best available, so unless the authors used them, you could get different results.

[Viviana Arrigoni]

As you can see, the system has really many parameters, which I am initializing in the function "syst2" with the values given in the table I attached. That's what all the a,b,c ,..., Yxs, Yco2e are. The system's got 4 equation, but looking at it, you may notice that dCO2/dt has a derivative inside, which is not explicity a function of t, but rather of S and X. That's why I have defined dydt(5): that's the second addend of dCO2/dt. Then dydt(1)=dX/dt, dydt(2)=dS/dt, dydt(3)=dCO2/dt, dydt(4)=dP/dt. Yet how you can see, I am not plotting dydt(5) at all (just giving it a coherent initial condition). I have sorted them differently in the function because in dX/dt, the term dS/dt appears, so I decided to define it first.

I have no idea of how the authors dealt with the second term of dCO2/dt, as they did not give the code of their simulations, which the made using Matlab.

[Star Strider]

If they used MATLAB, the last part of my original Answer is irrelevant.

I would email them and ask how they coded dCO2/dt. (Alternatively, one of their references might have it.) Please post their reply here, because I would like to know!

[Kevin Doherty]

As Star Strider said, the problem is probably with how the derivative terms are defined. At a particular time-step ode45 uses values from the previous time-step. But here you are using dydt values from the current time-step.

[Star Strider]

@Kevin — Thank you, I remain mystified just now. I thought about using persistent variables (and mlock) to track the previous values of the components of the derivative term in dCO2/dt and then numerically calculate the derivative from a vector-valued tspan with a known step, but I’m not certain that would work. (I tend not to use persistent variables, so I have no experience with them.)

@Viviana — Please post the PDF of the paper. It may have a clue or references to solving similar problems. I do not recognise the system you are simulating (other than it appears to be involving reactants, enzymes, and products, but then that is not a unique situation).