Past exam of the mathematics and the environment and mathematical physiology course of the University of Oxford / 2010 / B

Use $S,C_i,P$ for concentrations and define net complex-formation rates $R_i=k_i{C}_{i-1}S-(k_{-i}+k_i^{+})C_i$, with $R_{n+1}=0$. The law of mass action givesEach binding event removes one free substrate molecule; each unbinding returns one; each catalytic event releases one reaction product and reduces the number of occupied sites by one. Suitable initially unbound enzyme and substrate data are $S(0)=S_0$, $P(0)=0$, $C_0(0)=E_0$, $C_i(0)=0$ for $i\ge 1$. All rates are nonnegative on the boundary of the nonnegative concentration region in the direction required to preserve it. These equations also conserve $S+P+{\sum }_{i=1}^{n}i{C}_i=S_0$, accounting for bound substrate as well as free substrate.

Adding the ordinary differential equations for all enzyme states makes the net rates telescope:With initially unbound enzyme, this sum is $E_0$ for all time. HenceThis is enzyme conservation, not conservation of each complex individually. The reaction changes site occupancy and makes reaction product, but no event creates or destroys an enzyme molecule.

Besides the stated scales, set $P=S_{0}p$, $ϵ=E_0/S_0$, ${\alpha }_i=k_i/k_1$, ${\beta }_i=(k_{-i}+k_i^{+})/(k_1{S}_0)$, and ${\kappa }_i=k_i^{+}/(k_1{S}_0)$. Then $r_i={\alpha }_i{c}_{i-1}s-{\beta }_i{c}_i$, and direct substitution, remembering $d/dt=k_1{E}_{0}d/d\tau$, givesThe exact initial values are $s(0)=1$, $p(0)=0$, $c_0(0)=1$, $c_i(0)=0$ for $1\le i\le n$. The enzyme-state equations carry the small parameter because the enzyme concentration scale is small relative to the initial substrate concentration; the substrate evolves on the chosen slower depletion time scale.

For an outer solution with bounded derivatives as $ϵ\to 0$, the final complex equation gives $r_n=0$. Working backwards in $0=r_i-r_{i+1}$ then gives every $r_i=0$, including $r_1$; setting only their sum to zero would not suffice. This is the quasi-steady-state approximation, valid after the initial complex-formation layer with the relevant dimensionless rate parameters held fixed.

The recurrence is $C_i=S{C}_{i-1}/K_i$, where $K_i=(k_{-i}+k_i^{+})/k_i$. Thus $C_j={\varphi }_j{S}^j{C}_0$, ${\varphi }_j={\prod }_{i=1}^j{K}_i^{-1}$. Substitution into enzyme conservation yieldsEach state $j$ produces reaction product at rate $k_j^{+}{C}_j$, soThe source's uppercase $N$ in the index range is read as the number $n$ of binding sites. This approximation preserves the total enzyme population and the different catalytic rates of the occupied states.

The sequential-enzyme Hill limit uses the intended distinguished limit for $n\ge 2$ holds the other rates fixed while $k_1\to 0$, $k_n\to \infty$, and $k_1{k}_n$ tends to a positive finite constant. Then for every $1\le j<n$, ${\varphi }_j={\prod }_{i=1}^j{k}_i/(k_{-i}+k_i^{+})=O(k_1)\to 0$, whereas ${\varphi }_n$ remains finite and nonzero. Hence intermediate complexes have negligible occupancy relative to the unbound and fully bound states.

In the quasi-steady-state approximation rate, all intermediate numerator and denominator terms disappear. Multiplying by $1/{\varphi }_n={\prod }_{i=1}^n{K}_i$ givesThis is the Hill equation with maximal rate $E_0{k}_n^{+}$ and effective half-saturation concentration $(\prod K_i{)}^{1/n}$. Here its power $n$ follows from the sequential binding mechanism and the suppression of intermediate states.

The phrase “$k_1{k}_n$ finite” needs this scaling interpretation for a nontrivial limiting rate. Merely bounding that product, while permitting unrelated rates to vary, does not ensure that all intermediate weights vanish relative to the terminal one. The derivation explicitly supplies the usual fixed-other-rates, positive-product assumptions. For $n=1$, the rate in the box is already the exact one-state quasi-steady-state approximation rate, without this two-endpoint limiting argument.

Write $c=c_1$, $\beta =(k_{-1}+k_1^{+})/(k_1{S}_0)$ and $\kappa =k_1^{+}/(k_1{S}_0)$. Since $c_0=1-c$, the exact equations areAt the prescribed initial data $s=1,c=0$, $r_1=1$, rather than zero. The quasi-steady-state approximation would instead demand $c=s/(s+\beta )$, giving $c(0^{+})=1/(1+\beta )$. This mismatch signals an initial layer.

