Collective Behavior of Animals: Swarming and Complex Patterns

In this short note we review some of the individual based models of the collective motion of agents, called swarming. These models based on ODEs exhibit a complex rich asymptotic behavior in terms of patterns, that we show numerically. Moreover, we comment on how these particle models are connected to partial differential equations to describe the evolution of densities of individuals in a continuum manner. The mathematical questions behind the stability issues of these PDE models are questions of actual interest in mathematical biology research.


Introduction
Self-organization without the presence of a leader is a common feature in the collective behavior of certain animals: insects (locusts, ants, ...), fish, birds (in particular, starlings [8,35,2]) and some microorganisms such as myxobacteria [28].The variety and astonishing behavior of the aggregate patterns formed by groups of individuals is inspiring as well as attractive to scientists for deeper explanations.The biological reasons [35] of the grouping include, among others: protection against a predator, spawning migrations, energy savings due to flying or swimming in an ensemble, and food finding.It is well-known that effects such as pheromone trails for ants may produce very sophisticated behavior which can be very efficient, as in this case, at finding food, and is the result of the individual behavior of each ant, with no organization decided by any particular ant [6].
There have been some attempts to give mathematical models of collective motion, based on reasonable assumptions about the behavior of each individual, and which can reproduce the observed phenomena.Many of the models naturally involve some kind of evolution equation for the position, velocity, and possibly other characteristics, of individuals, be it in the form of a set of ordinary differential equations or a set of difference equations for the values of these quantities at certain discrete times [43,20,14,46].Such models are usually referred to as Individual-Based Models (IBMs).They follow many different strategies; for instance, one can try to give rules for the movement of each particular individual which are as realistic as possible.These rules then can become mathematically very complicated, and thus often the only way to understand the evolution of the model is to solve it numerically with a computer.One can then compare the result to the known behavior of real individuals, as has been done for example in [26] for starlings, using real data provided in [2], or in [24,25,4] for fishes, and possibly obtain a better understanding of the causes of the actual behavior of the animals under study.
One can also extract the basic features of the above models and try to understand the emergence of collective behavior from simpler models which lend themselves to a deeper mathematical study.While this is admittedly further away from the experimental explanation of the behavior of animals, it has the advantage that in some cases one is able to understand much better the behavior of the model, and identify the role that each effect plays.There are relatively simple models which pose very interesting mathematical questions, sometimes leading to problems in dynamical systems or kinetic equations which can benefit from existing techniques, which are interesting by themselves, and which may find applications in fields other than animal behavior (see below).Hence, the study of animal behavior suggests intriguing directions of research in mathematics; and conversely, a basic understanding of simple effects can give a valuable insight in the design of more complicated models that one cannot expect to be able to handle theoretically at present.
Here, we will be concerned with models which give the evolution of a set of individuals through a certain number of effects.Most of them include alignment or orientation effects, possibly involving some randomness.The self-organization of agents described by IBMs when self-propelling forces and pairwise attractive and repulsive interactions are considered was described in [34,31,18].The classification, in terms of particular strengths of interaction and propulsion, of the morphology of patterns obtained includes: translationally invariant flocks, rotating single and double mills, rings and clumps, patterns that we will show later.
Flocking patterns have also been shown in models of alignment or orientation averaging in [16,17].All these models share the objective of pinpointing the minimal effects or interactions leading to certain particular type of pattern or collective motion of the agents.These IBMs can be considered also as "particle" models in the optic of treating agents or animals as point particles in a physics-based description as in statistical mechanics.Let us also mention that this issue has also received attention from the control engineering viewpoint trying to reproduce these self-organization patterns with artificial robots or devices, see [12,16,36] and the references therein, with the aim of controlling unmanned vehicle operation.Finally, the combination of these minimal bricks in the modelling such as interaction, alignment and orientation with other effects is leading to rich complex behavior and detailed dynamical systems models for particular species, see for instance [5,4,25,2,26,6,46].
In some of the applications above, IBMs models are enough to describe the system under reasonable number of individuals/agents N .However, if the number is large as in migration of fish [47] or in myxobacteria [35,28], the use of continuum models for the evolution of a density of individuals is convenient for numerical simulation, and even necessary.Some continuum models in the literature [41,42,7] include attraction-repulsion mechanisms and spatial diffusion to deal with random effects.Other continuum models are based on hydrodynamic descriptions [13,10] derived from mean-field particle limits.In fact, as usually done in statistical physics, there is a middle ground in modelling between particle and hydrodynamic descriptions given by the mesoscopic kinetic equations describing the probability of finding particles in phase space.Kinetic models of swarming has recently been proposed [22,10,11] and the connection between the above IBMs and the continuum models via kinetic theory has been tackled very recently in [21,10,9].The interest of the kinetic theory models is to give a rigorous tool to connect IBMs and hydrodynamic descriptions as well as to interpret certain patterns as solutions of a given model, as in the case of the double mills [10].The analysis, numerical description of the complex behavior of these kinetic and hydrodynamic models and patterns stability are some of the open questions in this research direction.
Below we will focus on a review of certain IBMs for swarming proposed in the literature, including some of the basic effects above: attraction, repulsion and orientation.We will describe these particle models and explain some of their features and basic properties in terms of asymptotic patterns in Section 2. Subsection 2.3 is devoted to a nice short proof of flocking for the Cucker-Smale model and its variants, and finally Section 3 shows some of the numerical results and open questions regarding these asymptotic patterns for the continuum models.

