# Modeling collective motion: variations on the Vicsek model

Hugues Chaté, Francesco Ginelli, Guillaume Grégoire, Fernando Peruani,
Franck Raynaud

► To cite this version:

Hugues Chaté, Francesco Ginelli, Guillaume Grégoire, Fernando Peruani, Franck Raynaud. Modeling collective motion: variations on the Vicsek model. The European Physical Journal B: Condensed Matter and Complex Systems, 2008, 64, pp.451-456. 10.1140/epjb/e2008-00275-9. hal-01053505

HAL Id: hal-01053505

https://hal.science/hal-01053505

Submitted on 5 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling collective motion: variations on the Vicsek model

H. Chaté<fnref n="1" />,ª, F. Ginelli<fnref n="1" />,<fnref n="2" />, G. Grégoire<fnref n="3" />, F. Peruani<fnref n="1" />,<fnref n="2" />, and F. Raynaud<fnref n="1" />,<fnref n="3" />

Received 9 February 2008 / Received in final form 23 May 2008

Published online 11 July 2008 – © EDP Sciences, Società Italiana di Fisica, Springer-Verlag 2008

**Abstract.** We argue that the model introduced by Vicsek *et al*. in which self-propelled particles align locally with neighbors is, because of its simplicity, central to most studies of collective motion or “active” matter. After reviewing briefly its main properties, we show how it can be expanded into three main directions: changing the symmetry of the particles and/or of their interactions, adding local cohesion, and taking into account the fluid in which the particles move.

**PACS.** 64.70.qj Dynamics and criticality – 87.18.Nq Large-scale biological processes and integrative bio-physics

Collective motion is everywhere and at every scale, from herds of large mammals to amoeba and bacteria colonies, down to the cooperative behavior of molecular motors in the cell. The behavior of large fish schools and the dance of starling flocks at dusk are among the most spectacular examples. To the physicist, they are also highly non-trivial because they occur without any leader, external field, or geometrical constraint: collective motion can then be seen as the long-range orientational order following spontaneous symmetry-breaking.

Of course, collective motion has been studied in many fields long before statistical physicists got involved. Ethologists wonder about the signals exchanged between the moving animals, evolutionary biologists dissert on the benefit of moving in groups for individuals and for species. Robotics engineers strive to design robots which can accomplish a cooperative task without central control. Medical doctors try to understand tumor growth or wound healing, two situations in which cells move collectively. Because of these different motivations and viewpoints, the resulting modeling attempts come in many different kinds, have different goals, and often only aim at describing a particular situation in as much detail as possible, leading to over-parameterized models.

It is only recently (say fifteen years ago) that physicists, in their usual abrupt manner, approached the problem of collective motion by stripping it down to simple experiments or models having in mind the spontaneous symmetry-breaking picture mentioned above. On the experimental side, shaken asymmetric granular parti-

cles constitute the system of choice [1]. On the theoretical side, the simple model introduced by Vicsek and collaborators [2] —hereafter the “Vicsek model” (VM)— is central because of its “minimal” character. In short, in the VM, point particles move at fixed velocity, align locally with neighbors, while being submitted to some noise. In other words, an XY model in which the spins are actively moving.

In this paper, we briefly recall the essential properties of the Vicsek model, which have been comprehensively studied by us recently [3]. Then, we argue that the focal nature of the VM calls for it being extended, expanded, along three main directions in order to account for most of the different types of collective motion listed above. These are

- to vary the polarity of the particles and of their interaction

- to introduce an attraction/repulsion pairwise interaction to allow for cohesion

- to take into account the ambient fluid in which the particles move

## 1 The Vicsek model

Point particles move off-lattice in a space of dimension *d* with a velocity **v**ᵢ of fixed modulus *v*₀ = |**v**ᵢ| (driven-overdamped dynamics). The *direction* of motion of particle *i* depends on the average velocity of all particles (including *i*) in the spherical neighbourhood *S*ᵢ of radius *r*₀ centered

<footnote n="1">¹ CEA – Service de Physique de l’État Condensé, Centre d’Études de Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">² Institut des Systèmes Complexes Paris Île-de-France, 57-59 rue Lhomond, 75005 Paris, France</footnote>

<footnote n="3">³ Matière et Systèmes Complexes (MSC), UMR 7057, CNRS & Université Paris-Diderot, France</footnote>

