Appendix

The PREDPP Defined Data Items

Event Identification (EVID) Data Item

Value	Meaning
0	observation event
1	dose event
2	other-type event
3	reset event
4	reset-dose event

See Chapter V for a discussion of “extra” EVID’s: special names XVID1, XVID2, XVID3, XVID4, XVID5.

Time (TIME) Data Item

Time of event.

Dose Amount (AMT) Data Item

Value	Meaning
0	steady-state infusion
positive no.	amount of either bolus dose or infusion

Note: AMT data item must be 0 in a non-dose or non-reset-dose event record.

Dose Rate (RATE) Data Item

Value	Meaning
0	instantaneous bolus dose
-1	zero-order bolus dose, rate modeled in PK; or steady-state infusion, rate modeled in PK
-2	zero-order bolus dose, duration modeled in PK
positive no.	rate of infusion

Note: RATE data item must be 0 in a non-dose or non-reset-dose event record.

Steady-State (SS) Data Item

Value	Meaning
0	not a steady-state dose
1	a steady-state dose
2	a superimposed steady-state dose
3	Same as SS=1 but existing values of A are used as initial estimates (SS6, SS9, SS13)

Note: SS data item must be 0 in a non-dose or non-reset-dose event record.

Interdose Interval (II) Data Item

Value	Meaning
0	single non-steady-state dose or steady-state infusion
positive no.	interdose interval

Note: II data item must be 0 in a non-dose or non-reset-dose event record.

Compartment (CMT) Data Item

Value	Meaning
0	event (dose, reset-dose, observation) applies to the default compartment; has no effect with an other-type or reset event
integer	event (dose, reset-dose, observation) applies to the compartment with this number; for other-type or reset event, compartment with this number is turned on
-integer	for observation event, and if output compartment has this number, output compartment is observed and then turned off; for other-type or reset event, compartment with this number is turned off; not allowed for dose or reset-dose event

Note: With a reset event, a compartment is turned on or off after all compartments are reset to their initial status.

When the integer is 1000, it designates the default compartment for output.

When the integer is 100 and the number of compartments is <=99, it is another way of designating the default compartment for output.

Prediction Compartment (PCMT) Data Item

Value	Meaning
0	prediction is the one that would result for an observation of the default observation compartment
integer	prediction is the one that would result for an observation of the compartment with this number

Note: PCMT data item is ignored in an observation event record.

When the integer is 1000, it designates the default compartment for output.

When the integer is 100 and the number of compartments is <=99, it is another way of designating the default compartment for output.

Call (CALL) Data Item

Value	Meaning
0	has no effect
1	call to ERROR with event record is forced
2	call to PK with event record is forced
3	calls to both PK and ERROR with event record are forced
10	Call to ADVAN is forced. (Only with ADVAN9, ADVAN15, or ADVAN17)

Continuation (CONT) Data Item

Value	Meaning
0	data record ends event record
1	event record is continued on next data record

Additional Dose (ADDL) Data Item

Value	Meaning
0	no additional doses, or dose is a steady-state infusion
integer	number of additional doses of same type

Note: ADDL data item must be 0 in a non-dose or non-reset-dose event record. Also, when ADDL data item is nonzero, II data item must be nonzero.

Another Example: Theophylline

The population theophylline example of Guide I, section F, can be specified using PREDPP, as shown in this appendix. A NONMEM control stream, and PK and ERROR routines, are given on the following pages. A NM-TRAN control stream, using abbreviated PK and ERROR codes, is also given. The results are essentially identical to those given in Guide I, figures 76-83. Unlike the example of Guide I, the data set is not embedded in the control stream; it is found in a separate file, also given in this appendix. (This is unrelated to the fact that PREDPP is used.) A file stream, pointing to the data file, is also given.

This example uses ADVAN2, which implements an analytic computation for the one compartment linear model with first-order absorption. The same example can be handled using ADVAN6 and ADVAN7, although it is less efficient to use subroutines implementing “numerical” computations for general models when an analytic ADVAN for a specific model can be used. However, it may be instructive for the user to examine how all three ADVAN’s can be used for the same example, and to examine additional user-supplied codes to be used with these ADVAN’s. Therefore, this appendix also gives PK, MODEL, DES, and TOL routines to be used with ADVAN6, and PK and MODEL routines to be used with ADVAN7.

Finally, NM-TRAN control streams corresponding to the use of ADVAN6 and ADVAN7 with this example are also given. These use the very same data set as does the NONMEM control stream. However, the last two data item types, EVID and MDV, are ignored. Though NM-TRAN can always properly process these data item types, it can ignore them with this example.

FCON Stream (excerpt)

FILE FSTREAM
PROB THEOPHYLLINE POPULATION DATA
FMTN 3
LDFO 0
DATA        1 0 144 7 0 0 0 0   -1 0 0
ITEM        1 4 7 11 1 0 0 0 0 0 0 0 0 0 0 0
INDX        63200000000
LABL        ID,    DOSE
            TIME,  CP
            WT,    EVID
            MDV
