From: "Eleveld, DJ" d.j.eleveld@anest.umcg.nl
Subject: [NMusers] Problems with recursive PRED models...
Date: Thu, 22 Sep 2005 15:18:38 +0200

Hello Everyone,

I have a question about 'recursive PRED' models...

I have a model in which the bioavailability of a compartment depends on the
state of other compartments as well as a model parameter.  Since
bioavailability must be set in $PK but the compartmental amounts are not
available there, I used COMRES to provide $PK with information from $ERROR
using COM(1) following advice from the very helpful readers of this list.
The doses are directed to the correct compartment (CENTRAL or POTENT) with a
CMT field in the $DATA. 

In my searching of the mailing list archives regarding recursive PRED
routines I found the following quote from Alison Fri Sep 27 10:30:50 1996:

...any left hand quantity in abbreviated code that depends on random
variables (ETA and/or EPS) must be computed EXACTLY ONCE AND
UNCONDITIONALLY. 
If an either/or choice is involved, indicator variables (e.g., Q and 1-Q in
many examples) can be used, because this rule is obeyed. 
But if the computation involves random variables computed with a prior data
or event record, it cannot be done correctly.

As I understand it F5 (bioavailability) is a left hand quantity that is
calculated exactly once and unconditionally, so that is not a problem.
However, the calculation of F5 involves a random variable (POTR) as well as
prior data because it depends on the state of other compartments (A(4) via
COM(1)).  If this interpretation is correct as well as the comment from
Alison (I assume it is) then I must conclude that this model cannot be done
correctly in NONMEM.  I am of course reluctant to make such a sweeping
conclusion without expert advice.  

This model seems to simulate just fine and estimations on individual data
seem fine also.  Straightforward FO estimation leads to high condition
numbers with the data I have.   Setting some $OMEGA values to FIXED values
gets the condition number to less than 1000.  FOCE methods often abort with
a floating point overflow error or rounding errors, even those with a
corresponding FO estimation having a low condition number.  I have had only
one FOCE estimation run to completion.  The results and the prediction are
reasonable and parameters are in the expected range.  However, even small
changes to the initial estimates cause floating point overflow errors and it
is difficult to have confidence in the results when the estimation breaks so
easily.  Before diving deeper into the specifics of this problem I think it
is prudent to make sure that NONMEM can correctly implement the model.

Can anyone tell me if this model is indeed a recursive PRED type and whether
or not the model can be correctly implemented in NONMEM?

Thanks very much,

Doug Eleveld

The following control file aborts (crashes) after 2 iterations with a
floating point overflow error (Open Watcom compiler):

$PROB  Potentiation fitting
$DATA  potpd.prn IGNORE=C
$INPUT ID TIME DV MDV CFLG AMT RATE CMT V1 V2 V3 CL Q2 Q3
$SUBROUTINES ADVAN6 TOL=4
$ABBREVIATED COMRES=1
$MODEL COMP(CENTRAL DEFDOSE)
       COMP(PERIF1 NOOFF NODOSE) 
       COMP(PERIF2 NOOFF NODOSE) 
       COMP(EFFECT NOOFF NODOSE) 
       COMP(POTENT NOOFF)
$PK
    S1=V1
    KEO=THETA(1)*EXP(ETA(1))
    EC50=THETA(2)*EXP(ETA(2))
    GAMM=THETA(3)*EXP(ETA(3))
    POTR=ABS(THETA(4)+ETA(4))
    POTK=THETA(5)*EXP(ETA(5))
    SCAL=THETA(6)+ETA(6)
    PD1=COM(1)**GAMM             
    PD2=EC50**GAMM
    MNMB=1-PD1/(PD1+PD2)
    F5=POTR*MNMB/1000
$DES
    C1=A(1)/V1                  ; Conc in V1
    C2=A(2)/V2                  ; Conc in V2
    C3=A(3)/V3                  ; Conc in V3
    DADT(1)=-CL*C1+Q2*(C2-C1)+Q3*(C3-C1) ; Change in V1
    DADT(2)=Q2*(C1-C2)          ; Change in V2
    DADT(3)=Q3*(C1-C3)          ; Change in V3
    DADT(4)=(C1-A(4))*KEO       ; Effect compartment conc
    DADT(5)=-POTK*A(5)          ; Decay potentiation 
$ERROR
    CPRE=A(1)/V1
    CEFF=A(4)
    IF (CEFF.LT.0) CEFF=0
    DPOT=1+A(5)                 ; The degree of potentiation
    EPD1=CEFF**GAMM             
    EPD2=EC50**GAMM
    ENMB=1-EPD1/(EPD1+EPD2)
    RPRE=(SCAL+100)*DPOT*ENMB   ; Twitch prediction
    Y=RPRE+ERR(1)
    COM(1)=CEFF
$THETA (0,0.2,0.4)(1000,1800,3000)(3,5,8)
       (0,5,20)(0,0.05,1)(-1 FIXED)
$OMEGA (0.05 FIXED)(0.05)(0.05 FIXED)
       (4 FIXED)(0.1 FIXED)(4 FIXED)
$SIGMA (12)
$ESTIMATION MAX=9999 SIG=3 METHOD=1
       NOABORT POSTHOC PRINT=1
$COVARIANCE PRINT=E
$TABLE ID TIME CPRE CEFF RPRE MNMB DPOT Y DV
       NOPRINT NOHEADER NOAPPEND FILE=potent.txt
$TABLE ID V1 V2 V3 CL Q2 Q3 KEO EC50 GAMM POTR POTK SCAL
       NOPRINT NOHEADER NOAPPEND FILE=potent1.txt

 MONITORING OF SEARCH:

0ITERATION NO.:    0     OBJECTIVE VALUE:  0.1077E+05     NO. OF FUNC.
EVALS.: 8
 CUMULATIVE NO. OF FUNC. EVALS.:    8
 PARAMETER:  0.1000E+00  0.1000E+00  0.1000E+00  0.1000E+00  0.1000E+00
0.1000E+00  0.1000E+00
 GRADIENT:   0.2022E+04  0.1228E+04  0.6046E+03  0.1841E+04 -0.7099E+03
0.1724E+02  0.4353E+02
0ITERATION NO.:    1     OBJECTIVE VALUE:  0.1072E+05     NO. OF FUNC.
EVALS.:10
 CUMULATIVE NO. OF FUNC. EVALS.:   18
 PARAMETER:  0.8500E-01  0.9089E-01  0.9551E-01  0.8634E-01  0.1053E+00
0.9987E-01  0.9968E-01
 GRADIENT:   0.7187E+03 -0.1404E+03  0.3946E+03  0.8371E+03 -0.9275E+03
0.7030E+02 -0.1695E+04

_______________________________________________________