Negative Concentrations from Euler-Maruyama

Why a Positive Model Can Produce a Negative Grid Value

Euler–Maruyama is a direct way to simulate an SDE on a time grid: over each interval, it evaluates the drift and diffusion at the start of the interval and uses one Brownian increment. That simplicity has a cost. The random increment has a normal distribution, which has support over the entire real line. Consequently, a step can cross below zero even when the state at the start of the step is positive—and even when the continuous-time model is positive.

For a PK amount or concentration, a negative simulated value is usually not a plausible state. It is important, however, to diagnose why it occurred rather than immediately changing it. The probability of a negative step depends on the current state, the local drift, the local diffusion, and the time-step length. The calculation below isolates these effects. It concerns one step between dose events; a dose jump, if present, is applied separately according to the model.

Definition: Consider the scalar SDE $dX_t=a(X_t,t)\,dt+b(X_t,t)\,dW_t$, where $W_t$ is Brownian motion. On a grid with spacing $h>0$, the Euler–Maruyama update from $t_n$ to $t_{n+1}=t_n+h$ is $$ X_{n+1}=X_n+a(X_n,t_n)h+b(X_n,t_n)\sqrt{h}\,Z_n, $$ where $Z_n=(W_{t_{n+1}}-W_{t_n})/\sqrt{h}\sim N(0,1)$. The grid value $X_n$ is assumed known at the start of the step, and the Brownian increment is independent of the information available then.

The One-Step Probability of Negativity

Condition on a starting value $X_n=x>0$, and abbreviate $a=a(x,t_n)$, $b=b(x,t_n)$. The update is normal with mean $x+ah$ and variance $b^2h$. If $b\ne 0$, standardizing this normal variable gives the exact one-step probability under the Euler–Maruyama update. If $b=0$, the update is deterministic instead.

Result: For $b\ne 0$, $$ P(X_{n+1}<0\mid X_n=x) =\Phi\!\left(\frac{-x-ah}{|b|\sqrt{h}}\right), $$ where $\Phi$ is the standard normal distribution function. For $b=0$, this probability is $1$ if $x+ah<0$, and $0$ if $x+ah\ge 0$.

To verify the formula, write the update as $X_{n+1}=x+ah+b\sqrt{h}Z_n$. When $b\ne0$, its conditional standard deviation is $|b|\sqrt{h}$, so $$ P(X_{n+1}<0\mid X_n=x) =P\!\left(\frac{X_{n+1}-(x+ah)}{|b|\sqrt{h}}< \frac{-x-ah}{|b|\sqrt{h}}\,\middle|\,X_n=x\right) =\Phi\!\left(\frac{-x-ah}{|b|\sqrt{h}}\right). $$ The absolute value is needed because the standard deviation is nonnegative, regardless of the sign of $b$. This proves the result. $\square$

The formula separates two effects. The drift changes the conditional mean to $x+ah$; the diffusion spreads the update with standard deviation $|b|\sqrt h$. A negative mean makes a negative update more likely, while greater diffusion makes the normal distribution wider. The effects of changing $h$ can interact: increasing $h$ changes both the mean and the standard deviation.

Worked Examples

Worked Example: A Positive Geometric Brownian Model

Suppose concentration follows geometric Brownian motion, $dX_t=\mu X_t\,dt+\sigma X_t\,dW_t$, with $X_n=x>0$, $\mu=0.2$ per hour, and $\sigma=0.8$ per square root hour. The exact continuous-time process remains positive when it starts positive. Find the probability that an Euler–Maruyama step is negative for $h=1$ hour and for $h=0.25$ hour.

Here the update is $X_{n+1}=x(1+\mu h+\sigma\sqrt h Z_n)$. Since $x>0$, the update is negative exactly when the parenthetical factor is negative.

For $h=1$, the factor has mean $1+0.2=1.2$ and standard deviation $0.8$. Thus $$ P(X_{n+1}<0\mid X_n=x) =P\!\left(Z_n<\frac{-1.2}{0.8}\right) =\Phi(-1.5)\approx 0.0668. $$ So about $6.68\%$ of such steps are negative under this update.

