Abstract

Through combining theoretical models and empirical data, complexity science has increased our understanding of social behavior of animals, in particular of social insects, primates, and fish. What are missing are studies of collective behavior of huge swarms of birds. Recently detailed empirical data have been collected of the swarming maneuvers of large flocks of thousands of starlings (Sturnus vulgaris) at their communal sleeping site (roost). Their flocking maneuvers are of dazzling complexity in their changes in density and flock shape, but the processes underlying them are still a mystery. Recent models show that flocking may arise by self-organization from rules of co-ordination with nearby neighbors, but patterns in these models come nowhere near the complexity of those of the real starlings. The question of this paper, therefore, is whether such complex patterns can emerge by self-organization. In our computer model, called StarDisplay, we combine the usual rules of co-ordination based on separation, attraction, and alignment with specifics of starling behavior: 1) simplified aerodynamics of flight, especially rolling during turning, 2) movement above a “roosting area” (sleeping site), and 3) the low fixed number of interaction neighbors (i.e., the topological range). Our model generates patterns that resemble remarkably not only qualitative but also quantitative empirical data collected in Rome through video recordings and position measurements by stereo photography. Our results provide new insights into the mechanisms underlying complex flocking maneuvers of starlings and other birds.

More and more often, it has been shown that social behavior of groups of animals can be understood with the help of complexity science (Camazine et al. 2001; Hemelrijk 2002). Combining models of self-organization with empirical data has helped to explain diverse behavioral aspects of several taxa (Hemelrijk 2005). For instance, the foraging and grouping behavior of social insects (Gautrais et al. 2003; Dussutour et al. 2004; Ame et al. 2006) and of primates (te Boekhorst and Hogeweg 1994a, 1994b; Hemelrijk et al. 2008; Puga-Gonzalez et al. 2009) and the co-ordination in swarms of insects (Buhl et al. 2006), of fish (Becco et al. 2006; Ward et al. 2008; Hemelrijk et al. forthcoming), and of birds in a small flock (Biro et al. 2006). However, particularly, studies of self-organization of large flocks of birds are missing. For instance, what causes the complex flocking maneuvers of tens of thousands of European Starlings (Sturnus vulgaris) over the sleeping site at night has not been solved. These aerial flocking maneuvers take place in winter in the period before migration when starlings come together with tens of thousands of individuals at their sleeping sites. Before settling in the trees for the night, they circle above the roost for about half an hour, although their flocks split and merge and change in shape and density (Brodie 1976; Feare 1984; Beauchamp 1999). An unlikely explanation was given by Selous (1931): He attributed their wheeling and turning, splitting, and merging to “telepathy”. The aim of this paper is to study whether patterns similar to aerial displays can be generated in models of self-organization. We validate our model, called StarDisplay, by comparing its patterns to empirical data recently collected in Rome. These data consist in extensive video recordings of these spectacular patterns (Carere et al. 2009) and the 3-dimensional positions of all individuals of 10 huge flocks consisting of thousands of starlings that were obtained using stereo photography (Ballerini et al. 2008a, 2008b; Cavagna et al. 2008a, 2008b). Such data are very difficult to collect because of the huge number of individuals and the large distance of the flocks from the observer (Major and Dill 1978; Pomeroy and Heppner 1992; Budgey 1998). The study in Rome has revealed a number of interesting new traits (Ballerini et al. 2008a). For instance, that the flock banks as a whole in the direction of banking by the individuals, that nearest neighbors are usually located at the side, that the shape of flocks is diverse as regards its orientation in the movement direction (e.g. whether it is longer than wide or vice versa) and that flock density is similar between back and front (Ballerini et al. 2008b) and that flock density is independent of flock size (Ballerini et al. 2008b). The latter 3 traits differ from descriptions of fish schools in empirical data and also in models. In fish schools, we typically find that schools have 1) an oblong shape, 2) high frontal density (Bumann et al. 1997), and that 3) group density depends on the number of individuals (Kunz and Hemelrijk 2003; Hemelrijk and Kunz 2005; Hemelrijk and Hildenbrandt 2008; Hemelrijk et al. forthcoming).

The usual biologically oriented models of swarming are tuned toward fish, in that they reproduce single schools with these 3 traits, and they lack the complexity of the flocking maneuvers of starlings. These models are based on 3 rules: individuals avoid those that are too nearby, they align with those at intermediate distance, and are attracted to those further away (Reynolds 1986; Huth and Wissel 1992, 1994; Reuter and Breckling 1994; Romey 1996; Parrish and Edelstein-Keshet 1999; Hemelrijk 2002; Couzin et al. 2005; Helbing et al. 2005; Parrish and Viscido 2005).

To generate starling-like patterns in a model, we have in the present study extended our former model of fish schools (Hemelrijk and Hildenbrandt 2008) with aspects characteristic of starlings. For instance, in the new model, individuals, like in real starlings, interact with others in their so-called “topological range”, that is, with the nearest 6–7 neighbors, independently of the absolute distance to them (Ballerini et al. 2008a, 2008b). They fly at cruise speed (Videler 2005) from which they sometimes deviate to avoid collisions and to rejoin others. Their flocks remain above a sleeping site at a certain height (Ballerini et al. 2008a; Carere et al. 2009) because we give individuals, once they are outside the area above the sleeping site, a horizontal attraction to return to that area and a vertical attraction to a preferred height. As regards their flying behavior, individuals in StarDisplay experience lift, drag, and the force of gravity. They fly in the same way as supposed in some other studies of bird flight, following theories of fixed wing dynamics (Norberg 1990). When flying along a curve, like real birds they roll into the direction of the turn until they are at a certain angle to the horizontal plane, the so-called banking angle (Videler 2005). Their tendency to roll into the turn increases with the strength of their tendency to turn sideward, which is due to the urge to co-ordinate with their topological neighbors and to stay above the roost. Once an individual has banked in the model, its tendency to roll back is proportional to its actual banking angle.

In order to verify whether flocking maneuvers of real starlings may arise by self-organization, the aim of the present study is to investigate the resemblance between flocking patterns in the model and empirical data in Rome. We study this resemblance in 2 ways. First, we compare the model to qualitative empirical data (Ballerini et al. 2008b; Carere et al. 2009). Second, we compare it to quantitative data collected by stereo photography described in Table 1 of the paper by Ballerini et al. (2008a, 2008b). Finally, we discuss what may cause these patterns in the model and how these explanations could be studied through experiments in the model and collection of new empirical data in future.

Table 1

Default parameter valuesa

