bug/feature in $AES

Hi, Carl Kirkpatrick and I have been trying to use $AES to compute initial values for use with $DES. We can get $AES to produce the correct equilibrium values at times after time=0 but cannot access these values at time=0 or at time=0.001. Because the equilibrium (i.e. initial condition) values are not available at time=0 we cannot use the NMVI method for initalizing the amounts in the compartments using A_0(1)=GBASE and A_0(2)=IBASE. We have inserted dosing records using the old NMV bioavailability trick to try and initialize the compartments at time=0.002. This almost gets the correct solution at time=0.002 in both compartments. Can anyone explain why the algebraic equation results are not available at time=0 and time=0.001? Nick and Carl --- table file output ----- GLU=A(1) INS=A(2) ID TIME AMT CMT F1 F2 GBASE GLU IBASE INS 1 0 0 1 0 0 0 0 0 0 1 0 0 2 0 0 0 0 0 0 1 0.001 0 1 0 0 0 0.0009 0 0.0024 1 0.002 1 1 187.14 527.49 4.6785 4.6805 13.187 0.0049 1 0.002 1 2 187.14 527.49 4.6785 4.6805 13.187 13.192 1 100 0 1 0 0 4.6785 4.6785 13.187 13.187 1 100 0 2 0 0 4.6785 4.6785 13.187 13.187 ---- Control stream -------- $PROB GLUCOSE AND INSULIN INITIAL VALUE $DATA test.dat $INPUT ID TIME CMT AMT DV $SIM (200070927) ONLYSIM NSUB=1 $THETA 40 ; POP_RGLU MMOL/H/70KG 40 ; POP_VGLU L/70KG 4 ; POP_CLGLU L/H/70KG 1 ; POP_EMXGLU 10 ; POP_C50GLU MMOL/L 100 ; POP_RINS NMOL/H/70KG 40 ; POP_VINS L/70KG 10 ; POP_CLINS L/H/70KG 2 ; POP_EMXINS 10 ; POP_C50INS NMOL/L $OMEGA 0 FIX ;PPV_RGLU 0 FIX ;PPV_RINS $SIGMA 0.1 ;G_EXP_RUV 1 ;G_ADD_RUV 0.1 ;I_EXP_RUV 1 ;I_ADD_RUV $SUBR ADVAN9 TOL=3 $MODEL COMP (GLUCOSE) COMP (INSULIN) COMP (GLU EQUILIBRIUM) COMP (INS EQUILIBRIUM) $PK "FIRST " COMMON/PRCOMG/IDUM1,IDUM2,IMAX " INTEGER IDUM1,IDUM2,IMAX " IMAX=1000000 ; GLUCOSE RGLU=THETA(1)*EXP(ETA(1)) VGLU=THETA(2) CLGLU=THETA(3) EMXGLU=THETA(4) C50GLU=THETA(5) ;INSULIN RINS=THETA(6)*EXP(ETA(2)) VINS=THETA(7) CLINS=THETA(8) EMXINS=THETA(9) C50INS=THETA(10) GBASE=A(3)/VGLU ; GLUCOSE INITIAL IBASE=A(4)/VINS ; INSULIN INITIAL IF (AMT.EQ.1) THEN F1=GBASE F2=IBASE ELSE F1=0 F2=0 ENDIF S1=VGLU S2=VINS $AESINITIAL INIT=0 A(3)=5 ; MMOL/L A(4)=10 ; NMOL/L $AES EGLU=A(3)/VGLU EINS=A(4)/VINS EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) EIEFF=EMXINS*EINS/(C50INS+EINS) E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM $DES DGLU=A(1)/VGLU DINS=A(2)/VINS DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) DIEFF=EMXINS*DINS/(C50INS+DINS) 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 AMT CMT F1 F2 GBASE GLU IBASE INS NOAPPEND ONEHEADER NOPRINT FILE=test.fit ---- test.dat -------------------- #id,time,cmt,amount,dv 1,0,1,.,. 1,0,2,.,. 1,0.001,1,.,. 1,0.002,1,1,. 1,0.002,2,1,. 1,100,1,.,. 1,100,2,.,. -- 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

Hi Nick, It would help if you prepare a cleaner example that demonstrate the problem. What exactly are you trying to achieve (A(?)=?). A(1) and A(2) are not initialized in the code below, so they are zeros as they should. Is this about A3, A4? If yes, the problem could be in communication between $AES and $PK. Is it allowed to use A() in the PK block? > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > IBASE=A(4)/VINS ; INSULIN INITIAL From help: " Right-hand quantities (of PK block can contain - LG) in assignment statement and in conditions: Data item labels specified on the $INPUT statement. THETA(n). ETA(n) (Used if the data are population.) PK-defined items that appeared earlier as left-hand quantities. I think, PK and blocks below can communicate via COMMONs but you need to take special action (via COMRES) to make them available, and they have a 1-record delay (next record knows about COM() of the previous record). I could be wrong on that, did it just once, and it was related to ERROR, not AES, but we put extra record with the same time, just to pass variable to PK Leonid Nick Holford wrote: > Hi, > [Second attempt - there was an error in the copy of the control stream in our > first attempt] > Carl Kirkpatrick and I have been trying to use $AES to compute initial values > for use with $DES. > We can get $AES to produce the correct equilibrium values at times after time=0 > but cannot access these values at time=0 or at time=0.001. Because the > equilibrium (i.e. initial condition) values are not available at time=0 we > cannot use the NMVI method for initalizing the amounts in the compartments > using A_0(1)=GBASE and A_0(2)=IBASE. > > We have inserted dosing records using the old NMV bioavailability trick to try > and initialize the compartments at time=0.002. This almost gets the correct > solution at time=0.002 in both compartments. > > Can anyone explain why the algebraic equation results are not available at > time=0 and time=0.001? > > Nick and Carl > > --- table file output ----- GLU=A(1) INS=A(2) > ID TIME AMT CMT F1 F2 GBASE GLU IBASE INS > 1 0 0 1 0 0 0 0 0 0 > 1 0 0 2 0 0 0 0 0 0 > 1 0.001 0 1 0 0 0 0.0009 0 0.0024 > 1 0.002 1 1 187.14 527.49 4.6785 4.6805 13.187 0.0049 > 1 0.002 1 2 187.14 527.49 4.6785 4.6805 13.187 13.192 > 1 100 0 1 0 0 4.6785 4.6785 13.187 13.187 > 1 100 0 2 0 0 4.6785 4.6785 13.187 13.187 > ---- Control stream -------- > $PROB GLUCOSE AND INSULIN INITIAL VALUE > $DATA test.csv > $INPUT ID TIME CMT AMT DV > $SIM (200070927) ONLYSIM NSUB=1 > > $THETA > 40 ; POP_RGLU MMOL/H/70KG > > 40 ; POP_VGLU L/70KG 4 ; POP_CLGLU L/H/70KG > > 1 ; POP_EMXGLU > 10 ; POP_C50GLU MMOL/L > 100 ; POP_RINS NMOL/H/70KG > 40 ; POP_VINS L/70KG > 10 ; POP_CLINS L/H/70KG > 2 ; POP_EMXINS > 10 ; POP_C50INS NMOL/L > > $OMEGA > 0 FIX ;PPV_RGLU > 0 FIX ;PPV_RINS > > $SIGMA > > 0.1 ;G_EXP_RUV > 1 ;G_ADD_RUV > 0.1 ;I_EXP_RUV > 1 ;I_ADD_RUV > > $SUBR ADVAN9 TOL=3 > > $MODEL > COMP (GLUCOSE) > COMP (INSULIN) > COMP (GLU EQUILIBRIUM) > COMP (INS EQUILIBRIUM) > > $PK "FIRST " COMMON/PRCOMG/IDUM1,IDUM2,IMAX " INTEGER IDUM1,IDUM2,IMAX " IMAX=1000000 > > ; GLUCOSE > RGLU=THETA(1)*EXP(ETA(1)) > > VGLU=THETA(2) CLGLU=THETA(3) > > EMXGLU=THETA(4) > C50GLU=THETA(5) > > ;INSULIN RINS=THETA(6)*EXP(ETA(2)) > > VINS=THETA(7) > CLINS=THETA(8) > EMXINS=THETA(9) > C50INS=THETA(10) > > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > IBASE=A(4)/VINS ; INSULIN INITIAL > IF (AMT.GT.0) THEN > F1=GBASE*VGLU > F2=IBASE*VINS > ELSE > F1=0 > F2=0 > ENDIF > S1=VGLU > S2=VINS > > $AESINITIAL > INIT=0 > > A(3)=5 ; MMOL/L A(4)=10 ; NMOL/L > > $AES > EGLU=A(3)/VGLU > EINS=A(4)/VINS > EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) > EIEFF=EMXINS*EINS/(C50INS+EINS) > E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM > E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM > > $DES DGLU=A(1)/VGLU DINS=A(2)/VINS DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) > > DIEFF=EMXINS*DINS/(C50INS+DINS) > 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 AMT CMT F1 F2 GBASE GLU IBASE INS NOAPPEND ONEHEADER NOPRINT FILE=test.fit > > ---- test.dat -------------------- > #id,time,cmt,amount,dv > 1,0,1,.,. > 1,0,2,.,. > 1,0.001,1,.,. > 1,0.002,1,1,. > 1,0.002,2,1,. > 1,100,1,.,. > 1,100,2,.,. > > -- > 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

