Finding Unknown Parameters from Monod Equations

Hi, new to Matlab and still learning.

I am using a modified code provided by Star Strider.Thanks Star Strider! I am making this as a new post so I can accept this one as well =).

The code is shown below. I got the data from another Matlab post. The equation is as shown in the screenshot: I converted all constants into theta. When I run the code below, I get huge constants that doesn't makes sense. Anything I'm doing wrong? Is it my equations? C0 or theta0? Thank you!

Function Igor_Mour

function C=kinetics(theta,t)

c0=[0.1;9;0.1];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(3,1);

dcdt(1)=-theta(1).*c(1).*c(2)/(theta(2)+c(1));

dcdt(2)= -(1/theta(3)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));

dcdt(3)= (1/theta(4)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));

dC=dcdt;

end

C=Cv;

end

t = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]'

x = [0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093]';

s = [9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0]';

p = [0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869]';

%Note: x,s and p are combined (in separate columns) to form c.

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

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

fprintf(1, '\tRate Constants:\n')

for k1 = 1:length(theta)

fprintf(1, '9\t\tTheta(%d) = %8.5f\n', k1, theta(k1))

end

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

figure(1)

plot(t, c, 'p')

hold on

hlp = plot(tv, Cfit);

hold off

grid

xlabel(Time)

ylabel(Concentration)

legend(hlp, 'C_1(t)','C_2(t)','C_3(t)','C_4(t)', 'Location','N')

end

[Alex Sha]

Hi, 0.07662, one of data for p, should be 0.7662? if so, see the result below:

Sum Squared Error (SSE): 1.43925599360388
Root of Mean Square Error (RMSE): 0.208839215637285
Correlation Coef. (R): 0.982548927303979
R-Square: 0.9654023945462

Parameter	Best Estimate     
---------	-------------     
theta1   	0.0214068471653595
theta2   	-1.28632596117882 
theta3   	-0.113732223600425
theta4   	-1.55421915737194 

Answer 1

[Star Strider]

The estimated value for ‘theta(1)’ definitely does not appear to be as expected for a rate constant, however if it is than it could be appropriate. I am not certain what the other parameters are by comparing the code to the symbolic expressions It would help to have them clarified.

According to the posted symbolic differntial equations, the first differential equation should not be negated. Changing that gives a decent fit to the data. Beyond that, I suggest checking to be certain that the equations are coded correctly. Here, I have ‘x’ as ‘c(1)’, ‘s’ as ‘c(2)’, and ‘p’ as ‘c(3)’. Change those if necessary to be the correct columns.

       t = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]' ;
       x = [0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093]'; 
       s = [9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0]';
       p = [0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869]';
       %Note: x,s and p are combined (in separate columns) to form c.
       c = [x s p];
theta0=rand(7,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(1,numel(theta0)));

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

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

		Theta( 1) = 432.22494
		Theta( 2) = 11061.84689
		Theta( 3) =  0.13367
		Theta( 4) =  2.38808
		Theta( 5) =  0.02007
		Theta( 6) =  8.88929
		Theta( 7) =  0.05209

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

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

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

Please define the parameters in terms of the symbolic equations. With that information, it will be easier to determine if there are any inconsistencies.

.

[Learning]

Thanks for your responses! not sure of the symbolic equations you are referring to but they were written as shown below.

OR:

dX/dt = umax*a(1)*a(2)/(ks + a(2));

dS/dt = -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2));

dP/dt = (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];

I made umax as K(1), ks as K(2), Yxs as K(3), and Yps as K(4), s as c(1), and x as c(2).

I found a typo in my diff equations in my question. The equations should be

dcdt=zeros(3,1);

dcdt(1)=-theta(1).*c(1).*c(2)/(theta(2)+c(1));

dcdt(2)= -(1/(theta(3))).*theta(1).*c(1).*c(2)/(theta(2)+c(1));

dcdt(3)= (1/(theta(4))).*theta(1).*c(1).*c(2)/(theta(2)+c(1));

dC=dcdt;

Few questions

