Predator-prey (consumer-resource) interactions

class: center, middle, inverse, title-slide

.title[
# Predator-prey (consumer-resource) interactions
]
.author[
### Christopher Moore
]
.date[
### Week 06
]

---

class: inverse, center, middle

![Grey Grey Owl](https://cdn.thinglink.me/api/image/690547197309091840/1240/10/scaletowidth)

# Predator-prey (consumer-resource) interactions, example I

---
class: inverse, center, middle

![Oak](https://acornsonglendotcom.files.wordpress.com/2011/09/oak.jpg)

---
# Lotka and Volterra
- Alfred J. Lotka, American, first proposed predator-prey model in 1910 but more generally in 1925

- Vito Volterra, Italian, inspired by fisheries data, published 1926

*Lotka-Volterra* is used to describe (nearly) any two-species models, generally; but when used alone it refers to their independently-derived predator-prey model

---
# General two-species model

`$$\begin{aligned}
\frac{dR}{dt} &= f_1(R, C) \\
\frac{dC}{dt} &= f_2(R, C)
\end{aligned}$$`

Predators and prey (consumers and resources) also symbolically represented respectively as `$V$` and `$P$` (in book; victim and predator), `$N_1$` and `$N_2$`, and `$H$` and `$P$` among others.

---
# **T**he Lotka-Volterra model

We *minimally* have four componenets to add to a model of predator and prey dynamics:

1. Prey needs to grow
2. Prey get eaten
3. Predators grow by eating prey
4. Predators must die

Let's represent these processes using the following
--

1. `$+rR$`
--

2. `$-aRC$`
--

3. `$+bRC$`
--

4. `$-mC$`

---
# **T**he Lotka-Volterra model
Putting the pieces together:
1. Prey growth, `$+rR$`
2. Predators removing prey, `$-aRC$`
3. Predator growth from eating prey, `$+bRC$`
4. Predators die, `$-mC$`
--

`$$\begin{aligned}
\frac{dR}{dt} &= rR - aRC \\
\frac{dC}{dt} &= bRC - mC
\end{aligned}$$`

--
. . . before we analyze this model, let's consider a simpler, more familiar case between two species

---
# Neutralism
The case where two species do not affect one another.

Let's write down the coupled equations for two species, but without interaction, behaving as single species do (logistically). Write each species as `$N_1$` and `$N_2$`.

`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1 - \alpha _1N_1^2 \\
\frac{dN_2}{dt} &= r_2N_2 - \alpha _2 N_2^2
\end{aligned}$$`

`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1\left(\frac{K_1 - N_1}{K_1}\right) \\
\frac{dN_2}{dt} &= r_2N_2\left(\frac{K_2 - N_2}{K_2}\right)
\end{aligned}$$`

---
# Neutralism time series: code

```r
neutralism <- function(t, y, parameters) {
  N1 <- y[1]
  N2 <- y[2]
  r1 <- parameters[1]
  r2 <- parameters[2]
  a1 <- parameters[3]
  a2 <- parameters[4]
  dN1 <- r1*N1 - a1*N1^2
  dN2 <- r2*N2 - a2*N2^2
  return(list(c(dN1, dN2)))
}

n.parms <- c(r1 = 0.55, r2 = 0.5, a1 = 0.02, a2 = 0.01)

times <- seq(from = 0, to = 20, by = 0.1)
init <- c(N1 = 2, N2 = 0.5)
neut.out <- ode(y = init, times = times, parms = n.parms, func = neutralism)
```

---
# Neutralism time series: plot

```r
plot(x = neut.out[,1], y= neut.out[,2], lwd = 2, las = 1, xlab = "Time", ylab = "Density", type = "l", ylim = c(0, max(neut.out)), col = "#0000FFFF")
lines(x = neut.out[,1], y= neut.out[,3], lwd = 2, col = "red")
```

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-2-1.png)

---
# Neutralism in state space

```r
plot(x = neut.out[,2], y= neut.out[,3], lwd = 2, las = 1, xlab = bquote(N[1]), ylab = bquote(N[2]), type = "p", ylim = c(0, max(neut.out)), xlim = c(0, max(neut.out)), col = cp, pch = 16)
```

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-4-1.png)

*dark and light colors respectively correspond to early and later points in time*

---
# Neutralism: animated
I. Both `$N_1$` and `$N_2$` below equilibirum

---
# Neutralism: animated
II. `$N_1$` below and `$N_2$` above their equilibiria

---
# Neutralism: animated
III. `$N_1$` above and `$N_2$` below their equilibiria

