NONMEM ODE solver

Dear NMusers, 1. I have some difficulties to find information on the ODE solver implemented in NONMEM. Can anyone tell me about the numerical integrator(s) implemented when using ADVAN6 (recommended for non-stiff) and ADVAN8 (stiff problems)? In particular, is the integrator able to switch automatically between stiff and non-stiff resolution techniques? 2. Then, still related to the ODE solver, I was also interested in a message I¹ve picked up from the archive (from Mark Sale, posted January 2005, see below), for which I couldn¹t find any answer from the group. Additional information on these items would be greatly appreciated.

Benjamin, According the comments in the code in version VI (imagine that - comments in code) ADVAN6 (DVERK1.for) is originally from IMSL, and is RUNGE KUTTA-VERNER FIFTH AND SIXTH ORDER METHOD. ADVAN8, (DGEAR1.for) is VARIABLE ORDER ADAMS PREDICTOR CORRECTOR METHOD OR GEARS METHOD and ADVAN9 is: the Livermore solver, from the ODEPACK. and, no, I'm pretty sure there is no switching between algorithms. (does any software do this??) BTW, no one ever answered my question from 2 years ago. Mark Sale Next Level Solutions, LLC www.NextLevelSolns.com

Mark, The Adams and Gear methods are different methods. The ODEPACK implementation can switch between these two algorithms. But whether or not NONMEM uses this feature is unclear. The switching is controlled by an argument to the ODEPACK solver and I have not been able to follow the NONMEM code (surprise, surprise). Nick Mark Sale - Next Level Solutions wrote: > > Benjamin, > According the comments in the code in version VI (imagine that - > comments in code) > ADVAN6 (DVERK1.for) is originally from IMSL, and is RUNGE KUTTA-VERNER > FIFTH AND SIXTH ORDER METHOD. > ADVAN8, (DGEAR1.for) is VARIABLE ORDER ADAMS PREDICTOR CORRECTOR METHOD > OR GEARS METHOD > and > ADVAN9 is: the Livermore solver, from the ODEPACK. > and, no, I'm pretty sure there is no switching between algorithms. (does > any software do this??) > > BTW, no one ever answered my question from 2 years ago. > > Mark Sale > Next Level Solutions, LLC > www.NextLevelSolns.com > > > -------- Original Message -------- > > Subject: [NMusers] NONMEM ODE solver > > From: Benjamin Ribba <[EMAIL PROTECTED]> > > Date: Sun, May 27, 2007 12:53 pm > > To: <[email protected]> > > > > Dear NMusers, > > > > 1. I have some difficulties to find information on the ODE solver > > implemented in NONMEM. > > Can anyone tell me about the numerical integrator(s) implemented when > > using > > ADVAN6 (recommended for non-stiff) and ADVAN8 (stiff problems)? In > > particular, is the integrator able to switch automatically between > > stiff and > > non-stiff resolution techniques? > > > > 2. Then, still related to the ODE solver, I was also interested in a > > message > > I¹ve picked up from the archive (from Mark Sale, posted January 2005, > > see > > below), for which I couldn¹t find any answer from the group. Additional > > information on these items would be greatly appreciated. > > _______________________________________________________ > > > > From: [EMAIL PROTECTED] > > Subject: [NMusers] questions on ODE solver errors > > Date: Wed, January 26, 2005 4:40 pm > > > > If the NOABORT option is used for the ODE solvers, and the > > ODE solver has "numerical difficulties", giving the message: > > > > NUMERICAL DIFFICULTIES WITH INTEGRATION ROUTINE. > > MAXIMUM NO. OF EVALUATIONS OF DIFFERENTIAL EQUATIONS, 100000, EXCEEDED. > > 0END OF PRED EXIT CODE MESSAGES FROM THE ESTIMATION STEP > > > > But, NONMEM continues to plod along. > > Does anyone know what NONMEM does instead? I assume that the OBJ > > contribution for that subject cannot be calculated, so, the gradient > > for the dimension(s) requiring that OBJ cannot be calculated. Seems > > like an error to simply ignore the OBJ contribution from that person > > (resulting in a lower OBJ for that function evaluation, and a bad > > gradient and Hessian). Is that function evaluation discarded, and > > no gradient is generated for that dimension/iteration? But I don't > > see a 0 gradient anywhere in the output (I know, in NONMEM V, it > > isn't as simple as one function evaluation per dimension (+1)) > > > > Thanks > > > > > > Mark Sale M.D. > > Global Director, Research Modeling and Simulation > > GlaxoSmithKline > > 919-483-1808 > > Mobile > > 919-522-6668 > > _______________________________________________________ > > > > > > Thanks for your help in this matter. > > > > Benjamin Ribba > > > > Therapeutic Targeting in Oncology (EA3738) > > University Claude Bernard Lyon 1 > > Faculty of Medicine Lyon-Sud > > 165, chemin du grand revoyet > > 69921 Oullins cedex > > > > Tel: +33 4 26 23 59 57 / +33 6 66 06 15 21 -- Nick Holford, Dept Pharmacology & Clinical Pharmacology University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand email:[EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:373-7556 http://www.health.auckland.ac.nz/pharmacology/staff/nholford/

Dear Dr Sale, ADAPT (including Adapt II (Rel. 4), Adapt V, and S-Adapt) uses the LSODA (Livermore Solver for Ordinary Differential equations with Automatic method switching for stiff and nonstiff problems). According to the Adapt user guide: "The ability of LSODA to detect stiffness and switch between Adam's and Gear's methods, together with its capability to control the step size and order of the selected method, makes this program [ADAPT] the most robust we know of for this application." I have not yet run a head to head comparison of the stability of the ODE solver in Adapt and NONMEM on the same project myself. However, some of my colleagues told me that they found the LSODA in Adapt very robust. I would be happy, if someone could share his/her opinion on the robustness of ODE solvers e.g. implemented in PKPD programs. Best regards Juergen ----------------------------------------------- Juergen Bulitta, PhD, Post-doctoral Fellow Pharmacometrics, University at Buffalo, NY, USA Phone: +1 716 645 2855 ext. 281, [EMAIL PROTECTED]

Hi Juergen The LSODA solver is also implemented in the odesolve package in R and used together with nlme for mixed-effects modeling in nlmeODE. I don't think you will find many PKPD programs with more well-tuned ODE solvers for mixed-effects modeling compared to NONMEM. You might also want to take a look at programs in other scientific areas besides PKPD since I've heard there has been quite a bit of improvements in the ODE solver field especially tuned for estimation problems. Best Chris ________________________________ Fra: [EMAIL PROTECTED] på vegne af Jurgen Bulitta Sendt: sø 27-05-2007 18:19 Til: [email protected] Emne: Re: [NMusers] NONMEM ODE solver Dear Dr Sale, ADAPT (including Adapt II (Rel. 4), Adapt V, and S-Adapt) uses the LSODA (Livermore Solver for Ordinary Differential equations with Automatic method switching for stiff and nonstiff problems). According to the Adapt user guide: "The ability of LSODA to detect stiffness and switch between Adam's and Gear's methods, together with its capability to control the step size and order of the selected method, makes this program [ADAPT] the most robust we know of for this application." I have not yet run a head to head comparison of the stability of the ODE solver in Adapt and NONMEM on the same project myself. However, some of my colleagues told me that they found the LSODA in Adapt very robust. I would be happy, if someone could share his/her opinion on the robustness of ODE solvers e.g. implemented in PKPD programs. Best regards Juergen ----------------------------------------------- Juergen Bulitta, PhD, Post-doctoral Fellow Pharmacometrics, University at Buffalo, NY, USA Phone: +1 716 645 2855 ext. 281, [EMAIL PROTECTED]

My comments are attached. On Sun, 27 May 2007 18:53:21 +0200, "Benjamin Ribba" <[EMAIL PROTECTED]> said: > Dear NMusers, > > 1. I have some difficulties to find information on the ODE solver > implemented in NONMEM. > Can anyone tell me about the numerical integrator(s) implemented when > using > ADVAN6 (recommended for non-stiff) and ADVAN8 (stiff problems)? In > particular, is the integrator able to switch automatically between stiff > and > non-stiff resolution techniques? > > 2. Then, still related to the ODE solver, I was also interested in a > message > I¹ve picked up from the archive (from Mark Sale, posted January 2005, see > below), for which I couldn¹t find any answer from the group. Additional > information on these items would be greatly appreciated. > _______________________________________________________ > > From: [EMAIL PROTECTED] > Subject: [NMusers] questions on ODE solver errors > Date: Wed, January 26, 2005 4:40 pm > > If the NOABORT option is used for the ODE solvers, and the > ODE solver has "numerical difficulties", giving the message: > > NUMERICAL DIFFICULTIES WITH INTEGRATION ROUTINE. > MAXIMUM NO. OF EVALUATIONS OF DIFFERENTIAL EQUATIONS, 100000, EXCEEDED. > 0END OF PRED EXIT CODE MESSAGES FROM THE ESTIMATION STEP > > But, NONMEM continues to plod along. > Does anyone know what NONMEM does instead? I assume that the OBJ > contribution for that subject cannot be calculated, so, the gradient > for the dimension(s) requiring that OBJ cannot be calculated. Seems > like an error to simply ignore the OBJ contribution from that person > (resulting in a lower OBJ for that function evaluation, and a bad > gradient and Hessian). Is that function evaluation discarded, and > no gradient is generated for that dimension/iteration? But I don't > see a 0 gradient anywhere in the output (I know, in NONMEM V, it > isn't as simple as one function evaluation per dimension (+1)) > > Thanks > > > Mark Sale M.D. > Global Director, Research Modeling and Simulation > GlaxoSmithKline > 919-483-1808 > Mobile > 919-522-6668 > _______________________________________________________ > > > Thanks for your help in this matter. > > Benjamin Ribba > > Therapeutic Targeting in Oncology (EA3738) > University Claude Bernard Lyon 1 > Faculty of Medicine Lyon-Sud > 165, chemin du grand revoyet > 69921 Oullins cedex > > Tel: +33 4 26 23 59 57 / +33 6 66 06 15 21 > > > -- Alison Boeckmann [EMAIL PROTECTED] All emails to this topic have been appropriate. I can add a few comments. (1) By default, with ADVAN9, MF=21 (if there are no etas) or MF=22 (if there are etas). The significance of these values can be read in the comments in the ODEPACK routine LSODI1. By default, with ADVAN8, METH=2 and MITER=2. Again, the significance is documented in the comments in the IMSL routine DGEAR. Stuart Beal made a choice of the METH/MITER values a long time ago based on these comments and neither of us experimented with other choices. Routines from IMSL and ODEPACK are to some extent black boxes to me. Stuart may have studied their code when he designed the ADVAN routines. I never did. (2) The user has some control over the methods used by LSODI (with ADVAN9) and DGEAR (with ADVAN8). See the help/html entries for PRCOMG. I am almost 100% certain that no one has ever actually attempted to use non-default values of METH or MITER. (3) Why were there no comments in these routines till NONMEM VI? Actually, the commented versions were distributed with NONMEM III, in 1992. The organization of the distribution media (mag tape and floppies) was very different then from what it is today, there was less code, and there was room for commented code. Starting with NONMEM IV, the organization was changed (in part, to accommodate the SETUP and SETUP.BAT utilities, which were new), and the code grew. In order to fit the code on the floppies, we (Stuart and I) decided to strip the comments. Was this a good decision? I don't know. With NONMEM VI and the use of CDROM as the (sole) distribution medium, there is plenty of room for comments, and Tom Ludden and I decided to include them, if only so that you users could attempt to answer questions such as the current topic without the help of the original developers. (4) Re: The PREDPP error message NUMERICAL DIFFICULTIES WITH INTEGRATION ROUTINE. MAXIMUM NO. OF EVALUATIONS OF DIFFERENTIAL EQUATIONS, 100000, EXCEEDED. This message was not originally part of ADVANs 6 or 8 or 9. However, someone sent us at least one problem in which the integrating routine drove the step size to nearly 0, and the ADVAN routine got into a seemingly infinite loop. With NONMEM IV, the ADVANs were changed to count the number of calls to DES, and abort if it exceeded the above arbitrarily chosen limit. IMAX was added to PRCOMG to allow the user to supply a larger limit using verbatim code, as shown in the example in the help entry. I am certain that this feature has been used, to good effect. The value of 100000 is probably too small, especially with modern fast computers. (5) What does NONMEM do if the above error return occurs from PREDPP? If you do not use NOABORT option, then NONMEM quits. If NOABORT is present, then NONMEM attempts theta-recovery or eta-recovery. (see Guide V, chapter 12). As I recall, Stuart told me that NONMEM in effect puts a "spike", i.e., "infinity", in the objective function during the search procedure, for the eta and/or theta values that gave rise to the error condition. NONMEM tries other values for eta (for this subject, if doing conditional estimation), or other values for theta (for all subjects, if not conditional estimation). If it cannot find any values of eta and/or theta for which it can avoid the error condition, then it quits with a message from NONMEM. I think (but am not certain) that this would be a message such as MINIMIZATION TERMINATED DUE TO PROXIMITY OF NEXT ITERATION EST. TO A VALUE AT WHICH THE OBJ. FUNC. IS INFINITE If with NOABORT you see only a few error messages in file PRDERR, early in the estimation step, then you need not worry - a few unrealistic values of theta and/or eta were avoided. But if there are many, and they are still showing up near the end of the search, then NONMEM may have been hampered in the search because it was constrained to avoid what might have been an important area of search space. What Mark conjectured ("I assume that the OBJ contribution for that subject cannot be calculated, so, the gradient for the dimension(s) requiring that OBJ cannot be calculated") is not true. NONMEM always includes every subject's and every observation's contribution to the objective function. If it cannot, then it aborts.

