## Abstract

Nature is rich with many different examples of the cohesive motion of animals. Previous attempts to model collective motion have primarily focused on group behaviours of identical individuals. In contrast, we put our emphasis on modelling the contributions of different individual-level characteristics within such groups by using stochastic asynchronous updating of individual positions and orientations. Our model predicts that higher updating frequency, which we relate to perceived threat, leads to more synchronized group movement, with speed and nearest-neighbour distributions becoming more uniform. Experiments with three-spined sticklebacks (*Gasterosteus aculeatus*) that were exposed to different threat levels provide strong empirical support for our predictions. Our results suggest that the behaviour of fish (at different states of agitation) can be explained by a single parameter in our model: the updating frequency. We postulate a mechanism for collective behavioural changes in different environment-induced contexts, and explain our findings with reference to confusion and oddity effects.

## 1. Introduction

The ubiquitous features observed in animal collectives have inspired researchers from a range of disciplines to describe, model and reproduce these extraordinary displays of coordinated behaviour (Sumpter 2006). Most group-living animals are able to move coherently and collectively, preserving common features such as coordinated turns, maintenance of internal structures and apparently leaderless movement. Examples include the tightly bound tori exhibited by large shoals of sardines under predation pressure (Parrish *et al*. 2002) and the striking pre-roosting displays of starlings (Ballerini *et al*. 2008). Despite considerable research interest in group coordination there is still a significant gap between theory and experimental data. Attempts to bridge this gap are hindered by the emergent nature of collective motion (Viscido *et al*. 2004), and matching modelling studies to empirical data—as, for example, in Buhl *et al*. (2006) and Yates *et al*. (2009)—remains a challenging goal in this field.

Models have been central to understanding the mechanisms behind collective animal motion (Krause & Ruxton 2002). Individual-based models in particular have allowed researchers to examine the emergence of different collective behaviours resulting from simple mechanisms at the level of the individual (Krause & Ruxton 2002). These models typically assume that identical individuals react to the position and movement of their nearest conspecifics by a combination of alignment, attraction and repulsion (Krause & Ruxton 2002). Originally based on the extensive work by Aoki (1982), these simple ideas were also adopted in physics (Vicsek *et al*. 1995), computer science (Reynolds 1987) and control engineering (Liu *et al*. 2003). In biology, the connection between the metric inter-individual distance and the subsequent behaviour has given rise to a family of models investigated computationally (Couzin *et al*. 2002) and tested empirically (Tien *et al*. 2004). Models have been successful in shaping explanations and understanding mechanisms for different collective behaviours (Couzin *et al*. 2002, 2005; Hoare *et al*. 2004; Viscido *et al*. 2005; Hemelrijk & Hildenbrandt 2008), but only at a qualitative level.

One of the earliest empirical studies to quantify individual trajectories in collective motion was performed three decades ago (Aoki 1980). In his experiments, Aoki filmed shoals of tamoroko (*Gnathopogon elongatus*) and Japanese horse mackerel (*Trachurus japonicus*) under controlled conditions and extracted time series of the positions of individual fish from his films. Aoki assembled the distribution of speeds and nearest-neighbour distance distributions of individuals within fish shoals. It is surprising that so few models incorporate these findings and no model explains them. Most modelling studies, for example, use a constant and homogeneous speed (e.g. Couzin *et al*. 2002). Some studies have used Aoki's data by drawing an instantaneous speed at each time step from an appropriate distribution (Aoki 1982; Huth & Wissel 1992), but they do not explain the emergence of this distribution from first principles. In this article, we construct a model in which the speed distributions are emergent purely from the local interactions between the group members, and discuss its consequences in the context of experimental work on fish shoals.

## 2. Material and methods

### (a) Computational model

