Worked Diffusion Identifiability Analysis For A Dosed Compartment

From a Known Bolus to a Diffusion Estimate

A concentration trajectory after a dose combines two different features: a systematic decline governed by elimination and random fluctuations around that decline. Estimating the fluctuation scale requires separating it from the deterministic dose response, observation error, and uncertainty in the drift. The key is to write the transition distribution for the actual dosing schedule before calculating residuals.

Consider a one-compartment concentration $X_t$. Between dose events, suppose it follows the Ornstein–Uhlenbeck SDE

dX_t=-kX_t\,dt+\sqrt{q}\,dW_t,\qquad k>0,\quad q\geq 0,

where $k$ is the elimination rate, $W_t$ is standard Brownian motion, and $q$ is the diffusion variance per unit time. Thus the instantaneous diffusion coefficient is $\sqrt q$, not $q$. A known bolus of amount $D_j$ at time $s_j$ produces the jump $X_{s_j}=X_{s_j-}+D_j/V$, where $V$ is the known volume of distribution. We treat the concentration measurements in this tutorial as exact observations of the state unless stated otherwise.

Definition: A diffusion parameter $q$ is structurally identifiable under a specified observation scheme if two distinct admissible values of $q$, with the other model parameters fixed, cannot generate the same probability distribution for the observations. It is practically identifiable when the available finite data and design determine it with useful precision.

For observation times $t_i<t_{i+1}$, let $\Delta_i=t_{i+1}-t_i$. Integrating the linear SDE over this interval gives

X_{t_{i+1}}=e^{-k\Delta_i}X_{t_i} +\sum_{t_i<s_j\leq t_{i+1}}\frac{D_j}{V}e^{-k(t_{i+1}-s_j)} +\sqrt q\int_{t_i}^{t_{i+1}}e^{-k(t_{i+1}-u)}\,dW_u.

The known dose contributions enter the conditional mean. By the Itô isometry, the stochastic integral is normal with mean zero and variance $\int_{t_i}^{t_{i+1}}e^{-2k(t_{i+1}-u)}du$. Define $a_i=e^{-k\Delta_i}$, let $m_i$ denote the conditional mean displayed above, and define

g_i=\frac{1-e^{-2k\Delta_i}}{2k}.

Then $X_{t_{i+1}}\mid X_{t_i}\sim N(m_i,qg_i)$. Every $g_i$ is strictly positive for $k>0$ and $\Delta_i>0$. In particular, a known bolus changes the conditional mean but does not create process variance.

Identifiability and the Exact Transition Likelihood

Suppose $k$, $V$, dose times, and dose amounts are known, and the state is observed exactly at $n+1$ times. Given the observed state at each interval's start, the residual is $r_i=X_{t_{i+1}}-m_i$. The independent Brownian increments over disjoint intervals imply conditionally independent transition residuals, with $r_i\sim N(0,qg_i)$.

Theorem: If at least one observed interval has positive length, then $q$ is structurally identifiable from the exact-state transition distributions when the other parameters are fixed. For $n$ positive-length intervals and $q>0$, the maximum-likelihood estimator is $$ \hat q=\frac{1}{n}\sum_{i=0}^{n-1}\frac{r_i^2}{g_i}. $$ If the conditional means and $k$ are known, then $n\hat q/q\sim\chi^2_n$.

Proof of identifiability. For any positive-length interval, the conditional distribution is normal with variance $qg_i$, where $g_i>0$ is fixed by the known $k$ and interval length. If two values $q_1$ and $q_2$ produce the same conditional distribution, their variances must agree: $q_1g_i=q_2g_i$. Dividing by $g_i>0$ gives $q_1=q_2$. At $q=0$, the transition is degenerate at its conditional mean, which is distinct from the positive-variance transition for every $q>0$. This proves structural identifiability over $q\geq0$.

Proof of the estimator and its sampling distribution. For $q>0$, the conditional log-likelihood, omitting constants independent of $q$, is

