why PSI values are diffrent when nlmefit using ode45 solver compared to when it uses analytical solution?

nlmefit estimated parameters differently when it uses analytical solution compared to when it uses ode45 solver for the same dataset. PSI values are completely underestimated.

[beta,PSI,stats,b] = nlmefit(timeNew,DVnew,IDnew,dosePerID,model,logparameter0,...

'REParamsSelect',[1 2 3 ],'ErrorModel','Constant','ApproximationType','FOCE','CovPattern',P,'RefineBeta0','on',...

'ParamTransform',repelem(1,3),'OptimFun','fminunc');

function concentration=Getsolutionode(parameter,t,doseNew)

%% initial values

V=parameter(:,1);

ka=parameter(:,2);

CL=parameter(:,3);

predictor(:,1)=t;

minTime=min(t);

maxTime=max(t);

y0=[doseNew 0];

%% ode45

options=odeset('RelTol',10-9,'AbsTol',10-9);

solode=ode45(@odefunction,[minTime maxTime], [y0],options,V,ka,CL);

rspl=deval(solode,t);

concentration=rspl(2,:)'./V;

end

function dxdy = odefunction(~,y,V,ka,CL)

kel=CL/V;

dxdy= [-ka*y(1)

ka*y(1)-kel*y(2)];

end

[Arthur Goldsipe]

Is that the actual code? It looks like you've set the ODE solver tolerances to 10-9, which evaluates to 1. I assume you intended to set the tolerances to 1e-9 instead.

If possible, can you share complete reproduction steps, including data? I'd like to try to solve the problem using SimBiology.

[Sepi Sharif]

Hi I meant 10-9, you are right. so, the good news I did the dataset, in Simbiology (first order absorption, CL,v for elimination, one compartment)and I have a good estimations for the parameter and PSI aligned with nlmefit algabraic coding that I did.

If I knw how to do ode in Simbiology!? then I can see their coding output. but I havent done anything there. so I attached my dataset and my coding.

[Sepi Sharif]

Hi Arthur, any chance to try this, maybe?

[Star Strider]

@Sepi Sharif — Do you already have SimBiology? If so, please add it to the ‘Products’ section.

[Arthur Goldsipe]

@Sepi Sharif I haven't had time to look at this yet. Your replies arrived while I was sleeping. But I should have some time later today.

[Arthur Goldsipe]

@Sepi Sharif I just tried running your sample code. It seems to depend on a MAT file warfpk1 that was not attached. I can guess what it looks like based on the CSV file you attached, but I'd prefer to use the actual file if you can attach that, too.

[Sepi Sharif]

Hi Arthur,

Thanks for your time and help. Your reply arrived when I was sleeping :).

The policy here, says I cannot attach more files unless I waite 24 hours. Do you think I can have your email adress or shall I waite.

[Star Strider]

For the record, I am also interested in following this, so I would prefer that everything appear here.

[Sepi Sharif]

Hi Star Strider,

Sounds right. so, then I will attach my Simbiology codes as well.

[Sepi Sharif]

So, here are the mat. file dataset attached, plus, Simbiology codes and the saved projects.

I did sort out the problem, so ode 45 has to have tspan, the minimum of tspan always should be zero, then ka = absorption rate cant estimated wrongly, so for each subject, they have their own maxtime of tspan but minimum is always zero. also by choosing different type of error and type of approximations in nlmefit command, it is a trade off between speed and accuracy.

Might be useful for the future people :)

[Sepi Sharif]

Hi all,

a question, about nlmefit, approximation type: so, dose LME as the default mean that nlmefit also, can solve linear mixed effects as well? say a systme of linear ODEs could be solved by nlmefit?

Thanks for your help

Answer 1

[Star Strider]

From the description, it appears that the code is attempting to fit data to the solution of a differential equation, or a system of differential equations.

To do this with ode45 (or more appropriately ode15s for kinetic equations), see Coefficient estimation for a system of coupled ODEs for a useful illustration of a similar problem.

.

[Sepi Sharif]

Thanks for your reply. My problem is the estimated PSI values which are diffrenet from nlmefit for analytical solution compared to ode45 solver.

it is estimating parameters for a nonlinear mixed effect model, which makes different values of variance_co-variance matrix output, Psi, for ode and alnalytical solution.

[Star Strider]

As always, my pleasure!

The PSI values should be quite similar using the numeric ODE solver and the analytic solution. They may not be exactly the same, owing to inherent differences in the iterative and analytic solutions.

I never tried to use nlmefit with nonlinear parameter estimation, using the iterative differential equation solvers, so my only experience is with lsqcurvefit (that can fit matrix dependent variables) and ga (that is usually robust at finding the correct parameters). The approach in the Answer that I reference here is generally succesful in situations like the one presented. It also estimates the initial conditions of the differential equations as parameters, as well as the parameters in the differential equations, so the number of parameters will likely be different between the analytic and iterative solutions (unless tha analytic solution also estimates the initial conditions).