---
# Neutralism: animated
IV. `$N_1$` and `$N_2$` above their equilibiria

---
# Neutralism: animated
All four starting points (initial conditions) showing different behaviours

---
# Nullclines: the trajectory separators
**Nullclines** describe the curves (including lines) on a phase plane towards or away from which a populations will grow. This is analagous to the concept of an equilibirum, but on a plane (2 dimensions) populations will grow towards or away from a line instead of a point.

If we were to plot a logistic-growing population against something arbitrary, like the number of angels dancing on the head of a pin, then we could create a plane that looks something like this:

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-6-1.png)

If we were to plot a logistic-growing population against something arbitrary, like the number of angels dancing on the head of a pin, then we could create a plane that looks something like this:

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-7-1.png)

If we were to plot a logistic-growing population against something arbitrary, like the number of angels dancing on the head of a pin, then we could create a plane that looks something like this:

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-8-1.png)

If we were to plot a logistic-growing population against something arbitrary, like the number of angels dancing on the head of a pin, then we could create a plane that looks something like this, **or like this**:

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-9-1.png)

---
# Nullclines: the trajectory separators
Now, if we had **two** populations/variables (say, `$N_1$` and `$N_2$`), neither of which are interacting with one another, what would we expect to see?

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-10-1.png)

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-11-1.png)

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-12-1.png)

---
# Nullclines: calculation

Set `$\frac{dN_i}{dt} = 0$` for each equation, put into a graphical form (isolate the variable on the vertical axis), plot, . . . and crush it.

![CRUSH](https://c.tenor.com/s7hjKTx1nJYAAAAC/stranger-things-el.gif)

Find the nullclines for the neutralism model . . . and plot them (one for each species).

---
# Nullclines: calculation

The model
`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1 - \alpha _1N_1^2 \\
\frac{dN_2}{dt} &= r_2N_2 - \alpha _2N_2^2
\end{aligned}$$`

Setting each subequation to 0
`$$\begin{aligned}
0 &= r_1N_1 - \alpha _1N_1^2 \\
0 &= r_2N_2 - \alpha _2N_2^2
\end{aligned}$$`

.pull-left[
For `$N_1$`
`$$\begin{aligned}
N^*_1 &= 0 \\
N^*_1 &= \frac{r_1}{\alpha _1}
\end{aligned}$$`
]

.pull-right[
For `$N_2$`
`$$\begin{aligned}
N^*_2 &= 0 \\
N^*_2 &= \frac{r_2}{\alpha _2}
\end{aligned}$$`
]

---
# Nullclines: plot your results
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-13-1.png)

---
# Nullclines: plot your results
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-14-1.png)

---
# Neutralism: using R to graph dynamics
We will use a package called `phaseR` that is found in the **C**omprehensive **R** **A**rchive **N**etwork (CRAN).

You can download pacakges using the function `install.packages(pkgs = )`. Let's `install.packages(pkgs = "phaseR")`.

Now that it is installed and saved to your personal R library, let's load it: `library(package = "phaseR")`.

---
# Neutralism: using R to graph dynamics

