Understanding PK/PD Modeling in Pharmacometrics

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

PK/PD Modeling Procedure: 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 Cmax Concentration 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 dA = c k A 10 c dt Central = t / C A V Vc c c c ( ) = = 0 Dose A c Elimination k10 Dose = k t ( ) C t e Solution: 10 c V c 9

First-Order 1-Compartment Model (Intravenous injection) Parameterized with Clearance Input Another parameterization for the solution uses Clearance = Cl = k10 Vc Central Vc Clearance = Volume of drug eliminated per unit time Elimination k10 Dose Cl = / t V ( ) C t e Solution: c V c 10

First-Order 1-Compartment Model (Extravascular Administration) Input dA = a k A ka Absorption depot: a a dt dA Central Central compartment: = c k A k A 10 a a c Vc dt = t / C A V c c c ( ( ) ) = = 0 dose A Elimination a = ) = 0 0 A t k10 c )( k F Dose 10 = k t k t ( ) a C t e e 10 a F = Bioavailability (i.e., fraction absorbed) Solution: ( c V k k c a 11

First-Order 1-Compartment Model (Extravascular Administration) Parameterized with Clearance Input = c Cl V k ka 10 Central Vc ( ) k F Dose ( ) = / Cl V t k t ( ) a V C t e e c a Solution: c k Cl c a 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) dA ( ) = + c k A k k A 21 12 10 p c dt Input dA p = k A k A k12 12 21 c p dt Peripheral Central = / C A V (Vp) Vc c c c = / C A V k21 p p p ( ( ) ) = = 0 Bolus Dose A t Elimination c = = 0 0 A t k10 p ( ) t = + exp( ) exp( ) Cc A t B t General form of solution: 14

First-Order 2-Compartment Model (Intravenous Dose) Parameterized in terms of Micro constants Note that including Vp over- parameterizes the model since Input k12 Peripheral Central k (Vp) Vc = 12 V V p c k k21 21 Elimination Ac = Amount of drug in central compartment k10 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 is ij f th th the j response for the i subject y a is scalar function nonlinear in a is parameter 1 vector k = + ( , , ) y f t d ij ij i ij i th th ( ( , 0 ) ) is ij d the j time for the i subject t ~ , N D th is i the i subject' dose s n i ranges from 1 to j ~ N R i i i residual is ij k error a is D covariance n matrix k is i R an covariance matrix n i i Pharmacokineticists use the term population model when the model involves random effects. 18

For simplification at this stage, assume = + b i i ~ ( , 0 ) b N D i and ( , 0 ) 2 ~ N I i i 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( ) where is an m x 1 vector which can incorporate , D and other parameters, e.g. R( )= 2[f( )]2, = 22

Error Models used for PK modeling ( ) + , f ijt Additive error i )( )( ij ( ) ) , + 1 f ( t e Proportional error ij i ij + + , 1 f t e Additive and Proportional error ij ( i ij ij ) , f t e ij Exponential error ij i 23

For the 1-compartment model parameterized with Cl, V, ka Input ka logk Central a Vc = = , logCl Cl k V 10 logV Elimination 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 = i 1 N ( ) p iy 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: ( ) ( ) b N N ( ) y = = 1 = | p p y b d b i i i i i 1 i i Where p and are normal probability density functions. Maximization is in = 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 ) ( ( i t f f 1 exp ( ( ) ( ) 1 T = 1 | exp p y b B y f R y f 1 i i 2 i i i i i ) ) + , , f t b d 1 i i i + , , b d 2 i = i i ( ) + , , f t b d in i i i ( ) b T = 1 B b D b 2 i i i 2 = and B B b constant w.r.t. 1 2 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 Available in NONMEM First Order Approximate just the integrand First Order Conditional Estimation ( ) ( ) b M =1 | p y b d b Laplace s Approximation i i i i i Importance Sampling Gaussian Quadrature 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 ( t f ) ( ) = + , , , , d f t b d ij i i ij i i And Ri = 2 , 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 ( ( i ij d t f + , , ) + , , f t b d ij i i ) ( ) b / , , | x f t x d = ij i x i Let this approximation be f ( ) , Tay , , t b d ij i i You use this approximation in the integrand. 32

Substituting back in and simplifying ( i i y B b y p 2 ) ( ) ( ) 1 T 1 Tay Tay | exp f R y f 1 i i i i i ( ( ) ) Tay , , , f t b d 1 i i i Tay , , , f t b d 2 i i i = Tay f i ( ) Tay , , , f t b d in i i i 1 ( ) b T = 1 exp B b D b See slide 26. 2 i i i 2 = and constant w .r.t. b B B 1 2 i And now the exponent term is integrable and now we can maximize the likelihood. 33

Using Laplaces Approximation A second order approximation can be constructed by using Laplace s approximation / 2 k ( ) b 2 ( ) b i / 1 2 ( ) b n exp n db e i i i i i i i i n i where ( ) b b maximizes i i i In this manner, the whole integral is approximated so no integration is needed. ( ) i n ( ) ( ) b 1 = log | b p y b i i i i i 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

( ) y ( ) ( ) = | p p y b b db i i i i ( ( ) ) 1 = | E 1 p y b i i ( ) ih b ~ 1 exp y f ( ) f i / 2 N N 2 N 2 h ( ) b ( ) b ( ) b T = 1 where y y f R y f i i i i i i Where h is the index for the sampling from (bi). ~ b b and sampled the are ih i 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 Super- position principle 44

Super-position Principle Assume dosing every 24 hours Assume concentration for single dose is Then concentration, C(t) is ( , , ) f t d ij i i ( , , t ) 24 f t d t i i + d ( , , t ) ( 24 , , t ) 24 48 f t d f d t i i = ( , , ) C t d i i i + + ( , , ) ( 24 , , ) ( 48 , , ) 48 72 f t d f f d t i i i i i i i 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: = + ( , , ) + y f t d ij ij ( i ij i = , ) g x a b ij i i i ( ( , 0 ) ) ~ , 0 b N B i ~ N R i i 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