Transcript An Introduction of PK/PD Modeling - Part 1

Introduction to PK/PD Modeling for
Statisticians
Part 1
Alan Hartford, AbbVie
ASA Biopharm FDA-Industry Statistics Workshop
September 16, 2015
Objectives and Outline for Part 1
• Provide statisticians with the main concepts of PK/PD
modeling (a.k.a., Pharmacometrics) (see outline below)
• Encourage statisticians to support pharmacometricians in
their modeling efforts
• Provide an appreciation for Pharmacometrics
• We won’t be reviewing step by step modeling methods using
AIC or p-values
• Outline: PK Basic Concepts, PK Compartmental Models,
Technical Considerations, Software Considerations, Adding
covariates to PK models
2
Introduction – PK and PD
Pharmacokinetics is the study of what an organism does with a dose
of a drug
– kinetics = motion
– Absorbs, Distributes, Metabolizes, Excretes
PK Endpoints
– AUC, Cmax, Tmax, half-life (terminal), Cmin
The effect of the drug is assumed to be related to some measure of
exposure. (AUC, Cmax, Cmin)
Pharmacodynamics is the study of what the drug does to the body
(dynamics = change)
3
Procedure:
PK/PD Modeling
– Fit a model to exposure data and estimate exposure at the
same time points where we have PD data.
– Examine correlation between estimated exposure and PD (or
other endpoints, e.g., AE rates).
– Might need to fit a mechanistic model with exposure data as
explanatory variable and PD as response.
Purpose:
–
–
–
–
Estimate therapeutic window
Dose selection
Identify mechanism of action
Model probability of AE as function of exposure (and
covariates)
4
Concentration of Drug as a Function of Time
Model for Extra-vascular Absorption
Concentration
Cmax
AUC
Tmax
Time
5
Observed or Predicted PK?
• Exposure endpoints are not measured – only
modeled, i.e., estimated
• Concentration in blood or plasma is a biomarker for
concentration at site of action
• PK parameters are not directly measured
6
Compartmental Modeling
• A person’s body is modeled with a system of
differential equations, one equation for each
“compartment”
• If each equation represents a specific organ or set of
organs with similar perfusion rates, then called
Physiologically Based PK (PBPK) modeling.
• The mean function f is a solution of this system of
differential equations.
• Each equation in the system describes the flow of drug
into and out of a specific compartment.
7
Helpful Reminder …
𝑑 𝑡
𝑒 = 𝑒𝑡
𝑑𝑡
𝑑 𝑘𝑡
𝑒 = 𝑘𝑒 𝑘𝑡
𝑑𝑡
𝑑 −𝑡
𝑒 = −𝑒 −𝑡
𝑑𝑡
𝑑 −𝑘𝑡
𝑒
= −𝑘𝑒 −𝑘𝑡
𝑑𝑡
8
First-Order 1-Compartment
Model (Intravenous injection)
Input
dAc
 k10 Ac
dt
Cc  Ac / Vc
Central
Vc
Ac t  0  Dose
Elimination
k10
Solution:
Dose  k1 0t
Cc (t ) 
e
Vc
9
First-Order 1-Compartment
Model (Intravenous injection)
Parameterized with Clearance
Input
Central
Vc
Elimination
Another parameterization for the solution
uses Clearance = Cl = k10 Vc
Clearance = Volume of drug eliminated
per unit time
k10
Solution:
Dose t Cl / V
Cc (t ) 
e
Vc
10
First-Order 1-Compartment
Model (Extravascular Administration)
Input
ka
Central
Absorption depot:
Central compartment:
Vc
Aa t  0   dose
Elimination
Ac t  0   0
k10
Solution:
dAa
  k a Aa
dt
dAc
 k a Aa  k10 Ac
dt
Cc  Ac / Vc

k a  F  Dose  k1 0t
Cc (t ) 
e
 e kat
Vc k a  k10 

F = Bioavailability
(i.e., fraction absorbed)
11
First-Order 1-Compartment
Model (Extravascular Administration)
Parameterized with Clearance
Input
Cl  Vc  k10
ka
Central
Vc
Solution:

k a  F  Dose Cl / Vc t
Cc (t ) 
e
 e kat
Vc k a  Cl

Elimination
k10
F = Bioavailability
(i.e., amount absorbed)
12
Parameterization
• ka, k10, V
– Micro constant
• ka, Cl, V
– Macro constant
• Note that usually F, V, and Cl are not estimable
(unless you perform studies with both IV and
extravascular administration)
• Instead, apparent V (V/F) and apparent Cl (Cl/F) are
estimated when only extravascular data are available
13
First-Order 2-Compartment
Model (Intravenous Dose)
Input
k12
Peripheral
Central
(Vp)
Vc
k21
Elimination
k10
General form of solution:
dAc
 k 21 Ap  k12  k10 Ac
dt
dAp
 k12 Ac  k 21 Ap
