Teaching population dynamics with simulations in R

Graphic: Results of a discrete-time simulation with two competitors having a shared predator. Exercise for the reader: which trace is the predator?

Warning: wonkish. Of interest primarily to those who teach upper-level ecology courses.

I don’t have an important message today, or a big unresolved question to talk about. I just thought I’d share some teaching resources. If you teach ecology (past the introductory level), you may find this useful.

One of the major themes in my 3rd-year population ecology course is the diversity of population dynamics that can emerge even in fairly simple systems: extinction, stable equilibrium, damped oscillations, stable limit cycles, neutral cycles, chaos, and so on. We spend a lot of time on the kinds of ecology that tend to favour oscillations (things like time lags and enemy-victim interactions) as opposed to those that tend to favour stable equilibria (things like immediate density-dependence and, under some circumstances, interspecific competition).

The problem is that I can show examples of real population time series of various sorts, but the underlying mathematics (without which you can’t infer much from the time series) are rather abstract, and bridging the two is tough. I explore the dynamics using standard differential-equation and difference-equation models, starting with the geometric and exponential growth equations and building from there up to predator-prey and two-competitor models. As I argued in my widely-reviled post Do Biology Students Need Calculus?, there is little to be gained by asking student to integrate most of these (and most half-way interesting models are beyond 1st-year calculus integration anyway). So we examine the models graphically, which reveals a lot about their behaviour to those students who are comfortable with equations and modeling – but tends to freaks out those who are not.

So how can my students get some intuitive feel for how the simple population models are behaving, and what this means for how ecological factors determine the behaviour of real populations? I’ve been encouraged by the results of a lab exercise in which the students use simulations to explore two models – one very simple, and one more complex. I’m aware that there’s nothing original about this. Many courses have done the same, and there are many approaches – from simulation-friendly languages like MatLab and Maple to custom ecology-teaching software like EcoBeaker. But my lab is, I think, simple enough to give students a taste for how simulation works, but also rich enough to make the point that population dynamics can be really interesting and that any ecologist can use some coding skills to help understand them. I’ll post here both the population simulator lab handout and the two R scripts I have the students run; feel free to use them (non-commercially) if you like.

chaosThe first script is for a simple one-species model: the discrete-time logistic. By tuning parameters, students can see that the model generates population decline, population growth to equilibrium, damped oscillations, stable oscillations, and chaos; and they can explore the properties of chaotic dynamics (as seen in the plot at right). I even learned something about chaos myself from coding the simulations. I’d been teaching for years that chaotic dynamics are messy and easy to confuse with randomness. What I’d missed was the tendency for repeating “motifs” (look at the three similar runs beginning around generation 125). These occur in chaotic time series but not in random ones*, and it’s quite embarrassing that I’ve been teaching about chaos for years without recognizing that feature of it.

The second script is for a three-species model: two competitors being exploited by a shared predator. My students have seen both the Lotka-Volterra predator-prey model and the Lotka-Volterra competition model in lecture, and this lets them see those models in action – and also asks them to think about what happens as we move beyond 2-species interactions. I ask the students to explore how increasing a predator’s efficiency affects its population size (the answer isn’t what you probably think). I challenge them to find parameter combinations yielding all four possible outcomes of competition (competitive exclusion of 1 by 2, competitive exclusion of 2 by 1, exclusion dependent on starting conditions, and stable coexistence). With all three species in the model, I ask them to think about how predation can change the outcome of competition – something they haven’t seen in lecture. (The graphic at the top of this post shows predator-driven oscillations of the 3-species system).

A word about the approach: if you look at the materials, you’ll see I have the students tune parameters by directly editing the R scripts. This may seem clumsy – and I could easily have programmed something that offered them a dialog box or some kind of graphical control. There’s a reason I didn’t do that: I wanted the students to engage at least a little with what’s under the hood: I wanted them to understand that the simulations were happening because of computer code I had written, and that they could edit that code to change what we’re simulating. I did not require that they could write code themselves, or even that they understand my code in detail. I’d love to teach that course, but it would have to be a smaller-enrolment elective course. In addition, our students are gradually becoming exposed to R for statistics, and seeing it used in a simulation context gives some continuity and a broader understanding of what you can do with some coding experience.

