Modern PK/PD Modeling Tools

In its most basic presentation, pharmacokinetic (PK) and pharmacodynamic (PD) data is simply the two dimensional relationship between either drug concentration and time or effect and drug concentration, respectively.  The purpose of PK/PD modeling is to 1) to summarize this data with simple statistics or functions with a relatively small number of parameters and 2) to use statistical techniques to distinguish random variability within and among subjects from fixed, structural aspects of drug disposition or effect.  The tools available for this modeling are various techniques of statistical analysis, usually (but not always) embodied in specific software packages.  

Pharmacokinetics

                        To fully characterize drug disposition,  PK models are, if not necessary, very helpful.  It is essential to formulate both “structural” and “error” models.  Structural models characterize the deterministic, or fixed, aspects of drug disposition. The best known and most commonly used structural models are compartment models assuming instantaneous mixing to  simplify the mathematical details. However, the assumption of instantaneous mixing is clearly incorrect and standard compartment models tell us nothing about the early phases of drug distribution. An  exciting recent advance is the development of recirculatory models, which can accurately predict this mixing phase.  These types of models are particularly suited for understanding how hemodynamic variables may impact on PK behavior.

            Given a specific structural model,  determination of PK parameters requires the assumption of an error model. Whatever the structural model, the concentrations predicted by the model and the measured concentrations will not agree.  It is common to assume that the relationship between measured and predicted (by the structural model) concentrations are related by the following equation

            Cm = Cp +  M

where Cm is the measured concentration, Cp is the predicted concentration and is a function of the structural PK parameters, and ε is the “error” term.  The term “error model” refers to the assumed statistical distribution of ε.  For example it is common to assume that ε has a normal distribution. One then obtains estimates of  PK parameters using maximum likelihood.  The principle of maximum likelihood states that the “best” estimate of parameters are those which maximize the likelihood of the observations.  If ε has a normal distribution then the probability of any single observation is

                        Pi = (1/√(2πvar) )EXP(-(Cm,i-Cp,i)2/(2 var))

where var is the variance of the error term (and can be a function of the structural parameters).  The probability of a series of observations is equal to P= Πi Pi.  Maximizing the likelihood of the observations is equivalent to maximizing the natural logarithm of likelihood.  In this case the logarithm of P is equal to

      log likelihood = log P = Σi Pi =  -nlog(√(2πvar) – Σi (Cm,i-Cp,i)2/(2var)

By specifying the error model, i.e., by assuming a model for var, and then maximizing the log likelihood with respective to the parameters, consistent estimates of the same can be obtained.  In the most basic application of this method, it is assume that var is a constant and the equations reduce to simple least squares. However, the assumption of a constant variance for the model error is often unacceptable and a more realistic error model must be used.   In this case the approach is known as extended least squares. When the error model is relatively straightforward this method is also very useful for the analysis of “dense” data, i.e., when there is sufficient data from each subject for determination of  PK parameters for the individual patient.

             The most prominent focus of modern PK is population analysis in which data from multiple individuals is pooled for analysis. . The best known software for population analysis is NONMEM, although Pharsight Inc also now offers software for population analysis. The underlying mathematics is challenging since it requires distinguishing inter-patient from intra-patient variability. But population analysis offers great advantages, allowing the investigator to use sparse data (few data points per individual),  and to assess inter-patient variability.  The mathematical challenge can be appreciated by considering likelihoods.  For any individual, the probability of the observed results is a function of that individual’s structural parameters (denoted S), i.e., Pi = f(S).  However, when we pool data from multiple patients we have to also consider that the values of S are different for different patients. To calculate the likelihood we have to multiply the above probability of the observed results for an individual  given particular parameter values by the probability of those particular values (denoted h(S)) and then sum (integrate) over all possible values of S.  In mathematical terms

            Li = ∫ Pi (S) h (S) dS

And the likelihood of the all observations is

            L = ∏i   Li

In general, the integrals involved in the calculation of the likelihood are intractable and simplifying assumptions are required.  It is imperative that the investigator understands the nature of these assumptions.

Pharmacodynamics

            The mainstay of PD modeling is the Hill, or sigmoid Emax, equation, which postulates the following relationship between drug concentration (C) and drug effect (E)

            E = Emax (Cγ/(Cγ+C50γ))

where Emax is the maximal effect, C50 is the drug concentration that results 50% of the maximal effect, and γ is the slope parameter that determines the slope of the concentration-response curve. When the effect is a continuous variable, estimates of Emax, C50 and γ are usually obtained by extended least squares or iteratively reweighted least squares when there is sufficient data for analysis of individual subjects.   When sparse data is pooled from multiple patients, then a population analysis is a better approach.  However, the more interesting situation is dichotomous data.  For example, in the analysis of anesthetic effect we often have binary data, the patient is either responsive or nonresponsive.  In this situation the analysis is more complicated. The drug effect is the probability of response or nonresponse.  Emax by definition is unity.  C50 and γ are typically estimated by maximum likelihood.  The likelihood of the observed results is

            L = Πi (PiRi)(1-Pi)1-Ri

where Ri is one if there is a positive anesthetic effect and zero if there is not, and

            Pi = Ciγ/(Ciγ + C50γ)

and Ci is the concentration at the time of the observation.  At this time it is unclear how reliable the analytic tools when there is sparse data (minimal number of data points per patient).  Simulations reveal that C50 may be accurately estimated either by naïve pooling of data or population analysis.  However estimates of γ may be significantly biased.