Re: bug/feature in $AES

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 0 FIX ;PPV_RGLU 0 FIX ;PPV_RINS $SIGMA 0.1 ;G_EXP_RUV 1 ;G_ADD_RUV MMOL/L 0.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 0 0 1 3.4096 71.756 1 100 1 0 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