Influential groups for seeding and sustaining nonlinear contagion in heterogeneous hypergraphs

Contagion phenomena are often the results of multibody interactions—such as superspreading events or social reinforcement—describable as hypergraphs. We develop an approximate master equation framework to study contagions on hypergraphs with a heterogeneous structure in terms of group size (hyperedge cardinality) and of node membership (hyperdegree). By mapping multibody interactions to nonlinear infection rates, we demonstrate the influence of large groups in two ways. First, we characterize the phase transition, which can be continuous or discontinuous with a bistable regime. Our analytical expressions for the critical and tricritical points highlight the influence of the first three moments of the membership distribution. We also show that heterogeneous group sizes and nonlinear contagion promote a mesoscopic localization regime where contagion is sustained by the largest groups, thereby inhibiting bistability. Second, we formulate an optimal seeding problem for hypergraph contagion and compare two strategies: allocating seeds according to node or group properties. We find that, when the contagion is sufficiently nonlinear, groups are more effective seeds than individual hubs. Group interactions can dramatically alter social contagion dynamics and lead to the emergence of new phenomena like abrupt transitions and critical mass effects. The authors develop an approximate master equation framework to analytically describe contagion in heterogeneous hypergraphs and study the impact of large influential groups in seeding and sustaining epidemics.