<footnote>ª e-mail: hugues.chate@cea.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Main properties of the original Vicsek model ($d = 2, \rho = 2, v_0 = 0.5$, periodic boundary conditions). (a) Time-averaged order parameter $\langle\varphi\rangle_t$ vs noise strength $\eta$ for $L = 64$ (black circles) and $L = 256$ (red diamonds), time averages computed over $2 \times 10^7$ timesteps; Inset: corresponding Binder cumulants curves (a sharp minimum towards negative values is a signature of a first-order-like transition [3]). (b) Asymptotic phase diagram for the transition to collective motion; Inset: Log-log plot to compare the low density behavior with the mean field predicted behavior $\eta_c \sim \sqrt{\rho}$ (dashed red line). (c) Anomalous density fluctuations in the bandless regime (see text): $\Delta n$ scales approximately like $n^{0.8}$ ($L = 256, \eta = 0.25$) (The dashed line has slope 0.8). (d) Snapshot in a rectangular domain of size $1024 \times 256$ at $\eta = 0.42$. Only 50 000 particles are shown for clarity. The red arrow shows the mean direction of motion.</figcaption></figure>

on $i$. In a computer, at discrete timesteps $\Delta t$, we have:

$$ \mathbf{v}_i(t + \Delta t) = v_0 (\mathcal{R}_\eta \circ \vartheta) \left[ \sum_{j \in S_i} \mathbf{v}_j(t) \right] \quad (1) $$

where $\vartheta$ is a normalization operator ($\vartheta(\mathbf{w}) = \mathbf{w}/|\mathbf{w}|$) and $\mathcal{R}_\eta$ performs a random rotation uniformly distributed around the argument vector: in $d = 2$, for instance, $\mathcal{R}_\eta \mathbf{v}$ is uniformly distributed around $\mathbf{v}$ inside an arc of amplitude $2\pi \eta$. The positions $\mathbf{r}_i$ are then simply updated by streaming along the new direction as<fnref n="1" />

$$ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \Delta t \mathbf{v}_i(t + \Delta t). \quad (2) $$

Given the polar nature of each particle, the natural order parameter to monitor collective motion is just the (normalized) macroscopic velocity $\varphi(t) = \frac{1}{v_0} \langle \mathbf{v}_i(t) \rangle_i$. In most of the following, only its modulus $\varphi(t)$ will be considered.

The two main parameters of the VM are $\rho$, the density of particles, and $\eta$, the noise strength. At zero noise, perfect alignment eventually settles in the whole system (at least if the particles evolve in a domain with periodic

boundary conditions). At maximum noise ($\eta = 1$), particles are just non-interacting random walkers. Thus, a transition must occur in between.

Numerical results have shown that the onset of collective motion in the VM occurs at a finite noise level $\eta_c$, i.e. there exists a fluctuating ordered phase for $0 < \eta < \eta_c$. Extensive simulations going beyond the lengthscales considered by Vicsek et al. have shown that this transition is discontinuous (first-order like) (Fig. 1a). The transition line in the $(\rho, \eta)$ plane, follows, for small-enough $\rho$ values, the scaling law expected from a simple mean-field argument: $\eta_c \sim \rho^{1/d}$ (Fig. 1b).

The orientationally-ordered, fluctuating, collectively-moving phase possesses remarkable features. First of all, it shows true long-range order [4], even in two space dimensions, a remarkable departure from the equilibrium case of reference, the XY model, for which only quasi-long-range order arises.

In a large region of parameter space bordering the transition line, the density field is not homogeneous, but is organized in high-density high-order traveling objects spanning the dimensions transverse to the mean direction of motion (bands in $d = 2$, sheets in $d = 3$, see Fig. 1d). These objects, in spite of internal fluctuations, have rather well-defined, exponentially-localized profiles. They are solitary structures, and do not form regular wave-trains, and are separated by low-density disordered

<footnote n="1">¹ Note that the original updating scheme proposed by Vicsek et al. in [2] defined the speed as a backward difference, while we are using a forward difference, as most studies of Vicsek-style models.</footnote>
regions inside which particles essentially perform random walks.

It is only rather far away from the transition line that these objects disappear, leaving a spatially homogeneous moving phase. But if order is then high everywhere, the density field fluctuates anomalously strongly: Consider a local sub-system containing on average $n$ particles. Normal fluctuations of $n(t)$, the instantaneous number of particles in the box, would be characterized by the fact that $\Delta n \equiv \langle (n(t) - n)^2 \rangle^{1/2}$, the root mean square of the fluctuations scales like $\sqrt{n}$. Instead, here, we observed that $\Delta n$ scales like $n^\alpha$ with $\alpha > \frac{1}{2}$ (so-called “giant number fluctuations” (GNF), see Fig. 1c). Careful numerics in this region also reveal that, relative to the mean motion, individual particles undergo superdiffusive transverse behavior.

