Stochastic Reaction Jump Simulation and Coupling

Fix positive integers $d$ and $R$. A reaction propensity system consists of a state set $E\subseteq {\mathbb{Z}}^d$ that is finite or countable, increments ${\nu }_j\in {\mathbb{Z}}^d$ for $1\le j\le R$, and propensities $a_j:E\to [0,\mathrm{\infty })$. Here “countable” means finite or countably infinite, as in Finite, countably infinite, countable, uncountable. The total propensity $\lambda (x):={\sum }_{j=1}^R{a}_j(x)$ is finite for every $x\in E$, and $a_j(x)>0$ implies $x+{\nu }_j\in E$. Here $\mathbb{N}$ and $\mathbb{Z}$ have the meanings in The natural numbers $\mathbb{N}$ (von Neumann) and The integers as equivalence classes of pairs of naturals, and addition in ${\mathbb{Z}}^d$ is coordinatewise.

A reaction event of channel $j$ at state $x$ records the marked transition $x\mapsto x+{\nu }_j$. The counting process $N_j(t)$ records the number of channel-$j$ events by time $t$. A zero increment ${\nu }_j=0$ is allowed: it records a reaction event without changing the state.

For a time-homogeneous system, the reaction generator operator on bounded functions $f:E\to \mathbb{R}$ is $(Qf)(x):={\sum }_{j:a_j(x)>0}{a}_j(x)(f(x+{\nu }_j)-f(x))\mathrm{.}$ When $E$ is finite, set $\mathrm{\Lambda }:=\max (\{0\}\cup \{\lambda (x):x\in E\})<\mathrm{\infty }$; thus $\mathrm{\Lambda }=0$ for $E=\mathrm{\varnothing }$ and $\mathrm{\Lambda }={\max }_{x\in E}\lambda (x)$ otherwise. If no active reaction has zero increment, the corresponding matrix has off-diagonal entries $q_{xy}={\sum }_{j:x+{\nu }_j=y}{a}_j(x)$ and diagonal $q_{xx}=-\lambda (x)$. With zero-increment events, the operator formula remains valid and the matrix diagonal is instead minus the sum of state-changing rates; the zero-increment channels remain part of the marked event process.

For a finite horizon $T$, a time-dependent propensity system replaces $a_j(x)$ by functions $a_j:[0,T]\times E\to [0,\mathrm{\infty })$ whose sections $t\mapsto a_j(t,x)$ are measurable for the Borel sigma-algebra $\mathcal{B}([0,T])$ of the usual topology on $[0,T]$ (The Borel sigma-algebra of a topological space). Require finite total rate $\lambda (t,x):={\sum }_j{a}_j(t,x)$. Its time-dependent generator operator is $(Q_{t}f)(x):={\sum }_{j:a_j(t,x)>0}{a}_j(t,x)(f(x+{\nu }_j)-f(x))\mathrm{.}$ This definition specifies a stochastic model. It does not assert that a physical chemical system follows that model.

Assume Countable Choice. Let $R\ge 1$, let $r_1,\dots ,r_R\ge 0$, and set $r_0={\sum }_{j=1}^R{r}_j$. On a common probability space let $E_1,\dots ,E_R$ be independent random variables, each with the unit exponential law $\mathrm{Exp}(1)$, whose survival at $t\ge 0$ is $e^{-t}$. If $r_0>0$, define $T_{\ast }:={\min }_{j:r_j>0}\frac{E_j}{r_j},J_{\ast }:={\mathrm{argmin}}_{j:r_j>0}\frac{E_j}{r_j}\mathrm{.}$ Then ties have probability zero, $T_{\ast }\sim \mathrm{Exp}(r_0)$, $J_{\ast }$ is independent of $T_{\ast }$, and $\mathbb{P}(J_{\ast }=j)=\frac{r_j}{r_0}\mathrm{.}$ If $r_0=0$, every $r_j=0$ and no finite reaction clock is active.

