Simulating Concentration-Time Profiles By Formulation

Why Simulate Paired Profiles?

A fitted population pharmacokinetic model can be used to generate concentration-time curves under different formulations. The key design choice is what should be held constant across those curves. To compare test and reference formulations for the same simulated subject, the subject’s shared pharmacokinetic effects should ordinarily be the same under both formulations. The formulation-specific model terms then create the intended difference between the curves.

A curve calculated from fixed model parameters and subject effects is a structural prediction. A simulated observation profile additionally includes residual variation, the remaining variation specified by the observation model. Keeping these two outputs distinct makes a simulation interpretable: it shows which differences are attributed to the fitted structural model and which arise from residual error.

Definition: A paired simulated profile is a set of predicted or simulated concentrations for two or more formulations evaluated for the same simulated subject, using the same shared subject-level random effects. Formulation-specific effects differ only as specified by the fitted model. Residual errors are generated according to the fitted residual model and its dependence assumptions.

Structural Model and Simulation Assumptions

For a concrete illustration, use a one-compartment oral model with first-order absorption and elimination. Let \(D\) be the dose, \(F_f\) the formulation-specific relative availability, \(k_{a,f}\) the absorption rate constant for formulation \(f\), \(CL_i\) the subject’s clearance, and \(V_i\) the subject’s volume of distribution. Define \(k_{e,i}=CL_i/V_i\). Under this model, the predicted concentration is

C_{if}(t)=\frac{D F_f k_{a,f}}{V_i(k_{a,f}-k_{e,i})} \left[\exp(-k_{e,i}t)-\exp(-k_{a,f}t)\right].

This expression assumes a dose entering an absorption depot, first-order absorption and elimination, a well-mixed central compartment, and no lag time. It applies when \(k_{a,f}\ne k_{e,i}\); the equal-rate case is obtained by taking the continuous limit. A more elaborate fitted model requires using its own structural equation rather than substituting this one.

The model’s population parameters describe typical values and the distribution of subject effects. For example, a log-normal random-effects specification might be

\begin{pmatrix} \eta_{CL,i}\\ \eta_{V,i}\\ \eta_{F,i} \end{pmatrix} \sim N(\mathbf{0},\Omega), \qquad CL_i=CL_{\mathrm{typ}}\exp(\eta_{CL,i}),\quad V_i=V_{\mathrm{typ}}\exp(\eta_{V,i}),\quad F_{if}=F_f\exp(\eta_{F,i}).

Here \(\Omega\) is the random-effects covariance matrix. Its off-diagonal entries matter: correlated effects must be generated jointly, not as unrelated draws. The same \(\eta_{CL,i}\), \(\eta_{V,i}\), and shared \(\eta_{F,i}\) are used for both formulations when the model treats them as shared. If the fitted model includes a subject-by-formulation random effect, that term must also be generated and applied as specified; it should not be added merely because the simulation compares formulations.

A simulated observation can then be generated from an observation model. One possible proportional-error model on the log scale is

\log Y_{if}(t_j)=\log C_{if}(t_j)+\epsilon_{ifj},\qquad \epsilon_{ifj}\sim N(0,\sigma^2).

In this version \(C_{if}(t_j)\) is the conditional median of \(Y_{if}(t_j)\), and the error terms are independent across observations unless the fitted model specifies a covariance structure. Other residual models have different implications and must be simulated as fitted. In particular, do not add log-normal error to a model fitted with additive error simply for convenience.

From the Model to a Simulated Profile

The structural equation follows from the model’s differential equations. If the absorption depot contains \(A_g(t)\), then \(dA_g/dt=-k_aA_g\), so \(A_g(t)=D\exp(-k_at)\). With central amount \(A_c(t)\), zero initial central amount, and first-order elimination,

\frac{dA_c}{dt}=F k_aD\exp(-k_at)-k_eA_c.

Multiplying by the integrating factor \(\exp(k_et)\) and integrating from 0 to \(t\) gives

A_c(t)=\frac{F k_aD}{k_a-k_e} \left[\exp(-k_et)-\exp(-k_at)\right].

