Competition

Lotka-volterra Interspecific Continuous Logistic Growth Model

Today we will focus on “Interspecific Continuous Logistic Growth” which is common in nature. Let’s break this term into individual parts:

Inter: Between
Specific: Species
Continuous: Not discrete
Logistic: Density dependent

From lab 3, do you remember what the equation for continuous logistic growth?

\[ rN(1-\frac{N}{K}) \]

Why do we use \(r\) in this equations instead of \(\lambda\)?

Because we are modeling continuous population growth.

Why did we incorporate \(K\)?

Because we are modeling logistic population growth.

We can rearrange to help us convert this to relate to interspecific competition:

\[ rN(\frac{K}{K} - \frac{N}{K}) \]

\[ rN(\frac{K-N}{K}) \]

This continuous logistic growth equation is currently for a single species, but our model today should be interspecific. In order to do so, we can add the term \(\alpha\) which is the effect of species 2 on the population growth of species 1 (i.e. how many individuals of species 2 are equivalent to one individual of species 1 in regards to resource usage).

\[ \frac{dN_{1}}{dt} = r_{1}N_{1}(\frac{K_{1} - N_{1} - \alpha N_{2}}{K_{1}}) \]

We can simplify this equation by integrating carrying capacity \(K\) into the term \(\alpha_{11}\) representing conspecific/intraspecific density-dependence.

We can also combine terms into interspecific density dependence or “interspecific competition coefficient” \(\alpha _{12}\).

The result is the “Lotka-Volterra Interspecific Density-dependent Growth” model:

\[ \frac{dN_{1}}{dt} \ =\ r_{1} N_{1} \ ( \ 1\ -\ \alpha _{11} N_{1} \ -\ \alpha _{12} N_{2}) \]

\[ \frac{dN_{2}}{dt} \ =\ r_{2} N_{2} \ ( \ 1\ -\ \alpha _{22} N_{1} \ -\ \alpha _{21} N_{1}) \]

Now let’s create a function for this model in R. Once again, we can write differential equations in R using the ordinary differential equation (ODE) package deSolve(Soetaert, Petzoldt, and Setzer 2010) following the format:

Write a function to represent interspecific density-dependent growth between 2 species

interspec_dens_dep_growth <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1 * n[1] * (1 - a11 * n[1] - a12 * n[2])
    dn2dt <- r2 * n[2] * (1 - a22 * n[2] - a21 * n[1])
    list(c(dn1dt, dn2dt))
  })
}

Run your function with the following parameters:

Intrinsic rate of growth of species A is 1 and of species B is 0.1.
The interspecific negative density dependence of species A on B is 0.1. Vice versa is 0.01.
The intraspecific density dependence is 0.2 and 0.02 of species A and B respectively.
Initial population sizes of species A and B are 2 and 1 respectively
Time range will be 1:100

library(deSolve)

Warning: package 'deSolve' was built under R version 4.3.2

r1 <- 1
r2 <- 0.1
a11 <- 0.2
a21 <- 0.1
a22 <- 0.02
a12 <- 0.01

times <- c(1:100)

N0_1 <- 2
N0_2 <- 1

parms <- c(r1 = r1, r2 = r2, a11 = a11, a21 = a21, a22 = a22, a12 = a12)
initialN <- c(N0_1, N0_2)

out <- deSolve::ode(y = initialN, times = times, func = interspec_dens_dep_growth, parms = parms)

Now plot your results:

library(dplyr)

Warning: package 'dplyr' was built under R version 4.3.2


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.3.2

ode_df <- out %>% 
  as.data.frame() %>% # Convert to dataframe for ggplot2
  rename(time = 1, V1 = 2, V2 = 3)  # Rename because numbered column names don't work well with ggplot2


ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = V1, colour = "Species 1")) + 
  geom_line(aes(y = V2, colour = "Species 2")) + 
  labs(x = "Years", y = "Population Size", colour = "Species")

Describe what you observe in the output of this model with the given parameters. How are the populations of species A and B growing/declining? How are they influencing one another?

# A: Species 1 is strongly influenced by the abundance of species 2. At first species 1 increases following a logistic growth model, however once species 2 population is high enough it competes with resources with species 1 and species 1 population drops. If this continued, species 1 woud eventually go extinct

# Meanwhile, the population of species 2 continues to grow exponentially since it has not met its carrying capacity yet.

Instability

