Transcript Orion Corp

A case example: Building a population
pharmacokinetic model for theophylline
using NONMEM program.
Pyry Välitalo
Orion Pharma
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
1
Population pharmacokinetics:
Pharmacokinetics using nonlinear mixed-effects 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 between-subject
variability?
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
2
The NONMEM program
• NONMEM itself is a very general nonlinear mixed-effects
regression program, but it comes with packages:
– PREDPP is a package of subroutines for most common
problems.
– NM-TRAN 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.
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
3
Some 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 between-subject
variability.
– In general statistical nomenclature: bi
• Residual variability: Epsilon, ε. The standard deviation of the
residual variability is σ.
– The general statistical nomenclature is similar?
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
4
The 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
Apr;10(2):123-34.
• The dataset comes with NONMEM software.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
5
The dataset
• 12 individuals given a single dose of 3-6 mg/kg of oral
theophylline.
• 10 concentration measurements from 25 hours
– A total of 120 observations
• One continuous covariate: Weight.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
6
ADVAN2 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.
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
7
Oral pharmacokinetic data and ADVAN2
subroutine
• In its most basic form:
At t=0 A(1)=Dose and A(2)=0
dA(1)/dt=-KA*A(1)
dA(2)/dt=KA*A(1)-K*A(2)
Depot(1)
KA
Central(2)
K
• Define e.g. that
–
–
–
–
KA=THETA(1) * EXP(ETA(1))
K=THETA(2) * EXP(ETA(2))
V=THETA(3) * EXP(ETA(3))
;Volume of distribution
Y=A(2)/V + EPS(1)
;observation is the predicted
;concentration plus residual error
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
8
Run1: ”The zero model”
• This control stream came with the dataset, and likely
originates from year 1982:
$SUBROUTINES ADVAN2
KA=THETA(1)+ETA(1)
K=THETA(2)+ETA(2)
CL=THETA(3)*WT+ETA(3)
SC=CL/K/WT
;absorption rate constant
;elimination rate constant
;clearance is scaled by weight
;scaling of concentration observations is done
;by calculating V=CL/K and scaling by weight
;because dose is weight-adjusted.
$THETA (.1,3,5) (.008,.08,.5) (.004,.04,.9)
;Fixed effects estimates
$OMEGA BLOCK(3) 6 .005 .0002 .3 .006 .4
;Random effects variance-covariance
$ERROR
Y=F+EPS(1)
;additive error model
$SIGMA
;Residual error variance estimate
$EST
.4
MAXEVAL=450
PRINT=5 ;First Order method
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
9
Run1 results:
• The parameter estimates make sense.
• The parameter omega(K) near zero.
• High relative standard errors on some random effects:
–
–
–
–
Omega(Ka): 87%
Covariance(Ka-CL): 370%
Covariance(Ka-K): 270%
Omega(K): 51%
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
10
Run1 results: Plasma concentrations and
predictions.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
11
Final model: Run20 (back-up slides include
the steps taken in model-building process)
$SUBROUTINES
ADVAN2
$PK
D1=THETA(6)*EXP(ETA(3))
KA=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2)*THETA(5)) *WT/70
CL=THETA(3)*EXP(ETA(2)) *WT/70
K=CL/V
S2=V
$ERROR
IPRED=F
IRES=DV-IPRED
W=THETA(4)
;add error absorption phase
W2=THETA(7)
;add error postabs phase
IWRES=IRES/W
IF(TIME.GT.2) IWRES=IRES/W2
Y=IPRED+W*EPS(1)
IF(TIME.GT.2) Y=IPRED+W2*EPS(1)
$THETA (.1,2) (0,33) (0,3) (0,1) (0,.5) (0,.3) (0,.3)
$OMEGA .2 .2 .2
$SIGMA 1 FIX
$EST METHOD=1 INTER
$COV
MAXEVAL=9999
PRINT=5
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
12
Aspects of the improved model:
Residual error model
$SUBROUTINES
ADVAN2
$PK
D1=THETA(6)*EXP(ETA(3))
KA=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2)*THETA(5)) *WT/70
CL=THETA(3)*EXP(ETA(2)) *WT/70
K=CL/V
S2=V
$ERROR
IPRED=F
IRES=DV-IPRED
W=THETA(4)
;additive error defined with fixed effect. THETA(4) equals sd of res.var.
W2=THETA(7)
;additive error defined with fixed effect. THETA(7) equals sd of res.var.
IWRES=IRES/W
;IWRES are weighted with standard deviation of residual variability
IF(TIME.GT.2) IWRES=IRES/W2
;for post-absorption period
Y=IPRED+W*EPS(1)
;multiply by W because sd of residual variability is fixed to 1.
IF(TIME.GT.2) Y=IPRED+W2*EPS(1)
;for post-absorption period
$THETA (.1,2) (0,33) (0,3) (0,1) (0,.5) (0,.3) (0,.3)
$OMEGA .2 .2 .2
$SIGMA 1 FIX
;Variance of EPS(1) values is fixed to 1.
$EST METHOD=1 INTER
$COV
MAXEVAL=9999
PRINT=5
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
13
Benefits of coding residual error model this
way:
• We can obtain IWRES values as per definition:
IWRES = (observation – IPRED)/σ (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 post-absorption 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 Dec;23(6):651-72.
Chan PL, Weatherley B, McFadyen L. Br J Clin Pharmacol. 2008 Apr;65 Suppl 1:76-85.
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
14
Aspects of the improved model:
Distribution of random effects, estimation method, implementing weight
as a covariate
$SUBROUTINES
ADVAN2
;This is a predefined subroutine for 1-compartment model
;with first order absorption
$PK
D1=THETA(6)*EXP(ETA(3))
KA=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2)*THETA(5)) *WT/70
CL=THETA(3)*EXP(ETA(2)) *WT/70
K=CL/V
S2=V
$ERROR
IPRED=F
IRES=DV-IPRED
W=THETA(4)
;add error
W2=THETA(7)
;add err postabs
IWRES=IRES/W
IF(TIME.GT.2) IWRES=IRES/W2
Y=IPRED+W*EPS(1)
IF(TIME.GT.2) Y=IPRED+W2*EPS(1)
;log-normal
;log-normal
;log-normal
;log-normal
distribution
distribution
distribution, linear scaling by weight
distribution, linear scaling by weight
$THETA (.1,2) (0,33) (0,3) (0,1) (0,.5) (0,.3) (0,.3)
$OMEGA .2 .2 .2
$SIGMA 1 FIX
$EST METHOD=1 INTER
$COV
MAXEVAL=9999
PRINT=5
Orion Corporation
;method is First Order Conditional Estimation with
;Interaction
Pyry Välitalo SSL Presentation
1.1.2007
15
The practical difference between First Order
and First Order Conditional Estimation
• First Order (FO) method expands around ETA=0
– 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, FOCE-I takes into account the interaction
between random effects and residual error)
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
16
Implementing 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 55-86kg (median 70kg).
– Probably not wide enough range to estimate nonlinearity.
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
17
Aspects of the improved model:
Absorption model
$SUBROUTINES
ADVAN2
$PK
D1=THETA(6)*EXP(ETA(3))
;D1=duration of a hypothetical infusion into the depot
;compartment
KA=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2)*THETA(5)) *WT/70
CL=THETA(3)*EXP(ETA(2)) *WT/70
K=CL/V
S2=V
$ERROR
IPRED=F
IRES=DV-IPRED
W=THETA(4)
; add error
W2=THETA(7)
;add err postabs
IWRES=IRES/W
IF(TIME.GT.2) IWRES=IRES/W2
Y=IPRED+W*EPS(1)
IF(TIME.GT.2) Y=IPRED+W2*EPS(1)
$THETA (.1,2) (0,33) (0,3) (0,1) (0,.5) (0,.3) (0,.3)
$OMEGA .2 .2 .2
$SIGMA 1 FIX
$EST METHOD=1 INTER
$COV
MAXEVAL=9999
PRINT=5
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
18
The absorption model: Mixed 0- and 1storder absorption.
If (time<D1) dA(1)/dt=Dose/D1 –KA*A(1)
If (time>D1) dA(1)/dt=–KA*A(1)
• A lagtime model was also tried but the mixed 0- and 1st-order
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 Oct;34(5):711-26. Epub 2007 Jul 26.
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
19
What 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 post-absorption phases, the model has better
simulation properties.
– The model overall describes the data better (lower residual
error).
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
20
Original 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)
V=CL/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
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
21
Take-home messages
• If you have problems making the model converge:
– Realize that non-continuous 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”)
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
22
Resources
•
•
•
•
•
•
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.r-project.org/
Census. A helpful program for keeping record of NONMEM runs. Free.
– http://census.sourceforge.net/
PsN (Perl-speaks-Nonmem). 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/
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
23
(Back-up slides)
Run2: ”Modernize” the original run
• Dataset modifications: Input the exact dose that was given to
each patient rather than a weight-scaled 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 log-normal distributions very typical in
pharmacokinetics
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
24
(Back-up 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 (FOCE-I) method: In
$ESTIMATION block, ”METHOD=1 INTER” line is added
– First Order method would expand around ETA=0 (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).
Orion Corporation
Pyry Välitalo SSL Presentation
1.1.2007
25
(Back-up slides)
Run2 results:
• OFV 111.985
• Parameter estimate is near its lower boundary: omega(K)
• Left: individual time-concentration plots.
Right: Individual
random effects
KA versus CL. No
linear correlation
seems to be
present; thus, KACL covariance is
not implemented.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
26
(Back-up 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 additive+proportional 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 additive+proportional residual error might
improve the model a bit. However, more efforts will be spent on the residual
error model later on; that’s why additive+proportional residual error is
discontinued.
Orion Corporation
Author/Department
2.2.2010
27
(Back-up slides)
Run7: Based on run3. Estimate KA, V and
CL instead of KA, K and CL:
• BEFORE (run3):
KA=THETA(1)*EXP(ETA(1))
K=THETA(2)
CL=THETA(3)*EXP(ETA(2))
V=CL/K
• AFTER MODIFICATION (run7):
KA=THETA(1)*EXP(ETA(1))
V=THETA(2)*EXP(ETA(2))
CL=THETA(3)*EXP(ETA(2))
K=CL/V
• In run3, it was found that a random effect is not needed for K.
Since K=CL/V, in run7 the same random effect is included for
both V and CL. Results identical to run3 (OFV 111.983)
Orion Corporation
Author/Department
2.2.2010
28
(Back-up slides)
Runs 8-9: 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. -> Something fishy about
the dataset?
– r= covariance (ETA(CL)-ETA(V) )
(variance(ETA(CL))*variance(ETA(V))
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
29
(Back-up 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:
V=THETA(2)*EXP(ETA(3))
CL=THETA(3)*EXP(ETA(2))
$OMEGA BLOCK(2) .2 .2 .2
;separate random effect for each
;parameter.
;var(CL), covar(CL,V), var(V)
• And the code in run10:
V=THETA(2)*EXP(ETA(2)*THETA(5))
;scale by theta(5)
CL=THETA(3)*EXP(ETA(2))
$OMEGA .2
;var(CL), used also for V.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
30
(Back-up 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 time=lagtime, as absorption starts
instantaneously.
• How to implement lagtime in a model? Write
ALAG1=THETA(n)
;this would estimate lagtime for compartment 1.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
31
(Back-up slides)
Theophylline absorption modeling
• Run13: Mixed 0- and 1st-order absorption
– OFV 64.879
• Run14: Based on run13. Mixed 0- and 1st-order absorption; a
random effect is included for 0-order absorption rate.
– 31.140. The standard error of var(KA) becomes lower. This
seems the most reasonable solution.
• How to implement mixed 0- and 1st-order absorption? Add a
column ”RATE” into the datafile and set it to -2 for every
instance when AMT>0 (else ”.”). Write in the modelfile:
D1=THETA(n) ;this would estimate the duration of a hypothetical
;infusion into the absorption compartment
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
32
(Back-up 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 post-absorption phase.
• There are instances* when a time-dependent residual error
has been assigned to differentiate between absorption phase
and post-absorption phase, e.g.
IF(T.LE.2) Y=IPRED+EPS(1)
IF(T.GT.2) Y=IPRED+EPS(2)
;absorption phase, additive error
;post-abs phase, additive error
*For example:
•
Karlsson MO, Beal SL, Sheiner LB. J Pharmacokinet Biopharm. 1995 Dec;23(6):651-72.
•
Chan PL, Weatherley B, McFadyen L. Br J Clin Pharmacol. 2008 Apr;65 Suppl 1:76-85.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
33
(Back-up slides): Runs 15-16: Test using
different residual errors for absorption and
post-absorption phases
• Run15: Based on run10. Try setting the absorption-phase
residual error to before 2 hours and post-absorption phase
residual error to after 2 hours.
– Success: OFV 49.023, the residual error of post-absorption
phase drops to one-half of what it was previously.
• Run16: Based on run15. Try using proportional residual error
for post-absorption phase.
– OFV 67.746. Does not improve the model.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
34
(Back-up slides): Combine the time-dependent
residual error model with absorption model.
Implement weight as a covariate.
• Run17: Based on run14. Combine mixed 0- and 1st-order
absorption and time-dependent residual error
– OFV 19.804
• Run18: Based on run17. Check if a random effect is still
needed in the 0-order absorption rate.
– OFV 31.530 and the standard error of omega(KA) increases.
• Runs 19-22: 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.
Orion Corporation
Pyry Välitalo SSL Presentation
2.2.2010
35