Since concentration is amount divided by volume, \(C(t)=A_c(t)/V\), which yields the stated expression. Thus, for each simulated subject, first generate the random effects jointly, transform them into individual parameters, and evaluate the structural model at the desired sampling times for each formulation. Generate residuals only if simulated observations, rather than structural predictions alone, are wanted.

Specify the simulation target.
Choose whether to generate structural predictions for selected subjects or a population of simulated subjects, and whether to include residual error.

Generate subject effects.
Draw each subject’s random-effects vector from the fitted distribution using its covariance matrix, or use specified individual estimates when the target is a conditional profile.

Pair the formulations.
Use the same shared subject effects under test and reference; apply each formulation’s fitted structural terms.

Evaluate and, if needed, add residual error.
Calculate concentrations at the chosen times, then draw residuals according to the fitted error model and dependence structure.

Worked Calculations

Worked Example: Paired Structural Predictions for One Subject

Hypothetical data. Consider a 100 mg oral dose, \(CL_{\mathrm{typ}}=4\) L/h, and \(V_{\mathrm{typ}}=20\) L. For this subject, let \(\eta_{CL}=0.20\), \(\eta_V=0\), and \(\eta_F=0.05\). Set \(F_R=0.80\), \(F_T=0.90\), \(k_{a,R}=1.00\) h\(^{-1}\), and \(k_{a,T}=0.70\) h\(^{-1}\). These are illustrative model values, not data from a real study.

The subject-specific parameters are \(CL_i=4\exp(0.20)=4.886\) L/h, \(V_i=20\) L, and \(\exp(\eta_F)=\exp(0.05)=1.0513\). Therefore \(k_{e,i}=4.886/20=0.2443\) h\(^{-1}\). The effective availability values are \(0.80(1.0513)\) for reference and \(0.90(1.0513)\) for test. Evaluate both predictions at \(t=6\) h.

For reference,

C_{iR}(6)= \frac{100(0.80)(1.0513)(1.00)} {20(1.00-0.2443)} \left[\exp(-0.2443\cdot6)-\exp(-1.00\cdot6)\right] \approx 1.271\ {\rm mg/L}.

For test,

C_{iT}(6)= \frac{100(0.90)(1.0513)(0.70)} {20(0.70-0.2443)} \left[\exp(-0.2443\cdot6)-\exp(-0.70\cdot6)\right] \approx 1.569\ {\rm mg/L}.

Both calculations use the same \(CL_i\), \(V_i\), and shared availability multiplier. The different availability and absorption parameters are formulation-specific. These are structural predictions; no residual error has yet been added.

Worked Example: Adding Residual Variation to a Paired Prediction

Continue with the hypothetical structural predictions at 6 h above, and suppose the fitted observation model is \(\log Y=\log C+\epsilon\), with \(\epsilon\sim N(0,0.20^2)\), independent across the two formulation observations. For one illustrative simulation, draw \(\epsilon_R=0.10\) and \(\epsilon_T=-0.15\). The simulated observations are

Y_R=1.271\exp(0.10)=1.405\ {\rm mg/L},\qquad Y_T=1.569\exp(-0.15)=1.350\ {\rm mg/L}.

The noisy simulated observations happen to put reference above test at this time, even though the structural predictions put test above reference. This is not a contradiction: the simulated residuals are part of the observation model. The same subject effects remain paired, while the two residual draws differ. If the fitted model instead included correlated residuals, the two errors would need to be drawn jointly from that fitted structure.

Worked Example: Drawing Correlated Subject Effects

Suppose the hypothetical fitted random-effects distribution for \((\eta_{CL},\eta_V)\) has standard deviations 0.25 and 0.15, with correlation 0.40. One way to generate a draw is to take independent \(Z_1,Z_2\sim N(0,1)\) and set

\eta_{CL}=0.25Z_1,\qquad \eta_V=0.15\left(0.40Z_1+\sqrt{1-0.40^2}\,Z_2\right).