By modifying the parameters a12, a21, a11, and a22, can you find values where species A goes extinct? Can you find values where species B goes extinct?

r1 <- 1
r2 <- 0.2
a11 <- 0.2
a21 <- 0.1
a22 <- 0.02
a12 <- 0.03


parms <- c(r1 = r1, r2 = r2, a11 = a11, a21 = a21, a22 = a22, a12 = a12)
initialN <- c(2, 1)

out <- ode(y = initialN, times = 1:100, func = interspec_dens_dep_growth, parms = parms)

ode_df <- out %>% 
  as.data.frame() %>% # Convert to dataframe for ggplot2
  rename(time = 1, V1 = 2, V2 = 3)  # Rename because numbered column names don't work well with ggplot2


ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = V1, colour = "Species 1")) + 
  geom_line(aes(y = V2, colour = "Species 2")) + 
  labs(x = "Years", y = "Population Size", colour = "Species")

# A: We can use this function to find various combinations of values where one species goes extinct and the other reaches a carrying capacity (shown above). This further demonstrates how the populations of these two species influence one another. It is also possible to find combinations where neither species goes extinct.

Neutral Equilibria

Set your parameters to the following conditions:

r1 <- r2 <- 1
a11 <- a21 <- 0.01
a22 <- a12 <- 0.015
initialN <- c(2, 1)

What do you observe? Is this result stable or unstable? Or something more nuanced? Explain.

r1 <- r2 <- 1
a11 <- a21 <- 0.01
a22 <- a12 <- 0.015
initialN <- c(2, 1)

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

out <- deSolve::ode(y = initialN, times = 1:100, func = interspec_dens_dep_growth, parms = parms)

ode_df <- out %>% 
  as.data.frame() %>% # Convert to dataframe for ggplot2
  rename(time = 1, V1 = 2, V2 = 3)  # Rename because numbered column names don't work well with ggplot2


ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = V1, colour = "Species 1")) + 
  geom_line(aes(y = V2, colour = "Species 2")) + 
  labs(x = "Years", y = "Population Size", colour = "Species")

# A: We observe something called a "neutral saddle." The relationship is not considered stable because after perturbation although the ratio between the population sizes of the two species will remain constant, a new equilibrium point will be established for each species' population size.

Interspecific Continuous Logistic Growth With Real Data

From year 2000 to 2014, (Reif et al. 2019)(Reif et al. 2018) collected nightingale (n=552) distribution data from 959 plots, focusing on sibling species Common Nightingale (Luscinia megarhynchos) and the Thrush Nightingale (L. luscinia). This study was conducted in Poland where these species overlap in range:

Along with population counts, researchers noted the following relationships:

Allopatric: Two or more related species which have separate distributions geographically.
Sympatric: Two or more related species which have the same or overlapping distributions geographically.
- Syntopic: Two or more related species which have the same or overlapping distributions geographically, and share the same habitats within that distribution
- Allotopic: Two or more related species which have the same or overlapping distributions geographically, but do not share the same habitats within that distribution

Let’s load in the nightingales.csv data and assign to the variable nightingales_raw. (Hint: if you are having trouble loading the data, view the file in Rstudio and not the separator. You may need to define the sep paramater in your data loading function) Briefly explore the dataset:

nightingales_raw <- read.csv("nightingales.csv", sep = ";")

head(nightingales_raw)

  plot count  co.occ year  env.grad obs      lon      lat
1    1     1 CN.synt 2010 -132.3607 100 14.97798 52.35997
2    1     3 TN.synt 2010 -132.3607 100 14.97798 52.35997
3    2     0 TN.alot 2009  -87.9261 162 15.80042 52.38422
4    2     4 CN.alot 2009  -87.9261 162 15.80042 52.38422
5    3     1 CN.alot 2010 -125.6040 110 14.56317 52.74558
6    3     0 TN.alot 2010 -125.6040 110 14.56317 52.74558

summarize(nightingales_raw)

data frame with 0 columns and 1 row

Clean the data using the separate() function from the (tidyr)(Wickham and Wickham 2017) package to separate species and co_occurence_status.

library(tidyr)

Warning: package 'tidyr' was built under R version 4.3.2

# Cleaning data

nightingales_disag <- nightingales_raw %>%
  tidyr::separate(col = co.occ, into = c("species", "co_occurrence_status"), sep = "\\.")


head(nightingales_disag)

  plot count species co_occurrence_status year  env.grad obs      lon      lat