It could be worth beginning with lsqcurvefit to get acceptable estimates of the parameters, then use those as ‘beta0’ with nlmefit. The ‘model’ would be the ‘kinetics’ function.

I willl help as I can, however my experience with nlmefit is limited and not at all recent, so I will also be experimenting.

.

[Sepi Sharif]

Thanks alot for your time and help. Appreciate it

[Sepi Sharif]

I am developping one numerical solver in Matlab and my refrrence solutions are the analytical and then ode45 solver to see if the outpts are align together with NONMEM.

Actually, analytical solution has a very good out put align with NONMEM programm (the best robust one for this kinetic nonlinera mixed effect models).

I am new to nlmefit and as far as I have searched, there are similar questions but very limited answers or nothing at all.

[Star Strider]

As always, my pleasure!

I had not heard of NONMEM previously. I looked through the available online MATLAB documentation, specifically SimBiology, and discovered Nonlinear Mixed-Effects Modeling that may be worth exploring. (I do not have SimBiology, and so have no experience with it.)

If you want more details on it, Contact Support and ask.

.

[Sepi Sharif]

Thanks I am on it now. Appreciate it

[Star Strider]

As always, my pleasure!

I see that Arthur Goldsipe has responded, so I will leave you in his capable hands.

.

Answer 2

[John D'Errico]

There can sometimes be problems when you try to apply a nonlinear estimation scheme on top of an ODE solver. The nonlinear fitting tool presumes that the objective function is a differentiable function of the parameters. But is it? Remember that an ODE solver is a tool designed to work with a tolerance. Change the parameters slightly, and you will get a subtly different result, but only to within a tolerance. Effectively, your objective function is not a differentiable function of the parameters, because there is an arbitrary and unpredictable fuzz on top of your objective.

How is this different when you use an analytical solution to the problem? Now the objective will generally be a smooth, differentiable function of the parameters.

Effectively, this CAN be a problem anytime you layer one numerical solver on top of another. One thing to try is to crank down on the tolerance of the internal solver. But will that always be sufficient? Possibly not, since if the problem is a stiff one, then the solution can vary significantly as you cary the parameters.

Another trick may be to use a stiff solver. ODE45 is a tool that is not designed to solve stiff problems. So you may gain by simply moving to a stiff solver. That is, solvers where the name ends in s, like ODE15s, or ODE23s.

Can I know that this is the reason? Not without a concerted effort to dig into exactly what you did. Is it possible that your ODE may be a little stiff? Possibly. That would certainly exaggerate any issues in the solve. And of course, there may be a bug in the code you wrote. :) Surely not! That would never happen. ;)

[Sepi Sharif]

Thank you so much for your replying.

analytical solution with nlmefit has these estimated parameters

V=8.536; ka=0.585; CL=0.140, PSI=[0.062 0 0

0 0.348 0

0 0 0.083 ]

But;

ode45 solution has

V=8.827; ka=0.772; CL=0.171 , PSI=[0.048 0 0

0 0.0000 0

0 0 0.0001]

there is a discrepancy here!

So, the analytical solution is

concentration=doseNew.*ka./(V.*(ka-kel)).*(exp(-kel.*(t))-exp(-ka.*(t)));

where, kel=CL/V;

this is a nonlinear_nonstiff problem!

Answer 3

[Arthur Goldsipe]

The reason the results from nlmefit differ significantly is because your ODE solution is not equivalent to your analyitical solution. I think you need to update file Getsolutionode.m to set minTime to 0 rather than min(t). Otherwise, your initial dose condition gets applied at min(t) instead of at time 0, and the ODE solution essentially has a time-delay when compared to your analytical solution.

In addition, as I mention in my comments on the question, I recommend changing 10-9 to 1e-9 when setting the tolerances. Otherwise, you will be using extremely loose tolerances when solving the ODEs.

Finally, I see from your example code that you are looking at a standard problem from pharmacokinetics (the "warfarin" dataset). This is exactly the sort of problem that SimBiology is designed to solve. You can build your pharmacokinetic model and run nlmefit and other analyses using SimBiology's graphical user interface or using command line functions and scripts. Here are some advantages I see to using SimBiology to solve this sort of problem:

• You can build your models graphically.
• You can select many standard PK models (including this one) from a built-in library.
• You can represent your model in a more natural way, using compartments, species, and reactions.
• You don't have to worry about the "bookkeeping" required to write the underlying differential equations.
• You can use built-in programs/tasks in the graphical user interface to perform mixed effects modeling.
• Your simulations almost always run faster. SimBiology can solve linear ODE models like this one analytically, using highly optimzed C++ code. And models that can't be solved anlytically can be solved much faster by taking advatange of its acceleration feature.

If you have a SimBiology license, I can show you how I would solve this problem in SimBiology. Just let me know whether you prefer to use the graphical user interface, the command line, or both.