Propagating Plasma and Ke0 Uncertainty Into Effect Predictions

Why Propagate Uncertainty?

An effect prediction depends on more than a single estimated plasma concentration or a single fitted value of \(K_{e0}\). Both the plasma trajectory and the rate at which the effect site equilibrates with plasma are uncertain. Because the effect-site concentration depends on the history of plasma concentrations, uncertainty in either input can affect the predicted response at a given time and the predicted time and magnitude of a response peak.

The goal is to carry these uncertainties through the model, rather than report an effect prediction calculated only from point estimates. Simulation is useful because the effect-compartment calculation and the concentration–response relationship can both be nonlinear. It also allows each simulated trajectory to produce its own peak, so peak timing uncertainty is not mistaken for uncertainty in response at one fixed time.

Definition: Uncertainty propagation is the process of transforming draws from an uncertainty distribution for model inputs or parameters into corresponding draws of a model output. In this tutorial, each draw of a plasma concentration trajectory and \(K_{e0}\) produces an effect-site trajectory and, if specified, a response trajectory.

Model and Simulation Assumptions

For a one-compartment effect-site model, the effect-site concentration \(C_e(t)\) follows plasma concentration \(C_p(t)\) according to

\frac{dC_e(t)}{dt}=K_{e0}\bigl(C_p(t)-C_e(t)\bigr),

where \(K_{e0}>0\) has units of inverse time and \(C_p\) and \(C_e\) use the same concentration units. The model assumes that the effect site approaches the current plasma concentration at a rate proportional to their difference. It does not say that the two concentrations are equal at every time.

To propagate uncertainty, specify a joint distribution for the uncertain inputs. A draw might include a full plasma concentration trajectory generated from uncertain pharmacokinetic parameters, a draw of \(K_{e0}\), and any uncertain response parameters. The joint distribution matters: if the inputs are correlated, simulated draws should preserve that dependence. Each draw must also use a consistent set of units, dosing history, initial conditions, and model assumptions.

For a response that increases with effect-site concentration, one possible model is

E(t)=E_0+\frac{E_{\max}C_e(t)}{EC_{50}+C_e(t)},

where \(E_0\) is baseline response, \(E_{\max}\) is the maximum increase above baseline, and \(EC_{50}>0\) is the effect-site concentration producing half of that increase. These parameters may be fixed or included among the uncertain inputs. This tutorial treats the model as specified; simulation propagates its uncertainty but does not establish that the model is appropriate.

Calculating Effect-Site Draws

When plasma concentration is constant at \(C_p\) over an interval of length \(\Delta t\), the differential equation has an exact update:

C_e(t+\Delta t) = C_p+\bigl(C_e(t)-C_p\bigr)e^{-K_{e0}\Delta t}.

For a changing plasma curve, apply an appropriate numerical method or a piecewise solution using the plasma concentrations over successive intervals. The same method and time grid should be used consistently across simulation draws. A draw is one coherent model realization: use its plasma curve and its \(K_{e0}\) together to calculate its \(C_e(t)\), then calculate its response if needed.

Worked Example: Propagating Uncertainty at One Time

Suppose five hypothetical paired draws represent uncertainty in plasma concentration and \(K_{e0}\) over a one-hour interval. For this demonstration, plasma concentration is constant during the interval, \(C_e(0)=0\), and the response model has \(E_0=20\) response units, \(E_{\max}=80\) response units, and \(EC_{50}=5\) mg/L. The pairs \((C_p,K_{e0})\) are \((8,0.4)\), \((10,0.5)\), \((9,0.7)\), \((11,0.3)\), and \((10,0.8)\), with concentration in mg/L and \(K_{e0}\) in h\(^{-1}\).

For each draw, the exact update reduces to \(C_e(1)=C_p(1-e^{-K_{e0}})\), because \(C_e(0)=0\) and \(\Delta t=1\) h. For the first draw,

C_e(1)=8(1-e^{-0.4})=8(0.32968)=2.637\ \text{mg/L}.

Applying the same calculation to all five draws gives \(C_e(1)\) values of approximately 2.637, 3.935, 4.531, 2.851, and 5.507 mg/L. For the first draw, the response is

