Persistent homology of unweighted complex networks via discrete Morse theory.

Topological data analysis can reveal higher-order structure beyond pairwise connections between vertices in complex networks. We present a new method based on discrete Morse theory to study topological properties of unweighted and undirected networks using persistent homology. Leveraging on the features of discrete Morse theory, our method not only captures the topology of the clique complex of such graphs via the concept of critical simplices, but also achieves close to the theoretical minimum number of critical simplices in several analyzed model and real networks. This leads to a reduced filtration scheme based on the subsequence of the corresponding critical weights, thereby leading to a significant increase in computational efficiency. We have employed our filtration scheme to explore the persistent homology of several model and real-world networks. In particular, we show that our method can detect differences in the higher-order structure of networks, and the corresponding persistence diagrams can be used to distinguish between different model networks. In summary, our method based on discrete Morse theory further increases the applicability of persistent homology to investigate the global topology of complex networks.

Introduction

In recent years, the field of topological data analysis (TDA) has rapidly grown to provide a set of powerful tools to analyze various important features of data[1]. In this context, persistent homology has played a key role in bringing TDA to the fore of modern data analysis. It not only gives a way to visualize data efficiently, but also to extract relevant information from both structured and unstructured datasets. This crucial aspect has been used effectively in various applications from astrophysics (e.g., determination of inter-galactic filament structures[2]) to imaging analysis (e.g., feature detection in 3D gray-scale images[3]) to biology (e.g., detection of breast cancer type with high survival rates[4]). Informally, the essence of the theory is its power to extract the shape of data, as well as infer higher-order correlations between various parts of the data at hand which are missed by other classical techniques[1]. The basic mathematical theory used in this subject is that of algebraic topology, and in particular the study of homology, developed by the French mathematician Henri Poincaré at the turn of the 20th century. The origins of persistent homology lie in the ideas of Morse theory[5], which gives a powerful tool to detect the topological features of a given space through the computation of homology using real-valued functions on the space. We refer the reader to the survey article[6] for further details.

On the other hand, the discretized version of Morse theory developed by Robin Forman[7-9], gives a way to characterize the homology group of a simplicial complex in terms of a real-valued function with certain properties, known as a discrete Morse function. Examples of such simplicial complexes associated with discrete spaces are the Vietoris-Rips complex corresponding to a discrete metric space, or the clique complex of a graph. Forman[8,9] showed that given such a function, the so-called critical simplices completely determine the Euler characteristic of the space, which is a fundamental topological invariant.

The study of complex networks in the last few decades has also significantly raised our ability to understand various kinds of interactions arising in both natural and artificial realms[10-13]. Understanding how different parts of networks behave and influence each other is therefore an important problem[10-13]. However, for large networks, detecting higher-order structures remains a difficult task[14]. Moreover, recent studies[15-17] indicate that these higher-order correlations are not captured by usual network measures such as clustering coefficients. While a graph representation captures binary relationships among vertices of a network, simplicial complexes also reflect higher-order relationships in a complex network[15,16,18-25]. In this context, persistent homology has been employed to explore the topological properties of complex networks[18-21,23,26]. In this work, we present a systematic method to study the persistent homology of unweighted and undirected graphs or networks.

Previous work has investigated the persistent homology of weighted and undirected networks by creating a filtration of the clique complexes corresponding to threshold graphs obtained via decreasing sequence of edge weights[20,23]. However, the lack of edge weights in unweighted networks does not permit a filtration based on threshold graphs[20,23]. Thus for unweighted networks, Horak et al.[19] propose a filtration scheme based on the dimension of the simplices in the clique complex corresponding to the unweighted network. Horak et al.[19] do not assign weights to vertices, edges or higher-dimensional simplices in the clique complex corresponding to an unweighted graph. An unexplored filtration scheme involves transforming an unweighted network into a weighted network by assigning edge weights based on some network property, such as edge betweenness centrality[27,28] or discrete edge curvature[29,30], and then employing the filtration scheme based on threshold graphs[20,23]. As an alternative, we here use discrete Morse theory[7-9] to create a filtration scheme for unweighted networks by assigning weights to vertices, edges, triangles and higher-dimensional simplices in the clique complex of the graph. In our method, the weight of a simplex is chosen such that it reflects the degree of the vertices which constitute the simplex while simultaneously satisfying the conditions for the weighing function to be a discrete Morse function. Moreover, as explained in the Results section, an equally important intuition behind the choice of these weights is based on the goal of reducing the number of so-called critical simplices.

In the context of TDA, classical Morse theory, which involves smooth functions defined on topological spaces that admit a smooth structure has been used to compute persistent homology, e.g. in statistical topology[31], astrophysics[2]. Since the clique complex of a weighted or unweighted graph does not permit a smooth structure in general, applying classical Morse theory is not possible. However, discrete Morse theory[7-9] provides an efficient way of computing persistent homology. A discrete Morse function not only captures higher-order topological information of the underlying space, a “preprocessing” with respect to a suitable discrete Morse function leads to significant simplification of the topological structure. This makes computation of persistent homology groups or homology groups of filtered complexes much more efficient, see e.g.[32]. This is especially useful for large datasets where computationally efficient methods are key to compute their persistent homology.

We have applied discrete Morse theory to compute persistent homology of unweighted simple graphs. This is done by using the values given by the discrete Morse function to pass from an unweighted graph to a weighted simplicial complex (Fig. 1). This transformation automatically produces a filtration that is needed for the computation of persistent homology, through the so-called level subcomplexes associated with weights of critical simplices (See Theory section and Fig. 2). Moreover, this filtration is consistent with the topology of the underlying space and reveals finer topological features than the dimensional filtration scheme used in Horak et al.[19] The combination of these techniques have been used[3,33] with applications for image processing. However, to the best of our knowledge, this method has not been used for studying persistent homology in unweighted complex networks to date. Discrete Morse theory gives a theoretical lower bound on the number of critical simplices which can be attained by an optimal choice of the function on a simplicial complex. Interestingly, our method achieves close to the theoretical minimum number of critical simplices in several model and real networks analyzed here (See Results section). Furthermore, our algorithm for computing the discrete Morse function is easy to implement for complex networks.

Figure 1

