Local equilibrium in bird flocks.

The correlated motion of flocks is an instance of global order emerging from local interactions. An essential difference with analogous ferromagnetic systems is that flocks are active: animals move relative to each other, dynamically rearranging their interaction network. The effect of this off-equilibrium element is well studied theoretically, but its impact on actual biological groups deserves more experimental attention. Here, we introduce a novel dynamical inference technique, based on the principle of maximum entropy, which accodomates network rearrangements and overcomes the problem of slow experimental sampling rates. We use this method to infer the strength and range of alignment forces from data of starling flocks. We find that local bird alignment happens on a much faster timescale than neighbour rearrangement. Accordingly, equilibrium inference, which assumes a fixed interaction network, gives results consistent with dynamical inference. We conclude that bird orientations are in a state of local quasi-equilibrium over the interaction length scale, providing firm ground for the applicability of statistical physics in certain active systems.

Animal groups moving in concert such as mammal herds, fish schools, and bird flocks show that in biology, just as in physics, local coordination can result in large-scale order [1-3]. However flocks differ from classical statistical physics in that their constituents are active: they constantly move by self-propulsion, pumping energy into the system and keeping it out of equilibrium [4-7]. The key element is the rearrangement of the interaction network due to the active motion of individuals relative to each other, continuously changing their neighbours. Theoretical studies show that network rearrangement has major consequences, which include enhancing collective order, reducing from 3 to 2 the lower critical dimension, and affecting the critical exponents [4, 8].

However, the importance of activity must be assessed with respect to the relevant time scales of the system. The impact of network rearrangement depends on the interplay between its characteristic time scale, τnetwork, defined as the average time it takes an individual to change its interaction neighbours, and the local relaxation time scale, τrelax, defined as the time needed to relax locally the order parameter if the interaction network were fixed. If τnetwork ≤ τrelax, the interaction network rearranges at least as fast as the order parameter relaxes, and the system remains far from equilibrium. If on the other hand τrelax ≪ τnetwork, the relaxation of the order parameter is adiabatic, closely following the network as it slowly evolves. In this case, even though the system behaves in an out-of-equilibrium manner on the longest scales, it locally obeys a condition of equilibrium, and we expect some of the tools of equilibrium statistical physics to be applicable.

Here, we explicitly address the impact of network activity by developing a new inference method based on the exact integration of maximum-entropy dynamical equations, thus accounting for the reshuffling of the network. We apply the method to data of starling flocks of up to 600 individuals [9-12] (see Materials and Methods and Table S1 for data summary), inferring the relevant parameters of the interactions between individuals. We find that the alignment relaxation time, τrelax, is more than one order of magnitude shorter than the network rearrangement time, τnetwork. Consistently, we show that the parameters learned from the dynamics are consistent with those obtained by an equilibrium-like inference, which assumes a fixed network [13]. Our results suggest that natural flocks are in a state of local quasi-equilibrium over the interaction length scale, meaning that the relatively slow rearrangement of the local interaction network does not affect the ordering dynamics up to certain scales.

To compare the relevant time scales of the ordering process in flocks, we first need to learn the dynamical rules of their behaviour. Learning these rules usually relies on inferring the parameter of a chosen model directly from the data, as has been recently done in surf scoters [14] and fish schools [15-19]. Although in these studies the local rules of interaction were often learned using small groups, in some cases they could also be used to predict large-group behavior [17, 19]. Here, instead of assuming a model a priori, we apply the principle of maximum entropy to the trajectories of all birds in the group [20]. We look for a distribution of the stochastic process that is as random as possible, while agreeing with the data on a key set of experimental observables.

