Re: [NMusers] 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
>>> n.holford_at_auckland.ac.nz 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
> n.holford_at_auckland.ac.nz tel:+64(9)373-7599x86730 fax:+64(9)373-7090
> www.health.auckland.ac.nz/pharmacology/staff/nholford
>
>
Received on Fri Sep 28 2007 - 10:22:01 EDT