Positivity-Preserving Simulation of Stochastic Exchange

Why Simulate Exchange on the Amount Scale?

In a closed two-compartment exchange model, material can move from either compartment to the other, but it is neither created nor eliminated. A simulation should therefore respect two physical features: each compartment amount is nonnegative, and the total amount remains constant. A discretization that violates either feature can produce paths that are mathematically computable but pharmacologically uninterpretable.

We build the simulation around the fraction of total amount in compartment 1. This makes the state space explicit: the fraction must lie between zero and one. The resulting bounded transition is useful when an ordinary unconstrained Gaussian increment would risk leaving the allowable range.

The Exchange Diffusion and Its State Space

Let $X_1(t)$ and $X_2(t)$ be the amounts in two compartments, and let their total be $M>0$. Assume there is no input or elimination during the interval being modeled. Let $k_{12}>0$ and $k_{21}>0$ be the first-order transfer rates from compartment 1 to 2 and from 2 to 1, respectively. Define the compartment-1 fraction $P(t)=X_1(t)/M$, so $X_2(t)=M(1-P(t))$.

Definition: The two-compartment stochastic exchange diffusion is $$ dP(t)=\bigl[k_{21}(1-P(t))-k_{12}P(t)\bigr]\,dt +\eta\sqrt{P(t)(1-P(t))}\,dW(t), $$ where $W(t)$ is standard Brownian motion, $\eta\geq 0$, and $P(0)\in[0,1]$. The amount processes are $X_1(t)=MP(t)$ and $X_2(t)=M(1-P(t))$. Equivalently, $$ dX_1=(k_{21}X_2-k_{12}X_1)\,dt+\eta\sqrt{X_1X_2}\,dW, \qquad dX_2=-dX_1. $$

The drift moves the fraction toward the equilibrium share $k_{21}/(k_{12}+k_{21})$. The diffusion coefficient vanishes at both endpoints, where one compartment has zero amount. The standard Jacobi-diffusion boundary result applies: with positive $k_{12},k_{21}$, finite $\eta$, and an initial fraction in $[0,1]$, the solution exists uniquely and remains in $[0,1]$. The endpoints may be attainable for some parameter values, but the process does not cross them. Thus both amount processes remain nonnegative.

This continuous-time property does not automatically guarantee that an arbitrary numerical update will remain in the state space. The update itself must be designed to do so.

Two Properties to Preserve and Check

First, the total is conserved pathwise, not just on average. Adding the two amount equations cancels both the transfer drift and the shared stochastic increment:

d(X_1+X_2)=0, \qquad X_1(t)+X_2(t)=M.

Second, the mean follows a closed equation. Write $K=k_{12}+k_{21}$. Since $P(t)$ stays bounded, the stochastic integral has mean zero on finite time intervals. Taking expectations and writing $m(t)=E[P(t)]$ gives

\frac{dm(t)}{dt}=k_{21}-Km(t), \qquad m(t)=\frac{k_{21}}{K} +\left(P(0)-\frac{k_{21}}{K}\right)e^{-Kt}.

The mean amount in compartment 1 is $M m(t)$; the mean in compartment 2 is $M[1-m(t)]$. These expressions provide reference values for checking simulations. They do not imply that individual paths follow the mean trajectory.

Worked Example: Mean Amounts After One Time Unit

Take $M=1$ mg, $P(0)=0.7$, $k_{12}=0.4$ per time unit, and $k_{21}=0.2$ per time unit. The rates give $K=0.6$ and equilibrium fraction $k_{21}/K=1/3$. At $t=1$,

$$ E[P(1)]=\frac{1}{3}+\left(0.7-\frac{1}{3}\right)e^{-0.6} =\frac{1}{3}+\frac{11}{30}(0.5488116361) \approx 0.534564. $$

Therefore $E[X_1(1)]=1(0.534564)\approx 0.5346$ mg and $E[X_2(1)]=1-0.534564\approx 0.4654$ mg. Their sum is $1.0000$ mg, as required. The diffusion intensity $\eta$ affects path variability but not these means.

A Bounded Beta-Transition Simulation

A convenient one-step simulation uses a beta distribution, whose support is $(0,1)$. For a time step $h>0$, first propagate the drift exactly. With $p$ denoting the current fraction, set