In a flock of size N, we call the three-dimensional flight orientation of bird i at time t. The maximum entropy distribution over possible flock trajectories that is consistent with the correlation functions , as well as their derivatives , can be exactly mapped, in the limit of strong polarization , onto the following stochastic differential equation (see SI and [Ref. 20]): where is a random white noise, and where the projection onto the plane perpendicular to ensures that remains of norm 1. Equation (1) can be viewed as a generalization of the Vicsek model [21]: each bird modifies its flight direction according to a weighted average of the directions of its neighbours. The interaction matrix J encodes how much bird i is influenced by (i.e. interacts with) bird j. Given the experimentally measured correlation functions, entropy maximization yields equations that fix the values of the noise amplitude and the interaction matrix J. This matrix has too many parameters to be reliably determined from the data, but we can reduce its complexity by parametrising it. It was shown in [22] that the interaction decays exponentially with the topological distance k between birds, where k denotes the (time-dependent) rank of bird j among the neighbours of bird i ranked by distance. This interaction matrix has just two parameters: n is the topological interaction range, while J is the overall strength of the interaction. The noise is uncorrelated among birds and of uniform magnitude T, by analogy with physical temperature: , where d is the space dimension (d = 3 in the following).

In principle, to learn the parameters of Eq. 1 one needs actual continuous-time derivatives. In practice, we only have configurations separated by the finite experimental sampling time dt. A common solution is to use Euler’s approximation: where is a normally distributed vector of variance 1 in each direction. The conditional likelihood of the data given the model, , can be written in Gaussian form after expanding Eq. (3) in the spin-wave approximation (see Materials and Methods). Maximising this likelihood yields values for the alignment parameters n, J and T (see Ref. [20] and SI).

Euler’s approximation is used by virtually all methods that try to fit a dynamical equation to a discrete time series [15-17]. However, it is inappropriate when the experimental sampling time, dt, is larger than the intrinsic relaxation timescale, τrelax. In this case information spreads between subsequent frames beyond the directly interacting neighbours and Euler’s approximation overestimates the range of the interaction, as we shall see below. To overcome this issue, we rewrite Eq. 1 by formally subtracting from it: Bold symbols denote vectors and matrices over bird indices; the matrix , where is the connectivity matrix (2). Λ is analogous to a Laplacian defined on a lattice, and obeys the sum rule: . In the spin-wave approximation, where all orientations point in almost the same direction, this relation ensures that has almost no contribution along the common direction of flight, implying (see Materials and Methods and SI). Equation 4 is now linear and it can be integrated exactly: This result assumes a constant J in the interval dt, which is a good approximation if dt ≪ τnetwork. Fortunately, this same condition is necessary for the very possibility to collect data: tracking requires to follow each individual across time, which is only possible if individuals do not significantly change their neighbourhood between consecutive frames. The integrated noise in the right-hand side of (5) is Gaussian, of mean zero and covariance . Using the exact solution (5) we can write an explicit expression for the (Gaussian) conditional likelihood , which can then be maximised over the parameters of the model (see Materials and Methods).

We first tested our dynamical inference method on synthetic data simulated using the model of Eq. 1, with τrelax ≈ 0.7, for various values of the interaction range n (see Materials and Methods). We infer the parameters of the model using either Euler’s rule or the result of exact integration, for different values of the sampling time ranging from dt = 0.2 to dt = 0.8. The method based on exact integration predicts the interaction range n well, regardless of dt (Fig. 1A and B), while the method based on Euler’s approximation largely overestimates n at large dt (Fig. 1B). We can now apply our dynamical inference to real flocks and learn the model parameters. First, we used data of natural flocks to check the effect of changing the sampling time dt, from the real sampling time of our setup, dt = 0.2 s (see Materials and Methods), to 0.8 s. Although we cannot compare the inferred value of n to the ground truth as in simulations, we observe a similar trend as a function of dt (Fig. 1C), with the exact integration and Euler’s approximation methods agreeing only at small dt. This suggests that the sampling time of 0.2 s is of the same order as the orientation relaxation time τrelax, as we will confirm below. It also indicates that the inference method based on exact integration is extracting the parameters of alignment reliably.

Fig. 1

