1. Writing display and multiple lines of math in RMarkdown and \(\LaTeX\)

There are two ways of writing math in R, either inline or displayed. In-line is used when you want to write math that is simple and short, like if I wanted to write, \(C = \tau r\), as the circumference of a circle in relation to its radius. (See the Tau Manifesto.) ; For example, $C = \tau r$ gives you what I just typed. This is what you did on last week’s problem set. To write displayed, use two dollar signs and write the math in the middle. Displayed math create a line break, which you use for more complex or important math; e.g., \[a^2 + b^2 = c^2.\] To write math in display mode, you use two dollar signs; e.g., $$a^2 + b^2 = c^2$$. For different symbols, letters, etc. see the Resources page on the course website and follow the link to \(\LaTeX\).

1.1 Using the \(\LaTeX\) guide as a reference, take the solution from the continuous-time density-independent growth mode, \(N_t = N_0e^{rt}\), and algebraically solve for \(t\). For every algebraic step, please use a new line within \(\LaTeX\) display mode, like if I wanted to express \(E = mc^2\) in terms of the speed of light, \(c\). Separate lines using \\ and, to align the equations, use \begin{aligned} after the first pair of dollar signed and \end{aligned} before the second pair of dollar signed and put an & before each equals sign to align them, so they look something like this: \[ \begin{aligned} E &= mc^2 \\ \frac{E}{m} &= c^2 \\ c^2 &= \frac{E}{m} \\ c &= \sqrt[2]{\frac{E}{m}}. \end{aligned} \] 1.2. Do the same thing as above, but with the discrete-time density-independent solution, \(N_t = N_0\lambda ^ t\). Also, as a \(\LaTeX\) hint, use $\frac{}{}$ to write a fraction, where what you want in the numerator goes into the first set of braces, and the denominator in the second.

1.3. Finally, we’ll cover the concept of equilibrium later in the semester. But it’s when the system (e.g., population) is not changing; i.e., where the rate of change is 0, or \(0 = dN/dt\). For the logistic equation, \(dN/dt = rN(K - N)/K\), turn it into a per-capita rate of change by dividing both sides by \(N\), then find the equilibrium; i.e., set the per-capita rate of change to 0 and solve for \(N\). For every algebraic step, please use a new line within \(\LaTeX\) display mode, like above.

Answers

1.1. \[ \begin{aligned} N_t &= N_0e^{rt} \\ \frac{N_t}{N_0} &= e^{rt} \\ \log _e \left(\frac{N_t}{N_0}\right) &= rt \\ \frac{\log _e \left(\frac{N_t}{N_0}\right)}{r} &= t \\ t &= \frac{\log _e \left(\frac{N_t}{N_0}\right)}{r} \\ \end{aligned}\]

1.2. \[ \begin{aligned} N_t &= N_0\lambda ^t \\ \frac{N_t}{N_0} &= \lambda ^t \\ \log _\lambda \left(\frac{N_t}{N_0}\right) &= t \\ t &= \log _\lambda \left(\frac{N_t}{N_0}\right) \\ \end{aligned}\]

1.3. \[ \begin{aligned} \frac{dN}{dt} &= rN\left(\frac{K - N}{K}\right) \\ \frac{1}{N}\frac{dN}{dt} &= r\left(\frac{K - N}{K}\right) \\ 0 &= r\left(\frac{K - N}{K}\right) \\ 0 &= r - \frac{r}{K}N \\ \frac{r}{K}N &= r \\ rN &= rK \\ N^\ast = K \end{aligned}\]

2. Making line graphs in R

2.1. Before we proceed in making graphs, I want to note that you are free to use any method that you choose, but I would prefer that you use plot() in base R to make graphs in this course because (1) what you learn in base R is transferable but specialized packages like ggplot2 are not and (2) although there is no plotting in base R that I cannot do (feel free to test this out) I cannot say the same about ggplot2 and the rest of the tidy universe. (If you don’t know what any of this means, feel free to ignore it.)
Onward: making plots in R is fundamentally simple. The three steps you need for 2-dimensional graphs (i.e., have a horizontal and vertical axis) are: (1) create a vector of numbers representing your \(x\) variable, (2) plug \(x\) into an equation, and (3) plot it using plot(x = , y = ) where you plug in the names of your \(x\) and \(y\) variables. That’s really it.
Let’s try this by plotting a quadratic function, \(y = 0.5x^2 + 0.1x - 2\). To do this, I would create a sequence of values of a variable \(x\), like x_vec <- seq(from = -5, to = 5, length.out = 10) and make \(y\) a function of \(x\), like y_vec <- 0.5*x_vec^2 + 0.1*x_vec - 2. When I plot, plot(x = x_vec, y = y_vec, type = "l", las = 1), I get

