Construction of and efficient sampling from the simplicial configuration model

TL;DR

This paper introduces the simplicial configuration model as a null model for complex systems represented by simplicial complexes, along with an efficient MCMC sampler, enabling comparison with empirical data and revealing organizational structures.

Contribution

It presents the first efficient, uniform MCMC sampler for the simplicial configuration model, facilitating principled null model comparisons for multi-node interaction data.

Findings

The sampler is efficient and uniform.

Application to real systems shows non-random organization.

Rejection of null hypothesis in two systems.

Abstract

Simplicial complexes are now a popular alternative to networks when it comes to describing the structure of complex systems, primarily because they encode multi-node interactions explicitly. With this new description comes the need for principled null models that allow for easy comparison with empirical data. We propose a natural candidate, the simplicial configuration model. The core of our contribution is an efficient and uniform Markov chain Monte Carlo sampler for this model. We demonstrate its usefulness in a short case study by investigating the topology of three real systems and their randomized counterparts (using their Betti numbers). For two out of three systems, the model allows us to reject the hypothesis that there is no organization beyond the local scale.

Figures

Equations

$$\mathrm{Pr}(K;\mathbf{d},\mathbf{s})=1/|\Omega (\mathbf{d},\mathbf{s})|$$

$$\varphi =e^{-\frac{1}{2}\left(⟨d^2⟩/⟨d⟩-1\right)\left(⟨s^2⟩/⟨s⟩-1\right)},$$

$$L_{\max }^{\ast }\le 2\max \mathbf{s}.$$

Full text

Network science’s approach to complexity rests onto the tacit hypothesis that the structure of complex systems is reducible to the pairwise interaction of their constituents. It is often a valid premise and, as a result, network science has been extremely successful in, e.g., both predicting Pastor-Satorras et al. (2015) and controlling Liu and Barabási (2016) the behavior of complex systems, inferring their function from their structure Porter et al. (2009); Newman (2012), and so on. Networks, however, might not be as ubiquitous as previously thought. It has been shown recently that the structure of a number of complex systems, such as the brain Dabaghian et al. (2012); Giusti et al. (2016), protein interactions Xia and Wei (2014) and social systems Hébert-Dufresne et al. (2015); Stolz et al. (shed), cannot always be reduced to the sum of pairwise interactions. For these systems, it is now known that network representations can give an incomplete picture: When many-body interactions are broken down into multiple pairwise interactions (cliques), high-order information simply disappears Zuev et al. (2015).

Simplicial complexes generalize graphs by encoding many-body interactions explicitly; they have hence been proposed as a complementary description of the structure of complex systems Ellis and Klein (2014); Curto and Itskov (2008); Horak et al. (2009); Patania et al. (2017). Different from hypergraphs, they are equipped with an implicit notion of containment. If nodes $(v_{1},\ldots,v_{q+1})$ are involved in a $q$-dimensional interaction, then it is implicit that all possible lower dimension interactions involving the same nodes also exist [for example $(v_{1},\ldots v_{q})$ and $(v_{1},v_{3})$]. While it might appear constraining, this property actually arises in all systems where interactions are maximal, e.g., in scientific collaborations (largest cohesive group of collaborators) or gene activation pathways (largest group of collectively activated genes). Furthermore, it is found in many processed relational datasets, e.g., in clique complexes, obtained by mapping the cliques of networks to simplices Petri et al. (2013); Sizemore et al. (2017), or in filtered simplicial complexes Horak et al. (2009). Simplicial complexes thus offer a natural and compact description of the structure of complex systems, both when-high order structures are explicitly available, or when they are extracted from low-order information.

