Predator and Prey

The Lotka-Volterra equations describe two species of animals, a predator and its prey. They will provide us with an example of the use of phase-plane analysis of a $2 \times 2$ nonlinear system. The prey population is $x$, the predator is $y$, and the independent variable $t$ is time.

Without any predators, the prey would undergo exponential growth: $x' = a x$, where $a$ is a positive constant. But predation reduces the growth rate or even makes it negative. The simplest assumption is that the effect of predation is proportional to both $x$ and $y$. Thus $x' = a x - \alpha x y$. As for the predator, without any prey it would starve, and we assume this would happen exponentially: $y' = -b y$. But with prey they can grow. Again the simplest assumption is that the effect on $y'$ is proportional to both $x$ and $y$: $y' = -b y + \beta x y$. Thus we have the Lotka-Volterra equations, where $a$, $b$, $\alpha$ and $\beta$ are positive constants:

$$x' = x (a - \alpha y),\ \ y' = y (-b + \beta x)$$

For an equilibrium point of the equations, we need $x=0$ or $y = a/\alpha$, and $y=0$ or $x=b/\beta$. There are two equilibrium points: $(0,0)$ and $(b/\beta, a/\alpha)$. The Jacobian matrix is $$J(x,y) = \pmatrix{ a-\alpha y & -\alpha x\cr \beta y & -b+\beta x\cr}$$.

At $(0,0)$ the Jacobian matrix is $$\pmatrix{ a & 0\cr 0 & -b\cr}$$, with eigenvalues $a$ and $-b$, one positive and one negative. Thus the origin is a saddle point. The eigenvectors are $$\pmatrix{1\cr 0\cr}$$ for $a$ and $$\pmatrix{0\cr 1\cr}$$ for $-b$. Thus there will be two trajectories entering the equilibrium point vertically (one up and one down) and two leaving horizontally to the right and left.

At $(b/\beta, a/\alpha)$ the Jacobian matrix is $$\pmatrix{0 & -\alpha b/\beta \cr \beta a/\alpha & 0\cr}$$. The characteristic polynomial is $r^2 + ab$. Thus the eigenvalues are pure imaginary: $\pm \sqrt{ab} i$. In the linearization, this equilibrium point is a centre. Unfortunately, this does not tell us whether the equilibrium point in the nonlinear system is a centre or a stable or unstable spiral.

One way to resolve the question is to find an implicit equation that $x$ and $y$ satisfy. We find an equation in $x$ and $y$ by dividing the second equation of our system by the first (pretending that $y$ is a function of $x$):

$$\frac{dy}{dx} = \frac{y'}{x'} = \frac{y (-b + \beta x)}{x (a - \alpha y)}$$

This happens to be separable, so we can solve it:

$$\frac{a - \alpha y}{y}\, dy = \frac{-b+\beta x}{x} \, dx$$

$$a \ln \vert y\vert - \alpha y = - b \ln \vert x\vert + \beta x + C$$