I. INTRODUCTION
Mathematical models of contagion processes enhance our understanding of spreading dynamics and our predictive capabilities [1].To account for the interconnected nature of real-world systems, the last two decades of network science and computational epidemiology research have focused on modeling frameworks of increasing complexity [1][2][3].From the spreading of infectious diseases to rumors and innovations [4][5][6], a crucial as-pect remains the interplay between the contact patterns enabling transmission and the dynamics that unfolds on top.As a representation for these contacts, graphs have been widely exploited to better represent real-world patterns with increasing levels of accuracy [1,2,7].
While graphs remain a reference for the representation of complex systems, they come with a fundamental limitation: they can only encode pairwise interactions (represented by edges).Groups are instead the foundational units of many biological, ecological and social systems, whose processes can involve multi-body interactions between any number of elements.To overcome this limitation, higher-order network representations [8] are becoming more and more popular [9][10][11] to describe the structure of interacting systems [12,13], their growth [14][15][16][17] and dynamics in groups [18][19][20].Simplicial complexes and hypergraphs can encode relationships and interactions between any number of elements.They have been used to investigate the implications of higher-order interactions for landmark dynamical processes like synchronization [21][22][23], diffusion [24][25][26], and other social dynamics [27][28][29].Ref. [11] provides a comprehensive review of early efforts in this direction.
Recently, a model of simplicial contagion has been proposed [30].This is a standard Susceptible-Infected-Susceptible (SIS) compartmental model in which a susceptible individual can become infected through different transmission channels, beyond infectious edges.In models of simple contagion, the transition from susceptible to infected happens independently with each exposure to an infectious edge [1].In models of complex contagion instead [31], the transition requires multiple infectious edges or is reinforced by multiple exposures, thus accounting for the empirically observed mechanisms of social influence [32][33][34][35].In simplicial contagions, or more generally hypergraph contagions, a susceptible individual can become infected because of a multi-body process, i.e., through an exposure to an infectious group [30].In this way, the node is simultaneously exposed to the state of the entire group, whose effect can be interpreted as a mechanism of social reinforcement [32].In addition, the study of higher-order contagion models has applications as well in biological sciences: it has indeed recently been shown that the combination of multi-body interactions, heterogeneous temporal activity, and the concept of a minimal infective dose lead to nonlinear infection kernels in a model of biological contagions [36].
The analytical approaches derived so far have confirmed the rich phenomenology emerging from the contagion dynamics, characterized by discontinuous phase transitions, bistability and critical mass effects [30,[37][38][39][40][41].Most approaches follow a heterogeneous mean-field (HMF) framework in which nodes are divided into hyperdegree classes.These descriptions are analytically tractable, but do not consider the details of the structure and ignore the dynamical correlations within groups, which are especially important for hypergraph contagions since multi-body interactions naturally reinforce these correlations.
Other approaches like quenched mean-field theory (QMF) [42] and microscopic Markov-chains (MMCA) [43] can explicitly take the entire contact patterns into account.Along this line, the microscopic epidemic clique equations (MECLE) capture dynamical correlations as well, thereby highlighting the important impact these correlations have on critical points [44].The analytical tractability of these approaches is, however, sacrificed in favor of a more precise description of the structure.To fully understand the consequences of multi-body interactions in contagions on higherorder networks, we thus need a framework that is both analytically tractable and captures dynamical correlations.
In particular, we are interested in better understanding the notion of influence in hypergraph contagions.In classic contagion models on random networks, individual hubs are influential in the sense that they are the best seeds of contagions, but they are also the most apt at sustaining seeded contagions [45].However, in hypergraph contagions we must consider the influence of both individuals and groups, because dynamical correlations can allow groups to be more influential than sets of uncorrelated hubs.Regimes of bistability and hysteresis also imply a potential decoupling between the ability of nodes to seed a contagion and their ability to sustain it.We thus wish to determine (i) which set of groups can best maintain the stationary state of a hypergraph contagion, and when this becomes a dominant effect (ii) which set of groups and their configuration offer optimal initial conditions for a contagion, as compared to the classic notion of influential spreaders as individual hubs, and (iii) whether or not these two notions of group influence are aligned.
We use approximate master equations (AMEs) [19,20,[46][47][48][49] to study hypergraph contagions, capturing exactly the inner dynamics of groups.We apply this framework to a contagion model where the infection rate is a nonlinear function of the number of infected nodes in groups, and show that influential groups can drive both the stationary state of the contagion and its behavior in the transient state.
The paper is structured as follows.In Sec.II A, we map the simplicial contagion model onto a hypergraph contagion process and introduce the mathematical formalism of group-based AMEs [19,48,49].To illustrate the accuracy of our framework and assess its limitations, we compare the predictions of our AMEs with simulations on real-world hypergraphs and their randomized counterparts.
In Sec.II B, we characterize the importance of influential groups on the phase transition.We derive analytical expressions for the critical and tricritical points, and study the impact of a heterogeneous structure.We show that large values for the third moment of the membership (hyperdegree) distribution suppress the emergence of a discontinuous phase transition, a result consistent with HMF theories [38,39].Our formalism allows us to go further and capture the emergence of a mesoscopic localization phase, where infected nodes are concentrated in the largest groups [48,49].We show that localization effects can drive the onset of an endemic phase and, incidentally, inhibit bistability.
Finally, in Sec.II C, we showcase the importance of influential groups on the transient state by considering a problem of influence maximization.We define and solve the optimization problem based on two strategies, allocating seeds according to either node individual properties or according to group properties.When the contagion is sufficiently nonlinear, we find that the latter is FIG. 1. Mapping of the simplicial contagion model to a hypergraph contagion.We use a bipartite representation, where nodes (white circles) belong to groups (black circles).A facet of dimension n − 1 is mapped onto a group of size n.In the simplicial contagion model, contributions from higher-order interactions are taken into account by additional transmission rates (e.g., β2 for the 2-simplex) when all but one node of the simplex are infected [30].With our hypergraph representation, infections within a group are simply modeled by a general infection function β(n, i) that depends on both the size n of the group and the number i of infected nodes in the group (with i ≤ n).more effective.This suggests again that contagion processes driven by groups lead to a fundamentally different phenomenology than pairwise contagions, such that investigations and modeling of social contagions should carefully account for higher-order interactions.

A. Hypergraph contagion model
In the original version of the simplicial contagion model [30], the spreading process takes place on a simplicial complex [see Fig. 1] and obeys the following rule.If all nodes in a d-simplex are infected except a susceptible one, this remaining node gets infected at a rate β d , and also receives contributions from the lower-dimensional simplices included in the d-simplex with rates respectively equal to β d−1 , . . ., β 1 .For instance, in Fig. 1, two of the three nodes composing a 2-simplex are infected, hence the susceptible node gets infected at rate 2β 1 + β 2 , considering the contributions both from the two edges and from the "triangle".
A simplicial complex is a specific type of hypergraph, and thus we can always map the former on the latternote that the reverse direction is not always possible.To do so, each facet-a simplex that is not a face of any larger simplex-is represented by a single hyperedge (group).In this paper, we relax the requirement of the simplicial complex in favor of a more general hypergraph structure.We make use of the bipartite representation of hypergraphs [8,11], in which the two sets of nodes of the bipartite graph correspond respectively to the sets of nodes and groups of the original hypergraph, as illustrated in Fig. 1.
The size n of a group corresponds to the number of nodes belonging to this group, and it is therefore equivalent to the hyperedge cardinality.Note that a d-simplex consists of d + 1 nodes and is therefore mapped onto a group of size n = d + 1.Similarly, the membership m of a node, the node hyperdegree, corresponds to the number of groups to which it belongs, regardless of their size.
The hypergraph contagion model is defined as follows: for a group of size n, where i ≤ n members are infected, each of the n − i susceptible nodes gets infected at rate β(n, i).For susceptible nodes that belong to multiple groups, their total transition rate to the infected state is simply the sum of the infection rates associated with each group to which they belong-in other words, the infection processes are independent.All infected nodes transition back to the susceptible state at the same constant recovery rate µ.
Notice that the hypergraph contagion model allows to represent any type of simplicial contagion.For instance, in the simplicial contagion model case of Fig. 1, we would use β(3, 2) = 2β 1 + β 2 .In fact, the description offered by the infection rate function β(n, i) yields a variety of models more general than the original simplicial contagionin which a function β(i) would be sufficient to encode contributions from facets of any dimension.
In all our case studies, we will use an infection rate function of the form However, many results we derive hold true for a general infection rate function β(n, i).The parameter ν controls the nonlinearity of the contagion.A linear contagion is recovered by setting ν = 1, which is equivalent to a standard SIS model on networks, where each group is a clique [19].We intentionally chose the infection rate function independent of n to focus on the impact of a nonlinear dependence on i; it would be straightforward to generalize the results by considering The infection rate function in Eq. ( 1) is the simplest nonlinear generalization of standard epidemiological models, where β(n, i) ∝ i.Moreover, we can motivate the choice of exponents ν = 1 in the context of social contagions, by comparing our approach to the original formulation of the simplicial contagion model.A value of β 2 > 0 in Fig. 1 represents social reinforcement [30], and to correctly map the infection rate for a triangle, we need to use an exponent ν > 1 in our model.Similarly, a value β 2 < 0 represents social inhibition, and this case can be obtained with an exponent ν < 1.
Another motivation for the infection rate function at Eq. ( 1) is a recent study that shows this general form emerges in the occurrence of heterogeneous temporal patterns [36].More specifically, if you consider that the participation time of nodes-representing individualsto higher-order interactions is distributed according to a power law, and that individuals become infected according to a threshold mechanism based on the dose received in the interaction, then the probability for a node to get infected in a group is ∝ i ν , where ν is related to the temporal heterogeneity.In the continuous time limit, one recovers the infection rate function defined at Eq. ( 1).
In Ref. [36], the infection mechanism is motivated in the context of biological contagions, where the infective dose received could represent viral particles for instance, and the threshold would correspond to the minimal infective dose to develop a disease.While such type of complex contagions are rarely used in the context of biological contagions, they could help explain certain observed phenomena, such as super-exponential spread for certain diseases [50].Moreover, threshold models are very common in social contagions [32,[51][52][53], thus Eq. ( 1) could be interpreted as an effective mechanism of social spread accounting for heterogeneous temporal patterns.

Group-based Approximate Master Equations
To describe hypergraph contagions, we make use of group-based Approximate Master Equations [19,20,48,49].This means that we do not rely on specific hypergraph realizations.Instead, we assume that the structure is drawn from a random hypergraph ensemble described by the distributions p n , for the size n of a group, and g m , for the membership m of a node.Each of the m membership stubs of a node is assigned uniformly at random to a group available spot.Therefore, the membership m of a node and the sizes of the groups to which it belongs are uncorrelated.
To track the evolution of a contagion process on this ensemble of hypergraphs, we define two sets of quantities: s m (t), representing the fraction of nodes with membership m that are susceptible at time t and f n,i (t), the fraction of groups of size n having i infected members at time t.The last quantity can also be interpreted as a conditional probability (of observing i infected nodes in a group of size n) satisfying the normalization condition i f n,i = 1.We further define two mean-field quantities.First, let us take a random susceptible node.The mean-field infection rate resulting from a random group to which it belongs is defined as Indeed, the joint distribution for the size n and the number of infected nodes i in this group is proportional to (n − i)f n,i p n , and we just average β(n, i) over this distribution.
Second, let us randomly choose a susceptible node inside a group.The mean-field infection rate caused by all the external groups to which the susceptible node belongs (excluding the one from which we picked the node) can be written as To obtain ρ(t), we assume that infections coming from different groups are independent processes.We multiply r(t) with the mean excess membership of a susceptible node, i.e., if we pick a susceptible node in a group, it is the expected number of other groups to which it belongs.Since the membership distribution of a susceptible node picked in a group is proportional to ms m (t)g m , we simply average m − 1, its excess membership, over this distribution.
Using the definitions in Eqs. ( 2) and (3) we can write the following system of AMEs This system is composed of O(n 2 max +m max ) equations.It is approximate in the sense that the evolution of the fractions of infected nodes of membership m, s m , are treated in a mean-field fashion (still considering dynamic correlations between pairs node-group), while the evolution equations of the f n,i are treated as master equations.In the right-hand side of Eq. (4b), the first two terms are due to the recovery process, and the last two to the infection.The infection rate due to infected nodes inside the group is exact, while the infection rates associated to external groups (the terms making use of ρ) are treated in an approximate way.Without loss of generality (up to a change of time scale) we set µ ≡ 1 for the remainder of the document.
From Eq. ( 4), we can calculate the global prevalence and the group prevalence which correspond to the average fraction of infected nodes in the whole system and within groups of size n respectively.
In the stationary state, we obtain the following self-consistent relations: The latter equation can be simplified by noting that f * n,i must respect detailed balance in the stationary state, i.e., expressing that the probability to decrease the number of infectious nodes in a group of size n from i + 1 to i by a recovery process is equal to the probability of the opposite change (from i to i+1 infectious nodes) obtained through a contagion event.We thus finally obtain where f n,0 = 1 − i>0 f n,i .

Comparison with simulations
We provide a comparison of our approach with Monte Carlo simulations (see Methods Sec.IV A 1).To do this, we consider empirical higher-order structures constructed from real-world data, and their randomized counterparts.Details on the data collection and aggregation are provided in Methods Sec.IV A 2.
The motivation for this a priori validation is threefold: first, it allows us to illustrate the validity and accuracy of our analytical framework when our assumptions are met-namely when the structure is drawn from an ensemble of uncorrelated hypergraphs with fixed g m and p n .Second, because of the excellent agreement with simulations for random hypergraphs, we omit further comparison with Monte Carlo simulations for the many results we present in the following sections.Finally, it showcases the possible sources of discrepancies-and how our results could vary-when considering real datasets.
In Fig. 2, we show the phase diagram, i.e., the fraction I * of infectious nodes in the stationary state, of contagion dynamics on hypergraphs that encode higher-order social (face-to-face) interactions between individuals from a primary school in Lyon (see Methods Sec.IV A 2 for details).Both the membership and group size distributions are homogeneous for this dataset.We considered linear contagion (ν = 1), equivalent to the standard SIS model, and superlinear contagion (ν = 4).In both cases, our analytical formalism (continuous lines) agrees quite well with the Monte Carlo simulations (symbols) on the original (empirical) hypergraph [Figs.2(a) and 2(d)].The main source of errors for the observed discrepancy can be attributed to structural correlations, which do not appear to affect the threshold values, but reduce the stationary prevalence.Indeed, in Figs.2(b) and 2(e), the agreement improves by randomizing the hypergraph while preserving the memberships and the group sizes.The remaining discrepancies are due to the fact that simulations are affected by finite-size effects, while our formalism assumes an infinite size system: the agreement becomes almost perfect in Figs.2(c) and 2(f) by additionally increasing the size of the hypergraph by a factor 10.
Let us remark that the error is more important for ν = 4 on the original hypergraph [Fig.2(d)], which suggests that structural correlations have a greater effect on nonlinear contagions.In the Supplementary Information, we show how to generalize our framework to account for structural correlations.
We also considered a completely different dataset, which represents co-authorship relations in computer science publications, obtained from major journals and proceedings in the field.The resulting hypergraph is considerably larger than the previous one, and it also presents a very different structure (see Methods Sec.IV A 2).The results are shown in Fig. 3, where we plot the same phase diagram curves as in Fig. 2, using a superlinear contagion (ν = 2).In this case, however, the membership distribution is heterogeneous [Fig.3(a)], approximately of the form g m ∼ m −γm with γ m ≈ 2.3, while the group size distribution is more homogeneous [inset of Fig. 3(a)], but still extends to rather large values with n max = 25.By comparing the results for the original hypergraph [Fig.3(b)] against those for a randomized ensemble [Fig.3(c)], we see that structural correlations account for the major part of the discrepancies between simulations and theory, affecting both the invasion threshold and the stationary prevalence.However, Fig. 3(c) shows that structural correlations are not the only source of errors at high prevalence.
Let us recall that our formalism correctly captures the dynamical correlations within a group through a master equation description [Eq.(4b)] of the possible states, but it does not capture the dynamical correlations around nodes, since we use a heterogeneous mean-field description [Eq.(4a)].These correlations become especially important in the presence of hubs with a large membership, which is the case for the hypergraph considered in Fig. 3.In fact, when we use the same group size distribution as in Fig. 3, but with a more homogeneous membership distribution, the discrepancy at high prevalence disappears (see the Supplementary Information).

B. Phase transition
In this section, we unveil the important role of influential groups on the phase transition of hypergraph contagions.We first derive general expression for critical We compare the results of Monte Carlo simulations (symbols; see Methods Sec.IV A 1) with the predictions of our approach (solid and dashed lines for stable and unstable solutions respectively).We rescale λ with the invasion threshold λc, which is computed using Eq. ( 15) in Sec.II B. The symbols represent the average infected fraction measured over long runs and averaged over randomized hypergraphs when this applies.The error bars (sometimes too small to be seen) correspond to one standard deviation.The green circles are obtained with the quasistationary-state method, starting with a large fraction of infected nodes (I = 0.8) to sample the upper branch of the hysteresis loop when the phase transition is discontinuous.The blue triangles are obtained by ordinary simulations, starting with a small fraction of infected nodes (I = 0.02) to sample the lower branch of the hysteresis loop.We only do it for the expanded hypergraphs because finite-size effects makes unreliable the estimation of the lower branch with data of small size.
points, marking the limit of the domain of validity of a stationary solution.Secondly, we obtain an expression for tricritical points, indicating when the phase transition switches from continuous to discontinuous, and a bistable regime appears.These results are valid for any infection rate function such that β(n, i) > 0 for all i > 0 (see Supplementary Information for a consideration of threshold models).
We then concentrate our study on the infection rate function β(n, i) = λi ν .This allows us to define important threshold values, the invasion threshold λ c above which the disease-free solution I * = 0 is unstable, and the bistability threshold ν c , the smallest nonlinear exponent allowing for a discontinuous phase transition.In particular, we will show how heterogeneous membership and group size distributions alter these thresholds, especially in the presence of a mesoscopic localization driven by influential groups.
In the following, we assume that the system has reached the stationary state and we drop the asterisk to simplify the notation throughout this section.

Critical points and the invasion threshold
Equations (7a) and ( 8) imply that each s m and f n,i can be written in terms of r and ρ.In turn, r and ρ can be written in terms of s m and f n,i through Eqs. ( 2) and (3), which means the stationary solutions are determined by a set of self-consistent equations.
Satisfying all self-consistent equations means we can reexpress all quantities (s m , f n,i , ρ, r) as functions of a single mean-field quantity, which we choose to be r.s m (r) is given by Eq. (7a), and ρ(r) is given by Eq. ( 3), which is rewritten as f n,i is more simply written as a composite function, f n,i [ρ(r)], given by Eq. ( 8).Finally, r itself must sat- ) with the predictions of our approach (solid lines).We rescale λ with the invasion threshold λc, which is computed using Eq. ( 15) in Sec.II B. The infected fraction I * has been obtained as averages over time with long runs (and averaging over randomized hypergraphs when this applies).The error bars (too small to be seen) correspond to one standard deviation.We used ordinary simulations, starting with I = 0. isfy Eq. ( 2), which we can write as r = M[ρ(r)] where In Fig. 4, the stationary solutions thus correspond to the intersections of M[ρ(r)] (solid lines) and r (dashed line).We see that r = 0 is always a solution, while the solution r > 0 only exists for certain values of the parameter λ.This indicates the presence of a critical point.
Critical points mark the limit of the domain of validity of a solution to the equation r = M[ρ(r)].They arise when r is tangent to M[ρ(r)], which implies as can be seen in Fig. 4(a), where M[ρ(r)] is tangent to r at the point r = 0 for some value λ = λ c .When the tangent point r > 0, Eq. ( 11) needs to be solved numerically.However when the solution r → 0, we are able to obtain an analytical expression.In this limit, we expect to have s m → 1, ρ → 0 and f n,i → δ i,0 ∀n, i.e., all nodes are susceptible.Therefore, from Eq. ( 9) we have where • • • stands for the expectation value with respect either to g m or p n .From Eq. (10) (and using the fact that β(n, 0) ≡ 0) we also obtain To evaluate h n,i , we apply the derivative d/dρ to Eq. ( 8).We obtain h n,1 = n and Combining Eqs ( 12), ( 13) and ( 14) with Eq. ( 11), we obtain the following implicit expression for critical points in the limit r → 0 For the infection rate function β(n, i) = λi ν , Eq. ( 15) defines the invasion threshold λ c , i.e., the critical value of λ marking the limit of the validity for a solution r → 0. This solution is not always stable, but we always have that the trivial solution r = 0, corresponding to I = 0, is unstable for all λ > λ c .This is illustrated in Fig. 5(a on both the structure (g m , p n ) and the nonlinear exponent ν.The resulting phase diagram in the (λ, ν) plane is shown in Fig. 5(c).The invasion threshold spans the critical line (solid line) λ = λ c in the phase diagram.The dashed critical line λ = λ p is associated with the limit of validity of a solution where M[ρ(r)] is tangent to r for some r > 0-it is thus solved numerically.We call λ p the persistence threshold as it is the smallest value of λ such that a non-trivial solution is locally stable.For continuous phase transitions, λ c = λ p , but for discontinuous phase transitions, λ p < λ c .

Tricritical points and the bistability threshold
Depending on the structure and the form of β(n, i), we can have a continuous or a discontinuous phase transition, as can be seen in Fig. 5.When the phase transition is continuous, we have possibly two solutions for the stationary fraction of infected nodes, I 1 = 0 and I 2 > 0.
When I 2 exists (for instance when λ > λ c ), I 1 is unstable.
When the phase transition is discontinuous, we have typically three solutions, I 1 = 0 and 0 < I 2 < I 3 .In the bistable regime [for instance when λ ∈ (λ p , λ c )], all three solutions coexist, I 2 is unstable, and I 1 and I 3 are locally stable.In the endemic regime [for instance when λ ≥ λ c ], only I 1 and I 3 exist, and only I 3 is locally stable.
We are interested in knowing when the phase transition changes from continuous to discontinuous.In Fig. 5(c), the bistable regime starts at a tricritical point (star marker), where two critical lines meet.For the infection rate function β(n, i) = λi ν and a fixed hypergraph structure, the tricritical point happens at (λ, ν) = (λ c , ν c ), where ν c is what we call the bistability threshold, since a bistable regime only exists for ν > ν c .
To get some insights on the properties of tricritical points, we show in Fig. 4(b) the function M[ρ(r)] at λ = λ c and for values of ν below, at and above the bistability threshold.For ν < ν c , we have ∂ 2 r M| r→0 < 0 and I 3 does not exist.For ν > ν c , we have ∂ 2 r M| r→0 > 0 and there exists a solution I 3 > 0. At the tricritical point, Since a tricritical point is also a critical point, from Eq. ( 11) dM/dρ = (dρ/dr) −1 , so the condition can be rewritten as The derivatives on ρ with respect to r at a critical point where r → 0 can be easily evaluated, and the condition now becomes To evaluate the last term of Eq. ( 16), let us rewrite In the limit ρ → 0, f n,i → δ i,0 , which implies u(ρ) → 0 and v(ρ) → n , therefore First order derivatives can be evaluated using For the second-order derivative, let us define l n,i ≡ d 2 f n,i /dρ 2 | ρ→0 , so that we can write Finally, we can apply the second order derivative to Eq. (7b) to obtain the recurrence relation Again, l n,0 = − i>0 l n,i by definition.
Even though it is possible to express the l n,i in an explicit form, the expression does not give us more intuition, and it is simpler to calculate the l n,i using the recurrence equation just given.After rewriting d 2 M/dρ 2 | ρ→0 ≡ F [p n , β(n, i)], tricritical points are obtained by solving the equation Tricritical points result from an intricate relation between the structure (g m , p n ) and the infection rate β(n, i).Fig. 5 shows that either changing the structure [Fig.5(a)] or the shape of the infection rate function [Fig.5(b)-(c)] can lead to a change of behavior, from a continuous phase transition to a discontinuous one with a bistable regime.
A first hypothesis we can make from these simple examples is that more nonlinear infection rates (larger ν) and larger groups promote bistability.However, we will see that this intuition does not hold in general for heterogeneous structures due to the onset of mesoscopic localization.

Heterogeneous memberships
In this section we investigate the effects of a heterogeneous membership distribution g m while keeping p n = δ n,n0 homogeneous to disentangle the impact of the different types of heterogeneity.A first remark we can make about the invasion threshold [Eq.(15)] is that it is coherent with heterogeneous pair-approximation frameworks [54] on random networks when only dyadic interactions are considered, i.e., when p n = δ n,2 .In this case, we can set ν = 1 without loss of generality, thus recovering the standard SIS model.The associated threshold is , where g m can now be interpreted as the standard degree distribution of graphs.This threshold, although quite accurate for most structures, does not capture the hub reinfection mechanism [55], and thus could be inaccurate for graphs with hubs of a very large degree.More generally, for group interactions (p n = δ n,2 ) we can see that a larger average excess membership m(m − 1) / m always leads to a smaller invasion threshold λ c , akin to the standard SIS model, but the relationship is now nonlinear.To see this, let us rewrite Eq. ( 15) as Since β(n, i) is a monotonically increasing function of λ for all n, i, then the left-hand side of Eq. ( 22) is a monotonically increasing function (of λ) as well.Consequently, if the right-hand side decreases, λ c must decrease as well.7) for the corresponding membership distribution.Solid and dashed lines represent stable and unstable solutions, respectively.The infection rate function considered is β(n, i) = λi ν with ν = 2.8.The bistability threshold shows a maximum near γm = 3, suggesting that this is the maximum in the limit mmax → ∞.
Assessing the impact of membership heterogeneity on the bistability threshold ν c is more complicated.In fact, Eq. ( 21) explicitly depends on the first three moments of g m , but it also depends on the first two moments implicitly through λ c , at which F must be evaluated.
In order to build our intuition, let us assume that we are able to keep fixed the first two moments m and m 2 while increasing m 3 .This means that λ c would not change, hence the only dependence on g m would be explicit in Eq. (21).Since the term depending on m 3 is negative, increasing the third moment implies that F must increase if we want to balance Eq.(21).But since d 2 M/dr 2 | r→0 increases with ν [see Fig. 4], and thereby F as well, we can conclude that increasing m 3 while keeping the first two moments fixed leads to an increase of the bistability threshold ν c .This is validated in Fig. 6, where we consider two distributions g m sharing the same first two moments, but a different third moment [see Fig. 6(a)].By comparing the associated epidemic curves in Fig. 6(b), it is evident how the larger third moment suppresses the emergence of a bistable regime.
A corollary of this argument is that for certain structures, it is impossible to have bistability.To see this, let us consider a power-law membership distribution of the form g m ∼ m −γm .In this case, since the bistability threshold depends on the third moment of g m , while the invasion threshold only depends on the first two, by setting the exponent 3 < γ m < 4, the invasion threshold converges to a value λ c > 0, but ν c does not exist.In other words, it is impossible to have a discontinuous phase transition.
This second statement is validated in 7(a), where we show that that the bistability threshold ν c appears to grow without bound as m max → ∞ for γ m ≤ 4. Instead, for γ m ≥ 4, the bistability threshold appears to converge, as expected, since the first three moments of g m converge as well.What is more surprising is the non-monotonic behavior of ν c with respect to γ m , which we present in Fig. 7(b).The bistability threshold has a well defined maximum at a value of γ m that appears to converge to γ m = 3 for m max → ∞.In other words, γ m = 3 is the optimal value of membership exponent in suppressing the emergence of a discontinuous phase transition and the related bistability.This can be understood from Eq. ( 21): for γ m > 3, the invasion threshold does not vary much since the first two moments of g m are finite.Hence maximizing the third moment maximizes ν c , which corresponds to γ m → 3.One could still be surprised that the bistability threshold grows more slowly with m max in the range 2 < γ m < 3, since the invasion threshold λ c tends toward zero.In this case, the bistable regime exists, but its width (λ p , λ c ) simply vanishes as λ c → 0.

Heterogeneous group sizes
Let us now consider hypergraphs with heterogeneous group size distribution p n , and homogeneous membership distribution, namely g m = δ m,m0 .In this case, the invasion threshold, as defined by Eq. ( 15), depends on the whole distribution p n , which makes drawing general conclusions on the impact of a heterogeneous distribution p n more difficult.
To get some intuitions, let us consider the standard SIS model, i.e., the case ν = 1 in Eq. ( 1).With our AMEs, it was shown that [49] for power-law distributions p n ∼ n −γn with large cut-offs n max .The first term on the right-hand side of Eq. ( 23) suggests that more heterogeneous groups-size distributions p n (smaller values of γ n ) lead to smaller invasion thresholds.Intuition tells us that we should expect this behavior for ν = 1 as well.We have therefore investigated numerically in Fig. 8(a) the invasion threshold as a function of the group size exponent for different values of ν, confirming that more heterogeneous group sizes (smaller γ n ) do lead to a smaller invasion threshold, even for nonlinear infection functions (ν = 1).However, this effect is mitigated when larger values of ν are considered.For large ν and large n max , the value of the invasion threshold is dominated by the cut-off, and scales as λ c ∼ n −ν max , as illustrated in Fig. 8(b).
This behavior can be attributed to the onset of mesoscopic localization.In Refs.[48,49], it was shown analytically for ν = 1 that, for certain combinations of (γ m , γ n ), the epidemic near the invasion threshold [λ = λ c (1 + ) with 1] is dominated by the largest most influential groups.In these scenarios, the second term on the righthand side in Eq. ( 23) dominates the first one and, near λ c , the group prevalence I n grows exponentially with n, i.e., I nmax /I 2 = Ω (e anmax ) for some constant a.While an analytical characterization of mesoscopic localization in the general case of ν = 1 is out of the scope of this paper, we provide clear numerical evidence of localization phenomena in Fig. 8(c).The stationary distributions of the fraction of infected nodes in groups of increasing size n is concentrated in the largest group (n = 50) near the invasion threshold λ c .
Since mesoscopic localization was characterized using a linear contagion (ν = 1) and a continuous phase transition [48,49], two natural questions arise: How does ν = 1 affects localization?And what happens in the context of discontinuous phase transitions?In Fig. 9, we present the phase diagram of the group prevalence I n for different scenarios.Comparing Fig. 9(a) and Fig. 9(b), we see that increasing ν from 0.5 to 1.5 (while keeping g m = δ m,4 ) strengthens localization effects, which is expected since reinforcement effects are more important when the group prevalence is high.In Fig. 9(c), we show a similar diagram, but for a discontinuous phase transition.We see that the concentration of infected nodes in the largest groups is still possible, but the phenomenon is now associated with the unstable solution near the invasion threshold [λ = λ c (1 − ) with 1].Therefore, mesoscopic localization affects both continuous and discontinuous phase transition with a bistable regime, but the exponential growth of I n with n near λ c concerns the stable solution in the former and the unstable solution in the latter.
If we now reinterpret the results of Fig. 8 in light of these considerations, larger values of ν facilitate the onset of mesoscopic localization, where the largest groups drive the onset of the endemic phase, and make the invasion threshold scale as λ c ∼ n −ν max .This explains why λ c varies only slightly with γ n for ν = 2 in Fig. 8(a).
In Fig. 10, we finally investigate the role of heterogeneous group sizes on the bistability threshold by varying the group exponent γ n and the maximal group size n max .From Fig. 10(a), we see that a more heterogeneous group distribution, thereby increasing the fraction of larger groups, decreases the value of the bistability threshold ν c .This is consistent with our observation on regular structures [Fig.5], for which larger groups appear to promote bistability.However, Fig. 10(b) brings some nuance to this statement: for a fixed exponent γ n , there is an interesting non-monotonic relationship between ν c and the largest group n max .As such, the presence of larger groups does not always promote bistability.
We can again attribute this behavior to localization effects.In fact, we are able to illustrate this via a very simple example in Fig. 10(c).We look at the phase transition for a regular hypergraph with a fixed Impact of heterogeneous group sizes on the invasion threshold.We considered hypergraphs with power-law group size distributions pn ∼ n −γn with various exponents γn, a regular membership distribution gm = δm,4, and various exponents ν for the infection function in Eq. ( 1).We obtain the invasion threshold using Eq. ( 15).Lower values of γn implies a more heterogeneous group size distribution.(a) For a fixed nmax = 20, the invasion threshold increases with larger γn, but the effect is more limited for larger ν.(b) For a fixed γn = 3, the invasion threshold decreases like n −ν max for large nmax, indicating the onset of mesoscopic localization [48,49].(c) Stationary distribution fn,i for the starred case in (b), i.e., γn = 3 and ν = 1.5, with λ = 1.1λc, illustrating localization in the largest groups.n , the average stationary fraction of infected nodes in groups of various size n.We considered hypergraphs with power-law group size distributions pn ∼ n −γn with γn = 3, regular membership distributions of the form gm = δm,m 0 , and various exponents ν for the infection function in Eq. ( 1).Near the invasion threshold λc, I * n is larger for larger groups in (a) with ν = 0.5, but localization is much more pronouced for (b) with ν = 1.5.(c) For discontinuous phase transitions, mesoscopic localization is still possible, but near λc we must look at the unstable solution for I * n in the bistable regime.
group size, p n = δ n,4 , and a perturbed version of it, where we introduce a small proportion of larger groups, p n = (1 − )δ n,4 + δ n,15 with = 10 −3 .For the regular distribution, the phase transition is discontinuous, while for the perturbed distribution it is continuous, with the contagion localized in the largest groups near the invasion threshold.The bistability threshold ν c is larger for the perturbed distribution since mesoscopic localization reduces considerably the invasion threshold λ c .The largest most influential groups drive and self-sustain an endemic state for smaller values of λ, hence preventing a bistable regime.

C. Influence maximization
Influence maximization broadly refers to the problem of selecting a subset of nodes to initially spark a diffusion process in order to maximize the effect.The process could represent the spread of information, the diffusion of innovations or a viral marketing campaign [6,56].
The effectiveness of an influence maximization procedure is often measured by the fraction of affected nodes (in the limit t → ∞) for processes that terminate.How- FIG. 10.Impact of heterogeneous group sizes on the bistability threshold.We considered hypergraphs with power-law group size distributions pn ∼ n −γn with various exponents γn and a regular membership distribution gm = δm,4.We solve the bistability threshold using Eq. ( 21).Lower values of γn implies a more heterogeneous group size distribution.(a) For a fixed nmax, the bistability threshold increases with larger γn.(b) For a fixed γn, the bistability threshold has a non-monotonic relationship with nmax.(c) Phase transition using gm = δm,4 and ν = 2.3.We use pn = δn,4 for the regular case, and pn = (1 − )δn,4 + δn,15 with = 10 −3 for the perturbed case.Mesoscopic localization [48,49] inhibits bistability in the perturbed case.ever, because the final epidemic size in the SIS dynamics does not depend on the seeds (other than for stochastic extinction), we will consider the simpler task of maximizing İ(0), the initial spreading speed.This is often a straightforward task to solve for graphs.Considering the SIR model for instance, one just needs to maximize the number of outgoing edges from infected to susceptible nodes, which implies that nodes of maximal degree would be optimal influencers.However, we will show that additional considerations need to be accounted for in higher-order networks.More specifically, our goal is to use our formalism to answer the following question: Should we focus on finding influential nodes, or seed the spread from influential groups?
In this section, to simplify the notation, all dynamic quantities are evaluated at t = 0, e.g.I(0) ≡ I.
Let us assume that we are given a fixed hypergraph and an initial fraction of nodes that can be infected at the initial time I = 1 (the seeds of the contagion).Our task is to invade the system as fast as possible by maximizing İ for a hypergraph contagion, which is equivalent to maximizing the objective function where we define the initial node states S ≡ {s m } mmax m=1 and the initial group states The optimization problem is also constrained by While the first four constraints come from the definitions of the variables, the last one is less straightforward.Equation (25e) ensures the consistency between S and F, more specifically that the fraction of all memberships stubs belonging to susceptible nodes [left-hand side of Eq. (25e)] matches the fraction of susceptible nodes in groups [right-hand side of Eq. ( 25e)].By combining the constraint of Eq. (25e) with the definition of r as given by Eq. ( 2), the objective function can be simplified as Although it appears to be independent of S, it depends on it implicitly through Eq. (25e).
It is worth stressing that our formalism assumes that the membership stubs of nodes are assigned to groups uniformly at random, and thus we cannot engineer both S and F, i.e., choose at the same time the seeds according to their membership and the repartition of the seeds among the various group sizes.Indeed, if we decide for instance to infect only nodes of a certain membership m and we try to engineer F, there are no guarantees we can achieve such configuration in practice-e.g.we cannot infect a node in a group if none of its nodes have membership m .
We therefore compare two strategies to optimize the early spread: A) The influential spreaders strategy: we engineer S, i.e., we choose the fraction of seeds to assign to each membership class, and we assume a random configuration for the groups, i.e., all {f n,i } n i=0 are binomial distributions with probability q (to be determined).
B) The influential groups strategy: we engineer F, i.e., we assign a certain number of seeds in the groups depending on their sizes, and assume that nodes are infected at random through the group to which they belong.

