---
title: "Logistic Population Growth"
subtitle: "When space starts to matter: from the exponential to the logistic equation"
author: "Population Ecology Laboratory — Department of Biology"
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 3
    number_sections: true
    theme: flatly
    highlight: tango
    code_folding: show
    df_print: paged
  pdf_document:
    toc: true
    number_sections: true
    latex_engine: xelatex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 7,
  fig.height = 4.5,
  fig.align = "center"
)
library(ggplot2)
library(deSolve)
theme_set(theme_minimal(base_size = 12))
```

# Where we left off

In the previous laboratory we arrived at the exponential growth formula:

$$
\frac{dN}{dt} = r N \quad\Longrightarrow\quad N(t) = N_0\, e^{rt}
$$

With $N_0 = 100$ water lilies and $r = 0.1$/day, let's see what the model predicts over the long run:

```{r exp-projection}
N0 <- 100
r  <- 0.1

cat("In 100 days:", format(N0 * exp(r * 100),  big.mark = ","), "water lilies\n")
cat("In 200 days:", format(N0 * exp(r * 200),  big.mark = ","), "water lilies\n")
cat("In 365 days:", format(N0 * exp(r * 365),  big.mark = ","), "water lilies\n")
```

**It is absurd.** In less than a year, our small lake is supposed to have more water lilies than there are cells in the human body.

The problem is that the exponential model **ignores space, light, nutrients and competition**. In this laboratory we will fix this shortcoming.

---

# The idea of the carrying capacity $K$

Every real environment has a **limit** to the number of individuals it can sustain. We call this limit the **carrying capacity** and denote it by $K$.

For our lake, let's assume:

```{r set-K}
K <- 10000   # maximum number of water lilies the lake can hold
```

> *A note on the choice:* The value $K = 10{,}000$ is illustrative and was chosen for clarity in the plots. In a real study $K$ is estimated from data or from the physical features of the habitat (lake surface area, average water lily size, available light, etc.).

## $K$ is not just a "ceiling"

The naive idea would be: *"we let the population grow exponentially up to $K$, and then we stop it."* But that is not what happens in nature — growth **slows down gradually** as we approach $K$, *long before* we reach it.

Why?

---

# A mechanism: why $K$ reduces the rate

Imagine our lake with the water lilies. When the population is small, each water lily has plenty of space around it and can reproduce freely. But as the population grows, **the water lilies in the interior of the cluster are surrounded by other water lilies** and struggle to expand. Only the water lilies at the *periphery* have space. The fuller the lake, the smaller the fraction of water lilies that can actually reproduce.

This is **one example** of a mechanism. In reality several reasons can act together:

- **Space / light:** the water lilies in the interior are shaded.
- **Nutrients:** competition for root nutrients intensifies.
- **Pathogens:** in dense clusters diseases spread quickly.
- **Toxic by-products:** they accumulate as the population becomes denser.

All of these share the same **qualitative** feature: the **per capita reproduction rate decreases** as the population approaches $K$.

> So $K$ is not just a limit we bump into. $K$ **modifies the rate of change itself**.

---

# From physical intuition to mathematical form

## What we want from the "realized" per capita rate

We define the **realized** per capita rate of increase as a function $r_{\text{eff}}(N)$ of the current population. Two natural requirements:

1. When $N$ is **small**, there is plenty of space → $r_{\text{eff}}(N) \to r$ *(the intrinsic rate)*.
2. When $N \to K$, space runs out → $r_{\text{eff}}(N) \to 0$.

## The simplest function that satisfies these: linear

The simplest function that captures the two extremes is a **straight line** that falls from $r$ (at $N=0$) to $0$ (at $N=K$):

$$
\boxed{\;r_{\text{eff}}(N) \;=\; r\left(1 - \frac{N}{K}\right)\;}
$$

> **A note on the model choice (important):** The linear form is a *modeling assumption*, not a logical necessity. It is the simplest function that captures the correct qualitative behavior. Other possible forms (we will see them at the end):
>
> - **θ-logistic:** $r_{\text{eff}}(N) = r\bigl(1 - (N/K)^{\theta}\bigr)$
> - **Gompertz:** $r_{\text{eff}}(N) = r\,\ln(K/N)$
>
> All give similar sigmoidal behavior, but with different geometry.

```{r per-capita-comparison}
N_grid <- seq(0, K, length.out = 200)
df_pc <- data.frame(
  N        = rep(N_grid, 2),
  rate     = c(rep(r, length(N_grid)),                # constant r (exponential)
               r * (1 - N_grid/K)),                   # linearly decreasing (logistic)
  Model    = factor(rep(c("Exponential:  r constant",
                          "Logistic: r·(1-N/K)"),
                        each = length(N_grid)))
)