We have developed an individual-based model of group interactions, based on local rules, that replicates the speed distributions found in Aoki's and our experiments. The basis of our approach is to adopt stochastic asynchronous updating of individual fish positions and orientations; rather than using deterministic and sequential updating at each time step, fish can react to external stimuli with a stochastic rate. The average behaviour of individuals over short time intervals then varies probabilistically. Models with asynchronous updating have been previously introduced in simple one- and two-dimensional models (O'Loan & Evans 1999; Liu *et al*. 2003; Raymond & Evans 2006; Şamiloğlu *et al*. 2006), but we believe their potential to explain empirical observations in real animals, such as effectively modelling fish speed distributions, has been overlooked.

Individuals are represented by points on the plane moving continuously in a toroidal space (a square box of side length *L* with periodic boundary conditions). The notation and nomenclature below follows that of Couzin *et al*. (2002). Let *N* be the number of individuals indexed *i*, with positions *x*_{i} = (*x*_{j}, *y*_{j}) and direction of motion *θ*_{i}. We assume that each individual reacts with an identical stochastic rate, enabling us to exploit a particle-picking approach to exactly simulate the implicit underlying master equation of the system. The algorithmic implementation of our model contains aspects of earlier work (Tsitsiklis *et al*. 1986; O'Loan & Evans 1999; Couzin *et al*. 2002) and proceeds as follows:

Choose individual

*j*at random, where*j*∈ {1, …,*N*} (equal probabilities).Decide which one of the two behavioural rules

*j*will follow in this step (probabilities*p*and (1 −*p*), respectively).Update

*x*_{j}and*θ*_{j}according to the behavioural rule chosen in step (2).

*N* realizations of steps (1)–(3) constitute one update step of length Δ*t* seconds. The duration of this update step corresponds to the reciprocal value of the rate at which individuals update. The output of the model is obtained by recording the positions of all individuals every *T* = *λ*Δ*t* seconds, where *λ* ≥ 1. This is analogous to how data of animal motion are obtained empirically where individual positions and orientations are sampled according to the frame rate of video recordings (Aoki 1980; Buhl *et al*. 2006). In our simulations, we keep *T* fixed and only vary Δ*t* and therefore also *λ*.

We use purely metric behaviour rules in this implementation, based on those of Couzin *et al*. (2002). We will defer commenting on the appropriateness of this approach given recent findings (Ballerini *et al*. 2008) until the discussion. Each individual obtains information from interaction zones—zone of repulsion (zor), zone of orientation (zoo) and zone of attraction (zoa)—which are described by concentric circles, centred on the individual, of radius *r*_{R}, *r*_{O} and *r*_{A}, respectively. Both the zoo and the zoa are punctured by a ‘blind angle’, *α*, in which individuals cannot perceive other individuals. Suppose individual *j* has been chosen in the algorithm described above; our first behavioural rule, which is selected with probability *p*, implements either alignment or repulsion. The individual tries to move away from conspecifics within its zor or aligns to conspecifics in its zoo, where, in common with other models of collective motion, repulsive motion takes precedence over alignment. The distance of *j* to its nearest neighbour determines the behaviour. Let *R* ∈ {1, …, *N*} be the set of individuals within the zor of *j*, excluding *j*. If |*R*| ≥ 1, the desired direction of motion of *j* is given by
where angle (** y**) denotes the angle the vector

**makes with the horizontal axis and**

*y*

*r*_{ij}= (

*x*_{i}−

*x*_{j}) is the vector in the direction from

*j*to

*i*. However, if the distance from

*j*to its nearest neighbour is larger than

*r*

_{R}, then

*j*aligns with its conspecifics. Let

*O*∈ {1, …,

*N*} be the set of individuals within the zoo of

*j*, excluding

*j*. If |

*O*| ≥ 1, the desired direction of motion of

*j*is given by

If both *R* and *O* are empty, then , and the individual does not deviate its direction. It executes this move with an instantaneous speed *v* = *v*_{O}.

In the alternative case, we select our second behavioural rule with probability (1 − *p*). In this case, individual *j* gets attracted to conspecifics in its zoa and the distance *r*_{ij} once more determines its behaviour. Let *A* ∈ {1, …, *N*} be the set of individuals within the zoa of *j*, excluding *j*. If |*A*| ≥ 1, the desired direction of motion of *j* is given by

Once more, if *A* is empty, then , and no deviation occurs. Subsequent movement happens at instantaneous speed *v* = *v*_{A}. Throughout this study, we choose *p* = 0.5, in agreement with previous research in which an equal weight is assigned to orientation and attraction in individuals (Couzin *et al*. 2002). For both behavioural rules, once the desired direction of motion for *j* is calculated, the updated direction of motion, *new*(*θ*_{j}), is found by rotating the individual *j* by at most *β*Δ*t* from *θ*_{j} towards . Here, *β* denotes the maximum turning rate for individuals. Every time an individual *j* is updated, it is moved by *v* units in the updated direction
where *v* is selected to be either *v*_{O} (alignment or repulsion) or *v*_{A} (attraction), as described above. The average speed of an individual, over many update steps, is consequently given by *v*_{av} = *pv*_{O} + (1 − *p*)*v*_{A}. Parameters used in the model simulations are as follows: *N* = 8, *L* = 168.73 cm, *T* = 0.04 s, *p* = 0.5, *α* = 270°, *β* = 40°, *v*_{O} = 8.44 cm s^{−1}, *v*_{A} = 2*v*_{O}, *r*_{R} = 5.06 cm, *r*_{O} = 20.25 cm, *r*_{A} = 33.75 cm; values of Δ*t* are given in the figure legends and justified in the electronic supplementary material.

In a simple stochastic implementation, all individuals would have an identical instantaneous speed, independent of the rule they follow and the behaviour of their fellow individuals. This would produce an unskewed Poisson distribution for the individual speeds (when averaged over time) that is not supported by empirical data. The novelty of our implementation is that individuals adopt differing speeds according to the behavioural rule they follow. In such a way, we can obtain skewed distributions for the individual speeds as observed in empirical data (see below).

An inherent parameter in our model is the length of the update step Δ*t* (in seconds). This parameter reciprocally rescales the reaction rates in the system: small values of Δ*t* imply rapid updates, while large values of Δ*t* imply slow updates. It is important to stress at this stage that we are not explicitly relating the size of Δ*t* to biological or neurological reaction times of animals (but discuss the possibility of a connection later in this article). In addition, no direct physical meaning should be attached to the instantaneous positions on time scales close to Δ*t*. The fact that some individual might not move for one or more update steps does not imply they have stopped; rather that they are reacting more slowly to their surrounding than their fellows. At first glance it may appear that the same effect we obtain by varying Δ*t* could be obtained by varying the speed at which individuals move. This is not the case. In our model, different values of Δ*t* do not alter the average speed at which individuals move but they do alter the average rate at which individuals act upon information from within their sensory zones. Imposing different average speeds (i.e. changing *v*_{A} and *v*_{O}) changes the relationship between the average speed of individuals and the extent of their sensory zones.

### (b) Experimental methods

We extended Aoki's experiments (Aoki 1980), using small shoals of eight three-spined sticklebacks, *Gasterosteus aculeatus*, within an indoor circular tank of 1 m radius. From individual movement trajectories, we constructed the distributions of the individual speeds of fish within a shoal and the distribution of the individuals' distances to their nearest shoal mates (nearest-neighbour distances). To test predictions from our model, we designed a number of new experimental treatments (table 1) to produce varying levels of agitation in the fish, comprising most agitation (treatment 1), least agitation (treatment 4) and intermediate agitation levels. This allowed us to compare different model outcomes under different conditions to experimental data under parallel conditions. Evidence suggests that sticklebacks are in a greater state of agitation or excitement in higher light levels since, in experiments, fish of this species have preferred shaded regions in tanks (Ward *et al*. 2008). Fish that perceive a higher predation threat tend to use more shaded areas than fish that do not perceive the same level of threat (McCartt *et al*. 1997). This suggests that they perceive a lower predation threat in shaded areas than in well-lit areas. An explanation for this could be that fish can see approaching predators that are not in shaded areas better from the shade, and that fish are less likely to be seen in shade (Helfman 1981). We also varied the water depth in our experimental tank. Given the white background of the tank, the fish bodies are clearly visible and thus make fish potentially conspicuous to overhead predators, such as kingfishers and herons. In this situation sticklebacks show a strong tendency to move into deeper water (J. Krause 2008, unpublished data).

## 3. Results

### (a) Model output

Our model produces skewed speed distributions (in two dimensions) similar to empirical data as an emergent property of our novel update scheme (figure 1*c*,*d*). Individual speeds are approximated by calculating the distance covered by fish over a fixed time step, with each time step *T* comprising many multiples of Δ*t*. This is analogous to how speed distributions are determined empirically, where fish speeds are averaged over a range of video frames (see electronic supplementary material; see also Aoki 1980). The effect of varying Δ*t* is striking: large values of Δ*t* promote a strongly positively skewed distribution, and small values reduce the skewness and give rise to speed distributions that resemble normal, or Gaussian, distributions (figure 1*c*,*d*). We do not claim to reproduce speed distributions of real fish quantitatively as the influence of important factors on the speed distributions is unknown (e.g. interaction with environment). Rather, we show that our model is capable of producing similar speed distributions to the data without *a priori* assumptions or explicit addition of stochastic noise. Furthermore, our model suggests that the shape of the speed distributions can be varied by changing one parameter in our model.

From the speed and nearest-neighbour distance distributions of the simulated shoals, we extracted three summary statistics: the standard deviation of the speed distributions, skewness of the speed distributions and the median of the nearest-neighbour distance distributions. Substantial changes in the summary statistics for different values of Δ*t* are observed (figure 2*a–c*). Most prominent is the reduction in nearest-neighbour distance (for a given fish, this is the distance between this fish and the fish closest to it). Such a decrease in nearest-neighbour distances, or more compact group structures have previously been observed in empirical experiments for increasing threat or agitation levels (Krause 1993; Tien *et al*. 2004; Carere *et al*. 2009). The effect of Δ*t* on the summary statistics highlights that this parameter is important for biological interpretation and is not just an invisible model implementation parameter.

We propose that values of Δ*t* in our model correspond to states of agitation in animals. For example, low values of Δ*t* (that is, rapid updates) would correspond to high states of agitation, which might occur when animals feel threatened or at risk (Krause & Ruxton 2002). Our model predicts that increasing agitation makes the speed distribution of the shoal become more uniform and causes the nearest-neighbour distribution to contract. This allows us to form the following empirically testable hypothesis: the speed distributions and nearest-neighbour distributions of fish at different levels of agitation should qualitatively correspond to distributions in our model, where Δ*t* is varied appropriately. Specifically, our model predicts that:

— High states of agitation (low values of Δ

*t*) should result in strongly peaked, unskewed speed distributions and a contraction of nearest-neighbour distances.— Low states of agitation (high values of Δ

*t*) should result in well-spread distributions with positive skew and an increase in nearest-neighbour distances.

### (b) Empirical findings

Using our empirical system, we confirmed previous results (Aoki 1980) in finding long-tailed and positively skewed speed distributions (figure 1*a*,*b*). We then extracted the same three summary statistics from the empirical data as we did for the simulated data. To investigate the statistical significance of differences in the measurements of the summary statistics across treatments, we used a generalized linear mixed model (GLMM), taking into account the differences between shoals and the order in which the treatments were applied. In our analysis of the empirical data, we found that all three summary statistics were affected by one or more of the treatments in a statistically significant way (figure 2). Statistically significant differences between treatment 1 and treatments 2 and 3 (e.g. figure 2*d*,*f*) illustrate that water depth and light intensity can separately affect the animal's movement patterns. The lack of monotony in some of the trends is due to behavioural factors we cannot control, which are discussed in the electronic supplementary material. The fact that not all of the summary statistics show statistically significant differences between treatment 1 and treatments 2 and 3 is likely to be due to the fact that the contrast between these treatments is not large enough.

Overall, our experimental findings confirmed the predictions from our model that increasing agitation in fish makes the speed distribution of the shoal become more uniform (i.e. it decreases the distribution's standard deviation and skewness, and makes the nearest-neighbour distribution contract; figure 2).

## 4. Discussion

This study is a combined modelling and empirical effort that has successfully predicted and reproduced emergent empirical properties of coordinated group behaviour from a model based entirely on local interactions. Our model is relatively simple and therefore provides an ideal starting point for the inclusion of individual characteristics and excitement levels into models of collective motion based on stochasticity. Our model produces novel predictions as to how group properties will alter in different behavioural contexts, and our experiments provide supporting evidence for these predictions. This reveals the importance of threat or risk levels perceived by fish for the composition of their movement trajectories and coordination. It has been suggested that fish react more quickly to shoal mates in situations of higher perceived risk or threat levels (Ward *et al*. 2008). However, to our knowledge this is the first time that the updating frequency of individuals has been modelled and tested against empirical data explicitly.

Our parameter Δ*t* is the mean inter-update time that captures the relative frequency of updates within a given sampling time frame. It is important to emphasize that we are not making an explicit claim that this parameter is derivable from neurological information; we regard this parameter as controlling the dynamic averaging of the positional information from nearby conspecifics. The precise mechanism and quantities in our model provide an interesting avenue to be tested in further empirical research, focused on understanding the physiological interpretation of the stochastic rates we find emerging from our qualitative model comparisons. The simulations presented in figures 1 and 2 use values of Δ*t* that are significantly lower than the 0.1 s that has been previously recorded as the reaction time for fish (Partridge & Pitcher 1980), indicating that many multiples of Δ*t* make up a responsive reaction from the organism. The effect of reducing the step size in algorithms such as ours has previously been considered, but not in a biological context (Tsitsiklis *et al*. 1986).

It has recently been argued that the behavioural rules of collectively moving animals are based on the number of conspecifics each individual tracks (‘topological framework’) rather than on the distance between individuals as in our model (Ballerini *et al*. 2008). Some of the phenomena Ballerini and co-workers observed in their data have been reproduced in extensive simulation studies (Hildenbrandt *et al*. 2009) by assuming *a priori* that individuals only interact with a limited number of shoal mates. However, we suggest that this is not necessarily the only way in which the observations made in Ballerini *et al*. (2008) may arise in a model. Furthermore, a simple implementation of a topological framework in which individuals only interact with a fixed number of their nearest neighbours would not affect the emergence of speed distributions in our model. For these reasons we have continued to use a distance-based approach.

We suggest that by moving in a more coherent fashion with shoal members, an individual is able to reduce the risk of being targeted by predators as the ‘odd one out’, often termed the oddity effect (Krause & Ruxton 2002). The confusion effect—where predators find it more difficult to target an individual in a group than to target an isolated individual—is easily broken if one individual differs morphologically or behaviourally from others (Krause & Ruxton 2002). For example, in a threatened group where nearest-neighbour distances are generally low, an individual with a large nearest-neighbour distance will stand out from the crowd and probably be targeted by predators. This provides a mechanistic explanation for our findings: greater risk produces higher updating frequencies and higher updating frequencies produce lower oddity. Therefore, we suggest that the oddity effect could be the driving force for the behavioural changes in different contexts and the high degree of synchrony characterizing threat-induced collective behaviours.

Finally, our method of measuring the uniformity of speed distributions and nearest-neighbour distances could provide a simple way of empirically assessing stress levels of collectively grouping animals in a remotely collectable and non-obtrusive way.

## Acknowledgements

N.W.F.B.'s research is supported by the Natural Environment Research Council. J.J.F. acknowledges support from the Biotechnology and Biological Sciences Research Council. D.W.F. and A.J.W. are supported by RCUK Fellowships. D.W.F. acknowledges support from NERC grant no. NE/E016111/1. J.K. acknowledges support from NERC (NE/D011035/1).

## Footnotes

- Received April 22, 2010.
- Accepted May 6, 2010.

- © 2010 The Royal Society