\documentclass[12pt,a4paper,linespace=1.2]{article}
\usepackage{graphicx}
\usepackage{amsmath,amssymb}
\usepackage[small]{caption}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{hyperref}
\usepackage{enumerate}
\topmargin -1.5cm
\textwidth 16.59cm
\textheight 22.94cm
\oddsidemargin -0.04cm
\evensidemargin -0.04cm
\raggedbottom
\newcommand{\vsp}{\vspace{0.2cm}}
\newcommand{\tight}{\vspace{-0.1cm}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\nn}{\nonumber}
\title{The dynamics of two-species models}
\date{}
\author{Sam Dolan (Dec 2016)}
\begin{document}
\maketitle
\section{Introduction}
\noindent \textbf{Background}: Attempts to model population dynamics using mathematical methods can be traced back to the 18th century, and the writings of Robert Malthus \cite{Malthus}. The Malthusian model of exponential growth at a constant rate, $x(t) = x_0 e^{a t}$, arises as the solution to a simple differential equation, $\dot{x} = a x$, where $x$ is the population, $\dot{x} \equiv dx/dt$ is its rate of change with respect to time $t$, and $a$ is a constant parameter. Malthus' writings influenced the work of both the naturalist Charles Darwin and the economist Adam Smith.
In 1838, Verhulst \cite{Verhulst} refined Malthus' simple model to include the effect of resource limitations, introducing the logistic equation $\dot{x} = a x (1 - bx)$.
In the 1920s, attention turned to the dynamics of a \emph{pair} of interacting species. Independently, A.~J.~Lotka \cite{Lotka} and V.~Volterra \cite{Volterra} introduced the predator-prey (Lotka-Volterra) equations: an autonomous pair of first-order ordinary differential equations (ODEs)
\begin{eqnarray}
\dot{x} &=& \phantom{-} a x - b x y , \nn \\
\dot{y} &=& - c y + d x y , \label{eq:predprey1}
\end{eqnarray}
which purport to model the population dynamics of a prey species (e.g.~rabbits) interacting with a predator species (e.g.~foxes). Here, $x$ is the prey population; $y$ is the predator population; and $a,b,c,d$ are non-negative constant parameters, with $a$ the natural growth rate of the rabbit population in the absence of foxes; $c$ the natural decay rate of the fox population in the absence of rabbits (due to hunger); and $b$ and $d$ modelling the effects of predation, which boosts the fox population and diminishes the rabbit population.
The predator-prey equations (\ref{eq:predprey1}) are based on several idealised assumptions. For example, it is assumed that there is sufficient food available for the prey species, regardless of its population size; that predators have an unlimited appetite for the prey; and that the predators consume only one species. Weakening these assumptions leads to more complex models. In this report we will investigate both the original predator-prey equations (\ref{eq:predprey1}), and some modifications. \vsp
\noindent \textbf{(b) Mathematical results}:
\textbf{Fixed points}: A fixed point (or equilibrium) occurs where the populations are not changing, $\dot{x} = 0 =\dot{y}$. One fixed point is at $\bar{x} = 0 = \bar{y}$, with both species extinct. Another fixed point, at $\bar{y} = a / b$ and $\bar{x} = c / d$, can be found by first factorizing the right-hand sides of (\ref{eq:predprey1}) into $x (a - by)$ and $-y(c - dx)$ and then setting to zero. To classify the fixed points of a 2D system $\dot{x} = f(x,y)$, $\dot{y} = g(x,y)$, we examine the eigenvalues of the Jacobian matrix $J$. At the non-zero fixed point, the Jacobian matrix is
\begin{equation}
J = \left. \begin{bmatrix}
\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\
\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}
\end{bmatrix} \right._{x = \bar{x}, y = \bar{y}}
=
\begin{bmatrix}
0 & - b c / d \\
a d / b & 0
\end{bmatrix}.
\end{equation}
%
Its eigenvalues are found from the characteristic equation $\text{det}\| J - \lambda I \|= 0$. In this case, both eigenvalues are pure imaginary, $\lambda_1 = i \sqrt{a c}$, $\lambda_2 = \lambda_1^\ast$, and therefore fixed point is a centre. Equivalently, because $\text{Tr}(J) = 0$ and $\text{det} J > 0$ it follows that the fixed point is a centre. This implies that nearby orbits in phase space move around the fixed point in closed loops. \vsp
\textbf{A conservation law}. Let
\begin{equation}
C = a \ln y - b y + c \ln x - d x. \label{eq:C}
\end{equation}
Taking a derivative with respect to time $t$,
\begin{eqnarray}
\frac{dC}{dt} &=& \frac{a \dot{y}}{y} - b \dot{y} + \frac{c \dot{x}}{x} - d \dot{x} = \left(a - by\right) \frac{\dot{y}}{y} + \left(c - dx\right) \frac{\dot{x}}{x} \nonumber \\
&=& (a - by)(-c + dx) + (c - dx)(a - by) = 0 .
\end{eqnarray}
In the second line, I have used the differential equations (\ref{eq:predprey1}).
As $\dot{C} = 0$, it follows that $C$ is a constant along a given orbit in phase space (N.B.~different orbits will have different constants). \vsp
\textbf{Averaged quantities}. The mean average $\left< z \right>$ of a periodic quantity $z(t)$ is defined by $\left< z \right> = \frac{1}{T} \int_0^T z(t) dt$. Here $T$ is the period, such that $z(t+T) = z(t)$. Let us define $X(t) = x - \bar{x}$ and $Y(t) = y - \bar{y}$. The differential equations (\ref{eq:predprey1}) can be rewritten as $\dot{x}/x = -b Y$ and $\dot{y}/y = d X$. Thus
\begin{eqnarray}
\left< x \right> &=& \frac{1}{T} \int_0^T \left(\bar{x} + X\right) dt = \bar{x} + \frac{1}{Td} \int_0^T \frac{\dot{y}}{y} dt = \bar{x} + \frac{1}{Td} \int_{y_0}^{y_1} \frac{dy}{y} \nonumber \\
&=& \bar{x} + \frac{1}{T d} \left[ \ln y \right]^{y_1}_{y_0} = \bar{x} .
\end{eqnarray}
Here the integral is zero because the trajectory is periodic, $y_1 \equiv y(T) = y(0) \equiv y_0$. A similar argument follows for $\left$. Hence the temporal averages $\left$ and $\left$ are equal to the fixed-point values $\bar{x}$ and $\bar{y}$.
\vsp
\noindent \textbf{(c) Rescaled equations}:
\emph{Non}-\emph{dimensionalization} in mathematical modelling is the removal of (dimensionful) units from a set of equations by a suitable change of variables. This can be achieved in our case by rescaling the variables to eliminate as many constant parameters as possible. Let $x = \alpha \tilde{x}$, $y = \beta \tilde{y}$ and $t = \gamma \tilde{t}$, where $\alpha, \beta, \gamma$ are constants. After substituting into Eq.~(\ref{eq:predprey1}), and multiplying the upper equation by $\gamma / \alpha$ and the lower equation by $\gamma / \beta$, we obtain
\begin{equation}
\frac{d\tilde{x}}{d\tilde{t}} = (a \gamma) \, \tilde{x} - (b \beta \gamma) \, \tilde{x} \tilde{y} , \quad \quad
\frac{d\tilde{y}}{d\tilde{t}} = - (c \gamma) \, \tilde{y} + (d \alpha \gamma) \, \tilde{x} \tilde{y} .
\end{equation}
Any three of the quantities in parentheses can be made equal to unity by a suitable choice of the rescaling constants $\alpha$, $\beta$, $\gamma$. We will select $\gamma = 1/a$, $\beta = a / b$ and $\alpha = a / d$, so that
\begin{equation}
\frac{d\tilde{x}}{d\tilde{t}} = \tilde{x} - \, \tilde{x} \tilde{y} , \quad \quad
\frac{d\tilde{y}}{d\tilde{t}} = - f \, \tilde{y} + \, \tilde{x} \tilde{y} , \label{eq:predprey2}
\end{equation}
where $f = c \gamma = c / a$. Note that $f$ is dimensionless, unlike the original constants $a$ and $c$ whose units were inverse time (i.e.~frequency). Henceforth I will omit the tilde notation, for clarity.
\section{The Lotka-Volterra model}
I used Python code to numerically solve the system (\ref{eq:predprey2}) for given initial conditions $(x_0, y_0)$ and parameter $f$. My code used the {\tt odeint()} function in the module {\tt scipy.integrate}. Plots were created with {\tt matplotlib} and {\tt pyplot.plot()}.
Figure \ref{fig:part2a}(i) shows a typical evolution of (rescaled) prey and predators populations, as a function of (rescaled) time. The populations display regular oscillations with a certain period $T$ (with $T \approx 5$ in Fig.~\ref{fig:part2a}(i)). A trough in the fox population anticipates a peak in the rabbit population; similarly, a peak in the fox population leads to a trough in the rabbit population. The peak in the fox population lags the peak in the rabbit population by approximately one quarter-cycle. The `shape' of the cyclical response depends on the choice of initial conditions $(x_0, y_0)$ and the parameter $f$.
Figure \ref{fig:part2a}(ii) shows a sample of population evolutions in the phase space $(x,y)$. Each curve shows an evolution with different initial conditions $x_0, y_0$ but the same parameter $f=2$. The plot illustrates that the curves circle around a centre at $(2,1)$. All curves are closed, indicating that the evolution is always periodic.
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{../Code/plot1a.pdf}
\includegraphics[width=8cm]{../Code/plot1b.pdf}
\end{center}
\caption{Predator and prey dynamics for the basic Lotka-Volterra model, Eq.~(\ref{eq:predprey1}).
(i) \emph{Left:} a time-domain plot, showing predator and prey populations as a function of rescaled time. (ii) \emph{Right:} a phase portrait of predators versus prey. The closed loops indicate that the system is periodic. The fixed point (a centre) is at $(f,1)$, with $f=2$ in the plots above.
}
\label{fig:part2a}
\end{figure}
Figure \ref{fig:Cerror} shows numerical data for the constant of motion $C$, defined in Eq.~(\ref{eq:C}). Inevitably, numerical error in evolving the ODE system leads to a (spurious) secular drift in $C$ away from its initial value. However, the rate of this drift is small, and it can be reduced towards a negligible value by adjusting the tolerances in the numerical integrator ({\tt atol} and {\tt rtol}).
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{../Code/C_numerical_error.pdf}
\end{center}
\caption{Showing the effect of numerical error on the constant of motion, Eq.~(\ref{eq:C}).
The plot shows $\Delta C = C - \bar{C}$, where $C$ is the numerically-calculated value and $\bar{C}$ the exact value given by Eq.~(\ref{eq:C}). Here $\Delta C$ is shown as a function of time $t$, for a typical evolution with $x_0 = y_0 = 1$, $f=2$. The three curves show numerical data with the absolute and relative error tolerances ({\tt atol} and {\tt rtol}) in the function {\tt odeint()} set to (i) $\epsilon = 10^{-8}$, (ii) $\epsilon = 10^{-9}$ and $\epsilon = 10^{-10}$.
}
\label{fig:Cerror}
\end{figure}
\section{Refining the model}
Now we investigate a modified set of equations, \emph{viz.},
\begin{eqnarray}
\dot{x} &=& a(t) x (1 - gx) - \frac{xy}{1 + hx} , \nn \\
\dot{y} &=& -f y \quad + \frac{xy}{1 + hx} , \label{eq:refined}
\end{eqnarray}
We can give some biological interpretation to the modified equations. Broadly speaking, the term associated with the parameter $g$ models the effect of a finite food source for the prey species (cf.~the logistic equation), and $a(t)$ allow that food source to vary with time (for example, with the season). The factor of $1/(1 + hx)$ takes into account that the appetite of the predator is limited, or equivalently, that there are diminishing returns in increasing the predator's reproductive success as the prey population grows large. \vsp
\textbf{(a) Attractive fixed point}. In this part we consider Eq.~(\ref{eq:refined}) with $h=0$ and $a(t) = 1$. A 2D system of autonomous ODEs, $\dot{x} = F(x,y)$ and $\dot{y} = G(x,y)$, has a fixed point at $(\bar{x}, \bar{y})$ where $F(\bar{x},\bar{y}) = 0$ and $G(\bar{x}, \bar{y}) = 0$. Factorizing, we find $F(x,y) = x \left[ 1 - g x - y \right]$ and $G(x,y) = y \left[ -f + x \right]$. By setting $F$ and $G$ to zero, we find fixed points at $\bar{x} = \bar{y} = 0$ and at $\bar{x} = f$, $\bar{y} = 1-fg$. The latter is in the positive quadrant ($\bar{x} > 0$, $\bar{y} > 0$) provided that $f > 0$ and $0 < fg < 1$.
To classify the non-zero fixed point, we first compute its Jacobian matrix $J$, where
\begin{equation}
J =
\left. \begin{bmatrix}
\frac{\partial F}{\partial x} & \frac{\partial F}{\partial y} \\
\frac{\partial G}{\partial x} & \frac{\partial G}{\partial y}
\end{bmatrix} \right._{x = \bar{x}, y = \bar{y}}
=
\left. \begin{bmatrix}
1 - 2 g x - y & -x \\
y & x - f
\end{bmatrix} \right._{x = f, y = 1 - fg}
=
\begin{bmatrix}
-f g & -f \\
1 - fg & 0
\end{bmatrix} .
\end{equation}
It follows that $\text{det} J = f (1 - fg)$ and $\text{Tr} J = - f g$. For parameters satisfying $f>0$ and $0 0$ (right plot), all trajectories converge upon an attractive fixed point on the $x$-axis at $(1/g, 0)$. In other words, the fox population dies out, leaving behind a stable rabbit population. \vsp
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{../Code/attractive_fixed_point.pdf}
\includegraphics[width=8cm]{../Code/extinction_fixed_point.pdf}
\end{center}
\caption{
Examples of an attractive fixed point in the two-species model. \emph{Left:} A phase portrait for Eq.~(\ref{eq:refined}) with $a(t) =1$, $f = 1.5$, $g = 0.3$ and $h=0$. Any trajectory starting in the positive quadrant evolves towards the attractive fixed point at $(1.5, 0.55)$, representing a stable equilibrium of the rabbit and fox populations. \emph{Right:} A phase portrait with $f = 1.5$, $g = 0.8$. Trajectories converge on the fixed point at $(1.25, 0)$; in this scenario, the fox population dies out.
}
\label{fig:fixedpoint}
\end{figure}
\textbf{(b) Limit cycle}.
A limit cycle is a closed trajectory in phase space which is isolated, meaning that neighbouring trajectories are not closed. A stable limit cycle is one which attracts nearby trajectories. Figure \ref{fig:limitcycle} shows that, for the parameter choice $f = 8/3$, $g=3/50$ and $h=3/20$, there is a stable limit cycle, upon which the red and green lines converge.
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{../Code/limitcycle.pdf}
\end{center}
\caption{Showing a limit cycle in the phase plane for $f = 8/3$, $g=3/50$ and $h=3/20$.
The green trajectory starts outside the limit cycle, the red trajectory starts inside the limit cycle, and the blue trajectory starts near the repelling fixed point at $(40/9, 11/9) \approx (4.444, 1.222)$ and spirals outwards. All trajectories approach the limit cycle as $t \rightarrow \infty$.
}
\label{fig:limitcycle}
\end{figure}
Using the method of part (a), it is possible to show that there is also fixed point at $\bar{x} = f / (1 - fh)$, $\bar{y} = [1 - f(g+h)]/[1 - fh]^2$ provided that $0 < f(g+h) < 1$. For our parameter choice, the fixed point is at $(40/9, 11/9) \approx (4.444, 1.222)$. Further calculation shows that this fixed point is unstable, as the determinant and trace of the Jacobian are both positive. This can also be seen in Fig.~\ref{fig:limitcycle}.
Due to the stable limit cycle, in this case the fox and rabbit populations will eventually settle into a periodic cycle which is independent of the initial population numbers. \vsp
\textbf{(c) Existence of limit cycles}. A mathematician claims that a unique stable limit cycle exists if the parameters satisfy
\begin{equation}
h > g \frac{1 + f h}{1 - f h} . \label{eq:2c}
\end{equation}
I investigated this claim by looking at several particular cases. I have not been able to find a counter-example to the claim. Figure \ref{fig:2c} shows three cases. In cases (i) and (ii) the parameters chosen satisfy the bound and a stable limit cycle is observed. In case (iii), the parameters do not satisfy the bound, and there is no limit cycle; instead, there is an attractive fixed point.
\begin{figure}
\begin{center}
\includegraphics[width=5cm]{../Code/case1.pdf}
\includegraphics[width=5cm]{../Code/case2.pdf}
\includegraphics[width=5cm]{../Code/case3.pdf}
\end{center}
\caption{Example trajectories in phase space, for parameters $f = 1$ and (i) $g=4/50$, $h=3/20$, (ii) $g = 1/6$,$h =3/8$ and (iii) $g = 1/4$, $h = 1/10$. In cases (i) and (ii), Eq.~(\ref{eq:2c}) is satisfied and there is a stable limit cycle. In case (iii), Eq.~(\ref{eq:2c}) is not satisfied, and there is an attractive fixed point but no limit cycle.
}
\label{fig:2c}
\end{figure}
In case (i), the limit cycle passes close to $(0,0)$, suggesting that this biological system could be precarious, and at risk of an extinction event.
\section{A two-species model with harvesting}
In this section we consider Eq.~(\ref{eq:refined}) with a periodic harvesting term, $a(t) = A (1 + \sin \omega t)$, where $A$ and $\omega$ are real parameters. The period of the harvesting term is $T_0 = 2\pi / \omega$. We may rewrite Eq.~(\ref{eq:refined}) as a 3D autonomous system by introducing an auxilliary variable $z$ governed by $\dot{z} = \omega$. The phase space is then three-dimensional. \vsp
The predator-prey system can exhibit a periodic or aperiodic response to the periodic harvesting, depending on the choice of parameters. Figure \ref{fig:4a} shows examples of both behaviours, for $f=8/3, g=3/50, h=3/20$ and $ \omega=1$. For low $A$, we see evidence of \emph{frequency-locking}: the system exhibits a periodic response with the same period as the harvesting term [Fig.~\ref{fig:4a}(b)]. As the amplitude of harvesting increases, we see successive period-doubling [Fig.~\ref{fig:4a}(b)]. This leads on to aperiodic behaviour for $A \approx 0.63$ [Fig.~\ref{fig:4a}(c)].
\begin{figure}
\begin{center}
\includegraphics[width=5.3cm]{../Code/harvest_a.pdf}
\includegraphics[width=5.3cm]{../Code/harvest_b.pdf}
\includegraphics[width=5.3cm]{../Code/harvest_c.pdf}
\end{center}
\caption{Examples of response of predator-prey system to periodic harvesting.
Parameter choice $f=8/3, g=3/50, h=3/20, \omega=1$. (a) $A = 0.45$: the response is periodic with period $T_0 = 2\pi / \omega$; (b) $A = 0.61$: the response is periodic with period $4T_0$; (c) $A = 0.67$: the response is \emph{aperiodic}. See also Figs.~\ref{fig:4b}--\ref{fig:4d}.
}
\label{fig:4a}
\end{figure}
Figure \ref{fig:4b} shows the same data in the 3D phase space. I have plotted rabbits ($x$) and foxes ($y$) against $\sin \omega t$, so as to show the phase of the harvesting term on the $z$ axis. A phase space trajectory cannot not cross itself, because the 3D ODE system is autonomous. Closed curves [(a) and (b)] correspond to periodic trajectories. Curves that do not close can fill parts of phase space in an intricate way.
\begin{figure}
\begin{center}
\includegraphics[width=5.3cm]{../Code/3dphase_a.pdf}
\includegraphics[width=5.3cm]{../Code/3dphase_b.pdf}
\includegraphics[width=5.3cm]{../Code/3dphase_c.pdf}
\end{center}
\caption{Trajectories in 3D phase space for cases (a), (b) and (c) of Fig.~\ref{fig:4a}.
(a) A closed period-1 curve; (b) a closed period-4 curve; (c) an open curve which does not cross itself yet is tangled in an interesting way. The $z$-axis shows $\sin z = \sin \omega t$, i.e., the phase of the harvesting term.
}
\label{fig:4b}
\end{figure}
To understand this better, I have formed a `Poincar\'e section', whereby I plot a single point in the $x$-$y$ plane every time the condition $z = nT_0$ is met, where $T_0 = 2\pi / \omega$ is the period of the harvesting term. In periodic cases, we see a finite number of points, corresponding to the period of the response in units of $T_0$. Aperiodic cases will generate (in principle) an underending-yet-never-repeating sequence of points. These points trace out some underlying mathematical structure, which may have fractal properties (like the strange attractor in the Lorenz system); but I have not investigated this point in any detail.
\begin{figure}
\begin{center}
\includegraphics[width=10cm]{../Code/poincare.pdf}
\end{center}
\caption{Example of Poincar\'e section / `Surface of section'. This plot shows a 2D slice through the 3D phase space in Fig.~\ref{fig:4b}, determined by $t = nT_0$, $n \in \mathbb{Z}$ $\Leftrightarrow$ $\sin \omega t = 0, \cos \omega t = 1$ (i.e.~the $z=0$, $\dot{z} > 0$ slice of Fig.~\ref{fig:4b}). Cases (a), (b) and (c) of Fig.~\ref{fig:4a} are shown as 1 blue point, 4 red points and many green points, respectively.}
\label{fig:4c}
\end{figure}
Figure \ref{fig:4d} shows evidence that, in case (c), the long-term evolution of the system depends delicately on the initial conditions. I started two evolutions with nearly the same initial conditions. Initially, the response looks very similar in the two cases. However, after $t \approx 200$, the evolutions decouple, and do not appear to be correlated at late times. This behaviour is consistent with a system exhibiting \emph{deterministic chaos}, loosely, an exponential sensitivity to initial conditions in a deterministic system such as a set of ODEs.
\begin{figure}
\begin{center}
\includegraphics[width=9cm]{../Code/deterministic_chaos.pdf}
\end{center}
\caption{Sensitivity of aperiodic response to initial conditions. The plot compares two `nearby' evolutions, with initial conditions $x(0) = 1$, $y(0) = 1$ [blue] and $x(0) = 1 + \epsilon$, $y(0) = 1$ [red], where $\epsilon = 10^{-6}$. Note the decoupling after $t \approx 200$.
}
\label{fig:4d}
\end{figure}
Figure \ref{fig:4e} shows bifurcation diagrams. Here I have plotted the $x$ values (rabbits) in the Poincar\'e section of Fig.~\ref{fig:4d} for a range of harvesting amplitudes $A$. One can read off the period of the system from the number of points shown. In the case $\omega = 1.5$ (left), after several bifurcations in which the period doubles, the system starts to behave in an aperiodic fashion for large $A$. In case $\omega = 1.5$ (right) the harvesting frequency is initially quite different to the natural frequency of the (unharvested) system. For low amplitudes $A$, the frequency is not locked. As the amplitude of harvesting is increased, frequency-locking occurs through a $3 \rightarrow 1$ bifurcation. As $A$ is increased further, the system exhibits intervals of chaos and regularity.
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{../Code/bifurcation_w1_zoom.png}
\includegraphics[width=8cm]{../Code/bifurcation_w1pt5.png}
\end{center}
\caption{\emph{Left:} A bifurcation diagram for the harvested predator-prey system with parameter choice $f=8/3, g=3/50, h=3/20$ and $\omega=1$. For harvesting amplitudes $A \lesssim 0.5$, there is frequency-locking: the fox and rabbit populations oscillate with the same period $T_0$ as the harvester, i.e.~$T_0 = 2 \pi / \omega$. At $A \approx 0.505$, the period of the system doubles to $2T_0$, and it doubles again to $4T_0$ at $A \approx 0.595$. Successive period-doubling leads on to a chaotic (non-periodic) response beyond $A \gtrsim 0.63$. \emph{Right:} A bifurcation diagram for $\omega = 1.5$, with all other parameters as in Fig.~\ref{fig:4e}. See text.
}
\label{fig:4e}
\end{figure}
\textbf{Summary}: We have investigated some of the rich phenomenology of a family of predator-prey equations, showing that the populations may respond in a periodic or chaotic manner, depending on the choice of parameters in the model.
\begin{thebibliography}{20}
{\footnotesize
\bibitem{Malthus}
Malthus, R., \emph{An Essay on the Principle of Population} (London, 1798). \vspace{-0.3cm}
\bibitem{Verhulst}
Verhulst, P.~F., \emph{Notice sur la loi que la population poursuit dans son accroissement}, Correspondance math\'ematique et physique {\bf 10}, 113 (1838). \vspace{-0.3cm}
\bibitem{Lotka}
Lotka, A.J., \emph{Elements of Physical Biology}, Williams and Wilkins (1925). \vspace{-0.3cm}
\bibitem{Volterra}
Volterra, V., \emph{Variazioni e fluttuazioni del numero d'individui in specie animali conviventi}, Mem. Acad. Lincei Roma, {\bf 2}, 31 (1926). \vspace{-0.3cm}
\bibitem{Code} Example code available here: \url{http://sam-dolan.staff.shef.ac.uk/mas212/docs/assign2_code.ipynb}
}
\end{thebibliography}
\end{document}