Dear Bob and Alison, thank you for the prompt responses.
So it seems like the same TOL is used both for the precision of the differential equations and the stopping criterion of SS, correct? Bob I am not sure I feel brave/bold enough to recode my own version of NONMEM, but I can try. The question is, would this also decrease the TOL of the DE solver while calculating SS? Because that may have other consequences, while all I would like to do is to make the SS routine less picky, without decreasing the precision of the whole run. Alison, to get the values in the compartments, I have just asked NONMEM to print out the values while executing the $DES block. NONMEM spits out a huge table here, so I filtered out only the values at time 0, after each cycle of repeated dosing. It is true that at time 0 NONMEM puts the dose back in, but when looking at the tail of the curve (at around II hours) all the values in the absorption compartments are super small or also negative. If one uses transit compartments absorption, that's even worse. I have pasted code and dataset below, but I will also send you a separate email with attachments. Thanks to both for looking into this. I really think it would help to speed this up, as I have spoken to other people with the same problem. Ciao, Paolo Dataset: #ID TIME EVID AMT DV MDV SS II 1 0 1 450 . 1 1 24 1 1 0 . . 0 0 . 1 4 0 . . 0 0 . 1 11 0 . . 0 0 . 1 24 0 . . 0 0 . Model file: $PROBLEM SS investigation - normal 1st order NO LAG $INPUT ID TIME EVID AMT DV MDV SS II $DATA Sim_SS.csv IGNORE=# $SUBROUTINE ADVAN13 TRANS1 TOL=9 ; TOL=3 $MODEL NCOMPARTMENTS=2 COMP=(DEPOT,DEFDOSE) COMP=(CENTRAL) $THETA (0,8) ; 1. TVCL [L/h] $THETA (0,40) ; 2. TVV [L] $THETA (0,2) ; 3. TVKA [1/h] $THETA (0,0.15) ; 4. Sigma_prop $THETA (0,0.5) ; 5. Sigma_add [mg/L] $THETA (0,0,2) FIX ; 6. DUMMY $THETA (0,0,15) FIX ; 7. DUMMY $OMEGA 0 FIX ; 1. BSV_CL $OMEGA 0 FIX ; 2. BSV_Vd $OMEGA 0 FIX ; 3. BSV_Ka $OMEGA 0 FIX ; 4. BSV_BIO $OMEGA 0 FIX ; 5. BSV_MTT $SIGMA 1 FIX ; RESIDUAL $PK BSVCL= ETA(1) ; BSV_CL BSVV = ETA(2) ; BSV_V BSVKA= ETA(3) ; BSV_KA BSVBIO=ETA(4) ; BSV_BIO BSVMTT=ETA(5) ; BSV_MTT TVCL =THETA(1) ; Pop CL TVV =THETA(2) ; Pop V TVKA =THETA(3) ; Pop KA TVBIO=1 ; Pop BIO SIG_PROP=THETA(4) ; Sigma PROP SIG_ADD=THETA(5) ; Sigma ADD ;ALAG1 = THETA(6) CL=TVCL * EXP(BSVCL) V =TVV *EXP(BSVV) KA=TVKA*EXP(BSVKA) BIO=TVBIO*EXP(BSVBIO) K=CL/V F1 = BIO WRITE(50,2) 1, ID,TIME,A(1),A(2) ; TO OUTPUT INTERGRATION STEPS IN TEXTFILE $DES DADT(1) = -KA*A(1) DADT(2) = KA*A(1) -K*A(2) WRITE(50,2) 2, ID,T,A(1),A(2) ; TO OUTPUT INTERGRATION STEPS IN TEXTFILE $ERROR IPRED = A(2)/V IRES = DV - IPRED W=SQRT((SIG_PROP*IPRED)**2+SIG_ADD**2) ; Residual SD IWRES=IRES/W Y=IPRED+W*EPS(1) AA1=A(1) AA2=A(2) WRITE(50,2) 3, ID,TIME,A(1),A(2) ; TO OUTPUT INTERGRATION STEPS IN TEXTFILE $SIMULATION (123) ONLYSIM ;$ESTIMATION METHOD=1 INTER MAXEVAL=9999 PRINT=1 NOABORT NSIG=3; calculation method On 2016/05/19 16:23, Alison Boeckmann wrote: Dear Paolo, As Bob explained, ZSPOW uses the TOL defined by the user in the $SUBROUTINES record to determine precision. ADVAN9 and ADVAN13 also use a separate absolute tolerance ATOL, which is by default set to 12, but you can modify it with $EST … ATOL= , and see if it affects the precision of the SS evaluation. You could also try using ADVAN6/SS6. I can't understand why the SS option is taking so long with ADVAN13/SS13. The table of results you report do not correspond to the way SS is computed in PREDPP, and I can't understand where you got the numbers in the table. You say "the amount in the absorption compartment AA1 is literally numerical noise after II hours". But PREDPP never reports this amount, because SS amounts are reported *after* a dose has been administered. The way SS with multiple dosing is computed is described in the May 2015 version of Guide VI PREDPP. This is supplied with the latest alpha versions of NONMEM7.4. See Chapter V. Event Records, Section F.4 Implementation of Steady State. There are three calls to ADVAN with TIME=0 to obtain the initial estimate. After that, the root finder ZSPOW makes as many calls as it needs, but they are not necessarily sequential in TIME. Can you please send me your control stream and a fragment of the data so I can figure out what is going on? Alison PS. In case you don't have the May 2015 version of guide VI PREDPP, here is an extract: (The guide still needs to be updated for ADVAN13, but the same process is used.) ========================= When the model is defined by differential equations (SS6, SS8, SS9), PREDPP cannot use an ana- lytic solution, even if one happens to exist. Instead, the solutions are found by a numerical technique using ZSPOW, a root-finding subroutine from IMSL. With multiple dosing, ZSPOW searches for a vector of compartment amounts and their eta deriva- tives such that the initial and final (end of dosing interval) vector is the same. ADVAN6, ADVAN8, and ADVAN9 perform a numerical integration of the differential equations (and in the case of ADVAN9, algebraic equations). The ADVAN routine is used to advance (integrate) the state vector A from time 0 to time II with the appropriate dose pattern. Let AII represent the state vector after advance to time II with the appropriate dose. A search is made for Ass such that Ass - AII = 0 If SS=3, the state vector A contains user-supplied initial estimates. Otherwise, the SS routines sets A=0 and makes three calls to ADVAN. Each call advances the state vector from t=0 to t=II with the appropriate dosing pattern. The result is the initial estimate. Once the SS solution is found, the SS routine adds the final bolus dose to the state vector, or starts the final infusion. =================== Paolo:
From what I can tell in scanning the code of ..\pr\SS13.f90, the TOL setting on the $SUBROUTINES record appears to determine the precision of the steady state condition. You can change the TOL value, or, if you wish the SS precision to differ from the ODE evaluation precision, you can modify SS13.f90, as it is open source to users, by replacing the NRD argument (which contains the TOL value submitted by the $SUBROUTINES record) to all calls to ZSPOW1 with an alternative value of your choosing. Recompile SS13.f90, then execute NONMEM with nmfe73 script as usual.
Robert J. Bauer, Ph.D. Vice President, Pharmacometrics R&D ICON Early Phase 820 W. Diamond Avenue Suite 100 Gaithersburg, MD 20878 Office: (215) 616-6428 Mobile: (925) 286-0769 <mailto:robert.ba...@iconplc.com>robert.ba...@iconplc.com<mailto:robert.ba...@iconplc.com> <http://www.iconplc.com/>www.iconplc.com<http://www.iconplc.com> On Wed, May 18, 2016, at 03:15 AM, Paolo Denti wrote: Dear NMUsers (and developers), I am trying to speed up our run times when using SS in ADVAN13 models (user defined differential equations), and I would like to share some thoughts to get some feedback. We work on steady-state PK data, so ideally we would like to use the SS option, but the run times whenever we use ADVAN13 become unfeasibly long (even 20 times more), and this even for models that should reach steady state in 2-3 doses. As an alternative we end up using ad-hoc "patch-up" solutions, like initialising compartments, or just adding 4-5 doses "manually" in the dataset, but this is a bit tedious/tricky. I am writing to see if there is a way to speed up the SS feature. I decided to look into what was going on by asking NONMEM to simulate a SS PK profile and print out all the temporary iterations. It seems like NONMEM first opens a separate "SS session" used to work out the SS amounts in each compartment. The model is taken to SS by repeatedly giving doses every II, until the system is deemed to have reached steady-state. At this point, NONMEM goes back to the "main session", initialises all compartments to the amounts found in the "SS session", gives the final dose, and moves on with the analysis. Nothing surprising here, and it's understandable that things take longer with SS, cuz the calculation of these SS amounts implies solving the differential equations for all the extra time of the "SS session". What I found odd though, is that the number of doses that NONMEM uses to reach steady-state seems to me much higher than needed. In my test I used a simple 1-cmpt KA model coded in ADVAN13 and a drug with a half-life of 3.5 hours, so I was expecting 2 dosing intervals (48 hours) to be more than enough to get to steady state, as maths says the amounts should be 99.993% there. NONMEM instead used 13 doses in the example below. The amounts are reported in the table below. After iteration 3 the amount in AA2 just changes by tiny values, arguably comparable to numerical noise. #DOSE T TOT_TIME AA2 0 0 0 0 1 0 24 4.114874 2 0 48 4.148738 3 0 72 4.149017 4 0 96 4.149017 5 0 120 4.149017 6 0 144 4.149019 7 0 168 4.149019 8 0 192 4.149019 9 0 216 4.149019 10 0 240 4.149019 11 0 264 4.149019 12 0 288 4.149019 13 0 312 4.149019 Does anyone know what stopping criterion NONMEM uses to call it SS and move on? Is there a way to relax it? I think 0.1% would be fine in most practical cases - at least for preliminary runs - and in this example it would save 80% of the run time. Plus if one keeps in mind that this is a numerical solver, it is not clear how "real" the wee digits are, obviously depending on TOL. The other tricky thing I found is that the amount in the absorption compartment AA1 is literally numerical noise after II hours, meaning that for some entries it is even negative (e.g. -2.08038E-15) and jumps between positive and negative. Could it be that this is what is tricking NONMEM's stopping criterion to detect SS? Maybe if the stopping criterion only asks for a relative change <TOL, it will struggle to achieve that with values that change sign. Any settings that may help speed things up? Something like tinkering with TOL and ATOL? Or maybe is it possible to set a maximum to the number of doses that NONMEM tests to reach steady-state? Sorry for the nerdy topic, and thanks for any advice! Ciao, Paolo PS And big thanks to our PhD student Maxwell who ran all the tedious simulations! -- ------------------------------------------------ Paolo Denti, PhD Pharmacometrics Group Division of Clinical Pharmacology Department of Medicine University of Cape Town K45 Old Main Building Groote Schuur Hospital Observatory, Cape Town 7925 South Africa phone: +27 21 404 7719 fax: +27 21 448 1989 email: paolo.de...@uct.ac.za<mailto:paolo.de...@uct.ac.za> ------------------------------------------------ Disclaimer - University of Cape Town This e-mail is subject to UCT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 9111. If this e-mail is not related to the business of UCT, it is sent by the sender in an individual capacity. Please report security incidents or abuse via cs...@uct.ac.za<mailto:cs...@uct.ac.za> -- Alison Boeckmann alisonboeckm...@fastmail.fm<mailto:alisonboeckm...@fastmail.fm> -- ------------------------------------------------ Paolo Denti, PhD Pharmacometrics Group Division of Clinical Pharmacology Department of Medicine University of Cape Town K45 Old Main Building Groote Schuur Hospital Observatory, Cape Town 7925 South Africa phone: +27 21 404 7719 fax: +27 21 448 1989 email: paolo.de...@uct.ac.za<mailto:paolo.de...@uct.ac.za> ------------------------------------------------ Disclaimer - University of Cape Town This e-mail is subject to UCT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 9111. If this e-mail is not related to the business of UCT, it is sent by the sender in an individual capacity. Please report security incidents or abuse via cs...@uct.ac.za
