Re: bug/feature in $AES

Alison, I deliberately used the term bug/feature because it is often a question of perspective whether it is a bug or a feature based on one's expectations. The original design of NONMEM was based on a limited perspective of solving PK problems with the assumption that all compartments would start at zero. But the world is more interesting than that. When using $AES I had the expectation that the equilibrium state would be computed for the first event. I did not think that the $AES would wait for an arbitrary time (until after the first event) to 'wake up' and start to work. As you said in an earlier email the function of $AES is not necessarily dependent on time. Because the $AES compartments have no necessary connection with the $DES compartments there is no reason to force the $DES assumption that all compartments start at zero. Thus $AES should give the equilibrium answer for all events. For me this is a bug. For you and Stuart this seems to be a feature. The addition of a warning when $AES is used something like: WARNING: Equilibrium compartments are not computed for the first event record would at least warn people not to have expectations of reasonable behaviour. The work around mechanism of putting AMT and SS items in the data for the only purpose of getting the model to start as desired is an example of poor design. This failure to separate the model from the data was one of the points discussed at a NONMEM design thinktank (including Lewis and Tom) held at Globomax in December 2001. Allowing $DES compartment initialization without forcing the user to modify the data to implement a model choice has fortunately been addressed in NMVI. A better user interface for $AES would allow the user to control the model from the control stream instead of putting magic values in the data e.g. in $AESINITIAL a value of INIT=2 could mean $AES should be computed prior to the first event so that the equilibrium compartments truly reflect the initial equilibrium state. Thanks for the clue about ZSPOW1. It seems that this is hard to do even with the NMVI abbreviated function mechanism. But given that $AES can be forced to do the work then anybody who wants a general root finder can use ADVAN9. Best wishes, Nick Alison Boeckmann wrote: > > Nick, > > There is no bug. You yourself titled this thread "bug/feature". PREDPP > is designed to start every model with all compartments set to 0 prior to > the first dose. The fact that ADVAN9 is used does not and should not > imply that a system is automatically at steady state (equilibrium). > There are models in which a bolus dose is necessary before the D.E.'s > can be evaluated, e.g., compartment amounts are in the denominator of > the D.E.'s. > > Even if NM-TRAN could somehow infer from the D.E.'s that there is > endogenous drug, hence the system "must" be at SS, there is no mechanism > within PREDPP to make this happen with the first event record other than > via SS dose. If you don't want to include SS and RATE data items, it is > always possible to let a system advance from time 0 to a sufficiently > large time to initialize all the compartments with equilibrium values. > It just takes more computer time, hence the change to PREDPP to allow it > to be done with an SS dose record. (I recall that it was one of the > fellows who ran into this issue at a time when the computers at c255 > were very slow compared with current hardware.) > > NM-TRAN does try to warn in cases when surprising things happen, e.g., > when WRITE or PRINT is used in $DES: > > (WARNING 48) DES-DEFINED ITEMS ARE COMPUTED ONLY WHEN EVENT TIME > INCREASES. E.G., DISPLAYED VALUES ASSOCIATED WITH THE FIRST EVENT > RECORD OF AN INDIVIDUAL RECORD ARE COMPUTED WITH (THE LAST ADVANCE TO) > AN EVENT TIME OF THE PRIOR INDIVIDUAL RECORD. > > Would you like to see a similar warning, e.g., > > WARNING: ALL COMPARTMENT AMOUNTS ARE ZERO UNTIL AFTER THE FIRST ADVANCE > > (give this warning whenever ADVAN6 or 8 or 9 is used without a dose > event as the first event record ... or some such criterion ???) > > As for your comment about AMT=0/RATE=0/SS>0 being "a great find buried > in the documentation," I suspect that there are a number of nmusers who > have been quietly scratching their heads about this, thinking "but I've > been doing this for years. Why does he not know about it, and how could > she have forgotten?". What I will do is change Guide VIII and on-line > help so that "help endogenous" lists "ss dose" as a choice. This is a > mistake in the help system. (BTW: there are inevitably other errors and > omissions in the help system. Anyone who notices an error in the > documentation should please bring this to the attention of > [EMAIL PROTECTED]) > > You say "NONMEM must be doing some kind of computation to get the > initial steady state solution. This must require extra computation > compared to an explicit calculation in the code of initial conditions." > Indeed, with the analytic models SS1, SS2, etc., the steady state > solutions are computed explicitly. With ADVAN6 and 9, the SS condition > is computed by PREDPP using an IMSL routine, ZSPOW1, which is part of > the PREDPP package. As you surmise, it is a numerical root finder. With > steady state infusion, in SS6: > > C SS IS SOLUTION A OF: 0=DADT(A)+R > > Here, A is the state vector (compartment amounts) augmented by all the > eta derivatives dA/deta for each compartment A and eta. R is the vector > of infusions into each compartment (and their eta derivatives). In the > special case we are talking about, R is 0. DADT is computed by DES. > With SS9, the situation is a little more complicated, but is similar. > > I can think of no way to allow access to ZSPOW directly from > abbreviated code. It is a very complicated routine to use. But with > NONMEM VI there a new feature of abbreviated code, abbreviated > functions FUNCA, FUNCB and FUNCC. I suppose a person could study the > comments in ZSPOW1, study how it is used by PREDPP, and do something > similar in an abbreviated function. But the eta derivatives are always > going to be the hardest part. > > I hope this answers all your concerns. > > On Sat, 29 Sep 2007 09:50:37 +1200, "Nick Holford" > <[EMAIL PROTECTED]> said: > > Alison, Leonid, > > > > What a great find buried in the documentation. Well done Leonid for > > finding it and Alison for explaining how it can be done even more > > simply. Thank you! > > > > The AMT=0/RATE=0/SS>0 method not only solves the problem of how to get > > initial conditions for complex systems but can be used for any $DES > > problem that requires an initial condition. Thus it makes it > > unnecessary in NMV to add a special dosing record for each compartment > > in order to use the bioavailability trick. And for both NMV and NMVI > > it means that users no longer have to calculate the initial conditions > > in the code. > > > > For completeness I show below the full control stream and sample data > > file and output which shows how to use the AMT=0/RATE=0/SS>0 method > > for initialising DE systems. In this case the initial conditions do > > not have an algebraic solution but the method would work for any > > initial conditions. The glucose and insulin parameters are more > > realistic (Silber et al. 2007) than in my earlier code example. > > > > I see there are two prices that have to be paid for the > > AMT=0/RATE=0/SS>0 mechanism. > > > > 1. The data set must be modified to include RATE and SS data items and > > an extra record must be inserted for each subject at TIME=0 with > > AMT=0, RATE=0 and SS>0. Note that MDV and EVID are not required > > when using NM-TRAN. > > > > 2. NONMEM must be doing some kind of computation to get the initial > > steady state solution. This must require extra computation compared > > to an explicit calculation in the code of initial conditions. > > > > I wonder if Alison can explain what method is used to calculate the SS > > conditions for a DE system. I assume it uses a numerical root finder. > > Would be possible to make this root finder more generally available > > anywhere in abbreviated code? There are other modelling problems that > > could benefit from having a general purpose root finder available. > > > > Alison asks if there is any need to fix the bug in $AES that means it > > does not calculate the correct equilibrium solution on the first event > > record (if the AMT=0/RATE=0/SS>0 trick is not used). I dont think it > > is a good idea to leave known bugs to lurk around to trap unwary > > users. Presumably anyone who has used $AES in the past will have been > > bitten if their model relied upon having NONMEM find the correct > > equilibrium solution at the start. > > > > Best wishes, > > > > Nick > > > > > > --- sample control stream testssii.ctl ---- $PROB GLUCOSE AND INSULIN > > INITIAL VALUE $DATA testssii.csv $INPUT ID TIME CMT AMT RATE SS DV > > $SIM (200070927) ONLYSIM NSUB=1 ;Silber HE, Jauslin PM, Frey N, > > Gieschke R, Simonsson US, Karlsson MO. ;An integrated model for > > glucose and insulin regulation in healthy volunteers ;and type 2 > > diabetic patients following intravenous glucose provocations. ;J Clin > > Pharmacol. 2007 Sep;47(9):1159-71. $THETA 40 ; POP_RGLU MMOL/H/70KG 40 > > ; POP_VGLU L/70KG 5 ; POP_CLGLU L/H/70KG 1 ; POP_EMXGLU 10 ; > > POP_C50GLU MMOL/L 5 ; POP_HILGLU > > > > 5000 ; POP_RINS PMOL/H/70KG 5 ; POP_VINS L/70KG 70 ; POP_CLINS > > L/H/70KG 2 ; POP_EMXINS 50 ; POP_C50INS PMOL/L 2 ; POP_HILINS > > > > $OMEGA FIX ;PPV_RGLU FIX ;PPV_RINS > > > > $SIGMA .1 ;G_EXP_RUV 1 ;G_ADD_RUV MMOL/L .1 ;I_EXP_RUV 1 > > ;I_ADD_RUV PMOL/L > > > > $SUBR ADVAN6 TOL=3 > > > > $MODEL COMP (GLUCOSE) COMP (INSULIN) > > > > $PK > > > > ; GLUCOSE RGLU=THETA(1)*EXP(ETA(1)) VGLU=THETA(2) CLGLU=THETA(3) > > EMXGLU=THETA(4) C50GLU=THETA(5) HILGLU=THETA(6) > > > > ;INSULIN RINS=THETA(7)*EXP(ETA(2)) VINS=THETA(8) CLINS=THETA(9) > > EMXINS=THETA(10) C50INS=THETA(11) HILINS=THETA(12) > > > > S1=VGLU S2=VINS > > > > $DES DGLU=A(1)/VGLU DINS=A(2)/VINS DGLUH=DGLU**HILGLU > > DGEFF=EMXGLU*DGLUH/(C50GLU**HILGLU+DGLUH) DINSH=DINS**HILINS > > DIEFF=EMXINS*DINSH/(C50INS**HILINS+DINSH) DADT(1)=RGLU - > > CLGLU*(1+DIEFF)*DGLU DADT(2)=RINS*(1+DGEFF)- CLINS*DINS > > > > $ERROR GLU=A(1)/VGLU INS=A(2)/VINS IF (CMT.EQ.1)THEN > > Y=GLU*(1+ERR(1))+ERR(2) ENDIF > > > > IF (CMT.EQ.2)THEN Y=INS*(1+ERR(3))+ERR(4) ENDIF > > > > $TABLE ID TIME CMT AMT RATE SS GLU INS NOAPPEND ONEHEADER NOPRINT > > FILE=testssii.fit --- testssii.csv ---- #ID,TIME,CMT,AMT,RATE,SS,DV > > 1,0,1,0,0,1,. 1,100,1,.,.,.,. 1,100,2,.,.,.,. --- testssii.fit ---- ID > > TIME CMT AMT RATE SS GLU INS 1 0 1 1 > > 3.4096 71.756 1 100 1 0 0 > > 3.4096 71.756 1 100 2 0 0 0 > > 3.4096 71.756 > > > > Alison Boeckmann wrote: > > > > > > Leonid, this is brilliant. Thank you!! > > > > > > In fact, I'd completely forgotten that it can be done even more > > > simply. With SS dose, you don't even need to use ADVAN9 or $AES. > > > > > > With NONMEM V, the rules for the Steady state dose record were > > > changed. From the NONMEM V Supplemental Guide: > > > > > > Documentation change. Guide VI (pre 1998 editions) states (in > > > reference to the RATE data item): "A 0 indicates that the dose is > > > an instantaneous bolus dose." There is an exception: a record with > > > AMT=0, RATE=0, SS>0, and II=0 allows steady- state amounts to be > > > computed by PREDPP based on system input given by the DES routine. > > > This is useful when the differential equations provide for > > > endogenous drug production. > > > > > > Note that the ssdose help item says the same thing. > > > > > > With your problem, use ADVAN6 with the $DES as-is. Omit $AES, > > > $AESINT and all code for GBASE, IBASE, and F1. Add to $INPUT: SS > > > RATE The record at TIME 0 should have SS=1, RATE=0, AMT=0 > > > #id,time,cmt,amount,dv,ss,rate 1,0,1,.,.,1,0 Omit doses at .001 > > > and .002. > > > > > > And that does it. > > > A(1) and A(2) have their steady state amounts from time 0. (These > > > are the amounts produced by the differential equations from > > > their endogenous terms RGLU and RINS.) > > > > > > Nick, this is so simple that I can see no reason to change ADVAN9. > > > You suggested "If this can be fixed then I would suggest that an > > > extra state be added to the $AEINITIAL INIT variable. If INIT=-1 > > > then this would mean to obtain a solution at time=0 and not bother > > > to obtain solutions at subsequent times (for computational > > > efficiency)." Can you or anyone think of any reason to do this now?? > > > > > > Thanks again, Leonid, for suggesting a steady state dose record. > > > > > -- > > Nick Holford, Dept Pharmacology & Clinical Pharmacology University of > > Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand > > [EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090 > > www.health.auckland.ac.nz/pharmacology/staff/nholford > -- > Alison Boeckmann > [EMAIL PROTECTED] -- Nick Holford, Dept Pharmacology & Clinical Pharmacology University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand [EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090 www.health.auckland.ac.nz/pharmacology/staff/nholford