Computational Topology and Early Warning Signals of Population Collapse

By Laura S. Storch and Sarah L. Day

Understanding and predicting the abundance of populations in the natural world is a difficult task. The abundance of any given species is often strongly influenced by its environment and interactions with other individuals and species, including humans. The nonlinear nature of these relationships adds further complication, and population fluctuations can appear chaotic or random. We must continue to study these systems and improve our ability to predict a population’s reaction to perturbation in the context of global climate change and ever-increasing human resource usage. As important system parameters such as temperature and rainfall shift, populations may undergo drastic changes, sometimes resulting in population explosion or in local and global extinction events.

The study of critical transitions, or “tipping points,” has gained considerable traction over the last decade. A “tipping point” in ecology is the transition from one dynamical regime to another due to a change in the system. For example, an increase in temperature or a decrease in habitat size could destabilize a previously stable system steady state, possibly leading to a dramatic decrease in population abundance or even a global extinction event.  Tipping points often correspond to system hysteresis, making it difficult or impossible for the population to return to its pre-transition state, even when system parameters are returned to their pre-tipping point values. Seminal papers have proposed methods for the detection of impending critical transitions in the form of increased variance or autocorrelation and slowed recovery following a perturbation [1]. This “critical slowing down”—while largely theoretical—has been recently observed in real-world systems [2]. We propose an addition to the critical transitions toolbox by explicitly considering a population’s spatial distribution and the way in which distribution patterns change during a critical transition.

We use computational homology to measure and track topological changes in population distribution patterns over time. Cubical homology, in particular, lends itself well to (cubical) coupled-patch population model output as well as digital spatial image data. Our illustrative model employs a two-dimensional lattice to represent a regular grid of spatial “patches” on which we track population abundance. Black squares represent occupied patches—spatial patches with positive abundance values—and white squares represent unoccupied patches with patch abundance values of \(0\). Figures 1a through 1c show example population distribution patterns becoming more barren. We focus on the collection of (closed) occupied patches and use cubical homology to calculate the first and second Betti numbers, \(\beta_0\) and \(\beta_1\). These Betti numbers count the number of zero- and one-dimensional holes—the number of “connected components” and “holes,” respectively—in a digital (cubical) binary image. Figure 1 lists Betti numbers for each of the three example population distribution patterns.

Figure 1. Sample population distribution patterns. Black grid squares show occupied patches and white grid squares indicate an unoccupied patch. (Patch boundaries are shown for the purpose of visualization.) We compute Betti numbers for the occupied set. The boundary of a one-dimensional “hole” must consist of boundary edges of (closed) occupied patches in order to count as a hole. Because patches are closed, diagonally-connected occupied patches form part of the same connected component. 1a. One connected component and two holes (\(\beta_0=1\) and \(\beta_1=2\)). 1b. One connected component and four holes (\(\beta_0=3\) and \(\beta_1=4\)). 1c. Three connected components and zero holes (\(\beta_0=3\) and \(\beta_1=0\)). Figure courtesy of Laura Storch and Sarah Day.

We work with a simple density-dependent growth and coupled patch dispersal model on a two-dimensional, \(N\times N\) lattice to create population distribution patterns. An agent-based model with a maximum number of  \(M\) individuals avoids issues pertaining to numerical instability. Given an \(N\times N\) matrix \(X\) (where \(X_{ij}\) is the abundance on the \(ij\)th patch) as our input population distribution, the population in each patch is reset via the density-dependent, normalized Ricker map with an Allee effect: 

