Past exam of the mathematics and the environment and mathematical physiology course of the University of Oxford / 2007 (C8.1, DMAC 2481)

Original solutions to all eight questions and all eleven terminal parts. Source: released Mathematics papers, Part C 2007 11.pdf, examination of 28 May 2007. Modern examination statements are omitted.

In the non-scattering grey atmosphere approximation, $I$ is the frequency-integrated specific intensity along a ray, $s$ is distance along that ray, $\rho$ is the mass density, and $\kappa$ is the mass absorption coefficient. The local thermodynamic equilibrium source is the frequency-integrated Planck function $B={\sigma }_{\mathrm{SB}}{T}^4/\pi$. Across a thin element $ds$, absorption removes $\kappa \rho Ids$ and thermal emission adds $\kappa \rho Bds$. Their balance gives $dI/ds=-\kappa \rho (I-B)$.

Take the upward vertical as the polar axis, put $\mu =\cos \theta$, and define downward optical depth byAlong a ray, $dz/ds=\mu$ and $d\tau /ds=-\kappa \rho \mu$. The chain rule therefore gives $\boxed{\mu I_{\tau }=I-B}$. In a plane-parallel, azimuthally symmetric field, the solid angle element is $d\omega =d\phi d\mu$, soLocal radiative equilibrium requires zero net radiative heating of the material: absorption $4\pi \kappa \rho J$ balances emission $4\pi \kappa \rho B$. Where the absorption coefficient is nonzero this gives $B=J$, and substitution gives the required angular radiative transfer equation. It does not require the specific intensity itself to be isotropic.

At the top of an atmosphere with no external incoming long-wave radiation, impose $\boxed{I(0,\mu )=0(-1\le \mu <0)}$. These are inward directions; prescribing the same zero value for outward directions would remove the emergent radiation as well. The upward radiative flux and its derivative areThus $\boxed{\Phi \text{ is constant with height}}$: local radiative equilibrium expresses energy conservation through successive atmospheric layers.

The first Saint-Venant equation conserves water volume. The second balances material acceleration against downslope gravity, bed friction and the free-surface pressure gradient. With length scale $L$, depth scale $H$, slope $S$ and velocity scale $U$, the conventional gravity-friction scaling has $\epsilon =H/(SL)$ and Froude number $F=U/\sqrt{gH}$. The pressure-gradient/downslope-gravity ratio is $\epsilon$, and the inertia/downslope-gravity ratio is $\epsilon F^2$.

The uniform steady solution obeys $u^2=h$ and $hu=Q$, hence $\boxed{u_0=Q^{1/3},h_0=Q^{2/3}}$. In the constant-flux river momentum relaxation reduction, substitute $h=Q/u$ and $X=1-s$. Since $u_s=-u_X$ and $h_s=Q{u}_X/u^2$, the resulting quasilinear partial differential equation isA downstream boundary at $X=0$ supplies data into $X>0$ precisely when $c(u_1)>0$, equivalent to $\boxed{u_1<Q^{1/3}/F^{2/3}}$. Although water flows out there, the relevant characteristic curves carry information upstream in this subcritical regime.

Along a characteristic curve, $dX/dt=c(u)$ and $du/dt=(1-u^3/Q)/(\epsilon F^2)$. Positive $u$ relaxes toward $u_0$; near $u_0$, its displacement decays at rate $3/(\epsilon F^2{u}_0)$. Boundary-generated deviations therefore become small after a time of order $\epsilon$ and, while the characteristic speed stays bounded away from zero, within a distance of order $\epsilon$. The outer state is the initial uniform solution. A jump between initial and boundary values at their corner can produce a weak solution or a thin initial adjustment; a globally smooth corner solution requires compatible data.

For the moving boundary put $Y=X-vt$. The characteristic speed relative to it is $c(u)-v$, while the outer speed isFor the intended tidal slowing $0<u_1\le u_0$, $c$ is decreasing in $u$ and hence $c(u)\ge v_b$ throughout the relaxation. Thus $v<v_b$ leaves an incoming family and an attached thin boundary layer. The diagram shows numerical characteristics of this scalar reduction, with a slower boundary velocity. When $v>v_b$, the tidal boundary overtakes the unperturbed characteristic family; an attached smooth relaxation layer is no longer appropriate, and a tidal bore or shock wave can arise in a fuller river model.

