Joint Updating From Individual PK and PD Observations

Why Update PK and PD Together?

An individual’s concentration data inform how the drug is handled, while pharmacodynamic (PD) observations inform the resulting response. When the PK and PD models share parameters or depend on the same latent concentration trajectory, the two data streams can inform one another. A response observation may help distinguish PK-PD parameter combinations that concentration samples alone cannot, and PK data can constrain the exposure driving the response.

The central requirement is to combine the observations through a probability model that represents their dependencies. Multiplying a PK likelihood and a PD likelihood is appropriate only when they are conditionally independent given the model’s parameters and latent states. It is not generally appropriate to assume that the observations are independent without conditioning: both can be related through the individual’s underlying concentration.

Definition: Let $\theta$ denote the individual parameter vector, including the relevant PK and PD parameters, and let $C(t;\theta)$ be the modelled concentration trajectory. A joint Bayesian update uses the prior distribution $p(\theta)$ and the joint likelihood of the individual’s PK observations $Y$ and PD observations $R$: $$ p(\theta\mid Y,R) \propto p(\theta)\,p(Y,R\mid\theta). $$ If PK and PD observations are conditionally independent given $\theta$, then $p(Y,R\mid\theta)=p(Y\mid\theta)p(R\mid\theta)$. The posterior is normalized so that its density integrates to one.

A Joint PK-PD Model

For a basic model, suppose the PK observations $Y_j$ are measured at times $t_j$, and PD responses $R_k$ are measured at times $s_k$. Let $C(t;\theta)$ be the individual’s model-predicted concentration, and let $g(C;\theta)$ be the PD model’s expected response at concentration $C$. One possible observation model is

Y_j\mid\theta \sim N\!\left(C(t_j;\theta),\sigma_C^2\right), \qquad R_k\mid\theta \sim N\!\left(g(C(s_k;\theta);\theta),\sigma_R^2\right).

Here $\sigma_C$ and $\sigma_R$ are residual standard deviations, expressed in the units of concentration and response, respectively. For example, $g$ might be an increasing Emax relationship, $E_0+E_{\max}C/(EC_{50}+C)$. Other response models are possible; the likelihood must match the chosen response scale and clinical mechanism.

Under conditionally independent residual errors, the joint likelihood is the product of the observation densities. If there are $J$ concentration observations and $K$ responses, the normal-error version is

p(Y,R\mid\theta) = \prod_{j=1}^{J} \frac{1}{\sqrt{2\pi\sigma_C^2}} \exp\!\left[-\frac{\{Y_j-C(t_j;\theta)\}^2}{2\sigma_C^2}\right] \prod_{k=1}^{K} \frac{1}{\sqrt{2\pi\sigma_R^2}} \exp\!\left[-\frac{\{R_k-g(C(s_k;\theta);\theta)\}^2}{2\sigma_R^2}\right].

This factorization assumes that, conditional on the individual parameters and modelled concentration trajectory, the residual errors have the specified distributions and are independent across the observations and between PK and PD. If errors are correlated, or the PD model depends on an unobserved concentration that is not determined by $\theta$, the joint model must represent that dependence rather than use this factorization unchanged. A proper prior yields a proper posterior when the marginal likelihood is finite and positive; properness of the prior alone is not sufficient, and the likelihood and prior should assign positive support to plausible parameter values.

From Joint Likelihood to Posterior Forecast

The update follows directly from Bayes’ rule. Conditional independence permits the likelihood to be written as a PK contribution times a PD contribution; the same prior is then reweighted by both. The derivation is short but important: each observation changes the relative support for parameter values according to how well those values predict that observation.

p(\theta\mid Y,R) = \frac{p(\theta)p(Y,R\mid\theta)} {\int p(\vartheta)p(Y,R\mid\vartheta)\,d\vartheta} = \frac{p(\theta)p(Y\mid\theta)p(R\mid\theta)} {\int p(\vartheta)p(Y\mid\vartheta)p(R\mid\vartheta)\,d\vartheta}.

The integration variable $\vartheta$ is a dummy parameter vector. The denominator ensures that the posterior is normalized. In particle or grid calculations, the same result is implemented by multiplying each prior weight by the PK and PD likelihoods, then normalizing the resulting weights.

For a future exposure summary $A(\theta)$, such as AUC, and future response $R^{new}$, forecasts average over the updated parameter uncertainty. A response forecast also includes future observation variability if the target is an observed response rather than the model’s expected response.

