Collective predator evasion: Putting the criticality hypothesis to the test

According to the criticality hypothesis, collective biological systems should operate in a special parameter region, close to so-called critical points, where the collective behavior undergoes a qualitative change between different dynamical regimes. Critical systems exhibit unique properties, which may benefit collective information processing such as maximal responsiveness to external stimuli. Besides neuronal and gene-regulatory networks, recent empirical data suggests that also animal collectives may be examples of self-organized critical systems. However, open questions about self-organization mechanisms in animal groups remain: Evolutionary adaptation towards a group-level optimum (group-level selection), implicitly assumed in the “criticality hypothesis”, appears in general not reasonable for fission-fusion groups composed of non-related individuals. Furthermore, previous theoretical work relies on non-spatial models, which ignore potentially important self-organization and spatial sorting effects. Using a generic, spatially-explicit model of schooling prey being attacked by a predator, we show first that schools operating at criticality perform best. However, this is not due to optimal response of the prey to the predator, as suggested by the “criticality hypothesis”, but rather due to the spatial structure of the prey school at criticality. Secondly, by investigating individual-level evolution, we show that strong spatial self-sorting effects at the critical point lead to strong selection gradients, and make it an evolutionary unstable state. Our results demonstrate the decisive role of spatio-temporal phenomena in collective behavior, and that individual-level selection is in general not a viable mechanism for self-tuning of unrelated animal groups towards criticality.

Introduction Distributed processing of information is at the core for the function of many complex systems in biology, such as neuronal networks [1], genetic regulatory networks [2] or animal collectives [3,4]. Based on ideas initially developed in statistical physics and theoretical modeling it has been conjectured that such living systems operate in a special parameter region, in the vicinity of so-called critical points (phase transitions), where the system's macroscopic dynamics undergo 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 neuronal dynamics [10,12,13], gene regulatory networks [14,15], and collective behaviors of animals [16][17][18][19][20][21]. This evidence is often based on observation of characteristic features of critical behavior, such as power-law distribution or diverging correlation lengths in spatial systems. However these observations could in principle have different origins [12,[22][23][24]. Therefore, more convincing support for the "criticality hypothesis" can be obtained through additional identification of proximate mechanisms enabling biological systems to self-organize towards criticality. In neuronal systems, synaptic plasticity has been shown to provide such a mechanism [25][26][27]. For genetic regulatory networks, similar mechanisms based on network rewiring have been proposed [28,29]. Using an information-theoretic framework Hidalgo et al. [11] have shown that (coupled) binary networks evolve towards the critical state in heterogeneous environments. However, in their model already a single unit (individual) can exhibit a phase-transition and thus tunes itself individually to criticality. In addition, they assumed idealized random interaction networks between the agents. Thus, open questions remain whether evolutionary, individual-level adaptation is a possible self-tuning mechanism for (i) biological collectives, where phase transitions are purely macroscopic phenomena, and (ii) animal groups characterized by spatial, dynamic interaction networks. In general, if collective computation becomes optimal at a phase transition, a purely macroscopic phenomenon defined only at the group-level, then adaptation based on global fitness should be able to tune the system towards criticality. Therefore, at first glance Darwinian evolution appears a viable mechanism for emergence of self-organized criticality only for complex systems within a single individual, e.g. in the context of neuronal or genetic networks, or in collectives of closely related individuals such as eusocial insects [16]. In multi-agent systems group-level and individual-level evolutionary optima are often different [30,31], leading to socalled social dilemmas emerging in a broad range of multi-agent evolutionary game theoretic problems. In the context of animal groups consisting of non-related individuals, this questions individual-level adaptation as a proximate mechanism for self-tuning of collective systems to criticality as a potential group level optimum. Here multi-level selection has been proposed to address some related fundamental problems in the evolution of collective behavior [32,33]. However, it has been recently shown that even under strong group-level selection, as long as individual-level selection plays a non-negligible role, multi-level selection will also result in evolution of sub-optimal collective behaviors [34,35].
Whereas few empirical studies report signatures of criticality in collective animal behavior [18][19][20], most support for the criticality hypothesis in this context comes from mathematical models. For example, in agent-based simulation of fish schools it has been shown that at a critical point the collective state is influenced strongest by single or few individuals [36], or that collective response to external time-varying signals becomes maximal in idealized lattice models of flocks [37]. However, dynamical animal groups differ from lattice models [37,38] due to their dynamical neighborhood which may induce self-sorting of individuals according to their individual behavioral parameters [39][40][41]. This, in turn, has likely direct evolutionary consequences as for example predators may attack certain swarm regions more frequently [42][43][44].
Throughout this work, criticality or critical point will refer to the directional order-disorder transition, a prominent phase transition in statistical physics and collective behavior [45]: An initially disordered swarm, where the social coordination is weak compared to noise, shows spontaneous onset of orientational order, if the directional alignment (coupling strength) is increased beyond a critical parameter (critical point): The group starts to move collectively along a common "consensus" direction. A further increase of alignment results in highly ordered (polarized) schools [46]. This transition is characterized by a so-called spontaneous symmetry breaking: In disordered swarms there is no distinguished direction in space. In the ordered state, this symmetry is broken through the emergence of an average 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 coordinated collective behavior of the prey is believed to entail evolutionary benefits to individuals within the group [47]. In particular, we use an agent-based model of grouping and coordinating prey [39,40,[48][49][50][51], and analyze the role of the spatial structure of the group, its dynamical response and the individual-level selection by applying an evolutionary algorithm [52][53][54][55][56][57][58].
We show that the group-level behavior becomes optimal at criticality with respect to two measures: We observe i) 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, becomes minimal only due to the dynamical structure of the collective at criticality, independent of the direct response of individuals to the predator, and thus independent of information propagation within the school. Furthermore, through evolutionary simulations with individual-level selection, we show that the critical point is an evolutionary highly unstable state. This evolutionary instability can be linked to strong selection due to phenotypic-sorting with respect to the broken symmetry of the collective state. Finally, the observed evolutionary stable strategies (ESS) result from individual prey agents balancing the influence of social and private information on their movement response.

