Worked Oral Absorption Example at Periodic Steady State

From Repeated Oral Doses to a Periodic Distribution

A repeated oral dose does not enter the central compartment immediately: it first increases the amount in a depot, from which absorption feeds the central amount. With process noise, the state at a fixed phase of each dosing interval is random. At steady state, its distribution repeats from interval to interval, even though the process continues to evolve within each interval.

The aim here is to calculate that distribution for a linear two-compartment model and then find the concentration law at any time after a dose. The calculations depend on constant coefficients, additive Brownian noise, fixed periodic doses, and stable elimination and absorption dynamics. They do not imply that every model with random inputs has a Gaussian steady-state law.

Definition: Let $G(t)$ be the depot amount and $A(t)$ the central amount. Between doses, take $$ dX(t)=M X(t)\,dt+B\,dW(t),\qquad X(t)=\begin{pmatrix}G(t)\\A(t)\end{pmatrix},\qquad M=\begin{pmatrix}-k_a&0\\k_a&-k_e\end{pmatrix}. $$ Here $k_a>0$ is the absorption rate, $k_e>0$ is the elimination rate, and $W(t)$ is a two-dimensional standard Brownian motion. At each dose time $n\tau$, the state jumps by $q=(D,0)^\mathsf{T}$, where $D$ is the amount entering the depot and $\tau>0$ is the dosing interval. The observed concentration is $C(t)=A(t)/V$, with fixed volume $V>0$.

For the numerical examples, assume independent Brownian components and $B=\operatorname{diag}(\sigma_g,\sigma_a)$. More generally, the formulas below allow any constant matrix $B$, including one that gives correlated component noise. We use the convention that $X_n^-$ denotes the state immediately before dose $n$; the state immediately after it is $X_n^+=X_n^-+q$.

One-Dose-Interval Transition

The matrix exponential for the drift is $$ E(t)=e^{Mt}= \begin{pmatrix} e^{-k_a t}&0\\[2pt] \displaystyle k_a\frac{e^{-k_a t}-e^{-k_e t}}{k_e-k_a}&e^{-k_e t} \end{pmatrix},\qquad k_a\ne k_e. $$ If $k_a=k_e$, the lower-left entry is $k_a t e^{-k_a t}$, obtained by taking the limit. The linear SDE solution over an interval of length $t$ is $X(t)=E(t)X(0)+Z_t$, where $Z_t$ is a mean-zero Gaussian stochastic integral with covariance $$ Q_t=\int_0^t E(s)BB^\mathsf{T}E(s)^\mathsf{T}\,ds. $$ The covariance is positive semidefinite because, for every vector $v$, $v^\mathsf{T}Q_tv=\int_0^t\|B^\mathsf{T}E(s)^\mathsf{T}v\|^2ds\geq0$. The stochastic integral is independent of the state at the start of the interval when that state is determined by the past Brownian motion.

Worked Example: The Periodic Mean Under Repeated Doses

Take $k_a=1$, $k_e=2$, $\tau=\ln 2$, and $D=100$ mg. First calculate the one-interval transition matrix. Since $e^{-\tau}=1/2$, $e^{-2\tau}=1/4$, and $k_e-k_a=1$, $$ R=E(\tau)=\begin{pmatrix}1/2&0\\1/4&1/4\end{pmatrix}. $$

At periodic steady state, the pre-dose mean $\mu^-$ must equal the mean one interval later. The mean transition is $\mu^-=R(\mu^-+q)$, so $(I-R)\mu^-=Rq$. Here $Rq=(50,25)^\mathsf{T}$, giving $$ \begin{pmatrix}1/2&0\\-1/4&3/4\end{pmatrix} \begin{pmatrix}\mu_G^-\\\mu_A^-\end{pmatrix} = \begin{pmatrix}50\\25\end{pmatrix}. $$ The first row gives $\mu_G^-=100$ mg. The second gives $-25+(3/4)\mu_A^-=25$, hence $\mu_A^-=200/3$, or about $66.67$ mg. Immediately after the dose the mean depot amount is $200$ mg, while the central mean remains $66.67$ mg. This is a periodic mean, not a claim that the individual trajectory repeats exactly.

Steady-State Distribution at Dose Times

Writing $R=E(\tau)$, the pre-dose states satisfy the random affine recursion $$ X_{n+1}^-=R X_n^-+Rq+Z_n,\qquad Z_n\sim N(0,Q_\tau). $$ The $Z_n$ are independent across intervals and independent of the state at the beginning of their respective intervals. Since the eigenvalues of $M$ are $-k_a$ and $-k_e$, both negative, the eigenvalues of $R$ have absolute value less than one. This stability is what makes accumulated past doses and noise fade from the dose-time state.