Individual-Based Models
Most of the basic IBMs of collective behavior or swarming consider the so-called three-zone models, meaning that each individual is influenced in a different way by other individuals depending on their relative position, and distinguishing three different zones; see Figure 1.The inner region, repulsion zone, models the "vital space" in which any individual has a tendency to avoid the presence of any other in the swarm, due to collision avoidance or just for sociological reasons.The intermediate region, orientation or alignment zone, takes into account the process of mimicking other individual's behavior by trying to match the velocity of others, for instance with the others in that region.The outer region models the effect of inherent socialization of the animals since they want to be not too far from other individuals in that region.The strength, particularities and details of these different effects depends on the distance, and on the number of other individuals located at each moment in the different zones; hence the direction of the individual will be changed according to some weighted superposition of these effects.The "attraction", "repulsion" and "velocity mimicking" are loosely stated here and can take a wide range of forms.This modelling based on social forces finds its roots in the works of [1,27] for fish schooling and it has been widely used, expanded, and improved by theoretical biologists, physicists and applied mathematicians; see [15,34,29,24,30,44,45,3,31,32,19] and the references therein.This approach has been made much more specific for particular animals and species as in [5,4,47,25] for fish (studying migration patterns for the capelin around Iceland) and in [2,26] for birds (starlings grouping in Rome).
In this review, we concentrate mainly in these three effects as they constitute the main modelling "bricks", but the reader should be warned, as we mentioned before, that in order to reproduce realistic swarming behavior as in the starlings [2] or the migration of fish [4], one has to include much more involved interactions.For instance, real three dimensional effects due to the aerodynamics of birds, in which drag, lift and friction forces are included, or roosting forces trying to model their tendency to stay close to certain home area as in [25]; or changes in movement due to currents or temperature changes as in [4].Other improvements in the modelling include detailed vision zones for the individuals, assuming for instance that birds only see in a certain vision cone which depends on their position and direction of motion.
In this paper we address particle models for describing mathematically these different effects.Despite their simplicity, these models show how the combination of these simple rules can produce striking phenomena, such as pattern formations: flocks and single or double mills, resembling those observed in nature.

