Let’s set some parameters.

r <- 0.5 # growth rate
K <- 100 # carrying capacity
N0 <- 2 # starting population size
A <- K/2 # threshold for Allee effect
E <- 1 # effort

Equilibrium population size as a function of q, the catchability coefficient

# vector of timesteps
tset <- seq(from = 0, to = 100, length.out = 20000)

# vector of values for q (catchability coefficient)
qset <- seq(from = 0, to = 1, length = 100)

# holding vectors for 4 different scenarios:
# 1. harvest with logistic growth
N.star.log.q <- NaN*qset
# 2. harvest with depensating growth with a high starting population size
N.star.dep.up.q <- NaN*qset
# 3. harvest with depensating growth with a low starting population size
N.star.dep.down.q <- NaN*qset
# 4. harvest with overcompensating growth
N.star.over.q <- NaN*qset

for(j in 1:length(qset)) {
  
  q <- qset[j]
  
  N.simu.log.q <- NaN*tset
  N.simu.log.q[1] <- N0
  
  N.simu.dep.up.q <- NaN*tset
  N.simu.dep.up.q[1] <- 100
  
  N.simu.dep.down.q <- NaN*tset
  N.simu.dep.down.q[1] <- N0
  
  N.simu.over.q <- NaN*tset
  N.simu.over.q[1] <- N0
  
  for(i in 2:length(tset)) {
    
    dt <- tset[i] - tset[i-1]
    
    # 1. harvest with logistic growth
    dNdt.log <- (r * N.simu.log.q[i-1] * (1 - N.simu.log.q[i-1]/K) - q*E*N.simu.log.q[i-1]) * dt
    N.simu.log.q[i] <- N.simu.log.q[i-1] + dNdt.log
    
    # 2. harvest with depensating growth with a high starting population size
    dNdt.dep.up <- (r * N.simu.dep.up.q[i-1] * (1 - N.simu.dep.up.q[i-1]/K) * (N.simu.dep.up.q[i-1]/A) - q*E*N.simu.dep.up.q[i-1]) * dt
    N.simu.dep.up.q[i] <- N.simu.dep.up.q[i-1] + dNdt.dep.up
    
    # 3. harvest with depensating growth with a low starting population size
    dNdt.dep.down <- (r * N.simu.dep.down.q[i-1] * (1-N.simu.dep.down.q[i-1]/K) * (N.simu.dep.down.q[i-1]/A) - q*E*N.simu.dep.down.q[i-1]) * dt
    N.simu.dep.down.q[i] <- N.simu.dep.down.q[i-1] + dNdt.dep.down   
    
    # 4. harvest with overcompensating growth
    dNdt.over <- (r*N.simu.over.q[i-1] * (1 - N.simu.over.q[i-1]/K) * (2 - N.simu.over.q[i-1]/K) - q*E*N.simu.over.q[i-1]) * dt
    N.simu.over.q[i] <- N.simu.over.q[i-1] + dNdt.over
    
  }
  
  N.star.log.q[j] <- N.simu.log.q[length(N.simu.log.q)]
  N.star.dep.up.q[j] <- N.simu.dep.up.q[length(N.simu.dep.up.q)]
  N.star.dep.down.q[j] <- N.simu.dep.down.q[length(N.simu.dep.down.q)]
  N.star.over.q[j] <- N.simu.over.q[length(N.simu.over.q)]
  
}
plot(x = qset, y = N.star.log.q,
     type = 'l', 
     lwd = 2, 
     xlab = 'Catchability coefficient', ylab = 'Equilibrium population size, N*',
     col = "black")
lines(x = qset, y = N.star.dep.down.q,
      lwd = 2,
      col = "lightblue")
lines(x = qset, y = N.star.dep.up.q,
      lwd = 2,
      col = "coral")
lines(x = qset, y = N.star.over.q,
      lwd = 2,
      col = "tan")
legend(x = 0.55, y = 99,
       legend = c("logistic", "depensating (N0 = 100)", 
                  "depensating (N0 = 2)", "overcompensating"),
       col = c("black", "lightblue", "coral", "tan"),
       lwd = 2)