Influential spreaders
In this strategy, we are free to engineer S in order to maximize Φ, with respect to the constraints of Eq. ( 25).Let us assume that f n,i is a binomial distribution, Using Eq. (25e), we can identify An optimal solution S can be found by first finding the value q that maximizes the objective function Eq. ( 26), and then identifying any set S that satisfies the relation for q = q above.There are in general many optimal solutions possible, but they collapse into a single one when q is sufficiently small, which is reasonable for 1.In this case, we simply have that Φ ≈ q, and the optimal solution is intuitive: one needs to infect nodes of maximal membership first in order to maximize q.Interestingly, this is true irrespective of β(n, i), p n and g m .
The infection function and the structure affect the maximal value of such that this solution is unique and optimal.For example, in the simplest case of linear contagion, where β(n, i) ∝ i, it is possible to show that this strategy is optimal up to q = 1/2 for all g m and p n , and we expect even higher values for ν > 1.For all practical purposes, targeting nodes of highest membership is optimal, and this is the case in all experiments we considered (see Sec. II C 3).

Influential groups
In this second strategy, we want to engineer F in order to maximize Φ with respect to the constraints Eq. ( 25).Let us assume that we can do so by choosing a certain number of groups and infecting a certain portion of their nodes.Following this procedure, one can realize that not all sets F satisfying Eq. ( 25) are allowed.For instance, if we decide to infect i nodes in all groups of size n, the outcome is different from just having f n,i = 1.Indeed, nodes belong to more than one group, hence we need to account for spillover effects-groups of size n = n would have some infected nodes as well, and more than i nodes could be infected in some groups of size n.
To do so, let us first define fn,i as the fraction of all the groups of size n for which we infect i nodes at random.Note that if a node belongs to multiple groups, it can be chosen more than once for infection, but the duplicates have no effect.Spillovers are taken into account by considering that each of the n − i nodes that have not been chosen for infection in a group of size n could have been infected in another group, with probability u (to be determined).Therefore, we can write where Second, let us define η as the fraction of all spots in groups that have been chosen for infection, Since nodes within groups are chosen at random, a node of membership m is susceptible if it has not been chosen for infection in any of the groups to which it belongs, i.e., As a consequence, η is constrained by Eq. (25c), The probability u still needs to be obtained.It corresponds to the fraction of all memberships that are not matched with a spot chosen for infection in a group, but that are still associated with an infected node: With this formulation, we engineer Since the objective function is a linear function of each fn,i , the optimization problem can be solved using linear programming.However, there is an intuitive and more efficient way to solve this problem exactly.We just need to identify the most cost-effective fn,i by looking at the effect on Φ of increasing fn,i , versus the cost of increasing fn,i , i.e., the variation of η The most cost-effective fn,i maximizes the ratio Obviously, i = 0 is always the most cost-effective for all n (since it has zero cost), but to satisfy Eq. ( 28), we must also fill some fn,i with i > 0.
Optimal solutions tend to fill the fn,i with i > 0 that maximizes R(n, i), especially for sufficiently small .A general solution can be obtained using algorithm 1 presented in Methods Sec.IV B, building on this idea of cost effectiveness.In the worst case, the computational complexity to obtain an optimal solution F under the influential groups strategy is O(m max + n 3 max ) when using the method presented in Methods Sec.IV B, which is much more efficient than using a general-purpose linearprogramming method.
Equation ( 29) also gives us an intuition of what defines influential groups when trying to maximize the early spread.If 1, then u 1, hence we have when considering β(n, i) = λi ν .For simple contagions (ν = 1), picking the largest group with a single seed (i = 1) is always optimal.For hypergraph contagions with ν > 1, the largest groups are the most influential as well, but the optimal number of seeds is generally i > 1.
Hence, beyond its size, the initial configuration of a group determines whether or not it is influential.

