Past exam of the stochastic modelling of biological processes course of the University of Oxford / 2023 (B5.1, A15675W1)

Solutions to all three questions of the stochastic modelling of biological processes course of the University of Oxford Part B examination, sat on 8 June 2023. Source: original released PDF, bundle pages 272–275; paper pages 1–4. The solutions discuss the two readings of Question 1(c), correct the repeated noise-coordinate label in Question 2, and state the boundary qualification needed in Question 2(a).

Write $n=n_1$ and $m=n_2$. The three stoichiometric vectors are $(1,0)$, $(0,1)$ and $(-1,-1)$. The law of mass action gives the respective reaction propensity functionsIn particular, producing the second chemical species does not consume the first chemical species. Since the two reactants in the third chemical reaction are different chemical species, the number of possible pairs is $nm$.

For each Markov chain state, the chemical master equation adds incoming probability flux and subtracts outgoing probability flux:Here $p(n,m,t)=0$ if either index is negative, and the initial condition is $p(n,m,0)=1_{\{n=m=0\}}$. These conventions also give the correct chemical master equation on the coordinate axes: for example, the third chemical reaction has zero outgoing reaction propensity function when either chemical species is absent.

At the empty Markov chain state, only the constant-rate production chemical reaction is possible. Thus the first-passage time ${\tau }_1$ has an exponential distribution of rate $2$:Its expected value can be obtained by integrating this survival function, and its second moment by integrating $2t$ times the survival function:Consequently the requested expected value and variance are

Conditioning on the observed time. The native PDF defines ${\tau }_2$ as the first production time and then conditions on that production occurring at ${\tau }_2$. If the intended information is its numerical value ${\tau }_2=t$, the answer depends on $t$. We first compute this conditional probability using a conditional probability density; conditioning on an exact value of a continuous random variable means a density ratio, since $\mathbb{P}({\tau }_2=t)=0$.

Before ${\tau }_2$, no molecules of the second chemical species are present, so the joint-removal chemical reaction cannot occur. Molecules of the first chemical species arrive as a rate-$2$ Poisson process, and each independently starts a clock with a rate-$1/2$ exponential distribution for producing the first molecule of the second chemical species. Define the unnormalized probability mass functionThe chemical master equation before the first production, with production treated as killing, isMultiplying by $z^n$ and summing gives a partial differential equation for the probability generating function:Its solution, verifiable by direct differentiation, isThus $q_n(t)=e^{-2t}m(t{)}^n/n!$, the survival function is $S(t)=F(1,t)=e^{-2t+m(t)}$, and the conditional distribution of $A_1(t)$ given survival to $t$ is a Poisson distribution of mean $m(t)$.

Production in $(t,t+dt)$ from Markov chain state $n$ has joint probability $(n/2)q_n(t)dt+o(dt)$. Therefore the probability density function of ${\tau }_2$ isThe first production leaves $A_1$ unchanged. Taking the density ratio gives a size-biased distribution of the surviving count:In particular, under this reading the requested conditional probability isThe shift by one matters: an observation made exactly when production occurs preferentially selects larger counts, since their production reaction propensity function is larger.

Without observing the numerical time. If the conditioning merely repeats the definition of ${\tau }_2$, the requested probability is instead the unconditional count at that stopping time. Before production, $A_1$ increases by one at rate $2$, while production occurs at rate $n/2$ in Markov chain state $n$. The competing exponential clocks implyThe first factor requires a second arrival before production from count one; the second requires production before a third arrival. The first arrival itself is certain. Production eventually occurs with probability one, as is also clear from $S(t)\to 0$. As a consistency check, averaging the time-conditioned answer against $f_{{\tau }_2}$ givesThese are two distinct, explicitly specified interpretations of the printed conditioning.

A product-Poisson stationary distribution is suggested by balancing the three reaction fluxes around the directed cycle of reaction complexes. If their mean counts are $a,b$, this complex balance condition iswhich gives $a=4$ and $b=1/2$. This suggests the joint probability mass functionIt is normalized because it is the product of two Poisson distributions. We verify it directly in the chemical master equation, rather than assuming a product-Poisson stationary distribution theorem. The three incoming terms satisfyWith negative-index values set to zero, these identities also hold on the axes. Their sum equals $(2+n/2+nm)\pi (n,m)$, precisely the outgoing probability flux. The stationary average total exit rate is $2+{\mathbb{E}}_{\pi }n/2+{\mathbb{E}}_{\pi }(nm)=6<\infty$, so these stationary balance equations do define an invariant distribution of a continuous-time Markov chain for the nonexplosive model verified below. This is complex balance; it does not assert pairwise detailed balance for the irreversible chemical reactions.

To identify this invariant distribution with the limiting probability distribution in the question, first check that the model is a nonexplosive continuous-time Markov chain. On any finite time interval, $A_1$ is bounded by its initial count plus the finite number of arrivals in a rate-$2$ Poisson process. Conditional on that bound, the rate of producing $A_2$ is bounded, so only finitely many such productions occur. The number of joint removals cannot exceed the available molecules. Thus only finitely many chemical reactions occur in finite time.

