esitimate parameter and curve fitting kinetic model

4 views (last 30 days)
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

Answers (1)

Bjorn Gustavsson
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
  1 Comment
Alfiyya Isfahani
Alfiyya Isfahani on 29 Sep 2022
Thanks a lot for the corrections and suggestions @Bjorn Gustavsson, I'll try that first

Sign in to comment.

Categories

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

Community Treasure Hunt

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

Start Hunting!

Translated by