Question 1:
Create an R script that contains the function LogisticModel() from your text. Simulate the dynamics of a population growing logistically in a fluctuating environment for 1000 years, in the absence of harvesting, with R0=0.2, K=100, N0=100, and v=0.15.

  1. Using R, plot a graph of population size over time, and a frequency distribution of population sizes. Consider a population of fish described by the Logistic equation with R0 = 0.1, K=100.
# Define parameters
R0 <- 0.2
K <- 100
N0 <- 100
v <- 0.15
t <- 1000
time <- seq(0, t, by=1)

# Define logistic growth function
Logisticmodel <- function(N, R, K, v) {
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N
  return(dN)
}

# Calculate population sizes over time
population <- numeric(length(time))
population[1] <- N0

for (i in 2:length(time)) {
  population[i] <- population[i - 1] + Logisticmodel(population[i - 1], R0, K, v)
}

# Plot logistic growth vs time
plot(time, population, type = "l", col = "#2970a4",
     xlab = "Time (t)", ylab = expression("Population Size ( " * N[t] * "):"),
     main = "Population Size Over Time")


# Create a frequency distribution histogram
hist(population, breaks=20, col = "#2970a4", xlab = expression("Population Size (" * N[t] * "):"), ylab = "Frequency",
     main = "Frequency Distribution of Population Sizes")


  1. Interpret the shape of the frequency distribution of population sizes (N), and compare it with the histogram of R deviations from your text. How is the distribution of N values different? Why do you think these differences might be?

The histogram in the textbook has roughly the same frequencies for every bar, as expected, due to it being uniform. However, the frequency of population sizes somewhat resembles a positively skewed normal distribution. This is likely because of the change in population growth around K, i.e. the density dependent property. If the population is below K, it is positive and accelerates, but will slow down if it reaches or goes above it, leading to fluctuations around K. This is the self-correcting aspect of the model, causing the population to most frequently stay around K, as seen in the histogram.

Question 2:

  1. Use calculus to show that the function G(N) is maximized when N=K/2.

Calculating the derivative of G(N): \[ \begin{align*} G(N) &= \left(R_0 - \frac{R_0}{K}N\right)N \\ \frac{dG}{dN} &= \left(\frac{-R_0}{K}\right)N + \left(R_0 - \frac{R_0}{K}N\right) \\ &= R_0\left(\frac{-2}{K}N + 1\right) \end{align*} \]

Proving \(\frac{K}{2}\) is an extreme point: \[ \begin{align*} \frac{dG}{dN} &= 0 \\ R_0\left(\frac{-2}{K}N + 1\right) &= 0 \\ \frac{-2}{K}N + 1 &= 0 \\ \frac{2}{K}N &= 1 \\ N &= \frac{K}{2} \end{align*} \] Proving \(\frac{K}{2}\) is a maximum:

\[ \begin{align*} \frac{d^2G}{dN^2} &= \frac{-2R_0}{K} \\ \frac{d^2G}{dN^2} &< 0 \end{align*} \] Hence, as \(\frac{d^2G}{dN^2} < 0\), \(N = \frac{K}{2}\) is a maximum.

  1. Find maximum sustainable yield (MSY; the value of G(N) at its maximum) as a function of R0 and K. Show your workings below.

Plugging \(\frac{K}{2}\) into G(N): \[ \begin{align*} MSY &= G(\frac{K}{2}) \\ &= \left(R_0 - \frac{R_0}{K}*\frac{K}{2}\right)\frac{K}{2} \\ &= \left(\frac{R_0}{2}\right)\frac{K}{2} \\ &= \frac{R_0K}{4} \end{align*} \] Question 3:
Modify your simulation from question 1 to incorporate harvesting at MSY, set the initial population size (N0) to be equal to the population size that produces MSY, and reduce the length of the simulation to 100 years. Leave the other variables unchanged. Run this simulation 4 times, and create a graph of population size over time for each simulation. How are your graphs similar? How are they different?