Agent based model of predator-prey interactions
We consider a simple, yet generic agent-based model of schooling prey attacked by a predator. For simplicity we assume initially that the prey agents move with fixed speed v 0 and change their direction according to social forces ( Fig 1A): they tend to keep a preferred distance to, and align (alignment strength μ alg ) their velocities with the first shell of nearest neighbors, defined by a Voronoi tessellation ( Fig 1B) [4,59]. A distance regulating social force represents a continuous version of a two zone model, i.e. repulsion at short distances and an attraction zone at large distances with a "preferred" (equilibrium) inter-individual distance r d (Fig 1A). Randomness in the movement of individuals due to unresolved internal decisions or environmental noise is modeled as fluctuations in the heading of the agents (angular noise with intensity D). The prey responds with a flee-force with strength μ flee to a predator within its Voronoi neighbors (Fig 1C).
The predator moves with a fixed speed v p which is larger than the preys (here v p = 2v 0 ) and its direction changes towards the weighted mean direction of its frontal nearest prey, which represent possible targets (Fig 1C). The weight corresponds to the catch-probability of each target, which decreases linearly with distance until it equals zero at a distance larger than the catch-radius. If the predator launches an attack, with attack rate γ a , it selects equally likely among the possible targets and captures it according to the targets catch-probability. The predator is initiated outside the prey collective with a distance slightly above the capture-radius and a velocity vector oriented towards the center of mass of the prey school.
In evolutionary simulations for each generation we perform N r independent runs with different initial conditions for N agents, each with its behavioral phenotype defined by the evolvable social force parameter (alignment strength μ alg ). Fitness of a prey agent is defined through the negative number of deaths of this agent aggregated over the N r independent runs. The behavioral phenotypes, i.e. social force parameters, of the next generation are selected via fitness-proportionate selection (roulette-wheel-algorithm) [56,60,61] with mutations implemented through addition of Gaussian-distributed noise on the selected behavioral parameter. See Methods for model details. The social forceF soc acts on the focal agent (black triangle) and is a combination of alignmentF alg and distance regulationF d to its interaction partners. The alignment is proportional to the sum of the velocity differences F i;alg / P jṽ ji withṽ ji ¼ṽ j Àṽ i and thus not parallel to the neighbors mean velocity but tends to minimize the velocity difference. The distance regulating force F d is a continuous version of a two zone model, i.e. the focal agent is repelled from neighbors that are closer (red triangle=line) than the preferred distance r d (gray dashed circle=line) and attracted to those farther away (blue triangle=line). (B) A focal prey agent (yellow triangle) interact with it's nearest Voronoi neighbors (black triangles in yellow cells). (C) The predator (red point) pursues the weighted mean direction of the targets (small red triangles), which are the frontal Voronoi neighbors. Their weight is proportional to their probability of capture, which decreases linear with distance and is zero for r � r catch (magenta semicircle). All Voronoi neighbors of the predator flee with a repulsive forceF flee (red arrows). https://doi.org/10.1371/journal.pcbi.1008832.g001

