Importing packages, specifying some plotting “constants”
# Intro packages + define quantities
using Plots; pyplot()
using LaTeXStrings
using DifferentialEquations
#
LW1 = 2.5
LW2 = 1.5
LS1 = :solid
LS2 = :dot
LC1 = :red
LC2 = :blue;
Introducing LW1, etc. is to enforce consistent linewidths, colors, line styles, colors, etc., and make it easy to change these things in one location.
Population model
Comment on notation: in the description below, N is population number within some system boundary, \frac{dN}{dt} is the time derivative of N, while \dot{N} is a rate of N – either (i) the rate of N crossing the system boundary, or (ii) a generation rate within the system boundary. This notation is similar to what is common in some engineering disciplines. Specifically, \dot{N} is not the time derivative of N. In other words, I do not use Newton’s notation for (temporal) derivative.
Consider the population balance
\frac{dN_j}{dt} = \dot{N}_{\mathrm{i},j} - \dot{N}_{\mathrm{e},j} + \dot{N}_{\mathrm{g},j} - \dot{N}_{\mathrm{d},j} + \dot{N}_{\mathrm{m},j}
where N_j is the accumulated value (“population”) of species j, \dot{N}_{\mathrm{i},j} is the influent rate/immigration rate of species j, \dot{N}_{\mathrm{e},j} is the effluent rate/emigration rate of species j, \dot{N}_{\mathrm{g},j} is the growth rate of species j, \dot{N}_{\mathrm{d},j} is the death rate, and \dot{N}_{\mathrm{m},j} is the population management rate.
In the standard Lotka-Volterra model, an unmanaged “island” population is assumed, i.e., \dot{N}_{\mathrm{i},j}\equiv 0, \dot{N}_{\mathrm{e},j}\equiv 0, and \dot{N}_{\mathrm{m},j}\equiv 0. For the predator, growth is based on feeding off the prey, \dot{N}_{\mathrm{g},1} = \gamma_1 N_1 N_2, and death results from natural death/starvation \dot{N}_{\mathrm{d},1} = \delta_1 N_1. For the prey, growth without limiting factors (e.g., unlimited availability of grass) is assumed, \dot{N}_{\mathrm{g},2} = \gamma_2 N_2, and death comes from being consumed by the predator, \dot{N}_{\mathrm{d},1} = \delta_2 N_1 N_2. If prey growth is limited due to resource constraints, a modification is \dot{N}_{\mathrm{g},2} = \gamma_2 N_2\cdot \left( 1-N_2/N_{2,\max} \right); the logistic model.
Some standard parameters for the Lotka-Volterra model are, e.g., \gamma_1 = 0.1, \delta_1 = 0.4, \gamma_2 = 1.1, and \delta_2 = 0.4 (Wikipedia; cheetah and baboon). Typical initial conditions can be N_j(0) = 10.
The DifferentialEquations population model
#
# Population model
# -- population balance
function Nmod(dN,N,p,t)
dN .= Ndi(t,N,p) .- Nde(N,p) .+ Ndg(N,p) .- Ndd(N,p) .+ Ndm(t,N,p)
end
# -- immigration rate
function Ndi(t,N,p)
return zeros(size(N))
end
# -- emigration rate
function Nde(N,p)
return zeros(size(N))
end
# -- growth rate
function Ndg(N,p)
if p.limited_growth == false
return [p.g1*N[1]*N[2], p.g2*N[2]]
else
return [p.g1*N[1]*N[2], p.g2*N[2]*(1-N[2]/p.N2max)]
end
end
# -- death rate
function Ndd(N,p)
return [p.d1*N[1], p.d2*N[1]*N[2]]
end
# -- management rate
function Ndm(t,N,p)
if p.management_policy == 0
# zero management
return [0,0]
elseif p.management_policy == 1
# temporal management -- inserting prey for t>= 60
if t < 60
return zeros(size(N))
else
return [0,1]
end
else
# state based management
return [-0.5*max((N[1]-5),0),0]
end
end
;
Case: unlimited growth, no management
# Model parameters -- named tuple
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=false,management_policy=0)
# Initial population
N0 = [10.,10.]
#
tspan=(0.,100.)
#
prob = ODEProblem(Nmod,N0,tspan,p);
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="No management, unlimited prey growth")
The result is as follows:
Case: limited growth, no management
In a more realistic model, the growth rate of the prey is limited by the carrying capacity of the area — e.g., how much grass that is available before the rabbit population starts to eat all of it. The logistic model describes this. Let us set the capacity to N_{2,\max} = 400.
# -- modified growth rate
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=true,management_policy=0)
prob = ODEProblem(Nmod,N0,tspan,p)
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="No management, limited prey growth")
with result:
Case: unlimited growth, temporal management (sets in at t\ge 60)
# -- temporal management
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=false,management_policy=1)
prob = ODEProblem(Nmod,N0,tspan,p);
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="Temporal management")
Case: unlimited growth, state based management
Here, the simple policy is used of hunting the predator if its population exceeds 5.
# -- state based management
p = (g1=0.1,g2=1.1,d1=0.4,d2=0.4,N2max=400,limited_growth=false,management_policy=2)
prob = ODEProblem(Nmod,N0,tspan,p);
#
sol = solve(prob)
plot(sol,vars=(0,1),lc=LC1,lw=LW2,label="Predator")
plot!(sol,vars=(0,2),lc=LC2,lw=LW2,label="Prey")
plot!(title="State based management")
with result:
Remark: Here, the population number N is posed as an integer. The time derivative is meaningless for integers. An “engineering” approach to this is to assume that N is “large” and pretend that N “behaves” as a real number. This is hardly true for N in the order of 10 as above… and the plots show solution values that are non-integers. One rule of thumb seen in the literature is that N should be larger than 10^4 for this approximation to be valid. The requirement for the size of N for this approximation to be true (e.g., pretending that the integer N behaves like a real number) is probably problem dependent. The same difficulty shows up when considering molecules as particles and posing the “number” balance. In many problems, N is so large that this is hardly a problem (i.e., in the order of Avogadro’s number).