1    1     1      CN                 synt 2010 -132.3607 100 14.97798 52.35997
2    1     3      TN                 synt 2010 -132.3607 100 14.97798 52.35997
3    2     0      TN                 alot 2009  -87.9261 162 15.80042 52.38422
4    2     4      CN                 alot 2009  -87.9261 162 15.80042 52.38422
5    3     1      CN                 alot 2010 -125.6040 110 14.56317 52.74558
6    3     0      TN                 alot 2010 -125.6040 110 14.56317 52.74558

To help us visualize our dataset, let’s map our data and facet by co_occurence_status

# Map

ggplot(data = nightingales_disag, aes(x = lon, y = lat, color = species)) +
  geom_point() +
  facet_wrap(~co_occurrence_status)

Let’s aggregate our data by year:

nightingales_agg <- nightingales_raw %>%
  separate(col = co.occ, into = c("species", "co_occurrence_status"), sep = "\\.") %>%
  group_by(species, year, co_occurrence_status) %>%
  summarize(avg_count = mean(count)) %>%
  ungroup() %>%
  arrange(species, year)

`summarise()` has grouped output by 'species', 'year'. You can override using
the `.groups` argument.

# We are using summarize() instead of mutate because summarize() will make the group_by() permanent and make the `avg_count` permanent. If we used mutate() in this case, the column created would disappear once summarize() was called.

nightingales_agg

# A tibble: 90 × 4
   species  year co_occurrence_status avg_count
   <chr>   <int> <chr>                    <dbl>
 1 CN       2000 alop                     1    
 2 CN       2000 alot                     0.333
 3 CN       2000 synt                     0.688
 4 CN       2001 alop                     2.13 
 5 CN       2001 alot                     0.541
 6 CN       2001 synt                     0.630
 7 CN       2002 alop                     1.45 
 8 CN       2002 alot                     0.452
 9 CN       2002 synt                     0.686
10 CN       2003 alop                     1.12 
# ℹ 80 more rows

Inspect the dimensions, do they make sense?

(2 species) x (3 conditions) x (15 years) = 90 rows

group_by() and ungroup()

(dplyr)(Wickham et al. 2022)’s group_by() function is a powerful and convenient way to separate data into groups and perform operations on each group. However, it is important to remember to always use ungroup() after you are done with your group operations. This is because group_by() will continue to affect your data until you use ungroup(). This is even more dangerous because visually you will not see any difference in your data after using group_by().

Our goal is to find the mean of the count column for all rows. Let’s group_by() the column co.occ:

nightingales_grouped<- nightingales_raw %>%
  group_by(co.occ)

head(nightingales_grouped)

# A tibble: 6 × 8
# Groups:   co.occ [4]
   plot count co.occ   year env.grad   obs   lon   lat
  <int> <int> <chr>   <int>    <dbl> <int> <dbl> <dbl>
1     1     1 CN.synt  2010   -132.    100  15.0  52.4
2     1     3 TN.synt  2010   -132.    100  15.0  52.4
3     2     0 TN.alot  2009    -87.9   162  15.8  52.4
4     2     4 CN.alot  2009    -87.9   162  15.8  52.4
5     3     1 CN.alot  2010   -126.    110  14.6  52.7
6     3     0 TN.alot  2010   -126.    110  14.6  52.7

How does the dataset look different after grouping?

Besides the small “Groups: co.occ [4]” at the top of the output, the data does not appear any different after grouping

Now let’s try to meet our goal and find the mean of all count rows.

summarize(nightingales_grouped, mean(count))

# A tibble: 6 × 2
  co.occ  `mean(count)`
  <chr>           <dbl>
1 CN.alop         1.65 
2 CN.alot         0.491
3 CN.synt         0.723
4 TN.alop         1.25 
5 TN.alot         0.453
6 TN.synt         0.999

We have 6 mean values instead of the desired overall mean!

This is because we did not ungroup() like we always should. So even thought there is no visual difference after using group_by(), not using ungroup() affects your code downstream.

Now let’s try the same but use ungroup() appropriately:

nightingales_ungrouped <- nightingales_raw %>%
  group_by(co.occ) %>%
  ungroup()

head(nightingales_ungrouped)

# A tibble: 6 × 8
   plot count co.occ   year env.grad   obs   lon   lat
  <int> <int> <chr>   <int>    <dbl> <int> <dbl> <dbl>