p(R^{new}\mid Y,R) = \int p(R^{new}\mid\theta)\,p(\theta\mid Y,R)\,d\theta, \qquad E[A(\theta)\mid Y,R] = \int A(\theta)\,p(\theta\mid Y,R)\,d\theta.

Thus the forecast is not generally the response evaluated at a single “best” parameter estimate. A posterior predictive distribution carries uncertainty about the individual parameters forward into both exposure and response.

Worked Example: Conflicting PK and PD Evidence

Suppose two candidate parameter states, A and B, have equal prior probability. At a sampled time, their predicted concentrations are 4 and 6 mg/L, respectively, and the PK residual standard deviation is 1 mg/L. The measured concentration is 4.5 mg/L. At a response time, the candidate PD means are 10 and 14 response units, the response residual standard deviation is 2 units, and the observed response is 13.5 units. Assume conditional independence of the PK and PD residuals given the candidate state.

The normal density constants cancel when calculating likelihood ratios. For the PK observation, the ratio of the likelihood under A to that under B is

\frac{L_C(A)}{L_C(B)} = \exp\!\left[-\frac{(4.5-4)^2-(4.5-6)^2}{2(1)^2}\right] = \exp(1) \approx 2.718.

The concentration therefore favors A. For the response observation, the corresponding ratio is

\frac{L_R(A)}{L_R(B)} = \exp\!\left[-\frac{(13.5-10)^2-(13.5-14)^2}{2(2)^2}\right] = \exp(-1.5) \approx 0.2231.

The response favors B. Equal prior odds become posterior odds for A versus B of $2.718\times0.2231\approx0.6065$. Therefore,

P(A\mid Y,R)=\frac{0.6065}{1+0.6065}\approx0.3775, \qquad P(B\mid Y,R)\approx0.6225.

The observations pull in opposite directions, so neither should be treated as decisive alone. Suppose the future AUC predictions are 30 and 50 mg·h/L and the future mean responses are 9 and 15 units under A and B. Their posterior-averaged predictions are $0.3775(30)+0.6225(50)=42.45$ mg·h/L and $0.3775(9)+0.6225(15)=12.74$ response units. These are averages over the two states, not guarantees of the individual’s future outcomes.

Calculating Particle Weights and Forecasts

When the posterior cannot be obtained analytically, a finite grid or weighted particle set gives a transparent approximation. For prior particles $\theta^{(m)}$ with weights $w_0^{(m)}$, the unnormalized updated weight is the prior weight times the likelihood of all observations. The normalized weights sum to one. Any forecast that is a function of the particles is then averaged using those updated weights.

Worked Example: Updating Three Candidate Particles

Consider three candidate individual parameter particles with prior weights 0.50, 0.30, and 0.20. For the observed PK data, suppose their likelihood values are 0.20, 0.50, and 0.10; for the PD data, suppose the likelihood values are 0.10, 0.20, and 0.80. These likelihood values are on the same scale across particles for each data stream. Assume the conditional-independence factorization applies.

The unnormalized weights are $0.50(0.20)(0.10)=0.010$, $0.30(0.50)(0.20)=0.030$, and $0.20(0.10)(0.80)=0.016$. Their sum is 0.056, so the normalized posterior weights are

(w^{(1)},w^{(2)},w^{(3)}) = \left(\frac{0.010}{0.056},\frac{0.030}{0.056},\frac{0.016}{0.056}\right) \approx(0.1786,0.5357,0.2857).

The PD data substantially increase the weight of particle 3, while the PK data favor particle 2. Suppose their future concentration predictions are 3, 5, and 8 mg/L. The posterior mean concentration is $0.1786(3)+0.5357(5)+0.2857(8)=5.500$ mg/L. Suppose their future mean response predictions are 20, 35, and 50 units. The posterior mean response is $0.1786(20)+0.5357(35)+0.2857(50)\approx36.61$ units.

If future observed responses have a normal residual standard deviation of 5 units, the predictive variance is the posterior mean residual variance plus the variance between particle-specific response means. The latter is approximately 101.9 squared units; therefore, the predictive variance is $25+101.9=126.9$ squared units, and the predictive standard deviation is approximately 11.26 units. This wider uncertainty reflects both future residual variation and uncertainty about which parameter state best describes the individual.

Measurement Error and the Shared Concentration

A crucial modelling decision is whether a PD response is linked to the true underlying concentration or to a concentration reading treated as exact. If the concentration measurement is noisy, substituting its observed value into the PD model can understate uncertainty and misrepresent the PK-PD dependence. A joint model can instead relate both observations to a latent concentration, with the PK measurement model describing how the observed concentration arises.

