Forward Kolmogorov Equation for Pharmacokinetic Diffusions

From a Pharmacokinetic SDE to a Density Equation

A stochastic differential equation describes the evolution of a random pharmacokinetic state along a path. The forward Kolmogorov equation instead describes how the probability distribution of that state changes over time. Its coefficients come directly from the SDE: the drift determines the local directional change, while the diffusion determines how random variation spreads probability across states.

This tutorial uses Itô SDEs, consistent with the earlier material on stochastic calculus and pharmacokinetic SDE notation. The central task is to map the model’s coefficients to the probability density equation without confusing a diffusion amplitude with a variance rate. The equation is also called the Fokker-Planck equation.

Definition: Let \(X_t\) be a \(d\)-dimensional state satisfying the Itô SDE \(dX_t=b(t,X_t)\,dt+\sigma(t,X_t)\,dW_t\), where \(W_t\) is an \(m\)-dimensional standard Brownian motion, \(b\) is a \(d\)-vector, and \(\sigma\) is a \(d\times m\) matrix. The diffusion matrix is \(a(t,x)=\sigma(t,x)\sigma(t,x)^{\mathsf T}\). If \(X_t\) has density \(p(t,x)\) with respect to \(d\)-dimensional Lebesgue measure, its forward Kolmogorov equation is

\frac{\partial p(t,x)}{\partial t} = -\sum_{i=1}^{d}\frac{\partial}{\partial x_i}\bigl(b_i(t,x)p(t,x)\bigr) +\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d} \frac{\partial^2}{\partial x_i\partial x_j}\bigl(a_{ij}(t,x)p(t,x)\bigr).

Here \(t\) is elapsed time and \(x=(x_1,\ldots,x_d)\) is a possible state, such as a vector of compartment amounts or concentrations. The density is a probability density over states, not a path or a concentration value.

The factor \(1/2\) and the matrix \(a=\sigma\sigma^{\mathsf T}\) are essential. In one dimension, if \(dX_t=b(t,X_t)\,dt+s(t,X_t)\,dW_t\), then \(a=s^2\), and the equation is \(\partial_t p=-\partial_x(bp)+\tfrac12\partial_{xx}(s^2p)\). With several Brownian drivers, it is the sum of their variance contributions, including cross-covariances, that forms \(a\). It is generally not correct to substitute \(\sigma\) itself for \(a\).

Why This Equation Follows From the SDE

A precise derivation begins with a smooth test function \(\varphi:\mathbb R^d\to\mathbb R\) that has compact support and two continuous derivatives. Assume the SDE has a nonexplosive solution, the relevant expectations exist, and Itô’s formula applies. Its generator at time \(t\) is \(L_t\varphi(x)=\sum_i b_i(t,x)\partial_i\varphi(x)+\tfrac12\sum_{i,j}a_{ij}(t,x)\partial_{ij}\varphi(x)\). Itô’s formula and the zero expectation of the stochastic integral give the weak evolution of expectations.

Theorem: Suppose \(X_t\) satisfies the Itô SDE above and has density \(p(t,x)\). Assume the coefficients and density are sufficiently regular for the displayed integrals and derivatives to exist, and that the stochastic integral in Itô’s formula is a true martingale for each test function used. Then, for every compactly supported twice continuously differentiable \(\varphi\),

\frac{d}{dt}\int_{\mathbb R^d}\varphi(x)p(t,x)\,dx = \int_{\mathbb R^d} \left[ \sum_i b_i(t,x)\partial_i\varphi(x) +\frac12\sum_{i,j}a_{ij}(t,x)\partial_{ij}\varphi(x) \right]p(t,x)\,dx.

If \(p\) is continuously differentiable in \(t\), twice continuously differentiable in \(x\), and the coefficients are regular enough for integration by parts, this identity implies the forward Kolmogorov equation pointwise.

To justify the theorem, apply Itô’s formula between times \(0\) and \(t\): \(d\varphi(X_t)=L_t\varphi(X_t)\,dt+\nabla\varphi(X_t)^{\mathsf T}\sigma(t,X_t)\,dW_t\). Taking expectations removes the stochastic integral by the martingale assumption. The fundamental theorem of calculus then gives \(\frac{d}{dt}E[\varphi(X_t)]=E[L_t\varphi(X_t)]\) for almost every \(t\), and, when the integrands are continuous in time, for every \(t\). Since \(X_t\) has density \(p(t,\cdot)\), replacing each expectation by its density integral gives the displayed weak identity.

