A Worked Convergence Study for PK SDE Schemes

Why Refine the Simulation Step?

A numerical path can look plausible without being accurate. To assess discretization, simulate the same stochastic inputs on progressively finer grids and compare the resulting quantities. The comparison must be designed carefully: for a strong error, coarse and fine paths need to represent the same Brownian trajectory. A weak error, by contrast, compares expectations of a quantity of interest and does not require path-by-path closeness.

This study uses the one-compartment concentration model \(dC_t=k(\theta-C_t)\,dt+\sigma\,dW_t\), with \(C_0\) fixed, \(k>0\), and constant \(\theta\) and \(\sigma\). It is a convenient test case because its exact endpoint distribution and mean are known. It is not a claim that every PK SDE has additive noise or that every concentration simulation should use this model.

In applications the exact state is often unavailable. A fine-grid approximation \(C_T^{(\delta)}\) then acts as a reference, and comparisons estimate errors relative to that reference—not automatically errors relative to the exact SDE. Assume a fixed terminal time \(T\), grids that divide \(T\), and a common, correctly coupled Brownian path for every grid being compared.

Coupling the Grids and Estimating Error

Choose a finest increment size \(\delta\). For each replicate, generate independent Brownian increments of variance \(\delta\) on that grid. A coarser increment is the sum of all fine increments spanning its interval. This preserves the Brownian distribution and makes the coarse and fine paths use the same driving noise.

For \(M\) independent coupled replicates, write \(C_{T,i}^{(h)}\) for the coarse endpoint and \(C_{T,i}^{(\delta)}\) for the reference endpoint on replicate \(i\). The natural estimates are \(\widehat E_{\mathrm{strong}}(h)=\{M^{-1}\sum_{i=1}^M(C_{T,i}^{(h)}-C_{T,i}^{(\delta)})^2\}^{1/2}\) and \(\widehat E_{\mathrm{weak}}(h)=|M^{-1}\sum_{i=1}^M[g(C_{T,i}^{(h)})-g(C_{T,i}^{(\delta)})]|\). The second estimates the magnitude of a paired difference in expectations; it is not the average absolute path difference.

A Model with an Exact Weak-Error Benchmark

Set \(k=0.5\), \(\theta=0.6\), \(C_0=2\), \(\sigma=0.2\), and \(T=1\). The exact mean follows by solving the mean equation \(m'(t)=k(\theta-m(t))\), \(m(0)=C_0\): \(E[C_1]=0.6+1.4e^{-0.5}\approx1.449143\). Euler–Maruyama gives \(C_{j+1}=C_j+kh(\theta-C_j)+\sigma\Delta W_j\). Taking expectations, and using \(E[\Delta W_j]=0\), yields \(E[C_T^{(h)}]=0.6+1.4(1-0.5h)^{1/h}\) when \(1/h\) is an integer.

What Rate Should This Model Show?

For the smooth additive-noise model above, Euler–Maruyama has strong order one at a fixed terminal time under the stated assumptions \(k>0\), constant \(\theta,\sigma\), and finite \(T\). The usual generic strong-order-\(1/2\) guarantee for Euler–Maruyama is valid in broader settings, but it is not the expected asymptotic rate for this particular model. The reason is that its diffusion coefficient is constant: Euler–Maruyama captures the Brownian increment exactly, and the remaining within-step noise effect enters through the smooth drift.

The weak mean calculation is a useful independent check. Expanding \(\log(1-kh)=-kh+O(h^2)\) gives \((1-kh)^{T/h}=e^{-kT}(1+O(h))\), so the mean error is \(O(h)\). When simulations yield a strong rate near \(1/2\) in this model, do not immediately label that its asymptotic order: the generic guarantee is only a lower-order result here, and finite samples, reference-grid error, or an insufficiently fine range can obscure the model's order.

Diagnostics, Alternatives, and Limitations

A convergence table should include step size, replicate count, error estimate, and the rate computed from adjacent step sizes. A convincing rate is reasonably stable over several refinements, after transient coarse-grid behavior, and larger than the variation attributable to Monte Carlo error. For weak estimates, report the signed paired mean difference as well as its absolute value: cancellation can make the signed bias small even when replicate-level differences are variable. A standard error for the paired mean is the sample standard deviation of paired differences divided by \(\sqrt M\). The RMS strong estimate also has sampling variability; a bootstrap over independent paired replicates is one practical way to assess it.

The main alternative to a fine-grid reference is an exact solution or exact transition simulation, when available. That removes reference discretization error but may not be available for a nonlinear PK SDE. A fine-grid reference is broadly useful, but may be expensive and remains approximate. Independently simulated grids are adequate for comparing weak expectations in principle, but usually less efficient for paired comparisons because they discard the variance reduction from shared noise. In the convergence study, keep the Brownian coupling, terminal time, model parameters, and estimator fixed as \(h\) changes.

A rate can mislead if the finest reference is too coarse: both the reference and coarse approximations may share bias, or the reference error may dominate the difference. Check this by refining the reference itself and repeating the comparison. Also inspect path distributions, endpoint summaries, and the frequency and magnitude of negative concentrations. A positive-looking average does not rule out implausible individual paths; conversely, negative Euler–Maruyama values may be a model or scheme limitation rather than a convergence-rate calculation error. Do not silently truncate negative values, since truncation changes the numerical method and can alter both strong and weak errors.

A numerical convergence study measures the behavior of a scheme for a specified SDE; it does not establish that the SDE is an adequate pharmacokinetic model, that the simulated concentrations are clinically plausible, or that estimated parameters are reliable. It should be reported as a numerical diagnostic with its coupling, reference resolution, replicate count, and uncertainty made explicit.