Title: A case example: Building a population pharmacokinetic model for theophylline using NONMEM program.
1A case example Building a population
pharmacokinetic model for theophylline using
NONMEM program.
 Pyry Välitalo
 Orion Pharma
2Population pharmacokineticsPharmacokinetics
using nonlinear mixedeffects models
 Rather than model each individuals response
separately, combine all the data in one model.
The differences between individuals are described
with random effects.  Benefit Possible to use data that would not be
adequate for individual PK modeling (e.g. sparse
sampling).  Benefit More data per model yields more power.
 One model for each individual versus one model
for all data  Benefit Easier to locate sources of
betweensubject variability?
3The NONMEM program
 NONMEM itself is a very general nonlinear
mixedeffects regression program, but it comes
with packages  PREDPP is a package of subroutines for most
common problems.  NMTRAN makes NONMEM input and output more
userfriendly.  Main difference compared to SAS Designed
spesifically to be used in population
pharmacokinetics/pharmacodynamics.  Not as general as SAS.
 Lots of canned solutions to typical
pharmacokinetic problems.
4Some nomenclature
 Fixed effects Theta, ?. Mostly used to describe
the structural model (e.g. clearance) and
covariate relationships (e.g. effect of sex on
clearance)  In general statistical nomenclature b
 Random effects Eta, ?. The standard deviation of
the random effects is ?. In population
pharmacokinetics, random effects are most
commonly used to describe betweensubject
variability.  In general statistical nomenclature bi
 Residual variability Epsilon, e. The standard
deviation of the residual variability is s.  The general statistical nomenclature is similar?
5The dataset (theophylline)
 Part of the data used in
 Upton RA, Thiercelin JF, Guentert TW, Wallace SM,
Powell JR, Sansom L, Riegelman S. Intraindividual
variability in theophylline pharmacokinetics
statistical verification in 39 of 60 healthy
young adults. J Pharmacokinet Biopharm. 1982
Apr10(2)12334.  The dataset comes with NONMEM software.
6The dataset
 12 individuals given a single dose of 36 mg/kg
of oral theophylline.  10 concentration measurements from 25 hours
 A total of 120 observations
 One continuous covariate Weight.
7ADVAN2 subroutine The subroutine to be used
throughout this example.
 ADVAN2 subroutine is used throughout the example.
This subroutine defines an absorption compartment
and a central compartment.  ADVAN2 requires the following to be specified
 KA (absorption rate constant)
 K (elimination rate constant)
 How the observations are scaled (i.e. how do the
drug amounts relate to concentrations)?  How the observation and individual prediction
relate to each other, i.e. what type of residual
error is used?  Everything else is optional.
8Oral pharmacokinetic data and ADVAN2 subroutine
 In its most basic form
 At t0 A(1)Dose and A(2)0
 dA(1)/dtKAA(1) KA
 dA(2)/dtKAA(1)KA(2)

K  Define e.g. that
 KATHETA(1) EXP(ETA(1))
 KTHETA(2) EXP(ETA(2))
 VTHETA(3) EXP(ETA(3)) Volume of distribution
 YA(2)/V EPS(1) observation is the predicted
concentration plus residual error
Depot(1)
Central(2)
9Run1 The zero model
 This control stream came with the dataset, and
likely originates from year 1982  SUBROUTINES ADVAN2
 KATHETA(1)ETA(1) absorption rate constant
 KTHETA(2)ETA(2) elimination rate constant
 CLTHETA(3)WTETA(3) clearance is scaled by
weight  SCCL/K/WT scaling of concentration
observations is done  by calculating VCL/K and scaling by weight
 because dose is weightadjusted.
 THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9) Fixe
d effects estimates  OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4 Random
effects variancecovariance  ERROR
 YFEPS(1) additive error model
 SIGMA .4 Residual error variance estimate
 EST MAXEVAL450 PRINT5 First Order method
10Run1 results
 The parameter estimates make sense.
 The parameter omega(K) near zero.
 High relative standard errors on some random
effects  Omega(Ka) 87
 Covariance(KaCL) 370
 Covariance(KaK) 270
 Omega(K) 51
11Run1 results Plasma concentrations and
predictions.
12Final model Run20 (backup slides include the
steps taken in modelbuilding process)
 SUBROUTINES ADVAN2
 PK
 D1THETA(6)EXP(ETA(3))
 KATHETA(1)EXP(ETA(1))
 VTHETA(2)EXP(ETA(2)THETA(5)) WT/70
 CLTHETA(3)EXP(ETA(2)) WT/70
 KCL/V
 S2V
 ERROR
 IPREDF
 IRESDVIPRED
 WTHETA(4) add error absorption phase
 W2THETA(7) add error postabs phase
 IWRESIRES/W
 IF(TIME.GT.2) IWRESIRES/W2
 YIPREDWEPS(1)
 IF(TIME.GT.2) YIPREDW2EPS(1)
