A Worked Fit of a Bolus Effect-Compartment Model

From a Bolus to a Delayed Response

After an intravenous bolus, plasma concentration can decline while the measured response is still increasing. An effect compartment represents this delay by allowing effect-site concentration to approach plasma concentration gradually. A worked fit makes clear that the model must explain the whole response time course: both the rise and the decline, not merely the response at its peak.

The calculation below uses a simulated, hypothetical data set. The plasma profile is treated as known from a preceding pharmacokinetic analysis; the response observations are used to estimate the equilibration rate and the pharmacodynamic parameters. This separation is an assumption of the example, not a claim that plasma estimates are error-free in practice.

Definition: For a bolus with initial plasma concentration $C_0$ and monoexponential elimination rate $k$, let $C_p(t)=C_0e^{-kt}$. An effect compartment with equilibration rate $K_{e0}$ is defined by $$ \frac{dC_e(t)}{dt}=K_{e0}\bigl(C_p(t)-C_e(t)\bigr), \qquad C_e(0)=0. $$ Here $C_p(t)$ and $C_e(t)$ have concentration units, and $k$ and $K_{e0}$ have units of inverse time.

For the response, we use a baseline plus a stimulatory Emax relationship. The response approaches a maximum increment as effect-site concentration increases. The response errors are assumed independent, normally distributed, and of constant variance for the least-squares and uncertainty calculations discussed below.

Model: At response observation time $t_i$, $$ R_i=R_0+\frac{E_{\max}C_e(t_i)}{EC_{50}+C_e(t_i)}+\varepsilon_i, \qquad \varepsilon_i\sim N(0,\sigma_R^2). $$ The parameters are baseline response $R_0$, maximum response increment $E_{\max}$, effect-site concentration producing half of that increment $EC_{50}$, and equilibration rate $K_{e0}$. $EC_{50}$ has concentration units; $R_0$, $E_{\max}$, and $\sigma_R$ have response units.

Bolus Solution and Fitting Criterion

The effect-compartment equation is linear. Multiplying it by the integrating factor $e^{K_{e0}t}$, integrating from zero to $t$, and substituting $C_p(s)=C_0e^{-ks}$ gives

\begin{aligned} \frac{d}{dt}\left(e^{K_{e0}t}C_e(t)\right) &=K_{e0}C_0e^{(K_{e0}-k)t},\\ e^{K_{e0}t}C_e(t) &=K_{e0}C_0\int_0^t e^{(K_{e0}-k)s}\,ds,\\ C_e(t) &=\frac{C_0K_{e0}}{K_{e0}-k} \left(e^{-kt}-e^{-K_{e0}t}\right), \qquad K_{e0}\ne k. \end{aligned}

The factor $K_{e0}$ in the numerator is essential: it follows directly from the differential equation. If $K_{e0}=k$, taking the limit gives $C_e(t)=C_0K_{e0}t e^{-K_{e0}t}$. For positive rates and $C_0$, the unequal-rate expression is nonnegative, and $C_e(0)=0$.

For a fixed plasma profile, the four response parameters can be estimated by nonlinear least squares. At each trial $K_{e0}$, calculate $C_e(t_i)$, then calculate the predicted response. The fitting procedure chooses parameter values to minimize the sum of squared response residuals. Because $K_{e0}$ affects every predicted $C_e(t_i)$, it is estimated jointly with $R_0$, $E_{\max}$, and $EC_{50}$, rather than by fitting only one time point.

RSS=\sum_{i=1}^{n}\left(R_i-\widehat R_i\right)^2, \qquad \widehat R_i=R_0+\frac{E_{\max}C_e(t_i)}{EC_{50}+C_e(t_i)}.

The model assumes a single bolus, a monoexponential plasma profile over the time range, an initially drug-free effect compartment, and an instantaneous response to $C_e(t)$ through the stated Emax function. The fitting assumptions additionally include independent response errors with constant variance; normality is used for likelihood-based uncertainty calculations. These assumptions should be checked against the design and residuals rather than inferred from a successful optimizer run.

Worked Example: Calculating Effect-Site Concentration and Response

Suppose the hypothetical plasma profile is $C_p(t)=7.50e^{-0.200t}$ mg/L, with time in hours. A fit reports $K_{e0}=0.750$ h$^{-1}$, $R_0=20.2$ response units, $E_{\max}=82.0$ response units, and $EC_{50}=2.20$ mg/L. Calculate the prediction at $t=2$ h.

