Social contagion models on hypergraphs

Our understanding of the dynamics of complex networked systems has increased significantly in the last two decades. However, most of our knowledge is built upon assuming pairwise relations among the system's components. This is often an oversimplification, for instance, in social interactions that occur frequently within groups. To overcome this limitation, here we study the dynamics of social contagion on hypergraphs. We develop an analytical framework and provide numerical results for arbitrary hypergraphs, which we also support with Monte Carlo simulations. Our analyses show that the model has a vast parameter space, with first and second-order transitions, bi-stability, and hysteresis. Phenomenologically, we also extend the concept of latent heat to social contexts, which might help understanding oscillatory social behaviors. Our work unfolds the research line of higher-order models and the analytical treatment of hypergraphs, posing new questions and paving the way for modeling dynamical processes on these networks.

Indeed, modern messaging systems (e.g., WhatsApp, Telegram, and Facebook Messenger, among others) allow users to communicate in groups, which creates a direct channel among all members. In other words, modern information spreading is often a one-to-many process. Additionally, biological and team collaborations are also inherently group structured, similarly to some types of molecular interactions [20]. Complementary, there is evidence from social and biological studies indicating that higher-order structures have crucial dynamical effects [21][22][23][24]. Therefore, higher-order interactions are ubiquitous, and understanding their properties and impacts is of paramount importance. Notably, the sizes of such groups can span orders of magnitude. Thus, graph-projection-based approaches might not be sufficient to describe systems involving interactions over many different scales and orders.
Combinatorial higher-order models [19] offer a way to describe these systems as they overcome the limitations of lower-order network models. In a first attempt, Bodó et al. [21] proposed an SIS disease spreading in a hypergraph. Next, Iacopini et al. [22] presented a model of social contagion defined on simplicial complexes and provided approximate solutions for complexes of order three. Their model presented new phenomenological patterns associated with the critical properties of the dynamics. However, their proposed model is still very constrained, both structurally and dynamically.
Here we extend their model both structurally and dynamically. Structurally, we adopt hypergraphs, which generalize the concept of graphs, by allowing an edge to have an arbitrary number of nodes (see Fig. 1 for an illustration). Hypergraphs relax the structural restrictions required by simplicial complexes as they impose virtually no limitation on the type, size, and mutual inclusion of interactions, thus, representing more faithfully and naturally real systems. From the dynamical viewpoint, we incorporate explicit critical-mass dynamics (each hyperedge is an independent critical-mass process), generalizing the one modeled in Ref. [22]. The resulting model displays a rich complex phenomenology, remaining very flexible and able to cover a wide range of systems. We uncover the presence of discontinuous transitions and bistability led by higher-order interactions and critical-mass dynamics. Notably, these critical properties contrast with contagion models on graphs-models, which instead usually display continuous transitions, e.g., SIS or SIR disease spreading, regardless of the structural configuration. Therefore, assuming a graph projection of a hypergraph might lead to wrong results. Here we report analytical and numerical analyses of the theoretical framework that we introduce as well as results for several limiting cases and hypergraph structures. We round off the paper by discussing several implications of our study and, most notably, the role of critical mass dynamics in social contagion, providing insights that could help explain reported differences in experimental results [24][25][26][27].
Let us first introduce some formal definitions. A hypergraph is defined as a set of nodes, V = {v i }, where N = |V| is 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 |. If max (|e j |) = 2 we recover a graph. We remark that we do not explore this scenario since it recovers standard models. On the other hand, if for each hyperedge with |e j | > 2 its subsets are also contained in E, we recover a simplicial complex (for more on the hypergraph structure, see the Supplemental Material [28], Sec. I). Figure 1 shows an example of a hypergraph. In an arbitrary hypergraph, we associate with each individual v i a Bernoulli random variable Y i (complementary X i ). If the node v i is active Y i = 1 (X i = 0), otherwise Y i = 0 (X i = 1). To each active node, we associate a deactivation mechanism, modeled as a Poisson process with parameter δ i , For each hyperedge, j, we define a random variable T j = k∈e j X k , which is the number of active nodes in the hyperedge. If T j is equal or above a given threshold, j , we associate a Poisson process with parameter λ j , N In other words, the dynamics is given by a threshold process that becomes active only above a critical mass of activated nodes. Moreover, if |e j | = 2, we assume directed Poisson processes, implying that it is not a threshold process anymore. This definition allows for recovering traditional SIS contagion models. While the proposed model is general, allowing for arbitrary heterogeneity in parameters, we focus on more straightforward, but representative, cases. We assume that δ i = δ and λ j = f (|e j |), where f is an arbitrary function of the cardinality of the hyperedge. It is also convenient to define = * N , where the parameter * is a real number representing the fraction of active nodes.
The exact equation that describes the aforementioned dynamics can be written as where the first summation is over all hyperedges containing v i and the second over all the possible dynamical microstates inside the hyperedge e j , denoted by the set B. Furthermore, is an indicator function that is 1 if Y i = 0 and the critical mass in the hyperedge is reached and 0 otherwise. Naturally, the order parameter is defined as the expected fraction of active nodes, i.e., ρ = 1 (1) captures the exact process, it cannot be numerically solved. Thus, assuming that the random variables are independent and denoting y i = E(Y i ), we obtain the firstorder approximation, given as where we assume that the spreading rate is composed by the product of a free parameter and a function of the cardinality, i.e., λ j = λ × λ * (|e j |). In this formulation, the expectation of the indicator function follows a Poisson binomial distribution (for more on this approximation, see Ref. [28], Sec. II). Formally, where F l is the set of all subsets of k integers from {1, 2, . . . , n = |e j |}, A is one of those sets, and A c is its complementary. Intuitively, A accounts for the possibly active nodes and A c the possibly inactive ones. Thus, the summation over F l considers all possible nodal state configurations in a given hyperedge. Equation (4) is not numerically stable if |e j | is large [29]. Considering the discrete Fourier transform, we obtain a numeric stable solution as [29] P e j (K = k) (5) where C = exp ( 2iπ n+1 ), which then allows one to compute the solution for arbitrarily large hyperedges. Interestingly, although the whole argument is quite intricate, Eq. (5) is simple, allowing the numerical evaluation of Eq. (2).
Our main result is that a rich and diverse phase space, generally populated by continuous and discontinuous transitions and hysteretic behaviors, characterizes contagion on hypergraphs. In particular, we have analytically observed discontinuity and bistability in the order parameter on top of some regular structures. We provide full details of the calculations in Ref. [28] (see Secs. III and IV) for two limiting cases, namely, a hypergraph composed of a hyperedge containing all nodes in addition to (1) a random regular network (which we call a hyperblob) and (2) a star (referred to as a hyperstar).
For the sake of clarity, let us show the main results for the hyperblob. In this case, we can exploit the structural symmetries to solve ρ(λ, λ * , δ), obtaining two locally stable solutions. Specifically, consider a hypergraph built up as a homogeneous set of pairwise interactions with average degree k and a single additional hyperedge containing all nodes. In this case, the order parameter can be solved as (for the full derivations, see Ref. [28], Sec. III) where a second-order phase transition for ρ Lower is naturally obtained as λ δ 1 k [7]. Furthermore, the discontinuities can also be calculated as Phenomenologically, a discontinuity implies that our system possesses a "social latent heat" that is released or accumulated at a constant value of λ. More specifically, before the discontinuity, "energy" has been stored in the partial activation of the hyperedges. At the discontinuity this "energy" is absorbed (released) at once for a constant value of λ. In fact, the social latent heat can be expressed as where Q l (λ X c ) can be Q l (λ L c ) (energy absorbed) or Q l (λ U c ) (energy released). Therefore, for this structure, the latent heat is expressed as where λ X c can be (λ L c or λ U c ). In fact, this expression is true for any value of λ, but its physical interpretation is valid only near the discontinuity, which in turn depends on ρ c = * . We refer the reader to Ref. [28] for more details. Figure 2 shows the general phenomenology of the system obtained from the analytical solution, i.e., the first-order approximation, of the equations describing the contagion dynamics for the hyperblob. As can be seen in Fig. 2(a), there are two possible solutions, ρ Lower and ρ Upper . The solution depends on the initial conditions and the threshold, * , which, together with the structure, defines a value ρ c where the dynamics exhibit a discontinuity. If ρ(t = 0) ρ c , the solution is given by ρ = ρ Upper (forward phase diagram). On the other hand, if ρ(t = 0) < ρ c and ρ(t = 0) = 0, then ρ = ρ Lower (backward phase diagram). The arrows show these solutions as the size of the jump (i.e., the magnitude of the latent heat). Note that the lower solution can exhibit a secondorder phase transition. Next, in Fig. 2(b) we instead represent the corresponding parameter space, which is composed of five distinct regions, as explained in the figure caption. We assumed the most general case, where the lower solution has a transition from the absorbing state to an active state, here at λ c . We remark that, for some structures, the lower solution might have a vanishing critical point (e.g., the hyperstar), i.e., λ c → 0, thus slightly changing this picture. Finally, we have compared both analytic and numeric estimates of the latent heat for this hypergraph structure. The results show that the absolute error between analytical and numerical simulations is of order 10 −2 -10 −3 in hypergraphs with N = 10 4 (see Ref. [28], Sec. VI.E and Table I), indicating that the first-order approximation is accurate.
Thus, generically, the solutions for a social contagion dynamics on hypergraphs can be mathematically expressed as where ρ is obtained if ρ(t = 0) < ρ c , and ρ * if ρ(t = 0) ρ c , where ρ c is a global critical mass, i.e., the value of ρ at which the discontinuity appears. As before, we note that the lower solution (branch) might also exhibit a secondorder (continuous) phase transition, denoted by λ c in Fig. 2, from the absorbing state (ρ = 0) to the active state (ρ > 0). Furthermore, for a given hypergraph with fixed δ and λ * , the discontinuity points are formally defined as thus, also defining the bistable region, (λ U c , λ L c ). Although a closed solution for the general case is not possible, Monte Carlo simulations and numerical evaluation of Eq. (2) are reasonable alternatives to characterize our system (see Ref. [28], Sec. VI). Here we focus on a hypergraph with an exponential distribution of cardinalities (i.e., the number of nodes inside a hyperedge), P (|e|) ∼ μ exp (−μ|e|) with the constraint that |e| 2. Dynamically, we set λ * = log 2 (|e|). This choice is arbitrary, but we choose here the log 2 (|e|) function because it grows sublinearly. Note that, if a hyperedge cardinality goes to infinity, the average spreading value tends to zero: lim |e|→∞ log 2 (|e|) |e| = 0. The impact of such a function is yet unknown, and we leave this analysis for future work. Figures 3(a) and 3(b) show that the order parameter and the susceptibility follow the patterns expected for a first-order transition, i.e., both are discontinuous. Moreover, the order parameter is bistable, implying the presence of a hysteresis loop. This phenomenon is opposed to an SIS on a graph. The SIS has a second-order phase transition, characterized by a continuous behavior of the order parameter and a diverging susceptibility in the thermodynamic limit. Complementarily, Fig. 3(c) shows the distribution of active nodes in the upper and lower branches. In the former, we have a bell-shaped distribution, similar to the supercritical regime of an SIS process [7,30]. In the latter, we have a distribution peaked at one, similar to the subcritical regime (absorbing state) of an SIS process [7,30]. We emphasize that Fig. 3(c) displays the distribution of active nodes for the upper (left panel) and lower (right panel) branches and that the complete distribution for a given λ in a region where both solutions exist is bimodal. Intuitively, one expects that lower cardinality hyperedges are responsible for the lower branch as they are easier to activate than the higher cardinality ones. Note that, in regular cases, such as the hyperblob (see Ref. [28], Secs. III and IV), the giant pairwise component has N nodes. However, for the exponentially distributed cardinality scenario, this is not the case. Although P (|e| = 2) P (|e| = k), k = 3, 4 . . . , N (exponential distribution), the largest connected component, is very small, six nodes in the simulated hypergraph. Therefore, in this scenario, the lower branch is determined by hyperedges with cardinality lower but greater than two. The generality of the reported phenomenological behavior suggests that the observed dynamics are a consequence of group-group interactions. Additional experiments (reported in Ref. [28]) for a hypergraph with a power-law distribution of cardinality further corroborated these results. In all systems, we found similar qualitative behavior for finite networks.
Our results are relevant because they provide a theoretical foundation for, and a phenomenological explanation to, seemingly different experimental findings [24][25][26][27][31][32][33]. These works reported critical mass levels needed to change an established equilibrium of 10% in some experiments and 30%-40% in others, in apparent contradiction. The formalism here developed naturally brings forth plausible hypotheses for these observations and show that both ranges are possible. On the one hand, studies based on a single group suggest a threshold between 30%-40%, a situation that can be modeled as a single hyperedge in our formalism. On the other hand, a critical mass of 10% would correspond to a population that is composed of groups of diverse sizes, each one with a (larger) activation threshold. In other words, it is possible to have individual groups exhibiting a threshold * between 30%-40%, and at the same time, a global critical mass, ρ c , for the whole population of about 10% due to group intersections. A second reason that could explain the experimental findings is even more straightforward: admittedly, the fact that our model shows bistability also enables, for a given λ, two possible solutions for ρ corresponding to the lower and the upper branches. That is, the system might be operating in the region where both solutions are larger than zero and stable.
In summary, in this paper, we have developed a framework that allows extending the study of social contagion models when group interactions are relevant. This is achieved by considering hypergraphs as the substrates that capture such many-to-many interactions. First, our work opens the path to deal with new dynamical processes on top of higher-order models and specifically on hypergraphs. Second, we showed that simple dynamical processes could exhibit very rich dynamics, with different transitions, bistability, and hysteresis. Hypergraphs are ubiquitous, and our theory suggests that such a structure allows for the phase diagram reported here. Several findings support the relevance of this methodology. We remark that, depending on the structure, traditional graphprojected models may lead to wrong results. Ultimately, the uncovered phenomenology allows explaining seemingly contradictory experimental findings in which group interactions play a major role. We also mention that many interesting questions arise from our work. For instance, if one assumes that energy is proportional to ρ, our model might display phenomena reminiscent of a Carnot cycle for social contexts, which might help to understand abrupt changes and oscillatory patterns in social behaviors.
G.F.A. thanks E. Artiges, H. F. de Arruda, J. P. Rodriguez, L., Gallo and T. Peron for fruitful and inspiring discussions. G.P. acknowledges support from Compagnia San Paolo (ADnD project). Y.M. acknowledges partial support from the Government of Aragon, Spain, through Grant E36-17R (FENOL), and by MINECO and FEDER funds (FIS2017-87519-P). G.F.A., G.P., and Y.M. acknowledge support from Intesa Sanpaolo Innovation Center. Research carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (Grant 2013/07375-0). The funders had no role in study design, data collection, and analysis, decision to publish, or preparation of the manuscript.