\ell(q)=-\frac{n}{2}\log q-\frac{1}{2q} \sum_{i=0}^{n-1}\frac{r_i^2}{g_i}.

Writing $S=\sum_i r_i^2/g_i$, the derivative is $\ell'(q)=-n/(2q)+S/(2q^2)$. Setting this equal to zero gives $q=S/n$. When $S>0$, the second derivative at $S/n$ is $-n/(2\hat q^2)<0$, so this is the likelihood maximum. Under the stated model, $r_i/\sqrt{qg_i}$ are independent standard normal variables. Therefore $S/q$, the sum of their squares, has a chi-square distribution with $n$ degrees of freedom, proving $n\hat q/q\sim\chi^2_n$. If every residual is exactly zero, the likelihood instead increases without bound as $q$ approaches zero; the estimate is on the boundary. $\square$

Identifiability does not guarantee useful precision. A single transition identifies $q$ in the model-distribution sense, but provides only one random squared residual for estimation. Moreover, for short intervals $g_i$ is small, so a given observed residual can imply a large estimate of $q$, and finite data can yield substantial uncertainty.

Worked Analysis of a Bolus Trajectory

Worked Example: Estimating Diffusion After a Single Bolus

A bolus gives the known initial concentration $X_0=10$. Suppose $k=0.5$ per time unit, measurements are one time unit apart and exact, and the observed concentrations are $X_1=6.265307$, $X_2=3.700101$, and $X_3=2.544224$. There are no later doses. Estimate $q$ using the exact transition likelihood.

The one-step multiplier is $a=e^{-0.5}=0.6065306597$. The conditional predictions and residuals, calculated using this unrounded multiplier, are:

For the first interval, the prediction is $a(10)=6.065306597$, so $r_0=6.265307-6.065306597=0.199700403$.

For the second interval, the prediction is $a(6.265307)=3.8001008$ (approximately), so $r_1=3.700101-3.8001008\approx-0.0999998$.

For the third interval, the prediction is $a(3.700101)=2.2442247005\ldots$, which rounds to $2.244225$ at six decimal places. Thus $r_2=2.544224-2.2442247005\ldots=0.2999992995\ldots$, which rounds to $0.299999$, not $0.300000$.

All three intervals have $g=(1-e^{-1})/(2(0.5))=1-e^{-1}=0.6321205588$. Hence

\hat q= \frac{(0.199700403)^2+(-0.0999998)^2+(0.2999992995)^2} {3(0.6321205588)} \approx 0.0738.

The estimate is in concentration-squared per time unit. The calculation uses the exact OU transition variance, not the Euler approximation $q\Delta$. Here the known dose determines $X_0$ and the deterministic decay; it does not add a variance term.

Worked Example: An Exact Uncertainty Interval

Continue with the three residuals above and suppose the drift and conditional means are treated as known. Using the rounded estimate, $S=\sum_i r_i^2/g\approx 3(0.0738)=0.2214$; using the unrounded residuals gives $S$ approximately $0.2213$. For a 95% confidence interval, the exact pivot gives

\left[ \frac{S}{\chi^2_{3,0.975}}, \frac{S}{\chi^2_{3,0.025}} \right],

where $\chi^2_{3,p}$ is the $p$-quantile of a chi-square distribution with three degrees of freedom. Using $\chi^2_{3,0.975}\approx9.348$, $\chi^2_{3,0.025}\approx0.216$, and $S\approx0.2213$, the interval is approximately $[0.0237,1.02]$. Its broad range is a warning: structural identifiability does not make a three-transition estimate precise. This interval is exact only under the example's assumptions, including known drift and exact state observations.

Confounding, Diagnostics, and Alternatives

A concentration measurement is often noisy. Let $Y_i=X_{t_i}+\epsilon_i$, where the independent errors satisfy $\epsilon_i\sim N(0,R)$, and suppose $k$ is known. For a single interval, the contrast $Y_{i+1}-a_iY_i$, after subtracting any known dose contribution, contains process noise and measurement noise. Its variance is

