Mathematical biology

For almost a century the Hudson's Bay Company bought pelts from trappers across northern Canada, and its records show the numbers of lynx and snowshoe hare rising and falling in roughly ten-year cycles. In the Adriatic, the marine biologist Umberto D'Ancona noticed that during the First World War, when fishing was reduced, the proportion of predatory fish in the catch at ports such as Fiume (Rijeka) went up — and asked his father-in-law, the mathematician Vito Volterra, why. During the West African Ebola epidemic of 2014–16, teams of epidemiologists estimated how many people each case was infecting, to plan treatment centres and judge whether interventions were working. In each case a mathematical model turned a vague question into a precise one.

A model is a deliberately simplified description: we decide which quantities matter (population sizes, numbers infected), which processes change them (births, predation, infection), and write those processes as equations. Then mathematics tells us what the assumptions imply — equilibria, thresholds, oscillations, chaos — often things nobody expected. Finally the implications are compared with data, and the model is revised. The statistician George Box put the spirit memorably: "all models are wrong, but some are useful".

This lesson uses ordinary differential equations (solutions, equilibria, phase portraits) and linear algebra (eigenvalues, trace and determinant). Stochastic models, in which chance plays a role, connect to Markov chains.

A population with plentiful resources grows at a rate proportional to its size: $N' = rN$, so $N(t) = N_0e^{rt}$ — Malthusian, exponential growth. No real population can do this for long. Pierre-François Verhulst (1838) proposed that the per-capita growth rate falls linearly as the population approaches a carrying capacity $K$:

The picture already contains the whole qualitative story: two equilibria, $N = 0$ (unstable: a few individuals grow) and $N = K$ (stable: every positive population approaches it), and S-shaped growth in between. The same idea in discrete time — a population counted once per generation — produces something far richer, as we shall see.

Many insects, annual plants and fish breed once a year in non-overlapping generations, so populations are naturally modelled by maps. The logistic map $x_{n+1} = rx_n(1 - x_n)$ is the simplest non-linear one, and Robert May (1976) urged biologists to study it precisely because its behaviour is so rich.

For $f_r(x) = rx(1 - x)$: the fixed points are $x = 0$ and $x^* = 1 - \frac1r$ (in $(0, 1)$ for $r > 1$), with $f_r'(x) = r(1 - 2x)$, so $f_r'(0) = r$ and $f_r'(x^*) = 2 - r$. Hence extinction is stable for $r < 1$, and the positive equilibrium is stable exactly for $1 < r < 3$. At $r = 3$ the multiplier passes through $-1$, and a cycle of period two is born.

