Solutions_Computer_Lab_ECFM_Harvesting.Rmd

---
title: "Solutions for Computer Lab ECFM: Population dynamics and harvesting"
author: "Max Lindmark & Anna Gårdmark"
date: "24 augusti 2018"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Overview
In this practical you will use The Schaefer Model [1] to explore how harvesting affects population numbers and dynamics. You will be using R to explore the dynamics of the model. R is a statistical software that is very popular in especially biological sciences and can be downloaded for free. If you have not installed R and R-studio already, please see the document "First time R user". R-studio is just a interface for R with many helpful utilities. It is an advantage if you know some Basic R, but you should be able to follow along this lab without any prior knowledge if you consult the "First time R user". In this document R-code is text with grey background. You can copy these chunks of code to a new R-script in R-studio (instead of pasting it directly into the console because it quickly becomes hard to overview your work!). For some it might also be a repetition of calculus and algebra, and you will also practice some basic population dynamics and equilibrium theory.

There are hints for most questions **at the end of this document** if you need, and you can also look at the solutions or ask an instructor if you get stuck!
 
# The Schaefer Model
The Schaefer model is a harvest model based on the logistic growth equation (the Verhulst model). This is a continuous-time model and the following differential equation describes the rate of change in number of individuals, $N$: 
$$dN/dt = rN(1-N/K) - qEN$$ 
where $r$ is the intrinsic growth rate, $K$ is the carrying capacity, $q$ is catchability and $E$ is harvesting effort. In this lab, $N$ is referred to as "density" Below is some R-code to explore the population density over time in the absence of harvesting to recap the logistic growth model.

## Part 1: Recap the logistic growth model without harvesting
This part is mainly to become familiar with the growth model and read in parameters that you will use later. It is recommended you do not spend too much time on this part, so essentially just copy paste the chunks of code below until you reach Question #1!

```{r}
rm(list = ls()) # It is often a good practice to start your scrpits with this command, 
# which clears the workspace (console, "calculator") from all objects. 
# *Note though that if you want to start a completely clean R session you need to 
# quit the program the usual way.

# Two parameters govern the dynamics of the population: the growth rate (r) and  
# the carrying capacity (K).
# We can give them arbitrary values, it does not really matter here.
r <- 0.1
K <- 40

# We also need to set the initial population density and the time range 
N_ini <- 2
t <- 100

# Now we can calculate the density at the next time step!
N2 <- r * N_ini * (1-(N_ini/K)) + N_ini
N3 <- r * N2 * (1-(N2/K)) + N2
#... etc etc
```
We could repeat this iteratively, but it would be extremely tedious. This is where the loop comes in handy, which is used in programming to automate tasks. Basically, if you find yourself copy-pasting code like this, you can automate it. Here we will use a for loop. But do not worry, this lab is not about writing loops! They are just introduced so that you can get a feel of (1) how it is done (2) to explore the population dynamics.

This is essentially how you write a for loop:
```{r}
#for (each of these values) {
#    do the commands between the curly brackets
#}
```
Together with a feeling for indexing (basically that you can extract the "ith"-value of a vector using hard brackets): 
```{r}
vec <- seq(1, 20, 1)
vec
vec[3]
```
... we can build a for loop to simulate the time-dynamics in density:
```{r}
# Create an empty vector to store the values in
N <- rep(NA, t)
N[1] <- N_ini

for (i in 2:t) {
  N[i] <- r * N[i-1] * (1-(N[i - 1]/K)) + N[i-1]
}
```

```{r}
# Plot how population density change over time
plot(y = N, x = 1:t, xlab = "Time", ylab = "N", type = "l")
```

It is a bit tricky to understand exactly how density dependence works in this model. To better understand that, we can plot the population- and per capita growth rate as a function of $N$. 
```{r}
## We can start by creating a vector of N-values (ranging from 0 to carrying capacity):
N_vec <- seq(from = 0, to = K, by = 1)
N_vec

## Population growth rate
growth_pop <- r * N_vec * (1-(N_vec/K))

# Plot how population growth rate varies with population density
plot(growth_pop ~ N_vec, xlab = "N", ylab = "dN/dt", type = "l")

## Per capita growth rate
growth_cap <- growth_pop / N_vec

# Plot how per capita growth rate varies with population density
plot(growth_cap ~ N_vec, xlab = "N", ylab = "dN/dt * 1/N", type = "l")
```

* **Question 1:** a) How does density dependence operate in this model? b) When is the growth rate of the population highest? c) How could such density dependence occur in nature?

