Week 12

Resource competition: R* models

Lecture in a nutshell

  • Exploitative competition (consumer-resource dynamics)
    1. R* rule:
      • 2 consumers + 1 resource
        \(\begin{align}\frac {dN_{1}}{dt} = e_{1}a_{1}RN_{1}-dN_{1}\end{align}\\\)
        \(\begin{align}\frac {dN_{2}}{dt} = e_{2}a_{2}RN_{2}-dN_{2}\end{align}\\\)
        \(\begin{align}\frac {dR}{dt} = d(S_{0}-R)-a_{1}RN_{1}-a_{2}RN_{2}\end{align}\)
      • The two consumers CANNOT coexist: the one with lower R* (the resource level at the monoculture equilibrium) will win
      • Generalization: When multiple species are competing for a single resource, coexistence is not possible. The species with the lowest R* will win at equilibrium.
    2. Tilman’s resource ratio hypothesis:
      • 2 consumers + 2 resources
        \(\begin{align}\frac {dN_{1}}{dt} = e_{1}a_{1a}R_{a}N_{1}+e_{1}a_{1b}R_{b}N_{1}-dN_{1}\end{align}\\\)
        \(\begin{align}\frac {dN_{2}}{dt} = e_{2}a_{2a}R_{a}N_{2}+e_{2}a_{2b}R_{b}N_{2}-dN_{2}\end{align}\\\)
        \(\begin{align}\frac {dR_{a}}{dt} = d(S_{a}-R_{a})-a_{1a}R_{a}N_{1}-a_{2a}R_{a}N_{2}\end{align}\)
        \(\begin{align}\frac {dR_{b}}{dt} = d(S_{b}-R_{b})-a_{1b}R_{b}N_{1}-a_{2b}R_{b}N_{2}\end{align}\)
      • Outcomes:
        1. Consumer 1 wins: \(ZNGI_{2}\) lies above \(ZNGI_{1}\) in the \(R_{a}\)-\(R_{b}\) state space
        2. Consumer 2 wins: \(ZNGI_{1}\) lies above \(ZNGI_{2}\) in the \(R_{a}\)-\(R_{b}\) state space
        3. Coexist: \(ZNGI_{1}\) and \(ZNGI_{2}\) cross (trade-off in competitiveness); species should consume more of the resource that limits its growth (i.e., the resource of which the decline leads to the largest reduction in growth); resource supply should be balanced
        4. Priority effect
    3. Armstrong & McGehee relative nonlinearity
      • 2 consumers + 1 resource
        \(\begin{align}\frac {dN_{1}}{dt} = e_{1}a_{1}RN_{1}-d_{1}N_{1}\end{align}\\\)
        \(\begin{align}\frac {dN_{2}}{dt} = e_{2}\frac{a_{2}R}{k_{2}+R}N_{2}-d_{2}N_{2}\end{align}\\\)
        \(\begin{align}\frac {dR}{dt} = rR(1-\frac{R}{K})-a_{1}RN_{1}-\frac{a_{2}R}{k_{2}+R}N_{2}\end{align}\)
      • Logistically-growing resource + type-II functional response of one consumer
      • Fluctuating resource generates stable limit cycles (fluctuating-dependent coexistence): coexistence occurs if \(R^{**}\) < \(R^{*}\) < \(\bar{R}\)


Lab demonstration

Part 1 - R* rule

In this part of the lab section, we are going to analyze the consumer-resource dynamics with two consumers and one resource:

\(\begin{align}\frac {dN_{1}}{dt} = e_{1}a_{1}RN_{1}-dN_{1}\end{align}\\\)
\(\begin{align}\frac {dN_{2}}{dt} = e_{2}a_{2}RN_{2}-dN_{2}\end{align}\\\)
\(\begin{align}\frac {dR}{dt} = d(S_{0}-R)-a_{1}RN_{1}-a_{2}RN_{2}\end{align}\)

Please simulate the model using the parameter set (N1 = 2, N2 = 2, R = 0.1, a1 = 0.4, a2 = 0.6, e1 = 1, e2 = 1, d = 0.01, S0 = 0.1) and visualize the population trajectories of the two consumers as well as the resource.