Asymptotic speed with attraction-repulsion interaction model
The particle model proposed in [18] reads as: (2.1) where x i , v i ∈ R 3 represent the position and velocity of the i-th individual, for i varying from 1 to the total number N of individuals.Here, α and β are nonnegative parameters and U : R d −→ R is a given potential encoding the short-range repulsion and long-range attraction typical in these models.Here, the potential has been scaled depending on the number of particles as in [10], as it is convenient in order to study the limit of the system for large N .The term corresponding to α models the self-propulsion of individuals, whereas the term corresponding to β is the friction assumed to follow Rayleigh's law.The balance of these two terms imposes an asymptotic speed to the agent (if other effects are ignored), but does not influence the orientation vector.The main point here is that of having a "preferred" velocity α/β, and not the precise shape of the term involving α and β.On the other hand, a typical choice for U is the Morse potential, a radial potential given by where C A , C R and A , R are the strengths and the typical lengths of attraction and repulsion, respectively.This potential is only Lipschitz due to its singularity at the origin but the qualitative behavior of the particle system does not heavily depend on this fact [18].To analyze the limit of large number of particles N , and for simplicity, we scale the amplitude of the potential through a normalization which corresponds to assuming that all particles have mass 1/N .The interesting particularity of this system is that it exhibits very different behavior depending on the values of the parameters.Let us have a look at some possibilities.Taking the values C R = 50, R = 2, C A = 20, A = 100, β = .05,α = .07and solving numerically with some randomly chosen initial conditions, we find that individuals arrange into a sort of swirl or mill, as this type of arrangement has been called in the literature (see Fig. 3).After some hesitation, they seem to always agree to turn in only one direction.There are yet other possible patterns in which individuals tend to organize, depending on the value of the parameters.For instance, sometimes they organize in a kind of crystalline structure and they move translationally; this pattern is called flock.Identifying the tendency for each parameter seems difficult, and there are even some parameters for which one observes sometimes one type of organization, sometimes another one.
However, part of it can be analyzed, at least heuristically.In the case of the Morse potential, patterns of aggregation depend on the relative amplitudes C = C R /C A and = R / A , and the number of individuals N does not seem to affect the qualitative features of the observed patterns, so we may center the discussion on the values of C and .
In two dimensions, the different patterns were classified in [18] using the concept of Hstability of potentials [37,40].The most relevant set of parameters for biological applications concerns long-range attraction and short-range repulsion leading to C > 1 and < 1.For these potentials, there exists a unique minimum of the pairwise potential and a typical distance minimizing the potential energy.What they remarked is that the behavior of the system (2.1) depends on whether the parameters are in the so-called H-unstable or catastrophic region C d < 1, or in the H-stable region C d ≥ 1, where d stands for the space dimension, equal to 2 in the numerical solutions shown in figures.This terminology for potentials comes from the field of statistical mechanics.Complex behavior such as the formation of mills or double mills seems to happen only in the catastrophic region, while individuals seem to either form a swarm or just disperse in the stable region.There is also a difference between the catastrophic and stable regions in the size of "aggregates" (swarms or mills) as N grows.However, it does not seem easy to decide, in the catastrophic region, which parameters will produce swarms, mills, double mills, rings, or other patterns.
These shapes may contain some clue on the organization of animals in nature: while coherent flocks and single mill states are the most common patterns observed in biological swarms [35,38], double-mill patterns, as seen in Figure 2, are also reported in the biological literature; for instance M. xanthus cells show distinct cell subpopulations swarming in two opposite directions during part of their life cycle [28].

The Cucker-Smale model
In the Cucker-Smale model, introduced in [16,17], the only mechanism taken into account is the reorientation interaction between agents.Each agent in the swarm tries to mimick other individuals by adjusting/averaging their relative velocity with all the others.This averaging is weighted in such a way that closer individuals have a larger influence than further ones.For a system with N individuals the Cucker-Smale model, normalized as before, reads with the communication rate w(x) given by: γ for some γ ≥ 0. For this model, one would expect that individuals tend to adopt finally the same velocity and move translationally, as they change their velocity to adapt to that of others.This behavior is observed indeed by numerically solving the above equations, see Fig. 5, with γ = 0.45.
After several improvements, it has been shown in [16,17,22,21,11] that the asymptotic behavior of the system for γ ≤ 1/2 does not depend on N ; in this case, called unconditional non-universal flocking, the behavior of the population is perfectly specified: all the individuals tend to flock."Flocking" here means that they end up moving with the same velocity, this forming a group with fixed mutual distances, not necessarily in a crystalline-like pattern, but rather depending on the the initial positions and velocities (actually, any set of individuals moving initially at the same speed will continue to do so indefinitely regardless of their initial positions).On the other hand, in the regime γ > 1/2, flocking can be expected under certain conditions on the initial configuration but there are counterexamples to the generic flocking [16].We refer to [16,17,11] for further discussion about this model and its qualitative properties, but we do present now a short proof of the flocking behavior when γ ≤ 1/2.