This lab isn’t going to work miracles. Many of my students are still mathematically uncomfortable. I’m not teaching them to program in a single lab, and many of them will forget all the population dynamics they’ve learned within a week of the final exam. And no matter how often or how loudly I tell them that reviewing their lecture notes beforehand will make the lab much easier, to my enormous frustration, most of them just won’t**. But I do think the lab has made a difference: more students understand oscillations and what drives them, and some are even quite excited to see the simulations working. That’s a win, in my books; and if you can use my materials to similar effect, please go ahead.

© Stephen Heard (sheard@unb.ca) March 24, 2016

Related posts:

*^Roughly, this is because in chaotic dynamics small changes in initial populations lead to large changes in subsequent ones – but if the population sizes at times t1 and t2 are very similar, then the population size at t1+1 and t2+1 will be fairly similar. It takes a few generations for the two sequences to wander away from each other. If that doesn’t make intuitive sense, you can use the simulator to explore it!

**^I’m still looking for the magic elixir that will make students into professional learners, equipped with behaviours that let them take full advantage of the learning opportunities they’re offered. I’ll have more to say about this in a future post.

12 thoughts on “Teaching population dynamics with simulations in R

  1. seeddispersal

    We used to use Felsenstein’s PopG (http://evolution.gs.washington.edu/popgen/popg.html) for population genetic exercises. We’d give exercises asking students to identify key qualitative behaviours of the model and later assess what they learned by asking them to `reverse engineer’ simulation outputs, like what you’ve asked here. I felt like there was a great deal to be learned from making students sit and observe/dissect simulation outputs, and getting them asking questions like: do trajectories ever reach zero, are populations/alleles ever recovered, what are the relative magnitudes of the populations’ amplitudes, etc. I’ve not tried this with a simpler language that students could modify (R over Java for PopG) or for population ecology, but I like that idea from a heuristic perspective.
    Per the question, since I read competition first, I saw that blue and red were diverging before the lag, which suggested that black is the predator.
    Thanks for the thoughts!

    Liked by 1 person

  2. Amy Hurford

    Thanks for sharing this. I’ve been teaching Population biology and evolutionary ecology and been looking for ideas for labs. Qubes Hub is a resources where instructors can share materials pertaining to undergraduate education and quantitative biology. It might be worthwhile to post this material there: https://qubeshub.org/ I joined, but I’ve just started teaching Popn Biology and I haven’t yet developed resources that I’m fully happy with.


  3. Jeff Houlahan

    Hi Steve, thanks a lot for posting this. Coincidentally I was just about to go searching for some R code to demonstrate deterministic chaotic dynamics (I’m reading Case’s Intro to Theoretical Ecology and trying to play with some of the ideas) and popped in on your blog and found this. A very happy coincidence. best, Jeff H


  4. Pingback: Is this blog a “science blog”? If not, what is it? | Scientist Sees Squirrel

  5. Pingback: Squirrels, serendipity, and the reach of a blog | Scientist Sees Squirrel

  6. Pingback: Yes, I’m teaching a tropical field course. No, it’s not a vacation. | Scientist Sees Squirrel

  7. Arlin Stoltzfus

    Steve, thanks, this helped me to get pretty quickly to a presentation of a simple chaos example. The following code covers a series of values for the growth parameter and makes plots illustrating different regimes of behavior including cycles of period 2, 4 and 8.

    # parameters
    gens <- 100
    N0 <- 50 # initial population
    K <- 500 # capacity
    Rvals <- c(1, 1.9, 2.1, 2.5, 2.56, 2.7)
    Rdynamics <- c("asymptotic", "period 1", "period 2", "period 4", "period 8", "chaos")

    for (R in Rvals) {
    # initialize and run
    N <- N0
    result <- N
    for (i in 1:gens) {
    N <- N * (1 + R * (1 – N/K))
    if (N < 0) {
    N <- 0
    result <- c(result, N)
    # plot
    png(paste("R", R, "-plot.png", sep = ""), width = 600, height = 400)
    maintext <- paste("R = ", R, " (", Rdynamics[which(Rvals == R)], ")", sep = "")
    plot(0:gens, ylim = c(100, 700), result, xlab = "time (generation)", ylab = "Population size", type = "l", lwd = 2.5, main = maintext)



Comment on this post:

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.