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")
  1. Please use the MacArthur-Rosenzweig predator-prey model to create a time series where (1) prey exist at equilibrium, (2) prey and predators stably coexist, and (3) where prey and predators show unstable behavior. Choose your parameters arbitrarily, but make sure that your starting densities of predator and prey area positive (i.e., prey existing at equilibrium should come from predators going extinct and not because initial densities of predators were 0). Plot all three as a three-panel figure using: par(mfrow = c(1, 3), mar = c(4.5, 4.5, 0.5, 0.5)) (the first argument makes one row and three columns of graphs, the second sets the amount of white space respectively on the bottom, left, top, and right.)
MR <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    dVdt <- r*V*(K-V)/K - c*V*P/(g + V)
    dPdt <- b*V*P/(g + V) - m*P
    return(list(c(dVdt, dPdt)))
  })
}

parm_vec_eq <- c(r = 0.5, K = 15, c = 2, g = 10, b = 1.5, m = 1)
parm_vec_stab <- c(r = 0.5, K = 30, c = 2, g = 10, b = 1.5, m = 1)
parm_vec_unstab <- c(r = 0.5, K = 100, c = 2, g = 10, b = 1.5, m = 1)

time_vec <- seq(from = 0, to = 100, length.out = 1000)
N0 <- c("V" = 2, "P" = 1.5)

MRode_prey <- ode(y = N0, times = time_vec, func = MR, parms = parm_vec_eq)
MRode_stab <- ode(y = N0, times = time_vec, func = MR, parms = parm_vec_stab)
MRode_unstab <- ode(y = N0, times = time_vec, func = MR, parms = parm_vec_unstab)
y_max_prey <- max(MRode_prey[,2:3], na.rm = T)
y_max_stab <- max(MRode_stab[,2:3], na.rm = T)
y_max_unstab <- max(MRode_unstab[,2:3], na.rm = T)

par(mfrow = c(1, 3), mar = c(4.5, 4.5, 0.5, 0.5))
plot(x = MRode_prey[,1], y = MRode_prey[,2], type = "l", xlab = "Time", ylab = "Prey (blue) and predator (red) density", col = "red", ylim = c(0, y_max_prey), las = 1)
  lines(x = MRode_prey[,1], y = MRode_prey[,3], col = "blue")

plot(x = MRode_stab[,1], y = MRode_stab[,2], type = "l", xlab = "Time", ylab = "Prey (blue) and predator (red) density", col = "red", ylim = c(0, y_max_stab), las = 1)
  lines(x = MRode_stab[,1], y = MRode_stab[,3], col = "blue")

  plot(x = MRode_unstab[,1], y = MRode_unstab[,2], type = "l", xlab = "Time", ylab = "Prey (blue) and predator (red) density", col = "red", ylim = c(0, y_max_unstab), las = 1)
  lines(x = MRode_unstab[,1], y = MRode_unstab[,3], col = "blue")

  1. For each of the three scenarios you platted above, I would like for you to plot (1) the nullclines and (2) add the trajectory from above.
    As a bit of guidance, remember that for the equation for the nullcline you want to get it in a \(y = f(x)\) form. You then plot it as such (e.g., create a sequence of \(x\) and apply the function to get \(y\)). To add lines to a plot, use the function lines() where the arguments are the same as plot() (e.g., you give it an x =, a y =, col =, etc.) Make the prey nullclines blue and the predator nullclines red.
    To add trajectories, you want to do two things. First, add a point of what the initial population densities are using points() (takes teh same arguments as plot()). Second, use your numerical solution from question 1. and plot \(x\) against \(y\).
V_vec <- seq(from = 0, to = 100, length.out = 1000)

V_null_prey <- parm_vec_eq["r"]/parm_vec_eq["c"]*(parm_vec_eq["g"] + V_vec*(1 - (parm_vec_eq["g"]/parm_vec_eq["K"])) - 1/parm_vec_eq["K"]*V_vec^2)
P_null_prey <- parm_vec_eq["m"]*parm_vec_eq["g"]/(parm_vec_eq["b"] - parm_vec_eq["m"])

V_null_stab <- parm_vec_stab["r"]/parm_vec_stab["c"]*(parm_vec_stab["g"] + V_vec*(1 - (parm_vec_stab["g"]/parm_vec_stab["K"])) - 1/parm_vec_stab["K"]*V_vec^2)
P_null_stab <- parm_vec_stab["m"]*parm_vec_stab["g"]/(parm_vec_stab["b"] - parm_vec_stab["m"])

V_null_unstab <- parm_vec_unstab["r"]/parm_vec_unstab["c"]*(parm_vec_unstab["g"] + V_vec*(1 - (parm_vec_unstab["g"]/parm_vec_unstab["K"])) - 1/parm_vec_unstab["K"]*V_vec^2)
P_null_unstab <- parm_vec_unstab["m"]*parm_vec_unstab["g"]/(parm_vec_unstab["b"] - parm_vec_unstab["m"])

ylims <- c(0, 8)