Make this plot, but label the horizontal and vertical axes \(x\) and \(y\), respectively.

2.2. I want you to make a couple more plots. First, choose one to plot either of the equations from 1.1. or 1.2., with the \(x\) variable being \(t\) and the \(y\) variable being \(N_t\). Make sure you make the length.out argument within the sequence of \(x\) values longer than just 10: try whatever makes it sufficiently smooth like 100, 1,000, or 10,000. Second, I want you to add 3 additional curves onto a new graph of which ever equation you chose to plot. For each of the three curves, vary either the \(r\) or \(\lambda\) parameter. For instance, if I use the quadratic function from above, I would plot it and change, say, the coefficient of \(x^2\). To add this to your plot, use the lines() function. Change the color of the curve just to tell the difference between them; e.g.,

x_vec <- seq(from = -5, to = 5, length.out = 100)

y_vec <- 0.5*x_vec^2 + 0.1*x_vec - 2
y_vec1 <- 1*x_vec^2 + 0.1*x_vec - 2
y_vec2 <- 2*x_vec^2 + 0.1*x_vec - 2
y_vec3 <- 4*x_vec^2 + 0.1*x_vec - 2

par(mar = c(4.5, 4.5, 0.1, 0.1)) # sets the whitespace on the bottom, left, top, right
plot(x = x_vec, y = y_vec, type = "l", las = 1, xlab = "x", ylab = "y")
  lines(x = x_vec, y = y_vec1, col = "red")
  lines(x = x_vec, y = y_vec2, col = "#0000FF") # You can use hex codes
  lines(x = x_vec, y = y_vec3, col = "gold")

2.1

par(mar = c(4.5, 4.5, 0.1, 0.1))
x_vec <- seq(from = -5, to = 5, length.out = 10)
y_vec <- 0.5*x_vec^2 + 0.1*x_vec - 2
plot(x = x_vec, y = y_vec, type = "l", las = 1)

2.2.

par(mar = c(4.5, 4.5, 0.1, 0.1))
t_vec <- seq(from = 0, to = 10, length.out = 1e4)
r1 <- 0.1
r2 <- 0.2
r3 <- 0.3
N0 <- 0.5
y_vec1 <- N0*exp(r1*t_vec)
y_vec2 <- N0*exp(r2*t_vec)
y_vec3 <- N0*exp(r3*t_vec)
plot(x = t_vec, y = y_vec1, type = "l", las = 1, xlab = "Time (t)", ylab= "Density (N)", ylim = c(0, 10))
 lines(x = t_vec, y = y_vec2, col = "skyblue")
 lines(x = t_vec, y = y_vec3, col = "tomato")

3. Using ode() and some density-dependent models

In lab in week 2 I showed you the numerical solver you will be using for just about everything in the class, and in lab 3 you started to use it. It is the function called ode() in the deSolve package.
3.1. First, I want you to go to ode()’s documentation (use ?ode in the console) and copy the the code from example 1 through the line when out is created. Plot the predator and prey time series. Two ways to plot the columns of the ode() output saves as out could be, using the time variable as an example, out[,"time"] or out[,1]. Make sure to adjust the vertical axis in plot to show all the data.

3.2. Take the code above and modify it so that way you have the logistic equation, \(dN/dt = rN(K - N)/K\). This is what you did in lab in week 3. Give \(r\) and \(K\) whatever values you want and start at some low value of \(N\), say 0.1 Please plot the time series. Then, again numerically ingrate across a few different initial population sizes, say where \(N\) starts at 1/2 of the \(K\) you chose,where \(N\) starts just slightly about the \(K\) you chose, and where \(N\) starts at a value that’s 25% higher than the \(K\) you chose. Plot these curves all on the same graph.

