📂Dynamics

Lotka-Volterra Predator-Prey Model

Overview

The Lotka-Volterra predator-prey model models the interaction between species as a system, and the predator-prey model in particular represents the predatory relationship between two species. Since dealing with only two species can be extended indefinitely, it is sufficient to represent a food chain.

Model¹

$$ \begin{align*} \dot{x} =& a x - b y \cdot x \\ \dot{y} =& c x \cdot y - d y \end{align*} $$

Variables

• $x(t)$: Represents the number of the prey population $x$ at time $t$.
• $y(t)$: Represents the number of the predator population $y$ at time $t$.

Parameters

• $a>0$: The reproduction rate of $x$.
• $b>0$: The predation rate of $y$ on $x$.
• $c>0$: The reproduction rate of $y$ from preying on $x$.
• $d>0$: The mortality rate of $y$.

Derivation

$$ \dot{N} = rN $$

Let’s start from the Malthusian growth model. First, assume that the prey $x$ grows at a reproduction rate $a>0$ in the absence of predators, and predators $y$ decline at a mortality rate $d>0$ in the absence of prey. Though not truly a system of equations, it can be expressed as follows.

$$ \begin{align*} \dot{x} =& a x \\ \dot{y} =& - d y \end{align*} $$

It is assumed that predation by $y$ on $x$ is proportional to both $x$ and $y$, and it’s attempted to be expressed by some constant $h>0$ as term $hxy$. As $x$ decreases with predation activity and $y$ continues its existence and can reproduce through predation, the equation is modified as follows.

$$ \begin{align*} \dot{x} =& a x - hxy \\ \dot{y} =& hxy - d y \end{align*} $$

It’s important to note that predation activity does not necessarily result in 100% efficiency. In other words, not all consumption directly leads to reproduction, and first and foremost, $x$ and $y$ merely represent numbers and do not accurately reflect ’the amount consumed’. To modify the equation to apply $h$ separately for $x$ and $y$, it can be expressed as follows for some $b,c >0$.

$$ \begin{align*} \dot{x} =& a x - b y \cdot x \\ \dot{y} =& c x \cdot y - d y \end{align*} $$

Here, $by \cdot x$ is intentionally separately articulated as the meaning of $x$’s instantaneous mortality rate being proportional to the number of predators $by$ and $cx \cdot y$ as $y$’s instantaneous growth rate being proportional to the number of prey $cx$ for easier understanding.

■

Customization

In the case where the prey $x$ is a species like plankton or weeds, it might make sense to say that they grow Malthusian in the absence of predators. However, considering that rabbits or sheep can also compete among themselves over grass, Malthusian growth might not be a very realistic assumption.

$$ \dot{N} = {{ r } \over { K }} N ( K - N) $$

Then, simply converting the Malthusian growth to the logistic growth model could improve the model. If $x$ undergoes intraspecific competition, it can be modified for some $K > 0$ as follows.

$$ \begin{align*} \dot{x} =& {{ a } \over { K }} x (K - x) - b y \cdot x \\ \dot{y} =& c x \cdot y - d y \end{align*} $$

Fixed Points

$$ \overline{(x,y)} = \left({{ d } \over { c }} , {{ a } \over { b }}\right) $$

The non-trivial fixed points in the Lotka-Volterra predator-prey system are as above.

Existence

$$ \begin{align*} 0 =& a x - b y \cdot x \\ 0 =& c x \cdot y - d y \end{align*} $$

Let’s find a fixed point that satisfies. Since $(x,y) = (0,0)$ is a trivial solution, assuming $x \ne 0, y \ne 0$, eliminating $x$ from the first equation and $y$ from the second equation,

$$ \begin{align*} 0 =& a - b y \\ 0 =& c x - d \end{align*} $$

is obtained. Rearranging this for $x,y$ yields

$$ \overline{(x,y)} = \left({{ d } \over { c }} , {{ a } \over { b }}\right) $$

■

Stability

Let’s only check for non-trivial fixed points.

$$ \begin{align*} & \begin{bmatrix} a - by & -bx \\ cy & cx - d \end{bmatrix}_{\overline{(x,y)} = \left({{ d } \over { c }} , {{ a } \over { b }}\right)} \\ =& \begin{bmatrix} 0 & -b {{ d } \over { c }} \\ c{{ a } \over { b }} & 0 \end{bmatrix} \end{align*} $$

