Why Covariate Dependence Matters
A virtual cohort is not merely a collection of plausible values for individual physiological parameters. Its members should also have plausible combinations of values. Body size, lean and adipose mass, organ volumes, and regional blood flows are related in real populations. Sampling each quantity independently can therefore create individuals whose component values look reasonable in isolation but do not make sense together.
For PBPK simulations, this matters because model predictions depend on combinations of inputs. For example, an organ’s volume and its perfusion jointly affect how much drug can distribute to that organ and how quickly it receives drug. A cohort that breaks their dependence can produce misleading exposure distributions even if every input has a reasonable marginal range.
The objective is not to make every pair of inputs correlated. It is to represent the relationships justified for the target population, while preserving necessary identities and constraints. Correlations may differ by population: an association among body size and organ volume in an obesity cohort need not match the association in an older-adult or disease-state cohort.
Specify the Cohort Before Sampling
Before generating values, define the target population, the physiological inputs, and which quantities are directly sampled versus derived. A cohort intended to represent older adults with a disease should not be assembled by mixing unrelated marginal distributions from young healthy adults and then labeling the result as diseased. Specify the relevant age, obesity, or disease strata, and decide whether the cohort represents one stratum or a mixture.
For each input, specify a marginal distribution or a plausible range, units, and any conditioning variables. Then specify the dependence structure. For instance, an organ-volume model might depend on body size and age, while a disease-state flow model might depend on disease severity as well as body size. These are model choices that require suitable population information or defensible physiological assumptions; correlation alone does not establish a mechanism.
A common construction starts with a vector of correlated standard normal variables. Let \(R\) be the target correlation matrix, with ones on the diagonal. It must be positive semidefinite; otherwise it cannot be the correlation matrix of any real random vector. When \(R\) is positive definite, a Cholesky factorization gives \(R=LL^\mathsf{T}\). If \(\epsilon\) has independent standard normal components, then \(Z=L\epsilon\) has mean zero and covariance matrix \(R\). For normally distributed inputs, use \(X_j=\mu_j+\sigma_j Z_j\). For a nonnormal input with cumulative distribution function \(F_j\), a marginal transformation is \(X_j=F_j^{-1}(\Phi(Z_j))\), where \(\Phi\) is the standard normal cumulative distribution function.
The transformation gives the specified marginal distributions, but it does not generally preserve the original Pearson correlations. It preserves the latent ordering, while nonlinear transformations can change linear correlation. Therefore, check the correlations in the generated cohort on the physiological scale actually used by the PBPK model. If constraints such as positivity or mass balance require truncation or recalculation, check again after applying them.
A Practical Construction Sequence
Specify the population and any age, obesity, or disease strata. Avoid blending strata unless the intended cohort is explicitly a mixture.
For each individual, distinguish sampled covariates from derived quantities. Record identities such as total body mass being composed of component masses, and regional flows summing to total cardiac output.
Use a correlation structure or conditional models that represent plausible relationships for the target cohort. Check that the correlation matrix is valid.
Transform the correlated draws to the required scales, then derive constrained quantities rather than sampling them independently.
Check marginal summaries, pairwise relationships, physiological identities, boundary values, and combinations of inputs at the individual level.
Worked Example: Generating Correlated Body Size and Organ Volume
Suppose a hypothetical adult cohort uses body size \(X_1\) in kilograms and liver volume \(X_2\) in litres. For illustration only, assume normal marginal distributions with means 75 kg and 1.50 L, standard deviations 12 kg and 0.20 L, and target correlation 0.60. These are chosen teaching values, not population estimates.
The target correlation matrix is \(R=\begin{pmatrix}1&0.60\\0.60&1\end{pmatrix}\). A Cholesky factor is \(L=\begin{pmatrix}1&0\\0.60&0.80\end{pmatrix}\), since \(LL^\mathsf{T}=R\). For independent standard normal draws \(\epsilon_1=0.50\) and \(\epsilon_2=-1.00\), the correlated scores are \(Z_1=0.50\) and \(Z_2=0.60(0.50)+0.80(-1.00)=-0.50\).
Transforming to the stated normal marginals gives \(X_1=75+12(0.50)=81\) kg and \(X_2=1.50+0.20(-0.50)=1.40\) L. This individual has a larger-than-average body size and a smaller-than-average liver volume; a positive cohort-level correlation does not mean every individual must follow the average trend. The matrix encodes dependence across repeated draws, not a rule that every pair of values must move together.
The generated cohort should be checked for nonpositive or otherwise implausible values, since normal distributions have unbounded tails. If values are truncated, the correlation after truncation may differ from 0.60. In a real construction, marginal distributions and the target correlation would need support for the intended population.
Mass Balance and Derived Covariates
Some physiological inputs are linked by accounting identities rather than merely correlated. If body mass is represented as lean mass, adipose mass, and other mass, then the components must sum to body mass. Sampling all four independently can violate that identity. A safer approach is to sample a subset using a joint or conditional model, then derive the remaining quantity, while checking that every component is positive and plausible.
This distinction matters for correlations too. If total body mass is defined as a sum of components, correlations among total mass and its components are partly determined by that definition. They should not be treated as independent empirical relationships to impose arbitrarily. A physiological model can instead represent how component shares vary with the relevant covariates, and calculate total mass from those components or calculate a component from the total and the others.
Worked Example: Checking Composition in a Hypothetical Cohort
Consider a hypothetical virtual individual with body mass 90 kg, lean mass 55 kg, and adipose mass 28 kg. Suppose the model also represents “other” mass, so the accounting identity is body mass = lean mass + adipose mass + other mass.
The implied other mass is \(90-55-28=7\) kg. The three components sum to \(55+28+7=90\) kg, so the identity holds. The component shares are \(55/90=0.611\), \(28/90=0.311\), and \(7/90=0.0778\); they sum to approximately 1.000, with the slight rounding due to reporting.
If adipose mass had instead been sampled as 40 kg while body mass and lean mass remained unchanged, the implied other mass would be \(90-55-40=-5\) kg. That draw fails the nonnegative-component constraint even though each of the first three numbers might appear individually possible. The correction is not to ignore the negative value or to silently clip it to zero; the joint sampling or derivation must be revised so the identity and physiological bounds are satisfied together.
Regional Blood Flows Are Constrained
Regional blood flows require special care because they are components of a total. If a set of modeled systemic flows \(Q_1,\ldots,Q_k\) partitions cardiac output \(Q_{\mathrm{total}}\), then \(\sum_j Q_j=Q_{\mathrm{total}}\). Independently sampling each flow can violate this sum. Even if the sum is enforced afterward, the adjustment can change the intended relationships among flows, body size, and other covariates.
One option is to model total cardiac output and regional fractions separately. Positive fractions can be generated from latent scores using a softmax transformation: \(f_j=\exp(a_j)/\sum_k\exp(a_k)\). The flows are then \(Q_j=f_jQ_{\mathrm{total}}\), ensuring that they are positive and sum to the total. The latent scores can depend on covariates or share correlated random variation. Because fractions sum to one, an increase in one fraction must be balanced by decreases in one or more others. Thus, negative correlations among fractions can arise from the constraint itself; they do not necessarily indicate direct physiological antagonism.
Worked Example: Allocating Total Flow Without Breaking the Sum
For a hypothetical individual, let total cardiac output be 5.4 L/min. Suppose three regional latent weights, after exponentiation, are 1.2 for region A, 0.8 for region B, and 1.0 for region C. Their sum is 3.0, so the fractions are \(f_A=1.2/3.0=0.400\), \(f_B=0.8/3.0=0.267\), and \(f_C=1.0/3.0=0.333\), to three significant figures.
The regional flows are \(Q_A=0.400(5.4)=2.16\) L/min, \(Q_B=0.267(5.4)\approx1.44\) L/min, and \(Q_C=0.333(5.4)\approx1.80\) L/min. Using the unrounded fractions, the sum is \(2.16+1.44+1.80=5.40\) L/min. Rounding is not a change to the underlying identity.
This construction preserves the total by design. It does not, by itself, establish that these fractions or their relationships to body size are appropriate for any particular population. Those relationships must be specified and checked for the target cohort.
Diagnostics, Alternatives, and Limitations
Diagnostics should examine both the population-level structure and individual records. Compare simulated marginal summaries with their targets; inspect scatterplots or pairwise correlations for key inputs; and verify positivity, mass balance, and flow sums for every virtual individual. Also inspect extremes and combinations, such as high adipose mass paired with unusually small body mass, or high organ volumes paired with unexpectedly low total flow. A large cohort can still be systematically implausible.
A correlation-matrix approach is compact and useful when dependence is approximately summarized by pairwise relationships. Its limitations are that the matrix must be valid, linear correlation may not describe nonlinear or tail dependence, and Gaussian latent dependence can underrepresent unusual joint extremes. Conditional models can represent relationships such as organ volume given body size and age more directly, but require suitable conditional assumptions and parameters. Stratified cohort construction can retain different relationships across age or disease groups, whereas pooling strata can create misleading overall correlations. Each approach must be judged by whether it represents the intended population and preserves the required constraints.
A second pitfall is assuming that a target correlation matrix automatically becomes the observed correlation matrix after transformation, truncation, or constraint enforcement. It does not. The relevant diagnostic is the dependence in the final physiological inputs supplied to the PBPK model. If it differs materially from the intended structure, revise the construction and document the change.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Why can matching the marginal distributions of body size and organ volume still produce an implausible cohort?
- What condition must a proposed correlation matrix satisfy to be usable for correlated normal draws?
- A virtual individual has total body mass 80 kg, lean mass 50 kg, and adipose mass 25 kg. If these are the only components, does the mass balance hold? What is the discrepancy?
- Why must regional flows be checked after converting sampled fractions into flows?
- Does a positive correlation between body size and organ volume require every larger individual to have a larger organ volume? Explain.
Solutions. (1) The marginals do not specify how values occur together; independent pairing can break population dependence and yield implausible combinations. (2) It must be positive semidefinite; a Cholesky factor exists when it is positive definite. (3) The components sum to \(50+25=75\) kg, leaving a 5 kg discrepancy from the stated 80 kg; an additional component or revised inputs are needed. (4) They must satisfy positivity and sum to the modeled total cardiac output, and rounding or later adjustments can affect those properties. (5) No. Positive correlation is a population-level tendency; individual values can depart from it, and the association is not a deterministic rule.