Skip to content

Repository files navigation

Population Dynamics

Mathematical theory + computational simulations of population dynamics

A from-scratch mathematical and computational treatment of population dynamics. The project begins with the logistic growth model and builds out its discrete, continuous, stochastic, and multi-species generalizations. Each model is implemented in R and paired with a self-contained theoretical write-up.

Every folder has its own README with the complete story: derivations, assumptions, and simulation outputs (trajectories, phase planes, bifurcation diagrams).


Repository structure

Stage Model Folder
1 Discrete time logistic growth model /Discrete_time_logistic_model
2 Continuous time logistic growth model /Continuous_time_logistic_model
3 Harvesting and Holling functional response models /Holling_functional_harvesting_model
4 Two species competition and prey-predator dynamics /Two_species_population_models
5 Multi species dynamics work in progress

What's worth a closer look

  • Discrete models: an explicit assumptions section, and stochasticity layered on top of the mean field structure to make precise what a mean field is actually averaging over.
  • Stability and predictability: equilibrium stability is worked out through Jacobian and eigenvalue analysis throughout. Lyapunov exponents and their classification are discussed below for context, as the natural next step, even though they are not directly computed in the analysis.
  • Two species models: Holling type functional responses and predator prey dynamics. Multi species generalizations are under progress.

On the theory

The theory sections aim for mathematical rigor and are meant to be self-contained, though space constraints mean some notation is used loosely, and not every hypothesis is restated at the point it is invoked. I would suggest going through the simulations first to build intuition. The equations tend to click once you have watched them play out; the theory reads easier after that.

Some constructions depart from standard references and reflect my own framing. Where something reads as ambiguous, it probably is. Corrections are welcome.

All models here are deterministic mean field models. The one exception is a limited, deliberate use of randomness in the discrete models, introduced specifically to make the mean field structure explicit. This is not stochastic modelling in the general sense.

Acknowledgments and scope

This project began as an extension of Prof. Guttal's EC201 course on YouTube. The mathematical framing draws on Strogatz's Nonlinear Dynamics and Chaos and a handful of research papers, alongside constructions I developed myself, which are not guaranteed to be error free. AI tools were used to help catch errors and gaps in the theory. The simulation code and the mathematical write ups are otherwise written by hand.



General definitions and results

Dynamical systems

State Space (Phase Space) :

The set of all possible values of a system’s state variables.

Example: State variable := Population size.

  • $E \subset \mathbb{R}^n$ - Single species : $E := \mathbb{R}_{\ge 0}$
  • $n$ species: $E := \mathbb{R}_{\ge 0}^n$

Note: Why not $E = \mathbb{Z}_{\ge 0}$? refer last section in discrete_time_logistic_model.

Dynamical functions :

Function that tells you the next state given the current state.

(In population dynamics; $T = \mathbb{Z}_{\ge 0}$ for discrete time, and $T = [0, \infty)$ for continuous time.)

Evolution function (discrete time) :

$F : E \to E$, s.t. $p_{t+1} = F(p_t)$ where $p_t$ is the system’s state at time $t \in \mathbb{Z}_{\ge 0}$.

Flow (discrete time) :

$\phi: T \times E \to E$, $(t, p) \rightarrow \phi_t(p)$ s.t. $\phi_t(p_0) := F^t(p_0)$ is defined as flow starting at state $p_0$. where $F^t := F \circ F \circ ... F, \text{ t-times}.$

Vector field (continuous time) :

$f : E \to \mathbb{R}^n$ s.t. $f(p) = dp/dt$. (maps each point in $E$ to its “velocity”)

Flow (continuous time) :

$\phi: T \times E \to E$, $(t, p) \rightarrow \phi_t(p)$ is a flow generated by vector field $f$ on $E$ s. t.

$$ \begin{aligned} &\phi_0(p) = p & \forall p \in E\\ & \frac{\partial}{\partial t}\phi_t(p) = f(\phi_t(p)) & \forall t \in T \\ & \phi_{t+s}(p) = \phi_{t}(\phi_{s}(p)) = \phi_{s}(\phi_{t}(p)) \end{aligned} $$

  • ($\phi_t(p)$ is the state at time $t$ evolved under $\dot x = f(x)$ starting at $p$.)

Note: You can generalize time by taking $T := \mathbb{R}$ . For $T := \mathbb{R}_{\ge 0}$ we are only looking at the forward flow from given initial conditions.

Trajectory and orbit :

