Past exam of the mathematical physiology course of the University of Oxford / 2023 (C5.12, A12214W1)

Solutions to all three questions of the mathematical physiology course of the University of Oxford examination held on 31 May 2023. Original paper: public examination PDF, inclusive bundle pages 321–324; paper pages 1–4. All eighteen lettered parts are covered, including the third question even though only the best two answers counted. Question 1(c) contains a printed exponent error; Question 2(f) requires a signed recovery cutoff. Both are explained where they arise.

Write $a=[\mathrm{A}]$, $b=[\mathrm{B}]$ and $c=[\mathrm{C}]$ for the three concentrations. The letters denote chemical species; the forward arrow represents association of one molecule of each reactant into one complex, and the reverse arrow represents dissociation. The positive constants $k_{+}$ and $k_{-}$ are the respective reaction rate constants. For an elementary chemical reaction, the law of mass action gives forward and reverse reaction ratesEvery forward event removes one molecule of each reactant and adds one complex; a reverse event reverses those changes. Thus the reaction-rate equations areThe units of $k_{+}$ are inverse concentration per unit time, whereas $k_{-}$ has units of inverse time. Adding the equations gives the independent conservation lawsIn particular $a-b$ is constant, but this is a consequence of those two conservation laws. The total $a+b+c$ need not be constant: association changes the number of separate molecules.

For the general stoichiometry, take $p,q$ to be positive integers. The net reaction rate and stoichiometric vector areConsequentlyThe factors $p,q$ outside $R$ count molecules consumed per event; the powers inside $R$ express the law of mass action. The corresponding conservation laws are $a+pc=\text{constant}$ and $b+qc=\text{constant}$. No factorial is required in this deterministic convention for the reaction rate constant.

There are four reversible transitions of the ion channel, each changing exactly one gate:These are first-order reaction-rate equations: their reaction rate constants are per-gate transition rates. In the continuous-time Markov chain interpretation, simultaneous changes of both gates have probability of order $(dt{)}^2$, so there is no direct diagonal transition in its infinitesimal generator matrix.

Normalize the total channel density to one. For independent permissive gates in an ion channel, the probability that the $M$ gate is open is $m$, and the probability that the $H$ gate is open is $h$. The joint probability is the product of the appropriate open or closed probabilities, givingIn particularThe original PDF really prints $j$, rather than $1-j$, on the final factor. That expression fails normalization of a probability distribution and does not give the closed-$H$ probability at $j=0$; the displayed correction is necessary. If $S_{ij}$ are literal densities with conserved total density $N$, multiply each expression by $N$.

The product additionally assumes initially independent events and independent gate transitions. The transition scheme alone permits initially correlated gates and does not force this factorization at every time. The all-closed initial preparation used below satisfies the required product assumption.

For each state, add its incoming reaction rates and subtract its outgoing reaction rates. This gives the master equationSumming gives ${\sum }_{i,j}{\dot{S}}_{ij}=0$, verifying the conservation law for total channel density. Nonnegative initial densities stay nonnegative because every loss term vanishes when its own density reaches zero.

With total density one, the open probabilities for the individual gates are $m=S_{10}+S_{11}$ and $h=S_{01}+S_{11}$. Adding the relevant master equations cancels transitions of the other gate:These marginal ordinary differential equations hold even for correlated initial gates. Their product gives $S_{11}$ when the independence assumption in part (c) holds. Solving the first-order linear differential equations from the all-closed initial state givesPut $\lambda =\alpha +\beta =\gamma +\delta >0$. ThenFor positive opening rates this increases strictly from zero, so its half-equilibrium time is uniquely determined by $(1-e^{-\lambda t}{)}^2=1/2$:If an opening rate is zero, the all-open equilibrium density is zero and a nontrivial half-rise time is undefined.

