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.

The 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$

The 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}$$

or

$$\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

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

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")

Neutralism in state space

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)

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:

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:

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:

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, or like this:

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?

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, more simply?

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, more simply?

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.

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}$$

Nullclines: plot your results

Nullclines: plot your results

Neutralism: using R to graph dynamics

We will use a package called phaseR that is found in the Comprehensive R Archive Network (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)

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)))
}
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:

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

Neutralism: using R to graph dynamics

Neutralism: using R to graph dynamics

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

Neutralism: using R to graph dynamics

Plotting nullclines and direction field together

Neutralism: using R to graph dynamics

And if you want, you can add trajectories:

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

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

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

Amensalism

The phase plane

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}$$

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}$$

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:

Consumer-resource, both logistic

The model, initial densities, and parameter values

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

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}$$

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:

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

The model, initial densities, and parameter values

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

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

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}$$

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

The model, initial densities, and parameter values

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

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