## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Statistical mechanics for natural flocks of birds

Edited by Boris I. Shraiman, University of California, Santa Barbara, CA, and approved January 28, 2012 (received for review November 18, 2011)

## Abstract

Flocking is a typical example of emergent collective behavior, where interactions between individuals produce collective patterns on the large scale. Here we show how a quantitative microscopic theory for directional ordering in a flock can be derived directly from field data. We construct the minimally structured (maximum entropy) model consistent with experimental correlations in large flocks of starlings. The maximum entropy model shows that local, pairwise interactions between birds are sufficient to correctly predict the propagation of order throughout entire flocks of starlings, with no free parameters. We also find that the number of interacting neighbors is independent of flock density, confirming that interactions are ruled by topological rather than metric distance. Finally, by comparing flocks of different sizes, the model correctly accounts for the observed scale invariance of long-range correlations among the fluctuations in flight direction.

The collective behavior of large groups of animals is an imposing natural phenomenon, very hard to cast into a systematic theory (1). Physicists have long hoped that such collective behaviors in biological systems could be understood in the same way as we understand collective behavior in physics, where statistical mechanics provides a bridge between microscopic rules and macroscopic phenomena (2, 3). A natural test case for this approach is the emergence of order in a flock of birds: Out of a network of distributed interactions among the individuals, the entire flock spontaneously chooses a unique direction in which to fly (4), much as local interactions among individual spins in a ferromagnet lead to a spontaneous magnetization of the system as a whole (5). Despite detailed development of these ideas (6⇓⇓–9), there still is a gap between theory and experiment. Here we show how to bridge this gap by constructing a maximum entropy model (10) based on field data of large flocks of starlings (11⇓–13). We use this framework to show that the effective interactions among birds are local and that the number of interacting neighbors is independent of flock density, confirming that interactions are ruled by topological rather than metric distance (14). The statistical mechanics models that we derive in this way provide an essentially complete, parameter-free theory for the propagation of directional order throughout the flock.

We consider flocks of European starlings, *Sturnus vulgaris*, as in Fig. 1*A*. At any given instant of time, following refs. 11⇓–13, we can attach to each bird *i* a vector velocity and define the normalized velocity (Fig. 1*B*). On the hypothesis that flocks have statistically stationary states, we can think of all these normalized velocities as being drawn (jointly) from a probability distribution . It is not possible to infer this full distribution directly from experiments, because the space of states specified by is too large. However, what we can measure from field data is the matrix of correlations between the normalized velocities, . There are infinitely many distributions that are consistent with the measured correlations, but out of all these distributions, there is one that has minimal structure: It describes a system that is as random as it can be while still matching the experimental data. This distribution is the one with maximum entropy (10).

It should be emphasized that the maximum entropy principle is not a “modeling assumption;” rather it is the absence of assumptions. Any other model that accounts for the observed correlations will have more structure and hence (explicitly or implicitly) assumes something about the nature of the interactions in the flock beyond what is required to match the data. Of course the fact that the maximum entropy model is minimally structured does not make it correct. It could be, for example, that individual birds set their flight direction by computing a complicated nonlinear combination of the velocities from multiple neighbors, in which case correlations among pairs of birds would be insufficient to capture the essence of the ordering mechanism. We view the maximum entropy model as a powerful starting point, from which, as we will see, we can generate detailed and testable predictions.

The maximum entropy distribution consistent with the directional correlations *C*_{ij} is [1]where *Z*({*J*_{ij}}) is the appropriate normalization factor, or partition function; the derivation follows ref. 10, as explained in the *SI Appendix*. Notice that there is one parameter *J*_{ij} corresponding to each measured element *C*_{ij} of the correlation matrix. To finish the construction of the model, we have to adjust the values of the *J*_{ij} to match the experimentally observed *C*_{ij}, [2]where the symbol indicates an average using distribution *P* from Eq. **1**, whereas indicates an average over many experiments. This matching condition is equivalent to maximizing the likelihood that the model in Eq. **1** will generate the data from which the correlations were computed.