Alignment in the Cucker-Smale model
To see how the alignment comes up, instead of just considering the Cucker-Smale model, we may consider a more general model in which the averaging takes into account the strength of the relative speed, with p > 0. This model reduces to the original Cucker-Smale model for p = 2. First of all let us point out some facts and set some notation that will be useful in the following discussion.
Taking into account that these equations are invariant by translations, we can set the mean velocity of the system to zero without loss of generality.In this way, the center of mass, x c , will be preserved along the evolution.We fix R x 0 > 0 and R v 0 > 0 such that the particles are initially inside the ball B(x c , R x 0 ) and the initial velocities are in B(0, R v 0 ).Also, we define the function At any time t, we can choose an index i such that R v (t) = |v i (t)|.Note that since we are dealing with a finite number of particles and the curves are smooth, R v (t) is differentiable except possibly on a set of measure 0, and our argument below can be carried out rigorously.We will show the proof for values of γ < 1/2; improvements of this can be made following the steps in [11].
We prove that the concentration of the velocities around its mean value through a series of recurrent steps, in each of which the bounds we have for R x (t) and R v (t) are alternatively improved until we find an estimate for the concentration rate of the velocities.As a first step, we compute the derivative of R v (t) 2 with respect to time, which gives us due to the choice of the label i, which ensures that (v i − v j ) • v i ≥ 0, for all j.Hence, R v (t) is a non-increasing function.As a consequence, the distance between a particle and its original position can only grow linearly in time.Thus, where R = max(R x 0 , R v 0 ).In turn, this implies that we have a lower bound for the weight: Now, using the bound on w ij , we can extract more information from the computation of the time derivative of R v (t) 2 : At this point we need to be careful with the exponent p − 2. Due to the choice of i, since we have set the mean velocity to zero the second term in the right hand side of the equation disappears.Then we get where we have set w(t) := 1 + |2R(t + 1)| 2 −γ .In case 2 < p < 4, we can estimate it from below as j =i and expanding Summarizing, we have that for any 0 < p < 4 and integrating with respect to time we obtain, for p < 2 and 2 < p < 4, where W (t) := t 0 w(s) ds.This result can be extended by using a slightly more delicate argument to any p > 2, but since the next steps will require a stronger assumption on p than p < 4, used at this point, we omit it here.The case p = 2 needs a separate treatment to show that R v (t) is exponentially decreasing.These details can be found in [11] by the interested reader, here we will only report the details for the case p = 2. Coming again to (2.5), we notice that W (t) is an increasing function of time and for γ < 1  2 it diverges as t → ∞.Thus, for p < 2 there exists T < ∞ such that R v (T ) = 0, so the concentration in velocity holds in finite time.On the other hand, for p > 2, the exponent (2 − p) −1 is negative, R v (t) → 0 algebraically as t → ∞, and the concentration in velocity holds in infinite time.
Once a concentration in velocity has been established, we can improve the bound on the positions to show that |x i (t) − x j (t)| remains uniformly bounded in time.For 0 < p < 2, if we define f (t) ≡ 0 for t ≥ T we have that On the other hand, for 2 < p < 4, f (t) behaves like W (t) 1/(2−p) as t → ∞.In order for f (t) to be integrable in time up to infinity, we assume further that 2 < p < 3 − 2γ < 4. Proceeding as in (2.3) and (2.4), for 0 < p < 2 and 2 < p < 3 − 2γ, we have a uniform in time bound on |x i (t) − x j (t)| and an uniform bound from below in w ij ≥ W 0 > 0. Finally, we use again the computation of the derivative of R v (t) 2 to finally get for all 0 < p < 2 and 2 < p < 3 − 2γ.
Summarizing, we have shown that for this modified version of the Cucker-Smale model, the flocking behavior happens in a finite time for 0 < p < 2 and 0 < γ < 1/2 and with an algebraic speed if 2 < p < 3 − 2γ in contrast with the exponential speed obtained for the standard Cucker-Smale model with p = 2 and 0 < γ < 1/2.This can be checked numerically as seen in Fig. 6 for different values of p.It is an open problem to check if these limits for the parameter p and γ when p > 2 are sharp for the flocking pattern to appear for generic initial data.

