Here the problem numbers refer to the chapter numbers in the book.

4.1. The integrated form of the logistic equation is \[N(t) = \frac{KN(0)}{(K - N(0))\exp(-rt) + N(0))}\] Use this equation to generate a time series of \(N\) versus \(t\). Use parameters \(r = 1.5\), \(K = 100\) and construct two time series, one with \(N(0) = 10\) and one with \(N(0)=190\). (Use a time frame from 11 to 5 with intervals of 0.1.)

The “integrated form” is the same as a solution to the logistic model.

times_4.1 <- seq(from = 0, to = 5, by = 0.1)
r <- 1.5
K <- 100
N0_4.1_10 <- 10
Nt_4.1_10 <- (K*N0_4.1_10)/((K - N0_4.1_10)*exp(-r*times_4.1) + N0_4.1_10)

N0_4.1_190 <- 190
Nt_4.1_190 <- (K*N0_4.1_190)/((K - N0_4.1_190)*exp(-r*times_4.1) + N0_4.1_190)

par(mfrow = c(1, 2))
plot(x = times_4.1, y = Nt_4.1_10, las = 1, ylab = "Density (N)", xlab = "Time", type = "l", lwd = 2)
plot(x = times_4.1, y = Nt_4.1_190, las = 1, ylab = "Density (N)", xlab = "Time", type = "l", lwd = 2)

4.2. Repeat exercise 4.1 wth \(r= -0.1\). (Use a time frame of 0 to 100 with intervals of at least 1.0.)

Please comment on the difference between 4.1. and 4.2., and make sure to check in with your guts/intuition before and after.

times_4.2 <- seq(from = 0, to = 100, by = 0.1)
r <- -0.1
K <- 100
N0_4.2_10 <- 99
Nt_4.2_10 <- (K*N0_4.2_10)/((K - N0_4.2_10)*exp(-r*times_4.2) + N0_4.2_10)

N0_4.2_190 <- 90
Nt_4.2_190 <- (K*N0_4.2_190)/((K - N0_4.2_190)*exp(-r*times_4.2) + N0_4.2_190)

par(mfrow = c(1, 2))
plot(x = times_4.2, y = Nt_4.2_10, las = 1, ylab = "Density (N)", xlab = "Time", type = "l", lwd = 2)
plot(x = times_4.2, y = Nt_4.2_190, las = 1, ylab = "Density (N)", xlab = "Time", type = "l", lwd = 2)

When death rates are greater than birth rates and if the population is below its positive equilibrium (“carrying capacity”), then the population shrinks to 0. When death rates are greater than birth rates and if the population is above its positive equilibrium, then the population grows (strange—right?) then becomes negative (this is getting weirder), and lastly grows again until it asymptotes at 0.

4.3. Recall the population model of the Ricker equation from chapter 1, \[N_{t+1} = rN_t\exp(1 - bN_t)\]. Let \(r = 2.5\) and \(b = 2.5\) and project the population 45 times units, beginning with a starting population size of 0.01. What is the pattern of the time series Let \(r = 4\) and \(b=0.7\) and project the time series 45 times, beginning with a starting population of 3.5. Now what is the pattern of the time series?

No additional comments aside from use ode(..., method = "iteration") this time.

Ricker <- function(t, N, parms) {
 with(as.list(c(N, parms)), {
 Nt <- r*N*exp(1 - b*N)
 return(list(Nt))
 })
}

N0_4.3_a <- 0.01
parameters_a <- c(r = 2.5, b = 2.5)
times_4.3_a <- 0:45
out_4.3_a <- ode(y = N0_4.3_a, times = times_4.3_a, func = Ricker, parms = parameters_a, method = "iteration")
plot(x = out_4.3_a[,1], y = out_4.3_a[,2], las = 1, pch = 16, cex = 1.25, xlab = "Time steps", ylab = "Density (N)")
lines(x = out_4.3_a[,1], y = out_4.3_a[,2])