K=k_{12}+k_{21},\qquad \pi=\frac{k_{21}}{K},\qquad m=\pi+(p-\pi)e^{-Kh}.

This is the solution at time $h$ of the drift-only equation starting from $p$, and it lies in $[0,1]$ whenever $p$ does. For $\eta>0$, require $\eta^2h<1$ and define

\kappa=\frac{1}{\eta^2h}-1,\qquad P_{\mathrm{new}}\mid p\sim \operatorname{Beta}(m\kappa,(1-m)\kappa).

If $\eta=0$, set $P_{\mathrm{new}}=m$ deterministically. Given the new fraction, set $X_{1,\mathrm{new}}=MP_{\mathrm{new}}$ and $X_{2,\mathrm{new}}=M(1-P_{\mathrm{new}})$.

For a beta random variable with positive parameters $\alpha,\beta$, its mean is $\alpha/(\alpha+\beta)$ and its variance is $\alpha\beta/[(\alpha+\beta)^2(\alpha+\beta+1)]$. Substitution gives the update’s conditional moments:

E[P_{\mathrm{new}}\mid p]=m,\qquad \operatorname{Var}(P_{\mathrm{new}}\mid p)=\eta^2h\,m(1-m).

As $h$ decreases, $m=p+[k_{21}(1-p)-k_{12}p]h+O(h^2)$, and $m(1-m)=p(1-p)+O(h)$. Consequently, the conditional mean increment has the SDE drift to first order, and the conditional variance is $\eta^2p(1-p)h+O(h^2)$. This verifies local moment agreement; it is not, by itself, a proof of a global convergence rate.

Worked Example: One Beta Transition

Use the rates and initial fraction from the preceding example, with $h=0.1$ and $\eta^2=0.5$ per time unit. The drift-propagated mean is

$$ m=\frac{1}{3}+\left(0.7-\frac{1}{3}\right)e^{-0.06} \approx 0.678647. $$

Here $\eta^2h=0.05<1$, so $\kappa=1/0.05-1=19$. The beta parameters are $\alpha=19m\approx12.8943$ and $\beta=19(1-m)\approx6.1057$. Thus the update draws $P_{\mathrm{new}}$ from $\operatorname{Beta}(12.8943,6.1057)$, with conditional variance

$$ 0.05m(1-m)\approx 0.05(0.218085)=0.010904. $$

Every draw lies between zero and one, so both updated amounts are nonnegative. Their sum is exactly $M$, regardless of the particular draw. The expected compartment-1 amount after this step is $Mm$, because the transition’s conditional mean is $m$.

Worked Example: A Deterministic Exchange Step

Set $\eta=0$, $M=4$ mg, $P(0)=0.25$, $k_{12}=0.5$, $k_{21}=0.3$, and $h=2$. Then $K=0.8$, $\pi=0.3/0.8=0.375$, and

$$ P_{\mathrm{new}}=0.375+(0.25-0.375)e^{-1.6} =0.375-0.125(0.2018965) \approx0.349763. $$

The updated amounts are $X_{1,\mathrm{new}}=4(0.349763)\approx1.3991$ mg and $X_{2,\mathrm{new}}=4-1.3991\approx2.6009$ mg. The deterministic drift step preserves both nonnegativity and the total, even though the step is relatively long.

Why the Update Preserves Positivity

Assume the current fraction is in $[0,1]$, both transfer rates are positive, and $h>0$. The exact drift flow is a weighted average of $p$ and $\pi$, since $m=e^{-Kh}p+(1-e^{-Kh})\pi$. Both weights are nonnegative and sum to one, so $m\in[0,1]$. If $m$ is strictly between zero and one and $\eta^2h<1$, both beta parameters are positive and the draw is in $(0,1)$ with probability one. At an endpoint, the update can be interpreted by its limiting degenerate distribution; for positive rates and $h>0$, the drift flow moves an endpoint into the interior. Multiplication by $M>0$ then gives nonnegative amounts.

The total is preserved by construction at every step: $MP_{\mathrm{new}}+M(1-P_{\mathrm{new}})=M$. Repeating the update therefore preserves the total over the entire grid, up to floating-point rounding in an implementation. It also preserves the exact continuous-time mean at grid points. Indeed, the conditional mean is the affine drift flow $m$, so taking expectations at each step gives the same linear mean evolution over that step.