Parameter Description Default value 
Δt Integration time step 5 ms 
Δu Reaction time 50 ms 
v0 Cruise speed 10 m/s 
M Mass 0.08 kg 
CL/CD Lift-drag coefficient 3.3 
Lo Default lift 0.78 N 
D0, T0 Default drag, default thrust 0.24 N 
wβin Banking control 10 
wβout Banking control 
T Speed control 1 s 
Rmax Maximum perception radius 100 m 
nc Topological range 6.5 
S Interpolation factor 0.1 Δu 
rh Radius of maximum separation (hard sphere) 0.2 m 
rsep Separation radius 4 m 
Σ Parameter of the Gaussian g(x1.37 m 
ws Weighting factor separation force 1 N 
 Rear “blind angle” cohesion and alignment 90° 
wa Weighting factor alignment force 0.5 N 
wc Weighting factor cohesion force 1 N 
Cc Critical centrality below which an individual is assumed to be in the interior of a flock 0.35 
wξ Weighting factor random force 0.01 N 
RRoost Radius of roost 150 m 
wRoostH Weighting factor horizontal attraction to the roost 0.01 N/m 
wRoostV Weighting factor vertical attraction to the roost 0.2 N 
Parameter Description Default value 
Δt Integration time step 5 ms 
Δu Reaction time 50 ms 
v0 Cruise speed 10 m/s 
M Mass 0.08 kg 
CL/CD Lift-drag coefficient 3.3 
Lo Default lift 0.78 N 
D0, T0 Default drag, default thrust 0.24 N 
wβin Banking control 10 
wβout Banking control 
T Speed control 1 s 
Rmax Maximum perception radius 100 m 
nc Topological range 6.5 
S Interpolation factor 0.1 Δu 
rh Radius of maximum separation (hard sphere) 0.2 m 
rsep Separation radius 4 m 
Σ Parameter of the Gaussian g(x1.37 m 
ws Weighting factor separation force 1 N 
 Rear “blind angle” cohesion and alignment 90° 
wa Weighting factor alignment force 0.5 N 
wc Weighting factor cohesion force 1 N 
Cc Critical centrality below which an individual is assumed to be in the interior of a flock 0.35 
wξ Weighting factor random force 0.01 N 
RRoost Radius of roost 150 m 
wRoostH Weighting factor horizontal attraction to the roost 0.01 N/m 
wRoostV Weighting factor vertical attraction to the roost 0.2 N 
a

Note that D0 and T0 is calculated by Equation 15bc by inserting v0 for vi.

METHODS

The model

The behavior of each individual in StarDisplay is based on its cruise speed, its social environment (i.e., the position and heading of its nearby neighbors), its attraction to the roost, and the simplified aerodynamics of flight that includes banking while turning. Because flying implies movement in all directions, our model is 3 dimensional. We built the model in SI units and choose real parameter values where available (Table 1).

Even though social co-ordination depends on internal motivations, we model it in terms of (social) forces as has been done by others (Reynolds 1987; Couzin et al. 2002; Hemelrijk and Hildenbrandt, 2008).

Details of behavioral rules

Each individual is characterized by its mass, m, its speed, v, and its location, p. Its orientation in space is given by its local co-ordinate system (ex, ey, and ez). Following the model by Reynolds (1987), its orientation is indicated by its forward direction, ex, its sideward direction, ey, and its upward direction, ez, which it changes by rotating around these 3 principal axes, ex, ey, and ez (roll, pitch, and yaw) (Figure 1).

Figure 1

The local co-ordinate systems of 2 birds with different orientations in space and at different distances to the viewer. ex, its the bird's forward direction; ey, its sideward direction; and ez, its upward direction. It can change these by rotating around these 3 principal axes (roll, pitch and yaw).

Figure 1

The local co-ordinate systems of 2 birds with different orientations in space and at different distances to the viewer. ex, its the bird's forward direction; ey, its sideward direction; and ez, its upward direction. It can change these by rotating around these 3 principal axes (roll, pitch and yaw).

As to its speed, a force, fτi (Equation 1), brings an individual back to its cruise speed v0 after it has deviated from it (Hemelrijk and Hildenbrandt 2008). 

speed control (1)
graphic
where τ represents the relaxation time, m is the mass of the individual i and v0 its cruise speed, vi is its velocity, and exi its forward direction.

To represent topological interaction, each individual i in the model adapts its metric interaction range, Ri(t) (Hemelrijk and Hildenbrandt 2008) in such a way that individuals attempt to interact only with a constant number of their closest neighbors (Equations 2 and 3). 

adaptive interaction range (2)
graphic
 
neighborhood of an individual (3)
graphic
where Δu is the time span to react, s is an interpolation factor, Rmax is the maximal metric interaction range, Ni(t) is the neighborhood of individual i at time t, that is, the set of neighbors of an individual i which is composed of |Ni(t)| neighbors from the total flock of size N, nc is the fixed number of topological interaction partners it strives to have, and dij is the distance between individual i and j given by ‖pjpi‖, where pi gives the position of an individual i. Thus, the radius of interaction at the next step in reaction time, Ri (t + Δu), is increased if the number of interaction partners |Ni(t)| is smaller than the targeted number nc, and it is decreased if it is larger than that; it remains as before if |Ni(t)| equals nc. Here, Ri can neither decrease below the hard sphere in which individuals are maximally avoiding each other rh (Equation 4, Figure 2) nor increase beyond Rmax and s, the interpolation factor, which determines the step size of the changes and herewith, the variance of the number of actual influential neighbors.

Figure 2

Social interaction ranges for separation (A), cohesion (B), and alignment (C). Note that the values of the different radii in the figure are not to scale with the default values in Table 1.

Figure 2

Social interaction ranges for separation (A), cohesion (B), and alignment (C). Note that the values of the different radii in the figure are not to scale with the default values in Table 1.

As to separation, individual i is led by a force fsito move in the opposite direction of the average direction of the locations of the Ni(t) others in its neighborhood (Figure 2). Following other scientists (Couzin et al. 2002; Zheng et al. 2005), we have omitted the blind angle at the back (Equation 4 and see Parameterization, initial condition, and experiments). We gave individuals a hard sphere with radius rh as mentioned above in which they are not attracted to each other (Equation 5). Outside the hard sphere, but inside the radius of separation rsep, the degree of avoidance of others decreases with the distance to the neighbor following a halved Gaussian, g(dij), with σ the standard deviation of the Gaussian set such that at the border of the separation zone the force is almost zero, g(rsep) = 0.01 (Equation 4). 

separation (4)
graphic
Here, |Ni(t)| is the number of individuals in the neighborhood of interaction (Equation 3) and dij is the distance from individual i to individual j. The direction from individual i to individual j is specified by the unit vector dij = (pjpi)/‖pjpi‖ and ws is the fixed weighting factor for separation (Table 1). Note that if a flock is thinly spread out, some individuals will only be subject to alignment and cohesion because they are located outside the range of separation of others.

As to cohesion (Equations 5, 6 and 7), individual i is attracted by a force fcito center of mass (i.e., the average x, y, z position) of the group of Ni* (t) individuals located in its topological neighborhood. Following models of others (Huth and Wissel 1992; Reuter and Breckling 1994; Couzin et al. 2002; Kunz and Hemelrijk 2003; Hemelrijk and Kunz 2005), we excluded individuals located in its blind angle at its back. Here, wc is the weighing factor for cohesion (Equation 6, Table 1). Within the radius of the hard sphere rh, we ignore cohesion to others. We make individuals cohere more strongly when they are at the border of the flock than in its interior by multiplying the force of cohesion by a factor indicating the degree to which an individual is peripheral (Equations 5 and 7). To measure this, we use its degree of “centrality” in the group, Ci, calculated as the length of the average vector of the direction toward its neighbors NG(t) (Hemelrijk and Wantia 2005). The centrality of the individual Ci(t), we base on all its “neighboring” individuals NG(t) in a radius of twice its actual perceptual distance (Equation 7). 

cohesion (5)
graphic
 
reduced neighborhood (6)
graphic
 
centrality (7)
graphic

As regards its alignment behavior (Equation 8), individual i perceives a force, fai, to align with the average forward direction of its Ni* (t) interaction neighbors (the same neighbors as to whom it is attracted, Figure 2b,c). 

alignment (8)
graphic

Here, exi and exj are the vectors indicating the forward direction of individuals i and j and wa is the fixed weighting factor for alignment (Table 1).

The “social force” is the sum of these 3 forces (Equation 9). 

social force (9)
graphic

Individuals fly at a similar height above the sleeping site because we made them experience both in a horizontal and vertical direction a force of attraction fRoost to the “roosting area” (Equations 10, 11, and 12). The strength of the horizontal attraction, fRoostH, is greater, the more radial it moves away from the roost; it is weaker if it is already returning. The sign in Equation 11 is chosen to reduce the outward heading. The actual direction of the horizontal attraction force is given by eyi, which is the individual's lateral direction. Vertical attraction, fRoostV, is proportional to the individual's vertical distance from the preferred height above the roost (arbitrarily called the 0 level). Here z is the vertical unit vector. wRoostH and wRoostV are fixed weighting factors. 

roost  attraction  (10)
graphic
 
horizontal  (11)
graphic
 
vertical  (12)
graphic

The random force indicates unspecified stochastic influences (Equation 13) with ξ being a random unit vector from a uniform distribution and wξ being a fixed scaling factor. The sum of the social force, the forces that control speed and ranging, and the random force is labeled as “steering force” (Equation 14). 

random  force  (13)
graphic
 
steering  force  (14)
graphic

Physics of flight in the model follows the standard equations of fixed wing aerodynamics. These link the lift L, the drag D, and the thrust T produced by a bird to its current speed v (Equation 15): 

lift and drag (15a)
graphic
 
lift and drag at cruise speed v0(15b)
graphic
 
simplified lift and drag (15c)
graphic
where ρ is the density of the air and S represents the wing area of the bird. CL and CD are the dimensionless lift and drag coefficients. We assume that ρ is a constant and that the wing areas of all birds in the model are identical, and the ratio of the coefficients CL and CD is constant (Norberg 1990). When a bird is flying horizontally while maintaining a constant cruise speed v0, its lift balances its weight mg (i.e., mass × gravity) and its thrust balances its drag (Equation 15b). Division of L by L0 and of D by L in Equation 15a and 15b yields Equation 15c in which the lift and the drag depend on fixed values and only one variable, i.e. the actual speed.

By definition, gravity is directed toward the global “down” direction, g = (0, 0, −g), the lift toward the local “up” direction ez of the bird, and the drag is pointing in the direction opposite to the individual's actual “forward” direction ex. Thus, the flight forces are 

(16)
graphic

To represent banked turns, we first calculate the degree to which individuals want to turn, that is, their lateral acceleration, al, which is exerted by the steering force (Figure 3). Banking implies that the individual’s local frame rolls around the forward axis in the direction of the lateral acceleration, al. The lateral acceleration follows the first law of Newton (F = ma) 

lateral acceleration (17)
graphic
 
roll in (18)
graphic
 
roll  out  (19)
graphic
 
banking  angle  (20)
graphic
where βi is the actual banking angle, wβin and wβout, respectively, are the weights for rolling in and out the curve of turning, Δt is the update time and βin and βout are the angles over which an individual intends to move inwards and outwards. The degree of lateral acceleration determines the extent of inwards banking into the turn. The intended inward angle is given as the product of the weight for inward turning, the strength of the lateral acceleration, and the time interval over which it was updated (Equation 18, top row in Figure 3). To make the individual return again to level flight, it experiences also a tendency to roll out, which is larger when the actual banking angle is greater (Equation 19). The actual banking angle (Equation 20) is the sum of the current angle and the tendencies to roll in and to roll out (bottom row, Figure 3). The ratio of wβin and wβout determines the roll rate. Note that the acceleration, which was originally strictly lateral, has after the roll, obtained a component of magnitude ‖al‖ sin(βin) parallel to the lift (dotted lines in the last column, Figure 3).

Figure 3

Rotation of the individual around its forward direction during banked turns. The forward direction of the individual ex is pointing toward the reader. Top row, left: at the beginning of the time step (T = 0), the individual flies at level flight and it experiences a steering force to turn with a sideward horizontal component al. At the end of the time step (top row, right), it banks inwards over an angle β in order to take the turn. Bottom row, left: At the next time step, the individual experiences a tendency to both turn inwards βin and outwards βout. Its actual turning angle is given by the sum of its present angle β' and both tendencies to turn inwards β'in and outwards β'out. For illustration, we used an imaginary value for wβin * Δt ∼1 that differs from the default value. For further description see text.

Figure 3

Rotation of the individual around its forward direction during banked turns. The forward direction of the individual ex is pointing toward the reader. Top row, left: at the beginning of the time step (T = 0), the individual flies at level flight and it experiences a steering force to turn with a sideward horizontal component al. At the end of the time step (top row, right), it banks inwards over an angle β in order to take the turn. Bottom row, left: At the next time step, the individual experiences a tendency to both turn inwards βin and outwards βout. Its actual turning angle is given by the sum of its present angle β' and both tendencies to turn inwards β'in and outwards β'out. For illustration, we used an imaginary value for wβin * Δt ∼1 that differs from the default value. For further description see text.

After summing the forces of steering and flying, we use Euler integration to calculate the position and velocity at the end of each time-step Δt: 

(21)
graphic
 
(22)
graphic
where vi is the speed of individual i, m its mass, pi its location, and Δt is the update time.

Parameterization, Initial condition, and Experiments

As far as data are available, we parameterize StarDisplay to realistic values (Table 1). Values were available for starlings regarding body mass, lift/drag ratio, CL/CD, cruise speed (Ward et al. 2004), and reaction time (Pomeroy and Heppner 1977). We tune the interpolation factor s (Equation 2) to obtain the topological range |Ni(t)| of starlings of nc = 6.5 ± 1 neighbors (Ballerini et al. 2008a).

By tuning the weights for repulsion, alignment, and attraction, we adjust the steering force to the order of magnitude of the physical forces of flight, that is, around 1 Newton (L0 and D0, Table 1). The radius of the hard sphere, rh, equals half the starlings' average wingspan (Ward et al. 2004; Möller 2005; Ballerini et al. 2008b). We tune the weight of the separation force, ws, such that no collisions occur (Potts 1984). We adjust the separation radius, rsep, to obtain an average nearest neighbor distance similar to that of the starlings of Rome (Table 2, Supplementary Material Figure S1).

Table 2

Quantitative properties of (empirical and) modeled flocksa

Flock/flocking event Number of birds rsep(m) Cruise speed (m/s) Velocity(m/s) NND (m) Balance shift Volume (m3Thickness I1 (m) Aspect ratios
 
Orientation parameters
 
I2/I1 I3/I1 |I1·G|V·G|V·I1
M32-06 (E32-06) 781 1.6 10.0 9.8 (9.8) 0.69 (0.68) 0.04 (0.08) 532 (930) 4.27 (5.33) 2.02 (2.97) 4.60 (4.02) 0.85 (0.89) 0.01 (0.06) 0.49 (0.20) 
M28-10 (E28-10) 1246 1.75 11.0 10.8 (10.9) 0.70 (0.73) 0.1 (−0.06) 1500 (1840) 4.35 (5.21) 2.94 (3.44) 6.67 (6.93) 0.97 (0.80) 0.04 (0.09) 0.20 (0.41) 
M25-11 (E25-11) 1168 2.0 9.0 8.8 (8.8) 0.78 (0.79) 0.1 (−0.1) 1781 (2340) 5.87 (8.31) 2.15 (1.90) 5.40 (5.46) 0.99 (0.92) 0.03 (0.12) 0.08 (0.14) 
M25-10 (E25-10) 834 2.4 12.0 12.1 (12.0) 0.89 (0.87) −0.06 (0) 1609 (2057) 5.41 (6.73) 2.30 (2.65) 4.80 (4.98) 0.96 (0.99) 0.01 (0.18) 0.16 (0.18) 
M21-06 (E21-06) 617 3.0 12.0 11.9 (11.6) 1.10 (1.00) 0.04 (0) 2375 (2407) 6.42 (7.23) 2.18 (2.56) 4.41 (4.53) 0.99 (0.96) 0.04 (0.09) 0.01 (0.11) 
M29-03 (E29-03) 448 3.0 10.1 10.2 (10.1) 1.00 (1.09) 0.04 (0) 1218 (2552) 4.66 (6.21) 1.97 (3.58) 4.31 (5.96) 0.98 (0.97) 0.02 (0.27) 0.08 (0.06) 
M25-08 (E25-08) 1360 4.5 11.9 11.9 (11.9) 1.26 (1.25) 0.02 (0.16) 5171 (12 646) 6.00 (11.92) 3.95 (3.32) 5.65 (5.12) 0.98 (0.95) 0.01 (0.14) 0.12 (0.12) 
M17-06 (E17-06) 534 4.0 9.5 9.1 (9.2) 1.28 (1.30) 0.02 (0.5) 3483 (5465) 6.46 (9.12) 2.75 (2.76) 6.10 (6.94) 0.98 (0.91) 0.03 (0.09) 0.08 (0.32) 
M16-05 (E16-05) 2631 4.0 15.0 15.0 (15.0) 1.33 (1.31) 0.04 (0) 17388 (28 128) 13.04 (17.14) 2.36 (2.46) 4.13 (4.07) 0.98 (0.90) 0.01 (0.19) 0.20 (0.25) 
M31-01 (E31-01) 1856 5.4 6.9 7.0 (6.9) 1.54 (1.51) 0.05 (0.17) 11547 (33 487) 6.57 (19.00) 6.09 (2.44) 8.646 (4.07) 0.99 (0.95) 0.05 (0.09) 0.03 (0.13) 
Statistical tests model versus empirical data (Supplementary Table S1 ME1 ME8 ME7 ME2 ME3 ME6 ME4 ME5 
Correlation tests (Supplementary Table S1C7   C5, C6 C3, 4 C1, C2    
Flock/flocking event Number of birds rsep(m) Cruise speed (m/s) Velocity(m/s) NND (m) Balance shift Volume (m3Thickness I1 (m) Aspect ratios
 
Orientation parameters
 
I2/I1 I3/I1 |I1·G|V·G|V·I1
M32-06 (E32-06) 781 1.6 10.0 9.8 (9.8) 0.69 (0.68) 0.04 (0.08) 532 (930) 4.27 (5.33) 2.02 (2.97) 4.60 (4.02) 0.85 (0.89) 0.01 (0.06) 0.49 (0.20) 
M28-10 (E28-10) 1246 1.75 11.0 10.8 (10.9) 0.70 (0.73) 0.1 (−0.06) 1500 (1840) 4.35 (5.21) 2.94 (3.44) 6.67 (6.93) 0.97 (0.80) 0.04 (0.09) 0.20 (0.41) 
M25-11 (E25-11) 1168 2.0 9.0 8.8 (8.8) 0.78 (0.79) 0.1 (−0.1) 1781 (2340) 5.87 (8.31) 2.15 (1.90) 5.40 (5.46) 0.99 (0.92) 0.03 (0.12) 0.08 (0.14) 
M25-10 (E25-10) 834 2.4 12.0 12.1 (12.0) 0.89 (0.87) −0.06 (0) 1609 (2057) 5.41 (6.73) 2.30 (2.65) 4.80 (4.98) 0.96 (0.99) 0.01 (0.18) 0.16 (0.18) 
M21-06 (E21-06) 617 3.0 12.0 11.9 (11.6) 1.10 (1.00) 0.04 (0) 2375 (2407) 6.42 (7.23) 2.18 (2.56) 4.41 (4.53) 0.99 (0.96) 0.04 (0.09) 0.01 (0.11) 
M29-03 (E29-03) 448 3.0 10.1 10.2 (10.1) 1.00 (1.09) 0.04 (0) 1218 (2552) 4.66 (6.21) 1.97 (3.58) 4.31 (5.96) 0.98 (0.97) 0.02 (0.27) 0.08 (0.06) 
M25-08 (E25-08) 1360 4.5 11.9 11.9 (11.9) 1.26 (1.25) 0.02 (0.16) 5171 (12 646) 6.00 (11.92) 3.95 (3.32) 5.65 (5.12) 0.98 (0.95) 0.01 (0.14) 0.12 (0.12) 
M17-06 (E17-06) 534 4.0 9.5 9.1 (9.2) 1.28 (1.30) 0.02 (0.5) 3483 (5465) 6.46 (9.12) 2.75 (2.76) 6.10 (6.94) 0.98 (0.91) 0.03 (0.09) 0.08 (0.32) 
M16-05 (E16-05) 2631 4.0 15.0 15.0 (15.0) 1.33 (1.31) 0.04 (0) 17388 (28 128) 13.04 (17.14) 2.36 (2.46) 4.13 (4.07) 0.98 (0.90) 0.01 (0.19) 0.20 (0.25) 
M31-01 (E31-01) 1856 5.4 6.9 7.0 (6.9) 1.54 (1.51) 0.05 (0.17) 11547 (33 487) 6.57 (19.00) 6.09 (2.44) 8.646 (4.07) 0.99 (0.95) 0.05 (0.09) 0.03 (0.13) 
Statistical tests model versus empirical data (Supplementary Table S1 ME1 ME8 ME7 ME2 ME3 ME6 ME4 ME5 
Correlation tests (Supplementary Table S1C7   C5, C6 C3, 4 C1, C2    
a

Quantitative properties of flocks in the model and (between parentheses) in empirical data of flocks in Rome. Empirical data of flocks are labeled with an “E”, for example, E32-06, and corresponding flocks of our model start with “M”, e.g. M32-06. The empirical observations come from in Table 1 of Ballerini et al. (2008b). rsep: radius of separation; velocity is the speed of the center of mass; NND is the average nearest neighbor distance; Balance shift, positive values indicate higher frontal density, I1, I2, I3: shortest/medium/longest axis of the flock; G, unit vector parallel to gravity; V, unit vector of the normalized velocity. Left to balance shift parameters have been tuned to empirical data.

In order to force individuals to move mostly in a horizontal area above the sleeping site, we adjust the horizontal attraction to the roost (tuning wRoostH). The sleeping site is represented with radius, RRoost (Table 1), of the approximate size of the roost at Termini in Rome (Carere et al. 2009). We tune the vertical attraction of individuals to a specific altitude (calibrating wRoostV) in order to prohibit them from flying off vertically and to produce a “flat” flock shape such as has been observed in reality (Ballerini et al. 2008b). The modeled bird reacts to others at its back only for movements that are really necessary: Following Couzin et al. (2002), individuals avoid others in their blind angle, thus all around (Figure 2), because it is very dangerous to collide frontally with each other. Following models by others (Huth and Wissel 1992), there is a blind angle for the forces of alignment and cohesion because it is more difficult to react to others when the others are located at the back. We adjust the flying behavior by observing the behavior of a single individual in flocks of several sizes. We tune the relaxation parameter of the speed τ in combination with the vertical attraction to the roost wRoostV so that flocks stay at similar altitude when they are not disturbed. The reaction time of individuals is represented by the update time of the “steering” force (Δu seconds, Equation 2). This interval we make longer than the time step of integration of the location and speed (of Δt seconds, see Equations 21 and 22). Due to their “slow” reaction, minor errors in navigation result. We also add a random force scaled to 10% of the steering force. It reflects random errors, for example, in perception and control of movement (Equations 13 and 14).

We have made side-by-side comparisons of the visualization of the model and video recordings in Rome from Claudio Carere. To tune our model, we adjusted parameters of banking, wβin and wβout, so that the roll rates and bank angles match those observed in the videos of real flocks, and controlled the weight of the alignment force, wa, to obtain the strong polarization of real starlings.

At the beginning of each simulation, individuals start with random orientations at random positions inside a cylinder with radius, RRoost.

Measurements and statistical analysis

All variables are measured 10 times per second during 2 min (Table 2). Nearest neighbor distances are averaged per flock. To avoid effects of the flock border (Ballerini et al. 2008b), we consider only individuals in the interior of the flock. In the interior of the flock, individuals are surrounded at all sides by others as reflected by their values of centrality Ci(t) (Equation 7) being smaller than 0.35. The angular distribution of nearest neighbors of all individuals is measured as the distribution of cosines of the angle θ between the forward direction of an individual and the vector to its nearest neighbor. In a circle around the individuals, angles are labeled as follows: A nearest neighbor that is straight ahead is defined as having and angle θ of zero, thus, its cosine is +1. Exactly, orthogonally at either side is defined as being under an angle θ of 90° (thus, its cosine is 0), precisely behind implies an angle θ of 180° (thus a cosine of −1), “etcetera.”

The volume of a flock is measured by mapping the position of the individuals on a cubic lattice and counting the occupied lattice cells, so-called voxelization. For this, the cell size is adapted to the average distance to the nearest neighbors in the flock.

Balance shift is computed as the difference between center of mass and geometrical center divided by the length of the axis along the direction of motion. Positive values indicate higher frontal density.

The flock dimensions are measured by means of the minimum-volume bounding box of the flock. This bounding box is calculated by means of a principal component analysis of the co-ordinates of the flock members (Barequet and Har-Peled 2001) (Movie S3). The eigenvectors that are associated with the largest/medium/smallest eigenvalue of the covariance matrix provide a co-ordinate system oriented along the axes defining the longest/medium/shortest dimensions of the flock, respectively, I3, I2, I1.

Orientation parameters (Ballerini et al. 2008b) are given by 3 angles, the angle between the shortest dimension of the flock I1 and gravity G, between the shortest dimension I1 and the direction of movement of the flock’s center of mass V, and between gravity G and the direction of movement of the flock's center of mass V. These angles are represented by the absolute values of the dot products (which represent the cosines of these angles), |I1·G|, |V·I1| and |V·G|.

Flock banking is calculated every 0.1 s for 8 min. Given that I1 is the “up”-direction of the bounding box (because the flock is flat) and V is the direction of movement of the flock's center of mass, the banking of the flock is the angle between the horizontal plane and the vector orthogonal to the plane spanned by I1 and V (thus I1 × V, where “×” indicates the cross product).

We have used nonparametric statistics (the Wilcoxon matched-pairs signed-ranks test and the Kendall tau correlations) to avoid the difficulty of determining for small sample sizes whether data are distributed according to a normal distribution. For sample sizes above 10, we used the Pearson correlation because it is robust as regards the distribution of data (Sokal and Rohlf 1995). To be conservative, we have only considered 2-tailed probabilities.

Video recordings and other empirical data

The “quantitative” data of stereo photography were taken at the roost in the city center close to the central station “Termini” in Rome (Ballerini et al. 2008b). They were taken of flocks that were not threatened by attacks from predators. Predator attacks by falcons (Falco peregrinus) have been shown to influence the flocking patterns (Carere et al. 2009). For our “qualitative” empirical data, we used video recordings of flocks of the same roost taken between 19 January 2006 and 17 March 2006 and between 12 December 2006 and 2 March 2007 with a HD video camera (JY-HD10, JVC).

RESULTS

Comparison to qualitative empirical data

In our model, starling-like flocks emerge (Movies 1 and 2). They resemble the empirical data qualitatively as regards 1) the shapes of the flocks (Figure 4) and their splitting (Supplementary Material Figure S2), 2) the trajectories of groups and individuals during a turn (Figures 5A and 6), 3) their changes in shape and density over time (Movies 1 and 2 of the model vs. Movie S6 and Figure 6), and 4) the way complete flocks tilt sideward while they are turning (Figure 5B). The movement paths and changes of shape and density of flocks in the model can be compared with data of real starlings only qualitatively because quantitative data of time series do not yet exist. In the model, the degree of banking by individuals is correlated with the curvature κ of their paths (Figure 6A); the turn starts by the banking of the individuals (Figure 6A, I–II), leading to loss of effective lift (Figure 6B, I–II) and altitude and thus, to increase of vertical speed (Figure 6B, I–II). After passing the apex of the turn, banking diminishes (Figure 6A, II–III) and effective lift is gained (Figure 6B, II–III), whereas there is still a loss of altitude, which causes increase in horizontal speed (Figure 6C, II–III). Subsequently, individuals slowly return to level flight (Figure 6A, III–IV), to their original altitude and horizontal speed (Figure 6C, III–IV). While turning, the flock is compressed (Figure 6D), and although its thickness remains similar, its relative proportions change, particularly the largest over the shortest dimension, I3/I1 (Figure 6E). Thus, the variation of flock shape arises in the model mostly during turning, and due to this variability, the shape of the flock is often not oblong (Movie S3).

Figure 4

Resemblance between images of flocks of different sizes and densities from video recordings in Rome (Top row) and in the model (Bottom row). Estimated flock size (from left to right): N = 700, N = 2500–3000, and N = 2000. The outer flocks move straight forward, whereas the flock in the middle is turning. To match the sparse flock at the left, the separation radius in the model was set to rsep = 5. Otherwise, the parameters are kept at their default value.

Figure 4

Resemblance between images of flocks of different sizes and densities from video recordings in Rome (Top row) and in the model (Bottom row). Estimated flock size (from left to right): N = 700, N = 2500–3000, and N = 2000. The outer flocks move straight forward, whereas the flock in the middle is turning. To match the sparse flock at the left, the separation radius in the model was set to rsep = 5. Otherwise, the parameters are kept at their default value.

Figure 5

Properties of flocks turning behavior: (A) right turn of a small flock in the model (N = 50, default parameters). The trajectories of 3 individuals are marked at 3 points in time. (B) Flock banking versus average individual banking in the model (average and standard deviation for flock M32-06). Pearson correlation between tilting of flock and average degree of banking, N = 4800, r = 0.71, P < 0.001.

Figure 5

Properties of flocks turning behavior: (A) right turn of a small flock in the model (N = 50, default parameters). The trajectories of 3 individuals are marked at 3 points in time. (B) Flock banking versus average individual banking in the model (average and standard deviation for flock M32-06). Pearson correlation between tilting of flock and average degree of banking, N = 4800, r = 0.71, P < 0.001.

Figure 6

Time series of flight behavior and flock properties in the model. x axis in seconds. (A) average bank angle and path curvature, κ; (B) Effective lift, dashed line equals lift at cruise speed and at level flight; (C) horizontal speed and altitude (zero is the preferred level above the roost); At the right of (A), (B), and (C) a close-up is shown, in which I, II, III, and IV mark transitions between phases described in the text; (D) Volume and nearest neighbor distance (NND); (E) aspect ratio measured as I3:I1 (largest: shortest extent), I2:I1 (middle: shortest extent) and thickness (I1). (F) Trajectories (spatial top view). Numbers indicate points in time. Flock size N = 2000, default parameters.

Figure 6

Time series of flight behavior and flock properties in the model. x axis in seconds. (A) average bank angle and path curvature, κ; (B) Effective lift, dashed line equals lift at cruise speed and at level flight; (C) horizontal speed and altitude (zero is the preferred level above the roost); At the right of (A), (B), and (C) a close-up is shown, in which I, II, III, and IV mark transitions between phases described in the text; (D) Volume and nearest neighbor distance (NND); (E) aspect ratio measured as I3:I1 (largest: shortest extent), I2:I1 (middle: shortest extent) and thickness (I1). (F) Trajectories (spatial top view). Numbers indicate points in time. Flock size N = 2000, default parameters.

In small flocks in the model, the trajectories of all individuals during turns are similar in length and in radius (Figure 5A, Movie S5): individuals after the turn end up in another location inside the flock. For instance, during the right turn in Figure 5A, the “square” individual is initially located in the left of the flock (at T = 0 s), but after the turn (T = 3 s), it is located at the right side (grey line). Note that in this case, the shape of the flock is unchanged, but its orientation has changed relative to the movement direction because the point of the triangle (connecting 3 individuals) initially points backwards (T = 0 s), but after the turn (3, 6 s), it points forward.

During spontaneous turns of smaller flocks in the model, sometimes the whole flock banks in the direction of individual banking, so-called “flock banking” (e.g., M32-06, Figure 5B). This is apparent from the correlation between the degree of banking of the complete flock and the average degree by its members. It is found only in small flocks and during smooth turns.

Comparison to quantitative empirical data

In order to compare quantitative properties of 10 modeled flocks to the 10 corresponding flock events in Rome (Table 1 of Ballerini et al. 2008b), we first adjusted the following 3 traits to those of each of the 10 empirical flock events: flock size (N), speed (by tuning cruise speed v0), and average distance to the nearest neighbor (by tuning the separation radius rsep). Next, we showed that for each of the 10 empirical flock events, a flock in the model could be found that resembled its empirical counterpart in many respects (Table 2 and Figure 7A,B, for significance tests, see Supplementary Material Table S1).

Figure 7

(A) Distribution of nearest neighbor distance of 2 flocks that differ in their average density (Table 2). Closed circles: empirical data from flock E25-11, solid line: corresponding flock in the model M25-11. Open circles: empirical data from flock E32-06, dotted line: corresponding flock M32-06. (B) Distribution of directions of nearest neighbor (i.e., “bearing angles”). Open circles: average angular distribution of the corresponding empirical data E21-06, E32-06, E25-10, and E25-11 after (Ballerini et al. 2008b). Black rectangles: average angular distribution of the modeled flocks M21-06, M32-06, M25-10, and M25-11. Correlation with empirical distribution: Kendall tau, N = 9, τ = 0.72, P < 0.01. Absolute difference between distribution in model and empirical data: Wilcoxon matched-pairs signed-ranks test, N = 9, sr = 14, not signifcant. Dashed gray curve: averaged nearest neighbor distribution of 100 randomly chosen flocks in the model. Dashed line: expected value for a random isotropic system Prandom(cos(θ)) = 0.5. Error bars: standard error.

Figure 7

(A) Distribution of nearest neighbor distance of 2 flocks that differ in their average density (Table 2). Closed circles: empirical data from flock E25-11, solid line: corresponding flock in the model M25-11. Open circles: empirical data from flock E32-06, dotted line: corresponding flock M32-06. (B) Distribution of directions of nearest neighbor (i.e., “bearing angles”). Open circles: average angular distribution of the corresponding empirical data E21-06, E32-06, E25-10, and E25-11 after (Ballerini et al. 2008b). Black rectangles: average angular distribution of the modeled flocks M21-06, M32-06, M25-10, and M25-11. Correlation with empirical distribution: Kendall tau, N = 9, τ = 0.72, P < 0.01. Absolute difference between distribution in model and empirical data: Wilcoxon matched-pairs signed-ranks test, N = 9, sr = 14, not signifcant. Dashed gray curve: averaged nearest neighbor distribution of 100 randomly chosen flocks in the model. Dashed line: expected value for a random isotropic system Prandom(cos(θ)) = 0.5. Error bars: standard error.

Just as in the empirical data, in the modeled flocks, the balance shift was low, which indicates that the density is similar at front and back of the flock (ME1, Table 2).

The shape of all flocks as measured through the aspect ratios of the shorter dimensions, I2/I1, and the longest divided by the shortest dimension I3/I1 was similar and resembled the empirical values (ME2 and 3, Table 2). Just as in empirical data, I3/I1 did neither depend on flock size nor volume (C1 and 2, Table 2), but in contrast to empirical data, flock shape measured by the aspect ratio I2/I1 did depend on flock size and volume: Larger flocks (whether in number or in volume) are more asymmetric, and smaller ones are closer to being spherical (C3 and 4, Table 2).

As in empirical data, in StarDisplay, the plane of the flock was oriented parallel to the ground as measured by studying to what degree the shortest dimension of the flock I1 was parallel to the gravity G and, therefore, |I1·G| was close to 1. Despite the similarity of the average values in the model |I1·G| = 0.97 ± 0.04 and empirical data |I1·G| = 0.92 ± 0.06, the plane of the modeled flocks was significantly more parallel to the ground. Flocks moved mostly at a constant height above the ground, and therefore, both dot products were on average low, that is, |V·G| of the movement direction and the direction of gravity and |V·I1| of the angle between the velocity and the direction of the shortest dimension. However, the modeled flocks showed a smaller variation in altitude above the ground, that is, the average value of |V·G| = 0.03 ± 0.04 was lower than in empirical data |V·G| = 0.13 ± 0.06, even though the average values of |V·I1| were similar, that is, in StarDisplay |V·I1| = 0.15 ± 0.14 and in the empirical data |V·I1| = 0.19 ± 0.11. Thus, compared with the empirical data, our modeled flocks flew on average more horizontal and more parallel to the ground.

In the model, the thickness and volume are significantly smaller than in the empirical data (ME7 and 8, Table 2).

As in the empirical data, in our model, the nearest neighbor distance (NND) appears to be independent of group size for the 10 flocks (C7, Table 2) and also for 10 replicas of an even larger range of group sizes than the empirical data (viz. for a group size of 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, and 5000 individuals, default radius of separation rsep, Pearson correlation, N = 140, r = 0.32, P = 0.14).

For the following measures of internal structure, 1) the distribution of the distances to the nearest neighbors and 2) of the angles to the nearest neighbors, we compare our modeled flocks to their corresponding empirical flock. Of the empirical flock events E25-11 and E32-06, the frequency distribution of the distances to the nearest neighbors shows a peak around the hard sphere of 0.4 m, and the distribution is skewed to the right as measured by the second derivative of the shape of the curve (skewness for E25-11 and E32-06 are, respectively, +0.36 and +0.39). In the model of the corresponding flocks, the skew is also to the right (skewness for M25-11 and M32-06 are, respectively, +0.37 and +0.54, Figure 7A).

In empirical data of 4 flocks (E32-06, E17-06, E25-10, and E25-11), the frequency distribution of the angle between the heading of the individual and the location of the nearest neighbor (i.e., “bearing angle”) shows that the closest neighbors are lacking along the direction of the movement. In the corresponding 4 flocks of StarDisplay (M32-06, M17-06, M25-10, and M25-11), we find a similar distribution of bearing angles (compare model data in squares to empirical data in circles, Figure 7B).

DISCUSSION

Insights from the model

In line with earlier models of swarming (Couzin et al. 2002; Hemelrijk and Kunz 2005; Parrish and Viscido 2005; Buhl et al. 2006; Hemelrijk and Hildenbrandt 2008), our model StarDisplay shows that the flocking maneuvers of starlings may result from local interactions only. The interactions are local because they take place in a space (of <10 m3) that is much smaller than that of the complete flock (which ranges between 500 and 17.000 m3, Table 2). StarDisplay shows that neither perception of the complete flock nor any leadership or complex cognition is needed. With this result, the model reduces the need for perceptual and cognitive assumptions for the individual birds. As long as there is no evidence to the contrary, a model based on local perception and simple cognition is to be preferred to a model that assumes that the individual interacts with the whole group.

Density is independent of flock size, and it is similar in the front half and back half, both in our model and in empirical data of starlings (Ballerini et al. 2008b). This may be due to the individual's short range of interaction, which involves a low and constant number of (topological) interaction partners (Ballerini et al. 2008a). The greater density that is associated with larger groups of fish (Breder 1954; Keenleyside 1955; Nursall 1973; Partridge 1980; Partridge et al. 1980) is explained in our models by a greater attraction among a higher number of individuals in larger groups (Kunz and Hemelrijk 2003; Hemelrijk and Kunz 2005; Hemelrijk and Hildenbrandt 2008). This difference may indicate that fish interact with a higher number of influential neighbors than starlings do. To study this, modeling experiments with topological range should be done. Further, the independence of density and flock size and absence of high frontal density in starlings may also follow from the frequent turns the flocks take because in turning, groups are compressed, and this compression may neutralize the effect of the group size. To omit the effects of turning in empirical data and in our model StarDisplay, flocks should be studied when they are traveling in a straight line for a longer period (i.e., outside the roost).

In empirical data of starlings by Ballerini et al. (2008b), the aspect ratios of flocks were fixed and independent of flock size. However, in StarDisplay, aspect ratios vary with time and flock size. This maybe due to a bias in the selection of events of real flocks. Investigators may have selected only large flocks (i.e., with more than 447 individuals) that traveled in a straight line. This is also the kind of flocks we had to chose in our model in order to match their empirical data (Table 2). In our model, variability of relative dimensions of flocks is a consequence of abrupt turns, which lead to temporarily compression of the flock (Figure 6). To verify these suggestions, we have to measure in empirical data of starling flocks 1) whether aspect ratios of flock shape change over time and 2) whether turns are the major cause of these shape changes (in the absence of predation).

