Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 74 (Systems Biology)

Assume the decreasing functions are continuously differentiable so their local gains are defined. In the positive quadrant, the nullclines are $x=f(y)$ and $y=g(x)$, both decreasing. The phase plane arrows point right when $x<f(y)$ and left when $x>f(y)$; they point up when $y<g(x)$ and down when $y>g(x)$. Both axes point into the quadrant. Since $f,g$ are bounded, sufficiently large $x$ or $y$ decreases, giving a bounded forward-invariant rectangle containing all equilibria and eventually containing each trajectory. Thus mutual repression can have one attracting intersection or several intersections separating different basins of attraction.

An equilibrium obeys $x_{\ast }=f(y_{\ast })$, $y_{\ast }=g(x_{\ast })$. Its Jacobian matrix and eigenvalues areThe product is nonnegative. Hence an intersection is a stable equilibrium with two negative real eigenvalues when $f^{\prime }{g}^{\prime }<1$, a saddle equilibrium when $f^{\prime }{g}^{\prime }>1$, and nonhyperbolic when $f^{\prime }{g}^{\prime }=1$. This is the mutual repression stability criterion. To see why an intervening saddle normally accompanies multistability, eliminate $y$ and set $F(x)=f(g(x))-x$. Bounded positivity gives $F(0)>0$ and $F(x)<0$ for large $x$, whileTwo downward, stable crossings require an upward crossing between them. For simple roots that middle crossing has $f^{\prime }{g}^{\prime }>1$ and is a saddle equilibrium. Its stable manifold separates the attracting states. Thus robust multistability requires the loop gain to exceed one at an intervening equilibrium, with the stable equilibria themselves having gain below one. A tangency at gain one is the threshold where equilibria can merge; degenerate crossings require higher-order analysis. A gain exceeding one somewhere is necessary for multiple intersections, but by itself does not prove that a particular pair of nullclines has multiple stable intersections.

Define the logarithmic repression sensitivitiesBoth are negative. At equilibrium $f(y_{\ast })=x_{\ast }$ and $g(x_{\ast })=y_{\ast }$, so their product is exactly $f^{\prime }(y_{\ast })g^{\prime }(x_{\ast })$. The condition can therefore be expressed asEquivalently, the product of the magnitudes of the two negative logarithmic gains must exceed one at that saddle equilibrium. The divergence of the vector field is identically $-2$, so the Bendixson-Dulac criterion excludes periodic orbits in the positive quadrant. The relevant stable long-time states in this bounded phase plane are thus equilibria, rather than attracting limit cycles. The following sketches illustrate the single-attractor and two-attractor possibilities; the saddle in the second panel lies on the diagonal, with its stable diagonal separating the two basins.

For this equal-strength power-law mutual repression, the logarithmic gains and their product at equilibrium areThe mutual repression stability criterion immediately excludes multistability if $mn\le 1$, since $P<mn$ at every finite positive equilibrium. The common prefactor $\lambda$ gives a stronger exclusion: if either exponent is at most one, multistability is impossible for every $\lambda$, including some cases with $mn>1$.

To prove the additional cases, suppose $m\le 1<n$. At equilibrium the two expressions for $\lambda$ implyPut $a=x/(1+x)$ and $b=y/(1+y)$. Because $m\le 1$, the left side is at least $a$, so $a\le b(1-b{)}^{n-1}$. It follows thatThe maximum in the second inequality occurs at $b=2/(n+1)$. Its strict bound follows from $4n/(n+1{)}^2\le 1$, $m\le 1$, and $((n-1)/(n+1){)}^{n-1}<1$. Interchanging the components proves the case $n\le 1<m$; the case with both exponents at most one was already covered by $mn\le 1$. Every root of $F(x)=f(g(x))-x$ therefore has $F^{\prime }<0$. Since $F$ begins positive and ends negative, there is exactly one root: two downward crossings would require another crossing with nonnegative derivative between them. The unique equilibrium is stable.

For completeness, neither exponent can be relaxed above one in this exclusion. If $m>1$ and $n>1$, take $\lambda$ sufficiently large. On $x\in [\lambda /2,\lambda ]$, $g(x)=O({\lambda }^{1-m})\to 0$, so $f(g(x))=\lambda (1+o(1))$ maps that interval into itself. It has a fixed point withIt is a stable equilibrium. By exchanging the components, there is a second stable fixed point with $y\sim \lambda$ and $x\sim {\lambda }^{1-n}$. These two regions are disjoint for large $\lambda$, proving bistability for some production strength. Hence the full answer isFor example, $m=1/2,n=3$ has $mn>1$ but still a unique stable equilibrium for all $\lambda$. Thus the product-only necessary condition for strong feedback is not sufficient when both production prefactors are constrained to be equal.

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.

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.