Jupyter at Bryn Mawr College
Public notebooks: /services/public/dblank / BioCS115 Computing through Biology / 2016-Spring / Notebooks

1. Lab 9: Ecological Models

[In the cell below, put your name; then you can delete this cell]

YOUR ANSWER HERE

Goals:

• Explore ecological models through "individual/instance simulation"
• Observe complex interactions between agents
• Compare and contrast with mathematical models
• Practice using OOP
• Use widgets to run experiments

We now turn our attention to a different kind of biological modeling: ecological models. These models often involve spatial territory, and a possibly large number of individual agents woondering the landscape. This makes it perfect for using our new Object-Oriented Programming skills.

The model that we will develop today is designed to explore the stability of predator-prey ecosystems. Such a system is called unstable if it tends to result in extinction for one or more species involved. In contrast, a system is stable if it tends to maintain itself over time, despite fluctuations in population sizes.

This Python code and Jupyter notebook are based on a NetLogo model: http://netlogoweb.org/launch#http://netlogoweb.org/assets/modelslib/Sample%20Models/Biology/Wolf%20Sheep%20Predation.nlogo

1.1 Predator-Prey Model

There are two main experiments of exploration to this model.

In the first experiment, wolves and sheep wander randomly around the landscape, while the wolves look for sheep to prey on. Each step costs the wolves energy, and they must eat sheep in order to replenish their energy - when they run out of energy they die. To allow the population to continue, each wolf or sheep has a fixed probability of reproducing at each time step. This variation produces interesting population dynamics, but is ultimately unstable.

The second experiment includes grass (green) in addition to wolves and sheep. The behavior of the wolves is identical to the first point, however this time the sheep must eat grass in order to maintain their energy - when they run out of energy they die. Once grass is eaten it will only regrow after a fixed amount of time. This version is more complex than the first, but it is generally stable.

To put it another way: when the sheep have unlimited amount of grass, this leads to an unstable system. But putting some limits on the amount of grass generally produces a stable environment.

The construction of this model is described in two papers by Wilensky & Reisman referenced below.

1.2 How to use it

1. Set the grass parameter to True to include grass in the model, or to False to only include wolves (black) and sheep (white/yellow).
2. Adjust the parameters (see below), or use the default settings.
3. Instantiate the World class.
4. Call update to step the simulation.
5. Look at the history to see how the population sizes change over time.

Parameters:

• World.sheep: The initial size of sheep population
• World.wolves: The initial size of wolf population
• World.grass: Whether or not to include grass in the model
• World.grass_regrowth_time: How long it takes for grass to regrow once it is eaten
• Sheep.GAIN_FROM_FOOD: The amount of energy sheep get for every grass patch eaten
• Sheep.REPRODUCE: The probability of a sheep reproducing at each time step
• Wolf.GAIN_FROM_FOOD: The amount of energy wolves get for every sheep eaten
• Wolf.REPRODUCE: The probability of a wolf reproducing at each time step

Notes:

• one unit of energy is deducted for every step a wolf takes
• when grass is included, one unit of energy is deducted for every step a sheep takes

2. Wolves, Sheep, and Grass

There are an infinite number of ways to break the wolf, sheep, and grass model into Python objects. The following code explores just a single way.

2.1 World

The world is the class that contains all of the elements of the model. It contains a list of animals (animals that are either Wolf or Sheep) and also a discrete grid of Patches that contain the grass level for each fixed area. By default, the grid of grass is composed of 50 x 50 cells, but connected as a torus (wraps around).

Methods:

• World() - construct a new model
• World.setup() - call the setup methods, if parameters change
• World.get_stats() - get the populations of (sheep, wolves, grass/4)
• World.update() - update the world by one time tick
• World.draw() - create a nice graphical picture of the world

2.2 Patches

The world is made up of 50 x 50 patches. A Patch contains a level of grass. If World(grass=True), then Sheep can eat the grass, and it will eventually grow back.

Properties:

• Patch.level - grass level (states; -1, 0, 1, or infinity)
• Patch.time - counting time until regrowth

Methods:

• Patch.update() - update the grass (head toward regrowth, if grass enabled)

