The Lotka-Volterra predator-prey model

Concept

The Lotka-Volterra predator-prey model

T-179Home BU-104Threads evolution

Statement

Let a prey population of density \(N(t)\) and a specialist predator population of density \(P(t)\) obey the closed, unstructured, deterministic system \(\dot N = rN - aNP\), \(\dot P = \varepsilon aNP - mP\), with strictly positive intrinsic prey growth rate \(r\), mass-action attack rate \(a\), conversion efficiency \(\varepsilon\) and predator death rate \(m\), and with \(N(0)\gt 0\), \(P(0)\gt 0\). Then the open first quadrant contains exactly one equilibrium, \(\left(N^{*},P^{*}\right)=\left(m/\varepsilon a,\ r/a\right)\); the function \(H(N,P)=\varepsilon aN - m\ln N + aP - r\ln P\) is constant along every solution; every orbit other than the equilibrium is a closed curve, so every solution is periodic and neither population ever reaches zero; the equilibrium is Lyapunov stable but not asymptotically stable, and the period of small oscillations is \(2\pi/\sqrt{rm}\) with the prey peak leading the predator peak by a quarter cycle. Over one full period the time averages are exactly the equilibrium values, \(\langle N\rangle = N^{*}\) and \(\langle P\rangle = P^{*}\), from which Volterra’s principle follows: mortality applied to both species at a common rate \(h\lt r\) raises the average prey density to \((m+h)/\varepsilon a\) and lowers the average predator density to \((r-h)/a\).

Why it matters

Two species interacting is the smallest ecological system that is not just a population, and predation is the interaction that structures most food webs. The Lotka–Volterra pair is the minimal model of it: Logistic population growth asks what one population does when it limits itself, and this asks what two populations do when each is limited by the other. The answer is qualitatively unlike anything a one-species model can produce — not an approach to a steady state but a sustained oscillation, with the predator peak always trailing the prey peak, and with an amplitude fixed by history rather than by the parameters. That is the archetype every long-term abundance record is read against, from the boreal hare and lynx cycle to rodent–mustelid cycles at high latitudes.

Its importance to evolution, the home unit of this page, is that it supplies the demographic engine underneath an antagonistic arms race. Natural selection acts through differential survival and reproduction, and in a predator–prey pair the strength of selection on prey defence is set by the predation term \(aP\), which is itself a dynamical variable. A heritable improvement in prey defence lowers \(a\), and the model answers immediately: both averages rise, since \(\langle P\rangle=r/a\) and \(\langle N\rangle=m/\varepsilon a\) — the prey gains, and so does its enemy. A heritable improvement in predator efficiency raises \(\varepsilon\), which lowers the average prey density and leaves the average predator density exactly where it was, so the predator’s own improvement does nothing at all for its own numbers. Selection on one partner therefore alters the ecological context in which the other is selected, which is the formal content of Coevolution, and the reason predator–prey coevolution is treated as a feedback rather than as a sequence of one-sided improvements. The same equations, with the interaction sign flipped, give the competition model behind Competitive exclusion, and with the roles relabelled as susceptible and infectious hosts they become the mass-action core of The SIR model and of The basic reproduction number.

The applied consequences turn almost entirely on the time-average theorem below. It says that a non-selective agent of mortality — a broad-spectrum insecticide, a trawl fishery that catches predatory and prey fish alike, a cull — must on average increase the prey and decrease the predator. That is the standard model-based explanation for secondary pest resurgence after spraying, for the rise in the predatory fraction of Adriatic fish landings when fishing was suppressed during the First World War, and for why Keystone species arguments about top predators are quantitative rather than rhetorical.

Hypotheses

The prey grows exponentially at rate \(r\) in the absence of the predator: no self-limitation, no resource ceiling.This is the hypothesis that makes the closed orbits possible, and it is the one that is always false in nature. Restoring a carrying capacity, \(\dot N = rN(1-N/K)-aNP\), leaves the equilibrium prey density untouched but converts the centre into a stable spiral (Fails without, and Problem 4): the perpetual cycles are destroyed by an arbitrarily small amount of density dependence.

Predation follows mass action: encounters occur at rate proportional to \(NP\), and consumption per predator is \(aN\), unbounded in \(N\).A predator with a finite handling time cannot eat at a rate proportional to prey density for ever. Replacing \(aN\) by the Holling type II form \(aN/(1+a\tau_h N)\) saturates consumption at \(1/\tau_h\) and, once prey self-limitation is also present, produces a Hopf bifurcation to a genuine limit cycle as the environment is enriched — the paradox of enrichment (Problem 5).

The predator is an obligate specialist that dies exponentially at rate \(m\) with no prey.A generalist predator with alternative food does not crash when this prey becomes scarce, so the starvation phase that lets the prey recover is weakened or absent, and the predator equation loses its dependence on \(N\) at low prey density. Generalist predation typically stabilises rather than cycles.

Both populations are large enough, and reproduction continuous enough, that densities are differentiable functions of time and demographic noise is negligible.The orbits are neutrally stable, so a stochastic version has no restoring force at all: amplitude performs a random walk and the trajectory eventually hits an absorbing boundary. This is why Gause could not maintain Didinium and Paramecium together in a homogeneous culture, and it is a failure of the deterministic idealisation rather than of the biology.

The system is closed, spatially homogeneous and free of time lags: no immigration, no refuges, no age structure, no maturation delay.Space is the classic rescue. Huffaker’s orange-tray experiments with predatory and herbivorous mites persisted for many cycles only once dispersal barriers made the arena patchy, so that local extinctions were recolonised from elsewhere — the mechanism formalised in Metapopulation dynamics. A discrete-generation analogue with a full generation of lag, the Nicholson–Bailey host–parasitoid model, is not neutrally stable at all: its oscillations diverge.

All four parameters are constants: no seasonality, no evolution, no dependence on density or on age.The parameters are precisely the traits selection acts on, so on a timescale where the prey evolves escape (\(a\) falls) or the predator evolves efficiency (\(\varepsilon\) rises), the conserved quantity \(H\) is no longer conserved and the orbits drift between level curves. Seasonal forcing of \(r\) converts the neutral centre into a resonantly forced oscillator, which can be chaotic.

Proof

Steps 1–4 build the equations and strip them to a single dimensionless parameter. Steps 5–7 locate the equilibria and show that linearisation cannot settle the question at the coexistence point. Steps 8–11 supply the conserved quantity, prove from its convexity that every orbit is a closed curve, and identify the structural reason — the system is Hamiltonian in logarithmic coordinates. Steps 12–14 extract the period, the phase relation, the exact time averages and Volterra’s principle.

\[ [N]=[P]=\text{individuals per unit area},\qquad [r]=[m]=\text{time}^{-1},\qquad [a]=\frac{1}{[P]\,\text{time}},\qquad [\varepsilon]=\frac{[P]}{[N]} \]

Bookkeeping first, because every later formula must be checked against it. \(N\) and \(P\) are densities of a single closed population each; \(r\) and \(m\) are per-capita rates; \(a\) is fixed by requiring \(aNP\) to be a prey loss rate, so it is an area searched per predator per unit time; \(\varepsilon\) is the number of predators produced per prey consumed and therefore carries the ratio of the two density units. When the two species are counted in the same units \(\varepsilon\) is dimensionless. A

\[ \frac{1}{N}\frac{dN}{dt}=r-aP\qquad\Longrightarrow\qquad \dot N = rN - aNP \]