Through our model, we develop a hypothesis about the interconnection between 3 traits described empirically (Ballerini et al. 2008b): 1) while turning, the shape of the flock changes in relation to the direction of movement, 2) while turning, all individuals of a flock follow paths of the same length with equal radii in real starlings and rock doves (Pomeroy and Heppner 1992), and 3) in starlings, the neighbors that are nearest are rarely located behind or ahead. Our model confirms these 3 traits (Figures 5A, 6A,D,E, and 7B). They may be interconnected as follows. Because individuals follow paths of the same length while turning, this means that they turn without speeding up or slowing down, and therefore, their speed variance is low. Low speed variance in a turn causes the shape of the flock to change in relation to the movement direction; avoiding collisions by turning away rather than by speeding up or slowing down may bring their closest neighbors at their side more often (Figure 7B). This may be useful because collisions at an equal speed from the side may be less harmful than those caused by speeding up and slowing down. These suggested interrelations must be investigated in empirical and theoretical studies.

Limitations of the model

StarDisplay produces lower values than those in the empirical data in 4 aspects, namely, the volume and thickness of the starling flocks, the variation in altitude above the ground, and the deviation from the plane of the flock from being parallel to the ground. All 4 differences may be interrelated: They may be due to the absence of frequent disturbances in the model, such as predators, wind, rain, and change in lighting conditions. The relevance for real starlings of this hypothesis may be investigated by studying in real starlings whether attacks of a predator cause flocks to increase in volume on average (due to continuous reshuffling). This is a counter-intuitive hypothesis because on the short term, a predator attack may cause compression of the flock, for example, to a ball shape.