The probability distribution in Eq. **1** is mathematically identical to a model that is familiar from the physics of magnets—the Heisenberg model (5)—in which a collection of spins interact so that their energy (or Hamiltonian) is ; Eq. **1** then describes the thermal equilibrium or Boltzmann distribution at a temperature *k*_{B}*T* = 1. In this context, the constants *J*_{ij} are the strength of interaction between spins *i* and *j*, where *J* > 0 means that these elements tend to align. For many physical systems, once we know the Hamiltonian there is a plausible dynamics that allows the system to relax toward equilibrium, which is the Langevin dynamics [3]where is an independent white noise “force” driving each separate degree of freedom. Finding trajectories that solve Eq. **3** produces samples that are drawn out of the probability distribution in Eq. **1**. The interesting point is that this kind of dynamical model also is well known in biology: The direction of motion of an individual evolves in time according to “social forces” reflecting a weighted sum of inputs from neighboring individuals, plus noise (4). In this framework, *J*_{ij} measures the strength of the force that tries to align the velocity of bird *i* along the direction defined by bird *j*. We emphasize that we are doing an analogy but a mathematical equivalence.

In contrast to most networks, the connectivity in a flock of birds is intrinsically dynamic—birds move and change their neighbors. Thus, it may not make sense to talk about matrix of correlations *C*_{ij} or interactions *J*_{ij} between labeled individuals. On the other hand, the continuous and rapid change of neighbors induced by motion implies that the interaction *J*_{ij} between bird *i* and bird *j* cannot depend directly on their specific identities but only on some function of their relative positions.

The simplest form of interaction that is independent of the birds’ identity is one in which each bird interacts with the same strength, *J*, with the same number of neighbors, *n*_{c} (or with all birds within the same radius *r*_{c}; see below). If the interactions are of this form, then Eq. **1** simplifies to [4]where means that bird *j* belongs to the first *n*_{c} nearest neighbors of *i*. It is important to note that Eq. **4** can be also derived, without any assumption, as a maximum entropy distribution consistent with an experimental quantity simpler than the full correlation matrix. Instead of the correlation *C*_{ij} for each pair of individuals, we can measure its average value over birds within a neighborhood of size *n*_{c}, i.e., [5]It can be easily shown (see *SI Appendix*) that the maximum entropy distribution consistent with this scalar correlation is in fact the distribution [**4**]. As in the more general problem, finding the values of *J* and *n*_{c} that reproduce the observed correlation *C*_{int} is the same as maximizing the probability, or likelihood, that model Eq. **4** generates the observed configuration of flight directions in a single snapshot. Biologically, Eqs. **4** and **5** encapsulate the concept that the fundamental correlations are between birds and their directly interacting neighbors; all more distant correlations should be derivable from these correlations. If this concept is correct, a model as in Eq. **4** that appropriately reproduces the fundamental correlations up to the scale *n*_{c} must be able to describe correlations on all length scales.

Importantly, with large flocks we can estimate the correlations among interacting neighbors from a single snapshot of the birds’ flight directions , as indicated in the second step of Eq. **5**. In contrast, if we were trying to estimate the entire correlation matrix in Eq. **2**, we would need as many samples as we have birds in the flock (see *SI Appendix*), and we would have to treat explicitly the dynamic rearrangements of the interaction network during flight. This situation is an extreme version of the general observation that the sampling problems involved in the construction of maximum entropy models can be greatly reduced if we have prior expectations that constrain the structure of the interaction matrix (15, 16).

## Results

We now apply this analysis to data on real flocks of starlings. Given a snapshot of the flock, we have the configuration , and we need to evaluate the probability in Eq. **4** for any value of *J* and *n*_{c}, then maximize this probability with respect to these parameters (see *Materials and Methods* and *SI Appendix* for details of the computation). Special care must be devoted to birds on the outer edge, or border, of the flock, because these individuals have very asymmetric neighborhoods and may receive inputs from signals outside the flock. If we take the flight directions of these border birds as given, we can study how information propagates through the flock, without having to make assumptions about how the boundary is different from the interior. Technically, then, we describe the flock with Eq. **4** but with the flight directions of the border birds fixed (again, see *Materials and Methods* and *SI Appendix* for details).

