Simulating interacting populations
BIOL377
There are multiple ways that we can go about simulating species population dynamics in
order to explore the consequences of changes to those species’ abiotic and biotic
environment. For example, we might wish to explore what happens when changing
temperatures alter species metabolism and decrease or increase growth rates, how a
reduction in available habitat changes equilibrium population sizes, or how the introduction
of an exotic species might disrupt the dynamics of an undisturbed community. Each of these
disturbances, among many others, are increasingly common across modern-day
communities; however, most to all are impossible or unethical to test within real-world
conditions, hence the use of mathematical and computational models to investigate them in
silico . We will start the Research Workshop by interrogating three different (but related)
modelling scenarios that have been the subject of considerable study within the ecological
literature.
1. How to make an in silico model-based prediction
Before we can study any specific model, we need to understand the ground rules for how we
might go about using a population-dynamics model to make a prediction. By and large, there
are four critical ingredients required to make a model-based prediction about the species in
community of interest:
1. We need a mathematical description of how the size of the different populations are
expected to change over time. Typically, this mathematical model will depend on
parameters, which are things we normally treat as unvarying (e.g., temperature,
habitat availability, etc), and on variables, which are things that we allow to change
over time (e.g., the sizes of the different populations).
2. We need to decide what values to assign to the different parameters. Sometimes we
can directly relate these to measurable quantities to “ground” our simulation in a
specific realistic scenario. At other times, we will see that what the values actually
mean (and how they relate to each other) is a bit more nuanced.
3. We need to decide what the community state is at the “start” of our simulation. When
using population-dynamics models, we are effectively trying to make a prediction
about how populations will change over time. But to measure how they have
changed, we need to know where they changed from .
4. We need to decide how far into the future we wish to simulate and make our
prediction. Sometimes we are interested in making a short-term prediction, for
example if we want to know whether a new disturbance will make populations go up
or down. At other times, we are interested in making a long-term prediction because
we want to determine whether the community will recover from disturbance over a
long time horizon or whether or not it will end up in a completely different state due to
altered internal feedbacks.
2. Logistic growth of a single population
To make things as straightforward as possible, we will start by studying one of the
mathematically simplest models for a population varying over time: logistic growth. In this
model, we have a population of a single species whose abundance tends to increase when it
has a small population size relative to its carrying capacity and to decrease when it has a
large population size relative to its carrying capacity. Mathematically, the model looks like
= 1 ?
( )
where the equation describes the growth rate of the population density N at any given
moment of time. We can think of growth rate as the amount of change in population density
divided by the time over which that density changes. This implies that negative growth rates
lead to decreases in N over time and positive growth rates lead to increases in N over time.
In this model, the parameter r is the maximum growth rate of species N and K is the carrying
capacity of species N . Note that if N is much smaller than K , the total growth rate will be
positive as long as r is positive; if N = K , the total growth rate is zero; and if N > K the total
growth rate is negative. Code to explore this model in R can be found on Learn in the file
logistic_growth.R.
Though this model is simple mathematically, it sets an intuitively useful baseline for what
might happen if we add additional species or additional phenomena to the mix. This is what
is shown in the default scenarios outlined in the R code that has been made available to you.
What to do
1. Work through the R code in logistic_growth.R. Make sure that you can run it from the
top to bottom, that it all makes sense, and add comments to it in the places you think
you’d like reminders of what does what.
2. The different “experiments” in the code allow you to see what happens if you vary the
values of parameters like r and K . If you want to conduct further experiments, try
changing just one parameter at a time before you change both together. Before you
simulate what happens upon changing a parameter, try to make a guess about what
should happen based on your biological intuition. For example, if you increase K , do
you expect population sizes to go up or to go down?
3. The model as written makes predictions about what the population should do in the
absence of any external disturbance, but we know that global change is doing plenty
of things that impact species directly (and not necessarily by changing the underlying
growth rate or carrying capacity). One way to “disturb” the population is to introduce
some sort of perturbation and see what happens after (i.e. check its resilience). See
if you can change the R code you started with to explore the following simulated
disturbances:
a. Imagine that this is a threatened species. At time t=100 , we translocate 40
new individuals to the population. We then allow the population to go back to
its intrinsic dynamics. What happens to the population overall?
b. Imagine that this is a commercially viable species. At time t=100 , we remove
23% of the population. We then allow the population to go back to its intrinsic
dynamics. Can the population recover?
c. Imagine that this is a commercially viable species that we wish to exploit
continuously. Starting at time t=100 , we wish to remove a constant number of
individuals from the population at all times. What is the largest number of
individuals that we can exploit and maintain a viable population?
3. Dynamics of a predator and prey
We will now move to THE classic continuous-time model for interacting predator and prey
species: Lotka-Volterra. Mathematically, the model looks like
where the first equation describes the growth rate of N and the second the growth rate of P .
The parameter r is the maximum growth rate of N , K is the carrying capacity of N , a is the
attack rate of P attacking N , e is the efficiency of P eating N , and s is mortality rate of P . All
of these parameters should always be given positive values.
In this model, we have two different populations, N and P , who vary over time with respect to
each other. If the population of P is zero, the population of N shows logistic growth as
described above. If the population of P is not zero, then the realized growth rate of N
decreases because of the interaction with P . Likewise, if the population of N is zero, then the
population of P experiences mortality; but, if the population of N is not zero, then the realized
growth rate of P increases because of the interaction with N .
An important feature we can see in the model is that the equation for each species has both
positive and negative contributions to overall growth. Without these, we can end up with
models that don’t make biological sense. For example, the predator should die as time goes
on when there are no available prey. Likewise, the prey population should not grow forever,
even in the absence of the predator. Lastly, note that the terms that capture the interaction
between species (i) feature both densities being multiplied together (so that the interaction
“vanishes” whenever any species is absent) and (ii) are “paired” between the two models so
that loss of density of prey corresponds to a proportional gain of density for the predator.
Code to explore this model in R can be found on Learn in the file lotka_volterra.R. Similar to
the model explored above, changes to the different parameter values can change the
dynamic behavior and hence biological predictions that it makes. This is what is shown in the
default scenarios that have also been outlined in the R code.
In the R code, we start from a scenario where the prey population has an extremely large
carrying capacity and hence would grow exponentially without the predator. Here the
Lotka-Volterra model shows what’s called a fixed point equilibrium such that the populations
don’t vary over time. If we start the population away from this equilibrium, we see that the
two species will oscillate back toward this equilibrium. What does this mean in terms of what
the Lotka-Volterra model predicts for the resilience of a set of populations that might be
subject to some sort of external disturbance?
We next explore two related experiments in the code where we now have logistic growth
because intraspecific competition is no longer zero. In the first such scenario, we decrease
the prey’s carrying capacity which ends up driving the predator to extinction. (How small a
carrying capacity will allow the predator to avoid extinction?) In the second scenario, we
decrease the prey’s carrying capacity AND increase the predator’s attack rate. Somewhat
paradoxically, this helps prevent the predator from going extinct, even though the equilibrium
abundances of both predator and prey are much lower than where we started. Do these
results make sense to you? Why should competition between prey make it harder for the
predator to avoid extinction? Why would increased attack rates make things more stable?
What to do
1. Work through the R code in lotka_volterra.R. Make sure that you can run it from the
top to bottom, that it all makes sense, and add comments to it in the places you think
you’d like reminders of what does what.
2. See how the output (as captured in the plots) changes if you vary the parameters. Try
changing just one of them at a time before you change more than one. Before you
change a parameter, try to make a guess about what should happen based on your
biological intuition.
3. See if you can change the R code to explore the following scenarios:
a. The third and fourth experiments hint at the fact that the rate of intraspecific
competition and the predator attack rate combine to determine the fate of the
prey. Can you think of other parameters that you could have changed to get
similar qualitative predictions? See if you can change other parameters and
get figures that look similar to those in the middle and bottom row. Then work
out whether or not you can provide a biological explanation for why that
makes sense.
b. Explore what happens if you add another prey or predator to the system
described by the fourth experiment (since it gives stable coexistence). How
many new parameters do you need to add to describe the additional
biology/ecology? Think about how you might make simplifying assumptions
about the “behavior” of the new species. For example, can you assume that
some things can be identical to the species already present? If this is a new,
exotic species, what has to happen in order for it to establish? Does the
answer depend on if you add a new predator instead of a new prey?
4. Saturating predator functional response
All mathematical models are, at best, a caricature of reality. In the Lotka-Volterra model,
there is one assumption that has classically been regarded as less reasonable than all
others: a non-saturating functional response. If we look at the term that captures the
predator-prey interaction, we can see that the per capita rate of predation will always grow
regardless of the density of prey available. In real-world terms, this means that if you went to
a buffet the amount you would eat would always be proportional to the food available,
whether the buffet table was two meters long of two kilometers long. Plenty of empirical
evidence indicates that this is rarely the case, and instead it is much more likely that
consumption eventually saturates at some level. To build this into our model, we can change
the term for the interaction (which appears in both equations), as follows and get what is
known as the Rosensweig-MacArthur model
Note that, instead of a single attack rate a , we now have two different parameters a and h
which control the rate of consumption. In the literature, a is defined as the maximum capture
rate and h is the “handling time” of species P . If the number of prey is very very large, the
predators will never consume more than a prey per capita, but if the number of prey is small
then the consumption rate will increase proportional to N (just like in the original
Lotka-Volterra model!). Code to explore this model in R can be found on Learn in the file
rosenzweig_macarthur.R.
Similar to the models explored above, changes to the different parameter values can change
the dynamic behavior and hence biological predictions that it makes. This is what is shown in
the default scenarios outlined in the R code that has been made available to you. We start
from a scenario where the prey and predator populations can both co-exist with well-defined
equilibrium population sizes. If we increase the carrying capacity of the host, however, we
can shift to a situation where the populations cycle over time.
What to do
1. Work through the R code in rosenzweig_macarthur.R. Make sure that you can run it
from the top to bottom, that it all makes sense, and add comments to it in the places
you think you’d like reminders of what does what.
2. The difference between the experiments should remind you very strongly of
something you’ve seen in the previous examples. Can you figure out why the results
are so similar? See if you can figure out how to prove your hypothesis.
3. Explore how the output (as captured in the plots) changes if you vary the parameters.
Try changing just one of them at a time before you change both together. Before you
change a parameter, try to make a guess about what should happen based on your
biological intuition.
4. See if you can change the R code to explore the following scenarios:
a. Imagine that the habitat of the prey is undergoing land-use change. Estimate
how much habitat can be removed (via changes in the prey’s carrying
capacity) before either species is predicted to go extinct.
b. If this were a biological control situation in which a parasite was introduced to
control a host species’ population, what single process-related parameter
would be most critical to keep the host population under control and/or to
drive it extinct? How might you go about determining what the “ideal”
parasite’s behaviour should be and whether or not that combination is
realistic?
c. Many species exhibit an Allee effect, which implies that they show negative
density dependence at both low and high abundance. The observation at high
abundance corresponds to what we see with logistic growth, and the
observation at low abundance could arise because populations at very low
density have a harder time encountering each other and hence do not have
sufficient numbers to reproduce and maintain a viable population. We can
include an Allee effect in our model for the prey population by changing the
equation to be
where A determines the abundance above which the prey population can
grow (and hence should be less than K). Add an Allee effect to the R code for
the model and determine how the parameter A impacts the dynamic
predictions of the model. You can try whichever hypothetical scenarios you
wish, but you might find repeating (a) to be a useful starting point.
d. If you wish to boost the prey population by introducing a third species, will it
be better to add an alternative prey for the predator or to add a species that
attacks the predator (e.g., a top predator or parasitoid)? Consider carefully
biology/ecology, and try to start from the easiest case while working toward a
more “realistic” model. What assumptions do you need to make when adding
the new species?
5. Golden eagles, feral pigs, and insular carnivores
As discussed in the lectures, Roemer and colleagues made a number of model-based
predictions about the California Channel Islands that are a very good example of how
researchers employ mathematical models to understand global-change phenomena. The
primary article can be found at https://doi.org/10.1073/pnas.012422499 , and R code to
simulate and explore the author’s predictions can be found on Learn in the file roemer.R.