Collective information transfer and responsiveness
We first investigate whether operating at the order-disorder transition leads to optimal response of the prey school to the predator. Here, polarization F, i.e. the normalized average velocity of the group, is the relevant order parameter quantifying the amount of orientational order in the system: For large, disordered systems F is close to zero, while in completely ordered systems with all agents moving in the same direction it approaches 1 (see Methods). It increases with the strength of alignment μ alg and decreases with the intensity of angular noise D (S1 Movie) in a non-linear fashion: It remains small (F � 0) throughout most of the disordered regime, before showing the steepest increase in orientational order in the vicinity of the critical point, and finally asymptotically approaching F = 1. Both behavioral parameters, μ alg and D can be used as control parameters for crossing of the critical line (diagonal magenta line Fig 2A) between the disordered 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 average pair-wise correlation of velocity fluctuations C ij ¼ Cðdṽ i ; dṽ j Þ between interacting agents (see Methods). Here, dṽ i ¼ṽ i À hṽi is the deviation of the velocity of agent i from the average school velocity hṽi. C ij can be interpreted as a simple measure of directional information transfer between neighboring agents i and j: If agent i deviates from the average 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 quantifies the degree of responsiveness of the system to perturbations, and may become maximal at criticality. It can be defined analogous to magnetic susceptibility in physics [36,62] (see Methods).
Both measures, C ij and χ, show a peak at the transition between order and disorder (see Fig  2A and 2B) in line with predictions of the "criticality hypothesis" [13]. In terms of directional information transfer, i.e. the directional responsiveness to perturbation, it appears to be optimal for the collective to operate at criticality.

Fitness relevant performance measure
The validity of the above variables from a statistical physics point of view relies on the assumptions of homogeneity and temporal stationarity of the external field, which is not fulfilled in our predator-prey scenario: predator perturbation represents a strongly local, nonlinear perturbation. As a biologically relevant measure, independent of these assumptions, we use directly the predator capture rate γ c , computed as number of prey captured per time unit. In agreement with the previous response measures, we find that the capture rate also exhibits a distinct minimum at the critical point ( Fig 2C).
However, varying the behavioral parameters of the prey (alignment strength or noise) not only changes the polarization of the school and the information transfer capability but it also affects the spatial structure of the school (S1 and S2 Movies), e.g. the average inter-individual distance (IID) or the shape of the school. Our results show that structural properties of the prey school correlate strongly with the capture-rate, e.g. the inter-individual distance (inset Fig  2C) with C(γ c , IID) � −0.69. Thus, the reduced capture rate may be potentially related to changes in the structure of the school at criticality. To distinguish whether structure or information transfer is responsible for the optimal performance of the group at the critical point, we simulated for each predator attack a non-fleeing prey school (flee strength μ flee = 0) as a control. This non-responsive control school is identical to the responsive school in all the remaining parameters and in its positions and velocities at the time of predator appearance (see S3 Movie). The capture-rate of the non-fleeing prey γ c,NF depends only on the self-organized structure of the school. We compare the responsive and control school via two measures: (i) the simple difference between 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 2D) suggesting that the predator-response improves towards the ordered phase if we control for the differences in the self-organized spatial structure (compare column μ alg = 1 with μ alg = 2 in S2 Movie).
These results demonstrate that the direct cause of the optimal collective performance (minimal capture rate) is the dynamical structure, as a "passive" component, and surprisingly not the maximal responsiveness at criticality (see S1 Text, V.2 for theoretical reasoning on differences between susceptibility and predator response).