Below is a plot of resource level vs. consumers’ per capita growth rate. What would you predict about the system dynamics?

library(tidyverse)
library(deSolve)

### Parameters
a1 <- 0.4 
a2 <- 0.6
e1 <- 1
e2 <- 1
d <- 0.01 
S0 <- 0.1

### Resource level vs. consumers' per capita population growth
data.frame(R = seq(0, 0.1, 0.001)) %>%
  mutate(N1 = e1*a1*R-d,
         N2 = e2*a2*R-d) %>%
  pivot_longer(-R, names_to = "species", values_to = "N") %>%
  ggplot(aes(x = R, y = N, color = species)) +
  geom_line(size = 1.2) + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x = "Resource level", y = "Per capita growth rate") +
  scale_x_continuous(limits = c(0, 0.11), expand = c(0, 0)) +
  scale_y_continuous(limits = c(NA, 0.05), expand = c(0, 0)) +
  scale_color_brewer(name = NULL, palette = "Set1", direction = -1) + 
  theme_classic(base_size = 14)

### Model specification
CR_model_2C_1R <- function(times, state, parms){
  with(as.list(c(state, parms)), {
    dN1_dt = e1*a1*R*N1 - d*N1
    dN2_dt = e2*a2*R*N2 - d*N2
    dR_dt =  d*(S0-R) - a1*R*N1 - a2*R*N2
    return(list(c(dN1_dt, dN2_dt, dR_dt)))
  })
}

### Model parameters
times <- seq(0.1, 1000, by = 0.1)  
state <- c(N1 = 2, N2 = 2, R = 0.1)  
parms <- c(a1 = 0.4, a2 = 0.6, e1 = 1, e2 = 1, d = 0.01, S0 = 0.1)  # R and S0 should be the same in the chemostat

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

### Visualize the population dynamics
pop_size %>%
  as.data.frame() %>%
  pivot_longer(cols = -time, names_to = "species", values_to = "N") %>%
  mutate(trophic = case_when(species %in% c("N1", "N2") ~ "Consumer",
                             TRUE ~ "Resource")) %>%
  ggplot(aes(x = time, y = N, color = species)) + 
  geom_line(size = 1.5) +
  facet_wrap(~ trophic, 
             ncol = 2, 
             scales = "free_y",
             strip.position = "left") +
  theme_classic(base_size = 14) +
  theme(strip.background = element_blank(),
        strip.placement = "outside",
        legend.position = "top",
        legend.title = element_blank(),
        plot.margin = margin(r = 5)) + 
  labs(x = "Time", y = NULL) +
  scale_x_continuous(limits = c(0, 1050), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, NA), expand = c(0, 0)) +
  scale_color_manual(name = NULL, values = c("#377EB8", "#E41A1C", "#4DAF4A"))

We can see that N1 goes extinct, as predicted by the R* rule; the remaining consumer N2 and the resource R exhibit stable coexistence.


Part 2 - Tilman’s resource ratio hypothesis

Now we will extend the model in Part 1. into two consumers and two resources:

\(\begin{align}\frac {dN_{1}}{dt} = e_{1}a_{1a}R_{a}N_{1}+e_{1}a_{1b}R_{b}N_{1}-dN_{1}\end{align}\\\)
\(\begin{align}\frac {dN_{2}}{dt} = e_{2}a_{2a}R_{a}N_{2}+e_{2}a_{2b}R_{b}N_{2}-dN_{2}\end{align}\\\)
\(\begin{align}\frac {dR_{a}}{dt} = d(S_{a}-R_{a})-a_{1a}R_{a}N_{1}-a_{2a}R_{a}N_{2}\end{align}\)
\(\begin{align}\frac {dR_{b}}{dt} = d(S_{b}-R_{b})-a_{1b}R_{b}N_{1}-a_{2b}R_{b}N_{2}\end{align}\)

Here, the two consumers N1 and N2 compete for two substitutable resources R1 and R2. Under certain conditions, the two consumers are able to coexist:

library(tidyverse)
library(deSolve)

