# Flock2: A model for orientation-based social flocking

Rama Carl Hoetzlein
Quanta Sciences
Ithaca, NY 14850
<EMAIL_ADDRESS>

(April 26, 2024)

###### Abstract

The aerial flocking of birds, or murmurations, have fascinated observers while
presenting many challenges to behavioral study and simulation. We examine how
the periphery of murmurations remain well bounded and cohesive. We also
investigate agitation waves, which occur when a flock is disturbed, developing
a plausible model for how they might emerge spontaneously. To understand these
behaviors a new model is presented for orientation-based social flocking.
Previous methods model inter-bird dynamics by considering the neighborhood
around each bird, and introducing forces for avoidance, alignment, and
cohesion as three dimensional vectors that alter acceleration. Our method
introduces orientation-based social flocking that treats social influences
from neighbors more realistically as a desire to turn, indirectly controlling
the heading in an aerodynamic model. While our model can be applied to any
flocking social bird we simulate flocks of starlings, Sturnus vulgaris, and
demonstrate the possibility of orientation waves in the absence of predators.
Our model exhibits spherical and ovoidal flock shapes matching observation.
Comparisons of our model to Reynolds’ on energy consumption and frequency
analysis demonstrates more realistic motions, significantly less energy use in
turning, and a plausible mechanism for emergent orientation waves.

_Keywords_ flocking $\cdot$ murmuration $\cdot$ agitation waves $\cdot$
collective escape $\cdot$ bird flight $\cdot$ aerodynamics $\cdot$ boids
$\cdot$ modeling $\cdot$ computational biology

##  1 Introduction

Collective group formation is well known as many animals aggregate into dense
groups for protection. Fish form schools as a defense (Nikolsky, 1955)
(Breder, 1959) (Pitcher and Wyche, 1983) (Godin and Morgan, 1985) (Magurran,
1990) (Parrish, 2000) and grazing mammals such as oxen and sheep exhibit
strong inclinations to aggregate at all times for protection (Hamilton, 1971)
(Treves, 2000). The concept of the selfish herd, introduced by Hamilton,
provides one theoretical basis: many species aggregate out of self interest as
a form of protection (Hamilton, 1971). This attractive force to seek cover is
balanced by the competing need to maintain some distance between individuals
(Emlen, 1952).

Flocking is a form of collective animal behavior in birds that is strongly
influenced by the predator-prey relationship. Empirical studies have shown
that flocks are larger, flocking for longer times, and with more activity and
vigilance in monitoring when predators are present (Caraco et al., 1980)
(Carbone et al., 2003) (Procaccini et al., 2011) (Goodenough et al., 2017).
Birds such as quelea, juncos and starlings will defend from external predators
through aggregation into murmurations (Vine, 1971) (Caraco et al., 1980).
Other studies have noted that time budgets, the way a bird divides its time
between feeding, monitoring and flocking, are correlated with predators
(McNamara and Houston, 1992) (Sirot, 2006).

Aerial flocking is generally believed to be an airborne extension of predator
defense in combination with activities for foraging and roosting (Siegried and
Underhill, 1975) (Feare, 1984) (Cresswell, 1994) (Beauchamp, 2004) (Fischl and
Caccamise, 1987) (Tinbergen, 1981) (Caccamise, 1990a). Specific aerial
patterns such as escape waves, whereby oscillations are induced in a flock at
the point of attack, have been observed in relation to predators (Procaccini
et al., 2011). The geometric shape of an aerial flock is influenced by many
factors including season, temperature, food resources and predator presence
(Carere et al., 2009).

Computational modeling has become an important tool in the understanding of
aerial flocking. The work of Couzin et al. (2002) on group formation for fish
and birds in three dimensions is based in part on an earlier model by Reynolds
(1987). Reynolds defines avoidance as the tendency for a bird to avoid
collisions with neighbors, alignment as the tendency for a bird to match
direction and velocity with neighbors, and cohesion as the tendency to move
toward the local average position of neighbors. This approach, relying on
vector forces, continues to be the primary method employed in subsequent
models of aerial flocking dynamics (Couzin et al., 2002) (Couzin et al., 2005)
(Parrish and Viscido, 2005) (Buhl et al., 2006) (Hemelrijk and Hildenbrandt,
2008) (Hildenbrandt et al., 2010).

Despite the extensive use of Reynold’s model, there are discrepancies between
observational and computational biology that our work seeks to address. Some
models introduce a preferred direction (Couzin et al., 2005) or roosting site
(Hildenbrandt et al., 2010) to overcome the fact that Reynolds three basic
terms do not generally result in globular or ovoidal flocks, despite the fact
that these are among the most common observed real flock shapes (Carere et
al., 2009). This is generally unsatisfying as a continuous in-flight
adjustment since roosting is a resting behavior (Caccamise, 1990b).

Figure 1: Previous methods implement social factors as 3D vector forces in a)
Reynolds (1987) and b) with aerodynamic forces in Hildenbrandt et al. (2010).
We develop a flocking model c) in which bird motion is strictly controlled by
an aerodynamic model for which the social behaviors of avoidance, alignment
and cohesion are combined to act as a single actuator for a desired change in
orientation.

Hildenbrandt et al. (2010) model aerial flocking by combining the social
forces of Reynold’s model with aerodynamic forces of flight, even though the
only group of physical forces considered in the study of real bird flight are
aerodynamic since the Reynold’s social forces are not actual forces of
mechanical work (Pennycuick, 1968) (Withers, 1981). The social influences are
perceptual, which presumably influence the bird to change its wing shape,
direction and speed.

Agitation waves, whereby a predator causes the flock to deform or escape, have
been observed extensively in starlings (Carere et al., 2009) (Procaccini et
al., 2011) (Storms, 2023). Hemelrijk et al. (2015) develop a computational
model for agitation waves, yet this required the explicit coding of a
"transmission mechanism" for escape messages from one bird to the next, and
whereby a predator is strictly assumed to trigger the first message.

These discrepancies between observational and computational biology led us to
consider a new foundational model in which social forces are observational
rather than explicit 3D vectors doing mechanical work. This approach is more
realistic since a bird may sense its neighbors but is ultimately only able to
alter its wing shape, heading and speed.

##  2 Model and simulation

This work is motivated by the following observations:

  * •

Real birds are known to turn to a greater degree than altering their speed for
energy savings. "Speed is a stiffer mode than orientation, as it is more
costly for a bird to change its speed than its heading.” (Cavagna et al.,
2010)

  * •

Bird physiology provides that sensory information, including the vision of
neighbors, is integrated to control flight. (Altshuler et al., 2015)

  * •

Social forces in Reynold’s model are not actual forces doing mechanical work,
they are better thought of as perceptual factors that influence heading.

  * •

Computational models for the aerodynamics of birds must allow for some
articulation or heading control, since this is the only way a bird can alter
its desired direction.

These observations led to the consideration of an orientation-based social
model which controls an actuated aerodynamic model. An important aspect of
this design is that only the total aerodynamic forces are permitted to update
the bird velocity and position, whereas the social orientation-model is only
permitted to influence the desired heading.

###  2.1 Bird Flight

Understanding the flight of real birds should ideally provide the basis for
flocking models. As we introduce an aerodynamic model coupled to a flocking
model, the study of real bird flight is essential in model design and
parameter analysis. Real bird flight has been analyzed either as a flapping
wing (Pennycuick, 1968) (Rayner, 1993) (Rayner, 1999) or as a fixed wing
(Withers, 1981) (Norberg, 1985).

A distinction is made between mechanical power for the total external forces
of motion (Pennycuick, 1968) (Dial et al., 1997), and metabolic power which is
the total bodily energy produced including basal metabolism, respiration,
internal resistance and transformations of chemical energy (Daan et al., 1990)
(Rayner et al., 2015) (Ward et al., 2001) (Ward et al., 2004) (Tobalske,
2007).

We consider only mechanical power in our approach whereas a more complete
model would simulate both the internal and external functions of a bird. The
theory of flapping wings suggests that mechanical power output follow a
U-shaped curve as real birds use the dynamic stall effect to achieve better
lift at low speeds and with drag forces dominating at higher speeds
(Pennycuick, 1968). Our analaysis examines simulated birds in light of this
theory.

The starling is used as reference for our models with an average mass of 80
grams and typical mechanical power output of 10 watts flying at an average
speed of 10 meters/second (Withers, 1981) (Dial et al., 1997) (Nudds and
Bryant, 2000) (Schmidt-Wellenburg et al., 2008). The predominant energy use
corresponds to the lift force applied downward to move a volume of air with a
reactionary upward force on the body. The next most energetic expenditure is
to overcome drag including profile drag on the wings, with turbulence, and
parasitic drag on the body. Studies show that energy usage increases
proportionally with bird weight (Schmidt-Wellenburg et al., 2008). The
remaining mechanical energy is applied to accelerate, decelerate and turn. We
rely on this knowledge for both model design and analysis.

