---
title: "Exponential Population Growth"
subtitle: "From the lake's water lilies to the differential 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"
)
# Packages we will need
# install.packages(c("ggplot2", "deSolve")) # run once if you don't have them
library(ggplot2)
library(deSolve)
theme_set(theme_minimal(base_size = 12))
```

# The problem: a lake with water lilies

Imagine that we arrive at a lake and count the water lilies.

- **Day 0:** $N_0 = 100$ water lilies
- **Day 1:** $N_1 = 110$ water lilies

Returning the next day we observe a **10% increase**. Our goal is:

> How, starting from this simple observation, do we arrive at the *exponential growth formula* by way of a differential equation?

The journey has three stops:

1. **Observation → rate of change** (discrete thinking)
2. **Discrete → continuous** (we take the limit and write the differential equation $\frac{dN}{dt} = rN$)
3. **Solving the differential equation** (separation of variables → $N(t) = N_0 e^{rt}$)

---

# From the observation to the rate of change

## The change in the population

Over a time $\Delta t = 1$ day, the change in the population is:

$$
\Delta N \;=\; N_1 - N_0 \;=\; 110 - 100 \;=\; 10
$$

## The **per capita** rate $r$

We are not interested in the absolute change (10 water lilies), because it depends on how many water lilies we had at the start. We are interested in **how much the population changed per individual, per unit of time**:

$$
r \;=\; \frac{\Delta N}{N \cdot \Delta t} \;=\; \frac{10}{100 \cdot 1} \;=\; 0.1 \;=\; 10\%/\text{day}
$$

$r$ is called the **intrinsic rate of increase** of the population. It is the fundamental parameter of all classical population dynamics.

```{r per-capita-rate}
# Computing r from the data
N0 <- 100
N1 <- 110
dt <- 1   # days

r <- (N1 - N0) / (N0 * dt)
r
```

---

# From the discrete to the continuous

## Discrete equation

If we assume that this rate stays constant, then at every time step $\Delta t$:

$$
\Delta N \;=\; r \, N \, \Delta t
\quad\Longleftrightarrow\quad
\frac{\Delta N}{\Delta t} \;=\; r N
$$

### example values for r=0.1 (10%)

| N | ΔΝ/Δt |
| :-------- | :-------- |
| 100  | 10 |
| 200  | 20 |
| 1000 | 100 |
| 5000 | 500  |



This is a **discrete** approximation: we measure the population at discrete steps (every day).

## The limit $\Delta t \to 0$

But the water lilies do not wait for dawn to grow — they grow *continuously*. Taking the limit as the time step becomes infinitesimal:

$$
\lim_{\Delta t \to 0} \frac{\Delta N}{\Delta t} \;=\; \frac{dN}{dt}
$$

We arrive at the **differential equation of exponential growth**:

$$
\boxed{\;\frac{dN}{dt} \;=\; r\,N\;}
$$

> **Intuitive interpretation:** The rate of change of the population at every moment is **proportional** to the current size of the population. The more water lilies there are, the faster new ones appear.

---

# Solving the differential equation

We will solve $\frac{dN}{dt} = rN$ using the method of **separation of variables**.

## Step 1 — Separate the variables

$$
\frac{dN}{N} \;=\; r\,dt
$$

(All the $N$ on the left-hand side, all the $t$ on the right.)

## Step 2 — Integrate both sides

$$
\int \frac{dN}{N} \;=\; \int r\,dt
\quad\Longrightarrow\quad
\ln(N) \;=\; r\,t + C
$$

where $C$ is a constant of integration.

## Step 3 — Solve for $N$

$$
N(t) \;=\; e^{rt + C} \;=\; e^{C} \cdot e^{rt}
$$

## Step 4 — Find the constant from the initial condition

At $t = 0$ we have $N(0) = N_0$, so:

$$
N_0 \;=\; e^{C} \cdot e^{0} \;=\; e^{C}
\quad\Longrightarrow\quad
e^{C} = N_0
$$

## The final formula

$$
\boxed{\;N(t) \;=\; N_0\, e^{r t}\;}
$$

This is the **exponential growth formula**.

---

# Implementation in R

## Analytical solution

Let's plot $N(t) = 100 \cdot e^{0.1\,t}$ for 60 days:

```{r analytical-solution}
N0 <- 100
r  <- 0.1
t  <- seq(0, 60, by = 0.1)