3.3. Last, recreate parts of the Figure 5.12 from Stevens link. Choose a value of \(r_d\) that gives you asymptotic stability, oscillations (damped or perpetual), and chaos. Use the logistic equation he uses that is found 3 chunks above the plots where he creates dlogistic. Use the numerical solver like in 3.2, but in the ode() function, write the argument method = "iteration". Use the logistic equation he uses in dlogistic, give the initial abundance a value of 10, set alpha to 0.01, and create a time sequence from 0 to 20.

Please create your own plots using base plot(). To make a 2-row-by-3-column figure, in a line above plot() add a line for a function that adjusts plotting parameters, par(). If you use the argument mfrow (stands for matrix of figures by rows; there’s also an mfcol) equal to the rows and columns, c(2, 3), then the next 6 plots you make will be loaded in with the first 3 in the top row and the last 3 in the bottom row.

3.1.

library(package = "deSolve")

LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
    return(list(c(dPrey, dPredator)))
  })
}

pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
           rGrow  = 1.0,    # /day, growth rate of prey
           rMort  = 0.2 ,   # /day, mortality rate of predator
           assEff = 0.5,    # -, assimilation efficiency
           K      = 10)     # mmol/m3, carrying capacity
yini  <- c(Prey = 1, Predator = 2)
times <- seq(from = 0, to = 200, by = 1)

out   <- ode(y = yini, times = times, func = LVmod, parms = pars)

plot(x = out[,1], y = out[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 5))
  lines(x = out[,1], y = out[,3], col = "red")
  legend(x = "bottomright", bty = "n", legend = c("Prey", "Predator"), col = c("black", "red"), lty = 1)

3.2.

logistic <- function(time, state, pars) {
  with(as.list(c(state, pars)), {
    dN <- r*N*(K - N)/K
    return(list(c(dN)))
  })
}


yini1  <- c(N = 0.1)
yini2  <- c(N = 5)
yini3  <- c(N = 9)
yini4  <- c(N = 12.5)

pars  <- c(r = 0.5, K = 10)
times <- seq(from = 0, to = 25, by = 0.1)

out1   <- ode(y = yini1, times = times, func = logistic, parms = pars)
out2   <- ode(y = yini2, times = times, func = logistic, parms = pars)
out3   <- ode(y = yini3, times = times, func = logistic, parms = pars)
out4   <- ode(y = yini4, times = times, func = logistic, parms = pars)

plot(x = out1[,1], y = out1[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 14))
  lines(x = out2[,1], y = out2[,2], col = "gold")
  lines(x = out3[,1], y = out3[,2], col = "skyblue")
  lines(x = out4[,1], y = out4[,2], col = "tomato")

3.3

dlogistic <- function(time, state, pars) {
  with(as.list(c(state, pars)), {
    dN <- N + rd*N*(1 - alpha*N) # Sorry, using "iteration" means you need to put a + N in your model
    return(list(c(dN)))
  })
}

yini  <- c(N = 10)
times <- 0:20

rd_vec <- seq(from = 1.3, to = 2.8, by = 0.3)
pars1  <- c(rd = rd_vec[1], alpha = 0.01)
pars2  <- c(rd = rd_vec[2], alpha = 0.01)
pars3  <- c(rd = rd_vec[3], alpha = 0.01)
pars4  <- c(rd = rd_vec[4], alpha = 0.01)
pars5  <- c(rd = rd_vec[5], alpha = 0.01)
pars6  <- c(rd = rd_vec[6], alpha = 0.01)

out1   <- ode(y = yini, times = times, func = dlogistic, parms = pars1, method = "iteration")
out2   <- ode(y = yini, times = times, func = dlogistic, parms = pars2, method = "iteration")
out3   <- ode(y = yini, times = times, func = dlogistic, parms = pars3, method = "iteration")
out4   <- ode(y = yini, times = times, func = dlogistic, parms = pars4, method = "iteration")
out5   <- ode(y = yini, times = times, func = dlogistic, parms = pars5, method = "iteration")
out6   <- ode(y = yini, times = times, func = dlogistic, parms = pars6, method = "iteration")

