Effects of heterogeneous social interactions on flocking dynamics

Social relationships characterize the interactions that occur within social species and may have an important impact on collective animal motion. Here, we consider a variation of the standard Vicsek model for collective motion in which interactions are mediated by an empirically motivated scale-free topology that represents a heterogeneous pattern of social contacts. We observe that the degree of order of the model is strongly affected by network heterogeneity: more heterogeneous networks show a more resilient ordered state; while less heterogeneity leads to a more fragile ordered state that can be destroyed by sufficient external noise. Our results challenge the previously accepted equivalence between the {\em static} Vicsek model and the equilibrium XY model on the network of connections, and point towards a possible equivalence with models exhibiting a different symmetry.

Social relationships characterize the interactions that occur within social species and may have an important impact on collective animal motion. Here, we consider a variation of the standard Vicsek model for collective motion in which interactions are mediated by an empirically motivated scale-free topology that represents a heterogeneous pattern of social contacts. We observe that the degree of order of the model is strongly affected by network heterogeneity: more heterogeneous networks show a more resilient ordered state; while less heterogeneity leads to a more fragile ordered state that can be destroyed by sufficient external noise. Our results challenge the previously accepted equivalence between the static Vicsek model and the equilibrium XY model on the network of connections, and point towards a possible equivalence with models exhibiting a different symmetry.
Collective motion in living and complex systems [1], where simple interactions between constituent entities produce striking spatio-temporal patterns on scales larger than the entities themselves, are commonplace. Some of the examples that best highlight the emergence of such patterns are found in animal motion [2,3], where the animals collectively exhibit some of the most spectacular and fascinating sights in nature. These include flocks of birds turning in unison or migrating in wellordered formation, shoals of fish splitting and reforming as they outmaneuver a predator, seasonal migratory herds of large herbivores, etc.
The challenge of understanding how hundreds or thousands of organisms move together and give rise to such intriguing collective responses in the absence of any apparent leader or driving field has attracted the attention of the scientific community for a long time. Significant progress in understanding how some of these features come about has been achieved through the development of relatively simple models of self-propelled particles (SPP). In SPP models, the complex dynamics of individuals within a group are simplified to those of particles that move with given velocities and experience flocking interactions within a local interaction zone, combined with random fluctuations due to intrinsic or environmental factors. In the celebrated Vicsek model [4], these interactions consist in the alignment of the velocity of an SPP with the average velocity of some of its neighbors. Perfect alignment is, however, impeded by the addition of a noise term that mimics, for instance, the difficulties in gathering and processing the surrounding information. The success of the model lies in the production of a phase transition as a function of noise intensity, η, separating an ordered or polarized (flocking) phase at η ≤ η c , where particles travel in a common direction, from a disordered phase for η > η c , where particles behave as uncorrelated persistent random walkers [5,6]. This is particularly fruitful due to the analogies that can be drawn between the self-organization of herds of moving animals and standard phase transitions observed in condensed matter [1,7].
The main assumption of the Vicsek and other similar models of collective motion [8,9] is that particles tend to orient their velocity parallel to the average velocity in a local neighborhood, independently of their identity. This kind of interaction rule leaves aside, however, the important fact that real interactions between moving animals can be more intricate. One source of complication can be the presence of social interactions [10] between the group members, which can lead, in the framework of the Vicsek model, to a tendency to align one's velocity with that of individuals with which one has strong social ties, but that might be separated by a relatively long Euclidean distance. The presence of such social interactions, naturally represented in terms of social networks [11], has been observed in mammals [12][13][14] and fish [15,16], and has even been studied in the context of schooling fish [17].
The impact of social interactions given in terms of networks has already been considered in the context of collective motion and the Vicsek model [18][19][20][21][22][23]; but, to the best of our knowledge, an in-depth study is still lacking. Here, we focus on the effects of the topological heterogeneity observed in certain animal social networks [24,25], which can be represented by a degree distribution P (k), defined as the probability that a randomly chosen individual is connected to k other individuals, showing a scale-free signature [26] of the form: We study the behavior of the Vicsek model when applied to complex networks with varying heterogeneity (a varying degree exponent γ d ), generated using the uncorrelated configuration model (UCM) [27]. In this setting, each particle's neighbors always remain the same. As a consequence, it is usually assumed that, in this limit, the Vicsek model must be equivalent to the equilibrium XY model of ferromagnetism defined on the network of connections (see e.g. Refs. [4,5,28]). The XY model has arXiv:1801.03371v1 [physics.soc-ph] 10 Jan 2018 been theoretically and numerically analyzed in scale-free networks under various conditions [29][30][31]. By means of extensive numerical simulations, we show that on static scale-free networks the Vicsek and the XY models exhibit different critical behavior. Furthermore, our simulations are compatible with the behavior reported for the nonequilibrium majority-vote model with noise applied to complex networks [32].
We consider a version of the Vicsek model in which interactions are mediated by a static complex network, with links representing social interactions. A social network can be fully represented in terms of its adjacency matrix a ij [11], with value a ij = 1, if individuals i and j are socially connected, while a ij = 0 otherwise. By considering an ordering dynamics based on social interactions alone, we disregard spatial position, and thus the SPPs are uniquely specified in terms of their velocity v i (t), assumed to be normalized: |v i (t)| = v 0 . We fix v 0 = 1. We consider velocities in two dimension, fully determined by the angle θ i (t) they form with, say, the x axis, i.e. v i (t) = {cos θ i (t), sin θ i (t)}. With the original definition of the model [5], velocities are synchronously updated via the rule: where N is the network size, Θ[V] represents the angle described by vector V, ξ i (t) is random noise uniformly distributed within the interval [−π, π], and η ∈ [0, 1] is a parameter that reflects the noise strength. Note that η = 1 is the maximum possible noise, since it corresponds to a completely disordered system. The phase transition between ordered and disordered states in the Vicsek model is determined by the temporal evolution of an order parameter φ η (t), defined as [4]: From here, one defines the average order parameter , which close to the critical point behave as φ η ∼ (η c − η) β and χ η ∼ |η c − η| −γ , respectively, defining the critical exponents β and γ, in analogy with the ferromagnetic phase transition [7].
The model defined by Eq. (17) does not admit a feasible analytical treatment for general networks [33]. We can, however, solve it in the fully connected case. To proceed, it is convenient to write the order parameter in the alternative form: . Eq. (15) can be shown to be exactly equal to Eq. (19), see Appendix . For a fully connected network, the Vicsek model can be solved starting from Eq. (15) (see Appendix ), obtaining the result that the system is ordered for any η < 1. In the vicinity of this point, expansions of the solution lead to φ η ∼ 1 − η and χ η ∼ const., leading to the critical exponents β = 1, γ = 0.
In the case of sparse networks, it is usually assumed that, when the particles are immobile and the network of connections is sufficiently dense, the Vicsek model can be mapped to the equilibrium XY model [5], where the temperature T is a function of the noise intensity η, fulfilling the limits T → 0 for η → 0, and T → ∞ for η → 1. The XY model applied to networks can be solved within an annealed network approximation, obtaining a critical temperature T c = J k 2 /[2 k ] [29], where J is the coupling constant of the XY Hamiltonian. That is, for scale-free networks with γ d > 3, there is a true transition at a finite critical temperature; while for γ d ≤ 3, there is no transition and the system is always ordered for any finite T . These results have been confirmed by numerical simulations on heterogeneous [31] and homogeneous [34] networks.
In order to check the validity of the mapping to the XY model, we performed numerical simulations of the Vicsek model on UCM networks with different values of γ d and a minimum degree of m = 3 [27]. The order parameter, φ η , is computed by averaging over 50, 000 time steps, after letting the system initially relax for 10, 000 time steps. Fig. 1(inset) shows a plot of the average order parameter as a function of η, computed in networks of size N = 10 6 with different degree exponents. Fig. 1 (main) illustrates the effects of system size for two different values of γ d . From this figure it is apparent that for small γ d the effective threshold depends strongly on N . In order to explore size effects in greater detail, we proceed to compute the effective threshold by looking at the dynamic susceptibility: which is customarily used to detect phase transitions in complex networks [35,36]. The effective critical point, η c (N ), will be given, for a given network size N , by the position of the maximum of the dynamic susceptibility χ N (η). The critical point in the thermodynamic limit N → ∞ can be obtained by applying a finite-size scaling hypothesis [37] of the form: where ν is another characteristic critical exponent [35,38]. The height of the peak of the dynamic susceptibility, χ peak N , also scales with N , adopting the form [35,36]: In Fig. 2 we plot the dynamic susceptibility, χ N (η), for networks with different values of the degree exponent γ d . As can be seen from the figure, for γ d > 2.5, the location of the peak of the susceptibility appears to tend to a constant value smaller than 1. In contrast, for γ d < 2.5, this location shifts to larger values of η as N increases. We proceed to estimate the critical point in the thermodynamic limit by applying a non-linear fit to the position of the peak, η c (N ), as a function of N , according to Eq. (5); see Fig. 3 and Table I Table I. c) Scaling of the peak of susceptibility with network size for different values of γ d . The associated exponents χ peak n ∼ N δ , with δ = (β + γ)/ν, are given in Table I. than 1; while for γ d ≤ 2.50, the critical point tends to 1 in the thermodynamic limit. Therefore, in this latter case, the order-disorder transition characteristic of the model is suppressed, and the system is fully ordered for any η < 1 in sufficiently large networks. In contrast, for γ d > 2.50, there is a true order-disorder transition, which is preserved even in the limit of infinite network size. While the exponent ν is difficult to estimate due to statistical fluctuations in the non-linear fitting procedure, the exponent δ ≡ (β +γ)/ν, controlling the growth of the dynamic susceptibility peak, can be reliably computed; see Fig. 3 and Table I. The exponents obtained are again compatible with a radical difference in behavior between γ d ≥ 2.50, for which we obtain δ 0.75; and γ d < 2.50, where δ is an increasing function of γ d .
The numerical results obtained for heterogeneous scale-free networks provide a clear picture: when dealing with networks, the Vicsek model cannot be directly mapped to the XY model. The main evidence of this incompatibility comes from the behavior of the critical point. For the XY model, one expects a finite critical temperature (i.e., η c < 1) for γ d > 3; and an infinite critical temperature (i.e., η c = 1), or in other words, no phase transition, for γ d < 3. Meanwhile, when applied to networks, the Vicsek model only produces a true order-disorder transition for degree exponents larger than γ d = 5/2. Experimental characterizations of the degree exponent in groups of social animals (despite the conceivable difficulties associated with measuring it and the fact that it probably varies depending on the behavioral test chosen) provide γ d values within the range 1 − 3.5 [24,25]. It is also obvious that the thermodynamic limit cannot be achieved in experiments. Nevertheless, our results could have important consequences for the resilience of the ordered phase observed in differ-ent species, according to the heterogeneity of their social contact distribution. Strongly heterogeneous networks show a resilient ordered phase for the whole range of disorder values; while low heterogeneity leads to a more fragile ordered phase that can be destroyed by a sufficient amount of external noise.
In order to shed some light on the behavior observed, we put forward the following hypothesis: given that in networks, the dimensionality of the order parameter appears to be irrelevant (for example, the Ising and XY models share the same scaling of the critical point and the same critical exponents [29]), we conjecture that a model analogous to the Vicsek model, but with a scalar order parameter, might also share the same behavior as the Vicsek model in heterogeneous networks. In this way, we consider the majority-vote model [19], in which spin variables on the vertices of a network update their state taking the value of the majority of their nearest neighbors. This state is randomly flipped with a probability f , which plays a similar role to the noise strength, η, but takes a maximum value f max = 1/2 [39]. On a fully connected graph, the majority-vote model shows a critical point f c = 1/2, which in the Vicsek case translates to η c = 1 with exponents β = 1 and γ = 0, see Appendix . Meanwhile, in heterogeneous networks with a power-law degree distribution, a threshold f c = 1/2 − π 8 k k 3/2 has recently been reported [32]. This threshold shows a transition from f c < 1/2 for γ d > 5/2 to f c = 1/2 for γ d < 5/2 in the thermodynamic limit: in full agreement with the observations of the Vicsek model applied to networks. Moreover, above the threshold degree exponent, γ d = 5/2, the value of the exponent (β + γ)/ν 0.75 is also in agreement with the mean-field values of the majority-vote mode: β = 1/2, γ = 1, ν = 2 [40]. In order to confirm the equivalence of the majority-vote and Vicsek models on heterogeneous networks, we have performed additional extensive simulations of the latter for a range of different degree exponents on UCM networks. The results obtained are described in Appendix . From our simulations, we confirm the results in Ref. [32] regarding a threshold f c → 1/2 in the thermodynamic limit for γ d < 5/2, while f c < 1/2 for γ d > 5/2. The estimation of the exponent δ for the growth of the dynamical susceptibility peak with network size, Eq. (20), leads to the results δ = 0.57(1) for γ d = 2.10, δ = 0.61(1) for γ d = 2.25, δ = 0.67(2) for γ d = 2.40, and δ = 0.78(2) for γ d = 2.75. The excellent agreement of these exponents, compared with the ones for the Vicsek model reported in Table I, confirm our hypothesis regarding the equivalence of Vicsek and majority-vote model on complex networks.
We finally focus on the hierarchy of the order of the nodes of different degree in the Vicsek model, and compute a degree-restricted order parameter defined as where V k is the set of nodes with degree k, and N k is the number of such nodes. From this expression, a time-independent order parameter, φ η (k) , is defined by means of an appropriate time average over a large time window, T . In Fig. 4, we plot the restricted order parameter as a function of k. As can be seen, there is apparently a hierarchy in the order of the systems, with low-degree nodes being more disordered than high-degree nodes. This can be explained by the larger number of connections of high degree nodes, which average velocities over a larger ensemble that low-degree nodes do and are therefore less susceptible to the influence of the external noise. This effect can be interpreted as high degree nodes playing the role of leaders, that can keep the network ordered even close to the maximum possible value of disorder when they are large enough (i.e., for small values of γ d ).
In conclusion, we have studied numerically the Vicsek model applied to complex scale-free networks with a degree distribution P (k) ∼ k −γ d . By means of extensive numerical simulations, we observe that the nature of the possible order-disorder transition exhibited by the model depends on the level of heterogeneity of the network, as given by the value of the degree exponent γ d . For γ d > 2.5, there is a true transition, located at η c (γ d ) which increases with decreasing γ d . Meanwhile, for γ d < 2/5, we obtain a critical point in the thermodynamic limit equal to 1, indicating the lack of a true critical transition. These results indicate that flocking dynamics in scale-free social networks is more robust against noise effects in the case of high network heterogeneity (i.e. small γ d ). These numerical results are in disagreement with the validity of direct mapping of the Vicsek model to the equilibrium XY model on the network of connections, which is usually assumed to be valid. Nonetheless, our results do appear to be in agreement with those of the non-equilibrium majority-vote model on complex networks, which can be considered as a vari-ation of the Vicsek model with reduced symmetry of the order parameter. Our work highlights the role of the effects of social topology in flocking dynamics and open up intriguing questions regarding the role of symmetries in dynamical processes on networks. Deeper research effort in necessary to further our understanding of both questions.