$$ qg_i+R(1+a_i^2). $$

With only one such contrast, the variance identifies a combination of $q$ and $R$, not both separately. For example, with $k=0.5$, $\Delta=1$, and an observed contrast variance of $0.4$, $R=0$ corresponds to $q=0.4/0.6321206\approx0.633$. But $R=0.1$ corresponds to $q=[0.4-0.1(1+e^{-1})]/0.6321206\approx0.416$. Both parameter pairs produce the same contrast variance.

With repeated observations, the contrasts are not independent when measurement error is present: adjacent contrasts share a measurement error, and their covariance can carry information about $R$. A state-space likelihood can represent this structure; treating every noisy measurement as an exact state and applying the exact-state formula generally understates uncertainty or biases the diffusion estimate. If both $k$ and $q$ are unknown, $k$ affects the conditional mean and also $g_i$. They should be estimated jointly or assessed with a profile likelihood, rather than treating an uncertain estimate of $k$ as fixed without qualification.

The exact-transition likelihood is the natural approach when the model's transition distribution is available and the state data are exact. A common alternative is an Euler residual estimate, which uses $q\Delta_i$ as the transition variance. For finite intervals in this model, the exact variance is $q(1-e^{-2k\Delta_i})/(2k)$; replacing it with $q\Delta_i$ can misstate the scale when $k\Delta_i$ is not small. For noisy or latent states, filtering or a state-space likelihood accounts for the unobserved trajectory and observation process, but requires a correctly specified observation model and more computation.

Useful diagnostics begin with the model assumptions and the design. Check that dose amounts, times, units, and the initial concentration are represented correctly; a missed or mis-timed bolus changes the predicted mean and can masquerade as large residual variation. Inspect standardized residuals $r_i/\sqrt{\hat qg_i}$: unusually large values, systematic trends, or dependence across intervals may indicate a wrong drift, an unmodelled dose, outliers, or an inadequate diffusion model. These checks are not proof of model adequacy, especially with few observations.

Common pitfall: Subtracting a predicted concentration and calling the squared residual “the diffusion variance” ignores both elapsed time and the OU transition factor. Under the exact-state model, each squared residual is scaled by $g_i$, and the likelihood estimator averages $r_i^2/g_i$. If observations have error, even this correction is insufficient because the residuals also contain measurement noise.

The analysis has limits. It assumes additive, constant diffusion and linear elimination, so it can assign positive probability to negative concentrations; it may be a local approximation rather than a globally realistic concentration model. The confidence interval above is sensitive to fixed-parameter assumptions, and estimates near $q=0$ require boundary-aware interpretation. Sparse sampling, short observation windows, uncertain drift, and unknown observation noise can all make a structurally identifiable parameter practically weak.

Check Your Understanding

Try each question, then compare with the solutions below.

For $k=0.25$ and an interval of length $2$, calculate the exact transition variance as a multiple of $q$.
Why does a known bolus inside an observation interval change its conditional mean but not its process variance?
For exact-state observations with known drift and positive-length intervals, write the maximum-likelihood estimator of $q$.
With only one noisy transition contrast and both $q$ and $R$ unknown, explain why the two parameters are confounded.

Solutions. (1) Here $2k\Delta=1$, so $g=(1-e^{-1})/(0.5)=2(1-e^{-1})\approx1.26424$, and the variance is $1.26424q$. (2) The bolus contributes a known, deterministic jump. The Brownian integral still has variance $q\int_{t_i}^{t_{i+1}}e^{-2k(t_{i+1}-u)}du=qg_i$; adding a known constant changes the mean, not the variance. (3) With residuals $r_i=X_{t_{i+1}}-m_i$, the estimator is $\hat q=n^{-1}\sum_i r_i^2/g_i$. (4) The contrast variance is $qg_i+R(1+a_i^2)$, one equation involving two unknowns, so different pairs can yield the same distribution for that single contrast.

Worked Diffusion Identifiability Analysis For A Dosed Compartment

What You'll Learn