FORM
(5F10.0,2F2.0)
THCN        1000
THTA    3.000000000000000E+00, 8.000000000000000E-02, 4.000000000000000E-02
LOWR    1.000000000000000E-01, 8.000000000000000E-03, 4.000000000000000E-03
UPPR    5.000000000000000E+00, 5.000000000000000E-01, 9.000000000000000E-01
BLST    6.000000000000000E+00, 5.000000000000000E-03, 3.000000000000000E-01
        2.000000000000000E-04, 6.000000000000000E-03, 4.000000000000000E-01
DIAG    4.000000000000000E-01
ESTM  0 450 3 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
      0 0 0 0 -1 -1 0 0 0
COVR  000000100
TABL 1 4   10  20  50  30
SCAT  3   9   1  1  0   0  0  0  0
          0   0  0  -1
SCAT  3   10  1  1  0   0  0  0  0
          0   0  0  -1

Data Set (first records)

1  4.02         0.            79.6 1 1
1               0.   .74         00
1               0.25 2.84        00
1               0.57 6.57        00
1               1.12 10.5        00
1               2.02 9.66        00
1               3.82 8.58        00
1               5.1  8.36        00
1               7.03 7.47        00
1               9.05 6.89        00
1               12.12 5.94       00
1               24.37 3.28       00
2  4.4          0.            72.4 1 1
2               0.    0.         00
[…remaining records elided — see help Guide VIII for full theophylline data set…]

NM-TRAN Control Stream To Be Used With ADVAN2

$PROB THEOPHYLLINE POPULATION DATA
$INPUT     ID DOSE=AMT TIME CP=DV WT
$DATA      DATAFILE (5F10.0)
$SUBROUTINES ADVAN2
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
     CALLFL=1
     KA=THETA(1)+ETA(1)
     K=THETA(2)+ETA(2)
     CL=THETA(3)*WT+ETA(3)
     SC=CL/K/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$ERROR
     Y=F+EPS(1)
$SIGMA .4
$EST    MAXEVAL=450 PRINT=5
$COV
$TABLE              ID DOSE WT TIME
$SCAT               (RES WRES) VS TIME BY ID

PK Routine To Be Used With ADVAN6 and ADVAN7

SUBROUTINE PK(ICALL,IDEF,THETA,IREV,EVTREC,NVNT,INDXS,IRGG,GG,NETAS)
! PK SUBROUTINE FOR THE THEOPHYLLINE POPULATION DATA
! USED WITH ADVAN6 AND ADVAN7, AND TRANS1
! CLEARANCE PROPORTIONAL TO WEIGHT
! SCALING PARAMETER IS VOLUME/WEIGHT SINCE DOSE IS WEIGHT-ADJUSTED
USE NMPRD4P
USE SIZES,                  ONLY: DPSIZE,ISIZE
USE PRDIMS, ONLY: GPRD,HPRD,GERD,HERD,GPKD
USE NMPRD_REAL,ONLY: ETA,EPS
USE NMPRD_INT, ONLY: MSEC=>ISECDER,IFIRSTEM,FIRSTEM
USE NMPRD_INT, ONLY: IQUIT
USE PROCM_INT, ONLY: NEWIND=>PNEWIF
IMPLICIT REAL(KIND=DPSIZE) (A-Z)
REAL(KIND=DPSIZE) :: EVTREC
SAVE
INTEGER(KIND=ISIZE) :: ICALL,IDEF,IREV,NVNT,INDXS,IRGG,NETAS
DIMENSION :: IDEF(7,*),THETA(*),EVTREC(IREV,*),INDXS(*),GG(IRGG,GPKD+1,*)
FIRSTEM=IFIRSTEM
IF (ICALL <= 1) THEN
   CALL ASSOCNMPRD4
   ! SET UP IDEF ARRAY:
   IDEF( 1,0001)= -9
   ! CALL PK ONCE PER INDIV. REC.
   IDEF( 1,0002)= 1
   ! ROW INDEX FOR SCALING PARAMETER FOR COMPT. 2
   IDEF( 3,0002)= 3
   CALL GETETA(ETA)
   IF (IQUIT == 1) RETURN
   RETURN
ENDIF
IF (NEWIND /= 2) THEN
   IF (ICALL == 4) THEN
      CALL SIMETA(ETA)
   ELSE
      CALL GETETA(ETA)
   ENDIF
   IF (IQUIT == 1) RETURN
ENDIF
WT=EVTREC(NVNT,005)
KA=THETA(001)+ETA(001)
KE=THETA(002)+ETA(002)
CL=THETA(003)*WT+ETA(003)
S2=CL/KE/WT
P000001=KA
P000002=KE
! ABSORPTION RATE CONSTANT
GG(0001,1,1)=P000001
! ELIMINATION RATE CONSTANT
GG(0002,1,1)=P000002
! SCALING
GG(0003,1,1)=S2
! [eta and second-derivative blocks elided for brevity]
RETURN
END

MODEL Routine To Be Used With ADVAN6 and ADVAN7

