How to calculate the rate using a given data set?

Hey guys,

I need to estimate the reaction rate from given data in excel. To answer this I set up the differential equations and appointed the given data (for the sake of convenience I briefly reduced the data set here). My first thought was that the most simple way to estimate the rate was to use ODE45, but since it doesn't works I assume I made a mistake somewhere.

This was the code I made:

function [dxdt] = data(t)
A = [0.98 0.97 0.95 0.92 0.90]';
B = [0.99 0.97 0.95 0.93 0.90]';
C = [0.01 0.03 0.07 0.11 0.16]';
A_rad = [0.01 0.02 0.02 0.03 0.03]';
B_rad = [0.01 0.02 0.02 0.03 0.03]';

dAdt = -1 .* A - 100 .* B_rad .* A + 1 .* (A_rad).^2;
dBdt = -100 .* A_rad .* B + 1000 .* B_rad .* C;
dCdt = (20000 ((A).^0.5) .* B) ./(100 + 1000 .*(C./A)); %This is the rate I want to estimate

%dBr_radical = 0; Since equal to zero it can be left out
%dH_radical = 0; 
[dxdt] = [dAdt; dBdt; dCdt]

And next to call it:

init = [1 1 0 0 0];
tspan = [0 100];

[t,c] = ode45(@data, tspan, init, []);
plot(t,c)

Can someone help me improving my code?

Thanks in advance,

Kind regards,

Danny

[darova]

You don't need to use ode45 since you already have table data

This is a result you want (as you wrote above):

dCdt = (20000 ((A).^0.5) .* B) ./(100 + 1000 .*(C./A)); %This is the rate I want to estimate

Answer 1

[Star Strider]

You are not coding your differential equations and data correctly. See for example: Parameter Estimation for a System of Differential Equations.

[Danny Helwegen]

Hey Star Strider,

first of all thanks for you answer. I have looked at the link you gave and understand how to estimate a variable in the differential, but I have no idea how to estimate the outcome of the differential. Can you help me with setting up the right code?

Thanks in advance,

Danny

[Star Strider]

My pleasure.

I do not understand what you mean by ‘estimate the outcome of the differential.’ Note that the file I attached to that Answer has everything. You need to code your own version of the ‘kinetics’ function, include your data, and run the file.

I cannot make any sense out of the code you posted, so I cannot help you write an appropriate version of your differential equations and data to use with that Answer.

[Danny Helwegen]

Hey, I made a reading error hence I asked for the wrong outcome. You're right, the attached file answers almost anything. Only one question, in the attached file the fit is also dependend of a c with different values, only for me the c's have only one constant value. I tried rewriting the code, but the error I get is:

LSQCURVEFIT requires four input arguments.

So if i take the c out of the formula, what should be the fourth argument?

This is the code I tried:

function Solver
% 

function C=kinetics(theta,t)
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t)
    dcdt=zeros(4,1);
    dcdt(1)=-theta(1).*1-theta(2).*1;
    dcdt(2)= theta(1).*1+theta(4).*10-theta(3).*10-theta(5).*10;
    dcdt(3)= theta(2).*1+theta(3).*10-theta(4).*10+theta(6).*1;
    dcdt(4)= theta(5).*10-theta(6).*1;
    dC=dcdt;
    end
C=Cv;
end



t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];


theta0=[1;1;1;1;1;1];

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

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);


end

[Star Strider]

I am having a difficult time understand iing your version of the ‘kinetics’ function.

First, in this code, ‘theta’. is the parameter vector. The ‘c’ vector is the vector of the dependent variables for the differential equation. (That code models compartmental kinetics, so the ‘c’ vector is the time-varying concentration of the substance in various compartments. The differential equations describe the change in those concentrations in each compartment over time.)

It would be helpful if you could provide a symbolic description of your system of differential equations, preferably as a .pdf file. I can probably be of more help if I understand what you want to do.

[Danny Helwegen]

I will try to explain it now as clear as possible:

I have a reaction mechanism and from here I set up a number of differential equations that describe the rate of the diffrent species ( eg. A --> A_rad)

Next, I need to use this set of equations to simulate the concentrations of the species and compare them to expirimental results given in a data set.

My differential equations are:

dAdt = (20000 * (theta(1)).^0.5 .* theta(2))./(100+10.*(thetha(3))./(thetha(1)))
dBdt = sqrt(0.1.*theta(1))
dCdt = (10*sqrt(theta(1)).*theta(2))/(10.*theta(1) + 10.*theta(3))

