Sparse-Sampling Population Pharmacokinetics

Why Sparse Sampling Changes the Modeling Task

In a richly sampled pharmacokinetic profile, each participant may contribute enough measurements to show both the rapid distribution phase and the slower elimination phase. In a sparse trial, each participant contributes only a few concentrations, often at different scheduled times. The analysis must therefore distinguish what the population as a whole can inform from what the data say about any one participant.

A two-compartment model can be useful when concentration falls through an early distribution phase followed by a slower terminal phase. Sparse sampling does not make the model impossible to estimate, but it makes the sampling schedule, population assumptions, and predictive evaluation especially important. A good population fit does not imply that every participant's clearance, volume, or individual concentration curve has been precisely determined.

Definition: A sparse-sampling population PK analysis estimates typical pharmacokinetic parameters and between-participant variability from concentration data in which individual participants have relatively few observations. Population parameters describe the distribution of PK behavior; individual estimates are informed by both the participant's data and that population distribution.

Two-Compartment Model and Assumptions

For an IV bolus dose, let \(A_c(t)\) and \(A_p(t)\) denote drug amounts in the central and peripheral compartments. The central concentration is \(C(t)=A_c(t)/V_c\), where \(V_c\) is central volume. Let \(CL\) be clearance from the central compartment, \(Q\) intercompartmental clearance, and \(V_p\) peripheral volume. A linear model with first-order transfer and elimination is:

\frac{dA_c}{dt}=-(k_{10}+k_{12})A_c+k_{21}A_p,\qquad \frac{dA_p}{dt}=k_{12}A_c-k_{21}A_p

Here \(k_{10}=CL/V_c\), \(k_{12}=Q/V_c\), and \(k_{21}=Q/V_p\), each with units of inverse time. For a dose \(D\) placed initially in the central compartment, \(A_c(0)=D\) and \(A_p(0)=0\). The model assumes linear, time-invariant disposition over the modeled period; accurate dose and sample times; positive clearances and volumes; and no unmodeled input process such as absorption. If the administered drug is oral, the input model must also be specified, so the IV-bolus equations below should not be applied as though oral dosing were instantaneous.

The central concentration has a biexponential form. The two positive decay rates \(\alpha\) and \(\beta\), with \(\alpha>\beta\), are roots of the characteristic equation. With \(k_{10}, k_{12}, k_{21}>0\), the roots satisfy:

r^2-(k_{10}+k_{12}+k_{21})r+k_{10}k_{21}=0,\qquad \alpha,\beta=\frac{(k_{10}+k_{12}+k_{21})\pm \sqrt{(k_{10}+k_{12}+k_{21})^2-4k_{10}k_{21}}}{2}

The plus sign gives \(\alpha\), the faster rate; the minus sign gives \(\beta\), the slower rate. The concentration can then be written \(C(t)=A e^{-\alpha t}+B e^{-\beta t}\), where \(A\) and \(B\) are concentration coefficients, not amounts.

A=\frac{D}{V_c}\frac{\alpha-k_{21}}{\alpha-\beta},\qquad B=\frac{D}{V_c}\frac{k_{21}-\beta}{\alpha-\beta}

For a population model, a common assumption is that positive individual parameters vary lognormally: for example, \(CL_i=CL_{\mathrm{typ}}\exp(\eta_{CL,i})\), with \(\eta_{CL,i}\) drawn from a mean-zero normal distribution. A lognormal residual model can be written \(\log C_{ij}^{\mathrm{obs}}=\log C_{ij}^{\mathrm{pred}}+\epsilon_{ij}\), where the residuals are mean-zero normal on the log scale and conditionally independent given the individual parameters. These are modeling assumptions to assess, not properties guaranteed by sparse data.

Worked Example: Deriving a Hypothetical Two-Compartment Curve

Suppose a hypothetical IV-bolus study gives a dose of \(100\ \mathrm{mg}\). Consider parameters \(V_c=20\ \mathrm{L}\), \(CL=5\ \mathrm{L/h}\), \(Q=8\ \mathrm{L/h}\), and \(V_p=40\ \mathrm{L}\). These are illustrative model values, not results from a real study. Find the micro-rate constants, the biexponential rates, and several predicted concentrations.

