1. The Lotka-Volterra competition equations (6a and 6b in your text) are described: \[ \frac{\mathrm{d}N_1}{\mathrm{d}t} = r_1N_1\left(\frac{K_1 - N_1 - \alpha _{12}N_2}{K_1}\right) \\ \frac{\mathrm{d}N_2}{\mathrm{d}t} = r_2N_2\left(\frac{K_2 - N_2 - \alpha _{21}N_1}{K_2}\right). \] This is similar to what you’ve seen with mutualism and predator-prey, but with both of the interaction terms having negative signs. For these equations, first algebraically find the nullclines. Next, plot them all (not using phaseR), for the 4 different outcomes of competition; i.e., remake figure 8.4.

The nullclines are \(N_i^* = K_i - \alpha_{ij}N_j\).

# If you want, you can write functions to make code cleaner and clearer
comp_plot <- function(K1, K2, alpha12, alpha21) {
  plot(x = NA, type = "n", xlab = "N1", ylab = "N2", xlim = c(0, 1), ylim = c(0, 1), las = 1)
  segments(x0 = 0, y0 = 0, x1 = 0, y1 = 1, col = "red")
  segments(x0 = 0, y0 = 0, x1 = 1, y1 = 0, col = "blue")
  N1 <- seq(from = 0, to = 1, by = 0.1)
  N2_null <- K2 - alpha21*N1
  N1_null <- K1/alpha12 - 1/alpha12*N1
  lines(x = N1, y = N1_null, col = "red")
  lines(x = N1, y = N2_null, col = "blue")
}

par(mfrow = c(2, 2), mar = c(4.5, 4.5, 0.5, 0.5))
comp_plot(K1 = 0.9, K2 = 0.9, alpha12 = 0.25, alpha21 = 0.25)
comp_plot(K1 = 0.9, K2 = 0.9, alpha12 = 1.25, alpha21 = 0.25)
comp_plot(K1 = 0.9, K2 = 0.9, alpha12 = 0.25, alpha21 = 1.25)
comp_plot(K1 = 0.9, K2 = 0.9, alpha12 = 1.5, alpha21 = 1.5)

  1. Recreate figure 8.4. using phaseR. For each of the outcomes, plot two trajectories using trajctory() (see help for details): with one set of initial conditions with a larger \(N_1\) and a smaller \(N_2\) and one set with a smaller \(N_1\) and larger \(N_2\).
library(package = "deSolve")
library(package = "phaseR")
## -------------------------------------------------------------------------------
## phaseR: Phase plane analysis of one- and two-dimensional autonomous ODE systems
## -------------------------------------------------------------------------------
## 
## v.2.1: For an overview of the package's functionality enter: ?phaseR
## 
## For news on the latest updates enter: news(package = "phaseR")
comp <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
  dN1 <- r1*N1*(K1 - N1 - alpha12*N2)/K1
  dN2 <- r2*N2*(K2 - N2 - alpha21*N1)/K2
  return(list(c(dN1, dN2)))
})
}

time_vec <- seq(from = 0, to = 100, by = 0.01)

parm_vec_a <- parm_vec_b <- parm_vec_c <- parm_vec_d <- c(r1 = 1, r2 = 1, K1 = 0.9, K2 = 0.9, alpha12 = 0.25, alpha21 = 0.25)
parm_vec_b["alpha12"] <- 1.25
parm_vec_c["alpha21"] <- 1.25
parm_vec_d[c("alpha12", "alpha21")] <- 1.5

N01 <- c(N1 = 0.8, N2 = 0.1)
N02 <- c(N1 = 0.1, N2 = 0.8)

lims <- c(-0.1, 1)

par(mfrow = c(2, 2), mar = c(4.5, 4.5, 0.5, 0.5))
trash <- nullclines(deriv = comp, parameters = parm_vec_a, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = F, state.names = c("N1", "N2"), las = 1, add.legend = F, xlab = "N1", ylab = "N2")
trash <- trajectory(deriv = comp, y0 = N01, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_a, state.names = c("N1", "N2"))
trash <- trajectory(deriv = comp, y0 = N02, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_a, state.names = c("N1", "N2"))
points(x = c(N01[1], N02[1]), y = c(N01[2], N02[2]), pch = 19)

trash <- nullclines(deriv = comp, parameters = parm_vec_b, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = F, state.names = c("N1", "N2"), las = 1, add.legend = F, xlab = "N1", ylab = "N2")
trash <- trajectory(deriv = comp, y0 = N01, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_b, state.names = c("N1", "N2"))
trash <- trajectory(deriv = comp, y0 = N02, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_b, state.names = c("N1", "N2"))
points(x = c(N01[1], N02[1]), y = c(N01[2], N02[2]), pch = 19)

trash <- nullclines(deriv = comp, parameters = parm_vec_c, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = F, state.names = c("N1", "N2"), las = 1, add.legend = F, xlab = "N1", ylab = "N2")
trash <- trajectory(deriv = comp, y0 = N01, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_c, state.names = c("N1", "N2"))
trash <- trajectory(deriv = comp, y0 = N02, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_c, state.names = c("N1", "N2"))
points(x = c(N01[1], N02[1]), y = c(N01[2], N02[2]), pch = 19)

trash <- nullclines(deriv = comp, parameters = parm_vec_d, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = F, state.names = c("N1", "N2"), las = 1, add.legend = F, xlab = "N1", ylab = "N2")
trash <- trajectory(deriv = comp, y0 = N01, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_d, state.names = c("N1", "N2"))
trash <- trajectory(deriv = comp, y0 = N02, tlim = c(0, 100), tstep = 0.1, parameters = parm_vec_d, state.names = c("N1", "N2"))
points(x = c(N01[1], N02[1]), y = c(N01[2], N02[2]), pch = 19)