📘 Simulating Population Growth in Variable Environments

This tutorial simulates how a population might grow or decline over time when faced with varying environmental conditions — “good years” and “bad years”. It uses matrix population models (MPMs) and introduces the concept of a stochastic growth rate.

We’ll walk through each step of the code and explain what it’s doing in plain language. Don’t worry if some parts of the code seem confusing — just focus on the big ideas, try running the code, and ask your instructors for help if you get stuck!

🧬 Setup: Define Matrices for Good and Bad Years

We begin by defining two different projection matrices, representing a good year and a bad year for the population. Each matrix describes how individuals move between life stages and how many new individuals are produced.

# Set seed so results are the same each time we run
set.seed(42)

# Good year matrix
A_good <- matrix(c(
  0.0, 1.2, 1.5,   # Fecundity (new individuals from stage 2 and 3)
  0.5, 0.3, 0.0,   # Survival and growth between stages
  0.0, 0.6, 0.7
), nrow = 3, byrow = TRUE)

# Bad year matrix
A_bad <- matrix(c(
  0.0, 0.8, 0.75,
  0.3, 0.2, 0.0,
  0.0, 0.4, 0.5
), nrow = 3, byrow = TRUE)

Each matrix describes transitions between three life stages (say, juveniles, subadults, adults).

🧮 Initial Population and Simulation Parameters

Now we define the starting population and simulation settings.

# Starting population vector: 50 juveniles, 30 subadults, 20 adults
n0 <- c(50, 30, 20)

# Set parameters
timesteps <- 100     # Number of years to simulate
n_sim <- 1000        # Number of repeated simulations
p_good <- 0.5        # Probability that a given year is a good year

📦 Create Storage for Results

We’ll track total population size and stochastic growth rates (lambda_s) for each simulation.

# Vector to store growth rate estimates
lambda_s <- numeric(n_sim)

# Matrix to store population size over time for each simulation
trajectory_data <- matrix(NA, nrow = timesteps + 1, ncol = n_sim)

🔁 Run the Simulations

Here’s the heart of the analysis! We simulate population dynamics over time, choosing randomly between good and bad years.

for (i in 1:n_sim) {
  pop <- n0  # Start with initial population
  total_N <- numeric(timesteps + 1)
  total_N[1] <- sum(pop)
  
  for (t in 1:timesteps) {
    # Randomly choose either a good year or a bad year
    A_t <- if (runif(1) < p_good) A_good else A_bad
    pop <- A_t %*% pop  # Project population forward one year
    total_N[t + 1] <- sum(pop)  # Store total population size
  }
  
  # Estimate stochastic growth rate (slope of log population size over time)
  time <- 0:timesteps
  lambda_s[i] <- coef(lm(log(total_N) ~ time))[2]
  
  # Save trajectory
  trajectory_data[, i] <- total_N
}

💡 Note: Don’t worry if you don’t understand all of this code — try copying and running it, and focus on understanding what it’s modeling!

📈 Visualize Population Trajectories

Let’s plot the population trajectories on a log scale so we can compare them more easily.

# Log-transform the data
trajectory_data <- log(trajectory_data)

# Plot all trajectories
matplot(0:timesteps, trajectory_data, type = "l", lty = 1,
        col = adjustcolor("darkblue", alpha.f = 0.3),
        xlab = "Time (years)", ylab = "Total Population Size (log)",
        main = "Variation in Stochastic Trajectories (1000 Runs)")

# Add a red line showing the mean trajectory
lines(0:timesteps, rowMeans(trajectory_data, na.rm = TRUE), 
      col = "red", lwd = 2)

legend("topright", legend = c("Individual runs", "Mean"), 
       col = c("darkblue", "red"), lty = c(1, 1), lwd = c(1, 2), bty = "n")

📊 Estimate and Plot Stochastic Growth Rate (λs)

Now we calculate and visualize the stochastic growth rate (λs) — how fast the population grows on average across all simulations.

# Convert from log scale to multiplicative growth rate
lambda_s_exp <- exp(lambda_s)

# Histogram of growth rates
hist(lambda_s_exp, breaks = 30, col = "skyblue", border = "white",
     main = "Distribution of Stochastic Growth Rate (λs)",
     xlab = "λs (per time step)")

# Add mean and 95% confidence interval
mean_lambda_s <- mean(lambda_s_exp)
ci_lambda_s <- quantile(lambda_s_exp, probs = c(0.025, 0.975))

abline(v = mean_lambda_s, col = "red", lwd = 2)
abline(v = ci_lambda_s, col = "darkgreen", lty = 2)

legend("topright", legend = c("Mean", "95% CI"), 
       col = c("red", "darkgreen"), lty = c(1, 2), lwd = 2)

# Print results
cat("Estimated stochastic growth rate (λs):\n")

## Estimated stochastic growth rate (λs):

cat("Mean λs:", round(mean_lambda_s, 4), "\n")

## Mean λs: 1.0418

cat("95% CI:", round(ci_lambda_s[1], 4), "-", round(ci_lambda_s[2], 4), "\n")

## 95% CI: 0.992 - 1.0944

🧯 Risk of Decline or Extinction

We can now assess the risk of decline. If λs < 1, the population is shrinking.

# How many runs are showing decline?
table(lambda_s_exp < 1)

## 
## FALSE  TRUE 
##   951    49

🧪 Compare to Deterministic Growth Rates

Let’s compare the stochastic growth rate to the deterministic growth rates of the good and bad years (i.e., assuming only good or only bad years).

# Compare with deterministic growth rate from each matrix
(lambda_bad <- popdemo::eigs(A_bad, what = "lambda"))

## [1] 0.8255641

(lambda_good <- popdemo::eigs(A_good, what = "lambda"))

## [1] 1.314143

# Mean of the deterministic rates
mean(c(lambda_bad, lambda_good))

## [1] 1.069853

📌 You should notice that the stochastic growth rate is lower than the average of the good and bad year deterministic rates. This is a key insight: variation reduces growth!

🔍 Things to Explore

Try changing these parameters and see what happens:

What if bad years become more frequent? (p_good < 0.5)
What if the difference between good and bad years increases (e.g. make good better, bad worse)?
What happens to the shape of the histogram and the mean growth rate?

✅ Summary

This simulation models population dynamics under environmental stochasticity.
Even if good years occur half the time, variability can reduce long-term growth.
Stochastic growth rate (λs) helps us understand the long-term viability of a population under variable conditions.

Stochastic Matrix Models: Part 1

Owen Jones