For the illustrative draw \(Z_1=0.80\), \(Z_2=-0.50\), this gives \(\eta_{CL}=0.200\) and \(\eta_V=0.15[0.320-0.4583]=-0.02074\). With \(CL_{\mathrm{typ}}=4\) L/h and \(V_{\mathrm{typ}}=20\) L, the resulting values are \(CL_i=4\exp(0.200)=4.886\) L/h, \(V_i=20\exp(-0.02075)=19.59\) L, and \(k_{e,i}=4.886/19.59=0.2494\) h\(^{-1}\). The same generated effects are used for both formulations. For example, using the dose, availability values, and absorption rates in the previous example gives structural predictions at 6 h of approximately 1.267 mg/L for reference and 1.568 mg/L for test.

The correlation is preserved because both effects depend on \(Z_1\). Generating the two effects independently would instead impose zero correlation, changing the fitted population model and therefore the simulated profile distribution.

Diagnostics, Alternatives, and Limitations

A useful simulation check is to inspect both the generated parameters and the resulting profiles. Confirm that random-effects draws have approximately the intended means, standard deviations, and correlations; that formulation labels map to the correct parameters; and that concentrations are finite and physically interpretable over the simulated time range. Compare the simulated distributions of profiles with model predictions and observed data on the scales used for fitting. Persistent systematic discrepancies may indicate a coding error or model inadequacy; isolated unusual curves can be expected from a population distribution with substantial variability.

Several choices that can look similar answer different questions:

Typical-subject prediction: setting all random effects to zero gives a curve for typical parameter values. It is not a simulated population and generally does not represent the average curve or the range of subject profiles.
Conditional individual prediction: using estimated individual effects gives a profile conditional on information for that subject. This differs from drawing a new subject from the population distribution.
Population simulation: drawing random effects produces a distribution of subject profiles under the fitted population model. This is appropriate when the target is variability across simulated subjects.
Structural versus observation simulation: omitting residuals isolates the fitted concentration curves; adding residuals generates observations under the residual model. They should not be conflated when diagnosing model predictions.

For a paired formulation comparison, shared subject effects should be reused across formulations. However, residual errors should not automatically be made identical across formulations. Under a model with independent residuals, generate separate errors. A shared residual draw can be used as a computational coupling in some simulation comparisons, but it is not the same as simulating the fitted observation model’s independent errors and may change the distribution of paired observed differences. Use it only when the simulation’s explicit purpose and interpretation justify that coupling.

Finally, simulations using fitted parameter estimates as fixed values reflect random-effects and residual variability conditional on those estimates. They do not, by themselves, include uncertainty in the estimated population parameters or model choice. Simulations are also conditional on the structural and error assumptions: a smooth, plausible-looking curve is not proof that those assumptions are adequate. The model’s suitability and residual structure must be checked using the relevant data and diagnostics.

Common pitfall: Simulating test and reference curves with independent clearance and volume effects and then calling them a paired comparison. That creates two different simulated subjects, so their contrast mixes formulation differences with subject differences. For a paired profile, reuse shared subject effects across formulations; vary only model terms that are genuinely formulation-specific, and simulate residuals according to the fitted error structure.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For a paired test-reference simulation, which subject-level effects are reused, and what determines whether any effect differs by formulation?
If \(CL_i=5\) L/h and \(V_i=25\) L, calculate \(k_{e,i}\) and state its units.
In a log-residual model \(\log Y=\log C+\epsilon\), with \(\epsilon\sim N(0,\sigma^2)\), is \(C\) the conditional arithmetic mean or median of \(Y\)?
Why is drawing clearance and volume effects independently inappropriate when the fitted random-effects model specifies a nonzero correlation?
What does a simulation with all random effects set to zero omit?

Solutions. (1) Reuse all subject effects specified as shared in the fitted model; formulation-specific terms differ only when the model explicitly defines them as formulation-dependent. (2) \(k_{e,i}=CL_i/V_i=5/25=0.200\) h\(^{-1}\). (3) \(C\) is the conditional median, because the log error has median zero; the arithmetic mean is \(C\exp(\sigma^2/2)\). (4) Independent draws impose zero correlation and therefore do not reproduce the fitted joint distribution. (5) It omits between-subject random-effects variability, so the result is a typical-parameter prediction rather than a population of subject profiles.