1     1     1 CN.synt  2010   -132.    100  15.0  52.4
2     1     3 TN.synt  2010   -132.    100  15.0  52.4
3     2     0 TN.alot  2009    -87.9   162  15.8  52.4
4     2     4 CN.alot  2009    -87.9   162  15.8  52.4
5     3     1 CN.alot  2010   -126.    110  14.6  52.7
6     3     0 TN.alot  2010   -126.    110  14.6  52.7

summarize(nightingales_ungrouped, mean(count))

# A tibble: 1 × 1
  `mean(count)`
          <dbl>
1         0.940

Now plot avg_count across time for nightingales_agg using geom_point() and geom_smooth, faceted by co_occurrence_status

ggplot(data = nightingales_agg, aes(x = year, y = avg_count, color = species)) +
  geom_point() +
  geom_smooth() +
  facet_wrap(~co_occurrence_status)

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

It looks like there might be a difference between co-occurrence statuses. Let’s see if competition is occurring. In order to do so we will need to extract the following parameters from our data:

r1 = Intrinsic rate of growth of species CN
r2 = Intrinsic rate of growth of species TN
a21 = The interspecific negative density dependence of species CN on TN
a12 = The interspecific negative density dependence of species TN on CN
a11 = The intraspecific density dependence of species CN
a22 = The intraspecific density dependence of species TN
N0 = Initial population sizes of each species
times = Time range of dataset in years

Time Range

times <- c(1:(max(nightingales_agg$year) - min(nightingales_agg$year) + 1))

times

 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

Initial Population Sizes

N0_CN <- nightingales_agg %>%
  filter(species == "CN") %>%
  pull(avg_count) %>%
  first()

N0_CN

[1] 1

N0_TN <- nightingales_agg %>%
  filter(species == "TN") %>%
  pull(avg_count) %>%
  first()

N0_TN

[1] 1.236842

That will be sufficient for most datasets, however in this case we will now repeat for the other conditions:

N0_CN_alop <- nightingales_agg %>%
  filter(species == "CN", co_occurrence_status == "alop") %>%
  pull(avg_count) %>%
  first()

N0_CN_alot <- nightingales_agg %>%
  filter(species == "CN", co_occurrence_status == "alot") %>%
  pull(avg_count) %>%
  first()

N0_CN_synt <- nightingales_agg %>%
  filter(species == "CN", co_occurrence_status == "synt") %>%
  pull(avg_count) %>%
  first()

N0_TN_alop <- nightingales_agg %>%
  filter(species == "TN", co_occurrence_status == "alop") %>%
  pull(avg_count) %>%
  first()

N0_TN_alot <- nightingales_agg %>%
  filter(species == "TN", co_occurrence_status == "alot") %>%
  pull(avg_count) %>%
  first()

N0_TN_synt <- nightingales_agg %>%
  filter(species == "TN", co_occurrence_status == "synt") %>%
  pull(avg_count) %>%
  first()

Intrinsic Population Growth Rates

A simplified estimation of intrinsic population growth is to calculate the average value of intrinsic growth rate within a given interval:

\[ r = frac{log(\frac{\bar N_{t_{final}}}{\bar N_{0}})}{\bar N_{t_{final}}-{\bar N_{0}}} \]

r1 <- nightingales_agg %>%
  filter(co_occurrence_status == "alop", species == "CN") %>%
  summarise(r = log(last(avg_count) / first(avg_count)) / (last(year) - first(year))) %>%
  pull(r)


r2 <- nightingales_agg %>%
  filter(co_occurrence_status == "alop", species == "TN") %>%
  summarise(r = log(last(avg_count) / first(avg_count)) / (last(year) - first(year))) %>%
  pull(r)

r1

[1] 0.05923424

r2

[1] -0.008801167

Let’s plot to check to see if this is a decent approximation:

# Calculate predicted count using our intrinsic growth rates
nightingales_pred <- nightingales_agg %>%
  mutate(predicted_count = case_when(
    species == "CN" ~ N0_CN_alop * exp(r1 * (year - min(year))),
    species == "TN" ~ N0_TN_alop * exp(r2 * (year - min(year)))
  )) %>%
  filter(co_occurrence_status == "alop")