The time goes from 0 to 10 minutes in steps of 0.2 seconds.

This is pretty much what I need to do, I hope this is enough information to make the challenge clear.

[Star Strider]

Something is not correct.

Your differential equations each need to be functions of at least one of the variables (‘A’, ‘B’, ‘C’) and your parameter vector ‘theta’. They are not, so you cannot integrate them, and if you cannot integrate them, you cannot fit them to your data and estimate the parameters.

[Danny Helwegen]

I think that I get it, I have to few data to simulate?

Perhaps if I give the differentials exactly, maybe you know another way to simulate?

Real differentials:

dAdt = (20000 * (B).^0.5 .* C)./(100+10.*(A)./(B))
dBdt = sqrt(0.1.*B)
dCdt = (10*sqrt(B).*C)/(10.*B + 10.*A)

(For convenience I had already replaced the A's, B's and C's in the comment above)

[Star Strider]

That looks tractable.

You need to decide on a convention for the variable representation in a vector, so for example:

v(1) = A
v(2) = B
v(3) = C

your differential equations are then coded as:

dvdt(1,:) = (20000 * (v(2)).^0.5 .* v(3))./(100+10.*(v(1))./(v(2)))
dvdt(2,:) = sqrt(0.1.*v(2))
dvdt(3,:) = (10*sqrt(v(2)).*v(3))/(10.*v(2) + 10.*v(1))

Your constants apparently can vary over several orders of magnitude, so ode15s might be the best choice for a solver.

Where in those equations are the references to the parameters you want to fit? I can do my best to fit your differential equations to your data if I understand their structure and the parameters you want to estimate.

Also, how do ‘A_rad’ and ‘B_rad’ relate to your equations, including the parameters to be estimated?

[Danny Helwegen]

Where in those equations are the references to the parameters you want to fit?

The only reference is the data set I obtained, which are the concentrations of A, B, C, A_rad and B_rad

Also, how do ‘A_rad’ and ‘B_rad’ relate to your equations, including the parameters to be estimated?

I have given you the rewritten differentials were those 2 dropped out, when not rewritten this would be the set of equations:

dAdt = A^2 - A - 100 * A * B_rad
dBdt = 1000 * B_rad * C - 100 * A_rad * B
dCdt = 100 * A_rad * B + 100 * B_rad * A - 1000 * B_rad * C
dA_raddt = dB_raddt = 0 

[Star Strider]

I am lost at this point.

Are you estimating any parameters? (If so, what are they?)

Do you just want to run those differential equations with the values you posted? (If so, what result do you want?)

I am still not clear on what you want to do, and what ‘estimate the reaction rate’ means here.

[Danny Helwegen]

Okay, so I have these differentials and I have to estimate the outcome for A, B, C, A_rad and B_rad and these outcomes have to be compared to the dataset, so they have to be approximatly similar.

Given are the differentials:

dAdt = A^2 - A - 100 * A * B_rad
dBdt = 1000 * B_rad * C - 100 * A_rad * B
dCdt = 100 * A_rad * B + 100 * B_rad * A - 1000 * B_rad * C
dA_raddt = dB_raddt = 0 

A maximum time of 10, the begin values of A, B, C, A_rad and B_rad (=[1 1 0 0 0]) and the data set for comparison.

[Danny Helwegen]

I have found the anser, it was easier than I thought:

This was my final code:

function [dxdt] = Code(t, v)

A = v(1); 
B = v(2); 
C = v(3); 
D = v(4); %B_rad
E = v(5); %H_rad

dAdt = - A - 100 * E * A +  D^2 ;
dBdt = - 100 * D * B + 100 * E * C; 
dCdt = 100 * D * B + 100 * E * A - 1000 * E * C);
dDdt = 2 * A - 100 * D * B + 100 * E * A + E * C - 2000 * D^2 ; 
dEdt = 100 * D * B - 100 * E * A -  E * C; 

dxdt = [dAdt; dBdt; dCdt; dDdt; dEdt;];

And for calling:

init = [1 1 0 0 0 ];
tspan = [0:0.01:10];

[t,v] = ode45(@AJA, tspan, init, []);
plot(t,v)

Thanks for your help star strider.

[Star Strider]

As always, my pleasure!