Phase transitions and stability of dynamical processes on hypergraphs

Hypergraphs naturally represent higher-order interactions, which persistently appear in social interactions, neural networks, and other natural systems. Although their importance is well recognized, a theoretical framework to describe general dynamical processes on hypergraphs is not available yet. In this paper, we derive expressions for the stability of dynamical systems defined on an arbitrary hypergraph. The framework allows us to reveal that, near the fixed point, the relevant structure is a weighted graph-projection of the hypergraph and that it is possible to identify the role of each structural order for a given process. We analytically solve two dynamics of general interest, namely, social contagion and diffusion processes, and show that the stability conditions can be decoupled in structural and dynamical components. Our results show that in social contagion process, only pairwise interactions play a role in the stability of the absorbing state, while for the diffusion dynamics, the order of the interactions plays a differential role. Our work provides a general framework for further exploration of dynamical processes on hypergraphs.

Network science has been successful in describing the structure and dynamics of complex systems in many fields of science. Nevertheless, in the vast majority of studies, the mathematical and computational descriptions of these networks are limited to pairwise interactions. In most cases, this order of interaction is a good approximation, describing real phenomena [1][2][3][4]. Recently, it has become increasingly evident that simplifying higher-order interactions as a set of pairwise ones can be misleading in many situations [5,6], as it is the case, for instance, when comparing triangles and clustering [6]. An attempt to formally solve this problem is represented by the introduction of the simplicial complexes approach [7,8]. Although very powerful, as a topological space, simplicial complexes require mutual inclusion for the interactions [9], which is too restrictive for general purposes. To face this problem, one needs to resort to the use of hypergraphs, which by relaxing the assumption of mutual inclusion, allow the representation of a broader range of systems.
Admittedly, the attention to higher-order systems has proliferated lately [5,6,8,[10][11][12][13][14][15][16][17][18][19], with an increasing focus on hypergraphs in fields such as mathematics [6,12,13,15,16,20,21], physics [10,11,[17][18][19], and computer science [22][23][24][25][26][27][28]. Despite this interest, a general theory of dynamical processes on higher-order structures is still largely missing. To fill this important gap, here we build a mathematical framework that allows performing a linear stability analysis for general processes on arbitrary hypergraphs. This approach highlights the importance of the graph-projection −an underlying weighted graph representation of a hypergraph− for the dynamics. The proposed methodology makes it possible to decompose a hypergraph in uniform structures, hence enabling the characterization of their role in the system's dynamics. Finally, to show the usefulness of our approach, we analytically study and recover some results reported for social contagion [8,17], and diffusion processes, also providing new key insights into these paradigmatic dynamics. The framework * These two authors contributed equally discussed in this paper could be used to explore different processes when pairwise interactions are an oversimplification of the system, which has the potential to bring more insights into the understanding of higher order dynamics of interacting systems.
Hypergraph structure. A hypergraph, H = {V, E}, is defined as a set of nodes, V = {v i }, with N = |V| the number of nodes and a set of hyperedges E = {e j }, where e j is a subset of V with arbitrary cardinality |e j |. Note that if max (|e j |) = 2 we recover a graph, whereas one has a simplicial complex if for each hyperedge with |e j | > 2, its subsets are also contained in E. The adjacency matrix [15] can be defined as which can be interpreted as a weighted projected graph. An example of a hypergraph and its graph projection is shown in Fig. 1. The weighted counterpart of the adjacency matrix is defined as where w(|e j |) is a function of the cardinality of the hyperedge, e j , weighting differently the contribution of each hyperedge. For instance, if we consider w(|e j |) = |e j | − 1, we recover the adjacency matrix used in [19]. Furthermore, the Laplacian matrix is defined as L = D − A, where D = diag(k i ) and the degree is given as k i = N j=1 A ij . Note that, in order to define the Laplacian matrix on the weighted case we also need to redefine D accordingly.
We can also decompose the hypergraph into m-uniform hypergraphs, as motivated in [16]. Formally, therefore A = m A m as showed in Fig. 1. Consequently, we can also define D m , L m , and W m .
Stability analysis. A general dynamical process on the hypergraph can be written as where x i is the state of the node v i , f i (x i ) is a R → R function that depends only on the state of node v i and g j (x {ej } ) is a R |ej | → R function that takes all the states of the nodes on the hyperedge e j , here denoted as x {ej } , and compute its contribution to x i . While f i (x i ) represents the internal dynamics of the node, g j (x {ej } ) is the external interaction, expressed as the hyperedge contribution.
We can perform a linear stability analysis for Eq. 4 around a known fixed point, x i = x * i , similarly to [2]. The linearized equations are expressed as where x is a vector whose components are x i . We can express Eq. 5 in its matrix form decomposing a general hypergraph into uniform hypergraphs as where The fixed point x * is stable if where Λ i (M) is the i-th eigenvalue of M. The structures of cardinality m do not contribute to the stability of the fixed point We stress that, although we have a high-order structure (N max(|ej |) in the tensorial representation [13]), near the fixed point, the dynamics is described only by a N 2 matrix.
Let us first consider a diffusion-like process. In this case, g j (x {ej } ) is separable and we can express the interaction function as where the weight (|e j |−1) −1 is necessary due to conservation purposes. Assuming that f i (x i ) is the same for all nodes, the fixed point is symmetrical, i.e., x * i = x * j for any pair i, j. In this case, Eq. 7 reduces to where Therefore, Eq. 5 is expressed as which in matricial form reads where is the vector whose components are i . If µ i are the eigenvalues of L, the stability condition for the fixed point x * can be written as α + βµ i < 0, for all i. Next, considering that the Laplacian matrix is semi-positive definite, the former condition in terms of the largest eigenvalue, µ n , is Now, let us consider a second scenario for g j (x {ej } ) such that which depends, for a node v i , only on x k , where k ∈ e j and k = i, and the weight is the same as for Eq. 10. Assuming that we can write F and G in the same form as Eq. 11, we get or, equivalently, Denoting by λ i the eigenvalues of A, the stability condition reads In this two general cases, the projected adjacency matrix, Eq. 1, arises naturally from the linear stability analysis around the fixed point, which means that only the projected structure is relevant for the dynamics. Moreover, as it can be seen from Eqs. 15 and 19, we can decouple the structural and dynamical contributions to the stability conditions. Note that although the former derivations are specific for these two cases, Eq. (9) holds for a general dynamics in an arbitrary hypergraph. The framework outlined up to now is general and can be readily applied to many dynamical processes in which one is forced to go beyond pairwise interactions. In what follows, we discuss two of such applications, namely, social contagion [8,17] and diffusion processes.
Social contagion. In this dynamical process, the state x i represents the probability of an individual to be active. The deactivation mechanism is given by f i (x i ) = −δx i and the interaction functions g j (x {ej } ) are given by the product of the probability that v i is inactive and the other nodes in the hyperedge e j are active times a contact rate, β |ej | . Hence, Separating the pairwise and higher-order contributions to the state of the node, we have Assuming a symmetric fixed point, i.e., x i = x * , where we can apply Eq. 9 to evaluate the stability of the system's dynamics. The mean-field form of Eq. 21 is given by where ν = max j {|e j |} is the maximum cardinality and c(k) is the ratio between the number of hyperedges with cardinality k and the number of pairwise interactions, which characterizes the structure of the hypergraph. In the steady-state this is a polynomial equation whose solutions are the fixed points of the process and their stability can be evaluated using Eq. 22. A fundamental observation from Eq. 22 is that only pairwise interactions, given by the off-diagonal terms, are responsible for the stability of the fixed point x i = 0. Interestingly enough, in this scenario, the dynamics behaves like an SIS process on complex networks, and is characterized by a continuous phase transition [29]. Therefore, we claim that in a social contagion model, pairwise interactions are a necessary condition for a continuous phase transition. We can support this claim using an argument by contradiction as follows. Near the absorbing state the probability for a node to be active is very small, hence, assuming an arbitrary small number η ≈ 0, we can consider x i ∈ O(η). Therefore, from Eq. 21, the coefficients β |ej | are multiplying terms of order O(η (|ej |−1) ). At the steady state, δx i ∼ β |e|j x i , assuming δ ∈ O(1) without loss of generality. From this point, we can differentiate three cases depending on β |ej | : i If β |ej | ∈ O(1) for all |e j | ∈ {1, 2, ..., max(|e j |)}, then the leading coefficient is β 2 therefore β 2 x i ∈ O(η) as assumed. In this case we have a continuous phase transition with respect to the control parameter β 2 ; ii If β 2 ∈ O(1) and β |ej | ∈ O(η −(|ej |−2) ), for |e j | > 2, then the leading coefficient is β |ej | x i ∈ O(η), thus the transition can be either continuous or discontinuous; iii If β |ej | ∈ O(η −(|ej |−1) ) for all |e j | ∈ {1, 2, ..., max(|e j |)}, then the leading coefficient is β |ej | x i ∈ O(1), which contradicts the initial assumption of x i ∈ O(η). This suggests a discontinuity if β 2 ∈ O(1) and at least one of the higher-order coefficients is O(η −(|ej |−1) ).
To further extend the analysis in [8], in order to compute c(k), we consider higher-order contact patterns from two real datasets. The data were collected with proximity sensors in a high school [14,30], and a primary school [14,31]. Our results are summarized in Fig. 2 (a) and (b), where we consider the feasible fixed points and their stability. In this case, we assume β 2 = 1 and β k = e k , in order to emphasize the role of lower and higher-order structures. Notice that, given the same parameters, we can have either continuous or discontinuous phase transitions with hysteresis, depending only on the structure of the hypergraph. Considering the graph case, i.e., max j {|e j |} = 2, the social contagion model is an SIS epidemic spreading. Thus, the behavior near the fixed point Note that the stability of the absorbing state is given by λ max < δ β2 , as for the quenched mean-field [32]. Next, we analyze the simplicial complex case, where the structure is homogeneous and consists only of pairwise interactions and triangles, hence n = 3. In this case, Eq. 23 has the three solutions, x * 1 = 0, and . (25) From Eq. 22 we derive the following condition thus, x * 1 = 0 is stable for β 2 < δ. The feasibility of the solution x * ± 2 was shown in [8]. Complementarily, Fig. 2 (c) -(f) show the results of the corresponding stability analysis, illustrating how our approach could provide new insights. As noted before, we stress that Eq. 22 is not limited to the case of simplicial complexes, but it can be used to further explore higher-order structures.
Diffusion process. Another immediate application is the diffusion process [2]. In this case, g j (x {ej } ) is separable and g(x i ) = x i . The dynamics of the process is given by where D is the diffusion constant and the minus sign is a consequence of the stability condition in Eq. 15. We can arbitrarily decompose the Laplacian in m-uniform hypergraphs as done for the adjacency matrix in Eq. 3. In particular, we want to consider if the lower-order contributions are enough to describe the dynamics. To this end, we can separate the Laplacian in two components L = L O (θ) + L R (θ) given by where θ differentiates between higher and lower order Laplacian, L R (θ) and L O (θ), respectively. In order to underline the important role played by the structure on the dynamics, we next analyze a diffusion process on both a synthetic hypergraph generated according to a power-law cardinality distribution, P (|e|) ∼ |e| −2.05 , and a real hypergraph using data from email exchange at a European research institution [14,33,34]. In Fig 3 (a) and (b) we show the variance of x O (θ, t) = e −DL O (θ)t x(0), which is the solution of the diffusion process on the lower-order Laplacian L O . For a small enough θ, we have disconnected components and the process does not reach the final state of the complete Laplacian. On the other hand, when the final state is reached, different θ's correspond to different time-scales, and in particular, convergence is faster when increasing θ. In fact, for the synthetic case we need less information (structures) to represent the whole process, in the email dataset the convergence is only reached for θ − 1, as shown in Fig 3 (b). The relative error can be estimated as For the synthetic case, as shown in Fig 3 ( eneough times, different plateau are reached depending on how many components are considered. The same is true for the email case, Fig 3 (d), however we have a multiple peaked pattern, suggesting that we have various structures, each one with a different time scale. As expected, only for θ = 39 the error goes to zero, highlighting that real structures might have non-trivial behavior. Another important aspect is the role on the dynamics of considering different weights for the hyperedges [8,17,18]. In particular, by redefining the Laplacian with the weighted adjacency matrix, Eq. 2, different choices of w(|e j |) can be used, as shown in Fig. 3 (e) -(g). As expected, increasing the weights reduces the convergence time. Aside from the different time-scales in the two examples, we show that real data might also have a qualitatively different behavior, as shown in Fig. 3 (f) and (g). In particular, we focus on the powerlaw weighting functions, w(|e j |) = |e j | −γ , in (g). In this case, since higher-order structures are heavily penalized, we observe a greater time modulation by increasing γ.
Conclusions. In this paper, we have developed a framework that allows to describe the relationship between the structural organization of hypergraphs and general dynamical processes for which including more than pairwise interactions could be relevant. Within this framework, the graph projection appears naturally from a linear stability analysis around a given fixed point. Therefore, the dimensionality of the hypergraph, N max(|ej |) in its tensorial representation, is reduced to a N 2 space. Moreover, we showed that using a structural decomposition in m-uniform hypergraphs, given by Eq. 3, enables the individual analysis of each m-th order component. Altogether, the methodology provides the stability condition for any general dynamical process and also allows to identify the structural components that play a role in the stability of any known fixed point. Of particular interest are the special cases of separable functions and symmetric fixed points, as one is able to separate the structural and dynamical contributions to stability. We also applied the proposed framework to two relevant dynamical processes, namely, social contagion and diffusion dynamics. For the former, we not only recovered known results for graph-based and simplicial complexes, but we obtained a necessary condition for the stability of the absorbing state, showing that that only pairwise interactions play a role for such a state. Additionally, our analysis revealed the conditions for the continuity of the phase transition also for real systems, where different kinds of bifurcations might appear. Finally, in regard to diffusion processes, we characterized the role played by the structural decomposition of the hypergraph and different weighting functions for the hyperedges on the dynamics of the system. In concluding, we stress that the framework developed here is important beyond the interest in dynamical processes, e.g., potential applications can be found in optimization problems. We are confident that this investigation will motivate further research on the spectral properties of hypergraphs, their structural characterization and in providing relevant insights in different areas where higher-order interactions can not be neglected, most notably, in brain dynamics.