nightingales_agg %>%
  filter(co_occurrence_status == "alop") %>%
  ggplot(aes(x = year, y = avg_count, color = species)) + 
    geom_point() +
    geom_line(data = nightingales_pred, aes(x = year, y = predicted_count), linetype = "dashed") +
    geom_smooth() +
    labs(x = "Year", y = "Average Count", color = "Species", title = "Actual vs. Predicted Population Counts")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

It seems like our simplified estimate will be sufficient for this analysis, however note that the curves are not an excellent match.

Here are some examples of how a more accurate model of intrinsic growth rate could be made, by fixing to a non-linear model and using the continuous exponential population growth function we covered earlier this semester:

\[ Ne^{rt} \]

alop_data <- nightingales_disag %>%
  filter(co_occurrence_status == "alop") %>%
  mutate(min_year = min(year)) %>%
  mutate(year_since_min = year - min_year) %>%
  select(-min_year) %>%
  group_by(species) %>%
  mutate(N0 = if_else(year == min(year), count, NA_real_)) %>%
  mutate(N0 = if_else(is.na(N0), ave(N0, FUN = function(x) mean(x, na.rm = TRUE)), N0)) %>%
  ungroup()


# Setup the exponential growth model function for nls fitting
fit_growth_rate <- function(data, N0_col, time_col = "year_since_min", count_col = "count") {
  formula <- as.formula(sprintf("%s ~ %s * exp(r * %s)", count_col, N0_col, time_col))
  start_list <- list(r = 0.1) # Starting value for 'r', this is a "best guess"
  model <- nls(formula, data = data, start = start_list, algorithm = "port")
  return(model)
}

# Fit the model for each species
models <- lapply(split(alop_data, alop_data$species), function(species_data) {
  fit_growth_rate(species_data, N0_col = "N0")
})

# Print the models to see the fitted parameters
print(models)

$CN
Nonlinear regression model
  model: count ~ N0 * exp(r * year_since_min)
   data: data
     r 
0.0562 
 residual sum-of-squares: 2406

Algorithm "port", convergence message: relative convergence (4)

$TN
Nonlinear regression model
  model: count ~ N0 * exp(r * year_since_min)
   data: data
        r 
5.365e-05 
 residual sum-of-squares: 5673

Algorithm "port", convergence message: relative convergence (4)

alop_data <- nightingales_agg %>%
  filter(co_occurrence_status == "alop") %>%
  group_by(species) %>%
  mutate(N0 = avg_count[which.min(year)]) %>% 
  ungroup() %>%
  mutate(year_since_min = year - min(year))

# Define the exponential growth model function
exponential_growth_model <- function(t, r, N0) {
  N0 * exp(r * t)
}

# Setup the nls fitting function
fit_growth_rate <- function(data, time_col = "year_since_min", count_col = "avg_count") {
  formula <- as.formula(paste(count_col, "~", "N0", "* exp(r *", time_col, ")"))  # Using "N0" directly instead of N0_col
  start_list <- list(r = 0.1)  # Starting value for 'r'
  model <- nls(formula, data = data, start = start_list, algorithm = "port")
  return(summary(model))
}

# Fit the model for each species
models <- lapply(split(alop_data, alop_data$species), function(species_data) {
  fit_growth_rate(species_data)
})

# Print the summary of the models to see the fitted parameters
lapply(models, print)


Formula: avg_count ~ N0 * exp(r * year_since_min)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
r 0.055988   0.006361   8.802 4.43e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3807 on 14 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)


Formula: avg_count ~ N0 * exp(r * year_since_min)

Parameters:
  Estimate Std. Error t value Pr(>|t|)
r 0.001142   0.004100   0.278    0.785

Residual standard error: 0.1636 on 14 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)

$CN

Formula: avg_count ~ N0 * exp(r * year_since_min)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
r 0.055988   0.006361   8.802 4.43e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3807 on 14 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)


$TN

Formula: avg_count ~ N0 * exp(r * year_since_min)

Parameters:
  Estimate Std. Error t value Pr(>|t|)
r 0.001142   0.004100   0.278    0.785

Residual standard error: 0.1636 on 14 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)

Intraspecific Density Dependencies

For intraspecific density dependencies, we can perform a linear regression. Again, a non-linear regression would likely be better for this dataset, however this will provide a simplified estimate:

alop_data <- nightingales_agg %>%
  filter(co_occurrence_status == "alop") %>% # Only selecting "alop" so that interspecies dynamics at play. If your dataset is not separated into alopatric/sympatric samples, you can skip this step
  arrange(species, year) %>%
  group_by(species) %>%
  mutate(lag_avg_count = lag(avg_count)) %>% # Add a lagged population size column
  filter(!is.na(lag_avg_count)) %>% # Remove the first year since it won't have a lagged population size
  mutate(change_in_count = avg_count - lag_avg_count) # Calculate the change in population size
  

# Perform a linear regression to estimate the intraspecific density dependence
# For CN (a11)
a11_model <- lm(change_in_count ~ lag_avg_count, data = filter(alop_data, species == "CN"))
a11 <- unname(coef(a11_model)["lag_avg_count"])

# For TN (a22)
a22_model <- lm(change_in_count ~ lag_avg_count, data = filter(alop_data, species == "TN"))
a22 <- unname(coef(a22_model)["lag_avg_count"])

# Output the results
a11

[1] -0.4923871

a22

[1] -1.054365

Interspecific Negative Density Dependencies

As a simple way to estimate interspecific negative density dependencies if you do not have allopatric vs. sympatric data, you can still calculate a12 and a21 like we did in lecture, by taking the mean population size of each species and dividing them by one another:

\[ a_{12} = \frac{\bar N_{2}}{\bar N_{1}} \]

CN_avg <- nightingales_agg %>%
  filter(species == "CN", co_occurrence_status == "synt") %>% # Only selecting "synt" so that interspecies dynamics are strongest. If your dataset is not separated into alopatric/sympatric samples, you can skip this step
  summarise(avg = mean(avg_count)) %>%
  pull(avg)

TN_avg <- nightingales_agg %>%
  filter(species == "TN", co_occurrence_status == "synt") %>% # Only selecting "synt" so that interspecies dynamics are strongest. If your dataset is not separated into alopatric/sympatric samples, you can skip this step
  summarise(avg = mean(avg_count)) %>%
  pull(avg)

a12 <- TN_avg / CN_avg
a21 <- CN_avg / TN_avg

a12

[1] 1.522671

a21

[1] 0.6567406

However with this dataset we can take a more complex approach. We can calculate the mean population sizes under ‘alop’ and ‘synt’ conditions and then calculate the differences. We can then use these differences to estimate a12 and a21:

# Calculate the mean population sizes under 'alop' and 'synt' conditions
alop_data <- nightingales_agg %>%
  filter(co_occurrence_status == "alop") %>%
  group_by(species, year) %>%
  summarise(avg_count_alop = mean(avg_count), .groups = 'drop')

synt_data <- nightingales_agg %>%
  filter(co_occurrence_status == "synt") %>%
  group_by(species, year) %>%
  summarise(avg_count_synt = mean(avg_count), .groups = 'drop')

# Join the datasets and calculate the differences
comparison <- left_join(alop_data, synt_data, by = c("species", "year")) %>%
  mutate(difference = avg_count_synt - avg_count_alop)

# Estimate a12 and a21 by looking at the differences
# a21: effect of CN on TN
a21 <- comparison %>%
  filter(species == "TN") %>%
  summarise(a21 = mean(difference, na.rm = TRUE)) %>%
  pull(a21)

# a12: effect of TN on CN
a12 <- comparison %>%
  filter(species == "CN") %>%
  summarise(a12 = mean(difference, na.rm = TRUE)) %>%
  pull(a12)

# Output the results
a21

[1] -0.1859028

a12

[1] -0.875976

Now let’s do the same for each co-occurence status separately:

# Calculate mean population sizes for CN and TN under 'alop' condition
comparison_alop <- nightingales_agg %>%
  filter(co_occurrence_status == "alop") %>%
  group_by(year, species) %>%
  summarise(avg_count = mean(avg_count)) %>%
  pivot_wider(names_from = species, values_from = avg_count) %>%
  mutate(difference = TN - CN) %>%
  ungroup() %>%
  summarise(a21_alop = mean(difference),
            a12_alop = -mean(difference))  # Negative because it's the effect of TN on CN

`summarise()` has grouped output by 'year'. You can override using the
`.groups` argument.

a12_alop <- comparison_alop$a12_alop
a21_alop <- comparison_alop$a21_alop

a12_alop

[1] 0.3237545

a21_alop

[1] -0.3237545