SUBROUTINE MODEL (IDNO,NCM,NPAR,IR,IATT,LINK)
! DEFINES A 1 COMPARTMENT MODEL WITH DRUG DEPOT COMPARTMENT
! COMPT1: DRUG DEPOT COMPT2: CENTRAL COMPT
USE PRMOD_CHAR, ONLY: NAME
USE SIZES,                  ONLY: DPSIZE,ISIZE,SD
USE PRDIMS, ONLY: GPRD,HPRD,GERD,HERD,GPKD
INTEGER(KIND=ISIZE) :: IDNO,NCM,NPAR,IR,IATT,LINK,I,J
DIMENSION :: IATT(IR,*),LINK(IR,*)
SAVE
INTEGER(KIND=ISIZE), DIMENSION (2,7) :: MOD
CHARACTER(LEN=SD), DIMENSION(2) :: CMOD
DATA (MOD(I, 1),I= 1, 2)/&
! INITIAL STATUS: OFF ON
0,1 /
DATA (MOD(I, 2),I= 1, 2)/&
! ON/OFF ALLOWED: YES NO
1,0 /
DATA (MOD(I, 3),I= 1, 2)/&
! DOSE ALLOWED: YES YES
1,1 /
DATA (MOD(I, 4),I= 1, 2)/&
! DEFAULT FOR OBSERVATIONS: NO YES
0,1 /
DATA (MOD(I, 5),I= 1, 2)/&
! DEFAULT FOR DOSES: YES NO
1,0 /
DATA (MOD(I, 6),I= 1, 2)/&
0,0 /
DATA (MOD(I, 7),I= 1, 2)/&
0,0 /
DATA (CMOD(I),I= 1, 2) &
/'DEPOT','CENTRAL'/
FORALL (I=1:2) NAME(I)=CMOD(I)
FORALL (I=1:2,J=1:7) IATT(I,J)=MOD(I,J)
IDNO=9999
NCM= 2
NPAR=002
RETURN
END

DES Routine To Be Used With ADVAN6

SUBROUTINE DES (A,P,T,DADT,IR,DA,DP,DT)
USE NMPRD4P
USE SIZES,                  ONLY: DPSIZE,ISIZE
USE PRDIMS, ONLY: GPRD,HPRD,GERD,HERD,GPKD
USE NMPRD_INT, ONLY: MSEC=>ISECDER,IFIRSTEM,FIRSTEM
USE PRMOD_INT, ONLY: ICALL=>ICALLD,IDEFD
IMPLICIT REAL(KIND=DPSIZE) (A-Z)
SAVE
INTEGER(KIND=ISIZE) :: IR
DIMENSION :: A(*),P(*),DADT(*),DA(IR,*),DP(IR,*),DT(*)
FIRSTEM=IFIRSTEM
IF (ICALL == 1) THEN
   CALL ASSOCNMPRD4
   IDEFD(1)= 0
   IDEFD(2)= 0
   RETURN
ENDIF
DADT(1)=-P(001)*A(1)
DADT(2)=P(001)*A(1)-P(002)*A(2)
IF (FIRSTEM == 1) THEN
   E000004=-P(001)
   F000081=-A(1)
   E000006=P(001)
   E000007=-P(002)
   F000082=A(1)
   F000083=-A(2)
   DA(1,1)=E000004
   DA(2,1)=E000006
   DA(3,1)=E000007
   DP(1,1)=F000081
   DP(2,1)=F000082
   DP(3,1)=F000083
ENDIF
RETURN
END

TOL Routine To Be Used With ADVAN6

! TOL SUBROUTINE FOR THE THEOPHYLLINE POPULATION DATA
SUBROUTINE TOL(NRD,ANRD,NRDC,ANRDC)
USE SIZES,               ONLY: ISIZE
INTEGER(KIND=ISIZE) :: NRD(0:*), ANRD(0:*), NRDC(0:*), ANRDC(0:*)
NRD(1)=5
RETURN
END

NM-TRAN Control Stream To Be Used With ADVAN6

$PROB THEOPHYLLINE POPULATION DATA
$INPUT     ID DOSE=AMT TIME CP=DV WT
$DATA      DATAFILE (5F10.0)
$SUBROUTINES ADVAN6 TOL=5
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
     CALLFL=1
     KA=THETA(1)+ETA(1)
     KE=THETA(2)+ETA(2)
     CL=THETA(3)*WT+ETA(3)
     S2=CL/KE/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$DES
     DADT(1)=-KA*A(1)
     DADT(2)= KA*A(1)-KE*A(2)
$ERROR
     Y=F+EPS(1)
$SIGMA .4
$EST    MAXEVAL=450 PRINT=5
$COV
$TABLE              ID DOSE WT TIME
$SCAT               (RES WRES) VS TIME BY ID

NM-TRAN Control Stream To Be Used With ADVAN7

$PROB THEOPHYLLINE POPULATION DATA
$INPUT     ID DOSE=AMT TIME CP=DV WT
$DATA      DATAFILE (5F10.0)
$SUBROUTINES ADVAN7
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
$PK
;THETA(1)=MEAN ABSORPTION RATE CONSTANT (1/HR)
;THETA(2)=MEAN ELIMINATION RATE CONSTANT (1/HR)
;THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
;SCALING PARAMETER=VOLUME/WT SINCE DOSE IS WEIGHT-ADJUSTED
     CALLFL=1
     K12=THETA(1)+ETA(1)
     K20=THETA(2)+ETA(2)
     CL=THETA(3)*WT+ETA(3)
     S2=CL/K20/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$ERROR
     Y=F+EPS(1)
