Effects of heterogeneous social interactions on flocking dynamics

Effects of Heterogeneous Social Interactions on Flocking Dynamics M. Carmen Miguel,1,2 Jack T. Parley,3and Romualdo Pastor-Satorras3 1Departament de Física de la Mat`eria Condensada, Universitat de Barcelona, Martí i Franqu`es 1, 08028 Barcelona, Spain 2Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, 08028 Barcelona, Spain 3Departament de Física, Universitat Polit`ecnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain (Received 19 September 2017; revised manuscript received 30 November 2017; published 8 February 2018) 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. DOI: 10.1103/PhysRevLett.120.068303 Collective motion in living and complex systems [1], where simple interactions between constituent entities produce striking spatiotemporal 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 well-ordered 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 of 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–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–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 kother individuals, showing a scale-free signature [26] of the form PðkÞ∼k−γd. We study the behavior of the Vicsek model when applied to complex PHYSICAL REVIEW LETTERS 120, 068303 (2018) 0031-9007=18=120(6)=068303(5) 068303-1 © 2018 American Physical Society 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 been theoretically and numerically analyzed in scale-free networks under various conditions [29–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 aij [11], with value aij ¼1, if individuals iand jare socially connected, while aij ¼0otherwise. 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 viðtÞ, assumed to be normalized: jviðtÞj ¼ v0. We fix v0¼1. We consider velocities in two dimension, fully determined by the angle θiðtÞthey form with, say, the xaxis, i.e., viðtÞ¼ fcos θiðtÞ;sin θiðtÞg. With the original definition of the model [5], velocities are synchronously updated via the rule θiðtþ1Þ¼ΘviðtÞþX N j¼1 aijvjðtÞþηξiðtÞ;ð1Þ where Nis 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 η¼1is 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] ϕηðtÞ¼1 NX N i¼1 viðtÞ :ð2Þ From here, one defines the average order parameter hϕηi¼limT→∞ð1=TÞRT 0ϕηðtÞdt and the susceptibility χη¼ N½hϕ2 ηi−hϕηi2, which close to the critical point behave as hϕηi∼ðηc−ηÞβand χη∼jηc−ηj−γ, respectively, defining the critical exponents βand γ, in analogy with the ferromagnetic phase transition [7]. The model defined by Eq. (1) 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: ϕηðtÞ¼1 NX N i¼1 cos ½θiðtÞ−¯ θðtÞ;ð3Þ where ¯ θðtÞ¼Θ½PN i¼1viðtÞ. Equation (3) can be shown to be exactly equal to Eq. (2); see the Supplemental Material [34]. For a fully connected network, the Vicsek model can be solved starting from Eq. (3) (see Supplemental Material [34]), obtaining the result that the system is ordered for any η<1. In the vicinity of this point, expansions of the solution lead to hϕηi∼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 Tis a function of the noise intensity η, fulfilling the limits T→0for η→0,andT→∞for η→1.TheXY model applied to networks can be solved within an annealed network approximation, obtaining a critical temperature Tc¼Jhk2i=½2hki [29],whereJis 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 [35] 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 γdand a minimum degree of m¼3[27]. The order parameter, hϕηi, is computed by averaging over 50000 time steps, after letting the system initially relax for 10000 time steps. Figure. 1 (inset) shows a plot of the average order parameter as 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 γd = 2.10 γd = 2.40 γd = 2.75 γd = 3.50 0.5 0.6 0.7 0.8 0.9 1.0 η 0.0 0.1 0.2 0.3 <φη> N = 106 N = 105 N = 104 γd = 3.50 γd = 2.10 FIG. 1. Inset: Average order parameter as a function of the noise intensity ηfor different values of the degree exponent γdin UCM networks of size N¼106. Main: The order parameter as a function of ηfor different values of the network size N. The sets of plots correspond to γd¼2.1(right) and γd¼3.5(left). PHYSICAL REVIEW LETTERS 120, 068303 (2018) 068303-2 function of η, computed in networks of size N¼106with different degree exponents. Figure 1(main) illustrates the effects of system size for two different values of γd. From this figure it is apparent that for small γdthe 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: χNðηÞ¼Nhϕ2 ηi−hϕηi2 hϕηi;ð4Þ which is customarily used to detect phase transitions in complex networks [36,37]. The effective critical point, ηcðNÞ, will be given, for a given network size N,bythe 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 [38] of the form ηcðNÞ¼ηc−aN−1=ν;ð5Þ where νis another characteristic critical exponent [36,39]. The height of the peak of the dynamic susceptibility, χpeak N, also scales with N, adopting the form [36,37] χpeak N∼NðβþγÞ=ν:ð6Þ In Fig. 2we 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 Nincreases. We proceed to estimate the critical point in the thermodynamic limit by applying a nonlinear fit to the position of the peak, ηcðNÞ, as a function of N, according to Eq. (5);seeFig.3and Table I. From these results, it is apparent that, for γd>2.50, the critical point ηctends to a constant value of less 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 nonlinear 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 scalefree 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 0.6 0.7 0.8 0.9 η 102 104 χN N = 106 N = 3 x 105 N = 105 γd = 3.50 γd = 2.75 γd = 2.40 γd = 2.10 FIG. 2. Numerical dynamic susceptibility as a function of noise amplitude ηin UCM networks of different size. The groups of functions for different Ncorrespond to the values of γd(from left to right): 3.50, 2.75, 2.50, and 2.10. TABLE I. Critical point and exponent δfor the Vicsek model in scale-free networks with different degree exponent. γd2.10 2.25 2.40 2.50 2.75 3.50 ηc0.99(1) 1.00(5) 1.00(5) 1.00(1) 0.71(1) 0.61(1) δ0.574(3) 0.63(1) 0.71(1) 0.74(2) 0.77(3) 0.75(1) 104105106 10-1 ηc − ηc(N) 104105 10-4 10-2 104105106 N 102 103 104 χN peak γd = 2.10 γd = 2.25 γd = 2.40 γd = 2.50 γd = 2.75 γd = 3.50 (a) (b) (c) FIG. 3. Finite size scaling analysis of the position of the susceptibility peak, according to Eq. (5) for small (a) and large (b) values of the degree exponent γd. The critical points, in the thermodynamic limit, are given 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. PHYSICAL REVIEW LETTERS 120, 068303 (2018) 068303-3 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 γdvalues 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 different 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 fmax ¼1=2[40]. On a fully connected graph, the majority-vote model shows a critical point fc¼1=2, which in the Vicsek case translates to ηc¼1with exponents β¼1and γ¼0, see the Supplemental Material [34]. Meanwhile, in heterogeneous networks with a power-law degree distribution, a threshold fc¼1=2− ffiffiffiffiffiffiffiffiffiffiffiffi ðπ=8Þ p½hki=ðhk3=2iÞ has recently been reported [32]. This threshold shows a transition from fc<1=2for γd>5=2to fc¼1=2for γd<5=2in 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 [41]. In order to confirm the equivalence of the majorityvote 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 the Supplemental Material [34]. From our simulations, we confirm the results in Ref. [32] regarding a threshold fc→1=2in the thermodynamic limit for γd<5=2, while fc<1=2for γd>5=2. The estimation of the exponent δ for the growth of the dynamical susceptibility peak with network size, Eq. (6), 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 ϕηðt;kÞ¼ 1 NkX i∈Vk cos ½θiðtÞ−¯ θðtÞ;ð7Þ where Vkis the set of nodes with degree k, and Nkis the number of such nodes. From this expression, a timeindependent order parameter hϕηðkÞi 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 than 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, which 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, 101102 k 10-2 10-1 100 <φη(k)> γd = 2.10 γd = 2.40 γd = 2.75 γd = 3.50 FIG. 4. Restricted order parameter as a function of degree for UCM networks of size N¼3×105and different degree exponent values γd. Outline symbols correspond to a noise strength equal to the peak of the dynamic susceptibility; solid symbols correspond to a noise strength 3% above the peak. PHYSICAL REVIEW LETTERS 120, 068303 (2018) 068303-4 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 nonequilibrium majority-vote model on complex networks, which can be considered as a variation 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 opens up intriguing questions regarding the role of symmetries in dynamical processes on networks. Deeper research effort is necessary to further our understanding of both questions. We acknowledge financial support from the Spanish Government’s MINECO, under projects FIS2013-47282C2-1-P, FIS2013-47282-C2-2, FIS2016-76830-C2-1-P, and FIS2016-76830-C2-2-P. Also, R. P.-S. acknowledges additional financial support from ICREA Academia, funded by the Generalitat de Catalunya regional authorities. [1] T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012). [2] D. J. T. Sumpter, Phil. Trans. R. Soc. B 361, 5 (2006). [3] D. J. Sumpter, Collective Animal Behavior (Princeton University Press, Princeton, NJ, 2010). [4] T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995). [5] F. Ginelli, Eur. Phys. J. Spec. Top. 225, 2099 (2016). [6] V. M´endez, D. Campos, and F. Bartumeus, Stochastic Foundations in Movement Ecology (Springer Verlag, Berlin, Heidelberg, 2014). [7] J. M. Yeomans, Statistical Mechanics of Phase Transitions (Oxford University Press, Oxford, 1992). [8] I. Couzin, J. Krause, R. James, G. Ruxton, and N. Franks, J. Theor. Biol. 218, 1 (2002). [9] H. Chat´e, F. Ginelli, G. Gr´egoire, F. Peruani, and F. Raynaud, Eur. Phys. J. B 64, 451 (2008). [10] D. Croft, R. James, and J. Krause, Exploring Animal Social Networks (Princeton University Press, Princeton, NJ, 2008). [11] M. Newman, Networks: An Introduction (Oxford University Press, Inc., New York, 2010). [12] D. Lusseau and M. E. J. Newman, Proc. R. Soc. B 271, S477 (2004). [13] J. C. Flack, D. C. Krakauer, and F. B. M. de Waal, Proc. R. Soc. B 272, 1091 (2005). [14] M. Rhodes, G. W. Wardell-Johnson, M. P. Rhodes, and B. Raymond, Conservation biology : the Journal of the Society for Conservation Biology 20, 861 (2006). [15] D. P. Croft, J. Krause, and R. James, Proc. R. Soc. B 271, S516 (2004). [16] D. P. Croft, R. James, A. J. W. Ward, M. S. Botham, D. Mawdsley, and J. Krause, Oecologia 143, 211 (2005). [17] S. B. Rosenthal, C. R. Twomey, A. T. Hartnett, H. S. Wu, and I. D. Couzin, Proc. Natl. Acad. Sci. 112, 4690 (2015). [18] M. Aldana and C. Huepe, J. Stat. Phys. 112, 135 (2003). [19] M. Aldana, V. Dossetti, C. Huepe, V. M. Kenkre, and H. Larralde, Phys. Rev. Lett. 98, 095702 (2007). [20] J. A. Pimentel, M. Aldana, C. Huepe, and H. Larralde, Phys. Rev. E 77, 061138 (2008). [21] A. Sekunda, M. Komareji, and R. Bouffanais, Netw. Sci. 4, 244 (2016). [22] N. W. Bode, A. J. Wood, and D. W. Franks, Animal Behaviour 82, 29 (2011). [23] N. W. F. Bode, A. J. Wood, and D. W. Franks, Behav. Ecol. Sociobiol. 65, 117 (2011). [24] D. Lusseau, Proc. R. Soc. B Biol. Sci. 270, S186 (2003). [25] T. G. Manno, Animal Behaviour 75, 1221 (2008). [26] A.-L. Barabási and R. Albert, Science 286, 509 (1999). [27] M. Catanzaro, M. Boguñá, and R. Pastor-Satorras, Phys. Rev. E 71, 027103 (2005). [28] A. Czirók, H. E. Stanley, and T. Vicsek, J. Phys. A 30, 1375 (1997). [29] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008). [30] W. Kwak, J. Yang, J. I. Sohn, and I. M. Kim, Phys. Rev. E 75, 061130 (2007). [31] J.-S. Yang, W. Kwak, K.-I. Goh, and I.-m. Kim, Europhys. Lett. 84, 36004 (2008). [32] H. Chen, C. Shen, G. He, H. Zhang, and Z. Hou, Phys. Rev. E91, 022816 (2015). [33] See, however, Ref. [20] for a solution in the case of random regular (homogeneous) networks. [34] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.120.068303 for additional analytical calculations and numerical results. [35] K. Medvedyeva, P. Holme, P. Minnhagen, and B. J. Kim, Phys. Rev. E 67, 036118 (2003). [36] S. C. Ferreira, C. Castellano, and R. Pastor-Satorras, Phys. Rev. E 86, 041125 (2012). [37] C. Castellano and R. Pastor-Satorras, Eur. Phys. J. B 89, 243 (2016). [38] Finite Size Scaling, Current Physics-Sources and Comments, Vol. 2 edited by J. L. Cardy (North Holland, Amsterdam, 1988). [39] K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical Physics, 5th ed. (Springer-Verlag, Berlin, 2010). [40] T. Liggett, Interacting Particle Systems (Springer Verlag, Berlin, Heidelberg, 1997). [41] L. F. C. Pereira and F. G. Brady Moreira, Phys. Rev. E 71, 016123 (2005). PHYSICAL REVIEW LETTERS 120, 068303 (2018) 068303-5