### Inferring the Interaction Parameters from Data.

We proceed as follows. For a single flock, at a given instant of time, we compute the correlation *C*_{int} predicted by the model in Eq. **4** as a function of the coupling strength *J* and compare it with the experimental value of the correlation (Fig. 2*A*). The equation fixes *J*(*n*_{c}) for each value of *n*_{c}. Then we fix the interaction range by looking at the overall probability of the data as a function of *n*_{c}. In general there is a clear optimum (Fig. 2*B*), from which we finally infer the maximum entropy value of both parameters, *n*_{c} and *J*. We repeat this procedure for every snapshot of each flock and compute the mean and standard deviation of the interaction parameters for each flock over time. Alternatively, for a given flock we can average the log–likelihood over many snapshots, and then optimize, and this procedure gives equivalent results for *J* and *n*_{c} (see *SI Appendix*, section V).

In Fig. 2 *C* and *D* we report the value of the interaction strength *J* and of the interaction range *n*_{c} for all flocks, as a function of the flock’s spatial size, *L*. The inferred values of *J* and *n*_{c} are reproducible, although error bars are larger for smaller flocks. In particular, *J* and *n*_{c} do not show any significant trend with the flocks’ linear dimensions, with the number of birds, or with the density. This result is not obvious, nor is it in any way built in to our framework; for example, if the real interactions extended over long distances, then our fitting procedure would produce an increase of *n*_{c} and *J* with the size of the flock.

In Fig. 2*E* we also show that the interaction range *n*_{c} does not depend on the typical distance between neighboring birds, *r*_{1}, which is closely related to the flock’s density. Of course, we can run exactly the same method using a metric interaction range, *r*_{c}, rather than a topological range, *n*_{c}. We simply set *J*_{ij} = *J* if and only if birds *i* and *j* lie within *r*_{c} meters. In this way we find that the metric range *r*_{c} does depend on the nearest neighbor distance *r*_{1}, in contrast with the topological range *n*_{c} (Fig. 2*F*). This result provides strong support for the claim put forward in ref. 14 that birds interact with a fixed number of neighbors, rather than with all the birds within a fixed metric distance.

### Model Predictions.

Having fixed *J* and *n*_{c} by matching the scalar correlation in the flock, we have no free parameters—everything that we calculate now is a parameter-free prediction. We begin by computing the correlations between pairs of birds as a function of their distance, , as shown in Fig. 3*A*. There is extremely good agreement across the full range of distances. As we have seen, our maximum entropy calculation finds local interactions, i.e., a relatively small value of *n*_{c} (*n*_{c} ∼20 for flocks of up to thousands birds). This result implies that the scalar correlation *C*_{int}, used as an experimental input to the calculation, is the integral of *C*(*r*) only over a very small interval close to *r* = 0: Only the average value of pair correlations at very short distances is used as an input to the calculation, whereas all the long range part of *C*(*r*) is not. Nevertheless, we have very good agreement out to the overall extent of the flock itself. This finding confirms our expectation that a model for the local correlations is able to describe correlations on all length scales. We draw attention to the fact that the apparent correlation length, defined by *C*(*r* = *ξ*) = 0, is predicted to scale with the linear size of the flock (*ξ* ∝ *L*, Fig. 3*A, Inset*), as observed experimentally (17).

Correlations exist not just between pairs of birds, but among larger *n*-tuplets. In Fig. 3*B* we consider the correlations among quadruplets of birds. Although these correlations are small, their shape is nontrivial and quite noiseless. The model, which takes only local pairwise correlations as input, reproduces very accurately these four-body correlations, including a nonmonotonic dependence on distance, out to distances comparable to the full extent of the flock. Again, we are not making a fit but a parameter-free prediction.