Set $\sigma =\tau /ϵ$ and expand $s=s_0+ϵs_1+\cdots$, $c=c_0+\cdots$. At leading order $s_0^{\prime }=0$, so $s_0=1$, while $c_0^{\prime }=1-(1+\beta )c_0$ and $c_0(0)=0$. ThusFor a more informative substrate approximation, integration of $s_1^{\prime }=-1+(1+\beta -\kappa )c_0$ with $s_1(0)=0$ givesOn the slow scale $c_{\mathrm{out}}=s_{\mathrm{out}}/(s_{\mathrm{out}}+\beta )$ and $s_{\mathrm{out},\tau }=-\kappa s_{\mathrm{out}}/(s_{\mathrm{out}}+\beta )$, with matched leading value $s_{\mathrm{out}}(0)=1$. Integrating gives $s_{\mathrm{out}}+\beta \log s_{\mathrm{out}}=1-\kappa \tau$. In the overlap $1\ll \sigma \ll 1/ϵ$, the inner complex tends to $1/(1+\beta )$ and its substrate slope matches the expansion of this outer solution. A leading composite asymptotic expansion isIt satisfies both initial values at leading order and smoothly joins the initial complex build-up to substrate depletion. The displayed $O(ϵ)$ inner substrate term records the small binding-induced initial slip if greater accuracy is wanted.

The nullclines are the cubic $w=f(v)$ and the straight line $w=\gamma v$. At an equilibrium either $v=0$, or $\gamma =(a-v)(v-1)$. The latter quadratic has maximum $(1-a{)}^2/4$, so $\gamma >(1-a{)}^2/4$ leaves only the resting state $(0,0)$. Its Jacobian matrix isBoth eigenvalues have negative real parts; the rest state is linearly stable and asymptotically stable.

On the fast time $t/ϵ$, the recovery variable $w$ is approximately fixed. At $w\simeq 0$, the scalar fast equation has stable roots $v=0,1$ and an unstable root $v=a$. A small kick below $a$ decays towards rest; a kick above it moves rapidly towards the right attracting branch. For finite $ϵ$, the threshold is displaced from the frozen root by slow recovery, but the leading threshold mechanism is this unstable middle branch.

After the fast excitation, $w_t=\gamma v-w>0$ raises $w$ while $v$ follows the right attracting branch of the cubic. At its upper fold the branch ends and $v$ makes a fast downward jump to the left branch. Recovery then lowers $w$ slowly, with a negative undershoot in $v$, before return towards $(0,0)$. Thus the action potential time trace has rapid depolarization, a slower high-voltage excursion, rapid repolarization, and slow recovery. On the fast segments $w$ scarcely changes; on slow segments $w\simeq f(v)$ and $f^{\prime }(v)<0$. The figure shows a numerically integrated superthreshold excursion and its time trace, rather than just arrows on the nullclines.

Use $W=w-I$. The nullclines become $W=f(v)$ and $W=\gamma v-I$, while the slow equation becomes $W_t=\gamma v-W-I$. The cubic has foldswith $W_{\ast }=f(v_{\ell })$ its local minimum and $W^{\ast }=f(v_r)$ its local maximum. Let $V_{-}(W)$ and $V_{+}(W)$ be the left and right attracting roots of $f(v)=W$. A line whose intersection is on the repelling middle branch permits a relaxation oscillation: upward slow travel on the right branch from $W_{\ast }$ to $W^{\ast }$, a fast jump to the left branch, downward slow travel to $W_{\ast }$, and a return jump. For sufficiently large fixed $\gamma$ this occurs for currents between the fold values $\gamma v_{\ell }-f(v_{\ell })$ and $\gamma v_r-f(v_r)$, away from the small-$ϵ$ onset layers. The straight line then crosses each outer branch beyond the relevant slow excursion, so the slow velocities have the required signs.

Neglecting the fast-jump durations, integrate $dt=dW/(\gamma V_{\pm }(W)-W-I)$ along the two slow segments. Their opposite orientations giveThe printed period omits the current shift. With the source's definition that $V_{\pm }$ solve $f(v)=w$, both denominators must retain $-I$ in the shifted coordinate. Equivalently, in the physical recovery coordinate the bounds are $I+W_{\ast },I+W^{\ast }$ and the integrand uses $\gamma V_{\pm }(w-I)-w$. The printed formula is therefore not the leading slow period for the displayed $I>0$ equations unless its notation is changed or an additional approximation justifying the omitted term is supplied.

