Re: bug/feature in $AES

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