Imagine you're throwing a party in a small studio apartment. At first, it's great—people arrive quickly, the energy builds, and the room fills up fast. But as the space hits its limit, things change. It becomes harder for new guests to squeeze through the door, and some people start leaving because it's simply too crowded to breathe. The "growth" of your party doesn't just stop abruptly; it slows down as you approach the maximum capacity of the room.

In ecology, we call this the Logistic Growth model, and it's the foundation for most population modeling. Let's map that party scenario directly onto the math we use in R:

• The arrival rate is your intrinsic growth rate (r). This is how fast the population grows when resources are unlimited.
• The apartment size is the carrying capacity (K). This is the maximum population the environment can actually sustain.
• The "crowdedness" factor is the term (1 - N/K). When the population (N) is small, this value is close to 1 (full speed ahead). As N approaches K, this value drops toward 0, killing the growth rate.

Translating biology into differential equations

To do this right in R, we don't just want a static calculation; we want to see the population change over time. I usually reach for the deSolve package. It's the industry standard for handling ordinary differential equations (ODEs). If you haven't installed it yet, you'll need install.packages("deSolve") first.

Here is how I structure a basic logistic growth function. Notice how I keep the parameters in a separate vector—this makes it much easier to tweak the "apartment size" later without digging through the logic.

library(deSolve)

# The model function
logistic_growth <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    # N is our current population
    # r is growth rate, K is carrying capacity
    dN <- r * N * (1 - N / K)
    return(list(c(dN)))
  })
}

# Initial state: we start with 10 individuals
state <- c(N = 10)

# Parameters: 10% growth rate, max capacity of 1000
parameters <- c(r = 0.1, K = 1000)

# Time sequence: from day 0 to day 100
times <- seq(0, 100, by = 1)

# Solving the equation
out <- ode(y = state, times = times, func = logistic_growth, parms = parameters)
as.data.frame(out)

Letting the population run wild

The ode function is doing the heavy lifting here. It's essentially taking tiny steps through time, calculating the slope (the growth rate) at that exact moment, and moving the population forward. I've found that the biggest mistake people make is forgetting that state and parameters must be named vectors. If you just pass raw numbers, with(as.list(...)) will fail because it won't know which number is r and which is K.

If you plot the resulting data frame, you'll see that classic "S-curve." It starts exponential, then bends as it feels the pressure of the carrying capacity, and finally flattens out. It's a beautiful piece of math because it mirrors reality: nothing grows forever.

Handling environmental shocks

Real ecology is messier than a clean S-curve. What happens if a disease hits or a forest fire clears out half the population? In R, we handle this by breaking the simulation into chunks. You run the model for 50 days, manually change the state vector (e.g., state["N"] <- state["N"] / 2), and then use that new state as the starting point for the next 50 days.

📋 Practical Task

Exercise: Simulating a Rabbit Population Crash

You are modeling a rabbit population on a small island. The population starts at 50 rabbits, has a growth rate of 0.2, and a carrying capacity of 500. However, at day 30, a predator is introduced that instantly kills 60% of the population.

Write a script that:

• Uses the deSolve package to simulate the first 30 days of growth.
• Calculates the population at day 30 and reduces it by 60%.
• Uses that reduced population as the starting state to simulate another 30 days (from day 31 to 60).
• Combines both results into a single data frame to show the "crash" and the subsequent recovery.