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)

Nested for loop

n_lens <- 10
n_wids <- 10
len <-seq(from = 1, to = 10, length.out = n_lens) 
wid <-seq(from = 11, to = 20, length.out = n_wids)

area_mat <- matrix(data = NA, nrow = n_lens, ncol = n_wids)

for (l in 1:n_lens) {
  for (w in 1:n_wids) {
   # print(c(len[l], wid[w]))
  area <- len[l]*wid[w]
  area_mat[l, w] <- area
  }
}

My attempt to create a nested for loop with R-M

#[Add your attempt here]

Create a nested for loop with R-M

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)
n_iter <- 10000
time_vec <- seq(from = 0, to = 100, length.out = n_iter)
parm_vec <- c(r = 1, K = NA, c = 2, g = 1, b = 1.5, m = NA)

n_Ks <- 100
n_ms <- 100
K_vec <- seq(from = 2, to = 45, length.out = n_Ks)
m_vec <- seq(from = 1, to = 2, length.out = n_ms)

V_arr <- array(data = NA, dim = c(n_Ks, n_ms, n_iter))
P_arr <- array(data = NA, dim = c(n_Ks, n_ms, n_iter))

for (i in 1:n_Ks) {
  for (j in 1:n_ms) {
    parm_vec["K"] <- K_vec[i]
    parm_vec["m"] <- m_vec[j]
    
    RM1 <- ode(y = N0, times = time_vec, func = RM, parms = parm_vec)

    V_arr[i, j, ] <- RM1[,2]
    P_arr[i, j, ] <- RM1[,3]
  }
}

Create bifurcation diagram for P/A of predator

image(P_arr[,,n_iter])

P_bin <- P_arr[,,n_iter]
  P_bin[P_bin > 0.0001] <- 1
  P_bin[P_bin <= 0.0001] <- 0

image(x = K_vec, y = m_vec, z = P_bin, xlab = "K values", ylab = "m values", col = c("tomato","skyblue"), las = 1)

Create bifurcation diagram for 3 dynamical behaviors

P_bif_1 <- P_arr[,,(n_iter-8999):n_iter]
P_bin <- P_arr[,,n_iter]
  P_bif <- matrix(apply(X = P_bif_1, MARGIN = c(1,2), var), nrow = n_Ks, ncol = n_ms)
  P_bin[P_bif < 0.0001 & P_bin < 0.0001] <- 0 # is asymptotic, pred extinct
  P_bin[P_bif < 0.0001 & P_bin >= 0.0001] <- 1 # is asymptotic, stable coexistence
  P_bin[P_bif >= 0.0001] <- 2 # is unstsable

image(x = K_vec, y = m_vec, z = P_bin, xlab = "K values", ylab = "m values", col = c("gold", "tomato", "skyblue"), las = 1)