For $h=0.25$, the factor has mean $1+0.2(0.25)=1.05$ and standard deviation $0.8\sqrt{0.25}=0.4$. Therefore $$ P(X_{n+1}<0\mid X_n=x) =\Phi\!\left(\frac{-1.05}{0.4}\right) =\Phi(-2.625)\approx 0.00433. $$ The shorter step lowers the one-step probability, but does not make it zero. The negative values are a property of the numerical update, not of the exact geometric Brownian process.

Worked Example: A Square-Root PK State

Consider the square-root diffusion $dX_t=\kappa(\theta-X_t)\,dt+\sigma\sqrt{X_t}\,dW_t$, with $X_n=0.5$, $\kappa=1$ per hour, $\theta=0.5$, and $\sigma=0.8$ units of square root concentration per square root hour. The Feller quantity is $2\kappa\theta=1$, which exceeds $\sigma^2=0.64$; for this CIR model, a process started positive does not reach zero. Calculate the Euler–Maruyama negative-step probabilities for $h=1$ and $h=0.25$ hours.

At $x=0.5$, the drift is $\kappa(\theta-x)=0$, and the diffusion coefficient is $0.8\sqrt{0.5}$. For $h=1$, the update has mean $0.5$ and standard deviation $0.8\sqrt{0.5}\approx0.5657$. Hence $$ P(X_{n+1}<0\mid X_n=0.5) =\Phi\!\left(\frac{-0.5}{0.5657}\right) =\Phi(-0.8839)\approx0.1884. $$

For $h=0.25$, the drift at the starting state is still zero, so the mean remains $0.5$. The standard deviation is $0.8\sqrt{0.5}\sqrt{0.25}\approx0.2828$. Thus $$ P(X_{n+1}<0\mid X_n=0.5) =\Phi\!\left(\frac{-0.5}{0.2828}\right) =\Phi(-1.7678)\approx0.0386. $$ The process’s positivity does not imply positivity of this discretization. The step-size change reduces the one-step risk substantially, but a nonzero Gaussian tail remains.

Worked Example: Checking a Very Small Tail Probability

At a particular positive state $x=1.01$, suppose the locally evaluated drift is $a=0$, the diffusion coefficient is $b=1$, and the step is $h=0.01$. Calculate the one-step probability of a negative update.

The conditional mean is $x+ah=1.01$, and the standard deviation is $|b|\sqrt h=0.1$. Therefore the standardized cutoff for zero is $$ \frac{-x-ah}{|b|\sqrt h}=\frac{-1.01}{0.1}=-10.1. $$ It follows that $$ P(X_{n+1}<0\mid X_n=1.01)=\Phi(-10.1)\approx 2.76\times10^{-24}. $$ This is an extremely small one-step probability. It is not exactly zero, and it should not be confused with a general guarantee of positivity. This calculation describes the stated local update; whether the underlying continuous model is appropriate on the positive state space is a separate modeling question.

How Drift and Step Size Change the Risk

For fixed $x>0$ and constant local coefficients $a<0$, $b\ne0$, consider $z(h)=(x+ah)/(|b|\sqrt h)$, so that the negative-step probability is $\Phi(-z(h))$. Differentiating gives $$ z'(h)=\frac{ah-x}{2|b|h^{3/2}}. $$ Because $a<0$, $ah-x<0$, and therefore $z'(h)<0$ for every $h>0$. Since $\Phi$ is increasing, $\Phi(-z(h))$ increases as $h$ increases. This proves that, with fixed negative drift and diffusion at the starting state, increasing the step length increases the one-step negative probability. $\square$

This monotonicity conclusion is specific to negative drift. With positive drift, increasing $h$ can initially move the mean farther above zero; it is not valid to claim that a larger step always increases the negative probability. For state-dependent coefficients, the calculation is conditional on the current state: different values of $x$ can have different drift, diffusion, and risk.

A useful way to choose or assess a step size is to set a tolerated one-step risk $\varepsilon$, with $0<\varepsilon<0.5$. Let $z_{1-\varepsilon}$ be the $(1-\varepsilon)$-quantile of the standard normal distribution, defined by $\Phi(z_{1-\varepsilon})=1-\varepsilon$. When $b\ne0$, the formula above shows that the conditional negative probability is at most $\varepsilon$ exactly when $$ x+ah\ge z_{1-\varepsilon}|b|\sqrt h. $$ This condition is a local diagnostic, not a universal step-size guarantee for an entire path: later states, and therefore later coefficients, are random.

