Collective predator evasion: Putting the criticality hypothesis to the test

Abstract

[...]these observations could in principle have different origins [12, 22–24]. [...]more convincing support for the “criticality hypothesis” can be obtained through additional identification of proximate mechanisms enabling biological systems to self-organize towards criticality. [...]they assumed idealized random interaction networks between the agents. [...]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. [...]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]. [...]the observed evolutionary stable strategies (ESS) result from individual prey agents balancing the influence of social and private information on their movement response.

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–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–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–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–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 so-called 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–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–41]. This, in turn, has likely direct evolutionary consequences as for example predators may attack certain swarm regions more frequently [42–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–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–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.

Results

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 v0 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 rd (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).

[Figure omitted. See PDF.]

Fig 1. Social forces, interaction network and predator-prey interaction.

Schematic illustration of social forces (A), the Voronoi interaction network (B) and the predator prey interactions (C). (A) The social force acts on the focal agent (black triangle) and is a combination of alignment and distance regulation to its interaction partners. The alignment is proportional to the sum of the velocity differences with and thus not parallel to the neighbors mean velocity but tends to minimize the velocity difference. The distance regulating force Fd is a continuous version of a two zone model, i.e. the focal agent is repelled from neighbors that are closer () than the preferred distance rd () and attracted to those farther away (). (B) A focal prey agent (yellow triangle) interact with it’s nearest Voronoi neighbors (black triangles in yellow cells). (C) The predator () pursues the weighted mean direction of the targets (), 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 ≥ rcatch (). All Voronoi neighbors of the predator flee with a repulsive force ().

https://doi.org/10.1371/journal.pcbi.1008832.g001

The predator moves with a fixed speed vp which is larger than the preys (here vp = 2v0) 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 Nr 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 Nr 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.

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 Φ, 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 Φ 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 (Φ ≈ 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 Φ = 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).

[Figure omitted. See PDF.]

Fig 2. Group optimum.

Predation independent (A, B) and dependent (C, D) group measures. (A) Polarization Φ. The dashed vertical line marks the angular noise of D = 0.5 used in the evolutionary runs. (B) Directional information transfer , estimated via the correlation of velocity fluctuations between interacting agents, peaks at the transition. Inset: Susceptibility, estimated via polarization fluctuations. (C) Collective anti-predator performance quantified by the capture rate, which is strongly anti-correlated with the inter-individual distance R = −0.69 (IID, inset C). (D) Escape ratio Resc = 1 − γc/γc,NF. Inset: Difference between capture rates in schools of non-fleeing γc,NF and fleeing γc agents. In all panels: the disorder-order transition is indicated by the dash-dotted magenta line. Each parameter point corresponds to an average over Ns = 40 simulations, each with N = 400 agents attacked for Tsimu = 120 time units after an equilibration time of Teq = 200. For all insets(B, C, D): colorbars are shown separately in S7 Fig.

https://doi.org/10.1371/journal.pcbi.1008832.g002

A simple and intuitive measure of responsiveness of such a collective system to (local) perturbations is the average pair-wise correlation of velocity fluctuations between interacting agents (see Methods). Here, is the deviation of the velocity of agent i from the average school velocity . Cij 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 Cij indicates that agent j to a large degree is “copying” this velocity deviation or vice versa.

The velocity fluctuation correlation Cij 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, Cij 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 Resc, 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 (⟨μalg⟩ = 0), above (⟨μalg⟩ = 5) and far above (⟨μalg⟩ = 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.

[Figure omitted. See PDF.]

Fig 3. Evolution under predation.

(A) Overlay of three independent evolutionary runs starting at ⟨μalg⟩ = [0, 5, 10] over 1000 generations. The behavioral phenotype is determined only by the alignment strength as the evolving parameter. The predator attacks from random initial directions for Tsimu = 120. The inset shows the evolution of the population mean alignment parameter ⟨μalg⟩ of the three different evolutionary runs. (B) Assortativity coefficient () and smoothed fitness gradient ∇f (). The evolutionary stable state is defined by the zero crossing of the fitness gradient and represented as a vertical dashed black line. Black dots are the non-averaged fitness gradients for each generation (see Methods). (C) Self-sorting measured as correlation C(μalg, x) between the individual alignment strength μalg and variables quantifying its (spatial) location within the school: front-back position () and side-center position (black) and local density (). (D) Correlation C(f, x) of individual fitness with the average relative spatial positions. (E) Simulation snapshot illustrating the location variables: front-back position () and side-center position (black) and local density (). In all panels: the vertical dash-dotted magenta line marks the order-disorder transition and the vertical dashed black line the evolutionary stable state.

https://doi.org/10.1371/journal.pcbi.1008832.g003

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 Resc 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.

[Figure omitted. See PDF.]

Fig 4. Evolution for different flee strengths μflee.

(A) Sample evolutionary trajectories of the mean alignment strength μalg over 700 generations. (B) shows the dependence of evolutionary stable strategies (ESS) on the flee strength μflee. Solid diagonal line shows the theoretically predicted linear dependence of the ESS on μflee assuming balancing of social and private information (see SI Sec. VI). Dashed lines (A, B) connect the example individual evolutionary runs (A) to the corresponding ESSs (B) obtained as an average over multiple, longer evolutionary simulations. (C) Evolutionary stable states (circles) with respect to the group response, measured via the escape ratio Resc, for three selected flee-strengths indicated with dashed, solid and dotted lines for μflee = [2, 4, 8] respectively. In all panels: the dash-dotted magenta line marks the order-disorder transition and the different lines/points (red, black and blue) represent results for different flee strengths μflee = [2, 4, 8], respectively.

https://doi.org/10.1371/journal.pcbi.1008832.g004

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, model-independent 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–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–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 self-organized 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–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 = v0 with directional noise of intensity D [78] and responds to a combined force by adapting its position and heading φi as(1a)(1b)with as the combined force along the direction that is perpendicular to the agent’s heading direction and ξ(t) as Gaussian white noise. The alignment force () between a focal agent i and all its neighbors is the averaged velocity difference times the alignment strength μalg. The distance regulating force (see S1 Fig, panel A) is(2)with as direction from agent i to j, rd as preferred distance, μd as strength of the force and md as the slope of the change from repulsion (for rji < rd) to attraction (for rji > rd). If a predator p is a neighbor, the agent is repelled () from it with a flee strength μflee.

Predator-model

The predator moves with fixed speed vp = 2v0 according to(3)with as the pursuit force. It considers its frontal Voronoi-neighbors as targets and selects equally likely among them (). 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 rcatch:(4)In summary, the probability that a predator successfully catches a targeted agent within a small time window [t, t + δt] is(5)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 pcatch,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 Nf = 76 independent attack-simulations on the same prey population. For each simulation the γa⋅Ts agents with the largest cumulative pcatch are declared as dead. The fitness of agent i is fi = −Nk,i + max(Nk,j, j) with Nk,i as the number of simulations in which agent i was captured and max(Nk,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 pparent,i = fi/∑j fj.

(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 . The polarization Φ 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(6)(see S1 Text, Sec. V). It can be shown that Eq 6 is the same as the correlation of velocity fluctuations over all possible pairs (with , see S1 Text, Sec. V). However, in inset of Fig 2B we computed the correlation of velocity fluctuations only over neighboring pairs 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(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 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 ⟨ri,x⟩ 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 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).

Supporting information

S1 Text.

Sections: I. Model-Description. II. Model parameters. III. Evolutionary algorithm and ESS. IV. Measures of self-sorting. V. Susceptibility under a homogeneous global field. VI. Balancing social vs. direct predator information. VII. Robustness against modifications of the prey & predator dynamics and the selection mechanism.

https://doi.org/10.1371/journal.pcbi.1008832.s001

(PDF)

S1 Fig. Illustration of the distance regulating force and the relative polar coordinates.

A: Distance regulating force between agents i and j projected on the separation direction . The force equals zero at the preferred distance rd = 1 and is displayed for a distance regulating force steepness md = 2 (used in the simulations) and md = 4. B: Relative polar coordinates of an agent i with respect to the center of mass of the school () and to the average velocity of the school (). The angle αi,com () between the school velocity and the agents i current position () and the distance to the center of mass define the position in this relative coordinate system.

https://doi.org/10.1371/journal.pcbi.1008832.s002

(TIF)

S2 Fig. Fitness gradients for different flee-strength to estimate the ESSs.

Details on the estimation of evolutionary stable states of Fig 4 in the main text. A–G: Fitness gradient ∇f for evolution with different flee strength μflee. Black-dots indicate the estimated fitness gradients for each generation. Solid lines are averaged fitness gradients. Dashed vertical lines indicate where ∇f = 0 and thus mark the evolutionary stable states. H: All fitness gradients displayed together. Note that the peaks for μflee = 6 at μalg ≈ 3 and for μflee = 8 at μalg ≈ 4 are due to fluctuations in the standard-deviation of the population. If the standard-deviation is kept constant those peaks vanish (not shown).

https://doi.org/10.1371/journal.pcbi.1008832.s003

(TIF)

S3 Fig. Balancing social and private information via a directional compromise.

A: Illustration of angle-vector-relations for variables used in Eq S54 in S1 Text and the following. The black vector is the current velocity of agent i. The blue vector is the mean velocity of the neighbors of agent i. The red vector 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.

https://doi.org/10.1371/journal.pcbi.1008832.s004

(TIF)

S4 Fig. Robustness analysis of evolution results.

Evolutionary stable states of the alignment strength are estimated from the fitness gradient for different flee strength under slight variations of simulations parameters or predator attack implementation. The standard scenario of the main text () is compared to (A:) a prey population with varying speed which can avoid the predator additionally by acceleration (), a prey population with a angular diffusion coefficient which is doubled compared to the standard case (), a prey population with a continuous blind angle (), (B:) a less agile predator (“stiff”) which turns less quick () and a more agile predator which turns quicker () than the predator in the standard case. (C:) a non-binarized fitness estimate () 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 ().

https://doi.org/10.1371/journal.pcbi.1008832.s005

(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), self-sorting persists also in highly ordered regions. B: If prey agents can change their speed (variable speed), self-sorting vanishes for μflee ≤ 6.

https://doi.org/10.1371/journal.pcbi.1008832.s006

(TIF)

S6 Fig. 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 frontal individuals are declared as dead. The relative strength of the frontal-risk selection is defined by the ratio between agents killed at the front and by the predator, i.e. (Front Kills)/(Pred. Kills)∈[0, 0.05, 0.1, 0.2]. The evolutionary stable state (ESS) is defined by the zero-crossing of the fitness gradient with negative slope marked by a vertical dashed line. However, the lower bound is an additional ESS if the fitness gradient stays negative close to it which is marked by shaded points in the inset. Parameters are identical to the former simulations apart from the angular diffusion coefficient which is increased to D = 1 increasing the order-transition to μalg,c ≈ 1.6 marked by vertical dash-dotted magenta line. The flee strength is μflee = 4.

https://doi.org/10.1371/journal.pcbi.1008832.s007

(TIF)

S7 Fig. Enlarged insets from main text Fig 2.

The susceptibility χ (A), inter individual distance IID (B) and difference in capture rate between non-fleeing and fleeing individuals γc,NF − γc (C). All measures are shown with colorbars, which were omitted for clarity in Fig 2.

https://doi.org/10.1371/journal.pcbi.1008832.s008

(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).

https://doi.org/10.1371/journal.pcbi.1008832.s009

(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.

https://doi.org/10.1371/journal.pcbi.1008832.s010

(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).

https://doi.org/10.1371/journal.pcbi.1008832.s011

(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].

https://doi.org/10.1371/journal.pcbi.1008832.s012

(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.

https://doi.org/10.1371/journal.pcbi.1008832.s013

(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 v0 and preferred distance rd to 1. All length scales are thus measured in units of rd, and all time scales in terms of time needed to move the distance rd. Note that the flee strength μflee is strictly speaking a predator-prey parameter which reduces the prey-only parameters to four.

https://doi.org/10.1371/journal.pcbi.1008832.s014

(PDF)

Acknowledgments

We are grateful to Iain Couzin, Simon Levin, Jessica Flack, and David Krakauer for insightful and inspiring discussions leading to this work. We further acknowledge Ishan Levy for suggesting the investigation of non-responding prey dynamics.

Citation: Klamser PP, Romanczuk P (2021) Collective predator evasion: Putting the criticality hypothesis to the test. PLoS Comput Biol 17(3): e1008832. https://doi.org/10.1371/journal.pcbi.1008832

References

1. Rubinov M, Sporns O. Complex network measures of brain connectivity: Uses and interpretations. NeuroImage. 2010;52(3):1059–1069. pmid:19819337

2. Vijesh N, Chakrabarti SK, Sreekumar J. Modeling of gene regulatory networks: A review. Journal of Biomedical Science and Engineering. 2013;06(02):223–231.

3. Miller N, Garnier S, Hartnett AT, Couzin ID. Both information and social cohesion determine collective decisions in animal groups. Proceedings of the National Academy of Sciences. 2013;110(13):5263–5268. pmid:23440218

4. Strandburg-Peshkin A, Twomey CR, Bode NWFF, Kao AB, Katz Y, Ioannou CC, et al. Visual sensory networks and effective information transfer in animal groups. Current Biology. 2013;23(17):R709–R711. pmid:24028946

5. Packard NH. Adaptation Toward the Edge of Chaos. In: Kelso JAS, Mandell AJ, Shlesinger MF, editors. Dynamic Patterns in Complex Systems. Singapore, World Scientific; 1988.

6. Bak P, Chen K, Creutz M. Self-organized criticality in the’Game of Life’. Nature. 1989;342:780–782.

7. Langton CG. Computation at the edge of chaos: Phase transitions and emergent computation. Physica D. 1990;42:12–37.

8. Bak P, Sneppen K. Punctuated Equilibribum and Criticality in a simple model of evolution. Physical Review Letters. 1993;71(24):4083–4086.

9. Kinouchi O, Copelli M. Optimal dynamical range of excitable networks at criticality. Nature Physics. 2006;2(5):348–351.

10. Mora T, Bialek W. Are Biological Systems Poised at Criticality? Journal of Statistical Physics. 2011;144(2):268–302.

11. Hidalgo J, Grilli J, Suweis S, Munoz MA, Banavar JR, Maritan A. Information-based fitness and the emergence of criticality in living systems. Proceedings of the National Academy of Sciences. 2014;111(28):10095–10100. pmid:24982145

12. Beggs JM, Timme N. Being critical of criticality in the brain. Frontiers in Physiology. 2012;3 JUN(June):1–14.

13. Muñoz MA. Colloquium: Criticality and dynamical scaling in living systems. Reviews of Modern Physics. 2018;90(3):31001.

14. Krotov D, Dubuis JO, Gregor T, Bialek W. Morphogenesis at criticality. Proceedings of the National Academy of Sciences. 2014;111(10):3683–3688. pmid:24516161

15. Daniels BC, Kim H, Moore D, Zhou S, Smith HB, Karas B, et al. Criticality distinguishes the ensemble of biological regulatory networks. Physical review letters. 2018;121(13):138102. pmid:30312104

16. Gelblum A, Pinkoviezky I, Fonio E, Ghosh A, Gov N, Feinerman O. Ant groups optimally amplify the effect of transiently informed individuals. Nature Communications. 2015;6. pmid:26218613

17. Feinerman O, Pinkoviezky I, Gelblum A, Fonio E, Gov NS. The physics of cooperative transport in groups of ants. Nature Physics. 2018;14(7):1–11.

18. Attanasi A, Cavagna A, Del Castello L, Giardina I, Melillo S, Parisi L, et al. Finite-size scaling as a way to probe near-criticality in natural swarms. Physical Review Letters. 2014;113(23):238102. pmid:25526161

19. Bialek W, Cavagna A, Giardina I, Mora T, Pohl O, Silvestri E, et al. Social interactions dominate speed control in poising natural flocks near criticality. Proceedings of the National Academy of Sciences of the United States of America. 2014;111(20):7212–7217. pmid:24785504

20. Daniels BC, Krakauer DC, Flack JC. Control of finite critical behaviour in a small-scale social system. Nature Communications. 2017;8:1–8.

21. Ginelli F, Peruani F, Pillot MH, Chaté H, Theraulaz G, Bon R. Intermittent collective dynamics emerge from conflicting imperatives in sheep herds. Proc Natl Acad Sci. 2015;112(41):12729–12734. pmid:26417082

22. Lorimer T, Gomez F, Stoop R. Two universal physical principles shape the power-law statistics of real-world networks. Sci Rep. 2015;5:1–8. pmid:26202858

23. Reed WJ, Hughes BD. From gene families and genera to incomes and internet file sizes: Why power laws are so common in nature. Phys Rev E—Stat Physics, Plasmas, Fluids, Relat Interdiscip Top. 2002;66(6):4. pmid:12513446

24. Touboul J, Destexhe A. Power-law statistics and universal scaling in the absence of criticality. Physical Review E. 2017;95(1):012413. pmid:28208383

25. Bornholdt S, Rohlf T. Topological evolution of dynamical networks: Global criticality from local dynamics. Physical Review Letters. 2000;84(26):6114. pmid:10991137

26. Meisel C, Gross T. Adaptive self-organization in a realistic neural network model. Physical Review E. 2009;80(6):061917.

27. Ma Z, Turrigiano GG, Wessel R, Hengen KB. Cortical circuit dynamics are homeostatically tuned to criticality in vivo. Neuron. 2019;104(4):655–664. pmid:31601510

28. Liu M, Bassler KE. Emergent criticality from coevolution in random boolean networks. Physical Review E. 2006;74(4):041910. pmid:17155099

29. MacArthur BD, Sánchez-García RJ, Ma’ayan A. Microdynamics and criticality of adaptive regulatory networks. Physical review letters. 2010;104(16):168701. pmid:20482087

30. Torney CJ, Lorenzi T, Couzin ID, Levin SA. Social information use and the evolution of unresponsiveness in collective systems. Journal of the Royal Society Interface. 2015;12(103). pmid:25519991

31. Brush ER, Leonard NE, Levin SA. The content and availability of information affects the evolution of social-information gathering strategies. Theoretical Ecology. 2016;9(4):455–476.

32. Wilson DS. A theory of group selection. Proc Natl Acad Sci U S A. 1975;72(1):143–146. pmid:1054490

33. Wilson DS. Altruism And Organism: Disentangling The Themes Of Multilevel Selection Theory. Am Nat. 1997;150(S1):S122–S134. pmid:18811309

34. Cooney DB. The replicator dynamics for multilevel selection in evolutionary games. Journal of mathematical biology. 2019;79(1):101–154. pmid:30963211

35. Cooney DB. Analysis of multilevel replicator dynamics for general two-strategy social dilemma. Bulletin of Mathematical Biology. 2020;82:1–72. pmid:32474720

36. Calovi DS, Lopez U, Schuhmacher P, Chate H, Sire C, Theraulaz G. Collective response to perturbations in a data-driven fish school model. Journal of The Royal Society Interface. 2015;12(104):20141362–20141362.

37. Vanni F, Luković M, Grigolini P. Criticality and Transmission of Information in a Swarm of Cooperative Units. Physical Review Letters. 2011;107(7):078103. pmid:21902433

38. Chicoli A, Paley DA. Probabilistic information transmission in a network of coupled oscillators reveals speed-accuracy trade-off in responding to threats. Chaos An Interdiscip J Nonlinear Sci. 2016;26(11):116311.

39. Couzin ID, Krause J, James R, Ruxton GD, Franks NR. Collective memory and spatial sorting in animal groups. Journal of theoretical biology. 2002;218(1):1–11. pmid:12297066

40. Hemelrijk CK, Kunz H. Density distribution and size sorting in fish schools: An individual-based model. Behavioral Ecology. 2005;16(1):178–187.

41. Jamie Wood A. Strategy selection under predation; evolutionary analysis of the emergence of cohesive aggregations. Journal of Theoretical Biology. 2010;264(4):1102–1110. pmid:20346960

42. Krause J. Differential Fitness Returns in Relation To Spatial Position in Groups. Biological Reviews. 1994;69(2):187–206. pmid:8054444

43. Bumann D, Rubenstein D, Krause J. Mortality Risk of Spatial Positions in Animal Groups: the Danger of Being in the Front. Behaviour. 1997;134(13-14):1063–1076.

44. Handegard NO, Boswell KM, Ioannou CC, Leblanc SP, Tjøstheim DB, Couzin ID. The Dynamics of Coordinated Group Hunting and Collective Information Transfer among Schooling Prey. Curr Biol. 2012;22(13):1213–1217. pmid:22683262

45. Vicsek T, Czirk A, Ben-Jacob E, Cohen I, Shochet O. Novel type of phase transition in a system of self-driven particles. Physical Review Letters. 1995;75(6):1226–1229. pmid:10060237

46. Großmann R, Schimansky-Geier L, Romanczuk P. Active Brownian particles with velocity-alignment and active fluctuations. New Journal of Physics. 2012;14.

47. Krause J, Ruxton GD. Living in Groups. Oxford Univ. Press, Oxford; 2002.

48. Harpaz R, Tkačik G, Schneidman E. Discrete modes of social information processing predict individual behavior of fish in a group. Proceedings of the National Academy of Sciences. 2017;114(38):10149–10154.

49. Katz Y, Tunstrøm K, Ioannou CC, Huepe C, Couzin ID. Inferring the structure and dynamics of interactions in schooling fish. Proceedings of the National Academy of Sciences of the United States of America. 2011;108(46):18720–18725. pmid:21795604

50. Calovi DS, Litchinko A, Lecheval V, Lopez U, Pérez Escudero A, Chaté H, et al. Disentangling and modeling interactions in fish with burst-and-coast swimming reveal distinct alignment and attraction behaviors. PLOS Computational Biology. 2018;14(1):e1005933. pmid:29324853

51. Herbert-Read JE, Rosén E, Szorkovszky A, Ioannou CC, Rogell B, Perna A, et al. How predation shapes the social interaction rules of shoaling fish. Proceedings of the Royal Society B: Biological Sciences. 2017;284(1861):20171126. pmid:28855361

52. Wood AJ, Ackland GJ. Evolving the selfish herd: Emergence of distinct aggregating strategies in an individual-based model. Proceedings of the Royal Society B: Biological Sciences. 2007;274(1618):1637–1642. pmid:17472913

53. Olson RS, Hintze A, Dyer FC, Knoester DB, Adami C. Predator confusion is sufficient to evolve swarming behaviour. Journal of The Royal Society Interface. 2013;10(85):20130305–20130305. pmid:23740485

54. Olson RS, Knoester DB, Adami C. Evolution of Swarming Behavior Is Shaped by How Predators Attack. Artificial Life. 2016;22(3):299–318. pmid:27139941

55. Hein AM, Rosenthal SB, Hagstrom GI, Berdahl A, Torney CJ, Couzin ID. The evolution of distributed sensing and collective computation in animal populations. eLife. 2015;4(DECEMBER2015):1–43. pmid:26652003

56. Guttal V, Couzin ID. Social interactions, information use, and the evolution of collective migration. Proceedings of the National Academy of Sciences of the United States of America. 2010;107(37):16172–16177. pmid:20713700

57. Monk CT, Barbier M, Romanczuk P, Watson JR, Alós J, Nakayama S, et al. How ecology shapes exploitation: a framework to predict the behavioural response of human and animal foragers along exploration–exploitation trade-offs. Ecology Letters. 2018;21(6):779–793. pmid:29611278

58. Van Der Post DJ, Verbrugge R, Hemelrijk CK. The evolution of different forms of sociality: Behavioral mechanisms and eco-evolutionary feedback. PLoS ONE. 2015;10(1):1–19. pmid:25629313

59. Ballerini M, Cabibbo N, Candelier R, Cavagna A, Cisbani E, Giardina I, et al. Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study. Proceedings of the National Academy of Sciences. 2008;105(4):1232–1237. pmid:18227508

60. Guttal V, Romanczuk P, Simpson SJ, Sword GA, Couzin ID. Cannibalism can drive the evolution of behavioural phase polyphenism in locusts. Ecology letters. 2012;15(10):1158–1166. pmid:22882379

61. Lipowski A, Lipowska D. Roulette-wheel selection via stochastic acceptance. Physica A: Statistical Mechanics and its Applications. 2012;391(6):2193–2196.

62. Cavagna A, Giardina I, Grigera TS. The physics of flocking: Correlation as a compass from experiments to theory. Phys Rep. 2018;728:1–62.

63. James R, Bennett P, Krause J. Geometry for mutualistic and selfish herds: the limited domain of danger. Journal of Theoretical Biology. 2004;228(1):107–113. pmid:15064086

64. Nykter M, Price ND, Larjo A, Aho T, Kauffman SA, Yli-Harja O, et al. Critical networks exhibit maximal information diversity in structure-dynamics relationships. Phys Rev Lett. 2008;100(5):1–4. pmid:18352443

65. Berdahl A, Torney CJ, Ioannou CC, Faria JJ, Couzin ID. Emergent Sensing of Complex Environments by Mobile Animal Groups. Science. 2013;339(6119):574–576. pmid:23372013

66. Bastien R, Romanczuk P. A model of collective behavior based purely on vision. Sci Adv. 2020;6(6):1–10.

67. Romanczuk P, Couzin ID, Schimansky-Geier L. Collective Motion due to Individual Escape and Pursuit Response. Physical Review Letters. 2009;102(1):010602. pmid:19257176

68. Jhawar J, Morris RG, Amith-Kumar UR, Danny Raj M, Rogers T, Rajendran H, et al. Noise-induced schooling of fish. Nature Physics. 2020;16(4):488–493.

69. Lei L, Escobedo R, Sire C, Theraulaz G. Computational and robotic modeling reveal parsimonious combinations of interactions between individuals in schooling fish. PLOS Computational Biology. 2020;16(3):e1007194. pmid:32176680

70. Francisco FA, Nührenberg P, Jordan A. High-resolution, non-invasive animal tracking and reconstruction of local environment in aquatic ecosystems. Movement Ecology. 2020;8(1):27. pmid:32582448

71. Hansen MJ, Krause S, Breuker M, Kurvers RHJM, Dhellemmes F, Viblanc PE, et al. Linking hunting weaponry to attack strategies in sailfish and striped marlin. Proceedings of the Royal Society B: Biological Sciences. 2020;287 (1918). pmid:31937224

72. Graving JM, Chae D, Naik H, Li L, Koger B, Costelloe BR, et al. Deepposekit, a software toolkit for fast and robust animal pose estimation using deep learning. eLife. 2019;8:1–42. pmid:31570119

73. Tunstrøm K, Katz Y, Ioannou CC, Huepe C, Lutz MJ, Couzin ID. Collective States, Multistability and Transitional Behavior in Schooling Fish. PLoS Comput Biol. 2013;9(2):e1002915. pmid:23468605

74. Wickham DA, Russel GM. Evaluation of mid-water artificial structures. Fish Bull. 1974;72(1):181–191.

75. Sosna MMG, Twomey CR, Bak-Coleman J, Poel W, Daniels BC, Romanczuk P, et al. Individual and collective encoding of risk in animal groups. Proceedings of the National Academy of Sciences of the United States of America. 2019;116(41):20556–20561. pmid:31548427

76. Rosenthal SB, Twomey CR, Hartnett AT, Wu HS, Couzin ID. Revealing the 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. 2015;112(15):4690–4695. pmid:25825752

77. Herbert-Read JE, Romanczuk P, Krause S, Strömbom D, Couillaud P, Domenici P, et al. Proto-cooperation: group hunting sailfish improve hunting success by alternating attacks on grouping prey. Proc R Soc B Biol Sci. 2016;283(1842):20161671.

78. Romanczuk P, Schimansky-Geier L. Brownian Motion with Active Fluctuations. Physical Review Letters. 2011;106(23):230601. pmid:21770491

79. Newman MEJ. Mixing patterns in networks. Physical Review E. 2002;67(2):026126. pmid:12443515