Predator - prey interactions

Introduction

In the Lotka-Voterra common resource competition exercise, we saw that the dynamics of any given species can be influenced by those of other species which are competing with it for resources. This often takes place within a particular trophic level. Two species of ungulates could feed on grasslands within a common area, for instance. In nature, population dynamics are also regulated by species across trophic levels – the dynamics of predator populations are influenced by how many prey individuals are in the community is a common example of this. These predator-prey interactions are also described by Lotka-Voterra equations. In this case these are a pair of differential equations describing the population dynamics of one species at one trophic level as a function of the dynamics of species within another other trophic level. Although it is most commonly discussed in a “predator-prey” context, the framework here can me modified to suit other scenarios such as plant-herbivore, parasites-host interactions.

Exponentially growing prey

The simplest Lotka-Volterra predator-prey equation considers a prey or victim species \(V\) that grows exponentially at a rate \(r\), and shrinks as it gets consumed by predators \(P\), which attack the prey they encounter at a fixed, per-capita (i.e. dependent on the number of predators and prey) rate \(\alpha\):

\[\frac{dV}{dt} = rV - \alpha VP\]

The population of the predators \(P\) grows as it converts the prey which it encountered and consumed into new predator individuals at a fixed per-captia rate \(\beta\), and shrinks due to a constant per-capita mortality rate \(d\):

\[\frac{dP}{dt} = \beta \alpha VP - dP\]

Let’s take a look at what this model implies. When there are no predators in the community (i.e. \(P = 0\)), the prey species \(V\) grows due to simple exponential growth (\(rV\)). If there were limiting resources then the prey would grow until its carrying capacity but we will ignore that aspect in this example for simplicity. Conversely, when there are no prey individuals in the community (i.e. \(V = 0\)), the predator popultion has no way to sustain its growth grow and exponentially declines towards zero \(-dP\).

The rate of growth of the predator population increases when prey are abundant and predators are few. On the other hand, when there are plenty of predators, the prey population suffers and cannot grow very quickly. We can get some suggestion that this predator-prey relationship might lead to some sort of population cycling – as the prey population grows, so too does the predator popuplation, which in turn reduced the number of prey which in turn reduces the number of predators allowing the prey to increase once more.

As expected, cyclic population dynamics are common in predator-prey systems:

In the first plot above, the system begins with 25 prey and 10 predators – cyclic dynamics shortly ensue. In the middle panel, there are no predators in the beginning, so the prey population simply grows exponentially. In the third panel, there are no prey individuals in the system – so, the predator population shrinks exponentially to zero.

Logistically growing prey

In the previous model, prey experience exponential growth when predators are absent – in other words, the predators are the only factor regulating the prey population’s dynamics. Of course, this doesn’t generally happen in nature. Even when predators are absent a prey population will eventually level off at some carrying capacity dictated by the available resources within the habitat.

To incorporate this aspect into our equations we can modify the above prey equation by incorporating density-dependent growth simply by adding a term to represent the population’s carrying capacity:

\[\frac{dV}{dt} = rV - \alpha VP - cV^2\]

Where the value of the carrying capacity can be calculated as \(K = \frac{r}{c}\) which has been taken from logistic population growth. So if you want a very large \(K\) then make \(c\) very small. We keep the predator population dynamics the same as before:
\[\frac{dP}{dt} = baNP - dP\]

With these new equations, when predators are absent from the community, the prey species grows exponentially but it is then restricted by its own abundance in a similar manner to how the logistic growth equation would function. One possible outcome in this case is stabilised population cycles - in other words, the predator and prey populations can cycle up and down until they reach a joint equilibrium:

Note that in the middle graph, the prey population grows logistically to its carrying capacity of 10 when no predators are present. If we run this system out longer, we see that the populations eventually stop cycling as they reach an equilibrium when both predator and prey are present in the system:

Parameter exploration

Download Lotka-Voltera predator prey app here and answer the questions on the prac worksheet.