The prey equation, written as a per-capita statement because that is what the biology asserts. In the absence of predators each prey individual contributes a net \(r\) offspring per unit time (Hypotheses, no self-limitation). Predation removes prey at a rate proportional to the encounter rate between the two populations; assuming random mixing and independent movement, the rate of encounters per unit area is proportional to the product of the densities, exactly as a bimolecular reaction rate is proportional to the product of reactant concentrations. The per-capita prey death rate is therefore \(aP\), linear in predator density and independent of prey density. A

\[ \frac{1}{P}\frac{dP}{dt}=\varepsilon aN-m\qquad\Longrightarrow\qquad \dot P = \varepsilon aNP - mP \]

The predator equation. Every prey killed is converted into new predators with efficiency \(\varepsilon\), so the predator birth rate is \(\varepsilon\) times the same term \(aNP\) that appears with a minus sign in the prey equation — the two equations are coupled through one shared flux, which is what makes the pair a model of predation rather than two independent populations. The predator has no other food, so its per-capita rate falls to \(-m\) as \(N\to 0\). Note that \(\varepsilon\) collapses several distinct biological quantities (assimilation, respiration, the energy cost of reproduction) into one number; that \(\varepsilon\ll 1\) in real food chains is the demographic face of Trophic levels and energy flow. A

\[ x=\frac{N}{N^{*}},\quad y=\frac{P}{P^{*}},\quad \tau=rt,\quad \mu=\frac{m}{r}\qquad\Longrightarrow\qquad \frac{dx}{d\tau}=x\left(1-y\right),\qquad \frac{dy}{d\tau}=\mu\,y\left(x-1\right) \]

Nondimensionalisation, done before any analysis so that no conclusion can depend on a choice of units. Substituting \(N=N^{*}x\) into Step 2 gives \(\dot x = rx - aP^{*}xy = rx(1-y)\), because \(aP^{*}=r\); substituting \(P=P^{*}y\) into Step 3 gives \(\dot y = y(\varepsilon aN^{*}x-m)=my(x-1)\), because \(\varepsilon aN^{*}=m\). Measuring time in units of \(1/r\) leaves exactly one parameter, \(\mu=m/r\). Four biological constants therefore contain only one piece of dynamical information: \(a\) and \(\varepsilon\) fix the scales of the axes and nothing else, so no manipulation of attack rate or conversion efficiency can change the shape of the cycle, only where it sits. B

\[ \dot N=\dot P=0\iff\left\{N\left(r-aP\right)=0,\ P\left(\varepsilon aN-m\right)=0\right\}\iff \left(N,P\right)=\left(0,0\right)\ \text{or}\ \left(N^{*},P^{*}\right)=\left(\frac{m}{\varepsilon a},\ \frac{r}{a}\right) \]

The equilibria, obtained by factorising rather than by dividing, so that the boundary solutions are not lost. Each equation is a product of two factors; setting \(N=0\) forces \(P=0\), and \(P=0\) forces \(N=0\) unless \(r=0\). The prey nullclines are the line \(N=0\) and the horizontal line \(P=r/a\); the predator nullclines are \(P=0\) and the vertical line \(N=m/\varepsilon a\), and the only interior fixed point is where the last two cross. The result is the model’s most counterintuitive feature and it is exact: neither species’ equilibrium density contains its own demographic rate. The prey level \(m/\varepsilon a\) is assembled from predator traits alone and does not contain \(r\); the predator level \(r/a\) contains neither \(m\) nor \(\varepsilon\). A more fecund prey does not equilibrate at a higher density; it supports more predators. B

\[ J\left(N,P\right)=\begin{pmatrix} r-aP & -aN\\ \varepsilon aP & \varepsilon aN-m\end{pmatrix},\qquad J(0,0)=\begin{pmatrix} r&0\\0&-m\end{pmatrix},\qquad \lambda=r\gt 0,\ \ \lambda=-m\lt 0 \]

The origin is a hyperbolic saddle. Its unstable manifold is the \(N\) axis (prey alone grow exponentially) and its stable manifold is the \(P\) axis (predators alone starve exponentially). Both axes are invariant sets, since \(N=0\) makes \(\dot N=0\) and \(P=0\) makes \(\dot P=0\): a trajectory starting with both species present can never reach zero in finite time, so the deterministic model never predicts extinction. It is worth stating why that is not reassuring — the orbit can pass arbitrarily close to the axes, and in a finite population “arbitrarily close to zero” means extinct. A

\[ J\left(N^{*},P^{*}\right)=\begin{pmatrix} 0 & -\dfrac{m}{\varepsilon}\\[4pt] \varepsilon r & 0\end{pmatrix},\qquad \operatorname{tr}J=0,\qquad \det J=rm\gt 0,\qquad \lambda_{\pm}=\pm i\sqrt{rm} \]

Substituting \(N^{*}=m/\varepsilon a\) and \(P^{*}=r/a\) kills both diagonal entries: \(r-aP^{*}=0\) and \(\varepsilon aN^{*}-m=0\). The eigenvalues are purely imaginary, so the linearised system is a centre — and this is exactly the case in which linearisation proves nothing. The Hartman–Grobman theorem requires a hyperbolic fixed point, one with no eigenvalue on the imaginary axis, and this one is not hyperbolic. A linear centre may be a nonlinear centre, a stable spiral or an unstable spiral: the system \(\dot u=-v+u^{3}\), \(\dot v=u+v^{3}\) has the same linearisation and spirals outwards. The closed orbits must therefore be established globally, not locally, and Steps 8–10 do so. C

\[ \frac{dy}{dx}=\frac{\mu y\left(x-1\right)}{x\left(1-y\right)}\ \Longrightarrow\ \frac{1-y}{y}\,dy=\mu\,\frac{x-1}{x}\,dx\ \Longrightarrow\ \ln y-y=\mu\left(x-\ln x\right)+\text{const} \]

Eliminate time. Away from the nullclines the orbit satisfies a first-order equation in \(x\) and \(y\) alone, and it separates: every term on the left involves only \(y\) and every term on the right only \(x\). Integrating \(\int (1/y-1)\,dy\) and \(\mu\int(1-1/x)\,dx\) gives the displayed relation. Rearranged, it says that \(V(x,y)=\mu(x-\ln x)+(y-\ln y)\) takes the same value at every point of an orbit. This is the whole content of the model’s integrability, and it exists because the vector field, though nonlinear, is a product of a function of \(x\) and a function of \(y\) in each component. B

\[ H\left(N,P\right)=\varepsilon aN-m\ln N+aP-r\ln P,\qquad \frac{dH}{dt}=\left(\varepsilon a-\frac{m}{N}\right)\dot N+\left(a-\frac{r}{P}\right)\dot P=\left(\varepsilon aN-m\right)\left(r-aP\right)+\left(aP-r\right)\left(\varepsilon aN-m\right)=0 \]

The same invariant in dimensional variables, verified directly rather than inherited. Substituting \(\dot N=N(r-aP)\) and \(\dot P=P(\varepsilon aN-m)\) makes the two terms exact negatives of each other, so the derivative vanishes identically — on every solution, for every parameter set, with no approximation. \(H\) is a conserved quantity of an ecological model in the same sense that energy is a conserved quantity of a mechanical one, and it is dimensionally a rate: each term has units of inverse time. The scaled \(V\) of Step 8 is \(H/r\) up to an additive constant. B

\[ V=\mu g(x)+g(y),\quad g(u)=u-\ln u,\quad g''(u)=\frac{1}{u^{2}}\gt 0,\quad g(u)\ \xrightarrow[\ u\to 0^{+}\ \text{or}\ u\to\infty\ ]{}\ \infty,\qquad \nabla V=\left(\mu\left(1-\tfrac1x\right),\,1-\tfrac1y\right)=0\iff(x,y)=(1,1) \]

