Introduction to Lotka-Volterra Dynamics

Modeling real-world phenomena such as neural dynamics, population growth, and the spread of infections is of significant interest to researchers and practitioners. In this blog post, I will focus on phenomena that evolve over time. In other words, these phenomena result in time series data. Throughout the blog, I use predator–prey dynamics as an example to explain the key ideas and concepts that come along the way.

Mathematicians Vito Volterra and Alfred Lotka, in the 1920s, independently proposed a system of linear differential equations to characterize the oscillating nature of hare and lynx populations observed in Canada. These oscillations were observed and recorded between the years of 1821 and 1914 by a fur-trading company called Hudson Bay Company. The number of hare and lynx in the area were estimated by counting the number of pelts since it is reasonable to assume that the number of hare/lynx pelts is approximately proportional to the number of hare/lynx species. The resulting mathematical model proposed is the classic and well-known Lotka-Volterra model in mathematical biology.

Let $X(t)$ and $Y(t)$ denote, respectively, the number of predators (e.g. lynx) and prey (e.g. hare) at time $t$. Lotka-Volterra model characterizes the predator-prey interactions over time.

\[\dfrac{\mathrm{d}X(t)}{\mathrm{d}t} = \alpha X(t) - \beta X(t) Y(t) \text{ where } \alpha , \beta > 0\] \[\dfrac{\mathrm{d}Y(t)}{\mathrm{d}t} = -\gamma Y(t) + \delta X(t) Y(t) \text{ where } \gamma, \delta > 0\]

In the above equations $\alpha$, $\delta$ denote, respectively, prey and predator birth rates. The parameter $\beta$ denotes predation rate while $\gamma$ denotes the predator death rate. Note that a Lotka-Volterra initial value problem that satisfies $X(0) = x_0 > 0$ and $Y(0) = y_0 >0$ always has a unique solution.

A phase potrait is a common visualization technique to understand the dynamics associated with a system of differential equations. Consider the Lotka-Volterra dynamical system with $\alpha = \beta = \gamma = \delta = 1$. The points in the phase space that correspond to zero rate of change of hare and lynx populations are called equilibrium points. When all parameters are set to $1$, then the system has two such equilibrium points - $(0,0)$ and $(1,1)$. The trivial equilibrium point $(0,0)$ corresponds to extinction of both species while the non-trivial point $(1,1)$ corresponds to co-existance of hare and lynx. The locus of points where $\dfrac{\mathrm{d}X}{\mathrm{d}t} = 0$ and $\dfrac{\mathrm{d}Y}{\mathrm{d}t} = 0$ denote, respectively, $X$ and $Y$-nullclines of the phase space. Note that these nullclines partition the positive quadrant of the $XY$-plane into four regions, each with a unique lynx-hare behaviour.

For instance, in the lower-left region of the phase potraits, the hare population increases (arrows point rightward) while the lynx population decreases (arrows point downward). Hares grow because predation pressure is minimal with few lynx present. Simultaneously, lynx decline due to insufficient prey to sustain their population. Trajectories move rightward and downward, pushing the system toward the x-nullcline. Additionally, notice that the hare-lynx population trajectory is a closed loop so long as the initial conditions are positive. In other words, irrespective of the initial population, the hare-lynx populations experience a cyclic behaviour while always moving anti-clockwise in the phase space ¹.

If one is given a time series data like that of the hare-lynx dynamics, there are two common approaches used to understand underlying phenomena or predict future - mechanistic models and data-driven models. Mechanistic models often model the underlying process using a deterministic system of non-linear differential equations. Mechanistic modeling requires additional information or assumptions about the underlying system in addition to the time series data. On the other hand, a data-driven requires the time-series data alone to model the underlying system. In the case of predator-prey dynamics, a mechanistic approach involves using Lotka-Volterra equations to model the behaviour of the system. The goal then will be to estimate the unknown parameters $\alpha , \beta , \gamma , \delta$ using techniques like maximum likelihood estimation. On the other hand, a data-driven approach involves using methods like Multivariate Auto-Regressive (MAR) framework to model the hare-lynx system. For more details, go through the article by Olivença et al².

Temporal Mapper

Temporal Mapper, introduced by Zhang et al.³ in 2023, is based on Mapper - a popular topological data analysis approach used to obtain combinatorial summaries of very high-dimensional data. The authors of the article are primarily interested in utilizing fMRI time series to study brain dynamics. Brain dynamics, just like the simple hare-lynx system, can be modeled using a system of non-linear differential equations. In this article, the authors develop a data analysis method to represent time series data as a directed graph such that the nodes and edges reasonably map to the underlying “brain states” and phase transitions. Authors are trying to integrate two different approaches commonly used to understand brain dynamics - mechanistic and data-driven methods. Since brain dynamics are complicated, I will use the simple Lotka-Volterra dynamics to elaborate more about the temporal mapper pipeline.

