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

Let ${\mu }_i=⟨X_i⟩$. Applying the Markov jump-process generator to either coordinate gives the exact moment equationsIn each equation the joint reaction removes one molecule, at reaction propensity function $C{X}_1{X}_2$. Because the two reactants are different chemical species, the propensity uses their product; no same-species factorial correction is needed. The exact identity $⟨X_1{X}_2⟩={\mu }_1{\mu }_2+\mathrm{Cov}(X_1,X_2)$ shows explicitly where the equations fail to close. Neglecting this covariance relative to the product of means is the small-fluctuation moment closure. It gives the mean-field joint removal of two chemical species:The difference equation $d({\mu }_1-{\mu }_2)/dt=-\beta ({\mu }_1-{\mu }_2)$ actually holds exactly as well as after closure, because both the constant production and the joint removal terms cancel. The nonlinear moment closure will be used for the subsequent phase plane and stationary linear noise approximation calculations, as requested.

Use the approximate reaction-rate equations from part (a), with axes $u={\mu }_1$ and $v={\mu }_2$. In every case the positive quadrant is forward invariant, since the drift on either zero-concentration axis points inward. The nullclines specify the arrow signs: right if $\lambda >u(\beta +Cv)$, left if the inequality is reversed; up if $\lambda >v(\beta +Cu)$, down otherwise.

For $C=0$ and $\beta >0$, the nullclines are the vertical line $u=\lambda /\beta$ and the horizontal line $v=\lambda /\beta$. The two components relax independently:Every nonstationary trajectory is a straight ray towards the unique stable equilibrium $(\lambda /\beta ,\lambda /\beta )$, whose Jacobian matrix is $-\beta I$.

For $\beta =0$ and $C>0$, both derivatives equal $\lambda -Cuv$, so the two nullclines coincide onEvery point of this hyperbola is an equilibrium. The difference $u-v=d$ is conserved, and trajectories therefore run along straight lines of slope one. Below the hyperbola they point up and right; above it they point down and left. Each such line converges to its own intersection with the hyperbola,At this point the Jacobian matrix is $-C\left(\begin{array}{cc}{v}_{\ast }& u_{\ast }\\ v_{\ast }& u_{\ast }\end{array}\right)$, with eigenvalues $0$ and $-C(u_{\ast }+v_{\ast })$. Thus the equilibrium curve attracts transversely but has no restoring drift along its tangent. Its points are not individually asymptotically stable, because arbitrarily close initial states on different constant-difference lines converge to different equilibria.

For $\beta >0$ and $C>0$, the nullclines are the decreasing curves $u=\lambda /(\beta +Cv)$ and $v=\lambda /(\beta +Cu)$, reflected in the diagonal. At an intersection, subtracting the two stationary equations gives $u=v=\mu$. Hence there is exactly one intersection, whereThe Jacobian matrix at this point isIts symmetric eigenvector $(1,1)$ has eigenvalue $-(\beta +2C\mu )$, and the antisymmetric eigenvector $(1,-1)$ has eigenvalue $-\beta$. Both are negative, so the symmetric equilibrium is a stable node. All trajectories are bounded because each component above $\lambda /\beta$ decreases. The divergence $-2\beta -C(u+v)<0$ excludes periodic orbits by the Bendixson-Dulac criterion, and the unique interior equilibrium attracts the positive-quadrant trajectories. As $\beta$ becomes small, the drift restoring the difference becomes slow, so paths first approach the near-hyperbolic nullclines and then drift towards the symmetric intersection. These sketches show the three cases printed in the PDF; the middle case is $\beta =0$, not the different cases introduced by the converted TeX.

For $\beta >0$, symmetry of the unique stationary mean-field state gives $\mu =⟨X_1⟩=⟨X_2⟩$ and the balance $\lambda =\beta \mu +C{\mu }^2$. The fraction of the molecular removal flux using the joint reaction and the mean molecular lifetime are consequentlyIn particular $0\le E<1$, $C\mu =E/\tau$, and $\beta =(1-E)/\tau$.