This application of simplicial complexes has led to promising discoveries: We now better understand, for instance, how to detect large viral recombination events Chan et al. (2013), how brain networks reorganize under drugs Petri et al. (2014), and how the atomic structure of amorphous solids is hierarchically organized Hiraoka et al. (2016). It has become crucial to establish the statistical significance of these findings, a task for which random null models will be needed. There is already a rich and growing literature on random simplicial complexes and topology, ranging from simplicial generalization of Erdös-Rényi models, amendable to analytical treatment Kahle (2014); Costa and Farber (2016), to equilibrium formulations of simplicial complex ensembles Courtney and Bianconi (2016); Zuev et al. (2015), and growth models that reproduce various emergent patterns observed in real systems Bianconi and Rahmede (2016); Wu et al. (2015). However, null models—in the sense of network science—are still wanting Fosdick et al. (shed); Orsini et al. (2015).

We address this issue by refining a recently proposed generalization Courtney and Bianconi (2016) of the (simple) configuration model of network science Molloy and Reed (1995); Newman et al. (2001); Fosdick et al. (shed), which we dub the simplicial configuration model (SCM). Different from Ref. Courtney and Bianconi (2016), we think of our model as a null hypothesis for real systems; we therefore develop a numerical and statistical toolbox instead of focusing on closed ensemble averages. This entails a number of interesting results: One, we define the first simplicial configuration model able to describe arbitrary complexes, in line with our goal of obtaining a generic null model (Sec. I). Two, we propose and analyze an efficient and rigorous sampling algorithm for this model (Sec. II). Three, we use the model to investigate real datasets and show—now using sound statistical arguments—that the local structure of these systems does not always explain their mesoscale structure (Sec. III). We conclude by listing a few important open problems.

I Simplicial configuration model

Informally, a labeled simplicial complex $K$ is the high-order generalization of a network. Formally, it is a collection of simplices incident on a node set $V=\{v_{1},\ldots,v_{n}\}$ Hatcher (2000). A $q$–dimensional simplex—the generalization of an edge—is a tuple of $q+1$ distinct nodes $(v_{1},\ldots,v_{q+1})$; we say that this simplex is incident on $v_{1},\ldots,v_{q+1}$. All simplices not included in a larger simplex are called the facets of the complex, whereas a contained simplex is called a face; e.g., if $K$ comprises of $\sigma=(v_{1},v_{2})$ and $\tau=(v_{1},v_{2},v_{3})$, then $\sigma$ is a face of facet $\tau$. It is always assumed that if facet $\sigma=(v_{i},\ldots,v_{j})$ is in the simplicial complex $K$, all elements in the power set of $\sigma$ are also in $K$. Therefore, faces need not be enumerated: The structure of a simplicial complex is fully specified by the list of its facets.

Departing from other recent contributions Courtney and Bianconi (2016), we define the degree $d_{i}$ of a node $v_{i}$ as the number of facets incident on $v_{i}$ and the size $s_{i}$ of a facet $\sigma_{i}$ as the number of nodes it contains (its dimension plus one). This local information is summarized by the sequences $\bm{d}=(d_{1},\ldots,d_{n})$ and $\bm{s}=(s_{1},\ldots,s_{f})$, where $n$ is the number of nodes and $f$ is the number of facets.

With these notions in hand, we define the simplicial configuration model (SCM) as the uniform distribution over all labeled simplicial complexes with degree sequence $\bm{d}$ and facet size sequence $\bm{s}$. In other words, if $\Omega(\bm{d},\bm{s})$ is the set of all labeled simplicial complexes with joint sequences $(\bm{d},\bm{s})$, then the SCM places a probability

[TABLE]

on $K$ if it has sequences $(\bm{d},\bm{s})$, and 0 otherwise. The model of Ref. Courtney and Bianconi (2016) is recovered by setting the size of all facets to a constant $s$.

This particular choice of definition for the SCM is natural for three reasons. First, the SCM directly generalizes the simple CM of network science Fosdick et al. (shed); when $s_{i}=\{1,2\}$ for all facets, one recovers a graph ensemble with degree sequence $\bm{d}$. Second, the SCM does not include any correlation—the structure is maximally random beyond the local level. This is reminiscent of the equivalent network model. Third, the SCM can describe the local structure of any simplicial complex, since it allows for arbitrary degree and size sequences. This property is not common to all random models of simplicial complexes, for good reasons; many models are constructed with a focus on the calculation of closed-form expression for a few properties (e.g., the asymptotic entropy) Zuev et al. (2015); Courtney and Bianconi (2016); Kahle (2014). This commends simplifying assumptions, e.g., a regular facet size sequences Courtney and Bianconi (2016). Our definition of the SCM forgoes these simplifications to accommodate arbitrary local structures, at the expense of analytical tractability.