E(1) =20+\frac{80(2.637)}{5+2.637} \approx20+27.62 =47.62\ \text{response units}.

Transforming each effect-site value in the same way gives responses of approximately 47.62, 55.22, 58.03, 49.05, and 61.93 response units. In this small demonstration set, the median effect-site concentration is 3.935 mg/L and the median response is 55.22 response units. The response range is 49.05 to 61.93 response units. These summaries describe only the five illustrative draws; five draws are not enough to estimate stable uncertainty limits for an actual analysis.

Summarizing Peak Effect and Peak Timing

A pointwise summary describes predictions at a prespecified time. For example, the median of \(E(t)\) across draws at \(t=2\) h summarizes the predicted response at two hours. It does not generally summarize the distribution of each draw’s peak response. To summarize peaks, first find the maximum within each simulated response trajectory and record the time at which it occurs. Then summarize those peak responses and times across draws.

For an increasing \(E_{\max}\) response with fixed positive \(E_{\max}\) and \(EC_{50}\), the response is an increasing function of \(C_e\). Consequently, its maximum occurs at the time of maximum \(C_e\). If response parameters vary by draw, or the response is not monotone in \(C_e\), calculate the response trajectory for each draw and find its maximum directly.

Worked Example: Summarizing Simulated Peaks

Consider five hypothetical simulated effect-site trajectories. Their peak concentrations and peak times are, respectively: 2.8 mg/L at 1.5 h, 3.610 mg/L at 2.0 h, 4.0 mg/L at 2.5 h, 5.0 mg/L at 3.0 h, and 3.2 mg/L at 2.0 h. Use the response model \(E=20+80C_e/(5+C_e)\), with concentration in mg/L and response in response units.

For the trajectory with a peak \(C_e\) of 3.610 mg/L, the peak response is

E=20+\frac{80(3.610)}{5+3.610} \approx20+33.5424 =53.54\ \text{response units}.

Calculating the response at each trajectory’s peak concentration gives approximately 48.72, 53.54, 55.56, 60.00, and 51.22 response units. The sorted peak responses are 48.72, 51.22, 53.54, 55.56, and 60.00, so their median is 53.54 response units. The peak times are 1.5, 2.0, 2.5, 3.0, and 2.0 h; sorted, their median is 2.0 h. The ranges are 48.72–60.00 response units for peak response and 1.5–3.0 h for peak time. As in the previous example, these five draws illustrate the calculation rather than provide reliable interval estimates.

With many draws, empirical quantiles can summarize the distributions of peak response and peak time. For example, the 2.5th and 97.5th percentiles form a central 95% simulation interval under the chosen input distribution and model. This is a statement about the simulated distribution, not a claim that there is a 95% probability that a fixed, unknown parameter lies in that interval. The meaning of the interval depends on how the input uncertainty distribution was constructed.

Why Paired Draws Matter

The concentration at the effect site depends jointly on plasma concentration and \(K_{e0}\). Drawing their values independently when their uncertainty is correlated can change the resulting distribution of effect-site concentrations and responses. This issue applies to plasma trajectories as well: concentrations at nearby times are usually related, so a simulated curve should retain its time-course structure rather than be assembled from unrelated concentration draws at each time.

Worked Example: The Effect of Pairing Uncertain Inputs

Suppose the marginal draws for \(C_p\) are 8 and 12 mg/L, and those for \(K_{e0}\) are 0.3 and 0.8 h\(^{-1}\). With \(C_e(0)=0\), constant plasma for one hour, and positive pairing \((8,0.3)\), \((12,0.8)\), the effect-site values are

8(1-e^{-0.3})=2.073\ \text{mg/L}, \qquad 12(1-e^{-0.8})=6.608\ \text{mg/L}.

Their mean is \((2.073+6.608)/2=4.341\) mg/L. Pairing the same marginal values in the opposite way, \((8,0.8)\), \((12,0.3)\), gives

8(1-e^{-0.8})=4.405\ \text{mg/L}, \qquad 12(1-e^{-0.3})=3.110\ \text{mg/L},