N0_4.3_b <- 3.5
parameters_b <- c(r = 4, b = 0.7)
times_4.3_b <- 0:45
out_4.3_b <- ode(y = N0_4.3_b, times = times_4.3_b, func = Ricker, parms = parameters_b, method = "iteration")
plot(x = out_4.3_b[,1], y = out_4.3_b[,2], las = 1, pch = 16, cex = 1.25, xlab = "Time steps", ylab = "Density (N)")
lines(x = out_4.3_b[,1], y = out_4.3_b[,2])

The pattern of the time series for the first part is that the popluation fluxuates towards around a density just under 0.8, with the fluctuations dampening over time. In the second part, the population fluctuations increase over time, but eventually max and min out at just over 5 and 1, respectively.

4.4. A simple model of predator (\(P\)) and prey (\(V\), for victim) omteractomg om doscrete time is, for the prey, \[V_{t+1} = bV_t\left(\frac{K - V_t}{K}\right)\exp \left(-aP_t\right)\] and for the predator, \[P_{t+1} = cV_t\left(1 - \exp \left(-aP_t\right)\right).\] This model will be developed more fully in chapter 6. For now just examine the time series it generates. Use the parameters \(a = 0.1\), \(b = 1.5\), \(K = 40\), and \(c = 1.5\). Make a graph of the two species over time (from 0 to 200), and plot the two species on a graph of predator versus prey, showing the direction of change with arrows

I’ve done this problem for you, without arrows.

# Load package deSolve
    library(package = "deSolve")

# Function
    PV_disc <- function(t, state, parms) {
        with(as.list(c(state, parms)), {
          Vt1 <- b*V*((K - V) / K)*exp(-a*P)
          Pt1 <- c*V*(1 - exp(-a*P))
          return(list(c(Vt1, Pt1)))
          })
    }

# Parameters
    parm_vec <- c(b = 1.5, K = 40, a = 0.1, c = 1.5)

# Simulation conditions
    steps_vec <- seq(from = 0, to = 200, by = 1)
    init_Ns <- c(V = 12.18, P = 0.06)

# Simulation (discrete time, so method = iteration)
    odeOut_4.4 <- ode(y = init_Ns, times = steps_vec, func = PV_disc, parms = parm_vec, method = "iteration")

# Plot results
    par(mfrow = c(1, 2), mar = c(4.5, 4.5, 0.5, 0.5))

    # Time series
    plot(x = steps_vec, y = odeOut_4.4[,2], type = "l", col = "blue", xlab = "Step", ylab = "Abundance", ylim = c(0, max(odeOut_4.4[,2:3])), las = 1)
        points(x = steps_vec, y = odeOut_4.4[,2], pch = 16, col = "blue", cex = 0.5)
    lines(x = steps_vec, y = odeOut_4.4[,3], col = "red")
        points(x = steps_vec, y = odeOut_4.4[,3], pch = 16, col = "red", cex = 0.5)

    # State space
    plot(x = odeOut_4.4[,2], y = odeOut_4.4[,3], type = "l", xlab = "Prey abundance", ylab = "Predator abundance", las = 1)
        points(x = odeOut_4.4[,2], y = odeOut_4.4[,3], pch = 16, col = "black", cex = 0.5)

4.5.

A even simpler model (again, to be developed in chapter 6) is \[\begin{align} V_{t+1} &= bV_t - mV_tP_t \\ P_{t+1} &= b'V_tP_t - m'P_t. \end{align}\] Repeat exercise 4.4 with this model and the parameters \(b = 1.1\), \(m = 0.8\), \(b' = 0.5\), and \(m'=0.002\).

Same as 4.4, but change the formula and parameter values. For initial conditions, use densities of init.Ns <- c(V = 1.996, P = 0.125)

PV.disc.4.5 <- function(t, state, parms) {
 b <- 1.1
 m = 0.8
 b1 = 0.5
 m1 = 0.002
 V <- state["V"]
 P <- state["P"]
 Vt1 <- b*V - m*V*P
 Pt1 <- b1*V*P - m1*P
 return(list(c(Vt1, Pt1)))
}

