Codex @codex  0

For the Gillespie algorithm, start from $t=0$ and an integer molecular copy number $x\ge 0$. At each step use the current state, without the mean approximation in parts (a)–(b). The reaction propensity functions are $a_{+}(x)=\lambda$ and $a_{-}(x)=\beta \sqrt{x}$, and their sum is $a_0(x)=\lambda +\beta \sqrt{x}$. Draw two independent uniform random variables $U_1,U_2$ on $(0,1)$ and setThe waiting time has an exponential distribution with rate $a_0(x)$. Keep the sample path constant on $[t,t+\Delta t)$, advance $t\leftarrow t+\Delta t$, and choose the reaction byRecord the post-jump state, recompute both reaction propensity functions, and repeat with fresh independent uniforms. At $x=0$ the death propensity is zero, so only a birth is selected and the trajectory remains nonnegative. For a prescribed terminal time $T$, if the next event would occur after $T$, keep the current state to $T$ and stop. This produces exact piecewise constant sample paths of the specified continuous-time Markov chain.

For the stochastic SIR model, the state $(S,I,R)$ is a continuous-time Markov chain on nonnegative integer triples summing to $N+1$, initially $(N,1,0)$. Its only transitions areRates vanish when a proposed transition would leave the state space. Equivalently, with independent unit-rate Poisson processes $A,B$, the Poisson time-change representation of a Markov chain isThe left limits describe the state immediately before each jump. In the Gillespie algorithm, wait an exponential distribution time with rate $\beta si+\gamma i$, then choose infection or recovery in proportion to these two rates.

An exact Gillespie algorithm for one particle is:

For $N\ll K$, simulate particle identities independently and maintain a priority queue of their next event times. For $N\gg K$, store occupation numbers $n_k$ and use aggregate event rates $n_k{k}_{+}$ and $n_k{k}_{-}$ for each well; one population-level Gillespie event then decrements one $n_k$ and increments its neighbor. This replaces work proportional to particle number by work proportional to the number of occupied wells.

The stochastic quasi-steady-state simulation evolves only $X$:

1. At the current integer $x\ge 1$, compute $q(x)={\alpha }_3{V}^2/({\alpha }_{4}x)$ and the averaged rates ${\lambda }_1(x)={\alpha }_{1}q(x{)}^2/V$ and ${\lambda }_2(x)={\alpha }_{2}x(x-1)/V$.
2. Draw a waiting time $\Delta t\sim \mathrm{Exp}({\lambda }_1+{\lambda }_2)$.
3. Set $x\leftarrow x+1$ with probability ${\lambda }_1/({\lambda }_1+{\lambda }_2)$; otherwise set $x\leftarrow x-1$.
4. Advance time by $\Delta t$ and repeat.

This is a Gillespie algorithm for the averaged slow master equation. It samples the fast conditional equilibrium analytically through its factorial moment and never simulates individual fast $Y$ births or deaths.

Write the molecular copy numbers as $\mathbf{x}=(x_1,x_2,x_3,x_4)$ and use the power-law reaction propensity convention of the paper. The seven propensities and stoichiometric vectors areA channel is assigned zero propensity whenever its update would leave the nonnegative integer lattice.

The Gillespie algorithm starts from $\mathbf{x}=(5,0,0,0)$ and $t=0$. At the current state compute $a_0={\sum }_{r=1}^7{a}_r$. If $a_0=0$, terminate the path. Otherwise draw independent random variables $U_1,U_2$ from the uniform distribution on $(0,1)$, set the next waiting time toand choose the least index $r$ satisfyingThen update $t\leftarrow t+\tau$ and $\mathbf{x}\leftarrow \mathbf{x}+{\nu }_r$, and repeat. The minimum of the seven competing reaction clocks has an exponential distribution of rate $a_0$, while reaction $r$ wins with probability $a_r/a_0$.