5 Single species models

Learning outcomes

By the end of this week you should be able to:

• derive the exponential (Malthusian) population model from first principles;

• derive and solve the logistic model, including nondimensionalisation;

• find equilibrium points of autonomous ODEs and determine their stability using linearisation;

• analyse the effect of constant harvesting on a population model, including critical harvesting rates and bifurcation diagrams.

5.1 Introduction

In this part of the course we examine models that describe the behaviour of populations of a single species in terms of population size. In most cases the real population is measured in discrete units (it does not make much sense to have a third of a rabbit), but the equations we derive describe the behaviour in some average sense and we treat population size as a continuous variable. This is a subtlety worth bearing in mind when we construct mathematical models.

5.2 The exponential model

We discussed this simple model in the introductory chapter, now we will motivate it from first principals, for example consider a colony of rabbits with plenty of space and food. Rabbits do not migrate into or out of the colony and there are no predators.

5.2.0.0.1 Variables and parameters.

Let \(N(t)\) denote the number of rabbits at time \(t\), measured in appropriate units (weeks or months). In a time interval \([t,\, t+h]\), the change in population is entirely due to births and deaths: \[N(t+h) - N(t) = \text{births in } [t, t+h] - \text{deaths in } [t, t+h].\]

The birth rate \(\beta\) is the average number of births per unit time per member of the population. The death rate \(\delta\) is defined analogously. In any time interval, the number of births or deaths is (on average) \[\text{length of interval} \times \text{rate} \times \text{population size}.\]

5.2.0.0.2 Forming the model.

With birth rate \(\beta\) and death rate \(\delta\), the change in population between times \(t\) and \(t+h\) is \[N(t+h) - N(t) = \beta\,N(t)\,h - \delta\,N(t)\,h = \alpha\,N(t)\,h,\] where \(\alpha := \beta - \delta\) is the growth rate. Dividing by \(h\) and taking the limit \(h \to 0\) using the definition of the derivative gives the ODE \[\begin{equation} \frac{\dd N}{\dd t} = \alpha\,N. \label{eq:exponential} \end{equation}\] This is the exponential or Malthusian population model.¹

5.2.0.0.3 Solving the model.

As we saw earlier equation [eq:exponential] is separable. For \(N \neq 0\): \[\int \frac{\dd N}{N} = \int \alpha\,\dd t \qquad\Longrightarrow\qquad \log N = \alpha t + C.\] With the initial condition \(N(0) = N_0\), we find \(C = \log N_0\), giving \[\begin{equation} \boxed{N(t) = N_0\,e^{\alpha t}.} \label{eq:exp-solution} \end{equation}\] The solution depends on two parameters: the growth rate \(\alpha\) and the initial population \(N_0\).

5.3 Variations on the exponential model

A more general exponential model is \[\begin{equation} \frac{\dd N}{\dd t} = (\beta(N,t) - \delta(N,t))\,N, \label{eq:exp-general} \end{equation}\] where \(\beta\) and \(\delta\) can now be functions of population size and/or time. This allows us to model many interesting real-world phenomena that are not captured by constant birth/death coefficients.

5.3.1 Worked example: periodic growth

Consider the model \[\frac{\dd N}{\dd t} = -\alpha\cos(\omega t)\,N, \qquad N(0) = N_0.\] What could this model be used to describe? What are the maximum and minimum growth rates? What is the period of the growth rate? Solve the model using the initial condition \(N(0) = N_0\). What are the minimum and maximum values of the population size?

What is the value of \(\omega\) if we wish to model a seasonal phenomenon and \(t\) is measured in days?

A colony of wasps is reduced to size 250 in the depths of winter (when the growth rate of the colony is minimum), at what point in the year is the maximum growth rate? At what point in the year does the maximum population size occur? If the maximum growth rate of the colony is 1 per week, what is the maximum population size? ²

Solution: The model could be used to describe a seasonal growth rate in a population. The minimum growth rate is \(-\alpha\) and the maximum growth rate is \(\alpha\) . The period of the growth rate is \(T = 2\pi/\omega\). The differential equation is separable and the solution is:

\[\boxed{N(t) = N_0\exp\!\left(-\frac{\alpha}{\omega}\sin\omega t\right).}\]

The population oscillates between a minimum of \(N_0 e^{-\alpha/\omega}\) and a maximum of \(N_0 e^{\alpha/\omega}\).

For a yearly seasonal variation with \(t\) in days, \(\omega = 2\pi/365 \approx 0.0172\;\mathrm{day}^{-1}\). If a colony of wasps is reduced to \(N_0 = 250\) in the depths of winter (when the growth rate is minimum, at \(t = 0\)), the maximum growth rate occurs mid-year (June/July) and the maximum population size occurs in September/October (out of phase with the growth rate). If the maximum growth rate is \(\alpha = 1/7\;\mathrm{day}^{-1}\), the peak population is \[N_{\max} \approx 250\,e^{(1/7)/0.0172} \approx 10^6.\]