# Define parameters
R0 <- 0.2
K <- 100
N0 <- K/2
v <- 0.15
t <- 100
time <- seq(0, t, by=1)

# Define logistic growth function with harvest
Logisticmodel_h <- function(N, R, K, v){
  x <- runif(1,-v,v)
  MSY <- R0 * K / 4
  h <- MSY
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N - h
  return(dN)
}

# Calculate population sizes over time
population <- numeric(length(time))
population[1] <- N0

for (i in 2:length(time)) {
  population[i] <- population[i - 1] + Logisticmodel_h(population[i - 1], R0, K, v)
}

# Create a 2x2 grid of plots with adjusted margins
par(mfrow = c(2, 2))

for (i in 1:4) {
  # Calculate population sizes over time
  population <- numeric(length(time))
  population[1] <- N0
  
  for (j in 2:length(time)) {
    new_population <- population[j - 1] + Logisticmodel_h(population[j - 1], R0, K, v)
    
    # Ensure population doesn't go negative and set to 0 if negative
    population[j] <- max(new_population, 0)
  }
  
  # Adjust margins for the current plot
  par(mar = c(4, 4, 1, 2))  # Adjust top, right, bottom, left margins
  
  # Create and display plot for the current simulation
  plot(time, population, type = "l", col = "#2970a4",
       xlab = "Time (t)", ylab = expression("Population Size: " * N[t]),
       ylim = c(0, 100))  # Set y-axis limits
  
  # Adjust x-axis ticks
  x_ticks_small <- seq(0, max(time), by = 10)
  axis(1, at = x_ticks_small, labels = FALSE, tcl = -0.2)
  
  # Adjust y-axis ticks
  y_ticks_big <- seq(0, 100, by = 20)
  y_ticks_small <- seq(0, 100, by = 10)
  axis(2, at = y_ticks_big, labels = y_ticks_big)
  axis(2, at = y_ticks_small, labels = FALSE, tcl = -0.2)
}
Four Simulations of Logistic Growth with R0=0.2 ±0.15, K=100, N0=50 and h=5.


Two simulations show the population going extinct, while two it survives. One goes extinct pretty much immediately, while another holds on for slightly longer. One of the surviving populations (top right), despite not going extinct, has declined quite a bit, and there is a good likelihood it goes extinct if the simulation ran for a few more decades.

Question 4
(a) Run the simulation you designed in (3) for a total of 20 times, using the same parameter values and starting conditions as you did in question 3. Plot all 20 trajectories on a graph, and determine how many of these 20 populations went extinct in the 100 year period that you simulated. Show the plot and comment on the output.

# Define parameters
R0 <- 0.2
K <- 100
MSY <- R0 * K / 4
N0 <- K/2
v <- 0.15
t <- 100
time <- seq(0, t, by=1)

# Define logistic growth function with harvest
Logisticmodel_h <- function(N, R, K, v){
  x <- runif(1,-v,v)
  MSY <- R0 * K / 4
  h <- MSY
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N - h
  return(dN)
}

# Initialize plot
par(mfrow = c(1, 1))
plot(time, numeric(length(time)), type = "l", col = "white",
     xlab = "Time (t)", ylab = expression("Population Size ( " * N[t] * "):"),
     main = "Logistic Growth Simulated 20 Times",
     ylim = c(0,100))

# Simulate and plot 20 times
zero_counts <- numeric(20)  # To store zero counts
for (i in 1:20) {
  population <- numeric(length(time))
  population[1] <- N0
  hit_zero <- FALSE  # Flag to track if population hits zero, 
  # ensures if function only works the first time pop hits zero
  
  for (j in 2:length(time)) {
    new_population <- population[j - 1] + Logisticmodel_h(population[j - 1], R0, K, v)
    population[j] <- max(new_population, 0) # Replaces negative population values with zero, 
    # so the line trails along zero when it goes extinct
    
        if (!hit_zero && population[j] == 0) {
      zero_counts[i] <- zero_counts[i] + 1
      hit_zero <- TRUE
        }
  }
  
  lines(time, population, col = "#2970a4")
}
# Output zero counts below the plot
extinction_count <- sum(zero_counts)
paste("The population went extinct", extinction_count, "times.")


