NONMEM Users Network Archive

Hosted by Cognigen

Re: [NMusers] Steady-state (SS) and $DES veeeery slow

From: Paolo Denti <paolo.denti_at_uct.ac.za>
Date: Fri, 20 May 2016 17:14:48 +0200

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 differe=
ntial 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 ca=
lculating SS? Because that may have other consequences, while all I would l=
ike to do is to make the SS routine less picky, without decreasing the prec=
ision of the whole run.

Alison,
to get the values in the compartments, I have just asked NONMEM to print ou=
t the values while executing the $DES block.
NONMEM spits out a huge table here, so I filtered out only the values at ti=
me 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 aroun=
d II hours) all the values in the absorption compartments are super small o=
r also negative. If one uses transit compartments absorption, that's even w=
orse. I have pasted code and dataset below, but I will also send you a sepa=
rate 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 TEXTFIL=
E

$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 TEXTFIL=
E



$SIMULATION (123) ONLYSIM
;$ESTIMATION METHOD=1 INTER MAXEVAL=9999 PRINT=1 NOABORT NSIG=3; c=
alculation 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 $SUBROUTIN=
ES record to determine precision. ADVAN9 and ADVAN13 also use a separate a=
bsolute tolerance ATOL, which is by default set to 12, but you can modify i=
t 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 numerica=
l noise after II hours". But PREDPP never reports this amount, because SS a=
mounts are reported *after* a dose has been administered.

The way SS with multiple dosing is computed is described in the May 2015 ve=
rsion 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 i=
nitial 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 u=
sed.)

=========================

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-findin=
g subroutine from IMSL.

With multiple dosing, ZSPOW searches for a vector of compartment amounts an=
d their eta deriva- tives such that the initial and final (end of dosing in=
terval) vector is the same. ADVAN6, ADVAN8, and ADVAN9 perform a numerical =
integration of the differential equations (and in the case of ADVAN9, algeb=
raic 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 d=
ose. A search is made for Ass such that Ass - AII = 0

If SS=3, the state vector A contains user-supplied initial estimates. Oth=
erwise, the SS routines sets A=0 and makes three calls to ADVAN. Each cal=
l advances the state vector from t=0 to t=II with the appropriate dosin=
g 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 settin=
g on the $SUBROUTINES record appears to determine the precision of the stea=
dy state condition. You can change the TOL value, or, if you wish the SS p=
recision to differ from the ODE evaluation precision, you can modify SS13.f=
90, as it is open source to users, by replacing the NRD argument (which con=
tains the TOL value submitted by the $SUBROUTINES record) to all calls to Z=
SPOW1 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.Bauer_at_iconplc.com>Robert.Bauer_at_iconplc.com<mailto:Robert.Bau=
er_at_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 t=
o get some feedback.

We work on steady-state PK data, so ideally we would like to use the SS opt=
ion, 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 initial=
ising 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 t=
he SS amounts in each compartment. The model is taken to SS by repeatedly g=
iving doses every II, until the system is deemed to have reached steady-sta=
te. 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 wi=
th SS, cuz the calculation of these SS amounts implies solving the differen=
tial 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 e=
xpecting 2 dosing intervals (48 hours) to be more than enough to get to ste=
ady state, as maths says the amounts should be 99.993% there. NONMEM instea=
d used 13 doses in the example below.

The amounts are reported in the table below. After iteration 3 the amount i=
n 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 practic=
al cases - at least for preliminary runs - and in this example it would sav=
e 80% of the run time. Plus if one keeps in mind that this is a numerical s=
olver, it is not clear how "real" the wee digits are, obviously depending o=
n TOL.

The other tricky thing I found is that the amount in the absorption compart=
ment 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 rel=
ative 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 T=
OL 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 simula=
tions!



--
------------------------------------------------
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.denti_at_uct.ac.za<mailto:paolo.denti_at_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/abo=
ut/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 csirt_at_uc=
t.ac.za<mailto:csirt_at_uct.ac.za>

--
  Alison Boeckmann
  alisonboeckmann_at_fastmail.fm<mailto:alisonboeckmann_at_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.denti_at_uct.ac.za<mailto:paolo.denti_at_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/abo=
ut/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 csirt_at_uc=
t.ac.za

Received on Fri May 20 2016 - 11:14:48 EDT

The NONMEM Users Network is maintained by ICON plc. Requests to subscribe to the network should be sent to: nmusers-request@iconplc.com. Once subscribed, you may contribute to the discussion by emailing: nmusers@globomaxnm.com.