For a fixed initial condition $p_0 \in E$, trajectory through $p_0$ is a curve $\gamma_{p_0} : T \rightarrow E$, s.t. $$\gamma_{p_0} (t) = \phi_t(p_0)$$ And its Image ${ \phi_t(p_0) | t \in T } \subset E$ is the corresponding orbit.

Parameter of dynamic function :

A quantity $\mu$ (scalar or vector) that appears in the defining equation $\dot x=f(x;\mu)$ or $x_{n+1}=f (x_n;\mu)$ but is held fixed while the state $p$ evolves in time. e.g. max growth rate, max harvesting rate etc. set $M \subset \mathbb{R}^d$ of all $\mu$ is called parameter space.


Fixed point / equilibrium :

A state $p^\ast$ that doesn't change under the dynamics.

Discrete : $p^\ast \in E$ s.t. $F(p^\ast) = p^\ast$.

Continuous : $p^\ast \in E$ s.t. $f(p^\ast) = 0$.

Stability of the Equilibrium

(We have two definitions of stability, the latter being the stronger one.)

Lyapunov stable : Family ${ \phi_t }$ is Equicontinuous at point $p^\ast$.

$\forall \varepsilon \gt 0, \exists \delta \gt 0$ s.t. $\boxed{ |p_0 - p^\ast | \lt \delta \implies | p_t - p^\ast | \lt \varepsilon, \forall t \ge 0}$

Small perturbation to $p^\ast$ stays small over the time.

Asymptotically stable : Family ${ \phi_t }$ is Equicontinuous at point $p^\ast$ and converges pointwise to $p^\ast$.

$\boxed{\exists \delta_0 \ge 0 \text{ s.t. } |p_0 - p^\ast | \lt \delta_0 \implies p_t \to p^\ast}$

Small perturbation to $p^\ast$ stays small and eventually dies out over the time.

Note : Equilibrium is local attractor.

Linear stability test :

Discrete, $E \subset \mathbb{R}$

Let $F : E \to E$ be $C^1$ near fixed point $p^\ast$, then

  • $|F'(p^{\ast})| < 1 \implies p^{\ast} \text{ asymptotically stable}$
  • $|F'(p^{\ast})| > 1 \implies p^{\ast} \text{ unstable}$
  • $|F'(p^{\ast})| = 1 \implies \text{Inconclusive (possibly bifurcation point)}$

Continuous, $E \subset \mathbb{R}^n$

Let $J := Df(p^\ast)$ (Jacobian) and $\lambda_i$ be eigenvalues of $J$, then

  • $\text{Re}(\lambda_i) \lt 0, \ \ \text{for all } \lambda_i \implies p^\ast \text{ asymptotically stable}$
  • $\text{Re}(\lambda_i) \lt 0, \ \ \text{for any } \lambda_i \implies p^\ast \text{ unstable}$
  • $\text{Re}(\lambda_i) = 0, \ \ \text{for all } \lambda_i \implies \text{Inconclusive (possibly bifurcation point)}$

Attractor

$\omega -$limit set of $p_0$ :

Set of all subsequential limits of the forward orbit of $p_0$.

$\omega (p_0) = \{ p \in E \ : \ \lim_{t_k \rightarrow \infty} \omega_{t_k}(p_0)= p, \ \text{ for any subsequence } (t_k)_{k \ge 0} \subset T \}$

Attracting Set $A(\mu)$ :

A non empty closed set $A = A(\mu) \subset E$ is an attracting set if

  • $A \text{ is invariant under } \phi_t \text{ i.e. } \phi_t (A) = A, \forall t$.
  • $\exists \text{ open } U, A \subset U \subset E, \text{ s.t. } \text{dist}( \phi_t (p_0), A) \rightarrow 0 \text{ as } t \rightarrow 0, \ \forall p_0 \in U$.
  • $\forall \varepsilon \gt 0, \exists \delta \gt 0 \text{ s.t. } \text{dist}(p_0, A) \lt \delta \implies \text{dist}(\phi_t(p_0), A) \lt \varepsilon, \ \forall t \ge 0$.

$A(\mu)$ is an attractor if it is closed, invariant, Asymptotically stable (i.e. attracting and Lyapunov stable) and minimal set of such kind (no closed subset of A is attractor).

Here $\mu$ is a parameter value of the dynamic system $\dot x = f (x; \mu)$. For different $\mu$ we have may different $A$, $U$ structure even when the dynamic function has same form.

Basin of Attractor :

The entire region of starting points that will eventually get pulled into the attractor.

Attractor point:

$a \in A(\mu)$ is an attractor point. > Note: each $a \in \omega(p_0)$ for some $p_0 \in U$ and for each $p_0 \in U, \ \omega(p_0) \subset A$.