\[ f(X_{ij}) = \left\{ \begin{array}{cc} rX_{ij}\textrm{exp}(1-r\frac{X_{ij}}{M}) & {\textrm{if} \:\:\:\: X_{ij}\geq \epsilon M}\\ 0 & \mbox{otherwise.} \end{array} \right. \tag1 \]

with growth rate \(r>0\) and local extinction threshold \(0 \le \epsilon \le 1\) specifying the fraction of the maximum number of individuals (\(M\)) required for reproduction. The function \(f\) models the growth phase of the population. In what follows, we set the lattice size \(N=21\), growth rate \(r=10\), and carrying capacity \(M=10^9\). We also vary the local extinction threshold.

In the following dispersal phase, a set fraction \(0\leq d \leq 1\) of the population in each patch is dispersed symmetrically to the patch’s four nearest neighbors: up, down, left, and right.  Individuals that disperse beyond the lattice’s boundary are lost. Completion of the dispersal phase yields a new output population distribution that gives the abundance on each patch. We can use this model to calculate Betti numbers for the model output population distribution patterns, and track changes in these patterns through a critical transition by varying the local extinction threshold \(\epsilon\).

Figure 2 provides a summary of average, asymptotic Betti numbers over the \((d, \epsilon)\)-parameter space. We created the figure by running 100 trials of random initial conditions for each parameter combination until dynamical transients were removed (using 1000 model iterations). We then calculated the Betti numbers on the population distribution pattern from the final iteration and averaged the results across trials for each parameter combination. Grey coloring indicates that the parameter combination led to global extinction in all trials. By choosing two example dispersal parameters, we illustrate that the dispersal parameter \(d\) is critical to how the system shifts towards global extinction.

Figure 2. Average asymptotic Betti numbers for simulations at specified values of dispersal \(d\) and local extinction parameter \(\epsilon\) using 100 random initial conditions per parameter pair. For each initial condition, we first apply 1,000 model iterations to remove dynamical transients, and then calculate \(\beta_0\) and \(\beta_1\) for the population distribution pattern produced by the final model iteration. We then average the results for each parameter pair across the 100 trials to create a picture of parameter space in terms of average asymptotic Betti numbers. Grey coloring indicates that every trial went globally extinct. Figure courtesy of Laura Storch and Sarah Day.

Using \(d=0.15\) and \(d=0.40\), we cut vertical slices in the parameter space and observe how the Betti numbers change with increasing \(\epsilon\) en route to global extinction.  Figure 3 displays the results.  Both graphs show an increase in \(\beta_1\) prior to extinction, which indicates increasing topological complexity as small regions undergo local extinction events. For the \(d=0.15\) graph, a sharp increase in \(\beta_0\) is an obvious “warning signal” that indicates population fragmentation across the lattice. In contrast, no such spike in \(\beta_0\) occurs for the \(d=0.40\) graph.

These preliminary results provide some cautionary information and avenues for further study. For example, we obtained the average Betti numbers plotted in Figures 2 and 3 from the final iterations of model runs intended to remove dynamical transients. However, a parameter in the “real world” might constantly be drifting in time, meaning that the system will not “settle” before the parameter changes again. We are currently investigating how changing \(\epsilon\) in time will affect the shape of the Betti number curves in Figure 3. 

Figure 3. Changing Betti number averages with increasing local extinction threshold \(\epsilon\) for \(d = 0.15\) and \(d = 0.40\). The \(d = 0.15\) graph provides a “warning” via a sharp spike in \(\beta_0\) before extinction, while the \(d = 0.40\) graph offers no such warning. Figure courtesy of Laura Storch and Sarah Day.

Our current results indicate that the warnings of an impending transition can vary depending on the system’s parameters. For example, a \(\beta_0\) spike (which offers a measurement of spatial fragmentation of the population) does not occur in system measurements at the higher dispersal value of \(d = 0.40\). The \(\beta_1\) curves in both plots of Figure 3 depict an increase in the number of holes prior to a global extinction event; this means that trends in \(\beta_1\) might be a more reliable indicator of a population pattern on the cusp of collapse than trends in \(\beta_0\). Further study will illustrate the robustness of this hypothesis.

We ultimately want to use this computational topology “toolkit” on spatial time series measurements of physical systems to ascertain whether researchers can use the topological information (in satellite imagery, for example) to detect early warning signs of an impending critical transition. In a rapidly changing world, computational topology can play a role in developing new and better ways to monitor and predict the dynamics of complex natural systems.

References
[1] Scheffer, M., Bascompte, J., Brock, W.A., Brovkin, V., Carpenter, S.R., Dakos, V., …, Sugihara, G. (2009). Early-warning signals for critical transitions. Nature, 461, 53-59.
[2] Van Belzen, J., van de Koppel, J., Kirwan, M.L., van der Wal, D., Herman, P.M.J., Dakos, V., …, Bouma, T.J. (2017). Vegetation recovery in tidal marshes reveals critical slowing down under increased inundation. Nat. Comm., 8, 15811.

Laura Storch received her Ph.D. in applied mathematics at the University of New Hampshire. She is currently a postdoctoral scholar at Oregon State University and was previously a postdoctoral fellow at the College of William & Mary. She is interested in critical transitions and pattern formation in ecological systems. Sarah Day is a professor and director of the Computational Applied Mathematics and Statistics program at William & Mary.  Her research interests include topological data analysis and computational methods.