First, \(D/V_c=100/20=5\ \mathrm{mg/L}\), and the micro-rate constants are \(k_{10}=5/20=0.25\ \mathrm{h}^{-1}\), \(k_{12}=8/20=0.40\ \mathrm{h}^{-1}\), and \(k_{21}=8/40=0.20\ \mathrm{h}^{-1}\). Thus the sum in the characteristic equation is \(0.85\ \mathrm{h}^{-1}\), and the product is \(k_{10}k_{21}=0.0500\ \mathrm{h}^{-2}\).

The discriminant is \(\sqrt{0.85^2-4(0.0500)}=\sqrt{0.5225}=0.7228416\ \mathrm{h}^{-1}\). Therefore \(\alpha=(0.85+0.7228416)/2=0.7864208\ \mathrm{h}^{-1}\) and \(\beta=(0.85-0.7228416)/2=0.0635792\ \mathrm{h}^{-1}\). Substituting into the coefficient formulas gives \(A\approx4.0564\ \mathrm{mg/L}\) and \(B\approx0.9436\ \mathrm{mg/L}\); their sum is \(5.000\ \mathrm{mg/L}\), as required by \(C(0)=D/V_c\).

For example, at \(2\ \mathrm{h}\), \(C(2)=4.0564e^{-0.7864208(2)}+0.9436e^{-0.0635792(2)}\approx0.8415+0.8310=1.6725\ \mathrm{mg/L}\). At \(6\ \mathrm{h}\), the same calculation gives approximately \(0.6806\ \mathrm{mg/L}\); at \(12\ \mathrm{h}\), it gives approximately \(0.4403\ \mathrm{mg/L}\). The fast component has largely decayed by these later times, while the terminal component remains measurable.

What Sparse Data Can Identify

Identifiability has two related meanings. Structural identifiability asks whether distinct parameter values could produce exactly the same ideal, continuous concentration-time curve. Practical identifiability asks whether the actual finite, noisy, scheduled samples contain enough information to estimate parameters with useful precision. A model may be structurally identifiable yet practically weakly identified in a sparse study.

In the hypothetical curve above, early samples help characterize the rapid distribution decline, while sufficiently late samples help characterize the terminal slope. Samples concentrated only near the terminal phase may provide little direct information about the early distribution rate or its amplitude. Conversely, sampling only early can leave the terminal rate poorly determined. Even when the population distribution of sample times spans both phases, an individual with two concentrations generally cannot independently determine all four quantities \(CL\), \(V_c\), \(Q\), and \(V_p\). The population model pools information across participants under its assumptions.

Key distinction: A parameter can be estimable as a population parameter even when it is not precisely estimable for each participant. Individual empirical Bayes estimates are conditional estimates, not direct measurements of that participant's true parameters.

The design should therefore be evaluated against the parameters and summaries of interest. Check whether sample times cover the expected distribution and terminal phases, whether enough participants contribute observations at each informative time region, and whether sample timing is recorded accurately. If a phase is not sampled, adding complexity to the model does not create information about it. Sensitivity to starting values, wide uncertainty intervals, strong parameter correlations, or unstable estimates across reasonable model specifications can signal weak practical identifiability.

Worked Example: Planning and Reading a Sparse Schedule

Suppose 24 hypothetical participants each provide two samples after the same \(100\ \mathrm{mg}\) IV bolus: one sample at \(0.5\ \mathrm{h}\) and one at \(12\ \mathrm{h}\). Under the illustrative curve above, the typical-subject predictions at those times are approximately \(3.653\ \mathrm{mg/L}\) and \(0.4403\ \mathrm{mg/L}\), respectively.

The early observation is positioned to inform the distribution-rich part of the curve, and the late observation informs the terminal part. However, the schedule has only two observations per participant for four PK parameters. It cannot, by itself, provide four independent concentration constraints for each individual. Population estimates can draw on the repeated early-time and late-time observations across the 24 participants, but the reliability of each parameter still depends on the model, variability, residual error, and sample-time coverage.