Why the level curves are closed. \(g\) is strictly convex on \((0,\infty)\) with its unique minimum at \(u=1\), so \(V\) is strictly convex on the open quadrant with unique minimum \(V_{\min}=\mu+1\) at \((1,1)\), and \(V\) is coercive: it tends to \(+\infty\) at every boundary of the quadrant and at infinity. Hence for each \(c\gt V_{\min}\) the sublevel set \(\{V\le c\}\) is a compact convex body contained in the open quadrant, and \(\nabla V\ne 0\) on its boundary, so the level set \(\{V=c\}\) is a smooth simple closed curve. The vector field is tangent to that curve (Step 9) and vanishes nowhere on it (the only interior zero is the centre, where \(V=V_{\min}\)). A nowhere-zero flow on a compact one-dimensional manifold traverses the whole circle in finite time, so the orbit is the level curve and the solution is periodic. Every non-equilibrium orbit in the open quadrant is therefore a closed cycle, and \(V-V_{\min}\) is a positive-definite Lyapunov function with \(\dot V\equiv 0\): the equilibrium is Lyapunov stable, because every neighbourhood contains a whole sublevel set, but it is not asymptotically stable, because nothing approaches it. C

\[ u=\ln N,\quad v=\ln P:\qquad \dot u=r-ae^{v}=-\frac{\partial \mathcal{H}}{\partial v},\qquad \dot v=\varepsilon ae^{u}-m=\frac{\partial \mathcal{H}}{\partial u},\qquad \mathcal{H}(u,v)=\varepsilon ae^{u}-mu+ae^{v}-rv \]

The structural reason for everything above. In logarithmic coordinates the system is Hamiltonian, with \(v\) playing the role of a coordinate and \(u\) of its conjugate momentum, and \(\mathcal{H}\) is the invariant \(H\) of Step 9 rewritten. Two consequences follow immediately. The flow is area-preserving in the \((u,v)\) plane — its divergence \(\partial\dot u/\partial u+\partial\dot v/\partial v\) is identically zero — so no region of state space can contract or expand, and an asymptotically stable equilibrium or an attracting limit cycle would require contraction. Equivalently, the Dulac function \(1/NP\) makes the divergence of the rescaled field vanish identically — the degenerate case of Bendixson’s criterion, in which closed orbits are not forbidden but can never be isolated, so a limit cycle is impossible — that a whole continuum of cycles is present is the separate content of Step 10. And conservative systems are structurally unstable: an arbitrarily small perturbation that gives the divergence a definite sign near the equilibrium converts the whole family of closed orbits into a spiral. That is not a defect of this particular model but a property of the class it belongs to, and it is the single most important thing to know about the neutral cycles. C

\[ n=N-N^{*},\ p=P-P^{*}:\qquad \dot n=-\frac{m}{\varepsilon}p,\quad \dot p=\varepsilon r\,n\ \Longrightarrow\ \ddot n=-rm\,n,\qquad n=A\cos\omega t,\quad p=\frac{A\varepsilon\omega}{m}\sin\omega t,\qquad \omega=\sqrt{rm},\quad T_{0}=\frac{2\pi}{\sqrt{rm}} \]

Small oscillations. Eliminating \(p\) from the linearised pair of Step 7 gives the harmonic oscillator equation, whose angular frequency is the geometric mean of the prey growth rate and the predator death rate — a rate that belongs to neither species alone. The sine against the cosine says the predator deviation lags the prey deviation by exactly a quarter of a period, which is the model’s sharpest testable prediction and the reason a scatter plot of \(P\) against \(N\) runs anticlockwise. Dividing the two amplitudes by their own equilibria gives the relative amplitude ratio \(\sqrt{m/r}\), independent of \(a\) and \(\varepsilon\) as Step 4 requires. For large orbits the period is larger than \(T_{0}\) and is available only as the contour integral \(T=\left(1/r\right)\oint dx/\left[x(1-y)\right]\) over the scaled orbit; it grows without bound as the orbit approaches the axes, because such orbits pass close to the saddle at the origin, where the flow is slow. B

\[ \int_{0}^{T}\frac{d}{dt}\ln N\,dt=\ln\frac{N(T)}{N(0)}=0=\int_{0}^{T}\left(r-aP\right)dt\ \Longrightarrow\ \langle P\rangle\equiv\frac{1}{T}\int_{0}^{T}P\,dt=\frac{r}{a}=P^{*},\qquad \langle N\rangle=\frac{m}{\varepsilon a}=N^{*} \]

The time-average theorem, and the most useful line on the page. Because the solution is periodic (Step 10), \(\ln N\) returns to its starting value, so the integral of its derivative over one period vanishes; the derivative is \(r-aP\) by Step 2, and the result follows on rearranging. The same argument applied to \(\ln P\) gives \(\langle N\rangle=N^{*}\). This is exact for every orbit, however large the amplitude, and it is not an approximation about small oscillations: the equilibrium that the populations never sit at is nevertheless the mean they are obliged to average to. It also converts the model from an unfalsifiable statement about cycles into an arithmetic prediction about long-run means. B

\[ \dot N=rN-aNP-h_{N}N,\quad \dot P=\varepsilon aNP-mP-h_{P}P\ \Longrightarrow\ \langle N\rangle=\frac{m+h_{P}}{\varepsilon a},\qquad \langle P\rangle=\frac{r-h_{N}}{a}\quad\left(h_{N}\lt r\right) \]

Volterra’s principle. Constant per-capita harvesting of either species changes nothing but the parameters: \(r\mapsto r-h_{N}\) and \(m\mapsto m+h_{P}\), so the model retains its form and Step 13 applies unchanged. Three exact statements follow. Killing predators does not lower the average number of predators — \(\langle P\rangle\) does not contain \(h_{P}\) — it raises the average number of prey. Killing prey does not lower the average number of prey; it lowers the average number of predators. And a non-selective agent with \(h_{N}=h_{P}=h\) does both at once, raising \(\langle N\rangle\) and lowering \(\langle P\rangle\). The proviso \(h_{N}\lt r\) is essential: harvest the prey faster than it can grow and the interior equilibrium leaves the positive quadrant, the cycles disappear, and both species collapse. C

Result

\[ \dot N=rN-aNP,\quad \dot P=\varepsilon aNP-mP;\qquad \left(N^{*},P^{*}\right)=\left(\frac{m}{\varepsilon a},\frac{r}{a}\right),\qquad H=\varepsilon aN-m\ln N+aP-r\ln P=\text{const},\qquad \langle N\rangle=N^{*},\ \ \langle P\rangle=P^{*},\qquad T_{0}=\frac{2\pi}{\sqrt{rm}} \]

Reading. The two populations chase each other around a closed curve for ever, at an amplitude set by where they started and never forgotten, because the flow conserves \(H\) exactly. The centre of the curve is fixed by a crossed pair of parameters — the prey level is built from predator traits alone and the predator level contains no predator trait but the attack rate — and although neither population ever sits there, both average exactly to it. The prey peak leads the predator peak by a quarter of a cycle, and the natural frequency \(\sqrt{rm}\) is the geometric mean of the prey’s growth rate and the predator’s death rate.

Scope. Densities in individuals per unit area, rates per unit time, \(a\) in area per predator per unit time, \(\varepsilon\) in predators per prey. One prey, one specialist predator, closed and well mixed, continuous reproduction, large populations, constant parameters, no prey self-limitation and no saturation of predator intake. The closed orbits are exact within these hypotheses and structurally unstable outside them: any density dependence, saturation, generalist feeding, spatial structure or noise replaces neutral cycles by damped oscillation, by a limit cycle, or by extinction. The equilibrium values, the averages and the quarter-cycle lag are far more robust than the perpetual cycling and are what the model should be used to predict.