There are two qualifications to the printed reduction. First, without the physical slowing assumption, $v<v_b$ alone need not make the boundary incoming: for $u_1>u_0$, also require $v<c(u_1)$. Second, imposing $hu=Q$ exactly everywhere in the original continuity equation gives $h_t=0$, and consequently $u_t=0$. The time-dependent characteristics above solve the requested momentum-only approximation, rather than the complete Saint-Venant equation system with an exactly constant flux.

Bedload transport carries sediment by rolling, sliding and short hops near the bed. Suspended sediment is held in the water column by turbulent mixing. The competition between particle settling velocity and the turbulent vertical velocity scale chiefly determines which description applies; particle size and density influence that settling velocity.

Scale the depths and bed/free-surface elevations by $d$, both velocities by $U$, distance by $L$, and time by $T$. A choice producing exactly the requested coefficients isThe sediment continuity equation then has unit coefficients, while both material-acceleration time derivatives acquire the factor $L/(UT)=a/d$. Dividing water momentum by $gd/L$ and sediment momentum by $U^2/L$ givesFor example $d=1\mathrm{m}$, $a=10^{-3}\mathrm{m}$, $S=10^{-3}$, $f=10^{-2}$ and ${\rho }_s/{\rho }_w=2.65$ give $\epsilon =10^{-3}$ and $\delta =2.65\times 10^{-4}$. These are illustrative values: small effective mobile-layer thickness and small slope explain the two small parameters.

Set $\epsilon =\delta =0$ at leading order. Water continuity gives $hu=q$, with spatially constant flux; take a prescribed constant $q$. Water momentum integrates to $\eta +F^2{u}^2/2=C$. Since $s=\eta -h$,The bedload transport model with sediment inertia has remaining equations $s(u{)}_t+v_x=0$ and $v{v}_x=(u-v{)}^2-\kappa v$.

Choose a uniform state of this reduced model with $h_0=u_0=q=1$, constant $s_0=C-F^2/2-1$, andIt satisfies $(1-v_0{)}^2=\kappa v_0$. For small perturbations $U,V$ of the two velocities, put $a_1=2(1-v_0)>0$ and $b_1=a_1+\kappa >0$. Linearization givesFor a Fourier mode $e^{ikx+\lambda t}$ this yieldsEvery nonzero mode decays for $F<1$ and grows for $F>1$. The spatially uniform mode is neutral, and $F=1$ is a singular threshold where the neglected terms matter. In a general uniform state the controlling sign is $q/u_0^2-F^2{u}_0$: equivalently the actual mean-flow Froude number is less or greater than one. The displayed $F=1$ criterion uses the stated unit-depth, unit-velocity normalization of the reduced steady state.

Basal sliding allows ice to pass over bed obstacles by enhanced deformation and regelation: pressure-induced melting on the upstream side and refreezing downstream transfer ice around small obstacles. Cavitation opens water-filled gaps on the downstream side. Their size and the remaining ice-bed contact depend on effective pressure in subglacial drainage $N=p_{\mathrm{ice}}-p_{\mathrm{water}}$ as well as basal shear stress $\tau$. At fixed $\tau$, reducing $N$ generally enlarges cavities and reduces resistance, giving a sliding law $u(\tau ,N)$ with $u_N<0$ in this regime.

A linked cavity drainage system connects these cavities through passages along the bed. Under the prescribed decreasing relation $N(Q)$, the chain rule gives $\boxed{u_Q=u_N{N}^{\prime }(Q)>0}$. This conclusion uses that hydraulic relation; it need not apply to a different drainage regime.