Alternative form of the order parameter
In order to show that the proposed alternative form of the order parameter for the Vicsek model is identical to the original one, we start from the expression Using the cosine angle difference identity, we can write Now, using and substituting into Eq. (9), we obtain which is exactly the form of the temporal order parameter in its classical definition, Eq. (2) in the main paper.

Analytical solution for fully connected networks
In the case of fully connected networks, in which every node is connected to every other one, we have θ i (t) = θ(t) + ηξ i (t), an therefore As ξ i (t) is uncorrelated in both the particle index an time, we can write the average value, in the thermodynamic limit, This value of the average order parameter is different from zero for any η < 1, indicating η c = 1. In the vicinity of this critical point, a Taylor expansion of Eq. (14) leads to φ η ∼ 1 − η, defining the critical exponent β = 1.

Majority vote model on fully connected networks
In order to explore the behavior of the majority vote model on fully connected networks we have performed numerical simulations of the dynamics in which each node takes the majority state of the other N − 1 nodes with probability 1 − f , and the opposite state with probability f . The order parameter φ f is defined as the average absolute value of the magnetization in the steady state. The dynamic susceptibility is, in its turn, defined as In Fig. 5 we plot the results of numerical simulations, averaging over 50000 Monte Carlo time steps, after letting the system relax for 10000 steps. In Fig. 5a) we plot the average order parameter as a function of f . The clear linear behavior indicates a value f c = 1/2 with critical exponent β = 1. In Fig. 5b) we depict instead the dynamics susceptibility χ N (f ) as a function of the probability f . As we can see, it tends to diverge close to f c = 1/2, with rounding effects for small system sizes. In the inset in Fig. 5b) we plot χ N (f ) as a function of 1/( 1 2 − f ). The linear behavior for small values of this quantity leads to

