Week 7 - Lotka-Volterra competition model - Population dynamics

In this lab, we are going to analyze the two-species Lotka-Volterra competition model numerically and visualize the population dynamics under different parameter settings.

library(ggplot2)
library(tidyverse)
library(deSolve)

LV_model <- function(r1 = 1.4, r2 = 1.2, a11 = 1/200, a21 = 1/400, a22 = 1/200, a12 = 1/300, N1_0 = 10, N2_0 = 10) {

  ### Model specification
  LV <- function(times, state, parms) {
    with(as.list(c(state, parms)), {
      dN1_dt = N1 * (r1 - a11*N1 - a12*N2)
      dN2_dt = N2 * (r2 - a22*N2 - a21*N1)
      return(list(c(dN1_dt, dN2_dt)))
    })
  }

  ### Model parameters
  times <- seq(0, 100, by = 0.1)
  state <- c(N1 = N1_0, N2 = N2_0)
  parms <- c(r1 = r1, r2 = r2, a11 = a11, a21 = a21, a22 = a22, a12 = a12)

  ### Model application
  pop_size <- ode(func = LV, times = times, y = state, parms = parms)

  ### Visualize the population dynamics
  pop_size %>%
    as.data.frame() %>%
    gather(key = "Species", value = "pop_size", -time) %>%
    ggplot(aes(x = time, y = pop_size, color = Species)) +
    geom_line(size = 1.5)

}

  ### Different parameter settings
  ## N1_0 = 200 and N2_0 = 5
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/200, a21 = 1/100, a22 = 1/100, a12 = 1/200, N1_0 = 200, N2_0 = 5)  # N1 wins

  ## N1_0 = 5 and N2_0 = 200
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/200, a21 = 1/100, a22 = 1/100, a12 = 1/200, N1_0 = 10, N2_0 = 200)  # N1 wins

  ## N1_0 = 200 and N2_0 = 5
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/100, a21 = 1/200, a22 = 1/200, a12 = 1/100, N1_0 = 200, N2_0 = 5)  # N2 wins

  ## N1_0 = 5 and N2_0 = 200
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/100, a21 = 1/200, a22 = 1/200, a12 = 1/100, N1_0 = 5, N2_0 = 200)  # N2 wins

  ## N1_0 = 200 and N2_0 = 5
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/100, a21 = 1/200, a22 = 1/100, a12 = 1/300, N1_0 = 200, N2_0 = 5)  # stable coexistence

  ## N1_0 = 5 and N2_0 = 200
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/100, a21 = 1/200, a22 = 1/100, a12 = 1/300, N1_0 = 5, N2_0 = 200)  # stable coexistence

  ## N1_0 = 200 and N2_0 = 150
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/200, a21 = 1/100, a22 = 1/200, a12 = 1/100, N1_0 = 200, N2_0 = 150)  # priority effect (N1 wins)

  ## N1_0 = 150 and N2_0 = 200
  LV_model(r1 = 1.2, r2 = 1.2, a11 = 1/200, a21 = 1/100, a22 = 1/200, a12 = 1/100, N1_0 = 150, N2_0 = 200)  # priority effect (N2 wins)

#### phase diagram
phase_plane <- function(r1, r2, a11, a21, a22, a12, title, t){
    ### Vectors
  LV <- function(times, state, parms) {
    with(as.list(c(state, parms)), {
      dN1_dt = N1 * (r1 - a11*N1 - a12*N2)
      dN2_dt = N2 * (r2 - a22*N2 - a21*N1)
      return(list(c(dN1_dt, dN2_dt)))
    })
  }

  times <- c(0, t)
  parms <- c(r1 = r1, r2 = r2, a11 = a11, a21 = a21, a22 = a22, a12 = a12)

  x_inter<- max(c(r1/a11, r2/a21))
  y_inter <- max(c(r2/a22, r1/a12))

  # create position of arrows
  vector_grid <- expand.grid(seq(5, x_inter, length.out = 10),
                             seq(5, y_inter, length.out = 10))

  vector_data <- vector_grid %>%
    pmap(., function(Var1, Var2){
      state <- c(N1 = Var1, N2 = Var2)
      pop_size <- ode(func = LV, times = times, y = state, parms = parms)
      pop_size[2, 2:3]
    }) %>%
    bind_rows() %>%
    rename(xend = N1, yend = N2) %>%
    bind_cols(vector_grid) %>%
    rename(x = Var1, y = Var2)

    ### Phase plane
    ggplot() +
      geom_abline(slope = -a11/a12, intercept = r1/a12, color = "#E41A1C", size = 1.5) +
      geom_abline(slope = -a21/a22, intercept = r2/a22, color = "#377EB8", size = 1.5) +
      geom_segment(data = vector_data,
                   aes(x = x, y = y, xend = xend, yend = yend),
                   arrow = arrow(length = unit(0.1, "cm"))) +    
      scale_x_continuous(name = "N1", limits = c(0, x_inter), expand = c(0, 0)) +
      scale_y_continuous(name = "N2", limits = c(0, y_inter), expand = c(0, 0)) +
      theme_bw(base_size = 13) +
      theme(panel.grid = element_blank(),
            plot.title = element_text(hjust = 0.5)) +
      labs(title = title)
  }

  phase_plane(r1 = 1.2, r2 = 1.2, a11 = 1/100, a21 = 1/200, a22 = 1/100, a12 = 1/300, t = 0.3, title = "Stable coexistence")

  phase_plane(r1 = 1.2, r2 = 1.2, a11 = 1/200, a21 = 1/100, a22 = 1/200, a12 = 1/100, t = 0.3, title = "Unstable coexistence (saddle)")