2.3 Animals

Wolf and Sheep are both Animals. Animals has all of the properties and methods that they both share.

Properties:

• Animal.x, Animal.y - the location of the animal (floats, between 0 - 50)
• Animal.energy - the amount of energy this animal has

Methods:

• Animal.update() - calls: move(), eat(), hatch(), or dies
• Animal.move() - move in your random way
• Animal.eat() - try to eat
• Animal.hatch() - splits yourself
• Animal.forward(distance) - go forward (used in move)
• Animal.turnLeft(direction) - turn left in degrees (used in move)
• Animal.turnRight(direction) - turn right in degrees (used in move)
• Animal.eats(other_type) - can I eat this type?
• Animal.distance(a, b) - compute distance on torus

Now, we import all of the libraries that we will need:

import random
import math
from IPython.display import clear_output, display
from IPython.display import HTML as html
import time
import io
import matplotlib.pyplot as plt
import matplotlib.patches
from calysto.graphics import *
import threading

%matplotlib inline

And define the World class that will represent the entire model:

class World:
    def __init__(self, pwidth=50, pheight=50, grass=True, 
                sheep=100, wolves=50, grass_regrowth_time=30):
        # Width and Height (in patches):
        self.pwidth = pwidth
        self.pheight = pheight
        # -----------------------------------
        # Parameters:
        # -----------------------------------
        # grass:
        self.use_grass = grass
        self.grass_regrowth_time = grass_regrowth_time # 0 to 100
        # sheep:
        self.initial_number_sheep = sheep # 0 to 250
        # wolves:
        self.initial_number_wolves = wolves # 0 to 250
        self.state = "stopped"
        self.setup()
        
    def setup(self):
        self.history = []
        self.ticks = 0
        if self.use_grass:
            self.patches = [[Patch(self, random.choice([-1, 1])) 
                             for w in range(self.pheight)] for h in range(self.pwidth)]
        else:
            self.patches = [[Patch(self, math.inf) 
                             for w in range(self.pheight)] for h in range(self.pwidth)]
        self.animals = []
        for i in range(self.initial_number_sheep):
            self.animals.append(Sheep(self))
        for i in range(self.initial_number_wolves):
            self.animals.append(Wolf(self))
        
    def draw(self, pixels=10):
        canvas = Canvas(size=(self.pwidth * pixels, self.pheight * pixels))
        if self.use_grass:
            grass = 0
            non_grass = 0
            if self.use_grass:
                for h in range(self.pheight):
                    for w in range(self.pwidth):
                        if self.patches[w][h].level > 0:
                            grass += 1
                        else:
                            non_grass += 1
            if non_grass > grass:
                rectangle = Rectangle((0,0), (self.pwidth * pixels - 1, self.pheight * pixels - 1))
                rectangle.draw(canvas)
                rectangle.fill("brown")
                for h in range(self.pheight):
                    for w in range(self.pwidth):
                        if self.patches[w][h].level > 0:
                            rect = Rectangle((w * pixels, h * pixels), 
                                             (pixels, pixels))
                            rect.fill("green")
                            rect.stroke("green")
                            rect.draw(canvas)
            else:
                rectangle = Rectangle((0,0), (self.pwidth * pixels - 1, self.pheight * pixels - 1))
                rectangle.draw(canvas)
                rectangle.fill("green")
                for h in range(self.pheight):
                    for w in range(self.pwidth):
                        if self.patches[w][h].level <= 0:
                            rect = Rectangle((w * pixels, h * pixels), 
                                             (pixels, pixels))
                            rect.fill("brown")
                            rect.stroke("brown")
                            rect.draw(canvas)
        else:
            rectangle = Rectangle((0,0), (self.pwidth * pixels - 1, self.pheight * pixels - 1))
            rectangle.draw(canvas)
            rectangle.fill("green")
        for animal in self.animals:
            animal.draw(canvas, pixels)
        rectangle = Rectangle((0,0), (self.pwidth * pixels - 1, self.pheight * pixels - 1))
        rectangle.draw(canvas)
        rectangle.fill(None)
        return canvas
        
    def get_stats(self):
        sheep = 0;
        wolves = 0
        for animal in self.animals:
            if isinstance(animal, Sheep):
                sheep += 1
            elif isinstance(animal, Wolf):
                wolves += 1
        grass = 0
        if self.use_grass:
            for h in range(self.pheight):
                for w in range(self.pwidth):
                    if self.patches[w][h].level > 0:
                        grass += 1
        return (sheep, wolves, grass/4)
        
    def update(self):
        for animal in list(self.animals): ## copy of list in case changes
            animal.update()
        for h in range(self.pheight):
            for w in range(self.pwidth):
                self.patches[w][h].update()
                
    def stop(self):
        self.state = "stopped"
                
    def run(self, widgets, width=4.0, height=3.0, pixels=10):
        self.state = "running"
        while 0 < len(world.animals) and self.state == "running":
            self.ticks += 1
            self.update()
            stats = self.get_stats()
            self.history.append(stats)
            widgets["plot"].value = self.plot(width, height)
            widgets["canvas"].value = str(world.draw(pixels))
            widgets["sheep"].value = str(stats[0])
            widgets["wolves"].value = str(stats[1])
            widgets["grass"].value = str(stats[2])
            widgets["ticks"].value = str(self.ticks)
        self.state = "stopped"
            
    def plot(self, width=4.0, height=3.0):
        plt.rcParams['figure.figsize'] = (width, height)
        fig = plt.figure()
        plt.plot([h[0] for h in self.history], "y", label="Sheep")
        plt.plot([h[1] for h in self.history], "k", label="Wolves")
        plt.plot([h[2] for h in self.history], "g", label="Grass")
        plt.legend()
        plt.xlabel('time')
        plt.ylabel('count')
        bytes = io.BytesIO()
        plt.savefig(bytes, format='svg')
        svg = bytes.getvalue()
        plt.close(fig)
        return svg.decode()
            
    def __repr__(self):
        matrix = [[" " for w in range(self.pheight)] for h in range(self.pwidth)]
        for h in range(self.pheight):
            for w in range(self.pwidth):
                matrix[w][h] = repr(self.patches[w][h])
        for animal in self.animals:
            matrix[int(animal.x)][int(animal.y)] = animal.SYMBOL
        string = ""
        for h in range(self.pheight):
            for w in range(self.pwidth):
                string += matrix[w][h]
            string += "\n"
        return string
        