The two identical $M$ gates are independent Bernoulli random variables, so their open count has a binomial distribution. Multiplying by the independent $H$ factor gives the normalized channel densitiesExplicitly,The binomial coefficient $\left(\genfrac{}{}{0ex}{}{2}{1}\right)=2$ counts which of the two distinguishable $M$ gates is open. At open count $i$, the total $M$-opening rate is $(2-i)\alpha$ and the total closing rate is $i\beta$. Dividing the expected value of that count by two therefore recovers the same single-gate ordinary differential equation:For the specified initial conditions and common rate sum $\lambda$, the all-open fraction isThus, assuming positive opening rates,Three gates take longer to reach half their final all-open fraction than two gates: all three activation factors must build up simultaneously.

In the Hodgkin-Huxley model, $V$ is the intracellular-minus-extracellular membrane potential, $C_m$ is the membrane capacitance per unit area, and $I_{\mathrm{app}}$ is the applied inward electric current per unit area. A space clamp makes $V$ spatially uniform, leaving temporal ordinary differential equations. With outward ionic electric current positive, the capacitance balance is $C_m\dot{V}=I_{\mathrm{app}}-I_i$.

The variables $m,n,h$ lie between zero and one. They describe sodium activation, potassium activation and availability against sodium inactivation, respectively. Thus larger $h$ means less inactivation. Their voltage-dependent limiting values are $m_{\infty },n_{\infty },h_{\infty }$, and ${\tau }_m,{\tau }_n,{\tau }_h$ are their positive relaxation time constants.

Each ionic term has the Ohm's law form conductance times driving potential. The sodium term has instantaneous electrical conductance $g_{\mathrm{Na}}{m}^{3}h$, the potassium term has instantaneous electrical conductance $g_{\mathrm{K}}{n}^4$, and the leak term has constant electrical conductance $g_{\mathrm{L}}$. Here the $g$ constants are maximal electrical conductances per area, and $V_{\mathrm{Na}},V_{\mathrm{K}},V_{\mathrm{L}}$ are the corresponding reversal potentials. The factors $m^{3}h$ and $n^4$ arise from independent permissive gates in an ion channel. A conducting channel carries zero net ionic electric current at its reversal potential.

Use the three standard reduction assumptions: $m$ relaxes much faster than $V,n,h$ and can be replaced by $m_{\infty }(V)$; the two slower gates have the same time constant, ${\tau }_n(V)={\tau }_h(V)=\tau (V)$; and their limiting open fractions have a voltage-independent sum, $n_{\infty }(V)+h_{\infty }(V)=K$. The constant-sum assumption is essential: equality of two time constants by itself does not preserve a moving voltage-dependent sum.

For $q=n+h$, addition of the two gate ordinary differential equations givesThe specified initial condition is $q(0)=K$. Uniqueness for this first-order linear differential equation gives $q(t)=K$, hence $h=K-n$ for all time. The remaining two ordinary differential equations areThese assumptions give a planar reduction of the Hodgkin-Huxley model. The FitzHugh-Nagumo model keeps its fast activation and slow recovery variable geometry, approximating the voltage response by a cubic polynomial and the recovery response by a linear function after nondimensionalization. The polynomial equations in part (c) are a further model simplification; they are not an exact algebraic consequence of arbitrary original gate functions.

Assume $\epsilon >0$. At an equilibrium, $w=\gamma v=f(v)$, soThe nonzero roots areConsequently, with ${\gamma }_c=(1-a{)}^2/4$, there are three distinct equilibria for $0<\gamma <{\gamma }_c$, two distinct equilibria for $\gamma ={\gamma }_c$, and only the origin for $\gamma >{\gamma }_c$. At equality the two nonzero equilibria coincide; for $0<\gamma <{\gamma }_c$ they both have $a<v<1$. Geometrically, ${\gamma }_c$ is the maximum of the concave quadratic $(a-v)(v-1)$.