13Aspects of the improved modelResidual error
model
 SUBROUTINES ADVAN2
 PK
 D1THETA(6)EXP(ETA(3))
 KATHETA(1)EXP(ETA(1))
 VTHETA(2)EXP(ETA(2)THETA(5)) WT/70
 CLTHETA(3)EXP(ETA(2)) WT/70
 KCL/V
 S2V
 ERROR
 IPREDF
 IRESDVIPRED
 WTHETA(4) additive error defined
with fixed effect. THETA(4) equals sd of res.var.  W2THETA(7) additive error defined with fixed
effect. THETA(7) equals sd of res.var.  IWRESIRES/W IWRES are weighted with standard
deviation of residual variability  IF(TIME.GT.2) IWRESIRES/W2 for
postabsorption period  YIPREDWEPS(1) multiply by W because sd of
residual variability is fixed to 1.  IF(TIME.GT.2) YIPREDW2EPS(1) for
postabsorption period
14Benefits of coding residual error model this way
 We can obtain IWRES values as per definition
 IWRES (observation IPRED)/s (e.g. Karlsson
Savic 2005)  This is used to standardize the IWRES, so that
they should have a standard deviation of 1 (and
mean of 0).  Sometimes, using a fixed effect makes the model
more stable than estimating a lot of random
effects.  Regarding the use of separate residual errors for
absorption and postabsorption phases This makes
sense because the absorption phase typically has
more noise than postabsorption phase (and IV
data is yet more reliable).  Has been proposed in some articles, e.g.
 Karlsson MO, Beal SL, Sheiner LB. J Pharmacokinet
Biopharm. 1995 Dec23(6)65172.  Chan PL, Weatherley B, McFadyen L. Br J Clin
Pharmacol. 2008 Apr65 Suppl 17685.
15Aspects of the improved modelDistribution of
random effects, estimation method, implementing
weight as a covariate
 SUBROUTINES ADVAN2 This is a predefined
subroutine for 1compartment model  with first order absorption
 PK
 D1THETA(6)EXP(ETA(3)) lognormal
distribution  KATHETA(1)EXP(ETA(1)) lognormal
distribution  VTHETA(2)EXP(ETA(2)THETA(5))
WT/70 lognormal distribution, linear scaling
by weight  CLTHETA(3)EXP(ETA(2)) WT/70 lognormal
distribution, linear scaling by weight  KCL/V
 S2V
 ERROR
 IPREDF
 IRESDVIPRED
 WTHETA(4) add error
 W2THETA(7) add err postabs
 IWRESIRES/W
 IF(TIME.GT.2) IWRESIRES/W2
 YIPREDWEPS(1)
 IF(TIME.GT.2) YIPREDW2EPS(1)
16The practical difference between First Order and
First Order Conditional Estimation
 First Order (FO) method expands around ETA0
 Thus, the parameters are estimated for the mean
response and the estimation could be described
POPULATION AVERAGE.  First Order Conditional Estimation (FOCE) expands
around the expected values of each random effect.
 Thus, the parameters are estimated for the
median subject and the estimation could be
described SUBJECT SPESIFIC.  (On top of that, FOCEI takes into account the
interaction between random effects and residual
error)
17Implementing bodyweight in the model
 Linear scaling of both clearance and volume of
distribution by bodyweight.  Current trend, at least for pediatrics linear
scaling for volume of distribution, nonlinear
scaling for clearance (estimate an exponent for
nonlinear scaling).  Probably not applicable to obese patients, etc.
 Weight range in this data 5586kg (median 70kg).
 Probably not wide enough range to estimate
nonlinearity.
18Aspects of the improved modelAbsorption model
 SUBROUTINES ADVAN2
 PK
 D1THETA(6)EXP(ETA(3)) D1duration of a
hypothetical infusion into the depot
compartment  KATHETA(1)EXP(ETA(1))
 VTHETA(2)EXP(ETA(2)THETA(5)) WT/70
 CLTHETA(3)EXP(ETA(2)) WT/70
 KCL/V
 S2V
 ERROR
 IPREDF
 IRESDVIPRED
 WTHETA(4) add error
 W2THETA(7) add err postabs
 IWRESIRES/W
 IF(TIME.GT.2) IWRESIRES/W2
 YIPREDWEPS(1)
 IF(TIME.GT.2) YIPREDW2EPS(1)
19The absorption model Mixed 0 and 1storder
absorption.
 If (timeltD1) dA(1)/dtDose/D1 KAA(1)
 If (timegtD1) dA(1)/dtKAA(1)
 A lagtime model was also tried but the mixed 0