###  2.2 Mathematical background

In Reynold’s model the social terms for avoidance, alignment and cohesion are
summed to give the total force acting on the motion of the bird:

$F_{social}=F_{avoidance}+F_{alignment}+F_{cohesion}$ (1)

These social forces act directly on velocity as if they were doing mechanical
work with velocity updating position (Reynolds, 1987).

Aerodynamics have also been incorporated into social flocking in birds by
Hildenbrandt et al. (2010). In their approach, the Reynold’s social forces are
combined in the same way using vectors to give $F_{social}$. With a separate
term, they model the bird as a fixed wing where lift is the force resulting
from differential pressure flowing over the wing, and drag is the fluid
resistance of the air against forward motion. These are summed with gravity to
give $F_{flight}$, the total aerodynamic forces as a 3D vector. The social
forces are scaled by a constant factor $k$ so their magnitude is similar to
the flight forces, and added together to give the final bird motion:

$F_{total}=F_{flight}+kF_{social}$ (2)

Their work was the first to model a distinction between social forces, which
are based on bird behavior in response to neighbors or predators, and
aerodynamic forces, which are based on the in-flight forces of lift, drag,
thrust and gravity. Fig. 1b provides a summary of this approach.

Many recent models for both fish and birds rely on Reynold’s vector-based
model (Couzin et al., 2002) (Erra et al., 2004) (Hemelrijk and Hildenbrandt,
2008) (Hemelrijk and Hildenbrandt, 2011) (Hemelrijk and Hildenbrandt, 2012)
(Hemelrijk et al., 2015) (Silva et al., 2009) (Hildenbrandt et al., 2010)
(Costanzo et al., 2021) (Papadopoulou et al., 2023). Only a few models of
flocking take other approaches, such as fuzzy logic (Lebar Bajec et al., 2005)
(Demšar and Lebar Bajec, 2014), phase-transition physics (Vicsek and Zafeiris,
2012) or shape morphing (Ho et al., 2012).

Figure 2: Given bird at position $p_{i}$ with orientation expressed as a
quaternion $Q_{i}$, we use the inverse quaternion $Q^{-1}_{i}$ to transform
the direction of a neighboring bird $p_{j}$ from world space into the local
coordinate frame, allows us to find the yaw $\psi_{i}$ and pitch $\theta_{i}$
angles for the direction of a given target to attract or repel.

###  2.3 Model design

We introduce a new simulation model, Flock2 (or Floids), for social flocking
with flight. Our model consists of a high-level orientation model for social
factors which controls a low-level aerodynamic model. The aerodynamic model,
Flightsim, is fully integrated into Flock2 and provides actuated control over
the orientation of a fixed wing-body with lift, drag and gravity. Although
this is not a flapping model it is an improvement over previous methods for
aerodynamics in social flocking. An orientation model is implemented that more
realistically models the ability of a bird to turn rather than accelerate in a
specific direction.

###  2.4 Modeling social factors as orientation goals

The social factors of avoidance, alignment and cohesion in our system are
implemented as orientation targets, with the orientation of each bird
maintained as a quaternion $Q_{i}$. Given a bird $i$ at position $p_{i}$ the
avoidance term considers the direction to the nearest neighbor j. The
directional vector $p_{j}-p_{i}$ is transformed into the local reference frame
of the bird using the inverse quaternion $Q^{-1}_{i}$, as in Figure 2.

Avoidance

$\overrightarrow{v}_{avoid}=(\mathbf{p}_{j}-\mathbf{p}_{i})$ (3)
$T_{\psi_{i}}^{avoid}=k_{avoid}\cdot\tan^{-1}\left(\overrightarrow{v}_{avoid}\>Q^{-1}_{i}\right)\frac{1}{d},$
(4)
$T_{\theta_{i}}^{avoid}=k_{avoid}\cdot\sin^{-1}\left(\overrightarrow{v}_{avoid}\>Q^{-1}_{i}\right)\frac{1}{d},$
(5)

Within the local reference frame, the arctangent and arcsin give the local yaw
$\psi_{i}$ and pitch $\theta_{i}$ angles of the neighboring bird j with
respect to the forward direction of bird i. This can be understood as the
direction to bird j projected onto the sphere of vision of the current bird
(see Fig. 2). When searching for neighbors we only consider those within a
limited forward field-of-view.

Alignment considers the direction indicated by the average local velocity
$v_{j}$ for the $n_{i}$ neighbors of bird i.

$\overrightarrow{v}_{align}=\frac{\sum_{j}\mathbf{v}_{j}}{n_{i}}$ (6)
$T_{\psi_{i}}^{align}=k_{align}\cdot\tan^{-1}\left(\overrightarrow{v}_{align}\>Q^{-1}_{i}\right),$
(7)
$T_{\theta_{i}}^{align}=k_{align}\cdot\sin^{-1}\left(\overrightarrow{v}_{align}\>Q^{-1}_{i}\right),$
(8)

For cohesion we find the direction given by the average local position, or
centroid, of the $n_{i}$ neighbors of bird i.

$\overrightarrow{v}_{cohesion}=\left(\frac{\sum_{j}\mathbf{p}_{j}}{n_{i}}-\mathbf{p}_{i}\right)$
(9)
$T_{\psi_{i}}^{cohesion}=k_{coh}\cdot\tan^{-1}\left(\overrightarrow{v}_{cohesion}\>Q^{-1}_{i}\right),$
(10)
$T_{\theta_{i}}^{cohesion}=k_{coh}\cdot\sin^{-1}\left(\overrightarrow{v}_{cohesion}\>Q^{-1}_{i}\right),$
(11)

The constants $k_{align}$ and $k_{coh}$ express positive attraction for
alignment and cohesion respectively, whereas the parameter $k_{avoid}$ is
negative to express the strength of avoidance.

In each of these sets of equations the 3D vectors for the nearest bird
$\overrightarrow{v}_{avoid}$, the average neighborhood velocity
$\overrightarrow{v}_{align}$, and the neighborhood centroid
$\overrightarrow{v}_{cohesion}$, are identical to those used in the Reynold’s
force model. The key difference is that, rather than using these as forces
directly, we then project these vectors onto the visual sphere of the bird to
derive orientation targets $T_{\phi_{i}},T_{\theta_{i}}$ for each term.

Figure 3: With only the three terms for avoidance, alignment and cohesion, and
without roosting or peripheral bounds, simulations produce snake-like
formations rather than spherical or ovoid patterns. We introduce a peripheral
bound for birds on the margin to seek cover within the flock. This applies to
any bird which is more isolated and is distinct from the simulation domain
bounds. See section 2.5.

###  2.5 Peripheral boundary term

Boundary cohesion is also an important feature of murmurations that we
investigated further. In starlings the flock boundaries are highly structured
and coherent with few outlying birds. In previous models an additional term
for attraction to a preferred direction (Couzin et al., 2005) or to a roosting
site (Hildenbrandt et al., 2010) was often introduced to aid in grouping. This
is generally unsatisfying as a continuous in-flight model since roosting is a
ground behavior. Caccamise (1990b) found that most birds are attracted more to
feeding areas called diurnal activity centers (DACs). Birds nest for the night
at a communal roost, which may be near a DAC, yet European starlings are “less
faithful to roost sites than to their DACs.” They may commute 3 to 12
kilometers to 12 different roosts (Caccamise, 1990b). Neither feeding DACs nor
roosts are satisfactory explanations for the strong, near instantaneous in-
flight boundary cohesion of murmurations. Yet some peripheral cohesion is
required to present more spherical or ovoidal shaped flocks.

In models containing only the three classical Reynolds terms of avoidance,
alignment and cohesion - whether implemented as force vectors or as
orientation targets - we found these terms alone result in flocks that only
form into echelon or tubular snake-like patterns, Figure 3, while rarely
forming spherical or ovoidal shapes which are the most common (Carere et al.,
2009).

Marginal predation provides a behavioral basis for modeling peripheral birds
(Hirsch and Morrell, 2011) (Duffield and Ioannou, 2017). The theory of
marginal predation suggests that birds on the periphery are more exposed in-
flight, and therefore have a greater evolved tendency and inclination to
return to the flock. Hildenbrandt et al. (2010) explore peripheral coherence
by measuring the degree of centrality $C_{i}$ within the local topological
group, where the direction they use for this adjustment is still the Reynold’s
social force for cohesion.