Experiments
To compare the influential spreaders and the influential groups strategies, we measure the ratio where Φ F and Φ S are the values of the objective function for the optimal solution of the influential groups and influential spreader strategies, respectively.Therefore, ζ > 1 indicates that the influential groups strategy is better to maximize İ, and vice versa if ζ < 1.
In the Supplementary Information, we show that where (n , i ) is the pair that maximizes the ratio R(n, i), restricted to i > 0, in the limit → 0. For general , we need to solve numerically the optimization problem as discussed in the previous sections.Interestingly, with β(n, i) = λi ν , ζ is independent of λ, since Φ ∝ λ.As a consequence, ζ is agnostic to the underlying phase of the system (healthy, bistable, or endemic).Equation (31) In Fig. 11, we illustrate how ζ varies as we change ν, , and the underlying structure.For homogeneous memberships and group sizes [Fig.11(a)], we see that the influential groups strategy performs better as soon as the contagion process is sufficiently nonlinear (ν ≈ 2); for highly nonlinear contagions (ν ≈ 4), the influential group strategy is much more effective, with ζ up to 100.When considering heterogeneous memberships, but still homogeneous group sizes [Fig.11(b)], the influential spreaders strategy performs better for moderately nonlinear contagions (ν 2.8), otherwise the influential groups strategy is still a better choice.Finally, considering a heterogeneous p n as well [Fig.11(c)] helps the performance of the influential groups strategy, especially for larger .
When picking a pair ( , ν) far from ζ = 1 [dashed lines in Figs.11(a), 11(b) and 11(c)], the strategy that maximizes İ invades the system faster, as can be seen in Figs 11(d) and 11(f).However, sufficiently close to ζ = 1, maximizing İ does not necessarily imply that I(t) will be larger for all t > 0. For instance, in Fig. 11(e), ζ ≈ 1, but the influential spreader strategy is slightly better.Therefore, one must be careful when interpreting the results of Fig. 11.One way to improve on our approach would be to consider higher-order temporal derivatives of I to assess which strategy performs best, or refine the optimization procedure by trying to maximize these higher-order derivatives as well.Interestingly, Figs.11(d)-(f) suggest that the initial speed, İ(0), roughly correlates with the time taken by the disease to infect a given fraction of the population, a metric that has been used to measure influence for SI and SIR dynamics [67][68][69].
Figure 11(f) also illustrates a particular feature of highly nonlinear contagions: the time to reach the stationary state can be excessively long for suboptimal strategies, despite λ = 3λ c .In this regime, the initial conditions have a much more important impact on the capacity of the contagion to invade the system, especially considering the possibility of stochastic extinction in real systems due to finite size.
To understand why this is the case, we show in Fig. 12 the time evolution of f nmax,i (t), the distribution for the number of infected in the largest groups, when the nodes are initially infected at random.For up to t = 50, only a few nodes are typically infected in the largest groups.Then we see the apparition of a bimodal distribution around t = 80, where either almost all nodes of the largest groups are infected or almost all nodes are susceptible.So instead of having a homogeneous transition, where all the largest groups have a similar state for all t, they diffuse asynchronously to states where almost all nodes are infected, at which point equilibrium is reached.This transition is driven by stochastic fluctuations, which explains the slow take-off.Conversely, the influential group strategy is by design optimizing the group configurations to facilitate this transition, which explains its good performance in Fig. 11(f).
It is also worth mentioning that the transition in Fig. 12 hinges on the nonlinearity of the contagion process-as soon as the number of infected nodes is sufficiently large in a group, nonlinear effects boost infections, resulting in almost all nodes being infected.Thus, nonlinear effects can be important locally, even though the global prevalence is still pretty low.With mean-field approaches that ignore dynamical correlations, nonlinear effects-or contributions from higherorder interactions-kick in only for a sufficiently large global prevalence [70].
These results again highlight the importance of considering an accurate description of the inner dynamics of groups when studying hypergraph contagions.In the context of influence maximization, optimizing group configurations is a crucial component; one should not focus exclusively on identifying the most central nodes.Ultimately, an optimal strategy would capitalize on the synergy of these two important aspects.