a) As the population increases in numbers, the contribution of each individual to the population growth rate declines. The decline is linear with intercept $r$ and slope $r/K$. The expression $(1-N/K)$ reaches $0$ when $N$ approaches $K$, and at really low $N$, the population grows exponentially (proportional to $N$). b) K/2. You can get this value by inspecting the graph, or solve it analytically. You can do that by taking the derivative with respect to N of the logistic growth equation, and find the value of N that leads to that expression being equal to 0 (this method finds the maxima of the curve). $N$: 
$$dN/dt = rN(1-N/K)$$ 
$$dN/dt = rN-rN^2/K$$
$$(dN/dt)' = r-2rN/K$$
$$r-2rN/K = 0$$
$$2rN/K = r$$
$$2N/K = 1$$
$$N = K/2$$

c) This could arise as a result of increased competition for a limiting resource as density increases, such as food

## Part 2: Equilibrium, harvesting and bifurcations
Now let us explore the equilibrium dynamics and dynamical stability of the Schaefer model analytically and graphically. Here asymptotic stability is an important concept. An equilibrium is asymptotically stable if the system returns to it after small deviations. The density at an equilibrium is denoted $N$*.

* **Question 2:** Which are the two equilibria of the Schaefer model?
Rearrange the Schaefer model so that growth equals mortality, and then find values of N that satisfy the equation.

The trivial equilibrium is $N$*$=0$.

The non-trivial equilibrium is given by setting the growth rate to be equal fishing mortality:
$$dN/dt = 0$$
$$rN(1-N/K)-qEN = 0$$
$$rN(1-N/K)=qEN$$
$N=0$ is the solution for the trivial equilibrium. Continue to get the second solution:
$$r(1-N/K)=qE$$
$$1-(N/K)=qE/r$$
$$K-N=KqE/r$$
$$-N=(KqE/r)-K$$
$$N=(-KqE/r)+K$$
$$N=(1-qE/r)K$$

* **Question 3:** a) Can you show these equilibria graphically? Assume $q=0.2$ (you can do this for the rest of this lab!) and $E=0.1$. b) With this harvesting regime, what is the population density at the positive equilibrium?

```{r}
plot(growth_pop ~ N_vec, xlab = "N", ylab = "dN/dt", type = "l")
q <- 0.2
E <- 0.1
lines(q*E*N ~ N, col="red")
```
N is roughly ~32 at the non-zero equilibrium.

* **Question 4:** What can we say about the dynamical stability of this equilibrium?

Small perturbations to the right along the $N$-axis to population densities above the equilibrium leads to fishing mortality (red line) being higher than the growth rate, and $N$ declines back to the equilibrium. When $N$ decreases, the growth rate is higher than the fishing mortality, and $N$ increases.

* **Question 5**: What happens to the equilibria and their stability when we *gradually* increase the fishing effort?

```{r}
plot(growth_pop ~ N_vec, xlab = "N", ylab = "dN/dt", type = "l")
q <- 0.2
E <- 0.35
lines(q*E*N ~ N, col="blue") # Stability of the equilibrium is the same!

## Try even higher Effort:
E <- 0.6
lines(q*E*N ~ N, col="green")
```
With sufficiently high mortality (green line), it is always higher than the growth rate, in which case the trivial equilibrium ($N$*$=0$) is the only equilibrium and is then stable. That is, the population goes extinct when the fishing pressure is too high.

* **Question 6**: Plot $N$* as a function of fishing effort, $E$.

Refer to the expression in Question #2. Create a vector that holds the values of effort, $E$, that you want to explore. Plot $N$* as a function of $E$:
    
```{r}
E <- seq(from = 0, to = 0.6, by = 0.05)
N_eq <- (1 - q*E/r)*K

plot(N_eq ~ E, type = "l", xlab = "E", ylab = "N*")
abline(a = 0, b = 0, lty = 2)
```

* **Question 7**: a) Can you plot yield, $Y=qEN$*, as a function of fishing effort, $E$? b) At what effort and density is long term yield maximized? (You can give an approximate value based on the last two figures). c) How do the answers in a) and b) depend on $r$ and $K$? You can look at the analytical solutions and/or redo the last two plots.  

Replace $N$* in the equation above with the expression for $N$* found in Question #2:
    
```{r}
Y <- q*E*(K*(1 - q*E/r))
plot(Y ~ E, type = "l", xlab = "E", ylab = "Yield", ylim = c(0, 1))
```
a) see the figure b) $N=20$ and $E=25$
c)