Finally, instead of measuring correlations, we can ask the model to predict the actual flight directions of individual birds in the interior of the flock, given the directions chosen by birds on the border. This prediction can’t work perfectly, because within the model individual birds have an element of randomness in their choice of direction relative to their neighbors, but as shown in Fig. 3*D* the overlap between predicted and observed directions is very good, not just for birds close to the border but throughout the entire “thickness” of the flock. This result shows that the inference procedure and the model predictions work remarkably well for all individuals in the group. For a discussion of variability across different snapshots and flocks see *SI Appendix*.

### Testing the Mechanistic Interpretation.

The maximum entropy model has a mechanistic interpretation, from Eq. **3**, in terms of social forces driving the alignment of the flight directions. Given the success of the model in predicting the propagation of order throughout the flock, it is interesting to ask whether we can take this mechanistic interpretation seriously. As a test, we have simulated a population of self-propelled particles in three dimensions moving according to social forces that tend to align each particle with the average direction of its neighbors, as described by Eqs. **9** and **10** in *Materials and Methods*. We then compared the simulation parameters to the values obtained by applying the maximum entropy method to snapshots drawn from the simulation, just as we have analyzed the real data. Both the strength and the range of the interaction given by the maximum entropy analysis are proportional to the “microscopic” parameters used in the simulation (Fig. 4 *A* and *B*), although the maximum entropy interaction range is roughly 3× larger than the true number of interacting neighbors, . We believe that this overestimation is due to the fact that birds (unlike spins) move through the flock, encountering new neighbors before losing memory of the earlier flight directions and in so doing propagate information and correlation more effectively than if they were sitting on a fixed network. In other words, the maximum entropy model, where interactions are static by construction, compensates the dynamical nature of the true interaction network by giving a larger effective value of *n*_{c}. Hydrodynamic theories of flocking (7, 8) provide an analytic treatment of this effect, which is essential for collective motion of large two-dimensional groups. Indeed, in the limit of very large flocks, this ratio between the microscopic range of interactions and the effective range recovered by maximum entropy methods is predicted by the hydrodynamic theory (7, 8) to become arbitrarily large, but the flocks we study here seem not to be big enough for this effect to take over. If we use the “calibration” of the model from Fig. 4*B*, then the observation of *n*_{c} = 21.6 in the real flocks (Fig. 2) suggests that the true interactions extend over *n*_{c} = 7.8, in reasonable agreement with the result from (14, 18), *n*_{c} = 7.0 ± 0.6, using very different methods.

## Discussion

To summarize, we have constructed the minimal model that is consistent with a single number characterizing the interactions among birds in a flock, the average correlation between the flight directions of immediate neighbors. Perhaps surprisingly, this provides an essentially complete theory for the propagation of directional order throughout the flock, with no free parameters. The theory predicts major qualitative effects, such as the presence of long range, scale-free correlations among pairs of birds, as well as smaller, detailed effects such as the nonmonotonic distance dependence of (four-point) correlations among two pairs of birds. The structure of the model corresponds to pairwise interactions with a fixed number of (topological) neighbors, rather than with all neighbors that fall within a certain (metric) distance; the relevant number of neighbors and the strength of the interaction are remarkably robust across multiple flocking events.

For a long time, theoretical studies of collective animal behavior have relied on arbitrary (albeit reasonable) modeling assumptions. Recently, thanks to more sophisticated experimental techniques, it has finally been made possible to directly fit the models’ parameters to experimental data and perform model selection (19⇓–21). Our work, though, differs from current model fitting approaches in two respects. First, as mentioned in the introduction, we use a general guiding principle (maximum entropy) that spares us from the uneasiness of making assumptions. Second, we consider large groups of individuals, up to several thousands starlings, to be compared with a maximum of few tens (with a typical size of four) analyzed in the lab. Hence, we derive a quantitative theory for the large scale behavior, where the collective phenomenon is fully blown. At the moment, we do not know to what extent the interactions inferred at the level of few individuals are linearly “scalable” up to the thousands, so it seems important to formulate a theory directly at the collective level.

