Worked Oral-Absorption Model with Random Absorption Rates

Why random absorption rates change the profile

Two subjects given the same oral dose can differ in how quickly drug leaves the absorption depot. A mixed-effects model represents that between-subject variation separately from the stochastic fluctuations that occur within a subject. This distinction matters: a profile can differ from another subject’s profile because its absorption rate is different, because its central concentration fluctuates, or because of both.

We construct a first-order depot-to-central-compartment model. The depot amount follows its absorption dynamics, while central concentration follows an SDE with stochastic disposition. The random absorption rate is fixed for each subject over the observation period; conditional on that rate, the concentration process has explicit mean and variance dynamics.

Definition: For subject $i$, let $\eta_i\sim N(0,\omega_a^2)$ be an absorption-rate random effect and define $k_{a,i}=\exp(\mu_a+\eta_i)$. The rate $k_{a,i}$ is positive and constant within subject, but varies between subjects. Assume the $\eta_i$ are independent across subjects and independent of their Brownian motions $W_i(t)$.

Let $D$ be the oral dose, $F$ the bioavailability fraction, $V$ the apparent central volume, and $k_e>0$ the first-order elimination rate. At time zero, put the dose in the depot: $A_i(0)=D$, and take $C_i(0)=0$. The depot amount $A_i(t)$ and central concentration $C_i(t)$ are governed by

\begin{aligned} dA_i(t)&=-k_{a,i}A_i(t)\,dt,\\ dC_i(t)&=\left(\frac{F k_{a,i}}{V}A_i(t)-k_e C_i(t)\right)dt +\sigma C_i(t)\,dW_i(t), \end{aligned}

where $\sigma\geq 0$ is the multiplicative diffusion scale. The input term has units of concentration per time, and $\sigma$ has units of inverse square-root time. The model isolates random absorption rates from central stochasticity: $k_{a,i}$ varies across subjects, while $W_i(t)$ drives within-subject fluctuations.

Dynamics conditional on a subject’s absorption rate

Fix a subject and condition on $\eta_i$, so $k_{a,i}=k_a$ is a known positive constant. Solving the depot equation gives $A(t)=D e^{-k_a t}$. Write $S=FD/V$ and $q(t)=S k_a e^{-k_a t}$. The central equation is then a linear SDE with deterministic input $q(t)$.

Its solution, which also shows why the concentration remains nonnegative when the initial concentration and input are nonnegative, is

C(t)=\Phi(t)\left[ C(0)+\int_0^t \Phi(s)^{-1}q(s)\,ds \right], \qquad \Phi(t)=\exp\!\left[\left(-k_e-\frac{\sigma^2}{2}\right)t+\sigma W(t)\right].

This follows by applying the Itô product rule with the integrating factor $\Phi(t)^{-1}$. The representation is useful both for understanding the process and for justifying its moments. On any bounded time interval, the input $q$ is bounded. The ratio $\Phi(t)/\Phi(s)$ is an exponential of a Gaussian Brownian increment plus a deterministic term, so it has finite moments of every positive order. The integral representation and the integral form of Minkowski’s inequality therefore give finite moments of every fixed order for $C(t)$, bounded over bounded time intervals.

Theorem: Conditional on $k_a$, with $C(0)=0$, the mean $m(t)=E[C(t)\mid k_a]$ is $$ m(t)= \begin{cases} S\,\dfrac{k_a}{k_a-k_e}\left(e^{-k_e t}-e^{-k_a t}\right),& k_a\ne k_e,\\[6pt] S k_e t e^{-k_e t},& k_a=k_e. \end{cases} $$ For $v(t)=\operatorname{Var}(C(t)\mid k_a)$, the conditional variance satisfies $$ v'(t)=(-2k_e+\sigma^2)v(t)+\sigma^2m(t)^2,\qquad v(0)=0. $$

To justify the mean equation, apply Itô’s formula to the central SDE and take expectations. The stochastic integral has integrand $\sigma C(t)$; the finite second moments just established imply $E\int_0^T \sigma^2C(t)^2dt<\infty$ for every finite $T$. It is therefore a square-integrable martingale with expectation zero. Thus $m'(t)=q(t)-k_em(t)$, $m(0)=0$. Multiplying by $e^{k_et}$ and integrating gives

m(t)=e^{-k_e t}\int_0^t e^{k_e s}S k_a e^{-k_a s}\,ds.