Performance of the inference methods on the predicted interaction range n. A. Inferred versus real n obtained by applying our new inference method to simulated data generated with Eq. 1 at various interaction ranges. The method performs well for different values of the sampling rate dt. B. Dependence of the inferred n on the sampling time dt. On simulated data with n = 10 (dashed line), the inference method based on exact integration (red points) performs well regardless of the sampling time dt. By contrast, the inference method based on Euler’s integration method (green points) overestimates the true interaction range at large dt. C. A similar trend is observed when we apply the two inference procedures to real flocking data, as illustrated here on one flocking event. Note that in this case the true value is not known. Error bars represent standard errors over time frames.

Using the model parameters learned from the data, we evaluate the two time scales of interest for activity, namely relaxation of the orientations and network rearrangement. We estimated the network rearrangement time τnetwork experimentally for each flocking event as the characteristic decay time of its autocorrelation function , by fitting Cnetwork(t) ≈ C0 exp(−t/τnetwork) (Supplementary Fig. S1).

Working out the time scale of relaxation is more subtle. The relevant quantity is the product of the interaction strength J, which has units of inverse time, by the dimensionless connectivity matrix, Λ, as can be seen from Eq. (4). Since there are n neighbours acting on each individual, the total alignment force is of order Jn, suggesting that the characteristic time scale of relaxation of the orientations is τrelax ~ (Jn)−1. This result, however, seems at odds with the well-known fact that systems with spontaneously broken continuous symmetry - such as flocks - have correlation length and relaxation time that diverge with the system size L (Goldstone theorem [23]). On the other hand, we do not expect the large-scale modes responsible for this divergence to affect the local relaxation dynamics and its interplay with network reshuffling. To clarify this issue we calculate the dynamical autocorrelation function of the fluctuations of the order parameter, , where . We consider a fixed lattice, because we need to gauge relaxation in absence of network rearrangements, resulting in the autocorrelation function (see SI): where a is the lattice spacing. The infrared divergence at small k, which correspond to large-scale modes, makes the integral divergent in the L → ∞ limit for d = 2 (Mermin-Wagner theorem [24]). In d = 3 the integral is finite, but the correlation function is a power law, so that the relaxation time diverges with L (Goldstone theorem). The small k modes in (6) correspond to long wavelengths fluctuations spanning the entire flock, causing the local order parameter to relax slowly. However, these long wavelength fluctuations do not contribute to the disordering of the local interaction network: if the wavelength of a fluctuation is much larger than the interaction range, all directions of motion in the interaction neighbourhood fluctuate in unison, causing no change in the mutual positions of the birds. We conclude that the autocorrelation function that impacts on local network rearrangements only includes contributions from wavelengths up to the local interaction range (let us call it r). This amounts to restricting the integral in (6) to the modes , thus eliminating the infrared divergent modes k ~ 1/L. The resulting correlation function is exponentially decaying (see SI for the calculation of the integral), with finite relaxation time equal to τrelax = (Jn)−1, consistent with our initial guess. We note that, by considering wavelengths up to the interaction range, we are still dealing with a coarse-grained field theory, as in most biological systems the scale of interaction extends over tens of neighbours.

We can now proceed with the comparison of τnetwork and τrelax. Results are summarised in Fig. 2. The two time scales clearly separate, with local relaxation almost two orders of magnitude faster than network reshuffling. This separation of time scales suggests that flocks are in a state of local equilibrium. The network of interactions changes slowly enough for the dynamics of flight orientations to catch up before neighbours reshuffle. In other words, the orientation dynamics tracks network changes adiabatically. Note that this statement holds only locally, at the scale of the interaction range, as both τnetwork and τrelax are defined on that scale.

Fig. 2

Comparison between the two relevant time scales of active matter, as inferred in 14 natural flocks using our inference method based on exact integration. Histograms of the neighbour exchange time τnetwork versus the local alignment time τrelax = 1/Jn, show that the relaxation of orientations is much faster than the turnover of neighbours. Note that the experimental sampling time dt = 0.2 s (dashed line) is of the same order as the alignment time, justifying the use of exact integration. Inset: the scatter plot of τrelax versus τnetwork shows no correlation between the two quantities.

