Linear Gaussian State-Space Models for Pharmacokinetics

Why Separate the State from the Assay?

A pharmacokinetic concentration measured in a sample is not necessarily the underlying concentration process itself. The latent PK state evolves over time, while an assay reports a noisy observation of that state. A state-space model makes these two sources of variation explicit: process noise represents random evolution between observation times, and assay error represents measurement variation conditional on the state.

This distinction is central when a filtering method is used to infer an unobserved concentration or compartment amount from serial measurements. The model must specify what the state contains, how it changes between times, how the state maps to an observation, and the covariance of each noise source. The equations below define the linear Gaussian case; they do not imply that every PK process or assay is well described by Gaussian errors.

Definition: A discrete-time linear Gaussian state-space model consists of a latent state vector $x_k$ and an observation vector $y_k$ satisfying $$ x_k=F_kx_{k-1}+a_k+w_k,\qquad y_k=H_kx_k+d_k+v_k. $$ Here $F_k$ is the state-transition matrix, $a_k$ a known transition offset, $H_k$ the observation matrix, and $d_k$ a known observation offset. The process noise $w_k\sim N(0,Q_k)$, and the observation or assay error $v_k\sim N(0,R_k)$. The noises are independent across times and mutually independent, and are independent of the initial state. The initial state has a specified Gaussian distribution. Covariance matrices $Q_k$ and $R_k$ are symmetric positive semidefinite; when an ordinary observation density is required, $R_k$ is positive definite.

The state $x_k$ is a vector of latent quantities at time $t_k$, such as central and peripheral compartment amounts. Its dimension need not equal the number of assay measurements. The matrix $F_k$ describes how the state at the preceding time contributes to the current state. The covariance $Q_k$ describes random state changes not explained by that transition. In contrast, $H_k$ maps the current state to the scale and components of the observation, and $R_k$ describes residual assay variation after that mapping.

A state-space model is not fully specified by writing down these matrices without context. Their dimensions and units must agree. If $x_k$ has $r$ components and $y_k$ has $m$, then $F_k$ is $r\times r$, $Q_k$ is $r\times r$, $H_k$ is $m\times r$, and $R_k$ is $m\times m$. A covariance entry has units equal to the product of the units of the corresponding variables. A covariance between two state components need not be zero.

From a Linear Diffusion to a Discrete Transition

A useful connection to the earlier material on exact densities for time-inhomogeneous linear diffusions is the continuous-time linear model $$ dX(t)=\{AX(t)+c\}\,dt+L\,dW(t), $$ where $A,c,L$ are fixed matrices or vectors of compatible dimensions and $W(t)$ is a standard Brownian motion of suitable dimension. Over an interval of length $\Delta=t_k-t_{k-1}>0$, its exact transition has the discrete state-space form. No time-step approximation is needed for this linear model.

Theorem: If $X(t_{k-1})=x$ is fixed, then the solution of the linear diffusion above satisfies $$ X(t_k)\mid X(t_{k-1})=x\sim N\!\left(Fx+a,Q\right), $$ where $$ F=e^{A\Delta},\qquad a=\int_0^\Delta e^{A(\Delta-s)}c\,ds,\qquad Q=\int_0^\Delta e^{A(\Delta-s)}LL^\mathsf{T}e^{A^\mathsf{T}(\Delta-s)}\,ds. $$ The integral defining $Q$ is positive semidefinite.

Proof. Multiply the stochastic differential equation by the integrating factor $e^{-At}$ and apply the product rule. Integrating from $t_{k-1}$ to $t_k$ gives the variation-of-constants expression $$ X(t_k)=e^{A\Delta}x+\int_0^\Delta e^{A(\Delta-s)}c\,ds+ \int_0^\Delta e^{A(\Delta-s)}L\,d\widetilde W(s), $$ where $\widetilde W(s)=W(t_{k-1}+s)-W(t_{k-1})$ is Brownian motion on the interval. The first two terms are fixed given $x$. The last term is a stochastic integral with deterministic integrand, so it is a mean-zero Gaussian vector. By the Itô isometry for deterministic integrands, its covariance is the stated integral $Q$. For any vector $z$, $$ z^\mathsf{T}Qz=\int_0^\Delta \left\|L^\mathsf{T}e^{A^\mathsf{T}(\Delta-s)}z\right\|^2\,ds\geq 0, $$ so $Q$ is positive semidefinite. Adding the deterministic terms gives the claimed Gaussian transition. $\square$

This derivation identifies the process-noise covariance as accumulated diffusion over the interval, not as the assay-error covariance. In particular, $Q$ depends on the interval length and on the diffusion loading $L$. If the continuous model is time-varying, corresponding transition and covariance integrals use its time-varying fundamental solution; the constant-coefficient formulas above should not be applied unchanged.

