Re: bug/feature in $AES

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