The dynamics of network social contagion processes such as opinion formation and epidemic spreading are often mediated by interactions between multiple nodes. Previous results have shown that these higher-order interactions can profoundly modify the dynamics of contagion processes, resulting in bistability, hysteresis, and explosive transitions. In this paper, we present and analyze a hyperdegree-based mean-field description of the dynamics of the susceptible–infected–susceptible model on hypergraphs, i.e., networks with higher-order interactions, and illustrate its applicability with the example of a hypergraph where contagion is mediated by both links (pairwise interactions) and triangles (three-way interactions). We consider various models for the organization of link and triangle structures and different mechanisms of higher-order contagion and healing. We find that explosive transitions can be suppressed by heterogeneity in the link degree distribution when links and triangles are chosen independently or when link and triangle connections are positively correlated when compared to the uncorrelated case. We verify these results with microscopic simulations of the contagion process and with analytic predictions derived from the mean-field model. Our results show that the structure of higher-order interactions can have important effects on contagion processes on hypergraphs.

Including group interactions in network models of contagion can significantly affect epidemic behavior. By studying the susceptible–infected–susceptible epidemic model on networks with higher-order interactions, we observe that for certain parameters, there is a bistable regime, where above a critical number of infected individuals, the contagion spreads until it becomes an epidemic, and below this critical number, the epidemic dies out. We find that heterogeneity in the individual and group contact structure of a social network determines the existence of such tipping point events and derive conditions for their appearance. Last, we comment on how three group contagion mechanisms—collective contagion, infection by individuals, and the “hipster effect”—affect the onset of epidemics and the existence of bistability.

The study of contagion processes is a fundamental problem in network science, with applications including epidemics,1–7 social media,8 opinion formation,9 idea diffusion,10,11 sudden changes in social convention,12,13 and many more. Contagion processes can be of many types, ranging from discrete-state models such as the susceptible–infected–susceptible (SIS) model, to continuous models of opinion formation, to realistic models of disease such as those currently used to model the spread of COVID-19.14,15 Modeling the dynamics of such processes on pairwise interaction networks has been a hallmark of network science, providing many insights into the effect of network structures on the propagation of disease and information. Recently, the role of complex contagion mechanisms (i.e., contagion processes that cannot be described solely by pairwise interactions) has received much attention.16 It has been shown that higher-order interactions in networks (i.e., interactions involving multiple nodes) can have profound effects on dynamical network processes17 such as opinion formation,18 synchronization,19–21 and population dynamics.22 Efforts to map higher-order interactions in real-world networks have an uncovered rich structure,23 which is only now starting to be appreciated. In the context of contagion processes, it was recently shown24 that the addition of higher-order interactions to the SIS epidemic model on Erdös–Rényi networks results in bistability, hysteresis, and explosive transitions to an endemic disease state (see also Refs. 25–28). The simplicial SIS model has also been extended to scale-free uniform hypergraphs.29 The fact that the network SIS model with more general higher-order interactions results in bistability has been proven rigorously in Ref. 26. However, so far, there is no general theory explaining how heterogeneity and correlations in the structure of higher-order interactions affect the onset of bistability.

In this paper, we present and analyze a degree-based mean-field description of the dynamics of the SIS model in networks with higher-order interactions. To describe higher-order interactions, we consider the SIS model on a hypergraph, formed by a set of nodes and a set of edges of multiple sizes (so that edges of size larger than two represent higher-order interactions). Our formulation allows us to consider heterogeneous structures in the organization of the edges of a given size and correlations between the structure of edges of different sizes. Using the illustrative case of networks with edges of sizes 2 and 3, we derive conditions for the appearance of bistability and hysteresis in terms of moments of the degree distribution of the pairwise interaction network. We find that the onset of bistability and hysteresis can be suppressed by heterogeneity in the pairwise interaction network and promoted by positive correlations between the number of pairwise and higher-order interactions a node has. We also consider the effect of healing by higher-order interactions (a “hipster effect”).

The structure of the paper is as follows. In Sec. II, we present our hypergraph and contagion models. In Sec. III, we derive a mean-field description of the model and apply it to various illustrative cases. In Sec. IV, we study how model parameters affect the onset of bistability. Finally, we discuss our results and present our conclusions in Sec. V.

In this section, we present our hypergraph and contagion models. Our model consists of SIS contagion spreading on a hypergraph via pairwise and higher-order interactions. While we focus on the SIS epidemic model, we note that our formalism could be extended to other models. In the context of epidemic spreading, pairwise interactions could represent, for example, face-to-face interactions leading to contagion via viral droplets, while higher-order interactions could represent, for example, contagion via the shared spaces by a group. In the context of opinion dynamics, higher-order contagion could model, for example, a majority-vote process common in caucusing. In the following, we provide details about the hypergraph model representing the higher-order interactions and the contagion models that we consider.

We consider a population of N nodes labeled i=1,2,,N coupled via undirected hyperedges of sizes m=2,3,,M, where a hyperedge of size m is a set of m nodes, {i1,i2,,im}. We define the mth order degree of node i, ki(m), as the number of hyperedges of size m to which the node belongs, and its hyperdegree as the vector ki=[ki(2),ki(3),,ki(M)]. The 2nd order degree of a node corresponds to the number of pairwise connections of the node, while higher-order degrees measure the node’s participation in hyperedges of larger sizes. Figure 1 illustrates a hypergraph with hyperedges of sizes 2 and 3, which, for simplicity, we will henceforth denote as links and triangles, respectively.

FIG. 1.

Illustration of a hypergraph. Infected nodes (red) infect a healthy node (gray) via hyperedges of sizes 2 and 3 with rates β2 and β3, respectively.

FIG. 1.

Illustration of a hypergraph. Infected nodes (red) infect a healthy node (gray) via hyperedges of sizes 2 and 3 with rates β2 and β3, respectively.

Close modal

Extending degree-based descriptions of epidemic spreading on networks,30,31 we will develop a mean-field theory for the propagation of epidemics based on the assumption that nodes with the same hyperdegree have the same statistical properties. For this purpose, we assume that the number of nodes with the hyperdegree k, P(k), is given and that the probability that nodes with hyperdegrees k1, k2,…, km belong to a hyperedge of size m is given by fm(k1,k2,,km). This assumes that the statistical structure of the network is completely described by the hyperdegree distribution P(k) and the connection probabilities fm(k1,k2,,km). While this restriction rules out the possibility of assortative mixing by other node properties, it is straightforward to extend our formalism to include other node variables. Note that counting the number of hyperedges of size m in two different ways, the connection probabilities must be normalized such that

(1)