Diagnostics and Alternatives

A negative simulated value should prompt checks of both the numerical method and the model. First, compare the observed fraction of negative updates with the conditional probability calculated from the states at the start of those steps. A large discrepancy can indicate an incorrect implementation, such as using a variance $h$ rather than a standard deviation $\sqrt h$ for the Brownian increment. Second, inspect the risk by state and by step size. Negativity concentrated at small states or at large steps is consistent with the one-step formula; frequent negatives despite a very small predicted probability warrant further investigation.

For a sequence of steps, the one-step probability is not the probability that a whole simulated path ever becomes negative. The states on later steps depend on earlier increments, so the step events are generally dependent. If the conditional probability of a negative step is at most $q$ at every step before a negative value occurs, then over $N$ steps the probability of at least one negative value is at most $Nq$, by the union bound. This bound may be loose, but it warns that a tiny per-step risk is not automatically negligible over a very long grid.

Approach	What it changes	Main diagnostic or limitation
Reduce $h$	Uses more, shorter Euler–Maruyama steps.	Often reduces one-step risk, but does not generally eliminate it.
Transform the state	Simulates a coordinate such as log concentration and maps back to a positive scale.	Requires a valid transformed SDE and care near boundaries or at zero.
Use a positivity-preserving or boundary-aware scheme	Changes the update to respect the state space.	Check its assumptions and whether it targets the intended model.
Clip negative values to zero	Replaces a negative update after it has been generated.	Changes the simulated distribution; it is not the original Euler–Maruyama update.

Clipping is particularly easy to misinterpret. Replacing $X_{n+1}$ by $\max(0,X_{n+1})$ creates a point mass at zero whenever the original update is negative. It therefore changes transition probabilities and can affect quantities such as average exposure. A better response is to quantify the problem, check that the continuous model is suitable for a nonnegative state, and select a numerical method consistent with that model. The earlier material on positivity-preserving models and simulation develops boundary-aware alternatives; the focus here is diagnosing the Euler–Maruyama update itself.

Key takeaway: For a positive starting value, the one-step Euler–Maruyama risk is determined by the normal tail below zero with mean $x+ah$ and standard deviation $|b|\sqrt h$. Negative grid values can arise even for a positive continuous-time model, so assess them using the conditional formula rather than assuming that positivity of the model guarantees positivity of the numerical path.

Check Your Understanding

Answer each question, then compare with the solutions.

For $x=2$, $a=-0.5$, $b=0.4$, and $h=1$, calculate the one-step probability of a negative update.
For fixed $x>0$, $a<0$, and $b\ne0$, does increasing $h$ increase or decrease the one-step negative probability? State why.
Why can Euler–Maruyama produce a negative value for geometric Brownian motion even though its exact process stays positive?
What is the main statistical consequence of clipping every negative update to zero?

Solutions. (1) The mean is $2-0.5=1.5$, and the standard deviation is $0.4$. Thus the probability is $\Phi(-1.5/0.4)=\Phi(-3.75)\approx0.0000884$. (2) It increases: $z(h)=(x+ah)/(|b|\sqrt h)$ has derivative $(ah-x)/(2|b|h^{3/2})<0$, so $\Phi(-z(h))$ increases. (3) Its Euler–Maruyama update is $x(1+\mu h+\sigma\sqrt h Z_n)$; the normal random factor can be negative, even though the exact geometric Brownian solution is positive. (4) Clipping creates probability mass at zero and changes the simulated distribution, so it no longer represents the original update.

Approach	What it changes	Main diagnostic or limitation
Reduce \(h\)	Uses more, shorter Euler–Maruyama steps.	Often reduces one-step risk, but does not generally eliminate it.
Transform the state	Simulates a coordinate such as log concentration and maps back to a positive scale.	Requires a valid transformed SDE and care near boundaries or at zero.
Use a positivity-preserving or boundary-aware scheme	Changes the update to respect the state space.	Check its assumptions and whether it targets the intended model.
Clip negative values to zero	Replaces a negative update after it has been generated.	Changes the simulated distribution; it is not the original Euler–Maruyama update.