Majority vote model on scale-free networks
We consider the majority vote model on scale-free networks of degree distribution P (k) ∼ k −γ d , generated using the UCM network model. In this case, each node takes the majority state of its nearest neighbors with probability 1 − f , and the opposite state with probability f . In Fig. 6 we plot the dynamic susceptibility χ N (f ), Eq. (18) as a function of the probability f in networks of different size N and degree exponent γ d . This function is evaluated averaging over 50000 Monte Carlo time steps, after letting the system relax to its steady state. As we can see, the effective critical point for a given network size, f c (N ), defined as the position of the peak of the dynamic susceptibility, appears to tend to a constant for large γ d , while it approaches the limit of large noise f = 0.5 in the case of small γ d < 5/2 and large size. Indeed, the value of the critical point f c in the thermodynamic limit, extrapolated from the finite-size scaling ansatz leads to a result in agreement with the theoretical prediction in Ref. [32], namely a value approaching f c 0.5 for γ d < 5/2, and a value f c < 0.5 for γ d > 5/2. In Fig. 7 we plot the height of the peak of the dynamic susceptibility, which should scale with network size as with the exponent δ = (β + γ)/ν. By means of a linear regression in double logarithmic scale, we estimate the values of δ for the different degree exponents considered: γ d = 2.10, δ = 0.57(1); γ d = 2.25, δ = 0.61(1); γ d = 2.40, δ = 0.67(2); γ d = 2.75, δ = 0.78 (2). These exponents are in excellent agreement with the ones obtained for the Vicsek model in scale-free networks (see Table I in the main paper), and confirm the equivalence of behavior of the Vicsek model, with a continuous symmetry, and the majority vote model, with a discrete symmetry, on complex networks with a heterogeneous degree distribution.