Past exam of the further mathematical biology course of the University of Oxford

The quantities $n_i(\tau ,t)$ are cell-cycle age densities: $n_i(\tau ,t)d\tau$ counts biological cells whose time since entering phase $i$ lies in $[\tau ,\tau +d\tau ]$. Here $\tau$ is time spent in the current phase, rather than time since birth. In this two-stage approximation to the cell cycle, a biological cell passes from G1 phase to G2 phase, then completes cell division and produces two daughters entering G1 phase.

The time derivative measures change at fixed phase age. The age derivative transports biological cells towards larger $\tau$ at unit speed, because $d\tau /dt=1$. Thus each equation is an age-structured population equation. The nonnegative hazard function ${\Gamma }_i(\tau )$ means that, conditional on still being in phase $i$ at age $\tau$, the probability of leaving during the next $dt$ is ${\Gamma }_i(\tau )dt+o(dt)$. Its product with $n_i$ is the loss from that phase-age class; it is not an additional death term.

The boundary condition at $\tau =0$ specifies the incoming flux. Every G1 phase exit contributes one entrant to G2 phase, whereas every G2 phase exit contributes two entrants to G1 phase. The initial conditions specify the phase-age composition at $t=0$. Consequently integration over phase age gives the useful population balanceEach completed cell division increases the total number of biological cells by one.

The unit integrals of the profiles make $N_i(t)={\int }_0^{\infty }{n}_i(\tau ,t)d\tau$ the phase populations, and $f_i$ their normalized cell-cycle age densities. Assume the age tails vanish and the following integrals are finite. For a fixed profile define its mean exit hazard function byThese constants have units of inverse time. Integrating the age-structured population equations and using the boundary conditions givesThe incoming fluxes are $n_1(0,t)=2{\Lambda }_2{N}_2$ and $n_2(0,t)={\Lambda }_1{N}_1$, soThe separable-profile closure for an age-structured population has a qualification: arbitrary age-dependent hazard functions do not give this constant-coefficient ordinary differential equation for arbitrary initial conditions. Without fixed $f_i$, the effective exit rates are $\int {\Gamma }_i{n}_{i}d\tau /N_i$ and can vary with time. Moreover, an exact positive separation of variables with both profiles fixed forces $N_2/N_1$ to be constant through the incoming fluxes, hence selects a common exponential mode. The displayed ordinary differential equation is the necessary integrated system for such profiles; it is also exact for arbitrary populations when both hazard functions are age independent.

Write $A$ for the population matrix and assume ${\Lambda }_1,{\Lambda }_2>0$. Its characteristic polynomial isThe two eigenvalues have opposite signs. The dominant eigenvalue isand the other eigenvalue is $\nu =-({\Lambda }_1+{\Lambda }_2)-\mu <0$. A positive eigenvector for $\mu$ satisfiesFor any nonzero nonnegative population vector in this constant-coefficient system, the coefficient of the positive eigenvector is positive: a positive left eigenvector annihilates the other eigenvector and has positive pairing with that population vector. Thus $N_i(t)=C_i{e}^{\mu t}+O(e^{\nu t})$, with $C_i>0$, and $N_i^{\prime }/N_i\to \mu$. The long-time leading term therefore obeys $N_i^{\prime }=\mu N_i$.

Substituting that exponential mode into the age-structured population equation givesDefine the phase-residence survival function and its discounted integral byThen the normalized cell-cycle age densities areSince integration of $f_i^{\prime }=-(\mu +{\Gamma }_i)f_i$ gives $f_i(0)=\mu +{\Lambda }_i$, this is equivalentlyThe constants must be consistent with both incoming fluxes. PutThe boundary conditions for amplitudes $n_i=e^{\mu t}{\varphi }_i$ give ${\varphi }_1(0)=2H_2{\varphi }_2(0)$ and ${\varphi }_2(0)=H_1{\varphi }_1(0)$. Hence the Two-stage Euler-Lotka equation isIt determines $\mu$ directly from the hazard functions; the formula in terms of ${\Lambda }_i$ is not permission to choose those constants independently of the profiles. If each phase is eventually left with probability one, $H_i(0)=1$ and $H_i(\mu )$ decreases to zero, giving a unique positive root under the usual finite, non-instantaneous residence-time assumptions. More generally positive growth requires $2H_1(0)H_2(0)>1$. Nonnegative hazard functions alone do not ensure it: if one phase has zero exit hazard function, growth cannot have the asserted positive exponent. Convergence of the complete age distribution to these profiles needs the usual mixing assumptions on phase-residence times; the two-dimensional eigenvalue calculation establishes the long-time claim for the fixed-profile system.

