Libraries

library(plotly)

Discrete time dynamical models

Exponential growth

\[\begin{align}&N_{t+1} = N_t + bN_t - dN_t \\ &N_{t+1} = N_t + \left( b - d \right) N_t \\ &N_{t+1} = N_t + rN_t \\ &N_{t+1} = \left(1 + r \right) N_t \\ &N_{t+1} = \lambda N_t \end{align}\tag{1}\]

Here \(N_{t+1}\) is the new population, \(N_t\) is the old population, \(b\) is the per capita birth rate, \(d\) is the per capita death rate, and \(\lambda\) the growth rate such that \(\lambda = b - d + 1\).

We can solve this difference equation numerically by simulation with for loops and functions.

generations <- 10  # Discrete time steps
N <- numeric(generations)  # Vector to hold the calculated values
lambda <- 2.1
N[1] <- 3  # Initial condition using indexing, starting at 1 (which is t = 0)
str(generations)
##  num 10
str(N)
##  num [1:10] 3 0 0 0 0 0 0 0 0 0
str(lambda)
##  num 2.1
for (t in 1:(generations - 1)){
  N[t+1] <- lambda * N[t]
}
N
##  [1]    3.0000    6.3000   13.2300   27.7830   58.3443  122.5230  257.2984
##  [8]  540.3266 1134.6858 2382.8401

We can visualize this using Plotly.

exp_mod_plot <- plot_ly(x = seq(from = 0, to = 9, by = 1),
                        y = N,
                        type = "scatter",
                        mode = "lines+markers",
                        marker = list(size = 12),
                        name = "Population") %>% 
  layout(title = "Discrete exponetial growth",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Population"))
exp_mod_plot

A function can be created to automate the calculation.

pop.growth <- function(b, d, n0, gens){
  N <- c(n0, numeric(generations - 1))
  for (t in 1:(generations - 1)){
    N[t+1] <- N[t] + b * N[t] - d * N[t]
  }
  return(N)
}
population_vals <- pop.growth(2.0, 0.9, 3, 10)
population_vals
##  [1]    3.0000    6.3000   13.2300   27.7830   58.3443  122.5230  257.2984
##  [8]  540.3266 1134.6858 2382.8401

This problem can also be solved explicitely as shown in (2).

\[\begin{align}&N_1 = \lambda N_0 \\ &N_2 = \lambda N_1 = \lambda \left( \lambda N_0 \right) = \lambda^2 N_0 \\ \vdots \\ &N_t = \lambda^t N_0 \end{align}\tag{2}\]

Now for a quick use of (2).

lambda <- 2.1
t <- 0:9
N0 <- 3
Nt <- (lambda^t) * N0  # Broadcasting
Nt
##  [1]    3.0000    6.3000   13.2300   27.7830   58.3443  122.5230  257.2984
##  [8]  540.3266 1134.6858 2382.8401

We have used the object t differently.

str(t)
##  int [1:10] 0 1 2 3 4 5 6 7 8 9

This time we use the built-in plotting function to visualize our data.

plot(t, Nt, type = "o", main = "Discrete exponential population growth", xlab = "Time", ylab = "Population count", las = 1)

Logistic growth

This is a model of density-dependent growth, shown in equation (3).

\[N_{t+1} = N_t + rN_t \left( 1 - \frac{N_t}{K} \right)\tag{3}\]

Here \(N_i\) is as before, \(r\) is the difference between the per capita birth and death rates, and \(K\) is the carrying capacity. We can create a function with a for loop as before..

log_growth <- function(K, r, N0, gens){
  N <- c(N0, numeric(gens - 1))
  for (t in 1:(gens - 1)){
    N[t+1] <- N[t] * r * (1 - (N[t] / K))
  }
  return(N)
}

Let’s create a vector object to hold values over \(30\) generations given some input values and visualise the data.

log_pop_values <- log_growth(K = 100, r = 1.5, N0 = 10, gens = 30)
log_mod_plot <- plot_ly(x = seq(from = 0, to = 29, by = 1),
                        y = log_pop_values,
                        type = "scatter",
                        mode = "lines+markers",
                        marker = list(size = 12),
                        name = "Population") %>% 
  layout(title = "Discrete exponetial growth",
         xaxis = list(title = "Time"),
         yaxis = list(title = "Population"))
log_mod_plot

Other single-species discrete time-step population models

Ricker model

\[x_{n+1} = a x_n e^{-bx_n}\tag{4}\]

Hassel model

\[x_{n+1} = \frac{ax_n}{\left( 1 + bx_n \right)^c}\tag{5}\]

Beverton-Holt model

\[x_{n+1} = \frac{ax_n}{1 + bx_n}\tag{6}\]

Maynard Smith and Slatkin model

\[x_{n+1} = \frac{ax_n}{1+bx_n^c}\tag{7}\]

Skellam model

\[x_{n+1} = \frac{a}{b} \left( 1 - e^{-{b}{x}_{n}} \right)\tag{8}\]

Deterministic chaos

Here we look at modeling an intrinsic factor that introduces determined chaos. It can be modelled as shown in (9) (canonical form of the discrete-time logistic function).

\[X_{t+1} = rX_t \left( 1 - X_t \right)\tag{9}\]

Here \(X_{t+1}\) is the next ratio of the existing population to the maximum possible population, \(X_t\) is the previois ratio of the existing population over the maximum possible population, and \(r\) is the difference between the per capita birth and death rates.

With increasing rates of \(r\) we start at a 1-point cycle and move up through 2-point cycles, etc. until we reach chaos.

Let’s start by creating a function for the canocial form and then create some models with different valuesof \(r\).

can_log <- function(r, X0, gens){
  X <- c(X0, numeric(gens - 1))
  for (t in 1:(gens - 1)){
    X[t+1] <- r * X[t] * (1 - X[t])
  }
  return(X)
}
x0_val <- 0.9
gen_val <- 250
Xt_1 <- can_log(r = 2, X0 = x0_val, gens = gen_val)
Xt_2 <- can_log(r = 3.5, X0 = x0_val, gens = gen_val)
Xt_3 <- can_log(r = 3.9, X0 = x0_val, gens = gen_val) 
x_times <- seq(from = 0, to = 249, by = 1)
chaos_mod_plot_1 <- plot_ly(x = x_times,
                          y = Xt_1,
                          name = "r = 2",
                          type = "scatter",
                          mode = "markers",
                          marker = list(size = 12))
chaos_mod_plot_2 <- plot_ly(x = x_times,
                            y = Xt_2,
                            name = "r = 3.5",
                            mode = "markers",
                            marker = list(size = 12),
                            type = "scatter")
chaos_mod_plot_3 <- plot_ly(x = x_times,
                            y = Xt_3,
                            name = "r = 3.9",
                            mode = "markers",
                            marker = list(size = 12),
                            type = "scatter") 
subplot(chaos_mod_plot_1, chaos_mod_plot_2, chaos_mod_plot_3, nrows = 3, shareX = T)