Equilibrium population size as a function of E, harvesting effort

I’ll use \(q = 0.2\).

q <- 0.2
# vector of values for E (effort)
Eset <- seq(from = 0, to = 6, length.out = 200)

# holding vectors for 4 different scenarios:
# 1. harvest with logistic growth
N.star.log.E <- NaN*Eset
# 2. harvest with depensating growth with a high starting population size
N.star.dep.up.E <- NaN*Eset
# 3. harvest with depensating growth with a low starting population size
N.star.dep.down.E <- NaN*Eset
# 4. harvest with overcompensating growth
N.star.over.E <- NaN*Eset

for(j in 1:length(Eset)) {
  
  E <- Eset[j]
  
  N.simu.log.E <- NaN*tset
  N.simu.log.E[1] <- N0
  
  N.simu.dep.up.E <- NaN*tset
  N.simu.dep.up.E[1] <- 100
  
  N.simu.dep.down.E <- NaN*tset
  N.simu.dep.down.E[1] <- N0
  
  N.simu.over.E <- NaN*tset
  N.simu.over.E[1] <- N0
  
  for(i in 2:length(tset)) {
    
    dt <- tset[i] - tset[i-1]
    
    # 1. harvest with logistic growth
    dNdt.log <- (r * N.simu.log.E[i-1] * (1 - N.simu.log.E[i-1]/K) - q*E*N.simu.log.E[i-1]) * dt
    N.simu.log.E[i] <- N.simu.log.E[i-1] + dNdt.log
    
    # 2. harvest with depensating growth with a high starting population size
    dNdt.dep.up <- (r * N.simu.dep.up.E[i-1] * (1 - N.simu.dep.up.E[i-1]/K) * (N.simu.dep.up.E[i-1]/A) - q*E*N.simu.dep.up.E[i-1]) * dt
    N.simu.dep.up.E[i] <- N.simu.dep.up.E[i-1] + dNdt.dep.up
    
    # 3. harvest with depensating growth with a low starting population size
    dNdt.dep.down <- (r * N.simu.dep.down.E[i-1] * (1-N.simu.dep.down.E[i-1]/K) * (N.simu.dep.down.E[i-1]/A) - q*E*N.simu.dep.down.E[i-1]) * dt
    N.simu.dep.down.E[i] <- N.simu.dep.down.E[i-1] + dNdt.dep.down   
    
    # 4. harvest with overcompensating growth
    dNdt.over <- (r*N.simu.over.E[i-1] * (1 - N.simu.over.E[i-1]/K) * (2 - N.simu.over.E[i-1]/K) - q*E*N.simu.over.E[i-1]) * dt
    N.simu.over.E[i] <- N.simu.over.E[i-1] + dNdt.over
    
  }
  
  N.star.log.E[j] <- N.simu.log.E[length(N.simu.log.E)]
  N.star.dep.up.E[j] <- N.simu.dep.up.E[length(N.simu.dep.up.E)]
  N.star.dep.down.E[j] <- N.simu.dep.down.E[length(N.simu.dep.down.E)]
  N.star.over.E[j] <- N.simu.over.E[length(N.simu.over.E)]
  
}
plot(x = Eset, y = N.star.log.E,
     type = 'l', 
     lwd = 2, 
     xlab = 'Harvesting effort', ylab = 'Equilibrium population size, N*',
     col = "black")
lines(x = Eset, y = N.star.dep.down.E,
      lwd = 2,
      col = "lightblue")
lines(x = Eset, y = N.star.dep.up.E,
      lwd = 2,
      col = "coral")
lines(x = Eset, y = N.star.over.E,
      lwd = 2,
      col = "tan")
legend(x = 3.25, y = 99,
       legend = c("logistic", "depensating (N0 = 100)", 
                  "depensating (N0 = 2)", "overcompensating"),
       col = c("black", "lightblue", "coral", "tan"),
       lwd = 2)

Is the similarity of these two plots surprising to you? Why, or why not?