We developed an agent-based model of the emergence of coordinated movement in a group and implemented it in the simulation program P-Flock 2.0, which was written in C and Delphi and run in Windows. Both the model and the computer program are continuations of previous work on the modelling and simulation of spatial behaviour in groups (Beltran, Salas & Quera 2006; Quera, Beltran, Solanas, Salafranca & Herrando 2000; Quera, Solanas, Salafranca, Beltran & Herrando 2000). In the current model, agents adapt their movements locally through repeated interactions as they encounter other agents; adaptation is accomplished by applying local rules to reward the agents' movements that fulfil their predictions, and to penalize those that do not. As a result, under certain conditions, these dyadic adaptations lead to the emergence of coordinated movement in the group, a global phenomenon that is not specified in the local rules. A snapshot of a simulation with the P-Flock 2.0 program is shown in Figure 1. Figure 1. Snapshot of a P-Flock simulation with 50 agents, showing their positions and headings in a toroidal world. At time unit 1,600, a loose flock emerged; the program detected a possible leader (agent 39, shown in yellow). Agents with different perceptual features (see below) are displayed using different colours. A flocking index time series is displayed in the bottom section (see below). World and Agents

The agents move in a two-dimensional, discrete world that is either toroidal or closed. The world is composed of cells or patches, and one cell can be occupied by only one agent at a time. Each agent i is identified by its coordinates (x i (t), y i (t)) and heading α i (t) at time t. The heading of an agent at time t is defined as the vector connecting its location at t-1 with its location at t, and is expressed as the anti-clockwise angle between that vector and the X axis, α i (t) = tan-1 [ (y i (t) - y i (t-1)) / (x i (t) - x i (t-1))].

Agents have a perceptual field, which is defined as a circular sector whose centre is the agent's current location and which is bisected by its current heading vector; radius r i and angle θ i of the circular sector are the depth and scope of agent i's perceptual field, respectively, and are the model parameters that can be modified (see Figure 2). Agent movement is restricted to its current local neighbourhood of cells; Moore or Von Neumann neighbourhoods can be defined with several possible diameters (see Figure 3). Figure 2. An agent's coordinates at times t-1 and t, its heading α i (t), and its perceptual field (in white), defined by radius r i and angle θ i . Figure 3. Moore and Von Neumann neighbourhoods and possible diameters. The agent's current position is indicated by a black square, and its current movement is limited by its neighbourhood. Ideal Distances and Dissatisfaction

A matrix of ideal distances among the agents is also defined. The ideal distance from agent i to agent j at t, D ij (t), is the distance at which agent i wants to be from agent j at that time; ideal distances are not necessarily symmetrical. As agents encounter other agents, their ideal distances can change as a function of the outcome of their interaction. Agents' local goal is to minimize their dissatisfaction and to that end they can move within their current neighbourhood; an agent's dissatisfaction at time t, U i (t), is defined as a composite function of the absolute differences between its ideal distance from every other agent that it currently perceives and the real distance from them: (1) where d ij (t) is the Euclidean distance from agent i to agent j at t, d ij (t) = [(x i (t) - x j (t))2 + (y i (t) - y j (t))2]1/2 (2) and the sum is for all the agents currently perceived by agent i at t (subset Z i (t)). The sum is divided by the diameter of the world, or the maximum possible real distance (m), and by the cardinal number of Z i (t); dissatisfaction therefore has to be between 0 and 1 (for details, see Quera, Beltran, Solanas, Salafranca & Herrando 2000). The maximum possible real distance is a function of the size and shape of the world; for a rectangular world of size H × V, m = [H2 + V2]1/2 if it is non toroidal, and m = ([H2 + V2]1/2)/2 if it is toroidal.

At time t, agent i estimates its ideal or possible future dissatisfaction for time t+1 for all candidate positions p within its current neighbourhood, by computing the distances from them to the current positions of the agents it is perceiving: (3) where d ijp (t) is the real distance from position p within agent i's neighbourhood to agent j's position at time t (see Figure 4). Agent i then decides to move to that position within its neighbourhood for which it estimated U' ip (t+1) to be the lowest value. If several positions share the minimum ideal dissatisfaction, the agent decides on the move that requires the least change to its current heading. At time t+1, the other agents could also have moved, and therefore the ideal dissatisfaction that agent i estimated for time t+1 might not necessarily have been attained. At each time unit the agents make the decision to move within their respective neighbourhoods simultaneously and independently. However, if two or more neighbourhoods overlap, then some of the cells in those neighbourhoods can be candidates for more than one agent if they happen to provide the minimum dissatisfaction for all of them. At each time unit, agents' priorities to move are sorted randomly. The agents with lower priorities will then only be able to move to those cells not chosen by agents with higher priorities that provide less than the minimum dissatisfaction. Figure 4. How an agent estimates its possible dissatisfaction at a given time. For every candidate position within its current neighbourhood (in this example, a Moore neighbourhood with a diameter of 3 and nine candidate positions), agent 1 computes the real distances from itself to the current positions of the other agents it perceives. These real distances are compared with the ideal distances as per the Equation 3 and a dissatisfaction is obtained for each candidate position. Flock Synthesis Rules