Thus (at least for $x > 0$ and $y > 0$, which is the region we're interested in) the trajectories are level curves of the function $h(x,y) = a \ln y + b \ln x - \alpha y - \beta x$. Thus we should have a centre, not a spiral, at the equilibrium point $(b/\beta, a/\alpha)$.

The isoclines for $x'=0$ are $x=0$ and $y = a/\alpha$. Since $x=0$ is a vertical line, there are trajectories on it (the two trajectories that enter the equilibrium point $(0,0)$). The isoclines for $y'=0$ are $y=0$ and $x=b/\beta$, so again there are trajectories on the $x$ axis (the two that leave the $(0,0)$).

Here is the phase portrait of the system in the first quadrant, in the case $\alpha = \beta = a = b = 1$:

Every trajectory in the first quadrant is a closed curve. Thus the predator and prey populations oscillate go in cycles. Starting with a low population of both predator and prey (in the region $0 < x < b/\beta$, $0 < y < a/\alpha$), at first the prey increase while the predators decrease (but not all the way to 0, because trajectories don't meet) until $x=b/\beta$. As we cross into the region $x > b/\beta$, $0 < y < a/\alpha$ the predators start to increase, while the prey still increase until $y = a/\alpha$. Then we cross into $x > b/\beta$, $y > a/\alpha$, and the prey decrease while the predators increase until $x=b/\beta$. Next, in the region $0 < x < b/\beta$, $y > a/\alpha$, the predators and prey both decrease until $y = a/\alpha$. Finally, in the region $0 < x < b/\beta$, $0 < y < a/\alpha$, the prey increase while predators decrease, and we come back to the start of the cycle.

It is not clear whether the cycles predicted by the Lotka-Volterra equation are observed in nature. There are species, such as lynxes and snowshoe hares in the Canadian arctic, that go through population cycles, but the fit between the data and the equations is not very good, and there may be complicating factors such as disease.

The closed curves of the Lotka-Volterra equations are a very ``delicate'' phenomenon, a consequence of the fact that the equation for $dy/dx$ happened to be separable. Almost any modification to the equation would destroy this property. We'll consider two modifications: harvesting and self-limiting of the prey.

Harvesting

First consider harvesting. We harvest prey at a small but constant rate $\varepsilon > 0$. Thus the equations become

$$x' = x (a - \alpha y) - \varepsilon ,\ \ y' = y (-b + \beta x)$$

Where are the equilibrium points? For $y'=0$ we still have $y=0$ or $x=b/\beta$. If $y=0$, then $x' = a x - \varepsilon = 0$ when $x = \varepsilon /a$. Thus the equilibrium point that was at $(0,0)$ has shifted right to $(\varepsilon /a,0)$. This is easy to understand: $\varepsilon /a$ is just the population of prey for which, in the absence of predators, the natural increase will just balance our harvesting rate $\varepsilon$. On the other hand, with $x=b/\beta$, $x' = (a - \alpha y) b/\beta - \varepsilon = 0$ when $y = a/\alpha - \varepsilon \beta/(b \alpha)$. Thus the second equilibrium point has shifted down to $$(x^\star, y^\star) = \left( \frac{b}{\beta},\, \frac{a}{\alpha} - \frac{\varepsilon \beta}{b \alpha}\right)$$. We assume $\varepsilon$ is small enough that $y^\star > 0$.

The Jacobian matrix $J(x,y)$ is still unchanged (because we have just added a constant to $f(x,y)$), but the equilibrium points have shifted. For the first equilibrium point, $$J(\varepsilon /a, 0) = \pmatrix{a & -\alpha \varepsilon /a\cr 0 & -b + \beta \varepsilon /a \cr}$$. The eigenvalues are $a$, which is positive, and $-b + \beta \varepsilon /a$, which is negative (again, assuming $\varepsilon$ is small). So this equilibrium point is still a saddle. An eigenvector for the positive eigenvalue is still $$\pmatrix{1\cr 0\cr}$$, and in fact the trajectories leaving this saddle point are still on the $x$ axis (it is still true that if $y=0$ then $y'=0$). An eigenvector for the negative eigenvalue is $$\pmatrix{\alpha \varepsilon /a \cr a + b - \beta \varepsilon /a\cr}$$, which (when $\varepsilon$ is small) points up and slightly to the right.

At the equilibrium point $(x^\star,y^\star)$, the Jacobian matrix is $$\pmatrix{ \varepsilon \beta/b & -\alpha b/\beta\cr a\beta/\alpha - \beta^2 \varepsilon /(b\alpha) & 0\cr}$$. The characteristic polynomial is $r^2 - (\varepsilon \beta/b) r + a b - \varepsilon \beta$. We have $p < 0$ and (when $\varepsilon$ is small) $p^2 - 4 q < 0$, so this equilibrium point is an unstable spiral.

Here is the phase portrait, showing the trajectory that enters the saddle point. Note that the isocline for $x'=0$ is the curve $y = a/\alpha - \varepsilon /(\alpha x)$ (shown in blue), which has positive slope and goes through the two equilibrium points. The isoclines for $y'=0$ are still the same as before. Following backwards the trajectory that enters the saddle point from the first quadrant, we see it spirals out of the equilibrium point $(x^\star,y^\star)$. The picture below is for $a = \alpha = \beta = 1$, $b = 2$, $\varepsilon = 0.1$. About 20 turns of the spiral are shown. Note that the $y$ axis is not an isocline or a trajectory: trajectories can and do cross it. Of course the equation does not make sense for $x < 0$, so what this means is that the prey becomes extinct (and we then are forced to stop harvesting!). In fact, for any initial condition in the first quadrant that is not exactly on the trajectory that enters the saddle point or at the equilibrium point $(x^\star,y^\star)$, the solution will eventually hit $x=0$, perhaps after spiraling a number of times.

A Self-limiting Term

At this point the future of our predator-prey system looks bad: even a very small amount of harvesting causes the eventual extinction of the prey (followed, of course, by the demise of the predator). However, there are other factors that can stabilize the populations. One such factor is a self-limiting term for the prey, so that without predators or harvesting it would undergo logistic rather than exponential growth. Thus we now take

$$x' = x (a - \alpha y - c x) - \varepsilon ,\ \ y' = y(-b + \beta x)$$

where $c$ is a positive constant, assumed to be fairly small.

Again we begin the analysis by finding the equilibrium points. Again, either $y=0$ or $x=b/\beta$ to make $y'=0$. If $y=0$, then $x' = - c x^2 + a x - \varepsilon$, which is a quadratic in $x$. If $\varepsilon$ and $c$ are small, $a^2 - 4 c \varepsilon > 0$ so there are two real roots, and both are positive (their product is $\varepsilon /c$, and their sum is $a/c$). Thus we have two equilibrium points on the positive $x$ axis, $(x_1,0)$ and $(x_2, 0)$, where

$$x_1 = \frac{a - \sqrt{a^2-4c\varepsilon }}{2c}, \qquad x_2 = \frac{a + \sqrt{a^2- 4c\varepsilon }}{2c}$$

On the other hand if $x=b/\beta$ we have $x'=0$ for $y = a/\alpha - b c/(\alpha \beta) - \beta \varepsilon /(\alpha b)$. Thus we have a equilibrium point $(x^\star,y^\star)$, which is somewhat lower than it was without the harvesting and self-limiting terms, but is in the first quadrant if $c$ and $\varepsilon$ are small. The Jacobian matrix is $$J(x,y) = \pmatrix{ a - \alpha y - 2 c x & -\alpha x\cr \beta y & -b + \beta x\cr}$$.

At $(x_1,0)$ the Jacobian matrix is $$\pmatrix{\sqrt{a^2 - 4 c\varepsilon } & -\alpha x_1 \cr 0 & -b + \beta x_1\cr}$$. Since it is upper triangular, its eigenvalues are the diagonal entries $\sqrt{a^2 - 4 c\varepsilon }$ and $-b + \beta x_1$. The first is positive, and (since $x_1$ is small when $c$ and $\epsilon$ are small) the second is negative. So $(x_1,0)$ is a saddle point.

At $(x_2, 0)$ the Jacobian matrix is $$\pmatrix{-\sqrt{a^2 - 4 c\varepsilon } & -\alpha x_2 \cr 0 & -b + \beta x_2\cr}$$. Again the eigenvalues are on the diagonal, and this time the first is negative but (since $x_2$ is large when $c$ is small) the second is positive. So $(x_2, 0)$ is also a saddle point.

For both of these saddle points, an eigenvector for the first eigenvalue is $$\pmatrix{1\cr 0\cr}$$, and again the corresponding trajectories are on the $x$ axis, leaving $(x_1,0)$ and entering $(x_2, 0)$.

At $(x^\star,y^\star)$, the Jacobian matrix is $$\pmatrix{\varepsilon \beta/b-bc/\beta & -\alpha b/\beta\cr a\beta/\alpha - bc/\alpha - \beta^2\varepsilon /(b\alpha) & 0\cr}$$. The characteristic polynomial is

$$r^2 + (b^2 c - \beta^2 \varepsilon ) \frac{r}{b\beta} + ab - \frac{b^2 c}{\beta} - \beta\varepsilon$$

When $\varepsilon$ and $c$ are small, $p^2 - 4 q < 0$ so the eigenvalues are complex. Their real part has the same sign as $\beta^2 \varepsilon - b^2 c$. So, depending on the values of the parameters, we could have a stable spiral or an unstable spiral.

Consider the case $a = \alpha = \beta = 1$, $b = 2$, $\varepsilon = .1$, $c = .125$. This has $\beta^2 \varepsilon - b^2 c = -.4$, so $(x^\star,y^\star)$ should be a stable spiral. In fact, the Jacobian there is $$\pmatrix{-.2 & -2\cr .7 & 0\cr}$$, and its eigenvalues are $-.1 \pm 1.17898 i$. The equilibrium points are $(.101282,0)$, $(7.898718, 0)$ and $(2, .7)$. As we have seen, the trajectories entering and leaving the saddle points are often very important. In this case, it's particularly important to see the relationship between the trajectories in the first quadrant entering the first saddle point $(.101282,0)$ and leaving the second saddle point $(7.898718, 0)$. Both must hit the isocline $x = 2$ at some point above the equilibrium point $(2, .7)$. With the help of a computer (calculating solutions numerically) we can determine that the trajectory entering the first saddle point is above the trajectory leaving the second saddle point. In fact, it appears that the trajectory leaving the second saddle point spirals in toward $(2, .7)$, while the trajectory entering the first must continue along, with $x \to \infty$ and $y \to 0$ as $t \to -\infty$. This trajectory is a separatrix: from any initial point above it, the populations become extinct, while below it they spiral in toward the stable equilibrium $(x^\star,y^\star)$ as $t \to \infty$. So the populations can survive as long as they start out with enough prey and not too many predators. Here is the phase portrait.

If we increase the harvest rate $\epsilon$ to $.3$, keeping the other parameter values the same, something new happens. We still have three equilibrium points $(x_1,0) = (.312182, 0)$, $(x_2,0) = (7.687818,0)$ and $(x^\star, y^\star) = (2, .6)$, and the first two are saddles while the third is a stable spiral. However, when the trajectories entering $(x_1,0)$ and leaving $(x_2, 0)$ cross the isocline $x = 2$, it turns out that the one entering $(x_1,0)$ is lower. The result is that the trajectory leaving $(x_2, 0)$ continues on and hits $x=0$, while the one entering $(x_1,0)$ goes into a spiral. However, note that this one is spiraling outward, while near the stable equilibrium point $(x^\star,y^\star)$ the trajectories must be spiraling inward. How can these be reconciled? The answer is that somewhere between them is a closed curve trajectory, called a limit cycle. Just outside the limit cycle, trajectories spiral outwards, while just inside the limit cycle they spiral inwards. Thus the limit cycle is a separatrix: from any initial point inside the limit cycle, the populations approach the stable equilibrium $(x^\star,y^\star)$ as $t \to \infty$, while from any initial point outside the limit cycle (except those that are exactly on the trajectory that goes in to $(x_1,0)$) they eventually become extinct.

In the picture below, a trajectory spiraling inwards from the limit cycle is shown in red, the trajectory entering the saddle point $(x_1,0)$ in blue, and the trajectory leaving $(x_2, 0)$ in magenta. The limit cycle comes quite close to the saddle point $(x_1,0)$ and the $x$ axis.