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.

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$.