# Calculate mean population sizes for CN and TN under 'alot' condition
comparison_alot <- nightingales_agg %>%
  filter(co_occurrence_status == "alot") %>%
  group_by(year, species) %>%
  summarise(avg_count = mean(avg_count)) %>%
  pivot_wider(names_from = species, values_from = avg_count) %>%
  mutate(difference = TN - CN) %>%
  ungroup() %>%
  summarise(a21_alot = mean(difference),
            a12_alot = -mean(difference))  # Negative because it's the effect of TN on CN

`summarise()` has grouped output by 'year'. You can override using the
`.groups` argument.

a12_alot <- comparison_alot$a12_alot
a21_alot <- comparison_alot$a21_alot

a12_alot

[1] -0.01047506

a21_alot

[1] 0.01047506

# Calculate mean population sizes for CN and TN under 'synt' condition
comparison_synt <- nightingales_agg %>%
  filter(co_occurrence_status == "synt") %>%
  group_by(year, species) %>%
  summarise(avg_count = mean(avg_count)) %>%
  pivot_wider(names_from = species, values_from = avg_count) %>%
  mutate(difference = TN - CN) %>%
  ungroup() %>%
  summarise(a21_synt = mean(difference),
            a12_synt = -mean(difference))  # Negative because it's the effect of TN on CN

`summarise()` has grouped output by 'year'. You can override using the
`.groups` argument.

a12_synt <- comparison_synt$a12_synt
a21_synt <- comparison_synt$a21_synt

a12_synt

[1] -0.3663187

a21_synt

[1] 0.3663187

Lotka-Voltera Interspecific Competition Model

Now let’s input our values into our Lotka-Voltera interspecific competition model:

interspec_dens_dep_growth <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1 * n[1] * (1 - a11 * n[1] - a12 * n[2])
    dn2dt <- r2 * n[2] * (1 - a22 * n[2] - a21 * n[1])
    list(c(dn1dt, dn2dt))
  })
}

Set parameters:

# Set parameters
parms <- c(r1 = r1, r2 = r2, a11 = a11, a21 = a21, a22 = a22, a12 = a12)

# Set N0
N0 <- c(N0_CN, N0_TN)

# Set times

times <- times

Run the model while adding the parameter maxsteps = 1e5:

# Run the ODE model
out <- deSolve::ode(y = N0, times = times, func = interspec_dens_dep_growth, parms = parms, maxsteps = 1e5)

Plot your results:

# Convert to dataframe for plotting
ode_df <- out %>% 
  as.data.frame() %>% 
  rename(time = 1, V1 = 2, V2 = 3)

# Plot the results
ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = V1, colour = "CN")) + 
  geom_line(aes(y = V2, colour = "TN")) + 
  labs(x = "Years", y = "Population Size", colour = "Species")

In general, are these species suppressing the populations of one another? Is competition occurring?

From this analysis alone, it appears that as TN decreases, CN increases in a seemingly density independent manner. TN may be controlling CN.

Modify the parameters to see if the models agrees with your expectations:

N0 <- c(0, mean(N0_TN_alop, N0_TN_alot, N0_TN_synt))

times <- c(1:(max(nightingales_agg$year) - min(nightingales_agg$year) + 1))

out <- deSolve::ode(y = N0, times = times, func = interspec_dens_dep_growth, parms = parms, maxsteps = 1e5)

ode_df <- out %>% 
  as.data.frame() %>% 
  rename(time = 1, V1 = 2, V2 = 3)

ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = V1, colour = "CN")) + 
  geom_line(aes(y = V2, colour = "TN")) + 
  labs(x = "Years", y = "Population Size", colour = "Species")

N0 <- c(mean(N0_CN_alop,N0_CN_alot, N0_CN_synt), 0)

times <- c(1:(max(nightingales_agg$year) - min(nightingales_agg$year) + 1))

out <- deSolve::ode(y = N0, times = times, func = interspec_dens_dep_growth, parms = parms, maxsteps = 1e5)

ode_df <- out %>% 
  as.data.frame() %>% 
  rename(time = 1, V1 = 2, V2 = 3)

ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = V1, colour = "CN")) + 
  geom_line(aes(y = V2, colour = "TN")) + 
  labs(x = "Years", y = "Population Size", colour = "Species")

Without CN present, TN population levels continues to drop. When TN is removed, CN population levels continue to grow.

Now let’s look at the results across the different co-occurrence contexts to see if these results hold.

N0.s <- list(
  alop = c(N0_CN_alop, N0_TN_alop),
  alot = c(N0_CN_alot, N0_TN_alot),
  synt = c(N0_CN_synt, N0_TN_synt)
)

