Abstract
Many animal groups exhibit rapid, coordinated collective motion. Yet, the evolutionary forces that cause such collective responses to evolve are poorly understood. Here, we develop analytical methods and evolutionary simulations based on experimental data from schooling fish. We use these methods to investigate how populations evolve within unpredictable, timevarying resource environments. We show that populations evolve toward a distinctive regime in behavioral phenotype space, where small responses of individuals to local environmental cues cause spontaneous changes in the collective state of groups. These changes resemble phase transitions in physical systems. Through these transitions, individuals evolve the emergent capacity to sense and respond to resource gradients (i.e. individuals perceive gradients via social interactions, rather than sensing gradients directly), and to allocate themselves among distinct, distant resource patches. Our results yield new insight into how natural selection, acting on selfish individuals, results in the highly effective collective responses evident in nature.
https://doi.org/10.7554/eLife.10955.001eLife digest
In nature, we see many examples of highly coordinated movements of groups of individuals; think of a flock of birds turning swiftly in unison or a crowd of people filing through the exit of a building. A common feature of these behaviors is that they occur without any centralized control, and that they involve sudden and often dramatic changes in the 'collective state' of the group (i.e. speed, or the distances between individuals). In the past, researchers have likened these transitions in collective behavior to phase transitions in physical systems, for example, the transition between liquid water and water vapor. However, it is not clear how such collective responses could have evolved.
Natural selection is an evolutionary process whereby individuals with particularly 'fit' traits produce more offspring than others. Over many generations, these beneficial traits tend to become more common in the population. Hein, Rosenthal, Hagstrom et al. developed a mathematical model to investigate whether the capacity of a population to perform collective motions could evolve through natural selection.
The model shows that over many generations, populations consistently evolve a unique collective trait whereby small responses of individuals to an environmental cue can cause spontaneous changes in the collective state of the local population. These transitions in collective state greatly enhance the ability of individuals to locate and exploit resources. Hein, Rosenthal, Hagstrom et al.’s findings suggest that natural selection acting on the behavior of individuals can cause a population to evolve a distinctive, collective behavior.
The next challenge will be to identify a biological system in which the evolution of collective motion can be studied experimentally to test these predictions.
https://doi.org/10.7554/eLife.10955.002Introduction
In many highly coordinated animal groups, such as fish schools and bird flocks, the ability of individuals to locate resources and avoid predators depends on the collective behavior of the group. For example, when fish schools are attacked by predators, 'flash expansion' (Pitcher et al., 1993) and other coordinated collective motions, made possible above a certain group size, reduce individual risk (Handegard et al., 2012). Similarly, fish can track dynamic resource patches far more effectively when they are in a group (Berdahl et al., 2013). When an individual responds to a change in the environment (e.g., predator, resource cue), this response propagates swiftly through the group (Rosenthal et al., 2015), altering the group’s collective motion. How are such rapid, coordinated responses possible? These responses may occur, in part, because the nature of social interactions makes animal groups highly sensitive to small changes in the behavior of individual group members; theoretical (Couzin et al., 2002; D’Orsogna et al., 2006; Kolpas et al., 2007) and empirical (Tunstrøm et al., 2013; Buhl et al., 2006) studies of collective motion have revealed that minor changes in individual behavior, such as speed (Tunstrøm et al., 2013), can cause sudden transitions in group state, reminiscent of similarly sudden phase transitions between collective states in physical systems (such as the solidliquidgas transitions as a function of increasing temperature). It has been proposed that individuals may trigger such changes in collective state by responding to the environment, thereby initiating a coordinated response at the group level (e.g., Couzin et al. (2002); Kolpas et al. (2007); Couzin and Krause, 2003). This mechanism requires that the behavioral rules of individual animals within a population have evolved in a way that allows groups to transition adaptively among distinct collective states. The evolutionary processes that could lead to this populationlevel property, however, remain poorly understood.
The feedback between the behavioral phenotypes of individuals, the collective behaviors that these phenotypes produce, and individuallevel fitness consequences has made it challenging to study how complex collective behaviors evolve (Torney et al., 2011). Many species, including fish and birds, form groups in which members have low genetic relatedness, which implies that kin selection alone cannot explain the evolution of collective behavior. Moreover, while natural selection acts on the behavioral phenotypes of selfish individuals, collective behaviors are grouplevel, or perhaps even populationlevel, properties rather than heritable individual phenotypes. To understand how collective behaviors evolve, then, one must first understand the mapping between individual phenotypes and collective behavior, and between collective behavior and individual fitness.
Here, we take advantage of detailed studies of the social interaction rules and environmental response behaviors of schooling fish (Berdahl et al., 2013; Katz et al., 2011) to develop a biologicallymotivated evolutionary model of collective responses to the environment. Using analytical methods and evolutionary simulations, we study how individual behavioral rules produce collective behaviors, and how collective behaviors, in turn, govern the fitness and evolution of selfish individuals. To relate individual and collective behaviors to fitness, we consider a fundamental task faced by fish and other motile organisms: finding and exploiting dynamic resources (Stephens et al., 2007). In our model, individuals respond to the locations of near neighbors and also to local measurements of resource quality. Each individual achieves a fitness determined by the resource level it experiences over its lifetime. We use this framework to explore the evolution of complex collective responses to the environment, and how such responses are related to transitions in collective state.
Model development
Behavioral rules
We model the movement behaviors of each individual in a population of size $N$ using two experimentallymotivated (Berdahl et al., 2013; Katz et al., 2011) behavioral rules: a social response rule and an environmental response rule. The social response rule is motivated by experimental studies of pairwise interactions among golden shiners (Notemigonus crysoleucas) (Katz et al., 2011). Individual fish avoid others with whom they are in very close proximity. As the distance between individuals increases, however, interactions gradually change from repulsive to attractive, with maximum attraction occurring at a distance of twofour body lengths. For longer distances, individuals still attract one another but the strength of attraction decays in magnitude (Appendix section 1; Katz et al., 2011). As found in experimental studies of golden shiners (Katz et al., 2011) and mosquitofish (Gambusia holbrooki) (HerbertRead et al., 2011) there need not be an explicit alignment tendency; rather alignment can be an emergent property of motion combined with the tendencies for repulsion and attraction described above.
To capture these observed social interactions (or ‘social forces’), we model the acceleration of individuals using a forcebased method (Katz et al., 2011). The $i$th individual responds to its neighbors using the following rule:
where ${F}_{s\mathrm{,}i}$ is the social force on the $i$th individual, ${x}_{i}$ is the position of the $i$th individual, $\nabla $ is the twodimensional gradient operator, the term in brackets is a social potential, ${C}_{a}$, ${C}_{r}$, ${l}_{a}$, and ${l}_{r}$ are constants that dictate the relative strengths and length scales of social attraction and repulsion, and the set ${N}_{i}$ is a set of the $k$ nearest neighbors of the $i$th individual, where a neighbor is an individual within a distance of ${l}_{max}$ of the focal individual. Equation 1 does not include explicit alignment with neighbors. A similar model is discussed in D’Orsogna et al. (2006). In Equation 1, ${l}_{max}$ determines the length scale over which individuals are influenced by social interactions. If ${l}_{max}$ is greater than ${l}_{r}$ but less than ${l}_{a}$, individuals repel one another at short distances but do not attract one another. We refer to such individuals as asocial (Appendix section 1). If ${l}_{max}$ is greater than both ${l}_{r}$ and ${l}_{a}$, individuals repel one another at short distances and are attracted to one another at intermediate distances as observed by Katz et al. (2011). Finite $k$ ensures that individuals can only respond to a limited number of their neighbors in crowded regions of space and provides a simplified model of sensorybased social interactions (e.g., Rosenthal et al. (2015); StrandburgPeshkin et al. (2013)). Finite $k$ also ensures that individuals are limited to finite local density (Appendix section 3).
To model the response of individuals to the environment, we develop an environmental response rule based on experimentallyobserved environmental responses of golden shiners (Berdahl et al., 2013). In particular, in a dynamic, heterogeneous environment, individual golden shiners respond strongly to local sensory cues by slowing down in favorable regions of the environment, and speeding up in unfavorable regions. In contrast, fish respond only weakly to spatial gradients in environmental quality and instead adjust their headings primarily based on the positions of their near neighbors. Accordingly, we model the $i$th individual’s environmental response as a function of the level of an environmental cue (in this case, the level of a resource) at its current position:
where ${F}_{a\mathrm{,}i}$ is the autonomous force the $i$th individual generates by accelerating or decelerating in response to the environment, $\Psi (\cdot )$ is a monotonically decreasing function of the value of an environmental cue, $S\left({x}_{i}\right)$ is the cue value at the $i$th individual’s position, $\eta $ is a damping term that limits individuals to a finite speed, and ${v}_{i}$ is the $i$th individual’s velocity. In the absence of social interactions, individuals travel at preferred speed ${v}_{i}^{*}=\sqrt{{\Psi}_{i}/\eta}$ (for ${\Psi}_{i}>0$). Changes in speed are crucial in the schooling behavior of fish (Tunstrøm et al., 2013; Berdahl et al., 2013), and as we show below, are also responsible for generating effective collective response in our model. Following the experimental results in Berdahl et al. (2013) we assume that individuals do not change their headings in response to the cue. In what follows, we refer to 'cue' and 'resource' interchangeably as we model the case where the cue is the resource itself (see e.g., Torney et al. (2009); Hein and McKinley (2012) for cases where the cue is not a resource).
Combining social and environmental response rules yields two equations that govern each individual’s movement (in two dimensions):
and
where $m$ is mass. D’Orsogna et al. (2006) explores the behavior of a similar model with ${\Psi}_{i}=\Psi $ constant over the full parameter space. Here we focus on a parameter regime that yields behavioral rules that match the experimental observations of Katz et al. (2011) and Berdahl et al. (2013).
We simulate a discretized version of the system described by Equations 3 and 4. In particular, we choose a time step, $\tau $, within which the acceleration due to social influences (Equation 1) and resource value $S\left({x}_{i}\right)$ are assumed to be constant. Positions, speeds, and accelerations of all individuals at time $t+\tau $ are then given by the solutions to Equations 3 and 4 at time $t+\tau $, with the values of $S\left({x}_{i}\right)$ and ${x}_{i}{x}_{j}$ determined at time $t$. A navigational noise vector of small magnitude $\gamma $ and uniform heading 0 to 2$\pi $ is added to the velocity of each agent at each time step. Taking the limit as $\tau $ goes to zero means that individuals are constantly acquiring information and instantaneously altering their actions in response. In Appendix section 3−6, we analyze a continuum approximation of this limiting model and below we discuss results of this analysis alongside simulation results.
The social interaction rule allows us to build an interaction network for the entire population. Two individuals are socially connected if at least one of them influences the other through Equation 1. We define a 'group' as a set of individuals that belong to the same connected component in this network.
Evolutionary dynamics
The natural environments in which organisms live are often heterogeneous and dynamic (Stephens et al., 2007). Consequently, we simulate populations of individuals in dynamic landscapes, where individuals make decisions in response to local sensory cues (local measurements of a resource) and these decisions have fitness consequences for the individuals within the population (Guttal and Couzin, 2010; Torney et al., 2011). In keeping with experimental observations (Berdahl et al., 2013), we assume individuals follow a simple environmental response function: ${\Psi}_{i}={\psi}_{0}{\psi}_{1}S\left({x}_{i}\right)$, where ${\psi}_{0}$ dictates the $i$th individual’s preferred speed when the level of the environmental cue is zero and ${\psi}_{1}$ determines how sensitive the $i$th individual is to the cue value (Berdahl et al., 2013). Rather than prescribing values of ${\psi}_{0}$ and ${\psi}_{1}$, we use an evolutionary framework similar to that developed by Guttal and Couzin (2010) to allow these two behavioral traits to evolve along with the maximum interaction length ${l}_{max}$, which determines whether individuals are social (${l}_{max}>$ length scale of social attraction) or asocial (${l}_{max}<$ length scale of social attraction, Appendix section 1).
In each generation, $N$ individuals are located in a twodimensional environment in which each point in space is associated with a resource value that changes over time (see Materials and methods). Individuals move through the environment using the interaction rules described above, and each individual has its own value of the ${\psi}_{0}$, ${\psi}_{1}$, and ${l}_{max}$ parameters. At the end of each generation, we compute each individual’s fitness as the mean value of the resource it experienced during that generation. Each individual then reproduces with a probability proportional to its relative fitness within the population. $N$ offspring comprise the next generation where each offspring inherits the traits of its parent modified by a small mutation (Appendix section 2). For reference, we compare the evolution of populations in which ${\psi}_{0}$, ${\psi}_{1}$, and ${l}_{max}$ are allowed to evolve, to the evolution of populations of asocial individuals, for which ${l}_{max}$ is set to a constant (Appendix section 1).
Results
Evolution of behavioral rules
In populations of asocial individuals, the baseline speed parameter and environmental sensitivity increase consistently through evolutionary time (Figure 1A–B). Asocial individuals move through the environment, slowing down in regions where the resource value is high and speeding up when the resource value is low (Video 1). As one would expect from random walk theory (Schnitzer, 1993; Gurarie and Ovaskainen, 2013), individuals more rapidly encounter regions of the environment with high resource value when they travel at high preferred speeds (Equation A65; Gurarie and Ovaskainen, 2013), and the more they reduce speed in regions of the environment with high resource quality, the more time they spend in these regions (Schnitzer, 1993). Because of these two effects, the fittest asocial individuals have high baseline speeds (i.e., high ${\psi}_{0}$) and accelerate and decelerate rapidly in response to changes in the resource value (i.e., high ${\psi}_{1}$; Figure 1A–B, Appendix).
When populations are allowed to evolve sociality, the evolutionary process selects for very different behaviors (Figure 1C–E). Selection quickly favors sociality, and individuals evolve large maximum interaction lengths (Figure 1C). Over evolutionary time, selection removes individuals with high and low values of ${\psi}_{0}$ and ${\psi}_{1}$ from the population and an evolutionarily stable state (ESSt; Maynard Smith, 1982) emerges that is characterized by a single mode at the dominant value of each trait (Figure 1D–E; Appendix section 2). The ESSt resulting from selection on ${\psi}_{0}$, ${\psi}_{1}$, and ${l}_{max}$ is robust in that it is resistant to invasion by phenotypes near the ESSt, and by invaders with trait values far from the ESSt (Appendix section 2). Throughout evolution, populations of social individuals achieve mean fitness values that are approximately five times higher than those of asocial populations, and a coefficient of variation in fitness approximately four times lower than that of asocial individuals (Figure 1F).
Notably, a single individual drawn from a population at the ESSt can invade a resident population of asocial individuals and the social strategy quickly sweeps through the population (Appendix section 2). To understand why this invasion occurs, consider a population of asocial individuals that slow down in favorable regions of the environment. If the environment does not change too rapidly, such individuals will accumulate in regions where the resource level is high. This phenomenon has been studied mathematically in the context of positiondependent diffusion (Schnitzer, 1993), and will occur, in general, when individuals lower their speeds in response to the value of an environmental cue. A social mutant that responds to the environment, and to its neighbors, can take advantage of the correlation between density and resource quality by climbing the gradient in the density of its neighbors (Equation 1). In this case, the positions of neighbors contain information about the value of resources and social mutants quickly invade asocial populations leading to a rapid increase in mean fitness (Appendix section 2).
Evolved populations collectively compute properties of the environment
The high fitness of the evolved phenotype is due, in part, to a collective resource tracking ability, similar to that found in golden shiners (Berdahl et al., 2013). Evolved individuals can find and track resource peaks as they move through the environment (Figure 2A, Video 2; Materials and methods), whereas asocial individuals and social individuals with trait values far from the ESSt cannot (Videos 1, 3–4). Tracking occurs via a dynamic process. Individuals near the edge of the peak move rapidly, whereas individuals nearer to the peak center (where the resource value is high) move slowly (Equation 2). As in fish schools (Berdahl et al., 2013), individuals turn toward near neighbors (Equation 1) and travel toward the peak center. This collective tracking behavior is particularly important when the resource field changes rapidly over time. As a resource peak moves, individuals at its trailing edge experience a resource value that becomes weaker through time (Figure 2A). As the resource value becomes weaker, these individuals accelerate (Equation 2), but turn toward neighbors on the peak (Equation 1) and thus travel toward the moving peak (Figure 2A). When the environment contains multiple resource peaks, evolved populations fuse spontaneously to form groups whose sizes correspond to that of the peak they are tracking (Figure 2B), even though no individual is able to assess peak size, or know whether there are multiple peaks in the environment. This behavior is consistent with recent sonar observations of foraging marine fish showing that fish form shoals that match the sizes of dynamic resource patches (Bertrand et al., 2008; Bertrand et al., 2014). Our model demonstrates that collective tracking behavior similar to that observed in real fish schools can evolve through selection on the decision rules of individuals.
Evolved populations are poised near abrupt transitions in collective state
That individuals in evolutionarily stable populations have intermediate baseline speeds and intermediate environmental sensitivities (Figure 1D–E) raises a question: what determines the evolutionarily stable values of these traits? It is tempting to conclude that these trait values are determined by the nature of the environment alone. However, the fact that the evolutionary trajectories of social and asocial populations are so different (Figure 1), suggests that the collective behaviors discussed above strongly influence the outcome of evolution. Analysis of Equations 1–4 reveals that the preferred speed parameter divides the dynamical behavior of populations into distinct collective states (Figure 3; analysis in Appendix section 5). For $\Psi <0$, individuals have a preferred speed of zero and the interindividual distances are governed by initial conditions. In this state, individuals resist acceleration due to social interactions. For small $\Psi >0$, individuals form relatively dense groups that move through the environment as collectives, either milling, swarming, or translating (D’Orsogna et al., 2006), the collective motions exhibited by real schooling fish (Tunstrøm et al., 2013). Individual speeds are relatively low and interindividual distances are short. For large $\Psi $, interindividual distances are large, and individuals move through the environment quickly. Dynamic changes among theses states are evident in Video 2. These collective states are also clearly distinguishable in Figure 3 ($0<\Psi <1.6$ and $\Psi >2.9$) and Appendix Figure 9 ($\Psi <0$), and are separated by abrupt changes in the distances between near neighbors (the inverse of local density, Figure 3) or potential energy (Appendix Figure 9). The location of transitions between states depends on the parameters of the social response rule (e.g., number of neighbors an individual pays attention to $k$; Figure 4). The transitional regimes between these states are reminiscent of the firstorder phase transitions that occur in some physical systems, for example at the transition between liquid water and water vapor. As in the liquidvapor phase transition, transitions in collective state are characterized by strong hysteresis (Figure 3). If the population begins with large $\Psi $, mean distance to neighbors remains stable for decreasing $\Psi $ and then decreases abruptly (Figure 3, Appendix Figure 9 upper curve). If $\Psi $ is then increased, mean distance to neighbors increases but follows a different functional relationship with $\Psi $ (Figure 3, lower curve). We refer to the collective states as stationkeeping ($\Psi <0$; see Appendix Figure 9), cohesive (small $\Psi $), and dispersed (large $\Psi $). The analogy between transitions in collective state in our system and first order phase transitions in physical systems can be made more precise by analyzing the formation rate of groups when $\Psi $ is in the hysteresis region. In the hysteresis region, the rate at which groups of individuals form spontaneously (and therefore nucleate a transition from the dispersed to cohesive state) depends strongly on $\Psi $; when $\Psi $ is near the upper bound of the hysteresis region, the time required for a group to form spontaneously is very long (see Appendix section 5.4). From a thermodynamic perspective, this makes the spontaneous formation of groups extremely unlikely, which explains why populations that begin in the dispersed state follow the upper branch of the hysteresis curve shown in Figure 3.
For a wide variety environmental conditions (Appendix section 2) and social parameters (Figure 4), the evolutionarily stable trait values have a notable feature: the evolved values of the baseline speed parameter, ${\psi}_{0}$, place individuals in the population slightly above the transition between cohesive and dispersed states when $S=0$ (Figure 4, upper panels, Figure 5; points in both figures show mean ${\psi}_{0}$ values of population in the ESSt), and the evolved environmental sensitivity, ${\psi}_{1}$, is large enough that locally, groups of individuals cross from the dispersed state through the cohesive and stationkeeping states in regions of the environment where the resource value is high (Figure 2A, colors indicate instantaneous value of $\Psi $ for each individual). In other words, the evolved values of ${\psi}_{0}$ and ${\psi}_{1}$ allow local subpopulations to undergo sudden changes from one collective state to another in the proximity of favorable regions of the environment. Importantly, the approximate location of the transition between cohesive and dispersed states can be predicted by directly analyzing Equations 1–4 without considering details of the environment, or the mapping between behavior and fitness (Figure 4 compare upper panels [simulation] to lower panels [analytical prediction]). While the precise evolutionarily stable values of ${\psi}_{1}$ depend on the parameters of the environment (Appendix section 2), the evolutionarily stable values of ${\psi}_{0}$ place the population near the cohesivedispersed transition in many different kinds of environments (Appendix Figure 5). As we show below, being near this transition allows groups to respond quickly to changes in the environment. Our results demonstrate, that such locations in behavioral statespace are, in fact, evolutionary attractors.
The evolutionary results presented in Figure 1 assume that individuals do not appreciably deplete the resource. We can explore an alternative scenario in which resource peaks are depleted through consumption (Appendix section 2.8). In that case, the $i$th individual consumes resources at a rate $uS\left({x}_{i}\right)$ per time step. We repeated evolutionary simulations assuming either a high or low rate of resource consumption $u$. For high consumption rate (100 individuals can deplete a peak in roughly five time steps), ${l}_{max}$ still increases so that individuals are attracted to one another through social interactions, but selection for large ${l}_{max}$ is much weaker than the case shown in Figure 1C (see Appendix Figure 7). Moreover, ${\psi}_{0}$ and ${\psi}_{1}$ increase continually through evolutionary time. This result is intuitive because when resources are depleted rapidly, the locations of neighbors convey little information about the future location of resources and transitioning from the dispersed to cohesive state may actually be maladaptive. By contrast, when individuals consume the resource at a more moderate rate (Appendix Figure 7), evolutionary trajectories parallel the trajectory shown in Figure 1C–E; there is strong selection for high ${l}_{max}$, ${\psi}_{0}$ reaches a stable value that is situated directly above the hysteresis region shown in Figure 3, and ${\psi}_{1}$ evolves to a stable value that is large enough to allow individuals to cross from dispersed to cohesive, and stationkeeping states in regions of the environment where the resource value is high.
Changes in collective state allow for rapid collective computation of the resource distribution
Why do populations of selfish individuals evolve behavioral rules that place them near the transition between collective states? Dispersed, cohesive, and stationkeeping states are each associated with a characteristic density (low, intermediate, and high, respectively; Figure 3, Appendix Figure 9). If individuals enter the cohesive and stationkeeping states where the resource level is high, the density of individuals becomes strongly correlated with the resource distribution (Figure 6A). The similarity between the distribution of individuals and the distribution of the resource can be quantified by the KullbackLeibler divergence (KL divergence), an informationtheoretic concept that measures the distance between two distributions (Figure 6A inset). Though individuals cannot sense resource gradients, they can detect gradients in the density of their neighbors (Equation 1), and can therefore move up the resource gradient.
The abrupt transitions in the density of individuals between dispersed and cohesive states (Figure 3) mean that there is a strong density gradient in regions of the environment where individuals in the dispersed state border individuals in the cohesive state (e.g., Figure 2A, 6A, Video 2). This suggests that the behavior of an individual in this region can be approximated by considering only its interactions with individuals that are on the resource peak (i.e., where density is high). Using this assumption, we derive analytically the rate at which new individuals join (or rejoin) a group on the resource peak (Appendix section 6.5). Asocial individuals arrive at a resource peak at a rate ${\kappa}_{a}$, where ${\kappa}_{a}$ is a constant (Figure 6B, blue curves and points; Equation A65). However, social individuals initially arrive at a rate that increases as more individuals reach the peak, such that the number of individuals on the peak, ${N}_{s}$, increases exponentially with time: ${N}_{s}\approx {\kappa}_{s\mathrm{,1}}+\text{exp}\left({\kappa}_{s\mathrm{,2}}t\right)$, where ${\kappa}_{s\mathrm{,1}}$ and ${\kappa}_{s\mathrm{,2}}$ are positive constants (Figure 6B, red curves and points; Equation A68–A70). Analytical calculations (Figure 6B, solid lines) agree well with results of numerical simulations (Figure 6B, points and confidence bands). The rapid accumulation shown in Figure 6 is especially important when the environment changes quickly with time; it allows groups to respond swiftly to changes in the resource field and enables the emergent resource tracking behavior described above.
The form of Equations (3–4) implies that an individual’s behavioral response combines personal information about the environment (Equation 2) with social cues (Equation 1). In fact, under a time rescaling, our model is equivalent to one in which the relative strength of social forces varies across the environment (Appendix section 4). The tradeoff between using social information and personal information is inherent in social decisionmaking (Couzin et al., 2005; Couzin, et al., 2011). This tradeoff means that individuals with large ${\psi}_{0}$ and ${\psi}_{1}$ are, by default, less responsive to their neighbors. Perturbing the values of ${\psi}_{0}$ and ${\psi}_{1}$ of individuals in populations at the ESSt show that, in populations with high mean ${\psi}_{0}$, individuals fail to form large groups and are poor at tracking resource peaks (Appendix section 2.6, Appendix Figure 6). In populations with high mean values of ${\psi}_{1}$, individuals form groups (Appendix section 2.7), but fail to exploit regions with the highest resource quality. Individuals with low values of ${\psi}_{0}$ or ${\psi}_{1}$ form groups but do not effectively track dynamic resources (Appendix section 2.7).
Discussion
Our model demonstrates that selection on the behavioral phenotypes of selfish individuals can lead to the rapid evolution of distributed sensing and collective computation. The mechanism that promotes this evolution involves the use of public information: when individuals respond to the environment by slowing down in regions of high resource quality – a behavior that is adaptive even in the absence of social interactions (Appendix Figure 2) – their positions become correlated with the locations of resources. Social individuals can exploit this public information by climbing gradients in the density of their neighbors. As in simple, gametheoretic models of social foraging (e.g., Clark and Mangel, 1984), social individuals gain a fitness advantage by using information about the environment gleaned by observing neighbors. Because of this, asocial populations are readily invaded by social mutants and collective behaviors evolve (Appendix section 2).
Evolutionarily stable populations occupy a distinctive location in behavioral state space: one in which small changes in individual behavior cause large changes in collective state (Figures 4, 5). When individuals respond to local environmental cues by accelerating or decelerating, local populations transition between the collective states shown in Figure 3 (e.g. Figure 2A). This creates the strong spatial gradient in population density (Figure 6A) and allows groups to track dynamic features in the environment rapidly. Perturbations of this evolutionarily stable state cause individuals either to weigh social information too heavily (i.e., small ${\psi}_{0}$ and/or ${\psi}_{1}$), in which case groups fail to explore effectively (Video 3, Appendix Figure 7), or to weigh personal information too heavily (i.e., large ${\psi}_{0}$ and/or ${\psi}_{1}$), in which case individuals fail to exploit the social information that enables dynamic resource tracking (Video 4, Appendix Figure 7). Because of this, mutants with phenotypes far from the evolutionarily stable state are removed from the population by natural selection. The transitions we observe in collective state bear a resemblance to phase transitions in physical systems, and our results lend credence to the hypothesis that natural selection can result in the evolution of biological systems that are poised near such bifurcation points in parameter space. Importantly, we show that these highfitness regions of parameter space can be predicted a priori from the structure of individual decision rules, even without knowledge of the environment.
Collective computation is a notion that has strongly motivated research on animal groups (Berdahl et al., 2013; Couzin, 2007; Cvikel, et al., 2015). In our model, populations perform a collective computation through their social and environmental response rules. When individuals are exposed to a heterogeneous resource environment, their responses to the environment cause a modification of the local population density; individuals aggregate in regions where the resource cue is strong. The population performs a physical computation in the formal sense (Schnitzer, 2002): physical variables – the positions and relative densities of neighbors – represent mathematical ones – spatially resolved estimates of the quality of resources in the environment. The environments considered in our study bear a strong resemblance to those encountered in dynamic coverage problems in distributed control theory (Bachmayer and Leonard, 2002), dynamic optimization problems (Passino, 2002), and Monte Carlo parameter estimation (McKay, 2003). Combining an evolutionary approach to algorithm design with collective interactions may therefore be a useful starting point for optimization schemes or control algorithms for autonomous vehicles, particularly if the structure of social interactions leads to bifurcation points in behavioral parameter space as in the model studied here.
Understanding the feedback loop between individual behavior, collective behavior of populations, and selection on individual fitness is a major challenge in evolutionary theory (Guttal and Couzin, 2010; Torney et al., 2011; Pruitt and Goodnight, 2014). Our framework closes this loop and demonstrates how distributed sensing and collective computation can evolve through natural selection on the decision rules of selfish individuals.
Materials and methods
Resource environment
Request a detailed protocolOur model of the resource environment incorporates three salient features of the resource environments that schooling fish and other social foragers encounter in nature. These features are: 1) spatial variation in resource quality, 2) temporal variation in resource quality, and 3) characteristic length scales of resource patches (Stephens et al., 2007; Bertrand et al., 2008; Bertrand et al., 2014). Accordingly, we model a twodimensional environment in which the resource is distributed as a set of $M$ resource peaks. We assume the boundary of the environment is periodic such that individuals, interindividual potentials, and resource peaks are all projected onto a torus. Each of the $M$ peaks decays like a Gaussian with increasing distance to the peak center. The value of the resource in a single peak at a location, ${x}_{i}$, is given by
where ${\lambda}_{0}$ is a constant that determines the resource value at the peak center and ${\lambda}_{1}$ is a decay length parameter, and ${x}_{s}$ is the location of the centroid of the peak of interest. The total resource value the $i$th individual experiences $S\left({x}_{i}\right)$ is the sum over all peaks in the environment. Each peak moves according to Brownian motion with drift vector $\alpha $ and standard deviation $\beta $. At each time step, each peak has a probability $1/{\tau}_{p}$ of disappearing and reappearing at a new location, chosen at random from all locations in the environment.
Appendix
1 Social interaction rules
1.1 Model of social interactions
Past individualbased models that include social interactions have often depicted social interactions by assuming that individuals monitor metric 'zones'. Individuals avoid neighbors in a small zone of avoidance, and align and move toward neighbors within larger zones of social interactions (e.g., Guttal and Couzin, 2010; Couzin et al., 2002; Chou et al., 2012). Here, we use an alternative model that depicts social interactions as forces that act to modify individuals’ accelerations. This approach is closely related to force matching methods that have been applied to data to infer the strength of pairwise social interactions among individuals. We assume that social forces depend on distance in a way that creates shortrange repulsion among individuals, strong intermediate range attraction, and weak attraction for longer ranges in agreement with results of Katz et al. (2011). We model the social forces on a focal individual, $i$, by the following equation:
where, as described in the Main Text, ${x}_{i}$ and ${v}_{i}$ are the position and velocity of the $i$th individual, respectively, $\nabla $ is the twodimensional gradient operator, the term in brackets is a social potential, ${C}_{a}$, ${C}_{r}$, ${l}_{a}$, and ${l}_{r}$ are constants, and the set ${N}_{i}$ is a set of the $k$ nearest neighbors of the $i$th individual, where a neighbor is an individual within a distance of ${l}_{max}$ of the focal individual. Appendix Figure 1 shows the effective force exerted on a focal individual by a neighbor located along the focal individual’s trajectory, either behind ($x$axis $<0$) or in front of ($x$axis $>0$) the focal individual [compare to Appendix Figure 2 of Katz et al. (2011)]. Unlike many past models of interactions among individuals, we do not assume that individuals explicitly align with one another. However, because the r.h.s of Equation A1 is proportional to the gradient of a social potential, social interactions can cause the focal individual to turn. This turning toward neighbors causes the social gradient climbing behavior described in the Main Text and discussed in detail in Appendix section 6 below.
1.2 Definition of an asocial individual
To illustrate the collective behavior and evolution of social individuals it is useful to compare social individuals to individuals that are not influenced by social attraction. We refer to such individuals as 'asocial' and define them in terms of Equation A1 by setting ${l}_{max}$ to a value that corresponded to the distance at which the gradient of the social potential for a pairwise interaction is equal to zero (blue line in Appendix Figure 1: point at which potential crosses zero). We define asocial agents in this way because the shortrange repulsion included in the interagent potential shown in Appendix Figure 1 represents collisionavoidance–a behavior that should be common to all individuals, regardless of whether they are socially attracted to one another. While this definition of an 'asocial' individual is more biologically sensible, we have also tried modeling asocial individuals by assuming that the r.h.s. of Equation A1 is equal to zero for all individuals (this assumes, for instance, that these individuals are not limited to finite local density); this approach does not qualitatively change the results presented below and in the Main Text.
2 Evolutionary dynamics
2.1 Selection algorithm
To understand the connection between the evolution of collective behaviors and selection on the performance of individuals, we implement a simple evolutionary algorithm similar to that used in Guttal and Couzin (2010). In the first generation, $N$ agents with heterogeneous values of ${\psi}_{0}$ and ${\psi}_{1}$ and ${l}_{max}$ are initiated in an environment with $M$ resource peaks. The number of agents remains constant across generations, and generations are nonoverlapping. Each generation consists of a simulation run for 5,000 or 10,000 time steps over which we calculated the mean resource value experienced by each individual in the population. At the end of each generation, $N$ individuals are selected from the population (with replacement) to reproduce themselves, yielding a total of $N$ new offspring. An individual’s probability of being selected for reproduction is proportional to its mean resource value, normalized by sum of mean resource value over all individuals in the population. Individuals that perform well are more likely to be selected to reproduce and are likely to produce more offspring than individuals that perform poorly. The selection probability of the $i$th individual ${p}_{i}$ is defined as follows:
where ${S}_{i}\left(\mathrm{}t\right)\mathrm{}$ is the instantaneous resource value of individual $i$, and angular brackets represent timeaveraging over the particular generation under consideration. If an individual is selected for reproduction, a child is produced in the next generation with ${\psi}_{0}$ equal to that of the parent, with a small mutation. The ${\psi}_{0}$ value of an offspring is equal to the ${\psi}_{0}$ value of its parent, plus a normally distributed random number $\sigma $ with mean zero and variance ${\gamma}_{m}$:
where ${\psi}_{\mathrm{0,}{i}^{\prime}}$ is the ${\psi}_{0}$ value of an offspring of individual $i$. The ${\psi}_{1}$ and ${l}_{max}$ traits of offspring were determined in the same way.
2.2 Evolution of asocial populations
In general, populations of asocial individuals evolve to have increasing ${\psi}_{0}$ and ${\psi}_{1}$ values. While fitnesses of individuals in these populations are well below fitnesses of individuals in the evolutionarily stable states discussed below (see Main Text Figure 1F), selection on asocial populations still leads to an increase in mean fitness (Appendix Figure 2). This occurs because, as evolution progresses and ${\psi}_{0}$ and ${\psi}_{1}$ values evolve, asocial individuals spend more time in regions of the environment with high resource value.
2.3 Establishment of evolutionarily stable state (ESSt)
We allow populations to evolve according to the algorithm described above. Initial values of ${\psi}_{0}$ and ${\psi}_{1}$ phenotypes are drawn at random from uniform distributions between 0 and 6. Initial values of ${l}_{max}$ are drawn with uniform probability from the interval $\left(\mathrm{0,30}\right)$. The distribution of trait values quickly stabilizes for all three phenotypic traits as shown in Figure 1C–E of Main Text. We refer to this evolved state as an evolutionarily stable state (ESSt, Guttal and Couzin, 2010). The persistent variance in the distribution of ${\psi}_{0}$, ${\psi}_{1}$, and ${l}_{max}$ are partially due to mutations in the value of these traits, which are continually introduced into the population. We therefore expect such persistent of interindividual variation in phenotype as a result of mutationselection balance.
2.4 Robustness of ESSt
To evaluate the robustness of the evolutionarily stable state (ESSt) described in the Main Text, we performed evolution under invasion by phenotypes that are both near to, and far from the ESSt. We initiated the population with trait distributions from the ESSt (selected from the final generation of simulations used to establish ESSt). Then in each generation, we selected individuals to reproduce and applied ordinary mutations as described above. However, before initiating the next generation, a single individual was chosen to serve as an invader. That individual’s phenotype was replaced by values (${\psi}_{0}^{*}$, ${\psi}_{1}^{*}$, and ${I}_{max}^{*}$). ${\psi}_{0}^{*}$ and ${\psi}_{1}^{*}$ were chosen with uniform probability from the interval $\left(\mathrm{0,40}\right)$ and ${I}_{max}^{*}$ is chosen with uniform probability from the interval ${l}_{max}\in \left(\mathrm{0,30}\right)$. Though these intervals are somewhat arbitrary, we note that ${\psi}_{0}$ and ${l}_{max}$ must ultimately be bounded above by limits on the speed that individuals can sustain, and by limits on the distance over which individuals can perceive one another, respectively. ${\psi}_{1}$ should also be bounded above because it is limited by the rate at which individuals can accelerate (decelerate) in response to changes in the measured value of an environmental cue. Thus, all three traits are bounded above due to physical constraints. Applying higher bounds on these trait values did not qualitatively change our conclusions.
Appendix Figure 3 shows a typical evolutionary progression when a population at ESSt acquires mutations (i.e., small changes in phenotype) and receives an invader in each generation. Although invaders from across the phenotype space are introduced into the population (Appendix Figure 3 blue dots across phenotype space), none of these invaders establishes for more than a few generations (Appendix Figure 3 blue dots become extinct after few generations). The ESSt is resistant to invasion by both nearby phenotypes, introduced through ordinary mutation, and phenotypes far from the ESSt, introduced through invaders. We therefore refer to the ESSt as robust.
2.5 Invasion of asocial population by social strategy
To determine whether phenotypes from the ESSt could invade a population of purely asocial individuals, we performed another set of evolutionary simulations in which we initiated populations with $N1$ asocial individuals and a single individual, chosen at random from the ESSt. Appendix Figure 4 shows evolutionary progressions from this initial state. In panel A, the full trait distribution of the population is shown. The social invader increases in frequency and sweeps the population of asocials. Replicate invasions show a very similar progression (Appendix Figure 4B). The final distribution of trait values matches the ESSt. The change in phenotypes that occur when the ESSt phenotype invades the asocial population lead to a dramatic increase in mean fitness (Appendix Figure 4C) and a decrease in the range of fitnesses of different individuals in the population.
2.6 Dependence of evolutionary outcomes on the environment
One of the conclusions drawn in the Main Text is that the trait values of the evolved population at the EESt correspond to a location in behavioral state space where the population is in the dispersed state in regions of the environment with low resource quality, and that the population transitions from the dispersed to cohesive and stationkeeping states in regions of the environment with high resource quality. We evaluated whether this conclusion holds, more generally, by evolving populations in more complex environments in which environmental properties were selected at random. We initialized trait values of populations as described in Establishment of evolutionarily stable state (ESSt) above. However, to generate the environment, we chose the number of Gaussian resource peaks at random from 1 to 50 with uniform probability. The maximum resource value of each peak and the variance of the twodimensional Gaussian peak shape were also chosen at random. Maximum resource value was chosen with uniform probability from the interval $\left(\mathrm{0,10}\right]$ and variance was chosen with uniform probability from the interval $\left(\mathrm{0,30}\right]$. Finally, the variances of all peaks in a given simulation were rescaled so that the sum of the integral of all peaks over the environment was equal to $400\pi $. We enforce this latter condition to ensure that resource peaks are small relative to the size of the environment. All other parameter values were those listed in Figure 1 of the Main Text, except that $N$ was 300.
We allowed populations to evolve for 1500 generations and recorded values of ${\psi}_{0}$, ${\psi}_{1}$, and ${l}_{max}$ that evolved. Appendix Figure 5 show mean ${\psi}_{0}$ and ${\psi}_{1}$ trait values after 1500 generations for evolutionary simulations with different environmental conditions. The gray band in Appendix Figure 5 corresponds to the region of hysteresis between cohesive and dispersed states shown in Figure 3 in the Main Text. With the exception of a small number of simulated evolutions (Appendix Figure 5, points below gray band) populations in all environments had mean trait values of ${\psi}_{0}$ in or above the hysteresis region. In all cases, the combination of ${\psi}_{1}$ and ${\psi}_{0}$ caused individuals to exhibit values of $\Psi =\mathrm{}{\psi}_{1}{\psi}_{0}S\left(\mathrm{}{x}_{i}\right)\mathrm{}$ that were less than zero in the most favorable regions of the environment. Thus, for the large majority of random environmental conditions we generated, individuals transition from $\Psi $ values that correspond to the dispersed state, through $\Psi $ values that correspond to cohesive and station keeping states in favorable regions of the environment.
2.7 Perturbation of populations around the ESSt
To further understand how the evolutionarily stable trait values lead to high individual fitness we perturbed the entire populations at ESSt by shifting either ${\psi}_{0}$ or ${\psi}_{1}$ of all individuals in the population. This resulted in a change in the mean value of these traits over the entire population. We then simulated the dynamic behaviors of the new perturbed population in a simplified environment containing two resource peaks. Initially, all individuals were located in a single group near one of the peaks (the starting peak). Appendix Figure 6 shows that, for fixed ${\psi}_{1}$, group sizes and mean fitness vary strongly as a function of the mean value of ${\psi}_{0}$ of the population (both ${\psi}_{0}$ and ${\psi}_{1}$ are taken from a population at ESSt and values of ${\psi}_{0}$ are shifted to change $\u27e8{\psi}_{0}\u27e9$ of the population). For small values of ${\psi}_{0}$, individuals track the starting peak but do not find the second peak (Appendix Figure 6A, blue and red points, respectively). As ${\psi}_{0}$ reaches approximately 2.2, individuals begin to form a group on the second resource peak (Appendix Figure 6A, red points denoting size of group nearest second peak begin to increase). Mean performance of individuals in the group nearest the second peak rapidly increases (Appendix Figure 6A, red points rapidly increase for ${\psi}_{0}>2.2$). When performance is averaged over the entire population, there is a clear maximum at ${\psi}_{0}~3.6$ (Appendix Figure 6C), the value corresponding to mean (and modal) ${\psi}_{0}$ for the evolved population in the ESSt (orange point in Appendix Figure 6C). Selection on fitnesses of individuals and optimization for maximum fitness of the entire population lead to the same value of ${\psi}_{0}$. For larger values of ${\psi}_{0}$, the average performance over all individuals begins to decline (Appendix Figure 6C) because fewer individuals aggregate near peaks. Perturbations of ${\psi}_{1}$ also lead do decreases in mean fitness at the population level (Appendix Figure 6F). For mean ${\psi}_{1}$ below that of the ESSt, individuals form small groups near peaks (Appendix Figure 6D). For mean ${\psi}_{1}$ above that of ESSt, individuals form large groups, but individuals in the groups near peaks have low fitness because they do not effectively aggregate near peak centers (Appendix Figure 6E).
2.8 Evolution with resource consumption
As described in the Main Text, social interactions confer a fitness advantage to social individuals at least in part because the positions and local densities of a given individual’s neighbors contain information about the spatial distribution of resources. However, if individuals quickly consume resources, this may break down. For example, areas in which the density of neighbors is currently high may no longer contain resources in the near future if those neighbors consume the resources. To explore how resource consumption affects evolutionary dynamics, we repeated evolutionary simulations assuming individuals consume the resource. To model resource consumption, we assume each individual consumes resources at a rate given by the product of the resource value at its position $S\left(\mathrm{}{x}_{i}\right)\mathrm{}$ and a consumption rate constant $u$. At each time step, the height of peak $j$, ${\lambda}_{\mathrm{0,}\text{\hspace{0.05em}}j}$, is reduced by the sum ${\sum}_{i=1}^{N}uS\left(\mathrm{}{x}_{i}\mathrm{,}{x}_{s}\right)\mathrm{}$, where ${x}_{s}$ is the location of the peak and $N$ is the number of individuals in the population. We assume individuals abandon a resource when ${\lambda}_{\mathrm{0,}\text{\hspace{0.05em}}j}$ falls below ${s}^{*}$. To keep the number of resource peaks constant and the total amount of resource on the landscape from being completely depleted, we allow resource peaks that reach a height of ${\lambda}_{\mathrm{0,}\text{\hspace{0.05em}}j}={s}^{*}$ to regenerate at a new location chosen at random with equal probability from all points in the environment. The new peak has a peak height equal to the starting peak height, ${\lambda}_{0}$. Mean resource value for each agent is calculated in the same manner as in the case where peaks are not depleted.
Appendix Figure 7 shows the results of replicate evolutionary simulations with high (A) and low (B) rates of resource consumption. In the case of high consumption (Appendix Figure 7A), individuals evolve to have increasing mean values of ${\psi}_{0}$ and ${\psi}_{1}$, and ${\psi}_{0}$ values are well above the hysteresis regime between collective states. While values of ${l}_{max}$ still enable individuals to be attracted to one another at intermediate to large distances, the variation in ${l}_{max}$ values among replicate simulations suggests that there is not strong selection for large ${l}_{max}$. When individuals consume resources at a lower rate (Appendix Figure 7B), results parallel those shown in Figure 1 of the Main Text; populations evolve mean values of ${\psi}_{0}$ that are directly above the hysteresis region, ${\psi}_{1}$ approaches a stable value, and ${l}_{max}$ approaches the maximum allowable value of 30.
3 The cohesive state is characterized by a fixed, finite density
Agents obeying the equations described in the Main text exhibit several distinct collective states. One such state, which we call the cohesive state, is characterized by dense groups of agents occupying a fixed fraction of the environment. One of the salient properties of these groups is that they eventually reach a fixed density that becomes independent of group size. Using a small number of simple assumptions about the behavior of agents within a cohesive group, we are able to predict the density of agents directly from the model parameters. The motivation of our calculation comes from the structure of the equations, which include social potential terms and velocitydependent selfpropulsion terms. The social force on an agent is given by Equation A1. We will rewrite the social potential in Equation A1 (i.e., the term in brackets) as
The effect of the potential term in the equations is to exert a force on the entire system towards configurations where the potential energy is lower. The propulsive forces are nonconservative, causing phasespace volumes to contract, allowing the system to approach a potential energy minimum. We model the cohesive state as $N$ agents occupying a circular region of radius $l$. Further, we assume that the probability distribution of agents within this circular region is uniform, so that agent density is given by:
The density $\rho $ lets us define an interaction radius ${l}_{I}$ which is expected to contain $k$ individuals:
This expression is valid when $N>\mathrm{}k$, which is the case we are interested in. When $N<\mathrm{}k$ the interaction radius is simply the group radius, $l=\mathrm{}{l}_{I}$, and each agent interacts with every other agent. We calculate the expected potential by integrating over a circle of radius ${l}_{I}$:
This integral evaluates to the following expression:
It is illustrative to write this expression after the substitution $l=\mathrm{}\sqrt{\frac{N}{\pi \rho}}$:
The density that minimizes ${\Phi}_{i}$ will be the density of agents in the cohesive state. When written this way, it is clear that $N$ will not influence the location of the minimum of ${\Phi}_{i}$, as $N$ only appears in the expression as a constant multiplier. Thus, when $N>\mathrm{}k$ we expect cohesive groups to have a constant density, so that the radius of a group grows like $\sqrt{N}$. These predictions match the results of our simulation quite closely, which one can see from comparisons between Appendix Figure 8 to the lower branch of the hysteresis plot in Appendix Figure 9.
That density is necessarily constant with increasing $N$ is a hallmark of topological interaction laws which are repulsive at short range. When there is no restriction on the number of interaction neighbors, an interaction law of the type that we use can give rise to catastrophic behavior, where the group density increases without bound with increasing $N$ (D’Orsogna et al., 2006). One feature of the topological interaction is that it allows for biological realism for agent parameter values that would otherwise lead to catastrophic behavior.
4 Relationship between ${\psi}_{0}$ and the relative strength of social forces
In order to better understand how changing parameters affect our model, we ignore stochasticity and consider the equations for the acceleration and velocity in a homogeneous environment without resources (so that the background velocity is constant). First we define:
The equations become:
Let ${l}_{0}$ be a characteristic length scale in this problem, let ${v}_{0}$ be a characteristic velocity scale, and let ${t}_{0}=\mathrm{}\frac{{l}_{0}}{{v}_{0}}$ be a characteristic time scale. Then, let ${C}_{a}$, the attraction coefficient, be the scale of the potential. We nondimensionalize our equations by rewriting them in terms of the dimensionless variables:
The resulting dimensionless equations are:
The nondimensional number $\frac{{l}_{0}^{2}{C}_{a}}{{v}_{0}^{2}}=\mathrm{}\frac{\eta {l}_{0}^{2}{C}_{a}}{{\psi}_{0}}$ measures the relative strength of the social potential. When ${\psi}_{0}$ becomes large, social forces become negligible. The reason for this effect is that agents begin moving too quickly for the social forces to have any appreciable effect on their trajectories. Therefore, for constant ${l}_{0}$ and ${C}_{a}$, an alternative interpretation of ${\psi}_{0}$ is as a term that dictates the relative strengths of autonomous versus social forces.
5 Continuum description predicts cohesive and dispersed state
Figure 3 in the Main Text illustrates that populations exhibit distinct regimes, which we refer to as collective states, as a function of the preferred speed parameter. Two states are evident in Figure 3 in the Main Text: a state with short interindividual distances for small $\Psi $ and a state with large interindividual distances for large $\Psi $. A third state is evident if the mean potential energy of all individuals in the population is plotted as a function of $\Psi $ in a uniform environment, where potential energy is calculated from Equations 3 and 4 in the Main Text. For $\Psi <0$ there is a distinct drop in potential energy for decreasing $\Psi $ (Appendix Figure 9). We refer to the state that occurs for $\Psi <0$ as stationkeeping in the Main Text.
In order to better understand the behavior of agents in the context of our model, we have developed a continuum equation for the time evolution of agent density. In the context of a homogeneous environment this description can be used to predict the points in parameter space at which the uniform, purely solitary state becomes unstable, and to demonstrate that heterogeneous states cannot be stable at high enough background velocity. Although the continuum description is only an approximation, it is able to qualitatively predict many of the features of our multiagent simulations, which makes the mechanisms responsible for this behavior mathematically more explicit. In order to derive continuum equations, we begin with a Liouville equation for the probability density of all the $N$ particles within the full phase space, and derive a hierarchy of equations by taking moments (Flierl et al., 1999; Born and Green, 1946). This hierarchy can be closed by assuming that stochastic forces are sufficiently strong to ensure independence of the individual agents. For the analysis presented here, we assume that the agents travel at a constant velocity ${v}_{0}=\mathrm{}\sqrt{\Psi}$ (using the angular variable $\theta $ to describe the direction of the velocity), and that there is noise in the angular velocity driven by a Wiener process with variance $\sqrt{2\epsilon}$ per unit time. The assumption of a constant velocity implies that we have taken a limit where $\eta $, ${\psi}_{0}$, and ${\psi}_{1}$ go to infinity, though their ratios remain constant. We will let ${\psi}_{0}$ and ${\psi}_{1}$ stand for the limiting ratios of the original model. We will denote the space of positions by $V$, which will be a 2torus with length ${L}_{D}$. The assumption of a constant velocity and of angular noise lead only to small quantitative changes in agent behavior, and they make it possible to analyze the resulting equations.
Therefore we begin with the following set of stochastic differential equations:
Here ${F}_{i}$ is the social force on agent $i$, and ${\Gamma}_{i}$ is the angular direction of the social force on agent $i$. We assume that this force is produced through a topological interaction of the following form:
Here ${N}_{i}$ is the set of $k$ closest neighbors to agent ${x}_{i}$, and $w$ is an interaction kernel.
From these equations, we are able to write a Liouville equation by introducing a probability density on phase space, $P\left(\right\{\mathrm{}x\},\{\theta \left\}\right)\mathrm{}$, where $\left\{\mathrm{}x\right\}\mathrm{}$ is the set of $N$ agent positions and $\left\{\mathrm{}\theta \right\}\mathrm{}$ is the set of $N$ agent directions. The value of $P$ at a given set of positions and directions is the probability that the each of the agents have the specified positions and directions. The Liouville equation is:
In order to simplify this equation, we assume that the probability density function $P$ can be factorized into the product of $N$ identical single particle probability density functions:
This assumption is equivalent to assuming statistical independence of the positions and directions of the agents, a condition which could be reached either through large stochasticity $\epsilon $ or ergodic single particle trajectories. The assumption allows us to derive a closed equation for the single particle probability density function $p$, in a similar fashion to a closure of the usual BBGKY hierarchy in kinetic theory.
Then, we are able to write the following equation for the single agent probability density function $p\left(\mathrm{}x\mathrm{,}\theta \mathrm{,}t\right)\mathrm{}$ (where we have replaced a binomial distribution with a Poisson distribution):
Here, the expression $\lambda (\mathrm{}x\mathrm{,}x\text{'})\mathrm{}$ is the expected number of agents within a distance $\mathrm{}x\text{'}\mathrm{}$ from the point $x$, the function $w\left(\right\mathrm{}x\left\right)\mathrm{}$ is the social potential between two particles, and the expression $\Gamma (\mathrm{}x\mathrm{,}x\text{'})\mathrm{}$ is the angle of force from an agent located at $x\text{'}\mathrm{}$ to an agent located at $x$:
The presence of the terms involving $\lambda $ are a consequence of the topological interaction law between the agents. This equation is most accurate in the limit where ${\rho}_{c}{l}_{c}^{2}\gg 1$, where ${\rho}_{c}$ is a characteristic density and ${l}_{c}$ is the characteristic length scale of the interaction. In the examples which we considered in our study, this ratio is typically only slightly larger than 1 (see Peshkov et al., 2012; Chou et al., 2012) for derivations of continuum descriptions of topological interactions in a more collisional regime). Despite that, we find both quantitative and qualitative agreement between the continuum description and the agent based model. This kinetic description can be converted into a hierarchy of fluid equations by taking moments with respect to the angular direction variable $\theta $.
We introduce the phase space particle density $f$:
The Fourier series for $f$ gives the important macroscopic variables, for instance:
Here $u$ is the mean velocity of agents at each point in space. We will take moments of the kinetic equation through Fourier series:
The time evolution equation for the $m$th Fourier coefficient is:
Here, we have simplified this expression by introducing two new functionals of the density, $F\left(\mathrm{}x\mathrm{,}\rho \right)\mathrm{}$ and $\Theta \left(\mathrm{}x\mathrm{,}\rho \right)\mathrm{}$, which represent the social force exerted at the point $x$ due to the density $\rho $ and the direction of that social force. Explicitly these are given by:
The evolution of the $m$th moment depends on the value of the $m+1$st moment, so that we have an infinite hierarchy of equations. Moments with high values of $\left\mathrm{}m\right\mathrm{}$ experience strong damping, and we can use this to justify discarding all moments with $\left\mathrm{}m\right\mathrm{}$ above a given threshold. In the following treatment we will set to zero all Fourier coefficients with $\left\mathrm{}m\right>1$, which is the simplest truncation of the hierarchy that leads to nontrivial equations.
The right hand side of the momentum equation leads to rapid equilibration, and we can eliminate the time derivative in this equation. This allows us to find an expression for $\rho u$ in terms of $\rho $ only:
We can use this to write a single closed equation for $\rho $. In order to facilitate the analysis of this equation, we make one further approximation, replacing the sum over Poisson factors with a Heaviside function that is equal to $1$ if the expected number of agents between inside a ball of size $x\text{'}$ is less than $k$, and $0$ otherwise. This captures the dominant qualitative feature of the topological interaction in a simple way: the effective interaction radius is a function of the local density. This approximation is quantitatively consistent with the assumptions of the previous section and the results of our simulations. The resulting equation is:
The advectiondiffusion equation will be used to understand the behavior of our multiagent simulations. From its form one can see the effects of ${v}_{0}$: large ${v}_{0}$ enhances the diffusivity and reduces the effects of the potential.
5.1 Formation of the cohesive state and stability of the dispersed state
Any constant density function ${\rho}_{0}$ is an equilibrium solution of Equation A36. A crucial property of our multiagent model is that depending on the background velocity, the agents have the ability to spontaneously form a dense state which we call the cohesive state. In order for this to be possible, the uniform state must be unstable. One advantage of the continuum description is that it allows us to investigate such questions within a much simpler framework than in the original agent based model. In order to do so, we select a uniform state with value ${\rho}_{0}$ and linearize around it, neglecting terms second order or higher in the deviation away from ${\rho}_{0}$:
This is translationally invariant, and we have periodic boundary conditions, so we consider the Fourier coefficients of ${\rho}_{1}$:
The term $G\left(\right\mathrm{}j\left\right)\mathrm{}$ can be calculated by application of the convolution theorem and integration by parts, using the fact that the integrals are radially symmetric (which is true as long as ${L}_{0}<\mathrm{}{L}_{D}$):
Here ${J}_{1}$ stand for the corresponding Bessel functions of the first kind. Linear stability is determined by the sign of the coefficient on ${\stackrel{~}{\rho}}_{\mathrm{1,}j}$ on the right hand side of the following expression:
We can use our formula to determine the stability or instability of an arbitrary homogeneous equilibrium solution. We have plotted an example of this in Appendix Figure 10. A number of general features emerge from these diagrams. Increases in the background velocity ${v}_{0}$ always promotes stability of the dispersed state. The agents make use of this feature to enable themselves to transition from the dispersed state to the cohesive state in regions where $\Psi $ crosses below the stability threshold. Increases in the number of interaction neighbors $k$, the decay length of the attractive interaction ${l}_{a}$, and the strength of the attractive interaction ${C}_{a}$ all promote instability of the dispersed state, though at large $k$ further increases in $k$ have little effect. The background density of agents has a more complicated effect on the stability of the dispersed state, when ${\rho}_{0}$ is low the social forces are very weak because the distance between agents is large, so that a very small ${v}_{0}$ is required for formation of the dispersed state. When ${\rho}_{0}$ gets too large, the repulsive part of the interaction becomes more important and stability of the dispersed state is promoted.
5.2 Nonlinear stability of the dispersed state for high ${v}_{0}$
Equation A36 combines a diffusive term with effective diffusion coefficient $\frac{{v}_{0}^{2}}{2\epsilon}$, and a term due to the social forces, which is proportional to the magnitude of the social forces and $\frac{1}{2\epsilon}$. On the basis of this, we expect that as ${v}_{0}$ increases the diffusive terms become more important relative to the social forces. In the linear stability analysis, this manifested itself through instability of the homogeneous base state when ${v}_{0}$ was decreased below a threshold. When ${v}_{0}$ is large enough, we are able to prove using energy inequalities (Doering and Gibbon, 1995) that the homogeneous equilibrium state is a global attractor. The implication of this is that above a certain threshold the cohesive state can no longer exist, and all the agents enter the dispersed state. Combined with the results of the previous subsection, this provides an analytical demonstration of the hysteresis that we observe in our multiagent model.
In order to establish these results, we rewrite the dynamical equation for $\rho $ in terms of the deviation from the mean density $\overline{\rho}=\mathrm{}{\displaystyle {\int}_{V}}\frac{\rho}{{L}_{D}^{2}}{d}^{2}x$. We define ${\rho}_{1}=\mathrm{}\rho \overline{\rho}$. Then the equation for ${\rho}_{1}$ is:
We multiply by ${\rho}_{1}$ on both sides of the equation and integrate:
Here the expression $\parallel f{\parallel}_{p}={\left({\int}_{V}{d}^{2}x\leftf\right(x){}^{p}\right)}^{1/p}$ is the standard ${L}^{p}$ norm on the space $V$. The first term on the right hand side can be bounded through use of the Poincare inequality for a mean zero function on the torus:
The second term on the right hand side can be simplified by performing integration by parts in order to transfer the gradient operators onto the ${\rho}_{1}$ terms:
Here we have made use of Young’s Inequality, which states that $\left\rightf*g{}_{r}\le \leftf\right{}_{p}\xb7\left\rightg{}_{q}$ when $\frac{1}{r}=\mathrm{}\frac{1}{p}+\frac{1}{q}1$.
The third term on the right hand side is the most difficult to deal with because it contains three powers of ${\rho}_{1}$. To bound this term we make use of the fact that the density $\rho $ is always positive, which implies that ${\rho}_{1}>\mathrm{}\overline{\rho}$. Then, because ${\rho}_{1}$ has zero mean, we can bound $\left\right{\rho}_{1}{}_{1}$:
Using this expression (and Young’s Inequality), we find that:
Using these bounds, we can write a differential inequality for $\frac{d\left\right{\rho}_{1}{}_{2}^{2}}{dt}$:
If v_{0} satisfies the following inequality:
then the coefficient of $\left\right\nabla {\rho}_{1}{}_{2}^{2}$ in the differential inequality is negative, and we can apply the Poincare inequality to find:
This inequality allows us to use Gronwall’s lemma (Doering and Gibbon, 1995) to prove $\left\right{\rho}_{1}{}_{2}$ that converges to 0 as a function of time, which implies that the homogeneous state is globally attractive for sufficiently large v_{0}.
5.3 Conclusion from continuum model
The simulations of the agent based model indicated that our agents possessed two properties: for small ${v}_{0}$ a cohesive state forms spontaneously, for large ${v}_{0}$ only dispersed states are possible, and for moderate values of ${v}_{0}$ both cohesive and dispersed states are possible. We were able to create a continuum model that demonstrates the mechanisms behind these numerical observations. We showed that for small ${v}_{0}$, homogeneous background states are linearly unstable to the formation of clumped states. For larger ${v}_{0}$, the homogeneous background states are linearly stable. Further, we showed that for sufficiently large ${v}_{0}$, the homogeneous state is globally attractive, so that clumped states are not possible.
5.4 Some additional properties of the transition between collective states: nucleation rates and hysteresis
In the theory of first order phase transitions, hysteresis often arises because there is a free energy penalty for small droplets of the stable phase. This leads to extremely low probabilities of critical droplet formation near the transition temperature. In order to test whether this effect is responsible for the hysteresis in our model equations, we performed long time numerical simulations using values of $\Psi $ within the hysteresis region, allowing us to estimate the nucleation rate of the cohesive phase. We performed replicate simulations with ${10}^{3}$ agents, restarting the simulations each time the agents were able to form the cohesive state. The results of this simulation are shown in Appendix Figures 11 and 12, illustrating the superexponential growth in the mean nucleation time as $\Psi $ increases. This growth in nucleation times corresponds to an increase in the minimum radius of a stable group. When $\Psi $ increases above 1.7, the expected time for nucleating the cohesive state becomes extremely large, leading to strong hysteresis. We also computed $\frac{1}{\sqrt{log\overline{T}}}$ to illustrate the approximate scaling of the nucleation time (Appendix Figure 12). Because we use a topological interaction, we do not necessarily expect this scaling to hold for much larger values of $\Psi $, as groups with $N<25$ will have an increasing, rather than constant, potential energy per particle.
6 Social gradient climbing and aggregation on a resource peak
In this section we derive a model of collective exploration and exploitation that allows us to understand how the ability of agents to find the resource peak changes when model parameters like ${\psi}_{0}$ or sociality are varied. We model agents as being either in a cohesive state near a resource peak or in a dispersed state. Using the model parameters and some simple assumptions about the dynamics, we calculate the fraction of particles approaching the resource peak that are able to enter the cohesive state. Using this model, we quantitatively estimate the rate at which agents are able to find the peak and the advantage of social agents versus asocial agents. We begin by stating a number of assumptions, each of which arises from some feature of our multiagent model, that help make our theoretical analysis tractable.
6.1 Assumptions
The environmental response function is $\Psi \left(x\right)={\psi}_{0}{\psi}_{1}A{e}^{{r}^{2}/{\sigma}^{2}}$, indicating a single peak located at the origin and a preferred background velocity of ${v}_{0}=\sqrt{{\psi}_{0}}$ in regions far from the resource peak.
Agents travel at the speed dictated by the environmental response function Ψ, so that $v\left(x\right)=v\left(r\right)=\sqrt{\Psi}$. Only the direction of the velocity is allowed to vary.
Particles exist in one of two states, the cohesive state or the dispersed state. Particles in the cohesive state are close to the resource peak. Particles in the dispersed state have a uniform probability distribution in space and in direction. Particles in the cohesive state collectively produce a potential ${\Phi}_{\mathrm{peak}}\left(r\right)=\mathrm{min}(k,N){C}_{a}{e}^{\frac{r}{{l}_{a}}}$.
Agents in the dispersed state interact only with particles in the cohesive state, and this interaction is cutoff for distances $r{l}_{\mathrm{max}}={r}_{M}$. The potential force is projected normal to the velocity of the agents so that it can only effect the agent directions.
An agent enters the cohesive state if it has a trajectory that reaches the radius of zero velocity, ${r}_{t}={\sigma}_{}\sqrt{\mathrm{log}\left(\frac{{}_{A{\psi}_{1}}}{{}_{{\psi}_{0}}}\right)}$, which is assumed to mark the transition between the cohesive and dispersed states.
6.2 Critical angle for capture by the peak
Consider a particle reaching ${r}_{M}$, the radius where it begins to feel the influence of the agents on the environmental resource peak, as is depicted in Appendix Figure 13. If the angle of the agent’s trajectory is sufficiently directed towards the resource peak, the agent will reach the peak, and if the angle is directed sufficiently away, the agent will not reach the peak. There is a critical angle ${\Delta}_{i}$ at the boundary between these two scenarios. The size of ${\Delta}_{i}$ will determine the fraction of agents captured by the resource peak after crossing $r=\mathrm{}{r}_{M}$, and it will also determine the flux of agents onto the resource peak.
To derive an expression for ${\Delta}_{i}$, we write equations for an agent traveling with a velocity of magnitude $v\left(\mathrm{}r\right)=\mathrm{}\sqrt{{\psi}_{0}{\psi}_{1}A{e}^{{r}^{2}/\mathrm{}{\sigma}^{2}}}$, in direction $\omega $, experiencing the potential force $\nabla {\Phi}_{\text{p}\mathrm{eak}}\left(r\right)\mathrm{}$:
Our question is the following: given initial radius $r=\mathrm{}{r}_{M}$, initial angle ${\theta}_{0}$, and initial direction ${\omega}_{0}$, does the agent reach the zero velocity radius? The equations of motion are:
We define the angle $\Delta =\mathrm{}\pi \omega +\theta $, which is the angle between the velocity vector and the vector directed from the position of the agent to the origin. Then we can rewrite our equations in terms of the variables $r$ and $\Delta $ alone, leading to the following planar system:
The system of equations described here has the following properties:
If $\Delta =0$, then $\frac{d\Delta}{dt}=0$, so that $\Delta =0$ for all further times. Similarly, $\Delta =0$ implies $\frac{dr}{dt}<0$, so that any agent with $\Delta =0$ will reach the zero velocity radius.
If $\frac{\pi}{2}<\Delta <\mathrm{}\frac{\pi}{2}$, then $\frac{dr}{dt}<0$ and the agent will move closer to the peak.
If both $\left(\frac{{\Phi}_{peak}^{\text{'}}\left(\mathrm{}r\right)\mathrm{}}{v\left(\mathrm{}r\right)\mathrm{}}\frac{v\left(\mathrm{}r\right)\mathrm{}}{r}\right)>0$ and $\frac{\pi}{2}<\Delta <\mathrm{}\frac{\pi}{2}$, then $\frac{d\left\Delta \right\mathrm{}}{dt}<0$, and $\Delta $ becomes closer to $0$.
We make one additional assumption, which has been true in most practical cases, that allows us to make progress with the analysis.
Assumption: The function ${\Phi}_{\mathrm{peak}}^{\text{'}}\left(r\right)\frac{v(r{)}^{2}}{r}$ has exactly one sign change on the interval $({r}_{0},\infty )$. The location of the sign change occurs when ${\Phi}_{\mathrm{peak}}^{\text{'}}\left(r\right)=\frac{v(r{)}^{2}}{r}$, a point which we denote by r*. At r = r*, there is a balance between centrifugal and potential forces. We have the following inequalities:
This assumption allows us to divide the $(r,\u2206)$ plane into four regions:
 $\frac{\pi}{2}<\u2206<\frac{\pi}{2}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\text{r < r*}$
 $\u2206<\frac{\pi}{2}\phantom{\rule{1em}{0ex}}\text{or}\phantom{\rule{1em}{0ex}}\frac{\pi}{2}\text{<\u2206}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\text{rr*}$
 $\frac{\pi}{2}<\u2206<\frac{\pi}{2}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\phantom{\rule{1em}{0ex}}\text{rr*}$
 $\u2206<\frac{\pi}{2}\phantom{\rule{1em}{0ex}}\text{or}\phantom{\rule{1em}{0ex}}\frac{\pi}{2}\text{<\u2206}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\text{rr*}$
In region 1, the agents move towards the peak and the potential force is stronger than the centrifugal force. In region 2, the agents move away from the peak, and the potential force is weaker than the centrifugal force. In region 3, agents move towards the peak but the potential force is weaker than the centrifugal force. In region 4, the agents move away from the peak but the potential is stronger than the centrifugal force. We conclude that:
Any trajectory that enters region 1 will reach the zero velocity radius.
Any trajectory that enters region 2 will escape to $\infty $.
Consider again the hypothetical agent in region 3 and at radius r = r_{M}. The agent will eventually either reach region 1, region 2, or the boundary points between the two regions (which are unstable equilibria).
The points $\left(r*,\pm \frac{\pi}{2}\right)$ are unstable equilibrium points, each corresponding to a periodic orbit around the resource peak.
There are two values of Δ such that the solution of the initial value problem with initial condition (r_{M}, Δ) reach these equilibria. We call these angles $\pm {\u2206}_{i}$, where we define ${\u2206}_{i}>0$. Any trajectory with initial value $({r}_{M},\u2206)$, with $\u2206<{\u2206}_{i}$ will enter region I and be captured by the resource peak, and any trajectory with initial value satisfying $\u2206{\u2206}_{i}$ will enter region II and escape the resource peak.
The angle ${\Delta}_{i}$ is the critical angle that we seek.
6.3 Solving the reduced system
Instead of considering the time dependent differential equation, we search for an equation that describes the shape of a trajectory, that is, we assume $\Delta $ is a single valued function of $r$, and use the original equation in combination with the chain rule to write a differential equation for $\Delta \left(\mathrm{}r\right)\mathrm{}$. This method is valid in regions where $\Delta $ is actually a singlevalued function of $r$, and for this to be the case integration must be restricted to regions where ${\Phi}_{\mathrm{peak}}^{\text{'}}\left(\mathrm{}r\right)\mathrm{}\frac{{v}^{2}\left(\mathrm{}r\right)\mathrm{}}{r}$ has only one sign. The resulting equation will be valid only while a trajectory is in region $3$, and the coefficients of this equation will blow up at the border of region $3$.
The quotient of Equation A56 and A57 is the desired equation:
Equation A59 can be solved by integrating along a trajectory beginning at $\left(\mathrm{}{r}_{M}\mathrm{,}{\Delta}_{i}\right)\mathrm{}$, and ending at $\left({r}^{*}\mathrm{,}\frac{\pi}{2}\right)\mathrm{}$, leading to:
To simplify the resulting expressions, let $F\left(\mathrm{}r\right)\mathrm{}$ be the antiderivative of $\frac{{\Phi}_{\mathrm{peak}}^{\text{'}}\left(\mathrm{}r\right)\mathrm{}}{v{\left(\mathrm{}r\right)\mathrm{}}^{2}}$.
This leads to an integral equation:
This equation can be solved for the critical angle ${\Delta}_{i}$:
The critical angle ${\Delta}_{i}$ is a function of the parameters defining the agent behavior, such as ${\psi}_{0}$ and ${\psi}_{1}$, and the parameters defining the clump, such as the peak occupancy $N$. When ${\psi}_{0}$ is very small, trajectories spend much more time under the influence of the potential, and consequently it is much more likely that they are captured by the peak. Thus, for small ${\psi}_{0}$, ${\Delta}_{i}=\mathrm{}\pi /2$. As $N$ increases, the potential becomes stronger, and the values of ${\Delta}_{i}$ are increased for all ${\psi}_{0}$. When ${\psi}_{0}$ is too large ${\Delta}_{i}$ goes to $0$ and agents cannot find the peak. Appendix Figure 15 contains a plot demonstrating the aforementioned properties of ${\Delta}_{i}$.
Appendix Figure 16 contains a plot of the direction field Equation A59 and plots of trajectories that reach $\left(\mathrm{}{r}^{*}\mathrm{,}\pm \frac{\pi}{2}\right)\mathrm{}$, demonstrating the trapping of trajectories that enter region $I$, and providing numerical confirmation of our formula for the critical angle ${\Delta}_{i}$.
6.4 Equation for peak exploration
Using Equation A63 for ${\Delta}_{i}$, we can write an equation for the rate at which the number of agents occupying a resource peak increases.
We assume that there is a population of $P$ agents moving in a torus of width $L$, and that $N$ of these agents occupy a resource peak at the origin.
The spatial density of agents away from the peak is homogeneous and equal to $\frac{PN}{{L}^{2}}$.
The velocity of agents located at $r>\mathrm{}{r}_{M}$ has magnitude ${v}_{0}=\mathrm{}\sqrt{{\psi}_{0}}$ and is uniformly distributed in direction.
When an agent reaches $r=\mathrm{}{r}_{M}$, if it has ${\Delta}_{i}<\Delta <{\Delta}_{i}$, the agent will be captured by the peak. Otherwise it will escape.
This allows us to calculate the rate of capture of agents on the peak. The flux of agents to the radius ${r}_{M}$ and the angle $\Delta $ is equal to $\frac{\rho {v}_{0}}{2\pi}$. The flux of agents to a point on the circle with ${\Delta}_{i}<\Delta <{\Delta}_{i}$ is equal to $\frac{\rho {v}_{0}{\Delta}_{i}}{\pi}$. Then integrating over the circle with radius ${r}_{M}$ gives us the total flux to the peak, or the rate of change of the peak occupancy $N$:
6.5 Comparison of social versus asocial exploration
We can perform a simple calculation to demonstrate how sociality enhances the rate at which agents occupy a resource peak. In the context of this model, the difference between social and asocial agents is that the flux of asocial agents to a peak is not enhanced by the presence of agents on a peak. Thus the rate at which the number of asocial agents occupying a peak increases is linear in time. Indeed, if we assume that the total population is large in comparison with the number of agents on the peak, then we can approximate the arrival of asocial agents onto the peak with the following differential equation:
If the peak is unoccupied at time $t=0$, this equation has solution:
In contrast, the flux of social agents to a peak is enhanced by the presence of other agents on the peak. A similar approximation leads to the equations:
Appendix Figure 17 contains a plot of the function ${\Delta}_{i}\left(\mathrm{}N\right)\mathrm{}$ versus $N$. This plot motivates approximating ${\Delta}_{i}\left(\mathrm{}N\right)\mathrm{}$ as a piecewise linear function, linearly increasing from ${\Delta}_{i}\left(0\right)$ for small $N$ until the value ${N}_{M}$ at which ${\Delta}_{i}=\mathrm{}\frac{\pi}{2}$, at which point the flux becomes a constant function equal to $\pi \rho {r}_{M}{v}_{0}$. In the initial phase, we approximate the differential equation with:
When the peak is unoccupied at $t=0$, this has solution:
This solution is good up until $N=\mathrm{}{N}_{M}$, which happens at:
When $t>\mathrm{}{t}_{M}$, the solution is:
Appendix Figure 18 compares the function $N\left(\mathrm{}t\right)\mathrm{}$ for social and asocial agents.
7 Numerical methods
We used the CVODE subroutine of the SUNDIALS package to numerically solve the agentbased model (Hindmarsh et al., 2005). The resulting system of ODEs is stiff, so we utilized the variable order backwarddifferentiation methods provide by SUNDIALS. We found these implicit methods to be much more efficient than explicit methods for the particular problem that we considered. We also made use of the armadillo linear algebra library (Sanderson, 2010), the MATLAB statistics and machine learning toolbox for nearest neighbor searches, and the mex file libraries to interface all of these different tools (MATLAB, 2015).
References

1
Proceedings of 41st IEEE Conf. on Decision and Control112–117, Vehicle networks for gradient descent in a sampled environment, Proceedings of 41st IEEE Conf. on Decision and Control.
 2
 3
 4

5
A general kinetic theory of liquids. i. the molecular distribution functionsProceedings of the Royal Society of London 188:10–18.https://doi.org/10.1098/rspa.1946.0093

6
From disorder to order in marching locustsScience 312:1402–1406.https://doi.org/10.1126/science.1125142
 7

8
Foraging and flocking strategies: information in an uncertain environmentThe American Naturalist 123:626–641.https://doi.org/10.1086/284228
 9
 10
 11

12
Collective memory and spatial sorting in animal groupsJournal of Theoretical Biology 218:1–11.https://doi.org/10.1006/jtbi.2002.3065

13
Selforganization and collective behavior in vertebratesAdvances in the Study of Behavior 32:1–75.https://doi.org/10.1016/S00653454(03)010015
 14

15
Applied Analysis of the NavierStokes EquationsCambridge: Cambridge University Press.https://doi.org/10.1017/CBO9780511608803

16
Selfpropelled particles with softcore interactions: patterns, stability, and collapsePhysical Review Letters 96:104302.https://doi.org/10.1103/PhysRevLett.96.104302

17
From individuals to aggregations: the interplay between behavior and physicsJournal of Theoretical Biology 196:397–454.https://doi.org/10.1006/jtbi.1998.0842

18
Towards a general formalization of encounter rates in ecologyTheoretical Ecology 6:189–202.https://doi.org/10.1007/s1208001201704

19
Social interactions, information use, and the evolution of collective migrationProceedings of the National Academy of Sciences of the United States of America 107:16172–16177.https://doi.org/10.1073/pnas.1006874107
 20

21
Sensing and decisionmaking in random searchProceedings of the National Academy of Sciences of the United States of America 109:12070–12074.https://doi.org/10.1073/pnas.1202686109

22
Inferring the rules of interaction of shoaling fishProceedings of the National Academy of Sciences of the United States of America 108:18726–18731.https://doi.org/10.1073/pnas.1109355108

23
SUNDIALSACM Transactions on Mathematical Software 31:363–396.https://doi.org/10.1145/1089014.1089020

24
Inferring the structure and dynamics of interactions in schooling fishProceedings of the National Academy of Sciences of the United States of America 108:18720–18725.https://doi.org/10.1073/pnas.1107583108

25
Coarsegrained analysis of stochasticityinduced switching between collective motion statesProceedings of the National Academy of Sciences of the United States of America 104:5931–5935.https://doi.org/10.1073/pnas.0608270104
 26

27
Information Theory, Inference, and Learning AlgorithmsCambridge, UK: Cambridge University Press.

28
Biomimicry of bacterial foraging for distributed optimization and controlIEEE Control Systems Magazine 22:52–67.https://doi.org/10.1109/MCS.2002.1004010

29
Continuous theory of active matter systems with metricfree interactionsPhysical Review Letters 109:098101.https://doi.org/10.1103/PhysRevLett.109.098101

30
Functions of shoaling behaviour in teleostsIn: TJ Pitcher, editors. Behaviour of Teleost Fishes. Dordrecht: Springer Netherlands. pp. 363–439.https://doi.org/10.1007/9789401115780_12
 31

32
Revealing the hidden networks of interaction in mobile animal groups allows prediction of complex behavioral contagionProceedings of the National Academy of Sciences of the United States of America 112:4690–4695.https://doi.org/10.1073/pnas.1420068112

33
Armadillo: an open source c++ linear algebra library for fast prototyping and computationally intensive experiments technical reportNICTA.

34
Theory of continuum random walks and application to chemotaxisPhysical Review E 48:2553–2568.https://doi.org/10.1103/PhysRevE.48.2553
 35

36
Evolution and the Theory of GamesCambridge: Cambridge University Press.https://doi.org/10.1017/CBO9780511806292

37
Foraging: Behavior and EcologyUSA: University of Chicago Press.https://doi.org/10.7208/chicago/9780226772653.001.0001

38
Visual sensory networks and effective information transfer in animal groupsCurrent Biology : CB 23:R709–R711.https://doi.org/10.1016/j.cub.2013.07.059

39
Contextdependent interaction leads to emergent search behavior in social aggregatesProceedings of the National Academy of Sciences of the United States of America 106:22055–22060.https://doi.org/10.1073/pnas.0907929106

40
Signalling and the evolution of cooperative foraging in dynamic environmentsPLoS Computational Biology 7:e1002194.https://doi.org/10.1371/journal.pcbi.1002194

41
Collective states, multistability and transitional behavior in schooling fishPLoS Computational Biology 9:e1002915.https://doi.org/10.1371/journal.pcbi.1002915
Decision letter

Michael DoebeliReviewing Editor; University of British Columbia, Canada
eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.
Thank you for submitting your work entitled "The evolution of distributed sensing and collective computation in animal populations" for peer review at eLife. Your submission has been favorably evaluated by Ian Baldwin (Senior editor), a Reviewing editor, and two reviewers.
The reviewers are enthusiastic about your paper, and we would like to invite you to revise your manuscript according to the suggestions made by the reviewers.
Summary:
This paper considers a model for collective motion in the presence of a timevarying resource environment, in which the characteristics of future generations of the agents evolve according to how well the current generation obtained resources. The model is based on two behavioral rules: a social response rule and an environmental response rule. It is a forcebased model, in which the interactions with other agents and the environment are captured in terms of potentials that encode these rules. Parameters that are allowed to evolve include an individual's preferred speed when there is no environmental cue, and the sensitivity of an individual to an environmental cue.
It is found that the parameters evolve to robust evolutionarily stable states, and these are typically places in parameter space near transitions to other collective states. This allows groups to respond quickly to changes in the environment. It is argued that such changes in the collective state resemble a phase transition in physical systems.
The manuscript presents an important and substantial advance in understanding the evolution of collective motion of animals.
Reviewer #1:
This paper considers a model for collective motion in the presence of a timevarying resource environment, in which the characteristics of future generations of the agents evolve according to how well the current generation obtained resources. The model is based on two behavioral rules: a social response rule and an environmental response rule. It is a forcebased model, in which the interactions with other agents and the environment are captured in terms of potentials that encode these rules. Parameters that are allowed to evolve include an individual's preferred speed when there is no environmental cue, and the sensitivity of an individual to an environmental cue.
It is found that the parameters evolve to robust evolutionarily stable states, and these are typically places in parameter space near transitions to other collective states. This allows groups to respond quickly to changes in the environment. It is argued that such changes in the collective state resemble a phase transition in physical systems.
The paper has a nice combination of results from computations and theory, much of the details of which are described in the Appendix. Overall it is very well written and describes important properties of collective motion systems. I feel that it is a very strong contribution to the literature.
My only "major" comment is the following:
It is common for phase transitions to be characterized in terms of critical exponents that describe how quantities of interest scale near the phase transition. Is it possible to do this for the system in the paper? This would help to strengthen the somewhat vague claim that the changes in the collective state "resemble phase transitions".
Reviewer #2:
In my opinion, the manuscript "The evolution of distributed sensing and collective computation in animal populations" presents an important and substantial advance in understanding the evolution of collective motion of animals and could be published in eLife. The authors build up on the relatively wellstudied dynamics of nonevolving swarms of interacting particles, introduce important empiricallyderived features, and create a comprehensive evolving model which reproduces several common and intuitivelyappealing behavioural patterns of such animal groups. In addition to a thorough modeling effort, I especially appreciate the extensive analytic and simplified computational estimates of each important facet of the observed dynamics. Such highlyfocused "submodels" help a reader to verify that the contributions of basic mathematical and biological mechanisms to the rather complex observed phenomenology are understood correctly.
A major issue that naturally comes to one's mind while reading the manuscript is how the results would be affected by adding a realistic feature of a gradual depletion of welllocalized resources by swarming consumers. How a rate of such depletion would affect the evolutionary patterns? In the present form, the manuscript is already loaded with new insights, hence, I think, that a qualitative and perhaps even speculative discussion of this topic would suffice.
Reviewer #2 (Minor Comments):
The analogy to the 1order phase transition based on the density contrast and hysteresis curve could be developed a step further: In traditional firstorder phase transitions the hysteresis is caused by the instability of smallsize nuclei of an emerging phase. Could something similar be said about the transition between dispersed and cohesive states?
The word "criticality" seems to become too fashionable and abused in many not directly relevant contexts. In its true meaning, a critical state and critical point terminate the line of firstorder phase transition rather than sits "close to it". In addition, it is characterized by many unique properties (such as nonanalytic behaviour of potentials, scalefree correlations, etc.) which are not mentioned and apparently are irrelevant here. So I would recommend to simply refer to such evolutionary stable state as to one being close to the localizeddelocalized transition line.
Same with the word "computation", which, in my understanding, describes some mathematical actions. I think a more appropriate word to describe the collective behaviour would be "collective sensing".
Unless dictated by the journal format requirements, it would be easier on a reader to have the description of the resource dynamics in the Model Development section and to get rid of the Materials and methods one.
The captions to some figures, especially Figure 4 are very difficult to understand. Unless desperately pressed by the size limit, it would definitely help to expand the caption and even split the figure into several ones.
Would it be possible to embed the movies into the main text somewhere near the description of three typical dynamical regimes?
https://doi.org/10.7554/eLife.10955.031Author response
Reviewer #1:
[…] It is common for phase transitions to be characterized in terms of critical exponents that describe how quantities of interest scale near the phase transition. Is it possible to do this for the system in the paper? This would help to strengthen the somewhat vague claim that the changes in the collective state "resemble phase transitions".
We appreciate the reviewer’s point that the analogy between sharp transitions in density/potential energy we observe in collectives, and phase transitions in physical systems could be clarified further. We observe that agents undergo what appears to be a discontinuous change in local density as a function of $\text{\Psi}$, so the observed transition more closely resembles a first order phase transition than a second order phase transition. Because of this, we do not expect the kind of power law scaling near the transition point that is a characteristic of second order transitions, and the search for a critical exponent would not be meaningful. On the other hand, hysteresis is a hallmark first order phase transitions so that Figures 3 and Appendix Figure 9 provide evidence that the analogy to phase transitions is sound. Additionally, we have expanded discussion of the hysteresis behavior present in our model and added additional numerical results (Appendix section 5.4) to further clarify the connection between the behavior of our system and first order phase transitions in physical systems (subheading “Evolved populations are poised near abrupt transitions in collective state”, see also the response to reviewer #2’s comment).
Reviewer #2:
[…] A major issue that naturally comes to one's mind while reading the manuscript is how the results would be affected by adding a realistic feature of a gradual depletion of welllocalized resources by swarming consumers. How a rate of such depletion would affect the evolutionary patterns? In the present form, the manuscript is already loaded with new insights, hence, I think, that a qualitative and perhaps even speculative discussion of this topic would suffice.
To address the reviewer’s comment, we performed two additional sets of evolutionary simulations under alternative assumptions about resource depletion. In the first alternative scenario, we assumed that individuals deplete the resource very quickly. As one might expect, social cues provide little useful information about the locations of resources when resources are quickly depleted, and as a consequence, the evolutionary trajectories under fast resource depletion resemble those of asocial individuals. Under more moderate levels of resource depletion, however, the evolutionary trajectories shown in Figure 1 of the main text are largely unchanged; social individuals evolve to an evolutionarily stable state in which $\text{\Psi}$0 values are near the transition point in collective state when the resource is at a low value, and individuals cross the transition shown in Figure 3 in regions where the resource value is high. We have added a discussion of these new results to the manuscript (paragraph three, subheading “Evolved populations are poised near abrupt transitions in collective state”).
Reviewer #2 (Minor Comments):
The analogy to the 1order phase transition based on the density contrast and hysteresis curve could be developed a step further: In traditional firstorder phase transitions the hysteresis is caused by the instability of smallsize nuclei of an emerging phase. Could something similar be said about the transition between dispersed and cohesive states?
As the reviewer mentioned, when $\text{\Psi}$ is increased from a low value to a higher value, in the transitional regime, small nuclei of tightly clustered individuals coexist with individuals in the dispersed state. These nuclei eventually become completely unstable at large enough values of $\text{\Psi}$. We have added additional numerical results describing how the rate of nucleus formation depends on $\text{\Psi}$ in the hysteresis region and described these in subheading “Evolved populations are poised near abrupt transitions in collective state” of the revised manuscript.
The word "criticality" seems to become too fashionable and abused in many not directly relevant contexts. In its true meaning, a critical state and critical point terminate the line of firstorder phase transition rather than sits "close to it". In addition, it is characterized by many unique properties (such as nonanalytic behaviour of potentials, scalefree correlations, etc.) which are not mentioned and apparently are irrelevant here. So I would recommend to simply refer to such evolutionary stable state as to one being close to the localizeddelocalized transition line.
We have removed the term “criticality” throughout the manuscript to avoid confusion about terminology.
Same with the word "computation", which, in my understanding, describes some mathematical actions. I think a more appropriate word to describe the collective behaviour would be "collective sensing".
While we agree with the reviewer that the term “computation” is often misused or used casually, in this manuscript we have strived to define precisely what we mean when we use this term (subheading “Changes in collective state allow for rapid collective computation of the resource distribution” and paragraph three, Discussion). One of our motivations for doing this is to help develop a more rigorous definition of “collective computation” as it is a term that is frequently used in the literature on collective behavior.
Unless dictated by the journal format requirements, it would be easier on a reader to have the description of the resource dynamics in the Model Development section and to get rid of the Materials and methods one.
As far as we can tell, the Material and methods section is required by eLife. If it is allowed, we would be willing to move the details of the resource environment to the Model Development section and omit Materials and methods altogether.
The captions to some figures, especially Figure 4 are very difficult to understand. Unless desperately pressed by the size limit, it would definitely help to expand the caption and even split the figure into several ones.
We have divided Figure 4 into two figures and worked to clarify the other figure captions.
Would it be possible to embed the movies into the main text somewhere near the description of three typical dynamical regimes?
We have referenced Video 2 in the section of the text that describes dynamical regimes (subheading “Evolved populations are poised near abrupt transitions in collective state”).
https://doi.org/10.7554/eLife.10955.032Article and author information
Author details
Funding
James S. McDonnell Foundation
 Andrew M Hein
National Science Foundation (PHY0848755, IOS1355061, and EAGER IOS1251585)
 Iain D Couzin
Army Research Office (W911NG1110385 and W911NF1410431)
 Iain D Couzin
Office of Naval Research Global (N000140911074 and N000141410635)
 Iain D Couzin
Human Frontier Science Program (RGP0065/2012)
 Iain D Couzin
National Science Foundation (Dimensions of Biodiversity OCE1046001)
 George I Hagstrom
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was partially supported by National Science Foundation (NSF) Grants PHY0848755, IOS1355061, and EAGER IOS1251585; Office of Naval Research Grants N000140911074 and N000141410635; Army Research Office Grants W911NG1110385 and W911NF1410431; Human Frontier Science Program Grant RGP0065/2012 (to I.D.C.), NSF Dimensions of Biodiversity grant OCE1046001, and a James S McDonnell Foundation Fellowship (to A.M.H.).
Reviewing Editor
 Michael Doebeli, University of British Columbia, Canada
Publication history
 Received: August 20, 2015
 Accepted: November 1, 2015
 Accepted Manuscript published: December 10, 2015 (version 1)
 Accepted Manuscript updated: December 17, 2015 (version 2)
 Version of Record published: February 3, 2016 (version 3)
Copyright
© 2015, Hein et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,996
 Page views

 672
 Downloads

 38
 Citations
Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.