The Jacobian at $\left({{ d } \over { c }} , {{ a } \over { b }}\right)$ is calculated as above. Its eigenvalues are $$ \begin{align*} \det \left[ \begin{bmatrix} 0 & -{{ bd } \over { c }} \\ {{ ac } \over { b }} & 0 \end{bmatrix} - \lambda I \right] =& \det \begin{bmatrix} -\lambda & -{{ bd } \over { c }} \\ {{ ac } \over { b }} & -\lambda \end{bmatrix} \\ =& \lambda^{2} - \left( - {{ bd } \over { c }} \right)\left( {{ ac } \over { b }} \right) \\ =& \lambda^{2} + ad \\ =& 0 \end{align*} $$ the solutions for $\lambda_{1,2} = \pm i\sqrt{ad}$. Since all eigenvalues are purely imaginary, this non-trivial fixed point is a center.

■

Visual Understanding

$$ \begin{align*} \dot{x} =& {{ 2 } \over { 3 }} x - {{ 4 } \over { 3 }} y \cdot x \\ \dot{y} =& x \cdot y - y \end{align*} $$

Let’s visualize the above system. By randomly drawing 400 initial values and plotting their trajectories, it looks as follows.

It’s clear at a glance that there’s a fixed point $\displaystyle \left( 1, {{ 1 } \over { 2 }} \right)$ and all flows have a period without getting closer or farther.

Looking at the time evolution on each axis, as the predator population decreases, the prey population increases, and as the prey increases, the predator increases, repeats the cycle of the prey decreasing again.

Application

Volterra’s Principle²

Imagine that both predator and prey are harvested or removed by an external intervention by $\varepsilon > 0$. For example, assume both are fish affected by the same net, or both are pests affected by the same pesticide.

Before looking at the equation, one might think that the prey, even accounting for deaths from external intervention, would decrease in number due to the reduction in predators, and predators, already dying from outside, would find their prey reduced. However, the opposite could lead to an increase in prey and a decrease in predators, so intuitively predicting how this system would behave is not straightforward.

Let’s look back at the equation now. External intervention can simply be reflected by multiplying by the population and subtracting it as follows.

$$ \begin{align*} \dot{x} =& a x - b y \cdot x - \varepsilon x \\ \dot{y} =& c x \cdot y - d y - \varepsilon y \end{align*} $$

Then, the newly found fixed points in the changed system are as follows.

$$ \overline{(x,y)} = \left({{ d + \varepsilon } \over { c }} , {{ a - \varepsilon } \over { b }}\right) $$

The analysis of this fixed point is quite interesting, indicating that the number of prey actually increases when both predators and prey are killed equally. This is known as Volterra’s Principle, and, despite the simplicity of the dynamical analysis, it provides a clear insight into external intervention.

Example Code

• Source Code

Code

The following is the Julia code that created the images in the text.

cd(@__DIR__)

@time using Plots
@time using DifferentialEquations

function lotka_volterra_prey_pradator!(du,u,p,t)
    x, y = u
    a, b, c, d = p
    du[1] = dx = a*x - b*y*x
    du[2] = dy = c*x*y - d*y
end
u0 = [0.9, 0.9]
p = [2/3, 4/3, 1., 1.]
end_time = 20.

tspan = (0.,end_time)
prob = ODEProblem(lotka_volterra_prey_pradator!,u0,tspan,p)
sol0 = solve(prob; saveat=(0.0:0.01:end_time))

time_evolution = plot(sol0, vars=(0,1),
  linestyle = :solid,  color = :red,
  size = (400,300), legend=:bottomright, label = "prey");
plot!(time_evolution, sol0,vars=(0,2),
  linestyle = :solid,  color = :blue,
  size = (400,300), legend=:bottomright, label = "predator")
png(time_evolution, "lotka_volterra_prey_pradator.png")

Δt = 50
end_time = 50.
tspan = (0.,end_time)

senario = 400
J = 0.1:0.1:0.5
sol_ = []
for k in 1:senario
    u0 = [1.0rand()*cos(2π*rand())+1.,0.5rand()*sin(2π*rand())+0.5]
    prob = ODEProblem(lotka_volterra_prey_pradator!,u0,tspan,p)
    push!(sol_, solve(prob; saveat=(0.0:0.01:end_time)))
end

anim = @animate for t = 100:10:4950
    print('|')
    frame = plot();
    for j in 1:5
        u0 = [1.,J[j]]
        prob = ODEProblem(lotka_volterra_prey_pradator!,u0,tspan,p)
        frame = plot!(frame, solve(prob; saveat=(0.0:0.01:9)), vars=(1,2),
                        label = :none, color = colorant"#AAAAAA", linestyle = :solid)
    end
    for k in 1:senario
        frame = plot!(frame, sol_[k][t:(t+Δt)], vars=(1,2),
                    label = :none, color = :black, linealpha = (1:(1+Δt))/Δt)
    end
    frame = plot!(frame,
            size = (450,300), xlims = (0,3), ylims = (0,2), aspect_ratio=:equal,
            xlab = "prey", ylab = "predator")
end
gif(anim, "lotka_volterra_preypredator.gif")