Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 74 / 2

Let $\mu =⟨X⟩$. The Markov jump-process generator applied to $X$ gives the exact moment equation $\dot{\mu }=3\lambda -\beta ⟨\sqrt{X}⟩$. With the prescribed moment closure, $⟨\sqrt{X}⟩\approx \sqrt{\mu }$, stationarity givesThe production flux counts molecules, not reaction events: it is $v^{+}=3\lambda$, while the removal flux is $v^{-}=\beta \sqrt{\mu }=3\lambda$. The reaction-rate elasticity of removal relative to constant production isThe mean molecular lifetime is population divided by its stationary molecular removal flux. ThereforeThe noise-relevant average chemical event size is the molecular-flux-weighted reaction event size. For the stoichiometric vectors $+3$ and $-1$ it isThis weights each absolute jump by the molecular turnover it contributes. Equivalently, the balanced production and removal fluxes give equal weight to their sizes $3$ and $1$, yielding $(3+1)/2$. An event-count-weighted average would be $3/2$, since three quarters of firing events are deaths at stationarity; that different average does not determine the diffusion term in the linear noise approximation.

The normalized stationary fluctuation-dissipation relation for a reaction network isHere $M$ is the negative drift Jacobian matrix expressed in normalized population coordinates, and $D$ is the normalized noise matrix. More explicitly, with $S=\mathrm{diag}({\mu }_i)$, drift Jacobian $A$, stoichiometric vectors ${\nu }_r$, and stationary reaction propensity functions $a_r$, these matrices are $M=-S^{-1}AS$ and $D=S^{-1}B{S}^{-1}$, where $B={\sum }_r{a}_r{\nu }_r{\nu }_r^{\mathsf{T}}$. The linear noise approximation gives a covariance evolution $\dot{\eta }=-M\eta -\eta M^{\mathsf{T}}+D$; stationarity is the displayed Continuous Lyapunov equation.

In one dimension, the normalized restoring coefficient is $M=H/\tau =1/(2\tau )$. Squared chemical jumps determine the noise:The factor nine from a birth burst is essential; its three molecules arise in one correlated reaction event. The scalar Continuous Lyapunov equation is $2M\eta =D$, givingEquivalently, the normalized variance is $⟨s⟩/(H\mu )=2/((1/2)\mu )$. The Fano factor is approximately four. This calculation uses the prescribed mean closure and the linear noise approximation, so its small-relative-fluctuation regime is $\mu \gg 1$; it is not an exact assertion for arbitrarily small molecular copy numbers.

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.