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). \] These are the same equations that you studied in Ecology (BI/ES271). For these equations, (a) first, show your work and algebraically find the nullclines. Next, (b) plot the four nullclines all (not using phaseR), for the 4 different outcomes of competition; i.e., remake figure 8.4. Use segments() or lines(), but don’t use nullclines().
  1. For the first equation: \[ \begin{aligned} 0 &= \frac{\mathrm{d}N_1}{\mathrm{d}t} = r_1N_1\left(\frac{K_1 - N_1 - \alpha _{12}N_2}{K_1}\right) \\ N_1^* &= 0\text{ is one nullcline} \\ 0 &= r_1\left(\frac{K_1 - N_1 - \alpha _{12}N_2}{K_1}\right)\text{is the other nullcline}\\ &\text{I'm going to rearrange the eqn. above to more easily plot it:} \\ 0 &= r_1\left(\frac{K_1 - N_1}{K_1}\right) - r_1\frac{\alpha _{12}N_2}{K_1}\\ r_1\frac{\alpha _{12}N_2}{K_1} &= r_1\left(\frac{K_1 - N_1}{K_1}\right)\\ \frac{\alpha _{12}N_2}{K_1} &= \frac{K_1 - N_1}{K_1}\\ \alpha _{12}N_2 &= K_1 - N_1 \\ N_2 &= \frac{K_1}{\alpha _{12}} - \frac{1}{\alpha _{12}}N_1^* \\ &\text{I wrote the above in this form so it's } y = mx + b \end{aligned} \] For the second equation: \[ \begin{aligned} 0 &= \frac{\mathrm{d}N_2}{\mathrm{d}t} = r_2N_2\left(\frac{K_2 - N_2 - \alpha _{21}N_1}{K_2}\right) \\ N_2^* &= 0\text{ is one nullcline} \\ 0 &= r_2\left(\frac{K_2 - N_2 - \alpha _{21}N_1}{K_2}\right)\text{is the other nullcline}\\ &\text{I'm going to rearrange the eqn. above to more easily plot it:} \\ 0 &= r_2\left(\frac{K_2 - \alpha _{21}N_1}{K_2}\right) - r_2\frac{N_2}{K_2}\\ r_2\frac{N_2}{K_2} &= r_2\left(\frac{K_2 - \alpha _{21}N_1}{K_2}\right)\\ \frac{N_2}{K_2} &= \frac{K_2 - \alpha _{21}N_1}{K_2} \\ N_2^* &= K_2 - \alpha _{21}N_1 \\ \end{aligned} \]
# 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", lwd = 2)
 segments(x0 = 0, y0 = 0, x1 = 1, y1 = 0, col = "blue", lwd = 2)
 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", lwd = 2)
 lines(x = N1, y = N2_null, col = "blue", lwd = 2)
}

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’s functions flowField() and nullclines(). 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\).
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.05, 1)

par(mfrow = c(2, 2), mar = c(4.5, 4.5, 0.5, 0.5))
plot(x = NA, type = "n", xlim = lims, ylim = lims, xlab = expression("N"[1]), ylab = expression("N"[2]), las = 1)
trash <- flowField(deriv = comp, parameters = parm_vec_a, xlim = lims, ylim = lims, add = T, state.names = c("N1", "N2"), points = 22, arrow.type = "proportional")
trash <- nullclines(deriv = comp, parameters = parm_vec_a, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = T, state.names = c("N1", "N2"), add.legend = F, lwd = 2)
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)

plot(x = NA, type = "n", xlim = lims, ylim = lims, xlab = expression("N"[1]), ylab = expression("N"[2]), las = 1)
trash <- flowField(deriv = comp, parameters = parm_vec_b, xlim = lims, ylim = lims, add = T, state.names = c("N1", "N2"), points = 22, arrow.type = "proportional")
trash <- nullclines(deriv = comp, parameters = parm_vec_b, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = T, state.names = c("N1", "N2"), add.legend = F, lwd = 2)
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)

plot(x = NA, type = "n", xlim = lims, ylim = lims, xlab = expression("N"[1]), ylab = expression("N"[2]), las = 1)
trash <- flowField(deriv = comp, parameters = parm_vec_c, xlim = lims, ylim = lims, add = T, state.names = c("N1", "N2"), points = 22, arrow.type = "proportional")
trash <- nullclines(deriv = comp, parameters = parm_vec_c, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = T, state.names = c("N1", "N2"), add.legend = F, lwd = 2)
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)

plot(x = NA, type = "n", xlim = lims, ylim = lims, xlab = expression("N"[1]), ylab = expression("N"[2]), las = 1)
trash <- flowField(deriv = comp, parameters = parm_vec_d, xlim = lims, ylim = lims, add = T, state.names = c("N1", "N2"), points = 22, arrow.type = "proportional")
trash <- nullclines(deriv = comp, parameters = parm_vec_d, xlim = lims, ylim = lims, points = 201, col = c("red", "blue"), add = T, state.names = c("N1", "N2"), add.legend = F, lwd = 2)
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)