Our approach can be seen as part of a larger effort using maximum entropy methods to link the collective behavior of real biological systems to theories grounded in statistical mechanics (22⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–34). As with these other examples, we view the success of our theory as an encouraging first step. We have focused on the flight directions, taking the positions of the birds as given. A full theory must connect the velocities of the birds to their evolving positions, which requires more accurate measurements of trajectories over time, and we must consider the fluctuations in the speed as well as direction of flight. There are maximum entropy approaches to both of these problems, and the mapping from maximum entropy models to statistical mechanics suggests that the observation of scale-free correlations in speed fluctuations (17) will locate models of flocking at an especially interesting point in their parameter space (35).

## Materials and Methods

### Data.

Analyzed data were obtained from experiments on large flocks of starlings (*Sturnus vulgaris*), in the field. Using stereometric photography and innovative computer vision techniques (12, 13) the individual 3D coordinates and velocities were measured in cohesive groups of up to 4,268 individuals (11, 14, 17). The dataset comprises 21 distinct flocking events, with sizes ranging from 122 to 4,268 individuals and linear extensions from 9.1 m to 85.7 m (see *SI Appendix*, Table S1 for details). Each event consists of up to 40 consecutive 3D configurations (individual positions and velocities), at time intervals of 1/10 s. All events correspond to very polarized flocks (polarization between 0.844 and 0.992; see *SI Appendix*, Table S1). Given a flock at a given time, its exterior border can be computed using appropriate triangulation algorithms (see ref. 13 and *SI Appendix*). Boundary birds are then defined as those who belong to the exterior border.

### Analytic Approach to the Maximum Entropy Model.

To apply the maximum entropy analysis, we need to compute the expected values of correlation functions according to the measure defined in Eq. **1**. To this end we must first compute the partition function *Z*({*J*_{ij}}), which is, in general, a hard task. Flocks, however, are very ordered groups, in that birds’ velocities are neatly aligned to each other (17). In this case we can use the “spin wave” approximation (36), which exploits the strong ordering condition. Let us call the global order parameter, or polarization, measuring the degree of collective alignment, where is a unit vector. Individual orientations can be rewritten in terms of a longitudinal and a perpendicular component: . If the system is highly polarized, *S* ∼1, , and ; we verify this last condition for our data in *SI Appendix*, Fig. S1. The partition function can be written as an integral over the , and if *S* ∼1 the leading terms are (see *SI Appendix*, section II, for details): [6]where , , and the satisfy the constraint . If we consider the flight directions of birds on the border as given, integration must be performed with respect to internal variables only (see *SI Appendix*, section IV). After some algebra one gets [7]Here and represent the subsets of, respectively, internal and border individuals; is a ‘field’ describing the influence of birds on the border on internal bird *i*; and, now, . The integral in Eq. **7** can be carried out explicitly; see *SI Appendix*, Eqs. S35–S40. The reduced model in Eq. **4** corresponds to *J*_{ij} = *Jn*_{ij}, with *n*_{ij} = 1, 1/2, or 0 according to whether both individuals, just one, or none, belong to the local *n*_{c}-neighborhood of the other. Given the individual coordinates of birds in space, the matrix *A*_{ij} can be computed for any given snapshot, and (and correlation functions) can be calculated as a function of *J* and *n*_{c}. These two parameters must then be adjusted to maximize the log–likelihood of the data, [8]Maximizing with respect to *J* corresponds to equating expected and experimental correlations. In our case, this equation can be solved analytically, leading to an explicit expression of the optimal *J* vs. *n*_{c}; see *SI Appendix*, Eq S46. Maximization with respect to *n*_{c} can then be performed numerically. A graphical visualization of the solution can be found in Fig. 2.

### Self-Propelled Particle Model.