I have attached a reply to Nick's email. ---------- On Thu, 27 Sep 2007 23:08:26 +1200, "Nick Holford" <[EMAIL PROTECTED]> said: > Hi, > [Second attempt - there was an error in the copy of the control stream in > our first attempt] > Carl Kirkpatrick and I have been trying to use $AES to compute initial > values for use with $DES. > We can get $AES to produce the correct equilibrium values at times after > time=0 but cannot access these values at time=0 or at time=0.001. Because > the equilibrium (i.e. initial condition) values are not available at > time=0 we cannot use the NMVI method for initalizing the amounts in the > compartments using A_0(1)=GBASE and A_0(2)=IBASE. > > We have inserted dosing records using the old NMV bioavailability trick > to try and initialize the compartments at time=0.002. This almost gets > the correct solution at time=0.002 in both compartments. > > Can anyone explain why the algebraic equation results are not available > at time=0 and time=0.001? > > Nick and Carl > > --- table file output ----- GLU=A(1) > INS=A(2) > ID TIME AMT CMT F1 F2 GBASE GLU IBASE INS > 1 0 0 1 0 0 0 0 0 0 > 1 0 0 2 0 0 0 0 0 0 > 1 0.001 0 1 0 0 0 0.0009 0 0.0024 > 1 0.002 1 1 187.14 527.49 4.6785 4.6805 13.187 0.0049 > 1 0.002 1 2 187.14 527.49 4.6785 4.6805 13.187 13.192 > 1 100 0 1 0 0 4.6785 4.6785 13.187 13.187 > 1 100 0 2 0 0 4.6785 4.6785 13.187 13.187 > ---- Control stream -------- > $PROB GLUCOSE AND INSULIN INITIAL VALUE > $DATA test.csv > $INPUT ID TIME CMT AMT DV > $SIM (200070927) ONLYSIM NSUB=1 > > $THETA > 40 ; POP_RGLU MMOL/H/70KG > 40 ; POP_VGLU L/70KG > 4 ; POP_CLGLU L/H/70KG > 1 ; POP_EMXGLU > 10 ; POP_C50GLU MMOL/L > 100 ; POP_RINS NMOL/H/70KG > 40 ; POP_VINS L/70KG > 10 ; POP_CLINS L/H/70KG > 2 ; POP_EMXINS > 10 ; POP_C50INS NMOL/L > > $OMEGA > 0 FIX ;PPV_RGLU > 0 FIX ;PPV_RINS > > $SIGMA > 0.1 ;G_EXP_RUV > 1 ;G_ADD_RUV > 0.1 ;I_EXP_RUV > 1 ;I_ADD_RUV > > $SUBR ADVAN9 TOL=3 > > $MODEL > COMP (GLUCOSE) > COMP (INSULIN) > COMP (GLU EQUILIBRIUM) > COMP (INS EQUILIBRIUM) > > $PK > "FIRST > " COMMON/PRCOMG/IDUM1,IDUM2,IMAX > " INTEGER IDUM1,IDUM2,IMAX > " IMAX=1000000 > > ; GLUCOSE > RGLU=THETA(1)*EXP(ETA(1)) > VGLU=THETA(2) > CLGLU=THETA(3) > EMXGLU=THETA(4) > C50GLU=THETA(5) > > ;INSULIN > RINS=THETA(6)*EXP(ETA(2)) > VINS=THETA(7) > CLINS=THETA(8) > EMXINS=THETA(9) > C50INS=THETA(10) > > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > IBASE=A(4)/VINS ; INSULIN INITIAL > IF (AMT.GT.0) THEN > F1=GBASE*VGLU > F2=IBASE*VINS > ELSE > F1=0 > F2=0 > ENDIF > S1=VGLU > S2=VINS > > $AESINITIAL > INIT=0 > A(3)=5 ; MMOL/L > A(4)=10 ; NMOL/L > $AES > EGLU=A(3)/VGLU > EINS=A(4)/VINS > EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) > EIEFF=EMXINS*EINS/(C50INS+EINS) > E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM > E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM > $DES > DGLU=A(1)/VGLU > DINS=A(2)/VINS > DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) > DIEFF=EMXINS*DINS/(C50INS+DINS) > 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 AMT CMT F1 F2 GBASE GLU IBASE INS > NOAPPEND ONEHEADER NOPRINT FILE=test.fit > ---- test.dat -------------------- > #id,time,cmt,amount,dv > 1,0,1,.,. > 1,0,2,.,. > 1,0.001,1,.,. > 1,0.002,1,1,. > 1,0.002,2,1,. > 1,100,1,.,. > 1,100,2,.,. > > -- > 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, Thank you for sending this very interesting, challenging problem. I will try to explain what is going on here, partly because I may not have understood correctly myself (in which case you are encouraged to correct me), and also because not all NONMEM users may be familiar with ADVAN9. ADVAN9 solves a system of simultaneous algebraic and differential equations. Stuart Beal originally had in mind a system in which the equations are in some way coupled: the algebraic equations may involve amounts in all compartments, some of which are "equilibrium compartments" obtained by solving the algebraic equations, and some of which are obtained by solving differential equations. The differential equations themselves may also involve amounts from both kinds of compartments. With NONMEM V, he changed ADVAN9 so that it could be used to solve a system of only algebraic equations, with no differential equations at all. The TIME data item is irreleveant to such a system, and indeed cannot be used if there are only equilibrium compartments, i.e., no differential equations. You present a problem in which the two sets of compartments are (almost) completely separate. Comptartments 1 and 2 are obtained via differential equations involving only these compartments. Compartments 3 and 4 are obtained via algebraic equations involving only these compartments, and are the equilibrium solutions for compartments 1 and 2, i.e., the amounts that would be seen in 1 and 2 if the system were to run till equilibrium is obtained. TIME is irrelevant for compartments 3 and 4. You want to initialize 1 and 2 with their equilibrium amounts prior to the first dose. But ADVAN9 will not begin solving the system until it advances in TIME. There is no way to obtain the solution for compartments 3 and 4 until time increases. (Another way to achieve what you want is to use ADVAN6 with only compartments 1 and 2, and let the system run for a long enough time to achieve equilibrium, i.e., insert a new record at time 0 and add a value such as 100 to all subsequent TIME values.) Had Stuart seen your attempt to use ADVAN9, he might have considered yet another change: find some way to "force" a call to LSODI1 in ADVAN9 at TIME=0, to solve only the equilibrium part of the system. However, it is possible to obtain what you want using some existing features of PREDPP, and some features that are new to NONMEM VI. Below is a modified version of the data set and control stream that do this. There are three new records in the data set that are inserted before the original first record. The have TIME 0, 1, and 1, but any values would do, so long as time increases from the first to the second. Lets call them records a, b and c. Each results in a call to PK, which is the default ("CALL PK WITH EVERY EVENT RECORD"). There are 2 new data items: EVID and FLG. Record a is essentially irrelevant in this problem, but is needed. Record b provides the advance in time during which ADVAN9 solves A(3) and A(4). A(1) and A(2) are irrelevant. Record c has FLG=1. There will be no advance (because TIME is the same with records b and c), but with the call to PK with record c, $PK sees FLG=1 and saves the values of A(3) and A(4) that were computed during the previous advance. It saves them in variables called GBASE0 and IBASE0. These are recursive variables, a feature that is new to NONMEM VI. Previously, a random variable that is defined conditionally would default to 0 when the condition was not met, e.g., when FLG=0. With this special way of defining them, they retain both their values and their eta derivatives. Your original first data record (now the fourth) has EVID=3, reset. This resets all compartments to 0 so that the amounts introduced by the first 3 records are wiped out. PREDPP sets the comparment initialization flag A0_FLG to 1 with both the first event record and with the reset event record. When A0_FLG=1, $PK sets the amounts in compartments 1 and 2 to GBASE0 and IBASE0. When this is a reset record, then these are the amounts from comartments 3 and 4 that were saved with record c, which initializes them as desired. The new code in $PK: IF (FLG.EQ.1) THEN GBASE0=A(3) IBASE0=A(4) ELSE GBASE0=GBASE0 IBASE0=IBASE0 ENDIF IF (A_0FLG.EQ.1) THEN A_0(1)=GBASE0 A_0(2)=IBASE0 PRINT *,'Setting A_0 at TIME=',ICALL,TIME,GBASE0,IBASE0 ENDIF The PRINT statment is a debugging aid. The second $ERROR block and second table are also debugging aids. I did not bother to delete the extraneous dose records or F1 code, but they are no longer needed. ---------Revised data file testrev.csv #id,time,cmt,amount,dv,evid,flg 1,0,1,.,.,.,0 1,1,1,.,.,.,0 1,1,1,.,.,.,1 1,0,1,.,.,3,0 1,0,2,.,.,.,0 1,0.001,1,.,.,.,0 1,0.002,1,1,.,1,0 1,0.002,2,1,.,1,0 1,100,1,.,.,.,0 1,100,2,.,.,.,0 ---------Revised control stream $PROB GLUCOSE AND INSULIN INITIAL VALUE $DATA testrev.csv $INPUT ID TIME CMT AMT DV EVID FLG $SIM (200070927) ONLYSIM NSUB=1 $THETA 40 ; POP_RGLU MMOL/H/70KG 40 ; POP_VGLU L/70KG 4 ; POP_CLGLU L/H/70KG 1 ; POP_EMXGLU 10 ; POP_C50GLU MMOL/L 100 ; POP_RINS NMOL/H/70KG 40 ; POP_VINS L/70KG 10 ; POP_CLINS L/H/70KG 2 ; POP_EMXINS 10 ; POP_C50INS NMOL/L $OMEGA 0 FIX ;PPV_RGLU 0 FIX ;PPV_RINS $SIGMA 0.1 ;G_EXP_RUV 1 ;G_ADD_RUV 0.1 ;I_EXP_RUV 1 ;I_ADD_RUV $SUBR ADVAN9 TOL=3 $MODEL COMP (GLUCOSE) COMP (INSULIN) COMP (GLU EQUILIBRIUM) COMP (INS EQUILIBRIUM) $PK "FIRST " COMMON/PRCOMG/IDUM1,IDUM2,IMAX " INTEGER IDUM1,IDUM2,IMAX " IMAX=1000000 IF (FLG.EQ.1) THEN GBASE0=A(3) IBASE0=A(4) ELSE GBASE0=GBASE0 IBASE0=IBASE0 ENDIF IF (A_0FLG.EQ.1) THEN A_0(1)=GBASE0 A_0(2)=IBASE0 PRINT *,'Setting A_0 at TIME=',ICALL,TIME,GBASE0,IBASE0 ENDIF ; GLUCOSE RGLU=THETA(1)*EXP(ETA(1)) VGLU=THETA(2) CLGLU=THETA(3) EMXGLU=THETA(4) C50GLU=THETA(5) ;INSULIN RINS=THETA(6)*EXP(ETA(2)) VINS=THETA(7) CLINS=THETA(8) EMXINS=THETA(9) C50INS=THETA(10) GBASE=A(3)/VGLU ; GLUCOSE INITIAL IBASE=A(4)/VINS ; INSULIN INITIAL IF (AMT.GT.0) THEN F1=GBASE*VGLU F2=IBASE*VINS ELSE F1=0 F2=0 ENDIF S1=VGLU S2=VINS $AESINITIAL INIT=0 A(3)=5 ; MMOL/L A(4)=10 ; NMOL/L $AES EGLU=A(3)/VGLU EINS=A(4)/VINS EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) EIEFF=EMXINS*EINS/(C50INS+EINS) E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM $DES DGLU=A(1)/VGLU DINS=A(2)/VINS DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) DIEFF=EMXINS*DINS/(C50INS+DINS) 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 AMT CMT F1 F2 GBASE GLU IBASE INS NOAPPEND ONEHEADER NOPRINT FILE=testrev.fit $ERROR A1=A(1) A2=A(2) A3=A(3) A4=A(4) $TABLE TIME A1 A2 A3 A4 NOPRINT FILE=tabrev