5.4 The logistic model

The exponential model assumes unlimited resources. In reality, there is a maximum population size that the environment can support — the carrying capacity \(M\). As the population approaches \(M\), the growth rate must approach zero. The simplest way to model this is \[\alpha(N) = \alpha\!\left(1 - \frac{N}{M}\right),\] where \(\alpha\) is the linear growth rate (the growth rate when \(N \ll M\)). This gives the logistic model: \[\begin{equation} \boxed{\frac{\dd N}{\dd t} = \alpha\!\left(1 - \frac{N}{M}\right)N,} \label{eq:logistic} \end{equation}\] where \(\alpha > 0\) and \(M > 0\) are constants. This model was first studied by the Belgian mathematician Pierre François Verhulst in the mid-19th century.

Note two important limits:

• As \(M \to \infty\) (for fixed \(N\)), the growth rate \(\alpha(N) \to \alpha\), recovering the exponential model with infinite resources.

• For \(N \ll M\), the growth rate is close to \(\alpha\) — resources appear plentiful. The finiteness of resources only becomes noticeable when \(N\) is comparable to \(M\).

5.4.0.0.1 Solving the logistic model.

Rewrite [eq:logistic] as \(M\,\dd N/\dd t = \alpha(M - N)N\). The equilibrium solutions are \(N(t) = 0\) and \(N(t) = M\). For \(N \neq 0, M\), separate variables and use partial fractions: \[\int \frac{M\,\dd N}{(M-N)N} = \int \alpha\,\dd t \qquad\Longrightarrow\qquad \int \left(\frac{1}{N} + \frac{1}{M-N}\right)\dd N = \alpha t + C.\] Evaluating: \[\log\left|\frac{N}{M-N}\right| = \alpha t + C.\] Applying the initial condition \(N(0) = N_0\) and noting that \(N/(M-N)\) does not change sign along an orbit, we obtain \[\begin{equation} \boxed{N(t) = \frac{MN_0}{(M - N_0)\,e^{-\alpha t} + N_0}.} \label{eq:logistic-solution} \end{equation}\] This includes the equilibrium solutions (\(N_0 = 0\) gives \(N = 0\); \(N_0 = M\) gives \(N = M\)). The solutions depend on three parameters: \(N_0\), \(M\) and \(\alpha\). The resulting curves are called logistic curves.

5.4.0.0.2 Nondimensionalisation.

The logistic model has four dimensional quantities: \(\alpha\) (dimensions \([\mathrm{T}]^{-1}\)), \(N\) and \(M\) (dimensions \([\mathrm{N}]\)), and \(t\) (dimensions \([\mathrm{T}]\)). The natural time scale is \(\alpha^{-1}\) and the natural population scale is \(M\). Define \[\hat{t} = \alpha t, \qquad \hat{N} = \frac{N}{M}.\] The equation becomes \[\frac{\dd\hat{N}}{\dd\hat{t}} = \hat{N}(1 - \hat{N}), \qquad \hat{N}(0) = \hat{N}_0 = \frac{N_0}{M}.\] The only remaining parameter is \(\hat{N}_0\), the initial population as a fraction of the carrying capacity. All logistic solutions are essentially the same, the only differences are in scale.

5.5 Qualitative analysis of autonomous ODEs

In this section we will introduce some terms and ideas that are particularly useful when we move on to consider systems of differential equations that model populations of two or more species, or more generally systems of ODEs. Introducing these ideas in the simpler context of single species models will allow us to test the results against the previously obtained exact solutions. Consider a general autonomous first-order ODE \[\begin{equation} \frac{\dd y}{\dd t} = f(y). \label{eq:autonomous} \end{equation}\]

The value \(y(t)\) can be represented as a point on the real line, which we call the phase space. An orbit is the set of all points visited by the solution: \(\{y(t) : t \in \mathbb{R}\}\).

To determine the direction of motion along an orbit, note that for small \(h > 0\):

• \(\dot{y} > 0 \;\Longrightarrow\; y(t+h) > y(t)\) (moves to the right);

• \(\dot{y} = 0 \;\Longrightarrow\; y(t+h) = y(t)\) (stationary);

• \(\dot{y} < 0 \;\Longrightarrow\; y(t+h) < y(t)\) (moves to the left).

The sign of \(\dot{y}\) is determined directly from the ODE, since \(\dot{y} = f(y)\).

Any point \(y_e\) where \(f(y_e) = 0\) is called an equilibrium point (or fixed point), since a solution at \(y_e\) remains there for all time. The information about equilibrium points and the direction of orbits is displayed in a phase portrait: equilibrium points are marked by dots and arrows indicate the direction of motion.

5.5.0.0.1 The logistic model.

For \(\dd N/\dd t = \alpha(1 - N/M)N\), we have \(f(N) > 0\) for \(0 < N < M\) and \(f(N) < 0\) for \(N > M\). The equilibrium points are \(N = 0\) and \(N = M\). From the phase portrait:

5.6 Stability and linearisation

For the general autonomous ODE [eq:autonomous], let \(y_e\) be an equilibrium point. We say \(y_e\) is:

• stable if all solutions starting close to \(y_e\) remain close for all \(t\);

• asymptotically stable if all solutions starting close to \(y_e\) satisfy \(y(t) \to y_e\) as \(t \to \infty\);

• unstable otherwise.

For the logistic model, \(M\) is asymptotically stable and \(0\) is unstable.

5.6.0.0.1 Linearisation.

Let \(y_e\) be an equilibrium point with \(f(y_e) = 0\). Define the perturbation \(\tilde{y} = y - y_e\). Taylor-expanding \(f(y)\) about \(y_e\): \[f(y) = f(y_e) + (y - y_e)\,f'(y_e) + \tfrac{1}{2}(y - y_e)^2\,f''(y_e) + \cdots = f'(y_e)\,\tilde{y} + \text{h.o.t.}\] Since \(\dd\tilde{y}/\dd t = \dd y/\dd t\), for \(\tilde{y}\) small we obtain the linearisation of [eq:autonomous] about \(y_e\): \[\begin{equation} \frac{\dd\tilde{y}}{\dd t} = f'(y_e)\,\tilde{y}. \label{eq:linearisation} \end{equation}\] This is a linear ODE with solution \(\tilde{y}(t) = \tilde{y}(0)\,e^{f'(y_e)\,t}\). Provided \(f'(y_e) \neq 0\), the stability is determined entirely by the sign of \(f'(y_e)\):

5.7 Harvesting

To describe the effects of hunting, predation or migration out of a population, we add a harvesting term: \[\begin{equation} \frac{\dd N}{\dd t} = \alpha(N)\,N - H(N, t), \label{eq:harvesting} \end{equation}\] where \(H(N,t)\) represents the number of members removed per unit time. Common choices include:

• \(H = H_0\) (constant): a fixed number removed per unit time.

• \(H = H_1 N\) (proportional): a fixed fraction of the population removed.

• \(H = BN^2/(A^2 + N^2)\) (saturating): harvesting that “switches on” when \(N\) reaches \(A\) and saturates at rate \(B\).

5.8 Worked example: constant harvesting

5.9 Constant harvesting of the logistic model

Consider the logistic model with constant harvesting: \[\begin{equation} \frac{\dd N}{\dd t} = \alpha\!\left(1 - \frac{N}{M}\right)N - H, \qquad H \ge 0. \label{eq:logistic-harvest} \end{equation}\]

The equilibrium points are the roots of \(N^2 - MN + MH/\alpha = 0\), which has discriminant \(\Delta = M^2 - 4MH/\alpha\). Three cases arise:

5.9.0.0.1 Moderate harvesting (\(H < M\alpha/4\)).

Two equilibrium points exist: \[E_{\pm} = \frac{M \pm \sqrt{M^2 - 4MH/\alpha}}{2}.\] From the phase portrait, \(E_+\) is asymptotically stable and \(E_-\) is unstable. For all \(N_0 > E_-\), the population approaches \(E_+\). If \(N_0 < E_-\), the population goes extinct.

5.9.0.0.2 Critical harvesting (\(H = M\alpha/4\)).

There is a single equilibrium at \(E = M/2\). This is a borderline case and conclusions drawn from it are unreliable since the model is itself an approximation.

5.9.0.0.3 Excessive harvesting (\(H > M\alpha/4\)).

No equilibrium points exist. For all initial populations, \(N(t) \to 0\) — the population goes extinct.

5.10 Critical harvesting for general models

For a general autonomous model with constant harvesting, \[\begin{equation} \frac{\dd N}{\dd t} = f(N) - H, \label{eq:general-harvest} \end{equation}\] if \(f(N) - H < 0\) for all \(N \ge 0\) then every population goes extinct. The critical harvesting rate is \[H_c = \max_{N \ge 0} f(N).\]

5.11 Bifurcation diagrams

A bifurcation diagram displays the location and stability of equilibrium points as a function of a parameter. For the logistic model with constant harvesting, nondimensionalise by writing \(N = Mn\) and \(t = \tau/\alpha\): \[\frac{\dd n}{\dd\tau} = n(1 - n) - h, \qquad h = \frac{H}{M\alpha}.\] The single dimensionless parameter \(h\) controls the behaviour. The equilibrium points satisfy \(n(1-n) = h\), giving \[n_{\pm} = \frac{1 \pm \sqrt{1 - 4h}}{2}.\]

The bifurcation diagram plots \(n_{\pm}\) against \(h\). For \(h < 1/4\) there are two equilibria (\(n_+\) stable, \(n_-\) unstable). At \(h = 1/4\) they coalesce — this is a saddle-node bifurcation. For \(h > 1/4\) there are no equilibria and all populations go extinct.