Alison, Thanks for this, (the answer to my question from 2 years ago). This will help greatly with some automated diagnostics/corrections in work on automated model selection, especially with my favorite error message MINIMIZATION TERMINATED DUE TO PROXIMITY OF NEXT ITERATION EST. TO A VALUE AT WHICH THE OBJ. FUNC. which I never understood how that could happen (unless you have a proportional error and prediction of 0). The "spike", is an interesting solution to be problem (not that I have a better suggestion, since, in that situation you really can't calculate a meaningful OBJ). I'd like to ask some of the people out there who know about numerical method if there are other options for recovering from a THETA vector for which you cannot the ODE solution, in estimation. Perhaps an estimate of the OBJ based just on the first and second derivative (with some penalty, to keep it out of that range again)? Mark Mark Sale MD Next Level Solutions, LLC www.NextLevelSolns.com

Mark, just wondering: Why does it make a difference whether NONMEM runs into an error condition from PRED with ODE's versus any other method of computing the prediction F? There are error messages and error returns from other ADVAN routines. Also, a user $PRED has the option of using the EXIT statement if the values of THETA/ETA are such that it cannot compute the prediction. NONMEM does the same thing in all cases: back away from the problematic THETA/ETA because of the "spike", and search in a different direction for a better THETA/ETA. BTW, if PREDPP cannot return F, then it certainly cannot return the first and second eta derivatives of F, G and H. Sometimes a given error condition is raised in PREDPP not because F cannot be computed, but during the computation of the eta derivatives. There is no distinction made to NONMEM as to which was the problem. With the ODE's, for example, an augmented state vector (A and dA/deta and if appropriate d2A/detadeta) is integrated all together. (With ADVAN5 and 7, the ADVAN does know if the trouble was with a given eta and the text of the message itself indicates this, though nothing is said to NONMEM about it.) On Mon, 28 May 2007 13:04:22 -0700, "Mark Sale - Next Level Solutions" <[EMAIL PROTECTED]> said: > Alison, Thanks for this, (the answer to my question from 2 years > ago). This will help greatly with some automated > diagnostics/corrections in work on automated model selection, > especially with my favorite error message > > MINIMIZATION TERMINATED DUE TO PROXIMITY OF NEXT ITERATION EST. TO A > VALUE AT WHICH THE OBJ. FUNC. which I never understood how that could > happen (unless you have a proportional error and prediction of 0). > > The "spike", is an interesting solution to be problem (not that I > have a better suggestion, since, in that situation you really can't > calculate a meaningful OBJ). I'd like to ask some of the people out > there who know about numerical method if there are other options for > recovering from a THETA vector for which you cannot the ODE solution, > in estimation. Perhaps an estimate of the OBJ based just on the first > and second derivative (with some penalty, to keep it out of that > range again)? > > > > Mark > > Mark Sale MD Next Level Solutions, LLC www.NextLevelSolns.com > > > > -------- Original Message -------- Subject: Re: [NMusers] NONMEM ODE > > solver From: "Alison Boeckmann" <[EMAIL PROTECTED]> Date: > > Mon, May 28, 2007 1:01 pm To: "Benjamin Ribba" <[EMAIL PROTECTED] > > lyon1.fr>, [email protected] > > > > My comments are attached. > > > > > > On Sun, 27 May 2007 18:53:21 +0200, "Benjamin Ribba" <[EMAIL PROTECTED] > > lyon1.fr> said: > > > Dear NMusers, > > > > > > 1. I have some difficulties to find information on the ODE solver > > > implemented in NONMEM. Can anyone tell me about the numerical > > > integrator(s) implemented when using ADVAN6 (recommended for > > > non-stiff) and ADVAN8 (stiff problems)? In particular, is the > > > integrator able to switch automatically between > > stiff > > > and non-stiff resolution techniques? > > > > > > 2. Then, still related to the ODE solver, I was also interested in > > > a message I?ve picked up from the archive (from Mark Sale, > > > posted January > > 2005, see > > > below), for which I couldn?t find any answer from the group. > > Additional > > > information on these items would be greatly appreciated. > > > _______________________________________________________ > > > > > > From: [EMAIL PROTECTED] Subject: [NMusers] questions on ODE > > > solver errors Date: Wed, January 26, 2005 4:40 pm > > > > > > If the NOABORT option is used for the ODE solvers, and the ODE > > > solver has "numerical difficulties", giving the message: > > > > > > NUMERICAL DIFFICULTIES WITH INTEGRATION ROUTINE. MAXIMUM NO. OF > > > EVALUATIONS OF DIFFERENTIAL EQUATIONS, 100000, > > EXCEEDED. > > > 0END OF PRED EXIT CODE MESSAGES FROM THE ESTIMATION STEP > > > > > > But, NONMEM continues to plod along. Does anyone know what NONMEM > > > does instead? I assume that the OBJ contribution for that subject > > > cannot be calculated, so, the gradient for the dimension(s) > > > requiring that OBJ cannot be calculated. Seems like an error to > > > simply ignore the OBJ contribution from that person (resulting in > > > a lower OBJ for that function evaluation, and a bad gradient and > > > Hessian). Is that function evaluation discarded, and no gradient > > > is generated for that dimension/iteration? But I don't see a 0 > > > gradient anywhere in the output (I know, in NONMEM V, it isn't as > > > simple as one function evaluation per dimension (+1)) > > > > > > Thanks > > > > > > > > > Mark Sale M.D. Global Director, Research Modeling and Simulation > > > GlaxoSmithKline 919-483-1808 Mobile 919-522-6668 > > > _______________________________________________________ > > > > > > > > > Thanks for your help in this matter. > > > > > > Benjamin Ribba > > > > > > Therapeutic Targeting in Oncology (EA3738) University Claude > > > Bernard Lyon 1 Faculty of Medicine Lyon-Sud 165, chemin du grand > > > revoyet 69921 Oullins cedex > > > > > > Tel: +33 4 26 23 59 57 / +33 6 66 06 15 21 > > > > > > > > > > > -- > > Alison Boeckmann [EMAIL PROTECTED] > -- Alison Boeckmann [EMAIL PROTECTED]

Alison, I guess it doesn't make any difference why PRED returns an error, just happens so much more frequently with ODEs. When I was working on parallel NONMEM code, I walked though the THETA recovery (taking half the previous step size) because I had to make sure all the variable were handled correctly. The interesting thing is that, with the ODE solvers, and this estimation algorithm, the surface is so discontinuous (a spike to infinity, could be right next to a value for THETA that just barely completes the integration, but returns a reasonable value. But, I'm way out of my field here. WRT using the first and second derivative, I meant the derivative from the previous iteration, which, I think could give you and estimate of what the OBJ would be at the problem values of THETA, then perhaps add some (large) penalty to hopefully prevent the algorithm from going there again (at least right away). Mark Mark Sale MD Next Level Solutions, LLC www.NextLevelSolns.com