With $t=\epsilon \tau$, the ice continuity equation becomes $H_{\tau }+\epsilon (uH{)}_x=\epsilon s^{\prime }(x)$. For bounded spatial gradients over $\tau =O(1)$, $H$ is independent of the fast seasonal time to leading order. Meanwhile $S^{\prime }(Q)=\epsilon \kappa$ gives $S_t=\kappa Q_{\tau }$, soLet $J(\tau )={\int }_0^{\tau }M(r)dr$. Direct differentiation proves thatsatisfies the equation and $Q(0,\tau )=0$. This is the boundary-generated solution, or a periodic-forcing solution with the forcing continued into its past. An independently prescribed initial profile would also need matching in the region whose characteristic curves have not yet reached the upstream boundary.

If $M>0$, $J$ is invertible and rearranging gives $x=[\tau -J^{-1}(J(\tau )-\kappa Q)]/\kappa$. For $M=1-\delta \cos \tau$ with $|\delta |<1$,Substitution, followed by the sine-difference identity, givesFor fixed $Q$ and nonzero first-order amplitude, the extrema occur at $\tau =\kappa Q/2+n\pi +O(\delta )$. Since $x=Q+O(\delta )$, the seasonal phase lag in linked cavity drainage gives the stated phase relation $\boxed{\tau \approx \kappa x/2(\mathrm{mod}\pi )}$. Maxima and minima alternate; when $\sin (\kappa Q/2)=0$, the first-order approximation is flat and does not locate individual extrema.

For an experiment beginning by mixing unbound enzyme with substrate, take $\boxed{[E](0)=e_0>0,[S](0)=s_0>0,[X](0)=[Y](0)=[E_i](0)=[P](0)=0}$. The enzyme-substrate complexes, inactivated enzyme and product are initially absent. Nonzero pre-existing product could be allowed without changing the dynamics here, since it has no feedback reaction in the specified scheme.

Write concentrations as $E,S,X,Y,E_i,P$ in this calculation. The law of mass action assigns the forward bimolecular rate $k_{1}ES$ and the first-order rates $k_{-1}X,k_{2}X,k_{3}Y,k_{4}Y$. Each species gains the rates of reactions producing it and loses those consuming it:Adding the equations for all enzyme-bearing species cancels every reaction rate, so $\boxed{E+X+Y+E_i=e_0}$. This is conservation of the total enzyme-derived material, including the inactivated form. It does not assert conservation of active enzyme: mechanism-based inhibition transfers some of that conserved material irreversibly into $E_i$.

The conserved total is $e_0=[E](0)$ under the chosen initial conditions. Substitute $E=e_0-X-Y-E_i$ into the equations from the previous part. This closes the subsystem for $S,X,Y,E_i$; product can be recovered afterward by integrating $\dot{P}=k_{3}Y$.

Let $k_{\Sigma }=k_{-1}+k_2$ and $K_m=k_{\Sigma }/k_1$. The required nondimensionalization is completed byIndeed $X=e_0\sigma x/(1+\sigma )$ and the free-enzyme fraction is $1-\sigma x/(1+\sigma )-y-e_i$. Substitution and division by the indicated time scale giveFor example $k_2/k_{\Sigma }=1/(1+\rho )$ explains the coefficient in the third equation. The small multiplier on the intermediate-species derivatives identifies their potentially fast adjustment relative to substrate depletion.

For quasi-steady enzyme kinetics at low substrate saturation, the quasi-steady-state approximation is controlled here by $\epsilon =e_0/(s_0+K_m)$, rather than by $e_0/s_0$ alone. Even with $e_0\approx s_0$, taking $K_m\gg s_0$ makes $\epsilon \ll 1$ and $\sigma =s_0/K_m=O(\epsilon )$.

After the initial fast layer $T=O(\epsilon )$, the third equation givesSince $\varphi =O(\epsilon )$, the last equation gives $(e_i{)}_T=(\varphi /\epsilon )y=O(\epsilon )$. Starting from $e_i(0)=0$, this implies $e_i=O(\epsilon )$ on bounded intervals of the slow time. The quasi-steady second equation then gives $x=s+O(\epsilon )$. In the first equation,Matching to the initial substrate value $s(0)=1$ therefore yields $\boxed{s(T)\sim e^{-T/(1+\rho )}}$. The relation $x\sim s$ is an outer approximation: it does not replace the initial layer needed to reconcile $x(0)=0$ with $s(0)=1$. Likewise the smallness of accumulated inactivated enzyme has been justified for bounded slow times, not arbitrarily long times.