An illustration of the construction of a discrete Morse function f on a clique complex K corresponding to an unweighted and undirected graph G using our algorithm. (a) A simple example of an unweighted and undirected graph G containing 9 vertices and 11 edges. (b) The clique simplicial complex K corresponding to the simple graph G shown in (a). The clique complex K consists of 9 vertices or 0-simplices, 11 edges or 1-simplices and 2 triangles or 2-simplices. The figure also displays the orientation of the 1- and 2-simplices using arrows. (c) Generation of a discrete Morse function f on the clique complex K shown in (b) using our algorithm. The figure lists the state of the Flag variable in algorithm 1 and IsCritical variable in algorithm 2 (See SI Appendix) for each simplex in K. In this example, the clique complex has 4 critical simplices and their respective critical weights correspond to the filtration steps. The figure also lists the FiltrationWeight for each simplex in K obtained using algorithm 3 (See SI Appendix).

Figure 2

Filtration based on the entire sequence of weights satisfying discrete Morse function is equivalent to filtration based only on the subsequence of critical weights in terms of persistent homology. (a) Filtration of the network shown in Fig. 1 based on weights of the 4 critical simplices. There is a 0-hole (or connected component) that persists across the 4 stages of the filtration while another 0-hole is born at stage 2 on addition of critical vertex v2 but dies at the stage 3 which corresponds to the weight of the critical edge . Moreover, a 1-hole is born at the stage 4 on addition of the critical edge . (b) Five intermediate stages during the filtration between critical weights 1.1 (stage 2) and 2.35 (stage 3). (c) Four intermediate stages during the filtration between critical weights 2.35 (stage 3) and 3.48 (stage 4). It is seen that the homology of the clique complex remains unchanged during the intermediate stages of the filtration whereby the birth and death of holes occur only at stages which correspond to critical weights.

Our results underline the potence of persistent homology to detect inherent topological features of networks which are not directly captured by homology alone. For instance, the p-Betti numbers of the clique complexes corresponding to small-world[10] and scale-free[11] networks with similar size and average degree, respectively, are of comparable magnitude and thus, homology reveals no deep insight into the differences between the topological features of these two model networks. On the other hand, our observations on the persistent homology of these two networks indicate a clear demarkation with respect to the evolution of topological features in the clique complexes corresponding to these model networks during the filtration process. This dissimilarity in the evolution of topological characteristics that resonates across dimensions and the average degree of the underlying network, indicates an inherent disparity in the persistent homology of small-world and scale-free graphs. In addition to unravelling higher-order relationships in networks, this ability to capture inherent topological differences between two dissimilar networks thus motivates the application of our methods to study the persistent homology of real-world networks.

The remainder of the paper is organized as follows. We begin with a Theory section which gives a brief overview of concepts in persistent homology and discrete Morse theory. We then proceed to describe the model networks and real-world networks that have been studied in this work in the Network datasets section. In the subsequent section on Results and Discussion, we present our algorithm to construct a discrete Morse function on a simplical complex associated with a network. In the same section, we present our results for model networks and real-world networks. The final section on Conclusions gives a summary and outlook of our findings. In Supplementary Information (SI) Appendix, we give a brief review of the mathematical theory of homology groups of a simplicial complex. In SI Appendix, we also provide a rigorous proof of concept for our algorithm to construct a discrete Morse function. We then follow up with two algorithms both of which illustrate key procedures that are essential to construct the filtration of a simplical complex associated with the investigated networks. In SI Appendix, time complexity for our algorithms as well as some theoretical results on stability of persistent homology and persistence diagrams are also given.

Theory

Graphs and simplicial complexes

Consider a finite simple graph having vertex set and the edge set . Note that a simple graph does not contain self-loops or multi-edges[34]. Such a simple graph G can be viewed as a clique complex K[35]. A clique simplicial complex K is a collection of simplices where a p-dimensional simplex (or p-simplex) in K is a set of vertices that form a complete subgraph. In other words, vertices correspond to 0-simplices, edges to 1-simplices, and triangles to 2-simplices in the clique complex of a graph. Note that the dimension p of simplices contained in K is restricted to the range 0 to in the graph G. The dimension d of the clique complex K is given by the maximum dimension of its constituent simplices. A face γ of a p-simplex α is a subset of α with cardinality less than . Note that by definition, a face γ of a p-simplex α is a l-simplex where and this relationship is denoted as . Formally, the clique complex K corresponding to the simple graph G satisfies the following condition which defines an abstract simplicial complex, namely, K is a collection of non-empty finite sets or simplices such that if α is an element (simplex) of K then so is every non-empty subset of α. For additional details, the interested reader is referred to standard text in algebraic topology[36].

Figure 1 displays an example of the correspondence between a simple graph and its clique complex. The ordering of the vertex set of a p-simplex α determines its orientation. Moreover, two orderings of the vertex set of α are considered to be equivalent if and only if they differ by an even permutation. If the dimension of a p-simplex is greater than 1, then all possible orderings of its vertex set fall under two equivalence classes, with each class being assigned an orientation[36]. An exception is the 0-simplex with one vertex which has exactly one equivalence class and orientation. An oriented p-simplex α specifies the orientation of its vertices and is represented by [36]. In Fig. 1, the oriented 2-simplices and have opposite orientations, i.e., .

Persistent homology of a simplicial complex

In SI Appendix, we give a brief review of the mathematical theory of homology groups of a simplicial complex. In particular, we define p-chain group, p-boundary operator, p-boundary, p-cycle, p-hole, p-homology group and p-Betti number.

A subset K of a simplicial complex K is called a subcomplex of K if K by itself is an abstract simplicial complex. Then a filtration of a simplicial complex K is defined as a nested sequence of subcomplexes K of K where:

Note that each subcomplex K has an associated index i in the filtration. Moreover, each subcomplex K in the filtration has corresponding p-chain complexes , p-boundary operators , p-boundaries and -cycles .

The j-persistent p-homology group of K denoted as is defined as:

In the above equation, is the subgroup of which constitutes the p-boundaries of the subcomplex . The j-persistent p-Betti number of K denoted as is defined as:

An intuitive explanation of the above definitions of the j-persistent p-homology group and the corresponding Betti number is as follows. A p-hole of the subcomplex K can potentially become the boundary of a -chain of a later subcomplex with , and thus, no longer constitute a p-hole of K. The j-persistent p-Betti number of K represents the number of p-holes at the filtration index i that persist at the filtration index j + i. Therefore, each p-hole that appears across the filtration has a unique index that corresponds to its birth and death, and the persistence of such a p-hole can thus be characterized by its corresponding birth and death indices. Studying persistent homology allows us to quantify the longevity of such p-holes during filtration, and thus, measures the importance of these topological features which appear and disappear across the filtration.