Evolution of coordinated escape
The group-optimum at criticality with respect to prey-survival, does not need to coincide with the evolutionary stable state (ESS) with respect to evolutionary adaptations at the individual level. To explore whether the transition region is favored by individual-level adaptation, we let the individual alignment strength μ alg evolve over 500 generations, while keeping the angular noise constant (D = 0.5: vertical line Fig 2A). We repeat the evolutionary simulations from different initial conditions: below (hμ alg i = 0), above (hμ alg i = 5) and far above (hμ alg i = 10) the transition (μ c,alg � 0.9). To ensure that the evolution ends at the ESS we compute the fitness gradient which represent the strength of the selection pressure at a specific mean alignment strength (see Methods). Assuming a monomodal phenotype distribution, as observed in our evolutionary runs, a change in sign of the fitness gradient marks the location of the ESS. All three initiations end in the ordered region far above the critical point ( Fig 3A) and fluctuate around ESS(μ alg ) � 4.4 (vertical dashed line Fig 3B). Thus, the transition region is not an attractor of the evolutionary dynamics. On the contrary, it is a highly unstable point with fast evolutionary dynamics due to particularly strong selection pressure at criticality. The fitness gradient peaks shortly above the transition in the ordered phase (Fig 3B), with evolutionary dynamics pushing the system 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 to their behavioral parameters along specific spatial dimensions of the school, e.g. front-back or side-center, or in regions of higher or lower density (Fig 3C) [39]. We can quantify self-sorting through the Pearson correlation coefficient between the alignment strength (social phenotype) of an agent and variables quantifying its location within the school (see Methods). Another measure of self-sorting is the amount of assortative mixing in the school as quantified by the assortativity coefficient (see Methods). Assortativity (Fig 3B) as well as other self-sorting measures (Fig 3C) exhibit extrema which coincide with the fitness gradient peak. Note that a strong assortative mixing is equivalent to the formation of spatially coherent sub-groups within the school with similar behavioral parameter. In this context a peak in fitness gradient close to transition suggests that sub-groups with stronger alignment, thus better directional coordination, actively or passively perform better at avoiding capture. An increase in the escape ratio R esc with increasing alignment close to criticality (see Fig 2D) suggest an enhanced active avoidance. However, also passive effects appear to play an important role since the correlation between the fitness of a prey and its relative position becomes maximal in the same parameter region (Fig 3D). One specific mechanism of passive avoidance is the dilution effect [47] caused by local density differences correlating with behavioral phenotypes. Stronger aligning individuals form denser regions within the prey school (density-sorting Fig 3B). As a consequence they have a systematically smaller domain of danger [63] and are thus less frequently attacked by the predator.
It is possible to disentangle passive, structural effects from an active response, by setting the flee-strength to zero. This results in a significantly smaller, yet finite, fitness-gradient-peak at the transition (S2 Fig, panel H). This suggests that both, the structural, passive selection and the different active avoidance behavior of different phenotypes contribute to the strong selection pressure at criticality.
We note that the sudden increase in 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 phase every edge position is equivalent, but with the emergence of the common movement direction the sides and rear of the school become structurally different from the front. This can be clearly seen in the comparison of the correlations of individual alignment strength and 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 to deviate and show different behavior with increasing alignment strength (Fig 3C).