Since $f^{\prime }(0)=-a$, the Jacobian matrix at the origin isThus the origin is always asymptotically stable. Its eigenvalues areWriting ${\gamma }_N=(a-\epsilon {)}^2/(4\epsilon )$, the origin is a stable node for $0<\gamma <{\gamma }_N$ and a stable focus for $\gamma >{\gamma }_N$. At $\gamma ={\gamma }_N>0$ it is a degenerate node: the repeated negative eigenvalue has only one independent eigenvector, since the nonzero off-diagonal entry prevents $J_0$ being a scalar matrix. If $a=\epsilon$, then ${\gamma }_N=0$ and every permitted $\gamma >0$ gives a stable focus. This change from node to focus does not change stability and is not a Hopf bifurcation.

Take $\gamma >{\gamma }_c$ and let $\epsilon \to 0^{+}$ with the other parameters fixed. The nullclines are the cubic polynomial graph $w=f(v)$ and the line $w=\gamma v$. In the phase plane, $v$ increases below the cubic polynomial graph and decreases above it; $w$ increases below the line and decreases above it. The fast ordinary differential equation, in time $s=t/\epsilon$, freezes $w$ and gives $dv/ds=f(v)-w$.

At $w=0$, the fast roots are $0,a,1$. Their linear stability follows from $f^{\prime }(0)=-a<0$, $f^{\prime }(a)=a(1-a)>0$ and $f^{\prime }(1)=a-1<0$. Thus a fixed $0<v^{\ast }<a$ falls towards the resting branch, whereas a fixed $a<v^{\ast }<1$ rises towards $v=1$. If $v^{\ast }>1$, the first fast segment instead settles down towards $1$, followed by the same recovery excursion. The threshold $a$ is the leading singular perturbation threshold; finite-$\epsilon$ corrections matter for initial values approaching $a$ with $\epsilon$.

LetThese are the roots of $f^{\prime }(v)=0$. The outer branches $v<v_L$ and $v>v_R$ of the critical manifold are attracting under fast motion; the middle branch is repelling. The full excursion has four stages.

For this cubic polynomial, the landing point is explicit. Since $f(v)-f(v_R)=-(v-v_R{)}^2(v-v_J)$, comparison of the quadratic coefficients givesThus a graph of $v(t)$ has a rapid rise, a slowly falling elevated plateau, a rapid fall through zero to a negative undershoot, and a slow approach back to zero. The figure shows both the phase plane and a computed time trace. The smaller initial perturbation stays near the resting branch and produces no elevated plateau.

The four-stage description concerns the singular limit. Close to a simple fold point of a critical manifold, balancing the local quadratic activation term with the slow drift gives $\Delta v=O({\epsilon }^{1/3})$ and $\Delta t=O({\epsilon }^{2/3})$ for fixed parameters. These layers shrink compared with the slow phases and do not alter the leading time integrals below.

On an attracting branch of the critical manifold, $w=f(v)$ to leading order. The chain rule therefore givesThe initial fast segment ends at $v=1$, and the first slow phase ends at the right fold point of a critical manifold, $v_R$. HenceOn this interval $f^{\prime }(v)<0$ and $\gamma v-f(v)>0$. The integrand is negative, but the upper limit is smaller than the lower limit, so $T_1>0$. The numerator tends to zero at the fold while the denominator remains positive; the leading integral is finite. Fast segments and the fold delay give vanishing corrections as $\epsilon \to 0^{+}$.

The requested second phase is the second slow recovery phase, beginning after the downward fast jump. Its starting point is $v_3=v_J<0$, fixed by $f(v_J)=f(v_R)$, with the negative root chosen. Along it the same chain rule givesThe approach is to zero from below. Thus, if the printed cutoff $\beta$ is a signed value, it must satisfy $v_J<\beta <0$ and $v_4=\beta$. If $\beta$ denotes a small positive magnitude, the endpoint is instead $v_4=-\beta$. To express both conventions unambiguously, put $b=|\beta |>0$ and writeThere is no positive small endpoint on this leading slow branch. A positive cutoff would require a different event or finite-$\epsilon$ oscillatory recovery.

