DDynamics of Majority Rule on Hypergraphs
James Noonan ∗ and Renaud Lambiotte † Mathematical Institute, University of Oxford, Oxford, UK (Dated: January 12, 2021)A broad range of dynamical systems involve multi-body interactions, or group interactions, whichmay not be encoded in traditional graphical structures. In this work, we focus on a canonicalexample from opinion dynamics, the Majority Rule, and investigate the possibility to represent andanalyse the system by means of hypergraphs. We explore the formation of consensus and restrictour attention to interaction groups of size 3, in order to simplify our analysis from a combinatorialperspective. We propose different types of hypergraph models, incorporating modular structure ordegree heterogeneity, and recast the dynamics in terms of Fokker-Planck equations, which allowsus to predict the transient dynamics toward consensus. Numerical simulations show a very goodagreement between the stochastic dynamics and theoretical predictions for large population sizes.
Keywords: majority rule, consensus, higher-order, non-linear, networks, group dynamics, hypergraphs
I. INTRODUCTION
Opinion dynamics is concerned with the study of con-sensus formation in populations of interacting individu-als. An important evolution of the field has been to con-sider the impact of the structure of the underlying socialnetwork on the dynamics [1]. In contrast with mean-fieldapproaches, where all the agents are essentially connectedwith each other, network-based approaches assume thatthe agents are located at the nodes of a network and thatthey are sparsely related by direct, binary connectionsthrough edges. A combination of these edges allows forindirect connections through the notion of path. A widerange of models have been proposed so as to capture theways in which social interactions between agents affectopinion formation, through factors like peer pressure andconviction. Here, we will focus on the popular class ofmodels where each node may be in one of two states, rep-resented by the binary variables and . These allow usto model opinions relating to questions with a “yes”/“no”answer, as in a referendum for instance, or “left”/“right”political choices, and find interesting connections withstatistical physics models for spin dynamics.The Voter Model (VM) is a prime example of such bi-nary models [2] and is defined as follows. At discretetimes, an agent is chosen uniformly at random from thepopulation. This agent then adopts the opinion of oneof its randomly-chosen neighbours in the underlying net-work. This update is repeated ad infinitum , or until con-sensus is necessarily reached on a finite connected graph.Importantly, VM is linear and it is dyadic in nature giventhat pairwise interactions alone are sufficient to capturethe dynamics of the system. It is well known that VM issolvable on regular lattice structures in arbitrary spatialdimensions [4]. This is due to the fact that the averagenode state is conserved on degree-regular graphs. VMhas also been shown to be conservative on heterogeneous ∗ [email protected] † [email protected] networks [5], allowing for significant analytical progressto be made in the study of the associated dynamics on awide range of graph topologies.However, real-world dynamics often exhibit nonlin-ear and non-conservative behaviour. An important non-conservative generalisation of VM is the so-called Major-ity Rule (MR) model [6, 7] where, at each update event,we choose a group of G agents from the population, where G ≥
3. These agents form the interaction group, andmay be chosen uniformly at random (as in a mean-fieldscenario), or in a way that is constrained by the under-lying network, as we discuss below. After an interactiongroup is formed, all of its agents simultaneously adoptthe majority opinion in the group. When G is odd, themajority opinion is always well defined. When G is even,and a tie is observed between the opinions, the consensusis either decided randomly, or by introducing a bias forone opinion [8].An important aspect of MR is that interactions takeplace in groups, motivated by the mechanism of peerpressure, which questions the adequacy of networks toencode the interactions between agents. Indeed the dy-namical model is based on group interactions while theunderlying network only encodes pairwise interactions,and there is thus no simple way to choose a group of G agents, e.g. when G = 3, should they form a triangle, orsimply form a path of length 2? This type of question-ing has gained a lot of attention in recent years [3, 9], asit was observed that many systems exhibit multi-bodyor group interactions, such as in neural activity [10–12],robotics [13] or scientific collaborations [14], and that tra-ditional graphical structures are incapable of reflectingthe multi-body nature of such interactions. In order tocircumvent this issue, a natural choice is to adopt a moregeneral class of topological structures known as hyper-graphs where interactions may, as the dynamical system,be multi-body. While there are other topological struc-tures that are capable of encoding group interactions suchas simplicial complexes [15] and bipartite graphs [16], hy-pergraphs constitute a straightforward representation ofmulti-body interactions on networks, and appear to be an a r X i v : . [ phy s i c s . s o c - ph ] J a n ideal candidate for the representation and analysis of MRdynamics [17]. Importantly, as shown in [18], nonlinearmultibody interactions are not expressible as linear com-binations of pairwise interactions between adjacent nodeson a standard graph. MR dynamics on hypergraphs canthus not be reduced to a dynamics on standard networksin general.The main purpose of this article is to investigate MRdynamics on a broad range of hypergraph models. Forthe sake of simplicity, we will restrict our attention to in-teraction groups of size 3, referred to as triangles , even ifmost of our results could be generalised to general groupsize. In Section II, we start by reviewing results aboutthe MR model in the mean-field limit, revealing its non-conservative nature and deriving some of its propertiesthrough a Fokker-Planck equation. In Section III, wethen investigate MR dynamics on the so-called tripar-tite hypergraph, which is a natural generalisation of themean field case where the system is made of 3 types ofnodes. In Section IV, we consider a model of hyper-graphs with community structure, referred to as modularhypergraphs, and then extend our analysis to heteroge-neous hypergraphs in Section V. Section VI concludesour work. II. MEAN-FIELD ANALYSISA. Exact Analysis of the Exit Probability
The selection rule constitutes a key element of the MRmodel. We begin by considering the dynamics in themean-field. The model is defined as follows. Each nodeis endowed with a binary state variable, denoted by and . At each time step, 3 agents are chosen uniformlyat random and the Majority Rule is applied. This ran-dom selection can be formulated conveniently in the lan-guage of hypergraphs. Let us consider a fully connectedhypergraph of N agents. The hypergraph structure H consists of the node set V ( H ) = { , . . . , N } and, as werestrict ourselves to three-body interactions, the set ofall possible triangles given by T ( H ) = {{ i, j, k } : i, j, k ∈ V ( H ) , i (cid:54) = j (cid:54) = k } . The resulting object is a first gener-alisation of fully connected graphs to hypergraphs. Theselection of a group of 3 nodes is now defined as therandom selection of one hyperedge in the set T ( H ) ofavailable hyperedges.When studying the Majority Rule, as well as otheropinion dynamics models with discrete states, an impor-tant quantity is the exit probability defined as follows.Suppose that the system is initialised with n < N agentsin the state and N − n agents in the state. The exitprobability E n is the probability that the system reachesconsensus with all agents in the state given that n agents are initiated in the state. Krapivsky and Red-ner [7] adopted a combinatorial approach in deriving an FIG. 1: Exit probability for MR dynamics in the mean-field. exact expression for E n in the mean-field: E n = 12 N − n − (cid:88) j =1 Γ ( N − j ) Γ ( N − j − . (1)Figure 1 illustrates mean field exit probabilities using N = 20 ,
50 and 100 agents. The initial fraction of agents n/N in the state is plotted on the x -axis while the exitprobability E n is plotted on the y -axis. The exit proba-bility is sigmoidal in nature, rapidly changing in a nar-row interval centred at the pivotal value n/N = 0 .
5. As N → ∞ the exit probability converges to 0 for n/N < / n/N = 0 .
5. The complex functional form of E n is a man-ifestation of the nonlinear and non-conservative natureof MR interactions. B. The Fokker-Planck Approach
We now go beyond the workings presented in [7] and ap-proach the mean-field analysis from a different perspec-tive. In this section, we study the exit probability andthe dynamical approach towards consensus in the asymp-totic limit N → ∞ . Let ρ ( t ) denote the density of ’s inthe population at time t ; that is, the fraction of agentsin the state at time t . We shall write ρ ( t ) as ρ forconvenience, though the dependence in t is implicit. Fur-thermore, let δρ = N − denote the incremental changein ρ following an update event. We define the quantities ρ ± = ρ ± δρ to reflect the density of states in the popula-tion after an interaction takes place. Let R ( ρ ) and L ( ρ )denote the raising and lowering operators [5] that givethe transition probabilities for the update events ρ → ρ + and ρ → ρ − respectively. Finally, let p ( ρ, t ) be the prob-ability that the density of ’s in the population is ρ attime t . The probability density evolves according to thefollowing master equation over the incremental time pe-riod δt : p ( ρ, t + δt ) = R (cid:0) ρ − (cid:1) p (cid:0) ρ − , t (cid:1) + L (cid:0) ρ + (cid:1) p (cid:0) ρ + , t (cid:1) + [1 − R − L ] p ( ρ, t ) , (2)where R and L are assumed to denote R ( ρ ) and L ( ρ )respectively. A conventional choice of the quantity δt is N − [5, 19]. We proceed by Taylor expanding the lefthand side of equation (2) to first order in δt and the righthand side to second order in δp . This ultimately yieldsthe Fokker-Planck equation: ∂∂t p ( ρ, t ) = − δρδt (cid:20) ∂∂ρ ( R − L ) p ( ρ, t ) (cid:21) + ( δρ ) δt ∂ ∂ρ [( R + L ) p ( ρ, t )] . (3)In order to proceed, it is instructive to highlightthe connection between equation (3) and its associatedstochastic differential equation. Suppose that X = { X t ( ω ) } t ∈ I,ω ∈ Ω is a stochastic process on the samplespace Ω over the time interval I = [0 , T ] for some T > dX t = v ( X t , t ) dt + σ ( X t , t ) dW t (4)where W t is a standard Brownian motion. v ( X t , t ) and σ ( X t , t ) are referred to as the drift and diffusion coef-ficients respectively. We shall also refer to v ( X t , t ) asthe drift velocity. The Fokker-Planck equation associ-ated with the probability density p ( x, t ) of the randomvariable X t is given by ∂∂t p ( x, t ) = − ∂∂x [ v ( x, t ) p ( x, t )]+ ∂ ∂x [ D ( x, t ) p ( x, t )] (5)where D ( X t , t ) = σ ( X t , t ) /
2. Comparing equations (3)and (5) we observe that the density ρ may be describedas a stochastic process with drift velocity v ( ρ, t ) and dif-fusion coefficient D ( ρ, t ) defined as follows: v ( ρ, t ) = δρδt ( R − L ) = R − L, (6) D ( ρ, t ) = ( δρ ) δt ( R + L ) = 12 N ( R + L ) . (7)If N is sufficiently large so that ( N − /N ≈
1, equa-tions (6) and (7) simplify to give v ( ρ, t ) = 3 ρ (1 − ρ ) (2 ρ − , (8) D ( ρ, t ) = 32 N ρ (1 − ρ ) . (9)Equations (8) and (9) illustrate that the drift veloc-ity v ( ρ, t ) is O (1) whereas the diffusion term D ( ρ, t )is O (1 /N ). The drift velocities thus dominate the dy-namics of the system for sufficiently large N , as diffusive (a) ρ (0) = 0 . (b) ρ (0) = 0 . FIG. 2: Temporal density profiles for mean-field MR dynam-ics with N = 10 . contributions vanish in the asymptotic limit. This ob-servation will underpin much of the analysis presentedthroughout the course of this paper.In the absence of diffusive contributions, we can useequation (4) to write v = dρ/dt , and thus integrate equa-tion (8) directly to give the following result: ρ ( t ) = (cid:16) − (cid:112) − / (4 + κe t ) (cid:17) if ρ (0) < . , (cid:16) (cid:112) − / (4 + κe t ) (cid:17) if ρ (0) > . , where κ = (2 ρ (0) − / ( ρ (0) (1 − ρ (0))). This indi-cates that the system will rapidly reach consensus withall agents in the state if ρ (0) < .
5. Conversely, if ρ (0) > . state. As a validation of our results, wesimulate the mean-field MR dynamics in a population ofsize N = 10 . Figure 2 shows two representative resultswith ρ (0) > .
5. The density profiles are analogous for ρ (0) < . ρ = 0 .
5. The an-alytical trajectory is plotted in red, while the stochastictrajectory is plotted in blue. The time scale on the x -axisis measured in units of Monte Carlo steps per node, sothat δt = N − . The trajectories are in excellent agree-ment when the initial value ρ (0) is sufficiently far from0 .
5, and the system rapidly approaches consensus alongthe predicted trajectories. When ρ (0) = 0 . ± (cid:15) where0 < (cid:15) (cid:28)
1, however, diffusive fluctuations at early timescan lead to a temporal shift in the density profile alongthe x -axis, and deviations between the predictions andthe stochastic simulations can be observed, as illustratedin Figure 2a. Such diffusive effects become negligible as N → ∞ . III. THE TRIPARTITE HYPERGRAPH
In this Section we move beyond the mean-field andconsider hypergraphs with more sophisticated topologies.A natural starting point in this regard is the tripar-tite hypergraph, denoted by H . The tripartite hyper-graph consists of three distinct groups of nodes, whichwe shall refer to as G a , G b , and G c . We assume for thesake of simplicity that each group consists of N agents.Let V a ( H ) = { a , . . . , N a } , V b ( H ) = { b , . . . , N b } , and V c ( H ) = { c , . . . , N c } denote the sets of nodes in groups G a , G b , and G c respectively. Let us define the set oftriangles on H as T ( H ) = {{ i, j, k } : i ∈ V a ( H ) , j ∈ V b ( H ) , k ∈ V c ( H ) } . In other words, we consider all pos-sible triplets with one agent in each group, hence gen-eralising the notion of the complete bipartite graph. Asbefore, at each step, a group of 3 nodes is chosen byselecting one hyperedge in T ( H ) at random.Let us now adapt our treatment of the mean-field viaa Fokker-Planck equation to this setting. Let ρ a ( t ) de-note the density of nodes in state in G a at time t ,and let δρ a = N − denote the change in ρ a ( t ) result-ing from the changing of opinion of a single agent in G a .We then define ρ ± a = ρ a ± δρ a . Analogous expressionshold for agents in G b and G c where ρ b and ρ c denotethe associated densities of ’s respectively. We also de-fine the raising and lowering operators associated with G a , denoted by R a ( ρ a , ρ b , ρ c ) and L a ( ρ a , ρ b , ρ c ) respec-tively (written as R a and L a for convenience). R a is thetransition probability associated with the update event ρ a → ρ + a , whereas L a is the transition probability as-sociated with the update event ρ a → ρ − a . Once again,analogous expressions exist for the raising and loweringoperators associated with groups G b and G c . R a and L a are given as follows: R a = (1 − ρ a ) ρ b ρ c , (10) L a = ρ a (1 − ρ b ) (1 − ρ c ) . (11)Equation (10) follows from the fact that ρ a increases to ρ + a if we choose a node in state from G a as well as twonodes in state from G b and G c . Similarly, equation(11) is derived from the fact that ρ a decreases to ρ − a ifwe choose a node in state from G a as well as two nodesin state from G b and G c . The same logic applies forthe calculation of the transition probabilities for groups G b and G c .Using these transition probabilities we can deduce theprobabilistic master equation governing MR dynamicson the tripartite hypergraph, where p ( ρ a , ρ b , ρ c , t ) is theprobability of the system having densities ρ a , ρ b and ρ c at time t : p ( t + δt ) = R a (cid:0) ρ − a (cid:1) p (cid:0) ρ − a (cid:1) + L a (cid:0) ρ + a (cid:1) p (cid:0) ρ + a (cid:1) + R b (cid:0) ρ − b (cid:1) p (cid:0) ρ − b (cid:1) + L a (cid:0) ρ + b (cid:1) p (cid:0) ρ + b (cid:1) + R c (cid:0) ρ − c (cid:1) p (cid:0) ρ − c (cid:1) + L a (cid:0) ρ + c (cid:1) p (cid:0) ρ + c (cid:1) + [1 − R a − R b − R c − L a − L b − L c ] p. (12)We proceed by Taylor expanding the left hand side ofequation (12) to first order in δt and the right hand sideto second order in δρ a , δρ b and δρ c . We take δt to beequal to the reciprocal of the total number of nodes inthe hypergraph, in line with the convention introduced inSection II. This ultimately yields the three-dimensional Fokker-Planck equation: ∂p∂t = − δρ a δt ∂∂ρ a [( R a − L a ) p ] + ( δρ a ) δt ∂ ∂ρ a [( R a + L a ) p ] − δρ b δt ∂∂ρ b [( R b − L b ) p ] + ( δρ b ) δt ∂ ∂ρ b [( R b + L b ) p ] − δρ c δt ∂∂ρ c [( R c − L c ) p ] + ( δρ c ) δt ∂ ∂ρ c [( R c + L c ) p ] . (13)As in the mean field case, it is instructive to high-light the connection between equation (13) and its cor-responding stochastic differential equation. Consider an n -dimensional stochastic process X = { X t ( ω ) } t ∈ I,ω ∈ Ω on the time interval I = [0 , T ] and the sample space Ω,where n ≥ T >
0. Suppose X t ( ω ) = X t hasstochastic differential d X t = v ( X t , t ) dt + σ ( X t , t ) d W t , (14)where v ( X t , t ) is an n -dimensional random vector, σ ( X t , t ) is an n × m dimensional matrix and W t is astandard m -dimensional Weiner process where m ≥ X t by p ( x , t ), then p ( x , t ) obeys thefollowing Fokker-Planck equation: ∂p ( x , t ) ∂t = − n (cid:88) i =1 ∂∂x i [ v i ( x , t ) p ( x , t )]+ n (cid:88) i =1 n (cid:88) j =1 ∂ ∂x i ∂x j [ D ij ( x , t ) p ( x , t )] , (15)where v = [ v , ..., v n ] is the vector of drift velocities and D = σσ T is the diffusion tensor. Let us denote thedrift velocity associated with the density of ’s in G a by v a and the corresponding diffusion term by D a .By comparing equations (13) and (15) we can deducethe following expressions: v a = 3 [(1 − ρ a ) ρ b ρ c − ρ a (1 − ρ b ) (1 − ρ c )] , (16) D a = 32 N [(1 − ρ a ) ρ b ρ c + ρ a (1 − ρ b ) (1 − ρ c )] . (17)Note that in this instance the diffusion tensor is diago-nal, where D a , D b , and D c are equal to D , D , and D respectively. This is because we are considering in-teraction groups of size 3, meaning that only one of thequantities ρ a , ρ b , or ρ c is varied at each update event.This is inherent in the structure of equation (12). Anal-ogous expressions exist for the drift and diffusion termsfor ρ b and ρ c . As in the mean field case, we note thatthe drift velocity is O (1) whereas the diffusive term is O (1 /N ). This implies that the dynamics are dominatedby the drift velocities for large N , given that the diffusivecontributions become negligible as N → ∞ . Using thisfact, we posit that the dynamics of the system may bewell approximated by discarding the diffusive terms and (a) ( ρ a , ρ b ) (b) ( ρ b , ρ c ) (c) ( ρ a , ρ c ) FIG. 3: Simulation results on the tripartite hypergraph with N = 10 . The initial conditions are given by ρ = (0 . , . , . state. formulating a dynamical system using the drift velocitiesalone. This approach follows from the fact that in theabsence of diffusive contributions, equation (14) reducesto a system of ordinary differential equations: dρ a dt = v a = 3 [(1 − ρ a ) ρ b ρ c − ρ a (1 − ρ b ) (1 − ρ c )] ,dρ b dt = v b = 3 [(1 − ρ b ) ρ a ρ c − ρ b (1 − ρ a ) (1 − ρ c )] ,dρ c dt = v c = 3 [(1 − ρ c ) ρ a ρ b − ρ c (1 − ρ a ) (1 − ρ b )] . This system may be analysed via linear stabilityanalysis. The set of fixed points is given by { (0 , , , (1 / , / , / , (1 , , } . When evaluated atthe points (0 , ,
0) and (1 , ,
1) the eigenvalues of the Ja-cobian are − / , / , /
2) the eigenval-ues of the Jacobian are − / / , / , /
2) is asaddle point.We investigate the performance of our model byconducting simulations on a tripartite hypergraph of N = 10 nodes. Let us denote the initial conditions( ρ a (0) , ρ b (0) , ρ c (0)) by ρ . Figure 3 illustrates a sam-ple set of simulation results with ρ = (0 . , . , . ρ a , ρ b ) , ( ρ b , ρ c ) and( ρ a , ρ c ) planes. These projections are plotted from leftto right in Figure 3. The deterministic and stochastictrajectories are evidently in excellent agreement. Thisserves to validate our conjecture that the drift velocitiesdominate the dynamics of the system for large values of N .Our numerical simulations thus far have only focusedon a specific trajectory of the dynamics resulting from an arbitrary choice of initial condition. We now extendour analysis by approximating the exit probability E ( ρ )of the system for large N ; that is, the probability thatall of the agents reach consensus in the state if theinitial conditions are given by ρ . In order to proceed, weconsider the unstable fixed point at (1 / , / , / E u is spanned by the vector [1 , , E s spanned bythe remaining two eigenvectors. It is straightforward todeduce that the equation of the plane corresponding tothe stable linear subspace is given by ρ a + ρ b + ρ c = 3 / . (18)By symmetry we posit that in the limit of large N , theexit probability is given by the following piecewise func-tion: E ( ρ ) = (cid:40) ρ a (0) + ρ b (0) + ρ c (0) < / , ρ a (0) + ρ b (0) + ρ c (0) > / . (19)In order to confirm this prevision numerically, we vary ρ a (0) and ρ b (0) while keeping ρ c (0) fixed. Figure 4gives the simulation results for ρ c (0) = 0 . , . , and0 .
75, where N = 900. These plots were generated byinitialising the stochastic dynamics at each point on the( ρ a (0) , ρ b (0)) grid with uniform mesh size 0 .
01. Purpleindicates consensus in the state whereas yellow indi-cates consensus in the state. The red dashed line marksthe boundary predicted by equation (18) at t = 0. Onlyone simulation of the stochastic dynamics was conductedper coordinate, giving rise to diffusive effects at the in-terface separating the two domains. Nevertheless, thedeterministic exit probability is in excellent agreementwith the stochastic simulations of the system. (a) ρ c (0) = 0 . (b) ρ c (0) = 0 . (c) ρ c (0) = 0 . FIG. 4: Exit probability analysis on the tripartite hypergraph with N = 900. IV. THE MODULAR HYPERGRAPHA. The Symmetric Case
In this Section we consider MR dynamics on hyper-graphs with community structure, otherwise known asmodular hypergraphs [20, 21]. In our setting, the popu-lation is partitioned into a number of non-overlappingsets known as modules or communities [22, 23]. Themodular hypergraph H consists of two sets of N nodes,referred to as communities A and B respectively. Let V A ( H ) = { A , . . . , N A } and V B ( H ) = { B , . . . , N B } de-note the node sets in communities A and B respectively.The set of triangles on H is defined as T ( H ) = {{ i, j, k } : i, j, k ∈ V A ( H ) ∪ V B ( H ) , i (cid:54) = j (cid:54) = k } . In order to dis-tribute the triangles over the hypergraph, we define theparameter p ab which gives the probability of forming atriangle with a agents in A and b agents in B . It followsthat p ab ∈ { p , p , p , p } . These can be interpretedas hyperedge selection probabilities. Here, the hyper-graph is assumed to be symmetric, meaning p = p and p = p . The total number of triangles T = | T ( H ) | is given by T = 2 N (2 N −
1) (2 N − ≈ N . In order to vary the degree to which the two compo-nents of the modular hypergraph interact, we introduce aparametric dependence inspired by the study of VM dy-namics on the two-clique graph [5]. In graphical terms,communities A and B may be thought of as two fully-connected cliques of N nodes. We then define an inter-connectivity parameter C ∈ [0 , N ] such that each nodein clique A is, on average, connected to C nodes in clique B and vice versa. It follows that the average degree ofa node in the graphical representation is N + C , as eachclique is complete. Using this graphical representation itis straightforward to calculate hyperedge selection prob- abilities: p = p = 12 (cid:18) NN + C (cid:19) , (20) p = p = N C ( N + C ) + 12 (cid:18) CN + C (cid:19) . (21)These probabilities satisfy normalisation. When C = 0it follows that p = p = 1 / p = p = 0,hence the two modules evolve independently. Conversely,when C = N it follows that p = p = 1 / p = p = 3 /
8, corresponding to a mean-field of 2 N agents. We now define N ab as the number of hyperedgeswith a nodes in A and b nodes in B . Using the factthat p ab = N ab /T the hyperedge distribution may bedetermined as a function of C : N = N = 2 N N + C ) , (22) N = N = 4 N C N + C ) (cid:18) N + C (cid:19) . (23)Note that hyperedges may have multiplicity greater than1. This is due to the fact that the total number of hy-peredges T is fixed, regardless of the value of C . Forexample, when C = 0 we find that N = N = 2 N / N = N = 0. However, there are only N / N nodes, im-plying that all hyperedges have multiplicity 4 in this in-stance. This combinatorial consequence is immaterial tothe dynamics of the system.Once again, each agent in the population is assumed tooccupy one of two states, labelled and . Simulationsare conducted using the hyperedge selection probabili-ties given in equations (20) and (21). This is equiva-lent to choosing hyperedges uniformly at random fromthe hyperdegree distribution given by equations (22) and(23). Let us denote the density of nodes in state in A and B at time t by ρ A ( t ) and ρ B ( t ) respectively. Wewrite ρ A ( t ) as ρ A and ρ B ( t ) as ρ B for convenience. Let ρ ± A = ρ A ± δρ A and ρ ± B = ρ B ± δρ B where δρ A = N − and δρ B = N − . Let R A ( ρ A , ρ B ) and R B ( ρ A , ρ B ) bethe raising operators that give the transition probabili-ties from ρ A and ρ B to ρ + A and ρ + B respectively. Similarly,let L A ( ρ A , ρ B ) and L B ( ρ A , ρ B ) be the lowering opera-tors that give the transition probabilities from ρ A and ρ B to ρ − A and ρ − B respectively. R A and L A are given by R A = 32 (cid:18) NN + C (cid:19) ρ A (1 − ρ A )+ C (2 N + C )( N + C ) ρ A ρ B (1 − ρ A )+ C (2 N + C )2 ( N + C ) ρ B (1 − ρ A ) , (24) L A = 32 (cid:18) NN + C (cid:19) ρ A (1 − ρ A ) + C (2 N + C )( N + C ) ρ A (1 − ρ A ) (1 − ρ B )+ C (2 N + C )2 ( N + C ) ρ A (1 − ρ B ) . (25)Analogous expression exist for the operators R B and L B . Let p ( ρ A , ρ B , t ) be the probability that communities A and B have densities ρ A and ρ B respectively at time t . Its associated master equation is given by p ( t + δt ) = R A (cid:0) ρ − A (cid:1) p (cid:0) ρ − A (cid:1) + L A (cid:0) ρ + A (cid:1) p (cid:0) ρ + A (cid:1) + R B (cid:0) ρ − B (cid:1) p (cid:0) ρ − B (cid:1) + L B (cid:0) ρ + B (cid:1) p (cid:0) ρ + B (cid:1) + [1 − R A − L A − R B − L B ] p, (26)where δt = (2 N ) − . Taylor expanding the left hand sideof equation (26) to first order in δt and the right handside to second order in δρ A and δρ B yields the followingFokker-Planck equation: ∂p∂t = − δρ A δt ∂∂ρ A [( R A − L A ) p ] − δρ B δt ∂∂ρ B [( R B − L B ) p ]+ ( δρ A ) δt ∂ ∂ρ A [( R A + L A ) p ]+ ( δρ B ) δt ∂ ∂ρ B [( R B + L B ) p ] . (27)The drift velocities v A and v B of the clique densities FIG. 5: Phase diagram for the symmetric modular hyper-graph with N = 2500. The initial densities of (cid:48) s in A and B were taken to be 1 and 0 respectively, with a differentconsensus in each community. The dashed line indicates thetheoretical transition value C t . Data points were averagedover 50 simulations and 500 N update events were conductedper simulation to ensure distributional stationarity. Simula-tions confirm a transition between the coexistence of differentopinions in different communities to global consensus when C is increased. Note that global consensus is, for any value of C , an absorbing state. ρ A and ρ B may be identified using equation (27): v A = 3 (cid:18) NN + C (cid:19) ρ A (1 − ρ A ) (2 ρ A − C (2 N + C )( N + C ) ρ A (1 − ρ A ) (2 ρ B − C (2 N + C )( N + C ) (cid:16) ρ B (1 − ρ A ) − ρ A (1 − ρ B ) (cid:17) , (28) v B = 3 (cid:18) NN + C (cid:19) ρ B (1 − ρ B ) (2 ρ B − C (2 N + C )( N + C ) ρ B (1 − ρ B ) (2 ρ A − C (2 N + C )( N + C ) (cid:16) ρ A (1 − ρ B ) − ρ B (1 − ρ A ) (cid:17) . (29)Equations (28) and (29) imply that the drift velocitiesare O (1) for all C . However, as C tends to 0 the equa-tions decouple, implying that the two communities evolvealmost independently of one another. When C = 0 thedrift velocities decouple completely. Using equation (27)it is straightforward to show that the corresponding dif-fusion terms D A and D B are O (1 /N ) for all C . Thisimplies that the drift velocities dominate the dynamicsof the system for relatively high values of C , allowing fordiffusive contributions to be discarded. We approximatethe stochastic dynamics using a deterministic model bywriting dρ A /dt = v A and dρ B /dt = v B where v A and v B are given in equations (28) and (29) respectively. Theassociated linear stability analysis may be found in Ap-pendix VII. The fixed points at (0 ,
0) and (1 ,
1) are foundto be stable for all C . On the other hand, (1 / , /
2) isa saddle point for
C > C ∗ and a source for C < C ∗ where C ∗ ≈ . N . When C < C ∗ , two additionalfixed points are observed to appear at the coordinates (cid:0) ρ ± A , − ρ ± A (cid:1) where ρ ± A = 12 (cid:32) ± (cid:115) CC − N N + C (cid:33) . (30)Furthermore, they are found to be stable for C < C t andunstable otherwise, where C t ≈ . N . C t denotes the transition value, as it predicts the occurrence of a transi-tion in the hypergraph resulting in metastable state for-mation. Similar behaviour was observed by Lambiotteet. al [24] in their study of MR dynamics on modularnetworks. As C tends to 0 the metastable state coordi-nates converge to (0 ,
1) and (1 , C < C t , uni-form consensus will ultimately be reached. This is due tothe fact that diffusion governs the evolution of the systemonce these metastable states are reached.Figure 5 plots the phase diagram for the symmetricmodular hypergraph. The connectivity parameter is var-ied along the horizontal axis, and the absolute differencebetween the density of states at stationarity is plotted onthe vertical axis. Asymmetric initial conditions were as-sumed, with ρ A (0) = 1 and ρ B (0) = 0. The theoreticaltransition value C t is indicated by the dashed line, andis in very good agreement with the simulation results. B. The Asymmetric Case
We extend our analysis by considering the scenario inwhich p (cid:54) = p and p (cid:54) = p in general, hence leadingto an asymmetry in the hypergraph structure. The pa-rameter C AB ∈ [0 , N ] dictates the distribution of 2 N / N and N categories, whereasthe parameter C BA ∈ [0 , N ] distributes 2 N / N and N categories. The hyper-edge selection parameters are given by p = 12 (cid:18) NN + C AB (cid:19) , (31) p = N C AB ( N + C AB ) + 12 (cid:18) C AB N + C AB (cid:19) , (32) p = N C BA ( N + C BA ) + 12 (cid:18) C BA N + C BA (cid:19) , (33) p = 12 (cid:18) NN + C BA (cid:19) . (34)Using these probabilities, the hyperedge distribution maybe determined as a function of the interconnectivity pa- rameters: N = 2 N N + C AB ) , (35) N = 4 N C AB N + C AB ) (cid:18) N + C AB (cid:19) , (36) N = 4 N C BA N + C BA ) (cid:18) N + C BA (cid:19) , (37) N = 2 N N + C BA ) . (38)Simulations are conducted using the hyperedge selectionprobabilities given in equations (31)-(34), which is equiv-alent to sampling uniformly at random from the hyper-edge distribution given in equations (35)-(38). The driftvelocities v A and v B may be calculated in the usual way,and are found to be O (1) in all parameter regimes (ex-plicit expressions may be found in Appendix VIII). Thediffusion terms D A and D B are found to be O (1 /N ) forall C , implying that the drift velocities dominate the dy-namics of the system for large N .In this case the deterministic model consists of thesystem dρ A /dt = v A and dρ B /dt = v B , with v A and v B given in equations (58) and (59) respectively. Thefixed points (0 ,
0) and (1 ,
1) are asymptotically stable.Linear stability analysis at the fixed point (1 / , /
2) isconducted numerically. Visual inspection of the vectorfields associated with this system indicates that in cer-tain parameter regimes, stable fixed points emerge in thevicinity of the points (0 ,
1) and (1 ,
0) in the phase plane.In the symmetric case, these fixed points were found tolie on the line ρ B = 1 − ρ A . However, in the asymmetriccase this is no longer true. As C AB and C BA are in-creased from 0 these points are observed to move awayfrom the line ρ B = 1 − ρ A . Despite this, the absorbingstates are positioned close to this line, which allows usto determine their domain of existence heuristically. Thedrift velocities dρ A /dt and dρ B /dt are denoted by thefunctions f ( ρ A , ρ B ) and g ( ρ A , ρ B ) respectively.The fixed point coordinates are approximated by solv-ing the following system: f ( ρ A , − ρ A ) = 0 , (39) g ( ρ A , − ρ A ) = 0 . (40)Solving equation (39) yields two solutions for ρ A , referredto as ρ ± f . Similarly, solving equation (40) yields two solu-tions for ρ A , referred to as ρ ± g . Explicit expressions for ρ ± f and ρ ± g are given in Appendix VIII. When C AB = C BA , ρ ± f and ρ ± g reduce to the expressions for ρ ± A given byequation (30), where a correspondence between signs isunderstood. When C AB = C BA = 0, ρ + f = ρ + g = 1 and ρ − f = ρ − g = 0. By continuity, ρ + f and ρ + g are expectedto be approximately equal for C AB = (cid:15) and C BA = (cid:15) where (cid:15) and (cid:15) are small positive perturbation param-eters (and similarly for ρ − f and ρ − g ). However, once the Stable Unstable/Nonexistent
FIG. 6: Metastable state bifurcation plot for the asymmetricmodular hypergraph. approximation ρ ± f ≈ ρ ± g becomes invalid it follows thatthe system of equations (39) and (40) no longer has anyfixed points. We proceed by defining a distance thresh-old ε such that if (cid:107) (cid:16) ρ ± f , − ρ ± f (cid:17) − (cid:0) ρ ± g , − ρ ± g (cid:1) (cid:107) > ε then ρ ± f and ρ ± g are no longer deemed to be solutions toequations (39) and (40). Distance is calculated in theEuclidean norm.Using this threshold, a bifurcation diagram may beconstructed numerically. For a given set of connectivityparameters, metastable state existence is inferred if thedistance threshold criterion is satisfied, and if the eigen-values of the Jacobian when evaluated at (cid:16) ρ ± f , − ρ ± f (cid:17) and (cid:0) ρ ± g , − ρ ± g (cid:1) have negative real parts. A numer-ical exploration of the parameter space indicates that ε = 0 .
115 is a sensible value for the threshold parameter.The resulting bifurcation plot for the metastable states isgiven in Figure 6. The plot is symmetric about the line C AB = C BA . The red marker is located at ( C t /N, C t /N )where C t = 0 . N is the transition value on the sym-metric modular hypergraph. The marker lies on the bi-furcation boundary, as expected. Testing the nature ofparameter coordinates close to the boundary indicatesthat the region of stability in the bifurcation plane isgenerally a very good predictor of the true region of sta-bility in parameter space.Figure 7 provides simulation results on the asymmetricmodular hypergraph with ( C AB , C BA ) = (0 . N, . N ).Figure 6 confirms that metastable states do not emergein this instance. The red deterministic trajectory is anexcellent predictor of the blue stochastic path. This anal-ysis indicates that the exit probability is governed by thegeometry of the deterministic phase plane in the asymp-totic limit as N → ∞ for parameter regimes in whichmetastable states do not emerge.Diffusion becomes the dominating factor in dictatinghow the system reaches consensus for parameter coordi- FIG. 7: Simulations on the asymmetric modular hypergraphwith N = 10 , ρ = (0 , .
9) and ( C AB , C BA ) = (0 . N, . N ). nates in the stable region of Figure 6, due to metastablestate formation. These metastable states tend to thepoints (1 ,
0) and (0 ,
1) as C AB and C BA tend to zero, il-lustrating that local consensus is reached independentlyin each community. V. GENERALISED HETEROGENEOUSHYPERGRAPHSA. Heterogeneous Mean Field Analysis
In this Section, we consider the scenario where the hy-peredge distribution of a hypergraph H is given by aprescribed degree distribution, which we model with aheterogeneous mean field approximation [25], as is of-ten done in network science. Consider a system com-posed of N nodes with a given degree distribution. Wedefine the associated hypergraph H to consist of thenode set V ( H ) = { , . . . , N } and the set of triangles T ( H ) = {{ i, j, k } : i, j, k ∈ V ( H ) } . As noted byNeuhauser et al. [18], the structure of the hypergraphcan be encoded in the adjacency tensor A ∈ R N × N × N with entries A ijk = (cid:40) { i, j, k } ∈ T ( H ) , . (41)The adjacency tensor is symmetric in all dimensions. Inorder to prevent the same node from being chosen multi-ple times in the formation of a triangle, we assume that A ijk = 0 if any two of the indices i, j, k are the same.As usual, each node is assumed to occupy of two states,denoted by and . The degree distribution n k is givenby n k = N k N , (42)0where N k is the number of nodes of degree k . The mo-ments of the degree distribution are given by µ m = (cid:88) k k m n k . (43)Note that µ corresponds to the average node degree (cid:104) k (cid:105) .Let us denote the entire state of the system by η anddefine η ( x ) as the state of node x , which can take thevalues 0 and 1 (representative of the and states re-spectively). During each update event, exactly one nodecan change state. We represent the state of the systemafter changing the state of node x by η x , where η x = (cid:40) η ( y ) if y (cid:54) = x, − η ( x ) if y = x. (44)We first consider a fully-connected hypergraph. Thatis, any three nodes can be chosen to form an interactingtriangle. Suppose we choose a node uniformly at random,which we shall refer to as node x . The probability ofdoing so is 1 /N . To form a triangle, we can choose anytwo other nodes, of which there are (cid:0) N − (cid:1) ways of doingso. Therefore the transition probability at node x can becalculated as follows: P [ η → η x ] = 1 N (cid:0) N − (cid:1) (cid:88) yz A xyz ξ ( x, y, z ) , (45)where the indices y and z sum over all of the N − ξ ( x, y, z ) is anindicator-type function given as follows: ξ ( x, y, z ) = η ( y ) η ( z ) (1 − η ( x ))+ (1 − η ( y )) (1 − η ( z )) η ( x ) . (46)Note that ξ ( x, y, z ) is 1 if node x is in the minority opin-ion in the triangle { x, y, z } and 0 otherwise. The presenceof the adjacency tensor entries A xyz as multiplying fac-tors in the summation prevent the same node from beingchosen twice during triangle formation. The first termon the right hand side of equation (46) ensures that node x flips from state to state if x is in state and nodes y and z are in state . The second term on the righthand side ensures that node x flips from state to state if x is in state and nodes y and z are in state . Let ρ k ( t ) denote the density of ’s on nodes of fixed degree k at time t . We shall write ρ k ( t ) as ρ k for convenience.Furthermore let ρ ± k = ρ k ± δρ k where δρ k = N − k . Let R k [ { ρ k } ] and L k [ { ρ k } ] denote the raising and loweringoperators associated with nodes of degree k . They arecalculated as follows: R k [ { ρ k } ] = P (cid:2) ρ k → ρ + k (cid:3) = 1 N (cid:0) N − (cid:1) (cid:88) yz (cid:88) x (cid:48) A xyz η ( y ) η ( z ) (1 − η ( x )) , (47) L k [ { ρ k } ] = P (cid:2) ρ k → ρ − k (cid:3) = 1 N (cid:0) N − (cid:1) (cid:88) yz (cid:88) x (cid:48) A xyz (1 − η ( y )) (1 − η ( z )) η ( x ) , (48) where the primed summation indicates that the summa-tion is restricted to nodes x of degree k . We now appealto the heterogeneous mean field approximation [25] in as-suming that node states are independent and nodes of thesame degree behave similarly. We proceed by replacingthe entries of the adjacency tensor as follows: A xyz → k x k y k z N µ . (49)Equation (49) gives the interaction probability betweenthree nodes of degrees k x , k y and k z , where the scalingconstants are chosen for the purpose of normalization[25]. This results in triangles being clustered aroundnodes of higher degrees, which is more realistic than auniform random distribution of triangles. Finally, let usdefine the degree-weighted moments ω m [5] of the degreedistribution: ω m = 1 N µ m (cid:88) x k mx η ( x ) = 1 µ m (cid:88) k k m n k ρ k . (50)Using equations (49) and (50), equations (47) and (48)reduce to the following expressions for large N : R k [ { ρ k } ] ≈ kn k ω (1 − ρ k ) N , (51) L k [ { ρ k } ] ≈ kn k ρ k (1 − ω ) N . (52)Let p ( { ρ k } , t ) denote the probability that the density of ’s on nodes of degree k is ρ k at time t . We denote p ( { ρ k } , t ) by p for ease of notation, and similarly for R k and L k . The probability density obeys the followingmaster equation: p ( t + δt ) = (cid:88) k R k (cid:0) ρ − k (cid:1) p (cid:0) ρ − k , t (cid:1) + (cid:88) k L k (cid:0) ρ + k (cid:1) p (cid:0) ρ + k , t (cid:1) + (cid:34) − (cid:88) k ( R k + L k ) (cid:35) p. (53)Taylor expanding equation (53) yields the followingFokker-Planck equation: ∂p∂t = − (cid:88) k δρ k δt ∂∂ρ k (( R k − L k ) p )+ (cid:88) k ( δρ k ) δt ∂ ∂ρ k (( R k + L k ) p ) . (54)Using equation (54) the drift velocities { v k } and diffusionterms { D k } of the densities { ρ k } may be identified: v k = δρ k δt ( R k − L k )= 4 kN (cid:16) (1 − ρ k ) ω − ρ k (1 − ω ) (cid:17) , (55) D k = ( δρ k ) δt ( R k + L k )= 2 kN N k (cid:16) (1 − ρ k ) ω + ρ k (1 − ω ) (cid:17) . (56)1From equations (55) and (56) we observe D k is sup-pressed by a factor of 1 /N k relative to v k . This suggeststhat the dynamics of the densities { ρ k } are dominatedby the drift velocities { v k } when N k is large. Under thisassumption, we posit that the dynamics of the systemmay be well described by the set of equations dρ k dt = 4 kN (cid:16) (1 − ρ k ) ω − ρ k (1 − ω ) (cid:17) (57)where 1 ≤ k ≤ k max and k max is the maximum nodedegree in the network. The number of equations in thissystem is equal to the number of distinct node degrees.By symmetry, the system has fixed points when all ofthe ρ k ’s are equal. Suppose ρ k = ρ (cid:48) for all k , where ρ (cid:48) ∈ [0 , ω from equation(50) we observe that ω = ρ (cid:48) and therefore fixed pointsoccur for ρ (cid:48) ∈ { , / , } . The two states of uniformconsensus are asymptotically stable. When ρ k = 1 / k , linear stability analysis is conducted numerically. B. Numerical Implementation
In this illustrative example, we consider a systemwhere we impose that 1 / / / dρ dt = 4 N (cid:16) ω (1 − ρ ) − ρ (1 − ω ) (cid:17) ,dρ dt = 8 N (cid:16) ω (1 − ρ ) − ρ (1 − ω ) (cid:17) ,dρ dt = 12 N (cid:16) ω (1 − ρ ) − ρ (1 − ω ) (cid:17) . Numerical analysis of the Jacobian reveals that(1 / , / , /
2) is a saddle point. Figure 8 shows a sampleset of simulation results on such a network, where theinitial conditions are given by ( ρ (0) , ρ (0) , ρ (0)) =(0 . , . , . N becomes large. VI. CONCLUSION
In this paper, we have argued that hypergraphs providea natural, and efficient, framework to explore the rela-tions between structure and dynamics in situations wherebasic interaction units involve more than two nodes. Wehave conducted an in-depth study of a generic modelfrom opinion dynamics, the Majority Rule (MR), andhave proposed a number of hypergraph models to studyMR dynamics. Our analysis was achieved my recast-ing the dynamics in terms of Fokker-Planck equations.By deriving the Fokker-Planck equation associated withMR dynamics on a given hypergraph, the drift and dif-fusion terms governing the evolution of the system couldbe deduced. Interestingly, it was found that diffusivecontributions to the stochastic dynamics vanished as thepopulation size became increasingly large. This allowedfor the system to be modelled using deterministic non-linear dynamical systems which, in essence, reflected thenon-conservative nature of MR. On all of the hyper-graph topologies considered, the deterministic modellingapproach proved to be an excellent predictor of the be-haviour of the system, which allowed for the final stateof consensus to be predicted as a function of the initialconditions. This is markedly different to the behaviourobserved for VM dynamics on heterogeneous networks[5], where diffusion is non-negligible even in large popu-lations.There are a number of ways in which the analysis pre-sented in this paper could be extended. It would be de-sirable to find a way in which to accurately model MRdynamics on small hypergraphs. In this case, making an-alytical progress could be difficult as diffusion would benon-negligible, and as mean-field approximations wouldalso be expected to be less relevant. However, a numeri-cal study of the stochastic differential equations derivedfrom the associated Fokker-Planck equations could stillyield powerful insights. It would also be desirable to finda way in which to apply the analysis presented in SectionV to more sophisticated underlying network structures,as encoded by their adjacency tensor, or including inter-actions of arbitrary size. Finally, it would be interestingto find a way in which to characterise the consensus timeas a function of the population size. Krapivksy and Red-ner [7] succeeded in doing so for the mean field case,though it remains to be investigated on more complexhypergraph topologies.
ACKNOWLEDGMENTS
RL would like to thank Michael Schaub and LeonieNeuhauser for inspiring discussions related to thismanuscript.2 (a) ( ρ , ρ ) (b) ( ρ , ρ ) (c) ( ρ , ρ ) FIG. 8: Two-dimensional phase plane projections associated with the heterogeneous network described in Section V B. In thisinstance, N = 1 . × . VII. LINEAR STABILITY ANALYSIS ON THE SYMMETRIC MODULAR HYPERGRAPH
Here we provide additional details of the analysis presented in Section IV A. We conduct linear stability analysis ofthe system dρ A /dt = v A ,dρ B /dt = v B , where v A and v B are given in equations (28) and (29) respectively. Let J denote the associated Jacobian matrix.The fixed points of the system are given by (0 , , (1 / , /
2) and (1 , ,
0) and (1 ,
1) are asymptoticallystable. The analysis of the fixed point at (1 / , /
2) is slightly more complicated. The matrix J (1 / , /
2) haseigenvalues λ ± given by λ ± = 32 (cid:18) NN + C (cid:19) − C (2 N + C )2 ( N + C ) ± C (2 N + C )( N + C ) .λ + is always positive, therefore (1 / , /
2) is unstable for all C . For C = 0 we have λ − > C = N we have λ − <
0. In order to determine the value of C for which (1 / , /
2) transitions between a saddle point and a source,we solve the equation λ − ( C ) = 0. There are two possible values of C (denoted by C ± ) for which this equation issatisfied: C ± = (cid:32) − ± √ (cid:33) N.C − is negative, therefore it is nonphysical. Hence we conclude that the the fixed point at (1 / , /
2) changes in naturewhen C + ≡ C ∗ ≈ . N .Visual inspection of the ( ρ A , ρ B ) phase plane reveals that two new unstable fixed points emerge along the line ρ B = 1 − ρ A when C becomes sufficiently small. Furthermore, numerical simulations indicate that for some criticalvalue of C , which we shall refer to as C t , the two new fixed points transition between being stable and unstable. That isto say for C ∈ [0 , C t ] there exist four absorbing states in the phase plane. In order to determine the coordinates of thetwo new fixed points, and the range of values of C for which they exist, the drift velocities are evaluated along the line ρ B = 1 − ρ A . We do so by rewriting equations the dynamical system as dρ A /dt = f ( ρ A , ρ B ) and dρ B /dt = g ( ρ A , ρ B )respectively where f and g are functions of the clique densities. By symmetry, it suffices to consider only one ofthese equations given that we are interested in seeking solutions to the equation f ( ρ A , − ρ A ) = 0. This equationultimately reduces to a quadratic in ρ A , yielding the two solutions given in equation (30): ρ ± A = 12 (cid:32) ± (cid:115) CC − N N + C (cid:33) . ρ ± A in conjunction with the relation ρ B = 1 − ρ A gives the coordinates of the two new fixedpoints. As C tends to 0 equation (30) predicts that the fixed points tend to the coordinates (0 ,
1) and (1 , C for which these fixed points exist by looking at where the discriminant inequation (30) is positive. Solving the resultant quadratic equation reveals that the fixed points exist on the interval C ∈ [0 , C ∗ ]. This illustrates that as we approach C ∗ from above, the two new fixed points spontaneously emerge atthe exact instant when the saddle point at (1 / , /
2) becomes a source. In order to determine the stability of thesetwo fixed points, once again we proceed using linear stability analysis. When evaluating the Jacobian along the line ρ B = 1 − ρ A we find that J = J and J = J , where J ij denotes the Jacobian entry in row i and column j . Therefore the eigenvalues are given by λ ± = J ± J . Substituting ( ρ + A , − ρ + A ) into J yields the followingexpressions for the Jacobian entries: J = 7 C + 14 CN − N ( N + C ) , J = 2 C (2 N + C )( N + C ) . We note that J is always positive, therefore in order to determine the point at which these fixed points becomestable it suffices to consider solutions to the equation J + J = 0. This ultimately reduces to a quadratic equationin C , where the non-negative root C t is found to be C t = (cid:32) − √ (cid:33) N ≈ . N. The fixed points are therefore stable for C ∈ [0 , C t ]. In summary, C t is the parameter value below which metastablestate formation is predicted by the deterministic model on the symmetric modular hypergraph. VIII. EXPLICIT PARAMETER EXPRESSIONS ON THE ASYMMETRIC MODULAR HYPERGRAPH
Here we provide explicit expressions for ρ ± f and ρ ± g , the solutions to equations (39) and (40) respectively: ρ ± f = 12 ± (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) (cid:18) NC BA ( N + C BA ) + (cid:16) C AB N + C AB (cid:17) (cid:19) C AB (4 N − C AB )( N + C AB ) + 2 C BA ( C BA − N )( N + C BA ) − N ( N + C AB ) ρ ± g = 12 ± (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) (cid:18) NC AB ( N + C AB ) + (cid:16) C BA N + C BA (cid:17) (cid:19) C BA (4 N − C BA )( N + C BA ) + 2 C AB ( C AB − N )( N + C AB ) − N ( N + C BA ) . v A = 3 (cid:18) NN + C AB (cid:19) ρ A (1 − ρ A ) (2 ρ A − N C AB ( N + C AB ) ρ A (1 − ρ A ) (2 ρ B − (cid:18) C BA N + C BA (cid:19) ρ A (1 − ρ A ) (2 ρ B − N C BA ( N + C BA ) (cid:16) ρ B (1 − ρ A ) − ρ A (1 − ρ B ) (cid:17) + (cid:18) C AB N + C AB (cid:19) (cid:16) ρ B (1 − ρ A ) − ρ A (1 − ρ B ) (cid:17) , (58) v B = 3 (cid:18) NN + C BA (cid:19) ρ B (1 − ρ B ) (2 ρ B − N C BA ( N + C BA ) ρ B (1 − ρ B ) (2 ρ A − (cid:18) C AB N + C AB (cid:19) ρ B (1 − ρ B ) (2 ρ A − N C AB ( N + C AB ) (cid:16) ρ A (1 − ρ B ) − ρ B (1 − ρ A ) (cid:17) + (cid:18) C BA N + C BA (cid:19) (cid:16) ρ A (1 − ρ B ) − ρ B (1 − ρ A ) (cid:17) . (59)Equations (58) and (59) reduce to equations (28) and (29) respectively when C AB = C BA = C . [1] C. Castellano, S. Fortunato, and V. Loreto, Reviews ofmodern physics , 591 (2009).[2] P. Clifford and A. Sudbury, Biometrika , 581 (1973).[3] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lu-cas, A. Patania, J.-G. Young, and G. Petri, Physics Re-ports , 1 (2020).[4] S. Redner, A guide to first-passage processes (CambridgeUniversity Press, 2001).[5] V. Sood, T. Antal, and S. Redner, Physical Review E , 041121 (2008).[6] S. Galam, The European Physical Journal B-CondensedMatter and Complex Systems , 403 (2002).[7] P. L. Krapivsky and S. Redner, Physical Review Letters , 238701 (2003).[8] M. Friedman and R. D. Friedman, Tyranny of the statusquo (Penguin, 1985).[9] R. Lambiotte, M. Rosvall, and I. Scholtes, NaturePhysics , 313–320 (2019).[10] C. Giusti, E. Pastalkova, C. Curto, and V. Itskov, Pro-ceedings of the National Academy of Sciences , 13455(2015).[11] M. W. Reimann, M. Nolte, M. Scolamiero, K. Turner,R. Perin, G. Chindemi, P. D(cid:32)lotko, R. Levi, K. Hess, andH. Markram, Frontiers in Computational Neuroscience , 48 (2017).[12] F. A. N. Santos, E. P. Raposo, M. D. Coutinho-Filho,M. Copelli, C. J. Stam, and L. Douw, Phys. Rev. E , 032414 (2019). [13] R. Olfati-Saber, J. A. Fax, and R. M. Murray, Proceed-ings of the IEEE , 215 (2007).[14] A. Patania, G. Petri, and F. Vaccarino, EPJ Data Sci-ence , 18 (2017).[15] J. J. Torres and G. Bianconi, Journal of Physics: Com-plexity , 015002 (2020).[16] M. E. Newman, D. J. Watts, and S. H. Strogatz, Pro-ceedings of the National Academy of Sciences , 2566(2002).[17] N. Lanchier and J. Neufer, Journal of Statistical Physics , 21 (2013).[18] L. Neuh¨auser, M. T. Schaub, A. Mellor, and R. Lam-biotte, arXiv preprint arXiv:2004.00901 (2020).[19] P. Chen and S. Redner, Physical review E , 036101(2005).[20] T. Kumar, S. Vaidyanathan, H. Ananthapadmanabhan,S. Parthasarathy, and B. Ravindran, arXiv preprintarXiv:1812.10869 (2018).[21] P. Chodrow and A. Mellor, Applied Network Science ,9 (2020).[22] M. E. Newman, Proceedings of the national academy ofsciences , 8577 (2006).[23] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, andE. Lefebvre, Journal of statistical mechanics: theory andexperiment , P10008 (2008).[24] R. Lambiotte, M. Ausloos, and J. Ho(cid:32)lyst, Physical Re-view E , 030101 (2007).[25] N. Landry and J. G. Restrepo, arXiv preprint5