Theorem: Under $M>0$, $k_{12},k_{21}>0$, initial $P_0\in[0,1]$, and step size satisfying $\eta^2h<1$ when $\eta>0$, the beta-transition scheme has nonnegative compartment amounts and total $M$ at every grid point. Its one-step conditional mean equals the exact drift flow.

Proof. The drift-flow formula is a convex combination of $p$ and $\pi$, hence belongs to $[0,1]$. The beta draw is supported on $(0,1)$ when its parameters are positive; the zero-noise update is $m$ itself. Therefore $X_{1,\mathrm{new}}$ and $X_{2,\mathrm{new}}$ are nonnegative. Their sum is $M$ by their definitions. Finally, the beta mean formula gives $E[P_{\mathrm{new}}\mid p]=m$; in the zero-noise case the same identity holds deterministically. Iterating these facts proves the grid-point claims. $\square$

Diagnostics, Alternatives, and Limitations

For a simulated path, check each update rather than relying only on a final summary. The fractions should remain in $[0,1]$, the amounts should remain nonnegative, and $X_1+X_2$ should equal its initial value within numerical precision. Across many paths, compare the sample mean at selected times with $M\{\pi+(P_0-\pi)e^{-Kt}\}$. A systematic mean discrepancy suggests an implementation error, an incorrect rate direction, or an unintended change to the model.

The condition $\eta^2h<1$ is essential for this beta parameterization: if it fails, $\kappa$ is not positive and the specified beta distribution is invalid. Reduce the step size or use a different transition method; do not silently substitute negative beta parameters. Also check units: $k_{12}$ and $k_{21}$ have units of inverse time, while $\eta^2h$ must be dimensionless.

The beta update is not an exact sample from the SDE transition distribution. It matches the drift flow in its conditional mean and has the correct leading-order conditional variance, but finite-step distributions and higher moments can differ. Smaller steps reduce the local discrepancy, at additional computational cost. A Gaussian increment can be convenient but is unbounded and does not inherently respect the fraction’s interval. A transformed-coordinate approach, such as one based on Lamperti coordinates, handles boundaries differently and requires its own boundary-aware simulation choices. A discrete jump-exchange model is another alternative when individual transfer events, rather than a continuous diffusion approximation, are the intended mechanism. None of these methods removes the need to verify that the chosen model represents the pharmacological process of interest.

Common Pitfall

Clipping an invalid fraction after a Gaussian update. Replacing a value below zero by zero, or above one by one, may restore bounds, but it changes the transition distribution and generally changes its conditional mean and variance. The resulting path is not automatically a faithful simulation of the stated diffusion.
Checking only the mean total. A simulation can have the right average total while individual paths create or lose amount. Check $X_1+X_2=M$ path by path.
Confusing mean behavior with individual behavior. The mean follows a deterministic exchange equation, but stochastic paths fluctuate around it. Matching the mean alone does not verify the variance or the transition implementation.

Key takeaway: Simulate the exchange fraction with a transition supported on $[0,1]$, then define the second amount as the remainder of the conserved total. This construction makes nonnegativity and total conservation directly checkable at every grid point.

Check Your Understanding

Answer each question, then compare with the brief solutions.

For $k_{12}=0.3$ and $k_{21}=0.2$, what is the equilibrium fraction in compartment 1?
If $M=5$ mg and a beta-transition draw is $P_{\mathrm{new}}=0.6$, what are the two amounts and their total?
For $\eta^2=0.4$ and $h=0.5$, calculate $\kappa$ and check whether the beta step is admissible.
Why does the exchange diffusion have the same mean equation for every value of $\eta$?

Solutions. 1. $K=0.5$, so the equilibrium fraction is $0.2/0.5=0.4$. 2. The amounts are $5(0.6)=3$ mg and $5(0.4)=2$ mg; their total is $5$ mg. 3. $\eta^2h=0.2<1$, so $\kappa=1/0.2-1=4$, and the beta parameterization is admissible. 4. The diffusion term is a stochastic integral with mean zero because the fraction is bounded; taking expectations leaves only the linear drift equation.