Since flocks behave as if they were in local equilibrium, an equilibrium inference procedure, which takes as input the local spatial correlation computed from a snapshot of the birds’ flight orientation [13], should be consistent with the results of the dynamical inference. To check this prediction, we recall the equilibrium-like inference method of [13]. For symmetric J, Eq. 1 is the Langevin equation derived from the Hamiltonian of the Heisenberg model When J varies slowly in time, the fluctuations of are in quasi-equilibrium and distributed according to Boltzmann’s law: We recognise the maximum entropy distribution consistent with the local correlation index fitted in Ref. [13]. In practice, the equilibrium inference consists in maximising the likelihood of Eq. 8 over its parameters n and J/T (see Materials and Methods and SI). If the variations of n are slow compared to the dynamics of , τnetwork ≫ τrelax, this inference procedure should give an accurate estimate of the alignment parameters. If however the two time scales are comparable, we expect the equilibrium inference to overestimate the true n, as the frequent exchange of neighbours results in an effective number of interaction partners that is larger than the instantaneous one. We verified both these expectations on simulated data, by showing that the equilibrium inference is accurate for τnetwork ~ 100τrelax, but overestimates n for τnetwork ~ τrelax (see Supplementary Fig. S2). When applied to empirical data, the dynamical and equilibrium inferences give consistent results, and predict the same interaction range, n, and coupling-to-noise ratio, J/T (Fig. 3) Note that, while the dynamical inference provides the strength of the interaction, J, and the strength of the noise, T, separately, the equilibrium inference only gives the ratio J/T, which is the quantity to compare. To better appreciate this result, recall that the two inference procedures are based on independent pieces of information: the equilibrium inference uses instantaneous orientations, while the dynamical inference exploits how these orientations change in time. Their agreement confirms that the alignment dynamics of flocks are in an effective state of equilibrium over the range n.

Fig. 3

Inference on natural flocks. For each of the 14 flocking events, the parameters of the model were inferred using either the dynamical inference method presented here, with dt = 0.2 s, or an equilibrium inference method as in [13]. A. Both methods agree well on the predicted value of the alignment range n. B. While the dynamical method infers the alignment strength J and the noise amplitude T separately, the equilibrium method only infers their ratio J/T, the value of which is consistent between the two methods. Error bars represent standard errors over time frames.

Theoretical studies of active matter indicate that out-of-equilibrium effects induced by the rearrangement of the interaction network play a major role in the ordering of the system [4, 5]. In this light, any attempt to understand the properties of active biological systems based on equilibrium approaches may seem inappropriate. Does it mean that we should we always relinquish the methods of equilibrium statistical mechanics when dealing with active systems? Our results address this question by showing that bird flocks are in a state of local equilibrium, due to the rapid relaxation of orientations compared to the slow rearrangement of the network, over the local scale of interaction. As a consequence, an equilibrium inference method, which assumes a fixed interaction network, gives equivalent results to a full dynamical treatment.

Equilibrium inference seems to be justified in this system, not only as a formal mathematical equivalence allowing for useful insights and predictions, but as a tool to extract bona fide biological parameters. The equilibrium approach is mathematically simpler and computationally less expensive than the dynamical one in the limit of strong polarisation, making it easier to analyse larger groups. Although a dynamical approach such as the one presented here is still necessary for extracting the precise relaxation timescale of the ordering mechanism, there may be more straightforward ways to evaluate its order of magnitude and get a quick assessment of the local equilibrium hypothesis.

Our results do not mean that natural flocks are in global equilibrium and that network rearrangements play no role. The interaction network, far from being fixed as if individuals were linked by springs [25], completely reshuffles on long time scales [26]. The directions of motion relax on a faster time scale than the network over the local scale of interaction, but the network does move on longer time scales, and over larger length scales, with important consequences. To appreciate this point we must stress again the difference between local, short-wavelength modes, which set the balance between relaxation and network rearrangement, and long-wavelength modes, which govern the long time and long distance correlations. Capturing these large-scale properties requires to describe the active fluid using a hydrodynamic approach [4]. Equilibrium inference works despite the existence of these large-scale modes because it only uses information at the local scale of interaction, where relaxation is fast.