When $k_a\ne k_e$, evaluating this integral gives the first expression in the theorem. When $k_a=k_e$, the integrand inside the integral is the constant $S k_e$, giving $m(t)=S k_e t e^{-k_e t}$. This also verifies that the equal-rate expression is the continuous limit of the unequal-rate formula.

For the variance, set $r(t)=E[C(t)^2\mid k_a]$. Itô’s formula yields

d(C^2)=\left(2q(t)C+(-2k_e+\sigma^2)C^2\right)dt +2\sigma C^2\,dW.

The stochastic integrand here is $2\sigma C^2$, so its square-integrability requires a fourth moment, not just a second moment. The moment argument using the explicit solution establishes finite fourth moments bounded on bounded intervals; hence $E\int_0^T4\sigma^2C(t)^4dt<\infty$, and the stochastic integral has expectation zero. It follows that $r'=2qm+(-2k_e+\sigma^2)r$. Since $v=r-m^2$ and $m'=q-k_em$, differentiation gives $v'=(-2k_e+\sigma^2)v+\sigma^2m^2$. Solving this linear equation gives the equivalent form

v(t)=\sigma^2\int_0^t e^{(-2k_e+\sigma^2)(t-s)}m(s)^2\,ds.

This derivation makes the roles clear. The conditional mean is the usual first-order absorption and elimination profile; in this model, multiplicative diffusion does not change that mean. The diffusion scale does affect the conditional variance.

Worked profiles with between-subject variation

Worked Example: Conditional Mean and Peak Time

Suppose $D=100$ mg, $F=0.8$, $V=40$ L, $k_e=0.2\ \mathrm{h}^{-1}$, and one subject has $k_a=1.0\ \mathrm{h}^{-1}$. Find the conditional mean concentration at two hours and the time at which that mean peaks.

First, $S=FD/V=(0.8)(100)/40=2$ mg/L. Since $k_a\ne k_e$,

m(2)=2\frac{1.0}{1.0-0.2}\left(e^{-0.4}-e^{-2}\right) =2.5(0.670320-0.135335) \approx 1.337\ \mathrm{mg/L}.

For the peak, differentiate the mean expression. Its derivative is proportional to $-k_e e^{-k_e t}+k_a e^{-k_a t}$. Setting this equal to zero gives $k_a e^{-k_a t}=k_e e^{-k_e t}$, so

t_{\mathrm{peak}}=\frac{\log(k_a/k_e)}{k_a-k_e} =\frac{\log(5)}{0.8}\approx 2.012\ \mathrm{h}.

The peak time here describes the conditional mean curve, not the peak time of every stochastic realization. A realized path can fluctuate around that curve and may have a different local maximum.

Worked Example: Why the Population Mean Is Not a Plug-In Profile

For a transparent two-subject illustration, let $k_a$ equal $0.5\ \mathrm{h}^{-1}$ or $2.0\ \mathrm{h}^{-1}$, each with probability $1/2$. Use the dose and $F/V$ from the preceding example, but set $k_e=0.25\ \mathrm{h}^{-1}$, so $S=2$ mg/L. Calculate the population mean concentration at one hour and compare it with the profile evaluated at the average absorption rate.

For the slower subject,

m_{0.5}(1)=2\frac{0.5}{0.5-0.25} \left(e^{-0.25}-e^{-0.5}\right) =4(0.778801-0.606531) \approx0.6891\ \mathrm{mg/L}.

For the faster subject,

m_{2}(1)=2\frac{2}{2-0.25} \left(e^{-0.25}-e^{-2}\right) =\frac{4}{1.75}(0.778801-0.135335) \approx1.4708\ \mathrm{mg/L}.

Averaging over the two equally likely subjects gives $(0.6891+1.4708)/2\approx1.0799$ mg/L. Their average absorption rate is $(0.5+2)/2=1.25\ \mathrm{h}^{-1}$. Substitution of that rate into the profile instead gives

m_{1.25}(1)=2\frac{1.25}{1.25-0.25} \left(e^{-0.25}-e^{-1.25}\right) =2.5(0.778801-0.286505) \approx1.2307\ \mathrm{mg/L}.

The two answers differ because the concentration formula is nonlinear in $k_a$. In general, averaging subject-specific profiles is not the same operation as evaluating one profile at the average rate.

Worked Example: Reading a Short Simulated Path

Use $D=100$ mg, $F/V=0.02\ \mathrm{L}^{-1}$, $k_a=1\ \mathrm{h}^{-1}$, $k_e=0.2\ \mathrm{h}^{-1}$, and $\sigma=0.3\ \mathrm{h}^{-1/2}$. Starting from $A_0=100$ mg and $C_0=0$, take two Euler–Maruyama steps of size $h=0.1$ h, with Brownian increments $\Delta W_1=0.10$ and $\Delta W_2=-0.05$.