First calculate the effect-site concentration using the bolus solution:

\begin{aligned} C_e(2) &=\frac{(7.50)(0.750)}{0.750-0.200} \left(e^{-(0.200)(2)}-e^{-(0.750)(2)}\right)\\ &=10.2273(0.670320-0.223130)\\ &\approx 4.574\ \text{mg/L}. \end{aligned}

Then substitute this concentration into the response equation:

\widehat R(2)=20.2+\frac{(82.0)(4.574)}{2.20+4.574} \approx 20.2+55.37 =75.57\ \text{response units}.

If the observed response at 2 h is 74.8 units, the residual, defined as observed minus predicted, is $74.8-75.57=-0.77$ response units. A negative residual means that the model prediction is above this observation. The concentration calculation retains the $K_{e0}$ numerator; omitting it would not solve the stated effect-compartment equation.

A Simulated Fit Across the Time Course

The following seven hypothetical response observations were simulated around a delayed response curve. A nonlinear least-squares fit to this data set, with the plasma profile treated as fixed, gives the rounded parameter estimates used in the preceding example. The table shows predictions calculated from those estimates. Minor discrepancies in the last decimal can arise from displaying rounded parameters and predictions.

Time (h)	Observed response (units)	Predicted response (units)	Residual (units)
0	20.40	20.20	0.20
0.5	61.00	61.43	-0.43
1	71.30	70.78	0.52
2	74.80	75.57	-0.77
4	74.20	73.50	0.70
6	66.80	67.28	-0.48
8	60.20	59.65	0.55

Worked Example: Summarizing the Fit

Using the rounded residuals in the table, the residual sum of squares is

\begin{aligned} RSS &=(0.20)^2+(-0.43)^2+(0.52)^2+(-0.77)^2\\ &\quad +(0.70)^2+(-0.48)^2+(0.55)^2\\ &=0.0400+0.1849+0.2704+0.5929+0.4900+0.2304+0.3025\\ &=2.1111\ \text{response-units}^2. \end{aligned}

There are $n=7$ observations and four fitted response parameters, so the residual degrees of freedom are $df=7-4=3$. Under the independent, constant-variance error model, the residual standard deviation estimate is

s_R=\sqrt{\frac{RSS}{df}} =\sqrt{\frac{2.1111}{3}} \approx 0.839\ \text{response units}.

This summarizes the residual scale; it does not show by itself that the model is adequate. In particular, seven observations and three residual degrees of freedom provide little power to distinguish random scatter from a systematic pattern.

Interpreting the Time Course and Uncertainty

The fitted $K_{e0}=0.750$ h$^{-1}$ controls how quickly effect-site concentration responds to plasma concentration. With this bolus profile, $C_e(t)$ rises from zero, reaches a maximum, and then falls. For unequal positive rates, differentiating the bolus solution and setting its derivative to zero gives the peak time

t_{\mathrm{peak}}= \frac{\log(K_{e0}/k)}{K_{e0}-k}.

Here $t_{\mathrm{peak}}=\log(0.750/0.200)/(0.750-0.200)\approx 1.322/0.550=2.40$ h. Since the Emax response function increases with $C_e$ when $E_{\max}>0$ and $EC_{50}>0$, the predicted response has its maximum at the same time in this model. The fitted response is therefore still near its peak at 4 h and declines afterward, even though plasma concentration has been falling since the bolus.

An uncertainty assessment should account for the fact that changes in $K_{e0}$ can be partly offset by changes in $EC_{50}$ or $E_{\max}$. One local summary is the covariance matrix from the nonlinear fit, based on the curvature of the objective near its minimum. With only seven data points and four fitted parameters, such a local approximation may be unstable. Profiling $K_{e0}$ is an informative alternative: hold it at successive values, refit the other response parameters, and examine how much the fit deteriorates. For this simulated example, suppose the profile-likelihood interval is approximately 0.36 to 1.31 h$^{-1}$. It is broad relative to the point estimate, indicating meaningful uncertainty about the rate even though the fitted curve tracks the observations reasonably closely. This interval is a model-based uncertainty summary, not a 95% probability statement about the fixed parameter.

Worked Example: Reading the Profile Interval

