Interspecies Interactions: 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}}) \]

$$ = r_{1} N_{1} (1 - - )

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

\[ \alpha_{11} = \frac{1}{K_{1}} \]

We can also combine terms into interspecific density dependence or “interspecific competition coefficient” $\alpha _{12}$ which can be read as “the impact of species 2 on species 1.”

\[ \alpha_{12} = \frac{\alpha}{K_{1}} \]

For example:

If 1 aquatic turtle (species 1) has the same impact on resources as 4 fish (species 2) would, its presence affects species 2 as much as 4 individual fish would. We can write $\alpha_{21}=4.0$
Conversely, if the effect of species 2 on species 1 is reciprocal, i.e. one fish has the same resource impact as one-third of a turtle, then $\alpha_{12}=0.25$.

Substituting these in, 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_{2} \ -\ \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 n[1] and n[2].

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 out <- deSolve::ode(y = c(N0_1, N0_2), times = times, func = interspec_dens_dep_growth, parms = parms) with the following arguments:

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

parms <- c(
  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

out <- deSolve::ode(y = c(N0_1, N0_2), times = times, func = interspec_dens_dep_growth, parms = parms)

Now plot your results:

Hint: ggplot2 doesn’t like column names that are numeric, so include rename(times = 1, SP1 = 2, SP2 = 3) in your pipe

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


ggplot(data = ode_df, aes(x = times)) + 
  geom_line(aes(y = SP1, color = "Species 1")) + 
  geom_line(aes(y = SP2, color = "Species 2")) + 
  labs(x = "Years", y = "Population Size", color = "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?


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?

parms <- c(
  r1 <- 1,
  r2 <- 0.2,
  a11 <- 0.2,
  a21 <- 0.1,
  a22 <- 0.02,
  a12 <- 0.03
)

N0_1 <- 2
N0_2 <- 1

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

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


ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = SP1, color = "Species 1")) + 
  geom_line(aes(y = SP2, color = "Species 2")) + 
  labs(x = "Years", y = "Population Size", color = "Species")


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.

Which types of relationships are these different values reflecting?

Neutralism: \[ a_{12} = 0, \quad a_{21} = 0 \]
- No direct influence of either species on the other. They coexist without interaction.
Amensalism: \[ a_{12} < 0, \quad a_{21} = 0 \]
- One species is negatively affected while the other remains unaffected.
Commensalism: \[ a_{12} > 0, \quad a_{21} = 0 \]
- One species benefits without impacting the other (e.g., species 1 benefits while species 2 is unaffected).
Competition: \[ a_{12} < 0, \quad a_{21} < 0 \]
- Both species negatively impact each other, leading to a struggle for resources.
Antagonism: \[ a_{12} < 0, \quad a_{21} > 0 \]
- One species benefits at the expense of the other (e.g., predator-prey, host-parasite, consumer-resource interactions).
Mutualism: \[ a_{12} > 0, \quad a_{21} > 0 \]
- Both species benefit from their interaction (e.g., symbiosis, cooperative relationships).

Neutral Equilibria

Set your parameters to the following conditions:

r1 <- r2 <- 1
a11 <- a21 <- 0.01
a22 <- a12 <- 0.015
N0_1 <- 2
N0_2 <- 1

What do you observe? Is this result stable or unstable? Or something more nuanced? Explain. Hint: Try varying the initial population sizes and observe what happens.

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

N0_1 <- 2
N0_2 <- 1

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

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


ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = SP1, color = "Species 1")) + 
  geom_line(aes(y = SP2, color = "Species 2")) + 
  labs(x = "Years", y = "Population Size", color = "Species")


We observe something called a "neutral saddle." 
Even though the ratio between the population sizes of the two species will remain constant, 
the relationship is not considered stable,
because after perturbation a new equilibrium point will be established for each species' population size. 
It's not considered stable because the system doesn't return to the original state. 
However, it also doesn't lead to runaway growth or extinction unless further disturbed.

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. You may need to use the read_delim() function and define the delim paramater in your data loading function) Briefly explore the dataset:

nightingales_raw <- read_delim("Lab_5_data/nightingales.csv", delim = ";")

Rows: 5983 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ";"
chr (1): co.occ
dbl (7): plot, count, year, env.grad, obs, lon, lat

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(nightingales_raw)

