Starling ballet and circular statistics

The start of my PhD thesis was marked by the last few lockdowns imposed by the COVID-19 epidemic. Working from home, I spent many hours looking out of my window at a city I hadn’t yet discovered. As you can often read in the columns of La Dépêche, Toulouse was then a victim of the nuisance caused by starling colonies¹. Shameless language - at a time when mankind is methodically destroying biodiversity - but revealing of our inability to think about the habitability of our environment.

Figure 1: A flock of starlings, also called murmuration, above the Dôme de la Grave in Toulouse. Photo credits: Xavier de Fenoyl, La Dépêche.

Twice a day at sunrise and sunset, alerted by the concert of chirping birds, I watched in amazement as thousands of green-reflecting starlings burst into the city sky. These tiny birds repeated an unpredictable ballet, creating a succession of impressionistic tableaux whose vitality seemed to mock our poetically barren lifestyles. The starling murmuration hypnotized me, and I wanted to understand how these coherent, dynamic structures could emerge from the crowd.

Video: Starling murmuration captured by Marco Valk Nature & Wildlife Films in Geldermalsen, The Netherlands.

After some research, I learned that the coherence of starling flocks results from a series of simple decisions made at the individual level. Vicsek’s model² provides a simplified description of collective movements, based on the assumption that each individual #$k$ regularly updates its trajectory, taking as its new direction $\alpha_k$ the mean of the directions of its nearest neighbors:

$$\alpha_k = \langle \alpha_i \rangle_{i \neq k, | r_k - r_i |< r} + \varepsilon_k.$$

In this model, the mean direction is calculated from all individuals located within a given radius $r$. In practice, each starling defines its trajectory according to its six to seven nearest neighbors at most³. Finally, $\varepsilon_k$ is a random Gaussian noise that gives rise to the chaotic character of the collective motion.

Interested in observability issues, I decided to implement Vicsek’s model to see if the data assimilation algorithms I was developing for my thesis could describe these collective motions from a limited number of observations. However, my first simulations showed a behavior not predicted by Vicsek’s model: regardless of their initial trajectories, my starlings ended up heading west, as if attracted by an invisible force (Fig. 2).

Figure 2: First simulation of a starling murmuration with Vicsek’s model.

I tore my hair out over this strange behavior for a long time, questioning each line of code, function after function, parameter after parameter. Eventually, I realized that the problem was laying in the averaging used to update each starling’s trajectory from its neighbors. It turns out that the average calculation I’d learned in college, i.e. the arithmetic mean:

$$\langle \alpha \rangle_{artihmetic} = \frac{\alpha_1 + \alpha_2 + \dots + \alpha_{N}}{N}$$

the arithmetic mean is not suitable for angles. Take, for example, the angles $\alpha_1=\frac{\pi}{2}$ and $\alpha_2=\frac{11\pi}{6}$, their arithmetic mean is $\langle \alpha \rangle_{artihmetic} = \frac12 \times \frac{3\pi + 11\pi}{6} = \frac{7\pi}{6}$ and is indeed the only equidistant point between $\alpha_1$ and $\alpha_2$ on the real number axis (Fig. 3a). However, if we look at the trigonometric circle this averaged direction $\langle \alpha \rangle_{artihmetic} = \frac{7\pi}{6}$ is the opposite of the directions defined by $\alpha_1$ and $\alpha_2$ (Fig. 3b).

Figure 3: Example of arithmetic mean (a) and circular mean (b) of two angles.

Why doesn’t it work? The arithmetic mean applies to scalar numbers defined in $\mathbb{R}$, but the angles are measures of directions defined on the unit circle, which is periodic. This implies that the angles $\alpha_1$ and $\alpha_2$ are actually closer when viewed on the unit circle rather than on the real number axis (Fig. 3b). It also implies that there is another angle equidistant from $\alpha_1$ and $\alpha_2$ (Fig. 3b). This angle is oriented on the same side of the directions defined by $\alpha_1$ and $\alpha_2$ and therefore more in line with our intuition of what the average direction.

The starlings had shown me the way to directional statistics⁴, a theory born out of the need to adapt the fundamental concepts of statistics to cyclic quantities such as angles, or seasons. In circular statistics, the mean can be defined as the angle of the vector resulting from the sum of the $N$ unit vectors of directions $\alpha_k$:

$$\langle \alpha \rangle_{circular} = \mathrm{atan2} \left(\frac{\sin{\alpha_1} + \dots + \sin(\alpha_N)}{N}, \frac{\cos{\alpha_1} + \dots + \cos(\alpha_N)}{N} \right).$$

Note that in the world of circular quantities, we have to forget certain classical properties of the mean. In particular, the circular mean:

1. is not linear,
2. does not allow the data dispersion to be described as a second-order moment, i.e. $\sigma^2(\alpha) \neq \langle \left( \alpha -\langle \alpha \rangle_{circular} \right)^2 \rangle$,
3. is not defined if the sum of the vectors of directions $\alpha_k$ is the zero vector, i.e. $\sum e^{i \alpha_k} =0$.

Coming back to our example, we find that the circular mean of $\alpha_1=\frac{\pi}{2}$ and $\alpha_2=\frac{11\pi}{6}$, is $\langle \alpha \rangle_{circular} = \frac{\pi}{6}$, which fits well with our intuition of an average direction, as shown by the green vector in Fig. 3b.

Once this circular mean was implemented, my numerical starlings learned to lose their compass again and now mimic the same hypnotic figures that I like to observe in the sky over Toulouse (Fig. 4).

Figure 4: New simulation of a murmuration of starlings with Vicsek’s model using the circular mean operator to define the trajectory of the starlings.