10 Invasion waves: populations and epidemics
When a species is introduced into a new area, or a disease reaches a new population, it does not appear everywhere at once. It spreads as a wave that advances at a roughly constant speed. Grey squirrels, muskrats, cane toads, and rabies in European foxes all spread like this. In this final week we find out what sets the speed of such waves. We then combine the epidemic models of week 8 with the spatial models of week 9 to predict how fast an epidemic sweeps across a country, and we finish by looking back at the whole module.
By the end of this week you should be able to:
- reduce a reaction–diffusion equation to an ODE with the travelling-wave Ansatz and state the boundary conditions of an invasion wave;
- derive the minimum wave speed by linearising at the leading edge of the wave;
- formulate and non-dimensionalise a spatial SIR model;
- derive the fraction of susceptibles left behind an epidemic wave and its minimum speed;
- recognise the common toolkit (steady states, linearisation, thresholds) used throughout the module.
Open the week 10 notebook in Colab
10.1 Travelling waves
In the launch investigation you explored the Fisher–KPP equation \[ \partial_tu=ru(1-u)+D\,\partial_x^2u. \tag{10.1}\]
- Describe how the population spread. Did the front speed up, slow down or settle to a constant speed? Did the shape of the front change?
- Combine your tables of measured speeds. Agree on a formula for the speed in terms of \(r\) and \(D\). Check that it has the right units.
- The muskrat front advanced about 11 km per year, and muskrats have an intrinsic growth rate of about 1 per year. Estimate \(D\).
We look for solutions of fixed shape moving to the right with speed \(c\): \(u(x,t)=U(z)\) with \(z=x-ct\).
- Show that \(\partial_tu=-cU'\) and \(\partial_xu=U'\), and hence that Eq. 10.1 becomes the ODE \[DU''+cU'+f(U)=0,\qquad f(U)=rU(1-U).\]
- What should \(U(z)\) tend to as \(z\to-\infty\) (behind the front) and as \(z\to+\infty\) (ahead of it)?
- At the leading edge \(U\) is small. Linearise the ODE there and try \(U=e^{-\lambda z}\). Show that \(D\lambda^2-c\lambda+r=0\).
- If \(\lambda\) were complex, \(U\) would oscillate around \(0\) at the leading edge. Why is this impossible for a population density? Deduce that \[c\geq c_{\min}=2\sqrt{rD}.\]
- Compare \(c_{\min}\) with your measurements.
You will be asked to convert reaction–diffusion equations to travelling-wave form, state the boundary conditions and find minimum wave speeds.
With \(z=x-ct\), the reaction–diffusion equation \(\partial_tu=f(u)+D\partial_x^2u\) becomes \(DU''+cU'+f(U)=0\). An invasion wave connects the invaded state behind (\(U(-\infty)=1\)) with the empty state ahead (\(U(+\infty)=0\)), with \(U'(\pm\infty)=0\). At the leading edge, \(DU''+cU'+f'(0)U\approx0\) with solutions \(e^{-\lambda z}\), where \[ \lambda=\frac{c\pm\sqrt{c^2-4Df'(0)}}{2D}. \] A non-negative front needs real \(\lambda\), so \(c\geq2\sqrt{Df'(0)}\). Waves started from a localised population travel at the minimum speed \(c_{\min}=2\sqrt{rD}\). They approach it slowly from below, which is why the measured speeds were a little smaller. The speed is set entirely by the leading edge, where the population is small, grows at rate \(r\) and diffuses. It is like an epidemic seeded again and again just ahead of the front. For the muskrat, \(D\approx c^2/(4r)\approx30\) km² per year.
Without growth (\(f=0\)) the wave ODE \(DU''+cU'=0\) has solutions \(U=A+Be^{-cz/D}\). These blow up as \(z\to-\infty\), so there is no realistic travelling wave: diffusion alone spreads a population like \(\sqrt t\), not at a constant speed.
10.2 The spread of an epidemic
Rabies reached France from Eastern Europe in the 1960s and advanced at 30–60 km per year. Suppose it arrived in Dover. How soon would it reach York? Healthy foxes are territorial and stay put, but rabid foxes wander around randomly. So only the infected foxes diffuse: \[ \partial_tS=-\beta\frac{SI}N,\qquad \partial_tI=\beta\frac{SI}N-\gamma I+D\partial_x^2I, \tag{10.2}\] where \(S(x,t)\) and \(I(x,t)\) are densities and \(N\) is the (uniform) fox density. For rabies, the removed class consists of dead foxes.
Introduce \(u=S/N\), \(v=I/N\), \(\tilde t=\beta t\), \(\tilde x=\sqrt{\beta/D}\,x\) and \(r=\gamma/\beta\). Show that (dropping the tildes) \[ \partial_tu=-uv,\qquad \partial_tv=uv-rv+\partial_x^2v. \tag{10.3}\] How is \(r\) related to \(R_0\)?
- Watch the epidemic wave for \(r=0.5\). Describe the profiles of \(u\) and \(v\).
- Measure the speed \(c\) and the fraction \(a\) of susceptibles left behind the wave for several values of \(r\). What happens for \(r>1\)?
- Compare \(a\) with the final size of a (non-spatial) SIR epidemic with the same \(R_0\). What do you notice?
Make the Ansatz \(u=U(z)\), \(v=V(z)\), \(z=x-ct\).
- Show that \(-cU'=-UV\) and \(-cV'=UV-rV+V''\). What are the boundary conditions as \(z\to\pm\infty\)? Write \(a=U(-\infty)\).
- Use the first equation to write \(UV=cU'\) and \(V=c(\ln U)'\), and substitute into the second to show that every term is a derivative: \[-cV'=cU'-rc(\ln U)'+V''.\]
- Integrate from \(-\infty\) to \(\infty\) and use the boundary conditions to show that \[a-1=r\ln a.\] Compare with the final-size equation Eq. 8.3.
- Linearise the \(V\) equation at the leading edge (where \(U\approx1\), \(V\) small) and show that \(c\geq2\sqrt{1-r}\).
- Convert back to dimensional variables to show that the minimum speed is \(2\sqrt{D(\beta-\gamma)}\).
Take \(D=50\) km² per year, an infectious period of 5 days (\(\gamma\approx73\) per year) and \(\beta=80\) per year.
- Compute \(R_0\), the minimum wave speed, and the time the epidemic needs to travel the 350 km from Dover to York.
- What fraction of the fox population survives behind the wave?
- A proposal is to cull foxes ahead of the wave, reducing their density \(N\) by 20%. Assume that the contact rate \(\beta\) is proportional to the fox density (each fox meets more foxes where there are more of them). How does culling change \(\beta\), \(R_0\) and the wave speed? What reduction of the fox density would stop the wave altogether?
You will be asked to non-dimensionalise spatial epidemic models, find the fraction left behind a wave and compute minimum wave speeds.
The travelling-wave equations \(-cU'=-UV\) and \(-cV'=UV-rV+V''\), with \(U(\infty)=1\), \(U(-\infty)=a\) and \(V(\pm\infty)=0\), can be written as \(-cV'=cU'-rc(\ln U)'+V''\). Integrating over the whole line gives \[ a-1=r\ln a, \] which is exactly the final-size equation with \(S_0=N\) and \(R_0=1/r=\beta/\gamma\). Behind the wave, each location has experienced a complete local epidemic. At the leading edge \(V''+cV'+(1-r)V=0\), so real exponents need \[ c\geq2\sqrt{1-r}\quad\text{(dimensionless)},\qquad c\geq2\sqrt{D(\beta-\gamma)}\quad\text{(dimensional)}. \] A wave exists only if \(r<1\), i.e. \(R_0>1\). Just as for the Fisher wave, the speed is \(2\sqrt{D\times(\text{initial growth rate})}\), since \(\beta-\gamma\) is the initial growth rate of the epidemic.
For the group task, \(R_0=80/73\approx1.10\), \(c\approx2\sqrt{50\cdot7}\approx37\) km per year, so it takes about 9.4 years to reach York, and about 83% of foxes survive (\(a\) solves \(a-1=0.9125\ln a\)). If \(\beta\propto N\), reducing \(N\) by 20% gives \(\beta=64<\gamma\), so the wave is stopped. In general, culling must reduce \(N\) so that \(\beta<\gamma\), here by at least \(1-73/80\approx9\%\).
10.3 Looking back: one toolkit, many models
Over ten weeks we have met many models. Remarkably, a small set of ideas was enough to analyse all of them.
On a large sheet of paper, make a concept map of the module. Include at least: steady state, linearisation, eigenvalue, stability, bifurcation, threshold, reproduction number, \(R_0\), critical domain size, wave speed, separation of variables, non-dimensionalisation. Connect each concept with the models in which it appeared, and label the links. Then answer:
- In how many different models did a “threshold equals 1” condition decide between growth and decline? List them with their threshold quantity.
- In how many models did we linearise about the extinction state (or the disease-free state), and why is that such a powerful idea?
- Which models could you not have analysed without a computer? What did the computer add?
- Steady states and linearisation. For \(\dot N=f(N)\): sign of \(f'(N^*)\). For maps: \(|f'(N^*)|<1\). For planar systems: trace and determinant of the Jacobian. For PDEs: separation of variables or Fourier modes, followed by the same criteria for each mode.
- Thresholds. Growth versus decline is decided by a single number crossing \(1\) (or a rate crossing \(0\)): \(r\) in the exponential model, \(R\) in the Euler–Lotka equation, \(R_0\) for epidemics, \(\gamma_0>0\) for the critical reserve size, the existence of real exponents for the minimum wave speed. Very often this comes from linearising about the extinction or disease-free state: can a small population invade?
- Bifurcations happen where the linearisation has a zero eigenvalue, or \(-1\) for maps, or a purely imaginary pair (Hopf). They mark qualitative changes, such as collapse, tipping points, oscillations and patterns.
- Conserved quantities and exact trajectories: Lotka–Volterra, the SIR phase plane, the final-size equation, which reappears behind epidemic waves.
- Non-dimensionalisation reveals which parameter combinations matter.
The final exam will test exactly these skills: formulate a model, find and classify its steady states, find thresholds and bifurcations, and interpret the results ecologically, all by hand.
10.4 Exam-style questions
A population satisfies \(\partial_tu=f(u)+D\partial_x^2u\) with \(D>0\), \(f(0)=f(1)=0\) and \(f(u)>0\) for \(0<u<1\).
Convert the equation into travelling-wave form with \(z=x-ct\). [2 marks]
Assume a solution with \(U(-\infty)=1\) and \(U(\infty)=0\) exists. By linearising at the leading edge, show that a biologically realistic solution needs \(c\geq2\sqrt{Df'(0)}\). [5 marks]
Suppose instead \(f(u)=0\). Find the general solution of the travelling-wave equation. Is it realistic? [3 marks]
\(\partial_t=-c\frac{d}{dz}\) and \(\partial_x=\frac{d}{dz}\) give \(-cU'=f(U)+DU''\).
Near \(U=0\), \(f(U)\approx f'(0)U\), so \(DU''+cU'+f'(0)U=0\). With \(U=e^{-\lambda z}\), \(D\lambda^2-c\lambda+f'(0)=0\) and \(\lambda=\frac{c\pm\sqrt{c^2-4Df'(0)}}{2D}\). Complex \(\lambda\) would make \(U\) oscillate about \(0\), taking negative values. So we need \(c^2\geq4Df'(0)\).
\(DU''+cU'=0\) gives \(U=A+Be^{-cz/D}\). This is unbounded as \(z\to-\infty\) unless \(B=0\), so there is no front connecting \(1\) to \(0\): not realistic.
Consider \[ \partial_tu=-u^2v,\qquad \partial_tv=u^2v-\rho v+\partial_x^2v,\qquad \rho>0. \]
With \(u=U(z)\), \(v=V(z)\), \(z=x-ct\), show that \(-cV'=cU'-\frac{\rho c}{U^2}U'+V''\). [4 marks]
For a solution with \(U(\infty)=1\), \(U(-\infty)=a<1\) and \(V(\pm\infty)=0\), find \(a\). [4 marks]
Find a lower bound on \(c\). [3 marks]
The first equation gives \(-cU'=-U^2V\), so \(V=cU'/U^2\) and \(U^2V=cU'\). Substituting into \(-cV'=U^2V-\rho V+V''\) gives the result.
Integrate over \(\mathbb R\): \(0=c(1-a)+\rho c\left(1-\frac1a\right)\), using \(\int\frac{U'}{U^2}dz=\left[-\frac1U\right]\). Since \(a\neq1\), \(1=\rho/a\), so \(a=\rho\) (requires \(\rho<1\)).
At the leading edge \(U\approx1\): \(-cV'=(1-\rho)V+V''\). With \(V=e^{-\lambda z}\), \(\lambda^2-c\lambda+(1-\rho)=0\), and real \(\lambda\) needs \(c\geq2\sqrt{1-\rho}\).
Without rabies, foxes grow logistically. With rabies, \[ \partial_tS=rS\left(1-\frac SK\right)-\beta SI+D_1\partial_x^2S,\qquad \partial_tI=\beta SI-\gamma I+D_2\partial_x^2I. \]
Ignoring diffusion, show that with \(u=S/K\), \(v=I/K\) and a suitable \(\tilde t\) the model becomes \(\dot u=bu(1-u)-uv\), \(\dot v=uv-mv\), and give \(b\) and \(m\). [3 marks]
Find the steady states. When is there an endemic state? [3 marks]
With diffusion, find \(\tilde x\) and \(d\) such that \(\partial_tu=bu(1-u)-uv+\partial_x^2u\), \(\partial_tv=uv-mv+d\,\partial_x^2v\). [2 marks]
With \(u=A(z)\), \(v=B(z)\), \(z=x-ct\), write down the travelling-wave ODEs. [2 marks]
If \(A(\infty)=1\), what are \(A(-\infty)\), \(B(\infty)\) and \(B(-\infty)\)? Sketch \(A\) and \(B\). [3 marks]
Find a lower bound on \(c\). [3 marks]
\(\tilde t=\beta Kt\), \(b=r/(\beta K)\), \(m=\gamma/(\beta K)\).
\((0,0)\), \((1,0)\) and the endemic state \((m,\,b(1-m))\), which exists iff \(m<1\), i.e. \(\beta K>\gamma\) (\(R_0>1\)).
\(\tilde x=x\sqrt{\beta K/D_1}\) and \(d=D_2/D_1\).
\(-cA'=bA(1-A)-AB+A''\) and \(-cB'=AB-mB+dB''\).
Ahead of the wave there is no disease and the foxes are at carrying capacity: \(A(\infty)=1\), \(B(\infty)=0\). Behind it the endemic state is established: \(A(-\infty)=m\), \(B(-\infty)=b(1-m)\). \(A\) drops from \(1\) to \(m\) across the front (possibly with damped oscillations), while \(B\) rises from \(0\), overshoots and settles at \(b(1-m)\). The wave moves to the right.
At the leading edge \(A\approx1\), \(B\) small: \(dB''+cB'+(1-m)B=0\), so \(c\geq2\sqrt{d(1-m)}\).
Summary
- Travelling-wave Ansatz: \(u(x,t)=U(x-ct)\) turns a reaction–diffusion PDE into an ODE. Invasion waves connect the invaded and the empty states.
- Minimum wave speed from the leading edge: \(c_{\min}=2\sqrt{Df'(0)}\), so for the Fisher–KPP equation \(c=2\sqrt{rD}\).
- Epidemic waves (rabies): the fraction left behind solves the final-size equation \(a-1=r\ln a\), and the minimum speed is \(2\sqrt{D(\beta-\gamma)}\). A wave requires \(R_0>1\).
- A small toolkit (steady states, linearisation, eigenvalues, thresholds, bifurcations, non-dimensionalisation) underlies every model in this module.