ggplot(df_pc, aes(x = N, y = rate, colour = Model)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(
    title = "Per capita growth rate as a function of population",
    subtitle = "The essential difference between the exponential and logistic models",
    x = "Population N",
    y = expression(r[eff](N))
  ) +
  theme(legend.position = "bottom")
```

---

# The logistic equation

The total rate of change of the population is:

$$
\frac{dN}{dt} \;=\; r_{\text{eff}}(N) \cdot N \;=\; r\,N\left(1 - \frac{N}{K}\right)
$$

This is the **logistic equation** (Verhulst, 1838).

$$
\boxed{\;\frac{dN}{dt} \;=\; r\,N\left(1 - \frac{N}{K}\right)\;}
$$

---

# Qualitative analysis *before* we solve

One of the most useful tools in dynamical systems is to examine the behavior **before** solving analytically. Very often the solution does not even exist analytically — but the qualitative picture is enough for us.

## Equilibrium points

Where $\dfrac{dN}{dt} = 0$ the population stays stationary. We solve:

$$
r\,N\left(1 - \frac{N}{K}\right) = 0 \quad\Longrightarrow\quad N^* = 0 \;\text{ or }\; N^* = K
$$

- $N^* = 0$: empty lake (population extinct).
- $N^* = K$: full lake (saturation).

## Stability of the equilibrium points

- If we start **close** to $0$ (e.g. $N=1$): $\frac{dN}{dt} > 0$, the population *grows* — it moves away from 0. So $N^*=0$ is **unstable**.
- If we start **close** to $K$ (e.g. $N=0.9K$): $\frac{dN}{dt} > 0$, the population *approaches* $K$.
- If we start **above** $K$ (e.g. $N=1.5K$): $\frac{dN}{dt} < 0$, the population *decreases* toward $K$.

So $N^*=K$ is **stable** — it is an *attractor* of the system.

## Phase portrait: $\dfrac{dN}{dt}$ against $N$

```{r phase-portrait}
N_grid <- seq(-500, K * 1.4, length.out = 300)
dN_dt  <- r * N_grid * (1 - N_grid/K)
df_phase <- data.frame(N = N_grid, dN = dN_dt)

# Equilibrium points
eq_unstable <- data.frame(N = 0, dN = 0)
eq_stable   <- data.frame(N = K, dN = 0)

# Flow arrows on the N axis: the sign of (xend - x) gives the direction of motion
arrows_df <- data.frame(
  x    = c(K/4,        3*K/4,        1.2*K),
  xend = c(K/4 + 800,  3*K/4 + 800,  1.2*K - 800),
  y    = 0, yend = 0
)

ggplot(df_phase, aes(x = N, y = dN)) +
  geom_hline(yintercept = 0, colour = "grey50") +
  geom_line(linewidth = 1, colour = "darkgreen") +
  geom_point(data = eq_unstable, aes(x = N, y = dN),
             shape = 21, fill = "white", colour = "darkgreen",
             size = 4, stroke = 1.2) +
  geom_point(data = eq_stable, aes(x = N, y = dN),
             colour = "darkgreen", size = 4) +
  geom_segment(data = arrows_df,
               aes(x = x, xend = xend, y = y, yend = yend),
               arrow = arrow(length = unit(0.25, "cm")),
               colour = "tomato", linewidth = 0.9,
               inherit.aes = FALSE) +
  annotate("text", x = 0,    y = -150, label = "N* = 0\n(unstable)",  vjust = 1) +
  annotate("text", x = K,    y = -150, label = "N* = K\n(stable)", vjust = 1) +
  annotate("text", x = K/2,  y = max(dN_dt) + 30,
           label = "maximum rate\nat N = K/2", vjust = 0) +
  labs(
    title = "Phase portrait of the logistic equation",
    subtitle = "The arrows show the direction of change of N",
    x = "N", y = "dN/dt = r"
  )
```

The parabola represents **the whole dynamics of the system** in a single picture:

- Where the curve is **above** the axis → $\frac{dN}{dt} > 0$ → the population increases.
- Where the curve is **below** the axis → the population decreases.
- At the *roots* we have equilibrium points.
- The **peak** of the parabola (at $N = K/2$) is the maximum growth rate — that is where the population grows fastest. This is also the **inflection point** of the time trajectory.

---

# Analytical solution *(optional for the mathematically inclined)*

For those who want to see the solution: we start with separation of variables:

$$
\frac{dN}{N\,(1 - N/K)} = r\,dt
\quad\Longleftrightarrow\quad
\frac{K\,dN}{N\,(K - N)} = r\,dt
$$

**Partial fractions:**

$$
\frac{K}{N(K-N)} = \frac{1}{N} + \frac{1}{K-N}
$$

So:

$$
\int\!\left(\frac{1}{N} + \frac{1}{K-N}\right) dN = \int r\,dt
\;\Longleftrightarrow\;
\ln N - \ln(K-N) = rt + C
$$

$$
\ln\!\left(\frac{N}{K-N}\right) = rt + C
\;\Longrightarrow\;
\frac{N}{K-N} = A\,e^{rt}, \quad A = \frac{N_0}{K - N_0}
$$

Solving for $N$ we arrive at the **closed-form logistic solution**:

$$
\boxed{\;N(t) \;=\; \dfrac{K}{1 + \left(\dfrac{K - N_0}{N_0}\right) e^{-rt}}\;}
$$

**Sanity checks:**

- At $t = 0$: $N(0) = K \big/ \left(1 + \tfrac{K-N_0}{N_0}\right) = K \cdot \tfrac{N_0}{K} = N_0$. ✓
- At $t \to \infty$: $e^{-rt} \to 0$, so $N \to K$. ✓
- If $N_0 \ll K$, for small $t$ the term $\tfrac{K-N_0}{N_0} \approx \tfrac{K}{N_0}$ is huge, and the behavior is almost exponential (why?). ✓

---

# R: numerical solution with `deSolve`

Now that we know what to expect, let's simulate:

```{r logistic-ode}
logistic_model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N * (1 - N / K)
    list(dN)
  })
}