Types of attractors $A(\mu)$ :

Fixed point : $A = { p^\ast }$ i.e. equilibrium.

Periodic orbit : $A = \tau$ : closed curve with $\phi_T (p) = p, \ \forall p \in \tau$ : prey predator cycle etc.

Quasiperiodic attractor : oscillation with no exact period and dense orbit in A.

Chaotic attractor : bounded but unpredictable long-term behavior.

Bifurcation point :

Critical threshold in a parameter (growth rate, harvesting rate, infection rate etc.) where systems stability structure changes. (i.e. $A(\mu)$ changes in number, type or stability of attractors).

Let $M \subset \mathbb{R}^n$ be Parameter Space, $\mu_0 \in M$ is bifurcation point of $f(p; \mu)$ if every neighborhood $N$ of $\mu_0$ contains $\mu_1 , \mu_2$ s.t. no homeomorphism $h: E \rightarrow E$ maps orbits of $f(\cdot ; \mu_1)$ onto $f(\cdot; \mu_2)$ preserving the orientation (time-direction).

Bifurcation diagram :

A map of the system’s long-run states as a control parameter varies, marking where regime shifts occur. $$\mathcal{B} := { (\mu, p) \in M \times E : p \in A(\mu) }$$

Curves distinguished by stability : solid = attractor, dashed = unstable invariant set. Plotted with $\mu$ on the x-axis and $p$ on y-axis.


Determinism, Robustness, Predictability, and Attractor Structure

Deterministic system :

Initial condition $p_0$ uniquely determines entire future trajectory $\phi_t(p_0)$. No stochasticity involved in the system.

Robustness :

$\mu \in M$ is robust if it’s not a bifurcation point. Stability structure is robust against perturbations in $\mu$.

Predictability at $p_0$ :

If the measured initial condition has a small error $\delta_0$ can we determine how this error evolves over time? if small initial measurement error is forgiving i.e. it does not blow up into a large discrepancy, so the measured trajectory $\phi_t(x_0+\delta_0)$ remains a good approximation to the true trajectory $\phi_t(x_0)$. so we can still approximately predict future states of system. There are two distinct questions to focus. 1. Are our predictions valid in time window $[0, t^\ast]$ ? 2. Can we predict long run fate of system i.e. What happens to $\delta(t)$ as $t \to\infty$ ; does it eventually decay, drift, or grow without bound? This asymptotic question is answered by the Lyapunov exponent.

Let $p_0$ be true state and $p_0 + \delta_0$ be measured state, then we Define error at time $t$ as $$\delta(t) := \phi_t(p_0 + \delta_0) - \phi_t(p_0)$$

$(\varepsilon_{\text{tol}}, t^*)-$predictable at $p_0$ :

Let’s Fix tolerance $\varepsilon_\text{tol} \gt 0$ and predictability horizon $t^\ast \gt 0$. Define the system is $(\varepsilon_\text{tol}, t^\ast)-predictable$ at $p_0$ if, $$\exists \xi \gt 0 \text{ s.t. } \lVert \delta(0) \rVert \lt \xi \implies \lVert \delta(t) \rVert \lt \varepsilon_\text{tol}, \ \forall t \in [0, t^*]$$

Measurement error smaller than some threshold $\xi$ stays below the tolerance $\varepsilon_{\text{tol}}$​ throughout the entire window $[0,t^\ast]$.

Lyapunov exponent :

Lets assume $\lVert \delta(t) \rVert \approx e^{ \lambda t} \lVert \delta(0) \rVert$ (by linearization of first order taylor expansion (as $f$ is $C^1$) of $f$ about $\phi_t(p_0)$.) This assumption is only valid on $[0, T]$ for $T$ not too large relative to $\ln(1/ \lVert \delta(0) \rVert )$.

Here $\lambda$ is called Lyapunov constant.

Now taking logarithm of both side we define it as:

Maximal Lyapunov exponent :
$$\lambda(p_0) := \lim_{t \to \infty} \lim_{\lVert \delta(0) \rVert\to 0} {\dfrac{1} {t}} {\ln{\dfrac{\lVert \delta(t) \rVert}{\lVert \delta(0) \rVert}}}$$

  • The limit $\lVert \delta(0) \rVert \to 0$ ensures validity of approximation at any time.