Worked Example: Avoiding an Error-Free Concentration Assumption

At a response time, suppose two candidate states predict latent concentrations of 6 and 10 mg/L. A concentration measurement of 8 mg/L has residual standard deviation 2 mg/L. The two candidate states have equal concentration likelihood because each prediction is 2 mg/L from the observed value. Their predicted mean responses at their respective latent concentrations are 40 and 70 units. The observed response is 65 units, with residual standard deviation 10 units. Start with equal prior probabilities and assume conditional independence of measurement and response errors given the state.

The concentration observation leaves the odds unchanged. The response likelihood ratio for the low-concentration state versus the high-concentration state is

\frac{L_R(\text{low})}{L_R(\text{high})} = \exp\!\left[-\frac{(65-40)^2-(65-70)^2}{2(10)^2}\right] = \exp(-3) \approx0.04979.

With equal prior odds, the posterior probability of the low state is $0.04979/(1+0.04979)\approx0.04742$; the high state has probability approximately 0.9526. If future concentrations are predicted to be 5 and 9 mg/L, the posterior mean is $0.04742(5)+0.9526(9)\approx8.810$ mg/L. If the corresponding future mean responses are 35 and 68 units, the posterior mean response is approximately $0.04742(35)+0.9526(68)=66.44$ units.

Here the two states fit the measured concentration equally well, but the response favors the high state. The calculation uses each state’s latent concentration in its PD prediction; it does not treat 8 mg/L as an error-free exposure known to be the true driver of response.

Assumptions, Diagnostics, and Alternatives

A joint update is only as credible as its model. Before interpreting it, check that the prior represents plausible individual parameters, the PK model matches the dosing and sampling context, and the PD link has the correct direction, scale, and time dependence. Verify that residual variances are appropriate and that the assumed conditional independence is defensible. If responses are serially correlated, delayed relative to concentration, or measured with a different error distribution, the likelihood should reflect that structure.

Signs of a problematic update include posterior weights collapsing onto one or a few particles, parameter estimates moving to implausible boundaries, or PK observations favoring parameter combinations that imply implausible PD predictions. Such patterns can arise from sparse data, an overly narrow prior, a misspecified observation model, or incompatible PK and PD assumptions. They are warnings to inspect the assumptions and computation, not proof by themselves that the data are wrong.

PK-only updating uses $p(Y\mid\theta)$ and may be appropriate when PD data are unavailable or not informative for the parameters of interest, but it cannot use response evidence to refine the PK-PD relationship. PD-only updating can be weakly informative about exposure if the response model is flat or its parameters are confounded. A two-stage plug-in approach first estimates concentration or PK parameters and then treats them as fixed in the PD analysis; this is simpler, but generally omits uncertainty and dependence passed from PK to PD. Joint updating preserves that information when the model is correctly specified.

Joint updating does not guarantee precision or eliminate model dependence. A response may be uninformative over a saturated part of an Emax curve, and multiple PK-PD parameter combinations may explain the same observations. Forecasts can remain broad even with a carefully constructed posterior. They also depend on assumptions about future dosing, adherence, and the stability of the individual’s parameters. Use the posterior predictive distribution, not a point estimate alone, when uncertainty matters for exposure or clinical response.

Common pitfall: “The measured concentration is available, so I can insert it into the PD model as the true concentration and then multiply by the PK likelihood.” This treats a noisy measurement as error-free and can misstate uncertainty or count linked evidence incorrectly. Model the latent concentration and its measurement error, or explicitly justify an approximation that treats the concentration as known with negligible error.

Check Your Understanding

Answer each question, then compare with the solutions below.

Under what condition can the joint likelihood be factored into a PK likelihood times a PD likelihood?
A candidate particle has prior weight 0.25, PK likelihood 0.4, and PD likelihood 0.5. What is its unnormalized updated weight?
Why can a PD observation change the posterior for a PK parameter such as clearance?
When might a two-stage plug-in analysis give misleadingly narrow response forecasts?

Solutions. (1) The factorization is valid when the PK and PD observations are conditionally independent given the model parameters and any relevant latent states. (2) Its unnormalized weight is $0.25(0.4)(0.5)=0.050$; normalization requires the corresponding unnormalized weights for all particles. (3) Clearance affects the predicted concentration trajectory, which enters the PD model, so a response can favor parameter values that imply particular concentrations and hence particular clearance values. (4) It can omit uncertainty in estimated concentrations or PK parameters when those estimates are treated as fixed, and can omit dependence between PK and PD uncertainty. The result may understate uncertainty in predicted response.