Because the hazard functions are constant, normalization gives ${\Lambda }_1=2\alpha$ and ${\Lambda }_2=3\alpha$. The population characteristic polynomial factors asThusThese exponential distributions describe ages among currently present biological cells. Their rates differ from the phase-exit rates because exponential growth overrepresents younger entrants. The stable phase ratio is $N_2/N_1=1/2$, so the G1 phase and G2 phase fractions approach $2/3$ and $1/3$.

Increasing $\alpha$ rescales all times by $1/\alpha$: the population undergoes faster exponential growth, its doubling time is $\boxed{(\log 2)/\alpha }$, and the profile means $1/(3\alpha )$ and $1/(4\alpha )$ decrease. The phase fractions do not depend on $\alpha$. The population matrix has exact solutionso population-ratio transients decay relative to the growing term at rate $7\alpha$. This exact population formula does not require initially stationary cell-cycle age densities.

Put $k=\beta /\gamma$. The nullclines of the space-clamped FitzHugh-Nagumo model areThe first nullcline is a downward-leading cubic polynomial, shifted vertically by $I$; the second nullcline is a line through the origin with slope $k>0$. Their intersections satisfyWriting $M=(1-\alpha +{\alpha }^2)/3$, differentiation givesIf $k>M$, $h$ increases strictly, so there is exactly one equilibrium. If $k=M$, $h$ is still strictly increasing, although its derivative vanishes at one point, so the equilibrium is again unique. If $0<k<M$, defineHere $h$ has a local maximum at $u_{-}$ and a local minimum at $u_{+}$; thereforeAt either endpoint there are two distinct equilibria, one a double root. These are the saddle-node bifurcation boundaries omitted by the generic one-or-three wording. At $k=M$ and $I=h((1+\alpha )/3)$ there is one triple root. This is the Cubic FitzHugh-Nagumo nullcline criterion.

At an equilibrium $(v_{\ast },k{v}_{\ast })$, set $p=f^{\prime }(v_{\ast })$. The Jacobian matrix, its trace and its determinant areThe eigenvalues and their discriminant areThus the complete local linear stability test isIn the last two cases, $\mathcal{D}>0$ gives a node, while $\mathcal{D}<0$ gives a focus; the sign of $T$ distinguishes inward and outward trajectories. At $\mathcal{D}=0$, a nonzero repeated eigenvalue gives a degenerate node. At $p=\gamma <k$ the eigenvalues are purely imaginary, so linear stability analysis alone does not decide asymptotic stability; with the required crossing and nondegeneracy this is a Hopf bifurcation. At $p=k$ there is a zero eigenvalue, requiring a nonlinear test and generically marking the saddle-node bifurcation in the preceding answer.

For a unique simple equilibrium, $p<k$, so it can be a sink equilibrium or a source equilibrium. The number of equilibria alone does not decide stability. If $k\le \gamma$, the unique simple equilibrium is always a sink equilibrium. If $\gamma <k$, it can become a source equilibrium where $\gamma <p<k$. The same conclusion follows particularly simply if $\gamma >M$: the trace is negative at every nonsaddle equilibrium.

For three simple equilibria ordered by their $v$ coordinates, the middle one has $p>k$ and is a saddle equilibrium; the two outer ones have $p<k$. When both outer equilibria are sink equilibria, each has an attracting neighborhood, giving bistability. In the usual two-basin portrait, the stable manifold of the saddle equilibrium separates their basins. Each branch of its unstable manifold initially leaves towards a different side; its eventual destination requires a global argument. If $k\le \gamma$, both outer equilibria must be sink equilibria. If $\gamma <k$, either or both outer equilibria can instead be source equilibria; the same test $p\lessgtr \gamma$ applies separately to each. The three lower panels illustrate all three possibilities up to exchanging left and right.