For bounded negative $v$ and $\gamma \gg 1$,Integration therefore gives the useful more detailed approximationTaking also $b\ll 1$, the leading logarithm isIndeed near zero, $f(v)=-av+O(v^2)$ and the slow equation is $\dot{v}=-(\gamma +a)v/a+O(v^2)$, so at fixed $\gamma$ the exact leading logarithmic coefficient is $a/(\gamma +a)$. Its large-$\gamma$ limit is $a/\gamma$, agreeing with the integral.

The order of approximation matters: first take the singular limit at fixed $\gamma$, then the large-$\gamma$ and small-cutoff limits. For simultaneous limits, a sufficient separation near rest is $\epsilon \gamma \ll a^2$. Otherwise the fast and recovery time scales merge, and these slow-branch formulas need not describe the late-time motion.

In this two-pool calcium model, $c$ is the cytosolic calcium ion concentration, and $c_s$ measures calcium in an intracellular store, such as the endoplasmic reticulum. The variables in these balances are expressed using compatible effective volumes, so internal transfers enter their two equations with equal and opposite signs. If actual compartment concentrations use different reference volumes, the corresponding volume ratio must first be included in their balances.

The flux $J_{+}$ pumps calcium from the cytosol into the store, $J_{-}$ releases it back, and $k_s{c}_s$ is passive store leakage back to the cytosol. The external input $r$ supplies cytosolic calcium, while $kc$ exports it from the biological cell. Thus $k,k_s$ have units of inverse time and all five fluxes have units of concentration per time. Their directions can be sketched asAdding the two ordinary differential equations verifies the total calcium balanceInternal uptake, release and leakage redistribute calcium; only the exchanges across the cell membrane change this total. In particular any equilibrium must have $c=r/k$ when $k>0$.

Choose positive reference concentrations $C,C_s$ and a reference flux $J_0$, and defineOverdots now mean derivatives with respect to the dimensionless time $\tau$. The store balance becomesThe total balance becomes $\dot{u}+(C_s/C)\dot{v}=r/(kC)-u$. ThereforeandThe reference quantities are scale choices, not additional kinetic variables; typically the response functions supply natural saturation concentrations and a maximal uptake flux. The identity $\delta /\epsilon =k_s/k$ checks the nondimensionalization. A small $\epsilon$ means store exchange is fast compared with external cytosolic removal.

At any equilibrium, the first balance gives $u=\mu$ and the second gives $f(\mu ,v)=0$. For every fixed $u>0$,andThe intermediate value theorem and strict monotonicity give precisely one positive root $v=g(u)$ for each $u>0$. Since $f_v\ne 0$, the implicit function theorem makes this a smooth nullcline. For every $\mu >0$, therefore,is the unique positive equilibrium. Graphically, draw the vertical line $u=\mu$ through the given graph: it meets the graph once, even on its steep descending section. This vertical line is the condition for simultaneous stationarity; the full $u$ nullcline at finite $\epsilon$ is actually $f(u,v)=(\epsilon /\gamma )(\mu -u)$.

Eliminate $\dot{v}$ from the first equation before linearization:All partial derivatives below are evaluated at $(\mu ,g(\mu ))$. The Jacobian matrix isIts trace and determinant areBy linear stability of a planar equilibrium, a positive trace makes both real eigenvalues positive or gives a complex pair with positive real part. The equilibrium is consequently unstable precisely whenA negative trace gives asymptotic stability; equality leaves a purely imaginary pair and requires a separate nonlinear bifurcation calculation. In particular, the unique positive equilibrium is never a saddle equilibrium.

