65 views (last 30 days)

Show older comments

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.

Star Strider
on 24 Feb 2021

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
on 26 Feb 2021

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.)

Alan Stevens
on 24 Feb 2021

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!