## 2 From polar to apolar particles or interactions

### 2.1 The nematic equivalent of the VM

The particles in the VM can be considered as “polar” (they carry a velocity vector). Motivated by the prediction of GNF in active nematics by Ramaswamy et al. [7], we introduced, a few years ago, the (uniaxial) nematic equivalent of the VM [5] (in which each particle can be seen as carrying a rod).

The interaction rule (1) of the polar VM involves the local order parameter in the neighborhood surrounding each particle. A similar rule can be defined for apolar particles, where $\sum_{j \in S_i} \mathbf{v}_j(t)$ is replaced by the eigenvector of the largest eigenvalue of the nematic tensor calculated in the neighborhood. This eigenvalue is directly related to the local order parameter, which, for uniaxial nematics in two space dimensions, is $|\sum_{j \in S_i} \exp(2i\theta_j(t))|$ ($\theta_j$ is the angle defining the direction of $\mathbf{v}_j$, see [5] for details). To comply fully with the nematic symmetry, the streaming step (2) of the polar VM is replaced by

$$ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) \pm \Delta t \mathbf{v}_i(t + \Delta t). \quad (3) $$

(The $\pm$ sign means that either direction is chosen randomly with probability $\frac{1}{2}$.)

The change of local symmetry for both the particles and their interaction leads to drastic changes at the collective level. In particular, in two space dimensions, only quasi-long-range order arises and the transition is of the Kosterlitz-Thouless type, now in agreement with the XY model. In three dimensions, true long-range order is observed [6].

As predicted by Ramaswamy et al. [7], GNF are easily observed in the ordered phase in both two and three space dimensions (Fig. 2a). Spatial configurations in the ordered phase consists typically of a single, macroscopic high-density high-order elongated structure (Fig. 2b) with ill-defined, highly fluctuating interfaces. These fluctuations are strong enough to trigger the rare events of a complete rearrangement of the dense cluster, on typical time scales which diverge with system size.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Vicsek-like model with nematic particles and interactions (see text). (a) GNF in a two dimensional system of linear size $L = 256$ with periodic boundary conditions ($\rho = 0.5$, $\eta = 0.08$, $v_0 = 0.5$). The dashed line has slope 1, the prediction of Ramaswamy et al. [7]. (b) configuration in the ordered phase in three dimensions (cubic domain of linear $L = 50$, periodic boundary conditions, $\rho = 1/8$, $\eta = 0.05$, $v_0 = 0.5$, color code according to local density from blue (low) to yellow (high)).</figcaption></figure>

### 2.2 An interesting mixed case

Vicsek-style models can also be constructed in which the symmetry of the particles differ from that of their interactions. A case of particular interest is that of polar particles interacting nematically. Indeed, when thinking, say, of shaken elongated polar granular particles, it is likely that their inelastic collisions are rather nematic than purely polar: when colliding almost head-on, they will glide along each other rather than align. The same is true for self-propelled rods [8,9], or for molecular filaments displaced by motors [10].

The Vicsek-style version of this situation is as follows: particles carry a vector, so that their streaming is done according to (2), but their new velocity $\mathbf{v}_i(t + \Delta t)$ is given by the local nematic tensor in the neighborhood: the polar nature of the particles is disregarded and only their axis direction is retained to calculate the interaction. Among
the two opposite directions corresponding to the resulting nematic axis thus calculated, the one closest to $\mathbf{v}_i(t)$, the previous velocity of the particle, is chosen.

We are currently investigating this case in detail. Naturally, both the polar and the nematic order parameter have to be considered, and one main question is whether nematic order can set in before polar order (the reverse is trivially impossible). Some preliminary results at relatively low density ($\rho = \frac{1}{4}$) are shown in Figure 3. They indicate that the isotropic-nematic transition, at least for this density, is a discontinuous (first-order like) one. While the polar order parameter stays near zero (no polar order), the nematic order parameter, in the transition region, switches randomly between order and disorder (Fig. 3a), leading to a characteristic bimodal distribution (Fig. 3b). Thus this case is indeed intermediate between the pure polar (VM) and the pure nematic models. The ordered phase consists, at least in the transition region, of a typically single high-density band, with the particles circulating along its length in both directions: it has high nematic order, but essentially no polar order (Fig. 3c).