III. DISCUSSION
We have introduced group-based AMEs to describe hypergraph contagions.Our framework is analytically tractable, allowing us to obtain closed form implicit expressions for the critical and tricritical points.In addition, we have shown that it describes the dynamical process with remarkable accuracy when compared with Monte Carlo simulations.Our formulation in terms of an infection rate function β(n, i) makes it extremely flexible, allowing us to consider arbitrary group distribution with large group interactions, contrarily to existing HMF theories [30,38,39] which instead require specifying the rule for each different type of interaction separately.
Motivated by simplicity and recent results [36], we analyzed in depth the consequences of a nonlinear infection rate function β(n, i) = λi ν , highlighting the important role of influential groups in hypergraph contagions.
With our analytical results about the invasion and bistability thresholds, we were able to perform an exhaustive analysis of the phase transition and better understand the influence of a heterogeneous structure, both in terms of membership m and group size n.We found that the third moment of the membership distribution g m plays a crucial role, with large m 3 suppressing the onset of a discontinuous phase transition with a bistable regime, in line with other approaches [38,39].This is best exemplified for power-law membership distributions g m ∼ m −γm , where γ m = 3 most suppresses bistability, and in the limit m max → ∞, a discontinuous phase transition is only possible for γ m > 4.
The phenomenon of mesoscopic localization [48,49], driven by the most influential groups, also has important consequences on the phase diagram, with the effects be-ing enhanced by superlinear infection (ν > 1).In this case, the invasion threshold scales as λ c ∼ n −ν max , and for λ close to λ c , infected nodes are found almost exclusively in the largest groups.This localization of the contagion thereby inhibits bistability by enforcing an endemic state with a very small global fraction of infected nodes.
Our approach, furthermore, provided novel insights concerning the problem of influence maximization for hypergraph contagions.We focused on the problem of maximizing the early spread, and proposed two strategies: allocating seeds to the influential spreaders (engineering s m ), or to the influential groups (engineering f n,i ).For various types of structures, the latter strategy performs better for contagions that are sufficiently nonlinear, highlighting the key role of influential groups on the transient state of the system.
For the process we considered, the notion of influential groups to seed and sustain hypergraph contagions are mostly aligned-in both cases, the largest groups typically have a dominant role.In the case of influence maximization, however, we showed that a careful seed allocation is also essential to determine whether or not a group is influential.Moreover, a more realistic infection function β(n, i) that actually depends on n could affect which groups are most influential in both scenarios.
Our work constitutes a first step towards a better understanding of the role of higher-order interactions on the outreach of information spreading [6], and resonates with other recent theoretical findings on higher-order naming games, where big groups facilitate takeover of committed minorities in social convention [29].Interestingly, AMEs thus provide an analytical avenue to study recent empirical results showing how social contagions and movements defy classic influence maximization.As one example, networked counterpublics [71] are public spaces used by underrepresented groups to gather legitimacy and form tight-knit communities.Therein, non-dominant forms of knowledge can still spread and reach widespread attention through dense communities (influential groups) despite the limited connectivity of their members (noninfluential spreaders).These results provide one more addition to the mounting evidence that groups of elementary elements are the foundational unit of many complex systems.
Many avenues are now left open to explore and broaden the applicability of our group-based AMEs.While we restrained ourselves to a particular nonlinear infection rate function β(n, i) and a constant recovery rate, other dynamical processes could be considered, each having their own phenomenology and a rich dynamical behavior.In the Supplementary Information for instance, we briefly discuss how our framework can be applied to threshold models of the form β(n, i) = δ n−1,i , but one could consider other traditional dynamical processes, such as voter models [72].
In the Supplementary Information, we provide a roadmap to include structural two-point correlations, but a thorough characterization of the impact of correlation patterns on bistability, mesoscopic localization and influence maximization is still lacking.The inclusion of dynamical correlation around nodes is a more tedious task that would require a fusion between degreebased [46,47,73] and group-based [19,20,48,49] AMEs.This would allow to describe almost exactly short-range dynamical correlations, namely correlations between the states of nodes and their direct neighbors.Incorporating long-range correlations-beyond first neighbors-in AME frameworks, without a prohibitive computational time due to combinatorial explosion, is still an open problem.
Finally, many directions could be taken with regards to the influence maximization problem on hypergraphs.One avenue would be to analyze the notion of influential groups and influential spreaders from the perspective of centrality measures for hypergraphs [74].Another would be to investigate the closely related problem of targeted immunization [1,75].