# A tibble: 6 × 8
   plot count co.occ   year env.grad   obs   lon   lat
  <dbl> <dbl> <chr>   <dbl>    <dbl> <dbl> <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

summary(nightingales_raw)

      plot           count            co.occ               year     
 Min.   :  1.0   Min.   : 0.0000   Length:5983        Min.   :2000  
 1st Qu.:107.0   1st Qu.: 0.0000   Class :character   1st Qu.:2005  
 Median :256.0   Median : 0.0000   Mode  :character   Median :2009  
 Mean   :260.3   Mean   : 0.9402                      Mean   :2008  
 3rd Qu.:421.0   3rd Qu.: 1.0000                      3rd Qu.:2012  
 Max.   :546.0   Max.   :12.0000                      Max.   :2014  
    env.grad            obs             lon             lat       
 Min.   :-195.02   Min.   :  1.0   Min.   :14.23   Min.   :49.20  
 1st Qu.: -27.08   1st Qu.:223.0   1st Qu.:17.74   1st Qu.:51.12  
 Median :  60.42   Median :363.0   Median :19.24   Median :51.91  
 Mean   :  74.05   Mean   :324.3   Mean   :19.28   Mean   :52.01  
 3rd Qu.: 170.63   3rd Qu.:422.0   3rd Qu.:20.89   3rd Qu.:52.82  
 Max.   : 424.87   Max.   :502.0   Max.   :23.90   Max.   :54.81

Clean the data using the separate() function from the tidyr(Wickham and Wickham 2017) package to separate species and co_occurence_status. Assign to variable nightingales_disag as it is best practice to keep the raw data intact as its own variable.

# Cleaning data

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


head(nightingales_disag)

# A tibble: 6 × 9
   plot count species co_occurrence_status  year env.grad   obs   lon   lat
  <dbl> <dbl> <chr>   <chr>                <dbl>    <dbl> <dbl> <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

To help us visualize our dataset, let’s map our data and facet by co_occurence_status. We can make a simple map using ggplot2.

We can add the facet_wrap() function to our ggplot() workflow to rapidly generate plots for all orders at once following this format:

facet_wrap(~ <categorical variable>, scales = <>)

The scales parameter can be used to determine if y-axes should be independent (free_y), x-axes should be independent (free_x), or if all axes should be independent (free) between all plots. The default is fixed. Which of these do you think would be the best in t543his case?

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

Because we are more interested in inter-year trends than intra-year trends, 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), .groups = "drop") %>%
  ungroup()

# 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.

# Note that it is good practice to use `.groups = "drop"` with `summarize()` because it can cause lasting, unintentional grouping like `group_by()` does.

nightingales_agg

# A tibble: 90 × 4
   species  year co_occurrence_status avg_count
   <chr>   <dbl> <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?

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

min_time <- min(nightingales_agg$year)
max_time <- max(nightingales_agg$year)

times <- c(1:(max_time - min_time + 1))

times

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

Initial Population Sizes

# I am going to assume I can just average across the different groups to calculate N0

N0_CN <- nightingales_agg %>%
  filter(species == "CN") %>%
  summarize(N0 = mean(avg_count[year == min(year)]), .groups = "drop") %>%
  pull(N0)

N0_CN

[1] 0.6736111

# If your dataset has one value for the first time point, you would not need to average:
# N0_CN <- nightingales_agg %>%
#   arrange(year) %>%
#   summarize(N0 = first(avg_count), .groups = "drop") %>%
#   pull(N0)


N0_TN <- nightingales_agg %>%
  filter(species == "TN") %>%
  summarize(N0 = mean(avg_count[year == min(year)]), .groups = "drop") %>%
  pull(N0)

N0_TN

[1] 1.419225

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}})}{t_{final}-{t}} \]

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