We model peripheral birds as having a rapid, in-flight, evolved response based
on maintaining self-cover, whereby they not only cohere together more
strongly, but avoid the periphery by plunging into the overall flock. This is
modeled as a distinct social orientation term. Those in the center of the
flock are well protected all around so they only need avoid neighbors and stay
aligned. On the periphery, birds may sense that they are becoming isolated and
turn back toward the center. This evolved response is plausible based on an
evolutionary response to marginal predation, and does not require the explicit
presence of a predator.

The degree to which a bird is peripheral in the flock can be estimated from
the neighborhood density 0 < $n_{i}/B$ < 1, where B is an adjustable constant
defining the thickness of the border region. For these birds, we turn them
back toward the flock centroid $C$ with a strength based on density.

$\overrightarrow{v}_{bound}=C-\mathbf{p_{i}}$ (12)
$T_{\psi_{i}}^{bound}=k_{bound}\left(\frac{B-n_{i}}{B}\right)\cdot\tan^{-1}\left(\overrightarrow{v}_{bound}\>Q^{-1}_{i}\right),$
(13)
$T_{\theta_{i}}^{bound}=k_{bound}\left(\frac{B-n_{i}}{B}\right)\cdot\sin^{-1}\left(\overrightarrow{v}_{bound}\>Q^{-1}_{i}\right),$
(14)

When $n_{i}=0$ the bird has no neighbors in the field of vision and the
peripheral turning strength is 1. When $n_{i}=B$ the bird has maximal
neighbors and the peripheral turning strength becomes 0.

###  2.6 Total orientation target

The input to the aerodynamic model is a target orientation T that expresses
the direction that any bird wishes to fly in yaw and pitch, as shown in Figure
1c. To determine the total desired target orientation of bird i, we sum the
turning influence from all social factors.

$T_{\psi_{i},\theta_{i}}^{total}=T_{\psi,\theta}^{avoid}+T_{\psi,\theta}^{align}+T_{\psi,\theta}^{cohesion}+T_{\psi,\theta}^{bound}$
(15)

Care must be taken since the sum takes place on a sphere. Yaw $\psi$ values
must be permitted to wrap around +/- 180 degrees from forward, while pitch
$\theta$ values wrap around +/- 90 degrees.

###  2.7 Aerodynamic model

Real bird flight is analyzed with fixed wing models, or more commonly an
aerodynamic flapping model (Norberg, 1985) (Pennycuick, 1968). With respect to
simulating flocks of birds, the model by Hildenbrandt et al. (2010) was the
first to introduce aerodynamic flight based on fixed wings. In their model the
body and wings are always oriented along the forward velocity $v_{i}$, and the
social forces $F_{social}$ are not coupled to the aerodynamic forces
$F_{flight}$. Nonetheless, their model was the first to exhibit lift loss when
banking, whole flock turning, and vertical speed changes (climbing and
diving).

Typically the bird is oriented along the direction of forward velocity $v_{i}$
(Couzin et al., 2002) (Hildenbrandt et al., 2010) (Hemelrijk and Hildenbrandt,
2011) (Reynolds, 1987). On the other hand, a complete aerodynamic model for
fixed wing flight should include angular velocity, spin, rotational inertia,
and torque generated by each lift surface (Yechout et al., 2014). Our flocking
model relies on Flightsim, a fixed wing model that resides between these,
providing a level of control in which a target orientation deflects the
current heading, and the dynamic stability of the wing reorients the moving
body along the direction of motion with a delay.

At each time step the current orientation $Q_{i}$ defines the local axes for
the forward (F), up (U) and side (S) directions.

$F=<1,0,0>Q_{i},U=<0,1,0>Q_{i},S=<0,0,1>Q_{i}$ (16)

Deflection of the forward direction F’ is computed from the input orientation
target $T_{\psi_{i},\theta_{i}}$ by constructing quaternions that represent
the turning rate around the up (U) and side (S) axes respectively:

$Q_{\psi}=Q_{angleaxis}\left((dt/k_{r})T_{\psi_{i}}^{total},U\right)$ (17)
$Q_{\theta}=Q_{angleaxis}\left((dt/k_{r})T_{\theta_{i}}^{total},S\right)$ (18)
$F^{\prime}=FQ_{\psi}Q_{\theta}$ (19)

The operator $Q_{angleaxis}$ constructs a quaternion given an angle and axis.
F’ represents the modified forward direction vector toward the desired target,
whose angle with F depends on the reaction rate. At least two factors
influence the flight response rate of birds: 1) the reaction time of the
vision system which is measured to be around 76 ms (Pomeroy and Heppner,
1977), and 2) the rate at which a starling can turn to a desired orientation.
Since a measured value of the latter is unknown we simplify the turning
reaction speed $k_{r}$ as a single parameter: time to achieve a desired
orientation, with 250 ms used in our simulations.

Lift and drag are defined by the aerodynamic equations:

$F_{lift}=\frac{1}{2}\rho v^{2}AC_{L}$ (20) $F_{drag}=\frac{1}{2}\rho
v^{2}AC_{D}$ (21)

where $\rho$ is the air density, A is the wing area, $v$ is the airflow
velocity (bird velocity along F plus wind), and $C_{L}$ is the coefficient of
lift. We assume a constant weight, wing area and air density similar to
Hildenbrandt et al. (2010), with one notable exception: the angular difference
between the modified forward direction F’ and the body velocity $v$ of the
wing allows us to define the wing angle-of-attack in the classical sense.

$aoa=cos^{-1}(F^{\prime}\cdot v)$ (22)

The angle of attack is used to lookup a variable coefficient of lift $C_{L}$
defined by a standard wing profile with stalls occurring at higher angles if
power is not increased. Due to wing shape and flapping, birds are better able
to avoid stalls so this factor can be adjusted or eliminated (Dhawan, 1991).

Thrust is applied in the direction of motion. Realistically this would be
caused by flapping but is treated here as a constant force in the direction of
motion.

$F_{thrust}=\hat{v}\cdot k_{power}$ (23)

Forces for lift, drag, gravity and thrust are combined as 3D vectors to give
the total aerodynamic forces acting on the body.

$F_{flight}=F_{lift}+F_{drag}+F_{thrust}+F_{gravity}$ (24)

The principle of dynamic stability states that a deflected aerodynamic body
will naturally want to re-orient along its direction of motion, not
instantaneously, but with greater strength as the deflection increases. The
final step in the aerodynamic model applies this principle to fractionally re-
orient the bird $Q_{i}$ from the targeted direction F’ toward the body
velocity $v_{i}$ with the weighting factor $k_{s}$. In the following equation
the $Q_{fromto}$ operator constructs a quaternion rotating one vector to
another with fraction k.

$Q_{i}^{t+1}=Q_{i}^{t}\cdot Q_{fromto}\left(F^{\prime},v_{i},k_{s}\right)$
(25)

While the base aerodynamic model was intended for fixed wing aircraft we
assume that birds are capable of much greater control. In our model birds are
assumed to be able to adjust their roll, pitch and yaw rapidly at will merely
by changing wing shape.

####  2.7.1 Time integration

The complete system including both the social factors and aerodynamic model is
shown in Figure 1c. These are evaluated for each bird and integrated at each
time step. The neighbor search is performed first using efficient grid-based
acceleration techniques for every bird to find its seven topological
neighbors. For all birds, the current bird positions are used to find the
desired target orientations $T_{\psi_{i},\theta_{i}}$ based on the social
factors - this constitutes the bird’s perceptual activity. These orientation
targets update the aerodynamic model to give the total aerodynamic force
$F_{flight}$ and the new orientation $Q_{i}^{t+1}$ \- this constitutes the
bird’s mechanical activity. Lastly, position and velocity are integrated with
Euler’s method to update the bird for the next time step.

$v_{i}^{t+1}=v_{i}^{t}+\left(F_{flight}/mass\right)dt$ (26)
$p_{i}^{t+1}=p_{i}^{t}+v_{i}^{t}dt$ (27)

Note that only the $F_{flight}$ forces from the aerodynamic model are
permitted to update the bird velocities and positions.

##  3 Results

Our model produced emergent behaviors exhibiting many characteristics of real
flocks that were not explicitly programmed, including non-penetrating flock
collisions, vortex formations, plausible splitting and merging, and both low
and high frequency turn propagation. We define spontaneous orientation waves
as visible, low frequency waves that traverse a flock whether they occur with
a predator or not. Orientation waves were especially surprising to observe in
our model since they were initially not expected, not explicitly programmed
(emergent), and occurred in the absence of any predefined trigger such as a
predator.

Consistent with real flocks our introduction of the boundary term allowed for
the spherical and ovoidal pattern to be observed most frequently (Fig. 11a).
However, formations of miniflocks, singletons, and indefinite shapes could be
also observed within individual simulation runs (Fig. 11). This is elaborated
in the discussion section.