## [1] "The population went extinct 17 times."
  1. Based on the proportion of these simulations that went extinct, what would be your estimate of the probability that a population with these parameter values, harvested at MSY, would be extinct within 100 years?

\[ \begin{align*} P(population\ going\ extinct\ after\ 100\ years) &= \frac{no.\ of\ simulations\ where\ it\ goes\ extinct}{total\ simulations} \\ &= \frac{17}{20} \\ &= 0.85 \\ \end{align*} \]

Next, let’s explore what happens when we try to reduce the harvest, and manage the population at the stable and unstable equilibria, respectively.

Question 5:
Let’s try managing the population to the right of the peak in the production function, at N=3K/4. Set this to be the approximate population size, N0. Then set the harvest, h, at what it would need to be, in order to keep this population at equilibrium, if there were no randomness in the population dynamics. Leave R0=0.2, K=100, and v=0.15. Run another 20 simulations, plot them on a graph, and estimate the probability of extinction for this hypothetical population managed with this harvesting strategy.

The harvest was determined by evaluating \(G(\frac{3K}{4})\), giving a value of \(\frac{3R_0K}{16}\).

# Define parameters
R0 <- 0.2
K <- 100
N0 <- 3*K/4
v <- 0.15
t <- 100
h <- 3*R0*K/16
time <- seq(0, t, by=1)

# Define logistic growth function with harvest
Logisticmodel_h <- function(N, R, K, v){
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N - h
  return(dN)
}

# Initialize plot
par(mfrow = c(1, 1))
plot(time, numeric(length(time)), type = "l", col = "white",
     xlab = "Time (t)", ylab = expression("Population Size ( " * N[t] * "):"),
     main = "Logistic Growth with N=3K/4, h=3.75, Simulated 20 Times",
     ylim = c(0,120))

# Simulate and plot 20 times
zero_counts <- numeric(20)  # To store zero counts
for (i in 1:20) {
  population <- numeric(length(time))
  population[1] <- N0
  hit_zero <- FALSE  # Flag to track if population hits zero, 
  # ensures if function only works the first time pop hits zero
  
  for (j in 2:length(time)) {
    new_population <- population[j - 1] + Logisticmodel_h(population[j - 1], R0, K, v)
    population[j] <- max(new_population, 0) # Replaces negative population values with zero, 
    # so the line trails along zero when it goes extinct
    
        if (!hit_zero && population[j] == 0) {
      zero_counts[i] <- zero_counts[i] + 1
      hit_zero <- TRUE
        }
  }
  
  lines(time, population, col = "#2970a4")
}
# Output zero counts below the plot
extinction_count <- sum(zero_counts)
prob_extinct <- sum(zero_counts)/20
paste("The population went extinct", extinction_count, "times.", "The probability of extinction when N=3K/4 is", prob_extinct)


## [1] "The population went extinct 1 times. The probability of extinction when N=3K/4 is 0.05"

Question 6:
See what happens when we instead try to manage the population to the left of the production function, at N=1/4K. Set this as the initial population size, N0. Then set the harvest, h, at what it would need to be, in order to keep this population at equilibrium, if there were no randomness in the population dynamics. How does it compare with the harvest level in question 5? Again, leave R0=0.2, K=100, and v=0.15, run 20 simulations, plot them on a graph, and estimate the probability of extinction for this hypothetical population managed with this harvesting strategy.

