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

Parameter | Description | Default value |

Δt | Integration time step | 5 ms |

Δu | Reaction time | 50 ms |

v_{0} | Cruise speed | 10 m/s |

M | Mass | 0.08 kg |

C_{L}/C_{D} | Lift-drag coefficient | 3.3 |

L_{o} | Default lift | 0.78 N |

D_{0}, T_{0} | Default drag, default thrust | 0.24 N |

w_{βin} | Banking control | 10 |

w_{βout} | Banking control | 1 |

T | Speed control | 1 s |

R_{max} | Maximum perception radius | 100 m |

n_{c} | Topological range | 6.5 |

S | Interpolation factor | 0.1 Δu |

r_{h} | Radius of maximum separation (hard sphere) | 0.2 m |

r_{sep} | Separation radius | 4 m |

Σ | Parameter of the Gaussian g(x) | 1.37 m |

w_{s} | Weighting factor separation force | 1 N |

Rear “blind angle” cohesion and alignment | 90° | |

w_{a} | Weighting factor alignment force | 0.5 N |

w_{c} | Weighting factor cohesion force | 1 N |

C_{c} | 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 |

R_{Roost} | Radius of roost | 150 m |

w_{RoostH} | Weighting factor horizontal attraction to the roost | 0.01 N/m |

w_{RoostV} | 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 |

v_{0} | Cruise speed | 10 m/s |

M | Mass | 0.08 kg |

C_{L}/C_{D} | Lift-drag coefficient | 3.3 |

L_{o} | Default lift | 0.78 N |

D_{0}, T_{0} | Default drag, default thrust | 0.24 N |

w_{βin} | Banking control | 10 |

w_{βout} | Banking control | 1 |

T | Speed control | 1 s |

R_{max} | Maximum perception radius | 100 m |

n_{c} | Topological range | 6.5 |

S | Interpolation factor | 0.1 Δu |

r_{h} | Radius of maximum separation (hard sphere) | 0.2 m |

r_{sep} | Separation radius | 4 m |

Σ | Parameter of the Gaussian g(x) | 1.37 m |

w_{s} | Weighting factor separation force | 1 N |

Rear “blind angle” cohesion and alignment | 90° | |

w_{a} | Weighting factor alignment force | 0.5 N |

w_{c} | Weighting factor cohesion force | 1 N |

C_{c} | 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 |

R_{Roost} | Radius of roost | 150 m |

w_{RoostH} | Weighting factor horizontal attraction to the roost | 0.01 N/m |

w_{RoostV} | Weighting factor vertical attraction to the roost | 0.2 N |

Note that *D*_{0} and *T*_{0} is calculated by Equation 15bc by inserting *v*_{0} for *v*_{i}.

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

**,**

*e*_{x}**, and**

*e*_{y}**). Following the model by Reynolds (1987), its orientation is indicated by its forward direction,**

*e*_{z}**, its sideward direction,**

*e*_{x}**, and its upward direction,**

*e*_{y}**, which it changes by rotating around these 3 principal axes,**

*e*_{z}**,**

*e*_{x}**and**

*e*_{y}_{,}**(roll, pitch, and yaw) (Figure 1).**

*e*_{z}As to its speed, a force, $f\tau i$ (Equation 1), brings an individual back to its cruise speed *v*_{0} after it has deviated from it (Hemelrijk and Hildenbrandt 2008).

*τ*represents the relaxation time,

*m*is the mass of the individual

*i*and

*v*

_{0}its cruise speed,

*v*

_{i}is its velocity, and $exi$ its forward direction.