Worked Example: A Scalar Linear Concentration Process

Let the latent concentration $C(t)$, measured in mg/L, satisfy $$ dC(t)=-0.2C(t)\,dt+0.1\,dW(t). $$ Suppose observations are one day apart and the assay reports $Y_k=C(t_k)+\varepsilon_k$, with independent $\varepsilon_k\sim N(0,0.04)$ in squared mg/L.

For the transition, $A=-0.2$, $c=0$, $L=0.1$, and $\Delta=1$ day. Thus $F=e^{-0.2}=0.81873$, and the offset is $a=0$. The process variance is $$ Q=0.1^2\int_0^1 e^{-0.4(1-s)}\,ds =0.01\frac{1-e^{-0.4}}{0.4} \approx 0.008242. $$ Therefore $C(t_k)\mid C(t_{k-1})=c\sim N(0.81873c,0.008242)$. The observation equation has $H=1$, $d=0$, and $R=0.04$. The numbers $Q$ and $R$ are both variances in squared mg/L, but they describe different sources of variation: one-day latent evolution and assay error, respectively.

What the Observation Equation Means

The observation equation states that, conditional on the current state, the expected observation is $H_kx_k+d_k$, and its residual covariance is $R_k$. In particular, $$ y_k\mid x_k\sim N(H_kx_k+d_k,R_k). $$ This conditional statement is the measurement model; it does not say that observations are independent after the latent state has been marginalized out. Observations at different times can be associated because their latent states are associated through the transition model.

Result: If $x_k\sim N(m_k,P_k)$ and $v_k\sim N(0,R_k)$ is independent of $x_k$, then for $y_k=H_kx_k+d_k+v_k$, $$ E[y_k]=H_km_k+d_k,\qquad \operatorname{Var}(y_k)=H_kP_kH_k^\mathsf{T}+R_k. $$ Moreover, $y_k$ is Gaussian.

Derivation. The vector $(x_k^\mathsf{T},v_k^\mathsf{T})^\mathsf{T}$ is jointly Gaussian because its two components are independent Gaussian vectors. The observation is an affine transformation of that vector, so it is Gaussian. Linearity of expectation gives $E[y_k]=H_kE[x_k]+d_k+E[v_k]=H_km_k+d_k$. For the variance, independence makes the cross-covariance terms zero, hence $$ \operatorname{Var}(y_k)=H_kP_kH_k^\mathsf{T}+ \operatorname{Var}(v_k)=H_kP_kH_k^\mathsf{T}+R_k. $$ This separates variation in the state as seen through the assay from residual measurement variation. $\square$

Worked Example: Observing Central and Peripheral Amounts

Let $x_k=(A_{c,k},A_{p,k})^\mathsf{T}$ contain central and peripheral amounts in mg. Consider the transition and observation model $$ F=\begin{pmatrix}0.8&0.1\\0.2&0.7\end{pmatrix},\quad Q=\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix},\quad H=\begin{pmatrix}0.1&0\end{pmatrix},\quad R=(0.09). $$ The observation is central concentration in mg/L, with central volume $10$ L, so $H$ converts central amount to concentration and does not directly observe peripheral amount.

The dimensions are consistent: $F,Q$ are $2\times2$, $H$ is $1\times2$, and $R$ is $1\times1$. The process covariance is valid: its diagonal entries are positive and its determinant is $0.04(0.09)-0.01^2=0.0035>0$. The observation variance $0.09$ corresponds to assay-error standard deviation $0.3$ mg/L.

For a state mean $m=(100,40)^\mathsf{T}$ mg, the transition mean is $$ Fm=\begin{pmatrix}0.8(100)+0.1(40)\\0.2(100)+0.7(40)\end{pmatrix} =\begin{pmatrix}84\\48\end{pmatrix}\text{ mg}. $$ The expected observation at this state mean is $Hm=0.1(100)=10$ mg/L. If the current state covariance is $P=\begin{pmatrix}4&1\\1&9\end{pmatrix}$ in squared mg, the state contribution to observation variance is $$ HPH^\mathsf{T}=0.1^2(4)=0.04\;(\text{mg/L})^2. $$ Thus the marginal observation variance is $0.04+0.09=0.13\;(\text{mg/L})^2$. The $0.04$ is uncertainty in the central amount mapped through volume; the $0.09$ is assay error.

Worked Example: Replicate Assays and the Error Term