Inputs. A finite state set $E\subseteq {\mathbb{Z}}^d$, a finite-channel reaction propensity system from Finite-channel reaction propensity systems and their generator, an initial state $x_0\in E$, a horizon $T>0$, and an optional event budget $B\in {\mathbb{N}}_0\cup \{\mathrm{\infty }\}$, defaulting to $\mathrm{\infty }$, where ${\mathbb{N}}_0:=\{0,1,2,\dots \}$. The random source is a sequence $(E_n,U_n{)}_{n\ge 1}$ with $E_n\sim \mathrm{Exp}(1)$ and $U_n$ uniform on $[0,1]$, all coordinates mutually independent. The countable product construction for these input laws uses Countable Choice and Dependent Choice by Assuming countable and dependent choice, countable products of arbitrary probability spaces; the exponential law uses Countable Choice by Exponential and Weibull survival specifications define probability laws. For the uniform law, Assuming countable choice, $\mathcal{L}({\mathbb{R}}^n)$ is a sigma-algebra containing every elementary set and ${\lambda }_n$ is a complete measure extending elementary volume makes the set function in Lebesgue measurable sets, the family $\mathcal{L}({\mathbb{R}}^n)$, and the restricted set function ${\lambda }_n$ a measure, and A box in ${\mathbb{R}}^n$ with parameters $a_i\le b_i$ is Lebesgue measurable of measure ${\prod }_{i<n}(b_i-a_i)$, whichever of its faces are included gives its restriction to $[0,1]$ mass one and interval probabilities equal to interval lengths.

Procedure. Set $t=0$, $x=x_0$, $n=0$, and record the initial state. While $t<T$:

1. Compute $\lambda (x)={\sum }_{j=1}^R{a}_j(x)$. If $\lambda (x)=0$, keep the state equal to $x$ through $T$ and return status absorbed.
2. If $n=B$, return the path prefix through $(t,x)$ with status budget-exhausted; this status does not assert that the horizon was reached.
3. Set $\tau =E_{n+1}\mathrm{/}\lambda (x)$. If $t+\tau \ge T$, keep the state equal to $x$ through $T$ and return status horizon.
4. Set $J$ to the least index with $U_{n+1}<\lambda (x{)}^{-1}{\sum }_{j=1}^J{a}_j(x)$. If the null event $U_{n+1}=1$ occurs, use the largest index with $a_j(x)>0$. This rule never selects a zero-rate channel.
5. Set $t\leftarrow t+\tau$, $x\leftarrow x+{\nu }_J$, append the marked event $(t,J,x)$, and set $n\leftarrow n+1$.

If the loop reaches $t=T$, return the path through $T$ with status horizon. Each event uses a linear scan of the $R$ propensities and a linear cumulative channel search; no update-tree acceleration is part of this procedure. If $B<\mathrm{\infty }$ and the budget expires before $T$, the returned prefix is incomplete.

The procedure assumes an ideal independent random source with the displayed laws. It makes no claim that a finite-state pseudorandom generator supplies independent exact exponential or uniform variables.

Inputs. A finite state set $E\subseteq {\mathbb{Z}}^d$, a finite-channel time-homogeneous system from Finite-channel reaction propensity systems and their generator, initial state $x_0\in E$, horizon $T>0$, and optional reaction-event budget $B\in {\mathbb{N}}_0\cup \{\mathrm{\infty }\}$, defaulting to $\mathrm{\infty }$, where ${\mathbb{N}}_0:=\{0,1,2,\dots \}$. Supply independent unit exponential increments $E_{j,k}\sim \mathrm{Exp}(1)$ for $1\le j\le R$ and $k\ge 1$. The countable product construction uses Countable Choice and Dependent Choice by Assuming countable and dependent choice, countable products of arbitrary probability spaces; the exponential law assumes Countable Choice by Exponential and Weibull survival specifications define probability laws.

Initialization. Set $t=0$, $x=x_0$, $n=0$, internal times $A_j=0$, threshold indices $m_j=1$, and next thresholds $S_j=E_{j,1}$ for each channel.

Procedure. At each state $x$:

1. Compute $a_j(x)$ for every channel. If $\lambda (x)=0$, return the constant path through $T$ with status absorbed.
2. If $n=B$, return the path prefix with status budget-exhausted. For each $j$, set ${\mathrm{\Delta }}_j=(S_j-A_j)\mathrm{/}{a}_j(x)$ when $a_j(x)>0$ and ${\mathrm{\Delta }}_j=+\mathrm{\infty }$ when $a_j(x)=0$.
3. Let $\mathrm{\Delta }={\min }_j{\mathrm{\Delta }}_j$. If $t+\mathrm{\Delta }>T$, return the path held at $x$ through $T$ with status horizon.
4. Let $J$ be the least minimizing channel (ties have probability zero). For every channel set $A_j\leftarrow A_j+a_j(x)\mathrm{\Delta }$. Set $t\leftarrow t+\mathrm{\Delta }$, $x\leftarrow x+{\nu }_J$, append the marked event, then set $m_J\leftarrow m_J+1$ and $S_J\leftarrow S_J+E_{J,m_J}$; finally set $n\leftarrow n+1$.

