Spatial flocking: Control by speed, distance, noise and delay

Fish, birds, insects and robots frequently swim or fly in groups. During their three dimensional collective motion, these agents do not stop, they avoid collisions by strong short-range repulsion, and achieve group cohesion by weak long-range attraction. In a minimal model that is isotropic, and continuous in both space and time, we demonstrate that (i) adjusting speed to a preferred value, combined with (ii) radial repulsion and an (iii) effective long-range attraction are sufficient for the stable ordering of autonomously moving agents in space. Our results imply that beyond these three rules ordering in space requires no further rules, for example, explicit velocity alignment, anisotropy of the interactions or the frequent reversal of the direction of motion, friction, elastic interactions, sticky surfaces, a viscous medium, or vertical separation that prefers interactions within horizontal layers. Noise and delays are inherent to the communication and decisions of all moving agents. Thus, next we investigate their effects on ordering in the model. First, we find that the amount of noise necessary for preventing the ordering of agents is not sufficient for destroying order. In other words, for realistic noise amplitudes the transition between order and disorder is rapid. Second, we demonstrate that ordering is more sensitive to displacements caused by delayed interactions than to uncorrelated noise (random errors). Third, we find that with changing interaction delays the ordered state disappears at roughly the same rate, whereas it emerges with different rates. In summary, we find that the model discussed here is simple enough to allow a fair understanding of the modeled phenomena, yet sufficiently detailed for the description and management of large flocks with noisy and delayed interactions. Our code is available at http://github.com/fij/floc.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 small shoals [1], tracking the individual coordinates of up to 2700 birds in flocks [2], and obtaining GPS track logs of homing pigeons flying together [3]. Initially, experimental and modeling efforts were focused on planar (2 dimensional) motion. Due to these efforts it is now well known that in 2 dimensions bacteria, insects, horses, and also humans display collective motion patterns [4][5][6][7]. Compared to planar motion, an agent moving in space can be kept aligned by a higher number of nearest neighbor interactors. At the same time, it has also more directions to turn away from the consensus of those nearest neighbors.
The most straightforward local rule that can describe the alignment of moving agents with their neighbors is to set each agent's direction of motion explicitly to the average direction of its neighbors [8,9]. Turning continuously toward the average direction of the neighbors is also possible [10,11]. More detailed mechanisms of the alignment include anisotropic interactions caused by elongated shapes [12][13][14][15], also combined with a frequent reversal of the direction of motion [16], the preference for movements in the horizontal plane (as opposed to vertical movements) [17], a viscous medium [18], friction among the agents and inelastic collisions [19,20], and sticking together [21]. While these rules can set the direction of motion for the agents, collision avoidance and cohesion (staying together) are also necessary for flock formation. To avoid the collisions of moving and interacting agents the simplest solution is to let all agents (in the model) have zero size [8]. A more realistic solution is a strong short-range repulsive interaction, in which the magnitude of the repulsion force becomes very high when two agents come too close. Finally, for keeping the group together two commonly applied modeling tools are the spatial confinement of the group (e.g., periodic boundaries) and a weak attraction that is turned on when distances between the agents grow.
The model that we discuss here focuses on controlling the speed of the agents individually. The speed of an agent is adjusted to the preferred speed with a rate that is proportional to the difference from the preferred value (see Fig 1). This modeling approach is realistic, becauseaccording to recent experimental and modeling evidence-individual speed control plays a key role in the formation and the stability of bird flocks and fish shoals [1,22,23]. Also, experiments and models for vibrated self-propelled hard disks have shown that the binary collisions caused by maintaining speed can align velocity vectors first locally, and then also across the entire system [24].
We investigate the effect of time delay and noise, too. Time delay is a common phenomenon caused by latent communication between agents, information processing cost, and inertial reasons [25,26]. Noise at all levels is also inherent to the communication and decisions of all Spatial flocking: Control by speed, distance, noise and delay moving agents in a dynamic system and can lead to transitions between behavioral patterns [26][27][28][29]. Regarding the combination of time delay and noise, simulation results in [30] showed that a system with noise and delay displays bistability of several coherent patterns. Here we investigate both aspects and show their fundamental dissimilarities.

