Spatial Applications of Topological Data Analysis: Cities, Snowflakes, Random Structures, and Spiders Spinning Under the Influence

We now give a brief introduction to PH and tools for computing it. See [22, 49, 50] for more details. We begin by defining $k$-simplices and simplicial complexes. A $k$-simplex is $k$-dimensional polytope that is a convex hull of $k+1$ nodes. A face of a $k$-simplex is any subset (of dimension smaller than $k$) of the $k$-simplex that is itself a simplex. A simplicial complex $K$ is a set of simplices that satisfy the following requirements: (1) if $\sigma\in K$ is a $k$-simplex, then every face of $\sigma$ is in $K$; and (2) if $\sigma$ and $\tau$ are simplices in $K$, then $\sigma\cap\tau$ is a face of both $\sigma$ and $\tau$.

Given a data set $X$, we construct a sequence $X_{1}\subseteq X_{2}\subseteq\cdots\subseteq X_{l}$ of simplicial complexes of some fixed maximum dimension. We call the sequence $\{X_{i}\}$ a “filtered simplicial complex”, and we call each $X_{i}$ a “subcomplex” of the filtered simplicial complex. We equip each relation $X_{i}\subseteq X_{i+1}$ with an inclusion map. The filtered simplicial complex, along with its inclusion maps and the chain and boundary maps of each subcomplex, constitutes a “persistence complex”. The inclusion maps $X_{i}\hookrightarrow X_{j}$ induce a map $f_{i,j}\!:H_{m}(X_{i})\to H_{m}(X_{j})$ between homology groups. The map $f_{i,j}$ allows us to track an element of $H_{m}(X_{i})$ (the $m$th homology group of the subcomplex $X_{i}$) to an element of $H_{m}(X_{j})$. The $m$th homologies of the persistence complex are given by the pair

and we call them the “$m$th persistent homology” of $X$. We refer to the collection of all $m$th persistent homologies as the “persistent homology” (PH) of $X$.

Consider a generator $x\in H_{m}(X_{i})$ for some $m$ and $i$. If $x$ is not in the image of $f_{i-1,i}$, we say that $x$ is “born” at time $i$. Correspondingly, if $x\in H_{m}(X_{i})$ and $f_{i,i+1}(x)=0\in H_{m}(X_{i+1})$, we say that $x$ “dies” at time $i+1$. If for every $j\leq l$, we have that $f_{i,j}(x)\neq 0$, then we say that $x$ never dies, and we assign a death time of $\infty$ to the element $x$. For each element $x$ of the PH of $X$, there is a birth time $b_{x}$ and a death time $d_{x}$, and the collection of intervals $\{[b_{x},d_{x})\}$ is the “barcode” of $X$. Generators with longer associated half-open intervals $[b_{x},d_{x})$ are more persistent. It is traditional to construe more-persistent intervals as better indicators of a signal and construe less-persistent intervals as noise, although recent work (including our own [21, 51]) indicates that it is not always possible to interpret persistence in this way.

The collection of features in each $H_{m}(X_{i})$ describes the topological properties of the filtration $\{X_{i}\}$. Intuitively, each feature in $H_{m}$ corresponds to some $m$-dimensional void. In $H_{0}$, features are connected components; in $H_{1}$, features are loops. By considering the PHs of $\{X_{i}\}$, we can examine how the connectedness of $\{X_{i}\}$ changes for each step of the filtration in each dimension. For example, a short-persistence feature in $H_{0}$ is a connected component that appears and combines quickly with another component. PH records all features and their persistences, allowing us to take a global view of topological changes in each filtration step of $\{X_{i}\}$.

In the present paper, we use the software package Gudhi [52, 53] to compute PH of the filtered simplicial complex $\{X_{i}\}$. We construct $\{X_{i}\}$ from $X$ using two different constructions, which we developed recently in a paper on voting data [21].