The local equilibrium we have uncovered in natural flocks is not merely the consequence of the high degree of polarisation of this system. A high polarisation certainly implies slow network rearrangements, but it does not constrain the relaxation time, which could be even slower, as illustrated in our simulations (Supplementary Fig. S2). Conversely, there may be unpolarised systems where local relaxation is faster than network rearrangement – a limit easily obtained theoretically by considering weakly interacting, slowly moving individuals. Midge swarms may be such an example: they are not polarized, poised below the ordering transition [27], yet have been successfully analysed using standard equilibrium tools of critical phenomena [28]. In general, one must carefully quantify these two time scales to determine to what degree the tools of equilibrium statistical mechanics may be applied to a given active system.

Materials and Methods

Flocking data

The three-dimensional trajectories of all birds were reconstructed using imaging techniques. Stereoscopic experiments on natural flocks of European starlings were performed in the field in Rome using three high speed machine vision cameras shooting at 170 fps. The stereoscopic video acquisitions were then processed using a novel purpose-built three-dimensional tracking algorithm based on a recursive global optimization method [12]. This algorithm is extremely powerful, allowing for the reconstruction of full 3D trajectories of all individuals in groups of several hundreds individuals. We collected 3D data from 12 flocking events with sizes ranging from 50 to 600 individuals, and lasting from 2s to 6s (for details on the experiments and the dataset see Table S1 and [10, 29]). To avoid interference from birds flapping, which occurs at frequency ≈ 10 Hz, we subsampled all the 3D sequences so that two snapshots are separated by dt′ = 0.1 s. The instantaneous flight orientations were estimated by . To avoid overlap between two subsequent evaluations of , we used dt = 2dt′ = 0.2 s. The lower sampling rates of Fig. 1C, were obtained by taking dt′ = 0.2, 0.3, and 0.4 s.

Simulated data

Data were simulated in three dimensions with the continuous Vicsek model of Eq. 1 with the interaction matrix of Eq. 2. The positions of individuals are updated according to , with υ0 = 1. The simulations were set in a 8 × 8 × 8 box with periodic boundary conditions, and N = 512 birds, so that density is exactly 1. We set to obtain a polarization P ≈ 0.99 similar to natural flocks. Eq. 1 was integrated using Euler’s method with a simulation step dtsim = 0.01 that is much smaller than any other time scale in the system. The interaction range n varied from 7 to 25, and the interaction strength was picked so that Jn = 1.5, hence τrelax = (Jn)−1 ~ 0.7. The flocks were first brought to a steady state before taking snapshots for analysis.

Spin-wave approximation

The polarization P quantifies the level of order in the system. When P ≈ 1, we can expand each around the common direction of flight . This expansion gives , with . At leading order in , Eq. 4 becomes with . Similarly, the equilibrium distribution (Eq. 8) can be expanded into Since this distribution is Gaussian, Z can be calculated analytically and reads: , where λ are the eigenvalues of the matrix Λ.

Maximum likelihood Inference

The equilibrium inference is performed by maximising the likelihood of the data given by Eq. 10 over the parameters n and (J/T) (see SI for detailed formulas).

The dynamical inference based on Euler’s rule is implemented by maximising the likelihood calculated from Euler’s formula (Eq. 3). This likelihood reads The dynamical inference based on exact integration uses Eq. 5, rewritten as , where is a zero-mean Gaussian vector of covariance . The conditional likelihood now reads Depending on whether one uses Euler’s or exact integration rules, Eq. 11 or 12 is maximised over J, T and n (see SI for detailed formulas).

In all three inference procedures, the parameters are learned for each time t. Then the median and the associated standard error are calculated for each flocking event.

Data Availability

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon request.