Foxes & Chickens: A population dynamics model in five lines of Python

fox_chicken.jpg

As part of the work we do for our digital biology consulting firm Amber Biology, one of the activities that we enjoy the most is teaching computer programming to life science researchers. A couple of years back, we published our book Python For The LIfe Sciences, and have since adapted it to create a unique, workshop-based computer programming course specifically tailored to the needs of life scientists who want to be able to use computers in their own research. This course has already been taught at MIT and Harvard Medical School, and unlike most computer programming courses, almost all of the examples and tutorials are drawn directly from biology, and cover a broad spectrum of the kind of quantitative and computational problems that confront researchers in the life sciences. In this example taken from our course, we will implement the very simplest kind of deterministic population dynamics model imaginable, and provide a practical demonstration of the power of Python for modeling dynamic systems, .

First a little of the back story behind the model ...

Fur trappers working for the Hudon's Bay Company in Canada over an extended period from the 18th to the 20th century, kept copious records of the pelts they traded, of lynx and of the snowshoe hare that the lynx preyed upon. When the 20th century biologist Charles Gordon Hewitt examined those data closely, he observed a curious cyclical behavior. The populations of both species seemed to follow "boom and bust" cycles that were out of phase with one another - the spikes in the lynx population always lagging behind those in the hare population.

canadian_lynx.jpeg
snowshoe_hare.jpeg

This oscillatory behavior is an emergent property of predator-prey systems, and arises from the fact that the size of each population is dependent upon the size of the other. When hares are abundant the lynxes eat well and with plenty of resources available, their population booms. But once the lynx population reaches a certain size, its consumption of hares takes a toll on the hare population which starts to decline. As the hare population falls below the level at which a large lynx population can be supported, the lynx population similarly starts to decline and ... well you get the idea!

lv.png

In the 1920s, two scientists Alfred Lotke and Vito Volterra, formalized this oscillatory behavior in a famous set of equations that bears their names, and here we are going to use Python to build a very simple model of one of these predator-prey systems consisting of foxes and chickens. The actual model itself is essentially only 5 lines of Python code, yet despite its simplicity it is still able reproduce the 'boom and bust' behavior of the two interdependent populations of predator and prey.

First, we're going to set up two Python lists to store the population sizes of the foxes and chickens as they evolve over the course of our model simulation. As the simulation proceeds, the population size for the two species at each point in time, will be appended to their respective lists. We will be initializing the chicken population with 100 individuals, and the fox population with 10 individuals. 

chickens = [100]
foxes = [10]

Next, we define the rates for the events that can occur within our simulation - these are the birth rate of the prey (the chickens), the death rate of the prey (as a result of being eaten by foxes), the birth rate of the predator (the foxes), and the (natural) death rate of the predator. In the interests of keeping this model very simple and concise for this tutorial, we will make some simplifying assumptions, such as the one that all chickens succumb to predation before they have had a chance to live out their natural lifespan. We also do not adjust the lifespan of the foxes according to the availability of their food supply - so in effect, we are not accounting for any early deaths amongst the foxes due to starvation when their food supply is low.

chicken_birth_rate = 0.5
chicken_death_rate = 0.015
fox_birth_rate = 0.015
fox_death_rate = 0.5

We're almost ready for the Python code that defines the model itself, but first we need to define a time step for each cycle of the simulation and the number of cycles that we would like to run the simulation for.

delta_time = 0.01
cycles = 4500

And here's the model itself:

for t in range(0, cycles):
  updated_chickens = chickens[t] + delta_time * (chicken_birth_rate * chickens[t]  - chicken_death_rate * foxes[t] * chickens[t])    
  updated_foxes = foxes[t] + delta_time * (-fox_death_rate * foxes[t] + fox_birth_rate * foxes[t] * chickens[t])
  chickens.append(updated_chickens)
  foxes.append(updated_foxes)

Within a Python for loop that iterates over the numerical range from 0 to cycles, for each cycle we take the current values for the population size of the chickens and foxes, and update each of them in turn, according to their corresponding birth and death rates for the elapsed delta_time interval. You will also notice that the number of foxes in the current population appears as a multiplier to the death rate of the chickens, i.e. the more foxes there are, the more rapidly chickens die (through predation). Similarly, the chicken population size appears as a multiplier in the foxes birth rate, i.e. more foxes are born when their food is plentiful. At the end of the cycle, we append the updated values for each population to the corresponding list, such that these updated values can serve as the starting values for the subsequent cycle. And that's pretty much it!

Now we can run the model, but it would be nice to be able to visualize its output, and for that we will use the extremely useful matplotlib plotting library that comes bundled with the official Python distribution. Using the fox and chicken population lists that we generated using the model, we can plot their respective population sizes over the course of the simulation, like this:

import matplotlib.pyplot as plt
time_points = range(cycles + 1)
plt.plot(time_points, foxes)  
plt.plot(time_points, chickens) 
plt.xlabel('time')
plt.show()

Which gives us this:

fcplot.png

Using a few more of the many (many) features of matplotlib, we could have also labeled the two curves and included a small key in the graph, but since we have only two curves, it's pretty easy to see that the green curve represents the chickens and the blue curve, the foxes. Starting at 100 chickens and 10 foxes, both populations initially rise until the increasing fox population starts to thin the chicken population as a result of predation. The initially large chicken population starts to fall and eventually it reaches a point at which the growth in the fox population is no longer sustainable, and it also starts to fall. Once the struggling fox population has bottomed out as a result of the shortage of food, the chicken population, unburdened by the previous high levels of predation, bounces back, eventually driving a similar recovery in the fox population - and so the boom and bust cycles repeat, with the peaks and troughs of the two populations always slightly out of phase. This broadly recapitulates the patterns witnessed by Charles Gordon Hewitt in the Hudon's Bay Company's data for the lynx and snowshoe hare pelts they traded, and later formalized for a system of two interdependent predator-prey populations in the Lotke-Volterra equations.

In practice, most population dynamics models would be considerably more complex than the one presented here, but we like to use this introductory example in our Python For LIfe Scientists workshop course, as a simple and intuitive illustration of how useful Python can be for modeling biological systems. Beyond this kind of simple deterministic model, we also take the students through workshops in which we implement more advanced dynamic modeling algorithms such as the Gillespie Algorithm, which we use to generate stochastic, but kinetically accurate models of cell signaling pathways.

The course itself is very hands-on, with these kinds of programming tutorials regularly interspersed between the more traditional classroom sessions. The tutorials were also chosen to represent a range of biological research areas including genomics and next-generation sequencing, structural biology, cell signaling pathways, population dynamics, laboratory automation and many other examples of the application of Python to life science research. The course consists of 20 - 40 hours of classroom and workshop sessions and can be taught over several weeks, or as an intensive course over a few days. We also teach the course on-site at our clients' locations, so their researchers don't need to travel, and have the opportunity to learn Python programming right where they work. 

If you would be interested in having the course Python For LIfe Scientists  taught at your organization, please contact us and we would be happy to schedule a workshop at your location, for a class of 6 or more students.

© The Digital Biologist