The continuous-time Markov chain is also an irreducible Markov chain. From a Markov chain state with both counts positive, repeated joint removals reduce both. On the $A_2=0$ axis, producing one $A_2$ and then removing a pair reduces $A_1$ by one. On the $A_1=0$ axis, adding one $A_1$ and then removing a pair reduces $A_2$ by one. Each prescribed finite sequence has positive probability, so every Markov chain state reaches $(0,0)$. Conversely, arrivals followed by productions reach any $(n,m)$ with $n>0$; to reach $(0,m)$, reach $(1,m+1)$ and remove one pair. The continuous-time Markov chain is therefore positive recurrent, and the continuous-time Markov chain convergence theorem gives $\varphi =\pi$ from the empty initial condition.

The two Poisson factors happen to have cancelling numerical coefficients at $(2,2)$:

The native PDF repeats $d{W}_y$ in its list of independent noise terms, whereas its third stochastic differential equation uses $d{W}_z$. We use the intended three independent coordinate Brownian motions $W_x,W_y,W_z$.

Put $\mathbf{X}=(X,Y,Z)$, $r=|\mathbf{X}|$ and $S=r^2$. Before the origin is reached, Itô's formula givesThe $6$ is the sum of three quadratic variation terms, each $2dt$. The drift contribution is $-2\alpha dt$ because the numerator is $|\mathbf{X}{|}^2$. The radial stochastic integral can be written using a scalar Brownian motion $B$:Thus $S/2$ is a squared Bessel process of dimension $\delta =3-\alpha$ until it reaches zero. If this equation is taken globally with the usual nonnegative squared Bessel process continuation, taking expected values givesThe stochastic integral has zero expected value: stopping first at large levels bounds the integrands, yields a finite first moment for $S$, and then makes $\mathbb{E}{\int }_0^{t}S(u)du$ finite, so the unrestricted integral is a square-integrable martingale.

Boundary qualification. The punctured domain in the question does not itself specify that continuation. The boundary classification of a Bessel process says that the origin is inaccessible from $r_0>0$ when $\delta \ge 2$, equivalently $\alpha \le 1$. In this parameter range the original stochastic differential equation is global in the stated domain and the boxed result holds as written.

For $1<\alpha <3$, we have $0<\delta <2$: the origin is reached in finite time with probability one. One can also see accessibility directly from the radial generator of an SDE diffusionFor $r<b$, the harmonic function $1-(r/b{)}^{\alpha -1}$ is the probability of reaching zero before $b$. Finite-time access follows by using the nonnegative functionin the stopped Itô formula; it bounds the expected value of the time to reach zero or $b$. Letting $b\to \infty$ gives almost sure finite hitting of zero. The coefficients of the original vector stochastic differential equation are undefined there, so its stated punctured-domain version has no prescribed long-time evolution in this range.

The natural continuing radial squared Bessel process gives $6-2\alpha$, but stopping and holding the particle at the origin would give a different answer. Indeed, with $T_0$ the first hit of zero, Itô's formula for the stopped process givesSince $T_0<\infty$ almost surely, dominated convergence applied to $(t\wedge T_0)/t\le 1$ givesTherefore $6-2\alpha$ is the usual continued-radial answer; for $1<\alpha <3$ a boundary convention is necessary. It is not a uniquely defined limit for the punctured-domain model alone.

At $\alpha =0$ the transverse coordinates are independent Brownian motions with diffusion coefficient $D=1$. The cylinder is infinitely long, so the axial coordinate cannot affect its first contact time. By rotational invariance of Brownian motion, the mean first-passage time $u$ depends only on the transverse radius $\rho =\sqrt{x^2+y^2}$.

The Kolmogorov backward equation for a mean first-passage time is $Lu=-1$. Here $L={\partial }_x^2+{\partial }_y^2$, so the radial boundary value problem isThe absorbing boundary condition for a diffusion is imposed on the circular cross-section. At $\rho =0$, the solution must extend smoothly over the axis. Multiplying by $\rho$ and integrating givesSmoothness at the axis excludes the logarithm term, and $u(R)=0$ fixes $C_0=R^2/4$. ConsequentlyThis also verifies the finiteness assumption behind the Kolmogorov backward equation. Let $T$ be the first contact time and apply the stopped Itô formula to the bounded nonnegative function $u$ just obtained. Its bounded derivatives on the disk make the stopped stochastic integral a martingale, givingBy monotone convergence, $T$ has finite expected value; continuity and dominated convergence then make the final term vanish as $t\to \infty$. Hence $u({\hat{r}}_0)=\mathbb{E}T$. In particular, the mean first-passage time is $R^2/4$ on the axis and tends to zero at the adsorbing wall.