II Efficient sampling algorithm

II.1 Constraints on the support

For the SCM to be of any use, one needs to be able to sample from it. This is far from a trivial problem, because there are numerous constraints on the support of the model. It will be easier to see these constraints by first switching to the equivalent graphical representation of simplicial complexes.

In this representation, facets are replaced by nodes (we denote by $F$ this new node set, and by $V\cup F$ the complete node set), and an edge connects facet $\sigma_{i}\in F$ to node $v_{j}\in V$ in $B$ if and only if $\sigma_{i}$ is incident to $v_{j}$ in $K$ (see Fig. 1). Because $B$ encodes the structure of $K$ without ambiguity, we can think of the model in terms of either representations.

As such, one could be tempted to assume that sampling from the SCM of parameters ($\bm{d},\bm{s}$) is equivalent to uniformly sampling from all bipartite graphs with these degree sequences—a solved problem Miklós et al. (2013). But this would be wrong: The mapping is not bijective. This is, in fact, where the constraints on the support of the SCM become apparent Courtney and Bianconi (2016). Let us introduce the notion of sequence preserving bipartite graph to formalize these constraints. We say that a bipartite graph $B$ with joint degree sequences $(\bm{d},\bm{s})$ is sequence preserving if, upon interpretation of $B$ as a simplicial complex, one obtains a simplicial complex with facet size sequence $\bm{s}$ and generalized degree sequence $\bm{d}$.

Not all bipartite graphs are sequence preserving, and there are two reasons for this, both related to the fact that we think of the nodes in $F$ as facets. The first reason is the inclusion of at least one facet: If there is a $\sigma_{i}\in F$ such that the neighborhood of $\sigma_{i}\subseteq\sigma_{j}$ for some $j\neq i$, then $B$ is not sequence preserving [see Fig. 2 (a)]. When this occurs, $\sigma_{i}$ is included in $\sigma_{j}$; the corresponding simplicial complex is thus either ill specified (facets cannot contain other facets, by definition) or does not have the same degree and size sequences as $B$ (if we simply remove $\sigma_{i}$). For similar reasons, if two or more edges connect the same pair of nodes in $B$, then the graph is not sequence preserving [see Fig. 2 (b)].

The sampling space would not be too constrained if these non-sequence-preserving bipartite graphs were rare. Sampling would then be easy. Unfortunately, it is straightforward to show that-non sequence-preserving graphs are far more common than sequence–preserving ones, by adapting the calculations of Ref. Bender and Canfield (1978). We find that the fraction $\phi$ of bipartite graphs with degrees $(\bm{d},\bm{s})$ not featuring parallel edges rapidly tends to

[TABLE]

where $\langle x^{k}\rangle$ is the $k$th moment of the sequence $\bm{x}$, and where it is assumed that the elements of $\bm{d}$ and $\bm{s}$ do not grow with $n$ (i.e., $B$ is sparse). Thus, based on the presence of multi-edges alone, there is a stringent upper bound on the fraction of bipartite graphs that are actually in the support of the SCM. An even smaller fraction remains after the bipartite graph with included neighborhood are removed.

II.2 Markov chain Monte Carlo method