r2 <- nightingales_agg %>%
  arrange(species, year) %>%
  filter(co_occurrence_status == "alop", species == "TN") %>%
  summarize(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.

How do we plot intrinsic growth rate?

Let’s go back to our population growth lab! How did we model continuous exponential population growth?

\[ Ne^{rt} \]

Nt_CN <- N0_CN * exp(r1 * times)

Nt_CN_df <- data.frame(
  year = min_time + times - 1,  # Convert 'times' back to actual year values
  Population = Nt_CN,
  species = "CN" # This will just help the ggplot legends sync up between our datasets
)

Nt_TN <- N0_TN * exp(r2 * times)

Nt_TN_df <- data.frame(
  year = min_time + times - 1, 
  Population = Nt_TN,
  species = "TN"
)

nightingales_agg %>%
  filter(co_occurrence_status == "alop") %>%
  ggplot(aes(x = year, y = avg_count, color = species)) +
  geom_point() +
  geom_line(data = Nt_CN_df, aes(x = year, y = Population, color = species)) +
  geom_line(data = Nt_TN_df, aes(x = year, y = Population, color = species)) +
  labs(title = "Population Size", x = "Year", y = "Population Size")

Although not an excellent fit, for the purpose of this course we will assume these are good enough for our oversimplified analyses.

Intraspecific Density Dependencies

For intraspecific density dependencies ($a_{11}$ and $a_{22}$), what do you think a good estimate would be?

Back to population growth lab: density-dependent growth!

\[ rN(1-\frac{N}{K}) \] For this lab, we will just assumed carrying capacity is the maximum size of each population.

K_CN <- nightingales_agg %>%
  filter(species == "CN") %>%
  summarize(K = max(avg_count), .groups = "drop") %>%
  pull(K)

K_CN

[1] 2.608696

K_TN <- nightingales_agg %>%
  filter(species == "TN") %>%
  summarize(K = max(avg_count), .groups = "drop") %>%
  pull(K)

K_TN

[1] 1.9375

We can use one of the first equations ever used to express logistic growth to calculate our intraspecific competition terms:

\[ \alpha = \frac{r}{K} \]

a11 <- r1 / K_CN
a22 <- r2 / K_TN

a11

[1] 0.02270646

a22

[1] -0.004542538

Interspecific Negative Density Dependencies

As a very simplistic way to estimate interspecific negative density dependencies if you do not have allopatric vs. sympatric data, you can still calculate $a_{12}$ and $a_{21}$ 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
  summarize(avg = mean(avg_count)) %>%
  pull(avg)

TN_avg <- nightingales_agg %>%
  filter(species == "TN", co_occurrence_status == "synt") %>% 
  summarize(avg = mean(avg_count)) %>%
  pull(avg)

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

a12

[1] 1.522671

a21

[1] 0.6567406

Remember the above assumes that each species has a reciprocal effect on the other, which may not be the case in nature.

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 = 1e3:

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

Plot your results:

# Convert to dataframe for plotting
ode_df <- out %>% 
  as.data.frame() %>% 
  rename(time = 1, SP1 = 2, SP2 = 3)
  
# Plot the results
ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = SP1, color = "CN")) + 
  geom_line(aes(y = SP2, color = "TN")) + 
  labs(x = "Years", y = "Population Size", color = "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:

# Run with one species absent

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

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

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

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

# Run with the other species absent

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

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

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

# Plot the results
ggplot(data = ode_df, aes(x = time)) + 
  geom_line(aes(y = SP1, colour = "CN")) + 
  geom_line(aes(y = SP2, 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.

What if we run our model for 1000 years, what do you observe?

times <- c(1:5000)

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

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

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

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

We eventually observe stable equilibrium. 
Just like we saw with predator/prey models, CN is being controlled by TN. 
CN is going through more unstable fluctatuions and is the leading term, 
just like prey was for our predator/prey models. 
Eventually we see stable levels being reached.

Zero Growth Isoclines

K_CN <- nightingales_agg %>%
  filter(species == "CN") %>%
  summarize(K = max(avg_count), .groups = "drop") %>%
  pull(K)

K_CN

[1] 2.608696

K_CN/a12

[1] 1.713236

K_TN <- nightingales_agg %>%
  filter(species == "TN") %>%
  summarize(K = max(avg_count), .groups = "drop") %>%
  pull(K)

K_TN

[1] 1.9375

K_TN/a21

[1] 2.950176

ggplot() +
  geom_segment(aes(x = 0, xend = K_CN, y = (K_CN/a12), yend = 0, color = "red")) +
  geom_segment(aes(x = 0, xend = (K_TN/a21), y = K_TN, yend = 0, color = "blue")) +
  labs(x = "Common Nightingales", y = "Thrush Nightingale")

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


It appears like TN is indeed outcompeting CN.
This is in agreement with our previous models.

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, and Maintainer Hadley Wickham. 2017. “Package ‘Tidyr’.” Easily Tidy Data with’spread’and’gather ()’Functions.