par(mfrow = c(2, 3), mar = c(4.5, 4, 1, 0.5)) # mar changes the amount of space below, left, top, and right of each plot
plot(x = out1[,1], y = out1[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 125), lwd = 1.5)
plot(x = out2[,1], y = out2[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 125), lwd = 1.5)
plot(x = out3[,1], y = out3[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 125), lwd = 1.5)
plot(x = out4[,1], y = out4[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 125), lwd = 1.5)
plot(x = out5[,1], y = out5[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 125), lwd = 1.5)
plot(x = out6[,1], y = out6[,2], xlab = "Time", ylab = "Density", type = "l", las = 1, ylim = c(0, 125), lwd = 1.5)

4. Discrete-time density-dependent models

At the end of chapter 1 in Vandermeet and Goldberg you’re introduces to additional commonly-used discrete-time density-dependent models, the Beverton-Holt model (eqn. 28), the Hassell model (eqn. 29), and the Ricker model (eqn. 30). I’d like for you to visually/graphically analyze these models by plotting them and modifying their parameters.

4.1. First, let’s compare among the models. Let’s set all the growth terms the same, by setting \(\lambda = r = 1\), and do the same to the crowding terms, \(\alpha = b~\text{in Ricker} = 0.1\). The last remaining parameter is \(b\) in the Hassell model, and let’s just set it to 1. Create a time sequence, but remember that we are in discrete-time, so we are taking steps by 1. The easiest way to write this in R is 0:10 if you wanted to iterate/solve across 10 time steps. Last, start with \(N_0\) at 0.1. Plot all three using points on the same graph and make sure your axes are showing all the data.

4.2. Please look over the graphs you just made and comment on what you observe given what we did with the parameters.

4.3. Next, let’s vary the density-dependent terms in the model. On three graphs, please vary \(\alpha\) in the Beverton-Holt model in one, vary \(b\) in the Hassell model (even though \(\alpha\) is the density-dependent term, we vary it with the Beverton-Holt model that is analogous to Hassell), and \(b\) in the Ricker model. For each of those parameters, chose 3 values and plot them on each graph. Try to use a parameter value that gives you at least one different type of qualitative behavior (e.g., oscillations, chaos).

4.1.

# Write the functions
 bevholt <- function(t, y, p) {
   with(as.list(c(y, p)), {
     Nt1 <- lambda*N / (1 + alpha*N)
     return(list(c(Nt1)))
   })
 }
 
 hassell <- function(t, y, p) {
   with(as.list(c(y, p)), {
     Nt1 <- lambda*N / (1 + alpha*N)^b
     return(list(c(Nt1)))
   })
 }
 
 ricker <- function(t, y, p) {
   with(as.list(c(y, p)), {
     Nt1 <- N*exp(r*(1 - b*N))
     return(list(c(Nt1)))
   })
 }

# Set up parameters
 N0 <- c(N = 0.1)
 time_seq <- seq(from = 0, to = 25, by = 1)
 bevholt_parms <- c(lambda = 1, alpha = 0.1)
 hassell_parms <- c(lambda = 1, alpha = 0.1, b = 1)
 ricker_parms <- c(r = 1, b = 0.1)

# Run ode
 beveholt_out <- ode(y = N0, times = time_seq, func = bevholt, parms = bevholt_parms, method = "iteration")
 hassell_out <- ode(y = N0, times = time_seq, func = hassell, parms =hassell_parms, method = "iteration")
 ricker_out <- ode(y = N0, times = time_seq, func = ricker, parms = ricker_parms, method = "iteration")

# Plot
 plot(x = beveholt_out[,1], y = beveholt_out[,2], xlab = "Time (steps)", ylab = "Density", pch = 19, las = 1, ylim = c(0, 10))
   lines(x = beveholt_out[,1], y = beveholt_out[,2], col = "black")
   points(x = hassell_out[,1], y = hassell_out[, 2], pch = 1, col = "red")
   lines(x = hassell_out[,1], y = hassell_out[,2], col = "red", lty = 2)
   points(x = ricker_out[,1], y = ricker_out[, 2], pch = 19, col = "blue")
   lines(x = ricker_out[,1], y = ricker_out[,2], col = "blue")
   legend(x = "right", pch = c(19, 1, 19), col = c("black", "red", "blue"), legend = c("Beverton-Holt", "Hassell", "Ricker"), bty = "n")