The arrows in the phase plane follow directly from the nullclines: below $w=f(v)+I$ the horizontal component points right, and above it left; below $w=kv$ the vertical component points up, and above it down. Near a stable focus trajectories spiral inward; near an unstable focus they spiral outward. Near a saddle equilibrium they approach along its stable manifold and leave along its unstable manifold. These directions, together with the local tests, determine the sketches.

There is also a useful global conclusion when the unique equilibrium is a source equilibrium. Withthe leading term in $v$ is $-\beta v^4$, so sufficiently large ellipses are trapping regions. A small repelling neighborhood of the hyperbolic source equilibrium can be removed to form an invariant annulus. The Poincaré-Bendixson theorem then guarantees at least one periodic orbit in that annulus. The numerical sketch illustrates the resulting sustained oscillation; this argument does not assert uniqueness of a limit cycle. A stable equilibrium can instead represent a resting membrane potential, with perturbations producing damped relaxation; near a threshold, a sufficiently large perturbation can produce a large excursion before recovery. Global basin boundaries and coexistence of periodic orbits need more information than the number and local stability of equilibria.

There is a consistency issue in the printed model. If $w\equiv w_0$ is imposed while retaining its evolution equation, then $0=\beta v-\gamma w_0$ forces $v\equiv \gamma w_0/\beta$. Thus the full two-equation system has no nonconstant travelling wave with exactly constant $w$. Its constant solution exists only when $I=w_0-f(\gamma w_0/\beta )$. The intended wave calculation must treat the recovery variable as externally frozen, leaving the scalar reaction-diffusion equationFor this frozen-recovery reduction, letThe cubic polynomial $f$ has a local minimum at $a_{-}$ and a local maximum at $a_{+}$. Hence the conditiongives three distinct real roots $r_1<r_2<r_3$ of $g$, withThe outer roots are stable for the scalar spatially uniform dynamics and the middle root is unstable. Seek a travelling wave $v(x,t)=V(z)$, where $z=x-ct$. The travelling-wave reduction of a reaction-diffusion system givesLet $\delta =r_3-r_1$ and chooseIt has $V(-\infty )=r_3$, $V(+\infty )=r_1$ andSubstitution leavesConsequently the frozen-recovery cubic front isReflection $x\mapsto -x$ reverses the front and its speed. A nonzero speed requires $r_2\ne (r_1+r_3)/2$; equality gives a stationary front and is the Maxwell balance for a scalar reaction-diffusion front. Since the roots sum to $1+\alpha$, this stationary value is $I=w_0-f((1+\alpha )/3)$. For the simple admissible choice $I=w_0$, the roots are $0,\alpha ,1$ andThis particular front moves when $\alpha \ne 1/2$. At $\alpha =1/2$, other currents inside the displayed open interval give moving fronts.

The bacteria disperse by diffusion with coefficient $D_B$. Their additional chemotaxis flux isIt points up the nutrient concentration gradient, and its negative spatial derivative enters the balance equation. The positive $N_0$ prevents singular sensitivity at zero nutrient concentration; increasing the background nutrient concentration reduces the drift response to a fixed gradient. The logistic growth equation term $rB(1-B/K)$ represents net reproduction with low-population growth rate $r$ and carrying capacity $K$.

The nutrient undergoes diffusion with coefficient $D_N$. Its production term $a{B}^2/(B_0^2+B^2)$ is an increasing, bounded response to the bacteria, reaching half its maximum $a$ at $B=B_0$. It has the algebraic form of a Hill equation with exponent two; the equation alone does not specify a molecular mechanism. The term $dBN$ removes nutrient at a rate proportional to both bacteria abundance and nutrient concentration.

The full fluxes areThe homogeneous Neumann boundary conditions make both fluxes vanish at the endpoints, so neither bacteria nor nutrient crosses the boundaries. The initial conditions specify their starting spatial distributions. This model includes nutrient production by bacteria, rather than a separate external supply term.

Use $K$ as the population scale, $N_0$ as the nutrient concentration scale, $r^{-1}$ as the time scale, and $L$ as the length scale. The nondimensionalization isFor example, the chemotaxis term becomesDividing the population equation by $rK$ and the nutrient equation by $r{N}_0$ gives the requested system with the following six dimensionless quantities:The rescaled domain is $0<x<1$, $t>0$. Its boundary conditions and initial conditions are

