Estimating parameter values of a differential equation

Question

I am trying to fit the following ODE system which describes drug absorption and elimination.

If $A_e$ is the amount of drug in the absorption state and $A$ is the amount of drug in the body this can be described as:

$$\begin{aligned} {dA_e \over dt} &= -k_a A_e\\ {dA \over dt} &= k_a A_e - k A\end{aligned}$$

I am trying to estimate the parameter values the absorption rate ($k_a$) and the elimination rate ($k$).

Below is the Matlab code written to estimate the parameter values using lsqnolin. But I am getting warnings of RelTol has been increased to 2.2E-14. Can someone please tell me what is wrong with the way I have modelled it?

I want to know how to estimate using an ODE System (That is without writing the objective function as the solution of $A(t)$ as the actual model has many other states)

time=[0.5;1;1.5;2;3;4;5;6;8;12];%time in hours

Drugdata=[0.33;1.23;1.95;2.72;3.51;3.63;3.47;3.22;2.39;1.13];

%parameters estimated:k_a and k
initial=[0.5,0.5];%initial guesses of the estimates

ODEinit=[8.48,Drugdata(1)];%initial values to the ODE. 

lb=[0,0];
ub=[100,100];
[fittedVals,resnorm,~,~,~,~,~]=lsqnonlin(@(xEstimate)errorFun(xEstimate,Drugdata,ODEinit,...
    time),initial,lb,ub);

fittedValues=model2(fittedVals',time,ODEinit);


close all
figure
h2=scatter(time,Drugdata,20,'b','filled');
hold on
h1=plot(time,fittedValues,'color', 'r');
ylabel('Drug')
xlabel('time')

function err=errorFun(xEstimate,DrugMeas,ODEinit,thours)
    Drug_est=model2(xEstimate,thours,ODEinit);%estiamted values of drug using the ODE model
    err=DrugMeas-Drug_est;
end


function Drugout= model2(xEstimate,thours,ODEinit)


options = odeset('AbsTol',1e-25,'RelTol',1e-25);
[~,values] = ode45(@(t,y)Equations(t,y),thours,ODEinit,options);

function s=Equations(~,y)
    %parameters estimated:k_a, k

    k=xEstimate(2);%rate of elimination
    ka=xEstimate(1);%rate of drug absorption

    s=zeros(2,1);
    s(1)=-ka*y(1);%absorption compartment
    s(2)=ka*y(1)-k*y(2);%drug

end
Drugout=values(:,2);
end

Answer 1

You get a warning because in the usual 64bit floating point format a relative tolerance of 1e-25 is not achievable. Rounding real numbers to floating point already gives a relative error of up to 2e-16, add to that the accumulation over the computations of the solution method and a relative error of 2.2e-14 appears as justified for a minimal case for a positive value of it. Your data has 3 digits, thus a working precision of 6 digits for the integration seems sufficient. To be more than exact, set RelTol to 1e-9 and AbsTol to 1e-7.

Doing experiments with at translation of the code to python using scipy.optimize.least_squares, I get that

• Using the minimal (that is, largest) tolerances atol = 1e-4; rtol=1e-5; gives parameter estimates [0.25206566 0.30065899]
• Using moderate tolerances atol = 1e-5; rtol=1e-7; gives parameter estimates [0.25207333 0.30065716]
• Using extreme tolerances atol = 1e-12; rtol=1e-13; gives parameter estimates [0.25207395 0.30065723]

So as can be seen, the first 3 or even 4 digits of the estimated parameters do not depend on the integration tolerances.

The plot of the fitted model is

Setting the initial time at $0$ and making the initial values also variables of the optimization process results in the fitted parameter vector $$ [k, k_a, A_{e0}, A_0] = [ 0.281374,\, 0.28137101,\, 10.95249421,\, -1.34117042] $$ with the plot

Forcing $A_0$ to be positive gives a much worse fit.

Answer 2

This was not an answer as expected. This was a method to help in investigating the origin of the trouble.

Meanwhile Lutz Lehmann gave the clue which in fact was a matter of specification of the working precision for the software.

Nevertheless it is of interest to compare the results of both methods : Optimization of the parameters of the system of ODEs or optimization of the parameters of the function which is solution of the ODEs.

We will analytically solve the system of ODEs :

$${dA_e \over dt}=-k_a A_e\quad\implies\quad A_e=c\:e^{-k_a t}$$

$${dA \over dt}=k_a A_e-k\:A=k_a c\:e^{-k_a t}-k\:A$$ Solving the first order linear ODE $\quad {dA \over dt}+k\:A=k_a c\:e^{-k_a t}\quad$ leads to : $$A(t)=c\:\frac{k_a}{k-k_a}e^{-k_a t}+Ce^{-k\,t}\tag 1$$ $c$ and $C$ are two constants.

With $A(0)=0\quad\implies\quad C=c\:\frac{k_a}{k-k_a}$ $$A(t)=C\left(e^{-k\,t}-e^{-k_at} \right)\tag 2$$

Thus this function must approximately fit your data. Try a regression software to compute the approximate values of $k$ , $k_a$ and $C$.

If this method fails like the method of differential equations fitting, one would think that the model of ODEs is not fully convenient with regard to the data. If it is successful you obtains the approximate values of $k$ and $k_a$ which solves the problem.

EXAMPLE :

The mean least squares fitting for the four parameters of above Eq.$(1)$ which is equivalent to : $$A(t)=c_1e^{-k_a t}-c_2e^{-k\,t}$$ gives : $$k_a\simeq 0.251\quad;\quad k\simeq 0.335\quad;\quad c_1\simeq 37.961\quad;\quad c_2\simeq 39.450\quad;$$ RMSAE $\simeq 0.0879$

The results for $k_a$ and $k$ are slightly different to the Lutz Lehmann's results : $\quad k_a\simeq 0.252\quad;\quad k\simeq 0.301$

NOTE : The mean least squares fitting for the three parameters of the above Eq.$(2)$ is less good. So, assuming that $A(0)=0$ was not a good idea.