Model: A minimal continuous description of spatial flocking
Here we investigate the 3-dimensional version of the generic model that was suggested in [31] and analyzed for 2 dimensions in [32]. The model is continuous in both space and time, and it contains N agents. The ith agent's position isr i ðtÞ ¼ ðx i ðtÞ; y i ðtÞ; z i ðtÞÞ and its velocity is v i ðtÞ ¼ ðv x;i ðtÞ; v y;i ðtÞ; v z;i ðtÞÞ. During its motion each agent continuously adjusts its speed, v i (t), toward a constant preferred v 0 value with characteristic time τ. To model collision avoidance, we apply pairwise radial repulsion among the two agents with a magnitude of F(r) = cr −2 as a function of their distance, r (for correct dimensions we set c = 1Nm 2 ). For simplicity, we set the mass of each agent to m = 1kg and their time constants for adjusting speed to τ = 1s. Based on the above, the equations of motion are (for i = 1. . .N): Note that the first term on the right hand side of Eq 1 points toward the ith agent's own direction of motion,ṽ i =jṽ i j. In other words, this simplified model separates collision avoidance (second term) from keeping the preferred speed (first term). We solve the equations of motion numerically by applying the midpoint method with an integration time step of Δt = 10 −3 s in a cube that has side length L = 50m and periodic boundary conditions in all three directions. For two particles the precision of the forward Euler method is sufficient. To integrate the equations of motion of many particles we apply the midpoint method. When the simulation is started we place all agents at random positions-but no pair of them is allowed to be closer than 0.6LN −1/3 -and set all speeds to the preferred speed, v 0 = 5m/s. We start the system either from a disordered state or an ordered state. When starting the system from the disordered state, we set the directions of the velocity vectors to N different, randomly selected directions at simulation time t = 0. When starting from the ordered state, we set all directions to the same, randomly selected direction at time t 0 = −10 5 , then simulate the system until t = 0, and after that start to log positions and velocities.
The two additional aspects of the model that we investigate are its responses to noise and time delays. We include noise into the model by adding a randomx vector to the right hand side of Eq 1. This vector is uncorrelated both in time and among agents, its direction is distributed uniformly in space, and its magnitude is a constant, ξ. If ξ = 1, then the autocorrelation of x is 1, therefore, the noise vector added to the r.h.s. of Eq 1 during a simulation update of length Δt is x ffiffiffiffiffi Dt p . Similary, for the first half-step of the midpoint update the amplitude of the applied noise vector is x ffiffiffiffiffiffiffiffiffiffi Dt=2 p . In addition to the above, we compute interactions with a distance-based upper cutoff: in the simulations two agents interact only if their distance is below a fixed cutoff radius, R = 10m. Finally, for efficient computation, we apply a 3-dimensional grid with side length R and search for an agent's interactors-i.e., other agents closer than Ronly within its own grid cell and the neighboring 26 grid cells.

Results: Alignment and flocking
A necessary condition for the emergence of a single stable aligned group containing the majority of all agents is the ability of small groups to align. More specifically, the most frequently studied necessary condition is whether the alignment of two agents increases during their encounter, i.e., when they come close and depart as in the top panel of Fig 2.

Alignment of two agents in a symmetric encounter
As a simple-yet nontrivial-case of a two-agent encounter in space, we analyse the scenario when two identical agents move such that their velocity vectors are permanently mirror images of each other with respect to the two agents' center of mass (see Fig 3a). To parametrize this 3-dimensional motion of the two agents, we fit two parallel planes to the skew lines defined by their velocity vectors at t = 0. The distance of these two plains is d, and both agents have a speed of v 0 . When looking at the two agents from the direction perpendicular to the two planes, at t = 0 the distance of the two agents is 1,000m and the angle of their velocity vectors is φ. As in Eq 1, radial repulsion is F(r) = cr −2 (with c = 1Nm 2 for correct dimensions). We compute this pair encounter by integrating Eq 1 with the forward Euler method with an integration time step of dt = 10 −5 , and the parameters v 0 = 1m/s, τ = 1s and ξ = 0. We stop the integration when the horizontal distance of the agents reaches 1,100m. Finally, we obtain an estimate of the error of the numerical integration: we test how the integration time step parameter (dt)  influences the calculated change of this two-particle system's total momentum (ΔI) during the encounter. We find that the maximal difference of the final ΔI values between the dt = 10 −4 and dt = 10 −5 cases is 1.8 × 10 −4 . This difference is significantly smaller than the ΔI values in Fig 3 that are all above 10 −3 .
According to Fig 3b, we obtain the following result for the investigated symmetric collision. If the initial angle, φ, exceeds a threshold value-which is between 40 and 50 degrees depending on the distance, d, of the two particles' initial planes,-then the total momentum of the two particles increases. In other words, Eq 1 implies that two particles arriving symmetrically at a  large angle will depart at a smaller angle: they will become more parallel. With many interacting and moving particles the interactions are usually not pairwise and typically not symmetrical. Therefore, the results displayed in Fig 3b are not more than a significant microscopic result pointing at the possibility of the large-scale alignment of all moving agents. This possibility is tested in Sec below.