steps.vec <- seq(from = 0, to = 200, by = 1)
init.Ns <- c(V = 1.996, P = 0.125)

odeOut.2.15 <- ode(y = init.Ns, times = steps.vec, func = PV.disc.4.5, parms = NA, method = "iteration")

pos.steps <- 1:min(which(odeOut.2.15[,2:3] < 0, arr.ind = T)[,1])

plot(x = steps.vec[pos.steps], y = odeOut.2.15[pos.steps, 2], type = "l", col = "blue", ylab = "Step", xlab = "Abundance", ylim = , las = 1)
lines(x = steps.vec[pos.steps], y = odeOut.2.15[pos.steps, 3], col = "red")

plot(x = odeOut.2.15[pos.steps, 2], y = odeOut.2.15[pos.steps, 3], type = "l", xlab = "Prey abundance", ylab = "Predator abundance", las = 1)

4.9. Assume that a population is growing according to the logistic equation. To make things simple, presume that the value of both \(r\) and \(K\) is 1.0 (i.e., we represent the population as varying between 0 and 1). The equilibrium value of that population will be \[0 = N − N^2,\] which is a quadratic equation and has two roots, which are, by inspection, 0 and 1. Now suppose that a manager decides to impose a fixed harvesting rate on the population such that a constant number of individuals will be removed each year. A sensible model for this situation would be \[\frac{\mathrm{d}N}{\mathrm{d}t} = N(1 − N) − N_F\] where \(N_F\) is the fixed number (actually the proportion) of individuals removed. Plot the derivative versus N for the following values of \(N_F\): 0, 0.25, and 0.5. Also, directly solve for the roots (using the quadratic formula; remember, \(N_F\) is a constant). What do you conclude?

No additional comments.

N_vec <- seq(from = 0, to = 1.1, length.out = 1000)
NF1 <- 0
NF2 <- 0.25
NF3 <- 0.5
dNdt1 <- N_vec*(1 - N_vec) - NF1
dNdt2 <- N_vec*(1 - N_vec) -  NF2
dNdt3 <-N_vec*(1 - N_vec) - NF3

par(mfrow = c(1, 3), mar = c(4.5, 4.5, 0.25, 0.25))
plot(x = N_vec, y = dNdt1, las = 1, xlab = "Density (N)", ylab = "dN/dt", type = "n", ylim = c(-0.5, 0.25))
 # You can draw lines by specifying the beginning and end of the x and y cooridnates
 segments(x0 = -1, x1 = 2, y0 = 0, y1 = 0, col = "grey25")
 segments(x0 = 0, x1 = 0, y0 = -1, y1 = 1, col = "grey25")
 lines(x = N_vec, y = dNdt1, col = "red")
plot(x = N_vec, y = dNdt2, las = 1, xlab = "Density (N)", ylab = "dN/dt", type = "n", ylim = c(-0.5, 0.25))
 segments(x0 = -1, x1 = 2, y0 = 0, y1 = 0, col = "grey25")
 segments(x0 = 0, x1 = 0, y0 = -1, y1 = 1, col = "grey25")
 lines(x = N_vec, y = dNdt2, col = "red")
plot(x = N_vec, y = dNdt3, las = 1, xlab = "Density (N)", ylab = "dN/dt", type = "n", ylim = c(-0.5, 0.25))
 segments(x0 = -1, x1 = 2, y0 = 0, y1 = 0, col = "grey25")
 segments(x0 = 0, x1 = 0, y0 = -1, y1 = 1, col = "grey25")
 lines(x = N_vec, y = dNdt3, col = "red")

The roots are \[N^* = \frac{1}{2}\left(1 \pm \sqrt{4N_F + 1}\right).\] This means that there are roots—i.e., equilibria—that are dependent on the value of \(N_F\) we can see that in from the graphs: if \(N_F\) is large then there is no real root/poistive equilibirum.