The Patch class represents a single cell on the 50 x 50 matrix for keeping track of grass:

class Patch:
    def __init__(self, world, level):
        self.world = world
        self.level = level
        if level == -1:
            self.time = random.randint(0, self.world.grass_regrowth_time)
        else:
            self.time = 0            
        
    def update(self):
        if self.level == math.inf:
            pass
        elif self.level == 0:
            self.time = 0
            self.level = -1
        elif self.time == self.world.grass_regrowth_time:
            self.level = 1
        else:
            self.time += 1
            
    def __repr__(self):
        if self.level > 0:
            return "."
        return " "

The "base class" for our Sheep and Wolf classes:

class Animal:
    GAIN_FROM_FOOD = 1
    REPRODUCE = 0 ## percentage
    SYMBOL = "a"
    
    def __init__(self, world):
        self.world = world
        self.x = random.random() * self.world.pwidth
        self.y = random.random() * self.world.pheight
        self.direction = random.random() * 360
        
    def turnLeft(self, degrees):
        self.direction -= degrees
        self.direction = self.direction % 360
    
    def turnRight(self, degrees):
        self.direction += degrees
        self.direction = self.direction % 360

    def forward(self, distance):
        radians = self.direction * (math.pi/180.0)
        x = math.cos(radians) * distance
        y = math.sin(radians) * distance
        self.x = (self.x + x) % self.world.pwidth
        self.y = (self.y + y) % self.world.pheight
    
    def update(self):
        self.move()
        self.eat()
        if self.energy > 0:
            if random.random() < self.REPRODUCE:
                self.hatch()
        else: # die
            self.world.animals.remove(self)
    
    def move(self):
        self.turnRight(random.random() * 50)
        self.turnLeft(random.random() * 50)
        self.forward(1)
        if isinstance(self, Wolf) or (isinstance(self, Sheep) and self.world.use_grass):
            self.energy -= 1.0
        
    def hatch(self):
        child = self.__class__(self.world) ## hatch an offspring 
        child.x = self.x
        child.y = self.y
        self.energy /= 2                   ## divide energy between parent and offspring
        child.energy = self.energy
        child.forward(1)                   ## move offspring forward 1 step
        self.world.animals.append(child)   ## add offspring to list of animals

    def eats(self, other):
        return False
    
    def eat(self):
        pass
        
    def distance(self, a, b):
        # minimal distance on torus
        dx = abs(a.x - b.x)
        dy = abs(a.y - b.y)
        return math.sqrt(min(dx, self.world.pwidth - dx) ** 2 + 
                         min(dy, self.world.pheight - dy) ** 2)
            
    def __repr__(self):
        return self.SYMBOL
        