a12.s <- list(
  alop = a12_alop,
  alot = a12_alot,
  synt = a12_synt
)

a21.s <- list(
  alop = a21_alop,
  alot = a21_alot,
  synt = a21_synt
)

for (condition in c("alop", "alot", "synt")) {

  N0 <- N0.s[[condition]]
  a12 <- a12.s[[condition]]
  a21 <- a21.s[[condition]]
  
  times <- c(1:(max(nightingales_agg$year) - min(nightingales_agg$year) + 1))
  
  out <- deSolve::ode(y = N0, times = times, func = interspec_dens_dep_growth, parms = parms, maxsteps = 1e5)
  
  ode_df <- out %>% 
    as.data.frame() %>% 
    rename(time = 1, V1 = 2, V2 = 3)
  
  plot <- ggplot(data = ode_df, aes(x = time)) + 
    geom_line(aes(y = V1, colour = "CN")) + 
    geom_line(aes(y = V2, colour = "TN")) + 
    labs(title = condition, x = "Years", y = "Population Size", colour = "Species")
  
  print(plot)
}

What are your observations about interspecific competition in these populations from our model?

Based on our model, what do you think the status of each population will be in 2124?

for (condition in c("alop", "alot", "synt")) {

  N0 <- N0.s[[condition]]
  a12 <- a12.s[[condition]]
  a21 <- a21.s[[condition]]
  
  times <- c(1:100)
  
  out <- deSolve::ode(y = N0, times = times, func = interspec_dens_dep_growth, parms = parms, maxsteps = 1e5)
  
  ode_df <- out %>% 
    as.data.frame() %>% 
    rename(time = 1, V1 = 2, V2 = 3)
  
  plot <- ggplot(data = ode_df, aes(x = time)) + 
    geom_line(aes(y = V1, colour = "CN")) + 
    geom_line(aes(y = V2, colour = "TN")) + 
    labs(title = condition, x = "Years", y = "Population Size", colour = "Species")
  
  print(plot)
}

Zero Growth Isoclines

nightingales_isocline <- nightingales_disag %>%
  filter(co_occurrence_status == "synt") %>% # Only selecting "synt" so that interspecies dynamics are strongest. If your dataset is not separated into alopatric/sympatric samples, you can skip this step
  group_by(plot, year, species) %>%
  summarise(count = sum(count)) %>%
  spread(key = species, value = count, fill = 0) %>%
  ungroup()

`summarise()` has grouped output by 'plot', 'year'. You can override using the
`.groups` argument.

head(nightingales_isocline)

# A tibble: 6 × 4
   plot  year    CN    TN
  <int> <int> <dbl> <dbl>
1     1  2010     1     3
2    13  2003     2     1
3    15  2007     2     0
4    15  2008     2     1
5    15  2009     1     1
6    15  2010     3     0

ggplot(nightingales_isocline, aes(x = CN, y = TN)) +
  geom_point() +
  geom_smooth(method = "lm", color = "red", se = FALSE) +  # CN x TN
  geom_smooth(aes(x = TN, y = CN), method = "lm", color = "blue", se = FALSE) +  # TN x CN
  labs(title = "Zero Growth Isoclines for CN and TN",
       x = "CN",
       y = "TN")

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

From this model, which species do you think is the dominant competitor? Do you think there is an equilibrium point?

It appear like there is an equilibrium point around 5 CN: 1 TN. Interspecific dynamics would push populations towards this equilibrium point. From this plot we could also calculate estimated K values from the intercepts.

References

Reif, Jiřı́, Radka Reifová, Anna Skoracka, and Lechosław Kuczyński. 2018. “Competition-Driven Niche Segregation on a Landscape Scale: Evidence for Escaping from Syntopy Towards Allotopy in Two Coexisting Sibling Passerine Species.” Journal of Animal Ecology 87 (3): 774–89.

Soetaert, Karline, Thomas Petzoldt, and R. Woodrow Setzer. 2010. “Solving Differential Equations in r: Package desolve” 33. https://doi.org/10.18637/jss.v033.i09.

Stevens, Martin Henry Hoffman. 2009. A Primer of Ecology with r. Springer.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2022. Dplyr: A Grammar of Data Manipulation.

Wickham, Hadley, and Maintainer Hadley Wickham. 2017. “Package ‘Tidyr’.” Easily Tidy Data with’spread’and’gather ()’Functions.