To sample from the SCM, then, one needs to sample uniformly from a very constrained space, i.e., that of all sequence–preserving bipartite graphs with joint degree sequence $(\bm{d},\bm{s})$. Previously proposed approaches such as rejection sampling do not work well Courtney and Bianconi (2016), because natural proposal distributions (e.g., stub matching) give an appreciable weight to non-sequence-preserving bipartite graphs [see Eq. (2)]. Thus, we turn to the Markov chain Monte Carlo (MCMC) sampling strategy 111We provide a reference c++ implementation of the sampler as well as tutorials at https://www.github.com/jg-you/scm, which has been used with great success for the CM Fosdick et al. (shed); Miklós et al. (2013). The general idea is to construct a random chain of sequence–preserving bipartite graphs $B_{0},\ldots,B_{T}$, to sample from this chain at regular intervals, and to treat the samples as if they had been drawn identically and independently from the ensemble. The algorithm will be correct if the chain is ergodic (time averages equal ensemble averages) and uniform (all non isomorphic $B$ are represented equally). These properties are determined by the allowed transformations $B_{t}\to B_{t+1}$ and the resulting transition matrix $\bm{\pi}$, where $\pi_{ij}$ is the probability that $B_{j}$ follows $B_{i}$ in the chain. If the move set connects the space and the chain is aperiodic, then the chain will be ergodic. If the transition matrix is doubly stochastic (all rows and columns sum to 1), then the chain will be uniform.

We claim that the following set of moves satisfies all three conditions. Consider $L$, a random variable on the support $\mathcal{L}=\{2,3,\ldots,L_{\max}\}$, where $L_{\max}$ is a parameter and the distribution $\mathbb{P}[L=\ell]$ is arbitrary but non zero everywhere on $\mathcal{L}$ (for illustration purposes, we will use $\mathbb{P}[L=\ell;\lambda]\propto e^{-\lambda\ell}$). At each step of the chain, we pick $L$ edges in $B$ (uniformly at random). We cut these edges and randomly match the stubs stemming from facets to the stubs stemming from nodes. If this matching generates a sequence–preserving bipartite graph $B^{\prime}$, then we accept the move; otherwise we resample $B$. This set of moves is similar to the double-edge swap commonly used in graph MCMC Fosdick et al. (shed). The only difference is the variable number of rewired edges, added to help the sampler better navigate the constrained support Miklós et al. (2013). Much like its graphical counterpart, the resulting MCMC algorithm is efficient since drawing $L$ edges and checking for resampling can be done in polynomial times.

The chain is aperiodic because the above set of moves yields a doubly stochastic transition matrix for any distribution $\mathbb{P}$: The total number of possible transitions at each configuration is a constant independent from the configuration considered (resampling guarantees this) Fosdick et al. (shed). The chain is also aperiodic, because there exists orbits of period 1 (resampling steps) and 2 (all moves are reversible) for any nontrivial $(\bm{d},\bm{s})$.

This leaves open the question of whether the support of the SCM is connected by the set of moves or not. We argue that it is, for all $L_{\max}\geq L_{\max}^{*}$, where $L_{\max}^{*}$ is bounded by

[TABLE]

To prove this, one would have to show that given two sequence–preserving bipartite graphs $B_{1}$ and $B_{2}$, it is always possible to find a $B_{3}$ such that \big{|}\Delta_{+}\bigl{(}K(B_{1}),K(B_{2})\bigr{)}|\geq|\Delta_{+}\bigl{(}K(B_{1}),K(B_{3})\bigr{)}\big{|}, where $\Delta_{+}$ is the set of facets in $K(B_{2})$ that are not in $K(B_{1})$, and $\Delta_{-}$ is the set of facets in $K(B_{1})$ that are not in $K(B_{2})$ [$K(B)$ is the simplicial complex associated to the graph $B$]. Although a general proof remains elusive, we propose the following non-rigorous argument, valid for sparse simplicial complexes (simplicial complexes with bounded $\max\bm{d}$ and $\max\bm{s}$ in the limit $n\to\infty$).