In this model, ideal distances between agents vary as a consequence of the outcomes of their interactions, according to a set of reward rules, which we call Flock Synthesis Rules (FSRs) because a flock may emerge when they are applied locally, massively and repeatedly. Initially, agent i moves at random without considering any other agents; when it perceives agent j for the first time and their real distance apart is less than or equal to a critical value, P D , the FSRs for agent i are activated with regard to agent j, and the ideal distance from agent i to agent j is set so it is equal to a uniform random value between d ij , their current real distance apart, and m, the maximum possible real distance. From that time on, ideal distance D ij undergoes two different kinds of change: smooth change and abrupt change.

A smooth change is caused by agent i when it adapts to agent j's movements. As mentioned above, at each time step, agents i and j move to the positions in their respective current neighbourhoods that minimize their dissatisfaction. Before moving (time t-1), agent i predicted that its real distance from agent j at time t would be d' ij (t), assuming that agent j would not move. After moving, agent i evaluates the difference between its prediction, d' ij (t), and its previous real distance from agent j, d ij (t-1): c ij (t) = d' ij (t) - d ij (t-1) (4)

From the point of view of agent i, its movement was an "attempt to approach" agent j if c ij (t) < 0, or an "attempt to move away" from agent j if c ij (t) > 0. After moving, agent i also evaluates the discrepancies between the current real distance, d ij (t), the predicted distance, d' ij (t), and the ideal distance it wanted to keep from agent j, D ij (t-1): u ij (t) = |d ij (t) - D ij (t-1) | - | d' ij (t) - D ij (t-1) | (5) which is in fact the difference between the two partial dissatisfactions for agent i with respect to agent j: the actual one and the one it predicted it would have before moving. If u ij (t) = 0, we can say that agent i has momentarily adapted to agent j, because the dissatisfaction it expected is fulfilled; if u ij (t) < 0, agent i "overadapted" to agent j, because the dissatisfaction after moving is less than expected; and if u ij (t) > 0, agent i "underadapted" to agent j, because the dissatisfaction after moving is greater than expected.

The ideal distance from agent i to agent j is increased if agent i either attempted to approach and underadapted (c ij (t) < 0 and u ij (t) > 0) or attempted to move away and overadapted (c ij (t) > 0 and u ij (t) < 0); in the first case, increasing the ideal distance can be seen as a penalty, as the attempted approach was excessive, while in the second case, it can be seen as a reward, as the attempted move away was successful. On the other hand, the ideal distance is decreased if agent i either attempted to approach and overadapted (c ij (t) < 0 and u ij (t) < 0) or attempted to move away and underadapted (c ij (t) > 0 and u ij (t) > 0); in the first case, decreasing the ideal distance can be seen as a reward, as the attempted approach was successful, while in the second case, it can be seen as a penalty, as the attempted move away was excessive. Finally, the ideal distance either remains the same if agent i momentarily adapted (u ij (t) = 0), or it is randomly increased or decreased if agent i overadapted or underadapted but attempted neither to approach or move away (c ij (t) = 0 and u ij (t) ≠ 0). Change in the ideal distance at time t is thus: D ij (t) = (1 + kP C )D ij (t-1) (6) where P C is a parameter of the model that modulates the rate of change (0 < P C < 1), and k = +1, 0, or -1, for increase, no change and decrease, respectively (see Table 1 and Figure 5). Table 1: Decision Table for the Flock Synthesis Rules Underadapted

u ij (t) > 0 Adapted

u ij (t) = 0 Overadapted

u ij (t) < 0 Attempt to approach

c ij (t) < 0 Increase (k = 1)

Penalty No change (k = 0) Decrease (k = -1)

Reward No attempted move

c ij (t) = 0 Random

Increase or Decrease No change (k = 0) Random

Increase or Decrease Attempt to move away

