The effect of heterogeneity on hypergraph contagion models

The study of contagion processes is a fundamental problem in network science, with applications including epidemics,^1–7 social media,⁸ opinion formation,⁹ 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.¹⁶ It has been shown that higher-order interactions in networks (i.e., interactions involving multiple nodes) can have profound effects on dynamical network processes¹⁷ such as opinion formation,¹⁸ synchronization,^19–21 and population dynamics.²² Efforts to map higher-order interactions in real-world networks have an uncovered rich structure,²³ which is only now starting to be appreciated. In the context of contagion processes, it was recently shown²⁴ 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.²⁹ 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.

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

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′/(N⟨k⟩)$⁠, 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.²⁴ 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.¹ 

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 $t≥0$⁠, 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.

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

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 $m−1$ nodes participating in the hyperedge (⁠$k,k1,…,km−1$⁠), counting how many such combinations there are not counting permutations [$P(k1),…,P(km−1)/(m−1)!$], calculating what fraction of such combinations forms a hyperedge with the node in consideration [$fm(xk,xk1,…,xkm−1)$], multiplying by the probability that the hyperedge can transmit the infection [$G(xk1,…,xkm−1)$], 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,…,km−1)=fm(k,k1,…,km−1)$⁠. With these assumptions and using the collective contagion rule in Eqs. (3) and (2) becomes

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,

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

For the link connection probability $f2(k,k1)$⁠, we will take $f2(k,k1)=kk1/(N⟨k⟩)$⁠, 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/(N⟨k⟩)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)=2⟨k⟩/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.

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/(N⟨q⟩)2$⁠, which corresponds to three nodes being connected at random according to the configuration model for triangles.³² For the pairwise links, we define the degree-correlated and uncorrelated cases as before, where in the degree-correlated case, $f2(q,q1)=qq1/(N⟨q⟩)$ 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/(N⟨k⟩)2$⁠, Eq. (6) can be rewritten in terms of the fraction of infected links

as

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

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

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)/(N⟨k⟩)$⁠, and summing over $k$⁠, which yields the linearized equation for the evolution of the perturbation $δV$

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)=2⟨k⟩/N2$⁠. In this case, Eq. (6) can be rewritten in terms of the fraction of infected nodes,

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

As in the prior case, the equilibrium is

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

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

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)∝k−4$ 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.

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)∝k−4$ 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$⁠.

To quantify how the onset of bistability depends on the hypergraph parameters, we define the bistability index $B(β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)∝k−4$ for $67<k<1000$ and $0$ otherwise, and (c) $P(k)∝k−3$ 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.

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

For the correlated case, $f3(k,k1,k2)=2kk1k2/(N⟨k⟩)2$⁠, this can be rewritten as

with a fixed point

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

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

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.

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

with equilibrium

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

Linearizing, we obtain the system

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,

This relationship implies that there is a singularity when $β3=β3∗=γ⟨k2⟩/[2(⟨k2⟩−⟨k⟩2)⟨k⟩]$⁠. However, one can check that $β2$ is negative at $β3=β3∗$⁠, and therefore, the singularity is not physically relevant. Note that when $⟨k2⟩=⟨k⟩2$ 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.³⁷ 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,β3≥0$⁠.

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,β3≥0$⁠, 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.

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

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

for individual contagion and

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.²⁴ We focus, however, on the case of individual contagion. Linearizing, we find that the solution $xk=0$ becomes unstable when

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)∝k−r$ if $50≤k≤kmax$ 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.

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

First, note that $h(0,β2)=β2/β2c−1$ and that $h(1,β2)<0$⁠. Therefore, if $∂h∂V(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

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

which, when evaluated in

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

where

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 $a12−4a0a2=0$ must satisfy the inequality $0≤−a1/2a2≤1$⁠. In addition, we note that because of continuity, the sign of the $a2$ term must be negative because otherwise, $∂h∂V(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⟩/⟨k⟩4$⁠. Using these conditions, we can construct a piecewise definition of $β3c$

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.

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-19^14,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,²³ 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.

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

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, $ui∼Uniform(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$⁠.

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,⁴¹ 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.

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 $ϵ=10−5$ 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.)