1. From the equations above, I have only s and x and hence I'm thinking only 2 concentrations, c(1) and c(2)? where does the c(3) arises from? from dP/dt?
2. How do you choose the length of theta0, you have theta0 = [1;1;1;1;1;1;1]. How do you know you have to have 7 1s and not 4 1s?
3. When you were plotting the data after estimating theta's you made c0=theta(5:7); is it wrong if I made c0 be the same as the original initial concentrations?,

Thanks!

[Star Strider]

The first ‘dcdt(1)’ should not be negated, according to the symbolic equations. I now understand that ‘theta(1)’ is however ‘c(1)=s’ and ‘c(2)=x’ changes the code so that it now does not work.

Copying and pasting your corrected equations and changed compartments, I get this —

       t = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]' ;
       x = [0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093]'; 
       s = [9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0]';
       p = [0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869]';
       %Note: x,s and p are combined (in separate columns) to form c.
       c = [s x p];
theta0 = [rand(4,1); s(1); x(1); p(1)];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(1,numel(theta0)));

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

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

		Theta( 1) =  0.00014
		Theta( 2) = 973.54767
		Theta( 3) = 752.06527
		Theta( 4) =  0.00002
		Theta( 5) =  4.01626
		Theta( 6) =  0.70152
		Theta( 7) =  0.03684

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

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

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

I suspect that the new compartment assignments:

c = [s x p];

are incorrect.

Please experiment with it to get the differential equations to be correct, and the appropriate compartment assignments to be correct. The compartments in order are the column vectors in the ‘c’ matrix.

I will do my best to help you get this to work, however I am limited by not knowing the equations you want to code in the context of the compartments you want to model.

.

[Learning]

Final questions:

1. I think I now understand why you are using theta0 =rand(7,1). Is 7 due to the fact that I have 4 unknowns parameters (theta(1), theta(2), theta(3) and theta(4)) and 3 concentrations in c? In this case it appears the last 3 theta values in the estimated theta are the estimated initial concentrations for the provided data c? If I already know the initial concentrations c0, do I still have to estimate them? If I did not want to estimate them then does theta becomes rand(4,1)?

2. I see that when plotting the original data (t,c) on the estimated data (tv, Cfit), you are still using theta in the function C=kinetics(theta,t) and Cfit=kinetic(theta, tv). In this case if the estimated initial concentrations theta(5:7) are wrong, then the graph will be wrong. So can I just use the actual unknown parameters theta(1) to theta(4) in the kinetic function instead if I already know my c0 and not the estimated c0’s?

Thank you again.

[Star Strider]

As always, my pleasure!

‘Is 7 due to the fact that I have 4 unknowns parameters (theta(1), theta(2), theta(3) and theta(4)) and 3 concentrations in c? In this case it appears the last 3 theta values in the estimated theta are the estimated initial concentrations for the provided data c?’

Yes.

‘If I already know the initial concentrations c0, do I still have to estimate them? If I did not want to estimate them then does theta becomes rand(4,1)?’

For reasons not entirely clear to me, the code works best if the initial concentrations are parameters to be estimated. I generally use the initial concentrations as the initial estimates for those parameters. (I suspect that is because they are data like all the rest, and are subjkect to some measurement variation.)

‘So can I just use the actual unknown parameters theta(1) to theta(4) in the kinetic function instead if I already know my c0 and not the estimated c0’s?’

Yes, although I use the initial concentrations as the initial estimates for those parameters when I estimate them as parameters themselves. The estimated initial concentrations are generally close to the initial concentrations in the data.

.

[Learning]

Awesome!

1. I think I know why the plot in your second response and estimated theta and initial concentrations were did not look correct…the c compartment is kinda flipped…c should be [x s p] but you have it set to [s x p].

2. Do you suggest it’s best to estimate theta’s for each dataset or best to save the estimated theta’s and use it on new data? Assuming data is flowing into the code in real time.

3. Should Cfit then be written as Cfit=kinetics(theta(1:4),tv))? for the case where the first 4 values in the estimated theta are the actual unknown parameters being estimated.

3. Is it also possible to calculate the residuals or errors between the Cfit and the actual concentrations in c at each time point? (assuming linespace is not used).

Thank you!!

[Star Strider]

Thank you!

1. It originally was:

c = [x s p];

I changed it because in your Comment you wrote: ‘I made umax as K(1), ks as K(2), Yxs as K(3), and Yps as K(4), s as c(1), and x as c(2).’