ESS: Balancing benefits and costs of social information
Despite the importance of self-sorting for the maximal selection pressure at the transition, it does not provide an explanation for the observed location of the ESS. More specifically, it can not explain the negative fitness gradient for strong alignment μ alg > ESS(μ alg ) � 4.4. In this regime either the self-sorting is negligible, as for side-and density-sorting (Fig 3C), or the relative location has no effect on the individual fitness, as observed along the front-back dimension (Fig 3D). If the ESS is not determined by the structural self-organization of the school, it has to originate from individuals avoiding the predator better than others. Please note that avoidance does not only mean to escape if targeted but also to avoid becoming a target. In this case the ESS has to depend on the flee-strength μ flee as the main parameter tuning the strength of individual predator response.
We do find a clear dependence of the ESS on the flee-strength (Fig 4A). More specifically, the ESS exhibits a linear dependence on the flee-strength for μ flee � 2 (diagonal line in Fig 4B). The order transition acts as a lower bound since the non-fleeing agents (μ flee = 0) equilibrate closely above it. Thus, the ESS for non-responding agents matches the 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 information. Social information about the predator is beneficial if the prey is in the second neighbor shell of the predator, i.e. where its neighbors but not itself responses directly to the predator. Thus, by coordinating with its informed neighbors it gains distance to the predator. However, if a prey directly senses the predator, social information of uninformed neighbors conflicts with its private information and therefore may hinder evasion. Therefore, individual prey agents should continue to evolve towards stronger alignment strength until costs of the social inhibition of evasion counterbalance the benefits of social information. We find support for this conjecture by reproducing the observed linear dependence through a local mean-field approximation (see S1 Text, Sec. VI, S3 Fig) assuming the above balancing mechanism (Fig 4B). Interestingly, also the escape ratio, as a measure of group response while controlling against spatial effects, exhibits a maximum in the strongly ordered region away from criticality (Fig 2D).
This leads to the question whether the ESS coincides with the largest escape ratio. Indeed, the maximum of escape ratio shows the same trend as the ESS of moving towards higher alignment strengths with increasing flee strength (Fig 4C), but these maxima stay clearly below the corresponding ESSs (circles in Fig 4C). This suggests that the system does evolve towards unresponsiveness [30] by increasing the social responsiveness above the optimum (compare column μ alg = 2 with μ alg = 4 in S4 Movie). We propose that the evolution to unresponsiveness is due to only the targeted prey having a probability of being captured. It appears to be more beneficial for individuals to avoid becoming a target in the first place via a strong social response to fleeing neighbors, rather than being better at escaping once they end up as direct predator targets. Please note, if prey would ignore others during their escape, there would be no trade-off between social and private information about the predator and agents would remain responsive to the predator at the ESS.

Robustness analysis
The qualitative results are independent of model implementation details. We checked for robustness against the predator attack scheme (more and less agile predator), prey-modification (variable speed, persistence length, anisotropy of social interactions / blind angle), modifications in evolutionary algorithm (attack-rate, fitness-estimation) and importantly in a heterogeneous environment (see S1 Text, Sec. VII and S4 and S6 Figs). Note that we explicitly confirmed that considering prey with variable speeds, which enables them to accelerate away from the predator, does not change the qualitative results (S1 Text, Sec. VII.1 and S5 Fig). For strong flee forces corresponding accelerations resemble a typical startle response in fish (S5 Movie).
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 (S6 Fig). Note that our findings are expected to be robust because they are based on generic, modelindependent mechanisms: (i) the maximal self-sorting at the transition combined with the spatial explicit implementation of the predator avoidance (causing the transition to be evolutionary unstable) and (ii) the trade-off between social and personal information (causing the ESS to shift to larger social attention with increasing flee strength). It may be argued that the latter mechanism is biologically not plausible, because prey agents that detect the predator should just flee and ignore their conspecifics. However, this would correspond to a limiting case of a dominating flee-strength and would result in an ESS even further away from the critical point in the highly ordered state (Fig 4B).

