Social contagion models on hypergraphs

PHYSICAL REVIEW RESEARCH 2, 023032 (2020) Social contagion models on hypergraphs Guilherme Ferraz de Arruda ,1Giovanni Petri,1and Yamir Moreno 2,3,1 1ISI Foundation, Via Chisola 5, 10126 Torino, Italy 2Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, 50018 Zaragoza, Spain 3Department of Theoretical Physics, University of Zaragoza, 50018 Zaragoza, Spain (Received 4 October 2019; accepted 10 March 2020; published 10 April 2020) Our understanding of the dynamics of complex networked systems has increased significantly in the last two decades. However, most of our knowledge is built upon assuming pairwise relations among the system’s components. This is often an oversimplification, for instance, in social interactions that occur frequently within groups. To overcome this limitation, here we study the dynamics of social contagion on hypergraphs. We develop an analytical framework and provide numerical results for arbitrary hypergraphs, which we also support with Monte Carlo simulations. Our analyses show that the model has a vast parameter space, with firstand secondorder transitions, bistability, and hysteresis. Phenomenologically, we also extend the concept of latent heat to social contexts, which might help understanding oscillatory social behaviors. Our work unfolds the research line of higher-order models and the analytical treatment of hypergraphs, posing new questions and paving the way for modeling dynamical processes on higher-order structures. DOI: 10.1103/PhysRevResearch.2.023032 Network science has had a radical impact on our knowledge about critical dynamics in complex systems. In particular, new and relevant phenomena arise when investigating social and biological contagion processes [1–7]. For instance, while homogeneous spreading models predict finite critical points [2,6,7], heterogeneous networks often present vanishing transitions [1,2,6–8], supporting the predictions in real-world networks [9–12]. Contagion models cover many aspects, from different types to richer substrates underlying the process itself. A relevant development is the extension of contagion processes to multilayer networks, leading the way to combinatorial higher-order models. Indeed, multilayers’ structural [13–17], spreading [7,14], and diffusion properties [13,18] have new and richer phenomenology. Nevertheless, as recently argued in Ref. [19], real data are revealing that pairwise relationships—the fundamental interaction units of networks—do not capture complex dependencies. Indeed, modern messaging systems (e.g., WhatsApp, Telegram, and Facebook Messenger, among others) allow users to communicate in groups, which creates a direct channel among all members. In other words, modern information spreading is often a one-to-many process. Additionally, biological and team collaborations are also inherently group structured, similarly to some types of molecular interactions [20]. Complementary, there is evidence from social and biological studies indicating that higher-order structures have crucial dynamical effects [21–24]. Therefore, higher-order interactions are ubiquitous, and understanding their properties and impacts is of paramount importance. Notably, the sizes of such groups can Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. span orders of magnitude. Thus, graph-projection-based approaches might not be sufficient to describe systems involving interactions over many different scales and orders. Combinatorial higher-order models [19] offer a way to describe these systems as they overcome the limitations of lower-ordernetworkmodels.Inafirstattempt,Bodóetal.[21] proposed an SIS disease spreading in a hypergraph. Next, Iacopini et al. [22] presented a model of social contagion defined on simplicial complexes and provided approximate solutions for complexes of order three. Their model presented new phenomenological patterns associated with the critical properties of the dynamics. However, their proposed model is still very constrained, both structurally and dynamically. Here we extend their model both structurally and dynamically. Structurally, we adopt hypergraphs, which generalize the concept of graphs, by allowing an edge to have an arbitrary number of nodes (see Fig. 1for an illustration). Hypergraphs relax the structural restrictions required by simplicial complexes as they impose virtually no limitation on the type, size, and mutual inclusion of interactions, thus, representing more faithfully and naturally real systems. From the dynamical viewpoint, we incorporate explicit critical-mass dynamics (each hyperedge is an independent critical-mass process), generalizing the one modeled in Ref. [22]. The resulting model displays a rich complex phenomenology, remaining very flexible and able to cover a wide range of systems. We uncover the presence of discontinuous transitions and bistability led by higher-order interactions and critical-mass dynamics. Notably, these critical properties contrast with contagion models on graphs-models, which instead usually display continuous transitions, e.g., SIS or SIR disease spreading, regardless of the structural configuration. Therefore, assuming a graph projection of a hypergraph might lead to wrong results. Here we report analytical and numerical analyses of the theoretical framework that we introduce as well as results for several limiting cases and hypergraph structures. We round 2643-1564/2020/2(2)/023032(6) 023032-1 Published by the American Physical Society DE ARRUDA, PETRI, AND MORENO PHYSICAL REVIEW RESEARCH 2, 023032 (2020) FIG. 1. Graphical representation of a hypergraph. Mathematically, V={v1,v2,v3,v4,v5,v6,v7,v8},E={e1,e2,e3,e4},where the hyperedges are e1={v1,v2,v3},e2={v3,v4,v5,v6},e3= {v6,v7}and e4={v8}. off the paper by discussing several implications of our study and, most notably, the role of critical mass dynamics in social contagion, providing insights that could help explain reported differences in experimental results [24–27]. Let us first introduce some formal definitions. A hypergraph is defined as a set of nodes, V={vi}, where N=|V| is the number of nodes, and a set of hyperedges E={ej}, where ejis a subset of Vwith arbitrary cardinality |ej|.If max (|ej|)=2 we recover a graph. We remark that we do not explore this scenario since it recovers standard models. On the other hand, if for each hyperedge with |ej|>2 its subsets are also contained in E, we recover a simplicial complex (for more on the hypergraph structure, see the Supplemental Material [28], Sec. I). Figure 1shows an example of a hypergraph. In an arbitrary hypergraph, we associate with each individual via Bernoulli random variable Yi(complementary Xi). If the node viis active Yi=1(Xi=0), otherwise Yi=0 (Xi=1). To each active node, we associate a deactivation mechanism, modeled as a Poisson process with parameter δi, Nδi i(Yiδi −→ Xi). For each hyperedge, j, we define a random variable Tj=k∈ejXk, which is the number of active nodes in the hyperedge. If Tjis equal or above a given threshold, j, we associate a Poisson process with parameter λj,Nλj j (that is, if Tj⩾j, then Xk λj −→ Yk,∀k∈ej). In other words, the dynamics is given by a threshold process that becomes active only above a critical mass of activated nodes. Moreover, if |ej|=2, we assume directed Poisson processes, implying that it is not a threshold process anymore. This definition allows for recovering traditional SIS contagion models. While the proposed model is general, allowing for arbitrary heterogeneity in parameters, we focus on more straightforward, but representative, cases. We assume that δi=δand λj=f(|ej|), where fis an arbitrary function of the cardinality of the hyperedge. It is also convenient to define =∗N, where the parameter ∗is a real number representing the fraction of active nodes. The exact equation that describes the aforementioned dynamics can be written as dE(Yi) dt =E⎛ ⎝−δYi+[1−Yi] ej∩{vi}=∅ λj B 1{Yi=0,Tj⩾j}⎞ ⎠, (1) where the first summation is over all hyperedges containing viand the second over all the possible dynamical microstates inside the hyperedge ej, denoted by the set B. Furthermore, 1{Tj⩾j,Yi=0}is an indicator function that is 1 if Yi=0 and the critical mass in the hyperedge is reached and 0 otherwise. Naturally, the order parameter is defined as the expected fraction of active nodes, i.e., ρ=1 NiE(Yi). Although Eq. (1) captures the exact process, it cannot be numerically solved. Thus, assuming that the random variables are independent and denoting yi=E(Yi), we obtain the firstorder approximation, given as dyi dt =−δyi+λ(1−yi) ej∩{i}=∅ |ej|  k=j λ∗(|ej|)Pej(K=k), (2) where we assume that the spreading rate is composed by the product of a free parameter and a function of the cardinality, i.e., λj=λ×λ∗(|ej|). In this formulation, the expectation of the indicator function follows a Poisson binomial distribution (for more on this approximation, see Ref. [28], Sec. II). Formally, E(1{(Tj−Yk)⩾j})≈ |ej|  m=j Pej(K=l),(3) Pej(K=l)= A∈Fl i∈A yi j∈Ac (1 −yj),(4) where Flis the set of all subsets of kintegers from {1,2,...,n=|ej|},Ais one of those sets, and Acis its complementary. Intuitively, Aaccounts for the possibly active nodes and Acthe possibly inactive ones. Thus, the summation over Flconsiders all possible nodal state configurations in a given hyperedge. Equation (4) is not numerically stable if |ej| is large [29]. Considering the discrete Fourier transform, we obtain a numeric stable solution as [29] Pej(K=k)=1 n+1 n  l=0 C−lk n  m=1 [1 +(Cl−1)ym],(5) where C=exp ( 2iπ n+1), which then allows one to compute the solution for arbitrarily large hyperedges. Interestingly, although the whole argument is quite intricate, Eq. (5)is simple, allowing the numerical evaluation of Eq. (2). Our main result is that a rich and diverse phase space, generally populated by continuous and discontinuous transitions and hysteretic behaviors, characterizes contagion on hypergraphs. In particular, we have analytically observed discontinuity and bistability in the order parameter on top of some regular structures. We provide full details of the calculations in Ref. [28] (see Secs. III and IV) for two limiting cases, namely, a hypergraph composed of a hyperedge containing all nodes in addition to (1) a random regular network (which we call a hyperblob) and (2) a star (referred to as a hyperstar). For the sake of clarity, let us show the main results for the hyperblob. In this case, we can exploit the structural symmetries to solve ρ(λ, λ∗,δ), obtaining two locally stable solutions. Specifically, consider a hypergraph built up as a 023032-2 SOCIAL CONTAGION MODELS ON HYPERGRAPHS PHYSICAL REVIEW RESEARCH 2, 023032 (2020) homogeneous set of pairwise interactions with average degree kand a single additional hyperedge containing all nodes. In this case, the order parameter can be solved as (for the full derivations, see Ref. [28], Sec. III) ρLower =1−δ kλ,if λ δ⩾1 k 0,otherwise (6) ρUpper =−δ+kλ−λ∗λ+4kλ∗λ2+[δ+(−k+λ∗)λ]2 (2kλ),(7) where a second-order phase transition for ρLower is naturally obtained as λ δ⩾1 k[7]. Furthermore, the discontinuities can also be calculated as λL c=δ k−∗k,(8) λU c=− δ∗ λ∗∗−λ∗+(∗)2k−∗k.(9) Phenomenologically, a discontinuity implies that our system possesses a “social latent heat” that is released or accumulated at a constant value of λ. More specifically, before the discontinuity, “energy” has been stored in the partial activation of the hyperedges. At the discontinuity this “energy” is absorbed (released) at once for a constant value of λ. In fact, the social latent heat can be expressed as QlλX c=[ρUpper(λ, δ, λ∗,N)−ρLower(λ, δ, λ∗,N)]λ=λX c,(10) where Ql(λX c) can be Ql(λL c) (energy absorbed) or Ql(λU c) (energy released). Therefore, for this structure, the latent heat is expressed as QlλX c=δ−λ(λ∗+k)+[δ+λ(λ∗−k)]2+4λ∗kλ2 2kλλ=λX c ,(11) where λX ccan be (λL cor λU c). In fact, this expression is true for any value of λ, but its physical interpretation is valid only near the discontinuity, which in turn depends on ρc=∗. We refer the reader to Ref. [28] for more details. Figure 2shows the general phenomenology of the system obtained from the analytical solution, i.e., the first-order approximation, of the equations describing the contagion dynamics for the hyperblob. As can be seen in Fig. 2(a), 0.0 0.2 0.4 0.6 0.8 1.0 λ 0.0 0.2 0.4 0.6 0.8 1.0 ρ Ql(λL c) Ql(λU c) λc (a) ρLower ρUpper ρc(Θ∗,λ ∗(|e|)) ρ ρ∗ 0.0λcλ λ 0.0 0.2 0.4 0.6 0.8 1.0 Θ∗ (b) I II III IV V FIG. 2. Results for the hyperblob. (a) Possible solutions for a fixed ∗=0.5. In red and blue are the upper and lower solutions (branches), respectively. The transition from the lower to the upper solution (upper to lower) occurs at the intersection of the lower (upper) solution with a value of ρcin which the upper solution became stable (unstable). The discontinuity is characterized by the latent heat, Ql(λL c)orQl(λU c). At λc=0.2, the lower solution shows a second-order phase transition. (b) Schematic of the parameter space: Region I: the absorbing state for both the lower and upper solution; Region II: only the lower solution is stable (the global critical mass is not reached, ρ<ρ c); Region III: ρUpper is stable and ρLower =0 (below the critical point); Region IV: ρUpper >ρ Lower >0 and both are stable (bistable); Region V: only the upper solution is stable (the global critical mass was reached, ρ⩾ρc). 023032-3 DE ARRUDA, PETRI, AND MORENO PHYSICAL REVIEW RESEARCH 2, 023032 (2020) there are two possible solutions, ρLower and ρUpper.Thesolution depends on the initial conditions and the threshold, ∗, which, together with the structure, defines a value ρc where the dynamics exhibit a discontinuity. If ρ(t=0) ⩾ρc, the solution is given by ρ=ρUpper (forward phase diagram). On the other hand, if ρ(t=0) <ρ cand ρ(t=0) = 0, then ρ=ρLower (backward phase diagram). The arrows show these solutions as the size of the jump (i.e., the magnitude of the latent heat). Note that the lower solution can exhibit a secondorder phase transition. Next, in Fig. 2(b) we instead represent the corresponding parameter space, which is composed of five distinct regions, as explained in the figure caption. We assumed the most general case, where the lower solution has a transition from the absorbing state to an active state, here at λc. We remark that, for some structures, the lower solution might have a vanishing critical point (e.g., the hyperstar), i.e., λc→0, thus slightly changing this picture. Finally, we have compared both analytic and numeric estimates of the latent heat for this hypergraph structure. The results show that the absolute error between analytical and numerical simulations is of order 10−2–10−3in hypergraphs with N=104(see Ref. [28], Sec. VI.E and Table I), indicating that the first-order approximation is accurate. Thus, generically, the solutions for a social contagion dynamics on hypergraphs can be mathematically expressed as ρ=ρLower if ρLower <ρ c ρUpper if ρLower ⩾ρc,(12) ρ∗=ρUpper if ρUpper ⩾ρc ρLower if ρUpper <ρ c,(13) where ρis obtained if ρ(t=0) <ρ c, and ρ∗if ρ(t= 0) ⩾ρc, where ρcis a global critical mass, i.e., the value of ρat which the discontinuity appears. As before, we note that the lower solution (branch) might also exhibit a secondorder (continuous) phase transition, denoted by λcin Fig. 2, from the absorbing state (ρ=0) to the active state (ρ>0). Furthermore, for a given hypergraph with fixed δand λ∗,the discontinuity points are formally defined as λL c=argλ[ρLower(λ, δ, λ∗,N)=ρc],(14) λU c=argλ[ρUpper(λ, δ, λ∗,N)=ρc],(15) thus, also defining the bistable region, (λU c,λ L c). Although a closed solution for the general case is not possible, Monte Carlo simulations and numerical evaluation of Eq. (2) are reasonable alternatives to characterize our system (see Ref. [28], Sec. VI). Here we focus on a hypergraph with an exponential distribution of cardinalities (i.e., the number of nodes inside a hyperedge), P(|e|)∼μexp (−μ|e|) with the constraint that |e|⩾2. Dynamically, we set λ∗=log2(|e|). This choice is arbitrary, but we choose here the log2(|e|) function because it grows sublinearly. Note that, if a hyperedge cardinality goes to infinity, the average spreading value tends to zero: lim|e|→∞ log2(|e|) |e|=0. The impact of such a function is yet unknown, and we leave this analysis for future work. Figures 3(a) and 3(b) show that the order parameter and the susceptibility follow the patterns expected for a first-order transition, i.e., both are discontinuous. Moreover, the order FIG. 3. Estimation of ρand χusing the QS method in a hypergraph with an exponential distribution of hyperedge cardinalities and N=104. The dynamical parameter are: δ=1.0, λ∗=log2(|ej|), and ∗=0.1,0.2. (a) The susceptibility; (b) the order parameter. We considered two initial conditions for the QS method, ρQS(t=0) = 0.01, darker colors, and ρQS(t=0) =1.00 lighter colors. (c) The distribution of active node estimated using the QS method at λ= 0.086 and ∗=0.2 [the crossing between the two susceptibility curves, in Fig. 3(a)]. parameter is bistable, implying the presence of a hysteresis loop. This phenomenon is opposed to an SIS on a graph. The SIS has a second-order phase transition, characterized by a continuous behavior of the order parameter and a diverging susceptibility in the thermodynamic limit. Complementarily, Fig. 3(c) shows the distribution of active nodes in the upper and lower branches. In the former, we have a bell-shaped distribution, similar to the supercritical regime of an SIS process [7,30]. In the latter, we have a distribution peaked at one, similar to the subcritical regime (absorbing state) of an SIS process [7,30]. We emphasize that Fig. 3(c) displays the distribution of active nodes for the upper (left panel) and lower (right panel) branches and that the complete distribution for a given λin a region where both solutions exist is bimodal. Intuitively, one expects that lower cardinality hyperedges are responsible for the lower branch as they are easier to activate than the higher cardinality ones. Note that, in regular cases, such as the hyperblob (see Ref. [28], Secs. III and IV), 023032-4 SOCIAL CONTAGION MODELS ON HYPERGRAPHS PHYSICAL REVIEW RESEARCH 2, 023032 (2020) the giant pairwise component has Nnodes. However, for the exponentially distributed cardinality scenario, this is not the case. Although P(|e|=2) ⩾P(|e|=k), k=3,4...,N (exponential distribution), the largest connected component, is very small, six nodes in the simulated hypergraph. Therefore, in this scenario, the lower branch is determined by hyperedges with cardinality lower but greater than two. The generality of the reported phenomenological behavior suggests that the observed dynamics are a consequence of group-group interactions. Additional experiments (reported in Ref. [28]) for a hypergraph with a power-law distribution of cardinality further corroborated these results. In all systems, we found similar qualitative behavior for finite networks. Our results are relevant because they provide a theoretical foundation for, and a phenomenological explanation to, seemingly different experimental findings [24–27,31–33]. These works reported critical mass levels needed to change an established equilibrium of 10% in some experiments and 30%–40% in others, in apparent contradiction. The formalism here developed naturally brings forth plausible hypotheses for these observations and show that both ranges are possible. On the one hand, studies based on a single group suggest a threshold between 30%–40%, a situation that can be modeled as a single hyperedge in our formalism. On the other hand, a critical mass of 10% would correspond to a population that is composed of groups of diverse sizes, each one with a (larger) activation threshold. In other words, it is possible to have individual groups exhibiting a threshold ∗between 30%–40%, and at the same time, a global critical mass, ρc,for the whole population of about 10% due to group intersections. A second reason that could explain the experimental findings is even more straightforward: admittedly, the fact that our model shows bistability also enables, for a given λ,two possible solutions for ρcorresponding to the lower and the upper branches. That is, the system might be operating in the region where both solutions are larger than zero and stable. In summary, in this paper, we have developed a framework that allows extending the study of social contagion models when group interactions are relevant. This is achieved by considering hypergraphs as the substrates that capture such many-to-many interactions. First, our work opens the path to deal with new dynamical processes on top of higher-order models and specifically on hypergraphs. Second, we showed that simple dynamical processes could exhibit very rich dynamics, with different transitions, bistability, and hysteresis. Hypergraphs are ubiquitous, and our theory suggests that such a structure allows for the phase diagram reported here. Several findings support the relevance of this methodology. We remark that, depending on the structure, traditional graphprojected models may lead to wrong results. Ultimately, the uncovered phenomenology allows explaining seemingly contradictory experimental findings in which group interactions play a major role. We also mention that many interesting questions arise from our work. For instance, if one assumes that energy is proportional to ρ, our model might display phenomena reminiscent of a Carnot cycle for social contexts, which might help to understand abrupt changes and oscillatory patterns in social behaviors. G.F.A. thanks E. Artiges, H. F. de Arruda, J. P. Rodriguez, L., Gallo and T. Peron for fruitful and inspiring discussions. G.P. acknowledges support from Compagnia San Paolo (ADnD project). Y.M. acknowledges partial support from the Government of Aragon, Spain, through Grant E36-17R (FENOL), and by MINECO and FEDER funds (FIS201787519-P). G.F.A., G.P., and Y.M. acknowledge support from Intesa Sanpaolo Innovation Center. Research carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (Grant 2013/07375-0). The funders had no role in study design, data collection, and analysis, decision to publish, or preparation of the manuscript. [1] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86,3200 (2001). [2] A. Barrat, M. Barthlemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, New York, 2008). [3] M. Newman, Networks: An Introduction (OUP, Oxford, 2010). [4] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. Hwang, Phys. Rep. 424,175 (2006). [5] A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F. Mendes, Phys.Rev.Lett.109,128702 (2012). [6] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys. 87,925 (2015). [7] G. F. de Arruda, F. A. Rodrigues, and Y. Moreno, Phys. Rep. 756,1(2018). [8] S. Chatterjee and R. Durrett, Ann. Probab. 37,2332 (2009). [9] V. Colizza, A. Barrat, M. Barthelemy, A.-J. Valleron, and A. Vespignani, PLoS Med. 4,e13 (2007). [10] D. Balcan, H. Hu, B. Goncalves, P. Bajardi, C. Poletto, J. J. Ramasco, D. Paolotti, N. Perra, M. Tizzoni, W. Van den Broeck et al.,BMC Med. 7,45 (2009). [11] M. Tizzoni, P. Bajardi, C. Poletto, J. J. Ramasco, D. Balcan, B. Gonçalves, N. Perra, V. Colizza, and A. Vespignani, BMC Med. 10,165 (2012). [12] Q. Zhang, K. Sun, M. Chinazzi, A. Pastore y Piontti, N. E. Dean, D. P. Rojas, S. Merler, D. Mistry, P. Poletti, L. Rossi et al.,Proc. Natl. Acad. Sci. USA 114,E4334 (2017). [13] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, andM.A.Porter,J. Complex Netw. 2,203 (2014). [14] G. F. de Arruda, E. Cozzo, T. P. Peixoto, F. A. Rodrigues, and Y. Moreno, Phys. Rev. X 7,011014 (2017). [15] E. Cozzo, G. F. de Arruda, F. A. Rodrigues, and Y. Moreno, MultiplexNetworks:BasicFormalismandStructuralProperties (Springer International Publishing, Cham, 2018). [16] G. F. de Arruda, E. Cozzo, F. A. Rodrigues, and Y. Moreno, New J. Phys. 20,095004 (2018). 023032-5 DE ARRUDA, PETRI, AND MORENO PHYSICAL REVIEW RESEARCH 2, 023032 (2020) [17] A. Aleta and Y. Moreno, Annu. Rev. Condens. Matter Phys. 10, 45 (2019). [18] S. Gómez, A. Díaz-Guilera, J. Gómez-Gardeñes, C. J. PérezVicente, Y. Moreno, and A. Arenas, Phys. Rev. Lett. 110, 028701 (2013). [19] R. Lambiotte, M. Rosvall, and I. Scholtes, Nat. Phys. 15,313 (2019). [20] A. R. Benson, R. Abebe, M. T. Schaub, A. Jadbabaie, and J. Kleinberg, Proc. Natl. Acad. Sci. USA 115,E11221 (2018). [21] Á. Bodó, G. Y. Katona, and P. L. Simon, Bull. Math. Biol. 78, 713 (2016). [22] I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nat. Commun. 10,1(2019). [23] J. A. Caetano, G. Magno, M. Gonçalves, J. Almeida, H. T. Marques-Neto, and V. Almeida, in Proceedings of the 10th ACM Conference on Web Science (ACM, New York, 2019), pp. 27–36. [24] D. Centola, J. Becker, D. Brackbill, and A. Baronchelli, Science 360,1116 (2018). [25] R. M. Kanter, Am.J.Sociol.82,965 (1977). [26] D. Dahlerup, Scand. Polit. Stud. 11,275 (1988). [27] S. Grey, Politics Gender 2,492 (2006). [28] See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevResearch.2.023032 for a complete derivation of our model and additional experiments. [29] M. Fernandez and S. Williams, IEEE Trans. Aerosp. Electron. Syst. 46,803 (2010). [30] S. C. Ferreira, C. Castellano, and R. Pastor-Satorras, Phys. Rev. E86,041125 (2012). [31] J. Xie, S. Sreenivasan, G. Korniss, W. Zhang, C. Lim, and B. K. Szymanski, Phys.Rev.E84,011130 (2011). [32] X. Niu, C. Doyle, G. Korniss, and B. K. Szymanski, Sci. Rep. 7,41750 (2017). [33] A. Baronchelli, R. Soc. Open Sci. 5,172189 (2018). 023032-6