2. That is your call. If the new data matches the earlier data, then that would work. If not (for example if the temperature or other conditions were changed), the system parameters would likely need to be estimated from the new data.

3. That depends on whether you are estimating new data. If you are using the previously estimated parameters with the new data, and the system conditions are similar, yes.

4. The residuals can be calculated from the original data as:

resid = c - kinetics(theta, t)

That will produce a matrix of the residuals the same size as ‘c’ that you can then do some basic stitistics with.

You can use nlparci to get confidence intervals on the parameters of the‘kinetics’ function, however getting confidence intervals on the fit are (to the best of my knowledge) not possible with a function such as ‘kinetics’ because when I last looked (and just now checked the fitnlm documentation to be sure), the Statistics and Machine Learning Toolbox regression functions did not allow error estimation of functions with multiple outputs. I do not have the Curve Fitting Toolbox, however checking the documentation for fit just now, it apparently has the same restriction.

.

[Learning]

Thank you!!

[Star Strider]

As always, my pleasure!

[Learning]

Happy Thanksgiving! Quick question: I have used the same code in the previous discussions on another set of ODEs but it takes a VERY long time for results to be obtained…and even then the fit to the original data is not great. 1. What could be causing the slow execution of the code and 2. Could the poor fit be due to the fact that the ode equations are not correct or that the data simply do not fit the equations? Any other fit methods?

Thanks!

[Star Strider]

Happy Thanksgiving to you, too!

In all likelihood, the reason for the the code taking a long time is likely that the differential equations are ‘stiff’. One reason for this is the parameters differing by several orders of magnitude. Using ode15s (or other stiff solvers such as ode23s) should significantly shorten the time required to integrate them.

There could of course be other reasons (I’d have to have all the information — code and data — to be certain), however this is the most likely possibility in my experience.

.

[Learning]

Thanks for your response! My final question is how do I implement this in real time? Say I identify the parameters for the series of odes, do I need to have another set of data to estimate the parameters? Or I can fix the parameters if I know the process doesn’t change that much? Real time in the sense that assuming I have live stream data from a hardware that I want to compare with the concentrations estimated by the ode at specific times.

Thank you!

[Star Strider]

As always, my pleasure!

Updating the parameter estimates in real time is not something I’ve considered. I’ve not done real time parameter estimation in decades.

I’m not sure that a real time application of this is possible, however it depends on what ‘real time’ actually means. If there is time between samples to numerically integrate the differential equations, and good estimates of the existing paramerers are known, then the current approach could be used without alteration, using the existing estimates for the initial estimates for parameter estimates with the additional data.

Otherwise, a Kalman filter approach could be possible, however I’ve little experience with Kalman filters so I don’t know if one could be used in this instance. There are a large number of Kalman filter functions in MATLAB that might be possible to use. One such is kalman, and there are more than 50 others in the online documentation that come up with a search for ‘kalman’.

Answer 2

[Torsten]

Works for me.

t = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]' ;
x = [0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093]'; 
s = [9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0]';
p = [0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.7662 0.7869]';
c = [x s p];
theta0=[1;1;1;1;1;1];
options = optimset('MaxFunEvals',100000,'MaxIter',100000);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,[],[],options);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

theta

theta = 6×1

-3.5510 133.8605 -0.1503 -1.9775 1.0000 1.0000

C = kinetics(theta,t)

C = 12×3

0.0904 9.0115 0.0151 0.1811 8.4079 0.0610 0.2782 7.7620 0.1101 0.4889 6.3604 0.2166 0.6186 5.4979 0.2822 0.8008 4.2856 0.3743 1.0126 2.8768 0.4814 1.1491 1.9684 0.5505 1.2507 1.2924 0.6019 1.3212 0.8237 0.6375

hold on
plot(t,c,'o')
plot(t,C)
hold off

function C=kinetics(theta,t)
%c0=[0.1;9;0.1];
c0 = [0.0904;9.0115;0.0151];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
        dcdt=zeros(3,1);
        dcdt(1)=-theta(1).*c(1).*c(2)/(theta(2)+c(1));
        dcdt(2)= -(1/theta(3)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
        dcdt(3)= (1/theta(4)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
        dC=dcdt;
    end
C=Cv;
end