One of the most common interactions between species of animals is the predator-prey relationship, which can take many forms. These interactions are critical to the health of animal populations and to the flow of energy through an ecosystem, and altering these populations (by hunting, for example) can have significant consequences to an ecosystem.
The Michigan Department of Natural Resources (DNR) has been examining the population of deer on North Manitou Island (part of Sleeping Bear Dunes National Lakeshore). The main large predators are North Manitou Island are coyotes, and the main large prey animals are white-tailed deer. There are not enough coyotes on the island due to former residents of the island hunting them, and the deer population is getting too large. One consequence of this is that there is not enough readily available food for the deer population, which is destroying the vegetation on the island and also affecting the health of the deer. The DNR has decided that some number of deer should be removed from the population (by hunting or transportation to the mainland) in order for the remaining population to stay at a healthy level.
You have been asked by the Michigan Department of Natural Resources to make a model of the interactions between coyotes and deer to determine how many deer should be removed from North Manitou Island. You're going to do this by developing an agent-based model of the interactions between coyotes and deer, and vary the number and properties of both species to examine how the populations change in relation to each other. In the sections below, we will provide you with some instructions on how to set up the model and some specific questions that the Michigan DNR would like to know the answers to. Happy simulating!
Note: While it's fine to talk with others about this homework - and in fact we strongly encourage you to do so - you need to write your own code and turn in your own homework assignment!
SOLUTIONS
Here are the rules of this model:
Your task is to create a code that implements the predator-prey interaction model described above. No code has been provided to get you started - however, we strongly encourage you to look back at previous in-class and homework assignments, since there's a lot of code from previous assignments that you will find to be very helpful.
Some specific instructions are as follows.
Note that you will be graded both on the correctness of your solution and the quality of your code. Use functions when possible, and make sure that your code is clearly written and has comments explaining what everything does. Also make sure that all plots are clearly marked with axis labels and a title!
In [ ]:
# Put your code here, and add additional cells as necessary
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import random
from IPython.display import display, clear_output
import time
In [ ]:
def initialize_game_board(N_G,N_C,N_D,N_eat):
'''
Initializes game board. Takes in grid size, number of coyotes, number of dear, and how many turns
a coyote can live without eating before it dies. Returns a game board that is N_G * N_G * 2 cells
in size, with the extra slab of cells in the 3rd dimension keeping track of the number of turns a
coyote has left until it dies. Note: for first slab (index 0 in 3rd dimension), 0 = empty, 1 = deer,
2 = coyote
'''
# start out with empty board
game_board = np.zeros((N_G,N_G,2),dtype='int64')
# do a bit of error-checking
if N_C+N_D >= N_G*N_G:
print("board is too small for number of coyotes and deer!")
done = False
placed = 0
# place deer (if we want any)
while N_D > 0 and done==False:
i = random.randint(0,N_G-1)
j = random.randint(0,N_G-1)
if game_board[i,j,0] == 0:
game_board[i,j,0] = 1
placed += 1
if placed >= N_D:
done=True
done = False
placed = 0
# place coyotes (if we want any)
while N_C > 0 and done==False:
i = random.randint(0,N_G-1)
j = random.randint(0,N_G-1)
if game_board[i,j,0] == 0:
game_board[i,j,0] = 2
game_board[i,j,1] = N_eat
placed += 1
if placed >= N_C:
done=True
return game_board
In [ ]:
# take in game board and animal type requested, return total number of that
# type of animal
def number_of_animals(game_board,animal_type):
return (game_board[:,:,0]==animal_type).sum()
In [ ]:
def plot_board(game_board):
'''
This function takes in a 3D array and uses pyplot to make a matplotlib
plot. It returns no value.
'''
# first create two vectors based on the x and y sizes of the grid
x_range = np.linspace(0, game_board.shape[0], game_board.shape[0])
y_range = np.linspace(0, game_board.shape[1], game_board.shape[1])
# use the numpy meshgrid function to create two matrices
# of the same size as game_board with x and y indexes
x_indexes, y_indexes = np.meshgrid(x_range, y_range)
# make a list of all the x and y indexes that are either squares or triangles.
# the notation below is relatively new to us; it means that when game_board==(value),
# only record those values.
sq_x = x_indexes[game_board[:,:,0] == 1];
sq_y = y_indexes[game_board[:,:,0] == 1];
tr_x = x_indexes[game_board[:,:,0] == 2];
tr_y = y_indexes[game_board[:,:,0] == 2];
# plot the squares and triangles. make the size of the polygons
# larger than the default so they're easy to see!
plt.plot(sq_x,sq_y, 'bs',markersize=10)
plt.plot(tr_x,tr_y, '^g',markersize=10)
#Set the x and y limits to include half a space overlap so we don't cut off the shapes
plt.ylim([-0.5,game_board.shape[1]+0.5])
plt.xlim([-0.5,game_board.shape[0]+0.5])
In [ ]:
def move_animals(game_board):
'''
Move every animal on the board to an adjacent empty space. If no space is available,
don't move the animal. Game board is taken in; nothing is returned.
'''
Ngrid = game_board.shape[0]
# loop over the entire board
for i in range(0,Ngrid):
for j in range(0,Ngrid):
# Does this cell have an animal in it?
if game_board[i,j,0]>0:
counter = 0
# loop for a while - we're randomizing, so make it large enough we'll
# plausibly try every adjacent cell.
while counter < 20:
# calculate index to step to
newi = i + random.randint(-1,1)
newj = j + random.randint(-1,1)
# continue if we're trying to step off the board
if newi < 0 or newj < 0 or newi >= Ngrid or newj >= Ngrid:
counter += 1
continue
# if the cell is empty, we can move to it.
if game_board[newi,newj,0]==0:
# move animal to new cell
game_board[newi,newj,0] = game_board[i,j,0]
game_board[newi,newj,1] = game_board[i,j,1]
# set old cell to empty
game_board[i,j,0] = 0
game_board[i,j,1] = 0
# break the while loop
break
else:
counter+=1
In [ ]:
def animal_multiply(game_board,animal_type,prob_mult,N_eat):
'''
Loop through board and check to see if each animal of a given type multiplies.
If it DOES multiply, put the offspring in a random adjacent empty cell. If there
is no adjacent empty cell, don't create any offspring. Only do this for a single
species, because we might have different probabilities of each type of animal multiplying.
Takes in the game board, type of animal, probability of that animal multiplying, and the
number of turns a coyote can go without eating (for initialization of the new board).
Nothing is returned.
'''
Ngrid = game_board.shape[0]
# loop over board
for i in range(0,Ngrid):
for j in range(0,Ngrid):
# only do something if it's an animal of the appropriate type
if game_board[i,j,0]==animal_type:
il = i-1
jl = j-1
ir = i+2
jr = j+2
if il < 0:
il=0
if jl < 0:
jl=0
if ir > Ngrid:
ir = Ngrid
if jr > Ngrid:
jr = Ngrid
# how many of my neighbors are like me? (subtract 1 b/c I'm always like me)
like_me = (game_board[il:ir,jl:jr,0]==game_board[i,j,0]).sum()-1
# if one of my neighbors are like me, check to see if we breed.
if like_me > 0:
# breeding is determined by generating a random number in the interval [0,1)
# and then seeing if it's less than prob_mult.
if random.random() < prob_mult:
# if we breed, put the new animal in a random adjacent cell. If none
# is available, then don't put down a new animal. (This is not exactly
# right.)
counter = 0
#print("new animal",prob_mult,animal_type)
while counter < 20:
# calculate index of new cell
newi = i + random.randint(-1,1)
newj = j + random.randint(-1,1)
# is new cell outside of the grid?
if newi < 0 or newj < 0 or newi >= Ngrid or newj >= Ngrid:
counter += 1
continue
# if new cell is empty, put in our new animal.
if game_board[newi,newj,0]==0:
# clone original animal
game_board[newi,newj,0] = game_board[i,j,0]
# if it's a coyote, set N_eat
if animal_type==2:
game_board[newi,newj,1] = N_eat
#game_board[newi,newj,1]=game_board[i,j,1]
# break the while loop
break
else: # if cell isn't empty, just increment counter
counter+=1
In [ ]:
def eat_or_remove_deer(game_board,P_eat,P_rem,N_eat):
'''
If a coyote and deer are in adjacent cells, check to see if the coyote eats the deer
(only one deer per turn) with probability P_eat. If that happens, the coyote's eating
counter gets reset. There is also some probability P_rem that a deer will get removed
from the island each turn.
'''
Ngrid = game_board.shape[0]
# loop over grid
for i in range(0,Ngrid):
for j in range(0,Ngrid):
# if this cell contains a coyote, check to see if a deer is eaten.
if game_board[i,j,0]==2:
il = i-1
jl = j-1
ir = i+2
jr = j+2
if il < 0:
il=0
if jl < 0:
jl=0
if ir > Ngrid:
ir = Ngrid
if jr > Ngrid:
jr = Ngrid
# check for nearby deer
deer_near_me = (game_board[il:ir,jl:jr,0]==1).sum()-1
# if there ARE nearby deer loop over nearby deer and possibly eat one
if deer_near_me > 0:
#print("there are deer near me!")
have_i_eaten = 0
# loop over cells
for di in range(il,ir):
for dj in range(jl,jr):
# is it a deer?
if game_board[di,dj,0]==1:
# check to see if deer can get eaten and if the coyote CAN eat.
# If so, remove dear, reset N_eat counter, and also set have_i_eaten.
if random.random() < P_eat and have_i_eaten < 1:
game_board[di,dj,0]=0 # remove deer
game_board[i,j,1]=N_eat # reset N_eat counter
have_i_eaten = 1 # set have_i_eaten
#print("I ate a deer!")
# all remaining deer have some probability of being removed by a hunter.
if game_board[i,j,0]==1:
if random.random() < P_rem:
game_board[i,j,0]=0
In [ ]:
def check_remove_coyote(game_board):
'''
Checks to see if a coyote has starved to death, and if so, remove it. Takes in game
board. Returns nothing.
'''
Ngrid = game_board.shape[0]
# loop over grid
for i in range(0,Ngrid):
for j in range(0,Ngrid):
# is this cell a coyote?
if game_board[i,j,0]==2:
# remove one from the counter
game_board[i,j,1] -= 1
# if the counter has gone negative, coyote dies
if game_board[i,j,1] < 0:
#print("coyote died!")
game_board[i,j,0]=game_board[i,j,1]=0
In [ ]:
def evolve_system(N_G = 40, N_S=100, N_C=2, N_D=20, N_eat=10, P_D=0.25, P_C=0.25, P_eat = 1.0, P_rem = 0.0):
#fig = plt.figure(figsize=(N_G,N_G))
# arrays for record-keeping
step = []
N_coyotes = []
N_deer = []
# initialize game board
game_board = initialize_game_board(N_G, N_C, N_D, N_eat)
#plot_board(game_board)
this_step = 0
# loop for N_S steps.
while this_step <= N_S:
# move animals
move_animals(game_board)
# check to see if deer multiply
animal_multiply(game_board,1,P_D,N_eat)
# check to see if coyotes multiply
animal_multiply(game_board,2,P_C,N_eat)
# check to see if deer get eaten or removed
eat_or_remove_deer(game_board,P_eat,P_rem,N_eat)
# check to see if coyotes have died
check_remove_coyote(game_board)
#plot_board(game_board)
#time.sleep(1.0)
#clear_output(wait=True)
#display(fig)
#fig.clear()
# keep track of steps, number of coyotes, number of deer
N_coyotes.append(number_of_animals(game_board,2))
N_deer.append(number_of_animals(game_board,1))
step.append(this_step)
# increment step counter
this_step += 1
# return all of the requested info
return step, N_coyotes, N_deer, game_board
Note: make sure that all plots have appropriate x- and y-limits, as well as figure titles and axis labels. The questions may require both code and a written answer, so please make sure to do both!
For all models, assume a game board that is $N_G = 40$ cells on a side unless otherwise instructed!
Question 1: First, let's test the limits of our model. How does it behave when there are only Coyotes, and only Deer? Try setting $N_C = 20$ and $N_D = 0$, and then do a second model with $N_C = 0$ and $N_D = 20$. For both models, take $N_S = 100$ steps, set $N_{eat} = 10$, and give both deer and coyotes a probability of $P_D = P_C = 0.25$ of reproducing when they interact with others of their species. Set $P_{rem}=0$ (no deer are removed by humans). Use the code you have written to show how these two situations evolve over time (i.e., show the total population of coyotes and deer as a function of simulation time step), explain what you see, and explain why it makes sense.
In [ ]:
# Put your code and figures here - add additional cells if necessary
# N_G = grid size # grid size
# N_S = number of steps
# N_C = number of coyotes
# N_D = number of deer
# N_eat = number of steps before a coyote eats
# P_D = probability that two deer breed and produce a new deer
# P_C = probability that two coyotes breed and produce a new coyote
# P_eat = probability that a coyote eats a deer in an adjacent space
# P_rem = probability that a deer is removed from the island in any given timestep
step, N_coyotes, N_deer, game_board = evolve_system(N_G=20,
N_D=0, N_C=20,
N_S=100,
P_D=0.25, P_C=0.25,
N_eat=10, P_eat=0.5,
P_rem=0.0)
plt.plot(step,N_deer,'b-',step,N_coyotes,'r-')
#print(N_deer,N_coyotes)
In [ ]:
step, N_coyotes, N_deer, game_board = evolve_system(N_G=20,
N_D=20, N_C=0,
N_S=100,
P_D=0.25, P_C=0.25,
N_eat=10, P_eat=0.5,
P_rem=0.0)
plt.plot(step,N_deer,'b-',step,N_coyotes,'r-')
We expect to see exponential behavior in this circumstance - explosive growth of deer and coyote populations. The coyote behavior may be surprising because if there are lots of coyotes they breed before they die, and keep breeding.
Question 2: Now, let's try models with both species, but where they are far out of equilibrium. Consider a model with $N_C = 20$ and $N_D = 2$, and then do a second model with $N_C = 2$ and $N_D = 20$. For both models, take $N_S = 100$ steps, set $N_{eat} = 10$, give both deer and coyotes a probability of $P_D = P_C = 0.25$ of reproducing when they interact with others of their species, and give coyotes a probability $P_{eat} = 0.5$ of catching and eating a deer upon encountering it. Set $P_{rem}=0$ (no deer are removed by humans). Use the code you have written to show how these two situations evolve over time (i.e., show the total population of coyotes and deer as a function of simulation time step), explain what you see, and explain why it makes sense.
In [ ]:
# Put your code and figures here - add additional cells if necessary
step, N_coyotes, N_deer, game_board = evolve_system(N_G=20,
N_D=20, N_C=2,
N_S=100,
P_D=0.25, P_C=0.25,
N_eat=10, P_eat=0.5,
P_rem=0.0)
plt.plot(step,N_deer,'b-',step,N_coyotes,'r-')
In [ ]:
step, N_coyotes, N_deer, game_board = evolve_system(N_G=20,
N_D=2, N_C=20,
N_S=100,
P_D=0.25, P_C=0.25,
N_eat=10, P_eat=0.5,
P_rem=0.0)
plt.plot(step,N_deer,'b-',step,N_coyotes,'r-')
We expect to see exponential growth or decay and possibly some sinusoidal behavior. Out-of-equilibrium situations are fairly stochastic and dependent on initial conditions - be generous with grading.
Question 3: Now, let's modify the model above to try to create stable populations of coyotes and deer by removing deer from the island. Consider a model with $N_C = 2$ and $N_D = 40$. For both models, take $N_S = 100$ steps, set $N_{eat} = 10$, give both deer and coyotes a probability of $P_D = P_C = 0.25$ of reproducing when they interact with others of their species, and give coyotes a probability $P_{eat} = 0.5$ of catching and eating a deer upon encountering it. Start with $P_{rem}=0$ and do several runs where you increase it by steps of 0.1 to $P_{rem} = 1.0$, showing all of the results on a single plot with lines of different colors and/or types indicating different models. How does the situation change as deer are removed by humans? Over what range of $P_{rem}$ do you get populations that are relatively stable, meaning that they do not continuously increase or decrease?
In [ ]:
# Put your code here!
step, N_coyotes, N_deer, game_board = evolve_system(N_G=20,
N_D=20, N_C=2,
N_S=100,
P_D=0.25, P_C=0.25,
N_eat=10, P_eat=0.5,
P_rem=0.1)
plt.plot(step,N_deer,'b-',step,N_coyotes,'r-')
We expect to see: At really low values of P_rem you should get explosive growth of deer population; at higher values you should see the deer population plummet, followed by the coyote population. The difference between explosive deer population and hunting to extinction is quite small in terms of P_Prem - look for that.
Question 4: Run several models where you change the other parameters, and show two or three of the most interesting ones below. Explain in the text box below your figures what is going on, and why.
In [ ]:
# Put your code here!
This could be anything - basically we're looking for some interesting scenarios and analaysis.
Question 5: After you've run all of these different models, what general principles can you take away? In other words, what have you learned about the interactions of predator and prey species from this model?
We're looking for any thoughts that the students have on this subject. basically hoping that they get the idea that unchecked deer populations will grow explosively, and unchecked coyote populations will die off rapidly. Some hunting is necessary to maintain the population at a roughly constant level.
Question 6: In what ways could this model be modified to make it more realistic?
We're looking for pretty much anything sensible here. For example:
There are probably many others; I can't think of them now.
// Write your answer here
// Write your answer here