To represent topological interaction, each individual *i* in the model adapts its metric interaction range, *R _{i}*(

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

*u*is the time span to react,

*s*is an interpolation factor,

*R*

_{max}is the maximal metric interaction range,

*N*(

_{i}*t*) is the neighborhood of individual

*i*at time

*t*, that is, the set of neighbors of an individual

*i*which is composed of |

*N*(

_{i}*t*)| neighbors from the total flock of size

*N*,

*n*

_{c}is the fixed number of topological interaction partners it strives to have, and

*d*is the distance between individual

_{ij}*i*and

*j*given by ‖

*p**−*

_{j}

*p**‖, where*

_{i}

*p**gives the position of an individual*

_{i}*i*. Thus, the radius of interaction at the next step in reaction time,

*R*(

_{i}*t*+

*Δu*), is increased if the number of interaction partners |

*N*(

_{i}*t*)| is smaller than the targeted number

*n*

_{c}, and it is decreased if it is larger than that; it remains as before if |

*N*(

_{i}*t*)| equals

*n*

_{c}. Here,

*R*can neither decrease below the hard sphere in which individuals are maximally avoiding each other

_{i}*r*

_{h}(Equation 4, Figure 2) nor increase beyond

*R*

_{max}and

*s*, the interpolation factor, which determines the step size of the changes and herewith, the variance of the number of actual influential neighbors.

As to separation, individual *i* is led by a force $fsi$to move in the opposite direction of the average direction of the locations of the *N _{i}*(

*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

*r*

_{h}as mentioned above in which they are not attracted to each other (Equation 5). Outside the hard sphere, but inside the radius of separation

*r*

_{sep}, the degree of avoidance of others decreases with the distance to the neighbor following a halved Gaussian,

*g*(

*d*), with

_{ij}*σ*the standard deviation of the Gaussian set such that at the border of the separation zone the force is almost zero,

*g*(

*r*

_{sep}) = 0.01 (Equation 4).

*N*(

_{i}*t*)| is the number of individuals in the neighborhood of interaction (Equation 3) and

*d*is the distance from individual

_{ij}*i*to individual

*j*. The direction from individual

*i*to individual

*j*is specified by the unit vector

**= (**

*d*_{ij}**−**

*p*_{j}**)/‖**

*p*_{i}**−**

*p*_{j}**‖ and**

*p*_{i}*w*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.

_{s}As to cohesion (Equations 5, 6 and 7), individual *i* is attracted by a force $fci$to center of mass (i.e., the average *x*, *y*, *z* position) of the group of *N*_{i}^{*} (*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, *w*_{c} is the weighing factor for cohesion (Equation 6, Table 1). Within the radius of the hard sphere *r*_{h}, 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, *C _{i}*, calculated as the length of the average vector of the direction toward its neighbors

*N*(

_{G}*t*) (Hemelrijk and Wantia 2005). The centrality of the individual

*C*(

_{i}*t*), we base on all its “neighboring” individuals

*N*(

_{G}*t*) in a radius of twice its actual perceptual distance (Equation 7).

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

Here, $exi$ and $exj$ are the vectors indicating the forward direction of individuals *i* and *j* and *w*_{a} is the fixed weighting factor for alignment (Table 1).

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

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.

*w*

_{RoostH}and

*w*

_{RoostV}are fixed weighting factors.

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

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):

*S*represents the wing area of the bird.

*C*

_{L}and

*C*

_{D}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

*C*

_{L}and

*C*

_{D}is constant (Norberg 1990). When a bird is flying horizontally while maintaining a constant cruise speed

*v*

_{0}, its lift balances its weight

*mg*(i.e., mass × gravity) and its thrust balances its drag (Equation 15b). Division of

*L*by

*L*

_{0}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

**of the bird, and the drag is pointing in the direction opposite to the individual's actual “forward” direction**

*e*_{z}**. Thus, the flight forces are**

*e*_{x}To represent banked turns, we first calculate the degree to which individuals want to turn, that is, their lateral acceleration, ** a_{l}**, 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,

*a*. The lateral acceleration follows the first law of Newton (

_{l}**= m**

*F***)**

*a**is the actual banking angle,*

_{i}*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 ‖

**‖ sin(β**

*a*_{l}_{in}) parallel to the lift (dotted lines in the last column, Figure 3).

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*:

*v*

_{i}is the speed of individual

*i*,

*m*its mass,

*p**its location, and Δ*

_{i}*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, *C*_{L}/*C*_{D}, 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 |*N _{i}*(

*t*)| of starlings of

*n*

_{c}= 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 (*L*_{0} and *D*_{0}, Table 1). The radius of the hard sphere, *r*_{h}, 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, *w*_{s}, such that no collisions occur (Potts 1984). We adjust the separation radius, *r*_{sep}, to obtain an average nearest neighbor distance similar to that of the starlings of Rome (Table 2, Supplementary Material Figure S1).

Flock/flocking event | Number of birds | r_{sep}(m) | Cruise speed (m/s) | Velocity(m/s) | NND (m) | Balance shift | Volume (m^{3}) | Thickness I_{1} (m) | Aspect ratios | Orientation parameters | |||

I_{2}/I_{1} | I_{3}/I_{1} | |I_{1}·G| | |V·G| | |V·I_{1}| | |||||||||

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 S1) | C7 | C5, C6 | C3, 4 | C1, C2 |

Flock/flocking event | Number of birds | r_{sep}(m) | Cruise speed (m/s) | Velocity(m/s) | NND (m) | Balance shift | Volume (m^{3}) | Thickness I_{1} (m) | Aspect ratios | Orientation parameters | |||

I_{2}/I_{1} | I_{3}/I_{1} | |I_{1}·G| | |V·G| | |V·I_{1}| | |||||||||

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 S1) | C7 | C5, C6 | C3, 4 | C1, C2 |

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). *r*_{sep}: 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, *I*_{1}, *I*_{2}, *I*_{3}: 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 *w*_{RoostH}). The sleeping site is represented with radius, *R*_{Roost} (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 *w*_{RoostV}) 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 *w*_{RoostV} 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, *w*_{a}, 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, *R*_{Roost}.

### 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 *C _{i}*(

*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, *I*_{3}, *I*_{2}, *I*_{1}.

Orientation parameters (Ballerini et al. 2008b) are given by 3 angles, the angle between the shortest dimension of the flock *I*_{1} and gravity *G*, between the shortest dimension *I*_{1} 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), |*I*_{1}*·G*|, |*V·I*_{1}| and |*V·G*|.

Flock banking is calculated every 0.1 s for 8 min. Given that *I*_{1} 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 *I*_{1} and *V* (thus *I*_{1} × *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, *I*_{3}/*I*_{1} (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).

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 *v _{0}*), and average distance to the nearest neighbor (by tuning the separation radius

*r*

_{sep}). 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).

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, *I*_{2}/*I*_{1}, and the longest divided by the shortest dimension *I*_{3}/*I*_{1} was similar and resembled the empirical values (ME2 and 3, Table 2). Just as in empirical data, *I*_{3}/*I*_{1} 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 *I*_{2}/*I*_{1} 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 *I*_{1} was parallel to the gravity G and, therefore, **|***I*_{1}*·G***|** was close to 1. Despite the similarity of the average values in the model **|***I*_{1}*·G***|** = 0.97 ± 0.04 and empirical data **|***I*_{1}*·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·I*_{1}| 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·I*_{1}| were similar, that is, in StarDisplay |*V·I*_{1}| = 0.15 ± 0.14 and in the empirical data |*V·I*_{1}| = 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 *r*_{sep}, 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 m^{3}) that is much smaller than that of the complete flock (which ranges between 500 and 17.000 m^{3}, 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

*Sturnus vulgaris*, under different predation risk

*Sturnus vulgaris*) im Windkanal mit und ohne respiratorische Maske: Kenematik, Aerodynamic und Energetik

*Notropis hudsonius*), yellow perch (

*Perca flavescens*), and northern pike (

*Esox-Lucius*)

*Columba livia*) flocks

*Sturnus vulgaris*during flight in a wind tunnel, estimated from heat transfer modelling, doubly labelled water and mask respirometry