parameters <- c(r = 0.1, K = 10000)
state      <- c(N = 100)
times      <- seq(0, 200, by = 0.5)

out_log <- as.data.frame(
  ode(y = state, times = times, func = logistic_model, parms = parameters)
)

head(out_log)
```

## Comparison of the analytical and numerical solutions

```{r compare-analytical-numerical}
N0 <- 100
K  <- 10000
r  <- 0.1

logistic_analytical <- function(t, N0, r, K) {
  K / (1 + ((K - N0) / N0) * exp(-r * t))
}

df_log <- data.frame(
  t          = times,
  Analytical = logistic_analytical(times, N0, r, K),
  Numerical  = out_log$N
)

ggplot(df_log, aes(x = t)) +
  geom_line(aes(y = Analytical),  linewidth = 1.2, colour = "steelblue") +
  geom_point(aes(y = Numerical),
             data = df_log[seq(1, nrow(df_log), by = 25), ],
             colour = "tomato", size = 1.8) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95, y = K * 1.04,
           label = "K = 10,000", colour = "grey40", hjust = 1) +
  labs(
    title = "Logistic growth: analytical (line) vs numerical (dots)",
    x = "Time (days)",
    y = "Population N(t)"
  )

cat("Maximum error:", format(max(abs(df_log$Analytical - df_log$Numerical)),
                             scientific = TRUE), "\n")
```

---

# The sigmoidal curve and the three stages

Look carefully at the time trajectory. It has three distinct stages:

```{r three-stages}
stages <- data.frame(
  t = times,
  N = logistic_analytical(times, N0, r, K)
)

# Inflection point: N = K/2
t_inflection <- (1/r) * log((K - N0) / N0)

ggplot(stages, aes(x = t, y = N)) +
  geom_line(linewidth = 1.1, colour = "steelblue") +
  geom_hline(yintercept = K,   linetype = "dashed", colour = "grey50") +
  geom_hline(yintercept = K/2, linetype = "dotted", colour = "grey50") +
  geom_vline(xintercept = t_inflection, linetype = "dotted", colour = "tomato") +
  annotate("point", x = t_inflection, y = K/2, size = 3, colour = "tomato") +
  annotate("text",  x = t_inflection + 5, y = K/2,
           label = "Inflection point (N=K/2)\nmaximum rate", hjust = 0) +
  annotate("text", x =  20, y = 1500, label = "1. Almost exponential",     hjust = 0) +
  annotate("text", x =  60, y = 5500, label = "2. Transitional phase\n(acceleration → deceleration)", hjust = 0) +
  annotate("text", x = 130, y = 9300, label = "3. Saturation toward K",       hjust = 0) +
  labs(
    title = "The sigmoidal (S-shaped) logistic curve",
    x = "Time (days)",
    y = "N(t)"
  )