Let the matrix $\mathbf{X} \in \mathbb{R}^{n \times d}$ denote the $d$-dimensional time series vectors of interest. Let $\mathbf{D} \in \mathbb{R}^{n \times n}$ denote the pairwise distance matrix corresponding to the data in $\mathbf{X} = [ \vec{x}_1, \cdots , \vec{x}_n ]^T$. Let the data points in $\mathbf{X}$ correspond to nodes $v_1, v_2, \cdots, v_n$. The notation $v_1-v_2$ represents an undirected edge while $v_1 \rightarrow v_2$ denotes a directed edge between nodes $v_1$ and $v_2$.

Let $top(k, \vec{x}_i)$ represent the $k$ nearest neighbours of $\vec{x}_i$ (using $L_2$ norm). Then, there are two ways of defining “closeness” between points $\vec{x}_i$ and $\vec{x}_j$. Nodes corresponding to these “close” points have an edge between them. Let $G_t = (V, E)$ denotes the graph with $V = { v_1, v_2, \cdots, v_n}$. Then, one possibility is to define the edge set $E$ as follows:

\[E = \lbrace v_i-v_j : \vec{x}_i \in top(k, \vec{x}_j) \text{ OR } \vec{x}_j \in top(k, \vec{x}_i) \rbrace\]

The resulting graph $G_t$ is the “standard” $k$-nearest neighbor graph ($k$-NNG). On the other hand, if mutual “closeness” is required to include an edge, then the resulting graph $G_t$ is called the reciprocal $k$-nearest neighbour graph ($rk$-NNG).

\[E = \lbrace v_i-v_j : \vec{x}_i \in top(k, \vec{x}_j) \text{ AND } \vec{x}_j \in top(k, \vec{x}_i) \rbrace\]

Note that the authors consider the distance between temporal neighbours to be zero. On the other hand, distance between a node and itself is assumed $\infty$. In other words, for a given node $v$, the node itself cannot be one of the $k$ nearest neighbours.

Step 1: Add spatial and temporal links to generate graph $\tilde{G}$

Given the node set $V = { v_1, v_2, \cdots, v_n }$, construct an undirected $k$-NNG or $rk$-NNG $\tilde{G} = (V, \tilde{E})$ depending on how “closeness” is defined. Remove any undirected edges between consecutive temporal nodes. Now add directed edges between temporal neighbors $v_i$ and $v_{i+1}$. Therefore, the undirected edges encode spatial proximity while directed edges encode evolution in time. Notice that this step requires a user defined parameter $k$ i.e. the number of nearest neighbours to consider.

Step 2: Generate a simplified shape graph $G_s$

Using the graph $\tilde{G}$, construct a geodesic distance matrix $D_g \in \mathbb{R}^{n \times n}$ where $D_g[i,j]$ stores the geodesic distance between nodes $v_i$ and $v_j$ in the graph $\tilde{G}$. The graph $\tilde{G}$ is simplified further with a user defined parameter $d$. All nodes that are within $d$ units from each other are compressed to one node. This results in a simplified graph $G_s = (V_s, E_s)$ for which $V_s = \lbrace u_1, \cdots, u_m \rbrace$ for $m « n$. Each node $u_i$ in $G_s$ corresponds to $n_i$ nodes in $\tilde{G}$ such that $n_1 + n_2 + \cdots n_m = n$. Let $\lbrace C_1, C_2, \cdots, C_m \rbrace$ represent the set of nodes in $\tilde{G}$ that correspond to nodes in $G_s$. Then, there exists an edge between nodes $u_i$ and $u_j$ in $G_s$ if there is at least one directed temporal edge in $\tilde{G}$ between the nodes of $C_i$ and $C_j$.

\[E_s = \lbrace u_i \rightarrow u_j : \text{ } \exists \text{ nodes } x \in C_i \text{ and } y \in C_j \text{ s.t. } x \rightarrow y \in \tilde{G} \rbrace\]

Step 3: Add edge weights to $G_s$

Let $\tilde{A} \in \mathbb{R}^{n \times n}$ denote the adjacency matrix for the graph $\tilde{G}$. Then, the adjacency matrix for the graph $G_s$, denoted by $A_s \in \mathbb{R}^{m \times m}$, is defined as follows:

\[A_s[i,j] = \dfrac{1}{|C_i| \dot |C_j|}\sum_{x \in C_i} \sum_{y \in C_j} \tilde{A}[i,j]\]

Notice that the rows (and columns) of $A_s$ correspond to the nodes compressed using sets $C_1, \cdots , C_m$. A directed edge from node $u_i$ to $u_j$ (in $G_s$) exists iff $A_s[i,j] > 0$. The directed edge $u_i \rightarrow u_j$ is associated with the weight $A_s[i,j]$.

At the end of steps 1, 2 and 3, the directed graph $G_s$ that we obtain is called temporal mapper representation of the underlying dynamics.

Connection to Mapper

Step	Mapper (Classical)	Temporal Mapper (Extended)
Define Structure	Choose filter function f and partition range into overlapping intervals I₁, …, Iₛ	Construct k-nearest neighbor graph (spatial structure) + add directed temporal edges between consecutive time steps
Identify Components	Compute pre-images f⁻¹(Iⱼ) and cluster within each pre-image	Simplify graph by compressing nodes within geodesic distance d; identifies temporal clusters (states)
Generate Graph	Take nerve of clusters (connects clusters that intersect)	Add weighted directed edges based on transition frequencies between compressed nodes

References