Nick, Thank you for sending this very interesting, challenging problem. I will try to explain what is going on here, partly because I may not have understood correctly myself (in which case you are encouraged to correct me), and also because not all NONMEM users may be familiar with ADVAN9. ADVAN9 solves a system of simultaneous algebraic and differential equations. Stuart Beal originally had in mind a system in which the equations are in some way coupled: the algebraic equations may involve amounts in all compartments, some of which are "equilibrium compartments" obtained by solving the algebraic equations, and some of which are obtained by solving differential equations. The differential equations themselves may also involve amounts from both kinds of compartments. With NONMEM V, he changed ADVAN9 so that it could be used to solve a system of only algebraic equations, with no differential equations at all. The TIME data item is irreleveant to such a system, and indeed cannot be used if there are only equilibrium compartments, i.e., no differential equations. You present a problem in which the two sets of compartments are (almost) completely separate. Comptartments 1 and 2 are obtained via differential equations involving only these compartments. Compartments 3 and 4 are obtained via algebraic equations involving only these compartments, and are the equilibrium solutions for compartments 1 and 2, i.e., the amounts that would be seen in 1 and 2 if the system were to run till equilibrium is obtained. TIME is irrelevant for compartments 3 and 4. You want to initialize 1 and 2 with their equilibrium amounts prior to the first dose. But ADVAN9 will not begin solving the system until it advances in TIME. There is no way to obtain the solution for compartments 3 and 4 until time increases. (Another way to achieve what you want is to use ADVAN6 with only compartments 1 and 2, and let the system run for a long enough time to achieve equilibrium, i.e., insert a new record at time 0 and add a value such as 100 to all subsequent TIME values.) Had Stuart seen your attempt to use ADVAN9, he might have considered yet another change: find some way to "force" a call to LSODI1 in ADVAN9 at TIME=0, to solve only the equilibrium part of the system. However, it is possible to obtain what you want using some existing features of PREDPP, and some features that are new to NONMEM VI. Below is a modified version of the data set and control stream that do this. There are three new records in the data set that are inserted before the original first record. The have TIME 0, 1, and 1, but any values would do, so long as time increases from the first to the second. Lets call them records a, b and c. Each results in a call to PK, which is the default ("CALL PK WITH EVERY EVENT RECORD"). There are 2 new data items: EVID and FLG. Record a is essentially irrelevant in this problem, but is needed. Record b provides the advance in time during which ADVAN9 solves A(3) and A(4). A(1) and A(2) are irrelevant. Record c has FLG=1. There will be no advance (because TIME is the same with records b and c), but with the call to PK with record c, $PK sees FLG=1 and saves the values of A(3) and A(4) that were computed during the previous advance. It saves them in variables called GBASE0 and IBASE0. These are recursive variables, a feature that is new to NONMEM VI. Previously, a random variable that is defined conditionally would default to 0 when the condition was not met, e.g., when FLG=0. With this special way of defining them, they retain both their values and their eta derivatives. Your original first data record (now the fourth) has EVID=3, reset. This resets all compartments to 0 so that the amounts introduced by the first 3 records are wiped out. PREDPP sets the comparment initialization flag A0_FLG to 1 with both the first event record and with the reset event record. When A0_FLG=1, $PK sets the amounts in compartments 1 and 2 to GBASE0 and IBASE0. When this is a reset record, then these are the amounts from comartments 3 and 4 that were saved with record c, which initializes them as desired. The new code in $PK: IF (FLG.EQ.1) THEN GBASE0=A(3) IBASE0=A(4) ELSE GBASE0=GBASE0 IBASE0=IBASE0 ENDIF IF (A_0FLG.EQ.1) THEN A_0(1)=GBASE0 A_0(2)=IBASE0 PRINT *,'Setting A_0 at TIME=',ICALL,TIME,GBASE0,IBASE0 ENDIF The PRINT statment is a debugging aid. The second $ERROR block and second table are also debugging aids. I did not bother to delete the extraneous dose records or F1 code, but they are no longer needed. ---------Revised data file testrev.csv #id,time,cmt,amount,dv,evid,flg 1,0,1,.,.,.,0 1,1,1,.,.,.,0 1,1,1,.,.,.,1 1,0,1,.,.,3,0 1,0,2,.,.,.,0 1,0.001,1,.,.,.,0 1,0.002,1,1,.,1,0 1,0.002,2,1,.,1,0 1,100,1,.,.,.,0 1,100,2,.,.,.,0 ---------Revised control stream $PROB GLUCOSE AND INSULIN INITIAL VALUE $DATA testrev.csv $INPUT ID TIME CMT AMT DV EVID FLG $SIM (200070927) ONLYSIM NSUB=1 $THETA 40 ; POP_RGLU MMOL/H/70KG 40 ; POP_VGLU L/70KG 4 ; POP_CLGLU L/H/70KG 1 ; POP_EMXGLU 10 ; POP_C50GLU MMOL/L 100 ; POP_RINS NMOL/H/70KG 40 ; POP_VINS L/70KG 10 ; POP_CLINS L/H/70KG 2 ; POP_EMXINS 10 ; POP_C50INS NMOL/L $OMEGA 0 FIX ;PPV_RGLU 0 FIX ;PPV_RINS $SIGMA 0.1 ;G_EXP_RUV 1 ;G_ADD_RUV 0.1 ;I_EXP_RUV 1 ;I_ADD_RUV $SUBR ADVAN9 TOL=3 $MODEL COMP (GLUCOSE) COMP (INSULIN) COMP (GLU EQUILIBRIUM) COMP (INS EQUILIBRIUM) $PK "FIRST " COMMON/PRCOMG/IDUM1,IDUM2,IMAX " INTEGER IDUM1,IDUM2,IMAX " IMAX=1000000 IF (FLG.EQ.1) THEN GBASE0=A(3) IBASE0=A(4) ELSE GBASE0=GBASE0 IBASE0=IBASE0 ENDIF IF (A_0FLG.EQ.1) THEN A_0(1)=GBASE0 A_0(2)=IBASE0 PRINT *,'Setting A_0 at TIME=',ICALL,TIME,GBASE0,IBASE0 ENDIF ; GLUCOSE RGLU=THETA(1)*EXP(ETA(1)) VGLU=THETA(2) CLGLU=THETA(3) EMXGLU=THETA(4) C50GLU=THETA(5) ;INSULIN RINS=THETA(6)*EXP(ETA(2)) VINS=THETA(7) CLINS=THETA(8) EMXINS=THETA(9) C50INS=THETA(10) GBASE=A(3)/VGLU ; GLUCOSE INITIAL IBASE=A(4)/VINS ; INSULIN INITIAL IF (AMT.GT.0) THEN F1=GBASE*VGLU F2=IBASE*VINS ELSE F1=0 F2=0 ENDIF S1=VGLU S2=VINS $AESINITIAL INIT=0 A(3)=5 ; MMOL/L A(4)=10 ; NMOL/L $AES EGLU=A(3)/VGLU EINS=A(4)/VINS EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) EIEFF=EMXINS*EINS/(C50INS+EINS) E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM $DES DGLU=A(1)/VGLU DINS=A(2)/VINS DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) DIEFF=EMXINS*DINS/(C50INS+DINS) 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 AMT CMT F1 F2 GBASE GLU IBASE INS NOAPPEND ONEHEADER NOPRINT FILE=testrev.fit $ERROR A1=A(1) A2=A(2) A3=A(3) A4=A(4) $TABLE TIME A1 A2 A3 A4 NOPRINT FILE=tabrev