```{r}    
## Testing different r and K values for yield~effort plot
## Basically, below we recreate the last two figures while varying r and K by scaling 
## them with 0.5 and 2! ** Note, it is generally not very good practice to overwrite 
## objects like below, but at the same time creating new objects with slightly different 
## names is not perfect either. Both methods are quite error prone for different reasons, 
## and we can discuss other more safe options!
q <- 0.2
E <- seq(from = 0, to = 0.6, by = 0.05)
K <- 40
r <- 0.1

Y <- q*E*(K*(1 - q*E/r))
plot(Y ~ E, type = "l", xlab = "E", ylab = "Yield", ylim = c(0, 2.5))

Y <- q*E*(K*(1 - q*E/(r*0.5)))
lines(Y ~ E, type = "l", xlab = "E", ylab = "Yield", ylim = c(0, 1), col = "blue")

Y <- q*E*(K*(1 - q*E/(r*2)))
lines(Y ~ E, type = "l", xlab = "E", ylab = "Yield", ylim = c(0, 1), col = "blue", lty = 2)

Y <- q*E*((K*0.5)*(1 - q*E/r))
lines(Y ~ E, type = "l", xlab = "E", ylab = "Yield", ylim = c(0, 1), col = "red")

Y <- q*E*((K*2)*(1 - q*E/r))
lines(Y ~ E, type = "l", xlab = "E", ylab = "Yield", ylim = c(0, 1), col = "red", lty = 2)

legend("topleft", 
       c("default: r=0.1, K=40", 
         "r=0.05, K=40", 
         "r=0.2, K=40",
         "r=0.1, K=20",
         "r=0.1, K=80"),
       lty = c(1,1,2,1,2),
       col = c("black", "blue", "blue", "red", "red"),
       bty = "n")

## Testing different r and K values for N*~effort plot
E <- seq(from = 0, to = 0.6, by = 0.05)
N_eq <- (1 - q*E/r)*K
plot(N_eq ~ E, type = "l", xlab = "E", ylab = "N*")
abline(a = 0, b = 0, lty = 2)

N_eq <- (1 - q*E/(r*0.5))*K
lines(N_eq ~ E, type = "l", xlab = "E", ylab = "N*", col = "blue")

N_eq <- (1 - q*E/(r*2))*K
lines(N_eq ~ E, type = "l", xlab = "E", ylab = "N*", lty = 2, col = "blue")

N_eq <- (1 - q*E/(r))*K*0.5
lines(N_eq ~ E, type = "l", xlab = "E", ylab = "N*", col = "red")

N_eq <- (1 - q*E/(r))*K*2
lines(N_eq ~ E, type = "l", xlab = "E", ylab = "N*", lty = 2, col = "red")

legend("topright", 
       c("default: r=0.1, K=40", 
         "r=0.05, K=40", 
         "r=0.2, K=40",
         "r=0.1, K=20",
         "r=0.1, K=80"),
       lty = c(1,1,2,1,2),
       col = c("black", "blue", "blue", "red", "red"),
       bty = "n")
```
  
We can also derive the analytical expression for Effort at MSY:
$$rN(1-N/K) = qEN$$
$$r(1-N/K) = qE$$
$$E = r(1-N/K)/q$$
Insert $N=K/2$ ($E_M$ is effort at MSY!)
$$E_M = r/q*1/2$$
$$E_M = r/2q$$
$K$: affects $N$ at maximum growth, but not the actual growth rate at maximum growth. Thus, it affects yield for any given effort, and $N$* for any given effort, but not the effort at MSY (since there is no $K$ in the above expression!)
    
$r$: does not affect $N$ at maximum population growth rate (remember $K/2$), but does affect the actual value of the growth rate. It affects effort at MSY (see analytical solution above), therefore it affects yield (which has $E$) and $E_M$.
    
This also means that $r$ affects how much harvesting a population can sustain without going extinct, which $K$ does not.

## Part 3: Allee effects
In the logistic model, per capita growth rates decreases linearly with population density. However, in many cases there is a threshold density for positive growth [2], i.e. an Allee effect. In other words, there is a positive relationship between per capita population growth rate and population size at low density We can extend the Schaefer model to account for Allee-effects like this:
$$dN/dt = rN(N/K_0-1)(1-N/K)-qEN$$

* **Question 8**: What is the biological interpretation of the parameter $K_0$?

It is the threshold density for the population to exhibit positive growth rates. Or, when they fall below this threshold, the population collapses because then the per capita growth rate decrease with declining population density.
    
* **Question 9**: Plot the population- and per capita growth rate as a function of $N$ (you can set $K_0$=10). 

Repeat the analysis in the first part, but with the new growth model 
```{r}
## Introduce the parameter K_0
K_0 <- 10

## Population-level growth rate
G <- r * N_vec * ((N_vec/K_0)-1) * (1-(N_vec/K))

plot(G ~ N_vec, type = "l", xlab = "N", ylab = "dN/dt")
abline(a = 0, b = 0, lty = 2)

## Per capita growth
plot(G/N_vec ~ N_vec, type = "l", xlab = "N", ylab = "dN/dt*1/N", ylim = c(0, 0.06))
```