For a spatially uniform steady state with nonzero bacteria abundance, $b(1-b)=0$ forces $b_{\ast }=1$. The nutrient balance then givesTo test spatially uniform linear stability, writeAt the positive equilibrium, the Jacobian matrix isIts eigenvalues are $\boxed{-1,-\mu }$, both negative, so the positive equilibrium is asymptotically stable against spatially uniform perturbations. This remains true at $\mu =1$, where a repeated eigenvalue may produce a decaying $t{e}^{-t}$ contribution. The printed stability claim must be understood in this spatially uniform sense: sufficiently strong chemotaxis can destabilize nonuniform perturbations, as the next answer shows.

Write $b=1+u$, $n=n_{\ast }+v$ and defineThe first-order perturbation equations areThe homogeneous Neumann boundary conditions select Fourier modes proportional to $\cos (m\pi x)$, with squared wavenumber $q_m=(m\pi {)}^2$ for $m=0,1,\dots$. For a mode $e^{\sigma t}(U,V)\cos (m\pi x)$, the eigenvalue problem isIts trace is $-[1+\mu +(d_b+d_n)q]<0$. Therefore a growing mode exists exactly when its determinant is negative. SetThe dispersion relation isThe constant mode has $P(0)=\mu >0$; very short wavelengths also decay because $P(q)\sim d_b{d}_n{q}^2$. For a negative interval at positive $q$, the upward-opening quadratic polynomial must have its minimum below zero at positive $q$. This givesIn particular, $\boxed{b_0>1}$ is necessary, since otherwise $\eta \le 0$ and every coefficient of $P(q)$ is positive. When $b_0>1$, a locally increased population must increase net nutrient production, so attraction to that nutrient reinforces the population increase. This feedback must overcome both diffusion and local relaxation. The calculation is an instance of chemotactic instability with logistic growth and saturating production.

The displayed continuous-wavenumber inequality is necessary on the finite domain, but an allowed mode must actually lie in the unstable band:Equivalently, the exact discrete threshold isEquality gives a marginal zero eigenvalue rather than positive growth. At the unrestricted onset, $q_c=\sqrt{\mu /(d_b{d}_n)}$, so $k_c=(\mu /(d_b{d}_n){)}^{1/4}$. On the finite interval, Discrete Neumann mode selection first destabilizes the mode minimizing the discrete threshold; it need not attain $k_c$ exactly.

To predict the early form of a perturbation, use the larger growth eigenvalueFor $\eta >0$, the square root is real. Choose $m_{\mathrm{fast}}$ maximizing ${\sigma }_{+}(q_m)$ among allowed modes with positive growth. The corresponding perturbation has wavelength $2/m_{\mathrm{fast}}$ in $x$, or $2L/m_{\mathrm{fast}}$ in $X$. Its eigenvector satisfiesso bacteria and nutrient peaks initially coincide. A mode with larger $m$ has finer spatial structure. The two signs of $\cos (m\pi x)$ are both allowed, exchanging peaks and troughs.

The illustration holds $d_b=0.002$, $d_n=0.04$, $\mu =0.5$, $b_0=2$ and $\lambda =2.5$ fixed, giving $n_{\ast }=1$ and $\eta =0.3$, and varies ${\chi }_b$. Crossing the discrete threshold makes a nonuniform Fourier mode grow; stronger chemotaxis widens the unstable band and can favor finer structure. Increasing either diffusion coefficient raises the threshold of each fixed mode. Varying $\lambda$ can have a nonmonotone effect: with $s=1+b_0^2$,is maximal at $\lambda =\mu s$, where $n_{\ast }=1$. Additional production beyond that point weakens sensitivity enough to reduce this feedback. Changes in $L$ also change the allowed dimensional wavenumbers $m\pi /L$ and all three scaled transport coefficients. These are predictions of linear stability analysis: eventual pattern amplitude, competition between growing modes and possible merging require the nonlinear system. The threshold is a necessary condition for spontaneous pattern growth from arbitrarily small perturbations, rather than a claim that every finite-amplitude nonuniform state is excluded below it.