To construct $B_{3}$, we first select a facet $\sigma$ in $\Delta_{+}$ (incident on the set of nodes $\Sigma$ in $B_{2}$). The conservation of sizes and degrees guarantees that there exists a facet $\tau\in\Delta_{-}$ of the same size. The idea is then to start from $B_{1}$, cut all edges attached to $\tau$ and one edge from every node in $\Sigma$, match the stubs of $\sigma$ to those of $v\in\Sigma$, and finally match the remaining orphaned stubs. This algorithm ensures that $B_{3}$ is closer to $B_{2}$ than $B_{1}$ was, because it removes facets from $\Delta_{\pm}$ (and does not add new facets either: Each $v\in\Sigma$ has at least on facet in $\Delta_{-}$ by the conservation of degrees). In general, it is not guaranteed that the last step can be carried out without creating included faces. However, in sparse simplicial complexes, $\sigma$ is well separated from $\tau$ for almost all $(\sigma,\tau)$, since $B_{1}$ is locally treelike Newman et al. (2001). In such cases, no included faces are created at the last step, and the above algorithm can be carried through for some $(\sigma,\tau)$, generating $B_{3}$. Because this scheme involves at most $L_{\max}^{*}=2\max\bm{s}$ rewired edges (when $|\tau|=|\sigma|=\max\bm{s}$), we obtain the bound of Eq. (3) for infinite sparse SCM. In practice, $L_{\max}=2$ seems to always connect the space (we found no counterexamples), and sampling is more efficient when $L_{\max}\gg 2$ (see Fig. 3)—the value of $L_{\max}^{*}$ is thus more of theoretical than practical interest.

III Null model

We put our efficient MCMC algorithm to the test, by verifying the statistical significance of the structural patterns found in three relational datasets that can be represented as simplicial complexes (see caption of Fig. 4 for details).

Since an instance of the SCM is provided in each case (the real system), we use it as the initial condition for each independent run of the sampling algorithm. Ergodicity implies that the state of the sampler will be uncorrelated with the initial configuration after a sufficiently long burn-in period—the choice of initial condition is ultimately irrelevant. Extrapolating from the results of Fig. 3, we opt for the proposal distribution $\mathbb{P}[L=\ell]=e^{\lambda\ell}/Z$ with $\lambda=1$ and $L_{\max}$ set to $10\%$ of $m=\sum d_{i}=\sum s_{i}$. Non rigorous arguments from expander graph theory suggest $t_{f}=O(m\log m)$ as a good—if overzealous—choice of sampling interval 222We represent the support of the SCM as a graph $\mathcal{G}(L_{\max})=(V,E)$. If $\mathcal{G}$ is an expander, then the sampler yields uncorrelated configuration with high probability after $t_{f}=O(\log|V|)$ steps; the suggested $t_{f}$ follows from the loose upper bound $|V|\leq m!$. A proof that $\mathcal{G}(L_{\max})$ is in fact an expander will depend on the specifics of $(\bm{d},\bm{s},L_{\max})$; however, we note that $\mathcal{G}(L_{\max})$ shares the two fundamental properties with all expanders for sufficiently large $L_{\max}$: It is connected and not bipartite..

Significance results only make sense if they rely on a null model that embodies a natural null hypothesis for the problem at hand Fosdick et al. (shed). For example, the regular CM and its correlated variants usefully show that the network projection of datasets with high order interactions are abnormally clustered Newman (2003). Therefore, we use the sampler to investigate the distribution of a mesoscopic property only accessible when the datasets are encoded as simplicial complexes: The shape of the datasets, as captured by their homology, i.e., the pattern of holes, cavities and higher dimensional voids Hatcher (2000). The homology can be summarized by a series of Betti numbers $\bm{\beta}=(\beta_{0},\beta_{1},\beta_{2},\ldots)$, where $\beta_{k}$ counts the number of structural holes bounded by $k$-dimensional simplices. For example, $\beta_{0}$ counts the number of connected component, $\beta_{1}$ the number of homological cycles in $K$, $\beta_{2}$ the number of holes enclosed by facets of sizes $2$, etc. Since every instance of the SCM has the same fixed local structure but is otherwise maximally random, we expect significant differences between the Betti number $\bm{\beta}$ of an organized simplicial complex and the bulk of the distribution of $\bm{\beta}$ in the corresponding randomized ensembles.