The total loss flux of species one is $v_1^{-}={\mu }_1(\beta +C{\mu }_2)$, whereas its production flux $v_1^{+}=\lambda$ is constant. The reaction-rate elasticities are therefore $H_{11}=1$ and $H_{12}=C{\mu }_2/(\beta +C{\mu }_2)=E$ at the symmetric state. Species two gives $H_{22}=1$, $H_{21}=E$. All chemical jumps in either component have unit absolute size, so their molecular-flux-weighted reaction event sizes satisfy $⟨r_1⟩=⟨r_2⟩=1$. With $M_{ij}=H_{ij}/{\tau }_i$, this provesThe same matrix follows from $M=-S^{-1}AS$ with $S=\mu I$ and the drift Jacobian matrix in part (b); its positive eigenvalues are the restoring rates.

For the chemical diffusion matrix, each component receives noise from one unit birth, its own unit death, and the joint unit loss. The unnormalized diagonal noise intensity isOnly the joint reaction changes both components in one event. Its stoichiometric vector is $(-1,-1)$, so $B_{12}=B_{21}=(-1)(-1)C{\mu }^2=C{\mu }^2$. The positive sign reflects simultaneous decreases, rather than independent death events. Divide by the product of the means to normalize:Thus the chemical diffusion matrix in the normalized stationary fluctuation-dissipation relation for a reaction network isThe unequal roles of $M$ and $D$ matter: $M$ describes restoration of perturbations, while $D$ records the covariance of the immediate stochastic reaction increments.

By symmetry the normalized covariance matrix has the form $\eta =\left(\begin{array}{cc}V& W\\ W& V\end{array}\right)$. Substituting the matrices from part (c) into the normalized stationary fluctuation-dissipation relation for a reaction network, $M\eta +\eta M^{\mathsf{T}}=D$, givesHence $V+EW=1/\mu$ and $W+EV=E/(2\mu )$. Solving these two equations gives the linear noise of two species with joint removal:For $E=0$ the normalized variance is $1/\mu$ and the covariance is zero, as expected for independent immigration–death processes. For $E>0$, $V=1/\mu +E^2/(2\mu (1-E^2))$ exceeds this independent value and $W<0$. The correlation coefficient isThis negative stationary correlation is compatible with positive instantaneous joint-removal noise: the slow difference mode amplifies fluctuations of opposite signs in the two populations.

The symmetric and antisymmetric eigenvectors diagonalize all three matrices. Their normalized covariance eigenvalues areThe corresponding restoring rates are $(1+E)/\tau =\beta +2C\mu$ and $(1-E)/\tau =\beta$. Thus the sum mode retains finite fluctuations, with ${\eta }_{+}\to 3/(4\mu )$, whereas the difference mode becomes arbitrarily noisy. In particular, at a fixed nonzero mean, or along $\beta \to 0$ at fixed $\lambda ,C>0$ so $\mu \to \sqrt{\lambda /C}$,The Fano factor of either population is $\mu V=(2-E^2)/(2(1-E^2))$, which diverges as $E\to 1$ regardless of how the mean is varied. Unnormalized variances are ${\mu }^{2}V$ and covariance ${\mu }^{2}W$, so their absolute magnitudes also depend on the parameter path through $\mu$.

The phase plane in part (b) explains the singularity. At $\beta =0$, the mean-field system has a whole hyperbola of equilibria, rather than an isolated restoring state. A perturbation in the species difference selects another point on that curve. In the actual stochastic reaction network, joint complex formation still leaves $X_1-X_2$ unchanged, but the two constant-rate births change it by $+1$ and $-1$, each at rate $\lambda$. Therefore exclusive joint removal has a diffusing difference mode:Its continuous-time symmetric simple random walk on the integers has no normalizable stationary distribution, so there is no finite stationary covariance at $E=1$. The divergence as $E\uparrow 1$ signals both the loss of difference-mode restoration and eventual failure of the small-fluctuation moment closure; near the limit, the computed normalized variances need not remain small even when the individual mean populations are large.