cat("Inflection point at t =", round(t_inflection, 1), "days (where N = K/2 =", K/2, ")\n")
```

1. **Almost exponential phase** (small $N$): the term $1 - N/K$ is close to 1, so $\frac{dN}{dt} \approx rN$.
2. **Inflection point** ($N = K/2$): here the rate $\frac{dN}{dt}$ is maximized.
3. **Saturation phase** ($N \to K$): the term $1 - N/K \to 0$, so the rate decays exponentially toward zero.

---

# Comparison: exponential vs logistic

With the same $r$, the same $N_0$, but a different $K$:

```{r exp-vs-log}
df_compare <- data.frame(
  t           = times,
  Exponential = N0 * exp(r * times),
  Logistic    = logistic_analytical(times, N0, r, K)
)

# Long format for ggplot
df_long <- data.frame(
  t       = rep(times, 2),
  N       = c(df_compare$Exponential, df_compare$Logistic),
  Model   = factor(rep(c("Exponential", "Logistic"), each = length(times)))
)

ggplot(df_long, aes(x = t, y = N, colour = Model)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  scale_y_log10() +
  annotate("text", x = max(times) * 0.95, y = K * 1.3,
           label = "K = 10,000", colour = "grey40", hjust = 1) +
  labs(
    title = "Exponential vs logistic growth (logarithmic axis)",
    subtitle = "In the first days they are almost identical — then they diverge dramatically",
    x = "Time (days)",
    y = "N(t) [logarithmic scale]"
  ) +
  theme(legend.position = "bottom")
```

In the first $\approx 30$ days the two curves are **almost indistinguishable**. This gives us a warning: in the early stages of a population phenomenon, **we cannot easily tell** whether it is growing exponentially or logistically — both models fit equally well. The difference emerges only as we approach $K$.

---

# What happens if we start **above** $K$?

```{r overshoot}
times2 <- seq(0, 100, by = 0.5)

scenarios <- data.frame(
  t   = rep(times2, 3),
  N   = c(logistic_analytical(times2, N0 = 100,   r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 5000,  r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 15000, r = 0.1, K = 10000)),
  N0  = factor(rep(c("N0 = 100",  "N0 = 5,000",  "N0 = 15,000 (overshoot)"),
                   each = length(times2)))
)

ggplot(scenarios, aes(x = t, y = N, colour = N0)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  annotate("text", x = max(times2) * 0.95, y = K * 1.04,
           label = "K = 10,000", colour = "grey40", hjust = 1) +
  labs(
    title = "Logistic growth: three different initial scenarios",
    subtitle = "All converge to K, regardless of the initial condition",
    x = "Time (days)",
    y = "N(t)"
  ) +
  theme(legend.position = "bottom")
```

Notice that the case $N_0 = 15,000 > K$ **decreases** toward $K$. This follows automatically from the equation: if $N > K$, then $1 - N/K < 0$, so $\frac{dN}{dt} < 0$.

This models, for example, the case where we introduce an excessive number of water lilies into a lake — the population will **fall** to the carrying level.

---

# Estimating parameters from data

In practice biologists have data and want to **estimate** $r$ and $K$. This is done with non-linear regression (`nls()` in R).

```{r fit-logistic}
# We simulate data with noise
set.seed(123)
sample_days <- seq(0, 150, by = 10)
true_N      <- logistic_analytical(sample_days, N0 = 100, r = 0.1, K = 10000)
observed_N  <- true_N * exp(rnorm(length(sample_days), 0, 0.08))

# Fitting the logistic model
fit <- nls(observed_N ~ K / (1 + ((K - N0)/N0) * exp(-r * sample_days)),
           start = list(K = 8000, r = 0.05, N0 = 200))

summary(fit)$coefficients
```

```{r plot-fit}
fit_curve <- data.frame(
  t = times,
  N = predict(fit, newdata = data.frame(sample_days = times))
)

ggplot() +
  geom_point(data = data.frame(t = sample_days, N = observed_N),
             aes(x = t, y = N), size = 2.5, colour = "tomato") +
  geom_line(data = fit_curve, aes(x = t, y = N),
            linewidth = 1, colour = "steelblue") +
  geom_hline(yintercept = coef(fit)["K"],
             linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95,
           y = coef(fit)["K"] * 1.05,
           label = paste0("Estimate K = ", round(coef(fit)["K"])),
           colour = "grey30", hjust = 1) +
  labs(
    title = "Fitting the logistic model to data",
    subtitle = paste0("Estimates: r = ", round(coef(fit)["r"], 3),
                      ", K = ",  round(coef(fit)["K"]),
                      ", N0 = ", round(coef(fit)["N0"])),
    x = "Time (days)",
    y = "Population N(t)"
  )
```

> **Practical significance:** The `nls()` estimated $\hat{r}$ and $\hat{K}$ come with a standard error — so we have statistical confidence intervals for the parameters. This is crucial for comparing two populations or two experimental conditions.

---

# Beyond the logistic: other functional forms

As we mentioned, the linear form $r_{\text{eff}}(N) = r(1 - N/K)$ is **one choice**. Two important alternatives:

## θ-logistic

$$
\frac{dN}{dt} = r N \left(1 - \left(\frac{N}{K}\right)^{\theta}\right)
$$

- $\theta = 1$ → classical logistic.
- $\theta > 1$ → the population stays close to $K$ for a longer time before slowing down (e.g. $r$-strategists).
- $\theta < 1$ → the slowdown starts early (e.g. species sensitive to competition).

## Gompertz

$$
\frac{dN}{dt} = r N \ln\!\left(\frac{K}{N}\right)
$$

Used mainly for the growth of masses (e.g. cancerous tumors, organism development). It slows down **earlier** than the logistic.

```{r alt-models}
gompertz <- function(t, N0, r, K) {
  K * exp(log(N0/K) * exp(-r * t))
}
theta_logistic <- function(t, N0, r, K, theta) {
  # Numerical solution
  parms <- c(r = r, K = K, theta = theta)
  func  <- function(t, state, parms) {
    with(as.list(c(state, parms)),
         list(r * N * (1 - (N/K)^theta)))
  }
  out <- ode(y = c(N = N0), times = t, func = func, parms = parms)
  as.data.frame(out)$N
}

t_vec <- seq(0, 200, by = 0.5)
df_alt <- data.frame(
  t        = rep(t_vec, 4),
  N        = c(logistic_analytical(t_vec, 100, 0.1, 10000),
               gompertz(t_vec,             100, 0.1, 10000),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 0.5),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 3)),
  Model    = factor(rep(c("Logistic (θ=1)",
                          "Gompertz",
                          "θ-logistic (θ=0.5)",
                          "θ-logistic (θ=3)"),
                        each = length(t_vec)),
                    levels = c("Logistic (θ=1)",
                               "θ-logistic (θ=0.5)",
                               "θ-logistic (θ=3)",
                               "Gompertz"))
)

