parameter estimation in systems of ordinary differential equations

[MarceloLeite]

Hey guys. I'm trying to estimate parameters in a system of ordinary differential equations.
Here is the code I wrote (I know, is a disaster), with the error messages. Can anyone tell me how to fix the code?

title 'ODE BATALIM';
data BATALIM;
input     hour V   CELL    S     P;
datalines;
0.0 1.5 61.64 0. 30.39
1.0 2.2 34. 43. 31.18
2.0 2.9 36.34 56.06 35.95
3.0 3.6 21.92 66.29 38.87
4.0 4.3 17.75 69.12 42.31
5.0 5.0 17.01 70.14 44.17
;
proc model data=BATALIM;
dependent   V    CELL   S    P;

parms   Fa   mimax   Ks   Ki   Pmax   CELLmax   n    m;

/*Next I try to provide the initial conditions of the differential equations*/

if (hour=0) then
      do;
       V=1.5; CELL=61.64; S=0; P=30.39;
       end;
else
      do;
/*The system of differential equations is described below*/

      dert.V = Fa;
      dert.CELL = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLXmax)**m)*CELL - Fa*CELL/V;
     dert.S = Fa*(213-S)/V - (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL/0.017;
     dert.P = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL*24.94 - Fa*P/V;
end;

/*

Then I try to insert the initial values ​​of the parameters and get the fitted equations...

*/

fit V start=(Fa 0.5) / time=hour;
fit CELL start=(mimax 0.01 Ks=2 Ki=5 Pmax=40 n=2 CELLmax= 60 m=2) / time=hour;
fit S start=(Fa 0.5 mimax 0.01 Ks=2 Ki=5 Pmax=40 n=2 CELLmax= 60 m=2)/ time=hour;
fit P start=(mimax 0.01 Ks=2 Ki=5 Pmax=40 n=2 CELLmax= 60 m=2)/ time=hour;
run;

Error messages

WARNING: Missing value encountered for initial dert.CELL at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.S at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.P at time = 0. This variable is set to zero.

NOTE: At OLS Iteration 0 CONVERGE=0.001 Criteria Met.

WARNING: Missing value encountered for initial dert.CELL at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.S at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.P at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.CELL at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.S at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.P at time = 0. This variable is set to zero.

And more...

[SteveDenham]

I am probably answering the wrong question here, but I noticed that the derivative of CELL at hour=0 is, in turn, equal to zero.  That may be leading to the missing values. Looking at the example on fitting systems of ODEs in the PROC MODEL documentation, I see that it is set up to only calculate derivatives if 'time ne 0' and uses a quadratic solution for time=0.  I think that particular approach is designed to avoid the warnings that you are getting with your code.

The other issue is what you copied from the log - it appears that the process is "stuck" at hour=0, but I can't tell.  The repetitive nature makes me think that the same error is generated at each iteration.

Given these you ,might want to reparameterize your system of equations, if at all possible.  The derivatives look unwieldy, so I suspect the partials with respect to the parameters may be even worse.

SteveDenham

[MarceloLeite]

Hello Steve, how are you? Thanks for responding to my help request.

"I am probably answering the wrong question here, but I noticed that the derivative of CELL at hour=0 is, in turn, equal to zero.  That may be leading to the missing values... I see that it is set up to only calculate derivatives if 'time ne 0' ..."

This is a problem. In various physical phenomena dy/dt (x = 0) =0.

"The other issue is what you copied from the log - it appears that the process is "stuck" at hour=0, but I can't tell.  The repetitive nature makes me think that the same error is generated at each iteration."

I agree with you. The problem is that I can't fix it.

"Given these you ,might want to reparameterize your system of equations, if at all possible.  The derivatives look unwieldy, so I suspect the partials with respect to the parameters may be even worse."

I was also suspicious of the equations. However, I have consulted the literature and they are correct. I've tried everything here and I don't know how to solve it.

I will try to find a SAS code for solving EDOS systems that works and test it with my data.

[Rick_SAS]

> I was also suspicious of the equations. However, I have consulted the literature and they are correct.

Great. Please provide the reference to the literature that you consulted. The paper, article, or blog post will tell us what problem you are trying to solve, and we might be able to find the source of your problem.

[MarceloLeite]

Hi Rick!

The literature I am consulting is in Portuguese. In addition, it is very long, as it is a master's thesis. However, I have an article (unfortunately, also in Portuguese) with only 15 pages that portrays exactly what I want to do: simulate fermentation processes, both in batch and fed batch. If you want, I can send it. I'll understand if you don't want to.

[Rick_SAS]

Unfortunately, I would not be able to understand it. I wish you good luck in solving your problem.

[Rick_SAS]

I am not an expert on PROC MODEL, but I noticed two things:

1. You use the parameter CELLXmax, which is not defined. Presumably, that is a typo and you meant CellMax.

2. The parameter values for CellMax and PMax are less than the data values.  From the names, I assume that Cell < cellMax and P < PMax for all times? If so, change those initial parameter values:

parm   Fa 0.5  mimax 0.01  Ks 2  Ki 5  
       Pmax 50  CELLmax 65  n 2   m 2;

3. Do you need to add a BOUNDS statement? It seems like a lot of the parameters converge to large negative numbers. Check the model to see if it constrains any parameters.

4. I don't think you need the START= options on the FIT statements, but I don't really know what you are trying to accomplish.

Sometimes getting something that runs is the first step to debugging other problems. Here is some code that runs, but it does not provide a good fit to the data. I encourage you to check the ODE formulas (are all the signs correct?).  

title 'ODE BATALIM';
data BATALIM;
input     hour V   CELL    S     P;
datalines;
0.0 1.5 61.64 0. 30.39
1.0 2.2 34.   43. 31.18
2.0 2.9 36.34 56.06 35.95
3.0 3.6 21.92 66.29 38.87
4.0 4.3 17.75 69.12 42.31
5.0 5.0 17.01 70.14 44.17
;
proc model data=BATALIM plots=none;
dependent V 1.5 CELL 61.64 S 0. P 30.39;

parm   Fa 0.5  mimax 0.01  Ks 2  Ki 5  
       Pmax 50  CELLmax 65  n 2   m 2;
/* do you need to add BOUNDS????? */
Bounds 45 < PMax , 65 < Cellmax;
if (hour=0) then
      do;
       V=1.5; CELL=61.64; S=0; P=30.39;
       end;
else
      do;
/*The system of differential equations is described below*/
dert.V = Fa;
dert.CELL = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL - Fa*CELL/V;
dert.S = Fa*(213-S)/V - (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL/0.017;
dert.P = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL*24.94 - Fa*P/V;
end;

fit V CELL S P / time=hour out=ModelOut;
run;

/* visualize the results */
data results;
merge BATALIM(rename=(V=VObs CELL=CellObs S=SObs P=Pobs)) ModelOut;
by hour;
run;
ods graphics/width=400px height=300px;
proc  sgplot data=results;
scatter x=Hour y=VObs;
series x=Hour y=V;
run;
proc  sgplot data=results;
scatter x=Hour y=CellObs;
series x=Hour y=Cell;
run;
proc  sgplot data=results;
scatter x=Hour y=SObs;
series x=Hour y=S;
run;
proc  sgplot data=results;
scatter x=Hour y=PObs;
series x=Hour y=P;
run;

[MarceloLeite]

Hello how are you?

First of all, thank you very much for your help.

I tried to write a code to estimate the parameters of a system of ordinary differential equations. The result was that mess you saw.

The parameter values

Fa 0.5  mimax 0.01  Ks 2  Ki 5  Pmax 50  CELLmax 65  n 2   m 2

they are only initial estimates, since the estimation method is iterative. Pmax is expected to be close to the largest value of P in the dataset (slightly larger or slightly smaller). The same goes for CELLmax.

I used a system of equations taken from a master's thesis. Perhaps the author made a typo. I will check them by consulting the original source.

Thanks for writing and sending me code that works!

[clayt85]

I'd like to add a bit more perspective to the comments from @Rick_SAS that I think might assist as you debug this program.

If you take the program as supplied by @Rick_SAS but make the following modifications:

1. Remove the BOUNDS statement
2. Set the value of PMAX to 40 (as in the original code)
3. Set the value of CELLMAX to 60 (as in the original code)

Then you will receive exactly the same error ("Missing value encountered for intial dert.CELL at time = 0..." And so on.) I'm making an educated guess here (I do not know much about PROC MODEL) but I am quite certain this is the reason:

As @Rick_SAS notes, the values of CELL and P in the data exceed the values of the parameters CELLMAX and PMAX. In the most general case, this need not necessarily be a problem. (For example, a measured value might exceed a theoretical maximum because of random error during the measurement process.) BUT in your case this is almost certainly a problem.

To see why, notice the structure of the derivative equations (for example, dert.CELL). The equation contains the term "(1 - P/Pmax)**n". If P>Pmax, then the term inside the parentheses is negative. This is OK when n=2 exactly (as you specify for the initial value of the parameter). But this is not defined when, say, n=2.1. This latter fact is the reason for the error. For a nonlinear least squares problem, the derivatives of the objective function with respect to the parameters (these are computed internally by PROC MODEL) must be well-defined in an open region. For your case, the derivatives are defined only at a point (n=0).

Depending on the application, you might consider fixing the values of the parameters m and n.

[MarceloLeite]

Hello. I found your remarks interesting. 
However, I did what you suggested, but fitting the curves to the experimental data remained unacceptable...