BIO261: Coursework 2
B Tooby
2024-04-13
The impact of temperature and competition coefficients on the population dynamics between a native and invasive barnacle.
Load required packages:
library(deSolve)
library(viridis)## Loading required package: viridisLitelibrary(diagram)## Loading required package: shapeThe parameters used for each model with an explanation of their terms are:
# Parameters for all initial models:
r1 = 0.9 # Intrinsic growth rate of native sp
r2 = 1.5 # Intrinsic growth rate of invasive sp
K1 = 100 # Carrying capacity of native sp
K2 = 100 # Carrying capacity of invasive sp
alpha12 = 0.8 # Interaction coefficient -invasive sp on native sp - this is altered in Model 3
alpha21 = 0.1 # Interaction coefficient -native sp on invasive sp - this is altered in Model 3
T_opt1 = 11 # Optimal temperature of native sp
T_opt2 = 16 # Optimal temperature of invasive sp
sigma1 = 2.5 # Temperature sensitivity of native sp
sigma2 = 4.5 # Temperature sensitivity of invasive sp
temp_v = seq(10,26,1) # Temperature ranges
N0 <- c(10,1) # Initial population sizes (10 = native, 1 = invasive)The interaction coefficient in visual form:
forces <- matrix(c(
0, alpha12,
alpha21, 0),
nrow = 2, byrow = TRUE)
colnames(forces) <- c("Native", "Invasive")
plotmat(forces, relsize = 0.83, cex.txt = 0.7, dtext = 0.5,, box.cex = 0.7)
Both models used in this report have a temperature dependant growth rate:
\[ \begin{align*} \text{growth rate_1} &= r_1 *\exp\left(-\frac{{(T - T_{\text{opt1}})^2}}{{2 \cdot \sigma_1^2}}\right) \\ \text{growth rate_2} &= r_2 *\exp\left(-\frac{{(T - T_{\text{opt2}})^2}}{{2 \cdot \sigma_2^2}}\right) \end{align*} \]This means that the growth rate decreases as the temperature deviates from the optimal temperature, with a bell-shaped curve centred at the optimal temperature for the species (as shown below)
# Bell curve of growth rate of Native & Invasive species vs temperature:
# Growth rate equations for both species:
growth_rate1 <- function(T1, r1, T_opt1, sigma1) {
r1 * exp(-((T1 - T_opt1)^2) / (2 * sigma1^2))
}
growth_rate2 <- function(T1, r2, T_opt2, sigma2) {
r2 * exp(-((T1 - T_opt2)^2) / (2 * sigma2^2))
}
# Temperature range
T1 <- seq(0, 40, by = 0.1) # Temperature range
# Calculate growth rates for both species
growth_rates_native <- growth_rate1(T1, r1, T_opt1, sigma1)
growth_rates_invasive <- growth_rate2(T1, r2, T_opt2, sigma2)
# Plot the growth rates of both species on the same graph
plot(T1, growth_rates_native,
type = "l",
col = "dodgerblue",
xlab = "Temperature (°C)",
ylab = "Growth Rate (per capita)",
ylim = c(0,1.7),
lwd = 2)
lines(T1, growth_rates_invasive, col = "orange", lwd = 2)
legend("topright", legend = c("Native", "Invasive"),
col = c("dodgerblue", "orange"),
lty = 1,
lwd = 2)
As shown above, the invasive species has a higher growth rate and is able to still grow at a larger range than the native species. This can be then be implemented into each model.
Model 1
Density dependant abundance change with temperature dependant growth rate
The basis for this model was the Verhulst (1838) model for continuous density dependence:
\[ \begin{align*} \\\frac{{dN}}{{dt}}\ = N r \left(1-\frac{{N}}{{K}}\right) \end{align*} \]
with \(r\) substituted for the temperature dependant growth rate formula for each species which is defined in the previous section.
Define the model & solve differential equation
# Model 1
# Define the model
independent.model <- function(times, N0, parameters, T) {
N1 <- N0[1]
N2 <- N0 [2]
with(as.list(parameters), {
# Temperature dependant growth rate
growth_rate1 <- r1 * exp(-((T - T_opt1)^2) / (2 * sigma1^2))
growth_rate2 <- r2 * exp(-((T - T_opt2)^2) / (2 * sigma2^2))
# density dependant single species growth
dn1 <- N1 * growth_rate1 * (1 - (N1 / K1 ))
dn2 <- N2 * growth_rate2 * (1 - (N2 / K2 ))
return(list(c(dn1, dn2)))
})
}
# Assign parameters:
parameters <- list(r1, r2,
K1, K2,
T_opt1, T_opt2,
sigma1, sigma2,
alpha12, alpha21)
# Length of simulation
t.values <- seq(0,100, length.out = 501)
# Create empty lists to store results for both species
species1_results <- list()
species2_results <- list()
# Solve the differential equations for each temperature and store the results
for (temp in temp_v) {
independent.model.out <- ode(N0,
t.values,
independent.model,
parms = parameters,
T = temp)
species1_results[[paste("Temp =", temp, "°C")]] <- independent.model.out[, 2]
species2_results[[paste("Temp =", temp, "°C")]] <- independent.model.out[, 3]
}Investigate abundance vs time series
# Investigate time series
# Plot results for species 1 (native)
par(mfrow=c(1, 1)) # Arrange plots 1x1
matplot(t.values,
sapply(species1_results, function(x) x), # Extract species 1 abundance from the list
type = "l",
lty = 1,
lwd = 2,
ylab = "Abundance (N)",
xlab = "Time (t)",
cex.main = 0.9,
ylim = c(0, 110),
xlim = c(0, 110),
col = viridis(length(temp_v), option = "H"),
main = "Native")
legend("topright",
legend = paste(temp_v, "°C"),
lty = 1,
lwd = 2,
col = viridis(length(temp_v), option = "H"),
cex = 0.6,
title = "Temp")
# Plot results for species 2 (invasive)
matplot(t.values,
sapply(species2_results, function(x) x), # Extract species 2 abundance from the list
type = "l",
lty = 1,
lwd = 2,
ylab = "Abundance (N)",
xlab = "Time (t)",
cex.main = 0.9,
ylim = c(0, 110),
xlim = c(0, 110),
col = viridis(length(temp_v), option = "H"),
main = "Invasive")
legend("topright",
legend = paste(temp_v, "°C"),
lty = 1,
lwd = 2,
col = viridis(length(temp_v), option = "H"),
cex = 0.6,
title = "Temp")
Investigate changes in growth rate over time with temperature
par(mfrow=c(1,1))
# For Native Species
# Plot results for species 1 (native)
plot(NULL,
xlim = c(0, 105),
ylim = c(-2, 8),
main = "Native: Growth Rate change over time",
xlab = "Time (t)",
ylab = "Growth Rate (N(t+1)- N(t)")
# Plot results for species 1 (native) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
independent.model.out <- ode(N0,
t.values,
independent.model,
parms = parameters,
T = temp_v[i])
# Calculate growth rate for species 1 (native)
growth_rateN <- diff(independent.model.out[, 2])
# Add lines for species 1 (native)
lines(t.values[-1], growth_rateN,
col = viridis(length(temp_v), option = "H")[i],
lwd = 2,
lty = 1)
}
legend("topright", legend = paste(temp_v, "°C"),
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temperature",
cex = 0.6)
### Invasive Species
# Plot results for species 2 (invasive)
plot(NULL,
xlim = c(0, 105),
ylim = c(-2, 8),
main = "Invasive: Growth Rate change over time",
xlab = "Time (t)",
ylab = "Growth Rate (N(t+1)- N(t)")
# Plot results for species 2 (invasive) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
independent.model.out <- ode(N0,
t.values,
independent.model,
parms = parameters,
T = temp_v[i])
# Calculate growth rate for species 2 (invasive)
growth_rateI <- diff(independent.model.out[, 3])
# Add lines for species 2 (invasive)
lines(t.values[-1], growth_rateI,
col = viridis(length(temp_v), option = "H")[i],
lwd = 2,
lty = 1)
}
legend("topright", legend = paste(temp_v, "°C"),
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temperature",
cex = 0.6)
Investigate population growth maps and per capita growth rates at different temperatures
For native species:
# Create empty lists to store population growth functions and per-capita growth rates
f_N <- list()
pgrN <- list()
# Loop over temperature values
for (temp in temp_v) {
# Create population growth function for the current temperature
N1 <- seq(0, 1.5 * K1, length.out = 101)
growth_rate <- r1 * exp(-((temp - T_opt1)^2) / (2 * sigma1^2))
f_N[[temp]] <- exp(growth_rate * (1 - N1 / K1))
# Calculate per-capita growth rates
pgrN[[temp]] <- log(f_N[[temp]])
}
# Plot the results
par(mfrow = c(1, 1))
plot(NULL, xlim = c(0, 1.5 * K1), ylim = c(0, 1.5 * K1),
xlab = "Population size, N[t]",
ylab = "Population growth function, N[t+1] = N[t] f(N[t])",
main = "Population Growth Function Native")
# Add lines for each temperature
for (i in 1:length(temp_v)) {
lines(N1, N1 * f_N[[temp_v[i]]], col = viridis(length(temp_v), option = "H")[i])
}
legend("topright",
legend = paste(temp_v, "°C"),
lty = 1,
lwd = 2,
col = viridis(length(temp_v), option = "H"),
cex = 0.6,
title = "Temp °C")
grid()
# Plot per-capita growth rates for each temperature
plot(NULL, xlim = c(0, 1.5 * K1), ylim = c(-2, 2),
xlab = "Population size, N[t]",
ylab = "Per-capita growth rate",
main = "Per-capita Growth Rate Native")
for (i in 1:length(temp_v)) {
lines(N1, pgrN[[temp_v[i]]], col = viridis(length(temp_v), option = "H")[i])
}
legend("topright",
legend = paste(temp_v, "°C"),
lty = 1,
lwd = 2,
col = viridis(length(temp_v), option = "H"),
cex = 0.6,
title = "Temp °C")
grid()
For invasive species:
# Create empty lists to store population growth functions and per-capita growth rates
f_I <- list()
pgrI <- list()
# Loop over temperature values
for (temp in temp_v) {
# Create population growth function for the current temperature
N2 <- seq(0, 1.5 * K2, length.out = 101)
growth_rate2 <- r2 * exp(-((temp - T_opt2)^2) / (2 * sigma2^2))
f_I[[temp]] <- exp(growth_rate2 * (1 - N2 / K2))
# Calculate per-capita growth rates
pgrI[[temp]] <- log(f_I[[temp]])
}
# Plot the results
library(viridis)
par(mfrow = c(1, 1))
plot(NULL, xlim = c(0, 1.5 * K2), ylim = c(0, 1.5 * K2),
xlab = "Population size, N[t]",
ylab = "Population growth function, N[t+1] = N[t] f(N[t])",
main = "Population Growth Function Invasive")
# Add lines for each temperature
for (i in 1:length(temp_v)) {
lines(N2, N2 * f_I[[temp_v[i]]], col = viridis(length(temp_v), option = "H")[i])
}
legend("topright",
legend = paste(temp_v, "°C"),
lty = 1,
lwd = 2,
col = viridis(length(temp_v), option = "H"),
cex = 0.6,
title = "Temp °C")
grid()
# Plot per-capita growth rates for each temperature
plot(NULL, xlim = c(0, 1.5 * K2), ylim = c(-2, 2),
xlab = "Population size, N[t]",
ylab = "Per-capita growth rate",
main = "Per-capita Growth Rate Invasive")
for (i in 1:length(temp_v)) {
lines(N2, pgrI[[temp_v[i]]], col = viridis(length(temp_v), option = "H")[i])
}
legend("topright",
legend = paste(temp_v, "°C"),
lty = 1,
lwd = 2,
col = viridis(length(temp_v), option = "H"),
cex = 0.6,
title = "Temp °C")
grid()
Model 2
Inter-specific competition and the effect on abundance with temperature dependant growth rates
To determine the effect of competition on the populations of both species at different temperatures, Lotka Voltera inter-specific competition was uesed as the basis: \[ \begin{align*} \frac{dN_1}{dt} &= r_1 N_1 \left( \frac{K_1 - N_1 - \alpha_{12} N_2}{K_1} \right) \\ \frac{dN_2}{dt} &= r_2 N_2 \left( \frac{K_2- N_2 - \alpha_{21} N_1}{K_2} \right) \\ \end{align*} \]
with the additional growth dependent term used in the previous model:\[ \begin{align*} \text{growth rate_1} &= r_1 \cdot \exp\left(-\frac{{(T - T_{\text{opt1}})^2}}{{2 \cdot \sigma_1^2}}\right) \\ \text{growth rate_2} &= r_2 \cdot \exp\left(-\frac{{(T - T_{\text{opt2}})^2}}{{2 \cdot \sigma_2^2}}\right) \end{align*} \]
State Space graph
To determine the outcome of the competition: results show an equilibrium between the two species.
par(mfrow=c(1,1))
# Create an empty plot
plot(0, type = "n",
xlim = c(0, 1000),
ylim = c(0, 130),
xlab = "Native Abundance",
ylab = "Invasive Abundance")
# Determine intercept point
N_isocline <- K1/alpha12 # intercept of native population on the y axis
I_isocline <- K2/alpha21 # intercept of the invasive species on the x axis
# Add lines representing the nullcine
lines(c(0, K1), c(N_isocline, 0), lwd= 2,
col = "dodgerblue") # Line for species 1
lines(c(I_isocline, 0), c(0, K2), lwd = 2,
col = "orange") # Line for species 2
# Add legend
legend("topright",
legend = c("Native","Invasive "),
lty = 1,
lwd = 2,
col = c("dodgerblue", "orange"))
# Plot arrows showing direction
arrows(x0 = 90, y0 = 70, x1 = 40, y1 = 90, length = 0.1) # Bottom right arrow
arrows(x0 = -15, y0 = 115, x1 = 8, y1 = 105, length = 0.1) # Top left arrow
arrows(x0 = -10, y0 = 85, x1 = 10, y1 = 94, length = 0.1) # Bottom left arrow
arrows(x0 = 70, y0 = 110, x1 = 40, y1 = 100, length = 0.1) # Top right arrow
# Key features of each nullcline:
text(x= 85, y = 0, "K1", cex = 0.5, col = "gray29")
text(x= -5, y = 96, "K2", cex = 0.5, col = "grey29")
text(x = 1000, y = 20, expression(frac(K2, alpha_21)), cex = 0.4, col = "grey29")
text(x = 50, y = 123, expression(frac(K1, alpha_12)), cex = 0.4, col = "grey29")
# Equilibrium point
points(21, 98, pch = 21, col = "black", cex = 1.2)
Define Model 2
# Define model
competitive.lotka.volterra <- function(times, N0, parameters, T) {
N1 <- N0[1]
N2 <- N0[2]
with(as.list(parameters), {
# Temperature dependant growth rate (same as previous section)
growth_rate1 <- r1 * exp(-((T - T_opt1)^2) / (2 * sigma1^2))
growth_rate2 <- r2 * exp(-((T - T_opt2)^2) / (2 * sigma2^2))
# Lotka-Voltera inter-specific competition equation
dN1.dt <- growth_rate1 * N1 *((K1-N1-(alpha12*N2))/K1)
dN2.dt <- growth_rate2 * N2 *((K2-N2-(alpha21*N1))/K2)
return(list(c(dN1.dt, dN2.dt)))
})
}Confirm the initial conditions and parameters, including initial populations (10 for native species as they are already there and 1 for the invasive as they are ‘newer’), and the time vector.
# Set initial conditions and parameters
N0 <- c(10, 1) # Initial population for each species
parameters <- list(r1, r2,
K1, K2,
T_opt1, T_opt2,
sigma1, sigma2,
alpha12, alpha21)
# Time vector
t.values <- seq(0,100, length.out = 501)Then solve the differential equations for each temperature and plot them. First we plot with one temperature per graph to understand the relationship between the two species.
# Plot results for each temperature
par(mfrow=c(2, 2)) # Arrange plots in a 2x2 grid
for (temp in temp_v) {
# Solve the differential equations for each temperature
lotka.volterra.out <- ode(N0,
t.values,
competitive.lotka.volterra,
parms = parameters,
T = temp)
# Plot results
matplot(t.values, lotka.volterra.out[, c(2, 3)],
type = "l",
lty = 1,
lwd = 2,
ylab = "Abundance (N)",
xlab = "Time (t)",
cex.main = 0.9,
ylim = c(0,150),
col = c("dodgerblue", "orange"))
legend("topright",
c("Native", "Invasive"),
lty = 1,
lwd = 2,
col = c("dodgerblue", "orange"),
cex = 0.9,
title = paste("Temp:", as.character(temp), "°C"))
}
Then plot one graph for each species (with all temperatures on each one).
Abundance vs Time graphs:
par(mfrow = c(1, 1)) # Arrange graphs 1 x 1
# Plot results for species 1 (native)
plot(NULL, xlim = c(0, 105),
ylim = c(-10, 110),
xlab = "Time (t)",
ylab = "Abundance (N)",
main = "Native")
# Plot results for species 1 (native) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.out <- ode(N0,
t.values,
competitive.lotka.volterra,
parms = parameters,
T = temp_v[i])
# Add lines for species 1 (native)
lines(t.values, lotka.volterra.out[, 2],
col = viridis(length(temp_v), option = "H")[i],
lwd = 2,
lty = 1,
type = "l")
}
# Add legend for species 1 (native)
legend("topright", legend = temp_v,
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temp °C",
cex = 0.6)
# Species 2 (invasive)
# Plot results for species 2 (invasive)
plot(NULL, xlim = c(0, 105),
ylim = c(0, 110),
xlab = "Time (t)",
ylab = "Abundance (N)",
main = "Invasive")
# Plot results for species 2 (invasive) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.out <- ode(N0,
t.values,
competitive.lotka.volterra,
parms = parameters,
T = temp_v[i])
# Add lines for species 2 (invasive)
lines(t.values, lotka.volterra.out[, 3],
col = viridis(length(temp_v), option = "H")[i],
lwd = 2,
lty = 1,
type = "l")
}
# Add legend for species 2 (invasive)
legend("topright", legend = temp_v,
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temp °C",
cex = 0.6)
Interspecific competition population growth over time at different temperatures
Plotting growth rate change over time:
par(mfrow=c(1,1))
# For Native Species
# Plot results for species 1 (native)
plot(NULL,
xlim = c(0, 105),
ylim = c(-2, 8),
main = "Native: Growth Rate change over time",
xlab = "Time (t)",
ylab = "Growth Rate (N(t+1)- N(t)")
# Plot results for species 1 (native) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.out <- ode(N0,
t.values,
competitive.lotka.volterra,
parms = parameters,
T = temp_v[i])
# Calculate growth rate for species 1 (native)
growth_rate <- diff(lotka.volterra.out[, 2])
# Add lines for species 1 (native)
lines(t.values[-1], growth_rate,
col = viridis(length(temp_v), option = "H")[i],
lwd = 2,
lty = 1)
}
legend("topright", legend = temp_v,
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temp °C",
cex = 0.6)
### Invasive Species
# Plot results for species 2 (invasive)
plot(NULL,
xlim = c(0, 105),
ylim = c(-2, 8),
main = "Invasive: Growth Rate change over time",
xlab = "Time (t)",
ylab = "Growth Rate (N(t+1)- N(t)")
# Plot results for species 2 (invasive) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.out <- ode(N0,
t.values,
competitive.lotka.volterra,
parms = parameters,
T = temp_v[i])
# Calculate growth rate for species 2 (invasive)
growth_rate <- diff(lotka.volterra.out[, 3])
# Add lines for species 2 (invasive)
lines(t.values[-1], growth_rate,
col = viridis(length(temp_v), option = "H")[i],
lwd = 2,
lty = 1)
}
legend("topright", legend = temp_v,
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temp °C",
cex = 0.6)
Model 3:
Investigate changes in interaction coefficient:
First, cycle through model 2 with interaction coefficients and plot to determine extinction:
Then determine the alpha12 and alpha21 values where both species are extinct at different temperatures and look closer at them.
# Vary alpha12 and alpha21 sequence to work out when there is extinction for either species.
# Starting with:
alpha12a = seq (0.5, 2, 0.5)
alpha21a = seq (0.1, 2, 0.5)
par(mfrow=c(1,3))
# Define model
competitive.lotka.volterraC <- function(times, N0, parameters, T) {
N1 <- N0[1]
N2 <- N0[2]
with(as.list(parameters), {
# Temperature dependant growth rate (same as previous section)
growth_rate1 <- r1 * exp(-((T - T_opt1)^2) / (2 * sigma1^2))
growth_rate2 <- r2 * exp(-((T - T_opt2)^2) / (2 * sigma2^2))
# Lotka-Voltera inter-specific competition equation
dN1.dtC <- growth_rate1 * N1 *((K1-N1-(alpha12a*N2))/K1) # amended interaction coefficient
dN2.dtC <- growth_rate2 * N2 *((K2-N2-(alpha21a*N1))/K2) # amended interaction coefficient
return(list(c(dN1.dtC, dN2.dtC)))
})
}
###For native species 1:
# Loop through alpha12 and alpha21 combinations
for (a12 in alpha12a) {
for (a21 in alpha21a) {
# Set parameters for this combination
parameters <- list(r1, r2,
K1, K2,
T_opt1, T_opt2,
sigma1, sigma2,
alpha12a = a12, # new coefficients
alpha21a = a21) # new coefficients
# Plot results for species 1 (native)
plot(NULL, xlim = c(0, 105),
ylim = c(-10, 100),
xlab = "Time (t)",
ylab = "Abundance (N)",
main = paste("N a12=", a12, "a21=",a21))
# Plot results for species 1 (native) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.outC <- ode(N0,
t.values,
competitive.lotka.volterraC,
parms = parameters,
T = temp_v[i])
# Add lines for species 1 (native)
lines(t.values, lotka.volterra.outC[, 2],
col = viridis(length(temp_v), option = "H")[i],
lwd = 1,
lty = 1,
type = "l")
}
}
}
### Invasive species 2
# Loop through alpha12 and alpha21 combinations
for (a12 in alpha12a) {
for (a21 in alpha21a) {
# Set parameters for this combination
parameters <- list(r1, r2,
K1, K2,
T_opt1, T_opt2,
sigma1, sigma2,
alpha12a = a12, # new coefficients
alpha21a = a21) # new coefficients
# Plot results for species 2 (invasive)
plot(NULL, xlim = c(0, 120),
ylim = c(-10, 100),
xlab = "Time (t)",
ylab = "Abundance (N)",
main = paste("I a12=", a12, "a21=",a21))
# Plot results for species 2 (invasive) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.outC <- ode(N0,
t.values,
competitive.lotka.volterraC,
parms = parameters,
T = temp_v[i])
# Add lines for species 2 (invasive)
lines(t.values, lotka.volterra.outC[, 3],
col = viridis(length(temp_v), option = "H")[i],
lwd = 1,
lty = 1,
type = "l")
}
legend("topright", legend = temp_v,
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temp °C",
cex = 0.6)
}
}
extinction is recorded for species one (NATIVE) with the following interaction coefficients:
a12 = 1 a21 = 0.1 # all
a12 = 1 a21 = 0.6 # all
a12 = 1.5 a21 = 0.1 # all
a12 = 1.5 a21 = 0.6 # all
a12 = 1.5 a21 = 1.1 # not lowest 3 temps
a12 = 1.5 a21 = 1.6 # not lowest 5 temps
a12 = 2 a21 = 0.1 # all
a12 = 2 a21 = 0.6 # all
a12 = 2 a21 = 1.1 # not lowest 2 temps
a12 = 2 a21 = 1.6 # not lowest 5 temps
Extinction in species 2 (INVASIVE)
a12 = 0.5 a21 = 1.1 # only 7 lowest temps
a12 = 0.5 a21 = 1.6 # only lowest 6 temps
a12 = 1 a21 = 1.1 # only lowest 6 temps
a12 = 1 a21 = 1.6 # only lowest 5 temps
a12 = 1.5 a21 = 1.1 # only lowest 3 temps
a12 = 1.5 a21 = 1.6 # only lowest 4 temps
a12 = 2 a21 = 1.1 # only lowest 2 temps
a12 = 2 a21 = 1.6 # only lowest 4 temps
# Look at interesting dynamics where both are extinct at different temperatures (remove # from one variation at a time to produce graphs)
# Variation 1. a12 = 1.5, a21 = 1.1
#alpha12a = 1.5
#alpha21a = 1.1
# Variation 2. a12 = 1.5, a21 = 1.6
# alpha12a = 1.5
# alpha21a = 1.6
# Variation 3. a12 = 2, a21 = 1.1 - utilised in report
alpha12a= 2
alpha21a= 1.1
# Variation 4 a12 = 2, a21 = 1.6
# alpha12a = 2
# alpha21a = 1.6
par(mfrow= c(1,2))
# Define model
competitive.lotka.volterraC <- function(times, N0, parameters, T) {
N1 <- N0[1]
N2 <- N0[2]
with(as.list(parameters), {
# Temperature dependant growth rate (same as previous section)
growth_rate1 <- r1 * exp(-((T - T_opt1)^2) / (2 * sigma1^2))
growth_rate2 <- r2 * exp(-((T - T_opt2)^2) / (2 * sigma2^2))
# Lotka-Voltera inter-specific competition equation
dN1.dtC <- growth_rate1 * N1 *((K1-N1-(alpha12a*N2))/K1) # amended interaction coefficient
dN2.dtC <- growth_rate2 * N2 *((K2-N2-(alpha21a*N1))/K2) # amended interaction coefficient
return(list(c(dN1.dtC, dN2.dtC)))
})
}
###For native species 1:
# Loop through alpha12 and alpha21 combinations
for (a12 in alpha12a) {
for (a21 in alpha21a) {
# Set parameters for this combination
parameters <- list(r1, r2,
K1, K2,
T_opt1, T_opt2,
sigma1, sigma2,
alpha12a = a12, # new coefficients
alpha21a = a21) # new coefficients
# Plot results for species 1 (native)
plot(NULL, xlim = c(0, 105),
ylim = c(-10, 100),
xlab = "Time (t)",
ylab = "Abundance (N)",
main = paste("N a12=", a12, "a21=",a21))
# Plot results for species 1 (native) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.outC <- ode(N0,
t.values,
competitive.lotka.volterraC,
parms = parameters,
T = temp_v[i])
# Add lines for species 1 (native)
lines(t.values, lotka.volterra.outC[, 2],
col = viridis(length(temp_v), option = "H")[i],
lwd = 1,
lty = 1,
type = "l")
}
}
}
### Invasive species 2
# Loop through alpha12 and alpha21 combinations
for (a12 in alpha12a) {
for (a21 in alpha21a) {
# Set parameters for this combination
parameters <- list(r1, r2,
K1, K2,
T_opt1, T_opt2,
sigma1, sigma2,
alpha12a = a12, # new coefficients
alpha21a = a21) # new coefficients
# Plot results for species 2 (invasive)
plot(NULL, xlim = c(0, 120),
ylim = c(-10, 100),
xlab = "Time (t)",
ylab = "Abundance (N)",
main = paste("I a12=", a12, "a21=",a21))
# Plot results for species 2 (invasive) for each temperature
for (i in 1:length(temp_v)) {
# Solve the differential equations for each temperature
lotka.volterra.outC <- ode(N0,
t.values,
competitive.lotka.volterraC,
parms = parameters,
T = temp_v[i])
# Add lines for species 2 (invasive)
lines(t.values, lotka.volterra.outC[, 3],
col = viridis(length(temp_v), option = "H")[i],
lwd = 1,
lty = 1,
type = "l")
}
legend("topright", legend = temp_v,
col = viridis(length(temp_v), option = "H"),
lwd = 2,
title = "Temp °C",
cex = 0.6)
}
}
Phase plane graph for different competition coefficients:
alpha12a = 2
alpha21a = 1.1
Results in an unstable equilibrium.
par(mfrow=c(1,1))
alpha12a = 2
alpha21a = 1.1
# Create an empty plot
plot(0, type = "n",
xlim = c(0, 110),
ylim = c(0, 125),
xlab = "Native Abundance",
ylab = "Invasive Abundance")
N_isocline <- K1/alpha12a
I_isocline <- K2/alpha21a
# Add lines representing the nullcine
lines(c(0, K1), c(N_isocline, 0), lwd= 2,
col = "dodgerblue") # Line for species 1
lines(c(I_isocline, 0), c(0, K2), lwd = 2,
col = "orange") # Line for species 2
# Add legend
legend("topright",
legend = c("Native","Invasive "),
lty = 1,
lwd = 2,
col = c("dodgerblue", "orange"))
# Plot arrows showing direction
arrows(x0 = 90, y0 = 3, x1 = 98, y1 = -3, length = 0.1) # Bottom right arrow
arrows(x0 = 75, y0 = 15, x1 = 65, y1 = 23, length = 0.1) # Top left arrow
arrows(x0 = 70, y0 = -4, x1 = 82, y1 = 6, length = 0.1) # Bottom left arrow
arrows(x0 = 90, y0 = 40, x1 = 87, y1 = 20, length = 0.1) # Top right arrow
# Key features of each nullcline:
text(x= 105, y = 0, "K1", cex = 0.5, col = "gray29")
text(x= 0, y = 105, "K2", cex = 0.5, col = "grey29")
text(x = 0, y = 60, expression(frac(K1, alpha_12a)), cex = 0.4, col = "grey29")
# Equilibrium point
points(84, 8, pch = 21, col = "black", cex = 1.2)