$SIGMA .4
$EST    MAXEVAL=450 PRINT=5
$COV
$TABLE              ID DOSE WT TIME
$SCAT               (RES WRES) VS TIME BY ID

MTIME Details

III.A. General Description

Model event times are additional PK parameters defined in the PK routine or $PK block. A model event time is not associated with any compartment, but, like an absorption lag time, defines a time to which the system is advanced. When the time is reached, certain new indicator variables are set and a call to PK is made. At this call (and/or subsequent to this call) PK or DES or AES or ERROR can use the indicator variables to change some aspect of the system, e.g., a term in a differential equation, or the rate of an infusion. This feature may be used with any ADVAN routine and does not require that extra compartments be added to the model or that dummy doses be included in the data set. There may be up to PCT model event times, where PCT is defined in PSIZES and TSIZES. Its current value is 30.

Model times are independent of non-event dose times. See below (Rules for coincidence of times) when model times happen to coincide with other times.

In abbreviated code, the model event times have reserved names MTIME(i), where i is an integer from 1 to PCT.

The indicator variables currently implemented have reserved names MNEXT(i) and MPAST(i), each corresponding to the MTIME(i) variable with the same index i. An indicator variable may be used in the PK, ERROR, DES, or AES routine (as a right-hand quantity in abbreviated codes). An indicator variable may not be used on the right (other than in a write/print statement) unless the corresponding MTIME has already been defined. It is not required that one or more indicator variables be used.

There is no requirement that the MTIME(i) be in ascending order. That is, MTIME(i) may be less than, equal to, or greater than MTIME(i+1). Any MTIME(i) may be negative or have the value 0, in which case MPAST(i)=1 always and MNEXT(i)=0 always.

There may be gaps in the sequence of any length. E.g., PK may define MTIME(i) and MTIME(i+2) but not MTIME(i+1). This has the same effect as including the code MTIME(i+1)=0.

PK may redefine MTIME(i). See “Example 3” below.

An ETA may be used in the definition of MTIME(i).

PREDPP sets all indicator variables to 0 at start of an IR and at RESET events.

MNEXT(i)=1 when PREDPP calls PK prior to the advance to MTIME(i), and remains 1 during the advance from the previous time to MTIME(i). Otherwise, MNEXT(i)=0. (The previous time may be an event time, non-event time, or model event time).

MPAST(i)=0 until the call to PK subsequent to the one for which MNEXT(i)=1. At that call MPAST(i)=1. It then retains this value, unless MTIME is redefined, in which case MPAST will be appropriately redefined as another step function.

Another related variable is MNOW. MNOW=i if this is a call to PK such that MNEXT(i)=1 for some i. MNOW=0 otherwise.

Model event time MTIME is not a dose-related feature and has no effect on steady-state doses. See Note 3 of chapter V.

III.B. Examples

Example of indicator variable values

Suppose there are three MTIMES, and MTIME(1)<MTIME(2)<MTIME(3).

MTIME(1)                  MTIME(2)              MTIME(3)
............|(.............|(.............|(................time axis...

MNEXT(1)=1                MNEXT(2)=1            MNEXT(3)=1
                          MPAST(1)=1            MPAST(2)=1
                                                MPAST(3)=1

Another way to draw the same picture:

MTIME(1)                  MTIME(2)              MTIME(3)
............|..............|..............|..................

011111111111100000000000000000000000000000000000000 MNEXT(1)
000000000000011111111111111100000000000000000000000 MNEXT(2)
000000000000000000000000000011111111111111100000000 MNEXT(3)
000000000000011111111111111111111111111111111111111 MPAST(1)
000000000000000000000000000011111111111111111111111 MPAST(2)
000000000000000000000000000000000000000000011111111 MPAST(3)

If the relationship MTIME(1)<MTIME(2)<MTIME(3) does not hold, imagine the picture above with indices i, j, k (a permutation of 1, 2, 3) such that MTIME(i)<MTIME(j)<MTIME(k).

Event times (both dose and non-dose) can occur anywhere above. See below (Rules for coincidence of times).

If model event times depend on ETAS, so do the values of indicator variables. However, the effect of the model event times or of the indicator variables on the eta-derivatives of the state vector is handled by PREDPP in such a way that the indicator variables are not random variables, and their eta-derivatives are always 0.

Example 1. Flexible input

There is an infusion with modeled rate. The rate is 400*EXP(ETA(1)) from time 0 to 1.5. The rate is 300*EXP(ETA(2)) from time 1.5 to 2.5. The rate is 200*EXP(ETA(3)) from time 2.5 till end of infusion.

Code in $PK:

MTIME(1)=1.5
MTIME(2)=2.5
R1= 400*EXP(ETA(1))*(1-MPAST(1))
R1=R1+300*EXP(ETA(2))*(MPAST(1)-MPAST(2))
R1=R1+200*EXP(ETA(3))*MPAST(2)
WRITE (82,*) TIME,MNOW,MNEXT(1),MNEXT(2),MPAST(1),MPAST(2),R1,TSTATE

Events at times 0, 1, 3.

Values in fort.82 are:

TIME MNOW MNEXT1 MNEXT2 MPAST1 MPAST2     R1 TSTATE
0.00 0.00 0.00   0.00   0.00   0.00   400.00 0.00
1.00 0.00 0.00   0.00   0.00   0.00   400.00 0.00
3.00 1.00 1.00   0.00   0.00   0.00   400.00 1.00
3.00 2.00 0.00   1.00   1.00   0.00   300.00 1.50
3.00 0.00 0.00   0.00   1.00   1.00   200.00 2.50

Example 2. EHC control stream

Here are three versions of an EHC control stream (EHC stands for Enterohepatic Recycling). All give identical results. The transfer of drug from compt. 4 to 1 is controlled by “FLAG”:

DADT(1)=-KA*A(1)+K41*A(4)*FLAG
DADT(4)=K1G*A(2)-K41*A(4)*FLAG

1) No use of model event times.

Define 2 dummy compts. (5 and 6) and dummy doses with lag times:

CID DOSE TIME DV EVID CMT
101 10 0 0 1 2
101 1 0 0 1 5
101 1 0 0 1 6
101 0 2.0 0 0 2

In $PK:

ALAG5=THETA(8)*EXP(ETA(8))
DELTA=THETA(9)*EXP(ETA(9))
ALAG6=ALAG5+DELTA

1a) Define FLAG in $DES:

$DES
FLAG=0.0
IF (T.GE.ALAG5) FLAG=1
IF (T.GE.ALAG5.AND.DOSTIM.EQ.ALAG5) FLAG=0
IF (T.GE.ALAG5.AND.DOSTIM.EQ.0.AND.TIME.EQ.ALAG5) FLAG=0
IF (T.GE.ALAG6) FLAG=0
IF (T.GE.ALAG6.AND.DOSTIM.EQ.ALAG6) FLAG=1
IF (T.GE.ALAG6.AND.DOSTIM.EQ.0.AND.TIME.EQ.ALAG6) FLAG=1

1b) Define FLAG in $PK:

$PK
IF (NEWIND.LT.2) THEN
   FLAG=0.0
   OLDTIM=0
ENDIF
IF (OLDTIM.EQ.ALAG5) FLAG=1
IF (OLDTIM.EQ.ALAG6) FLAG=0
OLDTIM=DOSTIM

2) Use of model event times.

No dummy compartments, no dummy doses:

CID DOSE TIME DV EVID CMT
101 10 0 0 1 2
101 0 2.0 0 0 2

2a) Define MTIMES and FLAG in $PK:

$PK
MTIME(1)=THETA(8)*EXP(ETA(8))
DELTA=THETA(9)*EXP(ETA(9))
MTIME(2)=MTIME(1)+DELTA
FLAG=MPAST(1)-MPAST(2)

2b) Define FLAG in $DES. Might give better run time than 2a.

In 2a, FLAG is considered by NMTRAN to be an implicit basic PK parameter and hence its eta derivatives (which are always 0) are included in the chain rule calculations. In 2b, FLAG is known not to have derivatives because it is DES-defined.

$PK
MTIME(1)=THETA(8)*EXP(ETA(8))
DELTA=THETA(9)*EXP(ETA(9))
MTIME(2)=MTIME(1)+DELTA

$DES
FLAG=MPAST(1)-MPAST(2)

Example 3. MTIME’s may be changed

Suppose there are events at times 0 and 10 but one wants to advance in increments of 1 with stops at times 1, 2, 3, …, 9, 10. This code does not work:

IF (TIME.EQ.0) THEN
   MTIME(1)=0
ELSE
   MTIME(1)=MTIME(1)+1
ENDIF

Instead, use:

IF (TIME.EQ.0) THEN
   TEMP=0
ELSE
   TEMP=TEMP+1
ENDIF
MTIME(1)=TEMP

Reason: If a basic or additional PK parameter is set conditionally in the $PK block, NM-TRAN inserts statements setting it to 0 so that, if it is not set by the $PK block, its value is 0. In the code that does not work, MTIME(1) does not retain its value from one call to the next. In the code that does work, TEMP (being a PK-defined item that does not depend on etas) does retain its value.

Other Examples of the Use of Model Time Parameters

“Model Time Examples”. Some of the above examples are also found in Guide VIII and in on-line help entry for examples.

“Using MTIME to Model Periodic Discontinuities in $DES” contains important examples of modelling circadian rhythm. In Guide VIII and in on-line help entry for examples. This help entry also describes the following control streams.

• Indirect response model: idr_circexa.ctl in the NONMEM examples directory.
• Step Function model: step_circexa.ctl in the NONMEM examples directory.
• See ehc.exa (Enterhepatic circulation example). This gives a fully worked-out example for EHC. It cannot be used with Steady-State doses. This example also gives an alternative model, a smooth step model using Hill terms in a sigmoid emax model. This may be used with Steady-State doses.

III.C. Rules for Coincidence of Times

(1) MTIME(i) = MTIME(j)

There are two calls to PK. Assuming i<j, the calls are described by:

Call	MNEXT(i)	MNEXT(j)	MPAST(i)	MPAST(j)	MNOW
#1	1	0	0	0	1
#2	0	1	1	0	2

(2) MTIME(i) = ALAGj

MNEXT(i) and DOSTIM are set on separate calls to PK. First, MNEXT(i). Then, DOSTIM.

E.g., code in $PK:

ALAG1=.5
MTIME(1)=ALAG1
WRITE (82,*) TIME,MNOW,MNEXT(1),MPAST(1),DOSTIM,TSTATE

Events at times 0, 1, 3. Values in fort.82:

TIME MNOW MNEXT MPAST DOSTIM TSTATE
0.00 0.00 0.00  0.00  0.00   0
1.00 1.00 1.00  0.00  0.00   0
1.00 0.00 0.00  1.00  0.50   .5
1.00 0.00 0.00  1.00  0.00   .5
3.00 0.00 0.00  1.00  0.00   1.0

(3) MTIME(i) = event time

Two calls. First MTIME(1). Then, event time.

E.g., code in $PK:

MTIME(1)=1
WRITE (82,*) TIME,MNOW,MNEXT(1),MPAST(1),TSTATE

Events at times 0, 1, 3. Values in fort.82:

TIME MNOW MNEXT MPAST TSTATE
0.00 0.00 0.00  0.00  0
1.00 1.00 1.00  0.00  0
1.00 0.00 0.00  1.00  1
3.00 0.00 0.00  1.00  1

Compare with:

(3a) ALAGj = event time

Two calls. First, ALAG1. Then, event time.

E.g., code in $PK:

ALAG1=1
WRITE (72,*) TIME,DOSTIM,TSTATE

Events at times 0, 1, 3. Values in fort.72:

TIME DOSTIM TSTATE
0.00 0.00   0
1.00 1.00   0
1.00 0.00   1
3.00 0.00   1

(4) MTIME(i) = additional dose time

MNEXT(i) and DOSTIM are set on separate calls to PK. First, MNEXT(i). Then, DOSTIM.

E.g., code in $PK:

MTIME(1)=.5
WRITE (82,*) TIME,MNOW,MNEXT(1),MPAST(1),DOSTIM,TSTATE

Events at times 0, 1, 3. Dose at time 0 specifies ADDL=1, II=.5. Values in fort.82:

TIME MNOW MNEXT MPAST DOSTIM TSTATE
0.00 0.00 0.00  0.00  0.00   0.00
1.00 1.00 1.00  0.00  0.00   0.00
1.00 0.00 0.00  1.00  0.50   0.50
1.00 0.00 0.00  1.00  0.00   0.50
3.00 0.00 0.00  1.00  0.00   1.00

(5) MTIME(i) = additional dose time = event time

There are three calls to PK. MNEXT(i) and DOSTIM are set on separate calls to PK. First, MNEXT(i). Then, DOSTIM. Finally, a call at event time.

E.g., code in $PK:

MTIME(1)=1
WRITE (82,*) TIME,MNOW,MNEXT(1),MPAST(1),DOSTIM,TSTATE

Events at times 0, 1, 3. Dose at time 0 specifies ADDL=1, II=1. Values in fort.82:

TIME MNOW MNEXT MPAST DOSTIM TSTATE
0.00 0.00 0.00  0.00  0.00   0.00
1.00 1.00 1.00  0.00  0.00   0.00
1.00 0.00 0.00  1.00  1.00   1.00
1.00 0.00 0.00  1.00  0.00   1.00
3.00 0.00 0.00  1.00  0.00   1.00

Compare with:

(5a) MTIME(i) = lagged dose time = event time

There are three calls to PK. MNEXT(i) and DOSTIM are set on separate calls to PK. First, MNEXT(i). Then, DOSTIM. Finally, a call at event time.

E.g., code in $PK:

MTIME(1)=1
ALAG1=1
WRITE (82,*) TIME,MNOW,MNEXT(1),MPAST(1),DOSTIM,TSTATE

Events at times 0, 1, 3. Dose at time 0 specifies ADDL=1, II=1. Values in fort.82:

TIME MNOW MNEXT MPAST DOSTIM TSTATE
0.00 0.00 0.00  0.00  0.00   0.00
1.00 1.00 1.00  0.00  0.00   0.00
1.00 0.00 0.00  1.00  1.00   1.00
1.00 0.00 0.00  1.00  0.00   1.00
3.00 0.00 0.00  1.00  2.00   1.00
3.00 0.00 0.00  1.00  0.00   2.00

III.D. Implementation Details

MTIME parameters are defined in a MODULE. The declaration is:

MODULE PKERR_REAL
 USE SIZES, ONLY : DPSIZE,PCT
 REAL(KIND=DPSIZE) :: MTIME(PCT)
END MODULE PKERR_REAL

where the size is given in SIZES.f90 by PARAMETER (PCT=30) and PCT can be changed with the $SIZES record.

In FSUBS, the generated PK routine contains:

USE PKERR_REAL,ONLY: MTIME

PK contains code to copy MTIME parameters and their eta derivatives to GG. If there are gaps in the sequence of MTIME parameters, the corresponding positions in GG are left as-is. Since GG is initialized to 0 before calls to PK, the effective value of omitted MTIME parameters is 0.

The MTIME array is never listed in NMPRD4. Therefore, if an MTIME parameter is to be displayed in a table or scatter it must be copied to a PK-defined variable in NMPRD4. E.g.,

MTIME(1)=...
MT1=MTIME(1)
 ...
$TABLE MT1