A more informative sparse design could stagger additional samples across intermediate times, so that the population data cover more of the curve. This does not mean every participant must be sampled at every time. It does mean that conclusions about a parameter should be checked against where the data actually lie, rather than inferred solely from a successful optimizer run.

Estimation, Shrinkage, and Predictive Checks

Population estimation separates typical parameters and between-participant variability from residual measurement variability. In a nonlinear mixed-effects analysis, the likelihood accounts for the fact that individual random effects are not directly observed. Empirical Bayes estimates, sometimes called individual conditional estimates, combine an individual's sparse data with the estimated population distribution. When those data are weak, the estimates tend to move toward the population mean.

One commonly reported eta-shrinkage convention is \(1-\mathrm{SD}(\widehat{\eta})/\omega\), where \(\widehat{\eta}\) is the set of individual empirical Bayes estimates for a random effect and \(\omega\) is its estimated population standard deviation. This convention must be stated because variance-based shrinkage, \(1-\mathrm{Var}(\widehat{\eta})/\omega^2\), is a different calculation and can give a different percentage. High shrinkage means individual estimates contain limited information relative to the population variability; it does not by itself prove that the population model is wrong.

Shrinkage also affects diagnostics. If many empirical Bayes estimates are pulled toward zero, plots of those estimates against covariates or other variables can hide real patterns or show misleading patterns. Residual plots based on individual predictions can also look artificially reassuring. In sparse settings, prioritize checks that compare observed data with predictions from simulated populations using the same study design and analysis assumptions.

Reproduce the design.
Simulate participants, doses, sample times, and missingness to match the trial. Use the same residual and censoring assumptions where appropriate.

Simulate concentration data.
Draw individual parameters and residual errors from the fitted model, then calculate concentrations at the observed sampling times.

Compare distributions.
Compare observed concentration quantiles with corresponding simulated quantiles across time bins or scheduled sampling times. Investigate systematic departures, not isolated random differences alone.

Check the intended use.
Compare relevant exposure summaries, such as AUC or peak concentration, as well as concentration-time behavior. A fit to observed concentrations does not automatically validate every extrapolated summary.

Worked Example: Calculating Shrinkage and Checking Predictions

In a hypothetical analysis, the estimated population standard deviation for a clearance random effect is \(\omega=0.30\) on the log scale. The standard deviation of its empirical Bayes estimates across participants is \(0.12\). Using the stated SD-based convention, eta-shrinkage is:

1-\frac{\mathrm{SD}(\widehat{\eta}_{CL})}{\omega} =1-\frac{0.12}{0.30}=0.60

Thus the reported shrinkage is \(60\%\) under this convention. It indicates that the empirical Bayes estimates vary much less than the estimated underlying random effects; it does not mean that clearance estimates are “60% wrong.” If the variance-based convention were used instead, the calculation would be \(1-(0.12^2/0.30^2)=0.84\), or \(84\%\). Reporting the definition prevents these distinct metrics from being confused.

Now suppose a design-matched predictive check gives hypothetical observed median concentrations of \(3.8\), \(1.6\), \(0.73\), and \(0.39\ \mathrm{mg/L}\) at \(0.5\), \(2\), \(6\), and \(12\ \mathrm{h}\). The typical-subject reference predictions are approximately \(3.653\), \(1.672\), \(0.6806\), and \(0.4403\ \mathrm{mg/L}\), respectively. These small differences alone do not establish failure; uncertainty and simulated quantile bands are needed. If, for example, observed late-time quantiles repeatedly fall outside the simulated bands, investigate terminal-rate assumptions, sampling or assay issues, and residual variability before relying on terminal exposure predictions.

Exposure Summaries and Alternative Approaches

For a linear IV-bolus model with elimination from the central compartment, total exposure is \(AUC_{0-\infty}=D/CL\). This follows because total drug eliminated is the dose and the elimination rate at time \(t\) is \(CL\,C(t)\); integrating that rate over time gives \(D=CL\int_0^\infty C(t)\,dt\). Equivalently, integrating the biexponential curve gives \(A/\alpha+B/\beta\), which must agree with \(D/CL\) when the model and units are consistent.