Worksheet questions (30)

Coding questions (20)

We will be introducing a new aspect into our code in this practical. Because we have two separate equations to model it is easier for us to provide the computer with our rate of change equations and then for the computer to integrate these equations and calculate the population sizes automatically.

You may remember this being mentioned in earlier practicals but if you don’t then you’ll be glad to know that R has a fantastic package called deSolve (differential equation Solver) which is just right for the job! This package is very powerful but requires a bit of a learning curve. Understanding everything about this package is beyond the scope of these practicals but with what you have learned from the previous practicals together with a short tutorial on the package you should be able to obtain basic outputs from the package.

There are three main aspects needed for us to use deSolve successfully.

1. First we need to define a function() which contains all of the aspects relating to the differential equations. In the code below this function is called epg because the model it is working with is the exponential population growth model. This function takes three arguments; time, init, and parameters. The body of the code then contains with() which tells the next bit of code make use of time, init and parameters in the next code chunk. This next bit of code is where the differential equations are stored (in our case the equation is stored as dNdt. Each equation gets stored in an appropriately named variable and then the final call in this function is to return() these equations. Although this function is defined at the beginning of our script it does not actually run until we call epg() and supply arguments for it. All we are doing is just defining a function just as we would any other variable.

2. The next aspects of our code are to define the length and values of the predictor series (in our case this is time as a sequence of numbers beginning at 0, ending at 1000, and increasing by 1), the initial values of the population(s) we are dealing with and the parameters associated with each equation. These are the vectors that will be combined with with() and fed into the differential equations.

3. Our final next step is to generate out, a variable which stores the output of the ode function. ode() (full function name: General Solver for Ordinary Differential Equations) is the function within the deSolve package which will carry out the calculations necessary to determine the changes in population sizes between the time steps. This function needs four parameters fed to it - all of which we have already:

   • The first argument is y - the initial values for each of the differential equations - we would have stored these in init
     
   • Next we need times, the time values we want to use to integrate these differnetials at - we have these values stored in times
     
   • Then we need func which is the function that contains the equations - we have developed this function and called it equ()
     
   • And finally we need parms, the parameters which are fed into func which we have called params

Once we have all those inputted to ode and we have run it successfully we can simply convert the output to a data frame (using as.data.frame()) and then we can neatly plot our modified data frame using ggplot(). See the code below as an example - the three components listed as required above are also explicitly stated in the code required to produce the figure below:

# load required packages 
list.of.packages <- c("deSolve", "dplyr", "ggplot2")

#checking missing packages from list
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]

#install missing ones
if(length(new.packages)) install.packages(new.packages, dependencies = TRUE)

# load packages from library
library(deSolve) 
library(dplyr)
library(ggplot2) 

# define function which accepts arguments for the population parameter vectors (time, init, params), 
# parses these to the equation using with() 
# and then returns the output
epg <- function(time, init, params) {
  with (as.list(c(time, init, params)), {
    
    dNdt <- r*N

    return(list(c(dNdt)))
  })
}

# define values of the parameter vectors 
time <- seq(from = 0, to = 30, by = 0.1)
init <- c("N" = 1)
params <- c("r" = 0.2)

# populate a matrix called out in our case using ode() and supplying all the required arguments
out <- ode(y = init, 
           times = time, 
           func = epg, 
           parms = params)

# convert out to a data frame using dplyr
out <- out %>% as.data.frame() 

# plot the out dataframe using ggplot2
ggplot(out) + 
  geom_line(aes(x = time, 
                y = N), 
            size = 1.25, 
            alpha = 0.6) + 
  labs(title = "Exponentially growing population", 
       subtitle = "Parameters were as follows: N0 = 1, r = 0.2", 
       x = "Time", 
       y = "Population size") + 
  theme_bw() + 
  theme(axis.title = element_text(colour = "black", 
                                      size = 12), 
            axis.text = element_text(colour = "black", 
                                     size = 11)) 

Application questions (38)