N_analytical <- N0 * exp(r * t)

df_analytical <- data.frame(t = t, N = N_analytical)

ggplot(df_analytical, aes(x = t, y = N)) +
  geom_line(linewidth = 1, colour = "steelblue") +
  labs(
    title = "Exponential growth of the water lily population",
    subtitle = expression(paste(N(t) == N[0] * e^{r*t}, ",   ",
                                N[0] == 100, ",   ", r == 0.1)),
    x = "Time (days)",
    y = "Population N(t)"
  )
```

In 60 days the population reaches `r round(N0 * exp(r * 60))` water lilies — unrealistic for a finite lake! This is the **problem** of exponential growth that we will address later with the logistic model.

## Numerical solution with `deSolve`

Very often in biology differential equations **cannot be solved analytically**. In those cases we use numerical integrators. For this particular model, which we already know analytically, this gives us an opportunity to check whether the numerical method agrees.

```{r numerical-solution}
# 1. Define the differential equation as a function
exponential_model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N      # dN/dt = rN
    list(dN)
  })
}

# 2. Parameters, initial condition, time
parameters <- c(r = 0.1)
state      <- c(N = 100)
times      <- seq(0, 60, by = 0.1)

# 3. Solve
out <- ode(y = state, times = times, func = exponential_model, parms = parameters)
out_df <- as.data.frame(out)

head(out_df)
```

## Comparison of the analytical and numerical solutions

```{r comparison-plot}
df_compare <- data.frame(
  t          = times,
  Analytical = N0 * exp(r * times),
  Numerical  = out_df$N
)

ggplot(df_compare, aes(x = t)) +
  geom_line(aes(y = Analytical),  linewidth = 1.2, colour = "steelblue") +
  geom_point(aes(y = Numerical),
             data = df_compare[seq(1, nrow(df_compare), by = 50), ],
             colour = "tomato", size = 2) +
  labs(
    title = "Analytical (line) versus numerical (dots) solution",
    x = "Time (days)",
    y = "Population N(t)"
  )

# Maximum error
max(abs(df_compare$Analytical - df_compare$Numerical))
```

The error is negligible — the two solutions coincide.

---

# A subtle but **very important** observation

Notice something: we started by saying *"10% increase per day"*. That is a **discrete** statement. But then we wrote $r = 0.1$ inside a **continuous** differential equation. Are the two really equal?

## The discrete formula

If every day the population is multiplied by $1.10$:

$$
N(t) \;=\; N_0 \cdot (1+r)^{t} \;=\; 100 \cdot 1.1^{t}
$$

## Comparison of the two models

```{r discrete-vs-continuous}
days <- 0:30
N_discrete   <- 100 * (1 + 0.1)^days
N_continuous <- 100 * exp(0.1 * days)

df_dc <- data.frame(
  day        = rep(days, 2),
  N          = c(N_discrete, N_continuous),
  model      = factor(rep(c("Discrete  N0·(1+r)^t",
                            "Continuous  N0·exp(rt)"), each = length(days)))
)