Discrete time Lyapunov exponent :
$$\lambda(x_0) := \lim_{n \to \infty } \frac{1}{n} \sum_{i = 0}^{n-1} \ln{|F'(p_i)|} $$

Lyapunov spectrum :
In $n$-dimensional phase space, we have $n$-basis directions along which we can make the perturbation in initial condition $p_0$ so we have $n$ Lyapunov exponents ${ \lambda_i }_n$ for each basis direction. I will upload pdf regarding the choice of basis (which not even an fixed orthogonal basis as standard sense) .

Work in progress

Cases of $\lambda$ :

  • $\lambda \lt 0$ : exponential decay of error in long run; but still we can have transient spikes at some finite t. This is exactly why $\lambda \lt 0$ does not imply Lyapunov stability in general.

  • $\lambda = 0$ : linear/polynomial drift between trajectories i.e. phase/timing drift but no decay.

  • $\lambda \gt 0$ : exponential growth of $\delta(t)$ (until gets bounded by some attractor structure)

Predictability horizon for $\lambda \gt 0$ :

Let $\varepsilon_{\text{tol}}$ be fixed tolerance and $\lVert \delta(t) \rVert$ be the error, then there exist predictability horizon $t^\ast$ s.t. $\lVert \delta(t) \rVert \lt \varepsilon_\text{tol } \ \forall t \lt t^\ast$ ; given by $t^\ast = \frac{1}{\lambda} \ln \frac{\varepsilon_\text{tol} }{ \lVert \delta(0) \rVert }$.

  • Note: So, $\lambda \gt 0$ and $t^\ast$ tell you the average/eventual rate, they do not certify smooth monotonic growth of $\delta(t)$ up to $t^\ast$. So we may still have transient spikes even at $t \lt t^\ast$.

Attractor structure :

For fixed $\mu$, $E$ decomposes as $E = \left(\bigsqcup_i U_i\right)\sqcup\Sigma$, : disjoint open basins $U_i$ of distinct minimal attractors $A_i(\mu)$ (possibly of different types simultaneously i.e. multistability), plus a measure-zero exceptional set $\Sigma$ (basin boundaries, stable manifolds of saddles).

Note:

  • $A_i(\mu)$ need not be connected (e.g. period- $k$ points, multi-band chaotic attractors).

  • $\mu_0$ is a bifurcation point if the collection ${A_i(\mu)}$ changes in number, type, or stability as $\mu$ crosses $\mu_0$.

Attractor Structure flow chart:


Further reading:

* Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, third edition by Steven H. Strogatz*

Perko, L. Differential Equations and Dynamical Systems, 3rd ed., Springer, 2001. §2.9 — attracting sets, Lyapunov stability, asymptotic stability; Poincaré maps for periodic orbits. Primary source for the core topological definitions used throughout.

Guckenheimer, J. & Holmes, P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, 1983. Def. 1.1.2 and §1.2 — attracting sets with open neighborhoods, return/Poincaré maps, structural stability. Standard reference for bifurcation and limit-cycle stability arguments.

Kuznetsov, Y. A. Elements of Applied Bifurcation Theory, 3rd ed., Springer, 2004. Def. 2.1 — topological equivalence of vector fields; defines bifurcation as its failure. Source for the “no homeomorphism between orbit structures” formulation of a bifurcation point.

Milnor, J. “On the concept of attractor,” Communications in Mathematical Physics 99 (1985), 177–195; erratum 102 (1985), 517–519. Introduces the measure-theoretic (as opposed to topological/Lyapunov) definition of attractor — positive-measure, not-necessarily-open basin.

Alexander, J. C., Kan, I., Yorke, J. A., You, Z. “Riddled Basins,” International Journal of Bifurcation and Chaos 2 (1992), 795–813. Constructs basins of attraction that are positive-measure but nowhere dense/topologically riddled — the concrete counterexample motivating Milnor’s weaker definition.

Ashwin, P. & Timme, M. “Unstable attractors: existence and robustness in networks of oscillators with delayed pulse coupling,” Nonlinearity 18 (2005), 2035–2060. Source for “unstable attractors”: sets that are attracting ($\lim_{t\to\infty}$ sense) but not Lyapunov stable, due to large transient excursions from non-normal dynamics.

Oseledets, V. I. “A multiplicative ergodic theorem: Lyapunov characteristic numbers for dynamical systems,” Transactions of the Moscow Mathematical Society 19 (1968), 197–231. Establishes existence of the Lyapunov spectrum and the associated (Oseledets) eigenbasis almost everywhere — the theoretical basis for defining $\lambda_1,\dots,\lambda_n$ rigorously.

Katok, A. & Hasselblatt, B. Introduction to the Modern Theory of Dynamical Systems, Cambridge University Press, 1995. §4 — periodic orbits, recurrence, $\omega$-limit sets; general reference for the topological/ergodic framework connecting invariance, recurrence, and attractors.