Fitting Kinetic model to estimate kinetic parameter.

Dear matlab users.

I have the following problem of fitting this kinetics ecuations.

where i have these experimental data.

Time=[0 30 60 90 150 210 270 330];

c= [0.15780 0.12364 0.10460 0.05018 0.01702 0.00539 0.00146 0.00022

0.3111 0.2676 0.2490 0.1807 0.1214 0.1029 0.1079 0.0968

0.0000 0.0007 0.0000 0.0000 0.0002 0.0000 0.0002 0.0000

0.0000 0.0225 0.0591 0.1135 0.1368 0.1519 0.1773 0.1629

0.000000 0.000110 0.000390 0.000433 0.000400 0.000193 0.000262 0.000215

0.00000 0.00000 0.00482 0.00000 0.00000 0.00000 0.00000 0.00000

0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050];

So.

When i run the code i got something like this

but what i need its a better fit, so my question its how can I improve the fitting :c ?

Thanks in advise!!!.

Answer 1

[Torsten]

You treat the initial conditions at t=0 as unknowns k(8:14).

But you have the concentrations at t=0 from your experimental data and you include the difference between your experimental conditions at t=0 and the parameters in your function to be minimized.

Can you explain this contradiction ?

[Catalina Valdes Bustamante]

srry, you mean by treat k(8:14) like the initial concentrations.

when i do that i got this :c

with the initial conditions I got the same results. maybe I'm not doing it in the right place.

Answer 2

[Star Strider]

I was able to get a reasonably decent fit with the ga function:

Elapsed Time: 3.281440000000000E+02 00:05:28.144

Fitness value at convergence = 0.1428

Generations = 3144

Rate Constants:

Theta( 1) = 0.56724

Theta( 2) = 0.55739

Theta( 3) = 0.08265

Theta( 4) = 2.61614

Theta( 5) = 10.33273

Theta( 6) = 13.16327

Theta( 7) = 12.25866

Theta( 8) = 0.08987

Theta( 9) = 0.31539

Theta(10) = 0.00679

Theta(11) = 0.05491

Theta(12) = 0.00001

Theta(13) = 0.00208

Theta(14) = 0.09136

and:

Elapsed Time: 1.045420000000000E+02 00:01:44.542

Fitness value at convergence = 0.1467

Generations = 1099

Rate Constants:

Theta( 1) = 0.13283

Theta( 2) = 0.90345

Theta( 3) = 0.07339

Theta( 4) = 8.84979

Theta( 5) = 7.17460

Theta( 6) = 10.56757

Theta( 7) = 10.11484

Theta( 8) = 0.07805

Theta( 9) = 0.30280

Theta(10) = 0.00402

Theta(11) = 0.04414

Theta(12) = 0.00011

Theta(13) = 0.00162

Theta(14) = 0.09064

I usually expect better results than this, so be certain that the differential equations are coded correctly. Also, use the ode15s function, since the differential equations are ‘stiff’ (the parameters vary over a wide range).

.

[Catalina Valdes Bustamante]

Hi Star, could you get the source code to fully understand it :c I would love to take a look at it uwu

[Star Strider]

I deleted the file however I still had the code up and have not yet deleted it so I’m posting it here:

Time=[0 30 60 90 150 210 270 330];

c= [0.15780 0.12364 0.10460 0.05018 0.01702 0.00539 0.00146 0.00022

0.3111 0.2676 0.2490 0.1807 0.1214 0.1029 0.1079 0.0968

0.0000 0.0007 0.0000 0.0000 0.0002 0.0000 0.0002 0.0000

0.0000 0.0225 0.0591 0.1135 0.1368 0.1519 0.1773 0.1629

0.000000 0.000110 0.000390 0.000433 0.000400 0.000193 0.000262 0.000215

0.00000 0.00000 0.00482 0.00000 0.00000 0.00000 0.00000 0.00000

0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050];

t = Time;

ftns = @(theta) norm(c-kinetics(theta,t));

PopSz = 25;

Parms = 14;

% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-4, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',[randi(1E+4,PopSz,Parms/2)*(1E-3) repmat(c(:,1).',PopSz,1)], 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

% optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1); % With 'HybridFcn'

t0 = clock;

fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts);

t1 = clock;

fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

GA_Time = etime(t1,t0)

DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');

fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

fprintf(1,'\tRate Constants:\n')

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure

hd = plot(t, c, 'p');

for k1 = 1:size(c,1)

CV(k1,:) = hd(k1).Color;

hd(k1).MarkerFaceColor = CV(k1,:);

end

hold on

hlp = plot(tv, Cfit);

for k1 = 1:size(c,1)

hlp(k1).Color = CV(k1,:);

end

hold off

grid

xlabel('Time')

ylabel('Concentration')

legend(hlp, compose('C_%d',1:size(c,1)), 'Location','N')

function c = kinetics(k,t)

c0 = k(8:14);

[~,Cv] = ode15s(@DifEq,t,c0);

function dcdt=DifEq(t,c)

Cfenol=0.311083115088667;

dcdt = zeros(7,1);

r1=(k(1)*c(1)/((1+c(7)+k(6)*c(1)+k(7)*c(2))^5));

r2=(k(2)*c(6)*c(2)/((1+c(7)+k(6)*c(1)+k(7)*c(2))^2));

r3=(k(3)*c(2)^2/((1+c(7)+k(6)*c(1)+k(7)*c(2))^2));

r4=(k(4)*c(3)/((1+c(7)+k(6)*c(1)+k(7)*c(2))^3));

r5=(k(5)*c(3)/((1+c(7)+k(6)*c(1)+k(7)*c(2))));

dcdt(1)=(1/Cfenol)*-r1;

dcdt(2)=(1/Cfenol)*(-r2-(2*r3));

dcdt(3)=(1/Cfenol)*(r2-r4-r5);

dcdt(4)=(1/Cfenol)*(r3+r4);

dcdt(5)=(1/Cfenol)*r5;

dcdt(6)=(1/Cfenol)*(r1-r2);

dcdt(7)=(1/Cfenol)*((-2*r1)-r4+(2*r5));

end

c = Cv.';

end

It takes several attempts to get it to run continuously, since the initial parameters do not always work for all the individuals in the population. I expect a much better fit that was provided here, so be sure the differential equations are coded correctly. (I have no way of checking that.)

.