3-zone model
Let us consider now the model with attraction, repulsion, and Cucker-Smale effects included, which can then be properly called a "three-zone model" with three zones which are not at all disjoint: (2.7) Figure 6: Different convergence rates to the zero mean velocity, in log-scale, for nonlinear velocity-dependent Cucker-Smale models.
How does the asymptotic behavior of this model look like?Individuals tend to eventually adopt the same velocity due to the Cucker-Smale effect, as in section 2.2.They also try to arrange into some group, due to the potential interaction.The result is that they tend to an arrangement which should be a local minimum of the potential energy, but this does not seem obvious at all to prove.For instance, if we take the values C R = 500, R = 2, C A = 200, A = 100, β = .1,α = .2,γ = 0.45 and look at the system with 50 individuals with random initial conditions, we observe numerically that no matter how we start, they seem to end up arranging themselves in a group-like fashion in a crystalline structure, and moving in some direction at the preferred speed α/β, see Fig. 7.

Some open questions
Very few of the observed features of the above behavior can be actually proved in a rigorous way.Apart from the asymptotic convergence to the mean velocity of the Cucker-Smale model, we do not have a way to distinguish the values of the parameters in the system (2.1) for which one behavior or other takes place, and we do not know how to prove the convergence to a swarm, a mill or other kind of organization.A feasible way of studying these problems may be the following: one can consider, instead of the system (2.1) or (2.7), a partial differential equation which is obtained as its mean-field limit; this is an evolution equation whose solutions are approximated by the evolution of the density of particles of the system (2.1) (or (2.7)) when the number of particles is very large [9].Then, studying stationary states for this equation may be a simpler task than directly studying (2.7).For example, the partial differential equation obtained as a limit when N → ∞ of the system (2.7) is where f = f (t, x, v) represents the density of individuals at time t in an infinitesimal region dx dv, ρ = ρ(t, x) is the macroscopic density of individuals at time t obtained from f by integration on v, the first convolution is in the x variable, and the second one in both the x, v variables [9].The function H is given by "Flocking" solutions of the system (2.7) (solutions for which every individual moves at the same fixed velocity v 0 ) correspond to solutions of (3.8) of the form where δ is the Dirac delta function.These are solutions which have a constant mass profile ρ, and in which every point moves at the same velocity v 0 .Can we find solutions like this?Looking for such a profile ρ turns out not to be a simple problem, and very little is known about it.We may restrict ourselves to looking for solutions with velocity v 0 = 0, as all others are just translations in velocity of these.Then, if we look at the limiting shape of flocking solutions of the system (2.7), we see that the "flocks" tend to grow indefinitely with N when the potential is in the H-stable region; hence, in this case, the conjecture is that there are no nice (say, continuous and compactly supported) functions ρ for which (3.9) is a solution of (3.8).On the other hand, when the potential is not H-stable, or "catastrophic", the flocks seem numerically to converge to an asymptotic distribution as we add more and more particles.As shown in Table 1 and in the first plot of Fig. 8, the maximum radius of the particles with respect to the center of mass tends to stabilize as t → ∞ and as N gets larger to a fixed value.Moreover, we have computed the normalized cumulative distribution of particles in the radial direction starting from the center of mass as we increase the number of particles N .The second plot in Fig. 8 shows numerical evidence of the convergence towards a fixed continuous profile as N gets larger.The conjecture in this case is that there are continuous and compactly supported profiles ρ for which (3.9) is a solution of the kinetic equation, and that there are N -individual flocks whose density converges to ρ in the limit of N going to ∞.These problems are, to our knowledge, open for the moment, and they are also not particular to this field: they are more generally related to the shape of N -particle equilibria for a given potential U .

Figure 2 :
Figure 2: Mills in nature and in the IBMs!

Figure 3 :
Figure 3: Formation of a single mill.

Figure 4 :
Figure 4: Formation of a double mill.

Figure 5 :
Figure 5: Formation of a non-universal flock.

Figure 7 :
Figure 7: Formation of a flock.

Figure 8 :
Figure 8: Evolution with respect to time of the maximum radius and cumulative distribution of particles with respect to radius computed for 3 different number of particles.

Table 1 :
Maximum distance of the particles to the center of mass with respect to the number of particles in the simulation.