Multiple experiments were conducted to better understand our model. To
quantify our model in relation to others we performed an aggregate energy
analysis of bird flight, section 3.1. Since our primary contribution is an
orientation-based model, we perform a frequency and spectrum power analysis to
better understand the oscillations produced during orientation waves, section
3.2. We also perform sensitivity analysis with multiple simulations over
isolated model parameters, section 3.3. Finally, we conducted experiments
between two flocks to test the response rate of simulated birds to head-on
collisions, section 3.4.

We study the Flock2 model with the boundary term (F2) in comparison to the
Reynolds model (REY). Since existing models are all derived from Reynold’s
classical vector-based force model, we choose to compare against this original
work. The model of Hildenbrandt et al. (2010) includes aerodynamic forces, but
in their work these are separated from the social forces we wish to examine
with respect to energy. We simulated between 5,000 and 20,000 birds in our
starling model over a 200 $m^{2}$ region, with the number of birds and
simulation runs described in each experiment type.

Figure 4: Residual forces are found to analyze simulated energies in any model of flight. A certain portion of energy is applied to overcome drag, $-F_{drag}$. The residual $F_{residual}$ is the remaining amount of force which is applied to bird motion after lift, drag and gravity have been subtracted. Of this, the forward component of work along the previous velocity direction $v_{t}$ is applied to accelerate or decelerate the bird, $P_{accel}$. The side component, $P_{turn}$, is applied to turn the bird. In this figure we assume dt=1 and mass=1 for clarity. Output | F2 (ours) | REY | REY : F2 | REY
---|---|---|---|---
| watts | watts | % | % of $P_{mech}$
$P_{lift}$ | 7.80 | 8.30 | 106% | 76.54%
$P_{drag}$ | 2.62 | 2.51 | 95% | 23.15%
$P_{accel}$ | 0.0035 | 0.0020 | 56% | 0.02%
$P_{turn}$ | 0.00008 | 0.0320 | 39500% | 0.30%
$P_{mech}$ | 10.42 | 10.84 | 104% | 100.00%
Ave. Speed | 10.02 | 10.18 | 101% |
Table 1: Energy comparison of bird output in mechanical work averaged over
10,000 birds for the Flock2 and Reynolds models. Simulations are calibrated
using equilibrium parameters before analysis.

###  3.1 Energy analysis

Simulations were conducted to compare energy analysis between our model (F2)
and Reynold’s model (REY). As our model was designed in standard international
units (SI) we are able to compute mechanical energy in watts. To calibrate
both models we consider the equilibrum condition of steady horizontal flight.

Mechanical power $P_{mech}$ is the energy required to maintain flight, the
primary components of which are lift and drag. For steady flight, since we do
not model a flapping wing, our goal is to calibrate energy by finding the
coefficient of lift $C_{L}$ where the lift force balances that of gravity. We
assume a mass of 0.08 kg, wing area of 0.0224 $m^{2}$ and average speed of 10
m/s for starlings (Dial et al., 1997)(Bruderer and Boldt, 2001).

$\begin{split}|F_{lift}|=|F_{grav}|=\frac{1}{2}\rho v^{2}AC_{L}&=mg\\\
C_{L}=0.5714\end{split}$ (28)

Using a lift-to-drag ratio of 3.3 for starlings from [Wither 1981] we can find
$C_{D}=C_{L}/\frac{L}{D}=0.1731$. In equilibrium the forward thrust force
balances the total drag.

$\begin{split}|F_{thrust}|=|F_{drag}|=\frac{1}{2}\rho
v^{2}AC_{D}&=k_{power}\\\ k_{power}=0.2373\end{split}$ (29)

Applying these equilibrium parameters for $C_{L}$, $C_{D}$ and $k_{power}$ to
our simulation, with social factors eliminated, results in average energy for
induced power $P_{ind}$ of 7.8 watts, combined parasitic and profile power of
$P_{para}$ = 2.3 watts and total mechanical power of $P_{mech}$ = 10.0 watts.
These values match the typical range of power observed for Sturnus vulgaris
(Rayner, 1999) (Ward et al., 2004).

To calibrate the Reynold’s simulation we must make assumptions since there are
no aerodynamic forces in this model. Simulations typically converge to the
minimum bird speed parameter which we set equal to our model at 10 m/s. Since
there is no explicit orientation in Reynold’s model we assume the lift is
vertical against gravity, and drag is opposed to the forward velocity. These
forces are emulated to compute and analyze energy; they do not alter the
simulation itself.

We are interested in understanding the remaining energies used for forward
acceleration and for turning. Our hypothesis is that Reynold’s model greatly
exaggerates turning power since social forces applied in an arbitrary
direction are not something a bird is realistically capable of producing. We
define acceleration energy, $P_{accel}$, as any additional power applied after
overcoming drag to increase or decrease the bird’s speed, and turning energy,
$P_{turn}$, as any additional power applied to turn toward a given target, as
in Fig. 4. To analyze this, we compute these residual energies equivalently in
both models by subtracting the lift, drag and gravity forces from the model-
specific force that updates velocity ($F_{model}$ = $F_{social}$ in Reynolds,
$F_{flight}$ in ours).

$F_{residual}=F_{model}-F_{lift}-F_{drag}-F_{grav}$ (30)
${\Delta}v=F_{residual}\>dt/mass$ (31)
$P_{accel}=|F_{residual}|\>\left({\Delta}v\cdot\hat{v}\right)$ (32)
$P_{turn}=|F_{residual}|\>|{\Delta}v-\left({\Delta}v\cdot\hat{v}\right)\hat{v}|$
(33)

The component of residual force measured as work along the direction of motion
gives $P_{accel}$, while the perpendicular component gives $P_{turn}$, as
shown in Fig. 4. This framework enables the comparison of bird energy usage in
both models equivalently, in real units of watts.

Figure 5: Total mechanical energy curves for both Flock2 (F2) and Reynolds
(REY). Each point in this scatter plot for a simulation of 10,000 birds
represents a bird at a specific velocity and mechanical energy output,
$P_{mech}$. These curves exhibit the right-hand side of the theoretical
U-shaped curve predicted for real birds (Pennycuick, 1968). In the Reynolds
model the speed must be artificially limited to a minimum since bird
velocities will converge to the minimum. In our model, F2, thrust determines
the mean velocity. The F2 curve is composed of points yet appears solid in
this plot as the power and velocities are more evenly and symmetrically
distributed. Figure 6: Total turning energy curves for the Flock2 model (top)
and Reynolds model (bottom). Each point in this scatter plot for a simulation
of 10,00 birds represents a bird at a specific velocity and $P_{turn}$. Notice
that the F2 turn energies are two orders of magnitude smaller than the REY
turn energies. Birds in the Reynolds model have a narrower range of speeds
that converge to the minimum.

Simulations were done with 10,000 birds in both models with the critical
parameters $C_{L}$, $C_{D}$ and $k_{power}$ matching the equilibrium condition
for a bird traveling 10 m/s. In the Reynolds model birds converge to the
minimum speed parameter, therefore we artificially set this limit to 10 m/s
for energy comparisons. In our Flock2 model, the minimum and maximum speeds
are set to 5 m/s and 18 m/s respectively with the average bird velocity
determined by $k_{power}$. Table 1 shows the bird energy output averaged over
all birds. While the bird energies for lift, drag and acceleration are
similar, the bird turning energies in Reynolds model were 39.5x (39500%)
larger than our Flock2 model. Total mechanical energies and turning energies
for all 10,000 birds are shown as scatter plots relative to bird velocity in
Fig. 5 and 6, respectively. These energy results are compared to theories of
bird flight in the discussion section.

Figure 7: Detail of an orientation wave appearing in a Flock2 simulation of
10,000 birds. Birds are rendered as darts so that their aspect to the camera
causes them to appear more or less opaque based on banking angle. Gray arrows
highlight the crests of a moving wave where birds are presenting more wing
area to the camera (appearing darker). The wave travels backward relative to
bird motion at a speed much greater than the flock. See accompanying videos.
Figure 8: Frequency analysis for the new Flock2 model (left) and Reynolds
model (right). a) Angular accelerations for 200 birds are plotted here from a
simulation of 10,000 birds and analyzed using b) Short-Time Fourier Transforms
(STFT) to better understand the emergent low and high frequency oscillations
observed in bird turns. c) The spectral energy of frequencies below 50 hz are
accumulated to provide a picture of low frequency oscillations. Strong peaks
were observed in F2 which directly corresponded to visible orientation waves
in simulation videos.

###  3.2 Frequency analysis

Orientation waves are an emergent behavior in our simulation that occur as
oscillations in bird orientation and which appeared in videos as rippling
waves, shown in Fig. 7, that traverse from the point of impact of sub-flocks.

