Statistical models in pharmacokinetics and pharmacodynamics

Trabajo Fin de Grado Statistical Models in Pharmacokinetics and Pharmacodynamics Presented by: Sandra Ben´ıtez Pe˜na Supervisors: Dr. Rafael Blanquero Bravo, Universidad de Sevilla Dr. Emilio Carrizosa Priego, Universidad de Sevilla June 21, 2015 Contents 1 Introduction 4 2 Nonlinear Mixed Effects Models 5 2.1 NonlinearModels......................... 5 2.2 Nonlinear Fixed Effects Models . . . . . . . . . . . . . . . . . 7 2.3 Nonlinear Random and Mixed Effects Models . . . . . . . . . 8 2.3.1 The Two-Stage Hierarchy . . . . . . . . . . . . . . . . 9 2.3.2 Summary ......................... 15 2.4 Inference in Nonlinear Mixed Effects Models . . . . . . . . . . 16 2.4.1 Inferential Approaches . . . . . . . . . . . . . . . . . . 17 2.4.2 Inference.......................... 19 3 Applications in Pharmacokinetics and Pharmacodynamics 25 3.1 Pharmacokinetics and Pharmacodynamics . . . . . . . . . . . 25 3.1.1 Pharmacokinetics . . . . . . . . . . . . . . . . . . . . . 27 3.2 Some Pharmacokinetic Models . . . . . . . . . . . . . . . . . . 33 3.2.1 One compartment models . . . . . . . . . . . . . . . . 34 3.2.2 Two compartments models . . . . . . . . . . . . . . . . 41 3.3 Inference.............................. 45 3.4 Some Pharmacodynamic Models . . . . . . . . . . . . . . . . . 48 3.4.1 Inmediate response models . . . . . . . . . . . . . . . . 48 3.4.2 Turnover response models . . . . . . . . . . . . . . . . 50 2 Chapter 1 Introduction Population Pharmacokinetic / Pharmacodymanic analysis using Nonlinear Mixed Effects Modeling is an approach that has gained importance in pharmacometrics. One of the principal purposes of this analysis, as we will see later, is to keep the concentration of a drug high enough to produce a desirable response, but low enough to avoid toxicity. For this objective, studies in population are made. Nonlinear Mixed Effects Modeling is the primary tool employed on this analysis, so the first objects that are introduced in this work are Nonlinear Models. This work is structured in two chapters apart from this introduction. As mentioned before, Nonlinear Models are introduced in Chapter 2. Then, with this element defined, and after we define Nonlinear Fixed and Random Effects Models, Nonlinear Mixed Effects Models are presented. After that, and also in Chapter 2, inference on those models is made. We expose some methods of inference that can be made, as Methods Based on Individuals Estimates or Methods Based on Approximations of the Likelihood. Then, in Chapter 3 we begin explaining what Pharmacokinetics and Pharmacodynamics are. Once we have made that, some Pharmacokinetic and Pharmacodynamic models are presented. Also, we make inference using [10] in some of those models. 4 Chapter 2 Nonlinear Mixed Effects Models 2.1 Nonlinear Models In order to explain what a nonlinear model is, we must remember the definition of linear model. Let •Y∈Rbe a variable of responses (called dependent variable), •X∈Rpa vector of predictor variables (called independent variables or covariates), •θ∈Rra vector of regression parameters, •ε∈Ra variable of random intra-individual errors. A model Y=f(X;θ) + ε, is said to be linear if for any X, the function θ7→ f(X;θ) is affine. In particular, such fis differentiable and its partial derivatives with respect to any of the model parameters ∂f ∂θkk=1,··· ,p 5 are independent of the other parameters. Otherwise, a model Y=f(X;θ) + ε is nonlinear if any of the partial derivatives with respect to any of the model parameters ∂f ∂θi are dependent on any other model parameter θj(j6=i) or if any of the derivatives do not exist or are discontinuous. Example 2.1.1 An example of nonlinear model is the Emax pharmacodynamic model, see [2], [8] or [1]. It has a two dimensional parameter θ= (θ1, θ2) = (Emax, EC50), where •Emax is the maximal effect, •EC50 is the concentration that produces 50% of the maximal effect. Define: •Eas the observed effect, •Cas the concentration. The Emax model is given by: Y=E=EmaxC EC50 +C=f(X;θ) + ε, where C≡X. The partial derivatives are given by ∂Y ∂θ1 =∂E ∂Emax =C EC50 +C and ∂Y ∂θ2 =∂E ∂EC50 =−EmaxC (EC50 +C)2. We note that both partial derivatives are dependent on the another parameter, so this model is effectively a nonlinear model. 6 Example 2.1.2 Another example of nonlinear model is f(X;θ) = θ0+θ1X, X ≤X0, θ2+θ3X, X > X0, see [2]. We must note that in those models the regression parameters themselves are not dependent on any other model parameter, but the estimates may be dependent on the value of X0. For those models, the derivative is not generally continuous at X=X0, so they are nonlinear models. Models that are transformed in predictor variables are not, however, nonlinear. To illustrate that, we use the following example. Example 2.1.3 Let f(X;θ) = θ0+θ1pX1+θ2sin(X2). be a regression model. This model is not nonlinear because the partial derivatives are not dependent on any other model parameter and the derivative is not discontinuous (in fact, this is a linear model), see [2]. Now, we know how nonlinear models are, so we can advance and study the next section, the Nonlinear Fixed Effects Models. 2.2 Nonlinear Fixed Effects Models Before we explain what nonlinear fixed effects models are, we are going to define what a fixed effect variable is. Afixed effect variable is one where the researcher can choose the level(s) of the variable to represent the precise contrast of interest. Alternatively, fixed effect variables are those variables whose levels in a study represent an exhaustive set of all possible levels. Later, we will define those fixed effects variables more mathematically, and they will be represented as the (q×1) vector β. 7 Example 2.2.1 The doses of drug used in a study or the time points that blood samples are measured represent fixed effects variables, see [2]. Example 2.2.2 Using both males and females in a study exhaust the possible levels of this variable, thus making sex a fixed effect, see [2]. With fixed variables defined, we can now define nonlinear fixed effects models as those that only contain fixed effects variables. So, as we have already defined nonlinear fixed effects models, we can continue to define random and mixed effects models. 2.3 Nonlinear Random and Mixed Effects Models As we made before with fixed effects models, we will explain what random effects variables are, in order to define nonlinear random effects models, and then we will define what mixed effects models are. Random effects variables are variables whose levels do not exhaust the set of possible levels and each level is equally representative of other levels. As we have made with fixed effects variables, we also explain later and more mathematically random effects variables, which will be denoted as the (k×1) vector bi. Example 2.3.1 The most commonly seen random effect variables in clinical research are the subjects used in an experiment, since in most cases researchers are not specifically interested in the particular set of subjects that were used in a study. Now we can define what nonlinear random and nonlinear mixed effects models are. Any nonlinear model that only contains random effects variables is a nonlinear random effects model. Any nonlinear model that contains both fixed and random effects variables is a nonlinear mixed effects model. 8 one might regard the first term in (2.9) as negligible. This leads to a standard model in this application, where Ri(ui,θi,ξ) = Λi(ui,θi,λ), with Λi(ui,θi,λ) a diagonal matrix with diagonal elements σ2(f(ti,j,ui,θi))2λ for some λ= (σ2, λ)0; often λ= 1. 2.3.2 Summary We are now in a position to summarize the basic nonlinear mixed effects model. Let fi(ui;θi) =    f(Xi,1;θi) . . . f(Xi,ni;θi)   , and let zi= (u0 i,a0 i)0summarize all covariate information on subject i. Then, we may write the model in (2.1) and (2.2) succinctly as Stage 1: E[Yi|zi,bi] = fi(ui;θi) = fi(zi;β,bi) V ar[Yi|zi,bi] = Ri(ui,θi,ξ) = Ri(zi,β,bi,ξ) (2.10) Stage 2: θi=d(ai,β,bi),bi∼(0,G) (2.11) In (2.10), dependence of fiand Rion the covariates aiand fixed and random effects through θiis emphasized. This model represents individual behavior conditional on θiand hence on bi, the random component in (2.11). In (2.11), we assume that the distribution of bi|aidoes not depend on ai, so that all bihave common distribution with mean 0and covariance matrix G. Within-individual correlation: The nonlinear model (2.10)-(2.11) implies a model for the marginal mean and covariance matrix of Yigiven all covariates zi; i.e. averaged across population. 15 Letting Fb(bi) denote the cumulative distribution function of bi, we have E[Yi|zi] = Zfi(zi;β,bi)dFb(bi) V ar[Yi|zi] = E[Ri(zi,β,bi,ξ)|zi] + V ar[fi(zi,β,bi)|zi],(2.12) where expectation and variance are with respect to the distribution of bi. In (2.12), E[Yi|zi] characterizes the typical response profile among individuals with covariates zi. V ar[Yi|zi] involves two terms as we have seen in (2.12): –E[Ri(zi,β,bi,ξ)|zi], which averages realization and measurement variation that occur within individuals across individuals having covariates zi –V ar[fi(zi,β,bi)|zi], which describes how inherent trajectories vary among individuals which have the same zi. Note that E[Ri(zi,β,bi,ξ)|zi], is a diagonal matrix only if Γi(ρ) in (2.9), reflecting correlation due to within-individual realizations, is an identity matrix. However, V ar[fi(zi,β,bi)|zi], has non-zero offdiagonal elements in general due to common dependence of all elements of fion bi. Thus, correlation at the marginal level is always expected due to variation among individuals, while there is correlation from within-individual sources only if serial associations among intraindividual realizations are nonnegligible. In general, then, both terms contribute to the overall pattern of correlation among responses on the same individual represented in V ar[Yi|zi]. In many applications, the effect of within-individual serial correlation reflected in the first term of V ar[Yi|zi] is dominated by that from among individual variation in V ar[fi(zi,β,bi)|zi]. This explains why many published applications of nonlinear mixed models adopt simple, diagonal models for Ri(ui,θi,ξ) that emphasize measurement error. 2.4 Inference in Nonlinear Mixed Effects Models In this section we will introduce the theory behind nonlinear mixed effects models. Modeling and estimation of model parameters will be discussed, as 16 the relationship between the covariate and the dependent variable is modeled. 2.4.1 Inferential Approaches The main inferential objectives are: •Determine an appropriate population model d(ai,β,bi) and inference on elements of βas central interests. •Inference on Gand, in particular, on diagonal elements. •Inference on θiand f(t0,ui,θi) at some specific time point t0, for i∈ {1,2,· · · , m}. The nonlinear mixed effects model (2.10)-(2.11) is a so called subjectspecific model. The distinction between subject-specific and population-averaged (or marginal) models may not be important for linear mixed effects models, but it is critical under nonlinearity. A population-averaged model assumes that interest focuses on parameters that describe the marginal distribution of Yigiven covariates zi. From the discussion following (2.12), if E(Yi|zi) were modeled directly as a function of ziand a parameter β,βwould represent the parameter corresponding to the typical response profile among individuals with covariates zi. This is to be contrasted with the meaning of βin (2.11) as the typical value of individual-specific parameters θiin the population. Consider first linear models. A linear subject-specific model with second stage θi=Aiβ+Bibi as in (2.3), for design matrix Aidepending on ai, and first stage E[Yi,j|ui,θi] = Uiθi, where Uiis a design matrix depending on the ti,j and ui, leads to the linear mixed effects model E[Yi|zi,bi] = fi(zi,β,bi) = e Xiβ+Zibi for e Xi=UiAi 17 and Zi=UiBi where e Xithus depends on zi. From (2.12), this model implies that E[Yi|zi] = Z(e Xiβ+Zibi)dFb(bi) = e Xiβ, as E[bi] = 0. Thus, in a linear subject-specific model, βfully characterizes both the typical value os θiand the typical response profile, so that the interpretation is valid. Here, then, postulating the linear subject-specific model is equivalent to postulating a population-averaged model of the form E[Yi|zi] = e Xiβdirectly, in that both approaches yield the same representation of the marginal mean and hence allow the same interpretation of β. Consequently, the distinctions between subject-specific and population-averaged approaches have not generally been of concern in the literature on linear modeling. For nonlinear models, however, this is not longer the case. For instance, suppose that bi∼ N(0,G), and consider a subject-specific model of the form in (2.10) and (2.11) for some function fnonlinear in θiand hence bi. Then, from (2.12), the implied marginal mean is E[Yi|zi] = Zfi(zi,β,bi)p(bi;G)dbi,(2.13) where p(bi;G) is the N(0,G) density. For nonlinear f, this integral is clearly intractable and E[Yi|zi] is an expression that may not be available in a closed form and depends on both βand Gin general. As a result, if we start with a nonlinear subject-specific model, the implied population-averaged marginal mean model involves both the typical value of θi(β) and G. Accordingly, βdoes not fully characterize the typical response profile and thus cannot enjoy both interpretations. Conversely, if we were to take a populationaveraged approach and model the marginal mean directly as a function of ziand a parameter β,βwould indeed have the interpretation of describing the typical response profile. But it seems unlikely that it could also have the interpretation as the typical value of specific parameters θiin a subjectspecific model. Thus, for nonlinear models, the interpretation of βin subjectspecific and population-averaged models cannot be the same in general. The implication is that the modeling approach must be carefully considered to ensure that the interpretation of βcoincides with the questions of scientific interest. 18 2.4.2 Inference A number of inferential methods for the nonlinear mixed effects model are now in common use, see [4] and [7]. We provide a brief overview that can be seen below, beginning with Maximum Likelihood. Maximum Likelihood A natural starting point for inference is Maximum Likelihood. This is only a starting point here because the analytical intractability of likelihood inference has motivated many approaches on approximations. The individual model (2.10) along with an assumption on the distribution of Yigiven (zi,bi) yields a conditional density p(Yi|zi,bi;β,ξ), say; the ubiquitous choice is the normal. Under the popular assumption that Ri(zi,β,bi,ξ) is diagonal, the density may be written as the product of mcontributions p(Yi,j|zi,bi,β,ξ). Under this condition, the lognormal has also been used. At Stage 2, (2.11), adopting independence of biand ai, one assumes a k-variate density p(bi,G) for bi. As with other mixed models, normality is standard. With these specifications, the joint density of (Yi,bi) given ziis p(Yi,bi|zi;β,ξ,G) = p(Yi|zi,bi;β,ξ)p(bi;G) (2.14) A likelihood for β,ξ,Gmay be based on the joint density of the observed data Y1,· · · ,Ymgiven zi, m Y i=1 Zp(Yi,bi|zi;β,ξ,G)dbi= m Y i=1 Zp(Yi|zi,bi;β,ξ)p(bi;G)dbi(2.15) by independence across i. Nonlinearity implies that the m k-dimensional integrations in (2.15) generally cannot be done in a closed form; thus, iterative algorithms to maximize (2.15) in β,ξ,Grequire a way to handle these integrals. Although numerical techniques for evaluation of an integral are available, these can be computationally expensive when performed at each internal iteration of the algorithm. Hence, many approaches to fitting (2.10)- (2.11) are instead predicated on analytical approximations. Methods Based On Individual Estimates If the niare sufficiently large, an intuitive approach is to summarize the responses Yifor each ithrough individual-specific estimates ˆ θiand then 19 use these as the basis for inference on βand G. In particular, viewing the conditional moments in (2.10) as functions of θi, i.e., E[Yi|ui,θi] = f(Xi,j,θi), V ar[Yi|ui,θi] = Ri(ui,θi,ξ), fit the model specified by these moments for each individual. Usual large-sample theory implies that the individual estimators ˆ θiare asymptotically normal. Each individual is treated separately, so the theory may be viewed as applying conditionally on θifor each i, yielding ˆ θi|ui,θi∼ N(θi,Ci). Because of the nonlinearity of fin θi,Cidepends on θiin general, so Ciis replaced by ˆ Ciin practice, where ˆ θiis substituted. To see how this is exploited for inference on βand G, consider the linear second-stage model (2.3); the same developments apply to any general din (2.11). The asymptotic result may be expressed alternatively as θ∗ i≈θi+ε∗ i=Aiβ+Bibi+ε∗ i, ε∗ i|zi∼ N(0,ˆ Ci),bi∼ N(0,G) (2.16) where ˆ Ciis treated as known for each i, so that the ε∗ ido not depend on bi. This is the form of a linear mixed effects model with known error covariance matrix ˆ Ci, which suggest using standard techniques for fitting such models to estimate βand G. Some authors propose use of the EM algorithm in the case k=pand Bi=Ip. Alternatively, it is possible to use linear mixed model software to fit (2.16). EM Algorithm The EM algorithm is a procedure to compute the maximum likelihood estimates in the presence of missing or hidden data. Each iteration of the EM algorithm consists of two processes: 1. the E-step, in which the missing data are estimated given the observed data and current estimate of the model parameters. This is achieved using the conditional expectation, 20 2. the M-step, in which the likelihood function is maximized under the assumption that the missing data are known. The estimate of the missing data from the E-step are used instead of the actual missing data. Convergence is assured since the algorithm is guaranteed to increase the likelihood at each iteration. For more details, see [3]. In the case treated in this work, the EM-algorithm can be specified as follows, see [9]. At (k+ 1)th iteration, 1. E-step: Update the estimates of θi, i = 1,2,· · · , m as ˆ θi,(k+1) = ( ˆ C−1 i+ˆ G−1 (k))−1(ˆ C−1 iθ∗ i+ˆ G−1 (k)Aiˆ β(k)), i = 1,2,· · · , m. 2. M-step: Update the estimates of the population parameters as ˆ β(k+1) = m X i=1 A0 iˆ G−1 (k)Ai!−1m X i=1 A0 iˆ G−1 (k)ˆ θi,(k+1), ˆ G(k+1) =m−1 m X i=1 (ˆ θi,(k+1) −Aiˆ β(k+1))(ˆ θi,(k+1) −Aiˆ β(k+1))0+ +m−1 m X i=1 (ˆ C−1 i+ˆ G−1 (k))−1 Note that as starting values for the EM algorithm one can use the following estimates: • ˆ β(0) = m X i=1 A0 iAi!−1 m X i=1 A0 iθ∗ i!, • ˆ G(0) = (m−1)−1 m X i=1 (θ∗ i−Aiˆ β(0))(θ∗ i−Aiˆ β(0))0 This algorithm is iterated until the M-step converges. 21 Methods Based On Approximation of the Likelihood This last point is critical when the niare not large. Although this provides rich information for building population models d, there are insufficient data to fit the pharmacokinetic model fto any one subject. Implementation of (2.15) in principle imposes no requirements on the magnitude of the ni. Thus, an attractive strategy is instead to approximate (2.15) in a way that intractable integration is avoided. In particular, for each i, an approximation to p(Yi|zi;β,ξ,G) = Zp(Yi|zi,bi;β,ξ)p(bi;G)dbi is obtained. First Order Methods An approach is motivated by letting R1/2 ibe the Cholesky decomposition of Riand writing (2.10)-(2.11) as Yi=fi(zi,β,bi) + R1/2 i(zi,β,bi,ξ)i,i|zi,bi∼(0,Ini).(2.17) As nonlinearity in bicauses the difficulty for integration in (2.15), it is natural to consider a linear approximation. A Taylor series of (2.17) about bi=0 to linear terms, disregarding the term involving biias small and letting Zi(zi,β,b∗) = ∂ ∂bi (fi(zi,β,bi)) |bi=b∗ leads to Yi≈f(zi,β,0) + Zi(zi,β,0)bi+R1/2 i(zi,β,0,ξ)i(2.18) E[Yi|zi]≈fi(zi,β,0), V ar[Yi|zi]≈Zi(zi,β,0)GZ0 i(zi,β,0) + Ri(zi,β,0,ξ).(2.19) When p(Yi|zi,bi;β,ξ) in (2.15) is a normal density, (2.18) amounts to approximating it by another normal density whose mean and covariance matrix are linear in and free of bi, respectively. If p(bi;G) is also normal, the integral is analytically calculable analogous to a linear mixed model and yields a ni-variate normal density p(Yi|zi;β,ξ,G) for each iwith mean and covariance matrix (2.19). This suggests the proposal to estimate β,ξ,Gby jointly maximizing m Y i=1 p(Yi|zi;β,ξ,G), 22 which is equivalent to maximum likelihood under the assumption the marginal distribution Yi|ziis normal with moments (2.19). The advantage is that this approximate likelihood is available in a closed form. Standard errors are obtained from the information matrix assuming the approximation is exact. As (2.19) defines an approximate marginal mean and covariance matrix for Yi given zi, an alternative approach is to estimate β,ξ,Gby solving generalized estimating equations. An obvious drawback of all first order methods is that the approximation may be poor, as they essentially replace E[Yi|zi] = Zf(zi,β,bi)p(bi;G) by f(zi,β,0). This suggest that more refined approximation would be desirable, like First Order Conditional Methods. First Order Conditional Methods As p(Yi|zi,bi;β,ξ) and p(bi;G) are ordinarily taken as normal densities, a natural way to approximate integrals like those in (2.15) is to exploit Laplace’s method, a standard technique to approximate an integral of the form Re−l(b)dbthat follows from Taylor series expansion of −l(b) about the value ˆ b, say, maximizing l(b). The result is that p(Yi|zi;β,ξ,G) may be approximated by a normal density with E(Yi|zi)≈fi(zi,β,ˆ bi)−Zi(zi,β,ˆ bi)ˆ bi V ar[Yi|zi]≈Zi(zi,β,ˆ bi)GZ0 i(zi,β,ˆ bi) + Ri(zi,β,ξ) (2.20) ˆ bi=GZ0 i(zi,β,ˆ bi)Ri(zi,β,ξ)[Yi−fi(zi,β,ˆ bi)],(2.21) where Ziis defined as before, and ˆ bimaximizes l(bi)=[Yi−fi(zi,β,bi)]0R−1 i(zi,β,ξ)[Yi−fi(zi,β,bi)] + b0 iGbi in bi. In fact, ˆ bimaximizes in bithe posterior density for bi p(bi|Yi,zi;β,ξ,G) = p(Yi|zi,bi;β,ξ)p(bi;G) p(Yi|zi;β,ξ,G).(2.22) Equations (2.20) and (2.21) suggest an iterative scheme whose essential steps are 23 1. Given current estimates ˆ β,ˆ ξ,ˆ Gand ˆ bi, say, update ˆ biby substituting these in the right hand side of (2.21) 2. Holding ˆ bifixed, update estimation of β,ξ,Gbased on the moments in (2.20) It is well-documented by numerous authors that these first order conditional approximations work extremely well in general, even when niare not large or the assumptions of normality that dictate the form of (2.22) on which ˆ biis based are violated. These features and the availability of supported software have made this approach probably the most popular way to implement nonlinear mixed models in practice. Some final remarks follow. The methods in this section may be implemented for any ni. Although they involve closed-form expressions for p(Yi|zi;β,ξ,G) and moments (2.19) and (2.20), maximization or solution of likelihoods or estimating equations can still be computationally challenging, and selection of suitable starting values for the algorithms is essential. Results from first order methods may also be used as starting values for a more refined “conditional” fit. A common practical strategy is to first fit a simplified version of the model and use the results to suggest starting values for the intended analysis. For instance, one might take Gto be a diagonal matrix, which can often speed convergence of the algorithms; this implies the elements of θiare uncorrelated, which is usually highly unrealistic. 24 A(0) = 0 (3.2) Aa(0) = D, where D is the dose. It is easily seen that the Cauchy problem given by (3.1) under the initial value (3.2) has as solution A(t) = kaD ka−kee−ket−e−kat(3.3) Aa(t) = De−kat So concentration at time point t in the compartment is finally given by C(t) = A(t) V=kaD V(ka−ke)e−ket−e−kat Individual pharmacokinetics behavior is characterized by the three-dimensional pharmacokinetics parameter θ=  ka V Cl   Knowing the value of θallows one to determine concentrations achieved at any time point tunder different doses, because the pharmacokinetics model assumes that pharmacokinetics processes are dose-independent. So it can be used to develop dosing regimens. Now it can be studied its therapeutic window. It must be calculated the dose (D) which allows concentration to be bigger than a certain concentration, say Cmbut lower than another concentration, CM. First of all it must be calculated the global maxima of C(t). Then, at the time point t∗in which C(t) reaches its maximum is imposed to be Cm≤C(t∗)≤CM. After some simple calculations, it can be observed that such time point t∗is given by the expression t∗= ln ka ke ka−ke . 31 If it is desired to be C(t∗)≤CM, after some operations it is easily obtained that it must be D≤CMV(ka−ke) ka(e−ket∗−e−kat∗) Equivalently, if it is wanted to be C(t)≥Cm, it must be D≥CmV(ka−ke) ka(e−ket∗−e−kat∗) So, in summary, in order to be concentration inside the therapeutic window, dose must be D∈CmV(ka−ke) ka(e−ket∗−e−kat∗),CMV(ka−ke) ka(e−ket∗−e−kat∗). Example 3.1.2 Argatroban PK/PD study: Note: Argatroban is an anticoagulant that is a small molecule direct thrombin inhibitor. For the study described in [5], a sample of N= 37 subjects is used. Different constant infusion rates (doses) of 1to 5µ/Kg/min of argatroban are assigned. This infusion was administered by intravenous infusion for 4 hours (240 minutes). Blood samples for pharmacokinetics study were collected at {30,60,90,115,160,200,240,250,260,275,295,320}minutes. For pharmacodynamics study additional samples at 5-9 points, measuring activated partial thromboplastin time (aPTT, the response), were collected. The effect site of this example is blood. Argatroban pharmacokinetic model: Let the function (·)+be defined by x+= 0 if x≤0and x+=xif x > 0. Denoting by D,Cl,Tinf respectively the dose, clearance and duration of the infusion, the Argatroban concentration (C(t)) over time at the effect compartment is given by the equation C(t) = D Cl he−Cl V(t−tinf )+−e−Cl Vti, so in this model we have a two-dimensional parameter θ, namely, θ=Cl V 32 We must estimate typical values of θand how they vary in the population of subjects and understand the relationship between concentrations and response. Argatroban pharmacodynamic model: Argatroban pharmacodynamics model follows the Emax model presented in Example 2.1.1. E(t) = EmaxC(t) EC50 +C(t) The ultimate objective is to put pharmacokinetics and pharmacodynamics together to characterize the therapeutic window and how it varies across subjects and for develop dosing regimens targeting achieved concentrations leading to therapeutic response. 3.2 Some Pharmacokinetic Models In this section, we introduce some examples of pharmacokinetic models, see [1]. Some of them are models that we have already presented, but we also show other new models. We use the package PKfit in R, see [10], to simulate data of all this set of models, and then we are going to do the reverse: given some data, we try to fit the data into an specific model. Later in this work, some pharmacodynamic models are shown, but not in the same extension we present pharmacokinetic models. The presentation of the pharmacokinetic models is organised as follows: •First level: number of compartment –One compartment –Two compartments •Second level: route of administration –Intravenous (IV) bolus –Infusion –First order absorption –Zero order absorption 33 •Third level: elimination process –Linear –Michaelis-Menten •Fourth level: existence of a lag time for first and zero order absorption only •Last level: administration profile The equations presented here express the concentration C(t) in the central compartment at a time tafter a dose Dgiven at time tD(t≥tD). Note: For infusion, the duration of infusion is Tinf and Dis the total dose administrated. Note: For models with 1 and 2 compartments (as the models presented in this work), equations C(t) express concentration in the central compartment at a time tafter drug administration. PK/PD analysis, with intermediate response models, can use concentration C(t) in the central compartment but alternatively concentration Ce(t) in the effect compartment. There is an additional parameter to estimate, ke0, the equilibrium rate constant between central and effect compartment. For each model, the equation for Ce(t) is given after the corresponding one for C(t). 3.2.1 One compartment models Parameters •V= volume of distribution •ke= elimination rate constant •Cl = clearance of elimination •Vm= maximum elimination rate (in amount per time unit) •Km= Michaelis-Menten constant (in concentration unit) •ka= absorption rate constant •Tlag = lag time 34 •Tk0= absorption duration for zero order absorption Parameterisation. There are two parameterisations for one compartment models, (Vand ke) or (Vand Cl). The equations are given for the first parameterisation (V,ke). The equations for the second parameterisation (V,Cl) are derived using ke=Cl V. Example 3.2.1 IV bolus. Linear elimination. C(t) = D Ve−ke(t−tD) Ce(t) = D V ke0 (ke0−ke)(e−ke(t−tD)−e−ke0(t−tD)) R simulation data. Total subjects=6 D=300 ke=0.21 V=11.7 Time points: 0,0.1,0.2,0.3,0.4,0.6,0.8,1,2,4,6,8,12,14,16,18,24 35 0 5 10 15 20 0 5 10 15 20 25 30 Time(hours) Concentration (mg/L) 0 5 10 15 20 0 5 10 15 20 25 30 Time(hours) Concentration (mg/L) 0 5 10 15 20 0 5 10 15 20 25 30 Time(hours) Concentration (mg/L) Figure 3.3: Simulation IV Bolus. Linear elimination. Example 3.2.2 IV bolus. Michaelis-Menten elimination. Initial conditions:      C(t) = 0 ∀t<tD Ce(t) = 0 ∀t≤tD C(tD) = D V dC dt = Vm VC Km+C dCe dt =ke0(C−Ce) 36 R simulation data. Total subjects=6 D=300 Km=4.84 V=11.7 Vm=2.17 Time points: 0,0.1,0.2,0.3,0.4,0.6,0.8,1,2,4,6,8,12,14,16,18,24 0 5 10 15 20 0 5 10 15 20 25 30 Time(hours) Concentration (mg/L) Figure 3.4: Simulation IV Bolus. Michaelis-Menten elimination. Example 3.2.3 IV infusion. Linear elimination. C(t) =        D Tinf 1 keV(1 −e−ke(t−tD))if t−tD≥Tinf D Tinf 1 keV(1 −e−keTinf )e−ke(t−tD−Tinf )else 37 Ce(t) =            D Tinf 1 keV(ke0−ke)[ke0(1 −e−ke(t−tD))−ke(1 −e−ke0(t−tD))] if t−tD≤Tinf D Tinf 1 keV(ke0−ke)[ke0(1 −e−keTinf )e−ke(t−tD−Tinf )− −ke(1 −e−ke0Tinf )e−ke0(t−tD−Tinf )]else R simulation data. Total subjects: 6 Dose: 300 Tinf : 1 Ke: 0.12 V: 11.7 Time dose: 0,0.1,0.2,0.3,0.4,0.6,0.8,1,2,4,6,8,12,14,16,18,24,48,72 0 10 20 30 40 50 60 70 0 5 10 15 20 25 Time (hours) Concentration (mg/L) 0 10 20 30 40 50 60 70 0 5 10 15 20 25 0 10 20 30 40 50 60 70 0 5 10 15 20 25 0 10 20 30 40 50 60 70 0 5 10 15 20 25 0 10 20 30 40 50 60 70 0 5 10 15 20 25 0 10 20 30 40 50 60 70 0 5 10 15 20 25 Figure 3.5: Simulation IV Infusion. Linear elimination. Example 3.2.4 IV infusion. Michaelis-Menten elimination. 38 Initial condition: C(t) = 0 ∀t<tD Ce(t)=0∀t<tD dC dt = Vm VC Km+C+input dCe dt =ke0(C−Ce) input(t) =   D Tinf 1 Vif 0≤t−tD≤Tinf 0else Example 3.2.5 First order absorption. Linear elimination. In absence of a lag time. C(t) = D V ka ka−ke (e−ke(t−tD)−e−ka(t−tD)) Ce(t) = Dkake0 V(e−ka(t−tD) (ke−ka)(ke0−ka)+e−ke(t−tD) (ka−ke)(ke0−ke)+e−ke0(t−tD) (ka−ke0)(ke−ke0)) Example 3.2.6 First order absorption. Linear elimination. In presence of a lag time. C(t) =        0if t −tD≤Tlag D V ka ka−ke (e−ke(t−tD−Tlag)−e−ka(t−tD−T lag))else Ce(t) =                    0if t −tD≤Tlag Dkake0 Ve−ka(t−tD−Tlag) (ke−ka)(ke0−ka)+e−ke(t−tD−Tlag) (ka−ke)(ke0−ke)+ +e−ke0(t−tD−Tlag) (ka−ke0)(ke−ke0)else Example 3.2.7 First order absorption. Michaelis-Menten elimination. In absence of a lag time. 39 Initial condition: C(t)=0for t<tD Ce(t) = 0 for t < tD dC dt =− Vm VC Km+C+input dCe dt =ke0(C−Ce) input(t) = D Vkae−ka(t−tD) Example 3.2.8 First order absorption. Michaelis-Menten elimination. In presence of a lag time. Initial condition: C(t) = 0 for t<tD Ce(t) = 0 for t < tD dC dt =− Vm VC Km+C+input dCe dt =ke0(C−Ce) input(t) = (0if t −tD< Tlag D Vkae−ka(t−tD−Tlag)else Example 3.2.9 Zero order absorption. Linear elimination. In absence of a lag time. C(t) =          D Tk0 1 keV[1 −e−ke(t−tD)]if t −tD≤Tk0 D Tk0 1 keV[1 −e−keTk0]e−ke(t−tD−Tk0)else Ce(t) =                D Tk0 1 keV(ke0−ke)[ke0(1 −e−ke(t−tD))−ke(1 −e−ke0(t−tD))] if t −tD≤Tk0 D Tk0 1 keV(ke0−ke)[ke0(1 −e−keTk0)e−ke(t−tD−Tk0)− −ke(1 −e−ke0Tk0)e−ke0(t−tD−Tk0))] else 40 we obtain again a table that begins like that: Subject Time Concentration 1 0 0 1 0.1 5.016 1 0.2 9.964 1 0.3 14.885 1 0.4 19.791 . . .. . .. . . So we can consider again the data above as data in which we can make inference for determine Vm,Kmand V. We can choose as initial values •D= 300 (We consider the dose as known) •Tinf = 0.5 (We consider the infusion time as known too) •Vm= 2 •Km= 4 •V= 11, and R give us as a result the following: 47 --- initial values for parameters ---- Parameter Initial 1 D 300 2 Tinf 0.5 3 Vm2 4 Km4 5 V 11 --- wighting scheme: equal weight --- model selection: one-compartment, iv infusion PK model with M-M elim. << PK parameter obtained from Nelder-Mead Simplex algorithm >> Parameter Value 1 Vm2.183 2 Km4.792 3 V 11.670 As we see, we have get once more similar results as the original true values. 3.4 Some Pharmacodynamic Models Two different type of models are presented here: •Inmediate response models •Turnover models 3.4.1 Inmediate response models For these type of models, the effect (denoted as E(t)) is expressed as: E(t) = A(t) + S(t) (3.4) where A(t) represents the model of drug action and S(t) corresponds to the baseline/disease model. A(t) is a function of the concentration C(t) in the central compartment or of the concentration Ce(t) in the effect compartment. Parameters 48 •Alin = constant associated to C(t) •Aquad = constant associated to the square of C(t) •Alog = constant associated to the logarithm of C(t) •Emax = maximal agonistic response •Imax = maximal antagonistic response •EC50 = concentration to get half of the maximal response •γ= sigmoidicity factor •S0= baseline value of the studied effect •kprog = rate constant of disease progresion Drug action models •linear model A(t) = AlinC(t) •quadratic model A(t) = AlinC(t) + AquadC(t)2 •logarithmic model A(t) = Aloglog(C(t)) •Emax model A(t) = EmaxC(t) C(t) + EC50 •sigmo¨ıd Emax model A(t) = EmaxC(t)γ C(t)γ+ECγ 50 •Imax model A(t) = 1 −ImaxC(t) C(t) + EC50 •sigmo¨ıd Imax model A(t) = 1 −ImaxC(t)γ C(t)γ+ECγ 50 49 Baseline/disease models •null baseline S(t) = 0 •constant baseline with no disease progression S(t) = S0 •linear disease progression S(t) = S0+kprogt •exponential disease increase S(t) = S0e−kprogt •exponential disease decrease S(t) = S0(1 −e−kprogt) Note: Only, for the Imax and sigmo¨ıd Imax models A(t) is not added to S(t) but S0is multiplied by A(t) in the expression of S(t). For instance, for Imax model with linear baseline we have E(t) = S0∗A(t) + kprogt. 3.4.2 Turnover response models In these models, the drug is not acting on the effect Edirectly but rather on Rin (input (synthesis) rate) or kout (output (elimination) rate constant). Thus the system is described with differential equations, given dE dt as a function of Rin,kout and C(t), the drug concentration at time t. The initial condition is: while C(t) = 0, E(t) = Rin kout . Parameters: •Emax = maximal agonistic response •Imax = maximal antagonistic response 50 •EC50 = concentration to get half of the maximal response •γ= sigmoidicity factor •Rin = input (synthesis) rate •kout = output (elimination) rate constant Models with impact on the input (Rin) •Emax model dE dt =Rin 1 + EmaxC C+EC50 −koutE •sigmo¨ıd Emax model dE dt =Rin 1 + EmaxCγ Cγ+ECγ 50 −koutE •Imax model dE dt =Rin 1−ImaxC C+EC50 −koutE •sigmo¨ıd Imax model dE dt =Rin 1−ImaxCγ Cγ+ECγ 50 −koutE •full Imax model dE dt =Rin 1−C C+EC50 −koutE •sigmo¨ıd full Imax model dE dt =Rin 1−Cγ Cγ+ECγ 50 −koutE 51 Models with impact on the output (kout) •Emax model dE dt =Rin −kout 1 + EmaxC C+EC50 E •sigmo¨ıd Emax model dE dt =Rin −kout 1 + EmaxCγ Cγ+ECγ 50 E •Imax model dE dt =Rin −kout 1−ImaxC C+EC50 E •sigmo¨ıd Imax model dE dt =Rin −kout 1−ImaxCγ Cγ+ECγ 50 E •full Imax model dE dt =Rin −kout 1−C C+EC50 E •sigmo¨ıd full Imax model dE dt =Rin −kout 1−Cγ Cγ+ECγ 50 E 52 Bibliography [1] Julie Bertrand and France Mentr´e. Mathematical expressions of the pharmacokinetic and pharmacodynamic models implemented in the monolix software. Technical report, Technical 10 report, INSERM, 2008. [2] Peter L. Bonate. Pharmacokinetic-Pharmacodynamic Modeling and Simulation. Springer, 2nd edition, 2011. [3] Sean Borman. The expectation maximization algorithm-a short tutorial. 2004. [4] Marie Davidian. Nonlinear Mixed Effects Models: An Overview and Update. Journal of Agricultural, Biological, and Environmental Statistics, 2004. [5] Garrett Fitzmaurice, Marie Davidian, Geert Verbeke, and Geert Molenberghs. Longitudinal data analysis. CRC Press, 2008. [6] Mats O Karlsson, Stuart L Beal, and Lewis B Sheiner. Three new residual error models for population pk/pd analyses. Journal of pharmacokinetics and biopharmaceutics, 23(6):651–672, 1995. [7] E. Kuhn and M. Lavielle. Maximum likelihood estimation in nonlinear mixed effects models. Computational Statistics & Data Analysis, 49(0):1020–1038, 2005. [8] B. Meibohm and H. Derendorf. Basic concepts of pharmacokinetic/pharmacodynamic (PK/PD) modeling. International Journal of Clinical Pharmacology and Therapeutics, 35(0):401–413, 1997. [9] Hulin Wu and Jin-Ting Zhang. Nonparametric regression methods for longitudinal data analysis: mixed-effects modeling approaches, volume 515. John Wiley & Sons, 2006. [10] Chun ying Lee and Yung jin Lee. PKfit: A Data Analysis Tool for Pharmacokinetics, 2014. R package version 1.2.4. 53