### 3 Adding cohesion

The Vicsek model, as such, cannot maintain the cohesion of a moving group: if the particles evolve in an infinite domain, they will eventually fly apart. In other words, it is only at finite density that collective motion can arise. Moreover, particles in the VM are pointlike, and thus have no physical size. These two shortcomings are serious when modeling cohesive motion of a fish school, a bird flock, or, worse, cells in a tissue.

There are various ways of insuring cohesion without resorting to global interactions, but they all more or less amount to a classic pairwise attraction-repulsion mechanism. (This is in fact what most models of animal motion do, see in particular the seminal work of Huth and Wiszel [11].) The (polar) VM has been complemented along these lines [12]. The alignment interaction term in (1) then competes with a pairwise interaction and the noise strength, leading to:

$$ \mathbf{v}_i(t + \Delta t) = v_0 \vartheta \left[ \alpha \sum_{j \in S_i} \mathbf{v}_j(t) + \beta \sum_{j \in S_i} f_{ij} \mathbf{e}_{ij} + \eta n_i \mathbf{z} \right], \quad (4) $$

where $\mathbf{z}$ is a random unit vector, $n_i$ the number of particles interacting with particle $i$ including itself, $\alpha$ and $\beta$ measure the relative strengths of alignment and attraction/repulsion with respect to the noise, $\mathbf{e}_{ij} = \vartheta[\mathbf{r}_j^t - \mathbf{r}_i^t]$, and

$$ f_{ij} = \begin{cases} -\infty & \text{if } r_{ij} < r_c, \\ \frac{1}{4} \frac{r_{ij} - r_c}{r_a - r_c} & \text{if } r_c < r_{ij} < r_a, \\ 1 & \text{if } r_a < r_{ij} < r_0. \end{cases} \quad (5) $$

with $r_{ij}$ the distance between $i$ and $j$.

In addition, in [12], the neighbors have been further restricted to those being Voronoi neighbors: the distance criterion $r_{ij} < r_0$ is sued first, but among these neighbors only those belonging to the first Voronoi shell are retained. This confers a topological nature to the interactions, much alike what is advocated in [13], but with the difference that if a topological neighbor (strictly Voronoi) is much further away than the others, it will not interact with the central particle.

An extensive numerical study of this model in two space dimensions has allowed to sketch the asymptotic phase diagram in the $(\alpha, \beta)$ plane (at fixed $\eta$ and in the zero-density limit). At moderate $\alpha$ and $\beta$ values, a cohesive, orientationally-ordered phase exists. This “moving droplet” regime is the one of interest for real animal groups. It was also found that at the onset of collective motion (which is first-order, as in the VM), cohesion cannot be maintained.

Another finding of importance is the transverse superdiffusive nature of individual trajectories in the moving phase, relative to the center of mass motion. In $d=2$, the mean square displacement scales like $t^{\frac{4}{3}}$, in agreement with some calculations or Toner and Tu [4], whereas the exponent is close to 1.7 in three dimensions (Fig. 4a).

The flock shapes obtained in three dimensions are rather realistic when compared to typical starling flocks (Fig. 4b). Such superficial comparison can be misleading though: recently, the Roman team of the StarFlag project succeeded in reconstructing the three-dimensional coordinates of real starling flocks comprising up to 2000 birds [13]. This revealed that, often, the flocks are in fact rather flat. This flatness cannot emerge from the simple cohesive model described above as it is built on strictly polar particles and interactions. We are currently working on extensions on this model in which the complete symmetry of a “bird” is taken into account, namely the existence of a wingplane and a possible modulation of interactions depending on the positions of neighbors relative to it.

### 4 Role of the ambient fluid

