Matrix-Exponential Solutions for Linear Means

Why the Matrix Exponential Gives the Mean

A linear pharmacokinetic model can have several interacting states, such as amounts in a depot and a central compartment. When its drift is affine, the mean vector follows a deterministic linear differential equation. The matrix exponential solves that equation while retaining the effects of transfers between states. A constant infusion contributes an additional integral, which remains well-defined even when the drift matrix is singular.

The key distinction is between the stochastic state and its mean. The noise can change the variability around the mean without changing the mean trajectory, provided the stochastic integral has expectation zero and the drift is affine. The earlier material on affine linear SDEs for pharmacokinetic states supplies the model form; here the goal is to obtain and use its closed-form mean.

Definition: Let $X_t\in\mathbb{R}^d$ satisfy the Itô SDE $dX_t=(AX_t+b)\,dt+G_t\,dW_t$, where $A$ is a constant $d\times d$ matrix, $b\in\mathbb{R}^d$ is constant and deterministic, $W_t$ is a Brownian motion, and $G_t$ is a predictable matrix-valued process. Assume $E[\|X_0\|]<\infty$, $E[\int_0^T\|G_s\|_F^2\,ds]<\infty$, and $E[\int_0^T\|X_s\|\,ds]<\infty$ for every finite $T$. The mean vector is $m(t)=E[X_t]$.

These assumptions make the stochastic integral a square-integrable martingale with expectation zero, and ensure that expectations can be taken in the drift integral. A standard sufficient setting is a linear SDE with a finite-second-moment initial state and diffusion coefficients that ensure finite second moments on finite time intervals. The diffusion need not be constant, as long as the stated integrability conditions hold.

From the SDE to a Closed-Form Mean

Taking expectations in the integral form of the SDE gives $m(t)=m(0)+\int_0^t(Am(s)+b)\,ds$. Thus the mean solves a deterministic initial-value problem. The following result gives its solution without requiring $A$ to be invertible.

Theorem: Under the definition's assumptions, let $m_0=E[X_0]$. Then, for $t\geq 0$, $$ m(t)=e^{At}m_0+\int_0^t e^{A(t-s)}b\,ds. $$ This is the unique solution to $m'(t)=Am(t)+b$, $m(0)=m_0$.

Proof. The matrix exponential is defined by the convergent power series $e^{At}=\sum_{k=0}^{\infty}(At)^k/k!$, and satisfies $\frac{d}{dt}e^{At}=Ae^{At}=e^{At}A$. Define $q(t)$ by the displayed formula. Differentiating its integral, whose integrand is continuous, gives

q'(t)=Ae^{At}m_0+b+\int_0^t Ae^{A(t-s)}b\,ds =Aq(t)+b.

Also $q(0)=m_0$. To verify uniqueness, suppose $q_1$ and $q_2$ solve this initial-value problem and let $r=q_1-q_2$. Then $r'=Ar$, $r(0)=0$. Differentiating $e^{-At}r(t)$ gives $e^{-At}(r'-Ar)=0$, so $e^{-At}r(t)$ is constant and equals zero at $t=0$. Hence $r(t)=0$ for all $t$. The expectation equation obtained from the SDE has this same initial condition, so its mean is $q(t)$. $\square$

The integral is a convolution of the input with the system's matrix exponential. It accounts for the input at every time $s$, then propagates that contribution from $s$ to $t$. If $A$ is invertible, one may also write

m(t)=e^{At}m_0+A^{-1}(e^{At}-I)b.

This follows because $A\int_0^t e^{A(t-s)}\,ds=e^{At}-I$. It is a convenient shortcut in nonsingular cases, not the general solution: the integral formula remains valid when $A$ is singular.

Worked Example: One-Compartment Constant Infusion

Let $X_t$ be the amount in a one-compartment model with elimination rate $k=0.2$ per hour, constant infusion rate $R=10$ amount units per hour, and additive process noise: $dX_t=(-0.2X_t+10)\,dt+\sigma\,dW_t$. Take $X_0=5$ with probability one. Find the mean after three hours.

Here $A=-0.2$, $b=10$, and $m_0=5$. Since $A$ is a nonzero scalar, the inverse form gives

m(t)=e^{-0.2t}(5)+\frac{e^{-0.2t}-1}{-0.2}(10) =50-45e^{-0.2t}.

At $t=3$, $e^{-0.6}\approx0.548812$, so $m(3)=50-45(0.548812)\approx25.3035$ amount units. The steady value is $R/k=50$, and the mean has moved from 5 toward 50, as expected. The value of $\sigma$ does not appear in this calculation: under the stated additive-noise model it affects variability, not the mean.

Singular Matrices and Constant Inputs

In a model with a singular $A$, the expression $A^{-1}(e^{At}-I)b$ is undefined. This does not mean the mean lacks a solution; it means that the inverse shortcut is unavailable. Two reliable approaches are to keep the integral in the theorem or to include the constant input in an augmented matrix exponential.

Theorem: For constant $A$ and $b$, define $z(t)=\begin{pmatrix}m(t)\\1\end{pmatrix}$ and $C=\begin{pmatrix}A&b\\0&0\end{pmatrix}$. Then $$ z(t)=e^{Ct}z(0). $$ In particular, the upper-right block of $e^{Ct}$ is $\int_0^t e^{A(t-s)}b\,ds$, whether or not $A$ is invertible.