We show in Fig. 4 the distribution of $\beta_{0}$ and $\beta_{1}$ for the SCM associated to the real systems. Looking first at $\beta_{0}$, we find that the structure of the pollinator dataset is essentially random [Fig. 4(a)]. That is, the overwhelming majority of simplicial complexes with the same sequences have similar $\beta_{0}$. In contrast, the $\beta_{0}$ of the disease genome regulation (hereafter diseasome) and crime complexes are highly significant [Fig. 4 (b)–(c)]: A random instance of the SCM has fewer (diseasome) or more (crime) components than the real system with high probability. In one case (crime), the difference is a statistical signature of how the dataset was gathered, namely by looking up the ties of suspects, victims and witnesses already in the dataset, recursively Decker et al. (1991). Because this process creates much larger connected components than random sampling, the resulting $\beta_{0}$ is far from the ensemble average—an effect that we expect to find in any dataset constructed using a similar methodology. In the other case (diseasome), the real system has more components than one would typically expect from the local information alone. The construction procedure does not explain this disparity Goh et al. (2007), meaning that the system must self-organize in a fragmented way, likely for biological or evolutionary reasons.

Turning to $\beta_{1}$ we again find that the structure of the pollinator dataset is typical, and that the same cannot be said of the diseasome and crime datasets. Both simplicial complexes have significantly fewer cycles than expected; i.e., given a cycle, it is more likely to be filled by a simplex in the real system than in the randomized one, suggesting that some form of high order triadic closure is at play Zuev et al. (2015). The difference is, however, much more pronounced in the crime dataset; this could be due to the fact that it describes a social system, whose structure tend to be heavily driven by triadic closure Hidalgo (2016) (and potential high order analogs).

Finally, taking both distributions into account, we conclude that the shape of the pollinator dataset is completely determined by its local structure, while large–scale organizational principles influence the structure of the other datasets. This leads us to two final observations: One, care must be exerted in drawing conclusions about the shape of complex datasets—from the homology point of view there is nothing of note in the structure of the pollinator dataset. Two, some datasets—here the crime and diseasome datasets—are decidedly not random. This raises the question of just how much information must models account for, before they can capture such atypical Betti numbers. Would, for example, adding limited correlations among degrees be sufficient to capture the shape of most real datasets? Or do we need to embrace growth models, with their sophisticated rules, and clustered local structure Wu et al. (2015); Bianconi and Rahmede (2016); Hébert-Dufresne et al. (2015)?

IV Perspectives

As it stands, the SCM already establishes the analysis of simplicial complexes on firmer statistical ground. The next step will be to clarify a number of important open questions, e.g., what is the true value of $L_{\max}^{*}$ for arbitrary simplicial complexes, and what is optimal choice of proposal distribution $\mathbb{P}$ (cf. Fig. 3).

Beyond these obvious questions, the connection between the SCM and the simple CM lead us to a series of natural problems not addressed in this paper. These include the problem of the simpliciality of arbitrary pairs of sequences (i.e., is there a simplicial complex which realize a pair of sequences?) Courtney and Bianconi (2016), related to the problem of constructing initial conditions for the MCMC sampler, when no real system is available. We believe that the solution to such problems will require new insights, as the no–inclusion constraint appears to be a major obstacle to the application of classical methods developed for the analogous graphicality problem Havel (1955); Hakimi (1962).

In closing, we stress that all the above questions and challenges are of technical nature; the model and sampler can already be applied to practical problems Note (1). This could lead to improvements in persistent homology (e.g. statistically sound filtrations of weighted complexes) or a formulation of community detection of simplicial complexes (via modularity Newman and Girvan (2004)), and could provide a new glimpse into the emergence of homology and higher order structural properties in real complex systems.

Acknowledgments

We thank L. Hébert-Dufresne and G. Bianconi for helpful discussions and comments. The authors acknowledge the support of the ADnD project by Compagnia San Paolo (AP, GP), the Fonds de recherche du Québec-Nature et technologies (JGY), the Complex Systems Langrange Lab (FV), and the YRNCS Bridge Grant (AP, JGY). AP is grateful for the hospitality of L.J. Dubé at Université Laval, where parts of the research work was conducted.

AP and JGY contributed equally to this work.