Modelling Hierarchical Flocking

We present a general framework for modeling a wide selection of flocking scenarios under free boundary conditions. Several variants have been considered - including examples for the widely observed behavior of hierarchically interacting units. The models we have simulated correspond to classes of various realistic situations. Our primary goal was to investigate the stability of a flock in the presence of noise. Some of our findings are counter-intuitive in the first approximation, e.g., if the hierarchy is based purely on dominance (an uneven contribution of the neighbors to the decision about the direction of flight of a given individual) the flock is more prone to loose coherence due to perturbations even when a comparison with the standard egalitarian flock is made. Thus, we concentrated on building models based on leaderfollower relationships. And, indeed, our findings support the concept that hierarchical organization can be very efficient in important practical cases, especially if the leaderfollower interactions (corresponding to an underlying directed network of interactions) have several levels. Efficiency here is associated with remaining stable (coherent and cohesive) even in cases when collective motion is destroyed by random perturbations. The framework we present allows a the study of several further complex interactions among the members of flocking agents.


Introduction
Collective behaviour is a very important aspect through which small or large groups of organisms optimize their living [1]. It involves collective decision making an various contexts, such as such as searching for food [4], navigating towards a distant target [4,2,3,5] or deciding when and where to go [4,6,7]. Flocking is perhaps the most common and spectacular manifestation of collective behaviour not only in nature since recently has gained attention in the context of collective robotics as well [8,9,10]. Most of the experimental and modeling approaches aimed at describing flocking by assuming egalitarian interactions among the members of a flock.
However, just as flocking is a widespread behavioural pattern of a collective, the hierarchical structure of the interactions among the members of groups is also very much common [11]. Thus, starting with a trend-setting paper of Couzin et al [12] the question of leadership during flocking has attracted increasing interest. Early works assumed a two-level hierarchy while recent experimental observations involving some sophisticated animal groups such as pigeons or primates point towards the possibility of significantly more complex internal organization principles during group motion [2,7,13,14]. In socially highly organized groups beyond a given size (dozens or so) the roles related to leadership do not seem to be simply binary, but several levels of hierarchy can be identified. In particular the available few experimental results indicate that many comoving groups have an internal system of interactions which can be best interpreted in terms of pairwise hierarchical levels of interactions. Perhaps the best quantitatively evaluated hierarchical group motion was carried out for pigeon flocks [2,15] in which a hierarchically distributed set of interactions was demonstrated using GPS tracks and a velocity correlation-delay method [2].
On the other hand, in spite of its relevance, there have been no realistic models proposed and studied devoted to the understanding of the conditions under which hierarchical flocking can be optimal. The only vaguely related original model is in a beautiful work by Shen [16], but the assumptions are quite arbitrary and the effect of the hierarchy on the performance of the flock is not investigated. Therefore, our goal here is to design a general framework which allows the treatment of questions related to hierarchical leader-follower interactions in a general setting, when the flock is freely moving (in the present study in two dimensions). We also discuss the various variants allowed by our approach and characterize the efficiency of the behavioral rules by determining to what level a flock is stable against external perturbations. Hierarchy is thought to be prevalent because it can be shown to result in more efficient group performance [7,11].
The present work is the first one in which a flexible and plausible framework is introduced to model hierarchical flocking. It turns out that introducing pairwise interactions satisfying a realistic hierarchical dynamics is far from being trivial. After introducing the model the main aspect we investigate is the stability of a flock moving freely in a border-less two-dimensional space. We shall associate with stability the tendency of a flock not to break into parts under external perturbations.

Definitions used in our model
For an egalitarian system, every particle is believed to be the same, while for a hierarchical system, individual difference exists. Any complex multi-agent system can be classified as egalitarian or hierarchical according to the set of pairwise relationships between any two individuals. For each pair of particles in an egalitarian system, both of them have the same ability/contribution to the influence on the decision for the next time step. Meanwhile, for each pair of particles in a hierarchical system, when making a decision, they may have different level of contribution to the decision (weight) or have a directed information flow relationship, like the leader-follower mechanism (follower has no influence on the behavior of the leader).
Graphs/networks represent a useful tool for visualizing these pairwise interaction relationships of individuals belonging to one co-moving group. Therefore, we define several matrices to describe the internal properties of the hierarchical mechanisms from different points of view, in order to characterize the differences between the egalitarian and hierarchical systems more clearly.

