Convergence Rates to Periodic Steady State

Why convergence rate matters

A periodic steady-state law describes the distribution of a pharmacokinetic state at a fixed phase of a repeated dosing regimen. It does not, by itself, say how many doses are needed for a patient population starting from a different distribution to be close to that law. The rate depends on how strongly the dynamics retain information about the state at the start of each dosing interval. For a linear model, that dependence is summarized by a scalar transition multiplier or, for several compartments, by the eigenvalues of a cycle-transition matrix.

The rates below concern distributions observed at the same point in successive dosing cycles. A continuous-time path can behave differently within a cycle, and a small difference at dose times should not be confused with an identical difference at every time between doses.

Cycle transitions and their assumptions

Fix a dosing interval of length $\tau$, and observe the state just after each dose. For a linear diffusion with deterministic periodic coefficients, additive noise, and fixed affine dose maps, the state over one full cycle can be written as an affine random transition. Independent Brownian increments in different cycles give independent, identically distributed cycle-noise terms. In a $d$-dimensional model, write this transition as

X_{n+1}=M X_n+b+\eta_{n+1},\qquad n=0,1,2,\ldots,

Here $X_n\in\mathbb{R}^d$ is the state at the chosen dose-time phase, $M$ is the cycle-transition matrix, $b$ is the accumulated deterministic input over a cycle, and the $\eta_n$ are independent $N(0,Q)$ random vectors, independent of $X_0$. The covariance matrix $Q$ is allowed to be positive semidefinite. For a scalar model, $M$ is the multiplier $a$.

This representation applies when the between-dose state equation is linear with additive Gaussian noise and the dose-time maps are affine and repeat identically each cycle. The matrix $M$ includes both between-dose propagation and any linear parts of the dose maps. In a time-homogeneous linear SDE between doses, the propagation over an interval of length $\tau$ is a matrix exponential; the cycle map must also include the dose map. The analysis here does not assume that the state is nonnegative.

Definition: The cycle-boundary convergence rate describes how the law of $X_n$, observed at the same phase after $n$ cycles, approaches an invariant law for the cycle transition. For a matrix transition, the spectral radius is $\rho(M)=\max_i|\lambda_i|$, where the $\lambda_i$ are the eigenvalues of $M$.

Scalar convergence: the multiplier controls the gap

First consider a deterministic scalar recurrence $Z_{n+1}=aZ_n+c$, with a constant input $c$. If $|a|<1$, its fixed point is $z_*=c/(1-a)$. Subtracting the fixed-point equation from the recurrence gives $Z_{n+1}-z_*=a(Z_n-z_*)$. Repeating this identity yields the exact error after $n$ cycles:

Z_n-z_*=a^n(Z_0-z_*),\qquad |Z_n-z_*|=|a|^n|Z_0-z_*|.

Worked Example: A Dose-Time Concentration Gap

Suppose the dose-time concentration follows $Z_{n+1}=e^{-1.4}Z_n+2$, where concentrations are in mg/L and the cycle multiplier is $a=e^{-1.4}\approx0.246597$. Find the fixed point and the first two dose-time concentrations when $Z_0=0$.

The fixed point solves $z_*=az_*+2$, so $z_*=2/(1-a)=2/(1-e^{-1.4})\approx2.6546$ mg/L. The exact gap after $n$ cycles is $a^n z_*$, since $Z_0-z_*=-z_*$.

After one cycle, $Z_1=2$ mg/L, so the gap is $z_*-Z_1\approx0.6546$ mg/L. After two cycles, $Z_2=2(1+a)\approx2.4932$ mg/L, and the gap is $z_*-Z_2\approx0.1614$ mg/L. This agrees with $a^2z_*\approx0.1614$ mg/L. Each additional cycle multiplies the remaining gap by approximately $0.2466$.

A distributional rate for linear Gaussian cycles

For the stochastic recurrence, the state difference under two paths driven by the same cycle noise is especially informative. The noise cancels, leaving the deterministic difference $M$ times its previous value. To express a distributional bound, use the 2-Wasserstein distance. For probability distributions with finite second moments, $W_2(\mu,\nu)$ is the infimum, over all joint distributions with marginals $\mu$ and $\nu$, of the square root of $E[\|U-V\|^2]$.

Theorem: Suppose $X_{n+1}=MX_n+b+\eta_{n+1}$, where the $\eta_n$ are independent $N(0,Q)$ vectors, independent of $X_0$, and $M,b,Q$ are fixed. If $\rho(M)<1$, there is a unique invariant law with finite second moment. For any initial law with finite second moment, if $\pi$ denotes that invariant law, then $$ W_2\bigl(\mathcal{L}(X_n),\pi\bigr) \leq \|M^n\|_2\,W_2\bigl(\mathcal{L}(X_0),\pi\bigr). $$ Here $\|\cdot\|_2$ is the matrix operator norm induced by the Euclidean norm.

