Why Fit PK Before PD?
A pharmacodynamic response is often modeled as a function of drug concentration, but concentrations may not be observed at every time a response is measured. A two-stage procedure addresses this by first fitting a pharmacokinetic (PK) model separately for each subject, then using each fitted PK model to predict concentrations at the PD observation times. A population PD model is estimated from those predictions.
The procedure separates the computational tasks: subject-specific PK parameters are obtained from PK data, and common PD parameters are obtained from the relationship between predicted concentrations and observed responses. The resulting PD estimate is conditional on the chosen PK and PD models and on the predicted concentrations used in stage two. This makes the procedure straightforward to describe, but each stage requires careful specification and checking.
Model and Estimation Steps
Let \(Y^{PK}_{ij}\) be the measured concentration for subject \(i\) at time \(t_{ij}\). Suppose the PK model predicts concentration \(f(t_{ij};\theta_i)\), where \(\theta_i\) contains that subject’s PK parameters. Under a normal residual model, one possible stage-one specification is
For each subject, estimate \(\theta_i\) from that subject’s PK data. With constant residual variance and independent normal errors, maximum likelihood is equivalent to minimizing the sum of squared residuals for that subject. Other residual models, such as proportional error, require a different likelihood or fitting criterion.
For each PD observation at time \(s_{ik}\), calculate the concentration prediction \(\hat{C}_{ik}=f(s_{ik};\hat{\theta}_i)\). A PD model might be written \(Y^{PD}_{ik}=g(\hat{C}_{ik};\psi)+\varepsilon^{PD}_{ik}\), where \(\psi\) is a vector of population PD parameters. For example, \(g\) could be a linear concentration-response function or a nonlinear Emax function. Stage two estimates \(\psi\) using the observed PD responses and their corresponding predicted concentrations.
The individual PK fits and the PD fit rely on assumptions. The PK model must adequately describe each subject’s concentration-time profile, and each subject’s PK data must contain enough information to estimate the parameters being fitted. The chosen PD model must represent the response-concentration relationship adequately. If repeated PD observations from one subject have correlated errors, a fitting method that treats all those errors as independent may give inappropriate uncertainty estimates; the PD model or its error structure must account for that dependence.
Estimate \(\hat{\theta}_i\) using only subject \(i\)’s PK observations and the specified PK model.
Evaluate the subject-specific fitted PK model at the times of that subject’s PD measurements.
Relate the observed PD responses to their corresponding predicted concentrations and estimate the common parameters \(\psi\).
Worked PK Calculation
Worked Example: Predicting Concentration From an Individual Fit
Suppose hypothetical subject B has a measured concentration of \(10\) mg/L at time zero and \(5\) mg/L at hour 4. Assume a one-compartment model after an intravenous bolus, with concentration \(C(t)=C_0\exp(-kt)\). For this calculation, treat the two concentrations as exact model inputs to estimate the elimination rate and predict concentration at hour 2.
From \(C(4)=C_0\exp(-4k)\), the rate estimate is
The prediction at hour 2 must use the elapsed time of 2 hours:
Thus, the concentration supplied to the PD fit for a response measured at hour 2 would be approximately \(7.071\) mg/L. The units in the exponential cancel because the rate has units of inverse hours and time has units of hours.
Estimating a Population PD Relationship
To show the second stage explicitly, consider a hypothetical linear PD model with a common intercept and slope:
Here, \(\beta_0\) is the expected response at predicted concentration zero, and \(\beta_1\) is the expected change in response per 1 mg/L higher predicted concentration. These are population-level parameters in the fitted stage-two model; they are not individual PK parameters. For independent, constant-variance normal errors, the stage-two maximum-likelihood estimates equal the ordinary least-squares estimates.
Worked Example: Fitting the Stage-Two PD Model
Suppose a hypothetical analysis has one PD observation from each of four subjects. Their stage-one PK fits produce the predicted concentrations below at the response measurement times. Assume each response is an independent observation with constant-variance errors and that a linear mean response is appropriate for these concentrations.
| Subject | Predicted concentration, mg/L | Observed response, arbitrary units |
|---|---|---|
| A | 2 | 14 |
| B | 4 | 18 |
| C | 6 | 19 |
| D | 8 | 27 |
Let \(x_i=\hat{C}_i\) and \(y_i=Y^{PD}_i\). The sample means are \(\bar{x}=5\) mg/L and \(\bar{y}=19.5\) response units. The least-squares slope is
The intercept is \(\hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x}=19.5-(2.0)(5)=9.5\) response units. The fitted mean response is therefore \(9.5+2.0\hat{C}\). For example, at \(6\) mg/L, the fitted mean is \(9.5+(2.0)(6)=21.5\) response units. This describes the fitted association in this hypothetical data set; it does not establish that increasing concentration causes a response increase.
The residuals, observed minus fitted, are \(0.5,\ 0.5,\ -2.5,\ 1.5\) response units. Their squared sum is \(0.25+0.25+6.25+2.25=9.00\) response-units squared. With two fitted coefficients and four independent observations, the residual degrees of freedom are \(4-2=2\), giving a residual variance estimate of \(9.00/2=4.50\) response-units squared. With only four observations, this calculation does not establish that the error assumptions are satisfactory; it simply shows the fitting arithmetic.
Why the Least-Squares Slope Has This Form
For the linear stage-two model, least squares chooses \(\beta_0\) and \(\beta_1\) to minimize the sum of squared residuals:
At a minimum, the partial derivatives with respect to both coefficients are zero. Setting the derivative with respect to the intercept to zero gives \(\sum_i(y_i-\hat{\beta}_0-\hat{\beta}_1x_i)=0\), so \(\hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x}\). Substituting this relation into the derivative with respect to the slope yields
Provided the predicted concentrations are not all identical, the denominator is positive, and solving gives the slope formula used in the example. If all \(x_i\) values are identical, the data contain no concentration variation with which to estimate a linear slope separately from the intercept.
Individual Fits Versus a Common PK Fit
Worked Example: Why Subject-Specific Predictions Matter
Consider two hypothetical subjects, each with concentration \(10\) mg/L at time zero, but with different fitted elimination rates: \(k_A=0.10\ \text{h}^{-1}\) and \(k_B=0.30\ \text{h}^{-1}\). Under the one-compartment decay model, their hour-4 predictions are
If, instead, a common rate of \(0.20\ \text{h}^{-1}\) were imposed for illustration, both subjects would be assigned \(10\exp[-(0.20)(4)]=4.493\) mg/L. That common prediction is lower than subject A’s individual prediction by about \(2.210\) mg/L and higher than subject B’s by about \(1.481\) mg/L. This example illustrates why fitting PK separately can preserve observed subject differences when the data support individual parameter estimates. It does not show that every subject’s parameters can be estimated reliably from sparse data.
Diagnostics, Alternatives, and Limitations
Diagnostics should be considered at both stages. For each individual PK fit, compare observed concentrations with fitted concentration-time curves and inspect residuals against time and fitted concentration. Systematic curvature, time-dependent residual patterns, or a single observation dominating the fit can indicate a poor model, an unsuitable error assumption, or influential data. Check whether the fitted parameters are plausible and whether the subject has enough informative PK observations to estimate them. A numerical optimizer returning parameter values does not by itself demonstrate that those values are well determined.
For stage two, inspect the response-versus-predicted-concentration relationship and residuals against predicted concentration, fitted response, time, and subject where relevant. Curvature may indicate that a linear PD model is inadequate; changing residual spread may conflict with constant variance; and subject-patterned residuals may indicate dependence or subject-specific response differences not represented in the model. A nonlinear PD model, such as an Emax relationship, may be more appropriate when supported by the design and data. In all cases, the model should be checked against the actual observation structure rather than selected from a plot alone.
The two-stage approach is distinct from fitting one common PK model across subjects and then fitting PD: individual fits allow subject-specific PK parameters, while a common PK fit constrains them. It is also distinct from simultaneous PK/PD estimation, which estimates PK and PD parameters together using a joint model for both data types. The two-stage procedure can be useful when separate fits are practical or when the analysis deliberately defines the PD fit conditional on individual PK estimates. Its limitations include weakly identified individual PK fits, sensitivity to the chosen PK model, and dependence of the PD fit on the generated concentration values. A later treatment considers in detail how uncertainty in those fitted concentrations affects inference; here, the essential point is to state clearly that stage two uses the stage-one predictions as its inputs.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- In the two-stage procedure, which data are used to estimate subject \(i\)’s PK parameters, and what is passed to the stage-two PD fit?
- A subject has \(C_0=8\) mg/L and \(C(3)=4\) mg/L under exponential decline. Calculate \(k\) and the predicted concentration at hour 1.5.
- In a linear PD model, predicted concentrations are all identical across observations. Can the slope and intercept both be estimated separately? Explain.
- Name one diagnostic for the individual PK fits and one diagnostic for the stage-two PD fit.
- How does two-stage individual fitting differ from simultaneous PK/PD estimation?
Solutions. (1) The PK parameters are estimated from subject \(i\)’s PK observations; the fitted model then supplies predicted concentrations at that subject’s PD observation times. (2) \(k=\ln(8/4)/3=\ln(2)/3=0.2310\ \text{h}^{-1}\). Thus \(C(1.5)=8\exp[-(0.2310)(1.5)]=8\exp(-0.3465)\approx5.657\) mg/L. (3) No. With no concentration variation, the slope is not separately identified from the intercept. (4) For PK, inspect concentration residuals against time; for PD, inspect response residuals against predicted concentration for patterns such as curvature or changing spread. (5) In two-stage individual fitting, PK parameters are first estimated separately by subject and predictions are then used in a separate PD fit. Simultaneous estimation fits a joint model to PK and PD data together.