For the classical equation, write the left side as \(\int\varphi\,\partial_t p\,dx\). On the right, integrate each drift term once by parts and each diffusion term twice. Compact support of \(\varphi\) removes the boundary terms. The result is

\int_{\mathbb R^d}\varphi(x)\,\partial_t p(t,x)\,dx = \int_{\mathbb R^d}\varphi(x) \left[ -\sum_i\partial_i(b_i p) +\frac12\sum_{i,j}\partial_i\partial_j(a_{ij}p) \right]dx.

Both sides are integrals against every such test function. Under the stated continuity conditions, their integrands must agree at every \(x\); otherwise their continuous difference would have one strict sign on some neighborhood, and a nonnegative test function supported there would give a nonzero integral. This proves the pointwise equation. If the density or coefficients are not sufficiently smooth, the weak identity remains the appropriate formulation; a classical pointwise equation has not thereby been established. \(\square\)

Worked Pharmacokinetic Translations

Worked Example: One-Compartment Concentration With Additive Noise

Consider a one-compartment concentration \(C_t\) under constant infusion: \(dC_t=(R/V-kC_t)\,dt+\sigma\,dW_t\). Let \(R/V=2\) concentration units per hour, \(k=0.5\) per hour, and \(\sigma=0.3\) concentration units per square root hour. Identify the equation for the density \(p(t,c)\).

The state is \(c\), the time variable is \(t\), and the drift is \(b(c)=2-0.5c\). The diffusion amplitude is \(\sigma=0.3\), so the variance rate is \(a=\sigma^2=0.09\). Because this is a scalar model, the forward equation is

\frac{\partial p(t,c)}{\partial t} = -\frac{\partial}{\partial c}\bigl((2-0.5c)p(t,c)\bigr) +\frac{1}{2}(0.09)\frac{\partial^2p(t,c)}{\partial c^2} = -\frac{\partial}{\partial c}\bigl((2-0.5c)p(t,c)\bigr) +0.045\frac{\partial^2p(t,c)}{\partial c^2}.

The density describes uncertainty over possible concentrations at a fixed time; it is not itself the mean concentration. If the initial concentration is known to equal \(c_0\), the initial law is a point mass at \(c_0\), rather than an ordinary density function at time zero. This initial condition does not change the form of the equation for later times when a density exists.

Worked Example: Two Compartments With Independent Process Noise

Let \(X_1\) and \(X_2\) be central and peripheral amounts, with independent Brownian drivers: \(dX_1=[5-0.8X_1+0.2X_2]dt+0.4dW_1\) and \(dX_2=[0.6X_1-0.2X_2]dt+0.3dW_2\). Find the forward equation for \(p(t,x_1,x_2)\).

The drift components are \(b_1=5-0.8x_1+0.2x_2\) and \(b_2=0.6x_1-0.2x_2\). The diffusion matrix \(\sigma\) is diagonal with entries \(0.4\) and \(0.3\). Independence of the drivers means the off-diagonal entries of \(a=\sigma\sigma^{\mathsf T}\) are zero. Its diagonal entries are \(0.4^2=0.16\) and \(0.3^2=0.09\). Therefore,

\begin{aligned} \frac{\partial p}{\partial t} ={}&-\frac{\partial}{\partial x_1} \bigl[(5-0.8x_1+0.2x_2)p\bigr] -\frac{\partial}{\partial x_2} \bigl[(0.6x_1-0.2x_2)p\bigr]\\ &+\frac12(0.16)\frac{\partial^2p}{\partial x_1^2} +\frac12(0.09)\frac{\partial^2p}{\partial x_2^2}\\ ={}&-\partial_{x_1}\bigl[(5-0.8x_1+0.2x_2)p\bigr] -\partial_{x_2}\bigl[(0.6x_1-0.2x_2)p\bigr] +0.08\partial_{x_1x_1}p+0.045\partial_{x_2x_2}p. \end{aligned}

The cross-compartment terms in the drift describe exchange and appear in the first-derivative terms. They do not create a mixed second derivative here, because the two Brownian drivers are independent and the diffusion matrix has \(a_{12}=0\). Correlated Brownian drivers would require checking that entry rather than assuming it is zero.

Worked Example: A Concentration-Dependent Diffusion Amplitude

On the interior \(c>0\), consider \(dC_t=-0.4C_t\,dt+0.2C_t\,dW_t\). Write the formal forward equation and identify the diffusion term.

The drift is \(b(c)=-0.4c\), while the diffusion amplitude is \(s(c)=0.2c\). Squaring the amplitude gives \(a(c)=0.04c^2\), not \(0.2c\). Substitution into the one-dimensional equation gives

