Why Simulate Histories Rather Than Typical Schedules?
A regular-dose prediction can conceal how real or hypothetical dose histories interact with changing physiology. A dose given late may leave less time for elimination before a measurement; a covariate change can alter clearance partway through that interval. In an NLME simulation, each subject therefore needs an individual timeline: dose events, covariate values or changes, and observation times. The model is then propagated along that timeline rather than along a single average schedule.
The goal here is to generate such histories and compare their predicted concentration profiles. The examples use an intravenous bolus and a one-compartment elimination model so that the event-by-event calculations are transparent. The same simulation logic applies when a more complex NLME model requires numerical integration of additional pharmacokinetic states.
Specify the Model Before Generating Profiles
For the examples, let \(A_i(t)\) denote the amount of drug in the central compartment, \(V_i\) the volume of distribution, and \(CL_i(t)\) the time-varying clearance. Between dose events, the amount follows first-order elimination. An intravenous bolus adds its dose immediately to the amount. Let \(R_i(t)\) be renal function, with reference value \(R_{\mathrm{ref}}\). The individual parameters are defined by
Here \(CL_{\mathrm{pop}}\) and \(V_{\mathrm{pop}}\) are typical parameter values, \(\theta\) is the covariate effect, and \(\eta_{CL,i}\) and \(\eta_{V,i}\) are subject-specific random effects. In the examples, \(R_i(t)\) is constant within each interval between covariate events, and the covariate changes clearance immediately at its recorded event time. The concentration prediction is \(C_i(t)=A_i(t)/V_i\).
For a population simulation, random effects are drawn from their specified distribution, commonly a multivariate normal distribution with mean zero and a covariance matrix \(\Omega\). This tutorial’s numerical examples give the random-effect draws explicitly, so the arithmetic is reproducible. If residual variability is needed, one possible observation model is \(\log C_{ij}^{\mathrm{obs}}=\log C_i(t_{ij})+\epsilon_{ij}\), where \(\epsilon_{ij}\sim N(0,\sigma^2)\). This makes the model prediction the median observation conditional on the parameters; a simulation should state whether it reports predictions or includes simulated residual error.
Propagating the Amount Across Events
On an interval \([a,b)\) with constant clearance, define \(k=CL_i/V_i\). The state equation is \(dA_i/dt=-kA_i\). Solving the linear ODE directly gives the exact interval update, including when \(A_i(a^+)=0\):
The superscripts \(a^+\) and \(b^-\) mean immediately after and immediately before the endpoint, respectively. At a bolus event, \(A_i(t^+)=A_i(t^-)+D\). A covariate event changes the clearance used in the next interval; it does not itself add or remove drug. Applying these updates in chronological order gives the amount at any observation time. Equivalently, each past dose contributes its dose amount multiplied by the exponential survival factors for every interval since that dose.
Draw or specify random effects, calculate individual parameter values, and set the subject’s covariate history.
Combine dose times, covariate-change times, and observation times in chronological order, retaining the defined ordering for events that share a time.
Use the appropriate clearance to propagate the amount to the next event; then apply a dose jump or update the covariate-dependent clearance as required.
At each observation time, convert amount to concentration and compare scenarios using the same model, subject draws, and observation times.
Worked Example: One Subject With a Clearance Change
Consider a hypothetical subject with \(CL_{\mathrm{pop}}=5\) L/h, \(V_{\mathrm{pop}}=50\) L, \(R_{\mathrm{ref}}=90\), and \(\theta=1\). Specify \(\eta_{CL,i}=\log(1.2)\) and \(\eta_{V,i}=\log(0.9)\). The subject receives 400 mg at time 0 h and 250 mg at 5 h. Renal function is 90 until 6 h, then 60; calculate the predicted concentration at 10 h.
The individual volume is \(V_i=50\exp[\log(0.9)]=45\) L. Before 6 h, clearance is \(5(90/90)\exp[\log(1.2)]=6\) L/h. After 6 h, it is \(5(60/90)(1.2)=4\) L/h. Thus \(k\) is \(6/45=0.13333\) h\(^{-1}\) before 6 h and \(4/45=0.088889\) h\(^{-1}\) afterward.
The 400 mg dose has 6 h at the first rate and 4 h at the second rate before the observation. Its remaining amount is \(400\exp[-0.13333(6)-0.088889(4)]\approx125.95\) mg. The 250 mg dose has 1 h at the first rate and 4 h at the second: \(250\exp[-0.13333(1)-0.088889(4)]\approx153.33\) mg. Therefore \(A_i(10)\approx125.95+153.33=279.28\) mg, and \(C_i(10)=279.28/45\approx6.206\) mg/L.
This calculation assumes the clearance change occurs at 6 h and that the dose at 5 h is added before elimination over the subsequent interval. The separate dose contributions make the event history auditable: an incorrectly omitted covariate change would alter the elimination rate applied from 6 h onward.
Compare Covariate and Dose-History Scenarios
A useful simulation contrast changes one feature while keeping the rest of the setup fixed. The comparison below holds the individual, doses, and observation time constant and changes only whether renal function falls. This isolates the model-implied consequence of that covariate history; it does not establish a causal effect in observed patients.
Worked Example: Changing Renal Function Versus Constant Renal Function
For a hypothetical subject, let \(CL_{\mathrm{pop}}=5\) L/h, \(V_{\mathrm{pop}}=50\) L, \(\theta=1\), and both random effects be zero. Give 500 mg at 0 h and 300 mg at 7 h. Compare the 12 h concentration when renal function falls from 90 to 60 at 8 h with the concentration if it remains at 90.
In the changing-covariate scenario, clearance is 5 L/h before 8 h and \(5(60/90)=3.3333\) L/h afterward. For the first dose, the elapsed intervals are 8 h and 4 h, so its contribution at 12 h is \(500\exp[-(5/50)(8)-(3.3333/50)(4)]\approx172.08\) mg. For the 7 h dose, the corresponding intervals are 1 h and 4 h, giving \(300\exp[-(5/50)(1)-(3.3333/50)(4)]\approx207.91\) mg. Total amount is approximately 379.99 mg, so concentration is \(379.99/50\approx7.600\) mg/L.
With renal function held at 90, clearance remains 5 L/h. The 500 mg contribution is \(500\exp[-(5/50)(12)]\approx150.60\) mg, and the 300 mg contribution is \(300\exp[-(5/50)(5)]\approx181.96\) mg. Thus \(C(12)=(150.60+181.96)/50\approx6.651\) mg/L. Under this specified model and history, the changing-covariate prediction is higher by approximately \(7.600-6.651=0.9487\) mg/L because the lower clearance slows elimination after 8 h.
The same approach compares alternative dose schedules. A paired comparison uses the same random-effect draw and covariate history for each schedule; this prevents differences in simulated subject characteristics from obscuring the schedule contrast. Repeating the paired comparison across many subjects yields a distribution of individual contrasts, not just one deterministic difference.
Worked Example: Paired Comparison of Two Irregular Dose Times
Use a hypothetical subject with \(CL_{\mathrm{pop}}=5\) L/h, \(V_{\mathrm{pop}}=50\) L, \(\eta_{CL,i}=\log(1.2)\), \(\eta_{V,i}=0\), and \(\theta=1\). Renal function is 90 through 8 h and 60 afterward. Each schedule includes 500 mg at 0 h; the second 300 mg dose is given at either 6 h or 9 h. Compare concentrations at 12 h.
Individual clearance is 6 L/h before 8 h and 4 L/h afterward, giving rates \(0.12\) and \(0.08\) h\(^{-1}\). The initial dose contributes \(500\exp[-0.12(8)-0.08(4)]=500\exp(-1.28)\approx139.02\) mg at 12 h under either schedule.
For the schedule with the second dose at 6 h, its contribution is \(300\exp[-0.12(2)-0.08(4)]=300\exp(-0.56)\approx171.36\) mg. The total is approximately 310.38 mg and concentration is \(310.38/50\approx6.208\) mg/L. For the 9 h dose, the contribution is \(300\exp[-0.08(3)]=300\exp(-0.24)\approx235.99\) mg. The total is approximately 375.01 mg and concentration is \(375.01/50\approx7.500\) mg/L. The later dose schedule is higher by approximately 1.293 mg/L at 12 h for this subject and model.
This is one paired subject-level result, not a population estimate. To compare schedules for a population, draw each subject’s random effects once, apply the same draw to both schedules, and summarize the paired concentration differences at prespecified times or the paired profile-level quantities of interest.
Simulation Design, Diagnostics, and Alternatives
A population simulation requires explicit choices about what is random and what is held fixed. Random effects represent between-subject parameter variability; dose histories may be fixed scenarios or generated from a stated process; and covariate histories may be supplied as fixed trajectories or generated using a model. These choices answer different questions. For example, holding histories fixed compares model predictions conditional on those histories, whereas generating histories as well includes variation from the history-generation assumptions.
For a clear scenario comparison, use the same random-effect draws and observation times across alternatives. If the objective is to represent independent populations rather than make a paired contrast, independent draws may instead be appropriate. In either case, report the number of simulated subjects, the distributions used for random effects, how histories were generated, and whether residual error was added. A typical-subject profile sets random effects to zero; it is not generally the average of simulated subject profiles because the model is nonlinear.
Check simulated histories before interpreting concentrations. For every subject, verify that event times are ordered, dose amounts and units are valid, covariate values are attached to the intended intervals, and observation times are aligned with the correct pre-dose or post-dose state. Then inspect profiles for expected behavior: between boluses, amount should decrease under positive clearance; a bolus should produce the specified jump in amount; and a decrease in clearance should slow subsequent elimination when other parameters are unchanged. Compare exact interval updates with a numerical solver for selected cases if numerical integration is used.
Common failures include silently carrying an outdated covariate past its change time, applying a dose twice when event records are duplicated, and calculating every interval with a single average clearance. Another frequent mistake is to compare different schedules using unrelated random-effect draws and attribute all observed differences to dosing. The correction is to audit the event sequence, retain the correct interval-specific parameters, and use paired draws when the aim is a within-subject schedule contrast.
Exact exponential updates are efficient and accurate for this one-compartment model with piecewise-constant clearance. A numerical ODE solver is more flexible when the model includes multiple coupled states or continuously varying inputs, but it requires appropriate tolerances and event handling. A deterministic typical-subject calculation is useful for a compact illustration, but omits population variability. Monte Carlo simulation represents that variability and supports distributional summaries, but its conclusions depend on the random-effect and history-generating assumptions. No simulation alone establishes that the assumed model or generated histories accurately represent clinical behavior.
Check Your Understanding
Answer each question, then compare your reasoning with the concise solutions below.
- If clearance is constant over an interval, what fraction of the amount at the start remains after an interval of length \(\Delta t\)?
- In the one-compartment model, what happens to amount at a covariate-change event if there is no dose at that time?
- Why use the same random-effect draws for two schedules when estimating their paired concentration difference?
- What is the 12 h concentration for a 100 mg bolus at 0 h, \(CL=2\) L/h, and \(V=20\) L, with no other dose or covariate change?
- Why is the zero-random-effect profile not generally equal to the average profile in a population simulation?
Solutions. (1) The remaining fraction is \(\exp[-(CL/V)\Delta t]\). (2) Amount is continuous at that event; the covariate changes the clearance used for subsequent propagation. (3) Pairing controls subject-level parameter variation, so the within-draw difference more directly reflects the schedule change. (4) \(k=2/20=0.1\) h\(^{-1}\), so \(A(12)=100\exp(-1.2)\approx30.12\) mg and \(C(12)=30.12/20\approx1.506\) mg/L. (5) Nonlinear transformations mean averaging predictions over random effects need not equal prediction at their mean; for example, in general \(E[\exp(\eta)]\ne\exp(E[\eta])\).