Past exam of the mathematics course of the University of Cambridge / 2024 / iii / Paper 356 (Stochastic Processes in Biology)

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

For a function $h$ on the copy-number state space, define the shift operatorWith the propensities from part (a), the forward operator of a Markov jump process isEquivalently, the chemical master equation has the gain-minus-loss formwhere terms outside the state space vanish.

The adjoint Markov jump-process generator acts on observables:Indeed, shifting the summation index in the gain terms gives $⟨f,{\mathcal{L}}^{\ast }p⟩=⟨\mathcal{L}f,p⟩$.

Only reaction 1 changes $X_1$. Starting from five molecules, its successive states are $5\to 3\to 1$, with rates $5^2=25$ and $3^2=9$; state one is then absorbing for this channel. ThusThe probability of still being at three is the convolution of the first waiting time with an exponential distribution and survival of the second:Applying the law of total probability to the three possible states givesIt has $g(0)=0$ and tends to one as both reaction waiting times elapse.

For any observable $f$, the Markov jump-process generator identity givesSet $m(t)=⟨x_2⟩$ and $s(t)=⟨x_2^2⟩$. With ${\alpha }_7=0$, reaction 2 changes $x_2$ by $+1$ and reaction 3 changes it by $-1$. Thereforeand, using $(x_2+1{)}^2-x_2^2=2x_2+1$ and $(x_2-1{)}^2-x_2^2=-2x_2+1$,The equation for the second moment contains the third, whose equation contains the fourth, and so on. This is an infinite moment hierarchy.

Reaction 1 eventually leaves the odd initial copy number at $x_1^{\ast }=1$. PutThe stationary first-moment equation givesThe stationary second-moment equation then givesUnder the prescribed central-moment closure,Substitution yields the requested polynomial equationThis cubic comes from the closure approximation; it is not an exact equation for the stationary mean.

Reaction 4 creates $X_3$, reaction 5 removes one $X_3$ and one $X_4$, and reaction 6 creates $X_4$; reaction 7 leaves $X_4$ unchanged. The exact first-moment equations are thereforeAt the assumed unique stationary distribution, both left-hand sides vanish. Eliminating the common mixed moment givesNo moment closure or independence assumption is involved.

The odd $X_1$ copy number decreases by two until it reaches one, so $⟨x_1(t)⟩\to 1$. The exact $X_2$ first-moment equation isBecause copy numbers are nonnegative integers, $x_2^2\ge x_2$, and the last mixed moment is nonnegative. HenceThe integrating factor $e^{{\alpha }_{3}t}$, together with $⟨x_1(t)⟩\to 1$, gives the comparison bound

Subtracting the exact $X_4$ first-moment equation from the $X_3$ equation cancels reaction 5:The assumed inequality is preciselyChoose a positive $\epsilon$ smaller than this gap divided by ${\alpha }_6$. The limit superior bound implies that, for all sufficiently large $t$,Thus $⟨x_3⟩-⟨x_4⟩$ grows at least linearly. Since $⟨x_4⟩\ge 0$,so the required species index is $\boxed{i=3}$.

The two-dimensional Fokker-Planck equation iswith Fokker-Planck probability currentThe point initial condition iswhere $\delta$ is the Dirac delta function. The reflecting boundary condition for a diffusion iswith $\mathbf{n}$ the outward unit normal. This zero-flux condition conserves the integral of the probability density function over the square.

Let ${\xi }_{1,n},{\xi }_{2,n}$ be independent random variables with the standard normal distribution. The unconstrained Milstein method step isThe $X$ noise is additive, so its Milstein correction vanishes. The diagonal, independent noise fields commute, so no cross iterated stochastic integral is required.

Implement each reflecting boundary condition for a diffusion by folding the proposed coordinate back into $[-1,1]$. One formula that also handles multiple overshoots isThe complete update isUnder the usual smoothness assumptions, the Milstein discretization has strong order of convergence one and weak order of convergence one. The reflection enforces the boundary pathwise.

The adsorption mechanism is uniform along the two vertical sides, and $X$ evolves independently of $Y$. Consequently the mean first-passage time depends only on the initial $x$-coordinate. Its Kolmogorov backward equation for $a_1=1$ is

The forward partially absorbing boundary condition for a diffusion is $\mathbf{J}\cdot \mathbf{n}=\kappa p$. The boundary term in the adjoint relation isbecause the $X$ diffusion coefficient is one. It vanishes for every admissible $p$ precisely whenFor $\kappa =1$, the two Robin boundary conditions are therefore

The general solution of the ordinary differential equation isThe right condition gives $C=2$, and the left gives $D=-2/e$. At the prescribed initial position,

With an absorbing boundary condition for a diffusion on each vertical side, the splitting probability $u(x)$ of hitting the left side before the right is harmonic for the Kolmogorov backward equation:For $a_1\ne 0$,The particle starts at $x=0$, soThe continuous limit at zero is $g(0)=1/2$, as required by reflection symmetry. A large positive drift drives the particle toward the right and gives $g(a_1)\to 0$; a large negative drift drives it toward the left and gives $g(a_1)\to 1$.