This is the a working document from lab 06 as a first step of how to create a 2D bifurcation biagram. We did not come close to making the bifurcation diagram, but we took a few steps in that direction. We will finish this for the next lab. The steps we took and needed to take are: 1. Get a working model verified by a time series (done) 2. Create nullclines since we can do analysis (find nullclines) with this model (done) 3. Run the search algorithm as an alternative method to find nullclines using phaseR (done) 4. Use a loop to simulate over two parameters (e.g., \(K\) and \(m\)) 5. Convert time series data to qualitative dynamics (i.e., prey only, stable coexistence, oscillating) 6. Make bifurcation plot

Time series

RM <- function(t, y, parm) {
  with(as.list(c(y, parm)), {
    dV <- r*V*(K-V)/K - c*V*P/(g + V)
    dP <- b*V*P/(g + V) - m*P
    return(list(c(dV, dP)))
  })
}

N0 <- c(V = 0.1, P = 0.2)
time_vec <- seq(from = 0, to = 100, length.out = 10000)
parm_vec <- c(r = 1, K = 10, c = 2, g = 1, b = 1.5, m = 0.5)

RM1 <- ode(y = N0, times = time_vec, func = RM, parms = parm_vec)

plot(x = RM1[,1], y = RM1[,2], xlab = "Time", ylab = "Density", las = 1, type = "l")
  lines(x = RM1[,1], y = RM1[,3], col ="tomato")

Create a phase plane, manually

V_vec <- seq(from = 0, to = 12, length.out = 1000)
# V_NC <- r/c*(g + (1 - g/K)*V - V*V/K)
V_NC <- parm_vec["r"]/parm_vec["c"]*(parm_vec["g"] + (1 - parm_vec["g"]/parm_vec["K"])*V_vec - V_vec*V_vec/parm_vec["K"])
P_NC <- parm_vec["g"]*parm_vec["m"] / (parm_vec["b"] - parm_vec["m"])
  
plot(x = V_vec, y = V_NC, type = "l", col = "tomato", xlab = "Prey density (V)", ylab = "Predator density (P", las = 1)
segments(x0 = 0, x1 = 0, y0 = -2, y1 = 2, col = "tomato")
segments(x0 = 0, x1 = 12, y0 = 0, y1 = 0, col = "skyblue")
segments(x0 = P_NC, x1 = P_NC, y0 = -2, y1 = 2, col = "skyblue")

Create a phase plane, using phaseR

# Without assigning something to phaseR functions like `nullclines()`, then it prints out a lot of diagnostic information you don't need
trash <- nullclines(deriv = RM, xlim = c(-1, 12), ylim = c(-1, 2), parameters = parm_vec, state.names = c("V", "P"), col = c("tomato", "skyblue"), add = F, points = 200, las = 1)