In many models for collective motion (including the VM), the fluid in which the animals or cells move, is neglected. While this is not a problem for, say, herds of gnus, this is debated for fish schools and bird flocks, for which the high Reynolds number turbulent flow generated amounts to a modified effective medium [14]. More crucially, at the very low Reynolds numbers corresponding to bacteria swimming in a fluid, the long-range hydrodynamic interactions are probably dominant. (And for man-made micro swimmers or self-propelled nano rods they are the only interactions.) Recent experiments indicate that hydrodynamic effects influence the collective behavior of swimming bacteria at large enough particles concentrations [15]. In our view, the results of [15] show that they prevent the emergence of long-range order, as suggested earlier by Ramaswamy et al. [7]. In fact, freely swimming bacteria may be interacting only via hydrodynamic effects,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Vicsek-like polar particles with nematic interactions (see text). Two dimensional system of linear size L = 256 with periodic boundary conditions ($\rho = \frac{1}{4}$, $\eta = 0.2395$, $v_0 = 0.5$). (a) Time series of the nematic order parameter in the transition region. (b) Distribution function of the nematic order parameter from the time series of (a). (c) Snapshot in the ordered phase. Most particles move along the high-density band, in both directions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Vicsek-like model with cohesive interactions (see text). (a) root mean square displacement of individuals vs time in the “moving droplet” cohesive phase (no positional order) in three dimensions (groups of 4096 particles, $\alpha = 31.5$, $\beta = 1.5$, $\eta = 1$). The dashed line has slope 1.7. (b) snapshot of a three-dimensional cohesive flock of 16 384 particles in a similar regime.</figcaption></figure>

locally aligning in the quadrupolar field generated at short scales by their elongated shape, while long-range dipolar interactions may prevent large-scale collective motion.

The simplest way to introduce such interactions in a Vicsek-like model is to introduce a fluid velocity field $U(\mathbf{r})$ advecting the swimming particles which, in turn, generate the velocity field. In two spatial dimensions, a proxy for the experimentally relevant case of bacteria swimming in a fluid film, confinement acts as a sink for momentum perturbations generated by the swimmers, while only mass displacement plays a role in the far field [16]. Moreover, many-body screening effects cancel out, and long-range interactions can be computed as a linear superposition of elementary dipoles: $U(\mathbf{r}) = \sum_i \mathbf{u}_i(\mathbf{r})$ with

$$[\mathbf{u}_i(\mathbf{r})]_{\parallel} = v_0 \frac{x^2 - y^2}{(x^2 + y^2)^2}; \quad [\mathbf{u}_i(\mathbf{r})]_{\perp} = v_0 \frac{2xy}{(x^2 + y^2)^2} \quad (6)$$

where $[\mathbf{u}_i]_{\parallel}$ is the component along $\mathbf{v}_i$, $[\mathbf{u}_i]_{\perp}$ the transversal one, $x = [\mathbf{r} - \mathbf{r}_i]_{\parallel}$, and $y = [\mathbf{r} - \mathbf{r}_i]_{\perp}$. Note that $U(\mathbf{r})$ is irrotational, and thus only advective: it does not rotate the swimmer’s orientation. In this approximation, hydrodynamic interactions are indeed long-range and decay as $1/r^2$. In practice, the dipoles are screened-off below the (short) interaction range of the VM: When $|\mathbf{r} - \mathbf{r}_i| < r_0$, the denominators in (6) are replaced by $r_0^4$, and the VM alignment dominates (it can be thought of mimicking the quadrupolar component mentioned above in this context).

With this advecting field thus constructed, the streaming step of the VM is replaced by:

$$\mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \Delta t [\mathbf{v}_i(t + \Delta t) + \gamma \mathbf{U}(\mathbf{r}_i(t))]. \quad (7)$$

where $\gamma$ controls the strength of the long-range interactions.

Preliminary results indicate that long range hydrodynamic interactions can indeed destroy coherent motion (Fig. 5). The band structure can resist some amount of long-range interactions (Fig. 5b), but polar order is destroyed at higher $\gamma$ values (Fig. 5c). Surprisingly, at still larger values of $\gamma$ however, order is restored, but typical configurations reveal high-density high-order blobs rather than bands (Fig. 5d).
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Vicsek-like polar particles with hydrodynamic long range interactions (see text). Two dimensional system of linear size $L = 256$ with periodic boundary conditions ($\rho = 1/8$, $\eta = 0.3$). (a) Time-averaged order parameter as a function of hydrodynamic interaction strength. (b–d) Snapshots of particle positions and orientations for three different values of $\gamma$ (as indicated in (a)). The red arrows point to the mean direction of motion.</figcaption></figure>

Beyond, further work is clearly needed to clarify the behavior of this model. A “real” fluid need to be introduced in order to tackle the full problem. Several groups are currently studying the collective dynamics of “Stokeslets” [17], but this approach will remain numerically limited in the number of swimmers that can be considered.

## 5 Conclusion

We hope the results presented here may help put the Vicsek model in perspective. Clearly, much remains to be done to explore the avenues only sketched here. The variations on the Vicsek model presented here remain simple enough that one stays away from the dangers of overparameterization, and very large groups can still be studied (something desirable in view of the strong finite-size effects of the original VM). Yet, they allow to approach many interesting situations of current interest in physics and beyond.