dt
Cc  Ac / Vc
C p  Ap / V p
Ac t  0   Bolus Dose
Ap t  0   0
Cc t   A exp( t )  B exp( t )
14
First-Order 2-Compartment
Model (Intravenous Dose)
Parameterized in terms of
“Micro constants”
Note that including Vp overparameterizes the model since
Input
k12
Peripheral
Central
(Vp)
Vc
k21
k12
Vp 
Vc
k 21
Elimination
k10
Ac = Amount of drug in central compartment
Ap = Amount of drug in peripheral compartment
15
Web Demonstration
• http://vam.anest.ufl.edu/simulations/simulati
onportfolio.php
• http://vam.anest.ufl.edu/simulations/stochast
ictwocompartment.php
• (Requires installation of Adobe Shockwave
player.)
16
Technical Considerations
Outline
• General form of NLME
– Parameterization
– Error Models
• Model fitting
– (Approximate) Maximum Likelihood
– Fitting Algorithms
17
The Nonlinear Mixed Effects Model
yij is the jth response for the i th subject
f is a scalar function nonlinear in 
yij  f (tij ,  i , d i )   ij
 i ~ N  , D 
 i ~ N 0, Ri 
 is a k 1 parameter vector
tij is the jth time for the i th subject
d i is the i th subject' s dose
j ranges from 1 to ni
 ij is residual error
D is a k  k covariance matrix
Ri is an ni  ni covariance matrix
Pharmacokineticists use the term ”population”
model when the model involves random effects.
18
• For simplification at this stage, assume
 i    bi
b i ~ N 0, D 
and
 i ~ N 0,  I i 
2
19
Assay Variability
• Assays for measuring PK concentrations are validated
for specific concentration ranges.
• If the concentration is higher than the upper limit of
the validated assay, then the sample is diluted so the
resulting diluted sample has PK concentration within
the validated limits.
• If the concentration is lower than the lower limit of
the validated assay, then the concentration is
reported to be “below the limit of quantitation”
(“BLOQ” or “<LLOQ”).
20
Assay Variability (cont.)
• The result of the assay specifications and the needed
dilutions is that additional error is added into the
measurement system.
• These errors can be accounted for in the statistical
model.
21
Distribution of Error
• In each case, the errors are assumed to be normally
distributed with mean 0
• In PK literature, the variance is assumed to be
constant (2)
• Heteroscedastic variance is modeled using a
proportional error term
• Another option is to use the additive error model
assuming a variance function R(q) where q is an m x
1 vector which can incorporate , D and other
parameters, e.g. R(q)=2[f()]2, q[, T]T
22
Error Models used for PK modeling
f tij ,  i    ij
f tij ,  i 1  eij 
f tij ,  i 1  eij    ij
f tij ,  i e
 ij
Additive error
Proportional error
Additive and
Proportional error
Exponential error
23
For the 1-compartment model parameterized
with Cl, V, ka
Input
ka
Central
Vc
Elimination
logk a 
Cl  k10 V ,   logCl 
 logV 
k10
And cov(logCli, logVi) is assumed to be 0 by
definition of the pharmacokinetic parameters.
24
Maximum Likelihood Approach Is Used
We obtain the maximum likelihood estimate by
maximizing
N
 p y 
i
i 1
Where p(yi) is the probability distribution function
(pdf) of y where now we use the notation of yi as a
vector of all responses for the ith subject
The problem is that we don’t have this probability
density function for y directly.
25
• We use the following:
 p  y     p y
N
N
i
i 1
i 1



|
b

b
d
b
i
i
i
i
Where p and  are normal probability density
functions. Maximization is in
[T, vech(D), vech(R)] T
Notation: the vech function of a matrix is
equal to a vector of the unique elements of the
matrix.
26
Under Normal Assumptions






 1

T
1
p y i | b i  B1  exp  y i  f i Ri y i  f i 
 2

 f ti1 ,   b i , d i 
 f t , b ,d 
i2
i 
i

fi 





 f tini ,   b i , d i 






 1 T 1 
 b i   B2 exp  b i D b i 
 2

B1 and B2  constant w.r.t. b i
27
Maximum Likelihood
Given data yij, we use maximum likelihood to obtain
parameters estimates for , D, and 2.
Because the mean function, f, is assumed to be
nonlinear in i in pharmacokinetics, least squares
does not result in equivalent parameter estimates.
28
Approximate Methods
• Use numerical approaches to approximate the
integral and then maximize the approximation
• Options:
– Approximate the integrand by something we can
integrate
• First Order method (Taylor series)
– Approximate the whole integral
• Laplace’s approximation (second order approximation)
• Importance Sampling
• Gaussian Quadrature
– Use Bayesian methodology
29
Algorithms Used
First Order
Available in NONMEM
Approximate just the integrand
First Order Conditional
Estimation
Laplace’s Approximation
Importance Sampling
  p y
M
i 1
Gaussian Quadrature



|
b