To examine this aspect we perform Fourier analysis on bird accelerations with
the discrete-time Short-Time Fourier Transform (DT-STFT) to detect the
frequency range of these oscillations [Sejdic 2009]. In our case, the
magnitude of the angular acceleration of each bird, $|\alpha_{i}|$, is used as
sample input. These inputs were scaled up by 100x in the Reynolds model prior
to analysis to match the scale of Flock2 inputs, so that the STFT operates on
data of equivalent scale across models. Analysis is conducted on both our
Flock2 model and Reynolds model for a simulation of 10,000 birds.

Frequency anaylsis consists of the DT-STFT defined by a 1D real-to-complex
Discrete Fourier Transform (DFT), implemented with the FFTW library, and
applied over a moving window to isolate frequencies at each moment in time
(Jacobsen and Lyons, 2003) (Frigo and Johnson, 2005). This identifies the
temporal frequency spectrum for a single bird at a moment in time.

$X_{i}(f,t)=\sum_{n=-N/2}^{N/2}|\alpha_{i}(t+n)|\>H[n]e^{-j2\pi tk/N}$ (34)

The right side is a sliding window over the angular accelerations of a
specific bird, $\alpha_{i}$, at time t with a Hanning window filter $H_{n}$ of
width N, where N=512 is the number of frequency bins. The output X is a
complex number for each frequency bin $f$, at each time $t$, for a given bird
$i$. The magnitude squared of X gives the power of a particular frequency.

$F_{i}(f,t)=|X_{i}(f,t)|^{2}$ (35)

The STFT detects the frequencies for a 1D signal; one bird. To summarize the
entire flock we repeat the STFT for all birds (B) and accumulate the power at
each frequency. If we had accumulated the input accelerations before hand, the
differing phases of the birds would cause signals to interfere and cancel out.
By accumulating power after the STFT we can isolate and identify unique
frequencies throughout the flock. If many birds are turning at 10 hz, for
example, the power of the 10 hz frequency will be amplified regardless of
where they are in the wave.

$F_{total}(f,t)=\sum_{i=0}^{B}F_{i}(f,t)$ (36)

The spectrograms of $F_{total}$ for Flock2 and Reynolds simulations of 40
seconds are shown in Fig. 8b and 8f, respectively. The simulation rate is 5
milliseconds/frame therefore the sample rate is 200 hz. The Nyquist frequency
is 100 hz and the bin size is 0.39 hz/bin for N=512 bins. The spectrum
exhibited by our Flock2 model has strong, specific, low and high frequency
waves that are isolated in time.

Strong low frequency signals in Figure 8b were directly correlated with the
observations of orientation waves in simulation vidoes, occuring most often
when two sub-flocks interacted. This makes sense as the turning frequency must
be low enough to be visually observable as a wave. To quantify this, we
compute the total spectral energy over all frequencies below 50 hz as shown in
Fig. 8c and 8g.

$\text{SE}(t)=\sum_{f=0}^{50\>hz}F_{total}(f,t)$ (37)

Peaks in total spectral energy with positive curvature and a deviation greater
than 1% from a fit line were detected and marked. These peaks were found to
correspond very well with orientation waves observed in videos. We also
confirm this result by noting that spectural power peaks, Fig. 8c and 8g, also
correlate well with peaks in turning energy, $P_{turn}$, found in the earlier
force analysis and plotted in Fig. 8d and 8h. This provides a cross-validation
between the energy and frequency analyses that were performed.

###  3.3 Sensitivity analysis

Sensitivity analysis was conducted on various model parameters to understand
how they impacted flock activity. These results are reported in Appendix B.
Each simulation run was performed with 5000 birds, with a time step of 10
milliseconds/frame, for 40 seconds giving a total of 4000 frames/run. A
settling time of 5 seconds (500 frames) was allowed for the random initial
positions of birds to settle into a nominal state. Parameters for number of
birds, field of view, boundary size and reaction rate were independently
tested for spectral energy events (peaks >3%), total low frequency power (0 <
f <25 hz) and total high frequency power (75 < f < 100 hz).

Figure 9: A collision experiment between two flocks of 5,000 birds with a
timestep of 3 milliseconds. Purple represents birds with greater angular
acceleration, i.e. rapidly turning. Gray arrows indicate initial velocities.
a) Birds continue forward until b) the leading edge of birds responds to an
upcoming collision by rapidly adjusting orientation, which c) propagates
quickly to birds further from the edge, until d) reaching the entire flock.
Angular acceleration waves may propagate e) asymmetrically at different times.
On the order of one second the two flocks have established a new direction
that avoids collision.

###  3.4 Collision experiments

Collision experiments were inspired by a live murmuration event seen in real
starlings at 1:04 min in a public video by Clive and Smith (2011). At this
moment the boundaries of two sub-flocks of starlings collide head on and then
rebound without interpenetrating. We sought to replicate this. An experiment
with two flocks of 5,000 birds was conducted with a timestep of 3 milliseconds
over 1 second in the Flock2 and Reynolds models. Upon approaching another
flock in the F2 model, Fig. 9, leading birds turn quickly to avoid a collision
which propagates rapidly through the entire flock. This high frequency
response wave causes the two flock boundaries to avoid interpenetrating. Using
the same time step for the Reynolds model, shown in the accompanying video,
results in the two flocks merging and collapsing into a flattened disc rather
than avoiding one another. In our model, the two flocks detect a collision,
propagate turn information and resolve to a new direction without penetrating
peripheral boundaries on the order of one second.

Figure 10: A plausible mechanism for emergent orientation waves. They arise in
our model through the coupling of alignment (orange) and avoidance (purple)
within the bird’s field of vision, leading to an oscillatory turning as the
bird approaches and turns away from escaping birds in front. See discussion.

Figure 11: Flocking patterns exhibited by Flock2 (top) and Reynolds models
(bottom). Ours exhibited a) spherical, b) frontal, c) globular, d) mini-flock,
e) indistinct, and f) vortex patterns even within a single simulation run with
7 topological neighbors and the parameters given in Appendix A. Reynold’s
model shows patterns we describe as g) clustered, h) homogenous, i) large
amorphous, j) snake-like, k) rivers and l) tendrils, and has a more homogenous
fluid-like appearance that converges to a single large flock moving uniformly,
and required multiple runs with different settings.

##  4 Discussion

This work presents an improved simulation model, Flock2, for social flocking
that exhibits novel behaviors for boundary cohesion, orientation waves and
flocking patterns. We make this work available open source so that the
community may improve on and benefit from these computational experiments. To
our knowledge this is the first model to employ social orientation where the
behaviors of avoidance, alignment and cohesion are implemented as turn targets
rather than force vectors. Aerodynamic forces are computed as 3D vectors for
lift, drag, thrust and gravity, yet oriented solely through actuation of the
target heading.

Flocking patterns in our model are complex and appear to better mimic those in
nature, Figure 11. Carere et al. (2009) classifies the aerial flocking
patterns of European starlings into the categories of indefinite, globular,
miniflock, frontal, ovoid, snake-like, singleton, spherical, linear, V-shaped,
and planar. They found that the indeterminate, ovoidal and miniflock patterns
were most common over 1537 observations. The ovoidal pattern was also the most
common in our new simulation model. We were also able to observe spherical,
frontal, globular, mini-flock, indistinct and vortex patterns, similar to real
flocks, even within a single simulation run. In Reynold’s model, flocks appear
to have a more fluid-like nature that is either homogenous, clustered, or
tendril-like (Mototake and Ikegami, 2015). We observed that Reynold’s model
converges over a long simulation run into a single large flock moving
uniformly as confirmed by others (I. Lebar Bajec and Mraz, 2007), whereas ours
continues to split and merge indefinitely.

We distinguish between orientation waves, without predators, and agitation
waves which require predators as a trigger, and suggest that general
orientation waves may appear in real starling murmurations in the absence of
explicit predators. Our work does not refute the fact that agitation waves are
the result of an escape response to predator attacks (Procaccini et al.,
2011), yet it also validates that aerial displays of orientation waves may not
require an explicit attack by a predator (Carere et al., 2009). Rather,
orientation waves might also occur more generally when the flock self-
intersects or when the flock density suddenly triggers oscillating motions
within the flock.

Our method does not require a transmission mechanism (Hemelrijk et al., 2015)
as orientation waves are emergent in this system. The mechanism of action may
be explained as a feedback coupling of the competing turning forces of
alignment and avoidance within the bird’s field of vision (Fig. 10). When only
force-based avoidance is used, as in previous models, spontaneous orientation
waves do not occur. In the Flock2 model, a bird traveling with a flock (Fig.
10a) may approach too rapidly if birds ahead have reoriented, causing it to
turn in avoidance (Fig. 10b). Once a safer orientation is found, the bird is
compelled again to return to the flock (Fig. 10c), only needing to dodge again
(Fig. 10d). Sensitivity analysis shows this is coupled to the bird’s field of
view (Figure B.13).

