How to properly model a kinetic reaction system from experimental data using lsqcurvefit?

Hello,

I am currently trying to model chemical kinetics using MATLAB.

Currently, I am having a code that works but gives me a set on non-sense results. So there must be a mistake in the method, or in my way of thinking.

I will introduce the problem, the available data and my MATLAB integration. I hope you can help me find the mistakes and obtain a fit for my experimental data.

Thanks in advance.

Problem introduction:

This system consists of four reactions. Compound c1 reacts to compound c2. Compound c2 can react further to c3, although I am currently neglecting the reaction of c2 to c3. The challenge is that both c1 and c2 can undergo polymerization reactions to p1 and p2. This of course has an influence on their concentrations and therefore on kinetics. Fitting just reaction 1 (c1 to c2) and neglecting reactions 3 and 4 leads to extreme reaction orders. Also, during experiments the mass balance can never close properly if these reactions are to be neglected.

I have shared my code at the bottom of this question so you can see how I implemented the three reactions into MATLAB.

Available data:

Over a certain time period we have the concentration profiles of compounds c1 and c2. The amount of c3 is negligable, p1 and p2 can not be analyzed.

These profiles look like this:

This data is figurated however close to the experimental data. It should be well sufficient for this modeling challenge.

MATLAB Integration:

In MATLAB, I have used the lsqcurvefit method as described here: https://nl.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting. Thanks Star Strider for this starting point. I can now create a code that runs. However, the results from lsqcurvefit do not match the experimental data at all. Moreover, the results are different for every run.

My goal is to plot in figure 2 the experimental data and the model fit. As you will probably see, it is a mess although I am not sure why.

Your time and effort is very much apprieciated.

My code:

clear, clc

%Reaction system
%c1 --> c2 --> c3 (c2 --> c3 for now neglected)
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The last two reaction equations describe polymerization reactions, there
%is no data available for both p1 and p2.

%Data
time = [10 20 40 60 180];
c1 = [0.508 0.442 0.385 0.354 0.258];
c2 = [0.246 0.301 0.359 0.398 0.489];
c1_0 = 0.711;
c2_0 = 0;

%Plot
time_plot = [0,time];
c1_plot = [c1_0,c1];
c2_plot = [c2_0,c2];
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%Data processing
x_in = time';
y_in = [c1',c2'];

%Starting point for lsqcurvefit
p0 = [0.1,1,0.1,1,0.1,1];
lb=zeros(6,1);
ub=4*ones(6,1);

%Lsqcurvefit
[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@diff1,p0,x_in,y_in,lb,ub);
FitData = diff1(p,time,1);
c1_out = FitData(:,1);
c2_out = FitData(:,2);