For the travelling wave coordinate $\xi =x+st$, the chain rule givesFast-slow matching for a calcium pulse joins two fast transitions, $B\to C$ and $D\to A$, joined by slow portions $C\to D$ and $A\to B$. Write the nullcline as $V=g(U)$. The resting point has $U_B=\mu$, $V_B=g(\mu )$ in the leading matched description.

Put $\xi =\epsilon X$. Multiplying the first equation by $\epsilon$ gives $U_{XX}=s(U_X+\gamma V_X)+\epsilon (U-\mu )$. At leading order, integrate and use $U_X\to 0$, $(U,V)\to (U_B,V_B)$ at the start of the front. Together with the second equation this givesThe fast equilibria are intersections of $V=g(U)$ and $V=V_B-(U-U_B)/\gamma$. As the slope $-1/\gamma$ changes, there may be one or three simple intersections, with two distinct intersections at a tangency. At a simple point of the nullcline, $f_U=-f_V{g}^{\prime }$, and the sign convention above/below it gives $f_V<0$. The Jacobian matrix of the fast system hasFor the usual three-intersection configuration, the two outer points have $1+\gamma g^{\prime }>0$, hence negative determinant and one positive and one negative eigenvalue: they are saddle points. The middle point has positive determinant; it is linearly stable if $s^2<-f_V$ and unstable if $s^2>-f_V$. Zero trace is a possible oscillatory stability threshold, with a nonlinear classification requiring further information about $f$.

To select the front speed, follow the appropriate branch of the unstable manifold of $B$ and vary $s$ until it meets the stable manifold of $C$. This is a shooting method for a heteroclinic orbit. Follow the high-$U$ slow branch from $C$ to $D$; the falling front obeys the same fast equations with the integration constant anchored at $D$. Choose $D$ so that its descending saddle connection reaches a low-$U$ point $A$ at that same speed. The remaining slow branch returns toward the resting point, completing the matched spike. The construction requires the indicated connections; a nullcline sketch alone does not specify their numerical speed or prove their existence for every possible kinetic function $f$.

On the slow portions, $f(U,V)=0$ at leading order, so $V=g(U)$ and $V^{\prime }=g^{\prime }(U)U^{\prime }$. Neglecting the small diffusion term yieldsOn the outer rising branches the coefficient on the left is positive. Recovery from $C$ toward $D$ requires decreasing $U$, hence $\mu <U_D$; the low branch from $A$ back toward rest requires increasing $U$, hence $U_A<\mu$. Therefore $\boxed{U_A<\mu <U_D}$. These inequalities distinguish the recovery branches from a trajectory that would move away from the required matching points.

Systole is ventricular contraction and ejection; diastole includes ventricular relaxation and filling. Let $V_a,V_v,V_{LV}$ denote arterial, venous and left-ventricular blood volumes. Let $Q^{+}$ flow through the aortic valve into the arteries, $Q_c$ through the peripheral circulation, and $Q^{-}$ through the mitral valve into the ventricle. Conservation of volume givesThe same $Q_c$ leaves arteries and enters veins because the capillary compartment has no changing stored volume in this approximation. Adding the equations conserves total modeled blood volume.

Use constant vascular compliances $C_a,C_v$, with ${\dot{V}}_a=C_a{\dot{p}}_a$ and ${\dot{V}}_v=C_v{\dot{p}}_v$, and ventricular elastance $p_{LV}=E_{LV}(t)(V_{LV}-V_0)$. The hydraulic resistances and ideal one-way valves givewhere $[z{]}^{+}=\max (z,0)$. Substitution into the conservation equations gives the two stated arterial/venous pressure equations. Since $V_{LV}-V_0=p_{LV}/E_{LV}$, differentiating this identity gives the third equation. The positive parts encode valve closure against a reverse pressure difference.