With this growth model we have an equilibrium at $N$*$=0$ (trivial equilibrium). The non-zero equilibria are solutions to the quadratic equation:
$$r(N/K_0-1)(1-N/K)=qE$$
which we obtain by equaling mortality and growth and dividing by $N$. Rather than solving this equation, we can explore it graphically for different fishing efforts the same way we did it before. Or, we can plot $N$ as a function of $qE$:

* **Question 10**: From the above equation, plot the equilibrium density $N$* as a function of $qE$. Compare to Question #6 and explain how and why the equilibria in the two cases differ! 

```{r}
qE <- r*((N_vec/K_0)-1)*(1-N_vec/K)

plot(N_vec ~ qE, type = "l", xlab = "E", ylab = "N*", xlim = c(0, 0.06))
abline(a = 25, b = 0, col = "red")
```

They differ in that in the absence of Allee effects, the density at equilibrium declines monotonically with harvesting as $N$*$=(1-qE/r)K$, wheres here they collapse suddenly when $N=K_0$ (because the population density is below the threshold for positive growth). The equilibrium below the red line is unstable, and the unstable and stable equilibria collide where the red line intersects with the black line.

* **Question 11**: Plot yield as a function of fishing mortality (see Question #7). Locate the stable part of this curve. Which effort maximizes yield?

Yield is defined as: $qEN$, so we can go back to the quadratic equation for the non-zero equilibrium and multiply back $N$. In other words, redo the same plot as in the previous question but also multiply in $N$:
```{r}
## Yield ~ E
plot(qE*N_vec ~ qE, type = "l", xlab = "E", ylab = "Yield", xlim = c(0, 0.06))
abline(a = 0, b = 25, col = "red")
```

Here, increasing effort from 0 increases yield until the curve "folds" (this is a fold or a saddle node bifurcation). At the point where the stable branch starts to move towards lower $E$ the two equilibria (stable and unstable) collide and disappear.

* **Question 12**: What challenges to fisheries management do Allee effects pose?

When Allee effects are present, yield increases with Effort until it collapses, so unlike in the Schaefer model, you will not detect any decline in yield when you are suddenly over harvesting the stock.

## Bonus questions

* **B1**: Find the exact solution of the of the density at MSY in the Schaefer model.

MSY occurs at the population density where growth is fastest (compare your illustration of the growth and mortality parts of the Schaefer equation that you may have done to answer Question 3). We can find this point by taking the derivative of the logistic growth model with respect to $N$ and set it to 0:
$$dN/dt = rN(1-N/K)$$
$$dN/dt = rN-rN^2/K$$
Take the derivative of the growth rate with respect to $N$:
$$dN/dt = r-2rN/K$$
Set to equal 0, solve for N:
$$r-2rN/K = 0$$
$$2Nr/K = r$$
$$2N/K = 1$$
$$2N = K$$
$$N = K/2$$

* **B2**: Solve for the yield at MSY in the Schaefer model.

At MSY, the yield equals the maximum growth in the population. This occurs at $N=K/2$. The growth at $N=K/2$, and therefore the MSY, is $Y=r(K/2)(1-(K/2)/K)$.

Replace $N$ with $K/2$ and substitute $Y=qEN$.
$$qEN = Y = r(K/2)(1-(K/2)/K)$$
$$Y = r(K/2)(1-1/2)$$
$$Y = r(K/2)-r(K/4)$$
$$Y = 2r(K/4)-r(K/4)$$
$$Y = rK/4$$

## References
[1] Schaefer, M. B. 1954. Some aspects of the dynamics of populations important to the management of commercial marine fisheries. *Bulletin of the Inter-American Tropical Tuna Commission*, 1, 25-56.

[2] Perälä, T., & Kuparinen, A. (2017). Detection of Allee effects in marine fishes: Analytical biases generated by data availability and model selection. *Proceedings of the Royal Society B*, 284(1861), 20171284.

## Hints
Question 2: Find expressions of $N$ that leads to fishing mortality being equal to growth.

Question 4: What happens to the relationship between growth and mortality if the system is perturbed slightly to the left and right along the $N$-axis?

Question 5: Repeat Question #4 with higher $E$-values.

Question 6: Create a vector of different $E$-values ranging from 0-0.6 and look at the equation you found in Question #2.

Question 7: Insert the expression for the non-trivial equilibrium in the equation.

Question 9: See first part!

Question 10: Use the same N-vector you have used.

Bonus Q 1: Take the derivative of the logistic equation.

Bonus Q 2: Substitute N with K/2.