Simulation of contagions
We used a standard Gillespie algorithm for the simulation of contagions on hypergraphs.We decompose the whole process into events j ∈ J, that each happens at rate ω j .The next event to happen is chosen with probability and the time step between two events ∆t is distributed exponentially with mean ∆t = 1/ j∈J ω j .There are two types of events: infection and recovery.On the one hand, all susceptible nodes in a group can be considered equivalent with regards to infection.Consequently, each group is chosen for an infection event with rate Once a group is chosen for an infection event, one of the (n − i) susceptible nodes is chosen uniformly at random to become infected.On the other hand, all infected nodes perform a recovery event with rate ω rec = 1.
We store all possible events in an efficient data structure called a SamplableSet [76], where insertion, deletion and sampling of elements (events) all have a computational complexity O [log log (ω max /ω min )] [77], where ω max and ω min are respectively the maximal and minimal rate among {w j } j∈J .This makes the sampling and the updating of the data structure extremely fast, which is especially useful when {w j } j∈J spans multiple scales.
Once an event is performed-for instance a node recovers-we need to update the rate ω inf of all groups to which this node belongs.This is the most costly part of the algorithm, which unfortunately cannot be overcome.This essentially means the simulation procedure is slower for hypergraphs with large average excess membership m(m − 1) / m .In Fig. 2 and 3, we compare the stationary state solutions from our formalism with estimates from Monte Carlo simulations.To compute estimates, we let the system relax during a burnin period τ b ∈ [10 2 , 10 4 ] then we sample N ∈ [10, 10 4 ] states, both depending on the size of the hypergraph and if multiple randomized hypergraphs are being used.Sampled states are separated by a decorrelation period τ d = 1.
To simulate contagions in the stationary state, we used two approaches, ordinary simulations and the quasistationary-state method.
a. Ordinary simulation method.With this approach, we simply let the simulation run and do not intervene.This is usually not the method of choice, especially for small hypergraphs near the invasion threshold, because finite size systems all eventually reach the absorbing state where all nodes are susceptible.This is, however, more practical to obtain the lower branch in Fig. 2(f), or faster for large hypergraphs, as in Fig. 3.
b. Quasistationary-state method.This approach aims at sampling the quasistationary distribution of the contagion process [78], which is defined as the probability distribution for all states in the limit t → ∞, except for the absorbing state.We used the state-of-the-art method discussed in Ref. [78], where we keep a history of past states (in our case up to 50 states).We update the history by removing one uniformly at random and storing the current state after each decorrelation period τ d ∈ [0.1, 1].Each time the absorbing state is reached during the simulation, we pick a state from the history uniformly at random to replace the current one.This method is well suited for finite-size analysis, and especially useful for simulations on small hypergraphs, such as in Fig. 2.

