Critical-reaction binomial guarded tau-leaping

Inputs. A finite state reaction system Finite-channel reaction propensity systems and their generator with nonzero jump vectors, initial state $x_0\in E$, horizon $T>0$, maximum leap ${\tau }_{\max }>0$, critical threshold $n_c\in \mathbb{N}$, finite step budget $B$ that is a nonnegative integer, and independent ideal random inputs for the procedure. Use a fresh independent direct-method stream on every critical SSA call. Also use mutually independent unit-uniform streams $U_{q,j,r}\sim \mathrm{Unif}[0,1]$, independent of all SSA streams, for leap indices $q$, channels $j$, and positive trial indices $r=1,2,\dots$. A binomial draw with parameters $(m_j,p_j)$ is the number of the first $m_j$ entries in its stream that are less than $p_j$. The countable product of these input laws exists under Countable Choice and Dependent Choice by Assuming countable and dependent choice, countable products of arbitrary probability spaces. The first-event law of the critical SSA subroutine is the one in Finite-state Gillespie and next-reaction methods have the specified CTMC law. For a state $x$ and channel $j$, define its admissible firing capacity ${\mathbb{N}}_0:=\{0,1,2,\dots \},m_j(x):=\max \{m\in {\mathbb{N}}_0:x+r{\nu }_j\in E\text{ for every integer }r\in \{0,\dots ,m\}\}\mathrm{.}$ Finiteness of $E$ and ${\nu }_j\mathrm{\ne }0$ make this maximum finite. Let $C(x)=\{j:m_j(x)\le n_c\}$ be the critical set. When $C(x)\mathrm{\ne }\mathrm{\varnothing }$, hold $C(x)$ fixed and form the finite reaction system on $E$ using only those channels, with their original propensities. If $C(x)=\mathrm{\varnothing }$, do not form a restricted system and define the critical result to be no event. In the nonempty case, the critical SSA subroutine runs the direct Gillespie method on this restricted system with event budget one and horizon $\tau$; it returns its first event $(s,j_{\ast })$ if it occurs strictly before $\tau$, where $0<s<\tau$, and otherwise returns no event.

Output. The recorded states at accepted leap endpoints, with piecewise-constant interpolation, and a status in {completed, absorbed, invalid-state, inconclusive, step-size-infeasible, critical-interface-required}. A noncompleted step-size, interface, or budget failure returns the current path prefix. Absorption and invalid-state handling record the held state through $T$.

Procedure. Set $t=0$, $x=x_0$, step counter $q=0$, recorded path $P=[(0,x_0)]$, and status completed.

1. If $t=T$, return $P$ with status completed. If $x\mathrm{\notin }E$, append $(t,x)$ to $P$ if needed, append $(T,x)$ if $t<T$, and return $P$ with status invalid-state. If $q=B$ before $T$, return the current path prefix with status inconclusive. If $\lambda (x)=0$, hold $x$ fixed through $T$, append $(T,x)$ to $P$ if needed, and return $P$ with status absorbed.
2. Compute $C(x)$ and $\tau =\min \text{\hspace{0.17em}⁣}(T-t,{\tau }_{\max },{\min }_{\begin{array}{c}j\mathrm{\notin }C(x)\\ a_j(x)>0\end{array}}\frac{m_j(x)}{a_j(x)}),$ where the inner minimum over an empty set is $+\mathrm{\infty }$. If $\tau \le 0$, return step-size-infeasible with the current path prefix. If $C(x)=\mathrm{\varnothing }$, record that there is no critical event and skip the critical SSA call. If $C(x)\mathrm{\ne }\mathrm{\varnothing }$, call critical SSA for horizon $\tau$; it returns either no critical event before $\tau$, or a pair $(s,j_{\ast })$ with $0<s<\tau$ and $j_{\ast }\in C(x)$. Set $\delta$ to $s$ in the event case and to $\tau$ otherwise.
3. Starting from $\tilde{x}=x$, process the noncritical channels in increasing index order. Recompute $m_j(\tilde{x})$. If $a_j(\tilde{x})>0$ and $m_j(\tilde{x})=0$, return critical-interface-required with the current path prefix. Otherwise set $K_j=0$ when $m_j(\tilde{x})=0=a_j(\tilde{x})$. When $m_j(\tilde{x})>0$, set $p_j=a_j(\tilde{x})\delta \mathrm{/}{m}_j(\tilde{x})$; if $p_j>1$, return step-size-infeasible with the current path prefix and without rejecting and resampling. Set $K_j={\sum }_{r=1}^{m_j(\tilde{x})}{\mathbf{1}}_{\{U_{q,j,r}<p_j\}},\tilde{x}\leftarrow \tilde{x}+K_j{\nu }_j\mathrm{.}$
4. If critical SSA returned $(s,j_{\ast })$, return inconclusive with the current path prefix when $\tilde{x}+{\nu }_{j_{\ast }}\mathrm{\notin }E$; otherwise apply that critical reaction to $\tilde{x}$. If the call returned no event, leave $\tilde{x}$ unchanged. Set $x\leftarrow \tilde{x}$, set $t\leftarrow t+\delta$, append $(t,x)$ to $P$, set $q\leftarrow q+1$, and repeat.