The hypothetical profile interval is 0.36 to 1.31 h$^{-1}$, with point estimate 0.750 h$^{-1}$. Its width is $1.31-0.36=0.95$ h$^{-1}$, and the point estimate lies within it. Under the stated likelihood assumptions, the interval identifies values of $K_{e0}$ that are not excluded by the chosen profile-likelihood cutoff. It does not establish that all values in the interval are equally plausible, nor does it mean there is a 95% probability that the fixed $K_{e0}$ lies inside this particular interval.

If the interval were very wide or extended toward zero, the observations would provide weak information about equilibration speed. A useful follow-up is to examine the profile curve and the fitted response curves at its boundary values: if they nearly overlap at the sampled times, the design has limited ability to distinguish those rates.

Diagnostics, Alternatives, and Limitations

Inspect residuals against time and against predicted response, and compare observed and fitted trajectories. Residuals that remain mostly positive early and negative late, or show another ordered pattern, suggest that the model misses some feature of the time course. A single large residual may reflect an unusual observation, but it should not be removed without a scientific and data-quality justification. The assumed constant variance is questionable if residual spread systematically increases with predicted response; dependence between repeated measurements on one individual also violates the independent-error assumption unless it is modeled.

The sampling schedule matters. Observations concentrated only near the response peak may constrain the curve's height but provide little information about how quickly it rises or declines. Measurements on both limbs, including early times and a sufficiently informative decline, help separate equilibration delay from response parameters. In this example, the fixed plasma profile is another limitation: uncertainty in $C_0$ or $k$ is not included in the reported fit interval. Joint estimation or uncertainty propagation may be needed when those quantities are uncertain.

Approach	What it represents	Main limitation or distinction
Direct concentration-response model	Response is linked to contemporaneous $C_p(t)$.	Cannot represent a lag between plasma concentration and response without another mechanism.
Bolus effect-compartment model	Response depends on a filtered concentration $C_e(t)$ that approaches $C_p(t)$ at rate $K_{e0}$.	Requires a suitable plasma profile and sufficiently informative observations to estimate delay and response parameters.

A lower residual sum of squares alone does not establish that a delayed model is preferable; the earlier material on testing whether an effect compartment improves fit addresses formal comparison. Here the focus is fitting and interpreting the selected bolus model. Even a well-fitting $K_{e0}$ is a parameter of this model's time course. The estimate does not, by itself, prove a particular physiological mechanism.

Common Mistake

Wrong reasoning: “The formula is proportional to the difference of exponentials, so the coefficient is $C_0/(K_{e0}-k)$.” This drops the rate multiplying the concentration difference.
Correction: Solve the stated differential equation, which yields $C_0K_{e0}/(K_{e0}-k)$. Then calculate each $C_e(t_i)$ before calculating response predictions and residuals. An error in this coefficient changes the whole fitted trajectory and every subsequent fit statistic.
Another diagnostic error: Calling a model adequate because its predictions are close to the peak observation ignores the rising and falling limbs. Assess the complete residual pattern and the uncertainty in the fitted parameters.

Key takeaway: A bolus effect-compartment fit links plasma input, delayed effect-site concentration, and response through explicit equations. Check the bolus solution's rate factor, calculate predictions across the full time course, and interpret $K_{e0}$ together with its uncertainty and the model assumptions.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

For $C_p(t)=C_0e^{-kt}$ and $C_e(0)=0$, what factor appears in the numerator of the bolus solution for $C_e(t)$?
In the worked fit, what is the residual at 2 h, using observed minus predicted response?
Why can observations near the response peak alone give weak information about $K_{e0}$?
What does the profile interval 0.36 to 1.31 h$^{-1}$ say, and what does it not say?
Name one assumption that can fail when repeated response measurements are collected from the same individual.

Solutions: (1) The numerator is $C_0K_{e0}$. (2) It is $74.8-75.57=-0.77$ response units. (3) Different equilibration rates can produce similar peak responses while differing more on the rising and falling limbs. (4) Under the model and chosen cutoff, it is a range of $K_{e0}$ values not excluded by the profile analysis; it is not a 95% probability statement about the fixed parameter. (5) Independence of errors may fail because measurements from one individual can be correlated.

Approach	What it represents	Main limitation or distinction
Direct concentration-response model	Response is linked to contemporaneous \(C_p(t)\).	Cannot represent a lag between plasma concentration and response without another mechanism.
Bolus effect-compartment model	Response depends on a filtered concentration \(C_e(t)\) that approaches \(C_p(t)\) at rate \(K_{e0}\).	Requires a suitable plasma profile and sufficiently informative observations to estimate delay and response parameters.