Introduction of the peripheral boundary term relies on the fact that flocking
is an evolved response (Wood and Ackland, 2007) (Beauchamp, 2004) (Kwasnicka
et al., 2011). Bird flocks form a cohesive periphery in our model not because
a predator is immediately present, but because we assume they have evolved in
the predator-prey relationship to aggregate for protection. Proximity to the
periphery, or exposure, is the only variable required to model this. Since
Reynold’s original model did not account for peripheral bounds it is unable to
exhibit the most common ovoidal and spherical flocks our model is capable of.
This explanation is also more satisfying than one where attraction to a roost
is used to provide in-flight flock cohesion.

##  5 Conclusions & Future work

By decoupling aerodynamic and social factors, and requiring that social
orientation targets mimic the vision system to provide only heading goals for
an underlying aerodynamic model for the mechanics of flight, we believe this
system is more realistic with respect to both real bird flight and aerial
flocking. Nonetheless, an important limitation is that we do not yet model
flapping wings, which should require only a change to the aerodynamic model.

A more complete model would consider both the internal and external power in
flight including both mechanical and metabolic work. This requires a more
complete functional model of bird metabolism and would be an interesting
direction for future research.

Orientation waves are examined with energy and frequency analysis, yet a
deeper understanding of their arisal would be beneficial in both computational
models and real birds. An important question raised by this theoretical work
is the frequency of orientation waves in large murmurations of real birds, for
which there is still little experimental research since the focus of field
monitoring has been on predator-prey dynamics. This may require high
resolution, high speed cameras to provide a good control in ensuring there are
no predators actively flying in a very large flock.

In the future we hope to implement predators in order to compare different
triggers, types of waves, and mechanisms of propagation. Since predators are
known to trigger agitation waves we would expect them to emergently trigger
orientation waves.

We provide the Flock2 model as open source software in the hope that it
contributes to observational and computational biology by enabling further
research in aerial flocking, swarm dynamics and bird flight.

##  Appendix A Model Parameters

Parameter | Symbol | Value | Units
---|---|---|---
Social parameters |  |  |
Avoidance strength | $k_{avoid}$ | 0.02 |
Alignment strength | $k_{align}$ | 0.60 |
Cohesion strength | $k_{cohes}$ | 0.004 |
Boundary strength | $k_{bound}$ | 0.10 |
Boundary size | $B$ | 20 | # birds
Neighbor influence | $N_{topo}$ | 7 | # birds
Field of view | $fov$ | 290 | degrees
Flight parameters |  |  |
Mass | $m$ | 0.08 | kg
Wing area | $A$ | 0.0224 | $m^{2}$
Coeff. of lift | $C_{L}$ | 0.5714 |
Coeff. of drag | $C_{D}$ | 0.1731 |
Lift/drag ratio | L/D | 3.3 |
Reaction speed | $k_{r}$ | 250 | msec
Dynamic stability | $k_{s}$ | 70% |
Initial velocity | $v_{init}$ | 10 | $m/s$
Min velocity | $v_{min}$ | 5 | $m/s$
Max velocity | $v_{max}$ | 18 | $m/s$
Thrust | $k_{power}$ | 0.2373 | N
Environmental parameters |  |  |
Air density | $\rho$ | 1.225 | $kg/m^{3}$
Gravity | g | 9.8 | $m/s^{2}$
Simulation & testing |  |  |
Time step | $dt$ | 5 | msec/frame
Settle time | $t_{start}$ | 5 | sec
Run time | $t_{sim}$ | 20 | sec
FFT size, section 3.2 | N | 512 | bins
Sample rate | s | 200 | samples/sec
Nyquist freq. | $N_{f}$ | 100 | hz
Bin size | $B_{f}$ | 0.39 | hz/bin
Reynold’s parameters |  |  |
Avoidance strength | $R_{avoid}$ | 0.16 |
Alignment strength | $R_{align}$ | 0.60 |
Cohesion strength | $R_{cohes}$ | 0.08 |

##  Appendix B Sensitivity Analysis

Figure 12: Number of birds versus spectral events and energy for F2 and REY.
The spectral events, corresponding to orientation waves, are peaks in SE over
3% with average energy shown (squares). The total low frequency energy (<25
hz, circle) is also shown. The y-axis is in units of magnitude of angular
acceleration squared.  Figure 13: Field of view versus spectral events and
energy for F2 and REY. The spectral events, corresponding to orientation
waves, are peaks in SE over 3% with average energy shown (squares). The total
low frequency energy (<25 hz, cirlces) is also shown. The y-axis is in units
of magnitude of angular acceleration squared. The strongest event energy
corresponds very well to the actual field of view of Sturnus vulgaris, 296
degrees. Figure 14: Boundary count term (B), left, and reaction speed, right,
plotted against spectral events and energy in F2. Total turning energies in
low frequency (<25 hz, cirlces) and high frequency (>75 hz, triangles) are
shown. As the peripheral boundary term increases (width of boundary in number
of birds B) the total event energy, low frequency, and high frequency energy
all increase.

##  Appendix C Visual results

Figure 15: Visual output of Flock2. Birds are rendered as triangles with
orientation as required to exhibit orientation waves based on area of the
wings presented to the camera. Birds influenced by the peripheral boundary
term are colored orange. Extreme changes in angular acceleration are
visualized in green.

