Throughout this lab, we will model predator/prey interactions in increasingly complex and realistic ways. We will start with simple functional responses and move to Lotka-Volterra models and Rosenzweig-MacArthur models.
The term \(aN\) describes the functional response of the predator species, the rate at which a prey kills/consumes its prey (population size \(N\)) at an attack rate \(a\).
So \(aN\) can be used to describe the functional response of a single individual predator. We can plot the number of prey available by the number of prey killed per unit time to visualize the functional response.
library(ggplot2)Warning: package 'ggplot2' was built under R version 4.3.2a <- 0.1
df <- data.frame(N = c(0, 50))
ggplot(data = df, aes(x = N)) +
stat_function(fun = function(N) a * N) +
labs(x = "Prey Density", y = "Prey Killed Per Predator")
Following this example, the predator always consume a fixed proportion of prey no matter how many prey the predator has already consumed. This curve represents a Type I functional response.
In reality, predator consumption will often reach a saturation point when the predator has eaten enough prey.
Type II responses are more frequently seen in nature:
\[ \frac{aN}{1 + ahN} \] What are each of the terms in this equation?
When are Type II responses seen in nature?
Now plot a Type II functional response:
a <- 0.1
h <- 50
df <- data.frame(N = c(0, 1))
ggplot(data = df, aes(x = N)) +
stat_function(fun = function(N) a * N / (1 + a * h * N)) +
labs(x = "Prey Density", y = "Prey Killed Per Predator")
In contrast, when do we typically see Type III responses?
\[ \frac{\frac{1}{h}N^{2}}{(\frac{1}{ah})^{2} + N^{2}} = \frac{cN^{2}}{d^{2} + N^{2}} \]
Where \(c = \frac{1}{h}\) is the maximum consumption rate and \(d = \frac{1}{ah}\) is the half-saturation constant.
Plot a Type III functional response:
c <- 0.1
d <- 0.5
df <- data.frame(N = c(0, 1))
ggplot(data = df, aes(x = N)) +
stat_function(fun = function(N) c * N^2 / (d^2 + N^2)) +
labs(x = "Prey Density", y = "Prey Killed Per Predator")
As another way to represent functional response curves, represent the code below as a dataframe and plot it. What do you notice as you vary the values of the terms?
\[ F = \frac{bN^{1+q}}{1+bhN^{1+q}} \]
We find that all response curves can be modeled in the same equation:
If \(q>0\), we have a Type III response curve
If \(q = 0\), we have a Type II response curve
if \(q=0\) and \(h = 0\), then we have a Type I response curve
Explore the agent-based model below:
Compared to other simulation types that attempt to simulate systems as a whole, agent-based models are a type of simulation where the behavior of individual units (agents) are simulated, to observe the additive result of how individual behaviors and/or interactions might affect the system as a whole.
Select the model version sheep-wolves. Click setup to set-up the model based on your parameters. Use the go button to start and stop the simulation.
Change the model version to sheep-wolves-grass. Before running the simulation, what do you expect to happen?
Now run the simulation and modify the parameters to stabilize the system. What do you find?
In R, we can construct a simple Lotka-Volterra model to simulate single predator/prey interactions.
\[ \begin{array}{l} \frac{dN}{dt} \ =\ rN\ -\ aNP\\ \\ \frac{dP}{dt} \ =\ faNP\ -\ qP \end{array} \]
\[ \begin{array}{l} r\ =\ prey\ per\ capita\ growth\ rate\\ N\ =\ number\ of\ prey\ species\ individuals\\ a\ =\ per\ capita\ attack\ rate\\ P\ =\ number\ of\ predator\ species\ individuals\\ f\ =\ predator\ efficiency\\ q\ =\ predator\ per\ capita\ mortality\ rate \end{array} \]
Remember from last lab, \(r\) is generally used for continuously reproducing populations
Density independent, because it is just a growth rate and the population size like we saw last lab in our density independent formula:
\[ N_{t} \ =\ r N_{0} \]
Now let’s create a function for the Lotka-Volterra model in R. We can write differential equations in R using the ordinary differential equation (ODE) package deSolve(Soetaert, Petzoldt, and Setzer 2010) following the format:
function_name <- function(y, x, parameters) {
X <- x[1]
with(as.list(params), {
dX.dy <- X # complete formula here
return(list(c(dX.dy)))
})
}R’s ordinary differential equation (ODE) solver deSolve() can be used to analyze population growth. ODEs are commonly used to calculate rates of change in order to model dynamic systems.
If you have trouble installing, try the following:
install.packages("https://cran.r-project.org/src/contrib/deSolve_1.34.tar.gz", repos=NULL, type="source")If after restarting Rstudio/R that does not work, try the following:
library(devtools)
# Double check the most up to date version
install_version("deSolve", version = "1.34", repos = "http://cran.us.r-project.org")If neither option works and you are using a Mac, you can try the following in a Homebrew command line:
brew install wget
wget https://cran.r-project.org/src/contrib/deSolve_1.34.tar.gz
R CMD INSTALL deSolve_1.34.tar.gz::: callout-tip
Now you try to use this format to define an ODE function for the prey (dN.dt) and predator (dP.dt) Lotka-Volterra models named lotka_volterra
lotka_volterra <- function(t, x, params) {
N <- x[1]
P <- x[2]
with(as.list(params), {
dN.dt <- r * N - a * P * N
dP.dt <- f * a * P * N - q * P
return(list(c(dN.dt, dP.dt)))
})
}In the example below we call an example using this format package::function().
Namespaces are used to organize functions, data, and other objects within a package. The double colon operator allows us to explicitly access a namespace within a specified package.
We are implicitly accomplishing this same task whenever we run a package function or access a package dataset after loading that package using library() or require(). The double colon allows us to do this same tasks explicitly either to avoid conflicts with a function/dataset from another package with the same name (as above) or to be more explicit programmers.
Now let’s run our model and observe the outcome:
library(deSolve)Warning: package 'deSolve' was built under R version 4.3.2parms <- c(
r = 0.5,
a = 0.01,
f = 0.1,
q = 0.2
)
Time <- seq(0, 200, by = 0.1)
lotka_volterra.out <- deSolve::ode(c(N0 = 25, P0 = 5), times = Time, func = lotka_volterra, parms = parms)
df_lotka_volterra <- data.frame(Time, lotka_volterra.out)
ggplot(df_lotka_volterra, aes(x = Time)) +
geom_line(aes(y = P0), color = "red") +
geom_line(aes(y = N0), color = "black") +
labs(x = "Time", y = "Population Size")
Like in the agent-based models and the lynx-snowshoe hare cycle, we see a lag between predator and prey population peaks.
We also see prey population peaks much higher than predator population peaks. The model was constructed this way to reflect what is seen in nature.
What we are observing is called a “neutral cycle.” The populations oscillate indefinitely. The populations are unstable because population sizes never move towards a stable node.
The only way the population sizes converge to a node is through extinction which, as discussed in last lab, is an unstable equilibrium point.
Note that if you set the initial predator population size P0 to 0, the prey population grows exponentially in a density-independent manner.
Michael Rosenzweig (a graduate student at the time) and his adviser Robert MacArthur had issues with the applicability of the Lotka-Volterra in natural systems.
Especially the often inaccurate assumptions we have discusses in the past two labs: the Type I response of predators and the density-independence of prey populations. Because of these, Rosenzweig and MacArthur made the following changes to the model:
Predators would be unlikely to follow a Type I response due to:
Predators would become satiated and/or
Handling time (h) would physically limit the number of prey a predator could consume
In the absence of predators, prey would become self limiting (density-dependent)
These changes can be integrated into the predator portion of the Lotka-Volterra formula as:
\[ \frac{dP}{dt} \ =\ \frac{faNP}{1\ +\ ahN} \ -\ qP \]
As for the prey portion, we can add the handling time as follows:
\[ \frac{dN}{dt} \ =\ rN\ -\ \frac{aNP}{1\ +\ ahN} \]
We need to incorporate the carrying capacity, \(K\), into the model. From last we saw that we can represent density-dependent population growth as:
\[ \frac{dN}{dt} \ =\ rN\left( 1\ -\ \frac{N}{K}\right) \] Therefore, to incorporate into our model we would write:
\[ \frac{dN}{dt} \ =\ rN\left( 1\ -\ \frac{N}{K}\right) \ -\ \frac{aNP}{1\ +\ ahN} \]
rosenzweig_macarthur.
rosenzweig_macarthur <- function(t, x, parms) {
N <- x[1]
P <- x[2]
with(as.list(parms), {
dN.dt <- a * N * (1 - (N / K)) - (a * P * N)/(1 + a * h * N)
dP.dt <- ((f * a * P * N)/(1 + a * h * N)) - q * P
return(list(c(dN.dt, dP.dt)))
})
}deSolve::ode function like we did above with h <- 0.1 and K <- 1000 and plot with ggplot2.
parms <- c(
r = 0.5,
a = 0.01,
q = 0.2,
f = 0.1,
h = 0.1,
K = 1000
)
Time <- seq(0, 200, by = 0.1)
rosenzweig_macarthur.out <- deSolve::ode(c(N0 = 25, P0 = 5), times = Time, func = rosenzweig_macarthur, parms = parms)
df_rosenzweig_macarthur <- data.frame(Time, rosenzweig_macarthur.out)
ggplot(df_rosenzweig_macarthur, aes(x = Time)) +
geom_line(aes(y = P0), color = "red") +
geom_line(aes(y = N0), color = "black") +
labs(x = "Time", y = "Population Size")
Warning: package 'patchwork' was built under R version 4.3.2
Increasing carrying capacity increases the height of the prey population size each cycle before the predator interaction drives the prey population size downwards again. This destabilizes the system because the prey population is not trending towards a node. This concept is called the “paradox of enrichment.”
Decreasing the carrying capacity causes the prey population size to drop in a density-independent manner towards stabilization at K.
Decreasing the handling time causes prey population to stabilize towards a node. The predator is said to be “controlling” the prey population.
Increasing the handling time allows the prey population to “escape” the control of the predator. Note that this escape is self-limiting, the prey population size converges at K.
We see collapse of the entire system. The prey is driven to extinction which leads to the extinction of the predator.
f low, we see predator escapef at moderate level, we see predator controlf at high level destabilizes the systemf at a very high level drives to collapse (extinction)