We will examine the effect of tree density, tree health, and biodiversity on the spread-rate of the Emerald Ash Borer.
Given the variable behavior and life stages of the Emerald Ash Borer and the importance of spatial considerations when modeling their behavior, an ABM methodology is necessary. A Differential Equations approach would not be prudent in this case, given the relatively stochastic ways in which the Ash Borer is introduced and subsequently spreads.
In this model, I am primarily interested in the effects of tree density, tree health, and biodiversity on the spread-rate of the Emerald Ash Borer. To measure these effects, this model considers the length of time it takes for the entire tree population in the study area to be infected and killed by the Emerald Ash Borer.
In this model, the space will be a 2D square grid without wrapping edges. Each grid cell will contain zero or one tree.
Beetles are the main agent in this model which will spread outward from a randomly chosen focal point and degrade tree health to the point that the tree is killed. This model simulates their life cycle by stepping them through their life stages from egg to larvae to adult to dormant. Because reproduction is my main concern, I will only model female beetles.
For their step function, beetles actions will vary depending on life stage:
Trees are the other agent in this model. They will grow and die according to their interactions with beetles.
For their step function, the following will update:
At this time, this model has no institutions.
Based on the above description, we need the following:
In order to answer my question, I'll sweep tree density, tree health, and biodiversity (tree_population, min_/max_tree_heath, and prop_ash, respectively) to see how the Ash Borer spread rate is affected. Tree health will affect the length of time it takes for an ash borer infestation to kill a tree, and tree density/the proportion of ash trees will affect the ease with which the adult ash borers can find a suitable host. The model run for each set of parameters will stop once all ash trees in the study area (grid) have been killed. Then the number of timesteps that occurred will be counted and divided by four to quantify the number of years it took for the ash borer to decimate the ash population. By comparing these numbers, we can begin to see which factors affect Ash Borer spread rate, which could lead to management recommendations.
In [2]:
%matplotlib inline
# Standard imports
import copy
import itertools
# Scientific computing imports
import numpy
import matplotlib.pyplot as plt
import networkx
import pandas
import seaborn; seaborn.set()
# Import widget methods
from IPython.html.widgets import *
Below, the class constructor initializes the tree when we call Tree(). Also, we have the following methods:
In [ ]:
class Tree(object):
"""
Tree Class.
"""
def __init__():
"""
constructor for tree class. By default, all trees:
* are not ash
* have random dbh between specified min and max (basic uniform dist)
* have random health between specified min and max (basic uniform dist)
* have no beetles of any lifestage
"""
self.model = model
self.tree_id = tree_id
self.tree_species = tree_species
self.tree_dbh = tree_dbh
self.tree_health = tree_health
self.has_beetles = has_beetles
self.num_eggs = num_eggs
self.num_larvae = num_larvae
self.num_dormant = num_dormant
self.num_adult = num_adult
#because trees cannot move, we'll simplify later modeling by hard coding
#tree location into their properties
self.xloc = xloc
self.yloc = yloc
def tree_growth(self):
"""
Evaluates whether a tree will grow or not, given the presence
of beetles. If beetles are not present, tree dbh will increase
by a randomly generated number within a range.
"""
if self.has_beetles == "Yes":
self.tree_dbh += 0
else:
self.tree_dbh += numpy.random.uniform(self.min_growth,
self.max_growth)
def is_ash(self):
if self.tree_species == "ash":
return True
if self.tree_species == "nonash":
return False
def evaluate_larvae_num(self):
"""
This method correlates dbh with how many larvae a tree is able to host
"""
def is_infected(self):
if self.num_adults > 0 or self.num_dormant > 0 or self.num_larvae >0 or self.num_eggs >0:
return True
else:
return False
Below, the beetle class constructor creates a beetle every time we call Beetle(). It also houses the following methods:
In [ ]:
class Beetle(object):
"""
Beetle Class
"""
def ___init___():
"""
constructor for beetle class. By default, there is only one
adult beetle located on a randomly chosen ash tree.
"""
self.lifestage = lifestage
self.beetle_id = beetle_id
self.num_eggs_laid = num_eggs
self.egg_survival = egg_survival
self.winter_survival = winter_survival
self.wanderlust = wanderlust
def get_location(self):
"""
returns beetle position, calling through model
"""
return self.EABModel.get_beetle_location(self.beetle_id)
def egg_to_larva(self):
"""
determines whether an egg will become a larva
"""
if numpy.random.random() <= self.egg_survival:
return True
else:
return False
def dormant_to_adult(self):
"""
determines whether a dormant beetle will emerge
as an adult
"""
if numpy.random.random() <= self.winter_survive:
return True
else:
return False
Below is the constructor (init method), which initialized my model when I call EABModel(). It also has the following methods:
In [ ]:
class EABModel(object):
"""
Model Class.
"""
def __init(self,
initial_beetles = 1,
initial_lifestage = 3
tree_population = 40,
ash_prop = 0.2,
min_dbh = 12,
max_dbh = 25,
min_health = .75,
max_health = 1,
min_growth = .03 #inches per timestep
max_growth = .07 #inches per timestep
grid_size = 55,
tree_health = 0.9,
min_egg_survival = .75
max_egg_survival = .99
min_winter_survival = .7
max_winter_survival = .99
):
self.initial_beetles = initial_beetles
self.initial_lifestage = initial_lifestage
self.tree_population = tree_population
self.ash_prop = ash_prop
self.min_dbh = min_dbh
self.max_dbh = max_dbh
self.grid_size = grid_size
self.min_health = min_health
self.max_health = max_health
self.min_growth = min_growth
self.max_growth = max_growth
self.min_egg_survival = min_egg_survival
self.max_egg_survival = max_egg_survival
self.min_winter_survival = min_winter_survival
self.max_winter_survival = max_winter_survival
self.agent_list = []
# call setup methods to initialize model starting point.
self.setup_space()
self.create_agents()
self.infest_initial_trees()
def setup_space(self):
"""
creates an empty grid of tuples according to our specified grid_size
"""
self.space = numpy.full((self.grid_size, self.grid_size), numpy.nan)
"""
QUESTION: is it possible to create a grid of lists? I need to keep track of both tree_id
and various beetle_ids within each cell, and a given tree could have upwards of 300 beetles
over all possible lifestages at any given time
"""
def create_agents(self):
"""
create the agents for the model
"""
# create trees with default values, specified above
for i in range(self.tree_population):
self.agent_list.append(Tree(model = self, tree_id = i,
tree_species == "nonash",
tree_dbh = numpy.random.uniform(self.min_dbh, self.max_dbh),
tree_health = numpy.random.uniform(self.min_health, self.max_health),,
has_beetles = False,
num_eggs = 0,
num_larvae = 0,
num_dormant = 0,
num_adult = 0
xloc = 0
yloc = 0))
# randomly endow empty grid with trees, giving them a species
# based on our specified ash proportion.
for tree in self.agent_list:
#loop until unique
occupied = True
while occupied:
#find random position
row = numpy.random.randint(0, self.grid_size)
col = numpy.random.randint(0, self.grid_size)
#determine if space is empty
if numpy.isnan(self.space[row, col]):
occupied = False
else:
occupied = True
#determine tree species based on proportion of ash trees
if numpy.random.random <= self.ash_prop:
self.space[row, col] = tree.tree_id
tree.tree_species == "ash"
tree.xloc = row
tree.yloc = col
else:
self.space[row, col] = tree.tree_id
tree.tree_species == "nonash"
tree.xloc = row
tree.yloc = col
# create beetles based on model parameters
for agent_id in range(self.initial_beetles):
beetle_id = self.tree_population + agent_id
self.agent_list.append(Beetle(model = self, beetle_id,
lifestage = 4,
num_eggs_laid = numpy.random.uniform(25,45),
egg_survival = numpy.random.uniform(min_egg_survival, max_egg_survival),
winter_survival = numpy.random.uniform(min_winter_survival, max_winter_survival),
wanderlust))
def infest_initial_tree(self):
"""
choose an ash tree to be randomly infested with each initial beetle as specified in model
parameters. Keeps randomly choosing trees until an ash tree is found.
"""
for beetle in range(self.agent_list.Beetle):
# loop until true
ashtest = False
while ashtest:
random_infested = numpy.random.choice(range(self.agent_list.Tree))
if tree.is_ash(random_infested) is True:
ashtest = True
self.trees[random_infested].num_adults += 1
# add beetle to grid cell where tree is located. This record keeping method probably isn't
# valid. See other notes throughout model about making an array of lists of ID numbers
self.space[self.trees[random_infested].xloc,self.trees[random_infested].xloc,].append(beetle_id)
else:
ashtest = False
"""
QUESTION: How can I keep track of both beetles and trees within a single grid cell? I have total
beetles in each lifestage as an attribute of trees, but it is also important to keep track of
beetle_id so I can query an individual beetle's attributes in my step method.
QUESTION: I have included both the trees and the beetles in an overall agent_list, but called them
either Beetle or Tree when I appended them. Is self.agent_list.Tree the correct way to isolate
just the trees? Similar for self.agent_list.Beetle? I am following the methodology from our human/zombie
model for this piece of my ABM, but we didn't get far enough in that model to see how the different
classes get called later on.
"""
def find_neighborhood(self, x, y, distance = wanderlust):
"""
finds neighbors, accounting for grid edges, returns list of cells
locations within neighborhood (neighborhood type = extended Moore)
"""
neighbor_cells = []
try:
for x, y in itertools.product(xrange(x-wanderlust, x+wanderlust),
xrange(y-wanderlust, y+wanderlust)):
loc_tup = (x, y)
neigbor_cells.append(loc_tup)
except IndexError:
pass
return neighbor_cells
def find_hosts(self, x, y, distance = wanderlust):
"""
finds all ash trees within neighborhood, returns their locations
as a list. Based on beetle's current location and their flight
length ability (wanderlust)
"""
neighbor_cells = find_neighborhood(self, x, y, distance = wanderlust)
potential_hosts = []
for location in neighbor_cells:
if location[0] == x and location[1] == y:
continue
#we do not want to consider the tree on which the beetle
#is currently found
if not numpy.isnan(self.space[location[0], location[1]]):
#doesn't consider empty cells
if tree.is_ash is True:
#only ash trees can be hosts
potential_hosts.append(int(self.space[location[0], location[1]]))
def get_beetle_location(self, beetle_id):
"""
get the position of a beetle based on a query of their beetle_id
in self.space
"""
return numpy.reshape(numpy.where(m.space == beetle_id), (1,2))[0].tolist()
def get_beetle_neighbors(self, beetle_id):
"""
uses beetle position and flight ability to find all suitable hosts
to which the beetle is able to fly
"""
x, y = get_beetle_location(self, beetle_id)
return find_hosts(x, y, wanderlust)
def move_adults(self, beetle_id, x, y):
"""
finds all ash trees within range of beetle. Then, if ash tree is overbooked i.e. it cannot
host any more larvae based on its num_larvae and dbh, the beetle will try again. If the
beetle cannot find a new tree, it stays put. If the beetle moves, it is removed from
the old tree. Then, adults lay eggs.
"""
original_location = self.get_beetle_location(beetle_id)
potential_locations = get_beetle_neighbors(self, beetle_id)
for i in range(potential_locations):
overbooked = True
# loop until false
try:
while overbooked = True:
# random choice of potential hosts
new_host = numpy.random.choice(range(potential_locations))
# evaluates whether tree has room for more larvae or not
total_possible_larvae = tree.evaluate_larva_num(self)
if self.tree[new_host].num_larvae >= total_possible_larvae:
overbooked = True
else:
overbooked = False
self.space[original_location[0], original_location[1]].remove(beetle_id)
self.space[x,y].append(beetle_id)
self.tree[new_host].num_adults += 1
self.tree[new_host].num_eggs += self.beetle[beetle_id].num_eggs_laid
except IndexError:
continue
def step_update (self):
"""
updates tree properties based on presence/absence of beetles. The following are updated:
(1) if tree is uninfected, dbh increases
(2) if tree is infected, health decreases minimally according to num_adults (defoliation)
and health decreases more rapidly according to num_larvae.
(3) if tree health reaches zero, kill tree.
updates beetle lifestages:
(1) a proportion of eggs become larva based on egg_survival. Update tree.larva_num
(2) a proportion of larva become dormant according to tree.evaluate_larvae_num. Update tree.dormant_num
(3) a proportion of dormant beetles become adults. Update tree.adult_num
(4) if lifestage = 4 (beetle is adult), then call move_adults else, beetle stays put
(5) if tree dies, kill eggs, larva, and dormant. dormant beetles that become adults will search
for new tree to colonize.
"""
for beetle in range(self.agent_list.Beetle):
if self.agent_list.Beetle[beetle].lifestage = 4:
# delete dead beetles
self.agent_list.Beetle.pop([beetle])
elif self.agent_list.Beetle[beetle].lifestage = 1:
if numpy.random.random() <= self.egg_survival:
self.agentlist.Beetle[beetle].lifestage += 1
else:
self.agent_list.Beetle.pop([beetle])
elif self.agent_list.Beetle[beetle].lifestage = 2:
# check to see if tree can host more larva - if it can, grow them
# if it can't, then larva will die
elif self.agent_list.Beetle[beetle].lifestage = 3:
if numpy.random.random() <- self.winter_survival:
self.agent_list.Beetle[beetle].lifestage += 1
move_adults(self.agent_list[beetle].beetle_id)
else:
self.agent_list.Beetle.pop([beetle])
With each simulation, I’ll record the number of timesteps it took for all the ash trees in my simulated forest to be infected, which will tell me how well the initial parameters were able to curb EAB infestation rates. By sweeping model parameters over a large number of model runs, I can generalize spread rate for each parameter set and by comparing these spread rates, I can begin to draw conclusions about which spatial/environmental patterns are best able to slow EAB spread.
I think the most effective way to view these results might be a 3d graph in which the x- and y- axes represent the array of parameter values that were swept and the z-axis corresponds to the average length of time for the model to stop (stopping condition = all trees dead).
I believe that I will find that decreased tree density, increased biodiversity, and increased tree health will result in slower spread rates and therefore, slower tree mortality. However, I expect that these things will not decrease spread rate by a very large margin, given the extremely efficient reproduction and spread of the Emerald Ash Borer.