Density distributions and depth in flocks

TL;DR

This paper introduces a metric-free model for bird flocking that accounts for inhomogeneous density distributions and motional biases, successfully fitting real starling flock data and exploring how birds perceive their position within the flock.

Contribution

The paper presents a novel, metric-free collective behavior model incorporating distributed motional bias, explaining density variations and depth perception in bird flocks.

Findings

Inward bias at flock edges and outward bias inside the flock are necessary.

The model fits observed starling flock density data.

Simple visual analysis could enable birds to determine their depth.

Abstract

Recent experimental evidence suggests that interactions in flocks of birds do not involve a characteristic length scale. Bird flocks have also been revealed to have an inhomogeneous density distribution, with the density of birds near the border greater than near the centre. We introduce a strictly metric-free model for collective behaviour that incorporates a distributed motional bias, providing control of the density distribution. A simple version of this model is then able to provide a good fit to published data for the density variation across flocks of Starlings. We find that it is necessary for individuals on the edge of the flock to have an inward motional bias but that birds in the interior of the flock instead must have an outward bias. We discuss the ability of individuals to determine their depth within a flock and show how this might be achieved by relatively simple analysis of their visual environment.

Equations

$${\underline{r}}_i^{t+1}={\underline{r}}_i^t+v_0{\underline{\hat{v}}}_i^t$$

$${\underline{v}}_i^{t+1}=(1-{\varphi }_n){\underline{\hat{\mu }}}_i^t+{\varphi }_n{\underline{\hat{\eta }}}_i^t$$

$$\displaystyle\underline{\mu}_{i}^{t}=f_{i}\langle\underline{\hat{r}}_{ij}^{t}\rangle_{j\in S_{i}}+(1-f_{i})\vartheta\big{(}\langle\underline{\hat{v}}_{j}^{t}\rangle_{j\in B_{i}}\big{)}$$