We consider a model of self-propelled particles extensively studied in the literature (9). Each particle moves with vector velocity according to the following equations: [9][10]where Θ is a normalization operator that serves to keep the speed fixed at , and means that *j* belongs to the *n*_{c} interacting neighbors of *i*. The distance-dependent force acts along the direction connecting *i* and *j*; following ref. 9, if is the unit vector between *i* and *j*, we take [11][12][13]Finally, is a random unit vector, independent for each bird and at each moment of time. The parameters α and β tune the strength of the alignment and of the cohesion force, respectively; in particular, the strength of alignment is given by *J* = *v*_{0}*α*/*n*_{c}. To test the maximum entropy analysis, we modified the model in such a way that we could vary *n*_{c}. Specifically, we introduced an angular resolution *μ* such that only neighbors with mutual angles larger than *μ* were included in the neighborhood. When *μ* is of the order of the Voronoi angle the model is statistically equivalent to the original version (where Voronoi neighbors were considered), but increasing (decreasing) *μ* one can decrease (increase) the value of *n*_{c}. In this way both the number *n*_{c} of interacting neighbors and the strength of the interaction *J* can be arbitrarily tuned. Parameters were chosen as *r*_{0} = 1 (to set the scale of distance), *r*_{b} = 0.2, *r*_{e} = 0.5, *r*_{a} = 0.8, *α* = 35, *β* = 5, *v*_{0} = 0.05, and we simulated a flock of *N* = 512 birds. Additional simulations were run (see *SI Appendix*, section IX and Fig. S7) to check that the metric dependency of the distance-dependent force does not affect the relationship between inferred and real values of *J* and *n*_{c}.

## Acknowledgments

We thank I. Couzin, T. Grigera, N. Leonard, G. Parisi, G. Theraulaz, and J. Toner for helpful discussions. Work in Princeton was supported in part by National Science Foundation Grants IIS-0613435 and PHY-0957573; work in Rome was supported in part by Grants IIT-Seed Artswarm, ERC–StG n. 257126 and AFOSR-Z809101. Our collaboration was facilitated by the Initiative for the Theoretical Sciences at the Graduate Center of the City University of New York, supported in part by the Burroughs-Wellcome Fund.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: irene.giardina{at}roma1.infn.it.

Author contributions: W.B., A.C., I.G., T.M., and A.M.W. designed research; W.B., A.C., I.G., T.M., E.S., M.V., and A.M.W. performed research; A.C., I.G., and M.V. analyzed data; and W.B., A.C., I.G., T.M., and A.M.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1118633109/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- Camazine S,
- et al.

- ↵
- Anderson PW

- ↵
- Wilson KG

- ↵
- ↵
- Binney J,
- Dowrick NJ,
- Fisher AJ,
- Newman MEJ

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Ballerini M,
- et al.

- ↵
- ↵
- ↵
- Cavagna A,
- et al.

- ↵
- ↵
- Lukeman R,
- Li Y-X,
- Edelstein-Keshet L

- ↵
- Katz Y,
- Tunstrma K,
- Ioannou CC,
- Huepe C,
- Couzin ID

- ↵
- Herbert-Read JE,
- et al.

- ↵
- ↵
- Tkačik G,
- Schneidman E,
- Berry MJ II.,
- Bialek W

- ↵
- Shlens J,
- et al.

- ↵
- Lezon TR,
- Banavar JR,
- Cieplak M,
- Maritan A,
- Federoff NV

- ↵
- Tkačik G

- ↵
- Bialek W,
- Ranganathan R

- ↵
- Tang A,
- et al.

- ↵
- Yu S,
- Huang D,
- Singer W,
- Nikolić D

- ↵
- Tkačik G,
- Schneidman E,
- Berry ML II.,
- Bialek W

- ↵
- Weigt M,
- White RA,
- Szurmant H,
- Hoch JA,
- Hwa T

- ↵
- ↵
- Mora T,
- Walczak AM,
- Bialek W,
- Callan CG

- ↵
- ↵
- ↵