### Model specification
CR_model_2C_2R <- function(time, state, parms){
  with(as.list(c(state, parms)), {
    dN1_dt = e1*a1a*Ra*N1 + e1*a1b*Rb*N1 - d*N1
    dN2_dt = e2*a2a*Ra*N2 + e2*a2b*Rb*N2 - d*N2
    dRa_dt = d*(Sa-Ra) - (a1a*N1*Ra) - (a2a*N2*Ra) 
    dRb_dt = d*(Sb-Rb) - (a1b*N1*Rb) - (a2b*N2*Rb)
    return(list(c(dN1_dt, dN2_dt, dRa_dt, dRb_dt)))
  })
}

### Model parameters
times <- seq(0.1, 2000, by = 0.01)
state <- c(N1 = 0.05, N2 = 0.05, Ra = 0.3, Rb = 0.3)
parms <- c(a1a = 0.4, a1b = 0.8, a2a = 0.6, a2b = 0.5, e1 = 1, e2 = 1, d = 0.1, Sa = 0.3, Sb = 0.3)  # Ra/Rb and S0a/S0b should be the same in the chemostat

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

### Visualize the population dynamics
pop_size %>%
  as.data.frame() %>%
  pivot_longer(cols = -time, names_to = "species", values_to = "N") %>%
  mutate(trophic = case_when(species %in% c("N1", "N2") ~ "Consumer",
                             TRUE ~ "Resource")) %>%
  ggplot(aes(x = time, y = N, color = species)) + 
  geom_line(size = 1.5) +
  facet_wrap(~ trophic, 
             ncol = 2, 
             scales = "free_y",
             strip.position = "left") +
  theme_classic(base_size = 14) +
  theme(strip.background = element_blank(),
        strip.placement = "outside",
        legend.position = "top",
        legend.title = element_blank(),
        plot.margin = margin(r = 8)) + 
  labs(x = "Time", y = NULL) +
  scale_x_continuous(limits = c(0, 2050), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, NA), expand = c(0, 0)) +
  scale_color_manual(name = NULL, values = c("#377EB8", "#E41A1C", "#4DAF4A", "#984EA3"))


The phase diagram also illustrates stable coexistence of the two consumers:

### Parameters
a1a <- 0.4
a1b <- 0.8
a2a <- 0.6
a2b <- 0.5
e1 <- 1
e2 <- 1
d <- 0.1

### Slopes and intercepts of the ZNGI's
ZNGI_slope_N1 <- -a1a/a1b
ZNGI_intercept_N1 <- d/(e1*a1b)  
ZNGI_slope_N2 <- -a2a/a2b
ZNGI_intercept_N2 <- d/(e2*a2b)

### Consumption vectors
eqilibrium_Ra <- (d/e1)*((a1b-a2b)/(a2a*a1b-a2b*a1a))
eqilibrium_Rb <- (d/e2)*((a1a-a2a)/(a2b*a1a-a2a*a1b))

convec_df <- data.frame(x = c(eqilibrium_Ra + 6*a1a*eqilibrium_Ra, 
                              eqilibrium_Ra + 6*a2a*eqilibrium_Ra),
                        y = c(eqilibrium_Rb + 6*a1b*eqilibrium_Rb, 
                              eqilibrium_Rb + 6*a2b*eqilibrium_Rb),
                        xend = c(eqilibrium_Ra - a1a*eqilibrium_Ra, 
                              eqilibrium_Ra - a2a*eqilibrium_Ra),
                        yend = c(eqilibrium_Rb - a1b*eqilibrium_Rb, 
                              eqilibrium_Rb - a2b*eqilibrium_Rb),
                        species = c("N1", "N2"))