Proof. Since $\rho(M)<1$, the series of random vectors \[ X_*=\sum_{j=0}^{\infty}M^j(b+\eta_{-j}) \] converges in mean square when the noise is indexed by independent copies in the past. Indeed, $\|M^j\|_2$ is bounded by a constant times a polynomial in $j$ times $r^j$ for some $r<1$, as follows from the Jordan form of $M$. Thus both the deterministic series and the sum of the independent centered noise terms converge in mean square. Shifting the index in this series shows that $X_*$ has an invariant law for the recurrence. Its mean and covariance are finite.

For uniqueness among invariant laws with finite second moments, couple any two invariant initial states $U_0,V_0$ with finite second moments and drive them using the same future noises. Then $U_n-V_n=M^n(U_0-V_0)$. Because $\|M^n\|_2\to0$, the mean-square distance between the coupled states tends to zero. If both initial laws are invariant, their laws remain unchanged at every $n$, so their $W_2$ distance must be zero; the laws are identical.

For the bound, take a coupling of $X_0$ and an independent-of-future-noise $X_*$ that approaches the infimum defining $W_2(\mathcal{L}(X_0),\pi)$. Evolve both states with the same future noises. Their difference at step $n$ is $M^n(X_0-X_*)$; hence \[ E[\|X_n-X_{*,n}\|^2]\leq \|M^n\|_2^2 E[\|X_0-X_*\|^2]. \] Taking square roots and then the infimum over initial couplings gives the stated bound. $\square$

For this Gaussian recurrence, the invariant law is Gaussian. Its mean $\mu_*$ and covariance $\Sigma_*$ satisfy

\mu_*=M\mu_*+b,\qquad \Sigma_*=M\Sigma_*M^{\mathsf T}+Q.

The theorem gives a distance bound, not necessarily equality. It also shows why $\rho(M)<1$ is the stability condition: the powers of $M$ decay. The precise asymptotic exponential rate is determined by the largest eigenvalue magnitude. If $M$ is diagonalizable, its powers are bounded by a constant times $\rho(M)^n$ when $\rho(M)>0$. If $M$ has a Jordan block of size $k$ associated with an eigenvalue of magnitude $\rho(M)>0$, powers can instead contain a factor of order $n^{k-1}\rho(M)^n$. In either case, the $n$-th root of the power norm approaches $\rho(M)$.

Worked Example: A Stochastic Scalar Cycle

Consider $Z_{n+1}=aZ_n+0.4+\eta_{n+1}$, where $a=e^{-1}$, the independent cycle noises have distribution $N(0,0.09)$, and $E[Z_0]=0$, $\operatorname{Var}(Z_0)=0.16$. Calculate the stationary mean and variance, then the mean and variance after one cycle.

The stationary mean is $\mu_*=0.4/(1-a)=0.4/(1-e^{-1})\approx0.6328$. The stationary variance solves $v_*=a^2v_*+0.09$, giving $v_*=0.09/(1-e^{-2})\approx0.10409$.

The one-cycle mean is $E[Z_1]=a(0)+0.4=0.4$. Its difference from stationarity is $0.4-\mu_*\approx-0.2328$, which is also $a(0-\mu_*)$. The one-cycle variance is \[ \operatorname{Var}(Z_1)=a^2(0.16)+0.09 =e^{-2}(0.16)+0.09\approx0.11165. \] The variance difference from stationarity is $a^2(0.16-v_*)\approx0.00757$. Thus the mean gap is multiplied by $a$ per cycle, while the variance gap is multiplied by $a^2$. In this example the slowest distributional component is the mean, which contracts by approximately $0.3679$ each cycle.

Several compartments: eigenvalues and transient behavior

For multiple compartments, each eigenvalue describes a mode of the cycle transition. The mode with the largest magnitude generally determines the long-run exponential rate, but eigenvalues alone may conceal short-run behavior. In particular, a matrix that is not normal can temporarily amplify some errors even when all its eigenvalues lie inside the unit circle.

Worked Example: Two Decoupled Cycle Modes

Let $M=\begin{pmatrix}0.8&0\\0&0.3\end{pmatrix}$, and suppose the cycle noise has covariance $Q=\operatorname{diag}(0.04,0.01)$. The stationary covariance is diagonal, with entries $0.04/(1-0.8^2)=0.1111$ and $0.01/(1-0.3^2)\approx0.0110$.