For a radial function $f(r)$ in three dimensions, the Laplacian is $f^{\prime \prime }+2f^{\prime }/r$, while the imposed radial drift contributes $-\alpha f^{\prime }/r$. Thus the radial generator of an SDE diffusion isThe hitting probability for the inner sphere is harmonic for this diffusion generator in the shell, with value one on the inner sphere and zero on the outer sphere. At $\alpha =1$, the boundary value problem is thereforeIntegrating $(r{p}^{\prime }{)}^{\prime }=0$ gives $p(r)=A\log r+B$. The two boundary conditions determineBoth arguments of the logarithms are dimensionless. The result decreases from one to zero across the shell and lies strictly between zero and one for $R_1<r_0<R_2$.

For a direct probability interpretation, stop at the first contact $T$ with either sphere. The stopped process $p(r(t\wedge T))$ is a bounded martingale by Itô's formula. The shell has positive inner radius, so its drift is bounded and its diffusion is nondegenerate; exit occurs almost surely. Optional stopping and dominated convergence therefore giveThe logarithm appears because the inward drift cancels one of the two geometric $f^{\prime }/r$ terms in the three-dimensional radial Laplacian, leaving the same radial diffusion generator as a planar Brownian motion with diffusion coefficient one.

We use the equilibrium generalized Langevin equation and its generalized fluctuation-dissipation relation, with the factor $5$ in the noise autocorrelation fixing the equilibrium velocity variance. Let $g$ be the velocity relaxation function and let a tilde denote a one-sided Laplace transform:The equilibrium velocity autocorrelation function is $\chi (t)=5g(t)$, soFor completeness, this relation follows without assuming independence of the noise and velocity in a stationary infinite history. The two-sided noise power spectral density is $S_R(\omega )=10\mathrm{Re}\tilde{\kappa }(i\omega )$. The stable linear equation has velocity transfer function $H(i\omega )=[i\omega +\tilde{\kappa }(i\omega ){]}^{-1}$, soThis is precisely the power spectral density of the even extension of $5g(t)$. Stable relaxation removes initial transients with a finite second moment, giving the long-time autocorrelation function defined in the question. All the parameter choices used below have stable relaxation.

Integrating the velocity autocorrelation function twice gives the mean squared displacement from a fixed initial position:A decaying initial transient does not change its asymptotic slope. Therefore the Green-Kubo formula for diffusion givesFor the given friction memory kernel, direct integration yieldsAt $c_1=c_2=c_3=1$, the integrated friction is $2$, henceThe relaxation denominator here is $s^3+2s^2+2s+2$. Its roots have negative real parts by the cubic Routh-Hurwitz criterion, since all coefficients are positive and $2\cdot 2>2$. This justifies the decay and integrability used above. The diffusion coefficient depends on the area under the friction memory kernel, not just its value at zero.

For $c_1=2$, $c_2=4$, $c_3=5$, the friction memory kernel and its Laplace transform areUsing the generalized fluctuation-dissipation relation with equilibrium velocity variance $5$, the velocity autocorrelation function has Laplace transformAll three poles are negative, so the stationary relaxation is stable. A partial fraction decomposition givesTaking the inverse Laplace transform givesThe negative middle coefficient is compatible with an autocorrelation function; individual terms in a partial fraction decomposition need not themselves be separate physical modes with positive weights. Several useful checks follow directly:The Volterra integral equation for $\chi$ predicts exactly ${\chi }^{\prime }(0)=0$ and ${\chi }^{\prime \prime }(0)=-\kappa (0)\chi (0)=-10$. Its integral is also $8-10/9+1/18=125/18=5/\tilde{\kappa }(0)$, agreeing with the Green-Kubo formula for diffusion.

Write $c=c_1=c_3$, so the friction memory kernel is ${\kappa }_c(t)=c{e}^{-ct}$ and has area one for every $c>0$. Its Laplace transform is ${\tilde{\kappa }}_c(s)=c/(s+c)$. The generalized fluctuation-dissipation relation therefore givesTo verify the pointwise time-domain limit, rather than relying only on Laplace transform convergence, consider $c>4$. The two rootsare negative and satisfy ${\lambda }_{+}\to -1$, ${\lambda }_{-}\to -\infty$. Since ${\chi }_c(0)=5$ and ${\chi }_c^{\prime }(0)=0$, the inverse Laplace transform isThe coefficient of the slow exponential tends to $5$, and that of the fast exponential tends to zero. HenceThis is the velocity autocorrelation function of an Ornstein-Uhlenbeck process satisfying $dV=-Vdt+\sqrt{10}dW$. The friction memory kernel becomes local in time. Consistently, the even noise autocorrelation $5c{e}^{-c|\tau |}$ tends to $10\delta (\tau )$ in the sense of distributions, fixing the white noise normalization. The finite-$c$ right derivative ${\chi }_c^{\prime }(0)=0$ does not converge to the right derivative $-5$ of the limiting function: there is a short initial relaxation layer.

The integrated friction is exactlyFor every finite $c>0$, the roots of $s^2+cs+c$ have negative real parts. The Green-Kubo formula for diffusion therefore gives the diffusion coefficient directly:ConsequentlyIntegrating the limiting velocity autocorrelation function gives the same value, ${\int }_0^{\infty }5e^{-\tau }d\tau =5$. Computing $D_c$ before taking the limit also avoids having to justify an interchange of the infinite-time integral and the parameter limit.