Simulating Correlated PK and PD Variability

Why Correlation Matters in Virtual Participants

A simulated participant is more than a collection of independently sampled parameters. Clearance, volume, potency, and other subject-level characteristics may vary together. If the simulation ignores that dependence, it can generate plausible marginal distributions but implausible combinations of characteristics. Those combinations can alter predicted concentrations, responses, and the range of trajectories across participants.

The goal here is to specify and simulate correlated PK and PD variability, including repeated observations within each participant. The correlation structure is an assumption about the joint distribution of subject effects; it is not something that should be added merely to make a simulation look more realistic. Its scale, direction, and consequences should be explicit and checked.

Specify the Joint Distribution Before Simulating

For a simple PK/PD model, let \(\eta_{CL,i}\) and \(\eta_{EC50,i}\) be participant \(i\)'s random effects on clearance and the concentration producing half-maximal response. A common specification is a bivariate normal distribution with mean zero and covariance matrix \(\Omega\). The diagonal entries are variances; the off-diagonal entries are covariances.

Here \(\omega_{CL}\) and \(\omega_{EC50}\) are the standard deviations of the effects on their respective model scales, and \(\rho\) is their correlation, between \(-1\) and \(1\). The covariance matrix must be positive semidefinite. For this two-effect matrix, that condition requires \(|\rho|\leq 1\); in larger systems, all eigenvalues of the covariance matrix must be nonnegative. A Cholesky factorization is available when the matrix is positive definite.

For positive parameters, a lognormal parameterization keeps individual values positive:

In this parameterization, \(CL_{\mathrm{pop}}\) and \(EC50_{\mathrm{pop}}\) are typical values: with mean-zero normal effects, they are medians of the individual parameter distributions, not generally their arithmetic means. Correlation is specified between the normal effects on the log scale. The resulting correlation between the positive individual parameters need not equal \(\rho\).

To generate the effects, draw independent standard normal values \(Z_1,Z_2\), then form:

The shared \(Z_1\) term induces the intended correlation. For more effects, the same idea is applied using a matrix factor: if \(\Omega=LL^\mathsf{T}\) and \(Z\) is a vector of independent standard normal draws, then \(\eta=LZ\) has covariance \(\Omega\). Generate this vector once per participant and retain it for that participant’s entire trajectory.

From Subject Effects to Repeated PK and PD Observations

For illustration, consider a one-compartment model after an intravenous bolus dose \(D\), with fixed volume \(V\) and no additional structural complexity. The concentration at time \(t\) is:

Let the concentration-driven response follow an \(E_{\max}\) model, with maximum response \(E_{\max}\) and baseline \(E_0\):

Here \(m_i(t)\) is the model-predicted response, \(Y_{ij}\) is the observed response at visit \(j\), and \(\epsilon_{ij}\) is residual error. In the example below, residuals are assumed independent across visits and participants, normally distributed with mean zero and standard deviation \(\sigma\), conditional on the participant’s parameters. This is a modeling assumption, not a universal property of repeated measurements. If residuals are serially correlated or have a concentration-dependent variance, that structure must be specified rather than silently replaced by independent constant-variance error.

Why Observations Become Dependent

Even when residual errors are independent across visits, repeated observations from the same participant are generally dependent because they share the same random effects. The law of total covariance makes this distinction precise. Let \(Y_{i1}\) and \(Y_{i2}\) be two observations. Conditional on the subject effects, suppose their residual errors are independent. Then their conditional covariance is zero, and:

The first equality decomposes total covariance into covariance of conditional means plus expected conditional covariance. The second uses the stated conditional independence assumption. Since both predicted responses depend on the same \(\eta_i\), their means will often covary. The direction and size depend on the model, time points, parameter effects, and nonlinear response function. If residuals are themselves correlated, the expected conditional covariance term is not zero and must be included.

Correlation between PK and PD effects also does not imply that every PK measure and every PD outcome will have that same correlation. For example, the PD response depends on concentration, which depends on PK parameters; the nonlinear concentration-response function transforms the joint variability. The observed response correlation can therefore differ in magnitude, and sometimes in direction, from the correlation between any two underlying random effects.

Check the Dependence You Actually Induced

The covariance matrix describes a model input. A simulated data set is a finite random sample, so its empirical correlations will not equal the specified values exactly. Useful checks compare the intended structure with both generated subject parameters and repeated observations. Check on the scale where the model defines the effects, then also inspect the scales and outcomes relevant to the simulation.

For two jointly normal effects with standard deviations \(\omega_1,\omega_2\) and correlation \(\rho\), the corresponding lognormal parameters \(P_1=P_{1,\mathrm{pop}}e^{\eta_1}\) and \(P_2=P_{2,\mathrm{pop}}e^{\eta_2}\) have correlation:

This result follows from the normal moment-generating function: \(E[e^{\eta_k}]=e^{\omega_k^2/2}\), and \(E[e^{\eta_1+\eta_2}]=e^{(\omega_1^2+\omega_2^2+2\rho\omega_1\omega_2)/2}\). Subtracting the product of the means gives the covariance; dividing by the two lognormal standard deviations gives the stated correlation. It provides a reference for checking transformed parameters, not a substitute for checking the generated sample.

Diagnostics, Alternatives, and Limitations

A practical check is to simulate a large set of participants before running the full trial simulation. Examine the random-effects scatterplot and covariance or correlation matrix, then inspect the marginal distributions of individual parameters. For repeated data, plot trajectories by participant and compare within-participant variation with between-participant variation. If the model assumes independent residuals, residuals after accounting for the predicted trajectory should not show an unmodeled pattern across visits. These checks can reveal implementation errors; they cannot establish that the assumed correlation structure is empirically correct.

Common implementation failures include entering a correlation where a covariance is required, using a matrix that is not positive semidefinite, and drawing each random effect independently despite specifying a nonzero correlation. Another frequent error is generating new subject effects at every visit. That destroys the intended stable between-participant heterogeneity and produces trajectories that do not correspond to persistent individual parameters.

A further pitfall is to infer that the simulated response correlation must equal the input random-effect correlation. That is wrong because response is a nonlinear function of concentration and PD parameters, and residual error adds further variation. The correction is to state which scale the input correlation describes, derive or simulate the implications on the outcome scale, and check those implications directly. Correlation alone also does not establish a biological mechanism; it is a joint-distribution assumption in the model.