Why Longitudinal PK/PD Data Need a Joint Fit
Preclinical PK/PD experiments often measure drug concentration and pharmacological response repeatedly after a dose. These time courses carry complementary information: concentration observations constrain drug disposition, while response observations constrain how exposure drives the biological system. Fitting them together can therefore inform a mechanistic model more strongly than fitting either series alone.
The measurements are not automatically independent simply because they appear as separate rows in a data file. Multiple samples from one animal may share biological variation, and different assays may have different error scales. In addition, concentration and response measurements are often collected at different times. A sound fit represents these features explicitly rather than treating all rows as equally precise, independent observations on a common schedule.
Model and Assumptions
Let \(\theta\) denote the biological parameters of the mechanistic model. For an illustrative one-compartment PK model after an intravenous bolus dose, concentration \(C(t)\) follows
where \(D\) is dose, \(V\) is apparent volume, and \(CL\) is clearance. A possible PD component is an indirect-response model in which drug inhibits production of a response \(R(t)\):
This is an example structure, not a requirement for every application. The chosen PK/PD equations should represent the scientific system and the measured endpoints. The observation model connects the latent states to data: a concentration assay measures \(C(t)\), and a response assay measures \(R(t)\), possibly after a defined scale transformation or assay-specific offset.
For animal \(a\), assay \(k\), and observation \(i\), write the observed value as \(y_{aki}\), measured at time \(t_{aki}\). Let \(h_k(t_{aki};\theta)\) be the model prediction for that assay at that time. The calibration requires the following assumptions, or explicit alternatives to them:
- The model equations, dose history, time origin, and units are correctly specified for the data being fitted.
- The observation model gives a plausible distribution and scale for each assay’s measurement error. For example, additive normal error may suit an assay with roughly constant absolute error; lognormal error may suit positive measurements with roughly constant relative error.
- Any dependence among repeated measurements from an animal is represented in the likelihood or covariance structure. Measurements from different animals may be treated as independent only when that is defensible.
- Model predictions are evaluated at the recorded sampling times. An observation need not exist at every time point or for every assay.
These are modeling assumptions, not facts guaranteed by collecting time-course data. Earlier material on assay-specific observation error and within-study correlation provides the underlying concepts; here they enter directly into fitting the longitudinal model.
Constructing the Calibration Objective
If errors are independent and normally distributed, with known or specified standard deviation \(s_{aki}\), then \(y_{aki}\) has mean \(h_k(t_{aki};\theta)\) and variance \(s_{aki}^{2}\). Maximizing the likelihood is equivalent to minimizing the weighted sum of squared residuals
A residual is observed minus predicted. Dividing by the assay standard deviation expresses its size relative to the expected measurement error. Thus a residual of 2 concentration units should not count the same as a residual of 2 response units unless their error scales make those deviations comparably surprising. If error variances are estimated as part of the calibration, the likelihood must also include the corresponding variance terms; simply minimizing \(Q\) while allowing every error scale to grow would not define a valid fit.
With repeated measurements that are correlated, collect the residuals for an animal into a vector \(e_a\), and let \(\Sigma_a\) be their covariance matrix. Under a multivariate normal error model, the negative log-likelihood, up to a constant, is
The covariance matrix can represent different assay variances as well as correlation over time. If \(\Sigma_a\) is fixed, the determinant term does not affect which \(\theta\) minimizes the objective, but it remains part of the likelihood. If covariance parameters are estimated, the determinant term is essential: it prevents the fit from improving merely by inflating the variances.
Worked Fits and Calculations
Worked Example: Estimating a Terminal PK Rate
Suppose a hypothetical animal receives an intravenous bolus, and the concentration is measured as \(10.0\) mg/L at \(t=0\) h and \(2.50\) mg/L at \(t=4\) h. Assume the one-compartment model above and, for this calculation, treat the two measurements as exact. Estimate the elimination rate constant \(k=CL/V\) from these two points.
The model solution is \(C(t)=C(0)e^{-kt}\). Therefore,
The implied half-life is \(\log(2)/k=0.6931/0.3466=2.000\) h. The two observations are consistent with a fourfold concentration decline over two half-lives. This calculation identifies the rate constant under the assumed model; it does not separate \(CL\) from \(V\), and it does not account for measurement error. In an actual calibration, all available concentration observations would contribute through their observation model.
Worked Example: Comparing Joint PK/PD Parameter Candidates
Consider hypothetical concentration observations of \(8.0\), \(4.0\), and \(2.0\) mg/L at 0, 2, and 4 h, each with assay standard deviation \(0.8\) mg/L. Response observations are \(100\) and \(80\) normalized response units at 0 and 4 h, each with standard deviation \(5\) units. The unequal schedules are intentional: there is no response measurement at 2 h. Two candidate parameter sets give the following predictions:
Candidate A predicts concentrations \(7.2, 4.5, 2.8\) mg/L and responses \(100, 87\) units. Candidate B predicts concentrations \(8.4, 3.7, 1.8\) mg/L and responses \(100, 78\) units. Under independent normal errors with the stated standard deviations, calculate \(Q\) for each candidate using only the observations that exist.
For Candidate A, the concentration contribution is
Its response contribution is \(((100-100)/5)^2+((80-87)/5)^2=0+1.960=1.960\). Thus \(Q_A=2.391+1.960=4.351\).
For Candidate B, the concentration contribution is
Its response contribution is \(((100-100)/5)^2+((80-78)/5)^2=0+0.1600=0.1600\), so \(Q_B=0.4531+0.1600=0.6131\). Among these two candidates, B has the smaller weighted residual objective. This comparison does not establish that B is the global optimum or that the model is adequate; it only compares the specified candidates under the stated error assumptions. The calculation also shows how unequal schedules are handled: each available observation is evaluated at its own time, with no imputation at missing times.
Worked Example: Accounting for Correlated Repeated Responses
Suppose two response measurements from the same hypothetical animal have residuals \(e=(1,1)^T\) response units. Each has marginal standard deviation \(1\) unit, and their within-animal correlation is assumed to be \(0.5\). Then
The correlated quadratic residual contribution is
If the measurements were incorrectly treated as independent with the same marginal variances, the contribution would be \(1^2+1^2=2.000\). The positive correlation makes two same-direction residuals less informative as two independent pieces of evidence. The covariance determinant is \(|\Sigma|=0.75\), so the negative log-likelihood contribution, including the determinant and normalizing constant for these two observations, is
The independent-error value would be \(\frac{1}{2}[2+2\log(2\pi)]=2.838\). These likelihood values are meaningful only relative to the specified covariance models; correlation should be justified from the design or data, not added solely because it lowers an objective.
Practical Fitting, Diagnostics, and Alternatives
A typical calibration evaluates the PK/PD equations numerically at the actual times and doses, computes predicted measurements through the observation model, and optimizes the joint likelihood over \(\theta\) and any justified error parameters. A useful workflow is:
Confirm dose, units, time origin, animal identifiers, assay identity, and which samples belong to each animal.
Write the PK/PD equations and define how each assay measures or transforms its corresponding model state.
Choose assay-specific error distributions and represent repeated-measure dependence where supported by the design or data.
Minimize the resulting negative log-likelihood, then inspect predictions and residuals by assay, time, and animal.
Diagnostics should be performed on the scale used by the observation model. Plot observed values and predictions against time for each assay, and examine standardized residuals, such as \(e_{aki}/s_{aki}\), against time and predicted value. Persistent early-versus-late patterns may indicate a mismatch in absorption, distribution, turnover, or delay structure. Residuals that grow with predicted concentration can suggest that an additive constant-variance model is unsuitable. Residuals clustered by animal can signal unmodeled animal-level variability or an inadequate correlation structure. Numerical integration failures, estimates at parameter bounds, and strong sensitivity to starting values also warrant investigation.
A common alternative is to fit PK first and then fit PD using estimated exposure or a PK prediction as input. This staged approach can be easier to initialize and may be reasonable when PK is well determined and uncertainty propagation is handled. However, treating the fitted PK parameters as known can understate uncertainty in PD parameters. Simultaneous fitting propagates the shared dependence through the joint objective, but can be harder to optimize and may expose weakly identified parameter combinations. Independent weighted least squares is another option when errors are independent and variances are fixed; it is not a substitute for a correlated likelihood when observations from an animal share errors.
The joint fit has limits. A precise optimizer cannot compensate for an incorrect structural model, unrecorded dose timing, poor assay calibration, or too few informative sampling times. A large number of rows from a small number of animals does not necessarily provide the information of the same number of independent animals. Moreover, a good fit to preclinical time courses describes those data under the model assumptions; by itself it does not establish that the same parameters or response relationship apply in another species or setting.
Common Mistake
Check Your Understanding
Answer each question, then compare with the solutions.
- Why should a concentration measured at 3 h be compared with the model prediction at 3 h, even if most response samples were collected at 4 h?
- Two assays have residuals of 2 units, but their standard deviations are 1 unit and 4 units. What are their squared standardized residual contributions?
- For two repeated measurements with residuals \(1,1\), unit marginal variances, and correlation \(0.5\), what is the correlated quadratic residual contribution?
- Name one diagnostic pattern that could indicate that a constant-variance additive error model is inadequate.
Solutions. (1) Each observation contributes information about the model prediction at its recorded time; using another time would compare different quantities. (2) The contributions are \((2/1)^2=4\) and \((2/4)^2=0.25\), respectively, so equal raw residuals need not have equal weight. (3) The inverse covariance calculation gives \(4/3=1.333\), rather than the independent value \(2\). (4) For example, residual spread that systematically increases with predicted concentration suggests that a constant absolute error variance may be inappropriate.