We have to call (1) `parms`, `parameters` **and** (2) `state`, `y` in our function to play with `phaseR` (it's terrible, I know; I am sincerely sorry)

init <- c(N1 = 0.1, N2 = 0.05)
parm.vals <- c(r1 = 2/3, r2 = 2/3, a1 = 1, a2 = 1)
```
Then, let's plot the nullclines:

```r
nullclines(deriv = neutralism, parameters = parm.vals, xlim = c(-0.1, 1), ylim = c(-0.1, 1), add = F)
```

---
# Neutralism: using R to graph dynamics
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-17-1.png)

---
# Neutralism: using R to graph dynamics

```r
flowField(deriv = neutralism, xlim = c(-0.1, 1), ylim = c(-0.1, 1), parameters = parm.vals, add = F, points = 20)
```

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-19-1.png)

---
# Neutralism: using R to graph dynamics
Plotting nullclines and direction field together
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-20-1.png)

---
# Neutralism: using R to graph dynamics
And if you want, you can add trajectories:

```r
trajectory(deriv = neutralism, y0 = init, t.start = 0, t.end = 25, pch = 16)
```

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-22-1.png)

---
# A new model: amensalism
**Amensalism**: when one species has a negative effect on another, but that other has no effect on it. A couple biological examples include allelopathy or strong competition.

What might the model look like?

The model
`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1 - \alpha _1N_1^2 - \beta _1N_1N_2\\
\frac{dN_2}{dt} &= r_2N_2 - \alpha _2N_2^2
\end{aligned}$$`

Create a phase plane for an amensal relationship between two species.

---
# Amensalism
The model, initial densities, and parameter values

```r
amensalism <- function(t, y, parameters) {
  N1 <- y[1]
  N2 <- y[2]
  r1 <- parameters[1]
  r2 <- parameters[2]
  a1 <- parameters[3]
  a2 <- parameters[4]
  b1 <- parameters[5]
  dN1 <- r1*N1 - a1*N1^2 - b1*N1*N2
  dN2 <- r2*N2 - a2*N2^2
  return(list(c(dN1, dN2)))
}

init <- c(N1 = 0.1, N2 = 0.05)
parm.vals <- c(r1 = 2/3, r2 = 2/3, a1 = 1, a2 = 1, b1 = 0.2)
```

---
# Amensalism
The phase plane

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-24-1.png)

---
# Amensalism
The phase plane

.pull-left[
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-25-1.png)
]

.pull.right[
Does this make sense?

What are the equations for the nullclines?

]

---
# Amensalism
The model
`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1 - \alpha _1N_1^2 - \beta _1N_1N_2\\
\frac{dN_2}{dt} &= r_2N_2 - \alpha _2N_2^2
\end{aligned}$$`

Setting each subequation to 0
`$$\begin{aligned}
0 &= r_1N_1 - \alpha _1N_1^2 - \beta _1N_1N_2 \\
0 &= r_2N_2 - \alpha _2N_2^2
\end{aligned}$$`

.pull-left[
For `$N_1$`
`$$\begin{aligned}
N^*_1 &= 0 \\
N^*_1 &= \frac{r_1 - \beta _1N_2}{\alpha _1}
\end{aligned}$$`
]

.pull-right[
For `$N_2$`
`$$\begin{aligned}
N^*_2 &= 0 \\
N^*_2 &= \frac{r_2}{\alpha _2}
\end{aligned}$$`
]

---
# Amensalism: varying `$\beta$`

---
# Consumer-resource, both logistic
The model with both species growing logistically, but with one species benefiting at the other's expense (say, a consumer-resource interaction)

`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1 - \alpha _1N_1^2 - \beta _1N_1N_2\\
\frac{dN_2}{dt} &= r_2N_2 - \alpha _2N_2^2 + \beta_21N_1N_2
\end{aligned}$$`

---
# Consumer-resource, both logistic
The model
`$$\begin{aligned}
\frac{dN_1}{dt} &= r_1N_1 - \alpha _1N_1^2 - \beta _1N_1N_2\\
\frac{dN_2}{dt} &= r_2N_2 - \alpha _2N_2^2 + \beta_21N_1N_2
\end{aligned}$$`

Setting each subequation to 0
`$$\begin{aligned}
0 &= r_1N_1 - \alpha _1N_1^2 - \beta _1N_1N_2 \\
0 &= r_2N_2 - \alpha _2N_2^2 + \beta_21N_1N_2
\end{aligned}$$`

.pull-left[
For `$N_1$`
`$$\begin{aligned}
N^*_1 &= 0 \\
N^*_1 &= \frac{r_1 - \beta _1N_2}{\alpha _1}
\end{aligned}$$`
]

.pull-right[
For `$N_2$`
`$$\begin{aligned}
N^*_2 &= 0 \\
N^*_2 &= \frac{r_2 + \beta _2N_1}{\alpha _2}
\end{aligned}$$`
]

---
# Consumer-resource, both logistic
If we were to plot these on a graph, with red lines being the nullclines for `$N_1$` and the blue lines being the nullclines for `$N_2$`, then we'd draw something like this:

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-27-1.png)

---
# Consumer-resource, both logistic
The model, initial densities, and parameter values

```r
cr.log <- function(t, y, parameters) {
  N1 <- y[1]
  N2 <- y[2]
  r1 <- parameters[1]
  r2 <- parameters[2]
  a1 <- parameters[3]
  a2 <- parameters[4]
  b1 <- parameters[5]
  b2 <- parameters[6]
  dN1 <- r1*N1 - a1*N1^2 - b1*N1*N2
  dN2 <- r2*N2 - a2*N2^2 + b2*N2*N1
  return(list(c(dN1, dN2)))
}

init <- c(N1 = 0.1, N2 = 0.05)
parm.vals <- c(r1 = 2/3, r2 = 2/3, a1 = 1, a2 = 1, b1 = 0.2, b2 = 0.2)
```

---
# Consumer-resource, both logistic
The phase plane

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-29-1.png)

---
# Consumer-resource, both logistic

---
# C-R/L-V, R ex. growth and C cons. mort
This is the Lotka-Volterra predator-prey model:

`$$\begin{aligned}
\frac{dR}{dt} &= rR - \beta RC\\
\frac{dC}{dt} &= \gamma RC - \delta C
\end{aligned}$$`

---
# C-R/L-V, R ex. growth and C cons. mort
The model
`$$\begin{aligned}
F_1 = \frac{dR}{dt} &= rR - \beta RC\\
F_2 = \frac{dC}{dt} &= \gamma RC - \delta C
\end{aligned}$$`

Setting each subequation to 0
`$$\begin{aligned}
F_1 = 0 &= rR - \beta RC\\
F_2 = 0 &= \gamma RC - \delta C
\end{aligned}$$`

.pull-left[
For `$F_1$`
`$$\begin{aligned}
R^* &= 0 \\
C^* &= \frac{r}{\beta}
\end{aligned}$$`
]

.pull-right[
For `$F_2$`
`$$\begin{aligned}
C^* &= 0 \\
R^* &= \frac{\delta}{\gamma}
\end{aligned}$$`
]

---
# C-R/L-V, R ex. growth and C cons. mort
If we were to plot these on a graph, with red lines being the nullclines for `$R$` and the blue lines being the nullclines for `$C$`, then we'd draw something like this:

.pull-left[
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-30-1.png)
]

.pull-right[
![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-31-1.png)
]

---
# C-R/L-V, R ex. growth and C cons. mort
The model, initial densities, and parameter values

```r
cr.LV <- function(t, y, parameters) {
  R <- y[1]
  C <- y[2]
  r <- parameters[1]
  a <- parameters[2]
  b <- parameters[3]
  m <- parameters[4]
  dR <- r*R - a*R*C
  dC <- b*R*C - m*C
  return(list(c(dR, dC)))
}

init <- c(0.3, 0.5)
parm.vals <- c(r = 0.5, a = 0.5, b = 0.25, m = 0.5)
```

---
# C-R/L-V, R ex. growth and C cons. mort
The phase plane

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-33-1.png)

---
# C-R/L-V, R ex. growth and C cons. mort

---
class: inverse, center
background-image: url(http://www2.nau.edu/lrm22/lessons/predator_prey/hare_lynx_graph.png)
background-size: contain

# Snowshoe hare-lynx data from the Hudson Bay Trapping Company

---
# C-R, R log. growth and C cons. mort.
This is the Lotka-Volterra predator-prey model with generalist prey:

`$$\begin{aligned}
\frac{dR}{dt} &= rR - \alpha R^2 - \beta RC\\
\frac{dC}{dt} &= \gamma RC - \delta C
\end{aligned}$$`

---
# C-R, R log. growth and C cons. mort.
The model
`$$\begin{aligned}
F_1 = \frac{dR}{dt} &= rR - \alpha R^2 - \beta RC\\
F_2 = \frac{dC}{dt} &= \gamma RC - \delta C
\end{aligned}$$`

Setting each subequation to 0
`$$\begin{aligned}
F_1 = 0 &= rR - \alpha R^2 - \beta RC\\
F_2 = 0 &= \gamma RC - \delta C
\end{aligned}$$`

.pull-left[
For `$F_1$`
`$$\begin{aligned}
R &= 0 \\
R &= \frac{r - \beta C}{\alpha}
\end{aligned}$$`
]

.pull-right[
For `$F_2$`
`$$\begin{aligned}
C &= 0 \\
R &= \frac{\delta}{\gamma}
\end{aligned}$$`
]

---
# C-R, R log. growth and C cons. mort.
The model, initial densities, and parameter values

```r
cr.LVDD <- function(t, y, parameters) {
  R <- y[1]
  C <- y[2]
  r <- parameters[1]
  alpha <- parameters[2]
  a <- parameters[3]
  b <- parameters[4]
  m <- parameters[5]
  dR <- r*R - alpha*R*R - a*R*C
  dC <- b*R*C - m*C
  return(list(c(dR, dC)))
}

init <- c(0.3, 0.5)
parm.vals <- c(r = 0.5, alpha = 0.2, a = 0.25, b = 0.25, m = 0.4)
```

---
# C-R, R log. growth and C cons. mort.
The phase plane

![](Wk06_PredatorPrey_files/figure-html/unnamed-chunk-35-1.png)

---
# C-R, R log. growth and C cons. mort.