Corollaries & converses

Crossed control. \(N^{*}=m/\varepsilon a\) contains no prey parameter and \(P^{*}=r/a\) contains no predator death or efficiency parameter. Fertilising the prey’s resource, which raises \(r\), raises the predator average and leaves the prey average untouched: in this model, enrichment feeds the predator. Symmetrically, a predator that evolves a higher conversion efficiency \(\varepsilon\) depresses its prey and gains nothing itself, since \(\varepsilon\) is absent from \(P^{*}\).
Volterra’s principle. Uniform mortality \(h\lt r\) on both species gives \(\langle N\rangle=(m+h)/\varepsilon a\) and \(\langle P\rangle=(r-h)/a\), so the predator-to-prey ratio \(\varepsilon(r-h)/(m+h)\) falls monotonically in \(h\). Suppressing the harvest reverses it, which is the model’s account of the wartime Adriatic landings.
Quarter-cycle lag. To first order in the amplitude the predator lags the prey by \(T/4\), so a phase plot circulates anticlockwise. An observed lead of the predator over the prey falsifies the model outright; a lag near a quarter cycle is weak confirmation, since many oscillators produce one.
Geometric-mean frequency. \(T_{0}=2\pi/\sqrt{rm}\) uses one rate from each species, so a cycle period cannot be assigned to either alone; and because \(a\) and \(\varepsilon\) are absent, the period of small oscillations is unaffected by how efficient the predator is.
Amplitude is a constant of the motion, not a property of the system. The orbit is fixed by \(H\left(N(0),P(0)\right)\); two identical communities started differently cycle differently for ever. Nothing in the model selects an amplitude, which is precisely what makes it useless as an explanation of observed cycle amplitudes.
No extinction, and no attractor. The axes are invariant, so no interior orbit reaches them in finite time; the flow is area-preserving in \(\left(\ln N,\ln P\right)\), so no limit cycle and no attracting equilibrium can exist. Both statements fail under any perturbation that breaks the conservation law.
Converse fails. Out-of-phase cycles are not evidence of Lotka–Volterra dynamics. Seasonally forced populations, delayed density dependence in a single species, epidemiological cycles and consumer–resource models with entirely different functional forms all produce lagged oscillations; distinguishing them requires the quantitative predictions \(\langle N\rangle=m/\varepsilon a\), \(T_{0}=2\pi/\sqrt{rm}\) or the response to harvesting, not the mere existence of a cycle.
Kolmogorov’s generalisation. The qualitative behaviour of \(\dot N=Nf(N,P)\), \(\dot P=Pg(N,P)\) is determined by sign conditions on the partial derivatives of \(f\) and \(g\) rather than by their functional form; Lotka–Volterra is the degenerate member of that family in which \(\partial f/\partial N=0\), and it is the degeneracy that produces the conservation law.

Fails without

Drop the absence of prey self-limitation — the damped regime: with \(\dot N=rN(1-N/K)-aNP\) and \(K\gt N^{*}\), the equilibrium becomes \(\left(N^{*},\ (r/a)(1-N^{*}/K)\right)\) and the Jacobian acquires trace \(-rN^{*}/K\lt 0\) with determinant \(\varepsilon a^{2}N^{*}P^{*}\gt 0\), so the centre becomes a stable spiral and every orbit converges to it. The conserved quantity is destroyed by a term of any size whatever: the neutral cycles are not a robust prediction but a knife-edge one. Worked numerically in Problem 4.
Drop mass action for a saturating intake — the enrichment regime: with a Holling type II response and prey self-limitation (the Rosenzweig–MacArthur model), the prey nullcline becomes a hump-backed parabola and the equilibrium loses stability when it lies to the left of the hump, i.e. when \(N^{*}\lt \left(K-1/a\tau_h\right)/2\). Raising \(K\) — enriching the system — therefore destabilises it through a Hopf bifurcation into a large-amplitude limit cycle whose troughs push both species towards extinction. This is the paradox of enrichment, and it inverts the naive expectation that more resources make a community safer. Worked numerically in Problem 5.
Drop the deterministic idealisation — the stochastic regime: neutral stability means there is no restoring force acting on the amplitude, so demographic noise makes \(H\) execute a random walk with no reflecting boundary; the orbit inevitably wanders close enough to an axis for one population to be lost. Gause’s Didinium–Paramecium cultures ended in the predator eating out the prey and then starving, and persisted only once he added a sediment refuge for the prey or immigrated fresh individuals on a schedule — the deterministic model’s guarantee of persistence is an artefact of continuous densities.
Drop spatial homogeneity — the metapopulation regime: in a patchy arena local pairs go extinct while other patches are at high density, and recolonisation resets them; persistence then depends on dispersal rates rather than on the local dynamics at all. Huffaker’s mite microcosms persisted through many cycles only after barriers to dispersal were introduced, and asynchrony between patches, not stability within them, is what keeps the system alive (Metapopulation dynamics).
Drop continuous time for discrete generations — the divergent regime: the Nicholson–Bailey host–parasitoid model, \(N_{t+1}=\lambda N_{t}e^{-aP_{t}}\), \(P_{t+1}=cN_{t}\left(1-e^{-aP_{t}}\right)\), is the natural discrete analogue and its interior equilibrium is always unstable: oscillations grow until one species is lost. A whole generation of lag turns neutral cycles into divergent ones, so the continuous-time neutrality is itself special.

Common errors

“The eigenvalues are imaginary, so the equilibrium is a centre.” That establishes it for the linearisation only. The fixed point is not hyperbolic, so Hartman–Grobman does not apply and the nonlinear system could spiral either way; the closed orbits need the global argument of Steps 8–10.
“The cycles are stable, so the model explains observed predator–prey cycles.” They are neutrally stable. Amplitude is fixed by the initial condition and remembered for ever, and any perturbation moves the system permanently to a different orbit. Real cycles with a repeatable amplitude require a limit cycle, which this model cannot produce.
“Increasing the prey’s growth rate increases the prey.” \(\langle N\rangle=m/\varepsilon a\) does not contain \(r\). Raising \(r\) raises \(\langle P\rangle=r/a\) instead. The reflex that a species is controlled by its own parameters is exactly what the model exists to correct.
“Culling the predator reduces predator numbers.” Over a full cycle it does not: \(\langle P\rangle=r/a\) is independent of \(m\), so extra predator mortality is compensated by the higher average prey density it creates. Only the harvest term in the prey equation moves \(\langle P\rangle\).
“Spraying an insecticide reduces the pest.” If it kills the pest’s natural enemies too, Step 14 says the average pest density rises to \((m+h)/\varepsilon a\). Secondary pest resurgence is a prediction of the model, not an anomaly.
“\(\varepsilon\) and \(a\) determine the shape of the cycle.” After nondimensionalisation only \(\mu=m/r\) survives (Step 4); \(a\) and \(\varepsilon\) set the axis scales. Two systems with the same \(m/r\) have geometrically similar orbits.
“The populations settle down eventually.” There is no dissipation anywhere in the equations — the flow is area-preserving in logarithmic coordinates (Step 11) — so nothing can settle. If a simulation appears to converge, it is the numerical integrator losing or gaining the invariant \(H\); forward Euler spirals outwards on this system, and a symplectic or explicitly conservative scheme is required.
“Peak prey coincides with peak predation pressure.” Peak prey occurs when \(P=P^{*}\), on the way up for the predator; peak predator occurs a quarter cycle later, when \(N\) has already fallen back to \(N^{*}\). Reading the peaks as simultaneous destroys the mechanism.

Discussion