The updates are $A_{n+1}=A_n-k_aA_nh$ and $C_{n+1}=C_n+(F k_a A_n/V-k_eC_n)h+\sigma C_n\Delta W_{n+1}$. On the first step, $A_1=100-100(0.1)=90$ mg and $C_1=0+(0.02)(1)(100)(0.1)=0.200$ mg/L. There is no diffusion increment to concentration on this step because $C_0=0$.

On the second step, $A_2=90-90(0.1)=81$ mg. The input rate is $(0.02)(1)(90)=1.8$ mg/L/h, and the elimination rate at the start of the step is $0.2(0.200)=0.040$ mg/L/h. Therefore

C_2=0.200+(1.8-0.040)(0.1)+(0.3)(0.200)(-0.05) =0.200+0.176-0.003 =0.373\ \mathrm{mg/L}.

For comparison, the conditional mean at $0.2$ h is $2.5(e^{-0.04}-e^{-0.2})\approx0.355$ mg/L. The simulated value is above that mean, as a single realization can be. The step size is fairly coarse, so this short calculation illustrates the update and interpretation; it is not a precision simulation.

Diagnostics, alternatives, and limitations

A useful first diagnostic is to separate the profile features that the model assigns to absorption from those assigned to diffusion. Changing $k_a$ changes the timing and shape of the conditional mean: faster absorption generally moves the mean peak earlier. Changing $\sigma$ leaves that mean unchanged in this model but changes within-subject spread. If simulated profiles vary substantially in their mean peak timing across subjects, inspect the distribution of $k_a$; if paths fluctuate strongly around subject-specific mean curves, inspect $\sigma$.

These interpretations rely on the stated model assumptions: first-order depot loss, a fixed absorption rate within each subject, known dose and bioavailability-volume scaling, first-order elimination, and multiplicative central diffusion. A lognormal rate effect has median $e^{\mu_a}$ and mean $e^{\mu_a+\omega_a^2/2}$; confusing these quantities can misstate the typical rate. Also, a population-average concentration curve averages over both the random absorption rates and the conditional process variation. It is not generally the curve for a “typical” subject.

A fixed-effect model sets $\omega_a^2=0$, so every subject shares one absorption rate. It is simpler, but cannot represent systematic between-subject differences in absorption. A model with random effects on initial conditions, the subject of the earlier material on Random Effects on Initial Conditions, assigns heterogeneity to starting states instead; it does not by itself describe rate variation. Finally, omitting the diffusion term gives a deterministic conditional concentration curve. That alternative can represent absorption-rate heterogeneity but not the within-subject path variability represented here.

The model also has limitations. Multiplicative noise preserves nonnegativity, but does not guarantee that the assumed diffusion law is biologically appropriate at every concentration. A random absorption rate alone cannot represent delayed, transit, or zero-order absorption. Moreover, a noisy concentration path may have multiple local peaks, so the deterministic formula for $t_{\mathrm{peak}}$ should not be applied directly to every simulated path. Check that simulated concentrations remain nonnegative, that changes across subjects are consistent with the assigned rates, and that numerical results are stable to a smaller time step.

Key takeaway: Condition on the subject-specific absorption rate to obtain the subject’s mean curve and within-subject variance. Average those conditional distributions across the random-effects distribution to describe the population; do not substitute an average absorption rate into a nonlinear profile and call it the population mean.

Check Your Understanding

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

If $k_a=k_e$, what is the conditional mean concentration at time $t$?
In the stated model, which parameter changes the conditional mean profile: $k_a$, $\sigma$, or both?
Why can the profile evaluated at $E[k_a]$ differ from the population-average profile?
With $C(0)=0$, why is the first Euler–Maruyama diffusion increment zero when the first simulated concentration update begins at zero?

Solutions. (1) The equal-rate limit is $m(t)=S k_e t e^{-k_e t}$, obtained by integrating $S k_e$ in the integrating-factor formula. (2) $k_a$ changes the conditional mean; $\sigma$ enters the variance equation but not the mean equation. (3) The conditional mean is nonlinear in $k_a$, so averaging its values over the rate distribution need not equal evaluating it at the average rate. (4) The Euler diffusion increment is $\sigma C_0\Delta W_1$, which is zero when $C_0=0$, even though later diffusion increments can be nonzero.)

Worked Oral-Absorption Model with Random Absorption Rates

What You'll Learn