Further, in our model, robustness against parameter changes must be studied in greater depth. Because we were interested in modeling the flocking maneuvers of starlings in Rome, we have not systematically studied different parameters, apart from a range of empirically measured lift/drag ratios (Möller 2005). For these values of ratios of lift/drag, the results of the model showed no change.

A general often heard criticism of models that are based on many parameters is that they are lacking in explanatory value because they too easily fit single patterns. However, most of our parameters are “fixed” because they depend on data of starlings. Further, we have generated as emergent phenomena not merely a single pattern, but many of them. We have followed the so-called bottom-up procedure in which we trace the connection between simple behavioral rules and collective patterns (Hemelrijk 2002). Further, StarDisplay may be used as a “prediction” for empirical data. This concerns its flock dynamics over time because despite qualitative resemblance between the patterns of the model and those of flocks of rock doves (Pomeroy and Heppner 1992) and starlings (Figure 6), quantitative empirical data are still lacking.

CONCLUSIONS

In sum, we demonstrate in our model StarDisplay that local interactions and self-organization suffice to reproduce patterns of aerial display of starlings qualitatively and in many cases also quantitatively. In future, we will study the causes of these patterns in the model. Our model-based hypotheses may be useful in indicating suitable topics for empirical research of collective flocking maneuvers not only of starlings but also of other birds. This may be particularly useful because the empirical study of huge swarms is labor intensive.

