Re: bug/feature in $AES

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]