# esitimate parameter and curve fitting kinetic model

6 views (last 30 days)
Alfiyya Isfahani on 28 Sep 2022
Commented: Alfiyya Isfahani on 29 Sep 2022
Hi! I want to ask some question.
So, I want to estimate parameter from this kinetic model. The reaction is isomerization glucose to fructose with net rate constant (reversible reaction). The parameter that I want to estimate is 'k' (net rate constant) and 'K' (equilibrium constant) from experimental data.
I make this code with fminsearch and ODE45 for solve the problem, but when I change the initial guess, the result is change too. Is there anyone can help me with the code?
function estimpara
clc
clear all
function dcdt = subpro(t,C,p)
%Unknown Parameters
K = p(1); %equilibrium constant
k = p(2); %net rate constant
%Model
dcdt = zeros(2,1);
dcdt(1) = -k.*(C(1).*(1-(1./K)));
dcdt(2) = k.*(C(1).*(1-(1./K)));
dcdt = dcdt;
end
function obj = objective(p)
%initial concentration at t=0
C0 = [0.068;0.015];
%spesific time points
ts = [1, 2, 3, 5, 8, 10, 12, 15, 20, 30];
%function handle to pass p through substrate function and obtain numerical
%model by solving nonlinear differential equation
[t,C] = ode45(@(t,C)subpro(t,C,p),ts,C0);
%experimental data at each time point ts
Cmeasured = [0.068 0.015
0.065 0.016
0.062 0.017
0.06 0.018
0.059 0.02
0.058 0.0206
0.057 0.021
0.056 0.022
0.0553 0.023
0.055 0.024];
%objective function to minimise the square of the difference
%between the numerical model and experimental data
A = rms((Cmeasured(:)-C(:))./Cmeasured(:));
obj = sum(A);
end
%Parameter Initial Guess
K = 2;
k = 1;
p0 = [K;k];
fun = @objective;
options = optimset('Display','iter','PlotFcns',@optimplotfval);
[p,fval] = fminsearch(fun,p0,options);
%optimized parameter values (untuk katalis Amine)
K = p(1);
disp(['K : ' num2str(K)])
k = p(2);
disp(['k : ' num2str(k)])
%calculate model with updated parameter
p1 = K;
p2 = k;
ts = [1, 2, 3, 5, 8, 10, 12, 15, 20, 30];
C0 = [0.068;0.015];
[t,C] = ode45(@(t,C)subpro(t,C,p),ts,C0);
Cmeasured = [0.068 0.015
0.065 0.016
0.062 0.017
0.06 0.018
0.059 0.02
0.058 0.0206
0.057 0.021
0.056 0.022
0.0553 0.023
0.055 0.024];
%plot of updated parameter values
figure (1)
plot(ts,Cmeasured(:,1),'o',ts,C(:,1),'b')
hold on
plot(ts,Cmeasured(:,2),'o',ts,C(:,2),'b')
xlabel('time')
ylabel('concentration')
end
Star Strider on 28 Sep 2022

Bjorn Gustavsson on 28 Sep 2022
Is your model correct? To my (very much non-biology-expert) eyes the ODE for c(1) seem to reduce to:
dcdt(1) = C(1)*k*(1/K - 1);
That can be re-parameterized with only one reaction-konstant, lets say B:
B = k*(1/K-1);
Fminsearch can find you best values for B but any pair of k and K that produce that single-parameter B will give the same fit to your data. That's why you get different values from different starting-points.
Since this model leads to an exponential decrease of c(1) it doesn't make much sense to talk about equilibrium and reversible reaction - unless there should also be some transition from c(2) to c(1) as well.
HTH
Alfiyya Isfahani on 29 Sep 2022
Thanks a lot for the corrections and suggestions @Bjorn Gustavsson, I'll try that first

### Categories

Find more on Statistics and Machine Learning Toolbox in Help Center and File Exchange

R2016a

### Community Treasure Hunt

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

Start Hunting!