with mean \((4.405+3.110)/2=3.758\) mg/L. The two examples have identical marginal values but different pairings and therefore different propagated results. They are not interchangeable representations of uncertainty. In an analysis, use the joint draws implied by the parameter uncertainty calculation rather than pairing marginal draws arbitrarily.

Simulation Steps and Diagnostics

Specify the uncertain inputs.
Identify which plasma parameters or concentration trajectories, \(K_{e0}\), initial conditions, and response parameters are uncertain. State the source and meaning of the joint uncertainty distribution.

Generate coherent draws.
For each simulation replicate, draw the inputs jointly and retain dependence across parameters and across time in the plasma curve.

Propagate through the effect compartment.
Calculate a complete \(C_e(t)\) trajectory for each draw using the specified initial conditions and a sufficiently accurate numerical or piecewise solution.

Calculate response and peaks.
Apply the response model to each draw’s effect-site trajectory. Find its peak response and peak time before summarizing those quantities across draws.

Check stability and plausibility.
Increase the number of draws and refine the time grid to assess whether uncertainty summaries and peak locations are stable. Inspect trajectories for numerical artifacts or implausible predictions.

A credible propagation requires the model and simulation to behave sensibly. Check that concentrations and response parameters obey their intended constraints, that units agree, and that the initial effect-site concentration is appropriate to the dosing history. Compare solutions at successively finer time steps; a peak time that shifts substantially as the grid is refined may be a numerical artifact. If a peak falls at the start or end of the simulated window, extend the window or state that the maximum is only the maximum within the observed prediction interval.

Also check whether a few extreme draws dominate the summaries. Such draws may reflect genuine uncertainty, an unsuitable input distribution, or an implementation problem; they should not be discarded without justification. A narrow simulated interval is not evidence of low overall uncertainty if the input distribution omits important sources, such as uncertainty in the plasma trajectory or response parameters.

Alternatives and Limitations

A first-order approximation, such as a delta-method calculation, can be efficient when uncertainty is small, the model is smooth near the estimates, and a local approximation is adequate. It can be less reliable when the transformation is strongly nonlinear, uncertainty is broad, or the output is a maximum over time. Simulation handles those features more directly, but its summaries are only as appropriate as the joint input distribution and model supplied to it.

A point-estimate calculation is simpler but does not describe propagated uncertainty. Pointwise intervals are useful for response at specified times, but they are not intervals for a trajectory’s maximum: the time of the maximum can vary across draws. For peak uncertainty, calculate one peak per draw. The resulting intervals still do not address structural uncertainty in the effect-compartment model or prove that \(K_{e0}\) has a particular physiological interpretation. Those limitations require separate scientific and model assessment.

Common Mistake

Common mistake: “The upper pointwise response limit is the upper limit for peak response, and the peak occurs at the time of the median response curve’s maximum.” These conclusions do not generally follow. Pointwise quantiles at each time need not describe the distribution of per-draw maxima, and the maximum of a median curve need not equal the median of individual peak times. The correction is to find each draw’s peak response and peak time first, then summarize those draw-level quantities.

Check Your Understanding

Answer each question using the model assumptions stated above.

With constant \(C_p\), \(C_e(0)=0\), \(C_p=6\) mg/L, \(K_{e0}=0.5\) h\(^{-1}\), and \(\Delta t=1\) h, calculate \(C_e(1)\).
Why should correlated plasma and \(K_{e0}\) draws not automatically be sampled independently?
What should be calculated first to summarize uncertainty in peak response: pointwise response quantiles or the peak in each simulated trajectory?
What does a central 95% simulation interval describe, and what does it not establish by itself?

Solutions. (1) \(C_e(1)=6(1-e^{-0.5})=6(0.39347)=2.361\) mg/L. (2) Independent sampling can break the joint dependence in the uncertainty distribution and change the propagated effect-site and response distributions. (3) Find the peak in each simulated trajectory first, then summarize those peak values across draws. (4) It describes the central 95% of the simulated output distribution under the chosen inputs and model; by itself, it does not establish a 95% probability that a fixed unknown parameter lies in the interval or account for uncertainty omitted from the simulation.