Coupled rate ODEs with ode45

Dear all,

I have been having trouble trying to model the concentrations of 8 species in the attached file showing the reactions. I have the rate constant values, initial conditions and a time frame but for some reason the plot seems to display no changes in concentration. I have attached my function and script. Any pointers or tips would be greatly appreciated.

Thanks in advance.

Answer 1

[Star Strider]

The concentrations change appropriately, however they don’t change much and the concentrations are vanishingly small. That’s the reason they don’ appear to change.

Increasing the final time (by 15 times) demonstrates saturation kinetics:

[t,state]=ode15s(@reactions,[0 51.5*15],[2.71E-5;7.8E-6;9.2E-6;0;1.0614E-4;3.359E-3;3.9278E-4;4.521E-4]);

reactants = {'TAG','DAG','MAG','GLY','FAEE','ET','FFA','W'};

figure
for k = 1:8
    subplot(4,2,k)
    plot(t,state(:,k))
    xlabel('Time')
    ylabel('Concentration')
    title(sprintf('$C_{%s}(t)$',reactants{k}), 'Interpreter','latex')
    grid
end

producing:

Adding:

pos = get(gcf,'Position')
set(gcf, 'Position',pos+[0 -500 200 500])

after the loop enlarges the figure making the subplot plots easier to see. (I did not use that in the posted figure.)

[Star Strider]

I was in the process of responding when my desktop crashed again (third time today and the tenth time this week so far), wiping out everything I was working on. I so hate Windows10 at this point that no words will describe it. I’ll do my best to re-create what I was writing.

—————————————————————————————————————————————————

As always, my pleasure!

I have (an ancient) background in Chemistry, however not in Chemical Engineering. I’m not certain that it’s possible to calculate kinetic constants from equilibrium constants, since it’s not obvious to me how they’re related. I’ve posted several Answers on estimating rate constants from data (one such is Parameter Estimation for a System of Differential Equations that I’ve subsequently updated in other Answers), so if you have the data that created the plots you posted, estimating the kinetic constants should be straightforward.

If you have information on calculating kinetic constants from equilibrium constants, I’ll look at it and if possible suggest some example code to do the necessary calculations. That would be interesting for me as well. The only information I can find indicates that the equilibrium constants are derived from dividing the forward rate constant by the reverse rate constant, so for any specific reaction it’s likely possible to calculate the equilibrium constant from the rate constants, although not the other way round, unless at least one rate constant is known.

[Star Strider]

‘From what I can see they will be the only unknown values in the problem and it should be entirely possible to determine them from the 4 equations.’

Unfortunately, not. It seems to me that the ‘r’ values are important, and they’re not provided (or could be calculated from the existing data). With those and at least one of the concentration values and the ‘Keq’ values, it should be possible to calculate the ‘k’ values. In the absence of the ‘r’ values, I doubt that this is solvable.

Again, this is not my area of expertise and I’m drawing on what I remember from physical chemistry (back in the phlogiston era), so I could be missing something.

[Star Strider]

As always, my pleasure!

The problem with the latest approach is at least in part the problem of estimating four parameters with only two data points (‘Inlet’ and ‘Outlet’). If you have the time evolution of the various reactants and products over at least one more time instants than the number of parameters to be estimated (more is better), you can use the approach I used in the Answer I linked to. (I’ve updated that code since, so if you want to use it and if you post the necessary data here, as well as the appropriate differential equations, I’ll post back the updated code with an attempt at estimating the parameters, since lsqcurvefit may not be optimal and a genetic algorithm approach may be necessary.)

Answer 2

[Alan Stevens]

It can all be done in one script as follows. Because of the orders of magntude difference between the various concentrations they just look constant when plotted on a single graph. They do vary as can be seen by running the following:

tspan = [0 51.5];
Y0 = [2.71E-5;7.8E-6;9.2E-6;0;1.0614E-4;3.359E-3;3.9278E-4;4.521E-4];
[t,Y]=ode15s(@reactions,tspan,Y0);

TG = Y(:,1);   FAEE = Y(:,5);
DG = Y(:,2);   ET   = Y(:,6);
MF = Y(:,3);   FFA  = Y(:,7);
G  = Y(:,4);   W    = Y(:,8);

subplot(4,2,1)
plot (t,TG),grid
ylabel('TG')
subplot(4,2,3)
plot (t,DG),grid
ylabel('DG')
subplot(4,2,5)
plot (t,MF),grid
ylabel('MF')
subplot(4,2,7)
plot (t,G),grid
xlabel('t'),ylabel('G')
subplot(4,2,2)
plot (t,FAEE),grid
ylabel('FAEE')
subplot(4,2,4)
plot (t,ET),grid
ylabel('ET')
subplot(4,2,6)
plot (t,FFA),grid
ylabel('FFA')
subplot(4,2,8)
plot (t,W),grid
xlabel('t'),ylabel('W')

function dydt = reactions(~,y)
% y1=TG y5=FAEE
% y2=DG y6=ET
% y3=MF y7=FFA
% y4=G  y8=W

k1=0.6;
keq1=6.504E-6;
k2=0.8403;
keq2=0.02386;
k3=0.83524;
keq3=2.464;
k4=0.381;
keq4=12.182;

r1 = k1*((y(1)*y(6)*y(7))-((1/keq1)*y(2)*y(5)*y(7)));
r2 = k2*((y(2)*y(6)*y(7))-((1/keq2)*y(3)*y(5)*y(7)));
r3 = k3*((y(3)*y(6)*y(7))-((1/keq3)*y(4)*y(5)*y(7)));
r4 = k4*((y(7)*y(6)*y(7))-((1/keq4)*y(8)*y(5)*y(7)));

dy1dt=-r1;
dy2dt=r1-r2;
dy3dt=r2-r3;
dy4dt=r3;
dy5dt=r1+r2+r3+r4;
dy6dt=-(r1+r2+r3+r4);
dy7dt=-r4;
dy8dt=r4;

dydt=[dy1dt;dy2dt;dy3dt;dy4dt;dy5dt;dy6dt;dy7dt;dy8dt];

end

You can save and run this as a single script file.