Nick, I found another (simple !) trick how to do it: give steady state dose with F1=0. I tried it, and it is working. Leonid $DATA ../data.csv IGNORE=C $INPUT C,ID,TIME,CMT,AMT,DV,SS,II,MDV,EVID ; GLUCOSE RGLU=THETA(1)*EXP(ETA(1)) VGLU=THETA(2) CLGLU=THETA(3) EMXGLU=THETA(4) C50GLU=THETA(5) ;INSULIN RINS=THETA(6)*EXP(ETA(2)) VINS=THETA(7) CLINS=THETA(8) EMXINS=THETA(9) C50INS=THETA(10) GBASE=A(3)/VGLU ; GLUCOSE INITIAL IBASE=A(4)/VINS ; INSULIN INITIAL IF (TIME.EQ.0) F1=0 S1=VGLU S2=VINS output: ID TIME AMT CMT F1 F2 GBAS GLU IBAS INS 1.0000E+00 0.0000E+00 1.0000E+00 1.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 4.6785E+00 0.0000E+00 1.3187E+01 1.0000E+00 1.0000E-03 0.0000E+00 1.0000E+00 0.0000E+00 0.0000E+00 4.6785E+00 4.6785E+00 1.3187E+01 1.3187E+01 1.0000E+00 2.0000E-03 0.0000E+00 1.0000E+00 0.0000E+00 0.0000E+00 4.6785E+00 4.6785E+00 1.3187E+01 1.3187E+01 1.0000E+00 2.0000E-03 0.0000E+00 2.0000E+00 0.0000E+00 0.0000E+00 4.6785E+00 4.6785E+00 1.3187E+01 1.3187E+01 data: C,id,time,cmt,amount,dv,SS,II,MDV,EVID .,1,0,1,1,.,1,12,1,1 .,1,0.001,1,.,.,.,.,0,0 .,1,0.002,1,.,.,.,.,0,0 .,1,0.002,2,.,.,.,.,0,0 .,1,100,1,.,.,.,.,0,0 .,1,100,2,.,.,.,.,0,0 Nick Holford wrote: > Leonid and others, > > We are sorry we didn't really explain what this code and data are trying to do. > It was 11:30 pm when we wrote the email and it all seemed very clear to us over > a glass of excellent port (thanks Steve D!). > > The objective is to find a way to get the initial conditions for differential equations which do not have a simple algebraic solution. > > The example is for a glucose insulin model with feedback of glucose on insulin > (GEFF) and insulin on glucose (IEFF) using an Emax model in each case. Constant > endogenous input (RGLU and RINS) of glucose and insulin are used and we want to > start the system at its steady state value. We do not know how to solve this > system for the initial conditions by simple algebra. More complex feedback > functions e.g. sigmoid Emax have no explicit solution. > > We are using the $AES feature of ADVAN9 to numerically compute the initial > conditions by asking for the equilibrium solution based on the system shown in > $DES. The equilibrium solution in A(3) is the glucose steady state amount and > in A(4) is the insulin steady state amount. We use these solutions to compute > the steady state concs by dividing the equilibrium compartment amount by volume > (GBASE=A(3)/VGLU and IBASE=A(4)/VINS). > > The design of the data aims to illustrate that the system is correctly initialized > by predicting glucose and insulin at early times (time<=0.002) and at a later > time which is expected to be very close to steady state (Time=100). If the system > is correctly initialized then the values at time=0 and time=100 should be the same. > > The output shows that the system is not correctly initialized at time=0 (both GLU and INS are 0). It also shows that the equilibrium equation solution at is not available at time=0 or at time=0.001 (GBASE and IBASE are 0). > > However, after a nominal input AMT=1 to the glucose (CMT=1) and insulin (CMT=2) > compartments with F1=GBASE and F2=IBASE at time=0.002 we can show that $AES has > correctly computed the steady state value (compare IBASE with INS and GBASE > with GLU at time=100). We know from a Berkeley Madonna simulation of the same > system (using its root finder to get the initial conditions) that the values > predicted at time=100 are correct. It is almost correct at the second > time=0.002 event after the first two compartments have been loaded with the > eventual steady state amounts (F1=GBASE*VGLU and F2=IBASE*VINS). There is a > small error due to input from time=0 into both the GLU and INS compartments. > > So the way we have coded $AEINITIAL and $AES can predict the value we want but > not at time=0(see GBASE and IBASE at time=0). This seems to be either a feature > (i.e. Stuart had some higher purpose in mind to ensure that $AES was only > evaluated after the first non-zero time event -- but this purpose is not clear > to us) or a bug (i.e. $AES is not correctly evaluated at time=0). > > Leonid points out that in the code we sent that A(1) and A(2) are not > explicitly initialized at time=0. We did try doing this with this NMVI feature: > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > IBASE=A(4)/VINS ; INSULIN INITIAL > A_0(1)=GBASE ; set initial value for glucose in A(1) > A_0(2)=IBASE ; set initial value for insulin in A(2) > but it makes no difference to the output because GBASE and IBASE are not still > zero at time=0. > > Leonid asks if it allowed to access A(3) and A(4) in $PK. The answer is that > NMVI does not complain but NMV produces an NM-TRAN error. I had understood that > access to A() in $PK was a new feature in NMVI and had been used for other > purposes e.g. stochastic differential equations. > > It does seem that the value available in $PK may be delayed (which is why we don't get the correct IBASE and GBASE at time=0.001). There are tricky workarounds for this problem as Leonid indicates. > > But more fundamentally we want to know why doesn't $AES compute its solution so > that it is available at time=0? There is no computational reason we can think > of for this limitation. > > 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). > > Nick and Carl > > Leonid Gibiansky wrote: > > > Hi Nick, > > It would help if you prepare a cleaner example that demonstrate the > > problem. What exactly are you trying to achieve (A(?)=?). A(1) and A(2) > > are not initialized in the code below, so they are zeros as they should. > > Is this about A3, A4? > > > > If yes, the problem could be in communication between $AES and $PK. Is > > it allowed to use A() in the PK block? > > > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > > > IBASE=A(4)/VINS ; INSULIN INITIAL > > From help: > > " Right-hand quantities (of PK block can contain - LG) in assignment > > statement and in conditions: > > Data item labels specified on the $INPUT statement. > > THETA(n). > > ETA(n) (Used if the data are population.) > > PK-defined items that appeared earlier as left-hand quantities. > > > > I think, PK and blocks below can communicate via COMMONs but you need to > > take special action (via COMRES) to make them available, and they have a > > 1-record delay (next record knows about COM() of the previous record). I > > could be wrong on that, did it just once, and it was related to ERROR, > > not AES, but we put extra record with the same time, just to pass > > variable to PK > > > > Leonid > > > > Nick Holford wrote: > > > > > Hi, > > > [Second attempt - there was an error in the copy of the control stream in our > > > first attempt] > > > Carl Kirkpatrick and I have been trying to use $AES to compute initial values > > > for use with $DES. > > > We can get $AES to produce the correct equilibrium values at times after time=0 > > > but cannot access these values at time=0 or at time=0.001. Because the > > > equilibrium (i.e. initial condition) values are not available at time=0 we > > > cannot use the NMVI method for initalizing the amounts in the compartments > > > using A_0(1)=GBASE and A_0(2)=IBASE. > > > > > > We have inserted dosing records using the old NMV bioavailability trick to try > > > and initialize the compartments at time=0.002. This almost gets the correct > > > solution at time=0.002 in both compartments. > > > > > > Can anyone explain why the algebraic equation results are not available at > > > time=0 and time=0.001? > > > > > > Nick and Carl > > > > > > --- table file output ----- GLU=A(1) INS=A(2) > > > ID TIME AMT CMT F1 F2 GBASE GLU IBASE INS > > > 1 0 0 1 0 0 0 0 0 0 > > > 1 0 0 2 0 0 0 0 0 0 > > > 1 0.001 0 1 0 0 0 0.0009 0 0.0024 > > > 1 0.002 1 1 187.14 527.49 4.6785 4.6805 13.187 0.0049 > > > 1 0.002 1 2 187.14 527.49 4.6785 4.6805 13.187 13.192 > > > 1 100 0 1 0 0 4.6785 4.6785 13.187 13.187 > > > 1 100 0 2 0 0 4.6785 4.6785 13.187 13.187 > > > ---- Control stream -------- > > > $PROB GLUCOSE AND INSULIN INITIAL VALUE > > > $DATA test.csv > > > $INPUT ID TIME CMT AMT DV > > > $SIM (200070927) ONLYSIM NSUB=1 > > > > > > $THETA > > > 40 ; POP_RGLU MMOL/H/70KG > > > 40 ; POP_VGLU L/70KG > > > 4 ; POP_CLGLU L/H/70KG > > > 1 ; POP_EMXGLU > > > 10 ; POP_C50GLU MMOL/L > > > 100 ; POP_RINS NMOL/H/70KG > > > 40 ; POP_VINS L/70KG > > > 10 ; POP_CLINS L/H/70KG > > > 2 ; POP_EMXINS > > > 10 ; POP_C50INS NMOL/L > > > > > > $OMEGA > > > 0 FIX ;PPV_RGLU > > > 0 FIX ;PPV_RINS > > > > > > $SIGMA > > > 0.1 ;G_EXP_RUV > > > 1 ;G_ADD_RUV > > > 0.1 ;I_EXP_RUV > > > 1 ;I_ADD_RUV > > > > > > $SUBR ADVAN9 TOL=3 > > > > > > $MODEL > > > COMP (GLUCOSE) > > > COMP (INSULIN) > > > COMP (GLU EQUILIBRIUM) > > > COMP (INS EQUILIBRIUM) > > > > > > $PK > > > "FIRST > > > " COMMON/PRCOMG/IDUM1,IDUM2,IMAX > > > " INTEGER IDUM1,IDUM2,IMAX > > > " IMAX=1000000 > > > > > > ; GLUCOSE > > > RGLU=THETA(1)*EXP(ETA(1)) > > > VGLU=THETA(2) > > > CLGLU=THETA(3) > > > EMXGLU=THETA(4) > > > C50GLU=THETA(5) > > > > > > ;INSULIN > > > RINS=THETA(6)*EXP(ETA(2)) > > > VINS=THETA(7) > > > CLINS=THETA(8) > > > EMXINS=THETA(9) > > > C50INS=THETA(10) > > > > > > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > > > IBASE=A(4)/VINS ; INSULIN INITIAL > > > IF (AMT.GT.0) THEN > > > F1=GBASE*VGLU > > > F2=IBASE*VINS > > > ELSE > > > F1=0 > > > F2=0 > > > ENDIF > > > S1=VGLU > > > S2=VINS > > > > > > $AESINITIAL > > > INIT=0 > > > A(3)=5 ; MMOL/L > > > A(4)=10 ; NMOL/L > > > $AES > > > EGLU=A(3)/VGLU > > > EINS=A(4)/VINS > > > EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) > > > EIEFF=EMXINS*EINS/(C50INS+EINS) > > > E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM > > > E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM > > > $DES > > > DGLU=A(1)/VGLU > > > DINS=A(2)/VINS > > > DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) > > > DIEFF=EMXINS*DINS/(C50INS+DINS) > > > 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 AMT CMT F1 F2 GBASE GLU IBASE INS > > > NOAPPEND ONEHEADER NOPRINT FILE=test.fit > > > ---- test.dat -------------------- > > > #id,time,cmt,amount,dv > > > 1,0,1,.,. > > > 1,0,2,.,. > > > 1,0.001,1,.,. > > > 1,0.002,1,1,. > > > 1,0.002,2,1,. > > > 1,100,1,.,. > > > 1,100,2,.,. > > > > > > -- > > > 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 > > -- > 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

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. On Fri, 28 Sep 2007 10:22:01 -0400, "Leonid Gibiansky" <[EMAIL PROTECTED]> said: > Nick, > I found another (simple !) trick how to do it: give steady state dose > with F1=0. I tried it, and it is working. > Leonid > > $DATA ../data.csv IGNORE=C > $INPUT C,ID,TIME,CMT,AMT,DV,SS,II,MDV,EVID > > ; GLUCOSE > RGLU=THETA(1)*EXP(ETA(1)) > VGLU=THETA(2) > CLGLU=THETA(3) > EMXGLU=THETA(4) > C50GLU=THETA(5) > > ;INSULIN > RINS=THETA(6)*EXP(ETA(2)) > VINS=THETA(7) > CLINS=THETA(8) > EMXINS=THETA(9) > C50INS=THETA(10) > > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > IBASE=A(4)/VINS ; INSULIN INITIAL > IF (TIME.EQ.0) F1=0 > S1=VGLU > S2=VINS > > > output: > ID TIME AMT CMT F1 > F2 GBAS GLU IBAS INS > 1.0000E+00 0.0000E+00 1.0000E+00 1.0000E+00 0.0000E+00 > 0.0000E+00 0.0000E+00 4.6785E+00 0.0000E+00 1.3187E+01 > 1.0000E+00 1.0000E-03 0.0000E+00 1.0000E+00 0.0000E+00 > 0.0000E+00 4.6785E+00 4.6785E+00 1.3187E+01 1.3187E+01 > 1.0000E+00 2.0000E-03 0.0000E+00 1.0000E+00 0.0000E+00 > 0.0000E+00 4.6785E+00 4.6785E+00 1.3187E+01 1.3187E+01 > 1.0000E+00 2.0000E-03 0.0000E+00 2.0000E+00 0.0000E+00 > 0.0000E+00 4.6785E+00 4.6785E+00 1.3187E+01 1.3187E+01 > > data: > C,id,time,cmt,amount,dv,SS,II,MDV,EVID > .,1,0,1,1,.,1,12,1,1 > .,1,0.001,1,.,.,.,.,0,0 > .,1,0.002,1,.,.,.,.,0,0 > .,1,0.002,2,.,.,.,.,0,0 > .,1,100,1,.,.,.,.,0,0 > .,1,100,2,.,.,.,.,0,0 > > > Nick Holford wrote: > > Leonid and others, > > > > We are sorry we didn't really explain what this code and data are trying to > > do. It was 11:30 pm when we wrote the email and it all seemed very clear to > > us over a glass of excellent port (thanks Steve D!). > > > > The objective is to find a way to get the initial conditions for > > differential equations which do not have a simple algebraic solution. > > > > The example is for a glucose insulin model with feedback of glucose on > > insulin (GEFF) and insulin on glucose (IEFF) using an Emax model in each > > case. Constant endogenous input (RGLU and RINS) of glucose and insulin are > > used and we want to start the system at its steady state value. We do not > > know how to solve this system for the initial conditions by simple algebra. > > More complex feedback functions e.g. sigmoid Emax have no explicit solution. > > > > We are using the $AES feature of ADVAN9 to numerically compute the initial > > conditions by asking for the equilibrium solution based on the system shown > > in $DES. The equilibrium solution in A(3) is the glucose steady state > > amount and in A(4) is the insulin steady state amount. We use these > > solutions to compute the steady state concs by dividing the equilibrium > > compartment amount by volume (GBASE=A(3)/VGLU and IBASE=A(4)/VINS). > > > > The design of the data aims to illustrate that the system is correctly > > initialized by predicting glucose and insulin at early times (time<=0.002) > > and at a later time which is expected to be very close to steady state > > (Time=100). If the system is correctly initialized then the values at > > time=0 and time=100 should be the same. > > > > The output shows that the system is not correctly initialized at time=0 > > (both GLU and INS are 0). It also shows that the equilibrium equation > > solution at is not available at time=0 or at time=0.001 (GBASE and IBASE > > are 0). > > > > However, after a nominal input AMT=1 to the glucose (CMT=1) and insulin > > (CMT=2) compartments with F1=GBASE and F2=IBASE at time=0.002 we can show > > that $AES has correctly computed the steady state value (compare IBASE with > > INS and GBASE with GLU at time=100). We know from a Berkeley Madonna > > simulation of the same system (using its root finder to get the initial > > conditions) that the values predicted at time=100 are correct. It is almost > > correct at the second time=0.002 event after the first two compartments > > have been loaded with the eventual steady state amounts (F1=GBASE*VGLU and > > F2=IBASE*VINS). There is a small error due to input from time=0 into both > > the GLU and INS compartments. > > > > So the way we have coded $AEINITIAL and $AES can predict the value we want > > but not at time=0(see GBASE and IBASE at time=0). This seems to be either a > > feature (i.e. Stuart had some higher purpose in mind to ensure that $AES > > was only evaluated after the first non-zero time event -- but this purpose > > is not clear to us) or a bug (i.e. $AES is not correctly evaluated at > > time=0). > > > > Leonid points out that in the code we sent that A(1) and A(2) are not > > explicitly initialized at time=0. We did try doing this with this NMVI > > feature: > > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > > IBASE=A(4)/VINS ; INSULIN INITIAL > > A_0(1)=GBASE ; set initial value for glucose in A(1) > > A_0(2)=IBASE ; set initial value for insulin in A(2) > > but it makes no difference to the output because GBASE and IBASE are not > > still zero at time=0. > > > > Leonid asks if it allowed to access A(3) and A(4) in $PK. The answer is > > that NMVI does not complain but NMV produces an NM-TRAN error. I had > > understood that access to A() in $PK was a new feature in NMVI and had been > > used for other purposes e.g. stochastic differential equations. > > > > It does seem that the value available in $PK may be delayed (which is why > > we don't get the correct IBASE and GBASE at time=0.001). There are tricky > > workarounds for this problem as Leonid indicates. > > > > But more fundamentally we want to know why doesn't $AES compute its > > solution so that it is available at time=0? There is no computational > > reason we can think of for this limitation. > > > > 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). > > > > Nick and Carl > > > > Leonid Gibiansky wrote: > >> Hi Nick, > >> It would help if you prepare a cleaner example that demonstrate the > >> problem. What exactly are you trying to achieve (A(?)=?). A(1) and A(2) > >> are not initialized in the code below, so they are zeros as they should. > >> Is this about A3, A4? > >> > >> If yes, the problem could be in communication between $AES and $PK. Is > >> it allowed to use A() in the PK block? > >> > GBASE=A(3)/VGLU ; GLUCOSE INITIAL > >> > IBASE=A(4)/VINS ; INSULIN INITIAL > >> From help: > >> " Right-hand quantities (of PK block can contain - LG) in assignment > >> statement and in conditions: > >> Data item labels specified on the $INPUT statement. > >> THETA(n). > >> ETA(n) (Used if the data are population.) > >> PK-defined items that appeared earlier as left-hand quantities. > >> > >> I think, PK and blocks below can communicate via COMMONs but you need to > >> take special action (via COMRES) to make them available, and they have a > >> 1-record delay (next record knows about COM() of the previous record). I > >> could be wrong on that, did it just once, and it was related to ERROR, > >> not AES, but we put extra record with the same time, just to pass > >> variable to PK > >> > >> Leonid > >> > >> Nick Holford wrote: > >>> Hi, > >>> [Second attempt - there was an error in the copy of the control stream in > >>> our first attempt] > >>> Carl Kirkpatrick and I have been trying to use $AES to compute initial > >>> values for use with $DES. > >>> We can get $AES to produce the correct equilibrium values at times after > >>> time=0 but cannot access these values at time=0 or at time=0.001. Because > >>> the equilibrium (i.e. initial condition) values are not available at > >>> time=0 we cannot use the NMVI method for initalizing the amounts in the > >>> compartments using A_0(1)=GBASE and A_0(2)=IBASE. > >>> > >>> We have inserted dosing records using the old NMV bioavailability trick > >>> to try and initialize the compartments at time=0.002. This almost gets > >>> the correct solution at time=0.002 in both compartments. > >>> > >>> Can anyone explain why the algebraic equation results are not available > >>> at time=0 and time=0.001? > >>> > >>> Nick and Carl > >>> > >>> --- table file output ----- GLU=A(1) > >>> INS=A(2) > >>> ID TIME AMT CMT F1 F2 GBASE GLU IBASE INS > >>> 1 0 0 1 0 0 0 0 0 0 > >>> 1 0 0 2 0 0 0 0 0 0 > >>> 1 0.001 0 1 0 0 0 0.0009 0 > >>> 0.0024 > >>> 1 0.002 1 1 187.14 527.49 4.6785 4.6805 13.187 > >>> 0.0049 > >>> 1 0.002 1 2 187.14 527.49 4.6785 4.6805 13.187 > >>> 13.192 > >>> 1 100 0 1 0 0 4.6785 4.6785 13.187 > >>> 13.187 > >>> 1 100 0 2 0 0 4.6785 4.6785 13.187 > >>> 13.187 > >>> ---- Control stream -------- > >>> $PROB GLUCOSE AND INSULIN INITIAL VALUE > >>> $DATA test.csv > >>> $INPUT ID TIME CMT AMT DV > >>> $SIM (200070927) ONLYSIM NSUB=1 > >>> > >>> $THETA > >>> 40 ; POP_RGLU MMOL/H/70KG > >>> 40 ; POP_VGLU L/70KG > >>> 4 ; POP_CLGLU L/H/70KG > >>> 1 ; POP_EMXGLU > >>> 10 ; POP_C50GLU MMOL/L > >>> 100 ; POP_RINS NMOL/H/70KG > >>> 40 ; POP_VINS L/70KG > >>> 10 ; POP_CLINS L/H/70KG > >>> 2 ; POP_EMXINS > >>> 10 ; POP_C50INS NMOL/L > >>> > >>> $OMEGA > >>> 0 FIX ;PPV_RGLU > >>> 0 FIX ;PPV_RINS > >>> > >>> $SIGMA > >>> 0.1 ;G_EXP_RUV > >>> 1 ;G_ADD_RUV > >>> 0.1 ;I_EXP_RUV > >>> 1 ;I_ADD_RUV > >>> > >>> $SUBR ADVAN9 TOL=3 > >>> > >>> $MODEL > >>> COMP (GLUCOSE) > >>> COMP (INSULIN) > >>> COMP (GLU EQUILIBRIUM) > >>> COMP (INS EQUILIBRIUM) > >>> > >>> $PK > >>> "FIRST > >>> " COMMON/PRCOMG/IDUM1,IDUM2,IMAX > >>> " INTEGER IDUM1,IDUM2,IMAX > >>> " IMAX=1000000 > >>> > >>> ; GLUCOSE > >>> RGLU=THETA(1)*EXP(ETA(1)) > >>> VGLU=THETA(2) > >>> CLGLU=THETA(3) > >>> EMXGLU=THETA(4) > >>> C50GLU=THETA(5) > >>> > >>> ;INSULIN > >>> RINS=THETA(6)*EXP(ETA(2)) > >>> VINS=THETA(7) > >>> CLINS=THETA(8) > >>> EMXINS=THETA(9) > >>> C50INS=THETA(10) > >>> > >>> GBASE=A(3)/VGLU ; GLUCOSE INITIAL > >>> IBASE=A(4)/VINS ; INSULIN INITIAL > >>> IF (AMT.GT.0) THEN > >>> F1=GBASE*VGLU > >>> F2=IBASE*VINS > >>> ELSE > >>> F1=0 > >>> F2=0 > >>> ENDIF > >>> S1=VGLU > >>> S2=VINS > >>> > >>> $AESINITIAL > >>> INIT=0 > >>> A(3)=5 ; MMOL/L > >>> A(4)=10 ; NMOL/L > >>> $AES > >>> EGLU=A(3)/VGLU > >>> EINS=A(4)/VINS > >>> EGEFF=EMXGLU*EGLU/(C50GLU+EGLU) > >>> EIEFF=EMXINS*EINS/(C50INS+EINS) > >>> E(3)=RGLU - CLGLU*(1+EIEFF)*EGLU ; GLUCOSE EQUILIBRIUM > >>> E(4)=RINS*(1+EGEFF) - CLINS*EINS ; INSULIN EQUILIBRIUM > >>> $DES > >>> DGLU=A(1)/VGLU > >>> DINS=A(2)/VINS > >>> DGEFF=EMXGLU*DGLU/(C50GLU+DGLU) > >>> DIEFF=EMXINS*DINS/(C50INS+DINS) > >>> 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 AMT CMT F1 F2 GBASE GLU IBASE INS > >>> NOAPPEND ONEHEADER NOPRINT FILE=test.fit > >>> ---- test.dat -------------------- > >>> #id,time,cmt,amount,dv > >>> 1,0,1,.,. > >>> 1,0,2,.,. > >>> 1,0.001,1,.,. > >>> 1,0.002,1,1,. > >>> 1,0.002,2,1,. > >>> 1,100,1,.,. > >>> 1,100,2,.,. > >>> > >>> -- > >>> 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 > >>> > >>> > > > > -- > > 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]

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

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]

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

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]