Alignment of many collectively moving agents
To test whether Eq 1 can indeed lead to a stable coherent motion of all agents, we start the system from a fully disordered setup and check if it reaches a stable ordered state (see Fig 4). We simulate the motion of agents-i.e., perform the midpoint integration of Eq 1-within a cube that has side length L = 50 and periodic boundary conditions. Parameters are v 0 = 5 (preferred speed), τ = 1 (time constant for reaching the preferred speed), R = 10 (interaction cutoff radius), dt = 10 −3 (time step of integration) and ξ = 0 (no noise). At simulation time t = 0 the speed of each agent is v 0 and velocity vectors point to independently selected random directions. Starting the system with this scenario, we measure overall ordering through the quantity Spatial flocking: Control by speed, distance, noise and delay called efficiency: Note that in theory, Eq 1 does allow for E(t) to exceed 1, if many agents have jṽ i j > 1. However, in practice, even with 2 agents this is a very rare case.   numerically measured E values: We find that growing system size (and constant density) the transition from the ordered to the disordered state remains rapid. For practical applications this means that-in the absence of disturbances-once the transition to the ordered state starts, it is very hard to stop, and the resulting ordered state is very stable. To go into more detail, in the following we investigate the stability of the arising ordered state to two of the most common disturbances: noise (errors) and delay (also called time lag).

Stability of the ordered state to noise
Errors and small changes (also called fluctuations or noise) are ubiquitous in natural and social phenomena. In most cases, a low amount of noise still allows ordering, however, when noise amplitudes grow very large, order is destroyed and the system becomes disordered. A known exception to this rule is the phenomenon of "freezing by heating", when-compared to the case without noise-an intermediate amount of noise can lead to a new ordered structure [35].
In Fig 6 we test whether it is easy to destroy the ordered state (high E) of the agents of Eq 1 by Spatial flocking: Control by speed, distance, noise and delay increasing the initially low amount of noise to a high value. With the parameters selected in Sec we find that the noise level ξ = 50 can already keep an initially disordered system from reaching order, however, it is not sufficient to destroy order in an initially ordered system. We conclude that the system shows hysteresis, which is similar to the 2-dimensional version of the model [32]. Finally, note that this behavior stabilizes the ordered state, and at the same time makes transitions fast.
In this final paragraph of Sec, we add a note regarding how we start the system from the ordered state. For the curves labeled "Start: ordered" in Fig 6 we start each of the 30 simulations (with 30 different random seeds) at t start = −100,000s by setting the velocity vectors of the agents fully aligned and their coordinates randomly with a minimal distance. Here minimal initial distance means-similar to the case of the disordered initialization in Sec-a lower bound of 0.6LN −1/3 for the agent-agent distance. First, starting at t start , we run the simulation until t = 0s to reach spatial ordering. After the equilibration interval that starts at t start and ends at t 0 , we run the numerical integration of the equations of motion until t = 10 5 or t = 10 6 . Note that in Sec we do not apply the equilibration technique described here.

Stability of the ordered state to delay (time lag)
In most natural, social and technological phenomena the delays of interactions play a crucial role in shaping the emerging collective behavioural patterns. Stable formations of collective motion can usually emerge only if the distance traveled by an agent over the delay time interval with its characteristic relative speed (relative to nearby agents and obstacles) remains safely below the agent's distance to those other agents and obstacles. For the current model we set a single delay time interval, t d , that affects in an identical way how an agent responds to any internal or external effect. For the numerical integration, we implement delayed interactions by setting the time lag, t d , of all effects (forces) to a multiple of the numerical integration time step, Δt. Next, we replace all forces on the right hand side of Eq 1 by the values of the same forces measured t d time before the current time. Note that these forces include self-propelling, therefore, a particle's response to changes in its own speed are delayed as well.
We find that-similar to noise-interaction delays act differently on the ordered and disordered states of the system (see Fig 7). First, the two states respond to time delays at different time scales. Whereas the ordered state responds to a wide range of interaction delays on the same short time scale, the already slow disorder ! order transition is further significantly delayed even by small amounts of the interaction delay. Second, even large amounts (up to t d = 50ms) of the interaction delay cannot entirely destroy the ordered state. Third, based on our current results with t d = 0, 10, 20 and 50, we conclude that ordering is lost as a linear function of the amount of interaction time lag, t d .

Discussion and outlook
In this paper, we investigated a minimal continuous model of 3-dimensional collective motion. The model contains continuous adjustment of particle speed to a preferred value, pairwise radial repulsion for collision avoidance, and an effective weak attraction (periodic boundaries). We found that the combination of these three model components is sufficient for stable spatial ordering, and beyond these three no further model components are necessary. After investigating the model on the microscopic level, we found that for the majority of symmetric twoagent encounters the total momentum of the two agents increases. Regarding the macroscopic level, we found that the transition from disorder to order is fast for both small and large system sizes, which is in good agreement with previous results [27]. In the minimal continuous model that we investigated we found also that if the noise intensity is above a threshold value, then the system cannot reach the ordered state, similarly to results reported in [36]. Note that at the start of the simulation the efficiency drops slightly and then grows again to reach E % 1. To explain this effect, recall that this simulation is started with directionally ordered agents that are spatially not yet ordered. The slight drop of the efficiency shows as the particles attain their spatial order (close to crystalline ordering of coordinates) while temporarily losing some of their directional alignment. https://doi.org/10.1371/journal.pone.0191745.g007