%Plot
c1_out_plot = [c1_0,c1_out'];
c2_out_plot = [c2_0,c2_out'];
figure(2)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(time_plot,c1_out_plot,'b')
plot(time_plot,c2_out_plot,'r')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

function C = diff1(p,time,~)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%dc3/dt = k2*c2^n2 = 0
%variables: c1 = x(1), c2 = x(2)
%parameters: k1 = p(1), n1 = p(2), k3 = p(3), 
%n3 = p(4), k4 = p(5), n4 = p(6)
%for now we neglect reaction r2 of c2 to c3.

x0 = rand(2,1);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = -(p(1).*x(1).^p(2)+p(3).*x(1)^p(4));
          xdot(2) = p(1).*x(1).^p(2)-p(5).*x(2).^p(6);        
          dC = xdot;
      end
C = Cout;
end

Answer 1

[Star Strider]

Thank you for referring to my ’Monod Kinetics’ Answer! I’m glad you found it useful!

There are two problems with your code.

First, ‘xdot(1)’ should be:

xdot(1) = -(p(1).*x(1).^p(2)-p(3).*x(1)^p(4));
                            ↑

with the sign of the second term being negative (corrected here).

Second, it is a good idea to always include the initial values as parameters, as a rule, the last parameters. So here:

p0 = [0.1,1,0.1,1,0.1,1,c1(1),c2(1)];

and then incorporating those changes in ‘diff1’:

function C = diff1(p,time,~)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%dc3/dt = k2*c2^n2 = 0
%variables: c1 = x(1), c2 = x(2)
%parameters: k1 = p(1), n1 = p(2), k3 = p(3), 
%n3 = p(4), k4 = p(5), n4 = p(6)
%for now we neglect reaction r2 of c2 to c3.

x0 = p(7:8);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = -(p(1).*x(1).^p(2)-p(3).*x(1)^p(4));
          xdot(2) = p(1).*x(1).^p(2)-p(5).*x(2).^p(6);        
          dC = xdot;
      end
C = Cout;
end

the estimated parameters are:

p =
   0.105408507696630
   3.999984376449200
   0.000078501021266
   0.007282181170239
   0.002202264694465
   3.902023967475535
   0.502316713470111
   0.244896706878439

and your model fits your data almost exactly!

[John456]

Thanks for replying to my question on such short notice! After some thinking and evaluation of your remarks, I have made significant progress on my code.

To start off, your comment about the formule of 'xdot(1)'. It is true, when I change the sign, the code works and fits the data almost perfectly. However, I must disagree with your statement that the sign should be a minus. Mind the -(...) I am using. Using a minus sign would make the equation equal to:

xdot(1) = -p(1).*x(1).^p(2)+p(3).*x(1)^p(4);

This would, based on my reaction scheme, suggest that c1 is being formed from p1 which is definitely not the case, it's the other way around. Right?

I had to reconsider my data integration and found some made assumptions to be mistaken. After correcting my data sets, I got the system to fit nicely using the old set of equations, so using:

xdot(1) = -p(1).*x(1).^p(2)-p(3).*x(1)^p(4);

I still have a few questions though, cause the system is not perfect just yet. The code below is what I have at this point. I added some options, played around with initial guesses, trying to focus on the physical meaning of the outcome. In some cases, using a guess that is far off could result in lsqcurvefit not finding the right results, even if I use some extreme options.

My problems:

1. Using a corrected dataset (with a little more points than 5), I can properly fit the set of equations until t = 90. If I include t = 120 and t = 180, lsqcurvefit gives wrong results. I am not sure why, as these last datapoints do not seem far off compared to the others.
2. I would like to obtain a more smooth curve of the fit in figure 2. However, if I evaluate a different timespan time_2, the curve changes completely. Why is this and how can I solve this?
3. Eventually, I would like to determine for each fit a coefficient of determination. When changing my lower and upper boundaries, this will help to quickly see how good the fit is. In my experience with some simpler equations, I would integrate towards an analytical solution and use curve fitting to find parameters and the coefficient of determination. Maybe there is an easier, more convenient way to determine the coefficient of determination for a numerically obtained fit?

Many thanks for your time and help.

Current code

clear, clc

%Reaction system
%c1 --> c2 --> c3 (c2 --> c3 for now neglected)
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The last two reaction equations describe polymerization reactions, there
%is no data available for both p1 and p2.

%Data
time = [5 10 15 20 30 40 60 90];% 120 180];
c1 = [0.501 0.401 0.354 0.322 0.283 0.252 0.218 0.193];% 0.168 0.134];
c2 = [0.171 0.246 0.293 0.317 0.346 0.363 0.401 0.436];% 0.450 0.478];
c1_0 = 0.717;
c2_0 = 0;

%Plot
time_plot = [0,time];
c1_plot = [c1_0,c1];
c2_plot = [c2_0,c2];
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

%Data processing
x_in = time';
y_in = [c1',c2'];

%Starting point for lsqcurvefit
p0 = [0.1,1,0,1,0,1,c1(1),c2(1)];
lb = [0,0,0,0,0,0,-Inf,-Inf];
ub = [1,5,1,5,1,5,Inf,Inf];

%Lsqcurvefit
options = optimoptions('lsqcurvefit','StepTolerance',1e-100,'MaxFunctionEvaluations',10000,'MaxIterations',10000);
[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@diff1,p0,x_in,y_in,lb,ub,options);
time_2 = time;
%time_2 = 0:1:60;
FitData = diff1(p,time_2);
c1_out = FitData(:,1);
c2_out = FitData(:,2);