For example, for the configuration model for networks without higher-order interactions (i.e., only hyperedges of size 2, M=2), the hyperdegree of a node is just the number of links, k=k, connecting that node to other nodes and the connection probability is f2(k,k)=kk/(Nk), where k=i=1Nki/N=kkP(k)/N. For networks with hyperedges of sizes 2 and 3, f3(k1,k2,k3) is the probability that three nodes with degrees k1, k2, and k3 are connected by a hyperedge of size 3. The configuration model for hypergraphs and its associated statistical properties has been studied in Refs. 32 and 33.

This framework allows us to study networks with heterogeneously distributed higher-order interactions and correlations between nodal degrees of different orders. In addition, it allows us to treat the case in which nodes belonging to a triangle are not necessarily connected by links, as is assumed in simplicial complex models.24 We will study how the structure of higher-order interactions modifies some of the properties of epidemic spreading on networks with exclusively pairwise interactions (i.e., hyperedges of size 2 only), on which epidemic spreading has been studied extensively.1 

Now, we describe the contagion models we will study. As mentioned above, we will focus on the SIS model, but other epidemic models could be treated using the same formalism. We assume that at any given time t0, each node can be in either the susceptible (S) or infected (I) state. Infected nodes heal and become susceptible again at rate γ. Now, we specify how hyperedges mediate the contagion process. In general, the probability of contagion by a hyperedge could be a function of the number of infected nodes in the hyperedge (e.g., as in Ref. 27). Here, we will consider the two extreme cases where contagion occurs if all the other members of the hyperedge are infected or if at least one member of the hyperedge is infected. More precisely, in the collective contagion case, a susceptible node that belongs to a hyperedge of size m gets infected at rate βm if all the other members of the hyperedge are infected; in the individual contagion case, the node gets infected at rate βm if at least one member is infected. While we will analyze these two cases only, in principle, one could treat the case in which at least j other nodes of the hyperedge need to be infected for contagion to occur using the techniques presented below. This case corresponds to a quorum of size j, and there is evidence for such effects in collective behavior.34,35 For hyperedges of size 2, i.e., links, both cases reduce to the usual contagion via pairwise interactions. The social contagion model of Ref. 24 corresponds to the collective contagion case. The contagion processes are illustrated in Fig. 1 for hyperedges of sizes 2 and 3. Table I summarizes the notation and variables used.

TABLE I.

Relevant notation.

VariableDefinition
N Number of nodes 
k(m) Number of hyperedges of size m a node 
 belongs to 