Suppose three replicate assay readings of the same stable sample are $12.1,\ 11.7,\ 12.0$ mg/L. Under the model $Y_j=x+\varepsilon_j$, with independent errors $\varepsilon_j\sim N(0,R)$, the sample mean is $$ \bar y=(12.1+11.7+12.0)/3=11.9333\text{ mg/L}. $$ The deviations from this mean are $0.1667,-0.2333,0.0667$ mg/L. Their squared sum is approximately $0.08667\;(\text{mg/L})^2$, so the unbiased sample variance is $0.08667/(3-1)=0.04333\;(\text{mg/L})^2$. This gives an illustrative estimate of replicate assay variance $R$, not an estimate of process variance $Q$: the same underlying sample was measured repeatedly rather than allowed to evolve.

If $R=0.04333\;(\text{mg/L})^2$, the variance of the mean of three independent replicates is $R/3\approx0.01444\;(\text{mg/L})^2$, and its standard error is $\sqrt{0.01444}\approx0.1202$ mg/L. This calculation relies on independent, equal-variance replicate errors and a state that is effectively unchanged during the replicate measurements.

Assumptions, Diagnostics, and Alternatives

Before applying a linear Gaussian model, check the assumptions against the PK design and data. The transition is linear in the chosen state, and the process noise is additive Gaussian with the specified covariance. The observation is linear in that state, with additive Gaussian assay error. Conditional on the state, the measurement noise has mean zero and covariance $R_k$. If multiple assays are recorded together, off-diagonal entries of $R_k$ allow correlated errors; setting them to zero is an independence assumption, not a requirement of the model.

Several practical checks can reveal a mismatch. Covariance matrices must be symmetric and positive semidefinite; a negative variance or a nonzero asymmetry signals an invalid specification or numerical problem. Predicted concentrations or amounts that are frequently negative may indicate that an additive Gaussian state model is unsuitable for a quantity constrained to be nonnegative. Strongly skewed or concentration-dependent assay residuals can contradict constant-variance additive Gaussian error. Residual patterns over time can also signal an inadequate transition model, such as unrepresented changes in the latent dynamics. A good-looking fit alone does not establish that process noise and assay error have been separated correctly.

The linear Gaussian model is attractive because its state and observation distributions are described by means and covariances. Its main alternative depends on what assumption fails. A nonlinear state or observation relationship may motivate an extended Kalman filter or another nonlinear method; non-Gaussian transitions or observations may require a different filtering distribution. These alternatives do not automatically fix a poor biological model or an incorrect error structure. The linear Gaussian model is also not a substitute for a positivity-preserving process when nonnegative concentrations are essential.

Worked Example: A Negative Gaussian Concentration Is a Warning

Suppose a model assigns a latent concentration the distribution $C\sim N(0.10,0.20^2)$ mg/L. Its probability of being negative is $$ P(C<0)=\Phi\!\left(\frac{0-0.10}{0.20}\right) =\Phi(-0.5)\approx0.3085. $$ Thus this model assigns about $31\%$ probability to negative concentrations. A Gaussian approximation may still be used in a setting where that probability is negligible and the approximation is appropriate, but here the negative mass is substantial. It is a concrete warning to reconsider the state scale or model family, rather than to interpret the model as a physically positive concentration process.

Key takeaway: A filtering model needs a defined latent state, a transition law with process covariance $Q_k$, and an observation mapping with assay covariance $R_k$. Keep the units, dimensions, and scientific meanings of these quantities separate; $Q_k$ is not a substitute for $R_k$.

Check Your Understanding

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

A two-component state is observed by one assay. What are the dimensions of $F_k,Q_k,H_k,R_k$?
In a model with $y_k=H_kx_k+v_k$, which covariance describes replicate assay variability at an effectively fixed state?
For the scalar diffusion $dX=-0.1X\,dt+0.2\,dW$ over $\Delta=2$, write the exact $F$ and $Q$ formulas.
If $x_k\sim N(m_k,P_k)$, state the mean and variance of $y_k=H_kx_k+d_k+v_k$ when $v_k\sim N(0,R_k)$ is independent.

Solutions. (1) $F_k$ and $Q_k$ are $2\times2$, $H_k$ is $1\times2$, and $R_k$ is $1\times1$. (2) $R_k$, because replicate readings at a fixed state vary through the assay-error term. (3) $F=e^{-0.1(2)}=e^{-0.2}$ and $Q=0.2^2\int_0^2 e^{-0.2(2-s)}ds=0.04(1-e^{-0.4})/0.2\approx0.06594$. (4) $E[y_k]=H_km_k+d_k$ and $\operatorname{Var}(y_k)=H_kP_kH_k^\mathsf{T}+R_k$.