\frac{\partial p(t,c)}{\partial t} = -\frac{\partial}{\partial c}\bigl((-0.4c)p(t,c)\bigr) +\frac12\frac{\partial^2}{\partial c^2}\bigl(0.04c^2p(t,c)\bigr) = 0.4\frac{\partial}{\partial c}(cp) +0.02\frac{\partial^2}{\partial c^2}(c^2p).

The coefficient \(c^2\) remains inside the second derivative. Moving it outside would omit derivative terms and change the equation. This formal interior calculation does not specify boundary behavior at zero or establish that a density exists; those questions require separate analysis.

Probability Mass, Diagnostics, and Alternatives

Theorem: Suppose \(p\) is a classical solution on \(\mathbb R^d\) and the probability flux decays sufficiently fast at infinity that the boundary integrals below vanish. Then its total mass is constant in time. In particular, if \(\int_{\mathbb R^d}p(0,x)\,dx=1\), then \(\int_{\mathbb R^d}p(t,x)\,dx=1\) for all times in the interval of classical existence.

Proof. Integrate the forward equation over \(\mathbb R^d\). Each drift term is an integral of a first derivative, \(-\int\partial_i(b_i p)\,dx\), and each diffusion term is an integral of a second derivative, \(\tfrac12\int\partial_i\partial_j(a_{ij}p)\,dx\). By the assumed decay, integration by parts leaves no boundary contributions, so every term is zero. Thus \(\frac{d}{dt}\int p(t,x)\,dx=0\), and the integral equals its initial value. \(\square\)

This conservation check is useful when solving or discretizing the equation, but it is not a substitute for checking the model assumptions. For a bounded or half-line state space, boundary fluxes matter: probability can leave or enter through a boundary unless the model and boundary conditions rule that out. The theorem above is stated on \(\mathbb R^d\) with decay precisely to avoid silently imposing a boundary condition on a pharmacokinetic state.

Common warning signs include a density whose integral drifts away from one, a negative diffusion coefficient, or a second-derivative term that uses \(\sigma\) instead of \(\sigma^2\) in one dimension. Also check derivative placement: when \(a\) depends on state, the equation contains \(\partial_{ij}(a_{ij}p)\), not \(a_{ij}\partial_{ij}p\). These errors can be found by tracing the coefficients back to the SDE and checking the dimensions of drift (state per time) and variance rate (state squared per time).

The forward equation is not the only way to describe the process. Moment equations, such as equations for \(E[X_t]\), summarize selected averages and may close for special models, but do not in general determine the full density. The backward Kolmogorov equation instead acts on a function of the starting state to describe conditional future expectations or probabilities. These approaches answer different questions: use the forward equation for evolution of the state distribution, a moment equation for a specified summary, and the backward equation for a future quantity conditional on the current state.

The forward equation alone does not guarantee a smooth density, specify a valid boundary condition, or provide a closed-form solution. Those depend on coefficient regularity, state space, initial law, and the particular model. The next steps for specialized models can build on this coefficient-to-equation translation without changing its basic rule.

Check Your Understanding

For each item, identify the requested coefficient or equation before reading the solution.

For \(dX_t=(3-2X_t)dt+0.5dW_t\), what is the diffusion matrix in one dimension, and what is the coefficient multiplying \(\partial_{xx}p\)?
Two states have independent diffusion amplitudes \(0.2\) and \(0.4\). What are the diagonal entries of \(a\), and is there a mixed second-derivative term?
For \(dC_t=-kC_tdt+\gamma C_tdW_t\), write the diffusion contribution to the forward equation.
Why does integrating a classical forward equation over all states preserve total probability under suitable decay conditions?

Solutions. (1) The diffusion matrix is \(a=(0.5)^2=0.25\). The equation’s second-derivative coefficient is \(\tfrac12(0.25)=0.125\), so the term is \(0.125\partial_{xx}p\). (2) The entries are \(0.2^2=0.04\) and \(0.4^2=0.16\). Independence gives zero off-diagonal covariance, so there is no mixed second-derivative term. (3) Here \(a(c)=\gamma^2c^2\), so the contribution is \(\tfrac12\partial_{cc}(\gamma^2c^2p)\). (4) Each term in the forward equation is a spatial derivative. Integration over the whole state space reduces these terms to boundary contributions, which vanish under the stated decay condition; hence the derivative of total mass is zero.

Forward Kolmogorov Equation for Pharmacokinetic Diffusions

What You'll Learn