FUNDING

The 6th European framework under the STREP-project “StarFlag” (n12682) in the NEST-programme of “Tackling complexity in science.”

SUPPLEMENTARY MATERIAL

Supplementary material can be found at http://www.beheco.oxfordjournals.org/.

We would like to thank Andrea Cavagna, Volker Grimm, and Daan Reid for comments on a former draft. We are grateful to Dirk Visser for drawing Figures 2 and 5a.

References

Ame
JM
Halloy
J
Rivault
C
Detrain
C
Deneubourg
JL
Collegial decision making based on social amplification leads to optimal group formation
Proc Natl Acad Sci U S Am
 , 
2006
, vol. 
103
 (pg. 
5835
-
5840
)
Ballerini
M
Cabibbo
N
Candelier
R
Cavagna
A
Cisbani
E
Giardina
I
Lecomte
V
Orlandi
A
Parisi
G
Procaccini
A
, et al.  . 
Interaction ruling animal collective behaviour depends on topological rather than metric distance: evidence from a field study
Proc Natl Acad Sci
 , 
2008
, vol. 
105
 (pg. 
1232
-
1237
)
Ballerini
M
Cabibbo
N
Candelier
R
Cavagna
A
Cisbani
E
Giardina
I
Orlandi
A
Parisi
G
Procaccini
A
Viale
M
, et al.  . 
Empirical investigation of starling flocks: a benchmark study in collective animal behaviour
Anim Behav
 , 
