# how to model kinetic reaction system

22 views (last 30 days)
Anastasya Silitonga on 24 May 2020
Edited: darova on 26 May 2020
Research area: Chemical kinetics
problem description: So this is parallel kinetics. I want to determine rate of reaction (as dCdt) and constant rate (as k) from experimental data that I have, time (as t) and concentration (as C). There are 6 variable called as mag,dag,tag,gly,ffa, and air.
so there are 6 constant rate and rate of reaction that I need to know, using ode23 or ode45. and the equations are described below:
dCdt(mag)=k1*Ctag*Cgly+2*k2*Cdag*Cgly-k3*Ctag*Cmag+k5*Cdag*Cair+k6*Cgly*Cffa;
dCdt(dag)=k1*Ctag*Cgly-k2*Cdag*Cgly+2*k3*Ctag*Cmag+k4*Ctag*Cair-k5*Cdag*Cair;
dCdt(tag)=-k1*Ctag*Cgly-k3*Ctag*Cmag-k4*Ctag*Cair;
dCdt(ffa)= k4*Ctag*Cair+k5*Cdag*Cair-k6*Cgly*Cffa;
dCdt(gly)=-k1*Ctag*Cgly-k2*Cdag*Cgly-k6*Cgly*Cffa;
dCdt(air)=-k4*Ctag*Cair-k5*Cdag*Cair+k6*Cgly*Cffa;
all I want to do is fit my data to this function to solve for k1,k2,k3,k4,k5,k6 and find each dCdt. I've tried make a code by using fminsearch and ode23, but I dont know how to define each concentration for ode so as finding each k. Can anybody help me with the code?
function constantrate
clc
Tspan = [0.1 0.1]
constant=fminsearch(@rate, Tspan);
function constantrate = rate (k)
Data = [
60 90.5 5.2 1.1 10.5 2 3
120 90.6 5.1 1 10.1 1.8 3.3
180 90.7 5.1 1.4 12.5 1.9 4
240 90.7 5.4 2.1 10 1.7 4.5
];
t=Data(:,1); %time
Cmag=Data(:,2); %concentration of mag
Cdag=Data(:,3); %concentration of dag
Ctag=Data(:,4); %concentration of tag
Cgly=Data(:,5); %concentration of gly
Cair=Data(:,6); %concentration of air
Cffa=Data(:,7); %concentration of ffa
[tint, Cint]=ode23(@reaction,t,Cmag,[],k);
constantrate=sum((Cmag-Cint).^2);
plot (tint, Cint, 'r','o',t,Cmag,'b')
legend('calculation','experimental')
xlabel('time (min)')
ylabel('concentration (M)')
function dCdt=reaction(Ctag,Cmag,Cdag,Cgly,Cair,Cffa,k1,k2,k3,k4,k5,k6)
dCdt(mag)=k1*Ctag*Cgly+2*k2*Cdag*Cgly-k3*Ctag*Cmag+k5*Cdag*Cair+k6*Cgly*Cffa;
dCdt(dag)=k1*Ctag*Cgly-k2*Cdag*Cgly+2*k3*Ctag*Cmag+k4*Ctag*Cair-k5*Cdag*Cair;
dCdt(tag)=-k1*Ctag*Cgly-k3*Ctag*Cmag-k4*Ctag*Cair;
dCdt(ffa)= k4*Ctag*Cair+k5*Cdag*Cair-k6*Cgly*Cffa;
dCdt(gly)=-k1*Ctag*Cgly-k2*Cdag*Cgly-k6*Cgly*Cffa;
dCdt(air)=-k4*Ctag*Cair-k5*Cdag*Cair+k6*Cgly*Cffa;

Star Strider on 24 May 2020
There are several examples on fitting differential equations to data in Answers. See: ODE and Data fitting and Parameter Estimation for a System of Differential Equations. There are many others.

Star Strider on 25 May 2020
You need to have the Global Optimization Toolbox to use the ga function. I am running 2020a, and the code I posted (for reference) ran without error.
My objective otherwise was to present a set ot parameters (and slightly revised code) that you can use as initial estimates with fminsearch or other optimization functions to estimate the parameters. Those parameters should give a reasonable fit.
I still recommend that you be certain that the order and coding of the equations in ‘kinetics’ are correct. I generally expect a much better fit.
Anastasya Silitonga on 26 May 2020
Thanks again! this is new for me and I learn a lot! Do you mind if I ask again? With these equations, is it possible to find the best order and coefficient of rate at the same time? I understand that it could happen by using fminsearch and ode23 for single equation, but I don't understand since it contains so much parallel equations.
Star Strider on 26 May 2020
My pleasure!
This finds the best coefficients of the differential equation that minimise the norm or least-square difference between the differential equation and the data. I do not know the system you are modeling (it appears to be lipolysis or something similar), so I do not know what the parameters mean. As I remember (from undergraduate physical chemistry back when log tables and slide rules were all the rage, and computers were behemoths that ate punchcards), the order would be given by the exponents of the reactions in kinetic equations, and I do not see any exponents here (so I assume everything is first-order in the respective reactants).
So if you want to optimise or estimate the order and rate, and the differential equations are set up to do that, optimising the fit will estimate the parameters you want.