Theorem: For $k_a,k_e,\tau>0$, fixed $q$, and constant $B$, the dose-time recursion has a unique invariant distribution. It is Gaussian, $N(\mu^-,P)$, with $$ \mu^-=(I-R)^{-1}Rq,\qquad P=RPR^\mathsf{T}+Q_\tau =\sum_{j=0}^{\infty}R^jQ_\tau(R^\mathsf{T})^j. $$ The sum converges, and $P$ is positive semidefinite.

Proof. Because the eigenvalues of $R$ have absolute value less than one, $R^j$ tends to zero geometrically (a polynomial factor in $j$ in the repeated-eigenvalue case does not prevent convergence). Thus the displayed covariance series converges. Every term is positive semidefinite, so its limit $P$ is positive semidefinite. Multiplying the series on the left by $R$ and on the right by $R^\mathsf{T}$ removes its first term and yields $RPR^\mathsf{T}+Q_\tau=P$. The mean formula is well-defined since $I-R$ is invertible.

To establish existence, define a random vector as the convergent sum of the deterministic terms $\sum_{j=1}^{\infty}R^j q$ and the independent Gaussian terms $\sum_{j=0}^{\infty}R^j Z_{-1-j}$, with the indices representing past intervals. The deterministic series converges; the Gaussian series converges in mean square because its covariance series is finite. Its distribution is Gaussian with the stated mean and covariance, and shifting the interval index leaves its distribution unchanged under the recursion. For uniqueness, start the recursion from any finite-valued random state. After $m$ intervals, its contribution from that initial state is $R^mX_0^-$, which tends to zero in probability. The other terms converge to the same Gaussian series. Hence the distribution converges to the stated invariant law, and any invariant initial distribution must equal that limit. $\square$

Worked Example: Solving for the Pre-Dose Covariance

Continue with $k_a=1$, $k_e=2$, and $\tau=\ln2$, and set $\sigma_g=0.2$, $\sigma_a=0.1$. For $s\geq0$, $$ E(s)=\begin{pmatrix}e^{-s}&0\\e^{-s}-e^{-2s}&e^{-2s}\end{pmatrix}. $$ Using independent noise components, integration gives $$ Q_\tau= \begin{pmatrix} \sigma_g^2\int_0^\tau e^{-2s}ds& \sigma_g^2\int_0^\tau(e^{-2s}-e^{-3s})ds\\ \text{same}& \sigma_g^2\int_0^\tau(e^{-s}-e^{-2s})^2ds+ \sigma_a^2\int_0^\tau e^{-4s}ds \end{pmatrix}. $$

The four needed integrals are $3/8$, $1/12$, $25/960$, and $15/64$, respectively. For example, the third is $3/8-2(7/24)+15/64=25/960$. Therefore $$ Q_\tau=\begin{pmatrix}0.015&0.0033333\\0.0033333&0.0033854\end{pmatrix}. $$ Write $P=\begin{pmatrix}p&r\\r&s\end{pmatrix}$. Substitution of $R=\begin{pmatrix}1/2&0\\1/4&1/4\end{pmatrix}$ into $P=RPR^\mathsf{T}+Q_\tau$ gives $$ p=0.25p+0.015,\quad r=0.125p+0.125r+0.0033333,\quad s=0.0625p+0.125r+0.0625s+0.0033854. $$ Solving in order gives $p=0.0200$, $r=0.0066667$, and $s=0.0058333$, in squared amount units. As a check, $ps-r^2\approx0.0000722>0$, consistent with a positive definite covariance matrix.

Distribution Throughout the Dosing Interval

At elapsed time $t$, where $0\leq t<\tau$, the state is the transition from the post-dose state $X_n^-+q$. Applying the transition formula to a Gaussian pre-dose state gives a Gaussian state at every fixed phase of the interval.

Theorem: If the pre-dose state has its invariant distribution $N(\mu^-,P)$, then at elapsed time $t$ after the dose, $$ X(t)\sim N\!\left(E(t)(\mu^-+q),\ E(t)PE(t)^\mathsf{T}+Q_t\right). $$ Consequently, $C(t)=A(t)/V$ is normally distributed with mean $m_C(t)=e_2^\mathsf{T}E(t)(\mu^-+q)/V$ and variance $v_C(t)=e_2^\mathsf{T}[E(t)PE(t)^\mathsf{T}+Q_t]e_2/V^2$, where $e_2=(0,1)^\mathsf{T}$.