Discussion
We have shown, using a spatially-explicit agent-based model of predator-prey dynamics, that the group optimum with respect to predation avoidance is located in the vicinity of the critical point between disordered swarming and ordered schooling, in line with the so-called "criticality hypothesis". However, this optimality is not due to optimal transfer of social information but rather due to the highly dynamical structure of the group at the transition. Yet, this group optimum 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 animal groups when evaluating potential evolutionary benefits of grouping. It turns out that the mechanism responsible for the optimal collective performance (minimal capture rate) at the critical point, the highly dynamic and flexible structure of the collective, leads also to the steepest selection gradients in evolutionary dynamics, making the critical point evolutionary unstable. Evolution with random mutations enforces heterogeneity which in combination with the spatial symmetry breaking at the transition, results in maximal assortative mixing and self-sorting close to the transition. These effects of self-organized collective behavior play a decisive role for the evolutionary dynamics close to criticality and "drive" the ESS out of the transition region towards the aligned state. In our system the ESS is in the strongly ordered phase, which suggests the evolution towards external unresponsiveness by overestimating social information. Finally, we show that the ESS depends linearly on the flee strength, i.e. local perturbation strength, which can be explained by individual balancing of benefits of social information about the predators approach with the costs of 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 already evolve to its "individual" transition by tuning its own gene regulatory network. This appears to be essential for a critical point corresponding also to the evolutionary stable state in their information-based fitness framework. In our model, the disorder-order transition is a pure collective effect, i.e. individual agents cannot exhibit any transition behavior by themselves. Furthermore, at the disorder-order transition, small differences in behavioral parameters translate into systematic differences in the self-organized spatial positioning within the group, which in turn directly impacts the predation threat. This self-sorting [39][40][41] is maximal just above the transition and includes assortative mixing due to emergence of spatial "subgroups" with strong correlations between behavioral phenotype, spatial location and local school structure, which is potentially of interest in the broader context of collective task distribution and computation in spatially-explicit animal groups.
There is another consequence of the tight coupling between local school structure and individual dynamics: The extent of the collective is largest at the transition because the responsiveness to directional fluctuations is maximal, i.e. local fluctuations induce deviations in the movement of different parts of the school causing the school effectively to expand. In systems with a one-way influence from structure to dynamics (fixed networks) it is known that at the order-transition structural differences cause the largest dynamic variability [64]. We show here that in a system with additional feedback from the dynamics to the structure, also the structure has the highest variability at the transition, which may have important consequences for collective computations, as it may for example enhance collective gradient sensing [55,65]. It shows that interactions on fixed [31,37,38] or randomly rewiring [30] lattices might miss this functionally highly relevant features of collective behavior.
The general structure of the assumed social interactions (short ranged repulsion, alignment and long range attraction) is supported by experiments [49,50]. However, in different species the detailed dependence of social interactions on relative positions may differ (see e.g. [50]). Here, to be as general as possible, we used simple functional forms of social interactions. However, the fundamental mechanisms underlying our results such as self-sorting and the structure-dynamic feedback will not depend on a more complex, empirically derived, relative position dependence. Neither should alternative interaction mechanisms affect these findings [66][67][68][69].
Our finding suggests that evolutionary adaptations at individual level are not a general mechanism for self-organization towards criticality. In principle, one could consider the possibility of multi-level selection [32,33] as a potential mechanism which could make the system evolve towards the group-level optimum at criticality. However, recent theoretical investigations of models of multi-level selection have shown that social dilemma, i.e. differences between ESSs and group level optima, always emerge for non-negligible individual-level selection even in cases where group-level selection strongly dominates [34,35]. Thus even in this biologically implausible scenario for fission-fusion prey schools, multi-level selection by its own appears unable to enforce evolutionary stability of the critical point in predator-prey dynamics.
We do not exclude the general possibility that animal collectives may operate in the vicinity of phase transitions in order to optimize collective computations. However, our results clearly demonstrate the necessity for further research on biologically proximate mechanisms of selforganized criticality in animal groups. A general, fundamental difficulty is that besides predator evasion there are various ecological contexts and other dimensions of (collective) behavior which will affect individual fitness. Here, by focusing on a dominant selection pressure, namely predation, we neglect other mechanisms, as for example resource exploration and exploitation [31,52,55,57] whose ESS can also depend on the resource abundance [31,52,57]. This emphasizes the importance to study collective behavior in the wild [44,[70][71][72] to provide more empirical input on actual relevant behavioral mechanisms as well as variability of behavior across different contexts. However, we have shown that even by combining two opposing selection mechanisms (see S1 Text, Sec. VII.3), which on their own favor ordered or disordered state respectively, the critical point does not correspond to an evolutionary attractor, it remains an evolutionary highly unstable point.
We focused here on the prominent directional symmetry breaking transition between states which are commonly observed in natural systems of collective behavior (disordered swarm, polarized school). Another possible transition involves the milling state [36], however, the function of the milling state in natural systems is unclear. Experiments suggest that boundary effects are a main reason for emergence of milling behavior in the laboratory [73], while milling in predator-prey interactions appears only to occur in the final stages of the hunt when the prey school is confined by multiple predators [74].
Recently it was suggested that a transition in the speed relaxation coefficient may represent a functionally relevant critical point in flocking behavior [19]. Individuals with lower relaxation constants are less bound to their preferred speed and may gain fitness benefits due 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 [51].
Fish also exhibit a reflex-driven escape response, so-called startle, which was recently shown to spreads through fish schools as a behavioral contagion process [75,76]. This suggests that at least in the context of collective predator evasion in fish, another type of a critical point may be highly relevant, which is analogous to the critical threshold in epidemic models. It separates states of non-propagating startle response, with only small localized response of single or few individuals, from avalanche-like dynamics, where a single fish may cause a global startle cascade. Even if the prey escape behavior is more complex, the self-sorting that happens before or in between predator attacks is unaffected by it and therefore also our results. Additionally, if a school is continuously pursued by predators, as e.g. in pelagic fish [77], the individual prey are likely to swim at their speed limit at which no further acceleration is possible.
Overall, our study does not reject the general possibility that animal groups manifest critical behavior and that it may be adaptive. However, it highlights importance of 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 ecological niche.

Methods
All Model parameters are listed in S1 Table.

Prey model
A prey agent i moves in 2D with constant velocity v = v 0 with directional noise of intensity D [78] and responds to a combined forceF i ¼F i;alg þF i;d þF i;flee by adapting its positionr i and heading φ i as dr i ðtÞ dt ¼ṽ i ðtÞ ð1aÞ with F i;? ðtÞ ¼F i ðtÞ �ẽ i;? as the combined force along the directionẽ i;? ¼ ½À sin φ i ; cos φ i � that is perpendicular to the agent's heading direction and ξ(t) as Gaussian white noise. The alignment force (F i;alg ) between a focal agent i and all its neighbors j 2 N i is the averaged velocity differenceṽ ji ¼ṽ j Àṽ i times the alignment strength μ alg . The distance regulating force (see S1 withr ji ¼ ðr j Àr i Þ=jr j Àr i j as direction from agent i to j, r d as preferred distance, μ d as strength of the force and m d as the slope 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 (F i;flee ) from it with a flee strength μ flee .

Predator-model
The predator moves with fixed speed v p = 2v 0 according to withF 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=jN p j). It only attacks one prey at a time. If the predator launches an attack, with an attack rate γ a (also accounting for handling time), its success probability decreases linear with distance and is zero for distances larger than r catch : In summary, the probability that a predator successfully catches a targeted agent within a small time window [t, t + δt] is The pursuit force, with constant magnitude μ p , points to a weighted center of mass. Each prey position is weighted by its probability of 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 simulation the γ 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 all agents.
(ii) The new generation of N offspring is generated via fitness-proportionate-selection. Thus, a random offspring has the parameters of the parent i with probability p parent,i = f i /∑ j f j .
(iii) An offspring mutates with probability γ m (mutation rate), by adding a Gaussian random variable with zero mean and standard deviation σ 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 offspring population (without mutation to reduce noise) and define the fitness gradient as the offspring mean parameter from which the current mean parameter is subtracted. Thus, if the offspring have a larger mean parameter, the fitness gradient is positive and vice versa. The mean fitness gradient of a certain parameter region is the average of generations within it. For details see S1 Text, Sec. III.

