Exploring the Levins Metapopulation Model

Introduction

Many species live in landscapes made up of discrete habitat patches. Local populations may go extinct due to stochastic or deterministic events, but the species can persist at the regional level if colonisation of empty patches occurs frequently enough.

In 1969, Richard Levins proposed a simple and elegant model to capture these dynamics. His metapopulation model assumes that there are an infinite number of patches and that each patch is either:

Occupied, with probability \(p\)
Empty, with probability \(1 - p\)

The change in the proportion of occupied patches over time is given by the differential equation:

\[ \frac{dp}{dt} = c \cdot p \cdot (1 - p) - e \cdot p \]

Where:

\(c\) is the colonisation rate (probability an empty patch becomes occupied),
\(e\) is the extinction rate (probability an occupied patch becomes empty),
\(p\) is the proportion of occupied patches.

The model reaches equilibrium when:

\[ p^* = 1 - \frac{e}{c} \]

This exercise explores this model using R, helping you understand:

The effect of colonisation and extinction rates.
The role of initial occupancy.
The sensitivity of equilibrium to model parameters.

The Levins Model in R

We will use the following R function to simulate the model. Copy and paste it into R to make it available.

levins_model <- function(p0, c, e, t_max, dt = 0.01) {
  times <- seq(0, t_max, by = dt)
  p <- numeric(length(times))
  p[1] <- p0
  
  for (i in 2:length(times)) {
    dp <- c * p[i - 1] * (1 - p[i - 1]) - e * p[i - 1]
    p[i] <- p[i - 1] + dp * dt
  }
  
  data.frame(time = times, occupancy = p)
}

Exercise 1: Varying Colonisation and Extinction Rates

To simulate the model, choose values for the initial occupancy (p0), colonisation rate (c), extinction rate (e), and total simulation time (t_max). The function will return a time series of patch occupancy. You can then plot this to observe how the metapopulation changes over time.

# Define parameters
p0 <- 0.2     # initial proportion of occupied patches
c <- 0.3      # colonisation rate
e <- 0.2      # extinction rate
t_max <- 50   # total time to simulate

# Run the model
result <- levins_model(p0, c, e, t_max)

# Plot the results
plot(result$time, result$occupancy, type = 'l', col = 'blue', lwd = 2,
     xlab = "Time", ylab = "Proportion of Occupied Patches",
     main = paste("Levins Model: c =", c, ", e =", e))

# Add theoretical equilibrium line
abline(h = 1 - e / c, col = 'red', lty = 2)

This will give you a curve showing how the proportion of occupied patches changes over time, along with a dashed line indicating the expected equilibrium.

When that is working, try running the model using different combinations of colonisation and extinction rates. Predict what will happen to occupancy over time. Then run the model and check your predictions.

For example:

Scenario	Colonisation Rate (c)	Extinction Rate (e)	Expected Equilibrium \(p^* = 1 - \frac{e}{c}\)
A	0.1	0.2	Extinction (c < e)
B	0.3	0.3	0.00 (extinction)
C	0.6	0.2	0.67

Try some others…

Reflection:

Which systems reached equilibrium?
Which system went extinct?
Did the equilibrium match the formula \(1 - \frac{e}{c}\)?
Why do we see these patterns?

Exercise 2: Changing the Initial Occupancy

How does the starting occupancy affect the long-term dynamics?

Use the same code. This time, try setting colonisation rate (c) and extinction rate (e) to be 0.5 and 0.2 respectively. Then run the model with different starting values (‘p0’) (these can vary between just above 0 and 1).

Reflection:

How does the starting point influence the outcome?
Does the final equilibrium change?

Exercise 3: Sensitivity to Colonisation Rate

This exercise focuses on one of the theoretical predictions of the Levins model:

\[ p^* = 1 - \frac{e}{c} \]

You’ll explore how equilibrium occupancy changes as the colonisation rate \(c\) varies, while the extinction rate \(e\) stays fixed.

Set the extinction rate to:

e <- 0.2

Create a sequence of colonisation rates, starting at 0.05 and going to 1:

c_vals <- seq(0.05, 1, by = 0.01)

Use the formula to calculate the theoretical equilibrium for each \(c\). Make sure you handle values where \(c \leq e\) by setting equilibrium to 0 (extinction):

p_eq <- 1 - e / c_vals
p_eq <- ifelse(p_eq < 0, 0 ,p_eq)

Plot your results:

plot(c_vals, p_eq, type = "l", col = "blue", lwd = 2,
     xlab = "Colonisation rate (c)", ylab = "Equilibrium occupancy (p*)",
     main = "Theoretical Equilibrium Occupancy")

abline(v = e, col = "red", lty = 2)  # extinction threshold line

Try different values of e

What to Think About

What might happen to a species if its habitat becomes so fragmented that dispersal between patches is no longer possible? In other words, what happens to occupancy when \(c \leq e\)?
If we restore or connect habitat to make dispersal easier, how much improvement is needed before we start seeing real benefits for population survival? In other words, how quickly does occupancy increase as \(c\) rises, and does this vary along the curve?
Notice that the curve flattens as \(c\) gets large. Based on this, if a species already has good dispersal opportunities, should we keep investing in improving connectivity, or shift our focus elsewhere?

Summary

By now, you should understand:

The basic behaviour of the Levins metapopulation model.
The conditions for metapopulation persistence or extinction.
The relationship between colonisation/extinction rates and equilibrium occupancy.

This simple model offers a foundation for more complex spatial and stochastic models in conservation biology.