Proof. The mean equation and the constant last coordinate give $z'=Cz$. The matrix-exponential solution to this constant-coefficient initial-value problem is $z(t)=e^{Ct}z(0)$, as follows by differentiating the power series for $e^{Ct}$. To identify the input contribution, write $e^{Ct}$ in blocks as $\begin{pmatrix}P(t)&u(t)\\0&1\end{pmatrix}$. The equation $\frac{d}{dt}e^{Ct}=Ce^{Ct}$ implies $u'=Au+b$, $u(0)=0$. By the preceding theorem's uniqueness argument, $u(t)=\int_0^t e^{A(t-s)}b\,ds$. This proof does not use $A^{-1}$. $\square$

Worked Example: Singular Two-State Model with Infusion

Consider the mean dynamics $m'=Am+b$, with $A=\begin{pmatrix}-1&0\\1&0\end{pmatrix}$, $b=\begin{pmatrix}4\\0\end{pmatrix}$, and $m(0)=\begin{pmatrix}0\\0\end{pmatrix}$. Interpret the first state as a depot amount that transfers to a second state; the second state has no elimination in this simplified example. Since $\det(A)=0$, an inverse formula cannot be used.

The matrix is lower triangular and its exponential is

e^{At}= \begin{pmatrix} e^{-t}&0\\ 1-e^{-t}&1 \end{pmatrix}.

Multiplying by the input vector gives $e^{A(t-s)}b=\begin{pmatrix}4e^{-(t-s)}\\4(1-e^{-(t-s)})\end{pmatrix}$. Integrating from $s=0$ to $s=t$, or equivalently integrating over $u=t-s$, yields

m_1(t)=\int_0^t4e^{-u}\,du=4(1-e^{-t}),\qquad m_2(t)=\int_0^t4(1-e^{-u})\,du=4(t-1+e^{-t}).

At $t=2$, $e^{-2}\approx0.135335$, so $m_1(2)\approx3.45866$ and $m_2(2)\approx4.54134$. Their sum is 8, matching the total infused amount $4\times2=8$: in this model, transfer changes the state allocation but neither state eliminates drug.

Worked Example: Infusion in a Two-Compartment Linear Model

Take $A=\begin{pmatrix}-1&0\\1&-2\end{pmatrix}$, $b=\begin{pmatrix}8\\0\end{pmatrix}$, and $m(0)=\begin{pmatrix}0\\0\end{pmatrix}$. The input is constant in the first compartment. The component equations are $m_1'=-m_1+8$ and $m_2'=m_1-2m_2$.

Solving the first equation gives $m_1(t)=8(1-e^{-t})$. For the second, use its integrating factor $e^{2t}$:

m_2(t)=\int_0^t e^{-2(t-s)}\,8(1-e^{-s})\,ds =8\left(\frac{1-e^{-2t}}{2}-(e^{-t}-e^{-2t})\right) =4(1-e^{-t})^2.

At $t=\log 2$, $e^{-t}=1/2$, giving $m_1=8(1/2)=4$ and $m_2=4(1/2)^2=1$. The model's steady mean solves $Am+b=0$, which gives $m_1=8$, $m_2=4$; the values at $\log 2$ are below those steady values. This componentwise calculation is also a check of the matrix-exponential solution.

Choosing and Checking a Calculation

The convolution formula is the general analytic answer. The augmented matrix exponential is useful for numerical evaluation, especially when the matrix is singular or when a numerical library reliably computes exponentials of block matrices. Componentwise integration can be simpler when the equations are triangular, as in the two-state example. Alternatively, a numerical ODE solver can integrate $m'=Am+b$ directly; this is useful for checking an analytic formula or handling a time-varying input, but it is not itself a closed-form solution.

A few checks catch many errors. At $t=0$, the formula must return $m_0$. Its derivative at zero must equal $Am_0+b$. Units must also agree: $A$ has units of inverse time, $b$ has state units per time, and the integral contributes state units. In a stable model with invertible $A$, a constant-input steady mean must satisfy $Am_\infty+b=0$. If a computed trajectory moves away from a stable equilibrium under constant infusion, check the sign of the input term and the exponent.

Common pitfall: Writing $A^{-1}(e^{At}-I)b$ without checking whether $A$ is invertible. A zero eigenvalue makes this expression invalid, even though the mean trajectory may be finite and simple. Use the integral or augmented-exponential formula instead. Also check the order of matrix factors: in general, matrices cannot be rearranged as if they were scalars.

The matrix-exponential approach has limits. It relies on constant $A$ and $b$; a time-varying drift matrix generally requires a time-ordered fundamental solution rather than $e^{At}$. And a closed-form mean does not give the distribution or variability of the state. It answers a specific question: how the expectation evolves under affine drift and an integrable stochastic term.

Check Your Understanding

Try each calculation, then compare with the solutions below.

For $m'=-0.5m+3$, $m(0)=2$, find $m(t)$ and $m(2)$.
Why is the integral formula valid when $A$ is singular?
For $m'=Am+b$, what condition must a constant steady mean satisfy?
In the singular infusion example, explain why the sum of the two mean states at time $t$ equals $4t$.

Solutions. 1. $m(t)=6+(2-6)e^{-0.5t}=6-4e^{-0.5t}$, so $m(2)=6-4e^{-1}\approx4.5285$. 2. The matrix exponential and its integral are defined without an inverse; only the alternative inverse-based expression requires nonsingular $A$. 3. It must satisfy $Am_\infty+b=0$; existence and uniqueness of such a vector require additional conditions on $A$. 4. The displayed means sum to $4(1-e^{-t})+4(t-1+e^{-t})=4t$, because the input adds mass at rate 4 and there is no elimination.

Matrix-Exponential Solutions for Linear Means

What You'll Learn