k = [k(2), …, k(M)Hyperdegree 
P(kNumber of nodes with hyperdegree k 
γ Rate of healing 
βm Rate of infection by a hyperedge of size m 
fm(k1,k2,,km) Probability that m nodes form a hyperedge 
 of size m 
xk Fraction of nodes with hyperdegree k that 
 are infected 
VariableDefinition
N Number of nodes 
k(m) Number of hyperedges of size m a node 
 belongs to 
k = [k(2), …, k(M)Hyperdegree 
P(kNumber of nodes with hyperdegree k 
γ Rate of healing 
βm Rate of infection by a hyperedge of size m 
fm(k1,k2,,km) Probability that m nodes form a hyperedge 
 of size m 
xk Fraction of nodes with hyperdegree k that 
 are infected 

In this section, we present a mean-field analysis of the epidemic dynamics on a network specified by the hyperdegree distribution P(k)/N and the hyperedge connection probabilities fm(k1,k2,,km). Assuming that all nodes with the same hyperdegree behave similarly, we focus on xk, the fraction of nodes with hyperdegree k that are infected. The mean-field equation describing the evolution of xk is

(2)
(3)

The first term on the right-hand side of Eq. (2) corresponds to healing at rate γ and the second term accounts for infection by hyperedges. The number of hyperedges of size m that can pass an infection to a node with hyperdegree k is calculated by considering all the possible hyperdegrees of the other m1 nodes participating in the hyperedge (k,k1,,km1), counting how many such combinations there are not counting permutations [P(k1),,P(km1)/(m1)!], calculating what fraction of such combinations forms a hyperedge with the node in consideration [fm(xk,xk1,,xkm1)], multiplying by the probability that the hyperedge can transmit the infection [G(xk1,,xkm1)], and summing over all hyperdegree combinations. The probability that the hyperedge can transmit the infection, given by (3), depends on whether the collective contagion or individual contagion model is assumed. Note that the form for G taken above, and the mean-field treatment in general, assumes that the states of nodes are independent. A better approximation that includes correlations between connected nodes has been implemented in Ref. 28 for the case of unstructured hyperedges of sizes 2 and 3, leading to an improved quantitative agreement with the results of numerical simulations. Since our interest is in the effects of higher-order structures on qualitative aspects of the epidemic dynamics, we will use the mean-field approximation in Eq. (2). A similar mean-field equation for a node-based description of the contagion process was recently formulated in Ref. 26. In the following, we will apply the mean-field description to illustrative cases.

Here, we focus on the case where the hyperedge sizes are either 2 or 3; i.e., M=3. This corresponds to a network like in Fig. 1, with hyperedges of size 2 (links) and 3 (triangles). For simplicity, here, we denote the number of links per node as k, i.e., k=k(2), and the number of triangles a node belongs to by q, i.e., q=k(3). In addition, we will consider the case where the connection probabilities depend only on the node links, i.e., fm(k,k1,,km1)=fm(k,k1,,km1). With these assumptions and using the collective contagion rule in Eqs. (3) and (2) becomes

(4)

where the first term on the right-hand side represents healing, the second represents contagion by links, and the third represents contagion by triangles.

Since the connection probabilities do not depend on q, we can reduce the dynamics to the fraction of nodes with degree k that are infected,

(5)

where P(k)=qP(k,q) is the number of nodes with degree k. Multiplying Eq. (4) by P(k,q), summing over q and dividing by P(k), we obtain

(6)

For the link connection probability f2(k,k1), we will take f2(k,k1)=kk1/(Nk), which corresponds to nodes being connected completely at random according to their degree as in the configuration model. For the triangle connection probability f3, we will consider two cases: the uncorrelated case and the degree-correlated case. In the degree-correlated case, we assume that the connection probability is given by f3(k,k1,k2)=2kk1k2/(Nk)2 so that nodes that have a higher number of links also belong to more triangles. In the uncorrelated case, we assume instead that f3(k,k1,k2)=2k/N2 so that triangles are formed independent of the nodal degrees. The normalization is chosen using Eq. (1) so that the mean number of triangles per node, q=i=1Nki(3)/N, in each case is equal to k. We note that the model for triangle formation in Ref. 24 corresponds to the uncorrelated case. We can choose the mean triangle degree independent of the mean network degree by scaling f3(k,k1,k2) by q/k, but for simplicity, we assume q=k. Figure 2 illustrates the difference between the two cases in a small network, where in the degree-correlated case, the triangles cluster around nodes with a high pairwise degree and in the uncorrelated case, the triangles are distributed uniformly at random on the network.

FIG. 2.

Schematic illustration of the degree-correlated and uncorrelated cases. In the degree-correlated case (left), nodes with more links are more likely to belong to a triangle. In the uncorrelated case (right), triangles connect nodes with a probability independent of their degree.

FIG. 2.

Schematic illustration of the degree-correlated and uncorrelated cases. In the degree-correlated case (left), nodes with more links are more likely to belong to a triangle. In the uncorrelated case (right), triangles connect nodes with a probability independent of their degree.

Close modal

We can also specify the distribution of triangle degrees by defining f2(q,q1) and f3(q,q1,q2) and then reducing Eq. (4) by multiplying by P(k,q), dividing by P(q), and summing over k to reduce the dynamics to the fraction of infected nodes with triangle degree q. For the triangle connection probability, we take f3(q,q1,q2)=2qq1q2/(Nq)2, which corresponds to three nodes being connected at random according to the configuration model for triangles.32 For the pairwise links, we define the degree-correlated and uncorrelated cases as before, where in the degree-correlated case, f2(q,q1)=qq1/(Nq) and for the uncorrelated case, f2(q,q1)=q/N. From there, we can use the same formalism as our approach when specifying the pairwise degree.

Now, we consider separately the degree-correlated and uncorrelated cases. In the correlated case, where f3(k,k1,k2)=kk1k2/(Nk)2, Eq. (6) can be rewritten in terms of the fraction of infected links

(7)

as

(8)

In this case, the dynamics of nodes of degree k is determined by the global variable V. To study the qualitative characteristics of the dynamics, we find the steady-state solutions. The fixed point of Eq. (8) is

(9)

Inserting this in (7), we obtain a nonlinear equation that determines the fraction of infected links V,

(10)

The state with no infection, V=0, is a solution to (10). However, it is linearly unstable for β2>β2c=γk/k2, as can be seen by linearizing Eq. (8) about V=0, multiplying by kP(k)/(Nk), and summing over k, which yields the linearized equation for the evolution of the perturbation δV

(11)

The nonzero solutions of Eq. (10) represent states with a nonzero fraction of infected nodes.

Now, we study the uncorrelated case where f3(k,k1,k2)=2k/N2. In this case, Eq. (6) can be rewritten in terms of the fraction of infected nodes,

(12)

and the fraction of infected links V. In terms of these quantities, Eq. (6) reads

(13)

As in the prior case, the equilibrium is

(14)

Evaluating this expression in Eqs. (7) and (12), we obtain the coupled equations,

(15)
(16)

The state with no infection, U=0, V=0, is a solution of (15) and (16). By considering perturbations δU, δV from this solution, linearizing Eq. (13), and evaluating in Eq. (7) for the first equation and Eq. (12) for the second equation, we obtain the linear system

(17)
(18)

which shows that the no infection state is linearly unstable for β2>γk/k2, which is the same threshold we obtained for the correlated case.

In summary, nonzero solutions of Eq. (10) and Eqs. (15) and (16) for the degree-correlated and uncorrelated cases, respectively, represent states with a nonzero number of infected nodes. Figure 3 shows the fraction of infected nodes U for the uncorrelated case as a function of the normalized pairwise infectivity β2/β2c for three values of the triangle infectivity β3 obtained from numerical solution of Eqs. (15) and (16) with P(k)k4 for 67<k<1000 and 0 otherwise. Different solutions are plotted as solid and dashed lines to indicate stability or instability, respectively. The connected circles are obtained from numerical simulations of the full stochastic microscopic model. In these simulations, β2 was slowly increased in small steps up to a maximum value and subsequently decreased back to its initial value. For each β2, the average number of infected nodes after transient effects disappeared is shown as a filled circle. For more details about the simulations, see the  Appendix.

FIG. 3.

Fraction of infected nodes U vs link infectivity β2 obtained from the mean-field equations (15) and (16) (solid and dashed lines) and from microscopic simulations (connected circles) using P(k)k4 on [67,1000], γ=2, and N=10000 for β3=0.0194 (a), 0.0388 (b), and 0.05482 (c). Refer to the text for an explanation of the discrepancy between the mean-field equations and microscopic simulations.

FIG. 3.

Fraction of infected nodes U vs link infectivity β2 obtained from the mean-field equations (15) and (16) (solid and dashed lines) and from microscopic simulations (connected circles) using P(k)k4 on [67,1000], γ=2, and N=10000 for β3=0.0194 (a), 0.0388 (b), and 0.05482 (c). Refer to the text for an explanation of the discrepancy between the mean-field equations and microscopic simulations.

Close modal

The behavior of the microscopic simulations is captured qualitatively by the mean-field equations. The quantitative disagreement is likely due to the assumptions inherent to the mean-field approximation. In fact, Ref. 28 has shown that, for the particular case of uncorrelated triangles on an Erdös–Rényi network, the disagreement almost disappears when pair correlations are taken into account. Since our interest in this paper is on the qualitative dynamics, we use the mean-field theory, but note that the approaches proposed in Refs. 28 and 36 could be used to obtain better approximations. The qualitative aspects of interest, captured by the mean-field equations and the numerical solution of Eqs. (15) and (16), are the following. For small values of β3 [Fig. 3(a), β3=0.0194], the bifurcation from the state with no infection (U=0) to the infected state (U>0) is continuous. However, for larger values of β3 [Fig. 3(c), β3=0.0582], the transition is discontinuous: as β2 increases past a critical value β2c, the fraction of infected links increases explosively toward an epidemic equilibrium (upward arrow). If β2 is subsequently decreased, the fraction of infected links remains high until β2 decreases past the value at which the epidemic equilibrium solution disappears and then it decreases to zero (downward arrow). For such values of β3, there is hysteresis, bistability, and explosive transitions. At a critical value β3=β3c, which will be the focus of our interest, there is a transition from the type of bifurcation shown in Fig. 3(a) to the type of bifurcation shown in Fig. 3(c). Figure 3(b) shows U as a function of β2 for a value β3=0.0388β3c. We are interested in exploring how the presence of this bistable regime is affected by the degree distribution P(k) and other parameters of the model, in particular, the triangle infectivity, β3.

Figure 4 shows the phase diagram in the (β2,β3) plane for the degree-correlated, collective contagion model. The plot was obtained by counting the number of solutions of Eq. (10) as a function of β2 and β3 for γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise (all subsequent phase diagram plots are calculated using the same parameters). Light pink indicates that there is only the solution V=0 corresponding to a stable state with no contagion. Orange indicates two solutions, the unstable V=0 solution and another stable solution with V>0. Finally, dark red indicates a bistable regime with three solutions: the stable V=0 solution and a pair of stable and unstable solutions with positive V. As noted in Refs. 24 and 26, this regime is only present for large enough triangle infectivity, i.e., for β3>β3c. The phase space for the uncorrelated case (not shown) is qualitatively similar to the one in Fig. 4, but the transition to bistable behavior occurs at a larger value of β3.

FIG. 4.

Phase diagram for the degree-correlated, collective contagion model. The light pink region labeled “No infection” corresponds to 1 solution of Eq. (10), the orange region labeled “Infection, no bistability” to 2 solutions, and the region labeled “Bistability” to 3 solutions. The parameters are γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise.

FIG. 4.

Phase diagram for the degree-correlated, collective contagion model. The light pink region labeled “No infection” corresponds to 1 solution of Eq. (10), the orange region labeled “Infection, no bistability” to 2 solutions, and the region labeled “Bistability” to 3 solutions. The parameters are γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise.

Close modal

To quantify how the onset of bistability depends on the hypergraph parameters, we define the bistability indexB(β3) as the maximum separation, over all values of β2, between the largest and smallest stable solutions for the fraction of infected nodes U. The bistability index can be calculated from microscopic simulations of the contagion process such as those used to produce Fig. 3 or from the numerical solution of Eq. (10) for the correlated case and Eqs. (15) and (16) for the uncorrelated case. In Fig. 5, we plot the bistability index B as a function of β3 computed from microscopic simulations for three choices of the link degree distribution P(k), all with a mean degree of 100: (a) P(k) constant for 50<k<150 and 0 otherwise, (b) P(k)k4 for 67<k<1000 and 0 otherwise, and (c) P(k)k3 for 53<k<1000 and 0 otherwise. For each distribution, we considered the uncorrelated case (orange connected circles) and the degree correlated case (blue connected triangles). The dashed lines indicate the value β3c at which we expect the onset of bistability, obtained from the numerical solution of Eqs. (10) and (15)(16) for the degree-correlated and -uncorrelated cases, respectively (in Sec. IV, we provide analytical expressions for these values). As the degree distribution of the pairwise interaction network P(k) becomes more heterogeneous from (a) to (c), the value of β3 at which the onset of bistability occurs increases for the uncorrelated case, while it remains almost unchanged for the degree-correlated case. A heuristic interpretation of this phenomenon is the following: in the uncorrelated case, the triangle interactions do not depend on the heterogeneity of the link degree distribution. Therefore, as the link degree distribution P(k) becomes more heterogeneous, contagion becomes dominated by hubs of the pairwise interaction network, a mechanism which does not result in bistability. Therefore, bistability is suppressed in the uncorrelated case. On the other hand, for the degree-correlated case, both triangle and link contagion mechanisms increase their effectiveness in tandem as the heterogeneity of the link degree distribution is increased. It is important to note that the increase in β3c with heterogeneity, which is shown here in absolute terms, still occurs if one considers it relative to the value of β2c (i.e., β3c/β2c also increases with heterogeneity), as we will show later.

FIG. 5.

Bistability index B as a function of β3 for (a) P(k) constant for 50<k<150 and 0 otherwise, (b) P(k)k4 for 67<k<1000 and 0 otherwise, and (c) P(k)k3 for 53<k<1000 and 0 otherwise. For each distribution, we considered the uncorrelated case (orange connected circles) and the degree correlated case (blue connected triangles). The dashed lines indicate the value β3c at which we expect the onset of bistability, obtained from numerical solution of the mean-field equations (12) and (15)(16).

FIG. 5.

Bistability index B as a function of β3 for (a) P(k) constant for 50<k<150 and 0 otherwise, (b) P(k)k4 for 67<k<1000 and 0 otherwise, and (c) P(k)k3 for 53<k<1000 and 0 otherwise. For each distribution, we considered the uncorrelated case (orange connected circles) and the degree correlated case (blue connected triangles). The dashed lines indicate the value β3c at which we expect the onset of bistability, obtained from numerical solution of the mean-field equations (12) and (15)(16).

Close modal

Another interesting aspect seen in Fig. 5 is that the transition to bistable behavior seems sharper in the uncorrelated case for the more heterogeneous networks. As we will see in Sec. IV, the nature of the bifurcation is indeed different for the uncorrelated case and heterogeneous networks.

Finally, we have to point out that the numerical calculation of the bistability index from numerical simulations can be challenging. When the unstable solution is small, finite-size effects can cause transitions to the nonzero stable solution from the stable zero solution, making the numerical determination of the stable fixed points difficult and the bistability index plots noisy. Nevertheless, the mean-field theory predicts well the onset of bistability.

Now, we consider the case of individual contagion, in which an m-hyperedge infects a susceptible node with rate βm when at least one member of the hyperedge is infected. For simplicity, we will still consider only links and triangles (M=3) with infection rates of β2 and β3, respectively.

The analog to Eq. (6) for the individual contagion case is

(19)

For the correlated case, f3(k,k1,k2)=2kk1k2/(Nk)2, this can be rewritten as

(20)

with a fixed point

(21)

Inserting this into Eq. (7) like before, we obtain

(22)

Linearizing about the V=0 equilibrium, we find that the epidemic threshold is given by the condition

(23)

which defines a linear relationship between β2 and β3 for fixed γ, in contrast to the collective contagion mechanism that does not alter the epidemic threshold β2c=γk/k2. This relationship can be understood heuristically by noting that, close to the V=0 solution, the probability that two nodes in a hyperedge are simultaneously infected can be neglected. Under that assumption, infection of a susceptible node by a triangle when at least one other node is infected is equivalent to independent infection by either of the two other nodes in the triangle with rate β3. Since in the correlated case a node belongs, on average, to the same number of links and triangles, the individual contagion model reduces to the traditional SIS model with contagion rate β2eff=β2+2β3 in the linear regime (we emphasize, however, that the nonlinear behavior can be different).

In Fig. 6, we plot the (β2,β3) phase space for this scenario, with light pink indicating one solution (V=0) to Eq. (22) and orange indicating two solutions, the unstable V=0 solution and a stable V>0 solution.

FIG. 6.

Phase diagram for the degree-correlated, individual contagion model with parameters γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise.

FIG. 6.

Phase diagram for the degree-correlated, individual contagion model with parameters γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise.

Close modal

Considering the uncorrelated case where f(k,k1,k2)=2k/N2 and expressing Eq. (19) in terms of U and V, we obtain

(24)

with equilibrium

(25)

which has different first-order behavior than the degree-correlated case. Inserting this expression into Eqs. (7) and (12), we obtain

(26)
(27)

Linearizing, we obtain the system

(28)
(29)

Solving this system and canceling the zero solution, we find that the epidemic threshold is defined by a non-linear relationship between the three epidemic parameters,

(30)

This relationship implies that there is a singularity when β3=β3=γk2/[2(k2k2)k]. However, one can check that β2 is negative at β3=β3, and therefore, the singularity is not physically relevant. Note that when k2=k2 in the case of a k-regular network, the threshold reduces to that of the degree-correlated case.

Here, we consider the effect of higher-order healing for both collective and individual contagion. By higher-order healing, we refer to a situation where infected nodes that belong to a hyperedge of size m>2 with other infected nodes heal at rate βm. This can be thought of as a “hipster effect” where if an idea or trend is popular in groups, then this makes an individual less likely to adopt the trend, but the individual can be convinced to adopt the trend by their pairwise connections.37 For both the collective and individual contagion cases, we comment on the existence of bistability based on numerical phase plots.

When the contagion is collective, the model including higher-order healing can be written as Eq. (6) with the sign of the third term changed, and because the triangle healing mechanism is solely higher-order, there is no effect on the epidemic threshold, which is obtained by the linearization of the 0 solution. However, we find that explosive transitions do not occur for β2,β30.

Likewise, for the individual contagion model, higher-order healing can be written as Eq. (19) with the third term negative. In this case, the epidemic threshold for both the degree-correlated and uncorrelated case can be obtained by substituting β3 for β3 in Eqs. (23) and (30), respectively. Higher-order healing in individual contagion enables explosive transitions to occur for ranges of β2,β30, as can be seen in Fig. 7, which shows the phase space (β2,β3) for the degree-correlated case. As one might expect, for large enough higher-order healing β3, there is no infection, but there is a narrow band of bistable behavior separating the regions of no infection and monostable infection.

FIG. 7.

Phase diagram for the degree-correlated, higher-order healing with individual contagion with parameters γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise.

FIG. 7.

Phase diagram for the degree-correlated, higher-order healing with individual contagion with parameters γ=2 and P(k)k4 when 67<k<1000 and 0 otherwise.

Close modal

So far, we have considered hypergraphs with hyperedges of sizes 2 and 3 only. We now briefly discuss contagion in networks with hyperedges of all sizes; i.e., M=N. In the context of epidemic spreading, hyperedges could be interpreted as participation in social events such as parties, conferences, concerts, and sports events. For simplicity, we will focus on a hypergraph with degree-correlated hyperedges where

(31)

such that the average number of hyperedges of size m a node belongs to is k(m). In this case, by repeating the calculations of Sec. III B, the fraction of infected nodes of degree k evolves in terms of the fraction of infected edges V(7) as

(32)

for individual contagion and

(33)

for collective contagion. In the case of collective contagion, larger hyperedges can cause the emergence of new stable fixed points, which can lead to richer consensus dynamics.24 We focus, however, on the case of individual contagion. Linearizing, we find that the solution xk=0 becomes unstable when

(34)

If the sum yields a value larger than γk/k2, propagating social contagion will result. Social event restrictions implemented as a truncation of the series by prohibiting events larger than a certain size or practices that reduce contagion in social events and reduce βm (such as enforcing physical separation) can reduce the value of the sum so that contagion does not propagate.38,39

In Sec. III, we expressed the epidemic threshold β2c in terms of moments of the degree distribution of the underlying network structure. Similarly, we would like to express the critical value of β3 at which the explosive transitions appear, β3c, as a function of hypergraph structure. Explosive transitions and bistability occur when there are two stable steady-state solutions to Eq. (6). For the degree-correlated and uncorrelated cases, this occurs when there are two non-zero solutions to Eq. (10) and the coupled system of Eqs. (15) and (16), respectively. We can compute the critical value of β3 by finding the numerical solution of these mean-field equations and determining the value of β3c at which bistability appears. This method is much more efficient than using stochastic microscopic simulations of the contagion model to infer the onset of explosive transitions and to map the phase space. Figure 8 shows the predicted value of β3c normalized by β2c for the correlated (a) and uncorrelated (b) cases as a function of the power-law exponent r and the maximum degree kmax, where Eqs. (10) and (15)(16) were solved using P(k)kr if 50kkmax and P(k)=0 otherwise. Larger values of r and kmax correspond to larger heterogeneity of the degree distribution. We note that for the most homogeneous network—the k-regular network—β3c/β2c is 1, and we see in Figs. 8(a) and 8(b) that β3c increases relative to β2c as r or kmax increases except for small values of r and large values of kmax in the degree-correlated case. Thus, heterogeneity in the degree distribution of the pairwise interaction network appears to suppress explosive transitions. However, this effect is much more pronounced for the uncorrelated case (b) than that for the degree-correlated case (a), as we discussed previously. In the  Appendix, we describe in more detail the algorithm employed to find β3c from the mean-field equations.

FIG. 8.

β3c/β2c as a function of power-law distribution parameters for the degree-correlated case (a) and the uncorrelated case (b). β3c was calculated numerically from the mean-field equations (see the  Appendix) and β2c=γk/k2. The parameters are P(k)kr if 50kkmax and P(k)=0 otherwise and γ=2.

FIG. 8.

β3c/β2c as a function of power-law distribution parameters for the degree-correlated case (a) and the uncorrelated case (b). β3c was calculated numerically from the mean-field equations (see the  Appendix) and β2c=γk/k2. The parameters are P(k)kr if 50kkmax and P(k)=0 otherwise and γ=2.

Close modal

Although this method works well in predicting the value of β3c, it does not provide a direct relationship between the network structure and the onset of explosive transitions and is more computationally expensive than an analytical expression. For this reason, we present closed form approximations to β3c and describe the parameter regimes over which they are accurate. Starting with the degree-correlated case and canceling the zero solution of Eq. (10), we find conditions under which there are at least two solutions to

(35)

First, note that h(0,β2)=β2/β2c1 and that h(1,β2)<0. Therefore, if hV(0,β2c)>0, then by continuity, there will be at least two solutions for β2 less than, but sufficiently close to, β2c. This condition gives

(36)

which works well in predicting the onset of bistability for the degree-correlated case. The relative error with respect to the value obtained from directly solving Eq. (10) for all distributions tested is less than 2% (not shown).

The analysis for the degree-correlated case was based on the behavior of h(V,β2) near V=0. For the uncorrelated case, however, we find that a saddle-node bifurcation can occur at positive values of V, and it is necessary to expand Eqs. (15) and (16) to higher order.

Expanding Eqs. (15) and (16) to second order, setting β2=β2c=γk/k2, and subtracting the two equations yield

(37)

which, when evaluated in

(38)

and expanded to fourth order, again setting β2=β2c, yields

(39)

where

(40)
(41)
(42)

For continuous transitions to epidemics, there is only one equilibrium for V at β2=β2c, namely, V=0. The onset of bistability occurs when a second solution appears, which corresponds to the first appearance of a root of (39) in the interval (0,1). Such a root can appear at V=0 in a transcritical bifurcation or at V>0 as a pair of roots in a saddle-node bifurcation. A pair of roots appears when the discriminant of the quadratic equation a0+a1V+a2V2=0 is zero. However, this bifurcation is physically meaningless if it occurs for values of V outside the interval [0,1]. Therefore, we impose the constraint that the value of β3 found by solving a124a0a2=0 must satisfy the inequality 0a1/2a21. In addition, we note that because of continuity, the sign of the a2 term must be negative because otherwise, hV(0,β2c)>0 and the bifurcation has already occurred. The transcritical bifurcation occurs when a root crosses from a negative value to a positive value, which occurs when one root of a0+a1V+a3V2=0 is V=0, implying that a0=0 and β3c=γk3/k4. Using these conditions, we can construct a piecewise definition of β3c

(43)

The relative error in the value of β3c/β2c obtained from Eq. (43) compared with the numerically obtained value shown in Fig. 8(b) is shown in Fig. 9. In principle, one can expand to higher order to gain accuracy for the most heterogeneous of distributions. However, there is limited utility in increasing the order of the expansion further because the resulting conditions become extremely complicated.

FIG. 9.

Relative error in the value of β3c/β2c obtained from Eq. (43) compared with the numerically obtained value shown in Fig. 8(b).

FIG. 9.

Relative error in the value of β3c/β2c obtained from Eq. (43) compared with the numerically obtained value shown in Fig. 8(b).

Close modal

In this paper, we studied the SIS model of social contagion on hypergraphs with a heterogeneous structure. The mean-field description in Eq. (2) allowed us to explore the effects of hyperedge organization on the epidemic onset and the onset of bistability and explosive transitions. One of our main findings is that with increasing heterogeneity of the pairwise network degree distribution, the onset of explosive transitions is postponed when the pairwise and higher-order interactions have an independent structure. More generally, when considering a hypergraph contagion model, the group infection and the pairwise infection are competing mechanisms by which contagion spreads. Factors that promote contagion via pairwise infection, such as a heterogeneous degree distribution of the pairwise contact network, suppress discontinuous transitions. Conversely, heterogeneity in the degree distribution of hyperedges of higher order promotes such transitions.

We considered two ways in which the structure of hyperedges of different sizes could be organized: the uncorrelated case, in which they are independent, and the correlated case, in which hyperedges of different sizes connect preferentially to the same nodes. While the organization of hyperedges in real world networks is surely much more complicated, these cases can be considered null models against which the structure of real-world hypergraphs can be compared.

We studied various forms of higher-order contagion and healing: (i) collective contagion, in which all other members of the hyperedge need to be infected for contagion to occur, (ii) individual contagion, in which at least one member of the hyperdegree needs to be infected, and (iii) higher-order healing, in which pairwise interactions are infectious, while higher-order interactions heal. Other forms of higher-order contagion could in principle be studied with the same methodology, but we leave these studies for future research.

Now, we mention some of the limitations of our study. First, since we focused on the simplest contagion model, an important question left for future research is whether our results remain valid for more realistic epidemiological models (e.g., such as those used to model COVID-1914,40). Our model also does not apply to non-Markovian contagion dynamics, which are important when modeling real-world epidemics. From a technical standpoint, another limitation is that we used a mean-field description of the dynamics, and it is known that such a description is not quantitatively accurate for moderate values of the infected population value.28,41 Since we were mainly interested in the behavior close to the onset of epidemics, the mean-field approximation was enough for our purposes. However, more precise descriptions could be obtained as in Refs. 28 and 42. Another important limitation of our hypergraph model is that we assume that the probability that two nodes belong to the same hyperedge is a function of their hyperdegrees. While this assumption can be relaxed by considering additional nodal variables, it is possible that such a model might be inadequate to describe some real-world networks. Finally, we note that our model relies on knowledge of the functions fm, which encode the organization of hyperedges across different hyperedge sizes. These functions have not yet been estimated from real-world networks, but as progress is made toward understanding the organization of higher-order interactions,23 the determination of these functions could be a natural next step.

While in this paper we applied our hyperdegree-based mean-field equation to the SIS epidemic model, the same formalism could be applied to other dynamical processes on hypergraphs, such as synchronization, opinion formation, and other types of epidemic models. We believe that this methodology will be useful to study the effect of heterogeneity on these hypergraph dynamical processes.

Nicholas W. Landry acknowledges a helpful conversation with Danny Abrams and Lluís Arola-Fernández on the “hipster effect.”

The data that support the findings of this study are openly available in GitHub at https://github.com/nwlandry/SimplexSIS, Ref. 44.

1. Microscopic simulation of the hypergraph SIS model

We simulated the stochastic SIS model on a given hypergraph as a discrete-time Markov process on the nodes with transitions to infected and healthy states through the modalities described in Sec. II B, a variant of which was simulated in Ref. 44. We denote the binary states of the nodes at a time step t by a vector Xt=(X1t,X2t,,XNt), where Xit=0 if node i is healthy and Xit=1 if it is infected. In this model, we assume that the events that a hyperedge infects node i and that a pairwise connection infects node i are independent. Likewise, we assume that an infected neighbor, whether through a pairwise or group connection, infects a node independently of any other neighboring node. The probability that a single infected node infects its pairwise neighbor in the time interval [t,t+Δt] is β2Δt; therefore, the probability that no neighboring node infects a given node is

where A is the adjacency matrix with entries Aij=1 if nodes i and j are connected by a link and 0 otherwise.

In the collective contagion model, the probability that a triangle infects a node in the time interval [t,t+Δt] is β3Δt provided that the other two nodes are infected. Therefore, the probability of no triangles infecting node i can be written as

where the sum is over all triangles {i1,i2,i} with node i as a member. Last, the rate of healing is constant and independent of the infection status of any neighboring nodes; therefore, the probability that an infected node heals in a time interval [t,t+Δt] is γΔt.

The Markov process can then be described as

(A1)
(A2)

In our simulations, we updated the status of the nodes synchronously at times t=0,Δt,2Δt,,nΔt, where Δt=0.1.

Our specific implementation is described in what follows. We note that for all mechanisms of infection and healing described next, uiUniform(0,1) and this variable is drawn independently for each modality and each node i. At each time step, we iterate through every node and follow the following conditional logic. If a node i is already infected, it is healed if ui<γΔt and remains infected otherwise. Next, if the node i is currently healthy, it is infected by its pairwise neighbors if ui<1(1β2Δt)(AX)i and remains healthy otherwise. If node i still remains healthy after being subjected to a pairwise infection, the node is infected by its triangle neighbors if

and remains healthy otherwise. Note that each infection mechanism is only dependent on the prior time step; therefore, the order of these steps does not matter.

At t=0, the network is randomly and uniformly seeded with a small fraction (p=0.001) of infected nodes, and at each subsequent step, the current state is iterated as described above and the population average, xt=i=1NXit/N, is stored. To avoid the absorbing state Xt=0, we infect a single randomly chosen node if the population becomes completely healthy. To mitigate the effect of variability in the stochastic simulation, we average the time response of xt over a sufficient time window (determined from the average infected response curves) after it reached the steady-state. In this study, we ran the simulation for a fixed set of parameters {γ,β2,β3} for 1000 time steps and averaged over the last 300 time steps.

To find the bistability index, we initialize the simulation with a small fraction of infected agents for a fixed β3 value and incrementally increase β2 from a sufficiently small value (typically β2c/2) to a value above the critical value of β2 (typically 3β2c/2) and then incrementally decrease the value of β2 down to its original value. As described previously, if the equilibrium value while increasing the value of β2 is distinct from the equilibrium value while decreasing the value of β2 for the same β2 values, this indicates the presence of bistability. We simulated several equilibrium curves corresponding to different β3 values to observe the value of β3 at which the response curve starts to show bistability and thus infer the value of β3c.

2. Network models

We exclusively considered networks generated using the configuration model in order to isolate the effect of the degree distribution. Although the configuration model has the potential to contain both self-loops and multi-edges, in practice, the fraction of these types of edges is small,41 and in our numerical experiments, the number of self-loops was approximately 1% of the total number of nodes.

We used networks of size N=104 in the simulation of the hypergraph SIS model because this was sufficiently large enough to reduce the finite-size effects. Because the network realization was relatively large, we did not average over an ensemble of these random graphs as in Ref. 24. We have described in Sec. III A the particular distributions examined.

We generated the triangles in two different ways corresponding to the two separate cases: degree-correlated and uncorrelated. For the first case, we used the same degree sequence as used to generate the network using the configuration model and extended the configuration model to triangles as has been done in prior work.32,33 Because this is analogous to the construction of the network configuration model, there is also the possibility for self-loops and multi-edges, but this probability is low. For the independently distributed triangles, we drew with replacement a fixed number of triples (enforcing the mean triangle degree) containing node indices and assigned these nodes to a triangle. Again, as with the standard configuration model, there is the possibility for self-loops and multi-edges, but the probability of either occurring is small.

3. The numerical computation of β3c

In Sec. IV, we plotted the numerical solution of β3c for truncated power-law distributions as a function of the maximum degree and power-law exponent. In this section, we discuss the specific methodology in generating these results.

First, we describe the process for finding the bistability index accurately from the mean-field equations (10) and (15) and (16) for the correlated and uncorrelated cases, respectively. Since the V=0 solution becomes linearly unstable at β2=β2c and the stable V>0 solution is monotonically increasing with β2, the bistability index B(β3) coincides with the value of the largest root of Eq. (10) for the correlated case [or Eqs. (15) and(16) for the uncorrelated case] at β2=β2c, as shown schematically in Fig. 10. Therefore, using our analytical knowledge of β2c, we set B(β3)Vϵ, where Vϵ is the largest root at β2=β2cϵ with ϵ=105 being a small number added for numerical robustness. (We verified that this method gives numerically accurate results when compared with other methods that do not require knowledge of β2c but are more computationally intensive.)

FIG. 10.

Illustration of the bistability index with respect to the solutions to the mean-field equation in the bistable regime.

FIG. 10.

Illustration of the bistability index with respect to the solutions to the mean-field equation in the bistable regime.

Close modal

Being able to compute B(β3), we find β3c=sup{β3|B(β3)=0} by bisection: starting with an interval [β3min,0,β3max,0] such that B(β3min,0)=0, B(β3max,0)>0, we recursively define the interval [β3min,i+1,β3max,i+1] as [β3min,i,β~i] if B(β~i)>0 and [β~i,β3max,i] if B(β~i)=0, where β~i=(β3min,i+β3max,i)/2. When the length of the interval [β3min,i,β3max,i] is less than the tolerance 104, we set β3c=β3min,i.

1.
R.
Pastor-Satorras
,
C.
Castellano
,
P.
Van Mieghem
, and
A.
Vespignani
, “
Epidemic processes in complex networks
,”
Rev. Mod. Phys.
87
,
925
(
2015
).
2.
P.
Trapman
, “
On analytical approaches to epidemics on networks
,”
Theor. Popul. Biol.
71
,
160
173
(
2007
).
3.
I. Z.
Kiss
,
J. C.
Miller
,
P. L.
Simon
et al.,
Mathematics of Epidemics on Networks
(
Springer
,
Cham
,
2017
), Vol.
598
.
4.
T.
House
, “
Modelling epidemics on networks
,”
Contemp. Phys.
53
,
213
225
(
2012
).
5.
M. E.
Newman
, “
Spread of epidemic disease on networks
,”
Phys. Rev. E
66
,
016128
(
2002
).
6.
R.
Pastor-Satorras
and
A.
Vespignani
, “
Epidemic dynamics in finite size scale-free networks
,”
Phys. Rev. E
65
,
035108
(
2002
).
7.
Y.
Wang
,
D.
Chakrabarti
,
C.
Wang
, and
C.
Faloutsos
, “Epidemic spreading in real networks: An eigenvalue viewpoint,” in Proceedings of 22nd International Symposium on Reliable Distributed Systems, 2003 (IEEE, 2003), pp. 25–34.
8.
F.
Xiong
and
Y.
Liu
, “
Opinion formation on social media: An empirical approach
,”
Chaos
24
,
013130
(
2014
).
9.
D. J.
Watts
and
P. S.
Dodds
, “
Influentials, networks, and public opinion formation
,”
J. Consum. Res.
34
,
441
458
(
2007
).
10.
R.
Cowan
and
N.
Jonard
, “
Network structure and the diffusion of knowledge
,”
J. Econ. Dyn. Control
28
,
1557
1575
(
2004
).
11.
T. W.
Valente
,
Network Models of the Diffusion of Innovations
(
Hampton Press
,
1995
), 303.484 V3.
12.
M.
Gladwell
,
The Tipping Point: How Little Things Can Make a Big Difference
(
Little
,
Brown
,
2006
).
13.
D.
Centola
,
J.
Becker
,
D.
Brackbill
, and
A.
Baronchelli
, “
Experimental evidence for tipping points in social convention
,”
Science
360
,
1116
1119
(
2018
).
14.
A.
Arenas
,
W.
Cota
,
J.
Gomez-Gardenes
,
S.
Gómez
,
C.
Granell
,
J. T.
Matamalas
,
D.
Soriano-Panos
, and
B.
Steinegger
, “A mathematical model for the spatiotemporal epidemic spreading of COVID19,” medRxiv (2020).
15.
B.
Banerjee
,
P. K.
Pandey
, and
B.
Adhikari
, “A model for the spread of an epidemic from local to global: A case study of COVID-19 in India,” arXiv:2006.06404 (2020).
16.
F.
Battiston
,
G.
Cencetti
,
I.
Iacopini
,
V.
Latora
,
M.
Lucas
,
A.
Patania
,
J.-G.
Young
, and
G.
Petri
, “
Networks beyond pairwise interactions: Structure and dynamics
,”
Phys. Rep.
874
,
1
92
(
2020
).
17.
M. A.
Porter
, “Nonlinearity+ networks: A 2020 vision,” arXiv:1911.03805 (2019).
18.
L.
Horstmeyer
and
C.
Kuehn
, “
Adaptive voter model on simplicial complexes
,”
Phys. Rev. E
101
,
022305
(
2020
).
19.
P. S.
Skardal
and
A.
Arenas
, “
Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes
,”
Phys. Rev. Lett.
122
,
248301
(
2019
).
20.
C.
Xu
,
X.
Wang
, and
P. S.
Skardal
, “
Bifurcation analysis and structural stability of simplicial oscillator populations
,”
Phys. Rev. Res.
2
,
023281
(
2020
).
21.
A. P.
Millán
,
J. J.
Torres
, and
G.
Bianconi
, “
Explosive higher-order Kuramoto dynamics on simplicial complexes
,”
Phys. Rev. Lett.
124
,
218301
(
2020
).
22.
J.
Grilli
,
G.
Barabás
,
M. J.
Michalska-Smith
, and
S.
Allesina
, “
Higher-order interactions stabilize dynamics in competitive network models
,”
Nature
548
,
210
213
(
2017
).
23.
A. R.
Benson
,
R.
Abebe
,
M. T.
Schaub
,
A.
Jadbabaie
, and
J.
Kleinberg
, “
Simplicial closure and higher-order link prediction
,”
Proc. Natl. Acad. Sci. U.S.A.
115
,
E11221
E11230
(
2018
).
24.
I.
Iacopini
,
G.
Petri
,
A.
Barrat
, and
V.
Latora
, “
Simplicial models of social contagion
,”
Nat. Commun.
10
,
2485
(
2019
).
25.
G. F.
de Arruda
,
M.
Tizzani
, and
Y.
Moreno
, “Phase transitions and stability of dynamical processes on hypergraphs,” arXiv:2005.10891 (2020).
26.
F. B.
Pedro Cisneros-Velarde
, “Multi-group SIS epidemics with simplicial and higher-order interactions,” arXiv:2005.11404v1 (2020).
27.
G. F.
de Arruda
,
G.
Petri
, and
Y.
Moreno
, “
Social contagion models on hypergraphs
,”
Phys. Rev. Res.
2
,
023032
(
2020
).
28.
J. T.
Matamalas
,
S.
Gómez
, and
A.
Arenas
, “
Abrupt phase transition of epidemic spreading in simplicial complexes
,”
Phys. Rev. Res.
2
,
012049
(
2020
).
29.
B.
Jhun
,
M.
Jo
, and
B.
Kahng
, “
Simplicial SIS model in scale-free uniform hypergraph
,”
J. Stat. Mech.: Theory Exp.
2019
,
123207
).
30.
M.
Boguná
and
R.
Pastor-Satorras
, “
Epidemic spreading in correlated complex networks
,”
Phys. Rev. E
66
,
047104
(
2002
).
31.
J. C.
Miller
,
A. C.
Slim
, and
E. M.
Volz
, “
Edge-based compartmental modelling for infectious disease spread
,”
J. R. Soc. Interface
9
,
890
906
(
2012
).
32.
O. T.
Courtney
and
G.
Bianconi
, “
Generalized network structures: The configuration model and the canonical ensemble of simplicial complexes
,”
Phys. Rev. E
93
,
062311
(
2016
).
33.
J.-G.
Young
,
G.
Petri
,
F.
Vaccarino
, and
A.
Patania
, “
Construction of and efficient sampling from the simplicial configuration model
,”
Phys. Rev. E
96
,
032312
(
2017
).
34.
A. J.
Ward
,
D. J.
Sumpter
,
I. D.
Couzin
,
P. J.
Hart
, and
J.
Krause
, “
Quorum decision-making facilitates information transfer in fish shoals
,”
Proc. Natl. Acad. Sci. U.S.A.
105
,
6948
6953
(
2008
).
35.
S. C.
Pratt
,
E. B.
Mallon
,
D. J.
Sumpter
, and
N. R.
Franks
, “
Quorum sensing, recruitment, and collective decision-making during colony emigration by the ant Leptothorax albipennis
,”
Behav. Ecol. Sociobiol.
52
,
117
127
(
2002
).
36.
B.
Guerra
and
J.
Gómez-Gardenes
, “
Annealed and mean-field formulations of disease dynamics on static and adaptive networks
,”
Phys. Rev. E
82
,
035101
(
2010
).
37.
J. D.
Touboul
, “
The hipster effect: When anti-conformists all look the same
,”
Discrete Continuous Dyn. Syst. B
24
,
4379
4415
(
2019
).
38.
G.
St-Onge
,
V.
Thibeault
,
A.
Allard
,
L. J.
Dubé
, and
L.
Hébert-Dufresne
, “School closures, event cancellations, and the mesoscopic localization of epidemics in networks with higher-order structure,” arXiv:2003.05924 (2020).
39.
B. M.
Althouse
,
E. A.
Wenger
,
J. C.
Miller
,
S. V.
Scarpino
,
A.
Allard
,
L.
Hébert-Dufresne
, and
H.
Hu
, “Stochasticity and heterogeneity in the transmission dynamics of SARS-COV-2,” arXiv:2005.13689 (2020).
40.
D.
Balcan
,
B.
Gonçalves
,
H.
Hu
,
J. J.
Ramasco
,
V.
Colizza
, and
A.
Vespignani
, “
Modeling the spatial spread of infectious diseases: The global epidemic and mobility computational model
,”
J. Comput. Sci.
1
,
132
145
(
2010
).
41.
M.
Newman
,
Networks
(
Oxford University Press
,
2018
).
42.
R.
Albert
and
A.-L.
Barabási
, “
Statistical mechanics of complex networks
,”
Rev. Mod. Phys.
74
,
47
(
2002
).
43.
Á.
Bodó
,
G. Y.
Katona
, and
P. L.
Simon
, “
SIS epidemic propagation on hypergraphs
,”
Bull. Math. Biol.
78
,
713
735
(
2016
).
44.
N.
Landry
, (2020). “SimplexSIS,”
Zenodo
 V1.0.