On the returned time interval (through $T$ unless the budget is exhausted), the internal time is $A_j(t)={\int }_0^t{a}_j(X_{s-})ds$, and almost surely the random-time-change representation is $X_t=X_0+{\sum }_{j=1}^R{\nu }_j{Y}_j(A_j(t)),$ where $S_{j,k}={\sum }_{q=1}^k{E}_{j,q}$ and $Y_j(u)=\mathrm{\#}\{k:S_{j,k}\le u\}$. An event at exactly $T$ is recorded before returning. If the loop reaches $t=T$, return the path through $T$ with status horizon. With unlimited budget, Finite-state Gillespie and next-reaction methods have the specified CTMC law ↗ proves the Poisson threshold-clock laws and equality with the direct method's reaction path law. The implementation scans all $R$ channels at each event; no logarithmic update-tree cost is claimed.

Inputs. A finite or countable state set $E\subseteq {\mathbb{Z}}^d$, finitely many time-dependent propensities $a_j(t,x)$ from Finite-channel reaction propensity systems and their generator, a finite horizon $T>0$, initial state $x_0\in E$, optional candidate budget $B\in {\mathbb{N}}_0\cup \{\mathrm{\infty }\}$ defaulting to $\mathrm{\infty }$ (where ${\mathbb{N}}_0:=\{0,1,2,\dots \}$), and a deterministic envelope $\mathrm{\Lambda }<\mathrm{\infty }$ with an admissibility certificate establishing $0\le \lambda (t,x):={\sum }_{j=1}^R{a}_j(t,x)\le \mathrm{\Lambda }(0\le t\le T,x\in E)\mathrm{.}$ The global inequality is a precondition; the procedure does not verify the certificate exhaustively. The random source supplies independent $E_n\sim \mathrm{Exp}(1)$ and $U_n\sim \mathrm{Unif}[0,1]$ sequences. Their countable product model uses Countable Choice and Dependent Choice; the unit exponential law is as in Exponential and Weibull survival specifications define probability laws. For the uniform law, Assuming countable choice, $\mathcal{L}({\mathbb{R}}^n)$ is a sigma-algebra containing every elementary set and ${\lambda }_n$ is a complete measure extending elementary volume makes the set function in Lebesgue measurable sets, the family $\mathcal{L}({\mathbb{R}}^n)$, and the restricted set function ${\lambda }_n$ a measure, and A box in ${\mathbb{R}}^n$ with parameters $a_i\le b_i$ is Lebesgue measurable of measure ${\prod }_{i<n}(b_i-a_i)$, whichever of its faces are included gives its restriction to $[0,1]$ mass one and interval probabilities equal to interval lengths.

Procedure. Set $t=0$, $x=x_0$, and candidate counter $n=0$. If $\mathrm{\Lambda }=0$, return the constant path through $T$ with status absorbed. Otherwise, while $t<T$:

1. If $n=B$, return the path prefix with status budget-exhausted.
2. Set $\tau =E_{n+1}\mathrm{/}\mathrm{\Lambda }$. If $t+\tau \ge T$, return the state held through $T$ with status horizon.
3. Set $t\leftarrow t+\tau$ and evaluate $\lambda (t,x)$. If $\lambda (t,x)>\mathrm{\Lambda }$ at this candidate, stop with status invalid-envelope; this pointwise runtime check detects only violations at candidate times and states that are actually evaluated. The algorithm does not continue by clipping the rate.
4. If $U_{n+1}\ge \lambda (t,x)\mathrm{/}\mathrm{\Lambda }$, reject this candidate and leave $x$ unchanged; set $n\leftarrow n+1$ and continue. Otherwise choose the first $J$ such that $U_{n+1}\mathrm{\Lambda }<{\sum }_{j=1}^J{a}_j(t,x),$ set $x\leftarrow x+{\nu }_J$, append the marked event $(t,J,x)$, set $n\leftarrow n+1$, and continue.

A time-dependent propensity without a certified finite envelope is outside this algorithm; integrated-hazard inversion is a separate method and is not used here. If the loop reaches $t=T$, return status horizon. If the global precondition is false, the procedure need not detect that failure: status invalid-envelope is returned only when the pointwise check in step 3 encounters a violation. No path law is claimed under a false global envelope.

Assume Countable Choice and Dependent Choice. Let $E\subseteq {\mathbb{Z}}^d$ be countable and let a finite-channel system be as in Finite-channel reaction propensity systems and their generator. For an initial state $x_0\in E$, construct its minimal reaction process recursively: at state $x$, wait an exponential time with parameter $\lambda (x)$, choose channel $j$ with probability $a_j(x)\mathrm{/}\lambda (x)$, and move to $x+{\nu }_j$; if $\lambda (x)=0$, remain at $x$ forever. Use an independent stream of these exponential and uniform random inputs, whose countable product law exists by Countable Choice and Dependent Choice.

For $V:E\to [1,\mathrm{\infty })$, call $V$ proper when $\{x:V(x)\le n\}$ is finite for every integer $n\ge 1$. For this possibly unbounded $V$, define the finite local expression $(QV)(x):={\sum }_{j:a_j(x)>0}{a}_j(x)(V(x+{\nu }_j)-V(x));$ this does not assert that $V$ belongs to the bounded-function domain of the generator in [F1]. Suppose $QV\le cV$ for a constant $c\ge 0$. Define $T_m$ as the time of the $m$th reaction event, with $T_m=\mathrm{\infty }$ if fewer than $m$ events occur. Then the explosion time ${\tau }_{\mathrm{\infty }}:={\lim }_{m\to \mathrm{\infty }}{T}_m$ satisfies ${\mathbb{P}}_{x_0}({\tau }_{\mathrm{\infty }}=\mathrm{\infty })=1$. No full Axiom of Choice is used.

Assume Countable Choice and Dependent Choice. Let $E\subseteq {\mathbb{Z}}^d$ be finite or countable, with finitely many time-dependent propensities $a_j(t,x)$ on $[0,T]$ as in Finite-channel reaction propensity systems and their generator. Suppose a certified constant envelope $0\le \lambda (t,x):={\sum }_j{a}_j(t,x)\le \mathrm{\Lambda }<\mathrm{\infty }$ holds for all $t\in [0,T]$ and $x\in E$. Then the dominating-thinning algorithm in Dominating thinning for bounded time-varying propensities, with an unlimited candidate budget, has the following exact conditional law from time $t$ and state $x$:

Write ${\mathcal{F}}_t$ for the sigma-algebra generated by the candidate times and marks in $(0,t]$. Then, for $0\le t\le T$, $x\in E$, and $0\le u\le T-t$,

$\mathbb{P}(\text{no accepted reaction in }(t,t+u]\mid X_t=x,{\mathcal{F}}_t)=\exp \text{\hspace{0.17em}⁣}(-{\int }_t^{t+u}\lambda (s,x)ds),0\le u\le T-t,$

and for every Borel $B\subseteq (t,T]$ and channel $j$,

$\mathbb{P}(\text{first accepted reaction time lies in }B,J=j\mid X_t=x,{\mathcal{F}}_t)={\int }_B{a}_j(s,x)\exp \text{\hspace{0.17em}⁣}(-{\int }_t^s\lambda (v,x)dv)ds\mathrm{.}$

After a reaction, the same rule restarts at its new time and state. Candidate events on $[0,T]$ have a Poisson count of mean $\mathrm{\Lambda }T$; accepted reactions are a subset, so there are finitely many accepted reactions almost surely and $\mathbb{E}N(T)\le \mathrm{\Lambda }T$. The conclusion is conditional on the stated envelope; no law is asserted when it fails.

When $\mathrm{\Lambda }>0$, the same two conditional-law conclusions hold if, at each fixed $(s,x)$, channel $j$ is accepted on a Borel mark set $D_j(s,x)\subseteq [0,1]$, provided these sets are pairwise disjoint, have Lebesgue measure $a_j(s,x)\mathrm{/}\mathrm{\Lambda }$, and for each fixed state $x$ the region $\{(s,u)\in [0,T]\times [0,1]:u\in D_j(s,x)\}$ is measurable for $\mathcal{B}([0,T])\otimes \mathcal{B}([0,1])$. The standard cumulative intervals in the algorithm are one such choice: their endpoints are finite sums of the Borel time sections $a_k(\cdot ,x)\mathrm{/}\mathrm{\Lambda }$, so their time-mark regions are product-measurable. If $\mathrm{\Lambda }=0$, the step 5.1 constant-path case applies. This equivalent mark-set form is used by the minimum-intensity coupling construction.

Inputs. Two finite-state reaction systems $X$ and $Y$ with the same finite channel labels $j=1,\dots ,R$. Their state sets may differ, and channel $j$ has respective increments ${\nu }_j^X,{\nu }_j^Y$ and propensities $a_j^X(x),a_j^Y(y)$ as in Finite-channel reaction propensity systems and their generator. Supply an initial pair $(x_0,y_0)$, a horizon $T>0$, and an optional candidate budget $B\in {\mathbb{N}}_0\cup \{\mathrm{\infty }\}$, defaulting to $\mathrm{\infty }$, where ${\mathbb{N}}_0:=\{0,1,2,\dots \}$.

For each channel set $M_j:=\max \{{\max }_{x\in E_X}{a}_j^X(x),{\max }_{y\in E_Y}{a}_j^Y(y)\},M:={\sum }_{j=1}^R{M}_j\mathrm{.}$ These are finite because both state sets are finite. Supply mutually independent $E_n\sim \mathrm{Exp}(1)$ and $U_n\sim \mathrm{Unif}[0,1]$ for $n\ge 1$. The countable input product uses Countable Choice and Dependent Choice; the exponential law is that of Exponential and Weibull survival specifications define probability laws. For the uniform law, Assuming countable choice, $\mathcal{L}({\mathbb{R}}^n)$ is a sigma-algebra containing every elementary set and ${\lambda }_n$ is a complete measure extending elementary volume makes the set function in Lebesgue measurable sets, the family $\mathcal{L}({\mathbb{R}}^n)$, and the restricted set function ${\lambda }_n$ a measure, and A box in ${\mathbb{R}}^n$ with parameters $a_i\le b_i$ is Lebesgue measurable of measure ${\prod }_{i<n}(b_i-a_i)$, whichever of its faces are included gives its restriction to $[0,1]$ mass one and interval probabilities equal to lengths.

Procedure. Set $t=0$, $(x,y)=(x_0,y_0)$, and $n=0$. If $M=0$, return the two constant paths through $T$ with status absorbed. Otherwise define fixed mark blocks $I_j=[b_{j-1},b_j)$, where $b_0=0$ and $b_j=M^{-1}{\sum }_{k=1}^j{M}_k$; omit empty blocks when $M_j=0$, and include the endpoint $1$ in the last nonempty block. While $t<T$:

1. If $n=B$, return the pair of path prefixes with status budget-exhausted.
2. Set $\tau =E_{n+1}\mathrm{/}M$. If $t+\tau \ge T$, return both paths held through $T$ with status horizon.
3. Set $t\leftarrow t+\tau$. Find the unique nonempty block $I_j$ containing $U_{n+1}$ and set $v=(U_{n+1}-b_{j-1})\mathrm{/}(M_j\mathrm{/}M)$. Fire channel $j$ in $X$ if $v<a_j^X(x)\mathrm{/}{M}_j$, and fire it in $Y$ if $v<a_j^Y(y)\mathrm{/}{M}_j$. If both fire, record a shared event and update both states; if only one fires, record its residual event and update only that state; if neither fires, reject the candidate and leave both states fixed. Increment $n$ once for every candidate, including a rejected candidate. Endpoint marks have probability zero, so the last-block convention at $U_{n+1}=1$ does not affect the path law.

The candidate mark block for channel $j$ has length $M_j\mathrm{/}M$. Within it the two acceptance intervals start at the same endpoint and have respective lengths $a_j^X(x)\mathrm{/}M$ and $a_j^Y(y)\mathrm{/}M$. Thus their overlap has candidate intensity $\min (a_j^X(x),a_j^Y(y))$, and the two residual portions have intensities $(a_j^X(x)-a_j^Y(y){)}_{+}$ and $(a_j^Y(y)-a_j^X(x){)}_{+}$. The budget counts candidates, so an exhausted budget returns only an incomplete prefix. No logarithmic event-selection cost is claimed.

Assume Countable Choice and Dependent Choice, and let $E$ be finite with a reaction system from Finite-channel reaction propensity systems and their generator. Fix an initial state $x_0\in E$. Supply the independent random streams required by Gillespie direct method with horizon and budget statuses and Next-reaction method by internal Poisson clocks, and use unlimited reaction-event budgets. Define the reference marked process starting at $x_0$ recursively: at state $x$, if $\lambda (x)>0$, its next reaction time has law $\mathrm{Exp}(\lambda (x))$ and its channel is independent of that time with probabilities $a_j(x)\mathrm{/}\lambda (x)$; if $\lambda (x)=0$, it remains at $x$. Then both algorithms have the law of this reference process on every finite horizon. With a finite event budget, each algorithm instead returns the corresponding reference-process prefix, stopped according to its stated horizon, absorption and budget rules. The process has generator operator $(Qf)(x)={\sum }_{j:a_j(x)>0}{a}_j(x)(f(x+{\nu }_j)-f(x)),$ absorbing states are handled exactly, and for $\mathrm{\Lambda }={\max }_{x\in E}\lambda (x)$ the expected total reaction count on $[0,T]$ is at most $\mathrm{\Lambda }T$. A zero-increment channel is still a marked reaction event, but it contributes zero to the state generator.