R Script for a Basic Metapopulation Model

Explanation: Parameters:

n_patches: The number of habitat patches in the landscape. time_steps: The number of time steps over which to simulate the model. colonization_rate: The probability that an empty patch will be colonized at each time step. extinction_rate: The probability that an occupied patch will go extinct at each time step. Initial Setup:

The model starts with a random initial occupancy of patches. Metapopulation Dynamics:

For each time step, the script iterates over all patches. If a patch is occupied, there is a chance it will go extinct. If it is empty, there is a chance it will be colonized.

# Basic Metapopulation Model in R

# Parameters
n_patches <- 10          # Number of habitat patches
time_steps <- 50         # Number of time steps to simulate
colonization_rate <- 0.1 # Probability of colonization per time step
extinction_rate <- 0.05  # Probability of extinction per time step

# Initialize a vector to track the occupancy of patches (1 = occupied, 0 = empty)
occupancy <- numeric(n_patches)
occupancy[] <- sample(c(0, 1), n_patches, replace = TRUE) # Random initial occupancy

# Store results over time
occupancy_over_time <- matrix(0, nrow = time_steps, ncol = n_patches)
occupancy_over_time[1, ] <- occupancy

# Simulate metapopulation dynamics
for (t in 2:time_steps) {
  for (i in 1:n_patches) {
    if (occupancy[i] == 1) {
      # Patch is currently occupied, it may go extinct
      if (runif(1) < extinction_rate) {
        occupancy[i] <- 0
      }
    } else {
      # Patch is currently empty, it may be colonized
      if (runif(1) < colonization_rate) {
        occupancy[i] <- 1
      }
    }
  }
  # Record occupancy at this time step
  occupancy_over_time[t, ] <- occupancy
}

# Plot the occupancy over time
plot(1:time_steps, rowSums(occupancy_over_time), type = "l", col = "blue",
     xlab = "Time", ylab = "Number of Occupied Patches",
     main = "Metapopulation Dynamics Over Time")

To find the proportion of occupied patches at equilibrium, we can use a simple metapopulation model. At equilibrium, the rate of colonization equals the rate of extinction, and the proportion of occupied patches can be calculated directly.

Derivation: Let’s denote:

𝑝 p as the proportion of occupied patches, 𝑐 c as the colonization rate (probability of colonization per empty patch per unit time), 𝑒 e as the extinction rate (probability of extinction per occupied patch per unit time). At equilibrium, the proportion of patches colonized should equal the proportion of patches going extinct:

R Script to Compute and Simulate the Equilibrium

# Parameters
colonization_rate <- 0.1  # Probability of colonization per empty patch per unit time
extinction_rate <- 0.05  # Probability of extinction per occupied patch per unit time

# Calculate equilibrium proportion of occupied patches
equilibrium_proportion <- 1 - (extinction_rate / colonization_rate)

# Display the equilibrium proportion
equilibrium_proportion
## [1] 0.5
# Simulation to verify the equilibrium

# Parameters for simulation
n_patches <- 100          # Number of habitat patches
time_steps <- 200         # Number of time steps to simulate

# Initialize a vector to track the occupancy of patches (1 = occupied, 0 = empty)
occupancy <- numeric(n_patches)
occupancy[] <- sample(c(0, 1), n_patches, replace = TRUE) # Random initial occupancy

# Store results over time
occupancy_over_time <- numeric(time_steps)

# Simulate metapopulation dynamics
for (t in 1:time_steps) {
  for (i in 1:n_patches) {
    if (occupancy[i] == 1) {
      # Patch is currently occupied, it may go extinct
      if (runif(1) < extinction_rate) {
        occupancy[i] <- 0
      }
    } else {
      # Patch is currently empty, it may be colonized
      if (runif(1) < colonization_rate) {
        occupancy[i] <- 1
      }
    }
  }
  # Record the proportion of occupied patches at this time step
  occupancy_over_time[t] <- mean(occupancy)
}

# Plot the proportion of occupied patches over time
plot(1:time_steps, occupancy_over_time, type = "l", col = "blue",
     xlab = "Time", ylab = "Proportion of Occupied Patches",
     main = "Metapopulation Proportion of Occupied Patches Over Time")

# Add a horizontal line indicating the calculated equilibrium proportion
abline(h = equilibrium_proportion, col = "red", lty = 2)
legend("topright", legend = c("Simulation", "Equilibrium Proportion"),
       col = c("blue", "red"), lty = c(1, 2))

We can increase the number of occupied habitat patches by providing corridors between patches to increase the rate of colonization (c).

We can also increase the number of occupied patches by decreasing rates of extinction (e).

# Parameters
colonization_rate <- 0.7  # Probability of colonization per empty patch per unit time
extinction_rate <- 0.01   # Probability of extinction per occupied patch per unit time

# Calculate equilibrium proportion of occupied patches
equilibrium_proportion <- 1 - (extinction_rate / colonization_rate)

# Display the equilibrium proportion
equilibrium_proportion
## [1] 0.9857143
# Simulation to verify the equilibrium

# Parameters for simulation
n_patches <- 100          # Number of habitat patches
time_steps <- 200         # Number of time steps to simulate

# Initialize a vector to track the occupancy of patches (1 = occupied, 0 = empty)
occupancy <- numeric(n_patches)
occupancy[] <- sample(c(0, 1), n_patches, replace = TRUE) # Random initial occupancy

# Store results over time
occupancy_over_time <- numeric(time_steps)

# Simulate metapopulation dynamics
for (t in 1:time_steps) {
  for (i in 1:n_patches) {
    if (occupancy[i] == 1) {
      # Patch is currently occupied, it may go extinct
      if (runif(1) < extinction_rate) {
        occupancy[i] <- 0
      }
    } else {
      # Patch is currently empty, it may be colonized
      if (runif(1) < colonization_rate) {
        occupancy[i] <- 1
      }
    }
  }
  # Record the proportion of occupied patches at this time step
  occupancy_over_time[t] <- mean(occupancy)
}

# Plot the proportion of occupied patches over time
plot(1:time_steps, occupancy_over_time, type = "l", col = "blue",
     xlab = "Time", ylab = "Proportion of Occupied Patches",
     main = "Metapopulation Proportion of Occupied Patches Over Time")

# Add a horizontal line indicating the calculated equilibrium proportion
abline(h = equilibrium_proportion, col = "red", lty = 2)
legend("topright", legend = c("Simulation", "Equilibrium Proportion"),
       col = c("blue", "red"), lty = c(1, 2))