### Phase diagram
ggplot() + 
  geom_abline(slope = ZNGI_slope_N1, intercept = ZNGI_intercept_N1, color = "#377EB8", size = 1.2) + 
  geom_abline(slope = ZNGI_slope_N2, intercept = ZNGI_intercept_N2, color = "#E41A1C", size = 1.2) + 
  geom_segment(data = convec_df, aes(x = x, y = y, xend = xend, yend = yend, color = species), linetype = "blank") + 
  geom_segment(data = convec_df, aes(x = x, y = y, xend = xend, yend = yend, color = species), size = 0.5, linetype = "dashed", arrow = arrow(type = "closed", length = unit(0.1, "inches")), show.legend = F) +
  geom_path(data = as.data.frame(pop_size), aes(x = Ra, y = Rb)) +
  geom_point(data = as.data.frame(pop_size), aes(x = last(Ra), y = last(Rb)), size = 2.5) +
  theme_classic(base_size = 14) +
  labs(x = expression(italic(R[a])), y = expression(italic(R[b]))) +
  scale_x_continuous(limits = c(0, 0.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 0.5), expand = c(0, 0)) +
  scale_color_brewer(name = NULL, palette = "Set1", direction = -1,
                     guide = guide_legend(override.aes = list(
                       linetype = "solid", size = 1.2))) +
  coord_fixed(ratio = 1)

In this particular case, if the initial resource levels Ra0 and Rb0 fall within the upper area bounded by the two consumption vectors, then the two species will coexist.


Part 3 - Relative nonlinearity

We have seen in Part 1. that two species cannot coexist on one resource if the consumers consume it in a linear fashion. What if the consumers consume the resource in a non-linear fashion? Could that generate stable coexistence?

In this part of the lab section, we will introduce non-linear (type II) functional response into the consumer-resource model in Part 1. for one of the consumers. We also change the resource equation to logistic growth and allow the two consumers to have different conversion efficiencies and mortality rates:

\(\begin{align}\frac {dN_{1}}{dt} = e_{1}a_{1}RN_{1}-d_{1}N_{1}\end{align}\\\)
\(\begin{align}\frac {dN_{2}}{dt} = e_{2}\frac{a_{2}R}{k_{2}+R}N_{2}-d_{2}N_{2}\end{align}\\\)
\(\begin{align}\frac {dR}{dt} = rR(1-\frac{R}{K})-a_{1}RN_{1}-\frac{a_{2}R}{k_{2}+R}N_{2}\end{align}\)

library(tidyverse)
library(deSolve)

### Model specification
CR_model_nonlinear <- function(times, state, parms){
  with(as.list(c(state, parms)), {
         dN1_dt = e1*a1*R*N1 - d1*N1
         dN2_dt = e2*(a2*R/(k2+R))*N2 - d2*N2
         dR_dt =  r*R*(1-(R/K)) - a1*R*N1 - ((a2*R)/(k2+R))*N2
         return(list(c(dN1_dt, dN2_dt, dR_dt)))
       })
}

### Model parameters
times <- seq(0, 5000, by = 0.1)  
state <- c(N1 = 0.3, N2 = 19, R = 106) 
parms <- c(a1 = 0.003, a2 = 0.5, k2 = 50, e1 = 0.33, e2 = 0.3, d1 = 0.11, d2 = 0.1, r = 0.1, K = 300)

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

pop_size %>%
  as.data.frame() %>%
  pivot_longer(cols = -time, names_to = "species", values_to = "N") %>%
  mutate(trophic = case_when(species %in% c("N1", "N2") ~ "Consumer",
                             TRUE ~ "Resource")) %>%
  ggplot(aes(x = time, y = N, color = species)) + 
  geom_line(size = 1.5) +
  facet_wrap(~ trophic, 
             ncol = 2, 
             scales = "free_y",
             strip.position = "left") +
  theme_classic(base_size = 14) +
  theme(strip.background = element_blank(),
        strip.placement = "outside",
        legend.position = "top",
        legend.title = element_blank(),
        plot.margin = margin(r = 5)) + 
  labs(x = "Time", y = NULL) +
  scale_x_continuous(limits = c(0, 5050), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, NA), expand = c(0, 0)) +
  scale_color_manual(name = NULL, values = c("#377EB8", "#E41A1C", "#4DAF4A"))

As shown in the figure, the two consumers coexist in a fluctuating manner. In fact, the fluctuating properties generated by the combination of logistically growing resources and different nonlinear functions are critical for the coexistence of n species on k < n resources in certain parameter spaces.