We now describe a way to construct a filtered simplicial complex based on network adjacencies. We consider a network in the form of a graph $(V,E)$, with numerical data $f(v)$ associated with each node $v$. For a given filtration step $X_{i}$, let the $0$-simplices of $X_{i}$ be given by $v\in V$ such that $f(v)\leq\epsilon$ for some value $\epsilon$. For any edge $(u,v)\in E$, if $u\in X_{i}$ and $v\in X_{i}$, we add $(u,v)$ to $X_{i}$. Finally, to $X_{i}$, we add all triangles $(u,v,w)$ such that $(u,v)$, $(v,w)$, and $(u,w)$ are in $X_{i}$. We repeat this process for $X_{i+1}$, but now we use a larger value of $\epsilon$. By construction, each $X_{i}\subseteq X_{i+1}$, and we have a valid filtered simplicial complex. See Fig. 1 for an illustration of such a filtered simplicial complex.

This adjacency construction tracks topological changes in a network as it grows. The homology group $H_{m}(X_{0})$ characterizes the topology of the first filtration step. As one adds more nodes, edges, and faces to the simplicial complex, the topology changes and is recorded in $H_{m}(X_{i})$. By choosing $f$ carefully, we can control which subset of a network exists in the first filtration step, and we can also control how the network expands. For example, in [21], by attaching voting data to a network of precincts of a county, we used the adjacency construction to examine how the topology of a county changes as one considers precincts with an increasingly wide range of voting preferences.

In some of our applications, we use an alternate adjacency construction in which we associate data $g(u,v)$ to each edge $(u,v)$, instead of to the nodes. This construction differs from the one above only in that we define the function $\tilde{f}(v)=\min_{\{u:(u,v)\in E\}}g(u,v)$. We then proceed with the above adjacency construction, but we substitute $\tilde{f}$ for $f$. We recently introduced our main adjacency construction in [21], and we introduce this adaptation of it to edge-based data in the present work.

The other PH construction that we use (again see [21] for details) involves describing data as a manifold, rather than as a graph. Let $M$ denote a two-dimensional (2D) manifold, such as data in an image format. We consider the boundary $\Gamma$ of $M$ and construct a sequence

of manifolds, where at each time step, we evolve the boundary $\Gamma_{t}$ of $M_{t}$ outward according to the level-set equation. (See [54] for a thorough exposition of the level-set equation and level-set dynamics.) That is, for a manifold $M$ that is embedded in $\mathbb{R}^{2}$, we define a function $\phi(\vec{x},t)\!:\mathbb{R}^{2}\times\mathbb{R}\to\mathbb{R}$, where $\phi(\vec{x},t)$ is the signed distance function from $\vec{x}$ to $\Gamma_{t}$ at time $t$. We propagate $\Gamma_{t}$ outward at velocity $v$ using the equation

until all homological features have died. Because this evolution gives a signed distance function at each time step $t$, we take $M_{t}$ to be the set of points $\vec{x}$ such that $\phi(\vec{x},t)>0$. (This corresponds to points inside the boundary $\Gamma_{t}$.) We show an example of this evolution in Fig. 2. Throughout this paper, we use $v=1$. Different values of $v$ cause the level set to evolve faster (if $v>1$) or slower (if $v<1$), resulting in a different number of time steps (and hence a different number of filtration steps) in our evolution. However, we obtain the same homological features, although with different birth times and death times. If $v$ is sufficiently large, it is possible that all features would have the same birth and death time, such that no features would occur after the first filtration step. When $v=0$, there is no evolution.