2008
, vol. 
76
 (pg. 
201
-
215
)
Barequet
G
Har-Peled
S
Efficiently approximating the minimum-volume bounding box of a point set in three dimensions
J Algorithms
 , 
2001
, vol. 
38
 (pg. 
91
-
109
)
Beauchamp
G
The evolution of communal roosting in birds: origin and secondary losses
Behav Ecol
 , 
1999
, vol. 
10
 (pg. 
675
-
687
)
Becco
C
Vandewalle
N
Delcourt
J
Poncin
P
Experimental evidences of a structural and dynamical transition in fish school
Physica a-Stat Mech Its Appl
 , 
2006
, vol. 
367
 (pg. 
487
-
493
)
Biro
D
Sumpter
DJT
Meade
J
Guilford
T
From compromise to leadership in pigeon homing
Curr Biol
 , 
2006
, vol. 
16
 (pg. 
2123
-
2128
)
Breder
CM
Equations descriptive of fish schools and other animal aggregations
Ecology
 , 
1954
, vol. 
35
 (pg. 
361
-
370
)
Brodie
J
The flight behaviour of starlings at a winter roost
Br Birds
 , 
1976
, vol. 
69
 (pg. 
51
-
60
)
Budgey
R
The three dimensional structure of bird flocks and its implications for birdstrike tolerance in aircraft
Int Bird Strike Comm Proc
 , 
