SWEN90004
Modelling Complex Software Systems
Lecture Cx.03
ODE Models I: Predator-Prey
Artem Polyvyanyy, Nic Geard
artem.polyvyanyy@unimelb.edu.au;
nicholas.geard@unimelb.edu.au
Semester 1, 2021
SLIDE 1
Recap
complex systems
building a mathematical model assumptions
states
update rules
population growth models exponential
logistic
behaviour of dynamic systems fixed points
limit cycles chaos
SLIDE 2
Objectives
interpret mathematical equations describing a model
understand how these equations can be solved numerically and/or analytically to
give the model’s dynamic behaviour
describe the dynamic behaviour of the model, and how this depends on:
model parameters initial condition
SLIDE 3
The Lotka-Volterra model
SLIDE 4
Motivation
Lynxes and hares
Figure90 years of hare and lynx populations in the Hudson Bay, Canada
SLIDE 5
History
Predator-prey relationships
FigureRed fox
The Lotka-Volterra model
FigureRabbit
Independently formulated in 1925–26 by Alfred Lotka, a US mathematician and chemist, and Vito Volterra, an Italian mathematician and physicist. Both were inspired by the logistic equation of Verhulst.
SLIDE 6
Last week: dynamics of a single population (red foxes)
Today: dynamics of both a predator population (red foxes: F) *and* their prey (rabbits: R)
Assumptions
Predators: red foxes — Prey: rabbits
the rabbit population is sustainable (they have ample food)
the rabbit population grows at a rate proportional to its size if there are no red
foxes
the red foxes eat only rabbits
the red fox population declines at a rate proportional to its size if there are no
rabbits to eat
the environment is static and has no effect on the rates of population growth and
decay
SLIDE 7
Model formulation
Prey (rabbits): Predators (red foxes): Parameters:
Parameter
dR dt
dF dt
= αR − βRF
= δRF − γF
Meaning
α growth rate of the rabbit population
β rate at which foxes predate upon (eat) rabbits
δ growth rate of the fox population
γ decay rate of the fox population due to death and migration
SLIDE 8
The growth rate of the rabbit population and the decay rate of the fox population are both proportional to their size.
In the absence of any foxes (F = 0), the rabbit population will grow exponentially.
In the absence of any rabbits (B = 0), the fox population will decay exponentially.
The rate of predation of foxes upon rabbits and the growth rate of the fox population both derive from interactions between the two populations, and hence are dependent on the size of both populations. Note that the rate of rabbit depletion can be different from the rate of fox increase… we could think of this as capturing the relationship between how much ’energy’ a fox gains from eating a single rabbit, vs how much is required to produce a single new fox.
An aside: solving ODEs
SLIDE 9
An aside: solving ODEs using Euler method
Given a set of differential equations and an initial condition, how can we calculate the future states of the system?
In the system
where f defines the relationshp between variables y and their derivatives, the Euler
y′(t) = f(t,y(t))
forward difference approach computes the sequence of approximations:
yn+1 =yn +∆t.f(tn,yn) tn+1 = tn + ∆t
as a time-stepped simulation of the underlying behaviour
SLIDE 10
The simplest approach to solving this type of equation numerically is based on the principle that, from any point on a trajectory, you can find an approximation for a nearby point by taking a small step along a tangent to that trajectory.
However, to achieve good levels of accuracy, this approach needs small time intervals, which rapidly becomes computationally intensive.
An aside: using the mid-point to improve accuracy
Rather that basing our estimate of the next point on the slope at the start of the interval, we can improve accuracy by using the slope at the mid-point:
∆t ∆t yn+1 =yn +∆t.f tn + ,yn + f(tn,yn)
22
tn+1 = tn + ∆t
This approach will diverge from the ‘true’ situation more slowly, but requires more calculations per iteration
SLIDE 11
An aside: solving ODEs using the Runge-Kutta method
The Euler method is often not accurate enough. Accuracy can be increased by making use of more points in the period between tn and tn+1
The Runge-Kutta family of approaches take this approach, with the most commonly used being RK4, which computes:
yn+1 = yn +(
where the kn values are intermediate estimates of slope:
k1 =f(tn,yn)
k2 = f(tn + k3 = f(tn +
∆t ∆t ,yn + 22
∆t ∆t ,yn + 22
k1) k2)
∆t 6
).(k1 +2k2 +2k3 +k4)
k4 =f(tn +∆t,yn +∆t.k3)
SLIDE 12
The kn values correspond to estimates of slope taken at the beginning of the interval between tn and tn+1, at the midpoint (using the first slope), at the midpoint (using the second slope) and at the end (using the midpoint slope).
An aside: solving ODEs using Matlab
Matlab function for the Lotka-Volterra model:
function dx = lotka_volterra_function(t, x) dx = [0; 0];
alpha = 0.75; beta = 0.2; delta = 0.04; gamma = 0.5;
dx(1) = alpha * x(1) – beta * x(1) * x(2);
dx(2) = delta * x(1) * x(2) – gamma * x(2);
return
Matlab test script:
rabbits = 5; foxes = 2;
[t, x] = ode45(@lotka_volterra_function, [0 200], [rabbits foxes], options);
SLIDE 13
Matlab uses very general methods that employ adaptive time stepping — they detect when they need to speed up and slow down to keep the error within reasonable bounds.
the options parameter allows us to specify various details including tolerance—the error we are willing to accept—and constraints—such as the rabbit and fox population levels not becoming negative.
Type help ode45 or doc ode45 for more information.
Behaviour of the Lotka-Volterra model
SLIDE 14
Numerical solution of the Lotka-Volterra model
FigureRed fox (predator) and rabbit (prey) population levels over time
SLIDE 15
Note the lag between the peaks in the predator and prey populations: the rabbit population rises in response to decreased predation by foxes, which in turn supports growth in the fox population. Note also the differences between the dynamics here and what you observed in the ’Wolves Sheep’ NetLogo model (ie, these curves are smooth, compared to the ’noisy’ stochastic curves, and these oscillations are stable, compared to the 2-species stochastic model, which tended to result in one or both species going to extinction.
Numerical solution of the Lotka-Volterra model
FigureThe phase plane
SLIDE 16
This curve represents the trajectory (time series) from the previous slide in a phase plane, in which the value of one of the variables (the size of the fox) population is plotted against the size of the other variable (the size of the rabbit population).
Thus, the size of the two populations at a given moment in time is represented by a single point on this curve.
The trajectories run in a counter-clockwise direction.
Think about where the equilibria (next slide) are in this diagram, and the direction of movement of the system at various different points: on the axes, on the equilibria, near the equilibria, etc.
Analysing the Lotka-Volterra model
What is the long term behaviour of the system?
equilibria: points at which the system variables (the population levels R and F) do not change; ie, where dR =0 and dF =0
R(α − βF ) = 0 F (δR − γ) = 0
Two equilibria:
R = 0, F = 0
R=γ,F=α δβ
dt dt
Interpretation:
equilibrium 1: both populations are extinct
equilibrium 2: both populations sustained at non-zero levels indefinitely
SLIDE 17
Analysing the Lotka-Volterra model
How does system behaviour depend on its initial condition?
Stability analysis reveals that:
the first equilibrium (both populations going extinct) is unstable; it can only arise if
the prey level is set to zero, at which point the predators also die out
the second equilibrium is stable; specifically, it is neutrally stable, and surrounded by infinitely many periodic orbits
How realistic is this?
SLIDE 18
The first equilibrium is what is known as a ’saddle point’ (by comparison with a saddle in a mountain range). A system on the saddle point will remain there. Perturbing the system ’up’ the y-axis (increasing the number of foxes, while leaving the number of rabbits at zero) will cause the system to return to the saddle point (those foxes will die out in the absence of food). Perturbing the system ’right’ along the x-axis (increasing the number of rabbits, while leaving the number of foxes at zero) will cause the system to continue moving right (to infinity – the rabbits will reproduce unchecked). Perturbing the system along both axes simultaneously (adding both foxes and rabbits) will cause the system to move to a new limit cycle (periodic) attractor.
How realistic is this?
The amplitude of the periodic orbit depends on the initial conditions!
In the real world, the environment constantly changes, with the result that the systems location in the phase plane will constantly switch from one orbit to another. This is unrealistic – real cycles tend to have a ‘characteristic’ amplitude.
We would like a model in which the system behaviour stabilised to a single stable limit cycle. This type of system is said to be structurally stable.
Another limitation
FigureAnother Lotka-Volterra phase plane
What is wrong with this picture?
SLIDE 19
We have a fraction of a fox/rabbit at some points on the trajectory. In reality, foxes and rabbits are discrete, rather than continuous, units, and reducing their number below one will result in extinction. Compare to scenario in stochastic Netlogo predation model.
Extensions to the Lotka-Volterra model
SLIDE 20
Lotka-Volterra model with competition
There is no notion of carrying capacity for the prey species
The logistic equation captured intraspecific competition
The base Lotka-Volterra model captues interspecific competition
We can combine these ideas by adding an additional term to the Lotka-Volterra
model to account for competition among members of the same species:
dR = αR − βRF − aR2 dt
dF =δRF−γF−bF2 dt
SLIDE 21
Lotka-Volterra model with competition
FigureTime series FigurePhase plane
Other model refinements have been proposed that increase the ability to capture observed phenomena
SLIDE 22
Other possible issues:
The relation between the predator’s consumption rate and prey density. The efficiency with which extra food is transformed into extra predators. Unlikely to be linear as currently modelled.
The spatial distribution of species will play a role.
Lotka-Volterra model with multiple species
The two-species competitive model can also be written as:
dx1 dt dx2 dt
x1 − α12×2 K1
=r1x1 1−
=r2x2 1−
where αij is the effect that species j has on species i.
x2 − α21×1 K2
This can then be generalised to an n-species model:
dxi Nj=1 αijxj =rixi 1−
dt Ki
where xi represents the i-th species and αij is a matrix of interaction terms.
SLIDE 23
Lotka-Volterra model with multiple species
Another representation of the multi-species model is:
dt
This version pulls the carrying capacities Ki and interaction terms αij inside the interaction matrix A
For three species, growth rates r1 = 1.1; r2 = −0.5; r3 = α + 0.2, and the interaction
matrix:
0.5 0.5 0.1 A = −0.5 −0.1 0.1
α 0.1 0.1
can produce chaotic behaviour for values of α ≥ 1.5.
N dxi=xr+A (1−x)
ii ij j j=1
SLIDE 24
α in this model formulation plays a similar role to r in the logistic map from last week.
Lynxes and hares
Figure90 years of hare and lynx populations in the Hudson Bay, Canada
SLIDE 25
Very few real predator-prey systems exhibit such pure oscillations as the Lotka-Volterra model; one classic data set comes from records kept by the Hudson Bay Company on number of Canadian lynx and snowshoe hare pelts traded.
These exhibited a stable 10-year cycle over almost a century.
Trapping itself obviously impacted on this system (hence the records); however, the numbers taken correspond to only a small fraction of the total population and hence are unlikely to have had a large impact on dynamics.
Summary
SLIDE 26
Recap on dynamic behaviours
Dynamic systems can exhibit different types of long-term behaviour, including fixed points and limit cycles
Fixed points and limit cycles can be:
stable and attracting (often referred to as attractors) unstable and repelling
A single dynamic system can exhibit multiple stable and unstable fixed points
In a system with more than one stable fixed point, the long-term behaviour will
depend on the initial conditions
The set of initial conditions that approach an attractor are known as a basin of attraction
SLIDE 27
Broader applicability of Lotka-Volterra model
In the 1960s, Richard Goodwin, a mathematician and economist, proposed a model anal- ogous to the Lotka-Volterra model for the relationship between wages and employment in an economy:
wages = predators employment = prey
The system operates as follows:
at high levels of employment, employees can bargain for higher wages, which reduces profits
as profits shrink, fewer workers will be employed, which increases profits
once profit levels are higher, more workers will be employed and employment will
increase, leading to cyclic behaviour.
Complex systems often exhibit underlying patterns of behaviour across multiple domains
SLIDE 28
Whether these associations are meaningful is another question, but insights gained from the study of one system can certainly help stimulate thought about the dynamic mechanisms that operate in other systems.