Particularly intriguing is the emergence of “giant density fluctuations” in many of the cases considered, provided the spontaneous segregation between high-density/high-order regions and low-density/disordered regions takes place. This is similar to the clustering instability in granular particles colliding inelastically: indeed, the alignment in VM-style models can be seen as an effective description of two-body inelastic collisions. We stress that this non-equilibrium effect, in particular in the case shown in Figure 2a, is striking: the rms of the fluctuations scales then like the mean density: a local measurement of density is thus impossible!

At the theoretical level, reliable, faithful mesoscopic descriptions of the situations described here are needed before one attempts to apply field-theoretical techniques to them. Some progress has recently been recorded towards this [18], but we are still missing stochastic descriptions including the relevant noise terms (which we believe to be important, especially in view of the strong density fluctuations).

Finally, model experiments dealing with large numbers of self-propelled or swimming objects are needed. In this respect the recent findings on shaken rod-like objects [1] are encouraging.

Part of this work was funded by the European StarFlag and the French ANR Morphoscale projects.

## References

1. V. Narayan et al. Science **317**, 105 (2007)

2. T. Vicsek et al. Phys. Rev. Lett. **75**, 1226 (1995); A. Czirók, H.E. Stanley, T. Vicsek, J. Phys. A **30**, 1375 (1997)

3. H. Chaté, F. Ginelli, G. Grégoire, F. Raynaud, Phys. Rev. E **77**, 046113 (2008); G. Grégoire, H. Chaté, Phys. Rev. Lett. **92**, 025702 (2004)

4. J. Toner, Y. Tu, Phys. Rev. Lett. **75**, 4326 (1995); Phys. Rev. E **58**, 4828 (1998); J. Toner, Y. Tu, M. Ulm, Phys. Rev. Lett. **80**, 4819 (1998)

5. H. Chaté, F. Ginelli, R. Montagne, Phys. Rev. Lett. **96**, 180602 (2006)

6. H. Chaté et al. unpublished

7. R.A. Simha, S. Ramaswamy, Physica A **306**, 262 (2002); Phys. Rev. Lett. **89**, 058101 (2002); S. Ramaswamy, R.A. Simha, J. Toner, Europhys. Lett. **62**, 196 (2003)

8. F. Peruani, A. Deutsch, M. Bär, Phys. Rev. E **74**, 0309904(R) (2006)

9. A. Kudrolli, G. Lumay, D. Volfson, L.S. Tsimring, Phys. Rev. Lett. **100**, 058001 (2008)

10. P. Kraikivski, R. Lipowsky, J. Kierfeld, Biophysical J. **483** (2007); D. Karpeev et al. Phys. Rev. E **76**, 051905 (2007)

11. A. Huth, C. Wissel, J. Theor. Biol. **156**, 365 (1992)

12. G. Grégoire, H. Chaté, Y. Tu, Physica D **181**, 157 (2003); see also Phys. Rev. Lett. **86**, 556 (2001); Phys. Rev. E **64**, 011902 (2001)

13. M. Ballerini et al. Proc. Natl. Acad. Sci. USA **105**, 1232 (2008)

14. See for instance B.L. Partridge, T.J. Pitcher, Nature (London) **279**, 418 (1979); P.B.S. Lissaman, C.A. Shollenberger, Science **168**, 1003 (1970); D. Weihs, Nature **241**, 290 (1973); J.P. Badgerow, F.R. Hainsworth, J. Theor. Biol. **93**, 41 (1981)

15. A. Sokolov, I.S. Aranson, J.O. Kassler, R.E. Goldstein, Phys. Rev. Lett. **98**, 158102 (2007)

16. H. Diamant, B. Cui, B. Lin, S.A. Rice, J. Phys.: Condens. Matter **17**, S2787 (2005)

17. J.P. Hernández-Ortiz, J.J. de Pablo, M.D. Graham, Phys. Rev. Lett. **98**, 140602 (2007); D. Saintillan, M.J. Shelley Phys. Rev. Lett. **99**, 058102 (2007)

18. E. Bertin, M. Droz, G. Grégoire, Phys. Rev. E **74**, 022101 (2006); I.S. Aranson, A. Sokolov, J.O. Kessler, R.E. Goldstein, Phys. Rev. E **75**, 040901(R) (2007)