Datasets
The simulations shown in Fig. 2 and Fig. 3 run on two different empirical social structures that encode different types of social higher-order interactions.Here, we briefly describe the nature of these two datasets and the techniques used to construct the associated higher-order structures.
a. Face-to-face interactions in a French primary school.Originally collected in Ref. [79] as part of the SocioPatterns [80] collaboration, this dataset contains in-formation of face-to-face interactions between children of a primary school in Lyon recorded over 2 days.Participants are given wearable sensors (placed on their chests), and a contact is detected whenever two sensors are in close proximity (1.5m).The initial temporal resolution of this dataset is 20s, but contacts have been further preprocessed as in Ref. [30] in order to construct a static hypergraph from the temporal sequence of interactions.In particular, considering each child as a node, we aggregated different snapshots using a temporal window of 15 minutes and computed all the maximal cliques appearing in each window.Cliques were then aggregated across the entire time range, retaining only those that appeared at least 3 times, and finally "promoted" to groups.Some properties of the obtained structure are reported in the caption of Fig. 2.
b. Co-authorship in computer science.DBLP [81] is an online bibliography containing information on major computer science journals and proceedings.This dataset, already pre-processed in Ref. [82] (from the release 3, 2017), consists in a list of publications and respective authors that naturally calls for higher-order representations [83].In particular, each author corresponds to a node and any collaboration of n (co-)authors in a single publication corresponds to a group of size n.We constructed a hypergraph by aggregating all the resulting groups together, but without considering single author publications (these have been removed in order to avoid disconnected nodes).In addition, the original dataset contained 1 831 127 nodes and 2 954 518 groups, which is too large to perform simulations on a personal computer in a reasonable time.Therefore, we obtained a subhypergraph by performing a breadth-first search, starting from a random group, then visiting all groups at a maximum distance of 2 when considering the one-mode projection of the original hypergraph on the groups.This ensures that the resulting sub-hypergraph is connected.Some properties of the obtained structure are reported in the caption of Fig. 3.

Randomization and data augmentation
In Fig. 2 and 3, we make use of randomized versions of the original hypergraph.In Fig. 2, we also use expanded versions, where the size of the network is increased by a factor x. In all cases, we use the same procedure (x = 1 if the hypergraph is not expanded).
Let us first note m = [m 1 , m 2 , . . .] and n = [n 1 , n 1 , . . .] the membership sequence and the group size sequence of the original hypergraph, i.e., the list for the membership of each node and the list for the size of each group.From these sequences, we create two expanded sequences m and n , which are formed of x copies of m and n respectively.This can be seen as the membership and group size sequences for a hypergraph that is x times larger.
For each expanded sequence, we create a stub list.For instance, for the node j of the expanded hypergraph, we include m j copies of the label j in the stub list for the nodes.Similarly, we include n copies of the label in the stub list for the groups.By definition, these two stub lists are of the same length, M , which corresponds to the number of edges in the bipartite representation of the hypergraph.We can thus shuffle them and match the entries of both lists, thereby assigning nodes to groupsor equivalently creating edges between nodes and groups in the bipartite representation of the hypergraph.
We then remove multi-edges (nodes assigned multiple times to a same group) by performing edge swaps [84].We then perform M additional edge-swap attempts at random-picking two random edges, swapping the groups, and accepting the swap if it does not create multiedges.This ensures the uniformity of the generation process (see the Supplementary Information).The resulting hypergraph is a randomized version of the original hypergraph, expanded by a factor x.

B. Influential groups solutions
An intuitive approach to solve the problem would be to sort all pairs (n, i) in decreasing order of their R(n, i) values (for i > 0), then fill fn,i up to 1 following this order, until I = , or more directly until η reaches the value prescribed by .However, this approach does not account for the fact that one may encounter multiple times a same n value before the condition I = is reached.For instance, let us assume (n, i) is the next the pair with the highest value R(n, i), but there exists a pair (n, i ) with R(n, i ) ≥ R(n, i) and we have already assigned fn,i = 1.What is the best option? 1. Discard the (n, i) pair.
2. Fill the associated fn,i up to 1 and decrease the value of fn,i accordingly.
It turns out that an optimal solution is constructed by choosing one or the other depending on certain conditions.Option 1 is chosen whenever i < i , because it can only reduce the total contribution to Φ.If i > i , we assign a new cost-effective ratio to the pair (n, i), accounting for the fact that we need to decrease fn,i : R(n, i) = iR(n, i) − i R(n, i ) i − i .
This can be interpreted as the cost-effective ratio for the additional infected nodes (i − i ) that we add to the configuration.Note that R(n, i) can be negative, which is not a problem: this only means that infecting these nodes decreases the objective function Φ.If R(n, i) is still the highest ratio when compared with the ratios from other available pairs, option 2 is chosen.
The procedure to find an optimal set F is presented in algorithm 1, where we use a priority queue Q to keep sorted the pairs (n, i) in decreasing order of their costeffective ratio.Using a binary heap for Q, inserting any pair (n, i) or removing the pair with maximum costeffective ratio takes O(log |Q|) elementary operations in the worst case, where |Q| is the number of elements in the queue.Since |Q| ≤ n 2 max , insertions and removals have a complexity O(log n max ).In the worst case, we will perform O(n 2 max ) insertions and removals before the condition I = is reached, hence the worst-case complexity for algorithm 1 is O(n 2 max log n max ).We also need to precompute R(n, i) for all n, i > 0 and get F from F , both O(n 3 max ).Moreover, solving for η is O(m max ), so the overall computational complexity of the influential groups optimization procedure is O(m max + n 3 max ).

Algorithm 1 Optimization for influential groups
Input: pn, R(n, i), η Output: F 1: fn,i ← 1 if i = 0, else 0, ∀n, i 2: ψ ← η n Initial resources 3: iopt ← ∅ Associative array 4: Q ← ∅ Priority queue 5: for all (n, i) s.t n ≤ nmax, pn > 0 and 1 ≤ i ≤ n do  For the first term, we define the following quantities After some manipulations, one can write Combining all equations, we are able to compute the Jacobian J r in O(m max + n 2 max ) elementary operations.

FIG. 2 .
FIG.2.Fraction I * of infectious nodes in the stationary state of a contagion model on a hypergraph constructed from high-resolution face-to-face contact data from a primary school in Lyon (see Methods Sec.IV A 2).The hypergraph contains 242 nodes and 1188 groups.Both the membership and the group size distributions are homogeneous, with m ≈ 11.79, n ≈ 2.40, mmax = 32 and nmax = 5.We compare the results of Monte Carlo simulations (symbols; see Methods Sec.IV A 1) with the predictions of our approach (solid and dashed lines for stable and unstable solutions respectively).We rescale λ with the invasion threshold λc, which is computed using Eq.(15) in Sec.II B. The symbols represent the average infected fraction measured over long runs and averaged over randomized hypergraphs when this applies.The error bars (sometimes too small to be seen) correspond to one standard deviation.The green circles are obtained with the quasistationary-state method, starting with a large fraction of infected nodes (I = 0.8) to sample the upper branch of the hysteresis loop when the phase transition is discontinuous.(a),(d) We use the original hypergraph.(b),(e) We use 10 uniformly randomized versions of the original hypergraph.(c),(f) We use 10 uniformly randomized versions of the original hypergraph with its size expanded by a factor 10 (see Methods Sec.IV A 3). (f) The blue triangles are obtained by ordinary simulations, starting with a small fraction of infected nodes (I = 0.02) to sample the lower branch of the hysteresis loop.We only do it for the expanded hypergraphs because finite-size effects makes unreliable the estimation of the lower branch with data of small size.
FIG.2.Fraction I * of infectious nodes in the stationary state of a contagion model on a hypergraph constructed from high-resolution face-to-face contact data from a primary school in Lyon (see Methods Sec.IV A 2).The hypergraph contains 242 nodes and 1188 groups.Both the membership and the group size distributions are homogeneous, with m ≈ 11.79, n ≈ 2.40, mmax = 32 and nmax = 5.We compare the results of Monte Carlo simulations (symbols; see Methods Sec.IV A 1) with the predictions of our approach (solid and dashed lines for stable and unstable solutions respectively).We rescale λ with the invasion threshold λc, which is computed using Eq.(15) in Sec.II B. The symbols represent the average infected fraction measured over long runs and averaged over randomized hypergraphs when this applies.The error bars (sometimes too small to be seen) correspond to one standard deviation.The green circles are obtained with the quasistationary-state method, starting with a large fraction of infected nodes (I = 0.8) to sample the upper branch of the hysteresis loop when the phase transition is discontinuous.(a),(d) We use the original hypergraph.(b),(e) We use 10 uniformly randomized versions of the original hypergraph.(c),(f) We use 10 uniformly randomized versions of the original hypergraph with its size expanded by a factor 10 (see Methods Sec.IV A 3). (f) The blue triangles are obtained by ordinary simulations, starting with a small fraction of infected nodes (I = 0.02) to sample the lower branch of the hysteresis loop.We only do it for the expanded hypergraphs because finite-size effects makes unreliable the estimation of the lower branch with data of small size.