By imposing $\{M_{i}\}$ over a triangular grid of points, as described in [21], we obtain a corresponding simplicial complex $X_{i}$ for each $M_{i}$. In Fig. 3, we give a visualization of this simplicial complex. We construct this level-set complex using a polygon whose points we choose uniformly at random from $[0,1]\times[0,1]$ as an initial synthetic image. Because the level-set equation (2) evolves continually outward, we automatically satisfy that condition that $X_{i}\subseteq X_{i+1}$, so $\{X_{i}\}$ is a filtered simplicial complex. Our implementation of the level-set method works with any black-and-white image (or any image that one formulate as a piecewise-constant function $\mathbb{R}^{2}\to\{0,1\}$). We expect our level-set approach to capture information about $H_{0}$ and $H_{1}$ for any such image. The level-set approach also captures geometric information, which can be useful for some applications; however, this may make it difficult to capture information about holes that are visually irregular. Throughout this paper, we compare images that have roughly the same resolutions, where we take the image resolution from raw image data. Because image size should primarily affect the computation time of our level-set approach — but not the order in which features appear and disappear as an image evolves — we expect that it is possible to adapt our level-set construction when comparing images of different resolutions. Possible approaches for such an adaptation include normalizing image sizes or adjusting the resolution of the triangular grid that one uses for each image.

In this subsection, we discuss applications of our adjacency PH construction to a dynamical process on synthetic networks in which space plays an important role. For each network $(V,E)$, we run the Watts threshold model (WTM) [55] on it. Given a graph, we select a fraction $\rho_{0}=0.05$ of its nodes uniformly at random to be “infected” at time $0$. At each time step, we then we compute the fraction of each node’s neighbors that are infected. (That is, we synchronously update the states of the nodes [56].) If the fraction of a node’s neighbors that are infected meets or exceeds a threshold (in our case, the threshold is $\varpi=0.18$ for all nodes), the node becomes infected. We take this implementation of the WTM to be the generator of a function $f\!:V\to\mathbb{N}$ ¹¹1We use the convention that $\mathbb{N}$ includes $0$., where $f(v)$ is the time at which node $v$ becomes infected. We say that infected nodes are in the set $I$. If $v$ never becomes infected, we set $f(V)=\max_{v\in I}f(v)+1$, so that we eventually add all nodes to a filtered simplicial complex. The resulting filtered simplicial complex consists of the subgraphs that are generated by $I$ at each time step. See [58, 59, 60] for studies of the WTM on spatial networks.

We use parameter values of $\rho_{0}=0.05$ and $\varpi=0.18$ throughout this section. We expect changes in $\rho$ and $\varpi$ to affect the birth times and death times of features in a filtered simplicial complex. Using a different value of $\rho_{0}$ entails considering a different fraction of initially infected notes. Therefore, a larger value of $\rho_{0}$ yields a larger simplicial complex at the first filtration step, and smaller value of $\rho_{0}$ yields a smaller simplicial complex. Using a larger value of $\varpi$ results in fewer nodes becoming infected at each time step, and it thus takes more filtration steps before the simplicial complex stops growing. Because our underlying graph is the same for any choice of values of $\rho$ and $\varpi$, we do not expect changes in the homology of the last filtration step, unless $\rho$ or $\varpi$ are sufficiently small such that some nodes in a graph never become infected. However, one can certainly obtain a different PH for different values of $\rho$ or $\varpi$, as nodes and edges can join the filtered simplicial complex at different times and (more importantly) in different orders.

We examine topological changes in the infected subgraph of three different types of synthetic networks (see Fig. 4). We first examine random geometric graphics (RGGs) [61]. For each instance of an RGG, we pick $100$ nodes uniformly at random from the unit square. If the Euclidean distance between two nodes is less than or equal to $0.1$, we add an edge between them [see Fig. 4(a)]. Our second type of synthetic network is a square lattice with $100$ nodes. We arrange the $100$ nodes in a $10\times 10$ grid on the unit square, and we then connect the nodes along the grid lines [see Fig. 4(b)]. Our third type of synthetic network is a Watts–Strogatz (WS) small-world network [62, 63]. We begin with a ring of $100$ nodes, and we then connect each node to its $k=2$ nearest neighbors on each side. We then rewire each edge uniformly at random with a probability of $p=0.1$ using the implementation of the WS model in NetworkX [64]. In this version of the WS model, one removes each rewired edge before replacing it with a new edge. We show an example of a WS graph in Fig.. 4(c).

For each type of synthetic network, we consider $100$ instances, which we generate using NetworkX. For the RGG and WS networks, each instance is a different graph; the square lattice network is deterministic. For all three types of networks, each instance has a different initial set of infected nodes. We show visualizations of each of these types of networks (with WTM dynamics on it) in Fig. 4.