1998
, vol. 
24
 (pg. 
207
-
220
)
Buhl
J
Sumpter
DJT
Couzin
ID
Hale
JJ
Despland
E
Miller
ER
Simpson
SJ
From disorder to order in marching locusts
Science
 , 
2006
, vol. 
312
 (pg. 
1402
-
1406
)
Bumann
D
Krause
J
Rubenstein
D
Mortality risk of spatial positions in animal groups: the danger of being in the front
Behaviour
 , 
1997
, vol. 
134
 (pg. 
1063
-
1076
)
Camazine
S
Deneubourg
J-L
Franks
NR
Sneyd
J
Theraulaz
G
Bonabeau
E
Self-organization in biological systems
 , 
2001
Princeton and Oxford
Princeton University Press
Carere
C
Montanino
S
Moreschini
F
Zoratto
F
Chiarotti
F
Santucci
D
Alleva
E
Aerial flocking patterns of wintering starlings, Sturnus vulgaris, under different predation risk
Anim Behav
 , 
2009
, vol. 
77
 (pg. 
101
-
107
)
Cavagna
A
Giardina
I
Orlandi
A
Parisi
G
Procaccini
A
The STARFLAG handbook on collective animal behaviour: 2
Three-dimensional anal Anim Behav
 , 
2008
, vol. 
76
 (pg. 
237
-
248
)
Cavagna
A
Giardina
I
Orlandi
A
Parisi
G
Procaccini
A
Viale
M
Zdravkovic
V
The STARFLAG handbook on collective animal behaviour: 1. Empirical methods
Anim Behav
 , 
2008
, vol. 
76
 (pg. 
217
-
236
)
Couzin
ID
Krause
J
Franks
NR
Levin
SA
Effective leadership and decision-making in animal groups on the move
Nature
 , 
2005
, vol. 
433
 (pg. 
513
-
516
)
Couzin
ID
Krause
J
James
R
Ruxton
GD
Franks
NR
Collective memory and spatial sorting in animal groups
J Theor Biol
 , 