4.2. Ever with the same values for the growth parameters and crowding (i.e., density-dependence) parameters, the models yield different numeric/qunatitative results. Also, with \(b = 1\) in the Hassell model, it becomes the Beverton-Holt, and we see this equivalence graphically.

4.3.

 N0 <- c(N = 0.1)
 time_seq <- seq(from = 0, to = 150, by = 1)

# Beverton-holt, varying alpha
 bevholt_parms1 <- c(lambda = 1.1, alpha = 0.01)
 bevholt_parms2 <- c(lambda = 1.1, alpha = 0.05)
 bevholt_parms3 <- c(lambda = 1.1, alpha = 0.1)
# Hassel, varying b
 hassell_parms1 <- c(lambda = 1.1, alpha = 0.1, b = 0.25)
 hassell_parms2 <- c(lambda = 1.1, alpha = 0.1, b = 1)
 hassell_parms3 <- c(lambda = 1.1, alpha = 0.1, b = 2)
# Ricker, varying b
 ricker_parms1 <- c(r = 0.1, b = 0.1)
 ricker_parms2 <- c(r = 0.1, b = 0.5)
 ricker_parms3 <- c(r = 0.1, b = 1)

# Beverton-holt, varying alpha
 beveholt_out1 <- ode(y = N0, times = time_seq, func = bevholt, parms = bevholt_parms1, method = "iteration")
 beveholt_out2 <- ode(y = N0, times = time_seq, func = bevholt, parms = bevholt_parms2, method = "iteration")
 beveholt_out3 <- ode(y = N0, times = time_seq, func = bevholt, parms = bevholt_parms3, method = "iteration")
# Hassel, varying b
 hassell_out1 <- ode(y = N0, times = time_seq, func = hassell, parms =hassell_parms1, method = "iteration")
 hassell_out2 <- ode(y = N0, times = time_seq, func = hassell, parms =hassell_parms2, method = "iteration")
 hassell_out3 <- ode(y = N0, times = time_seq, func = hassell, parms =hassell_parms3, method = "iteration")
# Ricker, varying b
 ricker_out1 <- ode(y = N0, times = time_seq, func = ricker, parms = ricker_parms1, method = "iteration")
 ricker_out2 <- ode(y = N0, times = time_seq, func = ricker, parms = ricker_parms2, method = "iteration")
 ricker_out3 <- ode(y = N0, times = time_seq, func = ricker, parms = ricker_parms3, method = "iteration")

par(mfrow = c(1, 3), mar = c(4.5, 4, 1, 0.5))
# Beverton-holt, varying alpha
 plot(x = beveholt_out1[,1], y = beveholt_out1[,2], xlab = "Time (steps)", ylab = "Density", pch = 19, las = 1, ylim = c(0, 10), col = "red")
 points(x = beveholt_out2[,1], y = beveholt_out2[,2], pch = 19, col = "purple")
 points(x = beveholt_out3[,1], y = beveholt_out3[,2], pch = 19, col = "blue")
 legend(x = "right", bty = "n", pch = 19, col = c("red", "purple", "blue"), legend = c(0.01, 0.05, 0.1), title = expression(alpha))
# Hassel, varying b
 plot(x = hassell_out1[,1], y = hassell_out1[,2], xlab = "Time (steps)", ylab = "Density", pch = 19, las = 1, ylim = c(0, 10), col = "red")
 points(x = hassell_out2[,1], y = hassell_out2[,2], pch = 19, col = "purple")
 points(x = hassell_out3[,1], y = hassell_out3[,2], pch = 19, col = "blue")
 legend(x = "topright", bty = "n", pch = 19, col = c("red", "purple", "blue"), legend = c(0.25, 1, 2), title = expression(b))
# Ricker, varying b
 plot(x = ricker_out1[,1], y = ricker_out1[,2], xlab = "Time (steps)", ylab = "Density", pch = 19, las = 1, ylim = c(0, 10), col = "red")
 points(x = ricker_out2[,1], y = ricker_out2[,2], pch = 19, col = "purple")
 points(x = ricker_out3[,1], y = ricker_out3[,2], pch = 19, col = "blue")
 legend(x = "right", bty = "n", pch = 19, col = c("red", "purple", "blue"), legend = c(0.5, 0.1, 1.0), title = expression(b))