Our adjacency construction begins by selecting the initially infected nodes and the edges between them of a network to create an infected subgraph that we call an “infection network”. As an infection spreads, we add more nodes and edges to the infection network until eventually we have added all nodes and edges to it.

Examining the PHs of the RGGs (see Fig. 5), we see for our parameter values that an infection network tends to have several connected components, resulting in a large number of features in $H_{0}$. However, because of the spreading behavior of the WTM, new nodes can become infected only via their infected neighbors. Because features in $H_{0}$ record connected components of a graph, new infected nodes join existing connected components. Therefore, features can only be born at time $0$ or in the last step, which is when we add all remaining uninfected nodes to our filtered simplicial complex. By contrast, features in $H_{1}$ are relatively rare, as most cycles that occur in an RGG are filled because of the uniform probability distribution of the node locations.

For a square lattice network (see Fig. 6 for a PD of the WTM on such a network), we first note that there is only a single infinite-length feature in $H_{0}$, as the final infection network necessarily consists of a single connected component. Consequently, $H_{0}$ consists of a set of features that are born at time $0$ and eventually merge (and therefore die), resulting in a single infinite-length feature. Additionally, there are a constant number (81, to be precise) of features in $H_{1}$, because when we construct a simplicial complex, every grid cell of the lattice is a feature in $H_{1}$ at the last filtration step. However, these features can be born at a variety of times, as the filtration does not include every lattice cell until every node of the graph has entered the filtration.

From Fig. 7, we see that a WS small-world network also eventually has an infection network that consists of a single connected component. However, the WS networks consistently have more features in $H_{1}$ than the RGG networks, because the former’s (non-geometric) shortcut edges usually result in splitting an existing cycle (and hence a feature in $H_{1}$) into two cycles.

We summarize our observations about the various synthetic networks in Table 1, in which we give the means and standard deviations of the number of features during the temporal evolution of the WTM in each type of synthetic graph. Our counts include features that appear at any time during the WTM dynamics.

The field of urban analytics has grown rapidly in the last several years [1, 65, 11], Increasingly powerful computational tools have allowed researchers to characterize cities in terms of their street networks [66], and a variety of approaches from network analysis have been applied to the study of urban street networks [67, 68, 69, 8, 70, 71, 72]. In the present subsection, we use city street networks as base manifolds in our level-set construction, and we thereby characterize cities based on their PHs. We use these PHs to compare city morphologies both within a single city and across a variety of cities.

We use our level-set construction to obtain topological descriptors in the form of persistence for city street networks. We then use these city persistences to compare (1) different regions of the same city and (2) different cities to each other. We obtain all of our city street networks with the software package OSMnx [73] using latitude–longitude coordinates and taking a $1$ km block that is centered at specified coordinates. In each example, we indicate how we choose these coordinates.

The first filtration step of a filtered simplicial complex that results from our PH construction consists only of the streets in a network. As we increase the filtration time, we slowly add city blocks to the complex, and the topology changes as those blocks are filled in. More regular city blocks are more likely to be filled in without creating any new homological features, and larger blocks take longer to be filled in. Our construction is thereby able to capture information about the size and regularity of city blocks. The existence of dead ends tends to lead to the “pinching” of blocks into multiple homological features — as dead ends expand, they lengthen and eventually meet with nearby streets, cutting through blocks in the process — so our approach also yields information about dead ends.

As a second application that uses empirical data, we consider snowflake crystals [84]. We start with twelve different images (from [84]) of snowflakes with different crystalline structures. (See Fig. 20 in Appendix A.) Using the GNU Image Manipulation Program [85], we threshold these grayscale images (using a thresholding setting of 205) to create black-and-white images. From the black-and-white images, we compute level-set complexes and PHs, and we then perform average-linkage hierarchical clustering on the PDs to produce the dendrogram in Fig. 17.

