60 Examples
This appendix reproduces the control-stream examples from the original PDF (pages EX-1 through EX-159); they illustrate specific NONMEM features referenced from the Detailed Descriptions.
1 Abbreviated Function Examples
Here is an example of the use and construction of a function of both time and the first-order elimination rate constant of a simple PK one-compartment model, to compute the amount in the system at a given time.
(See Section 53.2)
$PRED
...
VECTRA(1) = THETA(1) * EXP(ETA(1))
VECTRA(2) = TIME
...
A = DOSE * FUNCA(VECTRA)
...With the $ABBREVIATED REPLACE feature, the function name can be changed to one that is perhaps easier to understand.
$ABBR REPLACE ONE_COMPARTMENT_LINEAR=FUNCA
...
$PRED
...
VECTRA(1) = THETA(1) * EXP(ETA(1))
VECTRA(2) = TIME
...
A = DOSE * ONE_COMPARTMENT_LINEAR(VECTRA)
...Either way, one needs to supply the function.
FUNCTION FUNCA(X,X1,X2)
! Implements the One Compartment Linear model
USE SIZES, ONLY: DPSIZE
REAL(KIND=DPSIZE), INTENT(IN) :: X
REAL(KIND=DPSIZE), INTENT(IN OUT) :: X1,X2
REAL(KIND=DPSIZE) :: FUNCA
DIMENSION :: X(9),X1(9),X2(9,9)
REAL(KIND=DPSIZE) :: EXPT
! THE FUNCTION ITSELF
EXPT=EXP(-X(1)*X(2))
FUNCA=EXPT
! 1ST. PARTIALS
X1(1)=-EXPT*X(2)
X1(2)=-EXPT*X(1)
! 2ND. PARTIALS
X2(1,1)=EXPT*X(2)*X(2)
X2(1,2)=EXPT*X(1)*X(2)
X2(2,1)=EXPT*X(1)*X(2)
X2(2,2)=EXPT*X(1)*X(1)
RETURN
ENDSuppose the above code is in file ONE_COMPARTMENT_LINEAR.f90. Then the control stream should also include
$SUBROUTINES OTHER=ONE_COMPARTMENT_LINEAR.f90Multiple functions can be used with the $ABBREVIATED REPLACE feature. Suppose two functions are needed, called (for example) RED and GREEN. The control stream might contain the following fragment:
$SUBR OTHER=red.f90 OTHER=green.f90
...
$ABBR REPLACE RED=FUNCA
$ABBR REPLACE GREEN=FUNCB
...
$PRED
...
VECTRA(1)=THETA(1)*EXP(ETA(1))
VECTRA(2)=TIME
A=DOSE*RED(VECTRA)
B=DOSE*GREEN(VECTRA)
...The file red.f90 contains code for FUNCA, and the file green.f90 contains code for FUNCB.
The $ABBREVIATED REPLACE feature could also be used to provide more meaningful names for reserved vectors. For example, suppose the arguments of the functions are different vectors.
$ABBR REPLACE RED=FUNCA, REDARG=VECTRA
$ABBR REPLACE GREEN=FUNCB, GREENARG=VECTRB
...
$PRED
...
REDARG(1)=...
REDARG(2)=...
GREENARG(1)=...
GREENARG(2)=...
A=DOSE*RED(REDARG)
B=DOSE*GREEN(GREENARG)
...2 Bayes Example 1
This is example1.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X SDIX SDSX
$DATA example1.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
; The thetas are MU modeled.
; Best that there is a linear relationship between THETAs and Mus
; The linear MU modeling of THETAS allows them to be efficiently
; Gibbs sampled.
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
; Initial values of THETA
$THETA
(0.001, 2.0) ; [LN(CL)]
(0.001, 2.0) ; [LN(V1)]
(0.001, 2.0) ; [LN(Q)]
(0.001, 2.0) ; [LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA
$SIGMA
( 0.6 ) ; [P]
;Prior information is important for MCMC Bayesian analysis,
;not necessary for maximization methods
;Note the syntax used for defining priors that is available
;as of NONMEM 7.3
$PRIOR NWPRI
; Prior information of THETAS
$THETAP (2.0 FIX)X4
; Variance to prior information of THETAS.
; Because variances are very large, this means that the prior
; information to the THETAS is highly uninformative.
$THETAPV BLOCK(4) FIX VALUES(10000,0.0)
; Prior information to the OMEGAS.
$OMEGAP BLOCK(4) FIX VALUES(0.2,0.0)
; Degrees of freedom to prior OMEGA matrix.
; Because degrees of freedom is very low, equal to the
; the dimension of the prior OMEGA, this means that the
; prior information to the OMEGAS is highly uninformative
$OMEGAPD (4 FIX)
; Prior information to the SIGMAS
$SIGMAP 0.06 FIX
; Degrees of freedom to prior SIGMA matrix.
; Because degrees of freedom is very low, equal to the
; the dimension of the prior SIGMA, this means that the
; prior information to the SIGMA is highly uninformative
$SIGMAPD (1 FIX)
; The first analysis is iterative two-stage,
; maximum of 500 iterations (NITER), iteration results
; are printed every 5 iterations, gradient precision (SIGL) is 4.
; Termination is tested on all of
; the population parameters (CTYPE=3),
; and for less then 2 significant digits change (NSIG).
; Prior information is not necessary for ITS, so NOPRIOR=1.
; The intermediate and final results of the ITS method will be
; recoded in row/column format in example1.ext
$EST METHOD=ITS MAPITER=0 INTERACTION FILE=example1.ext NITER=500
PRINT=5 NOABORT SIGL=4 CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=1 NSIG=2
; The results of ITS are used as the initial values for the
; SAEM method. A maximum of 3000 ; stochastic iterations (NBURN)
; is requested, but may end early if statistical test determines
; that variations in all parameters is stationary
; (note that any settings from the previous $EST
; carries over to the next $EST statement, within a $PROB).
; The SAEM is a Monte Carlo process,
; so setting the SEED assures repeatability of results.
; Each iteration obtains only 2 Monte Carlo samples ISAMPLE),
; so they are very fast.
; But many iterations are needed, so PRINT only
; every 100th iteration.
; After the stochastic phase, 500 accumulation iterations will be
; Performed (NITER), to obtain good parameters estimates with
; little stochastic noise.
; As a new FILE has not been given, the SAEM results will append to
; example1.ext.
$EST METHOD=SAEM INTERACTION NBURN=3000 NITER=500 PRINT=100
SEED=1556678 ISAMPLE=2
; After the SAEM method, obtain good estimates of the marginal
; density (objective function),
; along with good estimates of the standard errors.
; This is best done with importance sampling ; (IMP),
; performing the expectation step only (EONLY=1), so that
; final population parameters remain at the final SAEM result.
; Five iterations (NITER) should allow the importance sampling
; proposal density to become stationary.
; This is observed by the objective function settling
; to a particular value (with some stochastic noise).
; By using 3000 Monte Carlo samples
; (ISAMPLE), this assures a precise assessment of standard errors.
$EST METHOD=IMP INTERACTION EONLY=1 NITER=5 ISAMPLE=3000 PRINT=1
SIGL=8 NOPRIOR=1
; The Bayesian analysis is performed.
; While 10000 burn-in iterations are requested as a maximum,
; because the termination test is on (CTYPE<>0, set at the
; first $EST statement), and because the initial parameters are at
; the SAEM result, which is the maximum likelihood position,
; the analysis should settle down to a stationary distribution in
; several hundred iterations.
; Prior information is also used to facilitate Bayesian analysis.
; The individual Bayesian iteration results are important,
; and may be need for post-processing analysis.
; So specify a separate FILE for the Bayesian analysis.
$EST METHOD=BAYES INTERACTION FILE=example1.txt NBURN=10000
NITER=10000 PRINT=100 NOPRIOR=0
; Just for old-times sake, let's see what the traditional
; FOCE method will give us.
; And, remember to introduce a new FILE, so its results won't
; append to our Bayesian FILE.
; Appending to example1.ext with the EM methods is fine.
$EST METHOD=COND INTERACTION MAXEVAL=9999 NSIG=3 SIGL=10
PRINT=5 NOABORT NOPRIOR=1
FILE=example1.ext
; Time for the standard error results.
; You may request a more precise gradient precision (SIGL)
; that differed from that used during estimation.
$COV MATRIX=R PRINT=E UNCONDITIONAL SIGL=12
; Print out results in tables. Include some of the new weighted
; residual types
$TABLE ID TIME PRED RES WRES CPRED CWRES EPRED ERES EWRES NOAPPEND
ONEHEADER FILE=example1.TAB NOPRINT
$TABLE ID CL V1 Q V2 FIRSTONLY NOAPPEND NOPRINT FILE=example1.PAR
$TABLE ID ETA1 ETA2 ETA3 ETA4 FIRSTONLY NOAPPEND
NOPRINT FILE=example1.ETA3 Bayes Example 1B
This is example1B.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 1B (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X SDIX SDSX
$DATA example1b.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
; The thetas are MU modeled.
; Best that there is a linear relationship between THETAs and Mus
; The linear MU modeling of THETAS allows them to be efficiently
; Gibbs sampled.
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
; Initial values of THETA
$THETA
(0.001, 2.0) ; [LN(CL)]
(0.001, 2.0) ; [LN(V1)]
(0.001, 2.0) ; [LN(Q)]
(0.001, 2.0) ; [LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA
$SIGMA
( 0.6 ) ; [P]
;NTHETA=number of Thetas to be estimated
;NETA=number of Etas to be estimated
; (and to be described by NETAxBETA OMEGA matrix)
;NTHP=number of thetas which have a prior
;NETP=number of Omegas with prior
;Prior information is important for MCMC Bayesian analysis,
; not necessary for maximization methods
$PRIOR NWPRI NTHETA=4, NETA=4, NTHP=4, NETP=4
; Prior information of THETAS
$THETAP (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)
; Variance to prior information of THETAS.
; Because variances are very large, this
; means that prior information to the THETAS is highly uninformative.
$THETAPV BLOCK(4)
10000 FIX
0.00 10000
0.00 0.00 10000
0.00 0.00 0.0 10000
; Prior information to the OMEGAS.
$OMEGAP BLOCK(4)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
; Degrees of freedom to prior OMEGA matrix.
; Because degrees of freedom is very low, equal to the
; the dimension of the prior OMEGA, this means that the
; prior information to the OMEGAS is highly uninformative
$OMEGAPD (4 FIX)
; The first analysis is iterative two-stage,
; maximum of 500 iterations (NITER), iteration results
; are printed every 5 iterations, gradient precision (SIGL) is 4.
; Termination is tested on all of
; the population parameters (CTYPE=3),
; and for less then 2 significant digits change (NSIG).
; Prior information is not necessary for ITS, so NOPRIOR=1.
; The intermediate and final results
; of the ITS method will be recoded in row/column format in
; example1b.ext
$EST METHOD=ITS MAPITER=0 INTERACTION FILE=example1b.ext NITER=500
PRINT=5 CINTERVAL=5 NOABORT SIGL=4 CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=1 NSIG=2
; The results of ITS are used as the initial values for the SAEM
; method. A maximum of 3000 ; stochastic iterations (NBURN) is
; requested, but may end early if statistical test determines
; that variations in all parameters is stationary
; (note that any settings from the previous $EST
; carries over to the next $EST statement, within a $PROB).
; The SAEM is a Monte Carlo process,
; so setting the SEED assures repeatability of results.
; Each iteration obtains only 2 Monte Carlo samples ISAMPLE),
; so they are very fast.
; But many iterations are needed, so PRINT only
; every 100th iteration.
; After the stochastic phase, 500 accumulation iterations will be
; Performed (NITER), to obtain good parameters estimates with
; little stochastic noise.
; As a new FILE has not been given, the SAEM results will append to
; example1b.ext.
$EST METHOD=SAEM INTERACTION NBURN=3000 NITER=500
PRINT=250 CINTERVAL=100
SEED=1556678 ISAMPLE=2
; After the SAEM method, obtain good estimates of the marginal
; density (objective function),
; along with good estimates of the standard errors.
; This is best done with importance sampling ; (IMP),
; performing the expectation step only (EONLY=1), so that
; final population parameters remain at the final SAEM result.
; Five iterations (NITER) should allow the importance sampling
; proposal density to become stationary.
; This is observed by the objective function settling
; to a particular value (with some stochastic noise).
; By using 3000 Monte Carlo samples
; (ISAMPLE), this assures a precise assessment of standard errors.
$EST METHOD=IMP INTERACTION EONLY=1 MAPITER=0 NITER=5
ISAMPLE=3000 PRINT=1 SIGL=8 NOPRIOR=1
$EST METHOD=DIRECT INTERACTION EONLY=1 NITER=3 ISAMPLE=3000
PRINT=1 SIGL=8 NOPRIOR=1
; The Bayesian analysis is performed.
; While 10000 burn-in iterations are requested as a maximum,
; because the termination test is on (CTYPE<>0, set at the
; first $EST statement), and because the initial parameters are at
; the SAEM result, which is the maximum likelihood position,
; the analysis should settle down to a stationary distribution in
; several hundred iterations.
; Prior information is also used to facilitate Bayesian analysis.
; The individual Bayesian iteration results are important,
; and may be need for post-processing analysis.
; So specify a separate FILE for the Bayesian analysis.
$EST METHOD=BAYES INTERACTION FILE=example1b.txt NBURN=10000
NITER=10000 PRINT=1000 CINTERVAL=100 NOPRIOR=0
; Just for old-times sake, let's see what the traditional
; FOCE method will give us.
; And, remember to introduce a new FILE, so its results won't
; append to our Bayesian FILE.
; Appending to example1b.ext with the EM methods is fine.
$EST METHOD=COND INTERACTION MAXEVAL=9999 NSIG=3 SIGL=10
PRINT=5 NOABORT NOPRIOR=1
FILE=example1b.ext
; Time for the standard error results.
; You may request a more precise gradient precision (SIGL)
; that differed from that used during estimation.
$COV MATRIX=R PRINT=E UNCONDITIONAL SIGL=12
; Print out results in tables. Include some of the new weighted
; residual types
$TABLE ID TIME PRED RES WRES CPRED CWRES EPRED ERES EWRES NOAPPEND
ONEHEADER FILE=example1b.TAB NOPRINT
$TABLE ID CL V1 Q V2 FIRSTONLY NOAPPEND NOPRINT FILE=example1b.PAR
FORMAT=,1PE15.8
$TABLE ID ETA1 ETA2 ETA3 ETA4 FIRSTONLY NOAPPEND NOPRINT
FILE=example1b.ETA FORMAT=,1PE15.84 Bayes Example 2
This is example2.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two Compartment model with Clearance and
; central volume modeled with covariates age and gender
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example2 (from sampc)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT GNDR AGE
$DATA example2.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
; LCLM=log transformed clearance, male
LCLM=THETA(1)
;LCLF=log transformed clearance, female.
LCLF=THETA(2)
; CLAM=CL age slope, male
CLAM=THETA(3)
; CLAF=CL age slope, female
CLAF=THETA(4)
; LV1M=log transformed V1, male
LV1M=THETA(5)
; LV1F=log transformed V1, female
LV1F=THETA(6)
; V1AM=V1 age slope, male
V1AM=THETA(7)
; V1AF=V1 age slope, female
V1AF=THETA(8)
; LAGE=log transformed age
LAGE=DLOG(AGE)
;Mean of ETA1, the inter-subject deviation of Clearance,
; is ultimately modeled as linear function of THETA(1) to THETA(4).
; Relating thetas to Mus by linear functions is not essential for
; ITS, IMP, or IMPMAP methods, but is very helpful for MCMC methods
; such as SAEM and BAYES.
MU_1=(1.0-GNDR)*(LCLM+LAGE*CLAM) + GNDR*(LCLF+LAGE*CLAF)
; Mean of ETA2, the inter-subject deviation of V1,
; is ultimately modeled as linear function of THETA(5) to THETA(8)
MU_2=(1.0-GNDR)*(LV1M+LAGE*V1AM) + GNDR*(LV1F+LAGE*V1AF)
MU_3=THETA(9)
MU_4=THETA(10)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
CALLFL=0
; Option to model the residual error coefficient in THETA(11),
; rather than in SIGMA.
SDSL=THETA(11)
W=F*SDSL
Y = F + W*EPS(1)
IPRED=F
IWRES=(DV-F)/W
;Initial THETAs
$THETA
( 0.7 ) ; [LCLM]
( 0.7 ) ; [LCLF]
( 2 ) ; [CLAM]
( 2.0 ) ; [CLAF]
( 0.7 ) ; [LV1M]
( 0.7 ) ; [LV1F]
( 2.0 ) ; [V1AM]
( 2.0 ) ; [V1AF]
( 0.7 ) ; [MU_3]
( 0.7 ) ; [MU_4]
( 0.3 ) ; [SDSL]
;Initial OMEGAs
$OMEGA BLOCK(4)
0.5 ;[p]
0.001 ;[f]
0.5 ;[p]
0.001 ;[f]
0.001 ;[f]
0.5 ;[p]
0.001 ;[f]
0.001 ;[f]
0.001 ;[f]
0.5 ;[p]
; SIGMA is 1.0 fixed, serves as unscaled variance for EPS(1).
; THETA(11) takes up the residual error scaling.
$SIGMA
(1.0 FIXED)
;Prior information is important for MCMC Bayesian analysis,
; not necessary for maximization methods
; In this example, only the OMEGAs have a prior distribution,
; the THETAS do not.
; For Bayesian methods, it is most important for at least the
; OMEGAs to have a prior, even an uninformative one,
; to stabilize the analysis. Only if the number of subjects
; exceeds the OMEGA dimension number by at least 100,
; then you may get away without priors on OMEGA for BAYES analysis.
$PRIOR NWPRI
; Prior OMEGA matrix
$OMEGAP BLOCK(4) FIX VALUES(0.01,0.0)
; Degrees of freedom to OMEGA prior matrix:
$OMEGAPD 4 FIX
; The first analysis is iterative two-stage.
; Note that the GRD specification is THETA(11) is a
; Sigma-like parameter. This will allow NONMEM to make
; efficient gradient evaluations for THETA(11), which is useful
; for later IMP,IMPMAP, and SAEM methods, but has no impact on
; ITS and BAYES methods.
$EST METHOD=ITS INTERACTION FILE=example2.ext NITER=1000 NSIG=2
PRINT=5 NOABORT SIGL=8 NOPRIOR=1 CTYPE=3 GRD=TS(11)
; Results of ITS serve as initial parameters for the IMP method.
$EST METHOD=IMP INTERACTION EONLY=0 MAPITER=0 NITER=100 ISAMPLE=300
PRINT=1 SIGL=8
; The results of IMP are used as the initial values for the SAEM method.
$EST METHOD=SAEM NBURN=3000 NITER=2000 PRINT=10 ISAMPLE=2
CTYPE=3 CITER=10 CALPHA=0.05
; After the SAEM method, obtain good estimates of the marginal density
; (objective function),
; along with good estimates of the standard errors.
$EST METHOD=IMP INTERACTION EONLY=1 NITER=5 ISAMPLE=3000
PRINT=1 SIGL=8 SEED=123334
CTYPE=3 CITER=10 CALPHA=0.05
; The Bayesian analysis is performed.
$EST METHOD=BAYES INTERACTION FILE=example2.TXT NBURN=10000
NITER=3000 PRINT=100 NOPRIOR=0
CTYPE=3 CITER=10 CALPHA=0.05
; Just for old-times sake, lets see what the traditional
; FOCE method will give us.
; And, remember to introduce a new FILE, so its results wont
; append to our Bayesian FILE.
$EST METHOD=COND INTERACTION MAXEVAL=9999 FILE=example2.ext NSIG=2
SIGL=14 PRINT=5 NOABORT NOPRIOR=1
$COV MATRIX=R UNCONDITIONAL5 Bayes Example 3
This is example3.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Population Mixture Problem in 1 Compartment model,
; with Volume and rate constant parameters and their inter-subject
; variances modeled from two sub-populations
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example3 (from ad1tr1m2s)
$INPUT C SET ID JID TIME CONC=DV DOSE=AMT RATE EVID MDV CMT VC1 K101 VC2 K102 SIGZ PROB
$DATA example3.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS1
; The mixture model uses THETA(5) as the mixture proportion parameter,
; defining the proportion of subjects in sub-population 1 (P(1),
; and in sub-population 2 (P(2)
$MIX
P(1)=THETA(5)
P(2)=1.0-THETA(5)
NSPOP=2
$PK
; The MUs should always be unconditionally defined, that is,
; they should never be defined in IF/THEN blocks
; THETA(1) models the Volume of sub-population 1
MU_1=THETA(1)
; THETA(2) models the clearance of sub-population 1
MU_2=THETA(2)
; THETA(3) models the Volume of sub-population 2
MU_3=THETA(3)
; THETA(4) models the clearance of sub-population 2
MU_4=THETA(4)
VCM=DEXP(MU_1+ETA(1))
K10M=DEXP(MU_2+ETA(2))
VCF=DEXP(MU_3+ETA(3))
K10F=DEXP(MU_4+ETA(4))
Q=1
IF(MIXNUM.EQ.2) Q=0
V=Q*VCM+(1.0-Q)*VCF
K=Q*K10M+(1.0-Q)*K10F
S1=V
$ERROR
Y = F + F*EPS(1)
; Initial THETAs
$THETA
(-1000.0 4.3 1000.0) ; [MU_1]
(-1000.0 -2.9 1000.0) ; [MU_2]
(-1000.0 4.3 1000.0) ; [MU_3]
(-1000.0 -0.67 1000.0) ; [MU_4]
(0.0001 0.667 0.9999) ; [P(1)]
;Initial OMEGA block 1, for sub-population 1
$OMEGA BLOCK(2)
.04 ;[p]
.01 ; [f]
.027; [p]
;Initial OMEGA block 2, for sub-population 2
$OMEGA BLOCK(2)
.05; [p]
.01; [f]
.06; [p]
$SIGMA
0.01 ;[p]
; Prior information setup for OMEGAS only
$PRIOR NWPRI
; Prior OMEGA block 1. Note that because the OMEGA is separated
; into blocks, so their priors should have the same block design.
$OMEGAP BLOCK(2)
0.05 FIX
0.0 0.05
; Prior OMEGA block 2
$OMEGAP BLOCK(2)
0.05 FIX
0.0 0.05
; Degrees of Freedom defined for Priors.
; One for each OMEGA block defining each sub-popluation
$OMEGAPD (2 FIX) (2 FIX)
$EST METHOD=ITS INTERACTION NITER=20 PRINT=1 NOABORT SIGL=8
FILE=example3.ext CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=1
$EST NBURN=500 NITER=500 METHOD=SAEM INTERACTION PRINT=10 SIGL=6
ISAMPLE=2
$EST METHOD=IMP INTERACTION NITER=5 ISAMPLE=1000 PRINT=1 NOABORT
SIGL=6 EONLY=1 MAPITER=0
$EST METHOD=BAYES INTERACTION NBURN=2000 NITER=1000 PRINT=10
FILE=example3.txt SIGL=8 NOPRIOR=0
$EST MAXEVAL=9999 NSIG=3 SIGL=12 PRINT=1 FILE=example3.ext
METHOD=CONDITIONAL INTERACTION NOABORT
NOPRIOR=1
$COV MATRIX=R UNCONDITIONAL6 Bayes Example 4
This is example4.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Population Mixture Problem in 1 Compartment model,
; with rate constant parameter and its inter-subject variances
; modeled as coming from two sub-populations
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example4 (from ad1tr1m2t)
$INPUT C SET ID JID TIME CONC=DV DOSE=AMT RATE EVID MDV CMT VC1 K101 VC2 K102 SIGZ PROB
$DATA example4.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS1
$MIX
P(1)=THETA(4)
P(2)=1.0-THETA(4)
NSPOP=2
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
V=DEXP(MU_1+ETA(1))
K10M=DEXP(MU_2+ETA(2))
K10F=DEXP(MU_3+ETA(3))
Q=1
IF(MIXNUM.EQ.2) Q=0
K=Q*K10M+(1.0-Q)*K10F
S1=V
$ERROR
Y = F + F*EPS(1)
$THETA
(-1000.0 4.3 1000.0) ; [MU_1]
(-1000.0 -2.9 1000.0) ; [MU_2]
(-1000.0 -0.67 1000.0) ; [MU_3]
(0.0001 0.667 0.9999) ; [P(1)]
$OMEGA BLOCK(3)
.04 ;[p]
0.01 ;[f]
.027 ;[p]
0.01 ;[f]
0.001 ;[f]
.06 ;[p]
$SIGMA
0.01 ;[p]
; Prior information setup for OMEGAS only
$PRIOR NWPRI
; Prior OMEGA
$OMEGAP BLOCK(3)
0.05 FIX
0.0 0.05
0.0 0.0 0.05
; Degrees of Freedom defined for Priors.
$OMEGAPD (3 FIX)
$EST METHOD=ITS INTERACTION NITER=30 PRINT=5 NOABORT SIGL=6
FILE=example4.ext NOPRIOR=1 CTYPE=3 CITER=10 CALPHA=0.05
$EST METHOD=IMP INTERACTION NITER=20 ISAMPLE=300 PRINT=1
NOABORT SIGL=6 NOPRIOR=1
$EST NBURN=500 NITER=500 METHOD=SAEM INTERACTION PRINT=10 SIGL=6
ISAMPLE=2 NOPRIOR=1 MAPITER=0
$EST METHOD=IMP INTERACTION EONLY=1 NITER=20 ISAMPLE=3000 PRINT=1
NOABORT SIGL=6 NOPRIOR=1
$EST METHOD=BAYES INTERACTION NBURN=2000 NITER=5000 PRINT=10
FILE=example4.txt SIGL=6 NOPRIOR=0
$EST MAXEVAL=9999 NSIG=3 SIGL=12 PRINT=1
METHOD=CONDITIONAL INTERACTION
NOABORT FILE=example4.ext NOPRIOR=1
$COV MATRIX=R UNCONDITIONAL SIGL=107 Bayes Example 5
This is example5.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Population Mixture Problem in 1 Compartment model,
; with rate constant parameter mean modeled for two
; sub-populations, but its inter-subject variance is the same in
; both sub-populations Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example5 (from ad1tr1m4t)
$INPUT C SET ID JID TIME CONC=DV DOSE=AMT RATE EVID MDV CMT VC1 K101 VC2 K102 SIGZ PROB
$DATA example5.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS1
$MIX
P(1)=THETA(4)
P(2)=1.0-THETA(4)
NSPOP=2
$PK
Q=1
IF(MIXNUM.EQ.2) Q=0
MU_1=THETA(1)
; Note that MU_2 can be modeled as THETA(2) or THETA(3),
; depending on the MIXNUM value.
; Also, we are avoiding IF/THEN blocks.
MU_2=Q*THETA(2)+(1.0-Q)*THETA(3)
V=DEXP(MU_1+ETA(1))
K=DEXP(MU_2+ETA(2))
S1=V
$ERROR
Y = F + F*EPS(1)
$THETA
(-1000.0 4.3 1000.0) ; [MU_1]
(-1000.0 -2.9 1000.0) ; [MU_2-1]
(-1000.0 -0.67 1000.0) ; [MU_2-2]
(0.0001 0.667 0.9999) ; [P(1)]
$OMEGA BLOCK(2)
0.04 ;[p]
0.01 ;[f]
0.04 ;[p]
$SIGMA
0.01 ;[p]
$EST METHOD=ITS INTERACTION NITER=100 PRINT=1 NOABORT SIGL=8
FILE=example5.ext CTYPE=3
$EST METHOD=IMPMAP INTERACTION NITER=20 ISAMPLE=300 PRINT=1
NOABORT SIGL=8
$EST METHOD=IMP INTERACTION NITER=20 ISAMPLE=1000 PRINT=1
NOABORT SIGL=6
$EST NBURN=500 NITER=500 METHOD=SAEM INTERACTION PRINT=10
SIGL=6 ISAMPLE=2
$EST METHOD=IMP INTERACTION NITER=5 ISAMPLE=1000 PRINT=1
NOABORT SIGL=6 EONLY=1 MAPITER=0
$EST METHOD=BAYES INTERACTION NBURN=2000 NITER=5000 PRINT=10
FILE=example5.txt SIGL=8
$EST MAXEVAL=9999 NSIG=2 SIGL=8 PRINT=10 FILE=example5.ext
METHOD=CONDITIONAL INTERACTION NOABORT
$COV MATRIX=R8 Bayes Example 6
This is example6.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Receptor Mediated Clearance model with Dynamic Change
; in Receptors
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example6 (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C
; The new numerical integration solver is used, although ADVAN=9
; is also efficient for this problem.
$SUBROUTINES ADVAN13 TRANS1 TOL=4
$MODEL NCOMPARTMENTS=3
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
MU_6=THETA(6)
MU_7=THETA(7)
MU_8=THETA(8)
VC=EXP(MU_1+ETA(1))
K10=EXP(MU_2+ETA(2))
K12=EXP(MU_3+ETA(3))
K21=EXP(MU_4+ETA(4))
VM=EXP(MU_5+ETA(5))
KMC=EXP(MU_6+ETA(6))
K03=EXP(MU_7+ETA(7))
K30=EXP(MU_8+ETA(8))
S3=VC
S1=VC
KM=KMC*S1
F3=K03/K30
$DES
DADT(1) = -(K10+K12)*A(1) + K21*A(2) - VM*A(1)*A(3)/(A(1)+KM)
DADT(2) = K12*A(1) - K21*A(2)
DADT(3) = -VM*A(1)*A(3)/(A(1)+KM) - K30*A(3) + K03
$ERROR
CALLFL=0
ETYPE=1
IF(CMT.NE.1) ETYPE=0
IPRED=F
Y = F + F*ETYPE*EPS(1) + F*(1.0-ETYPE)*EPS(2)
$THETA
;Initial Thetas
(4.0) ; [MU_1]
(-2.1) ; [MU_2]
(0.7) ; [MU_3]
(-0.17) ; [MU_4]
(2.2) ; [MU_5]
(0.14) ; [MU_6]
(3.7) ; [MU_7]
(-0.7) ; [MU_8]
;Initial Omegas
$OMEGA BLOCK(8)
0.2 ;[p]
-0.0043 ;[f]
0.2 ;[p]
0.0048 ;[f]
-0.0023 ;[f]
0.2 ;[p]
0.0032 ;[f]
0.0059 ;[f]
-0.0014 ;[f]
0.2 ;[p]
0.0029 ;[f]
0.0027 ;[f]
-0.00026 ;[f]
-0.0032 ;[f]
0.2 ;[p]
-0.0025 ;[f]
0.00097 ;[f]
0.0024 ;[f]
0.00197 ;[f]
-0.0080 ;[f]
0.2 ;[p]
0.0031 ;[f]
-0.00571 ;[f]
0.0030 ;[f]
-0.0074 ;[f]
0.0025 ;[f]
0.0034 ;[f]
0.2 ;[p]
0.00973 ;[f]
0.00862 ;[f]
0.0041 ;[f]
0.0046 ;[f]
0.00061 ;[f]
-0.0056 ;[f]
0.0056 ;[f]
0.2 ;[p]
$SIGMA
0.1 ;[p]
0.1 ;[p]
$PRIOR NWPRI
; Omega prior
$OMEGAP BLOCK(8)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2
; degrees of freedom for OMEGA prior
$OMEGAPD
( 8 FIXED ) ; [dfo]
; Starting with a short iterative two stage analysis brings the
; results closer so less time needs to be spent during the
; burn-in of the BAYES analysis
$EST METHOD=ITS INTERACTION SIGL=4 NITER=15 PRINT=1
FILE=example6.ext NOABORT NOPRIOR=1
$EST METHOD=BAYES INTERACTION NBURN=4000 SIGL=4 NITER=10000
PRINT=10 CTYPE=3 FILE=example6.txt NOABORT NOPRIOR=0
; By default, ISAMPLE_M* are 2. Since there are many data points
; per subject, setting these to 1 is enough, and it reduces the
; time of the analysis
ISAMPLE_M1=1 ISAMPLE_M2=1 ISAMPLE_M3=1 IACCEPT=0.4
$COV MATRIX=R UNCONDITIONAL9 Bayes Example 6S
This is example6s.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory. Note the old style syntax for priors.
;Model Desc: Receptor Mediated Clearance model with Dynamic Change
; in Receptors
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example6s (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C
; The new numerical integration solver is used, although ADVAN=9
; is also efficient for this problem.
$SUBROUTINES ADVAN13 TRANS1 TOL=4
$MODEL NCOMPARTMENTS=3
$PRIOR NWPRI NTHETA=8, NETA=8, NTHP=0, NETP=8, NPEXP=1
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
MU_6=THETA(6)
MU_7=THETA(7)
MU_8=THETA(8)
VC=EXP(MU_1+ETA(1))
K10=EXP(MU_2+ETA(2))
K12=EXP(MU_3+ETA(3))
K21=EXP(MU_4+ETA(4))
VM=EXP(MU_5+ETA(5))
KMC=EXP(MU_6+ETA(6))
K03=EXP(MU_7+ETA(7))
K30=EXP(MU_8+ETA(8))
S3=VC
S1=VC
KM=KMC*S1
F3=K03/K30
$DES
DEL=1.0E-10
DADT(1) = -(K10+K12)*A(1) + K21*A(2) - VM*A(1)*A(3)/(A(1)+KM+DEL)
DADT(2) = K12*A(1) - K21*A(2)
DADT(3) = -VM*A(1)*A(3)/(A(1)+KM+DEL) - K30*A(3) + K03
$ERROR
CALLFL=0
ETYPE=1
IF(CMT.NE.1) ETYPE=0
IPRED=F
Y = F + F*ETYPE*EPS(1) + F*(1.0-ETYPE)*EPS(2)
$THETA
;Initial Thetas
(4.0) ; [MU_1]
(-2.1) ; [MU_2]
(0.7) ; [MU_3]
(-0.17) ; [MU_4]
(2.2) ; [MU_5]
(0.14) ; [MU_6]
(3.7) ; [MU_7]
(-0.7) ; [MU_8]
; degrees of freedom for OMEGA prior
( 8 FIXED ) ; [dfo]
;Initial Omegas
$OMEGA BLOCK(8)
0.2 ;[p]
-0.0043 ;[f]
0.2 ;[p]
0.0048 ;[f]
-0.0023 ;[f]
0.2 ;[p]
0.0032 ;[f]
0.0059 ;[f]
-0.0014 ;[f]
0.2 ;[p]
0.0029 ;[f]
0.0027 ;[f]
-0.00026 ;[f]
-0.0032 ;[f]
0.2 ;[p]
-0.0025 ;[f]
0.00097 ;[f]
0.0024 ;[f]
0.00197 ;[f]
-0.0080 ;[f]
0.2 ;[p]
0.0031 ;[f]
-0.00571 ;[f]
0.0030 ;[f]
-0.0074 ;[f]
0.0025 ;[f]
0.0034 ;[f]
0.2 ;[p]
0.00973 ;[f]
0.00862 ;[f]
0.0041 ;[f]
0.0046 ;[f]
0.00061 ;[f]
-0.0056 ;[f]
0.0056 ;[f]
0.2 ;[p]
; Omega prior
$OMEGA BLOCK(8)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2
$SIGMA BLOCK(1) 0.1 ;[p]
;$SIGMA BLOCK(1) 0.1 ;[p]
$SIGMA BLOCK(1) SAME
; Starting with a short iterative two stage analysis brings the
; results closer so less time needs to be spent during the burn-in
; of the BAYES analysis
$EST METHOD=ITS INTERACTION SIGL=4 NITER=30 PRINT=1
FILE=example6s.ext NOABORT NOPRIOR=1 CTYPE=3
$EST METHOD=IMP INTERACTION SIGL=4 NITER=30 PRINT=1
FILE=example6s.ext NOABORT NOPRIOR=1 ISAMPLE=100
$EST METHOD=BAYES INTERACTION NBURN=4000 SIGL=4 NITER=1000
PRINT=10 CTYPE=3 FILE=example6s.txt NOABORT NOPRIOR=0
PSAMPLE_M1=1 PSAMPLE_M2=1 PSAMPLE_M3=1
; By default, ISAMPLE_M* are 2. Since there are many data points
; per subject, setting these to 1 is enough, and it reduces the
; time of the analysis
ISAMPLE_M1=1 ISAMPLE_M2=1 ISAMPLE_M3=1 IACCEPT=0.4
$COV MATRIX=R UNCONDITIONAL10 Bayes Example 6SB
This is example6sb.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory. Note the old style syntax for priors.
;Model Desc: Receptor Mediated Clearance model with Dynamic Change
; in Receptors
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# example6sb (from r2compl)
$INPUT C SET ID JID TIME DV=CONC DOSE=AMT RATE EVID MDV CMT
$DATA example6.csv IGNORE=C
; The new numerical integration solver is used, although ADVAN=9
; is also efficient for this problem.
$SUBROUTINES ADVAN13 TRANS1 TOL=4
$MODEL NCOMPARTMENTS=3
$PRIOR NWPRI NTHETA=8, NETA=8, NTHP=0, NETP=8, NPEXP=1
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
MU_5=THETA(5)
MU_6=THETA(6)
MU_7=THETA(7)
MU_8=THETA(8)
VC=EXP(MU_1+ETA(1))
K10=EXP(MU_2+ETA(2))
K12=EXP(MU_3+ETA(3))
K21=EXP(MU_4+ETA(4))
VM=EXP(MU_5+ETA(5))
KMC=EXP(MU_6+ETA(6))
K03=EXP(MU_7+ETA(7))
K30=EXP(MU_8+ETA(8))
S3=VC
S1=VC
KM=KMC*S1
F3=K03/K30
$DES
DADT(1) = -(K10+K12)*A(1) + K21*A(2) - VM*A(1)*A(3)/(A(1)+KM)
DADT(2) = K12*A(1) - K21*A(2)
DADT(3) = -VM*A(1)*A(3)/(A(1)+KM) - K30*A(3) + K03
$ERROR
CALLFL=0
ETYPE=1
IF(CMT.NE.1) ETYPE=0
IPRED=F
Y = F + F*ETYPE*EPS(1) + F*(1.0-ETYPE)*EPS(1)
$THETA
;Initial Thetas
(4.0) ; [MU_1]
(-2.1) ; [MU_2]
(0.7) ; [MU_3]
(-0.17) ; [MU_4]
(2.2) ; [MU_5]
(0.14) ; [MU_6]
(3.7) ; [MU_7]
(-0.7) ; [MU_8]
; degrees of freedom for OMEGA prior
(8 FIXED) ; [dfo]
;Initial Omegas
$OMEGA BLOCK(8)
0.2 ;[p]
-0.0043 ;[f]
0.2 ;[p]
0.0048 ;[f]
-0.0023 ;[f]
0.2 ;[p]
0.0032 ;[f]
0.0059 ;[f]
-0.0014 ;[f]
0.2 ;[p]
0.0029 ;[f]
0.0027 ;[f]
-0.00026 ;[f]
-0.0032 ;[f]
0.2 ;[p]
-0.0025 ;[f]
0.00097 ;[f]
0.0024 ;[f]
0.00197 ;[f]
-0.0080 ;[f]
0.2 ;[p]
0.0031 ;[f]
-0.00571 ;[f]
0.0030 ;[f]
-0.0074 ;[f]
0.0025 ;[f]
0.0034 ;[f]
0.2 ;[p]
0.00973 ;[f]
0.00862 ;[f]
0.0041 ;[f]
0.0046 ;[f]
0.00061 ;[f]
-0.0056 ;[f]
0.0056 ;[f]
0.2 ;[p]
; Omega prior
$OMEGA BLOCK(8)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.2
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2
$SIGMA BLOCK(1) 0.1 ;[p]
; Starting with a short iterative two stage analysis brings the
; results closer so less time needs to be spent during the burn-in
; of the BAYES analysis
$EST METHOD=ITS INTERACTION SIGL=4 NITER=30 PRINT=1
FILE=example6sb.ext NOABORT NOPRIOR=1 CTYPE=3
$EST METHOD=IMP INTERACTION SIGL=4 NITER=30 PRINT=1
FILE=example6sb.ext NOABORT NOPRIOR=1 ISAMPLE=100
$EST METHOD=BAYES INTERACTION NBURN=4000 SIGL=4 NITER=1000 PRINT=10
CTYPE=3 FILE=example6sb.txt NOABORT NOPRIOR=0
; By default, ISAMPLE_M* are 2. Since there are many data points
; per subject, setting these to 1 is enough, and it reduces the
; time of the analysis
ISAMPLE_M1=1 ISAMPLE_M2=1 ISAMPLE_M3=1 IACCEPT=0.4
$COV MATRIX=R UNCONDITIONAL11 Bayes Example 7R
This is example7r.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Interoccasion Variability
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB run# example7r
$INPUT C SET ID TIME AMT RATE EVID MDV CMT DV OCC
$ABBR REPLACE ETA(OCC_CL)=ETA(3,4,5)
$DATA example7r.csv IGNORE=C
$SUBROUTINES ADVAN1 TRANS2
$PK
MU_1=THETA(1)
MU_2=THETA(2)
V=DEXP(MU_1+ETA(1))
S1=V
VC=V
CL=DEXP(MU_2+ETA(2))*EXP(ETA(OCC_CL))
$ERROR
IPRED=F
Y = F+F*EPS(1)
;Initial Thetas
$THETA
2.0 ;[MU_1]
2.0 ;[MU_2]
;Initial omegas
$OMEGA BLOCK(2)
.3 ;[p]
-.01 ;[f]
.3 ;[p]
$OMEGA BLOCK(1)
.1 ;[p]
$OMEGA BLOCK(1) SAME(2)
$SIGMA
0.1 ;[p]
$PRIOR NWPRI
; Degrees of freedom for Prior Omega blocks
$OMEGAPD (2.0 FIXED) (1.0 FIXED)
; Prior Omegas
$OMEGAP BLOCK(2)
.14 FIX
0.0 .125
$OMEGAP BLOCK(1) .0164 FIX
$OMEGAP BLOCK(1) SAME(2)
$EST METHOD=ITS INTERACTION FILE=example7r.ext NITER=10000
PRINT=5 NOABORT SIGL=8 CTYPE=3 CITER=10
NOPRIOR=1 CALPHA=0.05 NSIG=2
$EST METHOD=SAEM INTERACTION NBURN=30000 NITER=500 SIGL=8
ISAMPLE=2 PRINT=10 SEED=1556678 CTYPE=3
CITER=10 CALPHA=0.05 NOPRIOR=1
$EST METHOD=IMP INTERACTION EONLY=1 NITER=4 ISAMPLE=3000
PRINT=1 SIGL=10 NOPRIOR=1 MAPITER=0
$EST METHOD=BAYES INTERACTION FILE=example7r.txt NBURN=10000
NITER=10000 PRINT=100 CTYPE=3 CITER=10
CALPHA=0.05 NOPRIOR=0
$EST METHOD=COND INTERACTION MAXEVAL=9999 NSIG=3 SIGL=10 PRINT=5
NOABORT NOPRIOR=1
FILE=example7r.ext
$COV MATRIX=R PRINT=E UNCONDITIONAL12 Bayes Example 8
This is example8.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 8 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X SDIX SDSX
$DATA example8.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
include nonmem_reserved_general
; Request extra information for Bayesian analysis.
; An extra call will then be made for accepted samples
BAYES_EXTRA_REQUEST=1
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
; When Bayes_extra=1, then this particular set of individual
; parameters were "accepted" So you may record them if you wish
IF(BAYES_EXTRA==1 .AND. ITER_REPORT>=0 .AND. TIME==0.0) THEN
" WRITE(51,98) ITER_REPORT,ID,CL,V1,Q,V2
" 98 FORMAT(I12,1X,F14.0,4(1X,1PG12.5))
ENDIF
$ERROR
include nonmem_reserved_general
BAYES_EXTRA_REQUEST=1
Y = F + F*EPS(1)
IF(BAYES_EXTRA==1 .AND. ITER_REPORT>=0 ) THEN
" WRITE(52,97) ITER_REPORT,ID,TIME,F
" 97 FORMAT(I12,1X,F14.0,2(1X,1PG12.5))
ENDIF
; Initial values of THETA
$THETA
(2.0) ; [LN(CL)]
(2.0) ; [LN(V1)]
(2.0) ; [LN(Q)]
(2.0) ; [LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA
$SIGMA
(0.6) ; [P]
$PRIOR NWPRI
; Prior information to the Thetas.
$THETAP (2.0 FIX)x4
$THETAPV BLOCK(4) FIX VALUES(10000.0,0.0)
; Prior information to the OMEGAS.
$OMEGAP BLOCK(4)
0.2 FIX
0.0 0.2
0.0 0.0 0.2
0.0 0.0 0.0 0.2
$OMEGAPD (4 FIX)
$EST METHOD=BAYES INTERACTION FILE=example8.ext NBURN=10000
NITER=1000 PRINT=100 NOPRIOR=0 CTYPE=3 CINTERVAL=10013 Bayes Example 9
This is example9.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
;Model Desc: Two compartment Model, Using ADVAN3, TRANS4
;Project Name: nm7examples
;Project ID: NO PROJECT DESCRIPTION
$PROB RUN# Example 9 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X SDIX SDSX
$DATA example9.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4 OTHER=aneal.f90
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
; Initial values of THETA
$THETA
(0.001, 2.0) ; [LN(CL)]
(0.001, 2.0) ; [LN(V1)]
(0.001, 2.0) ; [LN(Q)]
(0.001, 2.0) ; [LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.05 ;[P]
0.01 ;[F]
0.05 ;[P]
0.01 ;[F]
0.01 ;[F]
0.05 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.05 ;[P]
;Initial value of SIGMA
$SIGMA
(0.6) ; [P]
$EST METHOD=SAEM INTERACTION FILE=example9.ext NBURN=5000
NITER=500 PRINT=10 NOABORT SIGL=6 CTYPE=3 CINTERVAL=100
CITER=10 CALPHA=0.0514 Bayes Example 10
This is example10.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
$PROB F_FLAG04est2a.ctl
$INPUT C ID DOSE=AMT TIME DV WT TYPE
$DATA example10.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$PK
CALLFL=1
MU_1=DLOG(THETA(1))
KA=DEXP(MU_1+ETA(1))
MU_2=DLOG(THETA(2))
V=DEXP(MU_2+ETA(2))
MU_3=DLOG(THETA(3))
CL=DEXP(MU_3+ETA(3))
SC=V/1000
$THETA 5.0 10.0 2.0 0.1 0.1
$OMEGA BLOCK (3)
0.5
0.01 0.5
0.01 0.01 0.5
; Because THETA(4) and THETA(5) have no inter-subject variability
; associated with them, the algorithm must use a more computationally
; expensive gradient evaluation for these two parameters
$SIGMA 0.1
$PRIOR NWPRI
; Priors to Omegas
$OMEGAP BLOCK (3)
0.09 FIX
0.0 0.09
0.0 0.0 0.09
$OMEGAPD (3 FIX)
$ERROR
EXPP=THETA(4)+F*THETA(5)
IF (TYPE.EQ.0) THEN
; PK Data
F_FLAG=0
Y=F+F*ERR(1) ; a prediction
ELSE
; Categorical data
F_FLAG=1
; Use protected exponent PEXP, to avoid numerical overflow
A=PEXP(EXPP)
B=1+A
Y=DV*A/B+(1-DV)/B ; a likelihood
ENDIF
$EST METHOD=ITS INTER LAP NITER=1000 PRINT=5 SIGL=6 NSIG=2
NOABORT NOPRIOR=1 CTYPE=3 CITER=10 CALPHA=0.05
FILE=example10.ext
; Because of categorical data, which can make conditional density highly
; non-normal, select a t-distribution with 4 degrees of freedom for
; importance sampling proposal density
$EST METHOD=IMP INTER LAP NITER=1000 PRINT=1 ISAMPLE=300 DF=4 IACCEPT=1.0
$EST METHOD=IMP EONLY=1 NITER=5 ISAMPLE=1000 PRINT=1 DF=4 IACCEPT=1.0 MAPITER=0
$EST METHOD=SAEM EONLY=0 INTER LAP NBURN=2000 NITER=1000 PRINT=50
DF=0 IACCEPT=0.4
$EST METHOD=IMP EONLY=1 NITER=5 ISAMPLE=1000 PRINT=1 DF=4
IACCEPT=1.0 MAPITER=0
$EST METHOD=BAYES NBURN=3000 NSAMPLE=3000 PRINT=100
FILE=example10.txt DF=0 IACCEPT=0.4 NOPRIOR=0
$EST METHOD=COND LAP INTER MAXEVAL=9999 PRINT=1 FILE=example10.ext
NOPRIOR=1 NOHABORT
$COV UNCONDITIONAL PRINT=E MATRIX=R SIGL=10
$TABLE ID DOSE WT TIME TYPE DV A NOPRINT FILE=example10.tab15 Bayes Example 10L
This is example10l.ctl from the NONMEM 7 distribution medium. It, along with the data file, can be found in the examples directory.
$PROB F_FLAG04est2a.ctl
$INPUT C ID DOSE=AMT TIME DV WT TYPE
$DATA example10l.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$PK
CALLFL = 1
MU_1 = THETA(1)
KA = DEXP(MU_1 + ETA(1))
MU_2 = THETA(2)
V = DEXP(MU_2 + ETA(2))
MU_3 = THETA(3)
CL = DEXP(MU_3 + ETA(3))
SC = V / 1000
$THETA 1.6 2.3 0.7 0.1 0.1
$OMEGA BLOCK (3)
0.5
0.01 0.5
0.01 0.01 0.5
; Because THETA(4) and THETA(5) have no inter-subject variability
; associated with them, the algorithm must use a more computationally
; expensive gradient evaluation for these two parameters
$SIGMA 0.1
$PRIOR NWPRI
; Priors to Omegas
$OMEGAP BLOCK (3)
0.09 FIX
0.0 0.09
0.0 0.0 0.09
$OMEGAPD (3 FIX)
$ERROR
EXPP = THETA(4) + F * THETA(5)
IF (TYPE.EQ.0) THEN
; PK Data
F_FLAG = 0
Y = F + F * ERR(1) ; A prediction
ELSE
; Categorical data
F_FLAG = 1
; Use protected exponent PEXP, to avoid numerical overflow
A = PEXP(EXPP)
B = 1 + A
Y = DV * A / B + (1 - DV) / B ; A likelihood
ENDIF
$EST METHOD=ITS INTER LAP NITER=1000 PRINT=5 SIGL=6 NSIG=2 NOABORT
NOPRIOR=1 CTYPE=3 CITER=10 CALPHA=0.05 FILE=example10l.ext
; Because of categorical data, which can make conditional density highly
; non-normal, select a t-distribution with 4 degrees of freedom for
; importance sampling proposal density
$EST METHOD=IMP INTER LAP NITER=1000 PRINT=1 ISAMPLE=300 DF=4 IACCEPT=1.0
$EST METHOD=IMP EONLY=1 NITER=5 ISAMPLE=1000 PRINT=1 DF=4 IACCEPT=1.0 MAPITER=0
$EST METHOD=SAEM EONLY=0 INTER LAP NBURN=2000 NITER=1000 PRINT=50
DF=0 IACCEPT=0.4
$EST METHOD=IMP EONLY=1 NITER=5 ISAMPLE=1000 PRINT=1 DF=4
IACCEPT=1.0 MAPITER=0
$EST METHOD=BAYES NBURN=3000 NSAMPLE=3000 PRINT=100
FILE=example10l.txt DF=0 IACCEPT=0.4 NOPRIOR=0
$EST METHOD=COND LAP INTER MAXEVAL=9999 PRINT=1 FILE=example10l.ext
NOPRIOR=1
$COV UNCONDITIONAL PRINT=E MATRIX=R SIGL=10
$TABLE ID DOSE WT TIME TYPE DV A NOPRINT FILE=example10l.tab16 Bind Example
EXAMPLES OF THE USE OF $BIND
Case 1.
Suppose a fragment of the dose records looks like this:
TIME AMT PREP X
0 100 4 .
10 100 1 90
20 0 0 80Suppose $INPUT and $PK include this code:
$INPUT ID TIME AMT DV PREP X
$PK
F1 = 1
IF (PREP .EQ. 4) F1 = THETA(3)
ALAG1 = THETA(5)
CL = THETA(1) * X * EXP(ETA(1))It is intended that theta(3) be the ratio of bioavailabilities of preparation 4 to preparation 1. However, calls to PK are made only at event times 0, 10 and 20. Suppose the value of theta(5) (lag time) is 5. After 5 time units, the first dose enters the system. The value of PREP (dose preparation type) is taken from the event record at time 10, which is 1. F1 is computed incorrectly. If the estimation step is run, the gradient with theta(3) will be 0, indicating the problem.
Case 2.
If lag time is not modeled (no ALAG1), but the doses are multiple (ADDL>0, II>0), a more subtle error can occur. Suppose a fragment of the data records looks like this:
TIME AMT PREP X II ADDL
0 100 4
.5 1
10 100 1
20 0 0 90 5 1
80 . .Suppose $INPUT and $PK include this code:
$INPUT ID TIME AMT DV PREP X II ADDL
$PK
F1 = 1
IF (PREP .EQ. 4) F1 = THETA(3)
CL = THETA(1) * X * EXP(ETA(1))F1 is computed correctly for the initiating dose records at times 0 and 10, and for the additional dose at time 15, but incorrectly for the additional dose at time 5. Theta(3) applies only to the first dose having PRED=4. (The gradient for theta(3) is not 0.) The same problem could have occurred with the dose having PREP=1, except that by an accident of coding, PREP=0 (on the record having TIME=20) is treated the same as PREP=1, so that the additional dose at time 15 has F1=1 as it should.
On the other hand, the values of X are presumably recorded so that 90 is valid for the time 0 through 10 and 80 is valid for the time 10 through 20. Therefore, the model for CL is correct in both examples.
In either of these two examples, the model can be corrected by the insertion of $BIND and CALLFL=-2:
$INPUT ID TIME AMT DV PREP X II ADDL ; II & ADDL only in case 2
$BIND - - - - DOSE NEXT
$PK
CALLFL=-2Because of CALLFL=-2 in the abbreviated code, additional calls to PK are made at the nonevent dose times (5 in case 1; 5 and 15 in case 2). Because of the $BIND record, PREP has the value at these calls from the initiating dose record (4 at time 5, 1 at time 15). F1 is computed appropriately for the dose preparation. X still has values from the next event record, so the computation of CL is unchanged, e.g., it is the same from time=0 to time=5 as it is from time=5 to time=10.
(See Section 50.5)
17 Bivariate Function
This is a fully-worked out example of the use of the $ABBREVIATED FUNCTION option to declare a user-defined function. In this example, the function BIVARIATE is used to compute the integral of the bivariate normal distribution of two correlated data points, DV1 and DV2, with correlation RHO.
This is bivariate.ctl from the NONMEM 7 distribution medium. This file, the data file bivariate.csv, and the function’s source code bivariate.f90 can be found in the examples directory.
$PROB BIVARIATE EXAMPLE
; THESE DECLARATIONS ALLOW ANY FUNCTION TO HAVE
; ALTERNATIVE DIMENSIONS FOR THEIR ARRAYS
; BUT, USER DEFINED DIMENSIONS ARE PASSED AS THE
; LAST ARGUMENT TO FUNC, SUCH AS:
; BV=BIVARIATE(VBI(1),FNC001_1(1,1),FNC001_2(1,1,1),5)
$ABBR FUNCTION BIVARIATE(VBI,5)
$INPUT SIM ID DOSE DV TIME
$DATA bivariate.csv IGNORE=@
$SUBROUTINES OTHER=bivariate.f90
$PRED
B1=THETA(1)
B2=THETA(2)
B3=THETA(3)
K =LOG(2)/EXP(THETA(4))
ED50=EXP(THETA(5))
U =(1-EXP(-K*TIME ))
MU_1=B1+B3*DOSE/(DOSE+ED50)
MU_2=B2
MXB=MU_1+ETA(1)
MXU=MU_2+ETA(2)
MX =MXB + MXU*U ;***Current model prediction***;
PHIMX=PHI(MX)
IF(NEWIND.NE.2) THEN
TIMEP=0
MXP=0
DVP=0
PHIMXP=0.5
ENDIF
RHOB=(2/(1+EXP(-THETA(8)))-1)
IF(RHOB>0.0) RHO=RHOB**(TIME-TIMEP)
IF(RHOB==0.0) RHO=0.0
IF(RHOB<0.0) RHO=-(-RHOB)**(TIME-TIMEP)
PC =(1-PHIMX) *(1-DV ) + PHIMX*DV
IF(PC.LE.0.0) EXIT
V=SQRT(1+OMEGA(1,1)+OMEGA(2,2)*U**2)
POPP = (B1+B2*U +B3*DOSE/(DOSE+ED50))/V ;*Population mean prediction*;
IF (TIME.EQ.1) THEN
JP=PC
PCP=1.0
ELSE
;***Pass information to bivariate normal***;
VBI(1)=RHO
VBI(2)=MX
VBI(3)=MXP
VBI(4)=1 ;*0 = Upper tail as in Drezner & Wesolowsky; 1 = Bottom tail*;
VBI(5)=1 ;*0 = 3 pt approximation; 1 = 5 point approximation*;
BV=BIVARIATE(VBI)
JP=((DV-1)*(DVP-1)+(DV-1)*(1-2*DVP)*PHIMXP+(DVP-1) &
*(1-2*DV)*PHIMX+(1-2*DV)*(1-2*DVP)*BV)
ENDIF
IF(JP.LE.0.0) EXIT
LOGL=LOG(JP/PCP)
Y = -2*LOGL
MXP=MX
PCP=PC
DVP=DV
TIMEP=TIME
PHIMXP=PHIMX
$THETA
-1.7 ; 1 B1
1.2 ; 2 B2
2.9 ; 3 B3
1.4 ; 4 LOG(B4)
1.2 ; 5 LOG(B5)
(0.0 FIXED) ; 6 LOG SQRT VAR(ETA1)
(0.0 FIXED) ; 7 LOG SQRT VAR(ETA2)
2.2 ; 8 RHO parameter
$OMEGA DIAGONAL(2)
0.8 ; V1
0.8 ; V2
;$EST METHOD=IMP LAPLACE -2LL PRINT=1 NITER=300 ISAMPLE=300
;SIGL=6 CTYPE=3 NOHABORT
$EST MAX=0 PRINT=1 METHOD=1 LAPLACE -2LL SIGL=10 NOHABORT
$COV COMPRESS MATRIX=R PRINT=E UNCONDITIONALHere is a fragment of the data file bivariate.csv:
SIM,ID,DOSE,DV,TIME
1,1,0,0,1
1,1,0,0,2
1,1,0,0,4
1,1,0,0,8
1,1,0,0,16
1,1,0,0,24
1,1,0,0,36This is the file bivariate.f90. (The extension .f90 in the file name is not actually required because the contents of the file is copied to FSUBS and compiled as part of FSUBS. It is not compiled independently.)
! Integral of Bivariate normal distribution of two correlated data points,
! DV1, DV2, correlation RHO.
! BV1, BV2 and GAUSS FROM Drezner and Wesolowsky,
! J. Statistical Computat. Simul. 35, pp. 101-107, 1990.
! BV2 is more accurate for extreme values of rho
! Partial derivatives worked out by Robert Bauer.
!
! How to use in NONMEM 7.4
! VBI is a suggested vector name. Could be any name.
! But the function name must match what is in bivariate.f90
! $ABBR FUNCTION BIVARIATE(VBI,5)
! $SUBROUTINES ... OTHER=bivariate.f90
! $PK
! ...
! DV1 AND DV2 SHOULD BE TWO SEPARATE DATA ITEMS ON EACH RECORD.
! VBI(1)=RHO
! VBI(2)=DV1 ! H
! VBI(3)=DV2 !K
! VBI(4)=DTYPE !=0 FOR H TO INF, K TO INF, OR 1 FOR -INF TO H, -INF TO K
! VBI(5)=BVTYPE ! BVTYPE=0 USES SIMPLE INTEGRATOR, BVTYPE=1 USES MORE
! ACCURATE ONE FOR LARGE VALUES OF RHO
! F_FLAG=1
! IT IS SAFER TO AVOID BV BEING EXACTLY 0, JUST MAKE SURE
! IT IS REALLY SMALL, TO AVOID LOG() ERRORS
! BV=BIVARIATE(VBI)+1.0D-30
! IF YOU USE -2LL FORMAT (DEFAULT):
! Y=-2.0D+00*LOG(BV)
FUNCTION BIVARIATE(X,X1,X2,NDIM)
! NDIM SHOULD EQUAL 5
! IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT NONE
INTEGER NDIM
REAL*8 X(NDIM),X1(NDIM),X2(NDIM,NDIM)
REAL*8 R,H,K,ITYPE,PI,PI2,SPI2,R1,RSQR,HRK,KRH,PHRK,PKRH,XH,XK
REAL*8 XHRK,XKRH,XC,BV,BV2,BIVARIATE,GAUSS,BV1
EXTERNAL GAUSS,BV1,BV2
INTEGER BVTYPE
! unload the arguments
R=X(1)
H=X(2)
K=X(3)
ITYPE=X(4)
BVTYPE=X(5)
! ITYPE=0 INTGRATE H TO INF, K TO INF
! ITYPE=1 INTGRATE -INF TO H -INF TO K
! BVTYPE=0 USES THE SIMPLER BV1
! BVTYPE=1 USES THE MORE ACCURATE BV2
! SET UP SOME USEFUL CALCUATIONS FOR FIRST
! AND SECOND DERIVATIVES OF THE ARGUMENTS
! DERIVATIVES NOT NEEDED IF YOU USE IMP, SAEM, BAYES, OR IF
! YOU USE LAPLACE OR ITS WITH OPTMAP=1 ETADER=3 (NONMEM 7.3), OR
! DERIVATIVES NEEDED IF YOU USE STANDARD LAPLACE
PI=3.141592653589793238D+00
PI2=2.0D+00*pi
spi2=sqrt(pi2)
R1=1.0D+00-R*R
RSQR=SQRT(ABS(R1))
IF(RSQR.LE.0.0D+00) RSQR=1.0D-100
IF(R1.LE.0.0D+00) R1=1.0D-100
HRK=(H-R*K)/RSQR
KRH=(K-R*H)/RSQR
PHRK=1.0D+00-GAUSS(HRK)
PKRH=1.0D+00-GAUSS(KRH)
XH=EXP(-H*H/2.0D+00)/SPI2
XK=EXP(-K*K/2.0D+00)/SPI2
XHRK=EXP(-HRK*HRK/2.0D+00)/SPI2/RSQR
XKRH=EXP(-KRH*KRH/2.0D+00)/SPI2/RSQR
XC=EXP(-(H*H-2.0D+00*R*H*K+K*K)/2.0D+00/R1)/RSQR/PI2
! parital F WRT RHO, from Drezner and Wesolowsky, equation (4)
X1(1)=XC
! parital F WRT H
X1(2)=XH*PKRH+(ITYPE-1.0D+00)*XH
! parital F WRT K
X1(3)=XK*PHRK+(ITYPE-1.0D+00)*XK
! 2ND parital F WRT RHO,RHO
X2(1,1)=R/R1*XC+XC*(H*K*(1.0D+00-R)*(1.0D+00-R)-R*(H-K)*(H-K))/R1/R1
! 2ND parital F WRT RHO,H
X2(1,2)=XC*(R*K-H)/R1
! 2ND parital F WRT RHO,K
X2(1,3)=XC*(R*H-K)/R1
X2(2,1)=X2(1,2)
X2(3,1)=X2(1,3)
! 2ND parital F WRT H,K
X2(2,3)=XC
X2(3,2)=X2(2,3)
! 2ND parital F WRT H,H
X2(2,2)=-X1(2)*H-XC*R
! 2ND parital F WRT K,K
X2(3,3)=-X1(3)*K-XC*R
! write(50,*) 'A ',x1(1),x1(2),x1(3)
! write(50,*) 'B ',x2(1,1),x2(1,2),x2(1,3),x2(2,1), &
! x2(2,2),x2(2,3),x2(3,1),x2(3,2),x2(3,3)
! bV2 is more accurate for exterme rho values.
IF(BVTYPE==0) THEN
BV=BV1(H,K,R)
ELSE
BV=BV2(H,K,R)
ENDIF
IF(ITYPE==1.0D+00) BV=BV-GAUSS(H)-GAUSS(K)+1.0D+00
BIVARIATE=BV
RETURN
END
FUNCTION BV1(H1,H2,R)
IMPLICIT NONE
INTEGER I
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 H1,H2,R,H12,BV,BV1,RR,RR2,H3,GAUSS
EXTERNAL GAUSS
REAL*8 X(5),W(5),DVAL
DATA X/.04691008,.23076534,.5,.76923466,.95308992/
DATA W/.018854042,.038088059,.0452707394,.038088059,.018854042/
H12 =(H1*H1+ H2*H2)/2.0D+00
H3 = H1*H2
BV = 0.0D+00
DO 1 I= 1,5
RR = R*X(I)
RR2= 1.0D+00-RR*RR
IF(RR2.LE.0.0D+00) RR2=1.0D-100
DVAL=(RR*H3 - H12)/RR2
IF(DVAL.GT.300.0D+00) DVAL=300.0D+00
BV = BV + W(I)*EXP(DVAL)/SQRT(RR2)
1 CONTINUE
BV = BV*R + GAUSS(H1)*GAUSS(H2)
BV1=BV
RETURN
END
FUNCTION BV2(H1,HK,R)
IMPLICIT NONE
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 X(5),W(5)
REAL*8 H1,HK,R,H2,H12,BV,BV2,R2,R3,H3,H7,H6,H5,AA,AB,R1,RR,GAUSS,RR2
REAL*8 DVAL,DVAL2
EXTERNAL GAUSS
INTEGER I
DATA X/.04691008,.23076534,.5,.76923466,.95308992/
DATA W/.018854042,.038088059,.0452707394,.038088059,.018854042/
H2 = HK
H12 = (H1*H1+ H2*H2)/2.0D+00
BV = 0.0D+00
IF(ABS(R).LT.0.7D+00)GO TO 4
R2= 1.0D+00-R*R
R3 = SQRT(R2)
IF(R.LT.0.0D+00)H2 = -H2
H3=H1*H2
H7 = EXP( -H3/2.0D+00)
IF(R2.EQ.0.0D+00) GO TO 3
H6 = ABS(H1 - H2)
H5 = H6*H6/2.0D+00
H6 = H6/R3
AA = 0.5D+00 - H3/8.0D+00
AB = 3.0D+00 - 2.0D+00*AA*H5
BV = .13298076D+00*H6*AB*GAUSS(H6) &
-EXP(-H5/R2)*(AB + AA*R2)*.053051647D+00
DO 2 I= 1,5
R1= R3*X(I)
RR = R1*R1
R2 = SQRT(DABS(1.0D+00 - RR))
IF(R2.EQ.0.0D+00) R2=1.0D-100
IF(RR.EQ.0.0D+00) RR=1.0D-100
DVAL=-H5/RR
IF(DVAL.GT.300.0D+00) DVAL=300.0D+00
DVAL2=-H3/(1.0D+00 + R2)
IF(DVAL2.GT.300.0D+00) DVAL2=300.0D+00
BV = BV - W(I)*EXP(DVAL)*(EXP(DVAL2)/R2/H7 - 1.0D+00 - AA*RR)
2 CONTINUE
3 IF(R.GT.0.0D+00)BV = BV*R3*H7 + GAUSS(MAX(H1,H2))
IF(R.LT.0.0D+00)BV = MAX(0.0D+00,GAUSS(H1) - GAUSS(H2)) - BV*R3*H7
BV2=BV
RETURN
4 H3=H1*H2
DO 1 I= 1,5
R1= R*X(I)
RR2= 1.0D+00-R1*R1
IF(RR2.EQ.0.0D+00) RR2=1.0D-100
1 BV = BV + W(I)*EXP((R1*H3 -H12)/RR2)/SQRT(RR2)
BV2 =GAUSS(H1)*GAUSS(H2) + R*BV
RETURN
END
FUNCTION GAUSS(Z)
IMPLICIT NONE
! IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 Z,GAUSS,X,G
REAL*8 A(4)
INTEGER I
DATA A/ -.72657601,.71070688,-.142248368,.127414796/
X= 1.0D+00/(1.0D+00 + .23164189D+00*ABS(Z))
G = .53070271D+00
DO 1 I= 1,4
1 G=G*X+A(I)
GAUSS = G*X*EXP( - Z*Z/2.0D+00)
IF(Z.LT.0.0D+00)GAUSS = 1.0D+00 -GAUSS
RETURN
END
! USE PROGRAM HEADER FOR STAND-ALONE EXECUTABLE TESTING
! PROGRAM BIVTEST
SUBROUTINE BIVTEST
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(9),A1(9),A2(9,9)
REAL*8 AH(9),AH1(9),AH2(9,9)
REAL*8 AJ(9),AJ1(9),AJ2(9,9)
REAL*8 B(9),B1(9),B2(9,9)
REAL*8 H,K
HDEL=1.0D-04
1 CONTINUE
WRITE(*,*) 'ENTER RHO,H,K,INTEGRAND-TYPE,BVTYPE'
READ(*,*) RHO,H,K,DTYPE,BVTYPE
A(1)=RHO
A(2)=H
A(3)=K
A(4)=DTYPE
A(5)=BVTYPE
BV=FUNCA(A,A1,A2)
WRITE(*,*) 'VALUE ',BV
WRITE(*,*) GAUSS(-1.96D+00),GAUSS(1.96D+00)
AH(1:5)=A(1:5)
AH(1)=RHO-HDEL
AJ(1:5)=A(1:5)
AJ(1)=RHO+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(1)=(BVJ-BVH)/2.0D+00/HDEL
B2(1,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(1,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(1,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG RHO',A1(1),A2(1,1),A2(1,2),A2(1,3)
WRITE(*,*) 'NUMER RHO',B1(1),B2(1,1),B2(1,2),B2(1,3)
AH(1:5)=A(1:5)
AH(2)=H-HDEL
AJ(1:5)=A(1:5)
AJ(2)=H+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(2)=(BVJ-BVH)/2.0D+00/HDEL
B2(2,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(2,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(2,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG H',A1(2),A2(2,1),A2(2,2),A2(2,3)
WRITE(*,*) 'NUMER H',B1(2),B2(2,1),B2(2,2),B2(2,3)
AH(1:5)=A(1:5)
AH(3)=K-HDEL
AJ(1:5)=A(1:5)
AJ(3)=K+HDEL
BVH=FUNCA(AH,AH1,AH2)
BVJ=FUNCA(AJ,AJ1,AJ2)
B1(3)=(BVJ-BVH)/2.0D+00/HDEL
B2(3,1)=(AJ1(1)-AH1(1))/2.0D+00/HDEL
B2(3,2)=(AJ1(2)-AH1(2))/2.0D+00/HDEL
B2(3,3)=(AJ1(3)-AH1(3))/2.0D+00/HDEL
WRITE(*,*)
WRITE(*,*) 'ANALOG K',A1(3),A2(3,1),A2(3,2),A2(3,3)
WRITE(*,*) 'NUMER K',B1(3),B2(3,1),B2(3,2),B2(3,3)
GO TO 1
END18 CCONTR Example
This example illustates the use of the CONTR and CCONTR stubs with NONMEM 7. This example was developed by Bill Frame.
$PROBLEM
$DATA NMDATA7.CSV
$INPUT ID DV
$SUB CONTR=CONTR.txt CCONTR=CCONTR.txt
$PRED
W = THETA(2) ;SD
CL = THETA(1)*EXP(ETA(1)) ;CL/F
PRE = 1/CL ;@ SS ASSUMING INPUT RATE = 1
LAM = THETA(3)
Y = (PRE**LAM-1)/LAM+EPS(1)*W
RES1 = (DV-PRE)/W
$THETA
(0,.1) ; CL/F
(0,2) ; SD ADDITIVE
(0,.251) ; BOX COX LAMBDA PARAMETER
$OMEGA
.003
$SIGMA
1 FIX
$EST MAXEVALS=9999 METH=1 PRINT=1
$COV PRINT=EFiles CONTR.txt and CCONTR.txt are as follows:
SUBROUTINE CONTR (ICALL,CNT,IER1,IER2)
DOUBLE PRECISION CNT
CALL NCONTR (CNT,IER1,IER2,L2R)
RETURN
END
SUBROUTINE CCONTR(I,CNT,P1,P2,IER1,IER2)
USE SIZES, ONLY: ISIZE,DPSIZE
USE ROCM_REAL, ONLY: THETA=>THETAC,Y=>DV_ITM2
USE NM_INTERFACE,ONLY: CELS
IMPLICIT NONE
INTEGER(KIND=ISIZE), INTENT(IN OUT) :: I,IER1,IER2
REAL(KIND=DPSIZE), INTENT(IN OUT) :: CNT,P1(:),P2(:,:)
REAL(KIND=DPSIZE) :: ONE,TWO,W
DATA ONE,TWO/1.00D+00,2.00D+00/
SAVE
IF (I.LE.1) RETURN
W=Y(1)
Y(1)=(Y(1)**THETA(3)-ONE)/THETA(3)
CALL CELS (CNT,P1,P2,IER1,IER2)
Y(1)=W
CNT=CNT-TWO*(THETA(3)-ONE)*LOG(Y(1))
RETURN
ENDThe first individual’s data follows.
id,dv
1,223.28
1,41.59
1,87.6
1,72.11
1,216.32
1,118.1
1,232.2919 Control3 Example
This is file CONTROL3 from the NONMEM distribution medium.
$PROBLEM THEOPHYLLINE SINGLE SUBJECT DATA
$INPUT DOSE=AMT TIME CP=DV
$DATA DATA3
$SUBROUTINES ADVAN2
$PK
CALLFL=1
KA=THETA(1)
K=THETA(2)
SC=THETA(3)
$ERROR
Y=F+ERR(1)
$THETA (0,1.7) (0,.102) (0,29)
$ESTIMATION MAXEVAL=240 PRINT=2
$COVR
$TABLE TIME
$SCAT CP VS TIME
$SCAT PRED VS TIME
$SCAT RES VS TIME
$SCAT PRED VS CP UNITThe data is file DATA3 from the NONMEM distribution media.
320 .0 .
1.71
. .27 7.91
8.31
. .52 8.33
6.85
. 1. 6.08
5.4
. 1.92 4.55
3.01
. 3.5
.90
. 5.02
. 7.03
. 9.
. 12.
. 24.320 Control4 Example
This is file CONTROL4 from NONMEM distribution media.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE TIME CP=DV WT
$DATA THEO
$PRED
;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)
;DOSE=WT-ADJUSTED DOSE (MG/KG)
;DS=NON-WT-ADJUSTED DOSE (MG)
IF (DOSE.NE.0) THEN
DS=DOSE*WT
W=WT
ENDIF
KA=THETA(1)+ETA(1)
KE=THETA(2)+ETA(2)
CL=THETA(3)*W+ETA(3)
D=EXP(-KE*TIME)-EXP(-KA*TIME)
E=CL*(KA-KE)
F=DS*KE*KA/E*D
Y=F+EPS(1)
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$SIGMA .4
$EST MAXEVAL=450 PRINT=5
$COV
$TABLE ID DOSE WT TIME
$SCAT (RES WRES) VS TIME BY IDThe data is file THEO from the NONMEM distribution media. The first individual’s data follows.
1 4.02 0. .74 79.6
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .21 Control5 Example
This is file CONTROL5 from NONMEM distribution media.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$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 IDThe data is file THEOPP from the NONMEM distribution media. The first individual’s data follows.
1 4.02 0. . 79.6
1 . 0. .74 .
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .22 Control6 Example
This is file CONTROL6 from NONMEM distribution media. It is discussed in Guide IV, Appendix VIII.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$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 IDThe data is file THEOPP from the NONMEM distribution media. The first individual’s data follows.
1 4.02 0. . 79.6
1 . 0. .74 .
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .23 Control7 Example
This is file CONTROL7 from NONMEM distribution media. It is discussed in Guide IV, Appendix VIII.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$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 IDThe data is file THEOPP from the NONMEM distribution media. The first individual’s data follows.
1 4.02 0. . 79.6
1 . 0. .74 .
1 . 0.25 2.84 .
1 . 0.57 6.57 .
1 . 1.12 10.5 .
1 . 2.02 9.66 .
1 . 3.82 8.58 .
1 . 5.1 8.36 .
1 . 7.03 7.47 .
1 . 9.05 6.89 .
1 . 12.12 5.94 .
1 . 24.37 3.28 .24 Data Average Example
A template record is a data record that contains a non-zero value for the raw-data (RAW_) data item. Only the DV values with those observation records (i.e., records with MDV=0) matching the template data record will be included in the average. This average replaces the DV item in the template record, and it may be displayed in a table or scatterplot.
An observation record matches a template if all its data items are identical to all those on the template, except that the values of NONMEM data items do not affect the match. (Any nonNONMEM data item type may in fact be excluded from the match (See Section 50.30)).
Example 1
Here is a fragment of a data set (which may or may not actually be for PREDPP):
#ID TIME AMT WT AGE EVID MDV DV RAW_ GEND
1 5 0 70 50 0 1 2 1 1 template record
4 5 0 70 50 0 0 4 0 1 observation recordThese two records match. The NONMEM items are the ID, MDV, DV, and RAW_ items. The DV item with the observation record will be included in the average. The DV item in the template record is ignored for the purpose of computing the average (MDV=1); the displayed DV item with this record will be the average value, and the displayed RES item will be the difference between this average and the PRED item. The average is taken over all records matching the TIME, AMT, WT, AGE, EVID, and GEND items of the template record.
If the omit record
$OMIT GENDappears, then the two records
#ID TIME AMT WT AGE EVID MDV DV RAW_ GEND
1 5 0 70 50 0 1 2 1 1 template record
4 5 0 70 50 0 0 4 0 2 observation recordalso match. In this case the average is taken over all records matching the TIME, AMT, WT, AGE, and EVID items of the template record.
Example 2
Here is an example of categorical dose-response with four different bolus doses per individual (separated by washouts), and data averages are formed for each of the four doses.
#ID DOSE DV MDV RAW_ set DV item to avr. for DOSE=5
100 5 . 1 1 set DV item to avr. for DOSE=10
100 10 . 1 1 set DV item to avr. for DOSE=20
100 20 . 1 1 set DV item to avr. for DOSE=40
100 40 . 1 1
1 5 10 0
1 10 1 0 0
1 20 2 0 0
1 40 2 0 0
2 5 10 0
2 10 1 0 0
2 20 1 0 0
2 40 1 0 0
3 5 10 0
3 10 2 0 0
3 20 2 0 0
3 40 3 0 0
.
.Example 3
An average of a PRED-defined item may also be displayed.
$ABBREV COMRES=1 COMSAV=1
...
$PRED
COM(1)=LOG(DV)
...
$TABLE COM(1)=LDVFor template records, the value of LDV displayed in the table is the average of LOG DV over all observation records matching the template.
(See Section 51.6, raw.)
25 Enterhepatic Circulation Examples
The help item “Model Time examples” gives a fragment of code for modelling EHC (Enterohepatic Circulation, which is also called Enterohepatic Recycling). It illustrates the use of MTIME parameters to model instantaneous changes in differential equations. This help item, “Enterhepatic circulation examples”, describes two fully-worked out control streams in the examples directory. They demonstrate how to generalize the fragment to multiple sequential doses.
Both examples use the same data. (hillss.dat and mtimess.dat are identical). There is a steady-state bolus dose at time 0. The interdose interval II is 12, and there are enough additional doses (ADDL=100) to continue the dosing pattern throughout the data set. There are “other” records every 4 units till time 140 to allow compartment amounts to be displayed and there is a final observation record at time 144.
mtimess.ctl
This example incorporates the fragment into a complete control stream. MTIME parameters are used to turn on and off the EHC terms in the differential equations. The variable FLAG is 1 between times MTIME(1) and MTIME(2) after each dose event and turns on the EHC terms. After time MTIME(2) is reached, a new set of MTIME’s is defined which affect the next dosing interval. MTIME parameters are not dose-related parameters and have no effect on steady-state dose events. Even if PK computes MTIME(i)< II, this produces future changes in the system, and does not apply retroactively to the preceding implied doses. (See Guide VI, Section V.F.4, Note 4). A steady-state dose record should not be used. Instead, the SS dose record is described as a transient dose with SS=DROP on the $INPUT record.
hillss.ctl
This control stream does not use MTIME. Instead, a smooth step model using Hill terms in a sigmoid emax model is used. The $DES code has to compute all the necessary variables. Flag1 and flag2 are continuous variables that change from 0 to 1 at the times corresponding to the MTIME’s. The FLAG variable is similar to FLAG in mtimess.ctl. The changes to the differential equations are not instantaneous, but they are continuous. If the exponent in the Hill term is made larger, the predictions approach those of the MTIME model. However, very large values of the exponent can lead to numerical difficulties in PREDPP. Smaller values of the exponent may be more realistic physiologically. A Steady-State dose event record is used with this model.
Note that SS dose records should only be used when the kinetics implemented in the model coincides with the II (interdose interval) of the SS record. Just as MTIME’s cannot affect the differential equations retroactively, changes to the differential equations that happen in the future cannot affect the Steady-state calculations. For example, with II=12, then the kinetics should not be different in the interval of time 0 to 12 vs. time 12 to 24 or time 24 to 36, etc. In each interval, changes occur at theta(8) and theta(8)+theta(9) after the start of the interval. $DES computes the first change time for each interval using the INT function:
mt1 = II * INT(T/II) + theta(8)NM-TRAN gives a warning about the use of the INT function in $DES:
(WARNING 68) THE INT, MOD, MIN, OR MAX FUNCTION IS BEING USED OUTSIDE
OF A
SIMULATION BLOCK. IF THE FUNCTION VALUE AFFECTS THE VALUE OF THE
OBJECTIVE FUNCTION, THEN AN ERROR WILL PROBABLY OCCUR.This warning may be disregarded. Discontinuties occur at the ends of the integration intervals, but the kinetics are unaffected. For example, when T<12, the value of 12INT(T/12) is 0. The values of FLAG1 and FLAG2 are initially 0 and FLAG is 0. As T approaches the end point T=12, the values of FLAG1 and FLAG2 both become 1 and FLAG is 0. At the end point when T=12, the value of 12INT(T/12) is 12 and both FLAG1 and FLAG2 are 0. The discontinuity in mt1 and mt2 does not affect the FLAG variable because FLAG1 and FLAG2 are both 0 or both 1 in the neighborhood of the discontinuity, and the kinetics are continuous.
In these examples, EHC is driven by the dose events. The EHC changes can also driven by the clock. Suppose every 12 hours, a new EHC cycle begins. (The value 12 is chosen so that the two versions will give the same predictions.) E.g., instead in mtimess.ctl, instead of
MTIME(1) = MTIME(1)+IIcompute
MTIME(1) = MTIME(1)+12In hillss.ctl, instead of
mt1=inter*INT(T/inter)+theta(8)compute
mt1=12*INT(T/12)+theta(8)There is no difference in the results.
25.1 Examples Using MTIME to Model Periodic Discontinuities in $DES
Some systems defined with differential equations can be discontinuous with respect to time. Discontinuities are typically introduced in the system by suddenly changing the value of one or more model variables at specific points of time. Such changes can be periodic due to e.g. a circadian rhythm. It is important to define discontinuous variables in $PK using model event time (MTIME) variables rather than updating them in $DES using IF-ELSE-ENDIF tests of T or using the discontinuous INT function. In fact, when integrating from time ta to time tb, routines such as DVERK (ADVAN6) may go slightly beyond time tb (i.e., it may happen that $DES is called at T>tb). The values of DADT should be continuous. If any element of DADT changes at T=b, this should happen at the next integration interval, integrating from time b to time c.
Two examples are provided for the implementation of periodic discontinuities using MTIME variables. The first (step_circexa.ctl) shows how to model the daily reset of a step function and illustrates how NONMEM update the MTIME and closely related MNOW variables at and in between event records. The second (idr_circexa.ctl) applies this step function in an indirect response model that describes a truncated sinusoidal baseline reponse. It also shows how values of T that range from 0 to infinity can be transformed into repeated values of 0-24h. All files may be found in the examples directory. These examples were suggested by Sebastien Bihorel, Luann Phillips and Jill Fiedler-Kelly.
The data set used in both examples is circadian.csv. This data set has 4 subjects with sparse data (one observation every 24 hours from time 0 to time 168) and 4 subjects with finely-spaced data (every 0.5 hours from time 0 to time 48). It is with the second group of subjects that scatters of the step function, variable T24, and indirect response vs. time can best be seen. There are no dose event records (i.e. AMT=0 for all records) because the step function and indirect response are generated endogenously. If this model is incorporated in a larger model, dose records can be used for other compartments.
step_circexa.ctl
This example illustrates how to turn on a step function, FLAG, at a certain time of the day, SHIFT, turn it off after a certain duration, DUR, and repeat this every 24 hours.
^
1 -| --------- --------- ---------
|| | | | | |
|| | | | | |
0 -|--- -------- -------- -------- FLAG
|
---------------------------------------------------> TIME
|| | |
SHIFT-24 0 24 48
<- SHIFT -><- DUR ->The MTIME(1) and MTIME(2) variables respectively define the times at which FLAG is set to 1 and 0, such as MTIME(2)-MTIME(1)=DUR. After MTIME(2) is reached, both variables are incremented by 24, and the FLAG update process perpetuates itself every 24 hours. Compartment 1 is a dummy compartment intended to monitor the proper update of the FLAG variable: the amount A1 starts at 0 and is incremented by DUR every day.
Because the FLAG variable is not necessarily updated at the time of event records, both $TABLE and WRITE statements within the abbreviated code are used in step_circexa.ctl to report the values of the variables of interest. WRITE statements are used because lines in a table file are produced only for event records, but not for nonevent times such as Model times (MTIMES) at which the variables computed in $PK are of special interest here.
Besides MTIME(1) and MTIME(2), the following variables are of interest in step_circexa.ctl:
- MPAST(1) and MPAST(2): variables automatically updated based on MTIME(1) and MTIME(2), and used to set FLAG in
$DES. MPAST(i)=0 until the call to PK subsequent to the one for which MNOW=i. At that call and until MTIME(i) is redefined, MPAST(i)=1. - INTMTIME (“interval MTIME”) is the value of MTIME(1) for the entire integration interval during which FLAG=1, ie from MTIME(1) up to and including the endpoint at MTIME(2). It is not specifically used in
step_circexa.ctl, but will be in the next example. INTMTIME is computed and displayed here so that it can be discussed with the other variables. - TSTATE is a reserved variable giving the time at which the current state vector (compartment amounts) was computed. It gives the time to which the system was most recently advanced.
The FLAG variable can be computed in either $PK or $DES, but the latter is preferred.
Background details:
Suppose that the values of TIME on the event records are t1 t2 t3 etc. and that PK is called with every event record and NONEVENT (ADDITIONAL AND LAGGED) DOSE TIMES AND AT MODEL TIMES.
Suppose there are k MTIME variables with values
t1 <= MTIME(1) <= MTIME(2) <= ... <= MTIME(k) <= t2The interval [t1,t2] is integrated by smaller intervals
[t1, MTIME(1)], [MTIME(1), MTIME(2)] ... [MTIME(k), t2].Calls to DES during any integration interval [ta, tb] will have ta <= T, but it may happen that T > tb. For calls to PK and ERROR, the sequence of calls for the record with TIME=t2 is:
- for each i from 1 to k: call PK with the record at MTIME(i) (MNOW=i, TSTATE=MTIME(i-1)) (if i=1, TSTATE=t1)
- call PK with the record at t2 (MNOW=0; TSTATE=MTIME(k))
- call ERROR with the record at t2 (MNOW is set to 9 to identify the WRITE in
$ERROR; TSTATE=t2)
The same calls to PK occur even if MTIME(k) happens to correspond to t2. In such a case, the call with MNOW=k precedes the call with MNOW=0, and the values in the table file are from the second call.
The lines of the file step_circexa.txt that corresponding to the second event record (TIME=24) for ID 1 (set with SHIFT=0.9088 and DUR=20.8385), are:
MNOW TIME MT1 MT2 MP1 MP2 TSTATE INTMTIME
(1) 1.0000 24.0000 0.9088 21.7473 0.0000 0.0000 0.0000 0.9088
(2) 2.0000 24.0000 24.9088 45.7473 1.0000 0.0000 0.9088 0.9088
(3) 0.0000 24.0000 24.9088 45.7473 0.0000 0.0000 21.7473 24.9088
(4) 9.0000 24.0000 24.9088 45.7473 0.0000 0.0000 24.0000 24.9088where MTi and MPi stand for MTIME(i) and MPAST(i).
At line 1, FLAG=0. This corresponds to the advance from time 0 to MTIME(1): no update of MTIME variables is performed. At line 2, FLAG=1. This corresponds to the advance from MTIME(1) to MTIME(2), at which point the MTIME(1) and MTIME(2) variables are updated, but INTMTIME retains its value. At line 3, FLAG=0. This corresponds to the advance from MTIME(2) to time 24, i.e. the call to PK with MNOW=0, in which INMTIME is updated to the new value of MTIME(1) defined at the previous call to PK. At line 4, FLAG=0. This corresponds to the call to ERROR.
The corresponding line of the table file step_circexa_debug.tab is:
TIME MT1 MT2 STS INTMTIME
24.0000 24.9088 45.7473 21.7473 24.9088Note that TSTATE values are:
- 0 (time of 1st. event record)
- 0.9088 after advance to MTIME(1)=0.9088
- 21.7473 after advance to MTIME(2)=21.7473
- 24.0000 after advance to t2=24.
There will be as many lines with MNOW=0.0 in step_circexa.txt, as there are lines in step_circexa_debug.tab. These lines will be consistent except for the values of MNOW.
INMTIME retains the value of MTIME(1)=0.9088 that pertains to the entire integration interval, up to and including the end point at MTIME(2). It is set in $PK before the MTIME values are changed when MNOW=2. The values in the table are those set in $PK when MNOW=0 (the final call to PK with this data record) by which time INTMTIME and the MTIME(i) have their new values for the next set of advances.
In this particular example, with simulated etas, it sometimes happens that MTIME(1)<t1. From a biological point of view, the step function (and the response) likely started well before the data collection. The beginning of the data collection may occur with FLAG=1 (in the next example, this is during the oscillation of the response) or with FLAG=0 (during the flat part of the response). When MTIME(1) < t1, MTIME(1) is ignored during the set of calls to PK described above.
idr_circexa.ctl
This example extends the concepts introduced in step_circexa.ctl and illustrates how to apply the daily reset of the step function to model a process whose baseline is characterized by a sinusoid function during part of the day and a flat line for the rest of the day. It also show how to create a 24h clock time variable T24 in $DES by transforming the T variable ranging from 0 to infinity into repeating intervals of continuous time between 0 and 24.
^
Rmin+Amp -| --- --- ---
||| || ||
|| | | | | |
|| | | | | |
Rmin -|-- -------- -------- ---- Response R
|
-------------------------------------------------> TIME
|| | |
SHIFT-24 0 24 48
<- SHIFT -><- DUR ->The baseline response depicted above can be described by:
R(t) = | Rmin + Amp*sin(z(t)) for t in [SHIFT+n*24,SHIFT+DUR+n*24] (Eq.1)
| Rmin for any other value of twhere
- Rmin is the minimum response
- Amp is the amplitude of the sinusoid function of time
- z(t) is a circadian function of time which scales the time intervals [SHIFT+n24, SHIFT+DUR+n24] to [0, PI], that is an interval of x such that sin(x)>=0: z(t) = PI * (t-(SHIFT+n*24))/DUR
- n is (-1,0,1,2,…) if (SHIFT+DUR)>24 is ( 0,1,2,…) if (SHIFT+DUR)<=24
The case when SHIFT+DUR>24 is when the flag function is already 1 at time 0, that is it started before the first event. Therefore, an offset of 24 is needed.
z(t) is the function that drives the circadian rhythm of R and is implemented using the INTMTIME variable described in the previous example and that represents SHIFT+n*24. MTIME(1) and MTIME(2) are defined as in step_circexa.ctl to periodically update INTMTIME and the FLAG variable used later in the $DES block.
The baseline function, R(t), can be fitted using Eq.1 in absence of drug effect. However, some biological processes which follow such circadian rhythms can be influenced by drugs (eg, metabolic effects of corticoids). If the effect of the drug is not direct, one can express the biological response and the drug effect in terms of an indirect response model (IDR). Assuming that the drug stimulates the formation of the biological response, one can parameterize the system using the following differential equation:
dR/dt = KIN(t)*(1+S(Cp)) - KOUT*R (Eq.2)where
- KOUT is the elimination rate of the biological response and is stationary
- KIN(t) is the formation rate of the biological response and is time-varying dependent on the step function FLAG: FLAG = MPAST(1)-MPAST(2)
- S(Cp) is a stimulation function of the drug concentration.
In absence of drug, there is no stimulation (i.e., S(Cp) = 0) and the baseline response R(t) described in Eq.1 can be reparameterized as follows:
dR/dt = KIN(t) - KOUT*R (Eq.3)The model is defined such as, when FLAG=0, KIN(t) is constant and equal to KOUT*Rmin, and such as, when FLAG=1, KIN(t) has a sinusoidal shape, increasing from its minimum value to its maximum, and decreasing back to its minimum within the duration DUR. One can obtain the explicit expression of KIN(t) within the time when FLAG=1, by determining the first derivative of R(t), as defined in the first part of Eq.1 and solving for KIN(t) using Eq.3. This expression of KIN(t) is used in $DES.
The $DES block also includes an independent line of code meant to define a continuous 24h clock time. This code is supported by the definition and circadian reset of the MTIME(3) variable in $PK. More specifically, this code is dependent on INTMTIME3, the value of MTIME(3) for the 24-hour intervals.
Note: Original idea for this example stems from Andreas Krause’s post on the NMUsers email list dated from 5/26/2011, subject “Coding INTEGER Function in NONMEM”. The original code was modified to use model event time variables.
(See mtime.)
26 Exogenous Supplementation Example
DISCUSSION:
In this example, an oral “drug” is given exogenously, and it also exists as an endogenous substance. There is an unknown dosing history prior to the observation period (i.e., prior to time zero). This example illustrates how three sources of drug can be modeled: pre-existing endogenous drug, pre-existing drug from an unknown prior dosing history, and drug from known doses. Any combination of the three could be modeled without the others.
The rate of endogenous drug production is assumed to be constant, with no feedback control of production. Thus endogenous drug is at steady-state, and, with linear kinetics, its effect is simply to add a constant increment to exogenous drug in the sampled compartment (the increment is modeled as theta(7)).
For the drug with unknown dosing history, it is assumed that the subject is at steady state with respect to this drug. This part of total drug is modeled by a steady state infusion dose into the depot compartment, ending at time 0, and having an unknown rate. The result of the SS dose is to introduce drug into all compartments of the system (not just the central compartment) because it is distributed throughout the system and is subject to elimination from the system. The unknown rate is modeled as theta(5). NONMEM will adjust theta(5) to best fit not only the “baseline” observation at time 0, but also the later observations.
Note that if samples are not taken sufficiently long after the time of the last dose ( > 4 half-lives), then theta(7) and theta(5) may not be separately identifiable. Note that the value of theta(7) may be determined by the residual concentration after all exogenous drug has disappeared.
A combined additive and ccv error model is used. Theta(8) is the ratio of the C.V. of the proportional component to the standard deviation of the additive component.
Any ADVAN/TRANS combination could be used. Population data could also be modeled in this manner, with eta variables in the $PK block.
$PROBLEM Example of pre-existing drug.
$INPUT ID TIME DV AMT SS II RATE
$DATA DATA1
$SUBROUTINES ADVAN4 TRANS5
$PK
AOB=THETA(1)
ALPHA=THETA(2)
BETA=THETA(3)
KA=THETA(4)
R1=THETA(5)
S2=THETA(6)
$ERROR
FP=THETA(7)+F ; adds endogenous component
W=(1+THETA(8)*THETA(8)*FP*FP)**.5
Y=FP+W*ERR(1)Note that, if there are other doses into the depot compartment with modeled rates, it is necessary to assign a value to R1 conditionally. E.g.,
IF (TIME.EQ.0) THEN
R1=THETA(5) ; rate for SS infusion record at time 0
ELSE
R1=....; rate of other kind of dose
ENDIFNote also that the combined additive and ccv error model can also be modeled using two random variables:
Y = F*(1+ERR(1)) + ERR(2)A fragment of the data follows. Record 1 specifies the SS infusion for the pre-existing drug, which ends at time 0. Record 2 gives the baseline observation. Record 3 specifies an oral bolus dose. Record 4 gives an observation.
1 0 . 0 1 0 -1
1 0 62.2 . . . .
1 0.01 . 95 . . .
1 0.50 235.93 . . . .27 Expectation Example
These are examples of the use of the Marginal (MRG_) data item.
Example 1
Suppose that the probability that a particular subject experiences a pain relief score of 2 is computed. Suppose also one wants to compute the (posterior population) expectation of the probability with each of 4 different bolus doses, not all of which are among those used to obtain observations. A fragment of the control file follows.
$INPUT ID DOSE MDV MRG_ ...
$PRED
....
Y = likelihood of observation given ETA
$EST METH=COND LAPLACE LIKELIHOOD
$TABLE DOSEA fragment of the data follows, with comments following “;”.
#ID DOSE DV MDV MRG_ ; PRED item is set to expectation for DOSE=5
100 5 2 1 1 ; PRED item is set to expectation for DOSE=10
100 10 2 1 1 ; PRED item is set to expectation for DOSE=20
100 20 2 1 1 ; PRED item is set to expectation for DOSE=40
100 40 2 1 1
1 3 10 0
1 10 1 0 0
1 25 2 0 0
1 30 2 0 0
2 3 10 0
2 10 1 0 0
2 25 1 0 0
2 30 1 0 0
3 3 10 0
3 10 2 0 0
3 25 2 0 0
3 30 3 0 0
... etc ...Example 2
This example produces a plot of four residuals, formed by the differences between the raw-data-averages and their (posterior population) expectations, versus the doses used to obtain the data.
$INPUT ID DOSE MDV MRG_ RAW_
$PRED
...
IF (ICALL.EQ.5) THEN
Y = expectation of observation given ETA
ELSE
Y = likelihood of observation given ETA
ENDIF
$EST METHOD=COND LAPLACE LIKELIHOOD
$SCAT RES VS DOSE#ID DOSE DV MDV MRG_ RAW_
100 5 . 1 1 1
100 10 . 1 1 1
100 20 . 1 1 1
100 40 . 1 1 1
1 5 10 0 0
1 10 1 0 0 0
1 20 2 0 0 0
1 40 2 0 0 0
... etc ...Example 3
This example produces a plot of four residuals, formed by the differences between the proportion of subjects in the data set with pain relief score 2 and the (posterior population) expectation of the probability that a subject experiences a pain score of 2, versus the doses used to obtain the data.
$INPUT ID DOSE MDV MRG_ RAW_
$PRED
...
Y = likelihood of observation given ETA
IF (ICALL.EQ.6) THEN
DVR=DV
DV=0
IF (DVR.EQ.2) DV=1
ENDIF
$EST METHOD=COND LAPLACE LIKELIHOOD
$SCAT RES VS DOSE#ID DOSE DV MDV MRG_ RAW_
100 5 2 1 1 1
100 10 2 1 1 1
100 20 2 1 1 1
100 40 2 1 1 1
1 5 10 0 0
1 10 1 0 0 0
1 20 2 0 0 0
1 40 2 0 0 0
2 5 10 0 0
2 10 1 0 0 0
2 20 1 0 0 0
2 40 3 0 0 0
3 5 10 0 0
3 10 2 0 0 0
3 20 2 0 0 0
3 40 3 0 0 0
... etc ...(See mrg, Section 53.40, Section 51.6, raw.)
28 Finalization Example ($PRED and $INFN)
This example contains abbreviated code which can be inserted in a $PRED or $INFN block. This code outputs final parameter estimates, standard errors, minimum value of the objective function, and conditional estimates of etas to various user files. The return codes from Estimation and Covariance steps (zero for normal termination) are also output.
IF (ICALL.EQ.3) THEN
DO WHILE(DATA)
IF (NEWIND.LE.1) WRITE (50,*) ETA
ENDDO
WRITE (51,*) OBJECT
WRITE (52,*) THETA
WRITE (52,*) THETAFR
WRITE (53,*) SETHET
WRITE (53,*) SETHETR
WRITE (54,*) OMEGA(BLOCK)
WRITE (55,*) SEOMEG(BLOCK)
WRITE (56,*) SIGMA(BLOCK)
WRITE (57,*) SESIGM(BLOCK)
WRITE (58,*) IERE,IERC
ENDIF29 Finedata Example
This example is new to NONMEM 7. It is an example of the use of the finedata utility program to fill in missing values of an independent variable (covariate).
Other examples (infn1, infn2) show the use of the INFN routines for interpolation. In data file THEOPP used for CONTROL5, where only one value of WT is present, the INFN routine is used to fill in missing values of WT. Presumably these examples would be part of a NONMEM run involving Estimation.
In this example, finedata is used to fill in the missing values of WT and create a new data set, THEOPPfull. This data set can be used in subsequent NONMEM runs for estimation.
USAGE:
finedata theofine.ctlContents of theofine.ctl are as follows:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$FINEDATA AXIS=TIME FILE=THEOPPfull WT=LIN MISSING=. NEVAL=-1(See control5.) (See infn1 example.) (See infn2 example.)
30 FOCE Parallel Example
The example control stream file foce_parallel.ctl is in the ..\examples directory.
; Used for comparing single versus parallel computing for FOCE method.
;$SIZES LVR=30
;$SIZES LTH=15
;$SIZES LIM1=100
$PROB RUN# Example 1 (from samp5l)
$INPUT C SET ID JID TIME DV=CONC AMT=DOSE RATE EVID MDV CMT CLX V1X QX V2X SDIX SDSX
$DATA example1.csv IGNORE=C
$SUBROUTINES ADVAN3 TRANS4
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_1+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
$ERROR
Y = F + F*EPS(1)
; Initial values of THETA
$THETA
(0.001, 2.0) ; [LN(CL)]
(0.001, 2.0) ; [LN(V1)]
(0.001, 2.0) ; [LN(Q)]
(0.001, 2.0) ; [LN(V2)]
;INITIAL values of OMEGA
$OMEGA BLOCK(4)
0.15 ;[P]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
0.01 ;[F]
0.01 ;[F]
0.01 ;[F]
0.15 ;[P]
;Initial value of SIGMA
$SIGMA
(0.6) ; [P]
;$EST METHOD=CHAIN FILE=foce_parallel.chn ISAMPLE=3 NSAMPLE=0
$EST METHOD=1 INTERACTION NSIG=3 PRINT=1 NOABORT FORMAT=s1PE23.16 MAXEVAL=9999 ; PARAFILE=parallel_file.pnm
$COV MATRIX=S PRINT=E UNCONDITIONAL ; PARAFILE=OFF31 Ignore Accept Example
The options IGNORE=(list) and ACCEPT=(list) of the $DATA record allow for simple selection criteria for records that are to be excluded or included, respectively, in the NONMEM data set. A list has the format:
cond,cond,cond,...where each cond is a test of a data item: “label operator value”. The commas are implied .OR. operators. The .AND. and .NOT. operators are not permitted.
Suppose an .AND. is needed. E.g.,
IGNORE=(GEN.EQ.1.AND.AGE.GT.60)If the expressions are simple, this may be done by negating both sides:
.NOT.IGNORE =.NOT.(GEN.EQ.1.AND.AGE.GT.60)
ACCEPT=.NOT.GEN.EQ.1.OR..NOT.AGE.GT.60Thus,
ACCEPT=(GEN.NE.1,AGE.LE.60)But logical operators inside parenthesis are not permitted, so that the following cannot be coded:
$DATA ACCEPT = ((A==1 .OR. A==2) .AND. B<100)Even if negated, this is:
IGNORE=(.NOT.(A==1.OR.A==2)).OR..NOT.B<100)which is:
IGNORE=((A.NE.1.AND.A.NE.2).OR.B>=100)There is still an .AND. inside parentheses.
There is a workaround that takes advantage of the fact that A takes numeric values and there are no values between 1 and 2. (A.NE.1.AND.A.NE.2) can be implemented with an “or”: (A.LT.1.OR.A.GT.2).
The following are equivalent:
IGNORE=(B.GE.100) IGNORE=(A.LT.1) IGNORE=(A.GT.2)
IGNORE=(B.GE.100,A.LT.1,A.GT.2)This workaround cannot handle continuous values of A or more complicated conditions.
With NONMEM through 7.3, two separate NONMEM problems can be used, in two separate runs, for a complicated condition. In the first problem, abbreviated code is used to append a new variable to the data with a value indicating ignore/accept; in the second run the variable is used in the $DATA statement to choose from the data set.
Suppose ACC is the variable that is to have values 0/1 for ignore/accept. The desired code is:
ACC=0
IF ((A == 1.OR.A == 2).AND.B<100) ACC=1The above can be coded using the FORTRAN language. But NM-TRAN abbreviated code is a subset of FORTRAN. Logical operators .NOT., .AND., .OR. may not be used within parentheses. There are two workarounds.
One is to clear the parentheses:
IF (A==1.AND.B<100.OR.A==2.AND.B<100) ACC=1(this is always possible, no matter how complicted the conditional expression).
The second is to use several statements. There are may ways to do this. In the following, the .AND. of multiple conditions is false if any of the conditions is false.
ACC=1
IF (A.NE.1.AND.A.NE.2) ACC=0
IF (B.GE.100) ACC=0In prob1.ctl, the first problem appends ACC to the (new) data set. In prob2.ctl, the $DATA record uses ACC to select records.
All required data items for NONMEM (and PREDPP, if used) must be present in the original data file testab.dat. The control stream for the first run:
$PROB test of accept Run #1
; Example of implementation of
; $DATA ACCEPT=((A == 1.OR.A == 2).AND.B<100)
$INPUT NO A B ID DV MDV
$DATA testab.dat IGNORE C
$PRED
; Implements:
; ACC=1 IF ((A == 1.OR.A == 2).AND.B < 100)
ACC=1
IF (A.NE.1.AND.A.NE.2) ACC=0
IF (B.GE.100) ACC=0
; The model for this run is unimportant. Keep it simple.
Y=THETA(1)+ETA(1)+EPS(1)
$THETA 1
$OMEGA 1
$SIGMA .4
$TABLE NO A B ID DV MDV ACC NOPRINT NOAPPEND NOHEADER FILE=run1.tabThe control stream for the second run:
$PROB test of accept Run #2
$INPUT NO A B ID DV MDV ACC
; Must do a numeric test for ACC==1, because it appears in
; table file run1.tab as 1.0000E+00
$DATA run1.tab ACCEPT=(ACC.EQN.1)
$PRED
; The model in this problem is the model used for simulation, analysis, etc.
Y=THETA(1)+ETA(1)+EPS(1)
$THETA 1
$OMEGA 1
$SIGMA .4
$TABLE NO A B ID DV MDV
NOPRINT NOAPPEND NOHEADER FILE=run2.tabWith NONMEM 7.4, it is possible to combine both problems into one NONMEM run using the $TABLE record option EXCLUDE_BY list. If any variable in list is non-0, the record is excluded from the table file probe.tab. The $DATA ACCEPT/IGNORE option is not needed. Because probe.tab does not exist when NM-TRAN processes the second problem in the control stream, $DATA .. NOOPEN is needed and a format specification is needed.
$PROB test of NONMEM 7.4 exclude_by feature
$INPUT NO A B ID DV MDV
$DATA testab.dat IGNORE C
$PRED
; sets ACC=0 if (a==1.or.a==2).and.b<100)
; uses the clear the parenthesis approach.
ACC=1
IF (A==1.AND.B<100.OR.A==2.AND.B<100) ACC=0
; The following is the model for both problems.
; It may be more complicated than this.
Y=THETA(1)+ETA(1)+EPS(1)
$THETA 1
$OMEGA 1
$SIGMA .4
; probe.tab has only those records that have ACC=0.
$TABLE NO A B ID DV MDV EXCLUDE_BY ACC NOPRINT NOAPPEND NOHEADER FILE=probe.tab
$PROB test of exclude_by #2
$INPUT NO A B ID DV MDV
$DATA probe.tab (10F12.0) NOOPEN
$THETA 1
$OMEGA 1
$SIGMA .4
; This problem should have $ESTIMATION and other tasks
$TABLE NO A B ID DV NOAPPEND NOHEADER NOPRINT FILE=probe2.tabHere is a fragment of the data testab.dat:
C NO A B ID DV MDV
1 1 100 1 0 0
2 2 100 1 0 0
3 3 100 1 0 0
4 1 101 1 0 0
5 2 101 1 0 0
6 3 101 1 0 0
7 1 99 1 0 0
8 2 99 1 0 0
9 3 99 1 0 0(See Section 50.9, Section 50.46)
32 INFN Interpolation Example 1
This example is adapted from Guide VI (PREDPP), Figure 37. It was revised for NONMEM 7.3. Earlier versions were untested and did not work correctly. This version allows both single-subject and population data. It requires a “missing independent variable” data item.
In this example, the INFN subroutine is used to fill in (interpolate) missing values of WT in the THEOPPMV data set. THEOPPMV is a version of the data set THEOPP used in CONTROL5. (See control5.)
If WT is recorded on every record, then the $PK block can be simpler. The $PK block can be called with every record (no need for CALLFL=1) and no special code is needed to save the value of WT for subsequent records.
Example infn2 shows how the INFN routine may be implemented in abbreviated code. Example fine1 shows how the finedata utility program can be used to fill in the values without the use of the INFN routine.
USAGE:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT MISS
$DATA THEOPPMV
$SUBROUTINES ADVAN2 INFN=infnsubThe file infnsub contains the following:
! THIS EXAMPLE IS ADAPTED FROM GUIDE VI (PREDPP), FIGURE 37.
! FOR HELP FILE infn1.exa
! MODIFIED BY AJB 9/2012
! INFN ROUTINE FOR COMPUTING LINEARLY INTERPOLATED VALUES
! OF AN INDEPENDENT VARIABLE V. ILLUSTRATES USE OF ROUTINE PASS.
! USE FOR SINGLE-SUBJECT OR POPULATION DATA
! WORKS WITH 1 OR MORE NON-MISSING VALUES.
! A SINGLE SUCH VALUE MAY BE ANYWHERE IN DATA SET.
! DATREC(UI)=TIME DATA ITEM
! DATREC(VI)=INDEPENDENT VARIABLE DATA ITEM
! DATREC(MI)=MISSING INDEPENDENT VARIABLE DATA ITEM
! =0 INDEP VAR NOT MISSING
! >0 IF THIS DATA RECORD HAS MISSING INDEP VAR AND:
! =1 THIS DATA RECORD PRECEDES FIRST NON-MISSING VALUE
! =3 THIS DATA RECORD FOLLOWS LAST NON-MISSING VALUE
! =2 THIS DATA RECORD IS BETWEEN NON-MISSING VALUES
!
SUBROUTINE INFN (ICALL,THETA,DATREC,INDXS,NEWIN)
USE SIZES, ONLY: ISIZE,DPSIZE,NO,MAXIDS
! NOTE THAT NEWIND IS NOT OBTAINED FROM THE SUBROUTINE ARGUMENT
USE NMPRD_INT, ONLY: NEWIND=>NWIND
INTEGER(KIND=ISIZE), INTENT(IN) :: ICALL,INDXS(*),NEWIN
REAL(KIND=DPSIZE), INTENT(IN) :: THETA(*)
REAL(KIND=DPSIZE), INTENT(IN OUT) :: DATREC(*)
REAL(KIND=DPSIZE) :: U,V
DIMENSION U(NO,MAXIDS),V(NO,MAXIDS),DEPVAR(MAXIDS),LASTI(MAXIDS)
INTEGER UI,VI,MI,IS,LASTI
! IMPORTANT: CHANGE UI, VI, MI ACCORDING TO LAYOUT OF DATA SET
UI=3
VI=5
MI=6
!
IF (ICALL.NE.1) RETURN
I=0
IS=0
! INITIALIZE PASS
MODE=0
CALL PASS (MODE)
MODE=2
! PASS THROUGH DATA
5 CALL PASS (MODE)
IF (MODE.EQ.0) GO TO 10
IF (NEWIND<2) THEN
I=0
IS=IS+1
ENDIF
! IF INDEP VAR IS PRESENT, STORE TIME AND VALUE
IF (DATREC(MI).EQ.0.) THEN
I=I+1
U(IS,I)=DATREC(UI)
V(IS,I)=DATREC(VI)
LASTI(IS)=I
ENDIF
IF (I == 1) THEN ! SAVE VALUE IN CASE ONLY RECORDED ONCE
DEPVAR(IS)=V(IS,1)
ENDIF
GO TO 5
! INITIALIZE PASS A SECOND TIME
10 I=0
IS=0
MODE=0
CALL PASS (MODE)
MODE=2
! PASS THROUGH DATA A SECOND TIME
15 CALL PASS (MODE)
IF (MODE.EQ.0) RETURN
IF (NEWIND<2) THEN
I=0
IS=IS+1
ENDIF
! IF INDEP VAR IS MISSING AND ONLY RECORDED ONCE, COPY IT
IF (DATREC(MI).NE.0. .AND. LASTI(IS) == 1 ) THEN
DATREC(VI)=DEPVAR(IS)
GO TO 15
ENDIF
! IF INDEP VAR IS MISSING, STORE INTERPOLATED VALUE
IF (DATREC(MI).EQ.0.) THEN
I=I+1
ELSE
IF (DATREC(MI).EQ.1.) THEN ! EXTRAPOL. FROM FIRST 2 VALUES
K=1
L=2
ELSEIF (DATREC(MI).EQ.2.) THEN ! INTERPOL. FROM BEFORE AND AFTER
K=I
L=I+1
ELSEIF (DATREC(MI).EQ.3.) THEN ! EXTRAPOL. FROM LAST 2 VALUES
K=I-1
L=I
ENDIF
B=(V(IS,K)-V(IS,L))/(U(IS,K)-U(IS,L))
DATREC(VI)=V(IS,K)+B*(DATREC(UI)-U(IS,K))
ENDIF
GO TO 15
ENDData file THEOPPMV is identical to THEOPP, but has an extra data item MISS (Missing Dependent Variable). The data for the first subject is:
1 4.02 0. . 79.6 0
1 . 0. .74 .3
1 . 0.25 2.84 .3
1 . 0.57 6.57 .3
1 . 1.12 10.5 .3
1 . 2.02 9.66 .3
1 . 3.82 8.58 .3
1 . 5.1 8.36 .3
1 . 7.03 7.47 .3
1 . 9.05 6.89 .3
1 . 12.12 5.94 .3
1 . 24.37 3.28 .3(See infn2 example.) (See fine1 example.)
33 INFN Interpolation Example 2
This example is new to NONMEM 7. It is an improved version of infn1.exa. In this example, instead of using a stand-alone INFN subroutine, $INFN abbreviated code is used to fill in missing values of WT in the THEOPP data set. (See control5.)
If WT is recorded on every record, then the $PK block can be simpler. The $PK block can be called with every record (no need for CALLFL=1) and no special code is needed to save the value of WT for subsequent records.
No “missing independent variable” data item is needed, so the THEOPP data file can be used as-is. The $ABBR abbreviated code is located in file infnabbr for convenient re-use. This is also an example of the $ABBR REPLACE feature, which allows easy replacement of character strings in the included code, and is also an example of declared subscripted variables.
Example fine1 shows how the finedata utility program can be used to fill in the values without the use of the INFN routine.
USAGE:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN2
$ABBR REPLACE INTVBL=WT
$ABBR REPLACE NULLVAL=0.0
INCLUDE 'infnabbr'The file infnabbr contains the following:
; INFNABBR: ABBREVIATED CODE FOR HELP infn2.exa
; THIS INFN BLOCK INTERPOLATES AND EXTRAPOLATES VALUES OF A SINGLE
; DATA ITEM CALLED INTVBL.
; BY DEFAULT, THE AXIS VARIABLE IS CALLED TIME.
; THE CHARACTER STRING NULLVAL DENOTES A NULL VALUE.
;
; IF THERE IS ONLY ONE NON-NULL VALUE, THIS VALUE IS COPIED TO
; ALL NULL VALUES. OTHERWISE:
; LINEAR INTERPOLATION IS PERFORMED FOR NULL VALUES THAT
; LIE BETWEEN TWO NON-NULL VALUES.
; LINEAR EXTRAPOLATION IS PERFORMED FOR VALUES BEFORE THE
; FIRST NON-NULL VALUE
; LINEAR EXTRAPOLATION IS PERFORMED FOR VALUES AFTER THE
; LAST NON-NULL VALUE
; INFNABBR CAN BE USED WITH BOTH SINGLE
; SUBJECT AND POPULATION DATA SET.
; WRITTEN BY AJB 9/2012
;
; TYPICAL USAGE:
; $ABBR REPLACE INTVBL=WT
; $ABBR REPLACE NULLVAL=0.0
; $ABBR REPLACE TIME=T (ONLY IF T IS THE AXIS VARIABLE RATHER THAN TIME)
; INCLUDE 'INFNINTERP'
;
;
$ABBR DECLARE U(MAXIDS,NO), V(MAXIDS,NO)
$ABBR DECLARE INTEGER IS
$ABBR DECLARE INTEGER I, INTEGER J, INTEGER L
$ABBR DECLARE IVAL(NO)
$INFN
; $INFN FOR COMPUTING INTERPOLATED VALUES OF AN
; INDEPENDENT VARIABLE NAMED INTVBL
; E.g., INTVBL MIGHT BE WEIGHT
; ASSUME THAT INTVBL = NULLVAL WHEN MISSING.
; IS=SUBJECT #
; IVAL(IS)= # OF NON-NULL VALUES OF THE INDEP VAR. FOR SUBJECT IS
; VI(IS,I) = ITH. NON-NULL VALUE OF THE INDEPENDENT VARIABLE
; UI(IS,I) = TIME OF THE ITH. NON-NULL VALUE OF THE INDEPENDENT VARIABLE
IS=0
IF (ICALL.EQ.1) THEN
; FIRST PASS. SAVE VALUES OF TIME & INDEP VAR. WHEN INDEP VAR. IS NON-NULL
DO WHILE (DATA)
IF (NEWIND < 2) THEN ; INITIALIZE NEW INDIVIDUAL
IS = IS + 1
IVAL(IS) = 0
ENDIF
IF (INTVBL /= NULLVAL) THEN ; SAVE NON-NULL VALUE
IVAL(IS) = IVAL(IS) + 1
I = IVAL(IS)
U(IS,I) = TIME
V(IS,I) = INTVBL
ENDIF
ENDDO
; SECOND PASS
IS=0
DO WHILE (DATA)
IF (NEWIND < 2) THEN ; INITIALIZE NEW INDIVIDUAL
IS=IS+1
I=0
ENDIF
; IF INDEP VAR IS MISSING AND WAS ONLY NON-NULL ONCE, COPY IT
IF (INTVBL == NULLVAL .AND. IVAL(IS) == 1) THEN
INTVBL=V(IS,1)
ELSE
IF (INTVBL /= NULLVAL) THEN ; COUNT ANOTHER NON-NULL VALUE
I = I + 1
ELSE ; CURRENT RECORD HAS NULL VALUE
IF (I == 0) THEN ; EXTRAPOLATE FROM FIRST 2 VALUES
J = 1
L = 2
ELSEIF (I > 0 .AND. I < IVAL(IS)) THEN ; INTERPOL. BEFORE AND AFTER
J = I
L = I+1
ELSEIF (I == IVAL(IS)) THEN ; EXTRAPOL. FROM LAST 2 VALUES
J = I-1
L =I
ENDIF
SLOPE = (V(IS,J)-V(IS,L))/(U(IS,J)-U(IS,L))
INTVBL = V(IS,J)+SLOPE*(TIME-U(IS,J))
ENDIF
ENDIF
ENDDO
ENDIF(See control5.) (See infn1 example.) (See fine1 example.)
34 Initial Steady State Example
This is an example of the use of the Initial Steady State (I_SS) reserved variable in $PK. It was provided by Nick Holford, Dept Pharmacology & Clinical Pharmacology, University of Auckland.
Glucose and insulin are both endogenous substances, whose rates of formation are given by terms RGLU and RINS in the differential equations. In the absence of dose events, both glucose and insulin are assumed to be at steady state levels. Exogenous doses of both glucose and insulin are also present in the first subject’s data.
It is possible to intialize the steady state condition with a dose record having AMT=0, RATE=0, SS=1, but this is artificial. The initial condition that the system is at steady state is a feature of the model, not of the data. In this example, the $PK block contains the statement I_SS=1. This causes the system to be initialized appropriately. The statement I_SS=1 could also have been placed among the $MODEL statements. However the initialization is performed, the result is the same: non-zero concentrations GLU and INS are computed at TIME=0.
The SS data item need not be present in the data set for the I_SS computation. It is present in this example in order to demonstrate that SS may also be computed using a dose record (See the third and fourth subjects’ data).
Note that SIGMA and OMEGA are fixed to 0 for purposes of illustration, so that the simulated values are exactly equal to those that arise from the model without the random variability due to ETA and EPS. NONMEM issues warning messages (“INITIAL ESTIMATE OF WITHIN INDIVIDUAL VARIANCE IS ZER0”). These may be ignored.
This example also illustrates the use of A(n) on the right in the $ERROR block to obtain the compartment amounts that are used for the computation of concentrations GLU and INS.
$PROB nicksiss6.ctl I_SS=1 in $PK ADVAN6
$DATA Nicksiss.dat
$INPUT ID TIME CMT AMT RATE SS DV
$SIM (200070927) ONLYSIM NSUB=1
;Silber HE, Jauslin PM, Frey N, Gieschke R, Simonsson US, Karlsson MO.
;An integrated model for glucose and insulin regulation in healthy
;volunteers and type 2 diabetic patients following intravenous glucose
;provocations. J Clin Pharmacol. 2007 Sep;47(9):1159-71.
$THETA
40 ; POP_RGLU MMOL/H/70KG - glucose input rate
40 ; POP_VGLU L/70KG - volume of distribution for glucose
5 ; POP_CLGLU L/H/70KG - glucose clearance
1 ; POP_EMXGLU - EMAX for effect of glucose on insulin input
10 ; POP_C50GLU MMOL/L - EC50 for effect of glucose on Insulin input
5 ; POP_HILGLU - Hill parameter for effect of glucose on Insulin input
5000 ; POP_RINS PMOL/H/70KG - insulin input rate
5 ; POP_VINS L/70KG - volume of distribution for insulin
70 ; POP_CLINS L/H/70KG - insulin clearance
2 ; POP_EMXINS - EMAX for effect of insulin on glucose clearance
50 ; POP_C50INS PMOL/L - EC50 for effect of insulin on glucose clearance
2 ; POP_HILINS - Hill parameter for effect of insulin on glucose clearance
$OMEGA
0 FIX ;PPV_RGLU
0 FIX ;PPV_RINS
$SIGMA ;G_EXP_RUV
0 FIXED ;0.1 ;G_ADD_RUV MMOL/L
0 FIXED ;1 ;I_EXP_RUV
0 FIXED ;0.1 ;I_ADD_RUV PMOL/L
0 FIXED ;1
$SUBR ADVAN6 TOL=3
$MODEL
COMP (GLUCOSE)
COMP (INSULIN)
$PK
I_SS=1
; GLUCOSE
RGLU=THETA(1)*EXP(ETA(1))
VGLU=THETA(2)
CLGLU=THETA(3)
EMXGLU=THETA(4)
C50GLU=THETA(5)
HILGLU=THETA(6)
;INSULIN
RINS=THETA(7)*EXP(ETA(2))
VINS=THETA(8)
CLINS=THETA(9)
EMXINS=THETA(10)
C50INS=THETA(11)
HILINS=THETA(12)
S1=VGLU
S2=VINS
$DES
DGLU=A(1)/VGLU
DINS=A(2)/VINS
DGLUH=DGLU**HILGLU
DGEFF=EMXGLU*DGLUH/(C50GLU**HILGLU+DGLUH) ; effect of glucose on insulin
DINSH=DINS**HILINS
DIEFF=EMXINS*DINSH/(C50INS**HILINS+DINSH) ; effect of insulin on glucose
DADT(1)=RGLU - CLGLU*(1+DIEFF)*DGLU ; glucose
DADT(2)=RINS*(1+DGEFF)- CLINS*DINS ; insulin
$ERROR
GLU=A(1)/VGLU ; glucose concentration
INS=A(2)/VINS ; insulin concentration
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 CMT AMT RATE SS GLU INS
ONEHEADER NOPRINT FILE=nicksiss6.tabThe data (Nicksiss.dat) follows.
#ID,TIME,CMT,AMT,RATE,SS,DV
1 0 1 0 0 0.
1 0 2 0 0 0.
1 .1 1 25 0 0 .
1 1.0 2 10 0 0 .
1 3 1 0 0 0.
1 3 2 0 0 0.
1 100 1 . . . .
1 100 2 . . . .
2 0 1 0 0 0.
2 0 2 0 0 0.
2 100 1 . . . .
2 100 2 . . . .
3 0 1 0 0 1.
3 100 1 . . . .
3 100 2 . . . .
4 0 1 0 0 1.
4 100 1 . . . .
4 100 2 . . . .(See Initial Steady State, I_SS, Section 53.45) (See i_ss, initial_condition.)
35 Interoccasion Variability Example
This example is based on work described in Karlsson and Sheiner, JPB 21(6):735-750 (1994). The model accounts for random variability in a subject’s parameters between study occasions.
It illustrates a data file in which EVID data items are present. In the fragment of data below, EVID has values 1 (dose event), 0 (observation event), and 3 (reset event).
The reset event is present at the start of the subject’s second study occasion. It causes PREDPP to re-initialize the kinetic system as if for a new individual, and permits time to be reset to 0. Note that the subject has the same dosage regimen prior to both study occasions (a bolus dose of AMT 1000 every 12 hours leading to steady state). However, it would have been possible to specify a different dosing regimen prior to the second occasion.
In the abbreviated code, EXIT statements are used to constrain TVCL and TVV (typical value of clearance and volume) to positive values. Because the models for TVCL and TVV include both thetas and covariates, it is difficult to bound the thetas so as to insure that TVCL and TVV are non-negative for all subjects. The NOABORT option is present on the $ESTIMATION record. Hence NONMEM will perform PRED error recovery, and will avoid values of theta (subsequent to the initial values) that result in non-positive TVCL or TVV for any subject.
Explanation of user data items:
- AGE: age
- RACE: race (1 or 2)
- HT: height
- HCTZ: hydrochlorothiazide (0 or 1)
- OCC: occasion (0 or 1)
;THIS IS A NMTRAN CONTROL STREAM For THE PRAZOSIN PK ANALYSIS
$PROB prazosin data
$INPUT ID AGE RACE HT HCTZ AMT TIME DV EVID SS II OCC
$DATA praz21
$SUBROUTINE ADVAN2 TRANS2
$PK
OCC2=1-OCC
R=0
IF (RACE.EQ.2) R=1
TMP=(HT-160)*THETA(1)
IF (TMP.LE.0.) TMP=0.0
TMP2=(AGE-60)*THETA(2)
IF (AGE.LE.60.) TMP2=0.0
TVCL=THETA(3)+TMP-TMP2
TVCL=TVCL+THETA(4)*R+THETA(5)*HCTZ
IF (TVCL.LE.0.) EXIT 1 100
TVV=THETA(6)+HCTZ*THETA(7)
IF (TVV.LE.0.) EXIT 1 200
TVKA=THETA(8)
CL=TVCL*EXP(ETA(3)*OCC+ETA(5)*OCC2+ETA(1))
V =TVV *EXP(ETA(4)*OCC+ETA(6)*OCC2+ETA(2))
KA=TVKA*EXP(ETA(8)*OCC+ETA(9)*OCC2+ETA(7))
S2=V
$THETA (0,0.4) (0.,0.3) 20. 4 -4.
$THETA 80 -25
$THETA (0,1.)
$OMEGA BLOCK(2) 0.2 0.1 0.2
$OMEGA BLOCK(2) 0.2 0.1 0.2
$OMEGA BLOCK(2) SAME
$OMEGA 1.2
$OMEGA BLOCK(1) .5
$OMEGA BLOCK(1) SAME
$ERROR
Y=F*(1+EPS(1))
$SIGMA 0.1
$EST NOABORT SIG=3 MAX=3000 PRINT=10 POSTHOC
$COVARIANCE
$TABLE ID OCC ETA1 ETA3 ETA5 ETA2 ETA4 ETA6 NOPRINT FILE=praz32.tabData for one subject follow.
21 57 1 180 1 1000 0 0.62 1 1 12 0
21 57 1 180 1 0.0 1 12.06 0 0 0 0
21 57 1 180 1 0.0 3 6.81 0 0 0 0
21 57 1 180 1 0.0 4 4.89 0 0 0 0
21 57 1 180 1 0.0 5 4.04 0 0 0 0
21 57 1 180 1 0.0 6 2.82 0 0 0 0
21 57 1 180 1 0.0 7 2.72 0 0 0 0
21 57 1 180 1 0.0 8 1.78 0 0 0 0
21 57 1 180 1 0.0 10 1.07 0 0 0 0
21 57 1 180 1 0.0 12 0.75 0 0 0 0
21 57 1 180 1 0.0 0 2.83 3 0 0 0
21 57 1 180 1 1000 0 2.83 1 1 12 1
21 57 1 180 1 0.0 1 6.27 0 0 0 1
21 57 1 180 1 0.0 2 12.72 0 0 0 1
21 57 1 180 1 0.0 3 9.99 0 0 0 1
21 57 1 180 1 0.0 4 6.90 0 0 0 1
21 57 1 180 1 0.0 5 6.59 0 0 0 1
21 57 1 180 1 0.0 6 4.07 0 0 0 1
21 57 1 180 1 0.0 7 3.11 0 0 0 1
21 57 1 180 1 0.0 8 2.79 0 0 0 1
21 57 1 180 1 0.0 10 1.72 0 0 0 1
21 57 1 180 1 0.0 12 1.17 0 0 0 1(See Section 50.15, Section 53.1, evid.)
36 Logistic Regression Example
Here is a typical NM-TRAN control stream for dichotomous response population data. Such data is referred to as “odd-type” data. Note that the eta is a population eta but there are no epsilons.
$PROB DICHOTOMOUS RESPONSE
$DATA data
$INPUT ID DOSE DV
$PRED
LOGIT=THETA(1)+THETA(2)*DOSE**THETA(3)+ETA(1)
A=EXP(LOGIT)
P=A/(1+A)
IF (DV.EQ.1) THEN
Y=P
ELSE
Y=1-P
ENDIF
$THETA .1 (25,100) (1,3,5)
$OMEGA 1
$ESTIMATION MET=COND LAPLACE LIKE MAX=500
$SCAT ETA(1) VS DOSEHere is a typical NM-TRAN control stream when, one wants to first simulate the data.
$PROB DICHOTOMOUS RESPONSE
$DATA data
$INPUT ID DOSE DV
$PRED
LOGIT=THETA(1)+THETA(2)*DOSE**THETA(3)+ETA(1)
A=EXP(LOGIT)
P=A/(1+A)
IF (ICALL.EQ.4) THEN
CALL RANDOM (2,R)
IF (R.LE.P) DV=1
IF (R.GT.P) DV=0
ENDIF
IF (DV.EQ.1) Y=P
IF (DV.EQ.0) Y=1-P
$THETA .1 (25,100) (1,3,5)
$OMEGA 1
$SIMULATION (7755399) (45211 UNIFORM)
$ESTIMATION MET=COND LAPLACE LIKE MAX=500
$SCAT ETA(1) VS DOSE37 Michaelis Menten Example
The kinetic system in this example is composed of depot, peripheral, and central compartments, and includes Michaelis-Menten Elimination from the central compartment. Either ADVAN6 or ADVAN8 could be used.
As with any ADVAN, one or more doses could be input to any of the compartments. Steady state doses are also possible.
$PROBLEM EXAMPLE OF M-M ELIMINATION WITH DEPOT AND PERIPH
$INPUT ID TIME AMT DV
$DATA datafile
$SUBROUTINES ADVAN6 TRANS1 TOL=4
$MODEL COMP=(DEPOT,DEFDOS),COMP=(CENTRAL,DEFOBS) COMP=(PERIPH)
$PK
VM = THETA(1)*EXP(ETA(1))
KM = THETA(2)*EXP(ETA(2))
S2 = THETA(3)*EXP(ETA(3))
K12 = THETA(4)*EXP(ETA(4))
K23 = THETA(5)*EXP(ETA(5))
K32 = THETA(6)*EXP(ETA(6))
$ERROR
Y = F + ERR(1)
$DES
C2 = A(2)/S2
DADT(1) = -K12*A(1)
DADT(2) = K12*A(1) -K23*A(2) +K32*A(3) -C2*VM/(KM+C2)
DADT(3) = K23*A(2) -K32*A(3)
; requires $THETA, $OMEGA, $SIGMA etc.If the units of KM is amount, rather than concentration, the $DES block need not compute C2. The change is as follows:
DADT(2) = K12*A(1) -K23*A(2) +K32*A(3) -A(2)*VM/(KM+A(2))(ADVAN10 assumes the units of KM is amount.)
38 Mixture Model Example
In this example of a mixture model, it is assumed that some unknown fraction of the population has one set of typical values of CL and Vd, and that the remaining fraction has another set of typical values. Reserved variables of interest:
MIXNUM-
NONMEM sets
MIXNUMto 1 or 2 according to whether PK is to compute parameter values for the 1st or 2nd subpopulations, respectively. MIXEST-
For each individual, NONMEM computes an estimate of the number of the subpopulation of which each individual is a member, and stores this estimate in the integer variable
MIXEST. EST-
In the
$PKblock,MIXESTis stored in the PK-defined variableEST, which enablesMIXESTto be displayed in tables and scatterplots.
$PROB PHENOBARB POPULATION DATA MIXTURE MODEL
$DATA DATA2 (6F7.0)
$INPUT ID TIME AMT WT APGR CP=DV
$SUBROUTINES ADVAN1 TRANS2
$PK
CALLFL=1
EST=MIXEST
IF (MIXNUM.EQ.2) THEN
CL=THETA(2)*THETA(1)*EXP(ETA(3))
V=THETA(4)*THETA(3)*EXP(ETA(4))
ELSE
CL=THETA(1)*EXP(ETA(1))
V=THETA(3)*EXP(ETA(2))
ENDIF
S1=V
$ERROR
Y=F*EXP(EPS(1))
$MIX
P(1)=THETA(5)
P(2)=1.-THETA(5)
NSPOP=2
$THETA (0,.0047) (0,1) (0,.99) (0,1) (0,.5,1)
$OMEGA BLOCK(2) .05 .01 .03
$OMEGA BLOCK(2) SAME
$SIGMA .02
$ESTM MAXEVAL=500 PRINT=5
$TABLE ID EST FIRSTONLY NOAPPEND(See displayed_values, Section 57.5, mixnum mixest).
REFERENCES: Guide VI, section III.L.2 , Figures 6, 23
39 Model Time Examples
These are examples of the use of model times in PREDPP.
Enterohepatic Recycling
This fragment of abbreviated code may be used to model EHC. The transfer of drug from compt. 4 to 1 is controlled by FLAG, which is 1 between the times specified by MTIME(1) and MTIME(2), and is 0 otherwise.
$PK
MTIME(1)=THETA(8)
MTIME(2)=MTIME(1)+THETA(9)
....
$DES
FLAG=MPAST(1)-MPAST(2)
DADT(1)=-KA*A(1)+K41*A(4)*FLAG
DADT(4)=K1G*A(2)-K41*A(4)*FLAG
....Flexible Input Rate
This fragment of abbreviated code may be used to model a flexible modeled infusion rate R1.
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.
In the three assignment statments for R1, exactly one of the right-side expressions is non-zero at each call to PK.
$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)Changing a Model Time
This fragment of abbreviated code shows how a model time parameter can 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.
IF (TIME.EQ.0) TEMP=0
TEMP=TEMP+1
MTIME(1)=TEMP
MTDIFF=1Changing Absorption Rate
This fragment of abbreviated code shows how the Absorption Rate KA can be changed. KA is XKA1 from time 0 to THETA(5). KA is XKA2 after time THETA(5).
$PK
MTIME(1)=THETA(5) ; change point for KA
XKA1=THETA(3)*EXP(ETA(3))
XKA2=THETA(4)*EXP(ETA(4))
KA=XKA1*(1-MPAST(1))+XKA2*(MPAST(1))(See mtime).
See also “Cirdadian example: Examples Using MTIME to Model Periodic Discontinuities in $DES”.
REFERENCES: None.
40 Moving Area Under The Curve Example
The following set of examples was developed in response to a question to the NONMEM Users about modeling a moving value of the area under the curve. The original question (from Pavel Belo) asked: “Efficacy is frequently considered a function of AUC. A disadvantage of this model of efficacy is that the effect is irreversable because AUC of concentration can only increase; it cannot decrease. In many cases, a more meaningful model is one where AUC is calculated from time t -a to t (kind of”moving average”), where t is time in the system of differential equations (variable T in NONMEM).”
DUPLICATE SYSTEM METHOD
Bob Bauer proposed a data set DELAYDATA and control file for computing a moving value of AUC. Bob wrote “In the following simple absorption model example developed by me and Alison Boeckmann for illustration purposes, compartments 1, 2, and 3 are the”real time” depot, central, auc, and compartments 4,5,6 are the “delayed time” depot, central, auc. So, the base model (non-time delay) system (compartments 1,2,3) is replicated (compartments 4,5,6) for the time delay portion. In addition, the data set duplicates the dose information of compartment 1 into compartment 4, and setting ALAG4 to a non-zero value in the control stream file provides a lag time to any doses inputted into compartment 4 (so this would take care of multiple dose problems as well). This allows for assessment and availability of AUC(t) and AUCT(t-time-delay) at any time t.”
DELAYDATA and aucdelay1.ctl
The data for one individual (DELAYDATA) and the control file (aucdelay1.ctl) from the email are provided in the examples directory of NONMEM 7.4 and higher. Note that aucdelay1.ctl and all the examples discussed below may be used with data from more than one individual and different numbers of event records per individual.
In aucdelay1.ctl, the model for time delay is
TDY=THETA(1)*EXP(ETA(1))and the variable DAUC (“difference in area under the curve”) gives the difference in the area under the curve
DAUC=A(3)-A(6) ; AUC(T)-AUC(T-TDY)One significant modification to the control file aucdelay1.ctl was made. The $ERROR model was changed from Y=F+EPS(1) to
Y=F+DAUC+EPS(1)This demonstrates that DAUC can modify Y (although the relationship would not usually be this simple).
SINGLE SYSTEM METHOD
In aucdelay1.ctl the value of DAUC is available continuously at every time value T during the integration, but in fact it is not used in the differential equations themselves. DAUC is used only in the $ERROR block. In this case, an alternate method is possible. Only the “real time” single system of differential equations is integrated, and the needed values of AUC are captured “on the fly” in the $PK block. This should result in faster run times because fewer differential equations need to be numerically integreated.
The model event time variable (MTIME) is used instead of absorption lag (ALAG) and is modelled with the same time delay variable TDY. For each event record, MTIME is used to interrupt the integration at the time of interest TIME-TDY. Variables called SAVEA are used to capture the value of A(3) at model event times for use in the $ERROR block.
The general algorithm is rather complicated, so a simple algorithm will be described first.
SINGLE SYSTEM METHOD: SIMPLE MODEL
A simple model is possible when TDY is always smaller than the time between any two event records.
aucdelay1S.ctl
This is identical to aucdelay1.ctl except that THETA(1) has initial estimate .5 and is constrained to be between 0 and 1, and TDY is constrained to be less than 1, which is the change in time between any two event records.
aucdelay2S.ctl
This implements the simple model. Only the original (single) set of compartments 1-3 are defined. aucdelay2S.ctl uses the same data file DELAYDATA, but ignores records with CMT=4 (doses into the delayed dose compartment). For each event record, it is necessary to have available NEXTT, the TIME of the next event record. This data item is not recorded in DELAYDATA. Instead, NEXTT is listed in the $INPUT record as an extra data item beyond the data items (columns) listed in the data file, and its entries and are created (“transgenerated”) during each run by code in the $INFN block. The control stream has the same constraints on THETA(1) and TDY as aucdelay1S.ctl.
Only one model event time variable, MTIME(1) is needed. Model event time MTIME(1) is modelled as NEXTT-TDY and changes with every event record. SAVEA is a variable that saves the value of the AUC compartment A(3) when $PK is evaluated at the model event time. At such times, MTIME(1) is reset for the next record. The important lines of code are relatively simple, and use TSTATE (“State Time”; the time at which the state-vector A was last computed.)
$PK
TDY=THETA(1)*EXP(ETA(1))
MTIME(1)=NEXTT-TDY
IF (TSTATE==TIME-TDY) SAVEA=A(3)
...
$ERROR
DAUC=A(3)-SAVEA
...The two control files generate the same predictions.
SINGLE SYSTEM METHOD: GENERAL MODEL
The model is more complicated if TDY may be greater than the time between event records, because each value of compartment A(3) must be saved until needed. For the ith. event record, with TIME(i), the value of SAVEAi at TIME(i)-TDY is needed, and is obtained by using a value of MTIME(i) equal to TIME(i)-TDY.
aucdelay2.ctl and aucdelay3.ctl
The general model is implemented in aucdelay2.ctl and aucdelay3.ctl, These versions integrate a single system of differential equations. File aucdelay3.ctl is a short-hand version of aucdelay2.ctl. After processing by the nmtemplate and doexpand utilites, it gives a control stream that is essentially identical to aucdelay2.ctl. It may be easier to look at the more compact code in aucdelay3.ctl.
The maximum number of records in any individual record is set using nmtemplate. DELAYDATA contains 16 records. Commands to run aucdelay3.ctl are:
nmtemplate aucdelay3.ctl temp1 maxrecs=16
doexpand temp1 temp2
nmfe75 temp2 aucdelay3.res(The most up-to-date doexpand and nmtemplate utilities are provided in https://nonmem.iconplc.com/utilities/)
Like aucdelay2S.ctl, aucdelay3.ctl uses DELAYDATA, and ignores records with CMT=4. An $INFN block is used to save all values of TIME in an array TIMES(i). There is no need for an explicit NEXTT variable because all values of TIME are available in this array. Note that the TIMES array contains all TIME values in the data set, which typically contains more than one subject’s data. In the abbreviated code, reserved variable IRECIDX gives the position in TIMES of the start of the current individual record. (LIREC gives the number of records in the current individual record and is needed if different subjects have different numbers of data records so that maxrecs is too large for some subjects.)
At the start of an individual record, a value of MTIME(i) is computed for every data record as TIMES(KREC)-TDY. (KREC uses reserved variable IRECIDX, the first record in the data set, and hence in TIMES, for the current individual record.) Some values of MTIME(i) may be <=0; such values are ignored by PREDPP.
When PK is called at a model event time MTIME(i) the value of A(3) is stored as SAVEAi. There may be several calls to PK for a given value of TSTATE and MTIME(i) but SAVEAi must be assigned on only the first of these calls because only at this call will A(3) have derivaties with respect to ETA (if MTIME depends on ETA). This is why MTIME(i)=0 is necessary after the value of A(3) (and its ETA derivatives, if any) has been saved.
REMARKS
- RECURSIVE CODE
Random variables are saved for future use by using with recursive code, such as in aucdelay2.ctl:
IF (NEWIND.LE.1.AND.1<=LIREC) THEN
MTIME(1)=TIMES(KREC)-TDY ; SET NEW VALUE OF MTIME (TDY MAY DEPEND ON ETA)
SAVEA1=0 ; ERASE OLD VALUE FROM PREVIOUS PASS
ELSE
MTIME(1)=MTIME(1)
SAVEA1=SAVEA1
ENDIFIt is instructive to note that there is one warning that may be ignored:
(WARNING 3) THERE MAY BE AN ERROR IN THE ABBREVIATED CODE. THE FOL-
LOWING
ONE OR MORE RANDOM VARIABLES ARE DEFINED WITH "IF" STATEMENTS THAT
DO NOT
PROVIDE DEFINITIONS FOR BOTH THE "THEN" AND "ELSE" CASES. IF ALL
CONDITIONS FAIL, THE VALUES OF THESE VARIABLES WILL BE ZERO.
SAThe code in $ERROR is such that SA is always given a value and SA never defaults to zero. Recursive code for SA could be used to avoid this warning if desired:
IF ([I]<=LIREC.AND.TIME==TIMES(KREC)) THEN
SA=SAVEA[I]
ELSE
SA=SA
ENDIF- SAVEA variables
The variables SAVEAi save the value of compartment 3 (AUC) at TIME(i)-TDY. They must be assigned individually because there is no true array feature in NMTRAN for random variables.
- SIMULATION
The control streams can be used to simulate values of DV. There will be a difference between simulated values from aucdelay1.ctl vs. aucdelay2.ctl (2 and 3 are the same). The extra record in aucdelay1.ctl’s FDATA causes an extra call to SIMEPS (even though it is a dose record), so the simulated etas and eps are not identical. To have identical simulated values, then do not use
$DATA DELAYDATA IGNORE=@ IGNORE=(CMT==4)Instead, use a version of DELAYDATA in which the second record has EVID=2 (“other-type event record”), so that it is physically present but does nothing. This adds another record to the data set. Instead of
nmtemplate aucdelay3.ctl temp1 maxrecs=16the value maxrecs=17 should be used.
All files may be found in the examples directory:
DELAYDATA
aucdelay1S.ctl
aucdelay2S.ctl
aucdelay1.ctl
aucdelay2.ctl
aucdelay3.ctlREFERENCES: Guide Introduction_7
41 Multiple Dose Example
This example is based on work described in Fattinger K. et al., Netilmicin in the neonate - population pharmacokinetic analysis and dosing recommendations, Clinical Pharmacology and Therapeutics, 1991, Jul, V50; N1:55-65. It illustrates a data file with multiple doses of different kinds.
A FORTRAN format statement is needed to describe the data file because each data record spans two FORTRAN records and also because blanks are sometimes used for null data items.
Explanation of data items:
TIME- hours
KG- bodyweight(kg)
LENG- length(cm)
SEX- 0 for female and 1 for male
GAGE- gestational age at birth
AGE- actual age (in days)
AGE1- age at start of therapy
IMIV- im or iv dose (1.0 is iv and 0.0 is im)
CREA- creatinine value
TAD- time after last dose (h)
All doses are given as zero order infusions. The RATE data item is -2, indicating that duration D1 is modeled in the $PK block. There are two types of doses: IV (intravenous) and IM (intramuscular). The duration of IV doses is 2 minutes. The duration of IM doses is unknown and is given by theta(7).
The Estimation Step produces a Model Specification File, which may be used in in a later run; (See Section 50.25, Section 52.2).
$PROBLEM NETILMICIN
$INPUT ID MDV TIME AMT RATE CP=DV KG LENG SEX GAGE AGE AGE1 IMIV CREA TAD
$DATA data (6F7.0/7X,9F7.0)
$SUBROUTINE ADVAN3,TRANS4
$PK
TVCL1=THETA(1)*KG
GAGEN=GAGE/39.0
PAN=(AGE+1.0)/7.0
IF(GAGE.LT.39.0) THEN
TVCL=TVCL1*(GAGEN**THETA(3))*(PAN**THETA(2))
ELSE
TVCL=TVCL1*(PAN**THETA(2))
ENDIF
CL=TVCL*EXP(ETA(1))
TVV1=THETA(4)*KG
V1=TVV1*EXP(ETA(2))
Q=THETA(5)
V2=THETA(6)*KG
S1=V1
IF(IMIV.GE.1.0) THEN
D1=.0333333 ; 2./60.
ELSE
D1=THETA(7)
ENDIF
$ERROR
Y=F*(1+EPS(1))
$THETA (0,0.07) (0,0.1) (0,1.1) (0,0.4) (0,0.07) (0,0.8) (0,.033)
$OMEGA 0.02 0.02
$SIGMA 0.02
$ESTIMATION PRINT=5 MSFO=msfo1
$COVARIANCE
$TABLE MDV ID KG LENG SEX GAGE AGE TAD BY MDV ID NOPRINT FILE=NTABLE.DAT
$SCATTER PRED VS CP UNITData for one subject follow. This subject received both IM and IV doses.
6.0 1.0 0.0 13.0 -2.0
6.0 2.590 47.0 0.0 36.0 0.0 0.0 1.0 56.0 0.0
6.0 1.0 24.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 1.0 0.0 1.0 56.0 0.0
6.0 0.0 26.0 0.0 0.0 6.95
6.0 2.590 47.0 0.0 36.0 1.0 0.0 1.0 56.0 2.0 0.0
6.0 0.0 32.0 0.0 0.0 3.45
6.0 2.590 47.0 0.0 36.0 1.0 0.0 1.0 56.0 8.0 0.0
6.0 0.0 40.0 0.0 0.0 1.8
6.0 2.590 47.0 0.0 36.0 2.0 0.0 1.0 56.0 16.0 0.0
6.0 1.0 48.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 2.0 0.0 0.0 56.0 0.0
6.0 1.0 72.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 3.0 0.0 0.0 56.0 0.0
6.0 1.0 96.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 4.0 0.0 0.0 56.0 0.0
6.0 1.0 120.0 10.4 -2.0
6.0 2.590 47.0 0.0 36.0 5.0 0.0 0.0 56.0 0.0
6.0 0.0 122.0 0.0 0.0 5.9
6.0 2.590 47.0 0.0 36.0 5.0 0.0 0.0 56.0 0.0
6.0 0.0 128.0 0.0 0.0 2.1
6.0 2.590 47.0 0.0 36.0 5.0 0.0 0.0 56.0 8.0 0.0
6.0 0.0 136.0 0.0 0.0 0.75
6.0 2.590 47.0 0.0 36.0 6.0 0.0 0.0 56.0 16.0 0.0(See duration, Section 50.9).
REFERENCES: Guide VI, section V.D REFERENCES: Guide V, section 6.8 , 6.10
42 Nonparametric Example
DISCUSSION:
Three examples are given of the NONPARAMETRIC step in NONMEM. In examples 1 and 2, it is assumed that there only two etas in the problem. With more than two etas, insert the appropriate additional code.
Example 1
This example illustrates how the cumulative density can be plotted versus the conditional etas.
The following code is inserted in the $PRED or $PK block.
DN1=CDEN_(1)
DN2=CDEN_(2)The following control records are used. Note that with the $NONPARAMETRIC record, the MARGINALS option is the default and need not be coded explicitly.
$ESTIMATION POSTHOC PRINT=1
$NONPARAMETRIC MARGINALS
$SCAT DN1 VS ETA1
$SCAT DN2 VS ETA2Example 2
This example illustrates how both the conditional etas and the conditional nonparametric etas (CNPE) can be obtained, and how they can be compared. (See Section 53.13).
The following code is inserted in the $PRED or $PK block:
;CONDITIONAL PARAMETRIC ETAS
IF (COMACT.EQ.2) ET1A=ETA(1)
IF (COMACT.EQ.2) ET2A=ETA(2)
;CONDITIONAL NONPARAMETRIC ETAS (CNPE)
IF (COMACT.EQ.3) ET1=ETA(1)
IF (COMACT.EQ.3) ET2=ETA(2)The following control statements may be used. The first two are required; the $SCATTER and $TABLE statements are examples of what can be used.
$ESTIMATION POSTHOC
$NONP ETAS
; Examine distribution of nonparametric etas
$SCAT ET1 VS ID
$SCAT ET2 VS ID
; Compare parametric and nonparametric etas
$SCAT ET1A VS ET1
$SCAT ET2A VS ET2
$TABLE ID ET1A ET2A ET1 ET2 NOPRINTExample 3
This example shows the specific statements needed to simulate from a nonparametric distribution using a Model Specification file generated from a previous analysis using $NONPARAMETRIC.
; Model specification file from a previous analysis using $NONPARAMETRIC
$MSFI nonpfit.msf
; Second seed will be used for simulation using the nonparametric density
; Must include TRUE=FINAL, otherwise initial instead of final parameter
; estimates will be used.
$SIMULATION (3432192) (6348371 NONPARAMETRIC) TRUE=FINALREFERENCES: None.
43 NWPRI Examples
These examples illustrate the use of the NONMEM utility routine NWPRI by the NM-TRAN $PRIOR record, and also by an equivalent user-supplied routine PRIOR. For a discussion of PRIOR: (See Section 57.6). For a discussion of the NWPRI routine and in particular, a discussion of the $THETA, $OMEGA, $SIGMA records: (See Section 59.12).
Example 1
Use of PRIOR and NWPRI and a control stream for single-subject data:
This example obtains parameter estimates from single subject data, using a two-compartment PK model, and it incorporates a prior of multivariate normal form for all of the THETA vector.
$PROB SINGLE SUBJECT DATA WITH PRIOR ON THETA
$INPUT ID DOSE TIME DV WT
$DATA data1
$PRIOR NWPRI NTHETA=3,NETA=1,NTHP=3,NPEXP=1
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=CLEARANCE (LITERS/HR)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.EQ.0) AMT=DOSE*WT
T1=THETA(1)
T2=THETA(2)
T3=THETA(3)
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+ETA(1)
$THETA (.4,1.7,7) (.025,.102,.4) (.3,3,30)
;mode for prior with THETA:
$THETA 2.77 FIX .0781 FIX 2.63 FIX
$OMEGA .388
;var-cov for prior with THETA:
$OMEGA BLOCK (3) 5.55 .00524 .00024 -.128 .00911 .515 FIX
$ESTInstead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=priorThe prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
NTHETA=3
NETA=1
NTHP=3
NPEXP=1
PLEV=0.
ITYP=0
NSAM=0
ISS=0
NEPS=0
NETP=0
NEPP=0
CALL NWPRI(NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,ITYP,PLEV, &
NSAM,ISS,CNT)
RETURN
ENDExample 2
First example of use of PRIOR and NWPRI and a control stream for population data:
This example obtains parameter estimates from population data, using a two-compartment PK model, and it incorporates a prior of multivariate normal/inverse Wishart form for all of the THETA vector and OMEGA matrix.
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data2
$PRIOR NWPRI NTHETA=3 NETA=3 NTHP=3 NETP=3 NPEXP=1
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(1)
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
;mode for prior on THETA:
$THETA 3 FIX .08 FIX .04 FIX
;df for prior on OMEGA:
$THETA 12 FIX
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
;var-cov for prior on THETA:
$OMEGA BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX
;mode for prior on OMEGA:
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08 FIX
$SIGMA .4
$ESTInstead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=priorThe prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
NTHETA=3
NETA=3
NTHP=3
NETP=3
NPEXP=1
PLEV=0.
ITYP=0
NSAM=0
ISS=0
NEPS=0
NEPP=0
CALL NWPRI(NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,ITYP,PLEV, &
NSAM,ISS,CNT)
RETURN
ENDExample 3
Second example of use of PRIOR and NWPRI and a control stream for population data:
This example simulates the THETA and OMEGA parameters of the same two-compartment population PK model, as well as the population data from this model. During problem finalization, the entire post-simulation theta vector and omega array in common NMPR16 are written to external files.
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data2
$PRIOR NWPRI NTHETA=3 NETA=3 NTHP=3 NETP=3 NPEXP=1 ITYP=0 PLEV=.9999
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(1)
IF (ICALL.EQ.3) THEN
WRITE (97,*) THSIMP,THSIMPR
WRITE (98,*) OMSIMP(BLOCK)
ENDIF
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
;mode for prior on THETA:
$THETA 3 FIX .08 FIX .04 FIX
;df for prior on OMEGA:
$THETA 12 FIX
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
;var-cov for prior on THETA:
$OMEGA BLOCK (3) .494 .00207 .0000847 .000692 .0000471 .0000292 FIX
;mode for prior on OMEGA:
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08 FIX
$SIGMA .4
$SIM (547676) ONLYSuppose a $TABLE record is added to the above control stream. To obtain values of DV computed with the simlulated values of eta, include the following code:
IF (ICALL.EQ.4) THEN
T0=THSIMP(1)*EXP(ETA(1))
T2=THSIMP(2)*EXP(ETA(2))
T3=THSIMP(3)*W*EXP(ETA(3))
ELSE
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T3=THETA(3)*W*EXP(ETA(3))
ENDIFInstead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=priorThe prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
PLEV=.9999
NTHETA=3
NETA=3
NTHP=3
NETP=3
NPEXP=1
NEPS=0
NEPP=0
ITYP=0
NSAM=0
ISS=0
CALL NWPRI(NTHETA,NETA,NEPS,NTHP,NETP,NEPP,NPEXP,ITYP,PLEV, &
NSAM,ISS,CNT)
RETURN
ENDREFERENCES: none.
44 Objective Function Value Individual Example
Note: With NONMEM 7, the additional output file root.phi contains the same information. (See additional_output_file).
This example illustrates how individual contributions to the objective function may be obtained. (See Section 53.69).
This example is a modification of CONTROL4 from the NONMEM distribution media. It is for NONMEM VI 2.0 and later.
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE TIME CP=DV WT
$DATA THEO
$PRED
;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)
;DOSE=WT-ADJUSTED DOSE (MG/KG)
;DS=NON-WT-ADJUSTED DOSE (MG)
IF (DOSE.NE.0) THEN
DS=DOSE*WT
W=WT
ENDIF
KA=THETA(1)+ETA(1)
KE=THETA(2)+ETA(2)
CL=THETA(3)*W+ETA(3)
D=EXP(-KE*TIME)-EXP(-KA*TIME)
E=CL*(KA-KE)
F=DS*KE*KA/E*D
Y=F+EPS(1)
IF (ICALL.EQ.3) THEN
WRITE (70,*) IIDX,CNTID
ENDIF
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
$SIGMA .4
$EST MAXEVAL=450 PRINT=5Pairs of values are written, one pair per line, one pair for each individual record. Each pair contains the ID data item and the contribution the objective function.
1.000000000000000E+00 1.328986294492194E+01
2.000000000000000E+00 2.071754196741930E+01
...They may also be displayed in a table, using
$ABBR COMRES=2and code such as the following in the $ERROR or $PK block:
IF (COMACT.EQ.1) THEN
COM(1)=IIDX(NIREC)
COM(2)=CNTID(NIREC)
ENDIFThe following, for example, will produced a separate table for the values:
$TABLE IID=COM(1) CNT=COM(2) FILE=comvals NOAPPEND NOPRINT FIRSTONLYNote: With earlier versions than NONMEM 7.3, verbatim code is needed:
" COM(1)=IIDX(NIREC)
" COM(2)=CNTID(NIREC)REFERENCES: None.
45 Phenobarb Example
This example illustrates simple kinetics.
$PROBLEM PHENOBARB SIMPLE MODEL
$INPUT ID TIME AMT WGT APGR DV
$DATA PHENO
$SUBR ADVAN1 TRANS2
$PK
TVCL=THETA(1)
TVV=THETA(2)
CL=TVCL*EXP(ETA(1))
V=TVV*EXP(ETA(2))
S1=V
$THETA (0,.0105) (0,1.05)
$OMEGA .25 .25
$ERROR
W=F
Y=F+W*EPS(1)
IPRED=F ; individual-specific prediction
IRES=DV-IPRED ; individual-specific residual
IWRES=IRES/W ; individual-specific weighted residual
$SIGMA .04
$ESTIMATION SIGDIGITS=4 POSTHOC
$COVARIANCE
$TABLE ID WGT APGR TIME AMT CL V ETA1 ETA2 IPRED IWRES NOPRINT ONEHEADER FILE=TABLE1
$SCATTER DV VS PRED UNIT
$SCATTER (RES,WRES) VS PRED
$SCATTER (CL,V) VS (WGT,APGR)
$SCATTER IWRES VS IPRED BY MDV ORDZEROData for the first individual follow.
1 0. 25.0 1.4 7 .
1 2.0 . 1.4 7 17.3
1 12.5 3.5 1.4 7 .
1 24.5 3.5 1.4 7 .
1 37.0 3.5 1.4 7 .
1 48.0 3.5 1.4 7 .
1 60.5 3.5 1.4 7 .
1 72.5 3.5 1.4 7 .
1 85.3 3.5 1.4 7 .
1 96.5 3.5 1.4 7 .
1 108.5 3.5 1.4 7 .
1 112.5 . 1.4 7 31.0REFERENCES: Guide III, section VI.2 REFERENCES: Guide V, section 2.3
46 PK PD Sequential 1 Example
This is an example for analysis of direct PD using previously estimated individuals’ PK parameters. (See PK_PD_sequential_2 example).
$PROBLEM EMAX MODEL APPLIED TO EFFECT SITE CONCENTRATION
; NOTE: INDIVIDUAL K, V ESTIMATED FROM PREVIOUS PK FIT
; AND RECORDED IN DATA
;
; THE DATA FILE CONTAINS *ONLY* EFFECT OBSERVATIONS,
; CMT = 2 NEEDED TO TELL PREDPP TO SET F = CE
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT KK V1
; |---|
; KK AND V1 ARE ESTIMATES OF K AND V
$SUBROUTINES ADVAN3
$PK
K = KK
K12= .001*K ; TRIVIAL LOSS TO EFFECT COMPT.
K21= THETA(1)*EXP(ETA(1))
EMAX=THETA(2)*EXP(ETA(2))
C50= THETA(3)*EXP(ETA(3))
S2= V1*K12/K21 ; SO THAT CESS = CPSS
$ERROR
Y =EMAX*F/(C50+F)*EXP(ERR(1))
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.If, moreover, an absorption lag (ALAG1) parameter is modeled in $PK with parameters estimated from the PD data only, then this introduces a lag in the dose - Ce relationship. Since this lag was not present in the dose - Cp relationship, it may be interpreted as a lag in the Cp - Ce relationship.
REFERENCES: Guide IV, section V.C.6 REFERENCES: Guide VI, section IV.B.2 , VII.C.3
47 PK PD Sequential 2 Example
This is an example for analysis of direct PD using previously estimated population PK parameters. (See PK_PD_sequential_1 example).
$PROBLEM EMAX MODEL APPLIED TO EFFECT SITE CONCENTRATION
; THE DATA FILE CONTAINS *ONLY* EFFECT OBSERVATIONS,
; POP PK PARAMS ESTIMATED IN PREVIOUS RUN & FIXED HERE.
; CMT = 2 NEEDED TO TELL PREDPP TO SET F = CE
;
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINES ADVAN3
$PK
K = THETA(1)*EXP(ETA(1))
K12= .001*K ; TRIVIAL LOSS TO EFFECT COMPT.
K21= THETA(2)*EXP(ETA(2))
V1= THETA(3)*EXP(ETA(3))
EMAX=THETA(4)*EXP(ETA(4))
C50= THETA(5)*EXP(ETA(5))
S2= V1*K12/K21 ; SO THAT CESS = CPSS
$ERROR
Y= EMAX*F/(C50+F)*EXP(ERR(1))
$THETA
.124 FIXED ; K - FIXED TO PREV EST
.1 ; KEO - ESTIMATED
28.6 FIXED ; V1 - FIXED TO PREV EST
100 ; EMAX- ESTIMATED
50 ; C50 - ESTIMATED
$OMEGA
.113 FIXED ; OMEGA(K) - FIXED TO PREV EST
.04
.627 FIXED ; OMEGA(V1)- FIXED TO PREV EST
.04
.04
$SIGMA .16
; ADD $ESTIMATION, ETC, AS DESIRED.If, moreover, an absorption lag (ALAG1) parameter is modeled in $PK with parameters estimated from the PD data only, then this introduces a lag in the dose - Ce relationship. Since this lag was not present in the dose - Cp relationship, it may be interpreted as a lag in the Cp - Ce relationship.
REFERENCES: Guide IV, section V.C.6 REFERENCES: Guide VI, section IV.B.2 , VII.C.3
48 PK PD Simultaneous 1 Example
This is an example for simultaneous analysis of PK and direct PD. Kinetics are one compartment. (See PK_PD_simultaneous_2, PK_PD_simultaneous_3).
$PROBLEM SIMULT PK/PD - ADVAN3
; THE DATA FILE CONTAINS BOTH CP AND EFFECT OBSERVATIONS.
; WHEN DV IS A CP OBSERVATION, CMT = 1 (OR 0)
; WHEN DV IS AN EFFECT OBSERVATION, CMT =2
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINES ADVAN3
$PK
K = THETA(1)*EXP(ETA(1))
K12= .001*K ; TRIVIAL LOSS TO EFFECT CMPT
K21= THETA(2)*EXP(ETA(2)) ; KEO
S1= THETA(3)*EXP(ETA(3)) ; V1 FOR DRUG
EMAX=THETA(4)*EXP(ETA(4))
C50= THETA(5)*EXP(ETA(5))
S2= S1*K12/K21 ; SO THAT CESS = CPSS
$ERROR
IF (CMT.EQ.2) THEN
Y = EMAX*F/(C50+F)*EXP(ERR(2)) ; DV = EFFECT IF CMT=2
ELSE
Y = F*EXP(ERR(1)) ; DV = CP IF CMT=1 (OR 0)
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.REFERENCES: Guide IV, section V.C.6 REFERENCES: Guide VI, section IV.B.2 , VII.C.3
49 PK PD Simultaneous 2 Example
This is an example for simultaneous analysis of PK and PD. Kinetics are one-compartment. (See PK_PD_simultaneous_1, PK_PD_simultaneous_3).
$PROBLEM SIMULTANEOUS PK/PD - ADVAN6
; NOTE: EXAMPLE USES 1-COMPT PK;
; MORE COMPLEX COMPARTMENTAL MODEL CAN BE USED
; CMT = 1 (OR 0) FOR DV = CP; CMT = 2 FOR DV = EFFECT
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINE ADVAN6 TOL=5
$MODEL
COMP=(CENTRAL,DEFDOSE,DEFOBS)
COMP=EFFECT
$PK
K10= THETA(1)*EXP(ETA(1))
V1= THETA(2)*EXP(ETA(2))
S1= V1
K20= THETA(3)*EXP(ETA(3))
C50= THETA(4)*EXP(ETA(4))
EMAX=THETA(5)*EXP(ETA(5))
$DES
DADT(1)=-K10*A(1) ;NOTE NO LOSS TO A(2)
DADT(2)= K20*(A(1)/V1-A(2)) ;A(2) IS CE
$ERROR
IF (CMT.EQ.2) THEN
Y = EMAX*F/(C50+F)*EXP(ERR(2)) ; DV = EFFECT IF CMT=2
ELSE
Y = F*EXP(ERR(1)) ; DV = CP IF CMT=1 (OR 0)
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.REFERENCES: Guide IV, section V.C.6 , V.C.7 REFERENCES: Guide VI, section IV.B.2 , VII.C.6
50 PK PD Simultaneous 3 Example
This is an example for simultaneous analysis of PK and direct PD. Kinetics are two-compartment. Either ADVAN5 or ADVAN7 could be used. (See PK_PD_simultaneous_1, PK_PD_simultaneous_2).
$PROBLEM SIMULTANEOUS PK/PD - ADVAN7
; NOTE: MODEL FOR CP CAN BE MORE COMPLEX BY ADDING CMPTS
;
; THE DATA FILE CONTAINS BOTH CP AND EFFECT OBSERVATIONS.
; WHEN DV IS A CP OBSERVATION, CMT = 1 (OR 0),
; WHEN DV IS AN EFFECT OBSERVATION, CMT =3.
$DATA data
$INPUT ID TIME DV AMT=DOSE CMT
$SUBROUTINES ADVAN7
$MODEL
COMP=(CENTRAL,DEFDOSE,DEFOBS)
COMP=PERIPH
COMP=EFFECT
$PK
K10= THETA(1)*EXP(ETA(1))
K12= THETA(2)*EXP(ETA(2))
K13= .001*K10 ; TRIVIAL LOSS TO EFFECT COMPT
K21= THETA(3)*EXP(ETA(3))
S1= THETA(4)*EXP(ETA(4)) ; V1 FOR DRUG
K30= THETA(5)*EXP(ETA(5)) ; KEO
EMAX=THETA(6)*EXP(ETA(6))
C50= THETA(7)*EXP(ETA(7))
S3=S1*K13/K30 ; SO THAT CESS = CPSS
$ERROR
IF (CMT.EQ.3) THEN
Y =EMAX*F/(C50+F)*EXP(ERR(2)) ; CMT = 3; DV = EFFECT
ELSE
Y = F*EXP(ERR(1)) ; CMT = 1 (OR 0); DV = CP
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.REFERENCES: Guide IV, section V.C.6 REFERENCES: Guide VI, section IV.B.2 , VII.C.7
51 PK PD Simultaneous 4 Example
This is an example for simultaneous analysis of PK and indirect PD, where the PREDPP PK prediction is modified by $ERROR code.
$PROB WARFARIN
$INPUT ID TIME DV AMT CMT
; DV = DRUG CONC OR EFFECT (PROTHROMBIN TIME)
$DATA data
$SUBROUTINE ADVAN6 TOL=5
$MODEL
COMP=(CENTRAL,DEFDOSE,DEFOBS)
; CAN ADD ADD'L COMPARTMENTS IF PK MULTICOMPARTMENTAL
COMP=(CLOTFAC)
$PK
K10= THETA(1)*EXP(ETA(1)) ;LOSS RATE DRUG
S1= THETA(2)*EXP(ETA(2)) ;VD DRUG
SYNTH=THETA(3)*EXP(ETA(3)) ;SYNTH RATE OF CLOTING FACTOR
LOSS= SYNTH
C50= THETA(4)*EXP(ETA(4))
$DES
DADT(1)=-K10*A(1)
EFF = C50/(C50+A(1)/S1)
DADT(2)= SYNTH*EFF - LOSS*A(2) ;CF LEVEL
$ERROR
PT = THETA(5)/(THETA(6)+F) ; HYPERBOLIC MODEL FOR PRO-TIME
IF (CMT.EQ.2) THEN
Y = PT*EXP(EPS(2)) ; OBS IS PRO-TIME
ELSE
Y = F*EXP(EPS(1)) ; OBS IS DRUG CONCENTRATION
ENDIF
; ADD APPROPRIATE $THETA, $OMEGA, $SIGMA, $ESTIMATION, ETC.The system should be initialized with a bolus dose of AMT=1 into compartment 2 (CLOTFAC) at time zero to establish the initial condition A(1)=0, A(2)=1 at time zero. This is appropriate because, in the absence of drug, the system is assumed to be at steady-state (DADT(2)=0).
REFERENCES: Guide IV, section V.C.6 , V.C.7 REFERENCES: Guide VI, section IV.B.2 , VII.C.6
52 Plasma Urine Example
Here are several examples of simultaneous modelling of plasma and urine data. The examples are adapted from a model for the analysis of Dextromethorphan data provided by Khaled Mohammed Abduljalil. The first involves data from the parent drug. The second involves both parent and metabolite data. In the control streams of these examples plasma and urine data are simulated. The same control streams could also be used for analysis of the data. MU_ modelling is used for the simulation to maintain consistency with subsequent analysis runs but MU_ modelling has no effect on simulation. If MU_ variables are renamed to avoid the character “_” (e.g., rename MU_1 as MU1), the same control streams may be run with NONMEM VI.
EXAMPLE 1 sim_parent_a2.ctl
The data dextroparent.dat contains data from only the parent drug. The parent drug is elimated directly to urine, and also via conversion to metabolite, but there are no measurements of metabolite in plasma or urine. The fraction F0 of parent drug in the urine is computed as in Guide V, Chapter 4, Section 3.1 equation (4.6) and Section 3.4.1 equation (4.7). (See Guide V, Chapter 4.3)
The L2 data item is not required, but is used because multiple observations at the same TIME are in effect multivariate observations. Values of SIGMA corresponding to unused elements EPS(3) and EPS(4) are defined so that the EPS vector has the same length in both examples. See Remarks below.
Although ADVAN2 is used, the same control stream could be implemented with ADVAN5 or ADVAN6 or ADVAN7, with minor changes.
$PROBLEM Parent drug, using ADVAN2
$INPUT ID TIME AMT UVOL DV CMT MDV EVID L2
$DATA dextroparent.dat IGNORE=#
$SUBROUTINES ADVAN2 TRANS1
$PK
K12=THETA(1)
KA=K12
MU_1=LOG(THETA(2))
V2=EXP(MU_1+ETA(1))
MU_2=LOG(THETA(3))
CLP=EXP(MU_2+ETA(2)) ; RENAL CL FOR PARENT
CLB=THETA(4) ; METABOLIC CL FOR METABOLIC
K23=CLP/V2+CLB/V2
K=K23
F0=CLP/(CLP+CLB)
S2=V2
S3=UVOL
$ERROR (EVERY EVENT)
ACMT=ABS(CMT) ; output compt. may have neg. value of CMT
IF(ACMT.EQ.2) Y=F*(1+EPS(1))
IF(ACMT.EQ.3) Y=F*(1+EPS(2))
$THETA
(0.01,0.8,6) ;KA
(0.01,43,1000);V2
(0.0001,20,190);CLP
(0.01,15,90);CLB
$OMEGA
0.05 0.05
$SIGMA
.01 .01
.01 .01 ; eps(3), eps(4) for consistent EPS with sim_metab_a6 example
$SIM (111111) ONLYSIM
$TABLE ID TIME AMT UVOL DV SIMP=PRED CMT MDV EVID L2 NOAPPEND FILE=sim_parent_a2.tabA fragment of the data file is as follows:
# ID TIME AMT UVOL DV CMT MDV EVID L2
2.000 0.000 30000.000 0.000 0.000 1.000 1.000 1.000 1
2.000 0.171 0.000 0.000 0.000 3.000 1.000 2.000 1
2.000 2.000 0.000 0.000 0.000 2.000 0.000 0.000 1
2.000 2.000 0.000 93.100 0.000 -3.000 0.000 0.000 1
2.000 2.000 0.000 0.000 0.000 3.000 1.000 2.000 2
2.000 3.150 0.000 0.000 0.000 2.000 0.000 0.000 2
2.000 3.150 0.000 133.600 0.000 -3.000 0.000 0.000 2There is a dose at TIME=0. At TIME=.171, the urine compartment (which is the default compartment for output) is turned on. At TIME=2, there are observations of both plasma and urine. The urine amount is reset to 0 and the compartment turned back on for the start of a new urine collection. At TIME 3.15 there are observations of plasma and urine.
In the model, different EPS variables are used for the two types of observations. In the table, values of PRED (SIMP) are computed with ETA and EPS equal to 0. The DV values are computed with simulated ETA and EPS, and would be used as the DV in subsequent analysis runs.
EXAMPLE 2 sim_metab_a6.ctl
The data dextrometab.dat includes measurements of metabolite in plasma and urine. The model also include predictions of the metabolite in both plasma (CMT=4) and urine (CMT=5; the default compartment for output). Different compartments are used for parent and metabolite in urine, but the same value of UVOL applies to both. The F0 parameter is not needed because the model for K23 now uses only CLP, the clearance of the parent drug. No analytic ADVAN model has two output compartments, so a general ADVAN (ADVAN6 or ADVAN5 or ADVAN7) must be used.
Note that Compartment 3 is no longer the default for output. Instead of being computed by PREDPP using mass balance, it is computed explicitly by the ADVAN routine itself (e.g., using a differential equation when ADVAN6 is used). If compartment 3 is defined in $MODEL with default attributes, e.g.,
COMP=(DEXURIN)there is an error messsage from PREDPP:
SPECIFIED COMPARTMENT MAY NOT BE TURNED OFF WITH AN OBSERVATION RECORDInstead, the compartment may be defined as an output-type compartment:
COMP=(DEXURIN INITIALOFF NODOSE)Now observations with CMT=-3 are permitted and no change to the original data records is needed.
$PROBLEM Parent drug and metabolite, using ADVAN6
$INPUT ID TIME AMT UVOL DV CMT MDV EVID L2
$DATA dextrometab.dat IGNORE=#
$SUBROUTINES ADVAN6 TRANS1 TOL=4
$MODEL
COMP=(DEPOT)
COMP=(PLASMA DEFOBS) ; PARENT IN PLASMA
COMP=(DEXURIN INITIALOFF NODOSE) ;PARENT IN URINE
COMP=METAB ;METABOLITE IN PLASMA
$PK
K12=THETA(1)
MU_1=LOG(THETA(2))
V2=EXP(MU_1+ETA(1))
MU_2=LOG(THETA(3))
CLP=EXP(MU_2+ETA(2)) ; RENAL CL FOR PARENT
CLB=THETA(4) ; METABOLIC CL FOR METABOLIC
CLMR=THETA(5) ; RENAL CL FOR METABOLITE
V4=1
K24=CLB/V2
K23=CLP/V2
; F0=CLP/(CLP+CLB) Omit F0 because parent and
; metab have different urine compts.
K45=CLMR/V4
S2=V2
S4=V4
S3=UVOL
S5=UVOL
$ERROR (EVERY EVENT)
ACMT=ABS(CMT) ; output compart. may have neg. value of CMT
IF(ACMT.EQ.2) Y=F*(1+EPS(1))
IF(ACMT.EQ.3) Y=F*(1+EPS(2))
IF(ACMT.EQ.4) Y=F*(1+EPS(3))
IF(ACMT.EQ.5) Y=F*(1+EPS(4))
$DES
DADT(1)=-K12*A(1)
DADT(2)=K12*A(1)-K23*A(2)-K24*A(2)
DADT(3)=K23*A(2)
DADT(4)=K24*A(2)-K45*A(4)
$THETA
(0.01,0.8,6) ;KA
(0.01,43,1000);V2
(0.0001,20,190);CLP
(0.01,15,90);CLB
(0.0001,5,90);CLMR
$OMEGA
0.05 0.05
$SIGMA
.01 .01 .01 .01
$SIM (111111) ONLYSIM
$TABLE ID TIME AMT UVOL DV SIMP=PRED CMT MDV EVID L2 NOAPPEND FILE=sim_metab_a6.tabA fragment of the data dextrometab.dat is:
# ID TIME AMT UVOL DV CMT MDV EVID L2
2.000
2.000 0.000 30000.000 0.000 0.000 1.000 1.000 1.000 1
2.000
2.000 0.171 0.000 0.000 0.000 3.000 1.000 2.000 1
2.000
2.000 0.171 0.000 0.000 0.000 5.000 1.000 2.000 1
2.000
2.000 2.000 0.000 0.000 0.000 2.000 0.000 0.000 1
2.000
2.000 2.000 0.000 0.000 0.000 4.000 0.000 0.000 1
2.000
2.000 2.000 0.000 93.100 0.000 -3.000 0.000 0.000 1
2.000
2.000 0.000 93.100 0.000 -5.000 0.000 0.000 1
2.000 0.000 0.000 0.000 3.000 1.000 2.000 2
2.000 0.000 0.000 0.000 5.000 1.000 2.000 2
3.150 0.000 0.000 0.000 2.000 0.000 0.000 2
3.150 0.000 0.000 0.000 4.000 0.000 0.000 2
3.150 0.000 133.600 0.000 -3.000 0.000 0.000 2
3.150 0.000 133.600 0.000 -5.000 0.000 0.000 2Observations of metabolite in plasma (CMT=4) and metabolite in urine (CMT=5) are present at the same values of TIME as observations of parent in plasma (CMT=2) and parent in urine (CMT=3), respectively, and have the same values of L2.
Remarks:
After the control streams are run, the predictions SIMP for compartments 2 and 3 are the same in the two table files sim_parent_a2.tab and sim_metab_a6.tab, as expected, because they are computed with EPS=0, i.e., they are not simulated values. Because of the use of L2 and the same-length EPS vector in both examples, the simulated DV’s for compartments 2 and 3 are also the same. This is useful for illustrative and debugging purposes, but is not usually the case when the EPS structure is different, or the L2 data item is not used.
To summarize: With L2, during simulation SIMEPS generates a new set of values for EPS only when L2 changes value. Without L2, during simulation SIMEPS generates a new set of values for EPS with every record, even when MDV is 1 or the PREDPP item EVID is not 0.
(See Section 59.18).
All control streams and data files are found in the NONMEM 7.4 examplesdirectory.
dextroparent.dat (also displayed here)
sim_parent_a2.ctl (also displayed here)
sim_parent_a5.ctl
sim_parent_a6.ctl
dextrometab.dat (also displayed here)
sim_metab_a5.ctl
sim_metab_a6.ctl (also displayed here)REFERENCES: None.
53 RECORDS=ID Example
These examples illustrate the use of the RECORDS=ID option of the $DATA record. (See Section 50.9).
Example 1 - Data checkout
During a checkout phase, one might want to see only the predictions for the first individual record, and therefore, use a data set comprised only of the data records therein. With RECS=ID, it is not necessary to count the number of records in the fist individual record.
$DATA filename RECS=IDExample 2 - Data Analysis
One might want to analyze each subject’s data independently of the others. One way this can be done is by creating a single control stream with as many problems as there are individual records in the data, and including one each of the following series of $DATA records in the series of problem specifications:
$DATA filename RECS=n1
$DATA filename RECS=n2 NOREWIND
$DATA filename RECS=n3 NOREWIND
etc.where “nk” is the number of records in the kth individual record. It would be necessary to count the exact number of records. Now, these records can simply be coded
$DATA filename RECS=ID
$DATA filename RECS=ID NOREWIND
$DATA filename RECS=ID NOREWIND
etc.Example 3 - Use of the INCLUDE record
With the INCLUDE record, the filename may now be followed by an integer n, whose default value is 1. Then NM-TRAN reads n copies of the named file.
In the example above if the problem specification for all subjects after the first are completely identical, a compact way of writing the control stream is possible. The control stream would contain the problem specification for the first individual, including the $SUBROUTINES record, abbreviated code, and the $DATA record for the first subject, and would end with (if there are e.g. 12 subjects in all)
INCLUDE ctlfile2 11The file ctlfile2 would contain one problem specification (including $PROBLEM) for one subject, with no $SUBROUTINES record or abbreviated code, but would include
$DATA filename RECS=ID NOREWINDA fully worked-out example is supplied as recid2.exa (See recid2.exa). An alternate approach for single-subject data is available with NONMEM 7. OMEGA diagonal values are fixed to a special value 1.0E+06. See Guide Introduction_7 “Single-Subject Analysis using Population with Unconstrained ETAs”.
REFERENCES: Guide Introduction_7
54 RECORDS=ID Example 2
This fully-worked out example shows how the RECORDS=ID option of $DATA and the INCLUDE statement may be used to analyze the Theophyline data set THEOPP as separate individuals. The NONMEM control file is as follows.
$PROB THEOPHYLLINE POPULATION DATA; Analysis of Individuals
; Modification of CONTROL5 control steam
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP RECS=ID
;RECS=ID: Data set will be read until ID changes or end-of-file
$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)
K=THETA(2)
CL=THETA(3)*WT
SC=CL/K/WT
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA .4
;For single subject data OMEGA is residual variance.
$ERROR
Y=F+ERR(1)
;ERR must be used instead of EPS.
$EST MAXEVAL=450 PRINT=5
$COV SPECIAL
;SPECIAL is required to obtain the variance-covariance matrix for
;single-subject data.
$TABLE ID DOSE WT TIME NOPRINT ONEHEADER FILE=indest.tab
$TABLE ID KA K CL SC NOPRINT FIRSTONLY ONEHEADER NOAPPEND
FILE=indest.par
INCLUDE indest.txt 11
; INCLUDE: Inserts copies of the file named indest.txt for each
; additional individual.The contents of file indest.txt are:
$PROB THEOPHYLLINE POPULATION DATA; Analysis of Individuals
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP RECS=ID NOREWIND
;NOREWIND: data set will be read starting after the previous individual
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
$OMEGA .4
;For single subject data OMEGA is residual variance
$EST MAXEVAL=450 PRINT=5
$COV SPECIAL
;SPECIAL is required to obtain the variance-covariance matrix for
;single-subject data
$TABLE ID DOSE WT TIME NOPRINT ONEHEADER FORWARD FILE=indest.tab
$TABLE ID KA K CL SC NOPRINT FIRSTONLY ONEHEADER FORWARD NOAPPEND
FILE=indest.parAn alternate approach for single-subject data is available with NONMEM 7. OMEGA diagonal values are fixed to a special value 1.0E+06. See Guide Introduction_7 “Single-Subject Analysis using Population with Unconstrained ETAs”.
REFERENCES: Guide Introduction_7
55 Repetition 1 Example
This example illustrates the use of the repetition feature, and in particular, the use of the variables RPTI and RPTO. For a discussion of these reserved variables: (See Section 53.94).
We suppose that a compartment amount A2 at time t1 is given by a convolution, i.e. an integral with respect to s, of the product D(t1-s) and r(s), where D is a unit distribution function, and r is an input function. The integral is taken from 0 to t1. The amount in question can be computed by integrating the differential equation
dA2/ds=D(t1-s)r(s)from time 0 to time t1. For A2 at a later time t2, one integrates a different differential equation:
dA2/ds=D(t2-s)r(s)from time 0 to time t2. This equation differs from the first in so far as the constant t2, rather than t1, is used in D.
Ordinarily, PREDPP will use a system of differential equations
dA2/ds=f(A2,t)to obtain the value of a vector A2 of compartment amounts at times t1 and t2, by first integrating the system to time t1, and then by integrating the same or different system from time t1 to time t2. This scheme cannot be used with the two equations described above, because the second equation - with the constant t2 - applies throughout the time interval 0 to t1, as well as throughout the interval t1 to t2.
PREDPP requires the event records to be ordered chronologically in the data set, and with its ordinary usage with differential equations (described above), these records are processed sequentially, starting with the first event record (which, for the purposes of this discussion, can be thought to have time 0). The records are passed sequentially to the PK routine.
The differential equations used for the convolution requires something special. During a regular sequential pass of the event records to the PK routine, when a record R with time t - at which compartment amounts are desired - is initially processed and passed to the PK routine for the first time, the routine should set the output variable RPTO to -1. This signals to NONMEM that the regular pass should be interrupted and a second pass of the records should be initiated, starting with the first event record once again and ending with record R (thus the pass of these records is “repeated”). Only at the end of the second pass, after all event records up to and including record R have once again been passed, will the value for A2(t) be established. The amounts at times on preceeding records have already been established, and they remain unchanged as the pass of these records is repeated. After R has been passed for the second time, the first pass will be continued, the next event record will be passed for the first time, and RPTO can again be set to -1.
Whenever an event record is passed for the first time, f(s) should be set to D(t-s)r(s), where t is the time on the record, so that during the second pass, the differential equation that is used is the desired one. When an event record is passed, if the input variable to the PK routine, RPTI, has the value 0, this signals that the record is not being passed as a part of a repeated pass, and so this can be used as a signal that the record is being passed for the first time.
Suppose the kinetic system is a two-compartment system where drug absorption from a depot compartment into a sampling compartment occurs in a first-order manner, and drug elimination from the sampling compartment occurs according to a Weibull unit disposition function D(s)=exp(-K(s**c)). Suppose too that the superposition principle is assumed to hold. So drug amount is given by a convolution. Suppose the event records for an individual are simply structured, e.g.
# TIME AMT DV EVID
1
.0 10 . 0 rec 1
0 rec 2
.25 . 2.44 0 rec 3
0 rec 4
.5 . 5.24 rec 5
.75 . 5.57
1 . 5.85Then the control stream should look like:
$PROB
$INPT TIME AMT DV EVID
$DATA data
$SUB ADVAN6 TOL=5
$MODEL COMP=(DEPOT DEFDOSE) COMP=(CENTRAL DEFOBS)
$INFN ;enables use of repetition feature
IF (ICALL.EQ.0) RPTO=1
$PK
KA=THETA(1)*EXP(ETA(1))
K =THETA(2)*EXP(ETA(2))
C =THETA(3)
V =THETA(4)*EXP(ETA(4))
S2=V
IF (RPTI.EQ.0) TI=TIME
IF (NEWIND.EQ.2) RPTO=-1
$DES
DADT(1)=-KA*A(1)
D=EXP(-K*(TI-T)**C)
DADT(2)=D*KA*A(1)
etc.Note that when TIME= 0.25, RPTO need not have been set to -1; after the second event record has been passed for the first time, the kinetics are advanced from time 0 to time 0.25 using TI=0.25.
- The above is supplied as repeat1.ctl and repeat1.dat in the examples directory
The code
IF (ICALL.EQ.0) RPTO=1 ;enables use of repetition featureis needed to inform NONMEM that the repetition feature will be used.
The code
IF (NEWIND.EQ.2) RPTO=-1initiates a pass starting with the first event record, which is the default with RPTO=-1.
- Another example, repeat1t.ctl, demonstates that a repeated pass may start with any designated record, not necessarily the first record. The statement
IF (NEWIND.NE.2) RPTO=1is used to indicate explicitly that the first record is a repetition base.
- Example repeat1s.ctl uses the
RPT_data item for the same purpose. This is repeat1s.ctl and repeat1s.dat in the examples directory. A non-zero value ofRPT_marks a record as a repetition base. A statement settingRPTO=-n(where n is the value of theRPT_data item) requests repetition starting with the marked record. WhenRPT_is used, it automatically enables use of the repetition feature and the$INFNblock may be omitted.
The data is augmented with the RPT_ item, as in this fragment:
# TIME AMT DV EVID RPT_
1 .0 10 . 11 rec 1
1 .25 . 2.44 0. rec 2and the control stream contains
$INPT ID TIME AMT DV EVID RPT_- A different example of the use of
RPT_is given in repeatf.ctl and repeatf.dat. Reserved variablesRPTO=nandRPTO=-nare intended to mark and repeat a sequence of more than 1 data records. If only one record is to be repeated,RPT_must be used. The example computes the factorial of a value supplied asKRPTin the data record. A fragment of the data is
ID TIME DV MDV RPT_ KRPT
11 1 0..
12 2 014The code sets
IF (RPT_.EQ.1) RPTO=-1
RPTON=KRPTNONMEM calls PRED KRPT+1 times for each record with RPT_=1. This allows the factorial to be computed. This is intended as an illustration of what can be done using RPT_ data item. With NONMEM 7.3, factorials can be computed more quickly using the DOWHILE feature and the GAMLN function.
REFERENCES: none.
56 Repetition 2 Example
This example illustrates the use of the repetition feature, and in particular, the use of the PRED output control flag PRDFL. For a discussion of the reserved variables RPTO, RPTI and PRDFL: (See Section 53.94).
This example computes a variable RS to be output in a table, defined as the difference between the observation in the data record and the average observation in the individual record. User-defined tabled quantities are computed during a copying pass through the data set - signaled by COMACT=1 (See Section 53.13). In this example, during the copying pass, the data records in each individual record are actually passed twice; the average is computed during a first pass through the individual record - signaled by RPTI=0, and the difference is computed during a second (repeated) pass - signaled by RPTO=1.
$PRED
Y=...
IF (COMACT.EQ.0) PRDFL=1
IF (ICALL.EQ.0) THEN
RPTO=1 ;initiates use of repetition feature
PRDFL=1 ;initiates use of the PRED output control flag
N=0
AV=0
ENDIF
IF (COMACT.EQ.1) THEN
IF (RPTI.EQ.0.AND.NEWIND.LE.1) THEN
N=0
AV=0
ENDIF
IF (RPTI.EQ.0.AND.MDV.EQ.0) THEN
N=N+1
AV=AV+DV
ENDIF
IF (RPTI.EQ.1) THEN
RS=0
PRDFL=1
ENDIF
IF (RPTI.EQ.1.AND.MDV.EQ.0) RS=DV-AV
IF (NDREC.EQ.LIREC) THEN
RPTO=-1 ;RPTO=-1 is a signal initiating
;a second pass. As RPTON isn't also set,
;there are no further passes.
AV=AV/N
ENDIF
ENDIF
$TABLE RS NOAPPENDThe variable NDREC contains the index of the current data record within the current individual record, and LIREC contains the index of the last data record within the current individual record. (See Section 53.92, Section 53.92) (See Section 53.101).
If a data record does not contain an observation (i.e. MDV=1), this code sets RS to 0.
If a variable is to be used only during a copying pass, but recursively across data records (such as N and AV), then it should be defined in an initialization block (e.g. at ICALL.EQ.0). Setting PRDFL=1 during a noncopying pass, assures that during such passes, PRED output will be accepted in the usual way: data record by data record.
Setting PRDFL=1 during a copying pass, but only with the second pass through the data records, assures that PRED output (which includes the RS variable) will be accepted data record by data record, but only during the repeated pass.
REFERENCES: none.
57 Simulation 1 Example
This is an example of simulation with population data. The presence of the $SIMULATION record causes NONMEM to generate DV data items. The data set in file SIMORIG may contain nulls in the DV column.
Table file SIMDATA1 is created. It is similar to SIMORIG, except that the DV column contains the simulated observations. Because of the option NOHEADER of $TABLE, table header lines are not present, and, with the option NOPRINT, the table is not included in the NONMEM output report. SIMDATA1 can be used in subsequent runs as an input file for NONMEM analysis without further modification. These runs can involve different initial estimates, different models, etc.
$PROBLEM Simulation of population data
$INPUT ID TIME WT AMT DV
$DATA SIMORIG IGNORE=#
$SUBROUTINE ADVAN1
$PK
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
K=CL/V
S1=V
$ERROR
Y=F+F*EPS(1)
$THETA .0625 10
$OMEGA .09 .05
$SIGMA .01
$SIMULATION (9215690) ; seed 1-7 digits
$TABLE ID TIME WT AMT NOPRINT FILE=SIMDATA1 NOHEADERContents of SIMORIG (for the first two individuals) follow.
# ID TIME WT AMT DV
1 0 80 100 .
1 1 80 . .
2 0 75 100 .
2 2 75 . .Contents of SIMDATA1 follow, with floating-point notation converted to fixed-point notation for the convenience of the reader. With versions of NONMEM prior to 7.4, this help item contained an incorrect set of values. The correct values are as follows.
1. 0. 80. 100. 0.000 10.000 0.000 0.000
1. 1. 80. 0. 13.646 9.938 3.708 1.531
2. 0. 75. 100. 0.000 10.000 0.000 0.000
2. 2. 75. 0. 8.379 9.876 -1.496 -0.625The last four columns are automatically included by NONMEM in every table (unless the NOAPPEND option is used). They are: DV (simulated values), PRED (calculated using zeros for etas and epsilons), RES, and WRES. The DV, RES, and WRES colunns are zero for dose event records. The last three columns may be ignored (or the DROP option of $INPUT may be used) when the file is subsequently used. Alternately, list DV on the $TABLE record, and include the option NOAPPEND.
It is possible to display simulated values of PRED-defined variables. With NM-TRAN, simulated etas are automatically displayable.
Add the following to the $ERROR block to capture simulated epsilons as ERROR-defined items:
EP1=EPS(1)One may change the $SIMULATION and $TABLE records as follows.
$SIMULATION (9215690) ONLYSIM
$TABLE ID TIME WT AMT CL V ETA1 ETA2 EP1 NOPRINT FILE=SIMDATA1 NOHEADERDisplayed values of CL, V, ETA1, ETA2, and EP1 are simulated values. Without the ONLYSIM option of $SIMULATION, the displayed values are typical values.
(See Simulation example 2, Simulation example 3).
REFERENCES: Guide IV, section III.B.13 , IV.I REFERENCES: Guide V, section 12.4.8 REFERENCES: Guide VI, section III.C , III.E , IV.B , IV.G.1
58 Simulation 2 Example
This is an example of simulation with population data. It is a continuation of Simulation Example 1.
First, during the Simulation Step, values of the covariate WT are simulated, as well as DV values, as PK-defined items SIMWT. Normally-distributed values of ETA(3) are used. Values of CL are calculated based on SIMWT. SIMWT replaces WT in the table file, and, because the ONLYSIM option is used, tabled values of SIMWT are simulated values.
Second, this process is repeated 10 times, using the same model, by means of the SUBPROBLEMS option. The table file contains the concatenated output of all the subproblems, and is the simulated data for 20 individuals.
$PROBLEM Simulation of population data and weight
$INPUT ID TIME WT AMT DV
$DATA SIMORIG IGNORE=#
$SUBROUTINE ADVAN1
$PK
SIMWT=70+70*ETA(3)
CL=THETA(1)*SIMWT*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
K=CL/V
S1=V
$ERROR
Y=F+F*EPS(1)
$THETA .0625 10
$OMEGA .09 .05 .04
$SIGMA .01
$SIMULATION (9215690) ONLYSIM SUBPROBLEMS=10
$TABLE ID TIME SIMWT AMT NOPRINT FILE=SIMDATA2 NOHEADERSIMORIG is the same as with Simulation Example 1. Contents of SIMORIG (for the first two individuals) follow.
# ID TIME WT AMT DV
1 0 80 100 .
1 1 80 . .
2 0 75 100 .
2 2 75 . .Contents of SIMDATA2 follow, with floating-point notation converted to fixed-point for the convenience of the reader.
1. 0. 77.486 100. 0.000 10.000 0.000 0.000
1. 1. 77.486 0. 6.825 6.457 0.369 0.000
2. 0. 72.882 100. 0.000 10.000 0.000 0.000
2. 2. 72.882 0. 3.515 4.169 -0.653 0.000REFERENCES: Guide IV, section III.B.13 , IV.I REFERENCES: Guide V, section 12.4.8 REFERENCES: Guide VI, section III.C , III.E , IV.B , IV.G.1
59 Simulation 3 Example
This is an example of simulation and estimation within the same NONMEM run. It is a continuation of Simulation Example 2.
As in Simulation Example 2, values of the covariate WT are simulated during the Simulation Step. They also replace the original WT values in NONMEM’s internal copy of the data set. This transgeneration makes them available in the Estimation and Table Steps.
There are two models for V, illustrating the use of diferent models in the Simulation and Estimation Steps:
(simulated V) VS = THETA(2) * EXP(ETA(3))
(non - simulated V) VN = THETA(3) + ETA(4)There are two epsilons in the $ERROR block, illustrating the use of an of using an initial estimate for estimation different from the value used to simulate.
Because theta(2), eta(2), eta(3) and eps(1) are used only in the Simulation Step, the FIXED attribute is specified for theta(2), omega(2,2), omega(3,3) and sigma(1,1) during the Estimation Step.
$PROBLEM Simulation of population data w/ covariates
$INPUT ID TIME WT AMT DV
$DATA SIMORIG IGNORE=#
$SUBROUTINE ADVAN1
$PK
IF (ICALL.EQ.4) THEN
WT=70+70*ETA(3)
ENDIF
CL=THETA(1)*WT*EXP(ETA(1))
IF (ICALL.EQ.4) THEN
V=THETA(2)*EXP(ETA(2)) ; Simulated V
ELSE
V=THETA(3)+ETA(4) ; Non-simulated V
ENDIF
K=CL/V
S1=V
$ERROR
IF (ICALL.EQ.4) THEN
Y=F+F*EPS(1)
ELSE
Y=F+F*EPS(2)
ENDIF
$THETA (0,.0625) (10 FIX) (0,12)
$OMEGA .09 (.05 FIX) (.04 FIX) 4
$SIGMA (.01 FIX) .03
$SIMULATION (9215690) SUBPROBLEMS=10
$ESTIMATION
$TABLE ID TIME WT AMT NOPRINT FILE=SIMDATA3 NOHEADERREFERENCES: Guide IV, section III.B.13 , IV.I REFERENCES: Guide V, section 12.4.8 REFERENCES: Guide VI, section III.C , III.E , IV.B , IV.G.1
60 Stieltjes Example
The following are examples of the specification of the Stieltjes method.
$EST PRINT=1 MET=1 STIELTJES
$EST MAX=9999 METH=1 LAPLACE LIKE PRINT=1 STIELTJES
$EST MAX=9999 METH=1 LAPLACE LIKE PRINT=1 STIELTJES GRID=(9,10,.1,.9)
$EST MAX=500 MET=COND LAPLACE LIKE STIELTJES GRID=(20,2,.1,.9)REFERENCES: Guide Introduction_7
61 Stochastic Differential Equation Examples
A number of examples are given of the use of Stochastic Differential Equations.
sde6- Example problem without sde
sde7-
Example problem with SDE, using expanded data set that have repeated lines to re-invoke
PREDseveral times for each data record. sde8-
Use of
XVIDto re-invokePREDseveral times per data record wthout requiring the redundancy of sde7.csv. sde9- Alternative method of sde programming, using add-on sde.f90
They, along with the data files, can be found in the examples directory of the NONMEM 7 distribution medium.
See INTRODUCTION TO NONMEM 7, Repeated Observation Records See INTRODUCTION TO NONMEM 7, Stochastic Differential Equation Plug-In
REFERENCES: Guide Introduction_7
62 Superproblem 1 Example
There are four problems:
Problem 1: The (true) probability that with an individual chosen at random, after a 200 mg dose, the true concentration at 5 hrs is above 5, is computed using simulation. Problem 2: A population data set is simulated. Problem 3: Parameter estimates are obtained. Problem 4: The probability in question is estimated, using simulation and the parameter estimates obtained in problem 3.
Problems 2-4 comprise a superproblem. The bias in the estimator defined by problem 4 is estimated by iterating this superproblem.
$PROB COMPUTE PERCENTAGE OF POP. WITH CP>5 AT DOSE=200, TIME=5
$INPUT ID TIME DOSE DV
$DATA dataB
$PRED
IF (ICALL.EQ.0) THEN
BIAS=0 ; initialize
TRUE=0
OPEN (50,FILE='BIAS') ; Give the output file a name
ENDIF
IF (ICALL.EQ.1) THEN
IF (IPROB.EQ.1.OR.IPROB.EQ.4) N=0 ; see rocm14
CALL SUPP (1,1)
ENDIF
KA=THETA(1)*EXP(ETA(1))
KE=THETA(2)*EXP(ETA(2))
V=THETA(3)*EXP(ETA(3))
A=EXP(-KE*TIME)
B=EXP(-KA*TIME)
C=KA-KE
D=A-B
E=KA*DOSE/(V*C)
F=E*D
Y=F+ERR(1)
IF (ICALL.EQ.4) THEN
IF (IPROB.EQ.1.OR.IPROB.EQ.4) THEN
IF (F.GT.5) N=N+1
IF (IREP.EQ.NREP) THEN ; see rocm10
PER=100.*N/NREP
IF (IPROB.EQ.1) TRUE=PER
IF (IPROB.EQ.4) THEN
BIAS=BIAS+(PER-TRUE)/TRUE
IF (S1IT.EQ.S1NIT) THEN
BIAS=BIAS/S1NIT
WRITE (50,*) BIAS
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
$THTA (.4,1.7,7) (.025,.102,.4) (10,29,80)
$OMEGA .04 .04 .04
$SIGMA 1.5
$SIM (5566898) ONLY SUB=1000
;
$SUPER SCOPE=3 ITERATIONS=10
$PROB SIMULATION
$INPUT ID TIME DOSE DV
$DATA dataA
$THTA (.4,1.7,7) (.025,.102,.4) (10,29,80)
$OMEGA .04 .04 .04
$SIGMA 1.5
$SIM (-1) ONLY
$TABLE ID TIME DOSE DV FILE=simulation NOHEADER NOPRINT NOFORWARD
;
$PROB ESTIMATION
$INPUT ID TIME DOSE DV
$DATA simulation (4E12.0) NRECS=500 NOOPEN
$THTA (.4,1.7,7) (.025,.102,.4) (10,29,80)
$OMEGA .04 .04 .04
$SIGMA 1.5
$EST PRINT=0 MSFO=msf
;
$PROB ESTIMATE PERCENTAGE OF POP. WITH CP>5 AT DOSE=200, TIME=5
$INPUT ID TIME DOSE DV
$DATA dataB
$MSFI msf
$SIM (-1) ONLY SUB=100 TRUE=FINALComments:
Variables whose values are to be retained across problems (e.g.
BIASandTRUE) should be defined in initialization blocks. In this instance these variables are defined atICALL=0, before any of the problems are implementedIREPandNREPare the number of the current replication and number of requested replications for the Simulation Step, respectively.S1ITandS1NITare the number of the current iteration and number of requested iterations for the (level 1) superproblem, respectively.Problems 2 & 3 can, of course, be combined, making the control stream much simpler. But separation can sometimes be useful, and here we illustrate how this can be done.
The data set of problem 3 can be taken to be the internal data set created with problem 2 simply by removing both the
$TABLErecord of problem 2 and the$DATArecord of problem 3. But these records are included here to illustrate the ability to use a table file as a data file in a subsequent problem of the same run. When so doing, all the items in the$INPUTrecord of the subsequent problem should be included in the$INPUTrecord of the problem with the$TABLErecord. (If aPRED-item in the table is to be used in the subsequent problem as a data item, thePRED-item name can be different from the data item name.) TheFORMATandNOOPENoptions should be included in the$DATArecord of the subsequent problem. TheNRECSoption is not needed, but if it is also included, it should be set to the number of records in the table file, i.e. to the number of data records in data set A (assumed here to be 500).
Contents of dataA (for individual 1 only) follow:
1 .27 320. .
1 .52 320. .
1 1. 320. .
1 1.92 320. .
1 3.5 320. .
1 5.02 320. .
1 7.03 320. .
1 9. 320. .
1 12. 320. .
1 24.3 320. .This data set cannot be run as-is. File dataA must include more subjects. If not, the estimation step in problem 3 will not terminate correctly. The final problem will produce the message:
THE SIMULATION TASK REQUIRES THAT A BETTER FINAL ESTIMATE
BE AVAILABLE IN THE MODEL SPECIFICATION FILE
OR
MULTIPLE SUB_PROBLEMS IS INCOMPATIBLE WITH OTHER INFORMATION
IN THE CONTROL STREAM AND INPUT MSFThis error message is related to an earlier message
PARAMETER ESTIMATE IS NEAR ITS BOUNDARY
THIS MUST BE ADDRESSED BEFORE THE COVARIANCE STEP CAN BE IMPLEMENTEDA larger version of dataA can be found in the examples directory of the NONMEM distribution medium. The content of dataA above is replicated 28 times. ID values 1 and 2 alternate so that there are 28 subjects.
Contents of dataB follow:
1 5.0 200. .Comments on NONMEM 7 changes:
In this example as distributed with earlier versions of NONMEM, the WRITE statement specified unit 42:
WRITE (42,*) BIASWith NONMEM 7, this must be changed to
WRITE (50,*) BIASThis avoids a possible NONMEM error message
FORTRAN UNIT BEING OPENED OR CLOSED IN PROBLEM 1 IS OUTSIDE
ACCEPTABLE RANGE (50-2000)A new statement is added:
OPEN (50,FILE='BIAS')This causes the output to be written to a file named BIAS rather than fort.50 or the like.
REFERENCES: None.
63 Superproblem 2 Example
Two problems, in order to effect a standard two-stage analysis:
Problem 1: First stage analysis Problem 2: Obtain simple population average of PK parameter.
Problem 2 comprises a superproblem. It is implemented once for CL and once for V, and can be implemented more times when there are more parameters of interest.
$PROBLEM TWO-STAGE STAGE 1
$DATA dataD
$INPUT ID DOSE TIME DV=CP
$PRED
IF (IPROB.EQ.1.AND.NEWIND.LE.1) D=DOSE
IF (IPROB.EQ.1) THEN
CL=THETA(1)*EXP(ETA(1))
V =THETA(2)*EXP(ETA(2))
K=CL/V
F=D*EXP(-K*TIME)/V
Y=F+F*EXP(ETA(3))*ERR(1)
ENDIF
IF (ICALL.EQ.4) THEN
; if S1IT >= 1 then problem 2 is being implemented.
; see comment 1 below.
IF (S1IT.GE.1) THEN
" IF (S1IT.EQ.1) DPV=DATREC(2)
" IF (S1IT.EQ.2) DPV=DATREC(3)
Y=DPV
RETURN
ENDIF
ENDIF
IF (S1IT.EQ.1) Y=THETA(1)+ERR(1)
IF (S1IT.EQ.2) Y=THETA(2)+ERR(2)
$THETA 5 50
$OMEGA 100 100 100
$SIGMA 1 FIX 0 FIX
$ESTIMATION MAXEVAL=0 POSTHOC METH=COND INTERACTION
$TABLE ID CL V FILE=stage2 NOHEADER NOPRINT FIRST NOFORWARD
;
$SUPER SCOPE=1 ITERATION=2
$PROBLEM TWO-STAGE STAGE 2
$INPT ID CLI VI DV
$DATA stage2 (4F12.0) NRECS=12 REWIND NOOPEN
$THETA 5 50
$OMEGA 0 FIX 0 FIX 0 FIX
$SIGMA 1 1
$SIMULATION (1)
$ESTIMATION ;population estimates are obtainedComment:
The simulation step in problem 2 is used to select the DV data item from one of the several possibilities (CL or V) output by problem 1.
Using generated code, data item names that do not appear in the
$INPUTrecord of problem 1, such asCLIandVI, cannot be used in abbreviated or verbatim code.DATRECis a vector-array of data items that can be referenced in verbatim code.
Contents of dataD (for individual 1 only) follow.
1.000 4.020 0.000 0.091
1.000 0.000 0.250 0.084
1.000 0.000 0.570 0.105
1.000 0.000 1.120 0.057
1.000 0.000 2.020 0.038
1.000 0.000 3.820 0.060
1.000 0.000 5.100 0.046
1.000 0.000 9.050 0.017
1.000 0.000 7.030 0.024
1.000 0.000 12.120 0.018
1.000 0.000 24.370 0.002REFERENCES: None.
64 TDM Example
This is an example of how F1 can be used to simplify the preparation of a data file when the dose is changed often, e.g., for Therapeutic drug monitoring (TDM).
Suppose the dosing pattern is consistent (e.g., every 12 hours) but the amount changes often. Suppose also there are observation events at various times.
In the following data set, there is only one dose record having AMT=1. ADDL and II are used to specify an (effectively infinite) number of additional doses, every 12 hours. When the dose amount changes, it is sufficient to supply a record giving the new dose amount in the DOSE data item. Such records have EVID=2. There is no need to interrupt the dosing pattern for an observation event, or to count the number of doses between observation events.
Suppose a fragment of the dose records looks like this:
ID DAT1=DROP TIME EVID DV DOSE AMT II ADDL
99999
1 01.01.2008 20 1 0 4 1 12
0
1 03.01.2008 8 0 6.1 0 00 0
0
1 04.01.2008 8 0 7.6 0 00 0
1 04.01.2008 20 2 0 3 00
1 06.01.2008 8 0 5.4 0 00Here is a fragment of the control stream:
$INPUT ID DAT1=DROP TIME EVID DV DOSE AMT II ADDL
...
$PK (ADDITIONAL OR LAGGED)
...
IF (DOSE> 0) SAVEDOSE=DOSE
F1=SAVEDOSEIn $PK, the value of DOSE is saved, and used on subsequent calls to PK as the value of F1. Although F1 is usually used to model the bio-availability of doses to compartment 1, in this example it is the actual dose amount.
REFERENCES: Guide V, section 7.4.4
65 Three Compartment Example
This example shows how to implement a three Compartment Linear Model with First-order Absorption using PREDPP, where elimination occurs from a central compartment and one of two peripheral compartments.
As with any ADVAN, one or more doses could be input to any of the compartments. Steady state doses are also possible.
The use of the DEFDOSE attribute on the following $MODEL record specifies that if the CMT data item is null on a given dose event record, then the dose(s) specified by that record are to be placed in the depot compartment. The use of the INITIALOFF attribute specifies that the depot compartment is to be off until a dose is placed in it.
The use of the DEFOBS attribute on the $MODEL record specifies that if the CMT data item is null on a given observation event record, then the quantity F which is available as a right-hand quantity in $ERROR is to be the scaled drug amount in the central compartment. (See Section 51.4).
$PROBLEM EXAMPLE OF A THREE COMPARTMENT MODEL WITH ABSORPTION
$INPUT ID TIME AMT DV CMT
$DATA datafile
$SUBROUTINES ADVAN7
$MODEL COMP = (DEPOT,DEFDOSE,INITIALOFF) COMP = (CENTRAL,DEFOBS) COMP = PERIPH1, COMP = PERIPH2
$PK
K12 = THETA(1)*EXP(ETA(1)) ; depot to central
K23 = THETA(2) ; central to periph1
K32 = THETA(3) ; periph1 to central
K30 = THETA(4)*EXP(ETA(2))
K24 = THETA(5) ; central to periph2
K42 = THETA(6) ; periph2 to central
K20 = THETA(7)*EXP(ETA(3)) ; elimination from central
S2 = THETA(8)*EXP(ETA(4)) ; scale for central
;There are other reasonable ways that etas can be assigned.
;Elaborate these expressions for the K's as is appropriate.
;Add $ERROR, parameter and task records as is necessary/desired.Other parameterizations are possible, as long as the above set of micro-constants is defined. E.g., the definitions of K20 and S2 can be replaced by the following.
CL = THETA(7)*EXP(ETA(3))
V = THETA(8)*EXP(ETA(4))
K20 = CL/V
S2 = VIf the depot compartment is not needed, it can be omitted. The numberings of the compartments, the thetas and etas must then be adjusted.
REFERENCES: Guide VI, section VII.C.5 , VII.C.7
66 Time After Dose (TAD) Example
This fully-worked out example shows how Time After Dose (TAD) may be computed in $PK abbreviated code. It works with transient and steady-state doses, and also with additional and lagged doses.
It is based on suggestions from the NONMEM tips directory but is not identical to any of them. See Note 8 for a discussion of the tips.
Two versions of the code are given.
- VERSION I
-
Code to compute
TADE(“TAD Effective”). The time that the dose enters the system is used, which is the time the dose was administered plusALAG. - VERSION II
-
Code to compute
TADA(“TAD Administered”). The time that the dose was administered is used, andALAGis ignored.
If PDxPop is used or any other software that looks for TAD in the table, specify TAD=TADA or TAD=TADE in the $TABLE record.
$PROB TIME AFTER DOSE (Based on RUN# 705 nonmem coding challenge #1)
; This example shows how Time After Dose can be computed.
; It contains two independent versions of the computation.
; Version I computes TADE (Time After Dose Effective)
; If there is absorption lag, TADE is the Time after the lagged dose,
; that is, the time the dose actually enters the sysem.
; Version II computes TADA (Time After Dose Administered)
; Absorption lag is ignored.
;
; As of Dec. 17, 2018
;
$INPUT C ID TIME DV AMT WT AGE CRCL SMK ADDL II EVID DROP
$DATA 706.csv IGNORE=C
$SUBROUTINE ADVAN2 TRANS2
$PK
;
ALAG1=0 ; ALAG1 is 0 for this example, but could be set to some other
; value. With ALAG1=0, either version can be used.
; However, VERSION I IS simpler and is preferred.
;
; VERSION I: CODE TO COMPUTE TADE = TAD EFFECTIVE (USE ALAG)
; PREDPP KEEPS TRACK INITIATING AND IMPLIED DOSES, AND COMPUTES DOSTIM
IF (NEWIND.LT.2.OR.EVID.EQ.3) THEN
TDOSE=-999 ; TIME OF MOST RECENT DOSE. -999 IF NO PREVIOUS DOSE.
TADE=0.0 ; TIME AFTER DOSE EFFECTIVE
ENDIF
IF (EVID.EQ.1.AND.ALAG1.EQ.0.OR.EVID.EQ.4.AND.ALAG1.EQ.0) TDOSE=TIME
IF (DOSTIM.GT.0) TDOSE=DOSTIM
IF (TDOSE.GT.-999) TADE=TIME-TDOSE
; If no ALAG and no ADDL doses, the above three lines become:
; IF (EVID.EQ.1.OR.EVID.EQ.4) TDOSE=TIME
; IF (TDOSE.GT.-999) TADE=TIME-TDOSE
;
; End of VERSION I
; VERSION II: CODE TO COMPUTE TADA = TAD ADMINISTERED (IGNORES ALAG)
; IF ALAG>0, PREDPP's DOSTIM CANNOT BE USED FOR TADA
; THE ABBREVIATED CODE MUST KEEP TRACK OF INITIATING AND IMPLIED DOSES.
; THIS CODE IS SIMILAR TO tip4-new-general-01-31-18.txt
IF (NEWIND.LT.2.OR.EVID.EQ.3) THEN
TDOSA=-999 ; TIME OF MOST RECENT DOSE. -999 IF NO PREVIOUS DOSE.
TADA=0.0 ; TIME AFTER DOSE ADMINISTERED
ENDIF
IF (DOSTIM==0) THEN ; IGNORE NON-EVENT DOSE TIMES (DOSTIM>0)
IF(EVID.EQ.1.OR.EVID.EQ.4) THEN ; NEW DOSE EVENT RECORD
TADA=0.0
DIV=II
TDOSA=TIME
TLAST=TDOSA+ADDL*II ; TLAST IS THE TIME OF THE FINAL IMPLIED DOSE
ENDIF
IF (TDOSA.GE.0.AND.EVID.NE.1.AND.EVID.NE.4) THEN ; THERE WAS AN EARLIER DOSE
IF (TIME>TLAST) THEN
TADA=TIME-TLAST ; CURRENT TIME IS PAST THE TIME OF THE FINAL IMPLIED DOS
ELSE
DIFF=TIME-TDOSA
TADA=MOD(DIFF,DIV) ; COMPUTES TIME OF THE MOST RECENT IMPLIED DOSE
ENDIF
ENDIF ; NOT A DOSE
ENDIF ; END OF DOSTIM==0
;
; End of VERSION II
CL=THETA(1)*EXP(ETA(1))
TVV=THETA(2)
V=TVV*EXP(ETA(2))
TVKA=THETA(3)
KA=TVKA*EXP(ETA(3))
S2=V
$THETA
(0, 10)
(0, 100)
(0, 0.5)
$ERROR
Y=F+ERR(1)
IPRED=F
$OMEGA
0.04 ;[P] INTERIND VAR IN CL
0.04 ;[P] INTERIND VAR IN V
0.04 ;[P] INTERIND VAR IN KA
$SIGMA
0.2 ;[A] ADDITIVE COMPONENT
$EST MAXEVAL=0 PRINT=2 NOABORT ;POSTHOC
;
; If PDxPop is used or any other software that looks for TAD in the table, specify
; TAD=TADA or TAD=TADE in the $TABLE record
$TABLE ID TIME EVID TDOSA TADA TDOSE TADE IPRED NOAPPEND NOPRINT ONEHEADER
FILE=tadexa.tab FORMAT=SF11.7The data for the first subject in 706.csv is:
C,Data Desc: test of TAD for ADDL,,,,,,,,,,
C,ID,TIME,DV,AMT,WT,AGE,CRCL,SMK,ADDL,II,EVID,TAD
0,1,0,0,0,58.4,51,4.49,1,0,0,2,0
0,1,0.1,0,500,58.4,51,4.49,1,1,24,1,0
0,1,1.1,1.67538,0,58.4,51,4.49,1,0,0,0,0
0,1,4.1,2.79283,0,58.4,51,4.49,1,0,0,0,0
0,1,10.1,0.583263,0,58.4,51,4.49,1,0,0,0,0
0,1,24.1,0.146255,0,58.4,51,4.49,1,0,0,0,0
0,1,25.1,1.67538,0,58.4,51,4.49,1,0,0,0,0The original data set 704.csv is embedded in tip3-2-25-02.txt. Data set 706.csv can be found in the tips directory. It was created from 704.csv by appending an extra column, TAD, whose values are zero. $INFN code in the tips assign values to the TAD data item (See Note 8). The code in the present example does not use the TAD data item (DROP is listed in $INPUT). Instead, TAD is a user-defined variable listed in MODULE NMPRD4 from which its values may be displayed in a table file.
- Note 1.
-
It is assumed that
TADis displayed in tables for the purpose of graphical display of the data. It is not a random variable, and is not intended to be part of the prediction model. - Note 2.
-
If there are non-dose records prior to the first dose record of the Individual record,
TADis set to zero for these records. This is an arbitrary convention; the user could use different values, especially if these conditions might signal an error in the data set. - Note 3.
-
ALAG1is used in the example, but a different dosing compartment could be used, e.g.,ALAG2for doses into compartment 2. The code would need modification to computeTADfor doses into two or more compartments. - Note 4.
-
This code can handle random
DOSTIM(Eta onALAGor ETA on modelled duration or rate.) - Note 5.
-
The code does not use
INFN. The information is developed with passes through the data set atICALL=2, in which NONMEM calls PREDPP. During a pass in whichINFNis called, PREDPP is not called. The abbreviated code would have to account for additional and lagged dose times. - Note 6.
-
TDOS=-999until the first dose.TDOS=0when the first dose is at TIME 0. This makes it possible to dispense with a flag variable whose only function is to say “a dose has occurred”. If there are negative time values in the data set and any time value is smaller than -999, then - Note 7.
-
The code is intended to compute
TADfor multiple individual doses or one initiating dose. (An initiating dose hasADDL>0,II>0). - Note 8.
- Changes were made to tip3, tip4, and tip4update in the NONMEM tips directory. See tips in
https://nonmem.iconplc.com/nonmem/tipsIn each case, the name of the revised tip contains 01-31-18 and supercedes the original. Although the tips were revised as part of the nonmem 7.4.2 release they apply to previous versions of nonmem.
tip3-revised-01-31-18.txt- The original version, tip3-2-25-02.txt, provides a data file 704.csv and a simple model. The revised version provides improved code, similar to the code in tadexa.exa.
tip4-revised-01-31-18.txt-
Uses Fortran suboutine
INFNand equivalent$INFNcode. The datafile is 706.csv. tip4update-revised-01-31-18.txt-
Uses
DOSTIM.
The original versions are still in the tips directory with new names that include the original dates.
tip3-2-25-02.txt
tip4-3-04-02.txt
tip4update-10-11-13.txtREFERENCES: None
67 TNPRI Examples
These examples illustrate the use of the NONMEM utility routine TNPRI by the NM-TRAN $PRIOR record, and also by an equivalent user-supplied routine PRIOR. For a discussion of PRIOR: (See Section 57.6). For a discussion of the TNPRI: (See Section 59.21).
Example 1
Use of PRIOR and TNPRI and a control stream for population data:
This example obtains parameter estimates from population data, using a two-compartment PK model, and it incorporates a frequentist prior of transformed normal form for all of THETA, OMEGA and SIGMA. The prior information is found in model specification file msf1.
$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID DOSE TIME DV WT
$PRIOR TNPRI (PROBLEM 2)
$MSFI msf1 ONLYREAD
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(1)
$PROB POPULATION DATA WITH PRIOR ON THETA, OMEGA AND SIGMA
$INPUT ID DOSE TIME DV WT
$DATA data REWIND
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
$SIGMA .4
$ESTInstead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=priorThe prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
USE NMPRD_INT, ONLY: IPROB
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
IF (IPROB==2) THEN
PLEV=0.
ITYP=0
NSAM=0
ISS=0
IFND=0
MODE=0
IVAR=0
CALL TNPRI(IFND,MODE,ITYP,PLEV,NSAM,ISS,IVAR,CNT)
ENDIF
RETURN
ENDExample 2
Use of PRIOR and TNPRI and a control stream for population data:
This example obtains parameter estimates from population data, using a two-compartment PK model, and it incorporates a frequentist prior of transformed normal form for all of THETA and OMEGA, but for only the first element of a diagonal SIGMA, which is a prior-specific parameter. Prior information for these parameters is found in model specification file msf1.
$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID DOSE TIME DV WT
$PRIOR TNPRI (PROBLEM 2)
$MSFI msf1 ONLYREAD
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(2)
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data REWIND
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
$SIGMA .4 FIX .5
$ESTInstead of a $PRIOR record, use the same $SUBROUTINE record and the same subroutine prior as in Example 1.
Example 3
Use of PRIOR and TNPRI and a control stream for population data:
This example simulates the THETA and OMEGA parameters of the same two-compartment population PK model, as well as the population data from this model. As with Example 2, the first element of SIGMA is prior-specific, and so its value need not be simulated. The value .5 is used for the value of the second element of SIGMA. The default variance-covariance matrix from the prior problem is used for the hyperparameter variance-covariance matrix (IVAR=1).
$PROB READ THE MODEL SPECIFICATION FILE
$DATA data
$INPUT ID DOSE TIME DV WT
$PRIOR TNPRI (PROBLEM 2) PLEV=.9999 IVAR=1
$MSFI msf1 ONLYREAD
$PRED
;
; THETA(1)=MEAN ABSORPTION RATE CONSTANT - MEAN ELIM. RATE CONSTANT (l/HR)
; THETA(2)=MEAN ELIM. RATE CONSTANT (1/HR)
; THETA(3)=SLOPE OF CLEARANCE VS WEIGHT RELATIONSHIP (LITERS/HR/KG)
; DOSE=WEIGHT-ADJUSTED DOSE (MG/KG)
;
IF (NEWIND.NE.2) THEN
AMT=DOSE*WT
W=WT
ENDIF
T0=THETA(1)*EXP(ETA(1))
T2=THETA(2)*EXP(ETA(2))
T1=T2+T0
T3=THETA(3)*W*EXP(ETA(3))
A=EXP(-T2*TIME)
B=EXP(-T1*TIME)
C=T1-T2
D=A-B
E=T3*C
Y=AMT*T1*T2/E*D+EPS(2)
$PROB POPULATION DATA WITH PRIOR ON THETA AND OMEGA
$INPUT ID DOSE TIME DV WT
$DATA data REWIND
$THETA (0,4,5) (0,.09,.5) (.004,.01,.9)
$OMEGA BLOCK (3) .7 .04 .05 .02 .06 .08
$SIGMA .4 FIX .5
$SIM (547676) ONLY TRUE=PRIORInstead of a $PRIOR record, the following may be used:
$SUBROUTINE PRIOR=priorThe prior routine is as follows:
SUBROUTINE PRIOR (ICALL,CNT)
USE SIZES, ONLY: DPSIZE,ISIZE
USE NMPRD_INT, ONLY: IPROB
REAL(KIND=DPSIZE) :: CNT
REAL(KIND=DPSIZE) :: PLEV
INTEGER(KIND=ISIZE) :: ICALL
IF (IPROB==2) THEN
PLEV=.9999
IVAR=1
ITYP=0
NSAM=0
ISS=0
IFND=0
MODE=0
CALL TNPRI(IFND,MODE,ITYP,PLEV,NSAM,ISS,IVAR,CNT)
ENDIF
RETURN
ENDREFERENCES: none.
68 Transit Compartment Examples
A number of examples are provided for the implementation of dose superimposition into a transit compartment.
BACKGROUND:
These examples are based on the following two papers:
Implementation of a transit compartment model for describing drug absorption in pharmacokinetic studies J Pharmacokinet Pharmacodyn (2007) 34:711-726 DOI 10.1007/s10928-007-9066-0 Radojka M. Savic, Daniel M. Jonker, Thomas Kerbusch, Mats. O. Karlsson
Implementation of dose superimposition to introduce multiple doses for a mathematical absorption model (transit compartment model) J Pharmacokinet Pharmacodyn (2012) 39:251-262 DOI 10.1007/s10928-012-9247 Jun Shen, Alison Boeckmann & Andrew Vick
The first paper defines what a transit compartment is, and gives a model for a single dose into such a compartment. The second paper gives several ways that multiple dosing with dose superimposition can be modelled with NONMEM 7.
Files on ftp site and NONMEM 7
The files for the second paper are provided on the IDS ftp site:
ftp://nonmem.iconplc.com/Public/nonmem/transit_compartment/The same files are also provided with NONMEM 7 in directory
examples\transit_compartmentThe general solution (Appendix 4; sumdosetf.ctl) is of particular interest in that it gives an example of the use of vectors (e.g., VECTRA(n), VECTRB(n)) and user-written functions (FUNCA, FUNCB) in abbreviated code.
NEW FILES FOR NONMEM 7:
The following files may be found in the examples directory.
sumdosetn.ctl, sumdosetf.dat- With NONMEM 7, it is also possible to code the general solution directly in abbreviated code, using the DO WHILE statement, user-declared arrays, and subscripted variables.
Note that data file sumdosetf.csv is identical to data file sall6.csv on the ftp site. As noted in the INTRODUCTION TO NONMEM 7, this includes abbreviated code for summation such as:
$abbr declare dosetime(100),dose(100)
$abbr declare dowhile i
$DES
INPT=0
I=1
DOWHILE (I<=NDOSE)
IPT=0
IF (T>=dosetime(I)) &
IPT=DOSE(I)*(T-dosetime(I))**NN*EXP(-KTR*(T-dosetime(I)))
INPT=INPT+IPT
I=I+1
ENDDOssmultidose.ctl, ssmultidose.dat- This is an example of a steady state reached using multiple transient bolus doses. The abbreviated code is the same as sumdosetn.ctl, but the data file ssmultidose.dat has nine identical bolus doses, at 12 hour intervals, leading to steady state (to the number of significant digits printed).
ssonedose.ctl, ssonedose.dat-
This is an example of a steady state computed with a single dose record having
SS=1andII=12. The data file is ssonedose.dat. The$PKblock saves the values ofAMTandIIfrom the dose record. Abbreviated code in the$DESblock contains a DO WHILE loop that adds drug from each implied dose, going back in time till the amount from the implied dose is negligible.
REFERENCES: Guide Introduction_7
69 User-Defined Reserved Function Example
If you wish to define your own function, and have the information about its proper use of arguments be conveyed upon its execution, so the compiler may detect errors, then one method is to package the definition of the function in a USE module. The function must be listed in an include file whose name starts with the characters
nonmem_reservedCaution: NMTRAN will permit the use of such a function but will not compute eta deriatives with respect to the arguments. Do not use such functions to compute any variable that affects the objective function!
Here is an example.
Myfuncmodule.f90 defines the functions mymin and mymax:
MODULE MYFUNCS
contains
function mymin(a,b,c,d,e)
integer mymin
integer a,b,c,d,e
mymin=min(a,b,c,d,e)
end function
function mymax(a,b,c,d,e)
integer mymax
integer a,b,c,d,e
mymax=max(a,b,c,d,e)
end function
END MODULE MYFUNCSNonmem_reserved_myfunc is the include file that declares its use:
" USE myfuncs, only: mymin,mymaxThe following control stream file uses the function:
$PROB THEOPHYLLINE POPULATION DATA
$INPUT ID DOSE=AMT TIME CP=DV WT
$DATA THEOPP
$SUBROUTINES ADVAN2 OTHER=myfuncmodule
$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
include "nonmem_reserved_myfunc"
CALLFL=1
KA=THETA(1)+ETA(1)
K=THETA(2)+ETA(2)
CL=THETA(3)*WT+ETA(3)
SC=CL/K/WT
I=mymin(1,2,3,4,5.0)
print *,'I ',I
$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 .4If you use the wrong argument type (real instead of integer), or perhaps use the wrong number of arguments, the compiler will readily flag this.
REFERENCES: Guide Introduction_7
70 YLO Example
In this example, blood concentrations below a known quantification limit (BQL) have been eliminated from the analysis. A one-sided interval with lower boundary equal to the quantification limit is defined. The conditional likelihood for an observation is “conditioned” on values of the population parameters and also on the condition that the observation is inside the one-sided interval. (See pr_y, Section 53.115) (See YLO YUP Probability: PR_Y).
Here is an NM-TRAN control stream.
$PROB YLO example
$INPUT ID TIME DV AMT EVID MDV
$DATA ylow.dat IGNORE=C
$SUBROUTINES ADVAN=1 TRANS=2
$PK
CL=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
K=CL/V
S1=V
$ERROR
; Limit of Quantification (LOQ) is 1
; DV's below LOQ are given MDV=1 in dataset
; PR_Y (PRB) is probability that DV is > LOQ
; PR_Y (PRB) is provided even if MDV is 1
YLO=1
PRB=PR_Y
Y = F+ERR(1)
$THETA
(0,0.693); CL
(0,24); V
$OMEGA
0.09
0.09
$SIGMA
.5
$ESTIMATION
MAXEVAL=9999
NSIG=4
METHOD=COND INTERACTION LAPLACIAN
PRINT=5
$TABLE ID TIME AMT MDV EVID PRB
NOPRINT
ONEHEADER
FILE=ylo.tabHere is a segment of the dataset for a subject with last DV < LOQ
ID TIME DV AMT EVID MDV
. . . . . .
. . . . . .
. . . . . .
2 0 0 100 1 1
2 0.333 2.47 0 0 0
2 2 3.26 0 0 0
2 6 2.32 0 0 0
2 15.917 2.46 0 0 0
2 24 0 100 1 1
2 48 0 100 1 1
2 72 0 100 1 1
2 95.917 4.04 0 0 0
2 96 0 100 1 1
2 96.08 8.85 0 0 0
2 120 3.16 0 0 0
C Last observation record (below) has DV < LOQ (LOQ = 1) and MDV = 1
C This observation does not contribute to the -2LL estimate. However,
C probability that DV is > LOQ can be tabled using PRB=PR_Y.
2 168 0.163 0 0 1REFERENCES: None.