ggplot(df_alt, aes(x = t, y = N, colour = Model)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 10000, linetype = "dashed", colour = "grey50") +
  labs(
    title    = "Different models of saturation",
    subtitle = "All models have the same r=0.1, K=10,000, N0=100",
    x = "Time (days)", y = "N(t)"
  ) +
  theme(legend.position = "bottom")
```

> **Pedagogical remark:** The fact that all these models give *different* curves is the most important message. On real data the choice of model **does not follow automatically** — it depends on the biology of the organism and must be checked statistically (e.g. AIC).

---

# Exercises

1. **Sensitivity to $K$.** Plot logistic curves with $r=0.1$, $N_0=100$ and $K \in \{1{,}000,\ 5{,}000,\ 10{,}000,\ 50{,}000\}$ on the same graph. What changes in the shape of the curve?

2. **Time to reach $90\%$ of $K$.** Solve analytically: for which $t$ does $N(t) = 0.9 K$ hold? Verify numerically.

3. **Inflection point.** Show that the rate $\frac{dN}{dt}$ is maximized at $N = K/2$. (*Hint:* differentiate $\frac{dN}{dt}$ with respect to $N$ and find the root.)

4. **From the data to $r$.** You are given measurements:
   ```r
   day <- c(0, 10, 20, 40, 60, 90, 120, 180)
   N   <- c(120, 280, 620, 2400, 5800, 9100, 9750, 9970)
   ```
   Fit the logistic model with `nls()` and compute the confidence intervals for $\hat{r}$ and $\hat{K}$ (`confint(fit)`).

5. **Conceptual.** In the discussion we mentioned *various* mechanisms that can justify the decrease of the per capita rate (space, light, nutrients, pathogens). For the water lilies in our lake, which mechanism do you consider **most important** and why? How would you verify it experimentally?

6. **Beyond the mean.** All the models here are deterministic — there is no chance involved. In real populations, however, there is stochasticity (births/deaths are random events). Try to simulate stochastic logistic growth: at each time step, add to $\frac{dN}{dt}$ a random term `rnorm(1, 0, σ)`. How does the behavior change for different $\sigma$?

---

# Next step

So far we have examined a *single* population. In reality, however, populations **interact** with each other:

- **Predator–prey**: Lotka–Volterra equations
- **Interspecific competition**: extended logistic with two or more equations
- **Symbiosis**: mutually beneficial interactions

These are **systems of differential equations** and will be the subject of the next laboratory.

---

# R session information

```{r session-info}
sessionInfo()
```