%Plot
time_2_plot = [0 time_2];
c1_out_plot = [c1_0,c1_out'];
c2_out_plot = [c2_0,c2_out'];
figure(2)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(time_2_plot,c1_out_plot,'b')
plot(time_2_plot,c2_out_plot,'r')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

function C = diff1(p,time)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%dc3/dt = k2*c2^n2 = 0
%variables: c1 = x(1), c2 = x(2)
%parameters: k1 = p(1), n1 = p(2), k3 = p(3), 
%n3 = p(4), k4 = p(5), n4 = p(6)
%for now we neglect reaction r2 of c2 to c3.

x0 = p(7:8);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = -p(1).*x(1).^p(2)-p(3).*x(1)^p(4);
          xdot(2) = p(1).*x(1).^p(2)-p(5).*x(2).^p(6); 
          dC = xdot;
      end
C = Cout;
end

[Star Strider]

My pleasure.

The sign I changed reflects the observation that everything is leaving ‘c1’, so all signs should be negative with respect to it.

I have no problems at all fitting your expanded data to my code, including the values at 120 and 180 minutes. The parameter values I get are:

p =
   0.208398498742923
   3.296912877018929
   0.000007461661372
   0.147667656984559
   0.005612972034412
   3.633461329118330
   0.477083265213336
   0.164431838417852

The ‘xdot(1)’ equation you mentioned using:

   xdot(1) = -p(1).*x(1).^p(2)-p(3).*x(1)^p(4);

is the same one I use in my code, with the change I suggested.

With respect to the problems you mentioned:

1. I got a good fit with the additional data, again using my code. You can’t expect an exact fit to experimental data, and I disagree that the two aded data are in any way ‘far off’. Note that you fit 5 minutes to 120 minutes, and then in your plot include the 0 minute points. I will let you ponder the wisdom of that approach. The parameters will only provide a good fit for the data you present to whatever optimisation function you use, however fitting your data with the 0 minute points does not work well (at least wioth lsqcurvefit). I’ve not experimented with using it with ga, so you may need to revise your model.
2. If you use linspace to generate a curve with more points, the curve is smoother and essentially matches you original data. I had no problems with it providing that I used the existing limits of your ‘time’ vector:

plottime = linspace(min(time), max(time), 150);

3. The coefficient of determination is easy enough to calculate:

SStotC1 = sum((c1 - mean(c1)).^2);
SSresC1 = sum((c1(:) - c1_out).^2);
R2C1 = 1 - SSresC1/SStotC1

SStotC2 = sum((c2 - mean(c2)).^2);
SSresC2 = sum((c2(:) - c2_out).^2);
R2C2 = 1 - SSresC2/SStotC2

producing (with my code):

R2C1 =
   0.991761664509532

R2C2 =
   0.984285221236630
   

It may be worthwhile to consider using the genetic algorithm (ga) function if you want a robust fit. It does not use gradient descent, so will generally find good parameter estimates for otherwise difficult-to-fit data.

Since my Answer apparently did not solve your problem, I will delete it in a few hours.

[John456]

Hi,

Thanks for elaborating on the points I mentioned.

So then we are both on the same page, using the equation:

xdot(1) = -p(1).*x(1).^p(2)-p(3).*x(1)^p(4);

That's good.

1. Yes, I mentioned. These two last datapoints do not seem off in any way, so I found it strange that the curvefit stopped telling me it reached its maximum amount of iterations or function evaluations. However, if I now use your results as my estimates, I can solve the system very quickly and find parameters for the system including the 120 and 180 points. Guess I need to pay more attention at what is going in there.
2. That worked perfectly at first try, thanks!
3. This as well.

What I want to do is step-wise change my estimates and boundaries to values that are more likely to occur in practise. We are not going to find extreme reaction orders, for example, at least not for the reaction c1 --> c2. When decreasing the upper limit for this reaction order, I obtain, of course, different parameters and a lower coefficient of determination. Like you mentioned, there is no use in trying to obtain a perfect fit from experimental data. On the other hand, I would like use the coefficient of determination to get a fast idea of how well a fit still fits the data when I start changing the boundaries.

I will look shortly into the ga method. For now I decided to include the 0 points because all the fits I would obtain for different sets of data and equations would never be able to properly predict the concentrations at time 0 within a reasonable range.

Lastly, what we have accomplished now with lsqcurvefit works. We can obtain a fit, read parameters, and re-evaluate our parameters by changing boundaries. I think hereby you have answered not only my question, but also helped me to quickly develop my code into something useful. Thanks a lot.

[Star Strider]

As always, my pleasure!

I had some problems yesterday getting ga to converge successfully. I intend to keep working on it, and if successful, will post the relevant code and results to a Comment here.