## --- Problem 1(a)---
x <- seq(15,75,5) #last argument is stepcount of 5
x
## [1] 15 20 25 30 35 40 45 50 55 60 65 70 75
## --- Problem 1(b)---
squared_elements <-x[3:6]^2 #use [] to subset from 3 to 6 elem.
squared_elements
## [1] 625 900 1225 1600
## --- Problem 1(c)---
x_no7element <- x[-7] #use [-] to remove 7 element
x_no7element
## [1] 15 20 25 30 35 40 50 55 60 65 70 75
## --- Problem 1(d)---
index <- which (x>42)[1] #[1] is for first match
#that exceeds 42
index
## [1] 7
## --- Problem 1(e)---
method1 <- x[c(1:5, 9:13)]
method2 <- x[c(seq(1, 5), seq(9, 13))]
method3 <- x[c(1,2,3,4,5,9,10,11,12,13)]
method1 #print results to verify methods work
## [1] 15 20 25 30 35 55 60 65 70 75
method2
## [1] 15 20 25 30 35 55 60 65 70 75
method3
## [1] 15 20 25 30 35 55 60 65 70 75
## --- Problem 1(f)---
x_over45 <- sort(x[x>45], decreasing=TRUE)
x_over45
## [1] 75 70 65 60 55 50
SCRIPT FOR DISCRETE TIME MODEL:
{r}
# Initial values for population sizes
init_vals <- c(25,75)
# Colors for each initial value/line
colors <- c("lightblue4","red")
# Time sequence for the simulation
times <- seq(0, 50, by = 1) # Discrete time steps
# Parameters for the density-dependent growth model
parms <- c(r=-0.5, K = 100) # r: growth rate, K: carrying capacity
# empty plot, set x and y axes
plot(
NULL,
xlim = c(0, max(times)),
ylim = c(0, 200),
xlab = "Time",
ylab = "Population Size (N)",
main = "Density-Dependent Growth with Different Initial Populations"
)
# Loop through the initial values
for (i in seq_along(init_vals)) {
# Initialize population size for this initial value
N <- numeric(length(times))
N[1] <- init_vals[i] # Set the initial population size
# Simulate the discrete model over time
for (t in 2:length(times)) {
N[t] <- N[t - 1] * exp(parms["r"] * (1 - (N[t - 1] / parms["K"])))
}
# Plot the results with distinct colors
lines(times, N, col = colors[i], lwd = 2)
}
# Add a legend to the plot
legend(
x = 35, y = 190,
legend = init_vals,
col = colors,
lty = 1,
title = "Initial Population (N)"
)
parameters: r=-0.5, initial population values = 25;75, K=100
To have a stable equilibrium at n*= 0, r must be a negative value ( r<0).
parameters: r=0.5, initial population values = 50;80, K=100
To have a stable equilibrium at n*= K approached without oscillations, r must be a small positive value ( 0>=r<1).
parameters: r=2, initial population values = 50;150, K=100
To have a stable equilibrium at n*= K approached via decaying oscillations, r must be a larger positive value. However, to maintain decay r must be less than or equal to 2 (1>x>=2).
parameters: r=4, initial population values = 50;200, K=100
To have a stable equilibrium at n*= K approached via persistent,regular oscillations, r must be a larger positive value greater than 2 and smaller than 3(2<x<3)
parameters: r=4, initial population values = 50;51, K=100
To demonstrate chaos, r must be a positive value greater than 3 (x>3).
\[ \frac{dp}{dt} = mp(1 - p) - ep \]
a. The state variables are:
The parameters are:
m : constant; relates to colonization rate
ep : rate of extinction
b. The colonization rate depends on both the number of occupied patches (p) and patches that are empty, or the term (1-p). When describing colonization we must consider what is available to colonize and what is already occupied.
c. simplifying assumptions:
There is a constant rate of colonization and extinction.
There is no interaction between the patches and they only have colonization and extinction rates as a commonality.
No spatial variation !
In my opinion the model is simplified as in spatial population ecology we must always consider variability. Extinction and colonization rates will not be the same for every patch due to external factors.
d.
\[ \frac{dp}{dt} = mp(1 - p) - ep \]
\[ \frac{dp}{dt} = 0 \]
\[ mp(1 - p) - ep = 0 \]
\[ mp-mp^2-ep=0 \]
\[ -mp^2 + (m - e)p = 0 \]
\[ p(-mp + (m - e)) = 0 \]
\[ p = 1 - \frac{e}{m} \]
The equilibrium points are:
\(p = 0\): all patches are empty!
\(p = 1 - \frac{e}{m}\): A positive equilibrium where a fraction of patches are occupied.
p>0 exists if m>e. P cannot exist otherwise; if e was greater than m, it would signify that extinction surpassed colonization, so all patches would be empty, otherwise stated as p=0. Stability can only be possible if colonization compensates extinction + growth.
e.
# Parameters
m <- 3 # Colonization rate
e <- 2 # Extinction rate
# Range of p values
p_values <- seq(0, 1, by = 0.01)
# Compute dp/dt for each p
dp_dt <- m * p_values * (1 - p_values) - e * p_values
# Plot dp/dt vs p
plot(
p_values, dp_dt,
type = "l",
col = "blue",
lwd = 2,
xlab = "Fraction of Occupied Patches (p)",
ylab = "Rate of Change (dp/dt)",
main = "Rate of Change vs Fraction Occupied"
)
abline(h = 0, col = "red", lty = 2) # Add a horizontal line at dp/dt = 0
# Mark equilibrium points
equilibria <- c(0, 1 - e / m) # p = 0 and p = 1 - e/m
points(equilibria, rep(0, length(equilibria)), col = "red", pch = 19)
text(equilibria, rep(0, length(equilibria)) + 0.02, labels = c("p = 0", paste("p =", round(1 - e / m, 2))), col = "red")
predicted stability: p=0.33 will be stable equilibrium.
f.
library(deSolve)
## Warning: package 'deSolve' was built under R version 4.3.3
# Define the Levins model
levinsModel <- function(time, state, parms) {
with(as.list(c(state, parms)), {
dp <- m * p * (1 - p) - e * p # Levins metapopulation model equation
list(c(dp))
})
}
# Parameters
parms <- c(m = 3, e = 2) # Colonization and extinction rates
times <- seq(0, 10, by = 0.1) # Time steps
init_vals <- c(0.1, 0.5, 0.9) # Initial values of p
# Set up the plot
colors <- c("red", "blue", "green")
plot(
NULL,
xlim = c(0, max(times)),
ylim = c(0, 1),
xlab = "Time",
ylab = "Fraction of Occupied Patches (p)",
main = "Levins Model Dynamics for Different Initial Conditions"
)
# Simulate for each initial condition
for (i in seq_along(init_vals)) {
init <- c(p = init_vals[i])
out <- lsoda(init, times, levinsModel, parms)
lines(out[, 1], out[, 2], col = colors[i], lwd = 2)
}
# Add legend
legend(
"right",
legend = paste("p0 =", init_vals),
col = colors,
lty = 1,
title = "Initial Conditions"
)
Like in (e), at p=0.33 the graph is at equilibrium.
g.
# Define new parameter sets
param_sets <- list(
list(m = 4, e = 2), # m > e
list(m = 2, e = 2), # m = e
list(m = 1, e = 2) # m < e
)
# Set up the plot
plot(
NULL,
xlim = c(0, 1),
ylim = c(-0.5, 1),
xlab = "Fraction of Occupied Patches (p)",
ylab = "Rate of Change (dp/dt)",
main = "Rate of Change vs Fraction Occupied for Different Parameters"
)
# Compute and plot dp/dt for each parameter set
colors <- c("blue", "green", "red")
for (i in seq_along(param_sets)) {
m <- param_sets[[i]]$m
e <- param_sets[[i]]$e
dp_dt <- m * p_values * (1 - p_values) - e * p_values
lines(p_values, dp_dt, col = colors[i], lwd = 2)
}
# Add a legend
legend(
"topright",
legend = c("m > e", "m = e", "m < e"),
col = colors,
lty = 1,
lwd = 2,
title = "Parameter Cases"
)
Part A: State Variables: P1, P2 (predator populations) and N (prey population)
parameters:
r: prey birth rate
a1: attack rate of Predator 1 on prey
a2:attack rate of Predator 2 on prey
d1,d2: death rates of predators
p: killing rate of predator 1 on predator 2
Part C: State Variables: P1, P2, P3 (plant density in patch #1,2,3 sequentially)
parameters:
g: seed germination rate
m: mortality rate
f: seed production rate
w: wind-dispersal probability (movement from one patch to another)
ODE Model:
\[ \begin{align*}\frac{dP_1}{dt} &= (1 - w)fgP_1 + wfg(P_2 + P_3) - mP_1 \\\frac{dP_2}{dt} &= (1 - w)fgP_2 + wfg(P_1 + P_3) - mP_2 \\\frac{dP_3}{dt} &= (1 - w)fgP_3 + wfg(P_1 + P_2) - mP_3\end{align*} \]
Explanation:
mPᵢ = death of plants in the selected patch (P1,P2,P3)
wfg(Pⱼ + Pₖ) = seeds dispersed into this patch from the other two patches (j ≠ i, k ≠ i).
(1 - w)fgPᵢ represents the seeds that stay in the same patch and successfully germinate.
I hope this makes sense!