Discrete Morse theory

Recalling from the preceeding section, to study the persistent homology corresponding to the clique complex K of a simple graph G, the primary requirement is a filtration of K. We here present a systematic method to study the persistent homology of unweighted and undirected networks by utilizing a refined filtration of the clique complex K based on discrete Morse theory[7-9]. It is important to note that the order in which the simplices of the clique complex K are added during the filtration affects the evolution of topological features which are observed by studying the persistent homology. Our proposed scheme which is based on discrete Morse theory developed by Forman[7-9] tackles this by assigning weights to 0-simplices (vertices), 1-simplices (edges), 2-simplices (triangles) and higher-dimensional simplices appearing in the clique complex corresponding to an unweighted and undirected network. Assigning weights to higher-dimensional simplices captures important higher-order correlations in addition to edges or 1-simplices. Moreover, we leverage the following important features of the framework of discrete Morse theory in our new scheme. Firstly, the framework enables assignment of weights to p-simplices which are concordant with weights of -simplices. Secondly, it captures the topology of a simplicial complex via the concept of critical simplices described below. Most importantly, the framework provides a natural way to create a filtration scheme to study persistent homology based upon the weights of the aforementioned critical simplices as will be described below.

We next provide the fundamental definitions in discrete Morse theory[8,9]. We remark that a p-dimensional simplex α in a simplicial complex K is denoted by . Also, if a p-simplex α in K is a face of a -simplex β in K then this is represented as in the sequel. Given a function , for each simplex , two sets and are defined as follows:

Simply stated, the set contains any -simplex β of which α is a face and the function value on β is less than or equal to the function value on α. The set contains any -simplex γ which is a face of α and the function value on α is less than or equal to the function value on γ. A function is a discrete Morse function[8,9] if and only if for each simplex :

Given a discrete Morse function f on the simplicial complex K, a simplex is critical[8,9] if and only if:

Simply stated, a p-simplex is critical if the following conditions are simultaneously satisfied. The first condition being that if β is any -simplex in K of which α is a face, then . The second condition being that if γ is any -simplex in K which is a face of α then . The concept of critical simplices in discrete Morse theory is in spirit a discrete analogue to the concept of critical points in classical Morse theory wherein the critical points corresponding to a smooth real valued function on are the points where the gradient of the function vanishes.

We remark that once a discrete Morse function f on a simplicial complex K is fixed, the sets and are denoted by U and V, respectively, to simplify the notation. A simple example for a discrete Morse function on a simplicial complex K is the dimension function used in Horak et al.[19]. The value of the dimension function on a given simplex α is the dimension of the simplex α. By a direct consequence of the definitions presented above in this section, for every simplex , the sets U and V corresponding to the dimension function are empty. Thus, the dimension function is indeed a discrete Morse function and every simplex in K is critical. In the results section, we present our new scheme and algorithm 1 to assign a discrete Morse function f to a clique complex K of an unweighted graph G.

We next describe the filtration of the clique simplicial complex K based on the discrete Morse function f. Given a discrete Morse function f on a simplicial complex K and a real number r, a level subcomplex K(r) is defined[8,9] as follows:

Simply stated, K(r) contains all simplices β in K with the value of the discrete Morse function or assigned weight along with any face α of β. Note that a face α of β is included in even if the discrete Morse function or assigned weight to a face α is greater than r.

Let denote the entire set of values assigned to simplices in K using the discrete Morse function f. Then, let denote the finite increasing sequence of the unique values in the set associated with the finite simplicial complex considered here. We now have a sequence of inclusions of level subcomplexes corresponding to this increasing sequence as follows:

This nested sequence gives a filtration of the simplicial complex K which enables the study of persistent homology in the context of unweighted networks.

According to Lemma 2.6 by Forman[9], if there are no critical simplices α with , then K(b) is homotopy equivalent to K(a).

The implications of this Lemma are as follows. Let denote the set of values assigned to critical simplices σ in K by the discrete Morse function f and denote the increasing sequence of the unique values in . We refer to the function values assigned to the critical simplices σ in K as critical weights. Note that the set defined for critical simplices is a subset of defined for all simplices in K and the increasing sequence is a subsequence of with . The above definition implies that there are no critical simplices α with . As homology is invariant under homotopy equivalence, Forman’s Lemma 2.6 gives us that for any x and y belonging to the real number interval , the homology groups of K(x) and K(y) are isomorphic. Thus, in order to observe the changes in homology as the filtration proceeds, it suffices to study the persistent homology of a filtration which corresponds to the subsequence of , where , and this results in a potential decrease in the required number of filtration steps. The new filtration sequence can be represented as:

Note that each simplex α in the clique complex K is first introduced as part of certain level subcomplex in the above nested filtration sequence. Therefore, each simplex α in K can be associated with a unique weight referred to as the filtration weight of α. In SI Appendix, we present algorithms 2 and 3 which depict the procedure to compute the filtration weights of simplices in the clique complex K of a graph G. In SI Appendix, we also give a sufficient condition for two discrete Morse functions to induce the same filtration and thus the persistent homology groups.

Using an example network in Fig. 2, we also show that the persistent homology observed using the filtration based on the entire sequence of weights satisfying discrete Morse function is equivalent to that observed using the filtration based on the subsequence of critical weights.

Let m represent the number of critical p-simplices in a simplicial complex K and let β denote the p-Betti number of K. Then Theorem 2.11 by Forman[9] can be stated as follows.

(i) For each p = 0, 1, 2, …, d (where d is the dimension of K), m ≥ β.

(ii) m0 − m1 + m2 −  + (−1)m = β0 − β1 + β2 −  + (−1)β.

In other words, the above theorem gives a lower bound of the number of critical p-simplices m for each dimension p as the p-Betti number β of K. In results section, we present our algorithm 1 to assign weights satisfying discrete Morse function to simplices in the clique complex K of a graph G. Our choice of the function in algorithm 1 to assign weights to simplices in the clique complex K tries to minimize the number of critical simplices (which has a lower bound given by Forman’s Theorem 2.11[9]), and thus, reduces the number of filtration steps required to compute the persistent homology without loss of information. In the results section, we will show that our algorithm achieves near-optimal number of critical weights in clique complexes corresponding to many model and real networks analyzed here.