The harvest would be the same, as the production function is symmetrical about K/2, and as 3K/4 and K/4 are the same distance from it, they have the same intersection required harvest for equilibrium. However, now we are harvesting the same amount from a smaller population, so we should expect a higher extinction rate.

# Define parameters
R0 <- 0.2
K <- 100
N0 <- K/4
v <- 0.15
t <- 100
h <- 3*R0*K/16
time <- seq(0, t, by=1)

# Define logistic growth function with harvest
Logisticmodel_h <- function(N, R, K, v){
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N - h
  return(dN)
}

# Initialize plot
par(mfrow = c(1, 1))
plot(time, numeric(length(time)), type = "l", col = "white",
     xlab = "Time (t)", ylab = expression("Population Size ( " * N[t] * "):"),
     main = "Logistic Growth with N=K/4, h=3.75, Simulated 20 Times",
     ylim = c(0,120))

# Simulate and plot 20 times
zero_counts <- numeric(20)  # To store zero counts
for (i in 1:20) {
  population <- numeric(length(time))
  population[1] <- N0
  hit_zero <- FALSE  # Flag to track if population hits zero, 
  # ensures if function only works the first time pop hits zero
  
  for (j in 2:length(time)) {
    new_population <- population[j - 1] + Logisticmodel_h(population[j - 1], R0, K, v)
    population[j] <- max(new_population, 0) # Replaces negative population values with zero, 
    # so the line trails along zero when it goes extinct
    
        if (!hit_zero && population[j] == 0) {
      zero_counts[i] <- zero_counts[i] + 1
      hit_zero <- TRUE
        }
  }
  
  lines(time, population, col = "#2970a4")
}
# Output zero counts below the plot
extinction_count <- sum(zero_counts)
prob_extinct <- sum(zero_counts)/20
paste("The population went extinct", extinction_count, "times.", "The probability of extinction when N=K/4 is", prob_extinct)


## [1] "The population went extinct 12 times. The probability of extinction when N=K/4 is 0.6"

As expected, the extinction rate went up, from 0.05 to 0.6.

Question 7:
Increase the magnitude of random variation in R to 0.25 (i.e., v=0.25), and repeat questions 5 and 6.

# Define parameters
R0 <- 0.2
K <- 100
N0 <- K/4
v <- 0.25
t <- 100
h <- 3*R0*K/16
time <- seq(0, t, by=1)

# Define logistic growth function with harvest
Logisticmodel_h <- function(N, R, K, v){
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N - h
  return(dN)
}

# Initialize plot
par(mfrow = c(1, 1))
plot(time, numeric(length(time)), type = "l", col = "white",
     xlab = "Time (t)", ylab = expression("Population Size ( " * N[t] * "):"),
     main = "Logistic Growth with N=3K/4, h=3.75, v=0.25, Simulated 20 Times",
     ylim = c(0,160))

# Simulate and plot 20 times
zero_counts <- numeric(20)  # To store zero counts
for (i in 1:20) {
  population <- numeric(length(time))
  population[1] <- N0
  hit_zero <- FALSE  # Flag to track if population hits zero, 
  # ensures if function only works the first time pop hits zero
  
  for (j in 2:length(time)) {
    new_population <- population[j - 1] + Logisticmodel_h(population[j - 1], R0, K, v)
    population[j] <- max(new_population, 0) # Replaces negative population values with zero, 
    # so the line trails along zero when it goes extinct
    
        if (!hit_zero && population[j] == 0) {
      zero_counts[i] <- zero_counts[i] + 1
      hit_zero <- TRUE
        }
  }
  
  lines(time, population, col = "#2970a4")
}
# Output zero counts below the plot
extinction_count <- sum(zero_counts)
prob_extinct <- sum(zero_counts)/20
paste("The population went extinct", extinction_count, "times.", "The probability of extinction when N=3K/4 is", prob_extinct)