and 1storder absorption proved better with the
following benefits  Describes the data better (lower OFV)
 Would make it possible to use proportional error
model?  Not as radical a change point as lagtime
possibly more stable than using lagtime?  Probably better still Transit compartmental
absorption (was not tried in this case because of
time limitations).  Savic RM, Jonker DM, Kerbusch T, Karlsson MO. J
Pharmacokinet Pharmacodyn. 2007 Oct34(5)71126.
Epub 2007 Jul 26.
20What was the additional effort worth?
 Compared to the original run we have
 More complicated absorption model.
 More complicated residual error model .
 (No random effect for elimination rate constant).
 What did we benefit from this?
 A better absorption model may also result in more
accurate estimates of V/F and CL/F.  Since different magnitudes of residual error are
identified for absorption and postabsorption
phases, the model has better simulation
properties.  The model overall describes the data better
(lower residual error).
21Original model parameter estimates versus final
model parameter estimates
Parameter Run3 Run20
OFV 111 (6 parameters) 1.84 (10 parameters)
KA (1/h) 1.6 (0.20) 2.8 (0.26)
VCL/K (l) 32 (0.045) 34 (0.038)
CL (l/h) 2.7 (0.079) 2.7 (0.074)
Omega(KA) 0.65 (0.53) 0.90 (0.44)
Omega(CL) 0.17 (0.30) 0.27 (0.56)
Sigma (mg/l) 0.71 (0.13) 0.66 (0.15) / 0.30 (0.11)
the numbers in parenthesis are relative standard
errors
22Takehome messages
 If you have problems making the model converge
 Realize that noncontinuous functions may be hard
to integrate (e.g. lagtime for absorption)  Use the appropriate residual error model.
 (In NONMEM start with a simple model and build
from ground up)
23Resources
 NONMEM. Currently the golden standard in
population pharmacokinetic modeling. Requires
license.  http//www.icondevsolutions.com/nonmem.htm
 Xpose An R package that helps in deciphering the
output of NONMEM. Free.  http//xpose.sourceforge.net/
 R A program needed by Xpose to operate. Free.
 http//www.rproject.org/
 Census. A helpful program for keeping record of
NONMEM runs. Free.  http//census.sourceforge.net/
 PsN (PerlspeaksNonmem). A collection of helpful
Perl scripts for NONMEM that make life easier in
a lot of ways. Free.  MONOLIX. Another population pharmacokinetic
program. Free.  Advantages Shorter runtimes than in NONMEM,
provides graphical output by itself.  Disadvantages Currently not as flexible as
NONMEM. May be hard to make sense of the error
messages encountered.  http//software.monolix.org/sdoms/software/
24(Backup slides)Run2 Modernize the original
run
 Dataset modifications Input the exact dose that
was given to each patient rather than a
weightscaled dose. Input weight as a covariate
for every observation.  Rather than model additive random effects, model
exponential random effects THETA(1)ETA(1)
becomes THETA(1)EXP(ETA(1))  Because lognormal distributions very typical in
pharmacokinetics
25(Backup slides)Run2 Modernize the original
run
 Take away covariance matrix for now, but leave
all the random effects parameters.  In NONMEM, unnecessary covariance structures can
slow down the estimation and add to model
instability.  Instead of First Order (FO) estimation method,
use First Order Conditional Estimation with
Interaction (FOCEI) method In ESTIMATION
block, METHOD1 INTER line is added  First Order method would expand around ETA0
(Fixed effects estimated to describe Population
Average)  First Order Conditional Estimation with
Interaction expands around ?E(?). Thus, the
fixed effects describe a typical subject rather
than population average (Subject Spesific
approximation).
26(Backup slides)Run2 results
 OFV 111.985
 Parameter estimate is near its lower boundary
omega(K)  Left individual timeconcentration plots.
Right Individual random effects KA versus CL.
No linear correlation seems to be present thus,
KACL covariance is not implemented.
27(Backup slides)Forward some runs
 Run3 Remove random effect on K
 OFV 111.983, large standard error found for
omega(KA)  Run4 Based on run3. Remove random effect on KA
 OFV 186.687. Even though there is difficulty
estimating the random effect, it seems to be an
important part of the model.  Run5 Based on run3. Implement additiveproportion
al error.  OFV 105.320 , large standard error found for
omega(KA)  Run6 Based on run3. Implement proportional
residual error.  OFV 140.536, fairly large standard error found
for omega(KA)  SUMMARY of this slide The random effect of K
(the elimination rate) is not essential to the
model. The additiveproportional residual error
might improve the model a bit. However, more
efforts will be spent on the residual error model
later on thats why additiveproportional
residual error is discontinued.
28(Backup slides)Run7 Based on run3. Estimate
KA, V and CL instead of KA, K and CL
 BEFORE (run3)
 KATHETA(1)EXP(ETA(1))
 KTHETA(2)
 CLTHETA(3)EXP(ETA(2))
 VCL/K
 AFTER MODIFICATION (run7)
 KATHETA(1)EXP(ETA(1))
 VTHETA(2)EXP(ETA(2))
 CLTHETA(3)EXP(ETA(2))
 KCL/V
 In run3, it was found that a random effect is not
