Re: bug/feature in $AES

Nick, I agree with everything you say. I'll put the new warning on the "todo" list. On Tue, 02 Oct 2007 09:31:45 +1300, "Nick Holford" <[EMAIL PROTECTED]> said: > 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 -- Alison Boeckmann [EMAIL PROTECTED]