c ij (t) > 0 Decrease (k = -1)

Penalty No change (k = 0) Increase (k = 1)

Reward Figure 5. An example of agent i adapting to agent j and updating its ideal distance. The agents' positions are indicated by green cells at the centre of their current diameter-3 Moore neighbourhoods. At time t-1, agent i decided to move to the upper left cell in its neighbourhood because it happened to provide minimum dissatisfaction; agent i computed its real distance to agent j, d ij (t-1) (shown in green), and the distance it predicted for time t, assuming that agent j would not move, d' ij (t) (shown in red). At time t, both agents moved (as indicated by the brown arrows, which are the agents' current headings) and the real distance between them is d ij (t) (shown in blue). Agent i's ideal distance to agent j at t-1, D ij (t-1), is shown in brown. As c ij (t) = d' ij (t) - d ij (t-1) > 0, agent i attempted to move away from agent j. As |d ij (t) - D ij (t-1)| < | d' ij (t) - D ij (t-1)|, u ij (t) < 0, and agent i (slightly) overadapted to agent j when it moved. Consequently, agent i's movement with respect to agent j is rewarded and the ideal distance is increased by factor P C .

An abrupt change in the ideal distance from agent i to agent j occurs when that distance remains constant or changes cyclically in small amounts at consecutive time steps for a certain period of time S ij (T), T being the number of times the abrupt change occurred. This reflects agent i's tolerance towards being 'stagnated' with respect to agent j, i.e. remaining stable or relatively stable with respect to it. Initially, S ij (1) = P S , which is a parameter of the model. When the tolerance limit is reached, the ideal distance is increased abruptly by a factor equal to P C · S ij (T), i.e. D ij (t) = D ij (t-1) + P C S ij (T) (7) and it remains constant for the next E ij (T) time units, where E ij (T) is agent i's tolerance for being 'exiled' from agent j and T is the number of times the exile occurred. Initially, E ij (1) = P E , which is a parameter of the model. By increasing agent i's ideal distance from agent j abruptly when it stagnates, agent i is given an opportunity to adapt to new situations while it is exiled from agent j. When the exile time is over, the ideal distance is decreased abruptly by a factor equal to P C · S ij (T) + P D , i.e. D ij (t) = D ij (t-1) - (P C S ij (T) + P D ) (8) and agent i resumes the smooth change phase with respect to agent j. Every time T that agent i is exiled with respect to agent j, its tolerance to stagnation is increased slightly by factor P C : S ij (T) = (1 + P C )S ij (T - 1) (9) and every time it resumes the smooth change phase, its tolerance to exile is decreased slightly by that same factor: E ij (T) = (1 - P C )E ij (T - 1) (10)

Thus, agent i progressively adapts to agent j by tolerating longer stagnation and shorter exile periods from it. The model parameters are listed in Table 2. Table 2: Parameters of the Flock Synthesis Rules P D Critical real distance at which the FSRs are activated P C Change rate (0 < P C < 1) for increasing and decreasing ideal distances, and stagnation and exile tolerance times P S Initial value of the stagnation tolerance time P E Initial value of the exile tolerance time

It should be noted that the model does not require every agent to interact with every other agent. That is, whenever agent i perceives agent j, only the FSRs for i towards j are activated. Agent j may never have perceived i, or certain pairs of agents may never have perceived each other. Thus, repeated application of the dyadic FSRs make the group move in a coordinated way, even if some, or many, of the agents never interacted with each other. We must stress that, even if an agent perceives several agents simultaneously, the FSRs are applied in a dyadic and independent way. That is, if at time t agent i perceives agents j 1 , j 2 , and j 3 , the FSRs are independently applied three times, for agent i towards j 1 , j 2 , and j 3 , respectively; the decision for agent i towards j 1 resulting from Table 1 has no effect on the decision for agent i towards j 2 , and so on.

The pseudocode of the main functions of the P-Flock program is shown in the Appendix. The C code includes specific functions (not shown there) for computing Euclidean distances on a toroidal surface and agent dissatisfaction, and for checking which agents perceive other agents at every simulation step. The program reads a parameters file, performs the simulations, and saves the agents' coordinates, headings, ideal distances, and flocking and leadership indices (see below) to a text file. The simulation can then be played, paused, slowed down and played backwards (e.g. Figure 1). P-Flock 2.0 runs on Windows and can be downloaded from http://www.ub.edu/gcai/ (go to download in the main menu). A demo video can also be downloaded from that web site.