Differentiate the identity $f(u,g(u))=0$ using the chain rule:Substitution into the trace givesBecause $f_v<0$, the exact instability condition isFor fixed $\delta >0$ and fixed equilibrium location, the correction tends to zero with $\epsilon$. To leading order this isThus instability requires the descending nullcline to be steeper than a line of slope $-1/\gamma$. The sign of the finite-$\epsilon$ correction makes the exact unstable interval slightly smaller. Since $\delta =0.001$ can make $|f_v|$ small, the approximation quantitatively requires $\epsilon /|f_v|\ll 1$ at the point in question; merely comparing $\epsilon$ with one can be insufficient.

On the descending part of the graph, mark the two tangency points whereChoose the crossings that enclose the portion with $g^{\prime }(u)<-1/\gamma$. For fixed $\mu$ strictly between these values and sufficiently small $\epsilon$, the unique equilibrium is a repelling equilibrium by part (e). The endpoints here are leading singular limit values; at finite $\epsilon$ the boundaries instead solve $f_v-\gamma f_u=\epsilon$ and approach these tangencies when the crossings are nondegenerate.

There is a necessary qualification: such an interval exists only when $\gamma$ is large enough that the descending graph actually attains slopes below $-1/\gamma$. If $M={\max }_{u>0}[-g^{\prime }(u)]$, one needs $\gamma M>1$. The supplied graph has a very steep descending section, but an arbitrarily small positive $\gamma$ need not satisfy this condition. For $\delta =0.001$, numerical evaluation gives $M\simeq 511.49$, so the leading threshold is $\gamma \simeq 0.0019551$. A graph alone does not remove this parameter restriction.

The fast motion is especially informative. In fast time $s=\tau /\epsilon$,At leading order $z=u+\gamma v$ is constant, so fast paths are straight lines with slope $dv/du=-1/\gamma$. Along such a line, the fast eigenvalue at the nullcline isIt is negative on the two outer attracting segments and positive on the intervening repelling segment. The tangencies defining ${\mu }_{\pm }$ are consequently the fold points of a critical manifold when it is expressed in the coordinates $(z,v)$. They are not generally the maximum and minimum of $g$ itself.

For ${\mu }_{-}<\mu <{\mu }_{+}$, slow motion on an attracting segment followsOn the left attracting segment $u<{\mu }_{-}$ this moves towards increasing $u$; on the right attracting segment $u>{\mu }_{+}$ it moves towards decreasing $u$. Each segment ends at a tangency and is followed by a fast jump along a line of slope $-1/\gamma$ to the other attracting segment. This constructs the limiting relaxation oscillation. The figure uses $\gamma =0.01$, for which ${\mu }_{-}\simeq 0.62970$ and ${\mu }_{+}\simeq 1.01157$; the marked endpoints change with $\gamma$.

To justify existence of a periodic orbit of the full planar system, rather than just draw a formal cycle, first establish a compact trapping region. The positive quadrant is forward invariant: at $u=0$, $\dot{u}=\mu +\gamma \delta v/\epsilon >0$, and at $v=0$, $\dot{v}=u^2/[\epsilon (1+u^2)]\ge 0$. MoreoverChoose $V>1/\delta$ above the initial store value. Since $f(u,v)\ge -1-\delta V$ when $0\le v\le V$,Thus the rectangle $0\le u\le U$, $0\le v\le V$ is a positively invariant set if $U>\mu +\gamma (1+\delta V)/\epsilon$ and the initial point lies inside it. Every forward solution is bounded. There are no equilibria on its axes for $\mu >0$, and its sole interior equilibrium is a hyperbolic equilibrium point with both eigenvalues having positive real parts in the strict unstable range. A sufficiently small neighborhood of this repelling equilibrium has outward flow; deleting it leaves a compact trapping region with no equilibria. The Poincaré-Bendixson theorem then supplies a periodic orbit.

Hence periodic solutions exist throughout the strict unstable interval, interpreted with the finite-$\epsilon$ correction above. This argument establishes existence, without asserting uniqueness, the criticality of possible endpoint Hopf bifurcations, or absence of coexisting periodic orbits outside the interval.