Quantification of collective behavior
The inter-individual distance is the distance between prey pairs averaged over all pairs IID ¼ hjr ij ji. The polarization F is the absolute value of the mean heading direction The susceptibility χ is the response of the polarization to an external field h and can be measured via polarization fluctuations (see S1 Text, Sec. V). It can be shown that Eq 6 is the same as the correlation of velocity fluctuations dṽ i ¼ṽ i À hṽi over all possible pairs (with hṽi ¼ P iṽi =N, see S1 Text, Sec. V). However, in inset of Fig 2B we computed the correlation of velocity fluctuations only over neighboring pairs Cðdṽ i ; dṽ j Þ ¼ P i;j2N i dṽ i � dṽ j because it is directly related to local 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 It is equal to the difference between the capture rates of non-fleeing and fleeing agents γ c,NF − γ c scaled by γ c,NF . The normalization of the capture difference by the baseline capture rate of non-fleeing prey γ c,NF accounts for potential differences in capture rates 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 and their mean relative location in the collective hr i,x i where x 2 {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 the relative locations see S1 Text, Sec. IV.1. Another, more general, quantification of self-sorting is how assortative the spatial arrangement of individuals with heterogeneous alignment is. We used the implementation of the assortativity coefficient [79] in igraph on the interaction network (Voronoi) with the values for each agent corresponding to their alignment strength (see S1 Text, Sec. IV for details). Illustration of angle-vector-relations for variables used in Eq S54 in S1 Text and the following. The black vectorṽ i is the current velocity of agent i. The blue vector hṽ j i N i is the mean velocity of the neighbors of agent i. The red vectorF flee represents the flee force experienced by agent i if the predator is sensed. The angle α is the angle between the mean velocity of neighbors and the velocity of agent i. The angle θ is the angle between the mean neighbor-velocity and the flee force. B: Numerical-results of the relative direction to neighbors α using Eq S54 in S1 Text. The initial condition is α = 0, i.e. the focal agent is perfectly aligned with its neighbors. The angle between mean neighbor velocity and flee force is θ = π/2. The different colors indicate that the effective flee direction, which is the compromise between the mean neighbor velocity and the flee-force, is faster the stronger the flee strength μ flee . We assumed, as discussed in S1 Text, Sec. VI, that the directional compromise represents the balance between benefits and costs of social information and is maintained, i.e. the prey evolve their alignment strength μ alg to keep the effective flee direction constant. population with a continuous blind angle (magenta dash dotted line), (B:) a less agile predator ("stiff") which turns less quick (black dotted line) and a more agile predator which turns quicker (red dashed line) than the predator in the standard case. (C:) a non-binarized fitness estimate (red dashed line) in which the preys fitness is not defined by captures but by the accumulated probability of capture, a fitness estimate based on captures during the simulation (black dotted line). (TIF) S5 Fig. Self-sorting with and without fixed speed. Self-sorting quantified via the Pearson correlation between the individual alignment parameter μ alg and the average relative position of the individuals (relative front-, side-or density-location as described in S1 Text, Sec. IV.1). A: If prey agents respond only by changing their direction but not their speed (fixed speed), selfsorting persists also in highly ordered regions. B: If prey agents can change their speed (variable speed), self-sorting vanishes for μ flee � 6. (TIF) S1 Movie. Animation of nine simulations without a predator below, at and above the phase transition. The red line are the past-and the empty red circle is the current center of mass of the collective. Animations in the same column are samples of the same parameter configuration. The columns differ in the alignment strength μ alg = [0, 1, 2] indicated at the top. The remaining parameters are identical to the ones used in the main text (listed in S1 Table). (MP4) S2 Movie. Animation of nine attack simulations below, at and above the phase transition. Same as S1 Movie but with a predator attacking the collective. (MP4) S3 Movie. Attack simulation on non-and fleeing prey. The left panel shows only the fleeing prey, the right the non-fleeing prey, and the center shows both. The color-code is black = fleeing prey, blue = non-fleeing prey, red = predator attacking fleeing prey, green = predator attacking non-fleeing prey. Parameters are identical to the ones used in the main text (listed in S1 Table). (MP4) S4 Movie. Animation of nine attack simulations above the phase transition. Same as S2 Movie but with other alignment parameters μ alg = [2,3,4]. (MP4) S5 Movie. Animation of nine attack simulations with variable prey speed. Same as S2 Movie but with preys that are able to accelerate according to the current force. The equations of motions for the prey with variable speed are defined in S1 Text, Sec. VII. (MP4) S1 Table. Default model parameters used. Time and space have been rescaled to dimensionless units by setting, without loss of generality, the prey speed v 0 and preferred distance r d to 1. All length scales are thus measured in units of r d , and all time scales in terms of time needed to move the distance r d . Note that the flee strength μ flee is strictly speaking a predator-prey parameter which reduces the prey-only parameters to four. (PDF)