Contribution matrix
is defined to describe the contribution strength (weight) of each particle during the decision making process regarding the new preferred directions of the particles.
(i) For an egalitarian system, every particle has the same contribution value, that is, c ij = q(q > 0), for i, j = 1, · · · , N. q is a constant.
(ii) For a hierarchical system, not all of these c ij (i, j = 1, · · · , N) have the same positive value. For example, c ij = q i (q i > 0), for i, j = 1, · · · , N, or c ij satisfies some probability distribution(such as log-normal), for i, j = 1, · · · , N. We name this kind of hierarchy as contribution driven hierarchical system.

Dominance Matrix
The dominance matrix B N = [b ij ] N ×N is defined to describe the direction of information flow between each pair individuals.
(i) For each pair of particles i and j in an egalitarian system, their behaviors can be influenced by each other, that is to say, the information between each pair of particles is transmitted bidirectionally (corresponding to an undirected graph). Thus, in an egalitarian system for each pairwise interaction b ij = b ji = 1, ∀i, j = 1, · · · , N.
(ii) If the information flow of paired particles is directional, that is, only one particle can obtain the other particle's information (directed graph), then we have a dominance driven mechanism, which is another kind of hierarchy. For paired particle i and j, if i is led by j, then we have b ij = 1 and b ji = 0, i, j = 1, · · · , N.
Leader-follower mechanism is a typical kind of dominance relation in hierarchical organizations. For each pair of particles, leader particle does not take into account the influence of the follower, but the follower considers the behavior of the leader particle when makes decision on its behavior at the next step. This feature is represented by the fact that the matrix B N is not symmetric. Without loss of generality, we number these particles according to the level of dominance from 1 to N. Particle 1 is the strongest one of the whole system, while particle N is the weakest one. Therefore, matrix B N is a complete and symmetric matrix for egalitarian systems, while matrix B N is a lowertriangular matrix for dominance hierarchical system.

Egalitarian versus hierarchical systems
Now we can give a formal definition of egalitarian and several kinds of hierarchical systems.

Egalitarian system
For each pair of individuals, if both of them will use each other's information with the same weight to decide the behavior at the next step, we say it is an egalitarian system. An egalitarian system satisfies the following two rules.
2.3.2. Contribution driven hierarchical system A contribution hierarchical system satisfies the following two rules.
(i) c ij follows some distribution for i, j = 1, · · · , N (thus, not all of these c ij have the same positive value); Single-layer leader-follower hierarchical system (dominance driven hierarchical mechanism) A single-layer leader-follower hierarchical system satisfies the following two rules.
(i) c ij = q, ∀i, j = 1, · · · , N, q is a positive constant; (ii) b ij = 1, i > j, i, j = 1, · · · , N, 2.3.4. Double-layer leader-follower hierarchical system (contribution driven dominance hierarchical mechanism) In case the weights in the contribution matrix are not equal, we associate the system with the presence of dominance (the contribution of agents having a larger weight dominate over the contribution by those with a smaller weight). We consider such systems whose behavior is determined through a contribution driven mechanism. Then, the particle with larger weight contribution is named as leader, while the other one is named as follower. A double-layer leader-follower hierarchical system satisfies the following two rules: (i) c ij meets some distribution for i, j = 1, · · · , N (but not all of these c ij have the same positive value); The above systems can be characterized by their intersection matrix c ij * b ij , i, j = 1, · · · , N. For an egalitarian system, it is a complete matrix with the same elements. For contribution hierarchical systems, it becomes a complete matrix with various elements. For single-layer leader-follower hierarchical system, it is a lowertriangular matrix with the same lower-triangular elements. For double-layer leaderfollower hierarchical system, it becomes a lower-triangular matrix with varying lowertriangular elements.
Besides, in a more general model, c ij and b ij can be time-dependent. If c ij and b ij is not time-dependent that means everyone in the system has fixed a set of relationships. If c ij (t) and b ij (t) are time-dependent, then the relationships among these particles changes with time. Other different variants can thus be defined according to the contribution matrix c ij (t) and dominance matrix b ij (t). Here, we only discuss the case that when c ij and b ij are constant.