And now our extended classes:

class Wolf(Animal):
    GAIN_FROM_FOOD = 20 # 0 to 50
    REPRODUCE = 0.05 ## percentage
    SYMBOL = "w"
    
    def eats(self, other):
        return isinstance(other, (Sheep))

    def eat(self):
        if self.energy >= 0: # if I am alive:
            for animal in list(self.world.animals):
                if (self.eats(animal) and 
                    self.distance(self, animal) < 1.0 and
                    animal.energy >= 0):
                    self.energy += self.GAIN_FROM_FOOD
                    animal.energy = -1 ## mark dead
                    break # don't eat any more!
    
    
    def __init__(self, world):
        super().__init__(world)
        self.energy = random.randint(5, 40)

    def draw(self, canvas, pixels):
        wolf = Ellipse((self.x * pixels, self.y * pixels), (pixels/2, pixels/3))
        wolf.fill("black")
        wolf.stroke("black")
        wolf.draw(canvas)

class Sheep(Animal):
    GAIN_FROM_FOOD = 4 # 0 to 50
    REPRODUCE = 0.04 ## percentage
    SYMBOL = "s"
    
    def __init__(self, world):
        super().__init__(world)
        self.energy = random.randint(1, 7)

    def eat(self):
        patch = self.world.patches[int(self.x)][int(self.y)]
        if patch.level >= 0 and patch.level < math.inf:
            self.energy += self.GAIN_FROM_FOOD
            patch.level -= 1

    def draw(self, canvas, pixels):
        sheep = Ellipse((self.x * pixels, self.y * pixels), (pixels/2, pixels/2))
        sheep.fill("white")
        sheep.stroke("white")
        sheep.draw(canvas)

3. Experiments

In this section we describe a method for running experiments, and then we will do some analysis and reflection.

3.1 Widgets

To run our experiments, we will use IPython widgets for setting parameters. Generally, you will:

1. Create a global variable world
2. Create the widgets
3. Change their values as appropriate
4. Click on "setup"
5. Click on "go"
   • Perhaps clicking on "pause" occasionally
   • Examine status
   • Display plot or world
   • Make notes
   • Continue by clicking "go" again

Make sure that your program is stopped when you stop.

First, we make a world. We keep the size small to make it run faster. Remember that the code is running on the server, but the displays (plot and world) must be computed and sent to your browser.

world = World(25, 25, grass=True, sheep=50, wolves=25)

To make our Graphical User Interface (or GUI, pronounced "gooey") we will import the IPython widget module:

from ipywidgets import *

We will use the following widgets: Text, HTML, Button, Checkbox, IntSlider, and FloatSlider. In addition, we will use HBox (horizontal box) and VBox (vertical box) for layout design.

Widgets are useful because they allow the controls in the browser to send commands to Python code, and vice versa.

Let's take a quick look at an example. First, we make a variable (called variable) and set it to an IntSlider (a slider control that sets integers). We can create it in the GUI:

variable = IntSlider(description="Setting Name:")
variable

And see its value, initially 0:

variable.value

0

If we change it in code, it changes in the GUI:

variable.value = 10

And if we change it in the GUI (slide slider to 31) then we see it in code:

variable.value

10

In the following code, note that:

• "px" stands for pixels in HTML (HyperText Markup Language)
• Our code will run in a "thread" which means it runs in the background. That makes it so that we can click a button while the program is running

# widgets:
ticks = Text(description="Ticks:")
sheep = Text(description="sheep:", margin=5, width="60px")
wolves = Text(description="wolves:", margin=5, width="60px")
grass = Text(description="grass/4:", margin=5, width="60px")
canvas = HTML(margin=10)
plot = HTML()
setup_button = Button(description="setup", width="47%", margin=5)
go_button = Button(description="go", width="47%", margin=5)
grass_checkbox = Checkbox(description="grass?", margin=5)
grass_regrowth_time = IntSlider(description="grass_regrowth_time:",
                               min=0, max=100, margin=5, width="300px")
sheep_slider = IntSlider(description="initial_number_sheep:", min=0, max=250, margin=5, width="100px")
wolves_slider = IntSlider(description="initial_number_wolves:", min=0, max=250, margin=5, width="100px")
sheep_gain_from_food_slider = IntSlider(description="sheep_gain_from_food:", min=0, max=50, margin=5, width="100px")
wolf_gain_from_food_slider = IntSlider(description="wolf_gain_from_food:", min=0, max=100, margin=5, width="100px")
sheep_reproduce_slider = FloatSlider(description="sheep_reproduce:", min=0, max=20, margin=5, width="100px")
wolf_reproduce_slider = FloatSlider(description="wolf_reproduce:", min=0, max=20, margin=5, width="100px")

# layout:
row1 = HBox([setup_button, go_button], background_color="lightblue", width="100%")
row2 = HBox([grass_checkbox, grass_regrowth_time], background_color="lightgreen", width="100%")
row3 = HBox([sheep_slider, wolves_slider], background_color="lightgreen", width="100%")
row4 = HBox([sheep_gain_from_food_slider, wolf_gain_from_food_slider, ], background_color="lightgreen", width="100%")
row5 = HBox([sheep_reproduce_slider, wolf_reproduce_slider], background_color="lightgreen", width="100%")
row61 = HBox([sheep], background_color="Khaki", width="100%")
row62 = HBox([wolves], background_color="Khaki", width="100%")
row63 = HBox([grass], background_color="Khaki", width="100%")
row7 = HBox([plot], width="100%")
widgets = HBox([VBox([row1, row2, row3, row4, row5, 
                      row61, row62, row63, row7], width="60%"), 
                VBox([ticks, canvas], background_color="Khaki")], width="100%")

# update:
ticks.value = str(world.ticks)
canvas.value = str(world.draw(15))
plot.value = world.plot(6.0, 3.0)
grass_checkbox.value = world.use_grass
grass_regrowth_time.value = world.grass_regrowth_time
sheep_slider.value = world.initial_number_sheep
wolves_slider.value = world.initial_number_wolves
sheep_gain_from_food_slider.value = Sheep.GAIN_FROM_FOOD
wolf_gain_from_food_slider.value = Wolf.GAIN_FROM_FOOD
sheep_reproduce_slider.value = Sheep.REPRODUCE
wolf_reproduce_slider.value = Wolf.REPRODUCE
stats = world.get_stats()
sheep.value = str(stats[0])
wolves.value = str(stats[1])
grass.value = str(stats[2])
ticks.value = str(world.ticks)

def update(args):
    if args["name"] != "value":
        return
    owner = args["owner"]
    name = args["name"]
    value = args["new"]
    if owner == grass_checkbox:
        world.use_grass = value
    elif owner == grass_regrowth_time:
        world.grass_regrowth_time = value
    elif owner == sheep_slider:
        world.initial_number_sheep = value
    elif owner == wolves_slider:
        world.initial_number_wolves = value
    elif owner == sheep_gain_from_food_slider:
        Sheep.GAIN_FROM_FOOD = value
    elif owner == wolf_gain_from_food_slider:
        Wolf.GAIN_FROM_FOOD = value
    elif owner == sheep_reproduce_slider:
        Sheep.REPRODUCE = value
    elif owner == wolf_reproduce_slider:
        Wolf.REPRODUCE = value