For $z=ct-x$ the chain rule gives ${\partial }_t=cd/dz$, ${\partial }_x=-d/dz$, ${\partial }_{xx}=d^2/d{z}^2$. ThusAs $z$ increases from the material ahead of the wave towards that behind it, the singular pulse construction is $A\to B\to C\to D\to A$. The front $AB$ raises $v$ at approximately fixed $w=0$ from the low root to the high root. On the plateau $BC$, $f(v)\simeq w$ on the right attracting branch while $w$ rises slowly. The back $CD$ decreases $v$ at approximately fixed $w=w_0$ to the left attracting branch. The recovery tail then lowers $w$ back to zero on that branch, returning to $A=(0,0)$. Both sharp transitions have thickness $O(ϵ)$ in $z$; the other segments have the slower recovery scale.

The figure labels all four points and the front/back orientations. This frozen-recovery construction needs $ϵ\gamma \ll 1$ as well as $ϵ\ll 1$ if $\gamma$ is large, so that the recovery variable remains nearly fixed in an inner transition. It also needs an admissible positive front speed for the chosen orientation, addressed below.

Put $X=z/ϵ$. The travelling equations become $c{v}_X=f(v)-w+v_{XX}$ and $c{w}_X=ϵ(\gamma v-w)$. Hence, for nonzero $O(1)$ speed and $ϵ\gamma \ll 1$, $w$ is constant to leading order across the front.

Writing $p=v_X$ gives the planar system $v_X=p$, $p_X=cp-f(v)+w$. At the outer cubic roots $A=(V_{-},0)$ and $B=(V_{+},0)$, $f^{\prime }(V_{\pm })<0$. The Jacobian determinant is $f^{\prime }(V_{\pm })<0$, so both are saddle points. The increasing front is a heteroclinic orbit from the unstable manifold of $A$ to the stable manifold of $B$. Varying $c$ supplies the shooting parameter required to match these manifolds. For this cubic one can also construct it explicitly, as in the next solution.

Multiply the second-order equation by $v_X$ and integrate across the front. Because $v_X\to 0$ at both ends, the term $\int v_{XX}{v}_{X}dX=[v_X^2/2{]}_{-\infty }^{\infty }$ vanishes. The remaining identity isThe denominator is positive for a nonconstant front, so its direction is determined by the signed area between the cubic and the horizontal recovery line: Maxwell balance for a scalar reaction-diffusion front gives zero speed when that area vanishes.

The derivative in the denominator must use the inner coordinate $X$. If one instead uses the source's original $z$ and denotes $J_z=\int v_z^{2}dz$, then $J_z={ϵ}^{-1}\int v_X^{2}dX$ and the identity is $c=\int (f-w)dv/(ϵJ_z)$. Thus the printed $v^{\prime }(z)$ is valid only if $z$ has silently been reused for the rescaled inner coordinate; with the original coordinate it misses a factor $ϵ$.

Let $c_{\uparrow }(w)$ denote the increasing-front speed from the preceding solution. Reversing that heteroclinic reverses its speed, so a decreasing back has speed $c_{\downarrow }(w)=-c_{\uparrow }(w)$. A rigidly travelling pulse has one speed throughout, requiringThis is the correct use of the integral formula: compare the oppositely oriented transitions, including their signs, rather than equating two increasing-front speeds.

For completeness, write $f(v)-w=-(v-r_1)(v-r_2)(v-r_3)$, $r_1<r_2<r_3$. The explicit increasing frozen-recovery cubic front isSubstitution, using $v_X=(v-r_1)(r_3-v)/\sqrt{2}$, verifies the equation. Since $r_1+r_2+r_3=1+a$, speed matching givesFor $0<a<1/2$ this selects the appropriate positive back level within the cubic's three-root range. The high branch connects $B$ to $C$ slowly, and the low branch supplies recovery after $D$.

The stated positive-speed, rising-front orientation imposes $a<1/2$ at leading order. The printed range $0<a<1$ by itself does not guarantee that pulse construction: at $a=1/2$ the leading front is stationary, while for $a>1/2$ the calculated increasing-front speed is negative. Existence and stability of a full finite-$ϵ$ pulse also require the slow connections; speed matching is its leading necessary construction, rather than a proof for every parameter in the printed range.

The common vector field $F$ represents an identical reference cell. The small cell-specific term $ϵG_i$ represents heterogeneity of cellular properties, such as differing ion-current parameters or intrinsic oscillation rates. It is not external random noise: the displayed model is deterministic.

The Kronecker delta factors restrict coupling to the first state component, the membrane potential. Substitution gives a coupling current in the voltage equation proportional to $ϵ{\sum }_j{d}_{ij}(u_j^1-u_i^1)$, while the other state components have no direct intercellular coupling. This is the linear voltage-difference response of gap junctions. Gating and chemical variables can still change indirectly through their own coupled intracellular dynamics. The zero diagonal removes self-coupling, which would contribute zero anyway. Physical reciprocal gap-junction coupling is often symmetric, but zero diagonal alone does not impose reciprocity; that distinction matters in the final part.