$$f_i=\{\begin{array}{ll}{\varphi }_e& {\underline{r}}_i^t\in C_0\\ 0& \text{otherwise}\end{array}$$

$$N=\rho {\int }_S{\tilde{r}}^{2}d\tilde{r}d\Omega$$

$$\frac{dN}{d\Omega }=\frac{\rho L(\theta {)}^3}{3}:=I(\theta )$$

$$\frac{dI}{d\theta}=-\rho L^{2}r_{p}\sin\theta\Bigg{(}1+\frac{r_{p}\cos\theta}{\sqrt{R^{2}-r_{p}^{2}\sin^{2}\theta}}\Bigg{)}=0$$

$$\Delta I=\frac{I_{\min }}{I_{\max }}=\frac{(1-D{)}^3}{(1+D{)}^3}$$

Full text

• Keywords: swarming, topological interaction, metric-free, distributed motional bias, inhomogeneous swarm density, flock depth, self-organisation.

• This Accepted Manuscript is available for reuse under a CC BY-NC-ND 3.0 licence after the 12 month embargo period provided that all the terms of the licence are adhered to.

• This is a peer-reviewed, un-copyedited version of an article accepted for publication/published in Journal of Physics D: Applied Physics. IOP Publishing Ltd is not responsible for any errors or omissions in this version of the manuscript or any version derived from it. The Version of Record is available online at 10.1088/1361-6463/aa942f.

\ioptwocol

1 Introduction

Swarming is the collective behaviour of animal aggregations, and can be observed in the flocking of birds [1, 2, 3], fish shoaling [4, 5, 6], mammal herding [7, 8], and insect swarming [9, 10]. Human crowds can also display this sort of collective trait [11, 12]. The emergence of global orientational order in groups of moving animals is arguably the most striking consequence of this type of social behaviour [10, 13]. In these systems collective animal behaviour is not thought to arise from centralised coordination but rather the system is believed to exhibit self-organisation due to the local rules of the interacting elements. This results in coherent motion with local rules manifesting global order [14, 15].

In recent years, a large number of theoretical models have been developed in which local interaction rules give rise to global ordering in animal systems [16, 17, 18] however empirical studies have been more rare [19, 20]. Testing models against data is essential if we are to determine which sorts of model give rise to specific characteristics: many models can generate some form of swarming, but which of these models give rise to swarms that resemble those seen in nature? It has been suggested that the specific interaction mechanism may vary with species and for some systems an interaction based on neighbour distance appears to be a good fit [21]. In contrast, recent field studies have reconstructed the internal dynamics of large flocks of Starlings and have determined that their nearest-neighbour interactions do not depend on interaction range [22, 23, 24].

Developing models with this metric-free characteristic is technically challenging as they typically support a zero density steady-state, such as described in the work of Ginelli and Chaté [25] in which diffusive expansion continues indefinitely. Pearce and Turner [26] describe a model that regulates swarm density using a motional bias on surface individuals and topological interaction rules, preserving the metric-free nature of the model and also generating a steady-state with finite spatial extent. This Strictly Metric-Free (SMF) model is therefore useful to compare with observations of bird flocks as it can produce bounded swarms in open boundary conditions. However we will show that, in its simplest form, it yields density distributions that are rather different to those observed.

In this work we propose a fully topological (metric-free) 3-dimensional model which includes a motional bias that is tunable throughout the swarm and not just on its surface. This bias has a topological character, preserving the fully topological nature of the model. Our aim is to explore the regulation of density across flocks of birds. We are motivated by findings from a field study [22] that reports a nonhomogeneous density variation across flocks of Starlings, specifically a higher density at the border of the flock than in the centre. This observation is counter to what has been observed in some other models of collective behaviour [27]. It is also counter-intuitive in relation to some theories of animal behaviour, such as the selfish herd hypothesis [28] in which the centre of the group would be the safest location and all individuals might therefore be expected to seek to occupy it. We show that our metric-free distributed motional bias model is able to support behaviour consistent with these empirical observations.

The model is introduced in section 2. The methods used to measure aggregate densities and fit the model to data are described in section 3. The resultant model and swarm density profiles are presented in section 4. Additionally, a biologically motivated basis for an individual determining their depth from within the flock is presented and discussed in section 5. Concluding remarks are in section 6.

2 Description of the model

The model we propose begins with the surface bounding effect introduced in the SMF model [26] and extends it to act on all individuals in the aggregate with strength prescribed by a function of the topological depth of the individual within the swarm. In contrast to classic models of self-propelled particles, such as those by Vicsek et al. [29], we identify two particles to be neighbours if they are directly connected to each other under a Voronoi tessellation [30, 31]. This is constructed for the particle positions at each time step, thus defining interacting neighbours as those in neighbouring Voronoi cells (i.e. particles which share an edge in the Delaunay triangulation of all particle locations).

We use this tessellation to determine topological depth for each of the particles in the dynamic aggregate (flock). We first identify a shell, or set, of particles as being those that occupy an infinite Voronoi cell. These are denoted as occupying shell [math] and correspond to particles that are on the convex hull of the system [31]. Particles that are connected to these shell [math] particles via Delaunay edges, but that are not themselves members of shell [math], are defined to lie in shell $1$. This process is repeated iteratively until all particles are assigned a shell number. This labelling encodes topological depth as it relates to the shortest path length from the border through the graph defined via the Delaunay triangulation. A driving term can then be included in the equation of motion that provides a motional bias on each particle. The direction of this bias (loosely “inwards” or “outwards”), is derived using the locations of its neighbours on the same shell.

The interaction rules governing all $N$ identical particles in the system are shown in equations 1–3.

[TABLE]

They involve the position $\underline{r}_{i}^{t}$ of particle $i$ at discrete time $t$, having direction of motion $\underline{\hat{v}}_{i}^{t}$ and constant speed $v_{0}$, which is set equal to unity in what follows. The “hat” symbol $\hat{\;\;}$ denotes a normalised (unit) vector and angled brackets $\langle\cdots\rangle$ indicate an average over the indicated particle subset. The operator $\vartheta(\;)$ performs normalisation via $\vartheta(\underline{w})=\nicefrac{{\underline{w}}}{{|\underline{w}|}}$ and $\underline{\hat{r}}_{ij}^{t}$ denotes the unit vector pointing from particle $i$ to particle $j$ at discrete time step $t$. The parameter $\phi_{n}$ encodes the strength of the (vectorial [32]) noise applied to each particle, multiplied by a random unit vector obeying $\langle\underline{\hat{\eta}}_{i}^{t}\rangle=0$ and $\langle\underline{\hat{\eta}}_{i}^{t}\cdot\underline{\hat{\eta}}_{j}^{t^{\prime}}\rangle=\delta_{i,j}\delta_{t,t^{\prime}}$. The neighbours of particle $i$ are denoted $B_{i}$ and particles which share the same shell number $\kappa(i)$ as particle $i$ form the set $C_{\kappa(i)}$. Therefore we denote the set of shell neighbours of $i$ as the intersection $S_{i}=B_{i}\cap C_{\kappa(i)}$. We average over the unit vectors pointing from particle $i$ to members of this set. In addition, figure 1 shows how shell 0 is defined as the members on the convex hull of the system, and also the procedure for identifying all other shells. Equation (1) represents a simple vectorial particle translation along the current velocity. Equation (2) encodes an update rule for the velocity that includes both some deterministic driving terms, weight $(1-\phi_{n})$, and some stochastic noise, weight $\phi_{n}$. Thus $\phi_{n}$, the degree of noise, is an important control parameter in what follows. Equation 3 defines the deterministic driving terms. It is comprised of two terms, the first, with weight $f_{i}$, encodes the motional bias constructed from the shell geometry, as described and the second term, with weight $(1-f_{i})$, provides co-alignment of each particle with its neighbours.

We denote $f_{i}$ as the “bounding function”, which encodes the relative strength of the bounding effect on each Voronoi shell. Changing this allows us to tune the bounding of the model across the aggregation as we wish. If we choose $f_{i}$ to have the form of equation 4, where $\phi_{e}$ is a parameter controlling the strength of the border shell effect, then we can recover the Strictly Metric-Free (SMF) model of Pearce and Turner [26] in its entirety. If instead we choose $f_{i}=0\;\,\forall i$ then we recover the unbounded metric-free model of Ginelli and Chaté [25].

[TABLE]

In our model, which uses topological shell depth, the value of $f_{i}$ is the same for all particles in the same shell and can therefore be mapped to a lower dimensional parameter set $f_{\kappa(i)}$. We believe that this generalisation of the SMF model is natural, allowing us to describe the motional bias, not as a specific characteristic for a subset of birds, but as a rule for all birds that has a strength that depends on the relative depth of an individual in the swarm.

3 Methodology

We are interested in measuring the density variation across our simulated swarms. As we wish to compare directly to the empirical study of Starling murmurations [22] we seek to compute this in a similar fashion. The type of flocks which were studied in [22] were non-columnar and compact, with sharp borders, containing on the order of hundreds to thousands of birds, and which moved nearly linearly for sufficiently long times so as to treat their behaviour as near steady-state. The type of density variation we are interested in here is the density profile across flocks in this steady state, which is observed to be higher near the edge and to decrease toward the centre: It is not the propagating density waves observed in response to specific events, such as turning or shock.

We determine the spatial extent of simulated swarms using the $\alpha$-shape method [22, 33], which allows for the presence of concavities within the swarm to the scale of $\alpha$. To measure density, individuals with distance less than $\delta$ from the border were removed and a new border of the reduced flock was computed. The reduced density was computed using this reduced volume and the number of internal birds. This process was repeated until the flock was empty (i.e. less than four members remaining such that no tetrahedra, and hence no volume, can be determined).

Simulated swarms typically have a non-negligible degree of concavity (as is also observed in the empirical study), therefore allowing for presence of a non-convex border is natural. Fixing the convexity scale $\alpha$ is non-trivial as we are not dealing with a few observations, but thousands of configurational snapshots per simulation, therefore we cannot do this manually (as is described in [34]). Instead we obtain a sensible estimate for $\alpha$ by selecting the smallest value possible that leaves the particle aggregation as a single connected component. This fixes the convexity scale throughout. We must also make a choice of the flock reduction parameter $\delta$ as this impacts on our measurement and ability to compare with the data. We select a value which on average provides a similar number of flock reduction iterations as the field study (which is 7).

In order to prevent this choice from impacting our measurements we scale the reduction so that shell number is mapped to the domain $[0,1]$ with [math] corresponding to the first reduction and $1$ the final reduced flock. This also allows for a much easier comparison with the observational data; we can map that data to the same domain and perform cubic splines interpolation to allow query of comparison points between simulated and empirical data. Additionally we normalise the density data such that the first flock density measurement is $1$, which makes our measurements and comparisons dimensionless, and allows us to look primarily at the density gradient across the aggregation. These transformations allow us to compare our simulation data more easily with the empirical data and minimises the impact of possible differences in choice of parameters.

Our primary goal is to identify a bounding function $f_{\kappa}$ that can produce simulations with density profiles that provide a good fit to the empirical data. There is some freedom in how one might parametrise $f_{\kappa}$. We choose $f_{\kappa}$ to be linear in shell depth (parametrised via gradient $a$ and intercept $b$). We allow the bounding strength on shell 0 individuals to be a separate parameter $\phi_{e}$ in order to include models in which individuals on the edge behave differently from the bulk.

We then use the Simultaneous Perturbation Stochastic Approximation (SPSA) algorithm [35, 36] for recursive optimisation of bounding function parameters ($\phi_{e}$, $a$, $b$) using gain sequences with suggested practical values from [37]. We used the mean-squared difference between simulated and empirical data, averaged over a specified number of density evaluations, as the cost function estimate. Using this method allows for a principled stochastic search of the parameter space and can be performed in parallel. Fresh simulations were performed at each parameter update, due to the presence of hysteresis in these types of systems [32, 38].

4 Results

In order to understand how the density across aggregation varies for swarms which interact in a metric-free fashion we generated simulations of our Distributed Motional Bias Strictly Metric-Free (DMBSMF) model, as described above. As we are interested in simulating real-world behaviour we choose the parameters for the model via stochastic optimisation using the previously described method, directly fitting to empirical data, obtaining fit parameters of $\phi_{e}=0.883$, $a=-0.944$, $b=0.056$. These parameters result in a bounding function $f_{\kappa}$ as displayed in figure 2. This translates to a strong surface effect generally pointing toward the centre of the flock, however the bulk of the flock has an outward motional bias of increasing strength as one approaches the centre.

In order to simulate a flock that is comparable to that observed in the field study, we note from the motivating empirical study [22] that the flock in question contains $1,360$ reconstructed birds. We also note details from a later study [24](S.I.) for the flock in question: $1,571$ reconstructed birds with a measured polarisation of $0.96\pm 0.03$ (i.e. observed flocks in high order regime). We therefore chose to simulate $1,500$ birds with noise parameter $\phi_{n}=0.22$, yielding a polarisation of $0.931\pm 0.003$, which is of similar magnitude to the observed flock. In each instance, we performed a simulation for $20,000$ time steps with the first $10,000$ steps discarded for equilibration of the system. The initial condition is a random (isotropic) orientation and a random location, uniformly distributed within a unit cube, for each individual. We measured the density variation across the flock (as described in section 3) every $10$ time steps after equilibration, resulting in $1,000$ measurements per simulation instance, which are then time-averaged. We combine the results from five independent simulation instances, with final values presented as the mean of these quantities and uncertainties corresponding to standard errors.

The simulated model matches closely to empirical data of Starling flocks, as can be seen in figure 3, and produces the observed effect that aggregation density is greater at the border and reduces in what appears to be a linear fashion. The rate of this decrease is also closely matched. This counter-intuitive observation appears to require a model with a surprising motional bias: whilst surface birds move toward the flock centre, ensuring global cohesion, the rest of the flock move toward the border with increasing strength the further from it they are, as determined by topological depth. Naturally then, the number of birds closer to the border of the flock increases and drops off toward the centre due to the strong gradient of the bulk bounding function.

Our model shares some similarity with another recently proposed flocking model, the “hybrid projection” model [39], that drives individuals to move towards features in their visual field, specifically the boundaries between light and dark regions, where light/dark encodes the absence/presence of a neighbour in each direction. This model effectively encourages the movement inwards of individuals near the flock border. This is because individuals at the border will experience featureless outward-directed visual fields, resulting in an inward bias. It will also generate a bias outwards from the bulk of the flock as there will typically be more features in the outward-pointing directions than toward the often opaque centre of the flock. It is notable then that the motional bias that fits data from real-world flocks is similar to the effective motional bias present in visual models of this type.

5 Determining topological depth

A key aspect of our model is the notion of topological depth within the flock. Individuals are assigned a shell number based on this quantity, encoding a non-metric measure of depth as the shortest path length from the individual to a member of the convex hull (shell 0). The motional bias experienced by this individual is a function of shell number, as shown in figure 2. It is therefore important to consider the accessibility of this quantity to the individual, from a biological/ sensory perspective - how might flock members determine their shell number? In this section we present a model for how this could be achieved using the degree of anisotropy in an individuals visual field as an indicator of their depth within the flock.

We analyse a simplified model of the system in which the density is homogenous, for simplicity. Consider the three-dimensional flock as a sphere $S$ of radius $R$ centred on the origin with particle mass distributed uniformly within this sphere. For a point $P$ on or inside the sphere we can define an axis $z$ along the vector from $P$ to the sphere centre at the origin, as seen in figure 4. In spherical polar coordinates $(r,\theta,\varphi)$ this necessarily has $\varphi$-rotational symmetry about the $z$ axis.

The number of particles $N=\int_{S}\rho(\underline{r})\,dV$ constrains the density, here assumed homogenous $\rho(\underline{r})=\rho$. If we transform to the frame in which $P$ as the origin, we can write:

[TABLE]

where $\tilde{r}$ is the radial component of a point in this frame and $d\Omega$ is the solid angle. Therefore,

[TABLE]

which is the particle mass per solid angle, where $L(\theta)$ is the distance from $P$ to the sphere surface. This quantity $I(\theta)$ is biologically accessible (i.e. can be sensed) via the visual field of an individual within the flock and is closely related (via a threshold function) to the fraction of sky occluded by individuals in the $\theta$ direction as observed from $P$.

For an individual at $P$ there are intuitively directions which have higher and lower particle mass per solid angle. The imprint of the flock on an individuals visual field is greater when looking through its centre than in the opposite direction, as can be seen in figure 4, panels A & B.

We are interested in the extrema of $I(\theta)$ and make use of the observation that $L(\theta)$ is the radial distance to the flock edge, see figure 4, with $P$ as the origin. This has the form $L(\theta)=r_{p}\cos\theta+\sqrt{R^{2}-r_{p}^{2}\sin^{2}\theta}$. To obtain the extrema of $I(\theta)$ we differentiate equation 6 which yields:

[TABLE]

For non-zero density $\rho$, there are a number of stationary points. First when $r_{p}=0$, from the perspective of an individual at the centre of the spherical flock, there is no variation in mass density in any direction, and $I(\theta)=\frac{\rho R^{3}}{3}$ is independent of $\theta$. More significantly, there is a maximum and minimum at $\theta=0$ and $\pi$ respectively. These correspond to $L(0)=R+r_{p}$, looking along a line from $P$ through the centre of the sphere (along $+z$), and $L(\pi)=R-r_{p}$, away from it (along $-z$). This also provides two features identifiable in the visual field of the individual at $P$: $I_{\max}=\frac{\rho(R+r_{p})^{3}}{3}$ and $I_{\min}=\frac{\rho(R-r_{p})^{3}}{3}$.

To obtain a quantity which captures the asymmetry of any individual’s visual field we take the ratio of the values of these two features to define the “visual anisotropy” as:

[TABLE]

where $D=r_{p}/R$ is the relative depth within the flock for an individual at $P$. Note how equation 6 does not explicitly feature $\rho$ and is “scale-free” by nature, being only a function of the dimensionless depth $D$, and is also monotonic on the interval $D\in[0,1]$.

This is useful if it can be linked to topological depth. To make this connection we now seek a relationship between relative depth $D$ and topological depth $\kappa$. For each time-step for our simulated, non-spherical flocks (example configurations can be seen in the Supplemental Materials), we determine the spatial extent of the flock as $R=\langle|\underline{r}_{i}-\underline{r}_{cm}|\rangle_{i\in C_{0}}$ the mean distance to centre of mass $\underline{r}_{cm}$ over all particles on the convex hull of the point set. Relative depth per individual is then determined as $D_{i}=|\underline{r}_{i}-\underline{r}_{cm}|/R$ which, on average, is one for individuals with zero topological depth. Figure 5 shows relative depth averaged over a thousand configurations from five simulations with parameters determined from the fit to empirical data, as shown in figure 2, compared with the corresponding topological depth. We observe a linear relationship with relative depth decreasing with increased topological depth: when an individual is closer to the centre ($|\underline{r}_{i}-\underline{r}_{cm}|$ is smaller) it has a higher topological depth and vice versa.

We can finally relate our biologically accessible quantity, the visual anisotropy $\Delta I$, to topological depth $\kappa$ and we show this for our model in Figure 6, providing a one-to-one map. An individual can therefore compare two features (the minimum and maximum projected density) from their visual field in order to determine their topological depth within the flock, and hence understand how they should adjust their motion. One could imagine such a relationship might be determined heuristically: an intuitive understanding of depth within the aggregation from visual observations. When $\Delta I$ is small, the ratio between minimum and maximum of particle mass per solid angle $I(\theta)$ is large, so there is a large distinction between the two directions these represent (away from and toward the bulk of the flock respectively). When $\Delta I$ is larger the curve has less extreme slope and presents distinct values for different topological depths suggesting an individual deep in the flock still has capacity to determine its depth. We intend to further develop this model, including the role of heterogeneity, in future work [40].

6 Conclusion

In conclusion, we have introduced a generalised topological model of collective behaviour with a tunable bounding function to distribute a metric-free motional bias across the swarm. This model was fitted to empirical data of Starling murmurations using stochastic optimisation to determine a suitable form of bounding function. Simulation data from this model was shown to match field study data and produce swarms which are more dense at the border than at the centre, which is a surprising characteristic of real-world Starling flocks. We compared this fitted model to a benchmark topological model with no motional bias on the bulk of the flock (only on the surface). This allowed us to understand the role of the specific form of distributed motional bias that we have identified, which is to produce the desired level of inter-individual exclusion across the swarm, and allow individuals to keep the necessary relative distance apart without directly enforcing what this should be. We also proposed how an individual might use the observed anisotropy of its visual field to determine its depth within the flock.

Models of swarming generally aim to obtain group cohesion and coalignment [41, 42]. Typically, these are explicitly included as rules imposed on the interacting agents in the system. Our model differs from current models in the literature. While it explicitly imposes coalignment in a familiar way, swarm cohesion (and density regulation) are controlled using a motional bias distributed across the flock, which is prescribed via metric-free interaction rules, consistent with experimental observations. We show that specific field observations of density variation in aggregations of Starlings can be reproduced using our model so that density is higher on the border of the flock than at the centre. This density profile may relate to the predator-evasion mechanisms of three-dimensional swarms and the evolutionary development of such behaviour.

This work was funded by the UK Engineering and Physical Sciences Research Council (EPSRC) through the Complexity Science Doctoral Training Centre (JML), Grant EP/E501311 and by EPSRC grant # EP/E501311/1, a Leadership fellowship (MST). Computing facilities were provided by the Centre for Scientific Computing of the University of Warwick with support from the Science Research Investment Fund. We also acknowledge the use of the Computational Geometry Algorithms Library (CGAL) [43]. The authors would like to thank the referees for encouraging us to think more about how topological depth might be sensed.