ggplot(df_dc, aes(x = day, y = N, colour = model)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
  labs(
    title = "Discrete versus continuous exponential growth",
    subtitle = "With the same value r = 0.1 the two models do **not** coincide",
    x = "Days", y = "Population", colour = ""
  ) +
  theme(legend.position = "bottom")
```

```{r dc-table}
data.frame(
  Day        = c(1, 5, 10, 20, 30),
  Discrete   = round(100 * 1.1^c(1, 5, 10, 20, 30), 1),
  Continuous = round(100 * exp(0.1 * c(1, 5, 10, 20, 30)), 1)
)
```

## Why do they differ?

In the **discrete** model, "10% increase/day" means: at the end of the day the population is 110% of the initial one.

In the **continuous** model with $r = 0.1$, the *instantaneous* rate is 10%, but because the newly born water lilies **immediately produce new water lilies themselves** (continuous compounding!), at the end of the day we have:

$$
N(1) \;=\; 100 \cdot e^{0.1} \;\approx\; 110.52
$$

that is, a little more than a 10% increase.

## How are they connected?

If we want the continuous model to give **exactly** a 10% increase per day, we need:

$$
e^{r_{\text{cont}}} = 1.1
\quad\Longrightarrow\quad
r_{\text{cont}} \;=\; \ln(1.1) \;\approx\; 0.0953
$$

```{r r-conversion}
log(1.1)   # in R, log() is the natural logarithm (ln)
```

> **In ecological practice**, most of the time we work with the continuous model and take $r$ to be the instantaneous rate. The difference is small for small $r$, but it is good to know that it **exists**.

---

# The logarithmic axis: the biologist's diagnostic tool

Taking the logarithm of the formula $N(t) = N_0 e^{rt}$:

$$
\ln N(t) \;=\; \ln N_0 \;+\; r\,t
$$

That is, $\ln N$ is a **linear** function of $t$ with slope equal to $r$. This is extremely useful in practice:

> **If your data on a semi-log plot ($\ln N$ against $t$) form a straight line, then the population is growing exponentially**, and the slope of the line gives you $r$.

```{r log-plot}
ggplot(df_analytical, aes(x = t, y = N)) +
  geom_line(linewidth = 1, colour = "steelblue") +
  scale_y_log10() +
  labs(
    title = "Exponential growth on a semi-log axis",
    subtitle = "The straight line indicates exponential growth",
    x = "Time (days)",
    y = "log10(N(t))"
  )
```

## Estimating $r$ from data with linear regression

```{r estimate-r}
# Simulated data with a little noise
set.seed(42)
sim_days <- 0:30
sim_N    <- 100 * exp(0.1 * sim_days) * exp(rnorm(length(sim_days), 0, 0.05))

# Linear regression log(N) ~ t
fit <- lm(log(sim_N) ~ sim_days)
summary(fit)$coefficients

cat("Estimated r =", round(coef(fit)[2], 4),
    "  (true r = 0.1)\n")
cat("Estimated N0 =", round(exp(coef(fit)[1]), 2),
    "  (true N0 = 100)\n")
```

---

# Doubling time

A useful quantity: how long does it take for the population to double?

$$
2 N_0 = N_0 e^{r T_2}
\;\Longleftrightarrow\;
T_2 = \frac{\ln 2}{r}
$$

```{r doubling-time}
T2 <- log(2) / r
cat("Doubling time of the water lilies:", round(T2, 2), "days\n")
```

So the water lilies double roughly every **7 days**.

---

# Exercises

1. **Play with $r$.** Plot on the same graph the exponential growth for $r = 0.05$, $r = 0.1$ and $r = 0.2$. How does the doubling time change?

2. **Negative $r$.** What happens if $r < 0$? Plot the case $r = -0.05$ and interpret it biologically.

3. **Estimation from real data.** You are given the following measurements:
   ```r
   day <- c(0, 3, 7, 14, 21, 28)
   N   <- c(50, 73, 109, 218, 437, 870)
   ```
   Estimate $r$ with linear regression on $\ln N$ and predict the population at 60 days.

4. **The mathematical danger.** With $r = 0.1$ and $N_0 = 100$ water lilies, in how many days will the water lilies exceed the number of humans on Earth ($\approx 8 \times 10^{9}$)? What does this tell us about the *limitations* of the exponential model?

---

# Next step: the law returns to reality

The exponential model **cannot** hold forever. In a real lake there are constraints: space, light, nutrients. At some point the population will *saturate*.

In the next lesson we will modify the differential equation so that it includes the **carrying capacity** $K$ of the environment:

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

This is the **logistic equation** — and it is the subject of the next laboratory.

---

# R session information

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