b
d
b
i
i
i
i
Or approximate whole integral
Bayesian (Gibb’s Sampler; Not
covered in this presentation)
(NONMEM is the gold standard software package for PKPD modeling.)
30
First Order Method
Approximate with a first order Taylor series expansion
If the model assumes
f tij ,  i , d i   f tij ,   bi , d i 
And Ri = 2I, then this is pretty straight-forward.
You use a Taylor series expansion about bi.
31
Taylor Series Expansion
With a first order Taylor series approximation
expanded about , the mean of the i
f tij ,   bi , di  
f tij ,  , di    / x f tij , x, di  |x   bi
Let this approximation be
f Tay tij ,  , bi , di 
You use this approximation in the integrand.
32
Substituting back in and simplifying …






 1

Tay T
p y i | b i  B1  exp  y i  f i Ri1 y i  f Tay i 
 2

 f Tay ti1 ,  , b i , d i 
 Tay



f
t
,

,
b
,
d
i2
i
i 
f Tay i  



 Tay

 f tini ,  , b i , d i 
 1 T 1 
 b i   B2 exp  b i D b i 
See slide 26.
 2

B1 and B2  constant w .r.t. b i


And now the exponent term is integrable and now we can
maximize the likelihood.
33
Using Laplace’s Approximation
A second order approximation can be constructed by
using Laplace’s approximation
 2
 exp[ni  i bi ]dbi   ni
where
bˆ maximizes  b 
i
i




k /2
 i bˆi
1 / 2
e
 
ni  i bˆi
i
In this manner, the whole integral is approximated so no
integration is needed.
1
 i bi   log p y i | b i  b i 
ni
34
[

]
Numerical Integration: Importance Sampling
Consider a function g(b), then
E ( g (b))   g (b) (b)db
To get a numerical solution to the integral simply use
a random number generator to sample many b from
the distribution (b) and change the expectation to a
sample mean.
35
p y    p yi | bi  bi dbi
 E  p yi | bi 
1

N
 
~ 
 1
h 2 N / 2  N exp 2 yi  f bih 
1
where yi  f bi   [ yi  f bi ] R 1 [ yi  f bi ]
T
Where h is the index for the sampling from (bi).
~
and bih are the sampled bi
36
Problem!
If each evaluation of the likelihood surface requires a
resampling, then you introduce a randomness to
your likelihood surface.
The likelihood surface would have small perturbations
which would affect your determination of a
maximum.
Solution: sample once and re-use this sample for each
evaluation of the likelihood.
37
It turns out that importance sampling is not very
efficient. To improve on this method, another method
takes advantage of the normal assumption of
distribution of bi.
This method is called Gaussian Quadrature. Instead of a
random sample, specific abscissas have been
determined to best evaluate the integral.
In particular, “adaptive Gaussian Quadrature” is a
preferred method (not covered here).
38
Review of Approximate Methods
First order: biased, only useful for getting starting
values for better methods; converges often even if
model is horrible. DON’T RELY ON THIS METHOD!
Laplacian: numerically “cheap”, reasonably good fit
Importance sampling: Need lots of abscissas, so not
useful
Gaussian Quadrature: GOLD STANDARD! But when
data set large, method is slow and difficult to get
convergence.
39
Additional Note
When your model does not converge, often it’s because
you have the wrong model.
Don’t switch algorithms just because of
nonconvergence. First plot data and scrutinize choice of
model.
40
Software
R – PKFIT package
NONMEM (industry standard, 1979, FORTRAN)
Monolix (Bayesian)
WinBugs (PKBugs)
Phoenix (windows program incorporating methods
from NONMEM)
SAS and R can be used to fit very simple PK models but,
in general, not very useful
41
Why is NONMEM the gold standard?
Software needs easy input of PK models.
Not many software packages allow for models written
in terms of ODEs instead of closed form solution.
More challenging for multiple dose settings.
Functional form dependent on data.
42
Multiple Dose Model
Daily Dose with Fast Elimination
43
Multiple Dose Model
Daily Dose with Slower Elimination
Superposition
principle
44
Super-position Principle
Assume dosing every 24 hours
Assume concentration for single dose is
Then concentration, C(t) is
f (tij ,  i , di )
f (t ,  i , d i )
t  24


f (t ,  i , d i )  f (t  24,  i , d i )
24  t  48

C (t ,  i , d i )  
 f (t ,  i , d i )  f (t  24,  i , d i )  f (t  48,  i , d i ) 48  t  72



45
Multiple Dose Model, Missed Third Dose
46
Dose Delayed by 3 Hours Every Other Day
47
Modeling Covariates
Assumed: PK parameters vary with respect to a
patient’s weight or age.
Covariates can be added to the model in a secondary
structure (hierarchical model).
“Population Pharmacokinetics” refers specifically to
these mixed effects models with covariates included in
the secondary, hierarchical structure
48
Nonlinear Mixed Effects Model
With secondary structure for covariates:
yij  f (tij ,  i , d i )   ij
 i  g ( x ij , ai )  b i
b i ~ N 0, B 
 i ~ N 0, Ri 
49
Part 1 Summary
• PK Basic Concepts (Cmax, Tmax, AUC …)
• PK Compartmental Models (derived mean function
from differential equations)
• Technical Considerations (approximate maximum
likelihood approach for fitting nonlinear mixed effects
model)
• Software Considerations (many complications!)
50