Derivation. The interval solution is $X(t)=E(t)X_n^++Z_t$, with $X_n^+=X_n^-+q$. The invariant pre-dose state is Gaussian and independent of the new interval noise $Z_t$. Thus the mean is $E(t)(\mu^-+q)$. The two independent random contributions have covariances $E(t)PE(t)^\mathsf{T}$ and $Q_t$, so their covariances add. A fixed linear function of a Gaussian vector is Gaussian; applying the linear map $e_2^\mathsf{T}/V$ gives the concentration formulas. As $t$ approaches $\tau$, these moments return to the pre-dose moments: the mean by $\mu^-=R(\mu^-+q)$, and the covariance by $P=RPR^\mathsf{T}+Q_\tau$.

Worked Example: Concentration at the Midpoint of the Interval

Use the parameters and pre-dose moments from the covariance example, take $V=10$ L, and examine $t=\tau/2=(\ln2)/2$. Then $e^{-t}=1/\sqrt2$, $e^{-2t}=1/2$, and the lower-left entry of $E(t)$ is $1/\sqrt2-1/2$. The post-dose mean is $(200,66.6667)^\mathsf{T}$ mg. Hence the central mean is $$ (1/\sqrt2-1/2)(200)+(1/2)(66.6667) =74.7547\ \text{mg}, $$ and the concentration mean is $74.7547/10=7.47547$ mg/L.

For the central variance, first calculate the central-component variance generated by noise during the half interval. Its two contributions are $$ 0.04\int_0^t(e^{-s}-e^{-2s})^2ds =0.04\left(\frac14-\frac{2(1-1/(2\sqrt2))}{3}+\frac{3}{16}\right) \approx0.0002614, $$ and $0.01\int_0^t e^{-4s}ds=0.001875$. Thus $(Q_t)_{22}\approx0.0021364$. The propagated pre-dose covariance contributes $$ \begin{pmatrix}1/\sqrt2-1/2&1/2\end{pmatrix} P \begin{pmatrix}1/\sqrt2-1/2\\1/2\end{pmatrix} \approx0.0036969. $$ The central amount variance is therefore approximately $0.0058333$ mg$^2$. Dividing by $V^2=100$ gives a concentration variance of $0.00005833$ (mg/L)$^2$. The midpoint concentration law is approximately $N(7.47547,0.00005833)$, with the second argument denoting variance.

Checks, Alternatives, and Limitations

A useful numerical check is to verify both fixed-point equations: $\mu^-=R(\mu^-+q)$ and $P=RPR^\mathsf{T}+Q_\tau$. Also verify that $P$ is symmetric positive semidefinite and that the interval-end moments approach the pre-dose moments. Material disagreement can indicate a timing error (using post-dose moments as pre-dose moments), a transposed matrix in the covariance update, or an incorrect diffusion covariance integral.

A deterministic repeated-dose calculation is the main alternative when process noise is absent: it gives the periodic mean trajectory but no distribution around that trajectory. Simulation can represent more complicated models, but finite simulation error and time-step error make it less direct than these exact linear-Gaussian formulas. The formulas here rely on a linear drift and additive noise; nonlinear drift, state-dependent diffusion, random doses, and random subject-specific parameters generally require different distributional calculations. In particular, additive Gaussian noise has unbounded support, so this model can assign positive probability to negative amounts or concentrations. That is a modeling limitation, not evidence that the computed Gaussian law is algebraically wrong.

Common pitfall: The pre-dose covariance is not simply $Q_\tau$. That would count noise from only the latest interval and discard the accumulated contribution from earlier intervals. The correct steady-state covariance includes every past interval, $P=\sum_{j=0}^{\infty}R^jQ_\tau(R^\mathsf{T})^j$, equivalently the fixed-point equation $P=RPR^\mathsf{T}+Q_\tau$.

Check Your Understanding

Use the definitions and formulas above. Solutions are included after the questions.

Why is the steady-state mean computed just before a dose using $Rq$, rather than $q$, in its fixed-point equation?
For the numerical example, what are the mean depot and central amounts immediately after the dose?
What equation should a candidate pre-dose covariance satisfy?
How is the central concentration variance obtained from the two-dimensional state covariance?

Solutions. (1) The state receives $q$ at the start of the interval, then evolves for $\tau$; its contribution at the next pre-dose time is $E(\tau)q=Rq$. (2) The pre-dose mean is $(100,66.6667)^\mathsf{T}$ mg, so the post-dose mean is $(200,66.6667)^\mathsf{T}$ mg. (3) It must satisfy $P=RPR^\mathsf{T}+Q_\tau$, with $Q_\tau$ the covariance accumulated over one interval. (4) If $e_2=(0,1)^\mathsf{T}$, the concentration variance is $e_2^\mathsf{T}Pe_2/V^2$ at pre-dose time, or $e_2^\mathsf{T}[E(t)PE(t)^\mathsf{T}+Q_t]e_2/V^2$ at elapsed time $t$.