needed for K. Since KCL/V, in run7 the same
random effect is included for both V and CL.
Results identical to run3 (OFV 111.983)
29(Backup slides)Runs 89 Still testing the
random effects
 Run8 Take away random effect on V
 OFV 130.861. Standard errors for parameter
estimates increase somewhat.  Run9 Based on run7. Include a separate random
effect for V and implement covariance for CL and
V.  dOFV 5.4
 According to the results, the correlation between
the two random effects would be absolute,
although the magnitude of these random effects
differs. gt Something fishy about  the dataset?
 r covariance (ETA(CL)ETA(V) )
 (variance(ETA(CL))variance(ETA(V))
30(Backup slides)The covariance between var(CL)
and var(V)?
 Run10 Reparameterize run9 Use the same random
effect for both CL and V but scale with a fixed
effect.  OFV 106.644. Results practically identical to
run9.  The code in run9
 VTHETA(2)EXP(ETA(3)) separate random
effect for each  CLTHETA(3)EXP(ETA(2)) parameter.
 OMEGA BLOCK(2) .2 .2 .2 var(CL), covar(CL,V),
var(V)  And the code in run10
 VTHETA(2)EXP(ETA(2)THETA(5)) scale by
theta(5)  CLTHETA(3)EXP(ETA(2))
 OMEGA .2 var(CL), used also for V.
31(Backup slides)Theophylline absorption modeling
 Run11 Based on run10 Lagtime on
 absorption? (1 parameter)
 OFV 64.447. Improvement on some individual
 plots.
 Run12 Based on run11 Lagtime for absorption
 with random effect? (1 parameter)
 OFV 47.179. Some difficulties to make this run
 work. A likely reason for difficulties is that
 it is hard to integrate at timelagtime, as
absorption starts instantaneously.  How to implement lagtime in a model? Write
 ALAG1THETA(n) this would estimate lagtime for
compartment 1.
32(Backup slides)Theophylline absorption modeling
 Run13 Mixed 0 and 1storder absorption
 OFV 64.879
 Run14 Based on run13. Mixed 0 and 1storder
absorption a random effect is included for
0order absorption rate.  31.140. The standard error of var(KA) becomes
lower. This seems the most reasonable solution.  How to implement mixed 0 and 1storder
absorption? Add a column RATE into the datafile
and set it to 2 for every instance when AMTgt0
(else .). Write in the modelfile  D1THETA(n) this would estimate the duration of
a hypothetical  infusion into the absorption compartment
33(Backup slides)Absorption phase and residual
errors?
 Oral dosing data is usually considered more
noisy than intravenous dosing data. One reason
for this is that despite the best absorption
modeling efforts, drug absorption from GI tract
can be erratic and/or unpredictable.  Thus, more variability can usually be observed in
the absorption phase than in postabsorption
phase.  There are instances when a timedependent
residual error has been assigned to differentiate
between absorption phase and postabsorption
phase, e.g.  IF(T.LE.2) YIPREDEPS(1) absorption phase,
additive error  IF(T.GT.2) YIPREDEPS(2) postabs phase,
additive error  For example
 Karlsson MO, Beal SL, Sheiner LB. J Pharmacokinet
Biopharm. 1995 Dec23(6)65172.  Chan PL, Weatherley B, McFadyen L. Br J Clin
Pharmacol. 2008 Apr65 Suppl 17685.
34(Backup slides) Runs 1516 Test using
different residual errors for absorption and
postabsorption phases
 Run15 Based on run10. Try setting the
absorptionphase residual error to before 2 hours
and postabsorption phase residual error to after
2 hours.  Success OFV 49.023, the residual error of
postabsorption phase drops to onehalf of what
it was previously.  Run16 Based on run15. Try using proportional
residual error for postabsorption phase.  OFV 67.746. Does not improve the model.
35(Backup slides) Combine the timedependent
residual error model with absorption model.
Implement weight as a covariate.
 Run17 Based on run14. Combine mixed 0 and
1storder absorption and timedependent residual
error  OFV 19.804
 Run18 Based on run17. Check if a random effect
is still needed in the 0order absorption rate.  OFV 31.530 and the standard error of omega(KA)
increases.  Runs 1922 Try different ways of scaling the
volume of distribution and clearance by weight.  Run20 with linear scaling seems to work best
 Conclusion Run20 shall be the final model.