def setup(widgets):
    world.setup()
    canvas.value = str(world.draw(15))
    plot.value = world.plot(6.0, 3.0)
    stats = world.get_stats()
    sheep.value = str(stats[0])
    wolves.value = str(stats[1])
    grass.value = str(stats[2])
    ticks.value = str(world.ticks)

def run(widgets):
    global thread
    if world.state == "stopped":
        go_button.description = "pause"
        w = {"canvas": canvas, "plot": plot, "ticks": ticks, 
             "sheep": sheep, "wolves": wolves, "grass": grass}
        thread = threading.Thread(target=lambda: world.run(w, 6.0, 3.0, 15))
        thread.start()
    else:
        go_button.description = "go"
        world.state = "stopped"
        thread.join()
        clear_output()
    
grass_checkbox.observe(update)
# FIX:
grass_regrowth_time.observe(update)
sheep_slider.observe(update)
wolves_slider.observe(update)
sheep_gain_from_food_slider.observe(update)
wolf_gain_from_food_slider.observe(update)
sheep_reproduce_slider.observe(update)
wolf_reproduce_slider.observe(update)
# END OF FIX

setup_button.on_click(setup)
go_button.on_click(run)

widgets

The above will look something like this when executed:

3.2 Things to Notice

When grass is not included, watch as the sheep and wolf populations fluctuate. Notice that increases and decreases in the sizes of each population are related. In what way are they related? What eventually happens?

Once grass is added, notice the green line added to the population plot representing fluctuations in the amount of grass. How do the sizes of the three populations appear to relate now? What is the explanation for this?

Why do you suppose that some variations of the model might be stable while others are not?

3.3 Things to Try

Try adjusting the parameters under various settings. How sensitive is the stability of the model to the particular parameters?

Can you find any parameters that generate a stable ecosystem that includes only wolves and sheep?

Try setting World.grass to True, but setting World.wolves to 0. This gives a stable ecosystem with only sheep and grass. Why might this be stable while the variation with only sheep and wolves is not?

Notice that under stable settings, the populations tend to fluctuate at a predictable pace. Can you find any parameters that will speed this up or slow it down?

Try changing the reproduction rules -- for example, what would happen if reproduction depended on energy rather than being determined by a fixed probability?

3.4 Analysis

At any point, you may want to pause and take some snapshots and make some notes:

display(html(world.plot(12.0, 6.0)))

display(html(str(world.draw(20))))

[In the cell below, enter your reflections on analysis]

In your reflections:

• Note items from the "things to notice" section
• Comment on analysis, variations
• Are there things that this model can capture that the mathematical model does not? Explain.

YOUR ANSWER HERE

3.5 Bonus Projects

3.5.1 Extend the Model

There are a number ways to alter the model so that it will be stable with only wolves and sheep (no grass). Some will require new elements to be coded in or existing behaviors to be changed. Can you develop such a version?

Can you modify the model so the sheep will flock?

Can you modify the model so that wolf actively chase sheep?

You can add cells here.

3.5.2 Further Reading

Read the paper below titled "Thinking like a Wolf, a Sheep or a Firefly: Learning Biology through Constructing and Testing Computational Theories" and provide insightful commentary. You can add cells here.

4. Reflections

[Enter your reflections in the following cell.]

YOUR ANSWER HERE

5. References

1. Wilensky, U. & Reisman, K. (1999). Connected Science: Learning Biology through Constructing and Testing Computational Theories -- an Embodied Modeling Approach. International Journal of Complex Systems, M. 234, pp. 1 - 12. (This model is a slightly extended version of the model described in the paper.)
2. Wilensky, U. & Reisman, K. (2006). Thinking like a Wolf, a Sheep or a Firefly: Learning Biology through Constructing and Testing Computational Theories -- an Embodied Modeling Approach. Cognition & Instruction, 24(2), pp. 171-209. http://ccl.northwestern.edu/papers/wolfsheep.pdf

5.1 COPYRIGHT AND LICENSE

• Copyright 1997 Uri Wilensky.
• Copyright 2016 Douglas Blank.

This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License. To view a copy of this license, visit https://creativecommons.org/licenses/by-nc-sa/3.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.