Predator-Prey#
Author:
Date:
Time spent on this assignment:
Remember to execute this cell to load numpy and pylab.
0. Background#
In this assignment we are going to consider the evolution of predators and preys within an ecosystem.
We will focus on a single population of prey \(B\) (bunnies) and a single population of predators \(F\) (foxes).
Three possible things can happen:
Prey reproduce. Each prey produces \(\alpha \Delta t\) new prey in the unit of time \(\Delta t\). Therefore, the total number of new prey generated at time \(t\) depends linearly on how many prey \(B(t)\) there are with \(\Delta B(t) = \alpha B(t) \Delta t\). We can write this as a differential equation
\(\frac{dB(t)}{dt} = \alpha B(t)\)
Predators die. Like reproduction, the total number of predators that die at time \(t\) depends linearly on the total number of predators \(F(t)\) such that the change \(\Delta F(t) = -\gamma F(t) \Delta t\). We can write this as a differential equation
\(\frac{dF(t)}{dt} \propto -\gamma F(t)\)
Predators eat prey.
The rate at which predators and prey run into each other is going to depend multiplicatively as \(F(t)B(t)\). One way to think about this is if I scatter predators and prey randomly over a grid, the probability that they will be at the same place goes as \(\propto F(t)B(t)\).
If a prey gets eaten, then it’s gone so the change in prey is \(\Delta B(t) = -\beta F(t)B(t)\).
If a predator eats a prey then it is no longer hungry and can reproduce and so \(\Delta F(t) = \delta F(t) B(t)\).
We can write these as the differential equations\(\frac{dB(t)}{dt} = -\beta F(t)B(t)\)
\(\frac{dF(t)}{dt} = \delta F(t)B(t)\)
Combining 1-3 gives us two coupled equations:
These are called the Lotka-Volterra equations
Exercise 1: Solving the predator-prey differential equations#
List of collaborators:
References you used in developing your code:
a. Euler integration of the Lotka-Volterra equations#
In this exercise we will solve the differential equations which correspond to the predator-prey equations.
🦉 Write a function that uses Euler integration to solve these differential equations. You can use any structure you want, but one possibility is:
@dataclass
class Volterra_parameters:
alpha: float = 1.0 # reproduction rate of bunnies
beta: float = 0.005 # rate at which foxes eat bunnies
delta: float = 0.003 # rate at which foxes eating bunnies results in more foxes
gamma: float = 0.6 # death rate of foxes
and
def run_predator_prey(p, dt, T, B0, F0):
"""
Parameters
----------
p: Volterra_parameters dataclass object
dt: timestep
T: total time
B0: Initial bunny population
F0: Initial fox population
Returns
-------
ts, B, F : np.ndarray (T//dt,) times, and population of bunnies and foxes over time.
"""
### ANSWER ME
b. Building Plots#
Let’s set \(\alpha=1.0\), \(\beta=0.005\), \(\gamma=0.6\) and \(\delta=0.003\). Use a timestep \(dt=0.001\)
We are going to start with \(400\) prey (bunnies) and \(200\) predators (foxes).
Set up your code to run for time \(T=50\)
🦉Plot two graphs
Population vs time
prey (Bunnies) vs time.
predator (Foxes) vs time.
predator vs prey (a phase space plot)
### ANSWER ME
The dynamics for the parameters above (predator vs prey) result in what is called a limit cycle. This means that if I start with a given number of predators and prey, the total number cycles in this periodic pattern forever. You can see this either by seeing periodicity in the prey vs. time plots, or by noting that the phase space plot has a repeating orbit. Make sure to answer all these questions in full sentences with your reasoning.
Q: For the parameters above, is the first derivative of bunnies or foxes positive (which initially increases)?
A:
Q: Suppose you were to reverse the initial conditions, so that there are more predators than prey at the beginning. How would the answer to the above change (you can try it!)?
A:
Q: What is the relationship between the maximum in bunnies vs foxes? In phase, out of phase by 180 degrees, something else?
A:
c. Predator Phase Space Plots#
🦉Generate the phase space plots for a \(x\) prey and \(x/2\) predators for \(x \in \{50,100,200,400,600,800,1000,1200,1600\}\). (Use the same parameters as above). Also generate a plot of the minimum prey seen during an ecosystem simulation for a given number of starting prey.
### ANSWER ME
Q: What is the minimum number of bunnies you ever see in your simulation?
A:
Q: As you increase \(x\), you will find that the populations start to crash near zero. Which population does that first?
A:
Q: Which value of \(x\) do you think is most likely to avoid extinction given these parameters ?
A:
d. Stationary points#
An interesting question we can ask: is there is some number of prey and predators we can start with so that the total number of prey and predators never change? This would be called a fixed point or stationary point.
We will determine the stationary point numerically. 🦉 To accomplish this, implement the following function:
def get_deviation(initial_population):
"""
Run the predator-prey simulation from the given initial conditions and
return the standard deviation of the prey population over time.
A standard deviation of zero means the prey population never changes,
i.e. the system is at a stationary (fixed) point.
Parameters
----------
initial_population : list
[preyPopulation, predatorPopulation]
Returns
-------
float
np.std of the prey-vs-time list produced by the simulation.
"""
Now, all we have to do is call get_deviation a bunch of times with different parameters and figure out when it returns zero. There are lots of rich and interesting approaches to optimization problems like this, but we don’t have time to explore them. Here we are going to use a package from python calling out=scipy.optimize.minimize(get_deviation,x0). When it’s finished running, if you print out.x you can see which parameters it found.
Let’s set \(\alpha=1.0\), \(\beta=0.005\), \(\gamma=0.6\) and \(\delta=0.001\). (note the different parameters)
🦉Run the optimization and report the number of predator and prey you find.
Then make a plot of
prey vs. time and predator vs. time for these initial conditions
a phase plot of predator vs prey using these initial conditions .
You also might find it useful to turn off matplotlib’s formatting as follows: matplotlib.rc('axes.formatter', useoffset=False) so you can see more clearly what numbers are changing.
### ANSWER ME
### ANSWER ME
Q: Qualitatively what do phase plots look like when a stationary point is found?
A:
Exercise 2: Stochasticity in Predator-Prey equations#
List of collaborators:
References you used in developing your code:
Problems with simulating the ecosystem using the Lotka-Volterra equations:
There are two qualitative problems with simulating ecosystems using the Lotka-Volterra equations:
(1) In our simulations so far, we have found either stationary points or limit cycles. Both of these are unrealistic. Ecosystems are more dynamic than this. It’s not the case that once you have a certain population of predator/prey it returns to this number over and over again forever. This occurs because we are simulating the ecosystem deterministically - there is no randomness involved. Deterministic equations essentially can only do a few things: they can converge to fixed points, limit cycles or be chaotic.
(2) millibunnies: You may have noticed in your simulations that you occasionally saw that your ecosystem had 1/1000 of a bunny and that this 1/1000 of a bunny recovered into a respectable population of 8000 bunnies. This is unrealistic, since there aren’t actually fractional animals (and also it prevents extinction events). To go beyond this we need to recognize the discrete integer-ness of animals, which will require us to also include stochasticity.
a. The first bunny#
Let’s think about what \(dB = \alpha B(t) dt\) means. Previously we thought that it meant that in a unit of time \(\Delta t\) we got \(\alpha B(t) \Delta t\) new bunnies. But if \(\alpha B(t)\) is 0.2 and \(\Delta t=0.1\) we don’t get 0.02 bunnies. Instead we can think of it as the probability density of getting one whole new bunny as \(20\%\) and so the probability of getting a new bunny in a small interval \(\Delta t\) is \(0.2 \Delta t\).
We now ask the following question: What is the probability of getting your first bunny after time \(T+\Delta t\)? Well, we’ve taken \(T/\Delta t + 1\) steps. You have to not have gotten a bunny the first \(T/\Delta t\) steps. That probability is \((1-0.2 \Delta t)^{T/\Delta t}\). The probability of getting a bunny at step \(T/\Delta t+1\) is \(0.2 \Delta t\). So the total probability density is \((1-0.2 \Delta t)^{(T/\Delta t)} \times 0.2\)
Suppose we take the limit where \(\Delta t \rightarrow 0\) The the probability we get a bunny in an interval \(\Delta t\) is
This follows because
This is called a Poisson process, and it shows up in many situations, from nuclear decay to phone call frequency. It is essentially saying that a bunny reproduces on average every \(1/\alpha\) units of time, and once it has reproduced, it starts again with no memory of when it last reproduced. You can imagine that such a process might not apply to human behavior but we will use it for bunnies.
Let’s test this.
🦉Implement the following function:
def time_of_first_bunny(alpha, deltaT):
"""
Simulate a Poisson process by stepping forward in increments of deltaT.
At each step, a new bunny appears with probability alpha * deltaT.
Parameters
----------
alpha : float
Reproduction rate (bunnies per unit time).
deltaT : float
Size of each time step.
Returns
-------
float
The time at which the first whole new bunny appears.
"""
🦉Run this function 10,000 times with \(\alpha=1\) and \(dt=0.001\) and plot it with plt.hist(ts,100,density=1) (100 sets the number of bins and density=1 makes sure it is a normalized histogram).
On the same plot, plot
from \(T\in[0,35]\) to compare them. (Recall that if you have a numpy array you can do math operations on them using functions like np.exp)
### ANSWER ME
b. A faster algorithm to find the first bunny.#
Waiting to determine the time of first of bunny using the method of a. above is pretty slow. There is another algorithm that’s much faster (in the limit of \(\Delta t\rightarrow 0\)). The time to first bunny is
where \(r\) is a uniform random number between 0 and 1. The idea here is that we are drawing a sample from the Poisson distribution. This method is sometimes called “continuous time” because \(t\) can be any real number, and is not decided by a discrete time step.
Instead of analytically verifying this, let’s check this numerically.
🦉Implement the following function and make a plot comparing it to \(\alpha \exp[-\alpha T]\):
def time_of_first_bunny2(alpha):
"""
Draw the time to the first bunny directly from the Poisson distribution
using the continuous-time formula t = -log(r) / alpha.
Parameters
----------
alpha : float
Reproduction rate (bunnies per unit time).
Returns
-------
float
A single sample of the waiting time until the first bunny appears.
"""
### ANSWER ME
c. Continuous time markov chain#
We will now use this to build a stochastic algorithm for predator-prey. Let’s translate everything we did above into this new continuous time language:
Prey reproduce at rate \(\alpha B(t)\)
\(\frac{dB(t)}{dt} = \alpha B(t)\)
At time \(t_1=-\log(r)/(\alpha B(t))\), the prey \(B \rightarrow B+1\)
Predators die at rate \(\gamma F(t)\)
\(\frac{dF(t)}{dt} \propto -\gamma F(t)\)
At time \(t_2=-\log(r)/(\gamma F(t))\), the predator \(F \rightarrow F-1\)
Predators eat prey at rate \(\beta F(t)B(t)\) which cause predators to reproduce at \(\delta F(t)B(t)\)
\(\frac{dB(t)}{dt} = -\beta F(t)B(t)\)
\(\frac{dF(t)}{dt} = \delta F(t)B(t)\)
At time \(t_{3a}=-\log(r)/(\beta F(t)B(t))\) the prey \(B \rightarrow B-1\)
At time \(t_{3b}=-\log(r)/(\delta F(t)B(t))\) the predator \(F \rightarrow F+1\).
This latter point is a bit weird but works. Both the fox eating the bunny, and the reproduction of new foxes is proportional to how much the fox eats, so these two processes in practice can be separated in this way.
At each step:
Compute times \(t_1, t_2, t_{3a}, t_{3b}\) using a different random number \(r\) for each.
Find the smallest \(t_i\) — that event happens first.
Apply event\(_i\) (e.g. add a bunny, remove a fox) and advance \(t \rightarrow t + t_i\).
Discard the other times (rates have changed). This is valid because Poisson processes are memoryless.
Stop if either population reaches 0.
🦉Implement the following function:
def run_ecosystem(p, T, B0, F0):
"""
Run the stochastic (continuous-time Markov chain) predator-prey simulation.
Parameters
----------
p : Volterra_parameters object
Parameters for the evolution.
T : float
Total simulation time.
B0 : int
Initial prey (bunny) population.
F0 : int
Initial predator (fox) population.
Returns
-------
ts : list of float
Times at which events occurred (not equally spaced).
Bs : list of int
Prey population immediately after each event.
Fs : list of int
Predator population immediately after each event.
"""
### ANSWER ME
d. Stochastic Populations#
🦉Plot
Population vs t
prey vs t
predators vs t
Phase plot (predator vs prey)
for \(T=50\), \(B(t=0)=100\), \(F(t=0)=50\), and \(\alpha\)=1.0, \(\beta\)=0.005, \(\gamma=0.6\), and \(\delta=0.003\)
On the same plots, show the results for the deterministic algorithms. Run it a few times to see the effect of stochasticity.
### ANSWER ME
e. Stationary again#
🦉Do the same thing as d. for the parameters you found above which give you a stationary point.
Remember to run using \(\alpha\)=1.0, \(\beta\)=0.005, \(\gamma=0.6\), and \(\delta=0.001\)
For this, go out to \(T=1000\)
### ANSWER ME
Q: Do the number of bunnies change as a function of time?
A:
Q: Qualitatively what is the effect of running the stochastic simulation at the stationary point?
A:
f. Extinction#
In this part we are going to determine the amount of time until one of the two species goes extinct given the initial number of prey/predators. Let’s call this an extinction event.
Use parameters
\(\alpha=1.0\)
\(\beta=0.008\)
\(\gamma=0.3\)
\(\delta=0.002\)
and assume we start with the same number of predators and prey.
🦉For \(B_0\in\) preyNum w/ preyNum=np.arange(20,600,1), run your ecosystem until there is one extinction event. Record the time \(t\) this event happened and graph \(t\) vs. initial prey number.
np.arange(20,60,1) first.
### ANSWER ME
g. Window Averaging (cleaning data)#
This data is going to look noisy. 🦉To make it look nicer do a running mean of your data with a window-size of 20 and plot it. This will be smoother.
What value is your largest time to extinction (post window-averaging)? I find my results are around 160-190
### ANSWER ME
Q: I would have expected that this curve would be more monotonic. The more prey you had, the less likely anyone would go extinct. Is this correct? Do a calculation with B0=600 and describe why the bunnies go extinct quickly.
A:
Exercise 3: Agent Based Predator-Prey Simulation (extra credit)#
List of collaborators:
References you used in developing your code:
In this exercise, we want to simulate an agent-based model of the predator-prey simulation. Instead of the previous exercises where you are using a differential equation, here we are going to simulate each individual bunny and fox separately as an agent. These bunnies and foxes are going to live on a \(s \times s\) grid. Foxes will keep track of how full they are. All foxes are born with a fullness amount \(F\). The dynamics of our simulation will be:
The bunnies and foxes will be allowed to move around the grid. At each time step, every bunny and every fox will decide either to stay where they are or try to move in one of four directions. If they decide to try to move, they will move if its empty of the same type of animal (i.e. foxes don’t move to where foxes are). If the grid space they wanted to go is occupied, they will stay where they are.
Every time step, every fox’s fullness decreases by \(f\).
When the foxes and bunnies are on the same grid point (after everyone moves) one fox on that grid (pick arbitrarily) eats all the bunnies on that grid. The foxes full-ness increases by the number of bunnies they eat.
A fraction \(\alpha\) fraction of bunnies reproduce (appearing at the same location) at each time-step.
All foxes whose fullness \(F\) is more than \(R\), reproduce (at the same location) at each time-step.
All foxes whose fullness \(F\) is less than 0 starve.
To do this, we are going to have a class to store each “fox” or “bunny”:
@dataclass
class fox:
location: list
full: int = 2
and
@dataclass
class bunny:
location: list
eaten: bool = False
The fox starts full with an amount 2 and the bunny starts not eaten. There is no default location for them.
To set up a fox you can do myFox=fox() and to access the location of the fox you can do myFox.location. If you want to push your fox back onto a list you can do myList.append(fox()).
To implement this project, you will use four data structures:
a list of foxes
a list of bunnies
a \(s\times s\) numpy array foxLocations that shows how many foxes are in each grid
a \(s\times s\) numpy array bunnyLocations that shows how many bunnies are in each grid.
To set things up, first set up the grid with a density of 10% bunnies and 1% foxes.
foxLocations=np.random.choice(2, size*size, p=[.99,0.01]).reshape(size,size)
bunnyLocations=np.random.choice(2, size*size, p=[.9,0.1]).reshape(size,size)
🦉Implement the following four functions:
def locations_to_list(locations):
"""
Convert a 2-D location array into a list of animal objects whose
.location attribute is set to the corresponding grid coordinate.
Parameters
----------
locations : np.ndarray, shape (s, s)
Grid where each non-zero entry indicates an animal at that cell.
Returns
-------
list
List of fox (or bunny) objects with .location = [row, col]
for every occupied cell.
"""
def get_locations(animalList):
"""
Convert a list of animal objects back into a 2-D location array.
Parameters
----------
animalList : list
List of fox or bunny objects, each with a .location attribute.
Returns
-------
np.ndarray, shape (s, s)
Grid where entry [i, j] equals the number of animals at cell (i, j).
"""
You can verify these are inverses by checking that
locations_to_list(get_locations(locations_to_list(myLocations))) returns the same animals you started with.
def move(animalList):
"""
Move every animal in animalList one step on the grid.
For each animal, flip a five-sided die: stay put, or attempt to move
north, south, east, or west. The animal moves only if the target cell
is unoccupied by the same species; otherwise it stays put.
Parameters
----------
animalList : list
List of fox or bunny objects with a .location attribute.
Returns
-------
animalList : list
Copy of the fox or bunny object list after moving
"""
Tip: After the function returns, np.sum(animalLocations) == len(animalList), np.max(animalLocations) == 1, and np.min(animalLocations) == 0 should all hold.
def eat(bunnies, foxes):
"""
Resolve predation: for every cell that contains both foxes and bunnies,
one fox eats all bunnies on that cell.
Creates a new list in which:
- eaten bunnies are removed from the grid.
- bunny.eaten: set to True for each eaten bunny.
- fox.full: increased by the number of bunnies eaten.
Eaten bunnies are removed from the bunnies list before the function returns.
Parameters
----------
bunnies : list
List of bunny objects
foxes : list
List of fox objects .
Returns
----------
bunnies : list
new list of bunny objects
foxes : list
new list of fox objects
"""
Finally, implement two reproduce functions separately for foxes and bunnies. The fox function should do the following:
def fox_reproduction(foxes, f, R):
"""
* The foxes' full state should decrease by $f$.
* If the foxes' full state is <=0 then the fox starves
* If the foxes' full state is >=R then the fox reproduces.
"""
The bunny function should do the following:
def bunny_reproduction(bunnies, alpha):
"""
* an alpha'th fraction of the bunnies should reproduce.
"""
Have them both return new lists.
Use
\(\alpha=0.01\) (reproduction rate of the bunnies)
\(R=5\) (reproduction fullness for the foxees)
\(s=80\) (size of the grid)
\(f=0.1\) (amount fullness decreases for each fox)
Run your code and plot both
bunnies and foxes vs. time
a phase plot of foxes vs. bunnies
an animation (see the code below). To do the animation you need to keep snapshots of the fox-locations and bunny-locations at each step.
You should see features reminiscent of your previous work.
### ANSWER ME
### ANSWER ME
Acknowledgements:
Bryan Clark (original)
© Copyright 2020
–