Another reason for defining MTIME in a MODULE is to make it available to the ERROR routine via the same USE statement. All other basic and additional parameters are available to ERROR only if listed in NMPRD4. If MTIME is defined in ERROR, then the values of the MTIME parameters are available as well. However, eta derivatives are not included in MTIME, only in GG. (This is done to save space and compute time.) Hence the ERROR routine should not use MTIME parameters in such a way as to influence the value of Y if those parameters have eta derivatives. NM-TRAN gives an appropriate warning if an MTIME variable is used in the $ERROR block.

Values of MTIME may be used in DES and AES, where they are obtained from the subroutine argument P (which is really GG) and may have eta derivatives.

Compact Arrays in DES

IV.A. Full Arrays and the Chain Rule

With general non-linear models ADVAN6, ADVAN8, ADVAN9, ADVAN13, ADVAN14, and ADVAN16, and ADVAN18 a DES subroutine is needed.

SUBROUTINE DES (A,P,T,DADT,IR,DA,DP,DT)
DIMENSION :: A(*),P(*),DADT(*),DA(IR,*),DP(IR,*),DT(*)

If DES sets IDEFD(2)=1, the arrays are full.

As explained in Chapter VI Section C:

• The partial of DADT(i) wrt A(k) is stored in DA(i,k).
• The partial of DADT(i) wrt P(k) is stored in DP(i,k).
• The partial of DADT(i) wrt T is stored in DT(i).

These arguments are used by PREDPP subroutines such as FCN1 as follows.

The state vector consists of compartment amounts \(A_i\) and their eta partials \(\partial A_i / \partial \eta_j\).

The ADVAN routines must integrate over time \(t\) both of the following, for each compartment \(i\) and eta \(j\):

\[\frac{dA_i}{dt} \quad \text{and} \quad \frac{d}{dt} \frac{\partial A_i}{\partial \eta_j}\]

The latter is not known to PREDPP, but can be computed by changing the order of differentiation:

\[\frac{d}{dt} \frac{\partial A_i}{\partial \eta_j} = \frac{\partial}{\partial \eta_j} \frac{dA_i}{dt}\]

The total derivative \(\partial/\partial \eta_j (dA_i/dt)\) can be obtained by the chain rule by adding the sum over all compartments \(A_k\) and parameters \(P_k\) and time \(T\):

\[ \begin{aligned} &\sum_k \frac{d}{dA_k} \left(\frac{dA_i}{dt}\right) \cdot \frac{\partial A_k}{\partial \eta_j} \\ &\quad + \sum_k \frac{d}{dP_k} \left(\frac{dA_i}{dt}\right) \cdot \frac{\partial P_k}{\partial \eta_j} \\ &\quad + \frac{d}{dT} \left(\frac{dA_i}{dt}\right) \cdot \frac{\partial T}{\partial \eta_j} \end{aligned} \]

But \(d/dA_k(dA_i/dt)\) is simply DA(i,k) and \(d/dP_k(dA_i/dt)\) is DP(i,k) and \(d/dT(dA_i/dt)\) is DT(i).

The derivative \(d/dt(\partial A_i/\partial \eta_j)\) can thus be computed as the sum of quantities that are available to the FCN routine:

\[ \sum_k \mathit{DA}(i, k) \cdot \frac{\partial A_k}{\partial \eta_j} + \sum_k \mathit{DP}(i, k) \cdot \frac{\partial P_k}{\partial \eta_j} + \mathit{DT}(i) \cdot \frac{\partial T}{\partial \eta_j} \]

This can be extended easily to second derivatives.

Note:

1. For every PK parameter (i.e., every variable defined in PK and used in DES) there is a loop over all the etas. If a given PK parameter does not have partials with respect to etas, the loop does a lot of multiplying by zero’s and adding of zero’s. This wastes time. It is better simply to define it in $DES.
2. There is no third dimension of DA or DP or DT for the second partials (e.g., partials of DADT with respect to A(k),A(n)), or cross-partials (e.g., partials of DADT with respect to A(k),P(n)). Therefore the Laplacian method cannot be used with the FULL arrays.
3. FCN1 always loops over all elements of DA(i,j), DP(i,j) and DT(i). With a large sparse system of differential equations, many of these elements are zero because they are never computed by DES. With compact arrays (below), the elements that are always zero are not included in the computation.

IV.B. Compact Arrays

If DES sets IDEFD(2)=0, the arrays are compact. (This is the default). FCN1 uses only those derivatives that are actually computed by DES. In addition, second partials (Laplace method) can be computed by DES.

DES computes all first and second partials and cross-partials of DADT. These may include (for all compartments and parameters defined in the model):

In DA:

• the first partials of DADT(i) wrt A(j)
• the second partials of DADT(i) wrt A(j) A(k)
• the cross partials of DADT(i) wrt A(j) P(k)

In DP:

• the first partials of DADT(i) wrt P(j)
• the second partials of DADT(i) wrt P(j) P(k)

In DT:

• the first partials of DADT(i) wrt T
• the second partials of DADT(i) wrt T T
• the cross partials of DADT(i) wrt T A(j)
• the cross partials of DADT(i) wrt T P(j)