Hierarchical model for flocking
We consider in this paper N particles moving continuously (off lattice) in a free area without any boundary limitation. As shown in Figure 1, the position and direction of N particles at the beginning are generated randomly, while over time these particles are expected to move coherently (ordered state). The figure is for the noiseless case. Suppose that the time interval between two updates of the directions and positions is ∆t = 1. This assumption can be made without loosing generality since ∆t occurs only in combination (being multiplied by it) with velocity terms. At t = 0, N particles were randomly distributed within an area of a given size and have the same absolute velocity υ as well as randomly distributed directions θ. At each time step, the position of the ith particle is updated according to In each time step, the velocity of a particle v i (t+1) is updated according to the following equation where v align i (t) is the alignment term, v rep i (t) is the repulsion term, and v adh i (t) is the attraction term.
The alignment term was constructed to have an absolute value υ and a direction given by the angle θ align i (t + 1). The angle was obtained from the expression where ∆θ(t) represents noise, which is a random number chosen with a uniform probability from the interval [−η/2, η/2]. < θ i (t) > denotes the average direction of the velocities of neighbors of the given particle i. The average direction is given by the angle The matrix L N (t) = [l ij (t)] N ×N describes the neighbor relationships of particles at time t, where The definition of adjacency matrix where N i (t) = {j| x i (t) − x j (t) < r. Here, r denotes the interaction radius. Using the above expressions the alignment term can be written as v align i (t + 1) = c align υe i (t) where c align is the coefficient of the alignment term. e i (t) is a unit vector with direction angle θ align i (t). The repulsion term exists only when the distance between any two particles is smaller than the repulsive radius r rep . And the repulsion term is defined as where x ij (t) < r rep and c rep is the coefficient of the repulsion term. The attraction term is only considered for the boundary particles [17] of the whole system when the distance between two particles is between r rep and r adh .
where r rep ≤ x ij (t) ≤ r adh and c adh is the coefficient of the attraction term. This term is introduced in order to prevent the flock spreading (or, in other words, "evaporating") due to perturbations. Figure. 2 shows the neighbor matrix of several variants of flocks mentioned in the last section. From Figure. 2, we can see more details on the difference among egalitarian system and other hierarchical system more clearly.

Simulations and discussion
The simulations were carried out in a free two-dimensional area. We considered groups of particles having various sizes ranging from 10 particles to 1280 particles. In order to keep the continuity and comparability of these simulation results of different group sizes, we chose the scale of the area for the random initial positions directly proportional to the scale of group size N, and generated the initial angle from the interval (−π, π]. In this simulation, we use υ = 0.1, r rep = 0.5, r adh = 2.2, and interaction radius r = 2.2. We aimed at comparing the stability of egalitarian systems versus various hierarchical systems for varying levels of external disturbance (noise) and for different group sizes. The alignment item of egalitarian flock model is much like the self propelled particle model proposed by Vicsek et. al. in 1995 [18], where the definition of neighbors of particle i includes itself and the contribution of the particles is the same. Thus, we call the egalitarian model as VEM (E for egalitarian, M for model). The contribution driven hierarchical model is called as CHM (C for contribution), while single-layer leaderfollower Hierarchical model is named by SHM (S for single-layer leader follower).
According to the definition of the dominance matrix B N , the neighbor matrix of a hierarchical system has zero-value diagonal elements, that is, we have l ii = 0, ∀i = 1, · · · , N for all hierarchical systems. l ii = 0 means that the neighbor set of particle i doesn't contain itself. In the following under hierarchy, we always mean hierarchical leader-follower kind of hierarchy. Therefore, we name the double leaderfollower hierarchical flock model with zero-s along the diagonal (l ii = 0, ∀i = 1, · · · , N) as 'Double Hierarchical Model with Zero-s (DHMZ)', while when l ii = 0, ∀i = 1, · · · , N, we call it 'Double Hierarchical Model(DHM)'.

Order Parameter
In this case we used the average normalized velocity is as the order parameter, defined as where T = 2000 is the simulation time for each experiment.
The error bar is defined as: whereφ is the standard deviation of the order parameter values, and n > 0 is the number of simulations for a given system size. The typical values of n were chosen as follows:

Results
We have compared the stability of VEM and some simple hierarchical flocking systems, such as CHM and SHM. According to our results (see Figure. 3 and Figure. 4), a simple dominance based system does not perform better than the much studied egalitarian flock.
In this paper we primarily report on our results concerning leaderfollower systems.
For more details about VEM versus CHM and SHM see hal.elte.hu/ vicsek/downloads/papers/Trieste-poster-JYN-TV-final.pdf Figure 3. Quantitative comparison on VEM and CHM (l ii = 0 for i ∈ 1, · · · , N ) in various of group sizes. c ij satisfies log-normal distribution (whose mean is 0 and standard deviation is 1) for the CHM system. Note that both Fig. 3 and 4. shows that the egalitarian system is more stable then the simple dominance-based hierarchical system. Figure 4. Quantitative comparison of the VEM and SHM (l ii = 0 for i ∈ 1, · · · , N ) for various of group sizes. c ij satisfies log-normal distribution(whose mean is 0 and standard deviation is 1) for the SHM system.

VEM vs DHM/DHMZ
4.2.1. Order Parameter As mentioned above, we shall associate the stability with the tendency of a flock not to break into parts under external perturbations. In order to measure the stability of the particle group, we used the following velocity correlation as the order parameter where j ∈ N i (t)/i means j ∈ N i (t) and j = i. This expression indicates the stability of the flock under different conditions. We did not use the average velocity as the order parameter, because the average velocity cannot give a right stability description when the system is divided into two or more coherently moving subgroups. Figure 5. quantitative comparison on VEM and DHM (l ii = 0 for i ∈ 1, · · · , N ) in various of group sizes. c ij satisfies log-normal distribution(whose mean is 0 and standard deviation is 1) for VDM system. Figure 6. quantitative comparison on VEM and DHMZ (l ii = 0 for i ∈ 1, · · · , N ) in various of group sizes. c ij satisfies log-normal distribution(whose mean is 0 and standard deviation is 1) for DHMZ system.

Results
We compare VEM with DHM and DHMZ separately. Figure. 5 shows the comparison results of the VEM and DHM, while Figure. 6 shows the comparison results of VEM and DHMZ. According to Figure. 5, we can see that VEM is a little bit better (better meaning that more ordered for the same amount of noise) than the DHM for different group size under lower noise, while it seems that DHM is a little better than VEM under larger noise. However, the advantage is so small that can be ignored. At the same time, it seems that there exists no relevant difference among the simulation results for various group sizes. Figure. 6 demonstrates that the DHMZ hierarchical model is significantly more stable than an egalitarian model for flocking, and the advantage is increasingly obvious as the size of the system increases. However, the stability gap between an egalitarian flock and a hierarchical flock stops increasing from 320 particles to 1280 particles. That is to say, the difference in the performances between the egalitarian system and the hierarchical system is increasing with system size, but only up to a given flock size and is not changing as the group size reaches a certain threshold. This is quite along the intuitive picture which suggests that hierarchy may play a relevant role only up to -at most -a few hundred of flock members.
Let us consider twenty-eight particles, for example, Figure. 7 shows some important scenarios during the flocking process (for a movie see Supplementary Material S1). The color bar indicates the weight of the contribution of the given particle. For example, the red particle is the strongest one, while the purple one is the weakest particle. The contribution values of these particles belonging to the DHMZ system satisfy log-normal distribution. That is, c ij , i, j = {1, · · · , N}, i > j satisfy log-normal distribution, with mean value 0 and standard deviation 1. At the same time, the contribution value of each particle belonging to VEM system is equal, and the sum of the contribution values of VEM system is equal to the sum of the contribution values of the DHMZ system. The first picture in Figure. 7 displays the initial state of all the particles at start. The second frame to the left depicts the structure of the VEM flock being less cohesive than the DHMZ (the right one). The third picture shows a key moment when VEM results in two separated groups while DHMZ results still in a single cluster. This can be taken as an important evidence to demonstrate that leader-follower hierarchical systems are more efficient regarding their stability against perturbations which the individual units are subject to.

Acknowledgments
This work was partly supported by the following grants: U.S. Air Force, Office of Scientific Research grant no.