par(mfrow = c(1, 3), mar = c(4.5, 4.5, 0.5, 0.5))

plot(x = V_vec, y = V_null_prey, type = "l", col = "blue", las = 1, xlab = "Prey density (V)", ylab = "Predator density (P)", ylim = ylims)
segments(col = "blue", x0 = 0, x1 = 0, y0 = 0, y1 = ylims[2])
segments(col = "red", x0 = P_null_prey, x1 = P_null_prey, y0 = 0, y1 = ylims[2])
segments(col = "red", x0 = 0, x1 = 100, y0 = 0, y1 = 0)

plot(x = V_vec, y = V_null_stab, type = "l", col = "blue", las = 1, xlab = "Prey density (V)", ylab = "Predator density (P)", ylim = ylims)
segments(col = "blue", x0 = 0, x1 = 0, y0 = 0, y1 = ylims[2])
segments(col = "red", x0 = P_null_prey, x1 = P_null_prey, y0 = 0, y1 = ylims[2])
segments(col = "red", x0 = 0, x1 = 100, y0 = 0, y1 = 0)

plot(x = V_vec, y = V_null_unstab, type = "l", col = "blue", las = 1, xlab = "Prey density (V)", ylab = "Predator density (P)", ylim = ylims)
segments(col = "blue", x0 = 0, x1 = 0, y0 = 0, y1 = ylims[2])
segments(col = "red", x0 = P_null_prey, x1 = P_null_prey, y0 = 0, y1 = ylims[2])
segments(col = "red", x0 = 0, x1 = 100, y0 = 0, y1 = 0)

  1. I next want to teach you how to create a “direction field,” like what I have shown in class (and what is in the discussion presentations) that effectively population trajectories on grid. To do so, you need to install a package in R called phaseR.
    Do this by, only one time ever, running the line of code in the console (not your script because you don’t need to download and install a package every time you run your code!): install.packages(pkgs = "phaseR"). Once you have done that, at the top of your script/markdown file, load the package from your library by writing the line: library(package = "phaseR"). Now, have a look at the function flowField() in the help section by typing into the console: ?flowField Notice the arguments: these are all the things you can modify! What you absolutely need to run the function are:
    • deriv =: this is where you put your function
    • xlim = and ylim =: your two-element vector of min and max values for both axes
    • parameters =: your vector of parameters that you create when you pass them to ode()
    • state.names =: a vector with two character elements identifying the names of your variables
    • add = F: add this because the default is for this function to add arrows onto an existing plot. Setting this argument to FALSE means it creates a new plot
    Just like the others above, make 3 graphs. Plot the direction field using flowField(), then add on the nullclines the same way you did for 2., and end by adding the trajectory just like you did for 2., too.
par(mfrow = c(1, 3), mar = c(4.5, 4.5, 0.5, 0.5))

plot(x = NA, type = "n", xlim = c(0,100), ylim = c(0, 8), xlab = "Prey density (V)", ylab = "Predator density (P)", las = 1)
trash <- flowField(deriv = MR, xlim = c(0,100), ylim = c(0, 8), parameters = parm_vec_eq, points = 10, add = T, state.names = c("V", "P"), las = 1)
trash <- nullclines(deriv = MR, xlim = c(-1,100), ylim = c(-0.1, 8), parameters = parm_vec_eq, add.legend = F, col = c("blue", "red"), state.names = c("V", "P"), points = 200)
points(x = MRode_prey[1, 2], y = MRode_prey[1, 3], pch = 19)
lines(x = MRode_prey[, 2], y = MRode_prey[, 3])

plot(x = NA, type = "n", xlim = c(0,100), ylim = c(0, 8), xlab = "Prey density (V)", ylab = "Predator density (P)", las = 1)
trash <- flowField(deriv = MR, xlim = c(0,100), ylim = c(0, 8), parameters = parm_vec_stab, points = 10, add = T, state.names = c("V", "P"))
trash <- nullclines(deriv = MR, xlim = c(-1,100), ylim = c(-0.1, 8), parameters = parm_vec_stab, add.legend = F, col = c("blue", "red"), state.names = c("V", "P"), points = 200)
points(x = MRode_stab[1, 2], y = MRode_stab[1, 3], pch = 19)
lines(x = MRode_stab[, 2], y = MRode_stab[, 3])

plot(x = NA, type = "n", xlim = c(0,100), ylim = c(0, 8), xlab = "Prey density (V)", ylab = "Predator density (P)", las = 1)
trash <- flowField(deriv = MR, xlim = c(0,100), ylim = c(0, 8), parameters = parm_vec_unstab, points = 10, add = T, state.names = c("V", "P"))
trash <- nullclines(deriv = MR, xlim = c(-1,100), ylim = c(-0.1, 8), parameters = parm_vec_unstab, add.legend = F, col = c("blue", "red"), state.names = c("V", "P"), points = 200)
points(x = MRode_unstab[1, 2], y = MRode_unstab[1, 3], pch = 19)
lines(x = MRode_unstab[, 2], y = MRode_unstab[, 3])