Worked Example: Computing AUC From the Population Model

For the hypothetical IV dose and clearance above, \(D=100\ \mathrm{mg}\) and \(CL=5\ \mathrm{L/h}\). Therefore:

AUC_{0-\infty}=\frac{D}{CL} =\frac{100\ \mathrm{mg}}{5\ \mathrm{L/h}} =20\ \mathrm{mg\cdot h/L}

This is the model-based total AUC for the stated dose and clearance. Integrating the concentration expression gives \(A/\alpha+B/\beta\), approximately \(20\ \mathrm{mg\cdot h/L}\), providing a useful conservation check. In a population analysis, individual or population AUC summaries should reflect the relevant parameter distribution; the AUC computed from typical clearance is not automatically the population mean AUC when clearance is variable.

A one-compartment model is a simpler alternative when the data do not support a distinct distribution phase; it should not be rejected merely because a more complex model can fit the sampled points. Noncompartmental analysis can summarize observed exposure with fewer structural assumptions, but sparse samples may make the area between observations and the terminal extrapolation uncertain. A population compartment model can use the full sampling design and quantify population variability, but its exposure summaries depend on the structural and statistical assumptions. The choice should follow the study question and available information.

Important limitations include sparse information about individual parameters, sensitivity to misspecified residual or random-effects distributions, and uncertainty in extrapolating beyond observed times. Concentrations below the assay's quantification limit should be handled with an appropriate censored-data method rather than automatically replaced by zero or by half the limit. A model can also fail because recorded dose times, sample times, or units are wrong; those operational checks are as important as examining parameter estimates.

Common Pitfalls and Practical Review

Claiming individual precision from a population fit. A stable population estimate does not imply that each participant's four PK parameters are individually identified. Examine information by parameter and acknowledge shrinkage.
Treating high shrinkage as evidence of model failure. Sparse data naturally produce shrinkage. Interpret its magnitude alongside the design, uncertainty, and predictive performance.
Using only empirical Bayes diagnostic plots. High shrinkage can mask patterns in individual estimates and residuals. Include design-matched population predictive checks.
Assuming a terminal-looking sample proves the terminal slope. A few late observations may be noisy or affected by quantification limits. Assess the consistency of late-time data and the uncertainty of the terminal rate.
Reporting exposure without its assumptions. \(D/CL\) applies to this linear IV setting. State dose, units, and whether the result is typical-subject, individual, or population-level.

Key takeaway: Sparse-sampling population PK relies on pooling across participants, but pooling cannot replace missing information without assumptions. Match the design to the phases and summaries of interest, state the parameterization and shrinkage convention, and test predictions against data simulated under the actual sampling scheme.

Check Your Understanding

Answer each question, then compare with the solution given directly below it.

A hypothetical IV-bolus model has \(D=60\ \mathrm{mg}\) and \(CL=3\ \mathrm{L/h}\). What is \(AUC_{0-\infty}\), and what assumptions support this calculation?
Solution: \(AUC_{0-\infty}=D/CL=60/3=20\ \mathrm{mg\cdot h/L}\). This uses linear disposition, IV administration, and clearance that is constant over the modeled period.
Why might a population distribution rate be estimable when most participants have only two samples, while their individual distribution rates remain uncertain?
Solution: Participants collectively contribute observations across sampling times, allowing population estimation under the shared model. Two concentrations per person do not independently determine all four two-compartment PK parameters, so individual estimates also rely strongly on the population distribution.
Using the SD-based convention, calculate eta-shrinkage when \(\omega=0.25\) and \(\mathrm{SD}(\widehat{\eta})=0.10\).
Solution: \(1-0.10/0.25=0.60\), so shrinkage is \(60\%\). This is the SD-based metric; the variance-based metric would be \(1-0.10^2/0.25^2=84\%\).
What is one reason to simulate the observed sampling schedule during a predictive check?
Solution: The simulated and observed distributions should be compared under the same design because the number and timing of samples determine which phases and concentration quantiles are represented. A check at a different schedule may not reveal problems relevant to the trial data.