The equations were written twice within a few years and for different reasons. Alfred Lotka arrived at them from physical chemistry, treating populations as reacting species and publishing the analysis in Elements of Physical Biology in 1925; Vito Volterra, a mathematical physicist, was asked by the marine biologist Umberto D’Ancona to explain why the proportion of predatory fish in Adriatic landings had risen during the First World War and fallen again afterwards, and produced the same system in 1926 together with the averaging theorem that answers the question. That history explains the model’s character: it was never an attempt to fit a time series, but a demonstration that a qualitative pattern in the data followed from an interaction structure alone. Volterra’s principle remains the best example in ecology of a counterintuitive result that is both exactly derivable and directly applicable.

What the model actually contributes to modern ecology is a set of statements about means and structure rather than about cycles. The averaging theorem holds for arbitrarily large amplitude and requires no measurement of the orbit; the crossed dependence of \(N^{*}\) on predator parameters is the seed of the top-down control arguments that dominate food-web ecology; and the failure of the cycles under perturbation is the standard illustration of structural instability in the biological literature. Kolmogorov showed in 1936 that the essential features of consumer–resource dynamics depend only on sign conditions on the per-capita growth functions, which both generalises the model and exposes how special its conservation law is. Rosenzweig and MacArthur’s graphical nullcline analysis in 1963, and Rosenzweig’s paradox of enrichment in 1971, are the direct descendants that ecologists actually use.

The deepest structural fact is that the system is Hamiltonian in the coordinates \(\left(\ln N,\ln P\right)\), which places it in the same class as the frictionless pendulum: a one-parameter family of periodic orbits, a conserved quantity, an area-preserving flow, and no attractors of any kind. Conservative systems in the plane are non-generic, in the precise sense that an arbitrarily small perturbation of the vector field removes the conservation law and with it the family of orbits, so no property that depends on exact conservation should ever be trusted as a prediction about a real community. The reliable content of the model is what survives perturbation: the location of the equilibrium, the time averages to first order, the sign of the response to harvesting, and the phase relation. Conversely, the same degeneracy is what makes the model so useful as a null hypothesis — because it produces cycles from nothing but the interaction, observing cycles in nature is not by itself evidence for any additional mechanism.

Common misconceptions. That the model predicts population cycles in the sense of predicting their amplitude or their repeatability; it predicts neither, and a data set showing a stable cycle amplitude is evidence against pure Lotka–Volterra dynamics and for a limit cycle. That its equilibrium is what populations tend towards; nothing tends towards it, yet everything averages to it. And that the model is a description of specific systems such as the hare and the lynx: the observed hare–lynx cycle involves plant–hare interactions, generalist predators, stress-mediated reproductive effects and a period far more regular than this model can generate, and the honest use of Lotka–Volterra there is as the first term of an explanation, not the explanation.

Worked examples

Example 1. A well-mixed laboratory microcosm contains the ciliate Paramecium as prey and the predatory ciliate Didinium, both counted as individuals per millilitre. Measurements give an intrinsic prey growth rate \(r=2.0\ \mathrm{d^{-1}}\), an attack rate \(a=0.050\ \mathrm{mL\ predator^{-1}\ d^{-1}}\), a conversion efficiency \(\varepsilon=0.40\) predators per prey, and a predator death rate \(m=1.0\ \mathrm{d^{-1}}\) in prey-free medium. The culture is started at \(N_{0}=80\ \mathrm{mL^{-1}}\) and \(P_{0}=40\ \mathrm{mL^{-1}}\). Find the coexistence equilibrium, the period of small oscillations, and the highest and lowest densities each species reaches on this particular orbit.

\[ N^{*}=\frac{m}{\varepsilon a},\qquad P^{*}=\frac{r}{a} \]

Symbols before numbers (Step 5). Check the dimensions first: \(\varepsilon a\) has units \(\mathrm{mL\ prey^{-1}\ d^{-1}}\), so \(m/\varepsilon a\) is a prey density in \(\mathrm{mL^{-1}}\); \(r/a\) has units \(\mathrm{d^{-1}}/\left(\mathrm{mL\ predator^{-1}\ d^{-1}}\right)=\mathrm{predators\ mL^{-1}}\). A

\[ N^{*}=\frac{1.0}{0.40\times0.050}=\frac{1.0}{0.020}=50\ \mathrm{mL^{-1}},\qquad P^{*}=\frac{2.0}{0.050}=40\ \mathrm{mL^{-1}} \]

Numbers. A useful sanity check on the parameter set: at equilibrium each predator consumes \(aN^{*}=0.050\times50=2.5\) prey per day and converts them at \(\varepsilon=0.40\) into \(1.0\) new predators per day, which is exactly \(m\) — the definition of the predator equilibrium, and a plausible intake for a ciliate that divides about once a day. A

\[ \omega=\sqrt{rm}=\sqrt{2.0\times1.0}=1.414\ \mathrm{d^{-1}},\qquad T_{0}=\frac{2\pi}{\omega}=\frac{6.2832}{1.4142}=4.44\ \mathrm{d} \]

The small-amplitude period from Step 12. The predator peak follows the prey peak by \(T_{0}/4=1.11\) days. Because the orbit through the given initial condition is not small, the true period is somewhat longer than \(4.44\) days; \(T_{0}\) is the limiting value as the orbit shrinks onto the equilibrium. A

\[ H=\varepsilon aN-m\ln N+aP-r\ln P,\qquad H_{0}=0.020\left(80\right)-\ln 80+0.050\left(40\right)-2\ln 40 \]

Evaluate the invariant on the initial condition (Step 9); every point of the orbit must return this same value, so it is the equation of the orbit. Note that \(\varepsilon a=0.020\) and that the logarithms are natural. B

\[ H_{0}=1.600-4.3820+2.000-7.3778=-8.1598\ \mathrm{d^{-1}},\qquad H^{*}=1.000-3.9120+2.000-7.3778=-8.2898\ \mathrm{d^{-1}} \]

Numbers, using \(\ln 80=4.3820\), \(\ln 50=3.9120\), \(\ln 40=3.6889\). The value at the equilibrium is smaller, as Step 10 requires since \(H\) has its strict minimum there; the excess \(H_{0}-H^{*}=0.1300\ \mathrm{d^{-1}}\) is the single number that labels this orbit among all the others. B

\[ \dot N=0\iff P=P^{*}\ \Longrightarrow\ \varepsilon aN-m\ln N=H_{0}-\left(aP^{*}-r\ln P^{*}\right)=-8.1598+5.3778=-2.7820 \]

Where the prey extremes are, found without integrating the equations. The prey density is stationary exactly on the line \(P=P^{*}\), so the largest and smallest prey densities on the orbit are the two roots of the invariant restricted to that line. The bracket is \(0.050(40)-2\ln 40=2.000-7.3778=-5.3778\). B

\[ N=50u:\quad u-\ln u=1.1300\ \Longrightarrow\ u=1.600\ \text{or}\ u=0.5729\ \Longrightarrow\ N_{\max}=80.0\ \mathrm{mL^{-1}},\quad N_{\min}=28.6\ \mathrm{mL^{-1}} \]

Solve the transcendental equation in scaled form. Writing \(N=N^{*}u\) turns \(0.020N-\ln N=-2.7820\) into \(u-\ln u=-2.7820+\ln 50=1.1300\), whose two roots straddle the minimum of \(u-\ln u\) at \(u=1\); the upper root is \(1.600\) because the starting point already lay on the line \(P=P^{*}\), and the lower is found by Newton iteration from \(u=0.6\). The prey therefore ranges over a factor of \(2.8\). B