## [1] "The population went extinct 8 times. The probability of extinction when N=3K/4 is 0.4"
# Define parameters
R0 <- 0.2
K <- 100
MSY <- R0 * K / 4
N0 <- K/4
v <- 0.25
t <- 100
h <- 3*R0*K/16
time <- seq(0, t, by=1)

# Define logistic growth function with harvest
Logisticmodel_h <- function(N, R, K, v){
  x <- runif(1,-v,v)
  dN <- (R * (1 - N / K) + x) * N - h
  return(dN)
}

# Initialize plot
par(mfrow = c(1, 1))
plot(time, numeric(length(time)), type = "l", col = "white",
     xlab = "Time (t)", ylab = expression("Population Size ( " * N[t] * "):"),
     main = "Logistic Growth with N=K/4, h=3.75, v=0.25, Simulated 20 Times",
     ylim = c(0,160))

# Simulate and plot 20 times
zero_counts <- numeric(20)  # To store zero counts
for (i in 1:20) {
  population <- numeric(length(time))
  population[1] <- N0
  hit_zero <- FALSE  # Flag to track if population hits zero, 
  # ensures if function only works the first time pop hits zero 
  
  for (j in 2:length(time)) {
    new_population <- population[j - 1] + Logisticmodel_h(population[j - 1], R0, K, v)
    population[j] <- max(new_population, 0) 
    
        if (!hit_zero && population[j] == 0) {
      zero_counts[i] <- zero_counts[i] + 1
      hit_zero <- TRUE
        }
  }
  
  lines(time, population, col = "#2970a4")
}
# Output zero counts below the plot
extinction_count <- sum(zero_counts)
prob_extinct <- sum(zero_counts)/20
paste("The population went extinct", extinction_count, "times.", "The probability of extinction when N=K/4 is", prob_extinct)


## [1] "The population went extinct 16 times. The probability of extinction when N=K/4 is 0.8"

The extinction rates for N=3K/4 and N=K/4 went up from 0.05 and 0.6 to 0.4 and 0.8 respectively with v set to 0.25. There were also high maximum population sizes. Both higher extinction rates and maximum population sizes would be expected to increase with higher variation.

Question 8:
Figure out what your target stock size would need to be (to within ±5), in order to reduce the probability of extinction for a population with this level of variation to <15%. When you test this out, set your initial population size, N0, to be equal to the target population size that you’re trying. Try increasing your number of simulations to 100, to get a better estimate. Explain how you made your determination of the necessary target stock size. What is the harvest associated with this target stock size?

Using the same parameters as previously, I incrementally increased the value of N0 by 1 with an initial value of \(\frac{K}{2}\) or 50. This continually decreased the probability of extinction until it was less than 0.15 (15%), whereby I then rounded the harvest to the nearest lowest integer. I repeated this several times, and the target stock size tended to fall in the range of 77-80. This range of values would be normal distributed so we can expect the mean to be within this range, and values above or below it would become increasingly more uncommon. In a real case, being conservative with this would be better, which is partially why harvest is always rounded down, the other reason being you cannot harvest part of a fish.

# Define parameters
R0 <- 0.2
K <- 100
N0 <- K / 2
v <- 0.25
t <- 100

time <- seq(0, t, by=1)

# Initialise variables
target_prob <- 0.15
results <- data.frame(N0 = numeric(0), 
                      ProbExtinct = numeric(0), 
                      Harvest = numeric(0))  # Empty data frame to store results

prob_extinct <- 1

