Collective predator evasion: Putting the criticality hypothesis to the test
CC OLLECTIVE PREDATOR EVASION : P
UTTING THE CRITICALITYHYPOTHESIS TO THE TEST
Pascal P. Klamser a,b and Pawel Romanczuk a,ba
Department of Biology, Institute for Theoretical Biology, Humboldt-Universit¨at zu Berlin, 10115 Berlin, Germany b Bernstein Center for Computational Neuroscience, 10115 Berlin, GermanySeptember 7, 2020 a r X i v : . [ q - b i o . P E ] S e p PREPRINT - S
EPTEMBER
7, 2020 A BSTRACT
According to the criticality hypothesis , collective biological systems should operate in a special pa-rameter region, close to so-called critical points, where the collective dynamics undergoes a qualita-tive change between different dynamical regimes. Critical systems exhibit unique properties, whichmay benefit collective information processing such as maximal responsiveness to external stimuli.Besides neuronal and gene-regulatory networks, recent empirical data supports that also animalcollectives may be examples of self-organized critical systems. However, open questions about self-organization mechanisms in animal groups remain: Evolutionary adaptation towards group-leveloptima (group-level selection), often implicitly assumed in the “criticality hypothesis”, is in generalnot a reasonable mechanism in fission-fusion groups composed of non-related individuals. Further-more, previous theoretical work relies on non-spatial models, which ignore potentially importantspatial self-organization effects.Here, using a generic, spatially-explicit model of schooling prey being attacked by a predator, weshow first that schools operating at criticality perform best. However, this is not due to optimal re-sponse of the prey to the predator, as suggested by the “criticality hypothesis”, but rather due to thespatial structure of the prey school at criticality. Secondly, by investigating individual-level evolu-tion, we show that the critical point is not an evolutionary stable state. On the contrary, strong spatialself-sorting effects increase the selection strength, and make it an evolutionary unstable state. Ourresults show that in collective behavior spatio-temporal dynamics are important, and that individual-level selection does not in general provide a robust proximate mechanism in support of the “critical-ity hypothesis” in animal groups.
Distributed processing of information is at the core for the function of many complex systems in biology, such as neu-ronal networks [1], genetic regulatory networks [2] or animal collectives [3, 4]. Based on ideas initially developed instatistical physics and theoretical modeling it has been conjectured that such living systems operate in a special param-eter region, in the vicinity of so-called critical points (phase transitions), where the system’s macroscopic dynamicsundergo a qualitative change, and various aspects of collective computation become optimal [5, 6, 7, 8, 9, 10, 11].In recent years some empirical support for the “criticality hypothesis” has been obtained from analysis of neuronaldynamics [12, 10, 13], gene regulatory networks [14, 15], and collective behaviors of animals [16, 17, 18, 19, 20].This evidence is often based on observation of characteristic features of critical behavior, such as power-law distribu-tion or diverging correlation lengths in spatial systems. However these observations could in principle have differentorigins [12, 21, 22, 23]. Therefore, more convincing support for the “criticality hypothesis” can be obtained throughadditional identification of proximate mechanisms enabling biological systems to self-organize towards criticality. Inneuronal systems, synaptic plasticity has been shown to provide such a mechanism [24, 25, 26]. For genetic regulatorynetworks, similar mechanisms based on network rewiring have been proposed [27, 28]. Using an information-theoreticframework Hidalgo et al [11] have shown that (coupled) binary networks evolve towards the critical state in hetero-geneous environments. However, in their model already a single unit (individual) can exhibit a phase-transition andthus tunes itself individually to criticality. In addition, they assumed idealized random interaction networks betweenthe agents. Thus, open questions remain whether evolutionary, individual-level adaptations is a possible self-tuningmechanism for (i) biological collectives, where phase transitions are purely macroscopic phenomena, and (ii) animalgroups characterized by spatial, dynamic interaction networks. In general, if collective computation becomes optimalat a phase transition, a purely macroscopic phenomenon defined only at the group-level, then adaptation based onglobal fitness should be able to tune the system towards criticality. Therefore Darwinian evolution appears a viablemechanism for emergence of self-organized criticality only for complex systems within a single individual, e.g. inthe context of neuronal or genetic networks, or in collectives of closely related individuals such as eusocial insects[16]. In general group-level and individual-level evolutionary optima may be different [29, 30]. Such social dilemmasemerge in broad range of multi-agent evolutionary game theoretic problems, which questions evolution as a generalmechanism for self-tuning of of such systems to criticality, in particular in the context of animal groups, consisting ofweakly or non related individuals.Whereas only few empirical studies report signatures of criticality in collective animal behavior [19, 18, 20], mostsupport for the criticality hypothesis in this context comes from mathematical models. For example, in agent-basedsimulation of fish schools it has been shown that at a critical point the collective state is influenced strongest by single orfew individuals [31], or that collective response to external time-varying signals becomes maximal in idealized latticemodels of flocks [32]. However, dynamical animal groups differ from lattice models [32, 33] due to their dynamicalneighborhood which may induce self-sorting of individuals according to their individual behavioral parameters[34,2
PREPRINT - S
EPTEMBER
7, 202035, 36]. This, in turn, has likely direct evolutionary consequences as for example predators may attack certain swarmregions more frequently[37, 38].Throughout this work, criticality or critical point will refer to the directional order-disorder transition, a prominentphase transition in statistical physics and collective behavior [39]: An initially disordered swarm, where the socialcoordination is weak compared to the noise, shows spontaneous onset of orientational order, where the group startsto move collectively along a common yet ”consensus” direction, if the directional alignment (coupling strength) isincreased beyond a critical parameter (critical point). A further increase of alignment results in highly ordered (polar-ized) schools. This transition is characterized by a so-called spontaneous symmetry breaking: In disordered swarmsthere is no distinguished direction in space. In the ordered state, this symmetry is broken through the emergence of anaverage heading direction of the school, which allows to distinguish front, back and sides of the group.We explore the criticality hypothesis in the context of spatially-explicit predator-prey dynamics, where coordinatedcollective behavior of the prey is believed to entail evolutionary benefits to individuals within the group [40]. Inparticular, we use an agent-based model of grouping and coordinating prey [41, 42, 43, 44, 34, 35], and analyze therole of the spatial structure of the group, its dynamical response and the individual-level selection by applying anevolutionary algorithm [45, 46, 47, 48, 49, 50, 51].We show that the group-level behavior becomes optimal at criticality with respect to two measures: We observei) maximal directional-information transfer between neighbors, and ii) minimal predator capture rates at criticality.However, a detailed analysis reveals that the capture rate, as a relevant measure of evolutionary fitness, becomesminimal only due to the dynamical structure of the collective at criticality, independent on the direct response ofindividuals to the predator, an thus independent of information propagation within the school. Furthermore, throughevolutionary simulations with individual-level selection, we show that the critical point is an evolutionary highlyunstable state. This evolutionary instability can be linked to strong selection due to phenotypic-sorting with respect tothe broken symmetry of the collective state. Finally, the observed evolutionary stable strategies (ESS) appear to resultfrom individual prey agents balancing the influence of social and private information on their movement behavior.
We consider a simple, yet generic agent-based model of schooling prey attacked by a predator. For simplicity weassume that prey agents move with fixed speed v and change their direction according to social forces which ensurethat prey keep a preferred distance to, and align (alignment strength µ alg ) their velocities with their nearest neighbors,defined by a Voronoi tessellation [4, 52]. Randomness in the movement of individuals due to unresolved internaldecisions or environmental noise is modeled as fluctuations in the heading of the agents (angular noise with intensity D ). If a predator is a nearest neighbor of a prey agent, the latter reacts by a repulsive flee-force with strength µ flee .The predator moves with a fixed speed v p which is larger than the preys (here v p = 2 v ) and its direction changestowards the weighted mean direction of its frontal nearest prey, which represent possible targets. The weight corre-sponds to the catch-probability of each target, which decreases linearly with distance until it equals zero at a distancelarger than the catch-radius. If the predator launches an attack, with attack rate γ a , it selects equally likely among thetargets and captures it according to the targets catch-probability. The predator is initiated outside the prey collectivewith a distance slightly above the capture-radius and a with velocity vector oriented towards the center of mass of theprey school.In evolutionary simulations for each generation we perform N r independent runs with different initial conditions for N agents each with a different phenotype (behavioral parameter). Fitness of a prey agent is defined through the negativenumber of deaths of this agent aggregated over N r independent realization. The behavioral phenotypes, i.e. socialforce parameters, of the next generation are selected via fitness-proportionate selection (roulette-wheel-algorithm)with mutations implemented through addition of Gaussian-distributed noise on the selected behavioral parameter. Seemethods for model details. We first investigate whether operating at the order-disorder transition leads to optimal response of the prey school tothe predator. Here, polarization, or normalized average velocity of the group, is the relevant order parameter, whichquantifies the amount of orientational order in the system: For large, fully disordered systems it is close to zero andapproaches 1 in completely ordered systems with all agents moving in the same direction (see methods). It increaseswith the strength of alignment µ alg and decreases with the intensity of angular noise D (Movie S1) in a non-linear3 PREPRINT - S
EPTEMBER
7, 2020Figure 1:
Group optimum:
Predation inde-pendent (
A, B ) and dependent (
C, D ) groupmeasures. ( A ) Polarization Φ . The dashedvertical line marks the angular noise of D =0 . used in the evolution-runs. ( B ) The di-rectional information transfer C ( δ(cid:126)v i , δ(cid:126)v j ) , es-timated via the correlation of velocity fluc-tuations between interacting agents, peaks atthe transition. Inset: Susceptibility, estimatedvia polarization fluctuations. ( C ) Collectiveanti-predator performance is quantified by thecapture rate which is strongly anti-correlatedwith the inter-individual distance R = − . (IID, inset C ). ( D ) Escape ratio R esc = 1 − γ c /γ c,NF . Inset: Difference between capturerates in schools of non-fleeing γ c, NF and fleeing γ c . In all panels: the disorder-order transitionis indicated by a the dash-dotted magenta line.Each parameter point corresponds to an aver-age over N s = 40 simulations, each with N =400 agents attacked for T simu = 120 time unitsafter an equilibration time of T eq = 200 .fashion: It remains small ( ≈ ) throughout most of the disordered regime, before showing the steepest increase inorder in the vicinity of the critical point, and finally asymptotically approaches 1. Both behavioral parameters, µ alg and D can be used as control parameters for crossing of the critical line (diagonal magenta line Fig. 1A) between thedisordered state (low µ alg , high D ) and the ordered state (high µ alg , low D ).A simple and intuitive measure of responsiveness of such a collective system to (local) perturbations is the averagepair-wise correlation of velocity fluctuations C ij = C ( δ(cid:126)v i , δ(cid:126)v j ) between interacting agents (see methods). Here, δ(cid:126)v i = (cid:126)v i − (cid:104) (cid:126)v (cid:105) is the deviation of the velocity of agent i from the average school velocity (cid:104) (cid:126)v (cid:105) . C ij can be interpreted asa simple measure of directional information transfer between neighboring agents i and j : If agent i deviates from theaverage group direction due to a perturbation, large values of C ij indicates that agent j to a large degree is ”copying”this velocity deviation or vice versa.The velocity fluctuation correlation C ij is closely related to the susceptibility χ , which in statistical physics quantifiesthe degree of responsiveness of the system to perturbations, and may become extremely large at criticality. It can bedefined here analogous to magnetic susceptibility in physics [31, 53] (see methods).Both measures, C ij and χ , show a peak at the transition between order and disorder (compare Fig. 1A with B) in linewith predictions of the “criticality hypothesis” [13]. In terms of directional information transfer, i.e. the directionalresponsiveness to perturbation, it appears to be optimal for the collective to operate at criticality. The validity of the above variables from a statistical physics point of view relies on the assumptions of homogeneityand temporal stationarity of the external field, which is not met in our predator-prey scenario: predator perturbationrepresents a strongly local, nonlinear perturbation. As a biologically relevant measure, independent of these assump-tions, we can use the predator capture rate γ c , computed as number of prey captured per time unit. In agreement withthe previous response measures, we find that the capture rate also exhibits a distinct minimum at the critical point (Fig.1C). 4 PREPRINT - S
EPTEMBER
7, 2020However, varying the behavioral parameters of the prey (alignment strength or noise) not only changes the polarizationof the school and the information transfer capability but it also affects the spatial structure of the school (Movies S1,S2), e.g. the average inter-individual distance (IID) or the shape of the school. Our results show that structuralproperties of the prey school correlate strongly with the capture-rate, e.g. the inter-individual distance (inset Fig. 1C)with C ( γ c , IID ) ≈ − . . Thus the reduced capture rate may be potentially related to changes in the structure of theschool at criticality. To distinguish whether structure or susceptibility is responsible for the optimal performance ofthe group at the critical point, we simulated for each predator attack a non-fleeing prey school (flee strength µ flee = 0 )as a control, which is otherwise identical to the responsive school in the remaining parameters and in its positions andvelocities at predator appearance (see Movie S3). The capture-rate of the non-fleeing prey γ c,NF depends only onthe structure of the school. We compare the responsive and control school via two measures: (i) the simple differencebetween both capture rates γ c,NF − γ c and (ii) the escape ratio R esc , which is more robust to fluctuations (see methods)and is defined as the fraction of surviving responsive prey, which would have been captured if they would not flee.Interestingly both measures show no peak at the transition but a continuous increase with alignment strength (Fig.1D) suggesting that the predator-response improves towards the ordered phase if we control for the differences in theself-organized spatial structure (compare column µ alg = 1 with µ alg = 2 in Movie S2).Thus the results suggest that the dynamical structural properties, and not the enhanced responsiveness, are the maincause of the minimum in capture rate, i.e. the group-optimum at criticality. But why is the structure at the transitionso different? The essent´ıal feature of a spatially-explicit system is a continuous co-adaptation between dynamics andstructure, i.e. individual preys movement decisions influences its relative position within the school, which in turn(via modifying the neighborhood and forces) affects its movement. This intrinsic reciprocal influence suggests thatthe highly dynamical structure of the school at the transition is to a large degree caused by the sensitivity to velocityfluctuations (inset in Fig. 1B). Concluding this section, we emphasize that the direct cause of the optimal collectiveperformance (minimal capture rate) is the dynamical spatial structure, as a ”passive” component, and surprisingly notthe maximal responsiveness at criticality. The group-optimum at criticality with respect to prey-survival, does not need to be conincide with the evolutionarystable state (ESS) with respect to evolutionary adaptations at the level individuals. To explore whether the transitionregion is favored by evolution, we let the individual alignment strength µ alg evolve over 500 generations, whilekeeping the angular noise constant ( D = 0 . : vertical line Fig. 1A). In this one-dimensional search-space the systemcan explore all collective states. We repeat the evolutionary simulations from different initial conditions: below( (cid:104) µ alg (cid:105) = 0 ), above ( (cid:104) µ alg (cid:105) = 5 ) and far above ( (cid:104) µ alg (cid:105) = 10 ) the transition ( µ c,alg ≈ . ). To ensure that theevolution ends at the ESS we compute the fitness gradient which represent the strength of the selection pressure at aspecific mean alignment strength (see methods). Assuming a monomodal phenotype distribution, as observed in ourevolutionary runs, a change in sign of the fitness gradient marks the location of the ESS. All three initiations end inthe ordered region far above the critical point (Fig. 2A) and fluctuate around ESS ( µ alg ) ≈ . (vertical dashed lineFig. 2D). Thus, the transition region is not an attractor of the evolutionary dynamics. On the contrary, it is an highlyunstable point with fast evolutionary dynamics due to particularly strong selection pressure at criticality. The fitnessgradient peaks shortly above the transition in the ordered phase (Fig. 2D), with evolutionary dynamics pushing thesystem out of the transition region towards stronger alignment.A possible driver of this maximal selection pressure is self-sorting, i.e. the tendency of individuals to sort according totheir behavioral parameters along specific spatial dimensions of the school, e.g. front-back or side-center, or in regionsof higher or lower density (Fig. 2C) [34]. We quantify this self-sorting through the Pearson correlation coefficientbetween the individuals alignment strength and its relative position/property (see methods). A more advanced measurefor self-sorting is assortative mixing of the school (assortativity coefficient, see methods). Assortativity (Fig. 2D) aswell as other self-sorting measures (Fig. 2B) exhibit extrema which coincide with the fitness gradient peak. Notethat a strong assortative mixing is equivalent to the formation of spatially coherent subpopulation within the schoolexhibiting similar behavioral parameter. In this context a peak in fitness gradient close to transition means that subpop-ulations with stronger alignment, thus better directional coordination, actively or passively perform better at avoidingcapture. An increase of the escape ratio R esc with increasing alignment close to criticality (see Fig. 2D) suggestan enhanced active avoidance. However, also passive effects appear to play an important role since the correlationbetween the fitness of a prey and its relative position becomes maximal in the same parameter region (Fig. 2E). Onespecific mechanism of passive avoidance is the dilution effect[40] caused by local density differences correlating withbehavioral phenotypes. Stronger aligning individuals form denser subgroups within the prey school (density-sortingFig. 2B). As a consequence they have systematically smaller domain of dangers [54] and are thus less frequentlyattacked by the predator. 5 PREPRINT - S
EPTEMBER
7, 2020Figure 2:
Evolution under predation: ( A ) Overlay of three evolutionary runs starting at (cid:104) µ alg (cid:105) = [0 , , over1000 generations. The behavioral phenotype is determined only by the alignment strength as the evolving parameter.The predator attacks from random initial directions for T simu = 120 . The inset shows the evolution of the populationmean alignment parameter (cid:104) µ alg (cid:105) of the three different evolutionary runs. ( B ) Assortativity coefficient (blue line)and smoothed fitness gradient ∇ f (red line). The evolutionary stable state is defined by the zero crossing of thefitness gradient and represented as a vertical dashed black line. Black dots are the non-averaged fitness gradients foreach generation (see methods). ( C ) Self-sorting measured as correlation C ( µ alg , x ) between the individual alignmentstrength µ alg and average spatial property of the indivdual as front-(red) and side-position(black) and local density(blue). ( D ) Correlation C ( f, x ) of individual fitness with the, latter mentioned, average relative spatial positions.( E ) Simulation snapshot illustrating the front-(red) and side-position(black) and local density (blue). In all panels: thevertical dash-dotted magenta line marks the order-disorder transition and the vertical dashed black line the evolutionarystable state.It is possible to disentangle passive, structural effects from an active response, by setting the flee-strength to zero. Thisresults in a significantly smaller, yet finite, fitness-gradient-peak at the transition (Fig. S.2H). This suggests that both,the structural, passive selection as well as the different active avoidance behavior of different phenotypes contribute tothe strong selection pressure at criticality.We note that the sudden increase of self-sorting at the transition is due to a coupled symmetry breaking. At the order-disorder transition the directional symmetry is broken and the school ”agrees” on a common movement direction.This also breaks the symmetry between relative locations within the school. For example in the disordered phaseevery edge position is equivalent, but with the emergence of the common movement direction the sides and rear ofthe school become structurally different from the front. This can be clearly seen in the comparison of the fitnesscorrelations with respect to specific relative spatial positions within the school (”side-sorting” versus ”front-sorting”):Below the transition the corresponding curves become indistinguishable, whereas above at the transition they start todeviate and show different behavior with increasing alignment strength (Fig. 2B). Despite the importance of self-sorting for the maximal selection pressure at the transition, it does not provide anyexplanation for the observed location of the ESS. The spatial properties can not explain the negative fitness gradientfor µ alg > ESS ( µ alg ) ≈ . . In this regime either the self-sorting is negligible, as for side- and density-sorting (Fig.2B), or the relative location has no effect on the individual fitness, as observed along the front-back dimension (Fig.2E). If the ESS is not determined by the structural self-organization of the school, it has to originate from individuallyoptimal predator evasion. In this case the ESS has to depend on the flee-strength µ flee as the main parameter tuningthe strength of individual predator response.We do find clear dependence of the ESS on the flee-strength (Fig. 3A), more specifically the ESS exhibits a lineardependence on the flee-strength for µ flee ≥ (diagonal line in Fig. 3B). The order transition acts as a lower bound6 PREPRINT - S
EPTEMBER
7, 2020Figure 3:
Evolution for different fleestrengths. ( A ) Evolution trajectory of the meanalignment strength µ alg over 700 generations.( B ) shows a dependence of evolutionary stablestrategies (ESS) on flee strength. Solid diago-nal line shows the theoretically predicted lin-ear dependence of the ESS on µ flee assum-ing balancing of social and private information(see SI Sect. C). Dashed lines ( A , B ) con-nect the individual evolutionary runs ( A ) to thecorresponding ESSs ( B ) obtained from multi-ple, longer evolutionary simulations. ( C ) Evo-lutionary stable states (circles) with respect tothe group response, measured via the escaperatio R esc , for three selected flee-strengths in-dicated with dashed, solid and dotted lines for µ flee = [2 , , respectively. In all panels:the dash-dotted magenta line marks the order-disorder transition and the different lines (red,black and blue) represent results for differentflee strengths µ flee = [2 , , , respectively.since the non-fleeing agents ( µ flee = 0 ) equilibrate closely above it. Thus the ESS for non-responding agents matcheshe group-level optimum due to the dynamical school structure at criticality.The linear dependence on the flee-strength may be explained by prey balancing social vs. personal predator informa-tion. Social information about the predator is beneficial if the prey is in the second neighbour shell of the predator,i.e. where its neighbors but not itself responses directly to the predator. Thus by coordinating with its informedneighbours it gains distance to the predator. However, if a prey directly senses the predator, social information ofuninformed neighbors conflicts with its private information and therefore may hinder optimal evasion.Therefore, individual prey agents should continue to evolve towards stronger alignment strength until costs of thesocial inhibition of evasion counterbalance the benefit of social information. We find support for this conjecture byreproducing the observed linear dependence through a local mean-field approximation (see SI Sect. A.3) assumingthe above balancing mechanism (Fig. 3B). Interestingly, also the escape ratio, as a measure of group response whilecontrolling against spatial effects, exhibits a maximum in the strongly ordered region away from criticality (Fig. 1D).This leads to the question whether the ESS coincide with the strongest group response. Indeed, the maximum ofescape ratio shows the same trend as the ESS of moving towards higher alignment strenghts with increasing fleestrangth (Fig. 3C), but these maxima stay clearly below the corresponding ESSs (circles in Fig. 3C). This suggeststhat the system does evolve towards unresponsiveness [29] by increasing the social responsiveness above the optimum(compare column µ alg = 2 with µ alg = 4 in Movie S4).The qualitative results are independent of implementation details. We checked for robustness against the predatorattack scheme (more and less agile predator), prey-modification (variable speed, persistence length), modifications inevolutionary algorithm (attack-rate, fitness-estimation) and importantly in a heterogeneous environment (see SI Sect.A.4 and Figs. S.4, S.6). Only by introducing an additional selection pressure, creating a heterogeneous environment,which favors disordered shoals and increasing its weight the ESS may be shifted into the disordered phase. However,even in this case the critical point acts as an unstable evolutionary point and an attractor (Fig. S.6). We have shown, using a spatially-explicit agent-based model of predator-prey dynamics, that the group optimum withrespect to predation avoidance is located in the vicinity of the critical point between disordered swarming and orderedschooling, in line with the so-called “criticality hypothesis”. However, this optimality is not due to optimal transfer7
PREPRINT - S
EPTEMBER
7, 2020of social information but rather due to the highly dynamical structure of the group at the transition. Yet, this groupoptimum at criticality does not represent an evolutionary stable state of individual-level selection.Our work demonstrates the crucial importance of taking into account the self-organized spatial dynamics of animalgroups when evaluating potential evolutionary benefits of grouping. It turns out that the mechanism responsiblefor the optimal collective performance (minimal capture rate) at the critical point, the highly dynamic and flexiblestructure of the collective, leads also to the steepest selection gradients in evolutionary dynamics, making the criticalpoint evolutionary unstable. Evolution with random mutations enforces heterogeneity which. in combination withthe spatial symmetry breaking at the transition, results in maximal assortative mixing and self-sorting close to thetransition. These effects of self-organized collective behavior play a decisive role for the evolutionary dynamics closeto criticality and “drive” the ESS out of the transition region towards the aligned state. In our system the ESS is inthe strongly ordered phase, which suggests the evolution towards external unresponsiveness by overestimating socialinformation. Finally, we show that the ESS depends linearly on the flee strength, i.e. local perturbation strength, whichcan be explained by individual balancing of benefits of social information about the predators approach with the costsof social interactions if the information is directly available.In contrast to Hidalgo et. al. [11], the critical state in our model is not evolutionary stable, despite the similar setup:evolving agents which respond to conspecifics and to a changing environment (here the appearance of a predator).This can be explained by crucial differences to our work. Most importantly, in [11] each agent in isolation can alreadyevolve to its “individual” transition by tuning its own gene regulatory network. This appears to be essential withintheir information-based fitness framework that the a critical point is also an evolutionary stable state. In our model,the disorder-order transition is a pure collective effect, i.e. individual agents cannot exhibit any transition behaviorby themselves. Furthermore, at the disorder-order transition, small differences in behavioral parameters translate intosystematic differences in the self-organized spatial positioning within the group, which in turn directly impacts thepredation threat. This self-sorting[34, 35, 36] is maximal just above the transition and includes assortative mixing andtherefore subpopulations which differ systematically in their behavioral parameters, possibly interesting for collectivetask distribution and computation.There is another consequence of the tight coupling between local school structure and individual dynamics: The extentof the collective is largest at the transition because the responsiveness to directional fluctuations is maximal, i.e. differ-ent parts of the school respond to uncorrelated local fluctuations causing the school effectively to expand. In systemswith a one-way influence from structure to dynamics (fixed networks) it is known that at the order-transition structuraldifferences cause the largest dynamic variability [55]. We show here that in a system with additional feedback fromthe dynamics to the structure, also the structure has the highest variability at the transition, which may have importantconsequences for collective computations, as it may for example enhance collective gradient sensing[56, 48]. It showsthat interactions on fixed [32, 30, 33] or randomly rewiring [29] lattices might miss a functionally highly relevantfeatures of collective behavior.For simplicity we used here an established social force model [34, 35, 41], which also received empirical support [42,43]. Alternative social interaction mechanism may also lead to highly ordered collective movement [57, 58, 59, 60],but the crucial spatial effects such as self-sorting and assortative mixing are generic and will be independent on thedetails of social interactions.Our findings, individual-level adaptation does not evolve the prey groups towards criticality, suggests that evolutionaryadaptations may not be viewed as a general mechanism for self-organization towards criticality. Whereas it does notexclude the possibility that animal collectives may operate in the vicinity of phase transitions in order to optimizecollective computations, it clearly demonstrates the necessity for further research on plausible proximate mechanimsof self-organized criticality in animal groups. A general, fundamental difficulty is that besides predator evasion thereare various ecological contexts and other dimensions of (collective) behavior which affect individual fitness. Here, byfocusing on a dominant selection pressure, namely predating, we neglect other mechanisms, as for example resourceexploration and exploitation [48, 30, 45] whose ESS can also depend on the resource abundance [30, 45]. Thisemphasizes the importance to study collective behavior in the wild [61, 62, 63, 64] to provide more empirical input onactual relevant behavioral mechanisms as well as variability of behavior across different contexts. However, we haveshown that even by combining two opposing selection mechanisms (see SI Sect. A.4.3), which on their own favorordered or disordered state respectively, the critical point does not correspond to an evolutionary attractor, it remainsan evolutionary highly unstable point.We focused here on the prominent directional symmetry breaking transition, however recently it was suggest that atransition in the speed relaxation coefficient may represent an functionally relevant critical point in flocking behavior[19]. Individuals with lower relaxation constants are less bound to their preferred speed and may gain fitness benefitsdue their ability to adapt faster to higher speeds of fleeing conspecifics. Consistent with this hypothesis, guppies(
Poecilia reticulata ) exhibit stronger accelerations in high-predation habitats[44]. However, fish exhibit a reflex-driven8
PREPRINT - S
EPTEMBER
7, 2020escape response, so-called startle, which was recently shown to spreads through fish schools as a behavioral contagionprocess [65, 66]. This suggests that at least in the context of collective predator evasion in fish, another type ofcritical points may be highly relevant, which is similar to the critical threshold in epidemic models. It separates statesof non-propagating startle response, with only small localized response of single or few individuals, and avalanche-like dynamics, where a single fish may cause a global startle cascade. Overall, our study does not reject the generalpossibility that animal groups manifest critical behavior and that it may be adaptive. However, it highlights importanceof identification of biologically plausible proximate mechanisms for self-organization towards - and maintenance of- critical dynamics in animal groups, which account for spatial self-organization and the corresponding ecologicalniche.
Methods
All Model parameters are listed in Tab.S1.
Prey model
A prey agent i moves in 2D with constant velocity v = v with directional noise, parametrized by the angular diffusion coefficient D [67] and responds to a combined force (cid:126)F i = (cid:126)F i,alg + (cid:126)F i,d + (cid:126)F i,flee by adapting its position (cid:126)r i and heading ϕ i as d(cid:126)r i ( t ) dt = (cid:126)v i ( t ) (1a) dϕ i ( t ) dt = 1 v (cid:16) F i,ϕ ( t ) + √ Dξ ( t ) (cid:17) (1b)with F i,ϕ ( t ) = (cid:126)F i ( t ) · (cid:126)e ϕ as the combined force along the agents angular direction (cid:126)e ϕ and ξ ( t ) as Gaussian white noise. Thealignment force ( (cid:126)F i,alg ) between a focal agent i and all its neighbors j ∈ N i is the averaged velocity difference (cid:126)v ji = (cid:126)v j − (cid:126)v i timesthe alignment strength µ alg . The distance regulating force (see Fig. S.1 A) is (cid:126)F i,d = 1 | N i | (cid:88) j ∈ N i µ d · tanh ( m d ( r ji − r d )) · ˆ r ji (2)with ˆ r ji = ( (cid:126)r j − (cid:126)r i ) / | (cid:126)r j − (cid:126)r i | as direction from agent i to j , r d as preferred distance, µ d as strength of the force and m d as theslope of the change from repulsion (for r ji < r d ) to attraction (for r ji > r d ). If a predator p is a neighbor, the agent is repelled( (cid:126)F i,flee ) from it with a flee strength µ flee . Predator-model
The predator moves with fixed speed v p = 2 v according to dϕ p dt = 1 v p (cid:126)e ϕ · (cid:126)F p (3)with (cid:126)F p as the pursuit force. It considers its frontal Voronoi-neighbors N p as targets and selects equally likely among them( p select,i = 1 / | N p | ). It only attacks one prey at a time. If the predator launches an attack, with an attack rate γ a (also accountingfor handling time), its success probability decreases linear with distance: p success,i = ( r catch − r ip ) /r catch (4)and is zero for distances larger than r catch . In summary, the probability that a predator successfully catches a targeted agent withina small time window [ t, t + δt ] is p catch,i ( t, δt ) = p success,i ( t ) p select,i ( t ) γ a δt. (5)The pursuit force, with constant magnitude µ p , points to a weighted center of mass. Each prey position is weighted by its probabilityof a successful catch p catch,i ( t, δt ) . Evolutionary algorithm
The algorithm consists of three components: fitness estimation, fitness-proportionate-selection and mutation.(i) The fitness is estimated by running N f = 76 independent attack-simulations on the same prey population. For each simulationthe γ a · T s agents with the largest cumulative p catch are declared as dead. The fitness of agent i is f i = − N k,i + max ( N k,j , j ) with N k,i as the number of simulations in which agent i was captured and max ( N k,j , j ) is the largest number of deaths among allagents.(ii) The new generation of N offspring is generated via fitness-proportionate-selection. Thus a random offspring has the parameters PREPRINT - S
EPTEMBER
7, 2020 of the parent i with probability p parent,i = f i / (cid:80) j f j .(iii) An offspring mutates with probability γ m (mutation rate), by adding a Gaussian random variable with zero mean and standarddeviation σ m to its alignment strength µ alg .Steps (i) till (iii) are repeated in each generation. To estimate the ESS we compute for each generation the expected offspringpopulation (without mutation to reduce noise) and define the fitness gradient as the offspring mean parameter from which thecurrent mean parameter is subtracted. Thus if the offspring have a larger mean parameter, the fitness gradient is positive and viceversa. The mean fitness gradient of a certain parameter region is the average of generations within it. For details see SI Sect. A.1.4. Quantification of collective behavior
The inter-individual distance is the distance between prey pairs averaged over all pairs
IID = (cid:104)| (cid:126)r ij |(cid:105) . The polarization Φ is theabsolute value of the mean heading direction Φ = | (cid:126) Φ | = | (cid:80) i (cid:126)u i /N | . The susceptibility χ is the response of the polarization to anexternal field h and can be measured via polarization fluctuations χ = ∂ Φ ∂h = N ( (cid:104) Φ (cid:105) − (cid:104) Φ (cid:105) ) (6)which is a form of the fluctuation dissipation theorem (see SI Sect. A.2.4). It can be shown that Eq. 6 is the same as the correlationof velocity fluctuations δ(cid:126)v i = (cid:126)v i − (cid:126) Φ over all possible pairs (see SI Sect. A.2.4). However, in inset of Fig. 1B we computed thecorrelation of velocity fluctuations only over neighboring pairs C ( δ(cid:126)v i , δ(cid:126)v j ) = (cid:80) i,j ∈ N i δ(cid:126)v i · δ(cid:126)v j because it is directly related tolocal transfer of social information than the correlation over all, including totally unrelated, prey pairs.We compare the performance of the fleeing prey to the non-fleeing prey (control) using escape ratio R esc = 1 − γ c γ c,NF . (7)It is equal to the difference between the capture rates of non-fleeing and fleeing agents γ c,NF − γ c scaled by γ c,NF . The normal-ization of the capture difference by the baseline capture rate of non-fleeing prey γ c,NF accounts for potential differences in capturerates due to differences in school structure for different parameters, which are unrelated to the fleeing response.The self-sorting is quantified via the Pearson correlation coefficient between the alignment parameter µ i,alg of individual agents andtheir mean relative location in the collective (cid:104) r i,x (cid:105) where x ∈ { f, s, d } which stands for front, side and local density respectively.Agents at the front (back) have the largest (smallest) front-location and at the side (center) have the largest (smallest) side-location.The local density sorting is the correlation of the agents local density and its alignment strength. For the detailed computation of therelative locations see SI Sect. A.2.1. Another, more general, quantification of self-sorting is how assortative the spatial arrangementof individuals with heterogeneous alignment is. We used the implementation of the assortativity coefficient[68] in igraph on theinteraction network (Voronoi) with the values for each agent corresponding to their alignment strength (see SI Sect. A.2 for details). Codes availability
The code to run the predator prey model is available at github ( https://github.com/PaPeK/PredatorPrey ). Acknowledgments
We are grateful to Iain Couzin, Simon Levin, Jessica Flack, and David Krakauer for insightful andinspiring discussions leading to this work. We further acknowledge Ishan Levy for suggesting the investigation of non-respondingprey dynamics. We acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) underGermany’s Excellence Strategy – EXC 2002/1 “Science of Intelligence” – project number 390523135 and DFG Emmy NoetherProgramme - RO 4766/2-1.
References [1] Mikail Rubinov and Olaf Sporns. Complex network measures of brain connectivity: Uses and interpretations.
NeuroImage , 52(3):1059–1069, sep 2010.[2] Nedumparambathmarath Vijesh, Swarup Kumar Chakrabarti, and Janardanan Sreekumar. Modeling of generegulatory networks: A review.
Journal of Biomedical Science and Engineering , 06(02):223–231, 2013.[3] N. Miller, S. Garnier, A. T. Hartnett, and I. D. Couzin. Both information and social cohesion determine collectivedecisions in animal groups.
Proceedings of the National Academy of Sciences , 110(13):5263–5268, mar 2013.[4] Ariana Strandburg-Peshkin, Colin R. Twomey, Nikolai W.F. F Bode, Albert B. Kao, Yael Katz, Christos C.Ioannou, Sara B. Rosenthal, Colin J. Torney, Hai Shan Wu, Simon A. Levin, and Iain D. Couzin. Visual sensorynetworks and effective information transfer in animal groups.
Current Biology , 23(17):R709–R711, 2013.[5] N. H. Packard. Adaptation Toward the Edge of Chaos. In J.A.S. Kelso, A.J. Mandell, and M.F. Shlesinger,editors,
Dynamic Patterns in Complex Systems . Singapore, World Scientific, 1988.10
PREPRINT - S
EPTEMBER
7, 2020[6] Per Bak, Kan Chen, and Michael Creutz. Self-organized criticality in the ’Game of Life’.
Nature , 342:780–782,1989.[7] C. G. Langton. Computation at the edge of chaos: Phase transitions and emergent computation.
Physica D ,42:12– 37, 1990.[8] Per Bak and Kim Sneppen. Punctuated Equilibribum and Criticality in a simple model of evolution.
PhysicalReview Letters , 71(24):4083–4086, 1993.[9] Osame Kinouchi and Mauro Copelli. Optimal dynamical range of excitable networks at criticality.
NaturePhysics , 2(5):348–351, may 2006.[10] Thierry Mora and William Bialek. Are Biological Systems Poised at Criticality?
Journal of Statistical Physics ,144(2):268–302, 2011.[11] Jorge Hidalgo, Jacopo Grilli, Samir Suweis, Miguel A. Munoz, Jayanth R. Banavar, and Amos Maritan.Information-based fitness and the emergence of criticality in living systems.
Proceedings of the NationalAcademy of Sciences , 111(28):10095–10100, jul 2014.[12] John M. Beggs and Nicholas Timme. Being critical of criticality in the brain.
Frontiers in Physiology , 3JUN(June):1–14, 2012.[13] Miguel A. Mu˜noz. Colloquium: Criticality and dynamical scaling in living systems.
Reviews of Modern Physics ,90(3):31001, 2018.[14] Dmitry Krotov, Julien O Dubuis, Thomas Gregor, and William Bialek. Morphogenesis at criticality.
Proceedingsof the National Academy of Sciences , 111(10):3683–3688, 2014.[15] Bryan C Daniels, Hyunju Kim, Douglas Moore, Siyu Zhou, Harrison B Smith, Bradley Karas, Stuart A Kauff-man, and Sara I Walker. Criticality distinguishes the ensemble of biological regulatory networks.
Physical reviewletters , 121(13):138102, 2018.[16] Aviram Gelblum, Itai Pinkoviezky, Ehud Fonio, Abhijit Ghosh, Nir Gov, and Ofer Feinerman. Ant groupsoptimally amplify the effect of transiently informed individuals.
Nature Communications , 6, 2015.[17] Ofer Feinerman, Itai Pinkoviezky, Aviram Gelblum, Ehud Fonio, and Nir S. Gov. The physics of cooperativetransport in groups of ants.
Nature Physics , 14(7):1–11, jul 2018.[18] Alessandro Attanasi, Andrea Cavagna, Lorenzo Del Castello, Irene Giardina, Stefania Melillo, Leonardo Parisi,Oliver Pohl, Bruno Rossaro, Edward Shen, Edmondo Silvestri, and Massimiliano Viale. Finite-size scaling as away to probe near-criticality in natural swarms.
Physical Review Letters , 113(23):238102, dec 2014.[19] William Bialek, Andrea Cavagna, Irene Giardina, Thierry Mora, Oliver Pohl, Edmondo Silvestri, MassimilianoViale, and Aleksandra M. Walczak. Social interactions dominate speed control in poising natural flocks nearcriticality.
Proceedings of the National Academy of Sciences of the United States of America , 111(20):7212–7217, may 2014.[20] Bryan C. Daniels, David C. Krakauer, and Jessica C. Flack. Control of finite critical behaviour in a small-scalesocial system.
Nature Communications , 8:1–8, 2017.[21] Tom Lorimer, Florian Gomez, and Ruedi Stoop. Two universal physical principles shape the power-law statisticsof real-world networks.
Sci. Rep. , 5:1–8, 2015.[22] William J. Reed and Barry D. Hughes. From gene families and genera to incomes and internet file sizes: Whypower laws are so common in nature.
Phys. Rev. E - Stat. Physics, Plasmas, Fluids, Relat. Interdiscip. Top. ,66(6):4, 2002.[23] Jonathan Touboul and Alain Destexhe. Power-law statistics and universal scaling in the absence of criticality.
Physical Review E , 95(1):012413, 2017.[24] Stefan Bornholdt and Thimo Rohlf. Topological evolution of dynamical networks: Global criticality from localdynamics.
Physical Review Letters , 84(26):6114, 2000.[25] Christian Meisel and Thilo Gross. Adaptive self-organization in a realistic neural network model.
PhysicalReview E , 80(6):061917, 2009.[26] Zhengyu Ma, Gina G Turrigiano, Ralf Wessel, and Keith B Hengen. Cortical circuit dynamics are homeostati-cally tuned to criticality in vivo.
Neuron , 104(4):655–664, 2019.[27] Min Liu and Kevin E Bassler. Emergent criticality from coevolution in random boolean networks.
PhysicalReview E , 74(4):041910, 2006.[28] Ben D MacArthur, Rub´en J S´anchez-Garc´ıa, and Avi Ma’ayan. Microdynamics and criticality of adaptive regu-latory networks.
Physical review letters , 104(16):168701, 2010.11
PREPRINT - S
EPTEMBER
7, 2020[29] Colin J. Torney, Tommaso Lorenzi, Iain D. Couzin, and Simon A. Levin. Social information use and the evolutionof unresponsiveness in collective systems.
Journal of the Royal Society Interface , 12(103), 2015.[30] Eleanor Redstart Brush, Naomi Ehrich Leonard, and Simon A. Levin. The content and availability of informationaffects the evolution of social-information gathering strategies.
Theoretical Ecology , 9(4):455–476, dec 2016.[31] Daniel S Calovi, Ugo Lopez, Paul Schuhmacher, Hugues Chate, C. Sire, and Guy Theraulaz. Collective responseto perturbations in a data-driven fish school model.
Journal of The Royal Society Interface , 12(104):20141362–20141362, jan 2015.[32] Fabio Vanni, Mirko Lukovi´c, and Paolo Grigolini. Criticality and Transmission of Information in a Swarm ofCooperative Units.
Physical Review Letters , 107(7):078103, aug 2011.[33] Amanda Chicoli and Derek A. Paley. Probabilistic information transmission in a network of coupled oscillatorsreveals speed-accuracy trade-off in responding to threats.
Chaos An Interdiscip. J. Nonlinear Sci. , 26(11):116311,2016.[34] Iain D Couzin, Jens Krause, Richard James, Graeme D Ruxton, and Nigel R Franks. Collective memory andspatial sorting in animal groups.
Journal of theoretical biology , 218(1):1–11, 2002.[35] Charlotte K. Hemelrijk and Hanspeter Kunz. Density distribution and size sorting in fish schools: An individual-based model.
Behavioral Ecology , 16(1):178–187, 2005.[36] A. Jamie Wood. Strategy selection under predation; evolutionary analysis of the emergence of cohesive aggre-gations.
Journal of Theoretical Biology , 264(4):1102–1110, 2010.[37] JENS KRAUSE. Differential Fitness Returns in Relation To Spatial Position in Groups.
Biological Reviews ,69(2):187–206, 1994.[38] Dirk Bumann, Dan Rubenstein, and Jens Krause. Mortality Risk of Spatial Positions in Animal Groups: theDanger of Being in the Front.
Behaviour , 134(13-14):1063–1076, 1997.[39] Tams Vicsek, Andrs Czirk, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet. Novel type of phase transition in asystem of self-driven particles.
Physical Review Letters , 75(6):1226–1229, aug 1995.[40] Jens Krause and Graeme D. Ruxton.
Living in Groups . Oxford Univ. Press, Oxford, 2002.[41] Roy Harpaz, Gaˇsper Tkaˇcik, and Elad Schneidman. Discrete modes of social information processing predictindividual behavior of fish in a group.
Proceedings of the National Academy of Sciences , 114(38):10149–10154,sep 2017.[42] Yael Katz, Kolbjørn Tunstrøm, Christos C. Ioannou, Cristi´an Huepe, and Iain D. Couzin. Inferring the structureand dynamics of interactions in schooling fish.
Proceedings of the National Academy of Sciences of the UnitedStates of America , 108(46):18720–18725, 2011.[43] Daniel S. Calovi, Alexandra Litchinko, Valentin Lecheval, Ugo Lopez, Alfonso P´erez Escudero, Hugues Chat´e,Cl´ement Sire, and Guy Theraulaz. Disentangling and modeling interactions in fish with burst-and-coast swim-ming reveal distinct alignment and attraction behaviors.
PLOS Computational Biology , 14(1):e1005933, jan2018.[44] James E. Herbert-Read, Emil Ros´en, Alex Szorkovszky, Christos C. Ioannou, Bj¨orn Rogell, Andrea Perna, In-dar W. Ramnarine, Alexander Kotrschal, Niclas Kolm, Jens Krause, and David J. T. Sumpter. How predationshapes the social interaction rules of shoaling fish.
Proceedings of the Royal Society B: Biological Sciences ,284(1861):20171126, aug 2017.[45] Andrew J. Wood and Graeme J. Ackland. Evolving the selfish herd: Emergence of distinct aggregating strategiesin an individual-based model.
Proceedings of the Royal Society B: Biological Sciences , 274(1618):1637–1642,jul 2007.[46] Randal S. Olson, Arend Hintze, Fred C. Dyer, David B. Knoester, and Christoph Adami. Predator confusion issufficient to evolve swarming behaviour.
Journal of The Royal Society Interface , 10(85):20130305–20130305,jun 2013.[47] Randal S. Olson, David B. Knoester, and Christoph Adami. Evolution of Swarming Behavior Is Shaped by HowPredators Attack.
Artificial Life , 22(3):299–318, aug 2016.[48] Andrew M. Hein, Sara Brin Rosenthal, George I. Hagstrom, Andrew Berdahl, Colin J. Torney, and Iain D.Couzin. The evolution of distributed sensing and collective computation in animal populations. eLife ,4(DECEMBER2015):1–43, 2015.[49] Vishwesha Guttal and Iain D. Couzin. Social interactions, information use, and the evolution of collectivemigration.
Proceedings of the National Academy of Sciences of the United States of America , 107(37):16172–16177, 2010. 12
PREPRINT - S
EPTEMBER
7, 2020[50] Christopher T. Monk, Matthieu Barbier, Pawel Romanczuk, James R. Watson, Josep Al´os, ShinnosukeNakayama, Daniel I. Rubenstein, Simon A. Levin, and Robert Arlinghaus. How ecology shapes exploitation:a framework to predict the behavioural response of human and animal foragers along exploration–exploitationtrade-offs.
Ecology Letters , 21(6):779–793, 2018.[51] Daniel J. Van Der Post, Rineke Verbrugge, and Charlotte K. Hemelrijk. The evolution of different forms ofsociality: Behavioral mechanisms and eco-evolutionary feedback.
PLoS ONE , 10(1):1–19, 2015.[52] M Ballerini, N Cabibbo, R Candelier, A Cavagna, E Cisbani, I Giardina, V Lecomte, A Orlandi, G Parisi,A Procaccini, M Viale, and V Zdravkovic. Interaction ruling animal collective behavior depends on topologicalrather than metric distance: Evidence from a field study.
Proceedings of the National Academy of Sciences ,105(4):1232–1237, 2008.[53] Andrea Cavagna, Irene Giardina, and Tom´as S. Grigera. The physics of flocking: Correlation as a compass fromexperiments to theory.
Phys. Rep. , 728:1–62, jan 2018.[54] R James, PG Bennett, and J Krause. Geometry for mutualistic and selfish herds: the limited domain of danger.
Journal of Theoretical Biology , 228(1):107–113, 2004.[55] Matti Nykter, Nathan D. Price, Antti Larjo, Tommi Aho, Stuart A. Kauffman, Olli Yli-Harja, and Ilya Shmule-vich. Critical networks exhibit maximal information diversity in structure-dynamics relationships.
Phys. Rev.Lett. , 100(5):1–4, 2008.[56] Andrew Berdahl, Colin J. Torney, Christos C. Ioannou, Jolyon J. Faria, and Iain D. Couzin. Emergent Sensingof Complex Environments by Mobile Animal Groups.
Science , 339(6119):574–576, feb 2013.[57] Renaud Bastien and Pawel Romanczuk. A model of collective behavior based purely on vision.
Sci. Adv. ,6(6):1–10, 2020.[58] Pawel Romanczuk, Iain D. Couzin, and Lutz Schimansky-Geier. Collective Motion due to Individual Escape andPursuit Response.
Physical Review Letters , 102(1):010602, jan 2009.[59] Jitesh Jhawar, Richard G. Morris, U. R. Amith-Kumar, M. Danny Raj, Tim Rogers, Harikrishnan Rajendran, andVishwesha Guttal. Noise-induced schooling of fish.
Nature Physics , 16(4):488–493, 2020.[60] Liu Lei, Ram´on Escobedo, Cl´ement Sire, and Guy Theraulaz. Computational and robotic modeling revealparsimonious combinations of interactions between individuals in schooling fish.
PLOS Computational Biology ,16(3):e1007194, mar 2020.[61] Nils Olav Handegard, Kevin M. Boswell, Christos C. Ioannou, Simon P. Leblanc, Dag B. Tjøstheim, and Iain D.Couzin. The Dynamics of Coordinated Group Hunting and Collective Information Transfer among SchoolingPrey.
Curr. Biol. , 22(13):1213–1217, jul 2012.[62] Fritz A Francisco, Paul N¨uhrenberg, and Alex Jordan. High-resolution, non-invasive animal trackingand reconstruction of local environment in aquatic ecosystems.
Movement Ecology (forthcoming) , page2020.02.25.963926, 2020.[63] M. J. Hansen, S. Krause, M. Breuker, R. H.J.M. Kurvers, F. Dhellemmes, P. E. Viblanc, J. M¨uller, C. Mahlow,K. Boswell, S. Marras, P. Domenici, A. D.M. Wilson, J. E. Herbert-Read, J. F. Steffensen, G. Fritsch, T. B.Hildebrandt, P. Zaslansky, P. Bach, P. S. Sabarros, and J. Krause. Linking hunting weaponry to attack strategiesin sailfish and striped marlin.
Proceedings of the Royal Society B: Biological Sciences , 287(1918), 2020.[64] Jacob M. Graving, Daniel Chae, Hemal Naik, Liang Li, Benjamin Koger, Blair R. Costelloe, and Iain D. Couzin.Deepposekit, a software toolkit for fast and robust animal pose estimation using deep learning. eLife , 8:1–42,2019.[65] Matthew M.G. Sosna, Colin R. Twomey, Joseph Bak-Coleman, Winnie Poel, Bryan C. Daniels, Pawel Ro-manczuk, and Iain D. Couzin. Individual and collective encoding of risk in animal groups.
Proceedings of theNational Academy of Sciences of the United States of America , 116(41):20556–20561, 2019.[66] Sara Brin Rosenthal, Colin R. Twomey, Andrew T. Hartnett, Hai Shan Wu, and Iain D. Couzin. Revealingthe hidden networks of interaction in mobile animal groups allows prediction of complex behavioral contagion.
Proceedings of the National Academy of Sciences of the United States of America , 112(15):4690–4695, 2015.[67] Pawel Romanczuk and Lutz Schimansky-Geier. Brownian Motion with Active Fluctuations.
Physical ReviewLetters , 106(23):230601, jun 2011.[68] M. E. J. Newman. Mixing patterns in networks.
Physical Review E , 67(2):026126, sep 2002.[69] Umberto Marini Bettolo Marconi, Andrea Puglisi, Lamberto Rondoni, and Angelo Vulpiani. Fluctuation-dissipation: Response theory in statistical physics.
Physics Reports , 461(4-6):111–195, 2008.13
PREPRINT - S
EPTEMBER
7, 2020[70] Robert Großmann, Lutz Schimansky-Geier, and Pawel Romanczuk. Active Brownian particles with velocity-alignment and active fluctuations.
New Journal of Physics , 14, apr 2012.[71] B. H. Lemasson, J. J. Anderson, and R. A. Goodwin. Collective motion in animal groups from a neurobiologicalperspective: The adaptive benefits of dynamic sensory loads and selective attention.
Journal of TheoreticalBiology , 261(4):501–510, 2009.[72] Bertrand H Lemasson, James J Anderson, and R A Goodwin. Motion-guided attention promotes adap-tive communications during social navigation.
Proceedings. Biological sciences / The Royal Society ,280(1754):20122003, 2013. 14
PREPRINT - S
EPTEMBER
7, 2020Figure S.1: ( A ) Distance regulating force (cid:126)F d ( r ij ) between agents i and j projected on the seperation direction ˆ r ji = (cid:126)r j − (cid:126)r i | (cid:126)r j − (cid:126)r i | . The force equals zero at the preferred distance r d = 1 and is displayed for a distance regulating forcesteepness m d = 2 (used in the simulations) and m d = 4 . ( B ) Folded polar coordinates of an agent i with respect tothe center of mass (cid:126)r com of the school (blue circle) and to the average velocity of the school (cid:126)v com (blue arrow). Theangle α i,com (magenta arc) between the school velocity and the agents i current position (cid:126)r i,com (magenta arrow) andthe distance to the center of mass | (cid:126)r i,com | define the position in this relative coordinate system. A Supporting Informations
A.1 Model-DescriptionA.1.1 Prey-Agents
The prey agents are modeled as active Brownian particles with constant speed v = v and angular noise [67]. Thestochastic equations of motion read: d(cid:126)r i ( t ) dt = (cid:126)v i ( t ) (S.1a) dϕ i ( t ) dt = 1 v (cid:16) F i,ϕ ( t ) + √ Dξ ( t ) (cid:17) , (S.1b)with F i,ϕ ( t ) = (cid:126)F i ( t ) · (cid:126)e ϕ being the force acting on agent i projected on the angular direction (cid:126)e ϕ , D being the angulardiffusion coefficient and ξ ( t ) as Gaussian white noise. For simplicity we omit in the following the explicit timedependence of positions, velocities and forces.Agents react to their environment by (i) coordinating their direction of motion with their neighbors through an align-ment interaction, (ii) by trying to maintain a preferred distance to conspecifics (long range attraction and short-rangerepulsion) and (iii) by a fleeing response (repulsion) from the predator. The alignment force between a focal agent i and all its neighbors j ∈ N i (cid:126)F i,a = 1 | N i | (cid:88) j ∈ N i µ alg · (cid:126)v ji . (S.2)acts towards minimizing the velocity difference (cid:126)v ji = (cid:126)v j − (cid:126)v i with the alignment strength µ alg .Individuals attempt to maintain a preferred distance r d to each other through a distance regulating force (cid:126)F i,d = 1 | N i | (cid:88) j ∈ N i µ d · tanh ( m d ( r ji − r d )) · ˆ r ji (S.3)with ˆ r ji = (cid:126)r j − (cid:126)r i | (cid:126)r j − (cid:126)r i | being the unit vector along the direction from agent i to j , µ d as strength of the force and m d as thesteepness of the change from repulsion (for r ji < r d ) to attraction (for r ji > r d ), as illustrated in Fig. S.1A. Finally ifa predator p is a neighbor of agent i , p ∈ N i , the agent is repelled with (cid:126)F i,f = − µ flee · ˆ r pi (S.4)otherwise (cid:126)F i,f = 0 . The total force governing the movement decision of agent i is defined as (cid:126)F i = (cid:126)F i,d + (cid:126)F i,alg + (cid:126)F i,flee . (S.5)15 PREPRINT - S
EPTEMBER
7, 2020
A.1.2 Predator-Agent
For simplicity the predator obeys a deterministic equation of motion for the heading angle, analogous to Eq. S.1b butwithout the angular noise term: dϕ p dt = 1 v p (cid:126)e ϕ · (cid:126)F p . (S.6)Here, v p is the fixed predator speed and (cid:126)F p is the predator pursuit force. In this study we consider a predator fasterthan the prey v p > v . We assume that the predator can only attack one prey at a time. It considers prey individualswhich are its frontal Voronoi-neighbors N p as targets and selects equally likely among them: p select,i = (cid:40) | N p | if i ∈ N p otherwise . (S.7)The limitation of potential targets to its frontal Voronoi-neighbors N p , is motivated by kinematic and sensory con-straints of the predator. If the predator launches an attack, with an attack rate γ a , which also accounts for potentialhandling time, it’s success probability is linearly dependent on and vanishes distances larger than r catch : p success,i = (cid:40) r catch − r ip r catch if r ip < r catch otherwise . (S.8)In summary, the probability that a predator successfully catches a targeted agent within a small time window [ t, t + δt ] is p catch,i ( t, δt ) = p success,i ( t ) · p select,i ( t ) · γ a δt . (S.9)The predators movement is biased towards the weighted center of mass of the prey school, where each prey positionis weighted by its probability of a successful catch p catch,i ( t, δt ) . Since p catch,i is non-zero only for the predator’sfrontal Voronoi-neighbors, the predator movement are governed by local, visually accessible information. The pursuitforce is thus F p = µ p · (cid:32)(cid:88) i p catch,i (cid:126)r ip (cid:33) . (S.10) A.1.3 Parameter
The default model parameters used are listed in Table S.1. Note that two parameters can be eliminated by renderingthe equations dimensionless. If, for instance, the preferred distance r d and the prey speed v are used to define thecharacteristic length L and time T : L = r d , T = r d v , (S.11)the Eq. S.1 can be reformulated to d(cid:126)r (cid:48) i dt (cid:48) = (cid:126)v (cid:48) i (S.12a) dϕ i dt (cid:48) = r d v (cid:18) F i,ϕ + √ D (cid:114) v r d ξ ( t (cid:48) ) (cid:19) (S.12b) = F (cid:48) i,ϕ + (cid:114) D r r d v ξ ( t (cid:48) ) . (S.12c)Here is D r = Dv the rotational diffusion coefficient (with the unit [ D ] = 1 /t ), the primed variables are the dimension-less counterparts t = r d v t (cid:48) , v i = v v (cid:48) i , r i = r d r (cid:48) i (S.13)and note that the Gaussian stochastic process is transformed according to ξ ( t ) = (cid:114) v r d ξ ( t (cid:48) ) . (S.14)16 PREPRINT - S
EPTEMBER
7, 2020 parameter symbol value p re y angular diffusion D µ alg evolvesdistance strength µ d m d r d v µ flee p re d a t o r speed v p µ p γ a r catch s i m u l . number of agents N dt T eq T simu γ m σ m v and preferred distance r d to 1. All length scales are thus measured in unitsof r d , and all time scales in terms of time needed to move the distance r d . Note that the flee strength µ flee is strictlyspeaking a predator-prey parameter which reduces the prey-only parameters to four.With this choice of characteristic length and time and setting v = 1 and r d = 1 , the dimensionless parameters keeptheir values listed in Table S.1.Since the flee strength µ flee is a predator-prey parameter, the prey system has effectively only four parameters fromwhich the alignment strength µ alg is evolving. The remaining prey-parameters are the angular-diffusion coefficient D which is set to D = 0 . resulting in a persistence time of τ p = v D = 2 , i.e. a solitary agents maintains it currentdirection of motion for approximately the distance of two body length. The distance regulating strength µ d = 2 ischosen to ensures that the prey group stays cohesive. The distance steepness m d = 2 regulates how quick the distanceregulating force saturates to its maximal/minimal values at distances below or above the preferred distance r d (Fig.S.1A).For the predator the speed must be larger than the prey-speed and is set to v p = 2 . Its pursuit strength µ p describestogether with the speed its turning ability and is set to µ p = 2 and therefore equals the preys distance regulating forcestrength. With an capture rate γ c = 1 / and a simulation time of T = 120 around forty prey are caught per roundwhich corresponds to ten percent. The catch radius is set to r catch = 3 and corresponds therefore to three body length.The simulation parameters, and in particular the shoal-size of N = 400 , have been chosen in order to simulatebiologically reasonable behavior, while at the same time limiting the computational costs. For each generation ofthe evolutionary simulations, 76 independent runs are performed, with each equilibrating for T eq = 200 before thepredator appears, and then running for T simu = 120 time units. The time-step is set to dt = 0 . which providessufficient numerical stability and efficient computation (see sectionA.1.5). A.1.4 Evolutionary algorithm and ESS
The evolutionary algorithm is designed to mimic a simplified natural selection at the level of behavioral phenotypes.Among others, the influence of fecundity selection or sexual selection is neglected and the fitness function is onlybased on how likely an individual is captured in a predator attack, which is a biologically reasonable simplificationin the context of predator-prey interactions. The algorithm consists of (i) a fitness estimation step, (ii) a fitness-proportionate-selection step and (iii) a mutation step.(i) The fitness is estimated by running N f = 76 independent attack-simulations on the same phenotype population.For each simulation the γ a · T simu agents with the highest cumulative probability of being caught (Eq. S.9) are declaredas dead. The fitness of agent i is: f i = − N c,i + max ( N c,j , j ) . (S.15)Here N c,i is the number of simulations in which agent i was captured and max ( N c,j , j ) is the largest number of deathsamong all agents. 17 PREPRINT - S
EPTEMBER
7, 2020(ii) The N offspring are generated via the fitness-proportionate-selection. Thereby has one offspring the parametersof the parent i with probability p parent,i = f i (cid:80) j f j . (S.16)(iii) An offspring agent mutates with a probability γ m , the mutation rate, by adding to its alignment strength µ alg aGaussian random variable with zero mean and standard deviation σ m , as the mutation strength.Steps (i) till (iii) are repeated in each generation.Note that instead of step (i) the agents could directly get captured during the simulation and removed from the groupduring the run. This however introduces an additional source of noise in the predation process and the resulting fitnessgradient of the prey would become more noisy. As a consequence the number of generations needed to reach an ESSincreases. Nevertheless, to ensure the robustness of our results we repeated the evolution with captures during theevolution, which did not change the final results (see Sect. A.4). Estimation of the evolutionary stable state (ESS)
In the evolutionary algorithm the finite mutation strength and the stochastic roulette-wheel selection introduce noiseon top of the intrinsic stochasticity of the the predator-prey dynamics (Eq. S.1). This stochasticity is essential for evo-lutionary adaptation and exploration of the phenotype space, but makes it challenging to identifying the evolutionarystable states (ESS) with high precision in evolutionary simulations.To circumvent this uncertainty about the exact optimum we estimate the evolutionary stable state based on the zero-crossing of fitness-gradient estimated from numerical simulations. For a system in generation g with agent parameters (cid:126)µ alg ( g ) ∈ R N + the estimated fitness gradient ∇ f ( g ) is computed by predicting the mean outcome of the fitness-proportionate selection (cid:104) µ alg (cid:105) predict ( g ) = (cid:126)p parent,i · (cid:126)µ alg (S.17a) = 1 (cid:80) Nj f j N (cid:88) i f i µ alg,i (S.17b)and subtracting from it the current mean-value: ∇ f ( g ) = (cid:104) µ alg (cid:105) predict − (cid:104) µ alg (cid:105) . (S.18)Note that, in sake of readability, we omitted for terms on the RHS of Eqs. S.17, S.18 the dependency on the generation g .The average fitness gradient corresponding to an alignment strength is ∇ f ( µ alg , ∆ µ ) = (cid:104)∇ f (cid:105) S µalg, ∆ µ = (cid:80) g ∈ S µalg, ∆ µ ∇ f ( g ) | S µ alg , ∆ µ | (S.19)where S µ alg , ∆ µ is the set of generations which fulfill the condition: µ alg − ∆ µ / ≤ (cid:104) µ alg (cid:105) ( g ) ≤ µ alg − ∆ µ / . (S.20)Therefore Eq. S.19 represents a simple binning of generations with a bin-width of ∆ µ . The maximum of the estimatedfitness landscape, i.e. the evolutionary stable state, is where the estimated fitness gradient is zero and where its slopeis negative. An detail illustration of all components needed to compute the ESS as proposed here is shown in Fig. S.2. A.1.5 Numeric stability
This section deals with the numerical stability of the Euler-Maruyama method used to simulate the stochastic differ-ential equations. The time-step dt should be much smaller than the persistence time τ p = 2 , smaller than the shortestcorrelation time and small enough to fulfill the stability criterion and to avoid oscillating behavior.A general linear stochastic differential equation (here in Langevin form) dxdt = µx + ση ( t ) (S.21)18 PREPRINT - S
EPTEMBER
7, 2020Figure S.2: Details on the estimation of evolutionary stable states of Fig. 3 in the main text. ( A - G ) Fitnessgradient ∇ f for evolution with different flee strength µ flee . Black-dots indicate the estimated fitness gradients foreach generation. Solid lines are averaged fitness gradients. Dashed vertical lines indicate where ∇ f = 0 and thusmark the evolutionary stable states. ( H ) Fitness gradients displayed altogether. Note that the peaks for µ flee = 6 at µ alg ≈ and for µ flee = 8 at µ alg ≈ are due to fluctuations in the standard-deviation of the population. If thestandard-deviation is kept constant those peaks vanish (not shown).which is simulated via the Euler-Mayurama method x n +1 = (1 + dtµ ) x n + σ √ dtη ( t ) (S.22)is algorithmic stable if | dtµ | < . (S.23)The above stability criterion is a Lyapunov-stability for discrete deterministic processes. A more strict criterion is thatthe process approaches its steady state continuously from above or below, i.e. it is not allowed to alternate/oscillatebetween the sides < dtµ < . (S.24)An even stricter criterion is that the time step should be much smaller than the correlation time of a process active inthe system, i.e. restricting the time step to be smaller than a tenth of the correlation time | µ | ≤ dt. (S.25)The above concepts are applicable to the discretization of linear SDEs. To apply the stability criteria to our high-dimensional non-linear stochastic processes (Eq. S.1) we consider the maximal angular change of a focal agent by19 PREPRINT - S
EPTEMBER
7, 2020an isolated social-force. The consideration of isolated forces is reasonable because the distance regulating force ispointing in general in a different direction that the alignment force, e.g. for strongly aligned agents the alignmentforce is parallel to its current velocity but the distance regulating force is pointing away from the closest neighbor.Without loss of generality we rotate the system such that the considered force is always pointing in the x-direction.We also assume that the force is constant.The distance regulating force ( (cid:126)F d ) and flee-force ( (cid:126)F f ) depend on the spatial position of the agent and therefore dependonly implicitly on the heading direction. This means that the dependence on the heading direction is introduced viathe projection of these forces on the direction of angular change ˆ e ϕ,i = [ − sin ϕ i , cos ϕ i ] , which results in dϕ i dt = µ x,i (cid:18) (cid:19) (cid:18) − sin ϕ i cos ϕ i (cid:19) + (cid:112) D r η ( t ) (S.26a) dϕ i dt = µ x,i sin ϕ i + (cid:112) D r η ( t ) . (S.26b)The above equation is non-linear but with respect to stability/convergence we can substitute the sine-function with alinear dependence dϕ i dt = µ x,i ϕ i (S.27)which will result in stronger forces, shorter correlation time and a stricter condition. Note that x ∈ { µ d , µ flee } .The alignment force F alg,i depends explicitly on the heading-direction and since large forces correspond to an orderedstates we assume that (i) the neighbors of the focal agent are strongly aligned with each other and (ii) the deviation inheading direction of the focal agent is small. Thus dϕ i dt = ( µ d (cid:88) j ∈ N i (cid:126)v j − (cid:126)v i ) · ˆ e ϕ,i + (cid:112) D r η ( t ) (S.28a) = µ d ( (cid:104) (cid:126)v (cid:105) N i − (cid:126)v i )) · ˆ e ϕ,i + (cid:112) D r η ( t ) (S.28b) = (cid:18) µ d (cid:18) (cid:19) − (cid:18) cos ϕ i sin ϕ i (cid:19)(cid:19) · (cid:18) − sin ϕ i cos ϕ i (cid:19) + (cid:112) D r η ( t ) (S.28c) = − µ d sin ϕ + (cid:112) D r η ( t ) (S.28d) ≈ − µ d ϕ + (cid:112) D r η ( t ) (S.28e)where we used in the last step the second assumption that deviations from the group-heading direction are small.Thus every force on a focal prey can be approximated by a linear SDE for which the stability consideration eq.S.25needs to hold. A.2 Measures
Here we explain in detail the susceptibility and the relative-positions in the swarm with respect to the front, side anddensity.
A.2.1 Relative positions
In order to define the relative positions with respect to the front and to the side we first represent every agent-positionby its distance to the center of mass of the collective r i,com = | (cid:126)r i,com | = (cid:126)r i − (cid:126)r com with (cid:126)r com = (cid:88) i (cid:126)r i /N (S.29)and the angle between its velocity and the mean velocity of the collective α i,com = ∠ ( (cid:126)r i,com , (cid:126)v com ) with (cid:126)v com = (cid:88) i (cid:126)v i /N . (S.30)We refer to this representation as the folded polar swarm-coordinates , illustrated in Fig. S.1B. Note that the x-axis isparallel to (cid:126)v com , the center of mass is at the origin and the quadrants IV and III are folded onto I and II respectively.The folding is reasonable if a left-right symmetry holds, which we assume. The relative front position is ˜ r i,f = r i,com cos α i,com (S.31)20 PREPRINT - S
EPTEMBER
7, 2020with its normalized version as r i,f = ˜ r i,f − min(˜ r j,f , j )max(˜ r j,f , j ) − min(˜ r j,f , j ) (S.32)which results in front positions in the interval r i,f ∈ [0 , .The relative side-position is ˜ r i,s = r i,com sin α i,com (S.33)with its normalized version as r i,s = ˜ r i,s / max(˜ r j,s , j ) . (S.34)We apply the normalization because we are interested if an individual is at the front and not how far the front is awayfrom the center of mass. In consequence are the normalized measures less noisy if we average over independentinitializations. The average normalized relative-position over S samples is (cid:104) r i,x (cid:105) = (cid:80) Sk =1 r i,x,k S (S.35)with r i,x,k as the normalized relative position of agent i in the k th sample run. Note that the normalized relativeposition is computed after the equilibration time T eq . A.2.2 Local density
The local density of agent i is computed with its distance to the k th nearest neighbour d i,kN to ρ i = k/A ( d i,kN , d i,e ) . (S.36)The term A ( d i,kN , d i,e ) represents the corrected area. If the agents distance to the edge of the collective d i,e is largeras d i,kN , no correction is needed and the area is the area of a circle with radius d i,kN . If the distance to the edge issmaller than d i,kN , the circle-area is corrected by subtracting the area of the circle segment with a sagitta (height) of h = d i,kN − d i,e . Therefore is the area computed as A ( d i,kN , d i,e ) = Φ d i,kN if d i,kN < d i,e Φ d i,kN − d i,kN (cid:18) d i,kN arccos d i,e d i,kN − d i,e (cid:114) − d i,e d i,kN (cid:19) otherwise. (S.37)This correction is good if the edge of the collective has a small local curvature compared to the curvature of the circlewith radius d i,kN . This should be fulfilled because a collective of N = 400 individuals with a preferred distance of r d = 1 and a spherical form has a radius of R ≈ while the distance to the k th nearest neighbor with k = 10 and aVoronoi-interaction network is between 1 and 2. A.2.3 Assortativity
The assortativity r is defined as r = 1 σ q (cid:88) j,k jk ( e j,k − q j q k ) (S.38)with e i,j as the joint probability that a randomly drawn edge connects vertices of type i and j , and q x is the probabilitythat a node of type x is at one end of a randomly drawn edge, i.e. it is the fraction of edges that have a vertex of type x at one end. The assortativity is the Pearson correlation coefficient over the values of the vertices connected by edges. A.2.4 Susceptibility under a homogeneous global field
From equilibrium systems which undergo a second order phase transition, we know that the susceptibility χ , the changeof an extensive macroscopic variable, diverges for an infinite system at the phase transition. If we use the Ising-modelas an example, the susceptibility would describe the change of the magnetization per spin m = M/N = 1 /N (cid:80) Ni =1 s i ,with s i as the spin at side i , under an external field hχ = ∂ (cid:104) m (cid:105) ∂h . (S.39)21 PREPRINT - S
EPTEMBER
7, 2020Interestingly this directly links the reaction of the system to fluctuations in the order parameter. This can be shown ifan energy H ( s i ) = ... − h (cid:80) i s i is assumed and the mean is expressed via the partition function (cid:104) m (cid:105) = 1 /N β ∂ ln Z∂h which leads to χ = 1 β ∂ ln Z∂h = βN [ (cid:10) M (cid:11) − (cid:104) M (cid:105) ] = βN [ (cid:10) m (cid:11) − (cid:104) m (cid:105) ] (S.40)The above fluctuation dissipation theorem connects formally only the response of the system to an infinitesimallysmall change of the external field h . The linear nature of this response to small changes can also be assessed by aTaylor-expansion to linear order of the canonical distribution around h = 0 (see for example Eq. 1.21 in [69]). Theresponse can be rephrased to create a link to the connected spin correlation function or spin pair correlation function χ = N β [ (cid:10) m (cid:11) − (cid:104) m (cid:105) ] = βN [ (cid:42)(cid:88) ij s i s j (cid:43) − (cid:42)(cid:88) i s i (cid:43) · (cid:42)(cid:88) j s j (cid:43) ] (S.41a) = βN (cid:88) ij [ (cid:104) s i s j (cid:105) − (cid:104) s i (cid:105) (cid:104) s j (cid:105) ] . (S.41b)In the following we establish an analog description for our model system with fixed speed. We first compute thesusceptibility for a global homogeneous force, analog to the one used in the Ising-model, and than point out differencesfor the force representing the fleeing behavior of prey, i.e. if it is a local and inhomogeneous force.From the change in heading due to the force F i,s acting on individual idϕdt = F i,s = (cid:126)h ˆ e ϕ v (S.42)part of the Hamiltonian can be computed via integration to H s,i = − (cid:126)h ˆ e r v . (S.43)The total Hamiltonian H = H m ( { v i } , { ϕ i } ) + (cid:80) i H s,i ( v i , ϕ i ,(cid:126)h ) is composed of the sum of isolated components H s,i , which only depend on the state of isolated agents, and of the part of the Hamiltonian which are influenced by theinteractions in between the particles H m . Only H s,i does directly depend on the external field (cid:126)h which means that theensemble average of the order-vector (cid:126) Φ = N (cid:80) Ni ˆ u i can be derived from the partition function as (cid:68) (cid:126) Φ (cid:69) = 1 N β (cid:126) ∇ (cid:126)h ln Z = 1 N β (cid:32) ∂∂h x ∂∂h y (cid:33) ln (cid:88) { r,ϕ } c H e β(cid:126)h · (cid:126)M , (S.44)with c H = e − βH m and (cid:126)M = N (cid:126) Φ . Thus the susceptibility is the gradient with respect to (cid:126)h of this ensemble averagewhich can be written more compact with the (cid:126)h -Laplace operator ∆ (cid:126)h = ∂ h x + ∂ h y to χ = (cid:126) ∇ (cid:126)h (cid:68) (cid:126) Φ (cid:69) = 1 N β ∆ (cid:126)h ln( Z ) (S.45a) = βN [ (cid:10) M x + M y (cid:11) − (cid:104) M x (cid:105) + (cid:104) M y (cid:105) ] (S.45b) = βN [ (cid:68) (cid:126)M · (cid:126)M (cid:69) − (cid:68) (cid:126)M (cid:69) · (cid:68) (cid:126)M (cid:69) ] (S.45c) = βN [ (cid:68) (cid:126) Φ · (cid:126) Φ (cid:69) − (cid:68) (cid:126) Φ (cid:69) · (cid:68) (cid:126) Φ (cid:69) ] = βN [ (cid:10) Φ (cid:11) − (cid:104) Φ (cid:105) ] . (S.45d)22 PREPRINT - S
EPTEMBER
7, 2020Which is analog to Eq. S.40 and the link to the pair-correlation between individual heading directions, analog to Eq.S.41, is also possible: χ = N β [ (cid:68) (cid:126) Φ · (cid:126) Φ (cid:69) − (cid:68) (cid:126) Φ (cid:69) · (cid:68) (cid:126) Φ (cid:69) ] (S.46a) = βN [ (cid:42)(cid:88) i ˆ u i · (cid:88) j ˆ u j (cid:43) − N (cid:68) (cid:126) Φ (cid:69) · (cid:68) (cid:126) Φ (cid:69) ] (S.46b) = βN [ (cid:42)(cid:88) ij ˆ u i · ˆ u j (cid:43) − (cid:88) ij (cid:68) (cid:126) Φ (cid:69) · (cid:68) (cid:126) Φ (cid:69) ] (S.46c) = βN (cid:88) ij [ (cid:104) ˆ u i · ˆ u j (cid:105) − (cid:68) (cid:126) Φ (cid:69) · (cid:68) (cid:126) Φ (cid:69) ] (S.46d) = βN (cid:88) ij (cid:68) (ˆ u i − (cid:68) (cid:126) Φ (cid:69) ) · (ˆ u j − (cid:68) (cid:126) Φ (cid:69) ) (cid:69) . (S.46e)At a continuous phase transition the susceptibility diverges at zero external field. Which means in equilibrium sta-tistical physics that a small change/perturbation in h leads to a new equilibrium state with an order parameter whichdiffers to the unperturbed one by χ . This difference in order parameter is at the order-transition largest. Alternative: maximum entropy method
It is not necessary to rely on results of equilibrium statistical physics if the maximum entropy approach is used. Theonly condition which needs to be fulfilled is the stationarity of the system, it should be in a stable state. The basicidea, as described e.g. by Bialek et. al.[19], is that the observables of the system are described by the most simpleprobability distribution which is quantified by the Shannon entropy of the distribution S ( P ) = − (cid:88) v P ( v ) ln P ( v ) . (S.47)Since the distribution also needs to result in the experimental observables (cid:104) O µ (cid:105) = (cid:88) { v } P ( v ) O µ = (cid:104) O µ (cid:105) exp (S.48)these has to be included in the generalized entropy with Lagrange multipliers h x and h y for the observables Φ x and Φ y respectively. There also exists the obvious constrain that the probability distribution sums up to 1 which is takeninto account by the lagrange multiplier λ resulting in S g ( P, h x , h y ) = S ( P ) − h x (cid:104) (cid:104) Φ x ( v ) (cid:105) P − (cid:104) Φ x ( v ) (cid:105) exp (cid:105) − h y (cid:104) (cid:104) Φ y ( v ) (cid:105) P − (cid:104) Φ y ( v ) (cid:105) exp (cid:105) − λ [ (cid:104) (cid:105) P − (S.49a) = S ( P ) − (cid:126)h · (cid:20)(cid:68) (cid:126) Φ (cid:69) P − (cid:68) (cid:126) Φ (cid:69) exp (cid:21) − λ [ (cid:104) (cid:105) P − . (S.49b)Note that it is possible to take the x and y component of the polarization, despite their obvious correlation, separatelyinto account. The observables are always computed from the detailed state of the system v and therefore always obeythe inequality (cid:113) Φ x + Φ y ≤ . The two constraints are also linear independent because only (cid:104) Φ x (cid:105) exp = 1 totallydefines (cid:104) Φ y (cid:105) exp to 0 and since we are dealing with a stochastic system this does not happen.To fulfill the simplicity criterion the generalized entropy Eq. S.49 needs to be maximized with respect to P ( v ) : ∂S g ∂P ( v ) = 0 (S.50a) = − (cid:88) v (cid:104) ln P ( v ) + 1 + (cid:126)h · (cid:126) Φ + λ (cid:105) . (S.50b)And reformulating with respect to P leads to P ( v ) = exp [ − (1 + λ )] exp (cid:104) − (cid:126)h · (cid:126) Φ (cid:105) . (S.51)23 PREPRINT - S
EPTEMBER
7, 2020The optimization, as done above, with respect to λ normalizes the distribution P ( v ) = exp (cid:104) − (cid:126)h · (cid:126) Φ (cid:105) Z ( { v } ) (S.52)with Z ( { v } ) = (cid:80) v exp (cid:104) − (cid:126)h · (cid:126) Φ (cid:105) .Remind that the stationarity is a necessary condition for applying the maximum entropy model. Therefore describesthe susceptibility as derivative of (cid:126) Φ with respect to (cid:126)h utilizing the stationary probability distribution P ( v ) the changeof the polarization-vector from one stationary state to another one. The transient regime is therefore not covered.longer time in a perturbed state and has therefore a longer collective memory. Theoretical distinction between susceptibility and predator response
For the derivation in Sect. A.2.4 we assumed that (i) the system is in thermodynamic equilibrium (ii) the changes ofthe external field are small and it is (iii) global and (iv) homogeneous. These four conditions might be violated for thereaction of a collective reacting to a predator. • equilibrium state : The system is an active system and therefore per definition a non-equilibrium system.The agents dissipate constantly energy (no conservation of momentum) but, due to an unspecified energysource, keep their preferred speed, i.e. the system is out of thermal equilibrium. However, there exist steadystates at which macroscopic measures which quantify the system do not change apart from inherent fluctua-tions. Therefore different steady states of the system can be compared instead of the equilibrium states (asdemonstrated in the section above). • small changes of an external field: In the context of an predator attack the perturbing force is the flee-forceof the agent. This flee-force is far from small and normally dominates all other forces. Therefore is the linearapproximation used for computing the susceptibility not justified at all. • global field : The global homogeneous field simplified the former analytical derivations of the susceptibility.The flee-force is however neither global nor homogeneous. The flee-force acts only on agents which directlysense the predator. If we assume visual interactions with occlusion by conspecifics, but also with metric-,Voronoi-interaction and other local interaction types, the sensing of the predator is per definition local. • homogeneous field : The flee-force is in the simplest case a repulsion force and therefore inhomogeneous.However, close individuals have similar relative position with respect to the predator and therefore also asimilar flee-force. Thus the homogeneity of the force is on the local scale given but not on the global one.We already demonstrated that the system does not rely on the first assumption because the steady states can be com-pared as done with the maximum entropy approach. However, the stationarity is a necessary condition for applyingthe maximum entropy model. Therefore describes the susceptibility as derivative of (cid:126) Φ with respect to (cid:126)h utilizing thestationary probability distribution P ( v ) the change of the polarization-vector from one stationary state to another one.However, phase transition are up to a certain degree analog to bifurcations, i.e. both mark the emergence of differentsteady states. And, as it is typical for bifurcations, also at phase transitions critical slowing down occurs. This meansthat the dynamic of the system slows down and the relaxation to the steady state takes longer the closer the systemis to the phase transition. The attack of a predator is fast and the predator does not wait for the collective to reach asteady state to continue. Critical slowing down might be beneficial in other context: the collective stays for a longertime in a perturbed state and has therefore a longer collective memory. However, in the predator-avoidance context itis an additional reason, with the other mentioned unmet assumptions, that the susceptibility can be different from anactual predator-response. A.3 Balancing social vs. direct predator information
We identified in the main text a possible explanation for the dependence of the evolutionary stable alignment strengthon the flee strength as observed in the main text Fig. 3B. A prey can benefit from stronger alignment if it has no personalinformation about the predators position. The benefit increases the quicker the alignment and therefore should increasewith alignment strength. But if the prey is fleeing already, i.e. it has personal information on the predator position,than alignment to uninformed neighbors can hinder an escape. Therefore we expect a balance of benefit and deficit.In the following we show its semi-analytical approximation which reveals the observed linear dependence. Note thatwe do not want to claim that it is the only explanation for the linear dependence but a reasonable one.24
PREPRINT - S
EPTEMBER
7, 2020Figure S.3: ( A ) Illustration of angle-vector-relations for variables used in Eq. S.54 and the following. The blackvector (cid:126)v i is the current velocity of agent i . The blue vector (cid:104) (cid:126)v j (cid:105) N i is the mean velocity of the neighbors of agent i .The red vector (cid:126)F flee represents the flee force experienced by agent i if the predator is sensed. The angle α is theangle between the mean velocity of neighbors and the velocity of agent i . The angle θ is the angle between the meanneighbor-velocity and the flee force. ( B ) Numerical-results of the relative direction to neighbors α using Eq. S.54.The initial conditions is α = 0 , i.e. the focal agent is perfectly aligned with its neighbors. The angle between meanneighbor velocity and flee force is θ = π/ . The different colors indicate that the effective flee direction, which is thecompromise between the mean neighbor velocity and the flee-force, is faster the stronger the flee strength µ flee . Weassumed, as discussed in Sect. A.3, that the directional compromise represents the balance between the benefit anddeficit of social information and is maintained, i.e. the prey evolve their alignment strength µ alg to keep the effectiveflee direction constant.The deficit to align with uninformed prey if the predator position is known can be viewed as a deviation from the fleedirection, i.e. the prey relaxes to an effective flee direction which is the compromise between the mean direction of itsneighbors and the flee direction Fig. S.3.We will use the following assumptions: • i) highly ordered: all neighbors are perfectly aligned with each other. • ii) strong forces: the acting forces are strong such that the agents equilibrate quickly in the direction of theforce. • iii) constant forces: the flee-angle and the heading of the neighbors are not changing. • iv) no noise: this will enable us to solve the problem analytically.Consequently the change of the direction-angle of Eq. S.1b can be reformulated to dϕ i dt = 1 v (cid:16) F i,ϕ + √ Dξ (cid:17) (S.53a) ≈ v ( F i,ϕ ) (S.53b) ≈ v (cid:16) µ flee ˆ f flee + µ alg [ (cid:104) (cid:126)v (cid:105) N i − ˆ e r,i ] (cid:17) · ˆ e ϕ,i . (S.53c)With (cid:104) (cid:126)v (cid:105) N i being the mean velocity of all neighbors of agent i and ˆ e r,i and ˆ e ϕ,i are its heading and angular direction,respectively. Without loss of generality we can permanently rotate the system such that ϕ = 0 , ∀ t which simplifies thevector products since ˆ e r,i = [1 ,
0] = ˆ e x and ˆ e ϕ,i = [0 ,
1] = ˆ e y . The angle α between (cid:126)v i and (cid:104) (cid:126)v (cid:105) N i behaves exactly25 PREPRINT - S
EPTEMBER
7, 2020opposite as ϕ (see Fig. S.3A) and we describe its dynamics instead: dαdt = − dϕdt (S.54a) ≈ − v (cid:16) µ flee ˆ f flee + µ alg [ (cid:104) (cid:126)v (cid:105) N i − ˆ e x ] (cid:17) · ˆ e y (S.54b) ≈ − v ( µ flee f flee,y + µ alg (cid:104) (cid:126)v (cid:105) N i ,y ) . (S.54c)With f flee,y = sin( θ − α ) and by assuming perfect order and unit speed the mean velocity of neighbors is (cid:104) (cid:126)v (cid:105) N i =1 (cid:18) cos( α )sin( α ) (cid:19) . Therefore the change of α simplifies to: dαdt ≈ − v ( µ flee sin( α − θ ) + µ alg sin α ) (S.55a) ≈ µ flee sin( θ − α ) − µ alg sin α. (S.55b)The fixed points are, as a sanity check, computed for the extreme cases µ alg (cid:29) µ flee and µ flee (cid:29) µ alg which are α (cid:63) = 0 and α (cid:63) = θ , respectively. There exist in general four fixed points from which only one fulfills the criteria α (cid:63) /θ ∈ [0 , ∀ ( µ flee > , µ alg > , < θ < π/ which is: α (cid:63) ( θ s , µ alg , µ flee ) = arccos µ alg + µ flee cos θ (cid:113) µ alg + µ flee + 2 µ alg µ flee cos θ . (S.56)Thus α (cid:63) is the effective flee angle with respect to the mean direction of the neighbors. The closer it is to the flee angle θ the smaller the deficit of being aligned given the knowledge of the predators position.Now we assume that individuals evolve such that they maintain α (cid:63) ( θ s ) with respect to a specific θ s . Thus, if we knowthe equilibration point µ (cid:63)alg,evo ( µ flee,evo ) for a specific flee strength which was used during the evolution µ flee,evo ,we can compute the effective flee angle α (cid:63) ( θ s , µ (cid:63)alg,evo , µ flee,evo ) = α (cid:63) ( θ s ) . If we assume that agents evolve suchthat the balance between alignment benefit and deficit, manifested in the effective flee angle, is kept constant, than wecan predict the evolutionary stable state µ (cid:63)alg for a given flee strength by reformulating Eq. S.56 to µ (cid:63)alg = sin( θ s − α (cid:63) ( θ s ))sin α (cid:63) ( θ s ) µ flee . (S.57)From this equation it becomes clear that the exact choice of θ s is irrelevant since sin( θ − α (cid:63) )sin α (cid:63) must be constant becauseit corresponds to the slope of the line through the origin and the one evolutionary stable state ( µ (cid:63)alg,evo , µ flee,evo ) usedto compute α (cid:63) ( θ s ) as shown by the blue line in Fig. 3B.Note that the equilibration alignment strength µ (cid:63)alg above but close to the order transition is systematically lower thanits predicted value, as seen for µ flee ∈ { , , } in Fig. 3B. This can be explained by a small signal due to the lowflee strength, beecause the system relaxes faster the greater the flee strength µ flee (see Fig. S.3B). An alternativeexplanation is that the spatial selection due to strong self-sorting dominates at the transition. This explanation is alsoin agreement with the ESS for low flee strength ( µ flee = 0 . ) being identical to the one with no flee strengt at all( µ flee = 0 ). A.4 Robustness against modifications in prey, predator and selection mechanism
To ensure that our results are robust we repeat the evolution (Fig. S.4) with (i) modified prey properties, i.e. changingthe angular diffusion coefficient and introducing variable speed, (ii) a changed predator behavior, i.e. its agility, and(iii) changes in the evolutionary selection mechanism, e.g. by an additional high-frontal-risk selection mechanism orby a prey capture during the simulation. Note that especially the additional high-frontal risk selection is of importance,because it introduces a heterogeneous environment which is assumed to be a general important condition for theevolution to criticality[11].
A.4.1 Prey modifications
The change in angular diffusion from D = 0 . to D = 1 shifts the order-transition to a larger mean alignment strengthof µ alg,c ≈ . and therefore also increases the lower bound for the ESS which is visible in larger ESS for small26 PREPRINT - S
EPTEMBER
7, 2020Figure S.4:
Robustness of evolution results:
Evolutionary stable states of the alignment strength are estimatedfrom the fitness gradient for different flee strength under slight variations of simulations parameters or predator attackimplementation. The standard scenario of the main text (blue line) is compared to ( A :) a prey population with varyingspeed which can avoid the predator additionally by acceleraton (black dotted line), a prey population with a angulardiffusion coeffient which is doubled compared to the standard case (red dashed line), ( B :)a stiff predator which turnsless quick (black dotted line) and an agile predator which turns quicker (red dashed line) than the predator in thestandard case. ( C :) a non-binarized fitness estimate (red dashed line) in which the preys fitness is not defined bycaptures but by the accumulated probability to get caught, a fitness estimate based on captures during the simulation(black dotted line),flee strength (compare dotted red with blue line in Fig. S.4A). For larger flee strength the results are nearly identicalsuggesting that the mechanism defining the ESS remains unchanged with respect to the standard scenario of the maintext.If the speed of the prey is not constant but can change according to social forces, the equations of motion (Eq. S.1)change to d(cid:126)r i dt = (cid:126)v i with (cid:126)v i = v i [cos ϕ i , sin ϕ i ] (S.58) dv i dt = β ( v − v i ) + F i,v ( t ) (S.59) dϕ i ( t ) dt = 1 v (cid:16) F i,ϕ ( t ) + √ Dξ ( t ) (cid:17) (S.60)with F i,v ( t ) = (cid:126)F i · ˆ e h,i as the projection of the social force of prey i on its heading direction ˆ e h,i and β as therelaxation coefficient which is set in the following to β = 4 . A value of β = 4 prevents the school to relax into a non-moving phase which exists for lower values of β [70]. In this non-moving state the speed of the prey fluctuates aroundzero. Additionally we set an upper bound for the preys speed corresponding to eighty percent of the predators speed v max = 0 . v p . Non-fleeing prey ( µ flee = 0 ) evolve to significant larger values compared to the standard scenariofrom the main text (compare dashed black with blue line in Fig. S.4A). The ESS for non-fleeing prey ( µ flee = 0 )coincides with the zero-crossing of the front-sorting (Fig. S.5). Not only is the ESS of the non-fleeing prey at largervalues due to a different self-sorting but also is the ESS much more sensitive to changes in the flee strength (compareslope of dashed black with blue line Fig. S.4A). This steeper increase is explainable with an additional social cue,the increased speed of fleeing-prey, which was not present in the standard-scenario and goes in hand with findings byLemmasson et. al.[71, 72]. A.4.2 Predator modifications
We repeated the simulations with (i) a stiffer predator which turns slower and (ii) a more agile predator which turnsfaster compared to the predator in the standard scenario of the main text. The different turning ability was implementedby modifying the pursuit strength µ p to µ p = 1 for the stiff and to µ p = 3 for the agile predator.The effect of using the stiff predator is negligible for low flee-strength, probably because the order-disorder transitionacts as lower bound for the ESS due to the explained maximum in assortative mixing and resulting subpopulationselection. However, for larger flee strength, e.g. µ flee ∈ { , } in Fig. S.4B, the ESSs are lowered compared to thestandard scenario in the main text. This can be explained by the missing feedback between the reaction of the prey andthe trajectory of the predator: in the standard scenario the predator heads for the closest prey, thus if certain prey aregood at avoiding the predator, they have an additional fitness benefit because the predator pursues less well-avoidingprey. 27 PREPRINT - S
EPTEMBER
7, 2020Figure S.5:
Self-sorting with and without fixed speed:
Self-sorting quantified via the Pearson correlation betweenthe individual alignment parameter µ alg and the average relative position of the individuals (relative front-, side- ordensity-location as described in Sect. A.2.1). ( A ) If prey agents respond only by changing their direction but not theirspeed (fixed speed), self-sorting persists also in highly ordered regions. ( B ) If prey agents can change their speed(variable speed), self-sorting vanishes for µ flee ≤ .Figure S.6: Evolution in heterogeneous environments:
Fitness gradients for different relative strength of the frontal-risk selection with respect to the simultaneously active predator-selection. In the frontal-risk selection the most frontalindividuals are decleared as dead. The relative strength of the frontal-risk selection is defined by the ratio betweenagents killed at the front and by the predator, i.e. ( Front Kills ) / ( Pred. Kills ) ∈ [0 , . , . , . . The evolutionarystable state (ESS) is defined by the zero-crossing of the fitness gradient with negative slope marked by a vertical dashedline. However, the lower bound is an additional ESS if the fitness gradient stays negative close to it which is markedby shaded points in the inset. Parameters are identical to the former simulations apart from the angular diffusioncoefficient which is increased to D = 1 increasing the order-transition to µ alg,c ≈ . marked by vertical dash-dottedmagenta line. The flee strength is µ flee = 4 .Consequently should the trajectory of the more agile predator have a stronger feedback with the reaction and lead tolarger evolutionary stable alignment strength. This is in fact observed (compare dotted black with blue line in Fig.S.4B).Despite the explainable variations due to the predator modifications the general finding in the main text, i.e. that theESSs are in the ordered phase and increase with increasing flee-strength, is robust. A.4.3 Selection modification: Evolution in a heterogeneous environment
In the simulations prey are not captured but a fixed fraction of them with the largest accumulated probabililty toget caught is declared as captured after the simulation. This means that no prey is removed during the simulationwhich reduces stochasticity of the fitness estimate but can be considered as unrealistic. If prey are removed during thesimulation based on their current probability to get caught and the predators attack rate, the evolution results remain28
PREPRINT - S
EPTEMBER
7, 2020unchanged (compare dotted black with blue line in Fig. S.4C). Hereby is the attack rate γ a adjusted at each generation g such that the mean capture rate (cid:104) γ c (cid:105) matches the initially set attack rate γ a ( g = 0) : γ a ( g + 1) = γ a ( g ) ∗ γ a (0) (cid:104) γ c ( g ) (cid:105) . (S.61)This ensures a constant evolutionary pressure.The attack rate parameter can be abandoned if the fitness is not estimated by the captures but by the negative accu-mulated probability to get caught. This modification does not alter the ESS identified in the main text at all (comparedashed red with blue line in Fig. S.4C).The chosen predator-prey interaction is set as general as possible, nevertheless reasonable alternatives exists and otherenvironmental interactions, e.g. exploration and exploitation of food-sources, might simultaneously impact the fitness.We introduce an additional selection mechanisms which favors a disordered phase and creates thus a heterogeneousenvironment. The self-sorting for this model predicts that a high mortality of front individuals leads to a disorderedstate which we implement by declaring the most frontal prey as dead. This extra selection is equivalent with theobserved high risk of being in the front in the presence of sit-and-wait predators[38]. Since the current transition isclose to the lower boundary of the alignment parameter ( min( µ alg ) = 0 ), we set the transition at larger values, i.e.at µ alg,c ≈ . , by increasing the angular diffusion to D = 1= 1