With $\sigma =\omega (ϵ)t$, $\tau =ϵt$, the multiscale chain rule gives $d/dt=(1+ϵ{\Omega }_1+\cdots ){\partial }_{\sigma }+ϵ{\partial }_{\tau }$. Expanding $F(u_{i,0}+ϵu_{i,1})=F(u_{i,0})+ϵDF(u_{i,0})u_{i,1}+\cdots$ gives at order oneAt order $ϵ$, moving the linearized reference dynamics to the left givesThe component sum runs over the $m$ state variables; the source uses $n$ in its first displayed sum and $m$ later. All functions for the solvability calculation are periodic in $\sigma$ with period one.

Differentiating the autonomous cycle equation $U^{\prime }=F(U)$ gives $U^{\prime \prime }=DF(U)U^{\prime }$, so $U^{\prime }$ lies in the periodic null space of $\hat{L}$. For a hyperbolic attracting limit cycle, its tangent direction has Floquet multiplier one and every transverse multiplier has modulus below one. Hence it is the only periodic direction, proving the required one-dimensional null space.

Strictly, mere nonlinear stability does not ensure this spectral assertion. In polar coordinates $\dot{r}=-(r-1{)}^3$, $\dot{\phi }=2\pi$ has an attracting unit-period cycle but neutral radial linearization at $r=1$, supplying a second periodic solution of the variational equation. Thus the intended phase reduction assumes the usual simple unit Floquet multiplier, or equivalently hyperbolic transverse attraction; this assumption is made explicit here.

Use the period-one inner product $⟨Y,V⟩={\int }_0^{1}Y(\sigma )\cdot V(\sigma )d\sigma$. Periodic integration by parts gives ${\hat{L}}^{\ast }Y=-Y^{\prime }-DF(U{)}^{\mathsf{T}}Y$. Let $Y$ span this adjoint null space and normalize ${\int }_0^{1}Y\cdot U^{\prime }d\sigma =1$. For oscillator $i$, write $s=\sigma +{\theta }_i(\tau )$. Then ${\partial }_{\tau }{u}_{i,0}={\theta }_i^{\prime }{U}^{\prime }(s)$ and ${\partial }_{\sigma }{u}_{i,0}=U^{\prime }(s)$.

The Fredholm solvability condition makes the forcing in the order-$ϵ$ equation orthogonal to $Y(s)$. Its intrinsic and frequency terms integrate to ${\xi }_i-{\Omega }_1-{\theta }_i^{\prime }$, whereVoltage-only coupling contributesThus withthe complete phase equation isPeriodicity makes these formulae independent of the absolute initial phase. The adjoint function is a phase response function; it projects a small perturbation of each state component onto its accumulated phase change. The subtraction of $h(0)$ ensures zero coupling between identical phases, as required by the original voltage differences. The source's reference to its equation (1) is read as the uncoupled cycle equation, numbered (3) on the actual page.

At a locked state choose the rotating frame with ${\theta }_i^{\prime }=0$. For small phase differences, Taylor expansion gives $h({\theta }_j-{\theta }_i)-h(0)=h^{\prime }(0)({\theta }_j-{\theta }_i)+O(|{\theta }_j-{\theta }_i{|}^2)$. Let the graph Laplacian have $L_{ii}={\sum }_j{d}_{ij}$ and $L_{ij}=-d_{ij}$ for $i\ne j$. The leading locked equations areIf coupling is balanced, ${\sum }_j{d}_{ij}={\sum }_j{d}_{ji}$ for every $i$, so $1^{\mathsf{T}}L=0$. Summing over cells cancels all linear coupling terms and givesIn particular symmetric reciprocal gap-junction conductances satisfy this condition. To leading order the collective physical frequency is $1+ϵN^{-1}{\sum }_i{\xi }_i$.

The arithmetic-mean conclusion requires this extra balance assumption. The printed condition $d_{ii}=0$ alone is insufficient. For a connected directed coupling network let ${\pi }^{\mathsf{T}}L=0$, normalized by ${\pi }^{\mathsf{T}}1=1$. The valid weighted collective oscillator frequency result isFor two cells with $d_{12}=p>0$, $d_{21}=q>0$, one has $\pi =(q,p)/(p+q)$ and ${\Omega }_1=(q{\xi }_1+p{\xi }_2)/(p+q)$ at linear order. For example $p=1,q=2,{\xi }_1=0,{\xi }_2=\delta$ gives ${\Omega }_1=\delta /3$, rather than $\delta /2$, while the locked phase difference is $O(\delta )$ when $h^{\prime }(0)\ne 0$. This supplies a counterexample with arbitrarily small phase differences. It also identifies exactly why reciprocal coupling recovers the requested mean: the left null space is then spanned by the uniform vector.