FIG. 3 .
FIG.3.Fraction I * of infectious nodes in the stationary state of a contagion model on a hypergraph constructed from coauthorship data in computer science (see Methods Sec.IV A 2).The hypergraph has been obtained by using a breadth-first search to extract a sub-hypergraph of 57 501 nodes and 55 204 groups; the original data contained more than 10 6 nodes and groups.(a) The membership distribution is heterogeneous, approximately of the form gm ∼ m −γm , while the group size distribution is more homogeneous, with m ≈ 3.75, n ≈ 3.90, mmax = 903 and nmax = 25.(b),(c) We compare the results of Monte Carlo simulations (circle markers; see Methods Sec.IV A 1) with the predictions of our approach (solid lines).We rescale λ with the invasion threshold λc, which is computed using Eq.(15) in Sec.II B. The infected fraction I * has been obtained as averages over time with long runs (and averaging over randomized hypergraphs when this applies).The error bars (too small to be seen) correspond to one standard deviation.We used ordinary simulations, starting with I = 0.8.(b) We use the original hypergraph.(c) We use 10 uniformly randomized versions of the original hypergraph.
FIG.3.Fraction I * of infectious nodes in the stationary state of a contagion model on a hypergraph constructed from coauthorship data in computer science (see Methods Sec.IV A 2).The hypergraph has been obtained by using a breadth-first search to extract a sub-hypergraph of 57 501 nodes and 55 204 groups; the original data contained more than 10 6 nodes and groups.(a) The membership distribution is heterogeneous, approximately of the form gm ∼ m −γm , while the group size distribution is more homogeneous, with m ≈ 3.75, n ≈ 3.90, mmax = 903 and nmax = 25.(b),(c) We compare the results of Monte Carlo simulations (circle markers; see Methods Sec.IV A 1) with the predictions of our approach (solid lines).We rescale λ with the invasion threshold λc, which is computed using Eq.(15) in Sec.II B. The infected fraction I * has been obtained as averages over time with long runs (and averaging over randomized hypergraphs when this applies).The error bars (too small to be seen) correspond to one standard deviation.We used ordinary simulations, starting with I = 0.8.(b) We use the original hypergraph.(c) We use 10 uniformly randomized versions of the original hypergraph.

FIG. 6 .
FIG. 6. Impact of the third moment of the membership distribution on the phase transition.The group size distribution is fixed to pn = δn,3.(a) Membership distributions with the same first two moments ( m ≈ 5.03 and m 2 ≈ 30.2) but with a different third moment m 3 .The distribution with higher m 3 is of the form gm ∼ (a + m) −b e −m/c with (a, b, c) = (−0.8,3.72, 40) and ranges from mmin = 4 to mmax = 100.The one with lower m 3 is Poisson distributed gm ∼ a m e −a /m! with a = 5 and ranges from mmin = 1 to mmax = 100.(b) Non-trivial solutions to Eq. (7) for the corresponding membership distribution.Solid and dashed lines represent stable and unstable solutions, respectively.The infection rate function considered is β(n, i) = λi ν with ν = 2.8.

4 FIG. 7 .
FIG.7.Impact of heterogeneous memberships on the bistability threshold.We considered hypergraphs with power-law membership distributions gm ∼ m −γm with various exponents γm, different maximal values mmax, a minimal value mmin = 2, and a regular group size distribution pn = δn,4.We obtain the bistability threshold by solving Eq.(21).Lower values of γm imply a more heterogeneous membership distribution.(a) The bistability threshold grows logarithmically or faster for γm ≤ 4, but converges for γm > 4. (b) The bistability threshold shows a maximum near γm = 3, suggesting that this is the maximum in the limit mmax → ∞.
FIG.8.Impact of heterogeneous group sizes on the invasion threshold.We considered hypergraphs with power-law group size distributions pn ∼ n −γn with various exponents γn, a regular membership distribution gm = δm,4, and various exponents ν for the infection function in Eq.(1).We obtain the invasion threshold using Eq.(15).Lower values of γn implies a more heterogeneous group size distribution.(a) For a fixed nmax = 20, the invasion threshold increases with larger γn, but the effect is more limited for larger ν.(b) For a fixed γn = 3, the invasion threshold decreases like n −ν max for large nmax, indicating the onset of mesoscopic localization[48,49].(c) Stationary distribution fn,i for the starred case in (b), i.e., γn = 3 and ν = 1.5, with λ = 1.1λc, illustrating localization in the largest groups.

10 FIG. 9 .
FIG.9.Mesoscopic localization in large groups, illustrated by I * n , the average stationary fraction of infected nodes in groups of various size n.We considered hypergraphs with power-law group size distributions pn ∼ n −γn with γn = 3, regular membership distributions of the form gm = δm,m 0 , and various exponents ν for the infection function in Eq. (1).Near the invasion threshold λc, I * n is larger for larger groups in (a) with ν = 0.5, but localization is much more pronouced for (b) with ν = 1.5.(c) For discontinuous phase transitions, mesoscopic localization is still possible, but near λc we must look at the unstable solution for I * n in the bistable regime.

FIG. 11 .
FIG. 11.Comparison of the influential spreaders and influential groups strategies.(a)-(c) If log 10 ζ > 0, this indicates that the influential groups strategy is better to maximize İ, and vice versa if log 10 ζ < 0. We use different combinations of homogeneous and heterogeneous distributions.More specifically, for homogeneous distributions we use gm ∝ a m e −a /m! and pn ∝ a n e −a /n!, with a = 5, mmin = 1, nmin = 2, and mmax = nmax = 20.For heterogeneous distributions we use gm ∝ m −γm and pn ∝ n γn , with γm = γn = 3, mmin = nmin = 2, and mmax = nmax = 100.The star markers indicate when ζ = 1 in the limit → 0, using Eq.(32).Irregularities of the level curves are due to the discrete nature of gm and pn, not to numerical errors.(d)-(f) Time evolution of the fraction of infected nodes for different strategies, with = 10 −2 and ν ∈ {1, 1.88, 3}, corresponding to the three empty markers in (a).We use λ = 3λc in each case.The random strategy corresponds to sm = 1 − for all m, and fn,i is a binomial distribution with probability .

50 FIG. 12 .
FIG.12.Time evolution of the distribution fn max ,i when using the random strategy.The structure, the dynamics and the parameters are identical to the ones of Fig.11(f).

7 :Q← i 16 : 3 IV. Threshold models 4 V 5 References 6 I
← (n, i) with priority R 8: end for 9: while Q is not empty and ψ > 0 do 10:else if i > iopt[n] then 17: i = iopt[n] 18: R ← [iR(n, i) − i R(n, i )]/(i − i ) 19: if R > R for all R ∈ Q then 20: fn,i ← min{1, ψ/[pn(i − i )]} 21:fn,i ← fn,i − fn,i 22: ψ ← ψ − fn,ipn(i − i ) assessment of the effect of dynamical correlation around nodes 2 III.Jacobian of the mean-field map .Uniformly sampling stub-labeled simple bipartite graphs 4 VI.Influence maximization: limit case .AME CORRECTION FOR CORRELATIONS It is straightforward to adapt our formalism to include two-point structural correlations between node memberships and group sizes.Let P (m, n) be the joint probability distribution for the membership of a node and the size of a group at the endpoint of an edge in the bipartite representation of the hypergraph.Marginalization results in m P (m, n) = np n n and n P (m, n) = mg m m .III. JACOBIAN OF THE MEAN-FIELD MAP All quantities are assumed to have reached the stationary sate.As shown in Sec.II B, the stationary state solutions are completely characterized by the mean field r.Therefore, one just needs to obtain solutions r that respect r = M[ρ(r)] .A useful quantity to find the fixed points numerically and analyze their stability is the Jacobian df n,i dρ = −f n,i j>0 f n,j ∆ n,j + f n,i ∆ n,i , , j) + ρ .For the second term we follow the same steps, starting with the definitionsu 2 (r) = m m(m − 1)s m g m , v 2 (r) = m ms m g m ,We insert the following in the previous definitions ds m dr = − m (1 + mr) 2 .