= Solution
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> are
$$
\begin{array}{c|c|c}
r&a_r(\mathbf x)&\nu_r\\ \hline
1&\alpha_1x_1^2/V&(-2,0,0,0)\\
2&\alpha_2x_1&(0,1,0,0)\\
3&\alpha_3x_2^2/V&(0,-1,0,0)\\
4&\alpha_4V&(0,0,1,0)\\
5&\alpha_5x_3x_4/V&(0,0,-1,-1)\\
6&\alpha_6x_2&(0,0,0,1)\\
7&\alpha_7x_2x_4/V&(0,-1,0,0).
\end{array}
$$
A 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}^7a_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 to
$$
\tau=-\frac{\log U_1}{a_0},
$$
and choose the least index $r$ satisfying
$$
\sum_{j=1}^ra_j\geq U_2a_0.
$$
Then 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$.
Back to article page