while(N0 <= 100){
 
  h <- N0 * R0 * (1 - N0 / K)  # Calculate h based on N0
  
  # Define logistic growth function with harvest
  Logisticmodel_h <- function(N, R, K, v){
    x <- runif(1,-v,v)
    dN <- (R * (1 - N / K) + x) * N - h
    return(dN)
  }
  
  zero_counts <- numeric(100)  # To store zero counts
  
    for (i in 1:100) {
    population <- numeric(length(time))
    population[1] <- N0
    hit_zero <- FALSE  # Flag to track if population hits zero, 
    # ensures if function only works the first time pop hits zero
  
      for (j in 2:length(time)) {
        new_population <- population[j - 1] + Logisticmodel_h(population[j - 1],
                                                              R0, K, v)
        population[j] <- max(new_population, 0)
      
        if (!hit_zero && population[j] == 0) {
          zero_counts[i] <- zero_counts[i] + 1
          hit_zero <- TRUE
        }
      }
    }
  
  
  # Calculate extinction probability
  prob_extinct <- sum(zero_counts) / 100
  
  # Append results to the data frame 
  results <- rbind(results, data.frame(N0 = N0, 
                                       ProbExtinct = prob_extinct, 
                                       Harvest = h))
  
  if (prob_extinct < target_prob) {
      break  # Exit loop when probability is below target
  }
  
  N0 <- N0 + 1  # Increase N0 by 1 for the next iteration
}

harvest <- results$Harvest[results$N0 == N0]
rounded_harvest <- floor(harvest)
print(results[c("N0", "ProbExtinct")])
paste("The target stock size is", N0, "with an associated harvest of", 
      harvest, "or roughly", rounded_harvest)
##    N0 ProbExtinct
## 1  50        0.94
## 2  51        0.90
## 3  52        0.92
## 4  53        0.94
## 5  54        0.91
## 6  55        0.88
## 7  56        0.82
## 8  57        0.88
## 9  58        0.85
## 10 59        0.85
## 11 60        0.80
## 12 61        0.87
## 13 62        0.78
## 14 63        0.77
## 15 64        0.81
## 16 65        0.80
## 17 66        0.82
## 18 67        0.68
## 19 68        0.66
## 20 69        0.65
## 21 70        0.53
## 22 71        0.38
## 23 72        0.51
## 24 73        0.47
## 25 74        0.38
## 26 75        0.35
## 27 76        0.26
## 28 77        0.19
## 29 78        0.16
## 30 79        0.12
## [1] "The target stock size is 79 with an associated harvest of 3.318 or roughly 3."

Question 9:
(Short essay). In 200-350 words, assess the robustness of our “rules of thumb” for sustainable harvesting from lecture. Do your results in questions 1 to 8 suggest that this needs to be modified? If so, how?

To assess the robustness of our model, we want to see the effects of variation on our set parameters, and if it significantly influences the outcome. In this case, we introduced variation to R, the net per capita rate of increase, and comparing the v=0.15 and v=0.25 cases, increasing variation increases the chances of extinction. This is somewhat expected, but is also problematic for applying the model to reality, as variation is likely much higher in reality, due to numerous unaccounted for factors. This suggests the model is not that robust, and likely needs modifications.

One flaw in the model is the isolation of years. It’s assumed population growth each year only depends on the previous year’s population, however, we’d expect various lag effects if population growth is lower one year. This is due to the age structure of populations, meaning only a proportion of the population are reproducing individuals. For example, overharvesting could reduce population size one year, but would not cause an immediate decline in births and instead be felt by the population several years later, potentially causing a delayed population collapse. This requires a more precise model representing the proper population dynamics of the Atlantic cod.

Another flaw is treating the h as constant. In reality, we would also expect it to vary, influencing the probability of extinction. In tandem with other unaccounted factors, like predator-prey dynamics, there could be serious consequences on the population. Introducing variation to harvest, as well as R, would allow for more realistic outputs.

In question 8, where a target probability of extinction determines the harvest is effective, but must account for two things. First, variation, and hence, the probability of extinction is likely higher than expected. Second, what should the target probability be? At 0.15, there’s roughly a 1 in 6 chance of extinction, probably even higher, so it should be lower, but how much lower? We also don’t want to underutilise the fishery. Overall, the model does not accurately represent reality, by definition. Still, accounting for higher variation and more conservative harvesting should be implemented if sustainable fishing is desirable.