The images of snowflakes consist of edges (the black lines in our images) and cells (the white spaces that are bounded by the edges). We refer to the outer areas that extend from the center of the snowflakes as their “points”. The twelve snowflakes have fairly regular crystalline structures, so our computation of PH predominantly records information about the distribution of cell sizes in a snowflake. The inherent hexagonal nature of snowflakes and the regularity of their crystalline structures largely overwhelms our ability to use PH to glean information about their spatial relationships and irregularities. Constructing a simplicial complex that is better suited to capturing information about images with a large number of regular structures may yield better results.

Examining the clusters (see Fig. 17) reveals that snowflake A and snowflake B each reside in their own cluster, and the remaining snowflakes split into two clusters. Snowflake A’s PD [see Fig. 16(a)] is dominated by a feature in $H_{1}$ with an early birth and a late death. (See the blue square in the top-left corner of the diagram. This arises from the large feature that is formed by the bold ring near the center of the snowflake. None of the other snowflakes have a bold central ring. More generally, we observe few features in the PD for snowflake A. By contrast, snowflake B’s features are largely close to the diagonal [see Fig. 16(b)] because the initial manifold of the snowflake does not have large holes. Notably, we do not observe any points in the top-left region of its PD. The cell sizes in snowflake B are smaller than those in most of the other snowflakes, and even its central ring structure includes a large number of small holes. The remaining snowflakes either have more large holes than snowflake B, or they do not have a bold central ring like the one in snowflake A. We also note the PD of snowflake B is much closer than that of snowflake A to those of the other snowflakes.

Our final application is to the topology of spiderwebs. In 1948, Peter Witt began research on the effects of drugs on spiders to test whether garden spiders would shift their web-building hours if they were administered drugs. Witt found that drugs affect the size and shape of the webs that are produced by spiders ²²2Interestingly, whiskey itself produces webs [90].. He also found that higher doses of most drugs (e.g., 100 $\mu$g per spider, as opposed to 10 $\mu$g per spider) tend to lead to larger changes in the shapes of webs, including yielding more irregular webs. Witt eventually published more than $100$ papers and several books on the behavior of spiders and spider webs. For more information on his experiments with psychotropic substances and spiders, see his 1971 review article [87]. In a 1995 technical briefing [88], NASA (which was inspired by Witt’s research) proposed that spiders who were administered more toxic substances produce webs that are more deformed (in comparison to a web that is spun by a drug-free spider) than less toxic substances. Additionally, using techniques from statistical crystallography, they concluded that spiders fail to complete more sides of their webs when they are under the influence of more toxic substances.

In our case study of PH in spiderwebs, we use five images from the NASA technical briefing [88] and two images from Witt [87] of webs that were spun by spiders under the influence of a variety of psychotropic substances, threshold grayscale images to turn them into black-and-white images (using a thresholding setting of 205 in the GNU Image Manipulation Program), apply our level-set construction to compute PH, and perform average-linkage hierarchical clustering to yield the dendrogram in Fig. 18. We show the images of the spiderwebs and their associated PDs in Fig. 19.

Our classification places the drug-free spider into its own cluster. The spiderweb of the drug-free spider is characterized by a clear central hole, threads radiating outward at approximately even intervals, and completed rings of threads that surround the center. We place the webs that were spun by spiders under the influence of marijuana, peyote, and LSD into the same cluster. In these webs, there is a clearly identifiable center; and most of the radial threads are evenly-spaced, straight, and radiate outward directly from the center. However, for the webs in this cluster, rings of threads are either difficult to see or are incomplete. The final cluster consists of webs that were spun by spiders under the influence of chloral hydrate, caffeine, and speed. In the caffeinated spider’s web, one cannot even clearly identify a center ³³3The web that was produced by the caffeinated spider is always fun to point out when giving presentations.. One can locate a center in the webs of the spiders that were under the influence of speed or chloral hydrate (a sedative that is used in sleeping pills), but many of the radial threads do not join the center and some of the radial threads are not straight. Almost no complete rings of thread are visible in any of the three webs in this cluster.