2002
, vol. 
218
 (pg. 
1
-
11
)
Dussutour
A
Fourcassié
V
Helbing
D
Deneubourg
JL
Optimal traffic organization in ants under crowded conditions
Nature
 , 
2004
, vol. 
428
 (pg. 
70
-
73
)
Feare
CJ
The starling
 , 
1984
Oxford
Oxford University Press
Gautrais
J
Theraulaz
G
Deneubourg
J-L
Anderson
C
Emergent polyethism as a consequence of increased colony size in insect societies
J Theor Biol
 , 
2003
, vol. 
215
 (pg. 
363
-
373
)
Helbing
D
Buzna
L
Johansson
A
Werner
T
Self-organised pedestrian crowd dynamics: experiments, simulations, and design solutions
Transp Sci
 , 
2005
, vol. 
39
 (pg. 
1
-
24
)
Hemelrijk
CK
Understanding social behaviour with the help of complexity science
Ethology
 , 
2002
, vol. 
108
 (pg. 
655
-
671
)
Hemelrijk
CK
Self-organisation and evolution of social systems
 , 
2005
Cambridge, UK
Cambridge University Press
Hemelrijk
CK
Hildenbrandt
H
Self-organized shape and frontal density of fish schools
Ethology
 , 
2008
, vol. 
114
 (pg. 
245
-
254
)
Hemelrijk
CK
Hildenbrandt
H
Reinders
J
Stamhuis
EJ
Emergence of oblong school shape: models and empirical data of fish
Ethology 116: doi: 10.1111/j.1439-0310.2010.01818.x
 , 
2010
Hemelrijk
CK
Kunz
H
Density distribution and size sorting in fish schools: an individual-based model
Behav Ecol
 , 
2005
, vol. 
16
 (pg. 
178
-
187
)
Hemelrijk
CK
Wantia
J
Individual variation by self-organisation: a model
Neurosci Biobehav Rev
 , 
2005
, vol. 
29
 (pg. 
125
-
136
)
Hemelrijk
CK
Wantia
J
Isler
K
Female dominance over males in primates: self-organisation and sexual dimorphism
PLoS ONE
 , 
2008
, vol. 
3
 pg. 
e2678
  
doi: 10.1371
Huth
A
Wissel
C
The simulation of the movement of fish schools
J Theor Biol
 , 
1992
, vol. 
156
 (pg. 
365
-
385
)
Huth
A
Wissel
C
The simulation of fish schools in comparison with experimental data
Ecol Modell
 , 
1994
, vol. 
75/76
 (pg. 
135
-
145
)
Keenleyside
MHA
Some aspects of schooling behaviour in fish
Behaviour
 , 
1955
, vol. 
8
 (pg. 
183
-
248
)
Kunz
H
Hemelrijk
CK
Artificial fish schools: collective effects of school size, body size, and body form
Artif Life
 , 
2003
, vol. 
9
 (pg. 
237
-
253
)
Major
PF
Dill
LM
3-Dimensional structure of airborne bird flocks
Behav Ecol Sociobiol
 , 
1978
, vol. 
4
 (pg. 
111
-
122
)
Möller
U
Schlagflug des Stars (Sturnus vulgaris) im Windkanal mit und ohne respiratorische Maske: Kenematik, Aerodynamic und Energetik
 , 
2005
Saarbrücken, Germany
Universität des Saarlandes
Norberg
UM
Vertebrate Flight: mechanics, physiology, morphology, ecology and evolution
 , 
1990
New York
Springer Verlag
Nursall
JR
Some behavioral interactions of spottail shiners (Notropis hudsonius), yellow perch (Perca flavescens), and northern pike (Esox-Lucius)
J Fisheries Res Board Canada
 , 
1973
, vol. 
30
 (pg. 
1161
-
1178
)
Parrish
JK
Edelstein-Keshet
L
Complexity, pattern, and evolutionary trade-offs in animal aggregation
Science
 , 
1999
, vol. 
284
 (pg. 
99
-
101
)
Parrish
JK
Viscido
SV
Hemelrijk
CK
Traffic rules of fish schools: a review of agent-based approaches
Self-Organisation and the Evolution of Social Behaviour
 , 
2005
Cambridge, UK
Cambridge University Press
Partridge
BL
Effect of school size on the structure and dynamics of minnow schools
Anim Behav
 , 
1980
, vol. 
28
 (pg. 
68
-
77
)
Partridge
BL
Pitcher
T
Cullen
JM
Wilson
J
The 3-dimensional structure of fish schools
Behav Ecol Sociobiol
 , 
1980
, vol. 
6
 (pg. 
277
-
288
)
Pomeroy
H
Heppner
F
Laboratory determination of startle reaction time of the starling (Sturnus vulgaris)
Anim Behav
 , 
1977
, vol. 
25
 (pg. 
720
-
725
)
Pomeroy
H
Heppner
F
Structure of turning in airborne Rock Dove (Columba livia) flocks
Auk
 , 
1992
, vol. 
109
 (pg. 
256
-
267
)
Potts
WK
The chorus-line hypothesis of manoeuvre coordination in avian flocks
Nature
 , 
1984
, vol. 
309
 (pg. 
344
-
345
)
Puga-Gonzalez
I
Hildenbrandt
H
Hemelrijk
CK
Emergent patterns of social affiliation in primates, a model
Plos Comput Biol
 , 
2009
, vol. 
5
 pg. 
e1000630
  
doi:1000610.1001371/journal.pcbi.1000630
Reuter
H
Breckling
B
Selforganization of fish schools—an object-oriented model
Ecol Modell
 , 
1994
, vol. 
75
 (pg. 
147
-
159
)
Reynolds
CW
Flocks, herds and schools: a distributed behavioral model
Comput Graph
 , 
1987
, vol. 
21
 (pg. 
25
-
36
)
Reynolds
RG
Flannery
KV
An adaptive computer model for the evolution of plant collecting and early agriculture in the Eastern Valley of Oaxaca, Mexico
Quila Naquitz: archaic foraging and early agriculture in Oaxaca, Mexico San Diego, California
 , 
1986
San Diego
Academic Press
Romey
WL
Individual differences make a difference in the trajectories of simulated schools of fish
Ecol Modell
 , 
1996
, vol. 
92
 (pg. 
65
-
77
)
Selous
E
Thought transference (or what?) in birds
 , 
1931
London
Constanble and Company Ltd
Sokal
RR
Rohlf
FJ
Biometry: the principles and practice of statistics in biological research
 , 
1995
New York
Freeman and company
te Boekhorst
IJA
Hogeweg
P
Brooks
RA
Maes
P
Effects of tree size on travelband formation in orang-utans: data-analysis suggested by a model
Artificial life
 , 
1994
Cambridge, MA
The MIT Press
(pg. 
119
-
129
)
te Boekhorst
IJA
Hogeweg
P
Selfstructuring in artificial ‘CHIMPS' offers new hypotheses for male grouping in chimpanzees
Behaviour
 , 
1994
, vol. 
130
 (pg. 
229
-
252
)
Videler
JJ
Avian flight
 , 
2005
Oxford
Oxford University Press
Ward
AJW
Sumpter
DJT
Couzin
ID
Hart
PJB
Krause
J
Quorum decision-making facilitates information transfer in fish schools
Proc Natl Acad Sci
 , 
2008
, vol. 
105
 (pg. 
6948
-
6953
)
Ward
S
Mueller
U
Rayner
JMV
Jackson
DM
Nachtigall
W
Speakman
JR
Metabolic power of European starlings Sturnus vulgaris during flight in a wind tunnel, estimated from heat transfer modelling, doubly labelled water and mask respirometry
J Exp Biol
 , 
2004
, vol. 
207
 (pg. 
4291
-
4298
)
Zheng
M
Kashimori
Y
Hoshino
O
Fujita
K
Kambara
T
Behavior pattern (innate action) of individuals in fish schools generating efficient collective evasion from predation
J Theor Biol
 , 
2005
, vol. 
235
 (pg. 
153
-
167
)

Supplementary data