At ICALL=1, DES uses the first row DA(*,1) to describe each position in DA that will contain a value and describes (“maps”) which A and eta it corresponds to (and which P, in the case of second partials). At other values of ICALL, DES puts the values themselves in DA(*,1). Similarly for DP and DT. Thus DA, DP, DT are linear (e.g., only the first row DA(*,1) is used for both the ICALL=1 “map” and the ICALL=2 value.)

The mapping uses two constants:

• PG = max. no. of basic and additional PK parameters; may be set by $SIZES PG
• PC = max. no. of compartments + 1; may be set by $SIZES PC

From these are computed:

MM=MAX(100,PC-1+PG+10)
MMSQ=MM*MM

For DA, the mapping is in position p is:

DA(p,1)=DA(i*MMSQ+j*MM+k,1)

Where:

• i tells which DADT(i)
• j tells which A(j)
• k=0 if this is a first partial; k>0 if this is a second partial.
  • if k < PC, k tells which A(k)
  • if k > PC, k tells which P(n), where n=k-PC.

Some derivatives may happen to be zero for some values of A’s and P’s and etas, but presumably they are non-zero at other values. Position p in DA will be used to store the value at ICALL=2. The DA values are terminated by DA(..)=0.

Similarly for partials of DADT(i) wrt P:

DP(i*MMSQ+j*MM+k,1)

Where:

• j tells which P(j)
• k=0 if this is a first partial; k>0 if this is a second partial.
• if k > 0, k tells which P(k).

Similarly for partials of DADT(i) wrt T:

DT(i*MMSQ+k,1)

Where:

• k=0 if this is a first partial; k>0 if this is a second partial.
• if 0 < k < PC, this is partial wrt T, A(k)
• if k > PC and k < PC+PG, this is partial wrt T, P(n) where n=k-PC
• if k = PC+PG+1, this is second partial wrt T T.

IV.C. Example

An example is given as app4example.ctl. This is not a realistic model, but serves to demonstrate the mapping. The only differential equation is:

DADT(1)=-T*T*A(1)*A(2)*P(1)*P(2)

The example contains:

$SIZES PC=5 PG=10

Then MM=100 and MMSQ=10000 and the mapping given by the right sides can be read easily.

For DA, right sides have the form 00i0j0k. The DA values in DES at ICALL=1 are:

DA( 1,1)=0000010100
DA( 2,1)=0000010200
DA( 3,1)=0000010201
DA( 4,1)=0000010102
DA( 5,1)=0000010106
DA( 6,1)=0000010107
DA( 7,1)=0000010206
DA( 8,1)=0000010207

• 1st partial of DADT(1) wrt A(1) will be stored in DA(1,1)
• 1st partial of DADT(1) wrt A(2) will be stored in DA(2,1)
• 2nd partial of DADT(1) wrt A(2) A(1) will be stored in DA(3,1)
• 2nd partial of DADT(1) wrt A(1) A(2) will be stored in DA(4,1)
• 2nd partial of DADT(1) wrt A(1) P(1) will be stored in DA(5,1)
• 2nd partial of DADT(1) wrt A(1) P(2) will be stored in DA(6,1)
• 2nd partial of DADT(1) wrt A(2) P(1) will be stored in DA(7,1)
• 2nd partial of DADT(1) wrt A(2) P(2) will be stored in DA(8,1)

For DP, right sides have the form 00i0j0k. The DP values in DES at ICALL=1 are:

DP( 1,1)=0000010100
DP( 2,1)=0000010200
DP( 3,1)=0000010201
DP( 4,1)=0000010102
DP( 5,1)=0

• 1st partial of DADT(1) wrt P(1) will be stored in DP(1,1)
• 1st partial of DADT(1) wrt P(2) will be stored in DP(2,1)
• 2nd partial of DADT(1) wrt P(2) P(1) will be stored in DP(3,1)
• 2nd partial of DADT(1) wrt P(1) P(2) will be stored in DP(4,1)

For DT, right sides have the form: 00i000k. The DT values in DES at ICALL=1 are:

DT( 1)=0000010000
DT( 2)=0000010001
DT( 3)=0000010002
DT( 4)=0000010016
DT( 5)=0000010006
DT( 6)=0000010007
DT( 7)=0

• 1st partial of DADT(1) wrt T will be stored in DT(1)
• 2nd partial of DADT(1) wrt T A(1) will be stored in DT(2)
• 2nd partial of DADT(1) wrt T A(2) will be stored in DT(3)
• 2nd partial of DADT(1) wrt T T will be stored in DT(4)
• 2nd partial of DADT(1) wrt T P(1) will be stored in DT(5)
• 2nd partial of DADT(1) wrt T P(2) will be stored in DT(6)

Here is the app4example.ctl:

$SIZES PC=5 PG=10

$PROB DEMONSTRATE COMPACT DES MAPPING

$INPUT      ID DOSE=AMT TIME CP=DV WT

$DATA       THEOPP RECS=1

$SUBROUTINES ADVAN6 TOL=5
$MODEL NCM=2

$PK
     P(1)=THETA(1)*EXP(ETA(1))
     P(2)=THETA(2)*EXP(ETA(2))

$DES
     DADT(1)=-T*T*A(1)*A(2)*P(1)*P(2)

$ERROR
     Y=F+EPS(1)

$THETA 1 2

$OMEGA 1 2

$SIGMA  1