## References

  * Nikolsky [1955] G.V. Nikolsky.  On the biological significance of schooling in fishes.  pages 104–107, Moscow, Russia, 1955. Nauk SSSR.
  * Breder [1959] Charles M. Breder.  Studies on social groupings in fishes.  volume 117, pages 393–482, New York, 1959. American Museum of Natural History.
  * Pitcher and Wyche [1983] Tony J. Pitcher and Christopher J. Wyche.  _Predator-avoidance behaviours of sand-eel schools: why schools seldom split_ , page 193–204.  Springer Netherlands, 1983.
  * Godin and Morgan [1985] JG.J Godin and M.J. Morgan.  Predator avoidance and school size in a cyprinodontid fish, the banded killifish (Fundulus diaphanus Lesueur).  _Behav Ecol.Sociobiol_ , (16):105–110, 1985.  doi:https://doi.org/10.1007/BF00295142.
  * Magurran [1990] Anne E. Magurran.  The adaptive significance of schooling as an anti-predator defence in fish.  _Annales Zoologici Fennici_ , 27(2):51–66, 1990\.
  * Parrish [2000] Julia K. Parrish.  Self-organized fish schools: an examination of emergent properties.  _The biological bulletin_ , 202:296–305, 2000.  doi:https://doi.org/10.2307/1543482.
  * Hamilton [1971] W.D. Hamilton.  Geometry for the selfish herd.  _Journal of Theoretical Biology_ , 31(2):295–311, 1971.  ISSN 0022-5193.  doi:https://doi.org/10.1016/0022-5193(71)90189-5.
  * Treves [2000] Adrian Treves.  Theory and method in studies of vigilance and aggregation.  _Animal Behaviour_ , 60(6):711–722, 2000.  ISSN 0003-3472.  doi:10.1006/anbe.2000.1528.
  * Emlen [1952] Jr. Emlen, J.T.  Flocking behavior in birds.  _The Auk_ , 69:160–170, 1952.  doi:https://psycnet.apa.org/doi/10.2307/4081266.
  * Caraco et al. [1980] Thomas Caraco, Steven Martindale, and H. Ronald Pulliam.  Avian flocking in the presence of a predator.  _Nature_ , 285(5764):400–401, 1980.  ISSN 1476-4687.  doi:10.1038/285400a0.
  * Carbone et al. [2003] C. Carbone, W. A. Thompson, L. Zadorina, and J. M. Rowcliffe.  Competition, predation risk and patterns of flock expansion in barnacle geese (branta leucopsis).  _Journal of Zoology_ , 259(3):301–308, March 2003.  ISSN 1469-7998.  doi:10.1017/s0952836902003278.
  * Procaccini et al. [2011] Andrea Procaccini, Alberto Orlandi, Andrea Cavagna, Irene Giardina, Francesca Zoratto, Daniela Santucci, Flavia Chiarotti, Charlotte K. Hemelrijk, Enrico Alleva, Giorgio Parisi, and Claudio Carere.  Propagating waves in starling, Sturnus vulgaris, flocks under predation.  _Animal Behaviour_ , 82(4):759–765, 2011.  ISSN 0003-3472.  doi:10.1016/j.anbehav.2011.07.006.
  * Goodenough et al. [2017] Anne E. Goodenough, N. Little, W.S. Carpenter, and A. Hart.  Birds of a feather flock together: Insights into starling murmuration behaviour revealed using citizen science.  _PLOS ONE_ , 12:1–18, 06 2017.  doi:10.1371/journal.pone.0179277.
  * Vine [1971] Ian Vine.  Risk of visual detection and pursuit by a predator and the selective advantage of flocking behaviour.  _Journal of Theoretical Biology_ , 30(2):405–422, 1971.  ISSN 0022-5193.  doi:https://doi.org/10.1016/0022-5193(71)90061-0.
  * McNamara and Houston [1992] J.M. McNamara and A.I. Houston.  Evolutionarily stable levels of vigilance as a function of group size.  _Animal Behaviour_ , 43(4):641–658, 1992.  ISSN 0003-3472.  doi:https://doi.org/10.1016/S0003-3472(05)81023-1.
  * Sirot [2006] Etienne Sirot.  Social information, antipredatory vigilance and flight in bird flocks.  _Animal Behaviour_ , 72(2):373–382, 2006.  ISSN 0003-3472.  doi:https://doi.org/10.1016/j.anbehav.2005.10.028.
  * Siegried and Underhill [1975] W.R. Siegried and L.G. Underhill.  Flocking as an anti-predator strategy in doves.  _Animal Behaviour_ , 23:504–508, 1975.  ISSN 0003-3472.  doi:https://doi.org/10.1016/0003-3472(75)90126-8.
  * Feare [1984] Christopher Feare.  _The Starling_.  Oxford University Press, London, 1984.
  * Cresswell [1994] Will Cresswell.  Flocking is an effective anti-predation strategy in redshanks, Tringa totanus.  _Animal Behaviour_ , 47(2):433–442, 1994.  ISSN 0003-3472.  doi:https://doi.org/10.1006/anbe.1994.1057.
  * Beauchamp [2004] Guy Beauchamp.  Reduced flocking by birds on islands with relaxed predation.  _Proceedings of the Royal Society of London. Series B: Biological Sciences_ , 271(1543):1039–1042, 2004.  doi:https://doi.org/10.1098/rspb.2004.2703.
  * Fischl and Caccamise [1987] Joseph Fischl and Donald F. Caccamise.  Relationships of diet and roosting behavior in the European starling.  _The American Midland Naturalist_ , 117(2):395–404, 1987.  doi:https://doi.org/10.2307/2425982.
  * Tinbergen [1981] Joost Tinbergen.  Foraging decisions in starlings (Sturnus vulgaris l.).  _Ardea_ , 69:1–67, 01 1981.  doi:https://doi.org/10.5253/arde.v69.p1.
  * Caccamise [1990a] Donald F. Caccamise.  European starling fidelity to diurnal activity centers: Role of foraging substrate quality.  _The Wilson Bulletin_ , 103(1):13–24, 1990a.  doi:http://www.jstor.org/stable/4162966.
  * Carere et al. [2009] Claudio Carere, Simona Montanino, Flavia Moreschini, Francesca Zoratto, Flavia Chiarotti, Daniela Santucci, and Enrico Alleva.  Aerial flocking patterns of wintering starlings, Sturnus vulgaris, under different predation risk.  _Animal Behaviour_ , 77(1):101–107, 2009.  ISSN 0003-3472.  doi:https://doi.org/10.1016/j.anbehav.2008.08.034.
  * Couzin et al. [2002] Iain Couzin, Jens Krause, Richard James, Graeme Ruxton, and Nigel Franks.  Collective memory and spatial sorting in animal groups.  _Journal of Theoretical Biology_ , 218(1):1–11, 2002.  ISSN 0022-5193.  doi:https://doi.org/10.1006/jtbi.2002.3065.
  * Reynolds [1987] Craig W. Reynolds.  Flocks, herds and schools: A distributed behavioral model.  _ACM SIGGRAPH Computer Graphics_ , 21(4):25–34, 1987.  ISSN 0097-8930.  doi:10.1145/37402.37406.
  * Couzin et al. [2005] I. Couzin, J. Krause, and N. Franks.  Effective leadership and decision-making in animal groups on the move.  _Nature_ , 433:513–516, 2005.  doi:https://doi.org/10.1038/nature03236.
  * Parrish and Viscido [2005] J.K. Parrish and Steven Viscido.  Traffic rules of fish schools: A review of agent-based approaches.  _Self-organisation and evolution of social systems_ , pages 50–80, 01 2005.  doi:10.1017/CBO9780511542275.006.
  * Buhl et al. [2006] J. Buhl, D. J. T. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson.  From disorder to order in marching locusts.  _Science_ , 312(5778):1402–1406, 2006.  doi:10.1126/science.1125142.
  * Hemelrijk and Hildenbrandt [2008] Charlotte K. Hemelrijk and Hanno Hildenbrandt.  Self-organized shape and frontal density of fish schools.  _Ethology_ , 114(3):245–254, 2008.  ISSN 1439-0310.  doi:10.1111/j.1439-0310.2007.01459.x.
  * Hildenbrandt et al. [2010] H. Hildenbrandt, C. Carere, and C.K. Hemelrijk.  Self-organized aerial displays of thousands of starlings: a model.  _Behavioral Ecology_ , 21(6):1349–1359, 10 2010\.  ISSN 1045-2249.  doi:10.1093/beheco/arq149.
  * Caccamise [1990b] Donald F. Caccamise.  Communal starling roosts: Implications for control.  volume 14. UC Agriculture & Natural Resources, 1990b.
  * Pennycuick [1968] Colin J. Pennycuick.  Power requirements for horizontal flight in the pigeon Columba livia.  _The Journal of Experimental Biology_ , 49(3):527–555, 1968.  doi:https://doi.org/10.1242/JEB.49.3.527.
  * Withers [1981] Philip C. Withers.  An aerodynamic analysis of bird wings as fixed aerofoils.  _Journal of Experimental Biology_ , 90(1):143–162, 1981.  ISSN 1477-9145.  doi:10.1242/jeb.90.1.143.
  * Storms [2023] Rolf F Storms.  _Collective escape in bird flocks_.  PhD thesis, University of Groningen Press, 2023.
  * Hemelrijk et al. [2015] Charlotte K. Hemelrijk, Lars van Zuidam, and Hanno Hildenbrandt.  What underlies waves of agitation in starling flocks.  _Behavioral Ecology and Sociobiology_ , 69(5):755–764, 2015.  ISSN 1432-0762.  doi:10.1007/s00265-015-1891-3.
  * Cavagna et al. [2010] Andrea Cavagna, Alessio Cimarelli, Irene Giardina, Giorgio Parisi, Raffaele Santagati, Fabio Stefanini, and Massimiliano Viale.  Scale-free correlations in starling flocks.  _Proceedings of the National Academy of Sciences_ , 107(26):11865–11870, 2010.  doi:10.1073/pnas.1005766107.
  * Altshuler et al. [2015] Douglas Altshuler, Joe Bahlman, Roslyn Dakin, Andrea Gaede, Benjamin Goller, David Lentink, Paolo Segre, and Dimitri Skandalis.  The biophysics of bird flight: Functional relationships integrate aerodynamics, morphology, kinematics, muscles, and sensors.  _Canadian Journal of Zoology_ , 93, 10 2015.  doi:10.1139/cjz-2015-0103.
  * Rayner [1993] J Rayner.  On aerodynamics and the energetics of vertebrate flapping flight.  _Contemporary mathematics_ , 141:351–400, 1993.
  * Rayner [1999] Jeremy M. V. Rayner.  Estimating power curves of flying vertebrates.  _Journal of Experimental Biology_ , 202(23):3449–3461, 1999.  ISSN 1477-9145.  doi:10.1242/jeb.202.23.3449.
  * Norberg [1985] Ulla M. Norberg.  Evolution of vertebrate flight: An aerodynamic model for the transition from gliding to active flight.  _The American Naturalist_ , 126(3):303–27, 1985\.
  * Dial et al. [1997] K. P. Dial, A. A. Biewener, B. W. Tobalske, and D. R. Warrick.  Mechanical power output of bird flight.  _Nature_ , 390(6655):67–70, 1997.  ISSN 1476-4687.  doi:10.1038/36330.
  * Daan et al. [1990] S. Daan, D. Masman, and A. Groenewold.  Avian basal metabolic rates: their association with body composition and energy expenditure in nature.  _American Journal of Physiology-Regulatory, Integrative and Comparative Physiology_ , 259(2):R333–R340, 1990.  ISSN 1522-1490.  doi:10.1152/ajpregu.1990.259.2.r333.
  * Rayner et al. [2015] Jeremy M. V. Rayner, Paolo W. Viscardi, Sally Ward, and John R. Speakman.  Aerodynamics and energetics of intermittent flight in birds.  _American Zoologist_ , 41(2):188–204, 08 2015\.  ISSN 0003-1569.  doi:https://doi.org/10.1093/icb/41.2.188.
  * Ward et al. [2001] S. Ward, U. Moller, J. M. V. Rayner, D. M. Jackson, D. Bilo, W. Nachtigall, and J. R. Speakman.  Metabolic power, mechanical power and efficiency during wind tunnel flight by the european starling, Sturnus vulgaris.  _Journal of Experimental Biology_ , 204(19):3311–3322, 2001.  ISSN 0022-0949.  doi:10.1242/jeb.204.19.3311.
  * Ward et al. [2004] S. Ward, U. Moller, J. M. V. Rayner, D. M. Jackson, W. Nachtigall, and J. R. Speakman.  Metabolic power of European starlings Sturnus vulgaris during flight in a wind tunnel, estimated from heat transfer modelling, doubly labelled water and mask respirometry.  _Journal of Experimental Biology_ , 207(24):4291–4298, 11 2004.  ISSN 0022-0949.  doi:10.1242/jeb.01281.
  * Tobalske [2007] Bret W. Tobalske.  Biomechanics of bird flight.  _Journal of Experimental Biology_ , 210(18):3135–3146, 2007.  ISSN 0022-0949.  doi:10.1242/jeb.000273.
  * Nudds and Bryant [2000] Robert L. Nudds and David M. Bryant.  The energetic cost of short flights in birds.  _Journal of Experimental Biology_ , 203(10):1561–1572, 2000.  ISSN 1477-9145.  doi:10.1242/jeb.203.10.1561.
  * Schmidt-Wellenburg et al. [2008] Carola A. Schmidt-Wellenburg, Sophia Engel, and G. Henk Visser.  Energy expenditure during flight in relation to body mass: effects of natural increases in mass and artificial load in Rose Coloured starlings.  _Journal of Comparative Physiology B_ , 178(6), 2008.  ISSN 1432-136X.  doi:10.1007/s00360-008-0267-7.
  * Erra et al. [2004] Ugo Erra, Rosario De Chiara, Vittorio Scarano, and Maurizio Tatafiore.  Massive simulation using gpu of a distributed behavioral model of a flock with obstacle avoidance.  pages 233–240, 01 2004.
  * Hemelrijk and Hildenbrandt [2011] Charlotte K. Hemelrijk and Hanno Hildenbrandt.  Some causes of the variable shape of flocks of birds.  _PLoS ONE_ , 6(8):e22479, 2011.  ISSN 1932-6203.  doi:10.1371/journal.pone.0022479.
  * Hemelrijk and Hildenbrandt [2012] Charlotte K. Hemelrijk and Hanno Hildenbrandt.  Schools of fish and flocks of birds: their shape and internal structure by self-organization.  _Interface Focus_ , 2(6):726–737, 2012.  doi:10.1098/rsfs.2012.0025.
  * Silva et al. [2009] Alessandro Ribeiro Da Silva, Wallace Santos Lages, and Luiz Chaimowicz.  Boids that see: Using self-occlusion for simulating large groups on gpus.  _Computers in Entertainment_ , 7(4):1–20, December 2009.  ISSN 1544-3574.  doi:10.1145/1658866.1658870.
  * Costanzo et al. [2021] A. Costanzo, H. Hildenbrandt, and C. K. Hemelrijk.  Causes of variation of darkness in flocks of starlings, a computational model.  _Swarm Intelligence_ , 16(2):91–105, 2021.  ISSN 1935-3820.  doi:10.1007/s11721-021-00207-4.
  * Papadopoulou et al. [2023] Marina Papadopoulou, Hanno Hildenbrandt, and Charlotte K. Hemelrijk.  Diffusion during collective turns in bird flocks under predation.  _Frontiers in Ecology and Evolution_ , 11, 2023.  ISSN 2296-701X.  doi:10.3389/fevo.2023.1198248.
  * Lebar Bajec et al. [2005] Iztok Lebar Bajec, Nikolaj Zimic, and Miha Mraz.  Simulating flocks on the wing: the fuzzy approach.  _Journal of Theoretical Biology_ , 233(2):199–220, 2005.  ISSN 0022-5193.  doi:https://doi.org/10.1016/j.jtbi.2004.10.003.
  * Demšar and Lebar Bajec [2014] Jure Demšar and Iztok Lebar Bajec.  Simulated predator attacks on flocks: A comparison of tactics.  _Artificial Life_ , 20(3):343–359, 2014.
  * Vicsek and Zafeiris [2012] Tamás Vicsek and Anna Zafeiris.  Collective motion.  _Physics Reports_ , 517(3):71–140, 2012.  ISSN 0370-1573.  doi:https://doi.org/10.1016/j.physrep.2012.03.004.  Collective motion.
  * Ho et al. [2012] Choon Sing Ho, Yew-Soon Ong, Xianshun Chen, and Ah-Hwee Tan.  _FAME, Soft Flock Formation Control for Collective Behavior Studies and Rapid Games Development_ , page 258–269.  Springer Berlin Heidelberg, 2012.
  * Hirsch and Morrell [2011] Ben T. Hirsch and Lesley J. Morrell.  Measuring marginal predation in animal groups.  _Behavioral Ecology_ , 22(3):648–656, 05 2011\.  ISSN 1045-2249.  doi:10.1093/beheco/arr026.
  * Duffield and Ioannou [2017] Callum Duffield and Christos C Ioannou.  Marginal predation: do encounter or confusion effects explain the targeting of prey group edges?  _Behavioral Ecology_ , 28(5):1283–1292, July 2017.  doi:10.1093/beheco/arx090.
  * Yechout et al. [2014] Thomas R. Yechout, Steven L. Morris, David E. Bossert, Wayne F. Hallgren, and James K. Hall.  _Introduction to Aircraft Flight Mechanics_.  AIAA American Institute of Aeronautics & Astronautics, 2014.
  * Pomeroy and Heppner [1977] Harold Pomeroy and Frank Heppner.  Laboratory determination of startle reaction time of the starling (Sturnus vulgaris).  _Animal Behaviour_ , 25:720–725, 1977.  doi:10.1016/0003-3472(77)90121-X.
  * Dhawan [1991] S Dhawan.  Bird flight.  _Sadhana_ , 16(4):275–352, December 1991.  doi:10.1007/bf02745345.
  * Bruderer and Boldt [2001] Bruno Bruderer and Andreas Boldt.  Flight characteristics of birds: I. radar measurements of speeds.  _Ibis_ , 143:178 – 204, 04 2001.  doi:https://doi.org/10.1111/j.1474-919X.2001.tb04475.x.
  * Jacobsen and Lyons [2003] E. Jacobsen and R. Lyons.  The sliding dft.  _IEEE Signal Processing Magazine_ , 20(2):74–80, 2003.  doi:10.1109/MSP.2003.1184347.
  * Frigo and Johnson [2005] M. Frigo and S.G. Johnson.  The design and implementation of fftw3.  _Proceedings of the IEEE_ , 93(2):216–231, 2005\.  doi:10.1109/jproc.2004.840301.
  * Clive and Smith [2011] Sophie Windsor Clive and Liberty Smith.  Murmuration (official video) by Sophie Windsor Clive & Liberty Smith, Dec 2011.  URL https://www.youtube.com/watch?v=iRNqhi2ka9k.
  * Mototake and Ikegami [2015] Yoh-Ichi Mototake and Takashi Ikegami.  A simulation study of large scale swarms.  Kyoto, Japan, 10 2015.
  * I. Lebar Bajec and Mraz [2007] N. Zimic I. Lebar Bajec and M. Mraz.  The computational beauty of flocking: boids revisited.  _Mathematical and Computer Modelling of Dynamical Systems_ , 13(4):331–347, 2007.  doi:10.1080/13873950600883485.
  * Wood and Ackland [2007] Andrew J Wood and Graeme J Ackland.  Evolving the selfish herd: emergence of distinct aggregating strategies in an individual-based model.  _Proceedings of the Royal Society B: Biological Sciences_ , 274(1618):1637–1642, 2007.  ISSN 1471-2954.  doi:10.1098/rspb.2007.0306.
  * Kwasnicka et al. [2011] Halina Kwasnicka, Urszula Markowska-Kaczmar, and Marcin Mikosik.  Flocking behaviour in simple ecosystems as a result of artificial evolution.  _Applied Soft Computing_ , 11(1):982–990, 2011\.  doi:https://doi.org/10.1016/j.asoc.2010.01.018.