Re: [NMusers] bug/feature in $AES

Nick and Carl, looks like our emails crossed in cyberspace.
I think I answered most of your questions, but will make some comments
on this latest email.
Look for lines marked *** at the start.

On Fri, 28 Sep 2007 07:49:07 +1200, "Nick Holford"
<n.holford_at_auckland.ac.nz> said:
> 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.

*** Just so. The call to PK at time .001 takes place prior to the
(first) advance.
*** The compartment amounts are all 0, hence so is F1. In the
*** revised data file and control stream that I sent, it is the third
call to PK
*** that takes place after the (first) advance, and now A(3) and A(4)
*** are properly computed.
>
> 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.

*** I think you are correct and could perhaps have persuaded Stuart of
this.
>
> 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).

*** It may well be possible to modify ADVAN9. I will have to think about
it.
*** I suspect that it will take some effort on my part.


>
> 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
--
  Alison Boeckmann
  alisonboeckmann_at_fastmail.fm

Received on Thu Sep 27 2007 - 17:52:11 EDT