The road to chaos. As $r$ increases, the 2-cycle loses stability at $r = 1 + \sqrt 6$ and a stable 4-cycle appears, then an 8-cycle, and so on: a period-doubling cascade at parameters $r_1 = 3 < r_2 \approx 3.449 < r_3 \approx 3.544 < \cdots$ accumulating at $r_\infty \approx 3.5699$. The gaps shrink geometrically: $\frac{r_n - r_{n-1}}{r_{n+1} - r_n} \to \delta = 4.6692\ldots$. Mitchell Feigenbaum (1978) discovered numerically that the same $\delta$ appears for every smooth unimodal map with a quadratic maximum — a universal constant — and Oscar Lanford (1982) proved his conjectures with a computer-assisted proof. Beyond $r_\infty$, chaotic parameters and periodic "windows" interleave; the largest window, containing a stable 3-cycle, opens at $r = 1 + \sqrt 8 \approx 3.828$. By the Li–Yorke theorem (1975), a continuous map of an interval with a point of period three has points of every period (a special case of Sharkovskii's theorem, 1964).

Chaos at $r = 4$. Substituting $x = \sin^2(\pi\theta)$ gives $f_4(x) = 4\sin^2(\pi\theta)\cos^2(\pi\theta) = \sin^2(2\pi\theta)$: the logistic map is conjugate to the doubling map $\theta \mapsto 2\theta \bmod 1$, which shifts the binary digits of $\theta$ one place to the left. Two starting points that agree to $40$ binary digits disagree completely after about $40$ steps: sensitive dependence on initial conditions, with errors doubling at each step (Lyapunov exponent $\ln 2$). A deterministic model with one parameter can thus be unpredictable in practice — a sobering lesson for population forecasting, since field data with noisy, chaotic-looking fluctuations may come from very simple deterministic rules.

The Hartman–Grobman theorem says more: if no eigenvalue has zero real part (a hyperbolic equilibrium), the phase portrait near $x^*$ is a continuous deformation of that of $y' = Jy$. When some eigenvalue has zero real part — for example a centre, with purely imaginary eigenvalues — linearisation is inconclusive, and the nonlinear terms decide.

In the plane, the eigenvalues of $J$ are the roots of $\lambda^2 - \tau\lambda + \Delta$ with $\tau = \operatorname{tr}J$ and $\Delta = \det J$, so both have negative real part if and only if $\tau < 0$ and $\Delta > 0$ (the trace–determinant test). The "if" direction is proved in Lean: if $\lambda = a + bi$ is a root, the imaginary part of the equation gives $b(2a - \tau) = 0$, so either $a = \tau/2 < 0$, or $b = 0$ and $a^2 - \tau a + \Delta = 0$, which is impossible for $a \ge 0$. A saddle has $\Delta < 0$.

In the Lotka–Volterra model, prey grow exponentially without predators ($\alpha x$), predators starve without prey ($-\gamma y$), and encounters, at a rate proportional to $xy$ (the law of mass action), remove prey and feed predators. The equilibria are $(0, 0)$ and $(x^*, y^*) = (\gamma/\delta, \alpha/\beta)$, with Jacobians

The origin is a saddle. At the coexistence equilibrium $\tau = 0$ and $\Delta = \alpha\gamma > 0$, so the eigenvalues are $\pm i\sqrt{\alpha\gamma}$: a linear centre, and the linearisation theorem says nothing. A conserved quantity settles the question.

This is Volterra's principle. Fishing that removes both species in proportion to their numbers changes $\alpha$ to $\alpha - \varepsilon$ and $\gamma$ to $\gamma + \varepsilon$, so the average prey $(\gamma + \varepsilon)/\delta$ rises and the average predators $(\alpha - \varepsilon)/\beta$ fall. Reduced fishing during the war did the opposite — D'Ancona's observation explained. The same principle warns that spraying an insecticide that kills both a pest and its natural enemy can increase the pest on average.

Structural instability. The closed orbits are fragile: almost any change to the model destroys them. With self-limiting prey, $x' = \alpha x(1 - x/K) - \beta xy$, the coexistence equilibrium (when it exists) has $\tau < 0$ and $\Delta > 0$ and becomes a stable spiral: oscillations damp out. With a saturating ("Holling type II", 1959) predator response $\beta xy/(1 + hx)$, as in the Rosenzweig–MacArthur model (1963), the equilibrium can lose stability through a Hopf bifurcation and give way to a stable limit cycle — robust oscillations of a fixed amplitude, closer to the lynx–hare records. Rosenzweig's "paradox of enrichment" (1971) is a consequence: increasing the prey's carrying capacity can destabilise the system. A model whose qualitative behaviour changes under arbitrarily small perturbations should not be trusted for quantitative predictions.

Two species competing for the same resources follow $x' = x(1 - x - a_{12}y)$, $y' = \rho y(1 - y - a_{21}x)$, where $a_{12}$ measures the effect of $y$ on $x$ relative to $x$'s effect on itself. Besides $(0, 0)$, $(1, 0)$ and $(0, 1)$, there may be a coexistence equilibrium

which is positive when $a_{12}, a_{21}$ are both less than $1$ or both greater than $1$.

Coexistence requires that each species limits itself more than it limits its competitor. Georgii Gause's laboratory experiments with Paramecium (1934) supported the resulting competitive exclusion principle: two species competing for exactly the same limiting resource cannot coexist indefinitely.

The SIR model assumes a closed population (no births, deaths or migration on the time scale of the epidemic), homogeneous mixing (everyone equally likely to meet everyone), infection at rate $\beta SI/N$ (mass action), recovery at constant rate $\gamma$ (an exponentially distributed infectious period with mean $1/\gamma$), and lasting immunity. Adding the three equations, $(S + I + R)' = 0$: the population $N$ is constant.

Herd immunity and overshoot. If a fraction $p$ is immune (by vaccination or previous infection), an introduced infection grows only if $R_0(1 - p) > 1$; the herd immunity threshold is $p_c = 1 - 1/R_0$ — about $67\%$ for $R_0 = 3$ and $92$–$94\%$ for measles, whose $R_0$ is classically estimated at $12$–$18$. An unchecked epidemic infects more than this — $94\%$ for $R_0 = 3$ — because infections already under way continue after $S/N$ has fallen below $1/R_0$. The peak prevalence, from the invariant at the moment $S = N/R_0$, is $I_{\max}/N = 1 - \frac{1 + \ln R_0}{R_0}$ (for $S_0 = N$).

What $R_0$ is, and is not. $R_0 = \beta/\gamma$ combines a contact rate and a duration: it depends on how people meet (household size, schools, markets, transport), on the environment (season, ventilation, mosquito density for malaria), on the pathogen, and on behaviour. It is a property of a pathogen in a population, not a biological constant. Estimates for the same disease vary widely between settings, and the effective reproduction number $R_t$, which accounts for immunity and interventions, is what epidemic monitoring tracks. In 2014 the WHO Ebola Response Team estimated reproduction numbers of roughly $1.5$–$2$ in Guinea, Liberia and Sierra Leone before control measures took hold.

Small numbers and chance. Deterministic equations describe large populations. When an infection is introduced by one person, chance dominates: each case might infect nobody. In the stochastic SIR model the early phase is a branching process, and with an exponentially distributed infectious period the probability that the introduction dies out by chance is $1/R_0$ (for $R_0 > 1$). Small populations also drift genetically and can go extinct by bad luck even when the deterministic model predicts persistence; the Wright–Fisher chain in Markov chains is the basic model.

Space. Populations spread: Fisher (1937) and Kolmogorov, Petrovsky and Piskunov (1937) showed that $u_t = Du_{xx} + ru(1 - u)$ has travelling waves moving at speed $2\sqrt{rD}$, a model for the spread of advantageous genes and invading species. Alan Turing (1952) showed that two diffusing, reacting chemicals can spontaneously form spatial patterns — a proposed mechanism for spots and stripes on animal coats, now supported in several developmental systems.

Models are not forecasts. A model's output is a conditional statement: if these mechanisms, with these parameters, then this. Its predictions are only as good as its assumptions, and structural choices (homogeneous mixing versus contact networks, fixed versus changing behaviour) can matter more than parameter values. Epidemic models are often used to compare scenarios — what happens with and without an intervention — rather than to forecast the future, and because people change their behaviour in response to an epidemic (and to model projections), a projection that leads to action may never come true, which does not make it wrong. Good practice reports parameter uncertainty, tests sensitivity to assumptions, compares several model structures, and confronts the model with data.

Mathlib has the calculus needed here — HasDerivAt, logarithms and exponentials, the fundamental theorem of calculus — but no library of biological models, and not the linearisation theorem, Hartman–Grobman or Poincaré–Bendixson. The module proves the lesson's computations as theorems about arbitrary solutions: a hypothesis HasDerivAt x (α * x t - β * x t * y t) t says that x satisfies the prey equation at time t.

Logistic growth: the explicit solution, verified by differentiation, and its limit:

/-- **The logistic equation** `N' = r N (1 − N/K)` is solved by
`N(t) = K / (1 + C e^{−rt})` (found by separation of variables; here verified by
differentiation). With `C = K/N₀ − 1` this is the solution with `N(0) = N₀`. -/
theorem logistic_ode_solves {r K C : ℝ} (hK : K ≠ 0) (hC : 0 ≤ C) (t : ℝ) :
    HasDerivAt (fun s => K / (1 + C * exp (-(r * s))))
      (r * (K / (1 + C * exp (-(r * t)))) * (1 - K / (1 + C * exp (-(r * t))) / K)) t := by
  have hpos : 0 < 1 + C * exp (-(r * t)) := by positivity
  have h0 := (hasDerivAt_id t).const_mul (-r)
  have hinner : HasDerivAt (fun s => -(r * s)) (-r) t := by
    convert h0 using 1
    · funext s
      simp only [id]
      ring
    · ring
  have hd : HasDerivAt (fun s => 1 + C * exp (-(r * s))) (C * (exp (-(r * t)) * (-r))) t :=
    ((hinner.exp).const_mul C).const_add 1
  have h := (hasDerivAt_const t K).div hd hpos.ne'
  convert h using 1
  field_simp
  ring

/-- Every positive logistic solution approaches the **carrying capacity**: if `r > 0` then
`K / (1 + C e^{−rt}) → K` as `t → ∞`, because `e^{−rt} → 0`. -/
theorem logistic_ode_limit {r K C : ℝ} (hr : 0 < r) :
    Tendsto (fun t => K / (1 + C * exp (-(r * t)))) atTop (𝓝 K) := by
  have hexp : Tendsto (fun t => exp (-(r * t))) atTop (𝓝 0) := by
    have : Tendsto (fun t => r * t) atTop atTop := tendsto_id.const_mul_atTop hr
    exact tendsto_exp_neg_atTop_nhds_zero.comp this |>.congr fun t => rfl
  have hden : Tendsto (fun t => 1 + C * exp (-(r * t))) atTop (𝓝 1) := by
    simpa using (hexp.const_mul C).const_add 1
  have h := (tendsto_const_nhds (x := K)).div hden one_ne_zero
  rw [div_one] at h
  exact h

The logistic map: fixed points, the multiplier $2 - r$, the stability range $1 < r < 3$, the 2-cycle and its multiplier; and Mathlib's Banach fixed-point theorem, the principle behind the stability criterion:

/-- The **logistic map** `f_r(x) = r x (1 − x)`, a discrete-time model of a population
(as a fraction of its maximum) in non-overlapping generations. -/
def logisticMap (r x : ℝ) : ℝ := r * x * (1 - x)

/-- **Fixed points of the logistic map.** For `r ≠ 0`, `f_r(x) = x` exactly when `x = 0`
(extinction) or `x = 1 − 1/r` (the non-trivial equilibrium, in `(0, 1)` when `r > 1`).
Proof: `r x (1 − x) − x = x (r − 1 − r x)`. -/
theorem logistic_fixed_points {r : ℝ} (hr : r ≠ 0) (x : ℝ) :
    logisticMap r x = x ↔ x = 0 ∨ x = 1 - 1 / r := by
  unfold logisticMap
  constructor
  · intro h
    have hfac : x * (r - 1 - r * x) = 0 := by linarith [h]
    rcases mul_eq_zero.mp hfac with h0 | h1
    · exact Or.inl h0
    · right
      field_simp
      linarith
  · rintro (rfl | rfl)
    · ring
    · field_simp
      ring

/-- The derivative of the logistic map is `f_r'(x) = r (1 − 2x)`; at the non-trivial fixed
point `x* = 1 − 1/r` it equals `2 − r`. -/
theorem logistic_map_deriv {r : ℝ} (hr : r ≠ 0) :
    HasDerivAt (logisticMap r) (2 - r) (1 - 1 / r) := by
  have h : HasDerivAt (fun x => r * x * (1 - x)) (r * 1 * (1 - (1 - 1 / r)) + r * (1 - 1 / r) * (-1))
      (1 - 1 / r) :=
    ((hasDerivAt_id' _).const_mul r).mul ((hasDerivAt_id' _).const_sub 1 |>.congr_deriv (by simp))
  have hv : r * 1 * (1 - (1 - 1 / r)) + r * (1 - 1 / r) * (-1) = 2 - r := by
    field_simp
    ring
  rw [hv] at h
  exact h

/-- **Stability of the non-trivial equilibrium.** A fixed point of a smooth map is
asymptotically stable when `|f'(x*)| < 1`; here `|2 − r| < 1` exactly when `1 < r < 3`.
At `r = 3` the multiplier passes through `−1` and the period-doubling cascade begins. -/
theorem logistic_stable_iff (r : ℝ) : |2 - r| < 1 ↔ 1 < r ∧ r < 3 := by
  rw [abs_lt]
  constructor <;> intro h <;> constructor <;> linarith [h.1, h.2]

/-- **The period-2 cycle.** For `r ≥ 3`, with `s = √((r − 3)(r + 1))`, the points
`x± = (r + 1 ± s)/(2r)` are swapped by the logistic map: `f_r(x₊) = x₋` and `f_r(x₋) = x₊`.
(They coincide with the fixed point `1 − 1/r` at `r = 3` and separate for `r > 3`.)
Proof: expand, using `s² = r² − 2r − 3`. -/
theorem logistic_period_two {r : ℝ} (hr : 3 ≤ r) :
    let s := Real.sqrt ((r - 3) * (r + 1))
    logisticMap r ((r + 1 + s) / (2 * r)) = (r + 1 - s) / (2 * r) ∧
      logisticMap r ((r + 1 - s) / (2 * r)) = (r + 1 + s) / (2 * r) := by
  intro s
  have hr0 : r ≠ 0 := by linarith
  have hs : s ^ 2 = (r - 3) * (r + 1) := Real.sq_sqrt (by nlinarith)
  unfold logisticMap
  constructor
  · field_simp
    linear_combination -hs
  · field_simp
    linear_combination -hs

/-- **Stability of the 2-cycle.** The multiplier of the cycle — the derivative of `f_r ∘ f_r`
along it, `f_r'(x₊) f_r'(x₋) = r²(1 − 2x₊)(1 − 2x₋)` — equals `−r² + 2r + 4`. It lies in
`(−1, 1)` exactly for `3 < r < 1 + √6 ≈ 3.449`, where the 2-cycle is stable and the next
doubling begins. -/
theorem period_two_multiplier {r : ℝ} (hr : 3 ≤ r) :
    let s := Real.sqrt ((r - 3) * (r + 1))
    (r * (1 - 2 * ((r + 1 + s) / (2 * r)))) * (r * (1 - 2 * ((r + 1 - s) / (2 * r)))) =
      -r ^ 2 + 2 * r + 4 := by
  intro s
  have hr0 : r ≠ 0 := by linarith
  have hs : s ^ 2 = (r - 3) * (r + 1) := Real.sq_sqrt (by nlinarith)
  field_simp
  linear_combination -hs

/-- **Why stable fixed points attract: the Banach fixed-point theorem.** If a map of a complete
metric space is a contraction (`dist (f x) (f y) ≤ K dist x y` with `K < 1`), its iterates
from any starting point converge to its unique fixed point. Near a fixed point with
`|f'(x*)| < 1` a smooth map is a contraction of a small interval, which is the proof of the
stability criterion. This is Mathlib's `ContractingWith.tendsto_iterate_fixedPoint`. -/
theorem contraction_converges {α : Type*} [MetricSpace α] [CompleteSpace α] [Nonempty α]
    {K : NNReal} {f : α → α} (hf : ContractingWith K f) (x : α) :
    Tendsto (fun n => f^[n] x) atTop (𝓝 (ContractingWith.fixedPoint f hf)) :=
  hf.tendsto_iterate_fixedPoint x

Predator and prey: Volterra's conserved quantity along any positive solution, and the prey average over a cycle, via the fundamental theorem of calculus:

/-- **The Lotka–Volterra conserved quantity.** Along any positive solution of
`x' = αx − βxy` (prey), `y' = δxy − γy` (predators), the function
`V(x, y) = δx − γ log x + βy − α log y` has zero derivative, so it is constant and orbits lie
on its level curves (which are closed curves around the equilibrium `(γ/δ, α/β)`).
Proof: `V' = (δ − γ/x) x' + (β − α/y) y' = (δx − γ)(α − βy) + (βy − α)(δx − γ) = 0`. -/
theorem lotka_volterra_conserved {α β γ δ : ℝ} {x y : ℝ → ℝ} {t : ℝ} (hx : 0 < x t)
    (hy : 0 < y t) (hx' : HasDerivAt x (α * x t - β * x t * y t) t)
    (hy' : HasDerivAt y (δ * x t * y t - γ * y t) t) :
    HasDerivAt (fun s => δ * x s - γ * log (x s) + β * y s - α * log (y s)) 0 t := by
  have h1 := hx'.const_mul δ
  have h2 := (hx'.log hx.ne').const_mul γ
  have h3 := hy'.const_mul β
  have h4 := (hy'.log hy.ne').const_mul α
  have h := ((h1.sub h2).add h3).sub h4
  convert h using 1
  field_simp
  ring

/-- **Averages over a cycle equal the equilibrium.** If a positive solution is periodic with
period `T > 0` (`y(T) = y(0)`), then the average prey population over a period is
`(1/T) ∫₀ᵀ x = γ/δ`. Proof: `(log y)' = y'/y = δx − γ`; integrate over a period, where
`log y` returns to its initial value. (This is Volterra's explanation of D'Ancona's
observation that reduced fishing during the First World War raised the share of predators.) -/
theorem lotka_volterra_average {α β γ δ T : ℝ} (hδ : δ ≠ 0) (hT : 0 < T) {x y : ℝ → ℝ}
    (hy : ∀ t, 0 < y t) (hx' : ∀ t, HasDerivAt x (α * x t - β * x t * y t) t)
    (hy' : ∀ t, HasDerivAt y (δ * x t * y t - γ * y t) t) (hper : y T = y 0) :
    (∫ t in (0 : ℝ)..T, x t) / T = γ / δ := by
  have hxc : Continuous x := continuous_iff_continuousAt.mpr fun t => (hx' t).continuousAt
  -- `(log y)' = δ x − γ`.
  have hlog : ∀ t, HasDerivAt (fun s => log (y s)) (δ * x t - γ) t := by
    intro t
    have := (hy' t).log (hy t).ne'
    convert this using 1
    field_simp [(hy t).ne']
  have hint : IntervalIntegrable (fun t => δ * x t - γ) MeasureTheory.volume 0 T :=
    ((hxc.const_mul δ).sub continuous_const).intervalIntegrable 0 T
  have hftc := intervalIntegral.integral_eq_sub_of_hasDerivAt (fun t _ => hlog t) hint
  rw [hper, sub_self, intervalIntegral.integral_sub
    ((hxc.const_mul δ).intervalIntegrable 0 T) intervalIntegrable_const,
    intervalIntegral.integral_const_mul, intervalIntegral.integral_const, smul_eq_mul,
    sub_zero] at hftc
  field_simp
  linarith

Epidemics: conservation of the population, the threshold, the invariant behind the final-size relation, and the herd immunity threshold:

/-- **The SIR model conserves the population.** If `S' = −βSI/N`, `I' = βSI/N − γI` and
`R' = γI`, then `(S + I + R)' = 0`: everyone is always susceptible, infected or removed. -/
theorem sir_total_constant {β γ N : ℝ} {S I R : ℝ → ℝ} {t : ℝ}
    (hS : HasDerivAt S (-(β * S t * I t / N)) t)
    (hI : HasDerivAt I (β * S t * I t / N - γ * I t) t) (hR : HasDerivAt R (γ * I t) t) :
    HasDerivAt (fun s => S s + I s + R s) 0 t := by
  convert (hS.add hI).add hR using 1
  ring

/-- **The epidemic threshold.** When infection is present (`I > 0`), the number infected is
increasing (`I' > 0`) exactly when `R₀ · S/N > 1`, where `R₀ = β/γ` (`β, γ, N > 0`). -/
theorem sir_threshold {β γ N S I : ℝ} (hγ : 0 < γ) (hN : 0 < N) (hI : 0 < I) :
    0 < β * S * I / N - γ * I ↔ 1 < β / γ * (S / N) := by
  have key : β * S * I / N - γ * I = γ * I * (β / γ * (S / N) - 1) := by
    field_simp
  rw [key]
  constructor
  · intro h
    have := (mul_pos_iff_of_pos_left (mul_pos hγ hI)).mp h
    linarith
  · intro h
    exact mul_pos (mul_pos hγ hI) (by linarith)

/-- **The invariant behind the final-size relation.** Along a solution of the SIR model with
`S > 0`, the quantity `S + I − (N/R₀) log S` (with `R₀ = β/γ`, written `γN/β`) is constant.
Proof: its derivative is `S' + I' − (γN/β) S'/S = −γI + γI = 0`. Comparing its values at the
start (`I ≈ 0`, `S ≈ N`) and the end (`I = 0`) of an epidemic gives the final-size relation
`log(S∞/S₀) = R₀ (S∞ − S₀)/N`. -/
theorem sir_invariant {β γ N : ℝ} (hβ : β ≠ 0) (hN : N ≠ 0) {S I : ℝ → ℝ} {t : ℝ}
    (hpos : 0 < S t) (hS : HasDerivAt S (-(β * S t * I t / N)) t)
    (hI : HasDerivAt I (β * S t * I t / N - γ * I t) t) :
    HasDerivAt (fun s => S s + I s - γ * N / β * log (S s)) 0 t := by
  convert (hS.add hI).sub ((hS.log hpos.ne').const_mul (γ * N / β)) using 1
  field_simp
  ring

/-- **Herd immunity threshold.** If a fraction `p` of the population is immune, an
introduced infection grows only if `R₀ (1 − p) > 1`. So it cannot grow exactly when
`p ≥ 1 − 1/R₀` (for `R₀ > 0`): about `92–94%` for measles (`R₀ ≈ 12–18`). -/
theorem herd_immunity {R₀ p : ℝ} (hR : 0 < R₀) : R₀ * (1 - p) ≤ 1 ↔ 1 - 1 / R₀ ≤ p := by
  constructor
  · intro h
    have : 1 - p ≤ 1 / R₀ := by
      rw [le_div_iff₀ hR]
      linarith
    linarith
  · intro h
    have : 1 - p ≤ 1 / R₀ := by linarith
    rw [le_div_iff₀ hR] at this
    linarith

Linear stability: the trace–determinant test, for complex roots of the characteristic polynomial:

/-- **The trace–determinant test.** If a real `2 × 2` Jacobian has trace `τ < 0` and
determinant `Δ > 0`, every root `λ` of its characteristic polynomial `λ² − τλ + Δ` has
negative real part, so the equilibrium is locally asymptotically stable (by the linearisation
theorem). Proof: write `λ = a + bi`. The imaginary part of the equation is `b(2a − τ) = 0`;
if `b ≠ 0` then `a = τ/2 < 0`, and if `b = 0` then `a² − τa + Δ = 0` forces `a < 0`. -/
theorem trace_det_stable {τ Δ : ℝ} (hτ : τ < 0) (hΔ : 0 < Δ) {z : ℂ}
    (hz : z ^ 2 - τ * z + Δ = 0) : z.re < 0 := by
  have hre := congrArg Complex.re hz
  have him := congrArg Complex.im hz
  simp [pow_two] at hre him
  by_cases hb : z.im = 0
  · rw [hb] at hre
    by_contra hnn
    push Not at hnn
    nlinarith
  · have h2 : 2 * z.re - τ = 0 := by
      have : z.im * (2 * z.re - τ) = 0 := by linarith
      rcases mul_eq_zero.mp this with h | h
      · exact absurd h hb
      · exact h
    linarith

What is and is not formalised. Proved here: all the identities and inequalities above — explicit solutions, fixed points and multipliers, the 2-cycle, the conserved quantities, the cycle average of the prey, the SIR threshold and invariant, herd immunity and the trace–determinant test (in the "if" direction). Cited from Mathlib: the Banach fixed-point theorem. Not formalised: existence and uniqueness of solutions for these models (Mathlib's Picard–Lindelöf theorem would provide them), the linearisation theorem, the closedness of the Lotka–Volterra level curves, the final-size limit argument ($I \to 0$), the competition analysis, period doubling beyond period two, Feigenbaum universality and chaos.

First, the logistic map: attractor periods, superstable parameters found by Newton's method with the ratios of their gaps approaching Feigenbaum's constant, and Lyapunov exponents.

"""The logistic map x -> r x (1 - x): period doubling, Feigenbaum's constant and chaos.

For each growth rate r we iterate the map past its transient and detect the
period of the attractor. Then we locate the "superstable" parameters R_n at
which the critical point 1/2 lies on a cycle of period 2^n, by Newton's
method, and form the ratios (R_n - R_(n-1)) / (R_(n+1) - R_n), which approach
Feigenbaum's constant 4.669... Finally the Lyapunov exponent - the average of
log|f'(x_k)| along the orbit - is negative on stable cycles and equals log 2
at r = 4, where nearby orbits separate exponentially.
"""

import math

def logistic(r, x):
    """One step of the logistic map."""
    return r * x * (1 - x)

def orbit(r, x0, n):
    """[x_0, ..., x_n] for the logistic map."""
    xs = [x0]
    for _ in range(n):
        xs.append(logistic(r, xs[-1]))
    return xs

def attractor_period(r, x0=0.2, transient=3000, max_period=64, tol=1e-9):
    """Smallest p <= max_period with x_(n+p) = x_n after the transient, else None."""
    x = x0
    for _ in range(transient):
        x = logistic(r, x)
    start = x
    for p in range(1, max_period + 1):
        x = logistic(r, x)
        if abs(x - start) < tol:
            return p
    return None

def superstable_gap(r, n):
    """f_r^(2^n)(1/2) - 1/2: zero exactly when 1/2 lies on a 2^n-cycle."""
    x = 0.5
    for _ in range(2 ** n):
        x = logistic(r, x)
    return x - 0.5

def superstable_parameters(count):
    """R_0 = 2, R_1 = 1 + sqrt 5, ... found by Newton's method with a numerical
    derivative, each started from an extrapolation of the previous gaps."""
    rs = [2.0, 1 + math.sqrt(5)]
    delta = 4.0
    while len(rs) < count:
        n = len(rs)
        r = rs[-1] + (rs[-1] - rs[-2]) / delta
        for _ in range(60):
            g = superstable_gap(r, n)
            h = 1e-10
            dg = (superstable_gap(r + h, n) - superstable_gap(r - h, n)) / (2 * h)
            step = g / dg
            r -= step
            if abs(step) < 1e-14:
                break
        rs.append(r)
        delta = (rs[-2] - rs[-3]) / (rs[-1] - rs[-2])
    return rs

def lyapunov(r, x0=0.2, transient=1000, n=20000):
    """Average of log|f_r'(x_k)| = log|r (1 - 2 x_k)| along the orbit."""
    x = x0
    for _ in range(transient):
        x = logistic(r, x)
    total = 0.0
    for _ in range(n):
        total += math.log(abs(r * (1 - 2 * x)) + 1e-300)
        x = logistic(r, x)
    return total / n

def main():
    print("Attractor period after 3000 iterations from x0 = 0.2:")
    for r in [2.8, 3.2, 3.5, 3.56, 3.7, 3.83]:
        p = attractor_period(r)
        print(f"  r = {r:<5} period {p if p else 'none found (chaotic?)'}")
    rs = superstable_parameters(8)
    print()
    print("Superstable parameters R_n and Feigenbaum ratios:")
    for n, r in enumerate(rs):
        line = f"  R_{n} = {r:.10f}"
        if n >= 2:
            ratio = (rs[n - 1] - rs[n - 2]) / (rs[n] - rs[n - 1])
            line += f"   ratio {ratio:.5f}"
        print(line)
    print("  (Feigenbaum's constant: 4.66920...)")
    print()
    print("Lyapunov exponents (negative: stable cycle; positive: chaos):")
    for r in [2.8, 3.2, 3.5, 3.7, 3.83, 4.0]:
        print(f"  r = {r:<5} {lyapunov(r):+.4f}")
    print(f"  log 2 = {math.log(2):.4f}")
    a, b = orbit(4.0, 0.3, 40), orbit(4.0, 0.3 + 1e-10, 40)
    gaps = [abs(x - y) for x, y in zip(a, b)]
    first = next(i for i, g in enumerate(gaps) if g > 0.1)
    print(f"At r = 4, starts 1e-10 apart differ by > 0.1 after {first} steps.")
    print()
    print("Epistemic status: floating-point iteration of finitely many steps.")
    print("Periods, ratios and exponents are numerical estimates; the")
    print("universality of Feigenbaum's constant was proved by Lanford (1982)")
    print("with computer assistance, and chaos at r = 4 is a theorem.")

if __name__ == "__main__":
    main()