Comparing persistence diagrams

Given a discrete Morse function f and its associated filtration of the clique complex K of a graph G (Eq. 10), each p-hole has a critical weight which corresponds to its birth index and which corresponds to its death index, with .

Persistence diagram D(f) for a d-dimensional simplicial complex K is the collection of points in whose first and second coordinates, x and y, respectively, correspond to the birth weight and death weight of a p-hole where [37]. Since two different holes can have the same birth and death weights, each point in the persistence diagram has a corresponding multiplicity, we refer the reader to the SI Appendix for more details. Thus, the persistence diagram is a multiset of points in . The persistence of a p-hole which has birth and death weights, and , respectively, is defined as . Thus, the persistence diagram for a clique complex K corresponding to a graph G is a compact representation of the persistent homology of a network.

Given two persistence diagrams X and Y (which may correspond to two different networks), the -Wasserstein distance between X and Y, also known as the bottleneck distance[37], is defined as follows:

Similarly, given two persistence diagrams X and Y, the q-Wasserstein distance[38] between X and Y is defined as follows:

In the above equations, ranges over all bijective maps from X to Y, and given , is the norm. In this work, we use Dionysus 2 package (http://www.mrzv.org/software/dionysus2/) to compute the Wasserstein distance between two persistence diagrams corresponding to two different model networks (See results section). Note that it is not generally true that two persistence diagrams X and Y have the same number of off-diagonal points, i.e., features with non-zero persistence, and we refer the readers to Kerber et al.[38] for details on circumventing this issue and further information regarding how the computation of the Wasserstein distance is reduced to a bipartite graph matching problem in Dionysus 2 package. We remark that the bottleneck distance between two persistence diagrams which are subsets of the unit square is in the range 0 to 1.

Stability of persistence diagrams with respect to small changes in the discrete Morse function is a key property of persistent homology. The first such stability theorem was given by Cohen-Steiner et al.[37], and later refined by Chazal et al.[39]. In the SI Appendix, using results from Chazal et al.[39], we give a stability result for persistence diagrams of discrete Morse functions with respect to the bottleneck distance.

Network Datasets

Model networks

We have investigated the following models of unweighted and undirected networks, namely, the Erdös-Rényi (ER)[40], the Watts-Strogatz (WS)[10], the Barabási-Albert (BA)[11] and the Hyperbolic Graph Generator (HGG)[41]. The ER model[40] is characterized by the property that the probability p of the existence of each possible edge between any two vertices among the n vertices in the graph G is constant. The existence of edges in the ER model are independent of each other, and thus, the model produces random graphs with average vertex degree . The WS model[10] produces small-world graphs as follows. The WS model starts with an initial regular graph with n vertices where each vertex is connected to its k nearest neighbours. Next, the endpoint of each edge in the initial regular graph of the WS model is randomly chosen for rewiring based on a fixed rewiring probability p and is rewired to another vertex in the graph which is chosen with uniform probability. The BA model[11] produces scale-free graphs which are characterized by a degree distribution that follows a power law decay. The BA model utilizes a preferential attachment scheme to produce scale-free graphs. The BA model generates an initial graph of m0 vertices, and then, at each successive iteration a new vertex is added with edges to m already existing vertices which are chosen with probability proportional to their degree at that particular iteration. The iterations in the BA model cease when the graph has attained the requisite number n of vertices. The HGG model[41,42] produces a random graph of n vertices by initially fixing n vertices to n points on a hyperbolic disk. In the HGG model, the probability of existence of an edge between two vertices is proportional to the hyperbolic distance between the two points on the hyperbolic disk that correspond to these two vertices. By tuning the input parameter γ, the HGG model can produce either a hyperbolic or a spherical random graph[41,42]. Specifically, the HGG model produces hyperbolic random graphs for whereas spherical random graphs for .

Real networks

We have also studied seven real-world networks which are represented as unweighted and undirected graphs. We have considered two biological networks, namely, the Yeast protein interaction network[43] with 1870 vertices and 2277 edges, and the Human protein interaction network[44] with 3133 vertices and 6726 edges. In both biological networks, each vertex represents a protein and an edge represents an interaction between the two proteins. We have considered two infrastructure networks, namely, the US Power Grid network[45] and the Euro road network[46]. In the US Power Grid network, the 4941 vertices represent the generators, transformers and substations in the Western states of USA and the 6594 edges represent power links between them. The 1174 vertices of the Euro road network correspond to cities in Europe and the 1417 edges correspond to roads linking the cities. We have also studied the Email network[47] of the University of Rovira i Virgili with 1133 vertices representing users and 5451 edges, each representing the existence of at least one Email communication between the two users corresponding to the vertices anchoring the edge. We have also studied the Route views network[45] which has 6474 autonomous systems as vertices and 13895 edges representing communication between the systems that are represented as vertices. We have considered a social network, the Hamsterster friendship network[48], containing 1858 vertices which represent the users and 12534 edges which represent friendships between the users. Note that we omit self-loops while constructing the clique complex K corresponding to the undirected graph G of a real-world network.

Results and Discussion

Algorithm to construct discrete Morse function on a simplicial complex

From an unweighted and undirected graph with vertex set and edge set , it is straightforward to construct a clique simplicial complex K with dimension d (See Theory section). Figure 1 shows the construction of a clique complex starting from an example network. Given a simplicial complex K, its dimension d and a non-negative real-valued function g on the 0-simplices of K, the algorithm 1 assigns weights to any simplex in K, producing a discrete Morse function f defined in Eq. 6. In the pseudocode of the algorithm 1, lines 2–6 initialize a variable Flag[α] for every simplex α in clique complex K with the value 0. We remark that the variable Flag[α] associated with a simplex α in K serves as a counter for the size of the set U defined in Eq. 4. Lines 7–9 assign weights to every 0-simplex in K based on the input non-negative function g. Lines 10–24 assign weights to 1- or higher-dimensional simplices in K in a manner which is consistent with the definition in Eq. 6 of a discrete Morse function. In summary, algorithm 1 outputs a discrete Morse function f on K, and in SI Appendix, we present a rigorous proof for the following theorem which states the same.

Theorem.

Algorithm 1 produces a discrete Morse function f on any simplicial complex K of finite dimension d.

Algorithm to construct a discrete Morse function on a d-dimensional simplicial complex K.

Given a simplicial complex K, its dimension d and a discrete Morse function f on K, the algorithm 2 in SI Appendix determines the weights of critical simplices in K. Given an unweighted and undirected graph G, we restrict the construction of clique complex K by including simplices up to a maximum dimension d. Then, the algorithm 3 in SI Appendix creates the filtration of clique complex K based on weights of critical simplices as described in the Theory section. In SI Table S1, we describe the role of key variables which appear in algorithms 1, 2 and 3. In SI Appendix, we also give a time complexity analysis for these algorithms.

Rationale for the choice of function on vertices

In order to construct a discrete Morse function f on clique complex K corresponding to a graph G using our algorithm 1, a real-valued function g has to be fixed on the 0-simplices of K (See lines 7–9 in algorithm 1). Let deg denote the maximum degree of a vertex in the graph . Our choice for the function value on the vertices or 0-simplices, , is as follows:where degree(v) is the degree of the vertex and corresponding to each vertex is a random number (noise) generated using the uniform distribution on the interval (0, 0.5).

In the Theory section, we had highlighted the Theorem 2.11 by Forman[9] which gives a lower bound on the number of critical p-simplices, m, in a simplicial complex K as the p-Betti number β. The choice of the real-valued function g in algorithm 1 plays a key role in determining if m is close to the theoretical minimum β stated above. In the Theory section, we have shown that the number of critical simplices determines the effective number of filtration weights to study the persistent homology of a clique complex (See Eq. 10). This motivated our choice for the real-valued function (Eq. 13) which determines the weights of 0-simplices, and the rationale for this choice is as follows.

Ignoring the noise term in Eq. 13, the reader can discern our intuition for choosing the function for any vertex v in G with the following example. Consider the simple example of the clique complex K corresponding to a graph G in Fig. 1. Here, we would like to obtain a discrete Morse function f on K such that the number of critical simplices is close to the theoretical minimum. This requirement applies to simplices of any dimension in K, and in the context of this example, we would like the number of critical 1-simplices (edges) to be as close as possible to the 1-Betti number β1 of K. Note that for the example clique complex K in Fig. 1.

Let us now examine lines 11–23 in algorithm 1. Consider any edge such that . While assigning the function value to the edge e in algorithm 1, the edge e and the vertex v are guaranteed to be not critical provided that the if condition in the line 16 is satisfied. This is a consequence of the definition of a critical simplex (See Eq. 7). Thus, we would like to force this if condition to be True for as many edges as possible. Moreover, once the function value of the 1-simplex e is set, we set the variable Flag[v] to 1 in line 18, and this subsequently forces the if condition in the line 16 to fail for all other edges in the graph that contain v and have function value .

Let us now examine the edge in Fig. 1 which is anchored by vertices v7 and v8 with degree 5 and 1, respectively. As the degree of a vertex gives the number of edges that contain the vertex, v7 is part of 4 other edges apart from e78 while v8 is part of only the edge e78. Suppose e78 is the first edge chosen for the function assignment in line 11 of algorithm 1 and both Flag[v7] and Flag[v8] for the anchoring vertices are 0. We would then prefer that the if condition in the line 16 is satisfied for e78, and as v7 is part of 4 other edges apart from e78 while v8 is part of only e78, ideally Flag[v8] is set to 1 instead of Flag[v7], in other words, we need the function value . We emphasize that this choice of v8 over v7 prevents the forced failure (described in the previous paragraph) of the if condition for the 4 other edges apart from e78 that contain v7.

The above example suggests a need for a function g on the vertices that has an inverse relationship with the degree of the vertices. Hence, our choice provides a simple and effective solution for the above requirement. As a consequence of this choice for the function g, the weights assigned to the simplices by algorithm 1 reflect the degree of its constituent vertices. The degree of a vertex can be thought of as a measure of its importance in the network. Hence, the intuition behind the assignment of weights to simplices by our method is to have an inverse or opposite relationship between ‘weight’ and ‘importance’ of a simplex while simultaneously satisfying the definition of a discrete Morse function. This inverse relationship instead of a proportional relationship between ‘importance’ and ‘weight’ plays a key role in our filtration scheme which is based on the sequence of level subcomplexes corresponding to the increasing sequence of critical weights (See Theory section and algorithm 3 in SI Appendix). Thus, within the constraint of ensuring that the definition for a simplicial complex is satisfied at each stage of the filtration, our scheme prioritizes the addition of simplices with higher ‘importance’ at earlier stages of the filtration due to their lower weights.

We now provide a rationale for the addition of a random noise in Eq. 13. As reasoned above, we would like to force the if condition in line 16 of algorithm 1 to be True for as many edges as possible. Consider an edge such that . The absence of a random noise in Eq. 13 forces . Thus, irrespective of the state of Flag[v] and Flag[w], the if condition fails. This implies that the set (See Eq. 5) corresponding to the edge e is forced to be empty since the only other possibility is which cannot be the case because algorithm 1 produces a discrete Morse function (See Theorem and Eq. 6). Thus, provided e is not a face of any higher dimensional simplex, e would be a critical simplex irrespective of the state of Flag[v] and Flag[w]. Hence, we would like while also retaining the inverse relationship of the function with the degree. Thus generating a small random noise in the range for each vertex as in Eq. 13 provides a simple resolution. We remark that the above argument can be generalized to higher-dimensional simplices, and thus, provides the intuition for the addition of noise in line 21 of algorithm 1.

We remind the readers that our initial motivation was not to develop a scheme to construct the optimal discrete Morse function on a clique complex corresponding to a graph. Rather, our main goal is to develop a systematic filtration scheme to study persistent homology in unweighted and undirected networks. In fact, constructing an optimal discrete Morse function in the general case has been shown to be MAX-SNP Hard[49]. The primary utility of our scheme is to create a filtration by assigning weights to simplices in the clique complex K of a graph G. However, we next report our empirical results from an exploration of model and real-world networks which underscore the following. Although our scheme is not optimal in the sense of minimizing the number of critical simplices, in practice, it achieves near-optimal results in several model and real-world networks (Table 1). Hence, our scheme based on discrete Morse theory reduces the number of filtration steps and increases the applicability of persistent homology to study complex networks.

Table 1

The table lists the number of p-simplices (n), the number of critical p-simplices (m) and the p-Betti number β for clique complexes corresponding to model and real networks.

Network	n 0	m 0	β 0	n 1	m 1	β 1	n 2	m 2	β 2	n 3	m 3	β 3
ER model with n = 1000 and p = 0.004	1000	90	21	2007	1090	1021	7	0	0	0	0	0
WS model with n = 1000, k = 4 and p = 0.5	1000	123	1	2000	991	864	137	5	0	0	0	0
BA model with n = 1000 and m = 2	1000	8	1	1996	949	942	55	0	0	0	0	0
Spherical model with n = 1000, T = 0, k = 4 and γ = ∞	1000	172	126	2028	118	0	2029	180	0	1321	554	446
Hyperbolic model with n = 1000, T = 0, k = 4 and γ = 2	1000	144	144	2593	20	0	5440	426	0	11456	8159	7753
US Power Grid	4941	573	1	6594	1671	1080	651	21	0	90	15	13
Email communication	1133	6	1	5451	1694	1186	5343	871	53	3419	1577	1262
Route views	6474	17	1	12572	2459	2157	6584	627	19	5636	3335	3013
Yeast protein interaction	1870	272	173	2203	424	318	222	12	0	41	12	7
Hamsterster friendship	1858	33	23	12534	4484	2970	16750	5324	1880	10015	4814	2874
Euro road	1174	213	26	1417	425	237	32	1	0	0	0	0
Human protein interaction	3133	269	210	6149	2454	2298	1047	109	1	142	35	24

Note that the dimension p of simplices ranges from 0 to 3. In case of model networks, the statistics is reported for a particular instance of ER graph with and , WS graph with , and , BA graph with and , Spherical random graph produced from HGG model with , , and , and Hyperbolic random graph produced from HGG model with , , and . For the statistics corresponding to a more comprehensive list of model networks with different parameters, we refer the readers to SI Table S2. Note that we omit self-loops in the real networks considered here.

Application to model and real networks

Given a model or real network G, we limit our study of persistent homology to the 3-dimensional clique simplicial complex K corresponding to G. In other words, during the construction of the clique complex, we only include p-simplices which have dimension (See Theory section). On this 3-dimensional clique complex K, we create the corresponding filtration based on the assigned weights to simplices using discrete Morse theory. In SI Appendix, we present algorithm 3 which outlines the procedure to compute the filtration weights of simplices. Thereafter, we make use of GUDHI[50], a C++ based library for Topological Data Analysis (http://gudhi.gforge.inria.fr/) to study the persistent homology of this filtration of K.

Each hole of any dimension in K has a corresponding birth filtration weight and a death filtration weight (See Theory section). We normalize the birth and death filtration weights of all holes in K by dividing with . In other words, w is 1 plus the maximum value among weights assigned to the simplices in K. We also make the convention that a normalized death filtration weight of 1 for a hole in K represents that the particular hole never dies.

A H-barcode diagram corresponding to a filtration of the clique complex K is a graphical representation containing horizontal line segments, each of which represents a p-hole in K, plotted against the x-axis ranging from 0 to 1 which corresponds to the normalized filtration weights of simplices in K[51]. A horizontal line in the H-barcode diagram of K is referred to as a barcode. Thus, a barcode in the H-barcode diagram of K which begins at a x-axis value of w1 and ends at a x-axis value of w2 represents a p-hole in K whose birth and death weights are w1 and w2, respectively. In Figs 3, 4 and 5, we display the barcode diagrams for model and real-world networks analyzed here.

Figure 3

Barcode diagrams for H0 and H1 in model networks. (a) ER model with and . (b) WS model with , and . (c) BA model with and . (d) Spherical random graphs produced from HGG model with , , and . (e) Hyperbolic random graphs produced from HGG model with , , and .

Figure 4

Barcode diagrams for H0 and H1 in real networks. (a) US Power Grid. (b) Email communication. (c) Route views. (d) Yeast protein interaction. (e) Hamsterster friendship.

Figure 5

Barcode diagrams for H2 in model and real networks. (a) Spherical random graphs produced from HGG model with , , and . (b) Hyperbolic random graphs produced from HGG model with , , and . (c) US Power Grid. (d) Email communication. (e) Route views. (f) Yeast protein interaction. (g) Hamsterster friendship.

Model networks

In this work, we have investigated the persistent homology of unweighted and undirected graphs corresponding to five model networks, namely, ER, WS, BA, hyperbolic random graphs and spherical random graphs. We have considered model networks with 1000 vertices and expected average degree 4, 6 and 8. In main text, we report results for model networks with expected average degree 4, and in SI, those with expected average degree 6 and 8.

The H0 barcode diagram of BA networks indicate a low number of 0-holes in BA networks across the entire filtration (Fig. 3 and SI Figs S1 and S2). A standard result in algebraic topology[36] gives that the 0-Betti number of a simplicial complex K is equal to the number of connected components in K. In other words, the above observation indicates that the scale-free BA network has a strong tendency to maintain a low number of connected components during filtration (Fig. 3 and SI Figs S1 and S2). In contrast, both the random ER network and small-world WS network have a relatively high number of connected components at initial phase of the filtration and then progress towards a more connected network at later stages of the filtration (Fig. 3 and SI Figs S1 and S2). This indicates that the simplices in the clique complex which are key to the connectivity of the model networks are introduced very early into the filtration for the scale-free BA network while this is not the case for the random ER network or small-world WS network. The H1 barcode diagram of BA networks also indicate late introduction of 1-holes during filtration in contrast to both ER and WS networks where 1-holes appear across a wider range of the filtration (Fig. 3 and SI Figs S1 and S2). Moreover, in ER and WS networks, it is interesting to observe that the 1-holes start to appear at roughly the same stage of filtration which corresponds to a sharp reduction in the number of connected components (Fig. 3 and SI Figs S1 and S2).

In contrast to ER, WS and BA networks, the spherical and hyperbolic networks are characterized by a relatively high 0-Betti number β0 and a low 1-Betti number β1 (Table 1 and SI Table S2). Simply stated, this observation on the magnitude of β0 indicates that both hyperbolic and spherical networks have a higher number of connected components in comparison to the ER, WS and BA networks of similar size, i.e., number of vertices, and average vertex degree (Table 1 and SI Table S2). Although, both spherical and hyperbolic networks exhibit a higher number of connected components, they differ from each other with respect to the evolution of these connected components during filtration (Fig. 3 and SI Figs S1 and S2). The hyperbolic model maintains a relatively low number of connected components until very late in the filtration wherein there is a sharp increase in the number of connected components (Fig. 3 and SI Figs S1 and S2). This is in contrast with the behavior of the H0 barcode diagram of the spherical model which exhibits a more distributed evolution of connected components during filtration (Fig. 3 and SI Figs S1 and S2). In addition, the low β1 for spherical and hyperbolic networks conveys the lack of 1-holes in these networks. A possible reason for this observation is the incidence of higher number of 2-simplices in the clique complex K of the spherical and hyperbolic networks in comparison to the ER, WS and BA networks (Table 1 and SI Table S2). Note that the formation of a 2-simplex can potentially fill in a 1-hole, and thus, result in a low value for β1 (See Theory section). Such a behaviour is also seen in the H2 barcode diagrams of spherical and hyperbolic networks wherein the 2-holes have very short persistence since the addition of 3-simplices successively fill in the 2-holes (Fig. 5 and SI Fig. S3). The H3 barcode diagrams of spherical and hyperbolic networks (See SI Figs S3 and S4) also indicate a clear difference in the evolution of their corresponding topological features during the filtration. The H2 and H3 barcode diagrams of ER, WS and BA networks do not provide any insight into network structure primarily due to a lack of higher-order correlations in these model networks that is essential for the formation of 2-holes and 3-holes.

A visual inspection of the barcode diagrams for the five model networks (Figs 3 and 5 and SI Figs S1–S4) suggests that the different models can be distinguished based on their persistent homology. In Theory section, we had introduced the bottleneck distance which can be employed to quantify the differences between persistence diagrams from filtration of clique complexes corresponding to different model networks. Recall that the persistence diagram of a d-dimensional simplicial complex K is a compact representation of the persistent homology of K which encompasses topological information across all d dimensions. Figure 6 and SI Table S3 give the bottleneck distance between different model networks with the same number of vertices and similar average vertex degree. For each of the five model networks, 10 random samples are generated by fixing the number of vertices n and other parameters of the model. In Fig. 6 and SI Table S3, we report the distance between two different models as the average of the distance between each of the possible pairs of the 10 sample networks corresponding to the two models along with the standard error. We find a relatively higher distance between a random instance of a BA network and a random instance of a ER network with the same number of vertices and similar average degree (Fig. 6 and SI Table S3). Similarly, we observe a relatively higher distance between a random instance of a BA network and a random instance of a WS network with similar size and average degree. In contrast, a relatively lower average distance is observed between a random instance of a ER network and a random instance of a WS network with similar size and average degree (Fig. 6 and SI Table S3). These observations indicate a similarity between networks generated by ER and WS models in terms of their corresponding persistence diagrams and also show that the BA model exhibits topological properties that are different from the other two model networks with the same number of vertices and similar average degree. Moreover, the average distance between a random instance of a spherical network and a random instance of a hyperbolic network with similar size and average degree is very high (Fig. 6 and SI Table S3). The last observation is a reflection of the differences in the persistent homology of the clique complexes corresponding to spherical and hyperbolic networks. Finally, the nature of differences observed between persistence diagrams of different model networks using bottleneck distance as shown in Fig. 6 remain consistent if the 1-Wasserstein or 2-Wasserstein distance is employed in place of the bottleneck distance (data not shown).

Figure 6

Bottleneck distance between persistence diagrams of model networks, namely, ER model with and , WS model with , and , BA model with and , Spherical random graphs produced from HGG model with , , and , and Hyperbolic random graphs produced from HGG model with , , and . For each of the five model networks, 10 random samples are generated by fixing the number of vertices n and other parameters of the model. We report the distance (rounded to two decimal places) between two different models as the average of the distance between each of the possible pairs of the 10 sample networks corresponding to the two models along with the standard error.

Table 1 and SI Table S2 list the empirical data on the number of p-simplices n, the number of critical p-simplices m, that our algorithm achieves and the p-Betti number β of the clique complexes corresponding to the five model networks. In Table 2 we report a value μ which indicates the optimality of our algorithm with respect to reducing the number of critical simplicies for each model network that has been analyzed. The definition for μ is as follows.

Table 2

The table lists the value of optimality indicator μ for various model and real networks analyzed here.

Network	Optimalitiy indicator μ
ER model with n = 1000, p = 0.004	0.924 ± 0.004
ER model with n = 1000, p = 0.006	0.947 ± 0.004
ER model with n = 1000, p = 0.008	0.959 ± 0.006
WS model with n = 1000, k = 4, p = 0.5	0.890 ± 0.003
WS model with n = 1000, k = 6, p = 0.5	0.917 ± 0.007
WS model with n = 1000, k = 8, p = 0.5	0.906 ± 0.003
BA model with n = 1000, m = 2	0.989 ± 0.003
BA model with n = 1000, m = 3	0.985 ± 0.004
BA model with n = 1000, m = 4	0.964 ± 0.006
Spherical model with n = 1000, γ = ∞, T = 0, k = 4	0.925 ± 0.007
Spherical model with n = 1000, γ = ∞, T = 0, k = 6	0.901 ± 0.003
Spherical model with n = 1000, γ = ∞, T = 0, k = 8	0.887 ± 0.003
Hyperbolic model with n = 1000, γ = 2, T = 0, k = 4	0.939 ± 0.013
Hyperbolic model with n = 1000, γ = 2, T = 0, k = 6	0.927 ± 0.007
Hyperbolic model with n = 1000, γ = 2, T = 0, k = 8	0.921 ± 0.005
US Power Grid	0.893937
Email communication	0.871847
Route views	0.952140
Yeast protein interaction	0.942157
Hamsterster friendship	0.793236
Euro road	0.840678
Human protein interaction	0.957924

Eq. 14 gives the definition of μ which is an indicator for the proximity of our algorithm to the optimal case. The value of μ ranges from 0 to 1 with indicating that our algorithm achieves exactly the theoretical minimum number of critical simplices while indicating that all simplices are critical. For model networks, the value reported in this table is the average of μ across 10 samples of each model network for a chosen set of parameter values along with the corresponding standard deviations.

Here d is the dimension of the corresponding clique complex. The value μ corresponding to a particular discrete Morse function on the clique complex of a given network is an indicator of its optimality with respect to minimizing the number of critical simplices. The value of μ ranges from 0 to 1. indicates that the discrete Morse function achieves exactly the minimum number of critical simplices thereby corresponding to the most optimal situation, while indicates that all simplices are critical thereby corresponding to the least optimal case. Note that the value of μ for a particular network increases linearly with a decrease in the number of critical simplices. In Table 2, for model networks, the value reported is the average of μ across 10 samples of each model network for a chosen set of parameter values along with the corresponding standard deviations. Based on the data presented in Table 2, we report that our algorithm achieves near-optimal results in terms of reducing the number of critical simplices for each of the five model networks analyzed. Moreover for these model networks, data presented in SI Table S2 underlines the close proximity of m to the theoretical lower bound β across all dimensions (See Theory section). We remark that the worst case scenario corresponds to all simplices in the clique complex being critical which is never attained by our algorithm in the analyzed networks.

Comparison of our method with dimension based filtration in model networks

In SI Figs S5 and S6, we show the H0, H1, H2 and H3 barcode diagrams for model networks with expected average degree 4 obtained using the dimensional filtration scheme used in Horak et al.[19]. We restrict our investigation to the three-dimensional clique complex while computing the barcode diagrams for model networks using the dimensional filtration scheme of Horak et al.[19]. In SI Figs S5 and S6, we normalize the filtration index to be in the range 0 to 1, and p-holes with normalized filtration index 1 indicate that they never die. Moreover, we also report in SI Table S4 and SI Fig. S7, the bottleneck distances between persistence diagrams of model networks obtained by the dimensional filtration scheme used in Horak et al.[19], and we find that the resultant barcode diagrams and the bottleneck distances between the persistence diagrams are inconclusive to distinguish between the five model networks which have similar size and average vertex degree.

The discrete Morse function produced by our algorithm is empirically near-optimal for the model networks analyzed here, in the sense of minimizing the number of critical simplices (See Table 2). On the other hand, it should be noted that the dimension function used in Horak et al.[19] is the most non-optimal in this regard. Moreover, in terms of persistent homology, our methods based on the discrete Morse function constructed using algorithms 1, 2 and 3 have the distinct feature of being able to distinguish between various model networks with an efficient filtration.

Although the dimension function used in Horak et al.[19] is algorithmically efficient in terms of having a low number of filtration steps, unlike our method it does not conclusively distinguish between different model networks based on the barcode diagrams and also the bottleneck distances between the corresponding persistence diagrams. Thus our discrete Morse function is a good candidate for applications in both, computational aspects of discrete Morse theory as well as in persistent homology of unweighted networks, and achieves a tradeoff between efficiency and applicability.

Real networks

In this work, we have investigated the persistent homology of unweighted and undirected graphs corresponding to seven real-world networks and the barcode diagrams for five of these real networks is shown in Figs 4, 5 and SI Fig. S4. Based on the H0 barcode diagrams, the behavior of most real networks considered here can be broadly classified into two categories. Real networks such as the Email communication, the Hamsterster friendship and the Route views exhibit a relatively low number of connected components across the entire range of filtration (Fig. 4). On the other hand, the two biological networks, namely, the Yeast protein interaction and the Human protein interaction, exhibit a sharp increase in the number of connected components at the later stages of filtration. The H0 barcode diagrams for the US Power Grid and Euro road do not conform with either of the above characterizations. The H0 barcode diagram of the US Power Grid network reveals that though there exists only a single connected component at the end of the filtration, there are a considerable number of non-persisting connected components that appear and subsequently disappear during filtration (Fig. 4). The H0 barcode diagram of the Euro road network shows a more distributed increase in the number of connected components (data not shown).

In the context of the H1 barcode diagrams, the real networks considered here exhibit similar properties with 1-holes appearing late in the filtration (Fig. 4). The H2 barcode diagrams reveal a lack of 2-holes with long persistence in both biological networks, as well as the Route views network, the Euro road network and the US Power Grid network (Fig. 5). In contrast, from the H2 and H3 barcode diagrams (Fig. 5, SI Fig. S4) we find that the social network, Hamsterster friendship, and the Email communication network exhibit a relatively high number of 2-holes and 3-holes with longer persistence.

Table 2 lists the value of the optimality indicator μ (See Eq. 14) for each of the seven real networks analyzed here. This data indicates near-optimal performance of our algorithm with respect to minimizing the number of critical simplices for each of these seven real networks. Table 1 also lists the empirical data on the number of critical p-simplices, m, that our algorithm achieves and the p-Betti number β of the clique complexes across each dimension p, corresponding to the seven real networks analyzed here.

Conclusions

To conclude, we have proposed a systematic scheme based on discrete Morse theory to study the persistent homology of unweighted and undirected networks. Our methods leverage the concept of critical simplices to permit a reduced filtration scheme while simultaneously admitting a finer inspection of the changes in topology across the filtration of a clique complex corresponding to an unweighted network. Moreover, our proposed algorithm to construct a discrete Morse function on the clique complex of a simple graph achieves close to optimal number of critical simplices for several model and real networks that have been studied here. Furthermore, based on visual representations of persistent homology such as the barcode diagrams as well as quantitative information in the form of distance between persistence diagrams, our methods successfully distinguish various model networks that exhibit inherently different properties. This motivates the application of our methods to real-world networks. We report the results obtained for seven real-world networks that are well studied in the network science community and observe certain patterns in the evolution of their topological features across the filtration. For instance both biological networks, namely the Yeast protein interaction network and the Human protein interaction network exhibit similar characteristics with respect to the H0, H1, H2 and H3 barcode diagrams. Similarly, both the Email network and the Hamsterster friendship network, exhibit shared features with respect to H0, H1, H2 and H3 barcode diagrams that vary from the characteristics of the two biological networks considered here. Our observations hint at the ability and possible applications of our methods to detect and classify real-world networks that are inherently different.

Future directions and ongoing work include examining the significance of critical simplices in the context of real-world networks. In other words, we aim to determine whether a critical edge in the context of discrete Morse theory holds any key significance when it is viewed as a link between two real entities in a real-world network. We also intend to explore the presence or absence of a correlation between the notion of critical simplices and network curvature. Since discrete Morse theory captures information about the Euler characteristic of the clique complex corresponding to a graph, the presence of such a correlation could potentially signify a close relationship between the discrete curvature of a graph and its topology, much like in the case of smooth, compact surfaces wherein the Gauss-Bonnet theorem relates the Gaussian curvature of a surface to its Euler characteristic.

Appendix

Supplementary Figures

Supplementary Tables