\[ \dot P=0\iff N=N^{*}\ \Longrightarrow\ aP-r\ln P=-8.1598+2.9120=-5.2478;\quad P=40v:\ v-\ln v=1.0650\ \Longrightarrow\ v=1.405,\ 0.6814 \]

The same construction on the other nullcline: the predator is stationary where \(N=N^{*}\), and the first bracket there is \(0.020(50)-\ln 50=-2.9120\). Dividing the resulting equation \(0.050P-2\ln P=-5.2478\) by \(2\) and substituting \(P=40v\) gives \(v-\ln v=1.0650\). B

\[ P_{\max}=40\left(1.405\right)=56.2\ \mathrm{mL^{-1}},\qquad P_{\min}=40\left(0.6814\right)=27.3\ \mathrm{mL^{-1}},\qquad \langle N\rangle=50\ \mathrm{mL^{-1}},\quad \langle P\rangle=40\ \mathrm{mL^{-1}} \]

The predator extremes, and the averages from Step 13, which need no information about the orbit at all. Observe that the arithmetic means of the extremes, \(54.3\) for prey and \(41.8\) for predators, are not the time averages: the orbit is not symmetric about its centre, and the population spends longer at low density than at high. C

\[ \left(N^{*},P^{*}\right)=\left(50,\,40\right)\ \mathrm{mL^{-1}},\quad T_{0}=4.44\ \mathrm{d},\quad N\in\left[28.6,\,80.0\right],\quad P\in\left[27.3,\,56.2\right]\ \mathrm{mL^{-1}},\quad \langle N\rangle=50,\ \langle P\rangle=40\ \mathrm{mL^{-1}} \]

Reading. Starting the culture at \(80\) prey and \(40\) predators per millilitre commits it to one particular closed loop for ever: prey between \(28.6\) and \(80.0\) per millilitre, predators between \(27.3\) and \(56.2\), circulating with a period a little over four and a half days and averaging exactly to the equilibrium the culture never occupies.

Scope. Deterministic, well-mixed, constant parameters. A real Didinium–Paramecium culture at these densities contains only thousands of individuals per millilitre and the predator saturates at high prey density, so the observed outcome is normally one or two large cycles ending in the extinction of the prey and then of the predator; the calculation above is what the idealisation predicts, and the discrepancy is the point of the Fails without section.

Example 2. An aphid pest, density \(N\) in aphids per square metre, is held in check by a predatory ladybird larva, density \(P\) per square metre, with \(r=0.80\ \mathrm{wk^{-1}}\), \(a=0.20\ \mathrm{m^{2}\ predator^{-1}\ wk^{-1}}\), \(\varepsilon=0.010\) predators per aphid and \(m=0.20\ \mathrm{wk^{-1}}\). A broad-spectrum insecticide is then applied on a schedule that imposes an extra per-capita death rate \(h=0.30\ \mathrm{wk^{-1}}\) on both species. Find the long-run average densities before and after spraying, the change in the predator-to-prey ratio, and the largest \(h\) for which the system persists. Then find what a perfectly selective aphicide, acting only on the pest at the same rate, would do.

\[ \langle N\rangle=N^{*}=\frac{m}{\varepsilon a},\qquad \langle P\rangle=P^{*}=\frac{r}{a} \]

The unsprayed averages, quoted from Step 13. No initial condition is needed and no cycle has to be observed: the theorem fixes the means for every orbit. A

\[ \langle N\rangle=\frac{0.20}{0.010\times0.20}=\frac{0.20}{0.0020}=100\ \mathrm{m^{-2}},\qquad \langle P\rangle=\frac{0.80}{0.20}=4.0\ \mathrm{m^{-2}} \]

Numbers. The check of Step 2 of Example 1 applies again: each larva takes \(aN^{*}=0.20\times100=20\) aphids per week, and \(\varepsilon\times20=0.20\ \mathrm{wk^{-1}}=m\). A

\[ h_{N}=h_{P}=h\ \Longrightarrow\ r\mapsto r-h=0.80-0.30=0.50\ \mathrm{wk^{-1}},\qquad m\mapsto m+h=0.20+0.30=0.50\ \mathrm{wk^{-1}} \]

Spraying does not change the structure of the model (Step 14): a constant per-capita death rate is absorbed into the linear term of each equation. The attack rate and the conversion efficiency are untouched, since the insecticide alters mortality rather than foraging. B

\[ \langle N\rangle_{h}=\frac{m+h}{\varepsilon a}=\frac{0.50}{0.0020}=250\ \mathrm{m^{-2}},\qquad \langle P\rangle_{h}=\frac{r-h}{a}=\frac{0.50}{0.20}=2.5\ \mathrm{m^{-2}} \]

The sprayed averages. The pest average rises by a factor of \(2.5\) and the predator average falls by a factor of \(1.6\), even though the insecticide kills aphids at \(0.30\) per week — a mortality that, applied to aphids alone in the absence of the predator, would look like effective control. The extra prey mortality is entirely absorbed by the released predator population, and then some. C

\[ \frac{\langle P\rangle}{\langle N\rangle}=\frac{\varepsilon\left(r-h\right)}{m+h}:\qquad \frac{4.0}{100}=0.040\ \longrightarrow\ \frac{2.5}{250}=0.010,\qquad \text{a fall by a factor of }4.0 \]

The ratio in symbols first, so the mechanism is visible: the numerator falls and the denominator rises, so the ratio is strictly decreasing in \(h\) and the effect compounds. This is Volterra’s principle in its testable form, and it is the same arithmetic that predicts a rise in the predatory fraction of a fish catch when fishing stops. B

\[ \langle P\rangle_{h}\gt 0\iff h\lt r=0.80\ \mathrm{wk^{-1}};\qquad h\ge 0.80\ \Longrightarrow\ \dot N\le -aNP\lt 0\ \text{and both species are lost} \]

The persistence condition. For \(h\ge r\) the prey cannot maintain itself even with no predators, the interior equilibrium leaves the positive quadrant, and the model collapses to extinction rather than to control — so the perverse result of Steps 4–5 is confined to the range \(0\lt h\lt 0.80\ \mathrm{wk^{-1}}\). Chemical control that succeeds does so by leaving this regime entirely, which requires killing the pest faster than it reproduces. B

\[ h_{N}=0.30,\ h_{P}=0:\qquad \langle N\rangle=\frac{m}{\varepsilon a}=100\ \mathrm{m^{-2}}\ \text{(unchanged)},\qquad \langle P\rangle=\frac{r-h_{N}}{a}=\frac{0.50}{0.20}=2.5\ \mathrm{m^{-2}} \]

The selective aphicide, from the general formulae of Step 14. It leaves the average pest density exactly where it was and cuts the natural enemy by nearly two fifths, from \(4.0\) to \(2.5\ \mathrm{m^{-2}}\): the aphids killed are precisely compensated by the reduced predation that follows, and the only lasting effect of the treatment is on the beneficial insect. The comparison is stark — a selective pesticide is useless here and a non-selective one is worse than useless — and it is the model’s argument for biological rather than chemical control. C

\[ \text{unsprayed: }\left(100,\,4.0\right)\ \mathrm{m^{-2}};\quad h=0.30:\ \left(250,\,2.5\right)\ \mathrm{m^{-2}};\quad \frac{\langle P\rangle}{\langle N\rangle}:\ 0.040\to 0.010;\quad \text{persistence for }h\lt 0.80\ \mathrm{wk^{-1}} \]

Reading. A spray that kills three tenths of both populations per week leaves two and a half times as many aphids as before, because it removes the predator’s numerical response faster than it removes the pest. A spray selective for the pest changes the average pest density not at all.