During isovolumetric contraction both valves are closed. Let $V^{+}$ be the end-diastolic ventricular volume, $V_0$ its unstressed reference volume, and $p_a^{+}$ arterial pressure before ejection. The rapidly raised elastance generates $E_s(V^{+}-V_0)$, so the aortic valve opens if $\boxed{E_s(V^{+}-V_0)>p_a^{+}}$.

During ejection, $R_a/R_c\ll 1$ makes the open-valve pressure drop small on the peripheral drainage time scale, so $p_{LV}\approx p_a$. The mitral valve is closed, and adding arterial and ventricular volume balances eliminates $Q^{+}$. With $p_v\ll p_a$ this givesFor approximately constant systolic elastance, set $C_s=1/E_s$ and reset time to the start of the slow ejection phase. ThenHere ${\hat{p}}_a$ is the matched pressure at that start. If peripheral outflow is negligible during the preceding rapid artery-ventricle pressure equilibration, conservation gives the more explicit ${\hat{p}}_a=[C_a{p}_a^{+}+V^{+}-V_0]/(C_a+C_s)$.

After aortic closure the ventricular volume is fixed at the end-systolic value $V^{-}$. Relaxation decreases $E_{LV}$ until the mitral valve opens at $p_{LV}=p_v$, hence $\boxed{E_{LV}\sim p_v/(V^{-}-V_0)}$, with $p_v$ the venous filling pressure at opening. During filling the aortic valve is closed, so $R_c{C}_a{\dot{p}}_a=-(p_a-p_v)\approx -p_a$. With $p_a^{-}$ the arterial pressure at the start of filling and time reset there,The change from $C_a+C_s$ to $C_a$ is physical: the ejecting ventricle shares the arterial pressure during systole, whereas a closed aortic valve disconnects it during filling.

In this delayed carbon-dioxide ventilation feedback model, $p$ is a representative carbon-dioxide partial pressure, $K$ an effective storage capacitance, $M>0$ metabolic production, and $pV$ removal by pulmonary ventilation. The delay represents transport and sensing between the gas-exchange compartment and the ventilatory response. The Hill equation response rises sigmoidally with pressure, reaches half its maximum $\alpha$ at $p=\theta$, and saturates at $\alpha$; $n>1$ permits a steep threshold-like response.

For positive parameters, $pV(p)$ increases strictly from zero to infinity because its derivative is $V+p{V}^{\prime }>0$. Its graph therefore intersects the horizontal line $M$ exactly once: $\boxed{M=p^{\ast }V(p^{\ast }),p^{\ast }>0}$. The figure shows this argument and the first imaginary-root crossing.

For a small perturbation $z=p-p^{\ast }$, linearization givesThe sigmoidal delayed ventilation threshold follows from the characteristic equation of a delay differential equation, $\lambda +A+B{e}^{-\lambda \tau }=0$. At an imaginary root $\lambda =i\omega$, separate real and imaginary parts:For $B>A>0$, the first positive crossing hasAt zero delay the scalar equation is stable. No root can cross the imaginary axis before this first crossing, and zero itself is not a root. To check the crossing direction, implicit differentiation givesAll subsequent imaginary crossings have the same direction. Thus the equilibrium is linearly stable below ${\tau }_c$ and unstable above it. This establishes the oscillatory linear threshold; the full nonlinear Hopf bifurcation additionally requires the usual nondegeneracy conditions. If $B\le A$, no finite imaginary-root crossing occurs and the equilibrium is stable for every finite delay. For the stated response, $B/A=n{\theta }^n/({\theta }^n+(p^{\ast }{)}^n)$, showing explicitly when the gain can exceed instantaneous damping.

If $A\ll B$, then ${\omega }_c\approx B$ and ${\omega }_c{\tau }_c\approx \pi /2$. Consequently the small-amplitude period near onset is $\boxed{2\pi /{\omega }_c\approx 4{\tau }_c}$. The comparison is to the critical delay near the first instability, rather than an arbitrary delay far beyond onset. A longer circulatory transit time can move the physiological system across this delayed-feedback threshold, providing a mechanism for periodic breathing associated with impaired cardiac output.