Start with a Gaussian law having this stationary covariance but a mean displaced from the stationary mean by $(1,2)^{\mathsf T}$. After $n$ cycles, the covariance remains stationary, while the mean displacement is \[ M^n(1,2)^{\mathsf T}=(0.8^n,2(0.3)^n)^{\mathsf T}. \] After three cycles it is $(0.512,0.054)^{\mathsf T}$. Since the covariance is unchanged, the $W_2$ distance is exactly the Euclidean length of this displacement: \[ \sqrt{0.512^2+0.054^2}\approx0.5148. \] For large $n$, the $0.8^n$ mode dominates, so the asymptotic rate is $0.8$, not $0.3$. The faster mode can still matter during the first few cycles.

A defective matrix can produce an additional polynomial factor. For example, for $M=\begin{pmatrix}0.8&1\\0&0.8\end{pmatrix}$, direct multiplication gives \[ M^n=\begin{pmatrix}0.8^n&n(0.8)^{n-1}\\0&0.8^n\end{pmatrix}. \] The eigenvalue magnitude is $0.8$, but the off-diagonal term initially grows before eventually decaying. Consequently, the spectral radius describes the asymptotic exponential rate, not necessarily the observed reduction from one early cycle to the next.

Rates at other phases, diagnostics, and limitations

To assess convergence at a time $s$ within a dosing interval, write the within-cycle transition from the dose-time state as $X(n\tau+s)=\Phi(s)X_n+c(s)+\xi_s$, where $\xi_s$ is the within-cycle noise. If two paths use the same within-cycle noise, their difference is $\Phi(s)$ times their dose-time difference. Thus their $W_2$ distance at phase $s$ is bounded by $\|\Phi(s)\|_2$ times the boundary distance. The phase-specific multiplier affects the size of the bound, while the cycle-to-cycle asymptotic rate still comes from $M$.

A useful diagnostic is to compare observed dose-time discrepancies over successive cycles with the predicted powers of $M$. For a scalar model, a nonzero discrepancy should change by the factor $|a|$ under the same deterministic transition. For a vector model, examine mode-specific errors or the singular values of $M^n$, not just individual eigenvalues. A discrepancy that grows initially does not by itself contradict eventual stability; it may reflect nonnormal transient amplification. Persistent failure to approach the predicted periodic law can instead indicate a misspecified cycle map, changing doses or parameters, non-repeating noise assumptions, or an unstable mode with magnitude at least one.

These rates have important limits. The $W_2$ theorem requires finite second moments and the stated affine, additive-noise recurrence; it does not establish the same rate in total variation or for arbitrary nonlinear, state-dependent diffusion models. If $\rho(M)\geq1$, the contraction argument fails, and an invariant law with finite second moment need not exist. Even with $\rho(M)<1$, a small spectral radius does not guarantee that every early cycle is close to steady state, particularly for a nonnormal matrix or a distant initial distribution.

The main alternatives answer different questions. A scalar multiplier gives an exact error ratio for a deterministic one-dimensional recurrence. Matrix powers give the detailed cycle-by-cycle change for a multicomponent linear model. Simulation can estimate convergence for a more complicated model, but its apparent rate can be obscured by Monte Carlo error; it should be checked against analytic transition structure where available. None of these cycle-boundary summaries replaces evaluation of the distribution at the specific within-cycle sampling times relevant to a study.

Key takeaway: For a stable affine linear cycle transition, the matrix powers $M^n$ quantify how quickly the initial-state effect fades. The spectral radius gives the asymptotic exponential rate, while the full matrix powers reveal transient amplification and phase-specific differences.

Check Your Understanding

Try each question, then check the solutions below.

A scalar recurrence has multiplier $a=0.5$. By what factor does its deterministic gap from the fixed point change after four cycles?
For a cycle matrix with eigenvalues $0.6$ and $-0.9$, what is the asymptotic exponential rate, and does the second mode alternate in sign?
Why can a matrix with all eigenvalue magnitudes below one still temporarily enlarge an error?
In the scalar stochastic recurrence $Z_{n+1}=0.4Z_n+1+\eta_{n+1}$, with independent noise of variance $0.12$, find the stationary mean and variance.

Solutions. (1) The factor is $a^4=0.5^4=0.0625$. (2) The spectral radius is $0.9$, so the asymptotic exponential rate is $0.9$; the eigenvalue $-0.9$ makes its mode alternate in sign. (3) Eigenvalues determine asymptotic behavior, but a nonnormal matrix can have powers whose operator norm initially grows; inspect $M^n$ or its singular values to diagnose this. (4) The stationary mean is $1/(1-0.4)=5/3$. The stationary variance is $0.12/(1-0.4^2)=0.12/0.84=1/7\approx0.1429$.