Predicting the kinetic constant of a reaction based on experimental data

I have the following optimised reaction based on experiments: A + 1.7B -> C (A = Rifamycin Oxazine; B = Piperazine; C = Rifampicin)

My postulated rate law is therefore:

r_A = -k*C_A*C_B^1.7

r _B= -k*C_A*C_B^1.7

r _C= k*C_A*C_B^1.7

Using my experimental data, I have used the following code to try and accurately find k (the goal is to get an R^2 > 90%):

function API2

function C=kinetics(theta,t)
c0=[0.752;1.278;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(3,1);
    dcdt(1)=-theta(1).*(c(1).^(1)).*(c(2).^(1.7));
    dcdt(2)=-theta(1).*(c(1).^(1)).*(c(2).^(1.7));
    dcdt(3)=theta(1).*(c(1).^(1)).*(c(2).^(1.7));
    dC=dcdt;
    end
C=Cv;
end

T = [0 1 2 5 10 15];
t = T';

%y values for a
a_ydata = [0.752 0.0596 0.0596 0.0596 0.0502 0.0424];
A_Ydata = a_ydata';

%y values for b
b_ydata = [1.278 0.378 0.101 0.101 0.085 0.072];
B_Ydata = b_ydata';

c_ydata = [0 0.692 0.692 0.692 0.702 0.71];
C_Ydata = c_ydata';

c = [A_Ydata B_Ydata C_Ydata];
theta0=[0.5];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
h = plot(t, c,'.');
set(h,{'Marker'},{'s';'d';'^'},{'MarkerFaceColor'},{'r';'b';'k'},{'MarkerEdgeColor'},{'r';'b';'k'});
hold on
hlp = plot(tv, Cfit,'LineWidth',1.5);
set(hlp,{'Color'},{'r';'b';'k'});
hold off
grid
xlabel('Time (min)')
ylabel('Concentration (M)')
legend(hlp, 'Rifamycin Oxazine', 'Piperazine', 'Rifampicin', 'Location','N')

end

However, when I run the code and produce a graph, there seems to be a clear discrepancy between my simulated curve and my experimental data (especially for piperazine). I have been messing around with the code but can't seem to get it to fit well with my data. I'm not sure what to do. Can anyone help me?

Answer 1

[Star Strider]

If you have the Global Optimization Toolbox, the code I attached will estimate the parameters with reasonable accuracy, although your model extimates only one parameter in three differential equations. (This is not something I have observde previusly, so it may be appropriate for you to reconsider the model.)

That aside, ga managed to converge, although with somewhat different estimates in each run.

Note — I set the initial conditions as parameters to be estimated in my code.

Typical parameters were:

	Rate Constants:
		Theta(1) =  2.55140
		Theta(2) =  0.92250
		Theta(3) =  0.99550
		Theta(4) =  0.01400

and

	Rate Constants:
		Theta(1) = 10.25306
		Theta(2) =  0.88801
		Theta(3) =  0.98375
		Theta(4) =  0.00308

with essentially identical fitness values, producing this plot:

with the second set of parameters.

That is likely as good as it gets for you model.

.

[BobbyJoe]

Thank you so much for your help star strider, I greatly appreciate this. It seems that it's possible my model could be incorrect, however, I'm not sure what other models would give me a more accurate fit. Do you have any suggestions of other differential rate equations I should try which would give me a more accurate fit?

[Star Strider]

My pleasure!

I have no idea what system you are modeling.

It would be straightforward for you to adapt different models to my code. (I am certainly willing to help you with that if necessary.)

[Image Analyst]

Like I said, you're going to get nothing good and useful out of that data. You need a lot more to get any sort of accuracy on the rate. Garbage in = garbage out.

[Star Strider]

Bobby Jo —

If my Answer helped you solve your problem, please Accept it!

.

Answer 2

[Image Analyst]

Not really sure I follow your code (like what is the function API2, and values theta and t), but for what it's worth, I'm attaching a program that fits data to the rate equation.

[BobbyJoe]

T is time, theta is the kinetic constant k. The function API2 is just simply for the whole document. I can't open the .csv file.

I'm not very proficient in MATLAB. If you were to use my experimental data results (A, B, C (these are the concentration results) and T) along with the rate equation, do you mind processing this on your matlab code and see if it works?

It would mean so much to me if that's possible.

[Image Analyst]

From the looks of the data in Star's answer you don't have nearly enough data points in the steep part of the curve to estimate a rate. Just can't do it. You just have an intial value and the final steady state value. You'd need way more data to determine accurately where the "knee" is and the rate equation coefficients. Please supply more data. Otherwise you know what they say : "garbage in garbage out".