esitimate parameter and curve fitting kinetic model

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?

Thanks in advance!

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]

Duplicate to: estimate parameter from kinetic model

Answer 1

[Bjorn Gustavsson]

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]

Thanks a lot for the corrections and suggestions @Bjorn Gustavsson, I'll try that first