Scope. Constant per-capita mortality averaged over the spray schedule, one prey and one specialist predator, and the averages taken over a whole number of cycles. Pulsed application, predator immigration from unsprayed refuges, insecticide resistance in either species (Antibiotic resistance is the same selection argument in microbes) and generalist predators all modify the numbers; the sign of the effect, which is what the principle asserts, is robust to all of them provided the predator remains dependent on this prey.

Problems

An island supports moose (prey, \(N\) individuals) and wolves (predator, \(P\) individuals) with \(r=0.30\ \mathrm{yr^{-1}}\), \(a=0.020\ \mathrm{wolf^{-1}\ yr^{-1}}\), \(\varepsilon=0.010\) wolves per moose and \(m=0.40\ \mathrm{yr^{-1}}\). Find the equilibrium, verify the dimensions of \(a\) and \(\varepsilon\), compute the period of small oscillations and the kill rate per wolf at equilibrium, and comment on whether that kill rate is biologically plausible.
Solution
Equilibrium, in symbols first: \(N^{*}=m/\varepsilon a\), \(P^{*}=r/a\). Numerically \(N^{*}=0.40/(0.010\times0.020)=0.40/2.0\times10^{-4}=2000\) moose and \(P^{*}=0.30/0.020=15\) wolves.
Dimensions. In \(aNP\) the product must be moose per year, and \(NP\) is moose \(\times\) wolves, so \([a]=\mathrm{wolf^{-1}\,yr^{-1}}\) as given. In \(\varepsilon aNP\) the product must be wolves per year, so \(\varepsilon\) carries wolves per moose and is here \(0.010\), i.e. one wolf produced per hundred moose consumed.
Period: \(\omega=\sqrt{rm}=\sqrt{0.30\times0.40}=\sqrt{0.12}=0.3464\ \mathrm{yr^{-1}}\), so \(T_{0}=2\pi/0.3464=18.1\) years, with the wolf peak trailing the moose peak by \(T_{0}/4=4.5\) years.
Kill rate: each wolf takes \(aN^{*}=0.020\times2000=40\) moose per year. That is roughly an order of magnitude above the kill rates wolves are actually observed to achieve, which are of order a few moose per wolf per year, and the reason is structural rather than a bad parameter estimate: with a linear functional response the intake per predator is proportional to prey density and cannot saturate, so any parameter set that reproduces a realistic predator equilibrium \(P^{*}=r/a\) forces an unrealistically high intake at a realistic prey equilibrium. Fixing this requires the Holling type II form of Problem 5.
For the microcosm of Example 1 (\(r=2.0\ \mathrm{d^{-1}}\), \(m=1.0\ \mathrm{d^{-1}}\), \(\varepsilon=0.40\), \(a=0.050\ \mathrm{mL\,d^{-1}}\)), work out the small-amplitude solution explicitly: write \(N=N^{*}+n\), \(P=P^{*}+p\), solve the linearised system with \(n(0)=A\) and \(p(0)=0\), and find (a) the lag in days between the prey and predator peaks, (b) the ratio of the fractional amplitudes \(\left(p_{\max}/P^{*}\right)\big/\left(n_{\max}/N^{*}\right)\), and (c) the direction of circulation in the \((N,P)\) plane.
Solution
Linearising about the equilibrium (Step 7) gives \(\dot n=-aN^{*}p=-(m/\varepsilon)p\) and \(\dot p=\varepsilon aP^{*}n=\varepsilon r\,n\). Differentiating the first and substituting the second: \(\ddot n=-(m/\varepsilon)(\varepsilon r)n=-rm\,n\), the harmonic oscillator with \(\omega=\sqrt{rm}=1.4142\ \mathrm{d^{-1}}\).
With \(n(0)=A\) and \(p(0)=0\) the solution is \(n=A\cos\omega t\) and, from \(p=-\dot n\varepsilon/m\), \(p=(A\varepsilon\omega/m)\sin\omega t\). Numerically \(\varepsilon\omega/m=0.40\times1.4142/1.0=0.5657\), so \(p=0.5657A\sin\omega t\).
(a) The sine lags the cosine by a quarter period. \(T_{0}=2\pi/1.4142=4.443\) d, so the lag is \(1.111\) days — about \(27\) hours.
(b) \(p_{\max}/P^{*}=0.5657A/40=0.014142A\) and \(n_{\max}/N^{*}=A/50=0.020A\); the ratio is \(0.7071\). In symbols this is \(\left(\varepsilon\omega/m\right)\left(N^{*}/P^{*}\right)=\omega/r=\sqrt{m/r}=\sqrt{0.5}=0.7071\), independent of \(a\) and \(\varepsilon\) as the nondimensionalisation of Step 4 demands. The predator oscillates relatively less than the prey whenever \(m\lt r\).
(c) At \(t=0\) the prey is at its maximum and the predator at its equilibrium; a quarter period later the predator is at its maximum and the prey at its equilibrium. Plotting \(P\) against \(N\), the point moves from (right, centre) to (centre, top), which is anticlockwise. A data set circulating clockwise — predator leading prey — is inconsistent with the model at any parameter values.
Prove the time-average theorem \(\langle N\rangle=N^{*}\), \(\langle P\rangle=P^{*}\) from the equations, then apply it to a fishery. A trawl fishery takes a prey fish and its predator non-selectively at a common per-capita rate \(h\). With \(r=1.00\ \mathrm{yr^{-1}}\) and \(m=0.50\ \mathrm{yr^{-1}}\), compare pre-war fishing at \(h_{1}=0.40\ \mathrm{yr^{-1}}\) with wartime fishing at \(h_{2}=0.10\ \mathrm{yr^{-1}}\), and give the factor by which the predator-to-prey ratio in the sea changes.
Solution
Proof. Divide the prey equation by \(N\), which is legitimate because \(N\gt 0\) for all time (the axes are invariant, Step 6): \(\dfrac{d}{dt}\ln N=r-aP\). Integrate over one full period \(T\). The left side gives \(\ln N(T)-\ln N(0)=0\) because the solution is periodic (Step 10). Hence \(0=rT-a\displaystyle\int_{0}^{T}P\,dt\), so \(\langle P\rangle=\frac{1}{T}\int_{0}^{T}P\,dt=\frac{r}{a}=P^{*}\). Applying the identical argument to \(\dfrac{d}{dt}\ln P=\varepsilon aN-m\) gives \(0=\varepsilon a T\langle N\rangle-mT\), so \(\langle N\rangle=m/\varepsilon a=N^{*}\). No approximation and no restriction on amplitude is used.
Fishery. Non-selective harvesting maps \(r\mapsto r-h\) and \(m\mapsto m+h\) (Step 14), so \(\dfrac{\langle P\rangle}{\langle N\rangle}=\dfrac{(r-h)/a}{(m+h)/\varepsilon a}=\dfrac{\varepsilon\left(r-h\right)}{m+h}\).
Pre-war: \(\varepsilon(1.00-0.40)/(0.50+0.40)=\varepsilon(0.60)/(0.90)=0.667\varepsilon\). Wartime: \(\varepsilon(1.00-0.10)/(0.50+0.10)=\varepsilon(0.90)/(0.60)=1.500\varepsilon\).
The ratio rises by a factor \(1.500/0.667=2.25\) when fishing is cut from \(0.40\) to \(0.10\ \mathrm{yr^{-1}}\), and falls back by the same factor when fishing resumes. Both averages move: \(\langle N\rangle\) falls by a factor \(0.90/0.60=1.5\) and \(\langle P\rangle\) rises by \(0.90/0.60=1.5\). This is the qualitative pattern D’Ancona reported from the Adriatic landings, and the point of the calculation is that it follows from the interaction structure alone, with no assumption whatever about which species the war favoured. Note also that \(\varepsilon\) cancels out of the factor, so the prediction is testable without knowing the conversion efficiency.
Add prey self-limitation: \(\dot N=rN\left(1-N/K\right)-aNP\), \(\dot P=\varepsilon aNP-mP\), with the Example 1 parameters and \(K=200\ \mathrm{mL^{-1}}\). (a) Find the interior equilibrium and the condition on \(K\) for it to exist. (b) Compute the Jacobian there, its trace and determinant, and its eigenvalues. (c) Give the quasi-period and the time for the amplitude to fall by a factor \(e\), and say by what factor the amplitude decays over one cycle. (d) State the structural conclusion.
Solution
(a) The predator equation is unchanged, so \(\varepsilon aN^{*}=m\) still gives \(N^{*}=m/\varepsilon a=50\ \mathrm{mL^{-1}}\). Setting \(\dot N=0\) with \(N=N^{*}\): \(r(1-N^{*}/K)=aP^{*}\), so \(P^{*}=\dfrac{r}{a}\left(1-\dfrac{N^{*}}{K}\right)=40\left(1-\dfrac{50}{200}\right)=40\times0.75=30\ \mathrm{mL^{-1}}\). The equilibrium is interior only if \(K\gt N^{*}=50\ \mathrm{mL^{-1}}\); a prey ceiling below the predator’s break-even density starves the predator out.
(b) With \(f=rN(1-N/K)-aNP\): \(\partial f/\partial N=r-2rN/K-aP\), which at the equilibrium equals \(r-2rN^{*}/K-r(1-N^{*}/K)=-rN^{*}/K\), and \(\partial f/\partial P=-aN^{*}\). With \(g=\varepsilon aNP-mP\): \(\partial g/\partial N=\varepsilon aP^{*}\) and \(\partial g/\partial P=\varepsilon aN^{*}-m=0\).
Numerically \(-rN^{*}/K=-2.0\times50/200=-0.500\ \mathrm{d^{-1}}\), \(-aN^{*}=-2.50\), \(\varepsilon aP^{*}=0.020\times30=0.600\). So \(\operatorname{tr}J=-0.500\ \mathrm{d^{-1}}\) and \(\det J=0-(-2.50)(0.600)=1.500\ \mathrm{d^{-2}}\).
Eigenvalues: \(\lambda=\tfrac12\left(\operatorname{tr}J\pm\sqrt{\operatorname{tr}^{2}J-4\det J}\right)=\tfrac12\left(-0.500\pm\sqrt{0.250-6.000}\right)=-0.250\pm 1.199i\ \mathrm{d^{-1}}\). Both real parts are negative, so the equilibrium is a stable spiral; in general \(\operatorname{tr}J=-rN^{*}/K\lt0\) and \(\det J=\varepsilon a^{2}N^{*}P^{*}\gt0\) for any \(K\gt N^{*}\), so it is always stable.
(c) Quasi-period \(T=2\pi/1.199=5.24\) d, longer than the undamped \(4.44\) d of Example 1. The amplitude envelope is \(e^{-0.250t}\), so it falls by \(e\) in \(1/0.250=4.00\) d, and over one quasi-period by \(e^{-0.250\times5.24}=e^{-1.31}=0.27\) — roughly a quarter per cycle, so the oscillation is visually gone after four or five cycles.
(d) The undamped cycles of the pure model are destroyed by a self-limitation term of any strength: as \(K\to\infty\) the damping rate \(rN^{*}/2K\to0\) continuously, so there is no threshold below which the neutral cycles survive. The Lotka–Volterra centre is structurally unstable, and its perpetual oscillation should never be quoted as a prediction about a real community.
Rosenzweig–MacArthur: replace mass action by a Holling type II response, \(\dot N=rN\left(1-N/K\right)-\dfrac{aNP}{1+a\tau_h N}\), \(\dot P=\dfrac{\varepsilon aNP}{1+a\tau_h N}-mP\), with \(r=2.0\ \mathrm{d^{-1}}\), \(a=0.050\ \mathrm{mL\,d^{-1}}\), \(\varepsilon=0.40\), \(m=1.0\ \mathrm{d^{-1}}\) and handling time \(\tau_h=0.10\ \mathrm{d}\) per prey. (a) Find the maximum intake per predator and the condition on \(\varepsilon\) for the predator to persist at all. (b) Find \(N^{*}\). (c) The prey nullcline is \(P=\dfrac{r}{a}\left(1-\dfrac{N}{K}\right)\left(1+a\tau_h N\right)\); locate its maximum and hence find the carrying capacity \(K_{c}\) at which the equilibrium loses stability. (d) Say what happens at \(K=500\ \mathrm{mL^{-1}}\) and name the phenomenon.
Solution
(a) As \(N\to\infty\) the intake \(aN/(1+a\tau_h N)\to 1/\tau_h=10\) prey per predator per day, so handling time caps consumption however dense the prey. The predator’s per-capita growth rate is bounded above by \(\varepsilon/\tau_h-m\), so persistence requires \(\varepsilon\gt m\tau_h=1.0\times0.10=0.10\); here \(\varepsilon=0.40\), comfortably above.
(b) Set the predator per-capita rate to zero: \(\varepsilon aN^{*}=m\left(1+a\tau_h N^{*}\right)\), so \(N^{*}\left(\varepsilon a-ma\tau_h\right)=m\) and \(N^{*}=\dfrac{m}{a\left(\varepsilon-m\tau_h\right)}=\dfrac{1.0}{0.050\left(0.40-0.10\right)}=\dfrac{1.0}{0.0150}=66.7\ \mathrm{mL^{-1}}\). Saturation has pushed the predator’s break-even prey density up from \(50\) to \(66.7\ \mathrm{mL^{-1}}\).
(c) Differentiate the nullcline: \(P(N)=\dfrac{r}{a}\left[1+a\tau_h N-\dfrac{N}{K}-\dfrac{a\tau_h N^{2}}{K}\right]\), so \(\dfrac{dP}{dN}=\dfrac{r}{a}\left[a\tau_h-\dfrac{1}{K}-\dfrac{2a\tau_h N}{K}\right]=0\) at \(N_{h}=\dfrac{K}{2}-\dfrac{1}{2a\tau_h}=\dfrac{K-1/a\tau_h}{2}\). Here \(1/a\tau_h=1/(0.050\times0.10)=200\ \mathrm{mL^{-1}}\), so \(N_{h}=(K-200)/2\).
The standard nullcline criterion is that the equilibrium is stable when it lies on the descending arm of the hump (\(N^{*}\gt N_{h}\)) and unstable on the ascending arm (\(N^{*}\lt N_{h}\)); the Hopf bifurcation is at \(N^{*}=N_{h}\), i.e. \(K_{c}=2N^{*}+\dfrac{1}{a\tau_h}=2(66.7)+200=333\ \mathrm{mL^{-1}}\).
(d) At \(K=500\ \mathrm{mL^{-1}}\) the hump sits at \(N_{h}=(500-200)/2=150\ \mathrm{mL^{-1}}\), well to the right of \(N^{*}=66.7\), so the equilibrium is unstable and the system settles onto a stable limit cycle of large amplitude whose prey trough is far below \(N^{*}\). Enriching the environment — raising \(K\) from \(333\) to \(500\) — has therefore destabilised a stable community and driven both species close to extinction at the troughs. This is Rosenzweig’s paradox of enrichment. Note the contrast with Problem 4: with a linear response, self-limitation always stabilises; with a saturating response, it stabilises only while \(K\) is small enough, and the sign of the effect of enrichment reverses.