Collective Motion in a Network of Self-Propelled Agent Systems

Collective motions of animals that move towards the same direction is a conspicuous feature in nature. Such groups of animals are called a self-propelled agent (SPA) systems. Many studies have been focused on the synchronization of isolated SPA systems. In real scenarios, different SPA systems are coupled with each other forming a network of SPA systems. For example, a flock of birds and a school of fish show predator-prey relationships and different groups of birds may compete for food. In this work, we propose a general framework to study the collective motion of coupled self-propelled agent systems. Especially, we study how three different connections between SPA systems: symbiosis, predator-prey, and competition influence the synchronization of the network of SPA systems. We find that a network of SPA systems coupled with symbiosis relationship arrive at a complete synchronization as all its subsystems showing a complete synchronization; a network of SPA systems coupled by predator-prey relationship can not reach a complete synchronization and its subsystems converges to different synchronized directions; and the competitive relationship between SPA systems could increase the synchronization of each SPA systems, while the network of SPA systems coupled by competitive relationships shows an optimal synchronization for small coupling strength, indicating that small competition promotes the synchronization of the entire system.


Introduction
Groups of animals sometimes display fascinating collective motions [1][2][3][4][5][6] in which animals move in the same direction, such as schools of fish can move in a rather orderly fashion or change direction amazingly abruptly [2], and flocks of birds can fly as a uniformly moving group. Moreover, collective motion are common both in living and non-living worlds, ranging from biology [7], ecology [8], climate [9], society to technology [10,11] and even art [12], the studies of which can help us understand the nature world and improving infrastructure systems in the man-made world. A simple model proposed by Vicsek et al. [13], unveiling the collective motion and phase transition of self-propelled agent systems (SPA), can be potentially applied to man-made systems, such as distributed sensor networks [14], unmanned aerial vehicles [15], underwater vehicles [16], altitude alignment of satellite clusters [17], and many more. Underlying the behavior of collective motions, synchronization process is the cause of such fascinating phenomenon. Various models have been proposed to mimic synchronization processes and many strategies has presented to improve the properties of synchronization [18,19]. Most of these studies are focus on isolated networks.
Increasing evidence shows that one system may interact or couple with other systems, such as: different social networks (e.g., Facebook, Twitter) are interconnected with each other because the nodes in different networks share the same actors [20]; transportation networks (e.g., buses, airplanes) are coupled with each other since the nodes in each network are in the same geographic locations [21]; the infrastructure systems (e.g. communication networks and power grid) are interdependent because the nodes in one network support the nodes in another network [22,23]. Indeed, some realistic systems such as neuronal system could also be self-propelled or self-adjusting due to the autapse connection to neuron, that the collective behaviors of neurons could be regulated by autapse driving when continuous pulse or traveling wave is induced [24,25]. All these examples unveil that real systems usually interact with each other, leading to the emerging new field in network science, interdependent networks, interconnected networks, a network of networks, multi-layered networks, multiplex networks and many more [26][27][28][29][30][31][32]. In real scenarios, different self-propelled agent systems also coupled with each other and form a network of SPA systems. For example, unmanned aerial vehicles may work with underwater vehicles to achieve some tasks, exhibiting a symbiotic relationship between these SPA systems; a flock of birds might have predator-prey relationships with a school of fish; and schools of fish may compete with each other for sharing the same food. Such coupling relations between different SPA systems could influence the synchronization of these systems. However, no such model exists for showing how the interdependence between different SPA systems influence the synchronization of a network of SPA systems. Studies of synchronization of network of SPA systems enable us to design high efficient systems of different coupled unmanned vehicles or robots.
In this work, we propose a model to show the collective motion in a network of SPA systems, generalizing the Vicsek model [13] (VM). In the model, we introduce a coupling strength β 2 [0, 1] denoting the ratio of a node's the radius in its coupled system to that in its current subsystem. Furthermore, we construct three networks of SPA systems coupled by three types of interaction: symbiosis, predator-prey, and competition. Besides, we find that (1) there exits an optimal coupling strength in the coupled systems with symbiotic and predator-prey relationship to achieve optimal synchronization for each subsystem and the entire system; (2) furthermore, increasing the radius and absolute velocity and decreasing of system size could increase the optimal synchronization of every subsystem and the entire systems; (3) in the systems coupled by competitive relationships, increasing the coupling strength could increase the synchronization of each subsystem but decrease the synchronization of the entire system.

Model
Our model contains N coupled subsystems of self-propelled agents with size n k in subsystem k, where the agents in each subsystem move in a L × L square with a same constant speed towards different directions. Initially, the agents are randomly distributed on the L × L square plane, and their directions are also uniformly randomly distributed in the interval (0, 2π). At each time step, the direction of each agent is determined by the average directions of all the agent within a circle centered at the given agent with a radius R (Fig 1 System A). At each step t, the position of a specific agent is updated according to where x w (t) is the two dimensional vector of position of agent w at time t, v 0 is the absolute speed of each agent and θ w (t) is the direction of agent w at time t. Then its direction is updated following e iy w ðtþ1Þ ¼ e iDy w ðtÞ P j2G w ðtþ1Þ e iy j ðtÞ k P j2G w ðtþ1Þ e iy j ðtÞ k ; ð2Þ In system A, the position of node w is (x A,w ,y A,w ), and it has 7 neighboring agents within a circle with a radius being R [Eq (3)], whose position projecting on system B is also (x A,w ,y A,w ), and it has 4 neighboring agents in system B [Eq (4)]. Accordingly, we could identify the neighborhoods of all nodes in both systems. where kÁk is the standard norm [33] defined by k(z 1 , z 2 , . . ., z l )k = (|z 1 | 2 + |z 2 | 2 + . . . + |z l | 2 ) 1/2 , 4θ w 2 [−η, η] representing the white noise, e iy w ðtÞ is a unit directional vector, and Γ w (t+1) is the set of neighbors of agent i at time step t+1, defined as G Aw ðtÞ ¼ fjjðx A;w ðtÞ À x A;j ðtÞÞ 2 þ ðy A;w ðtÞ À y A;j ðtÞÞ 2 R 2 g; ð3Þ and G Bw ðtÞ ¼ fjjðx B;w ðtÞ À x B;j ðtÞÞ 2 þ ðy B;w ðtÞ À y B;j ðtÞÞ 2 R 2 g ð 4Þ for system A and B respectively. In order to measure the synchronization of the system, an order parameter is introduced as [13,34]: A larger value of V αk indicates that the subsystem k shows a better synchronization, and when V αk = 1 all the agent are moving towards the same direction. A system of two coupled SPA systems (A and B) is shown in Fig 1, where each SPA system forms a subsystem. An agent A w in system A moves inside a L × L square and its position is (x A,w , y A,w ). The neighborhood of node A w in system B is defined as and the neighborhood of node B w in system A is defined as where r = β × R. Especially, (i) when β = 0, the two SPA systems are separated from each other corresponding to the case of two isolated SPA systems; (ii) when β = 1, the entire system fully mixed with each other forming a new subsystem with the density being the sum of these two SPA subsystems. When β > 0, an agent w in system A updates its direction according to both the agents in Γ Aw and O AB,w together, and updates its location according to Eq (1). Furthermore, different impact of an agent from one system on an agent from another system exists according to different relationships between these two SPA systems, which will be discussed as follows.
(1) Symbiotic relationship A symbiotic relationship between two coupled systems benefit the agents of both systems [35]. The agents of one system have positive impact to the synchronization of the agents in another system, leading a result that an agent w in system A changes its direction according to the According to Eq (8), one can obtain the direction changes of agents in system B symmetrically. Note that this model can be potentially applied to the systems where both subsystems intend to achieve a global alignment, such as the system of unmanned aerial vehicles and underwater vehicles. Furthermore, this model is similar to the classic model of interdependent networks [21] where two nodes from different networks depending on each other because they support each other showing a symbiotic relationship.
(2) Predator-prey relationship In the model where SPA systems coupled by predator-prey relationship [36], without lose of generality, we assume that system A is composed by the predators, and system B is composed by the preys. Then each agent of system A updates its direction by the average directions of all its neighbors in both systems A and B, because predators intends to synchronize with the preys; and each agent in system B updates its direction using the average directions of all its neighbors in system B and average of all its neighbors' opposite directions in system A, since preys avoid to synchronize with predators. respectively. Note that this model can be potentially applied to coupled systems such as a flocking of birds/ducks interacting with a schooling of fish, where the birds/ducks eat the fish.

(3) Competitive relationship
The competitive relationship between two SPA systems has negative effects on both since the food, space or other resources that they competing for are limited [37]. Thus, the direction of each agent is obtained by the average directions of the neighboring agents of its own subsystem and opposite directions of the neighboring agents in the other subsystem, whose mathematical expression is: According to Eq (11), one can obtain the direction changes of agents in system B symmetrically.

(4) Network of self-propelled agent systems
In real world, usually more than two SPA systems are coupled with each other, which is also addressed in this work. As shown in Fig 2, there are 6 SPA subsystems coupled together by different relationships between each pair of subsystems. In such a network of SPA systems, different subsystem may converge to different level of synchronization, thus we define a vector of order parameter V s as where V α is the synchronization degree for the entire system, which denotes the synchronization of all the agents of n subsystems. V αw (w = 1, 2, . . ., N) represents the synchronization of system w. The order parameter not only represents the degree of each subsystem, but also contains the synchronization of the entire system. Considering the situation without noise (i.e. noise amplitude η = 0), we simulate the synchronization of networks of SPA systems coupled by three different relationships, and it can  help us understand the behavior of a network of n SPA systems. Moreover, we also investigate the effect of the four important parameters on the synchronization of a network of SPA systems, namely the coupling strength β, the radius R, the absolute speed v 0 , and the system size n.
In the simulations, we use the periodic boundary condition [13] and the fixed density [1] and the density of the system is ρ = n k /L 2 = 1. If there is no special statement, we use the parameters as R = 0.3, v 0 = 0.1, n k = n = 200 for each subsystem and all the results are averaged over 400 realizations.

Network of SPA systems coupled by symbiosis relationship
In this section, we show that how these three factors: the system size n, the absolute velocity v 0 , and the radius R, affect the synchronization of two SPA subsystems coupled by symbiotic relationships.
In Fig 3(a), the synchronization V α is a function of the coupling strength β under three different values of system size n. The results demonstrate that V α increases and then decreases with coupling strength β increasing, indicating that there exists an optimal value of β o and when β = β o the system achieves optimal synchronization V so . As shown in Fig 3(b), the optimal synchronization decreases with the system size for both the subsystems and the entire system.
As shown in Fig 4(a), V s increases and then decreases showing a peak as the coupling coefficient β increasing. As shown in Fig 4(b), the optimal synchronization V so is an monotonously increasing function of R and each subsystem and entire system reach fully synchronized when R > 0.7. We also find that the optimal coupling strength β o decreases as R increasing in each subsystem and entire system [Fig 4(c)]. Furthermore, as shown in Fig 4(d), the optimal synchronization V so increases as v 0 increasing, indicating that the absolute velocity v 0 could increase the synchronization of each subsystem and entire system. In summary, as we know, when β = 0, the system is equivalent to two isolated systems with the same density ρ = 1; and for the special case that β = 1, the network is equivalent to each subsystem with density ρ = 2. However, the case β = 1 does not correspond to the optimal synchronization. Counter-intuitively, we find that the coupled system with symbiotic relationship shows an optimal value of β, indicating that we can design optimal coupled system by choosing suitable value of coupling strength.

Network of SPA systems coupled by predator-prey relationship
Different SPA systems may coupled by the predator-prey relationships, and its synchronization might be influenced by the properties of each SPA system and their relationships. Next we present how the system size n, the absolute velocity v 0 and the radius R affect the synchronization of the network of SPA systems coupled by predator-prey relationships. We denote V α1 and V α2 as the synchronization of the subsystem of predators and the subsystem of preys, the behaviors of which under different system size n as the coupling strength α increasing is showing in Fig 5(a). Both the degree of synchronization of these two subsystems present a peak as the coupling strength α increasing, as shown in Fig 5. The peak value V so decrease as the system size n increasing.
We next study the effect of absolute velocity on synchronization for the paired SPA system with predator-prey relationships. In Fig 6(a), the optimal synchronization V so for each subsystem and the entire system varies with the absolute velocity v 0 . The results demonstrate that, V so  monotonously increase with the velocity v 0 increasing, which implies that the synchronization is significantly improved when the velocity v 0 increases. In particular, the synchronization converges to full synchronization when velocity v 0 is large. Besides, for a fixed velocity v 0 , the synchronization of each subsystem is better than the case of the entire system. To further study the difference between the synchronization of each subsystem and the entire system, we treat the corresponding optimal value β o as a function of the absolute velocity v 0 [Fig 6(b)] and find that optimal value β o decreases quickly as the absolute velocity v 0 increases for each subsystem and the entire system. Especially, for a fixed velocity v 0 , the optimal value β o of each subsystem is greater than that of the entire system.
In Fig 7(a), the optimal synchronization V so can be achieved with the radius R for each subsystem and the entire system with various radius β. For a fixed β, V so is an increasing function of radius R, which implies that the synchronization is significantly improved when the radius R increases. In particular, when R is large enough each subsystem achieve a full synchronization but the entire system goes to a constant level of synchronization (but less than 1), indicating that the predator subsystem and the prey subsystem converges to different directions as illustrated in Fig 1 as well. To further study the difference between the synchronization of each subsystem and the entire system, we regard the corresponding optimal value β o as a function the radius R in Fig 7(b) and find that optimal value β o decreases as the radius R increases for each subsystem and the entire system.

Network of SPA systems coupled by competitive relationship
We study the synchronization of the paired SPA systems coupled by competitive under different system size n, radius R, and different values of velocity v. As a function of the coupling strength for different system size n, the synchronization of subsystem A V α1 monotonically increases as the coupling strength β increases, as shown in Fig 8(a). On the contrary, when we consider synchronization of the entire system V α as a function of the coupling strength β, we can see that V α shows a peak at small β, indicating that tiny competition can increase the synchronization of the entire system. This result is as surprising as the result presented in Fig 8(b) that small noise can improve the synchronization of a system. Besides, we show the synchronization V s as a function of the absolute velocity β for each subsystem and the entire system [Fig 8(c)]. Our simulation results indicate that, the difference between each subsystem and the entire system synchronization changes sharply as the coupling strength β increases. In particular, the total synchronization of the entire system is significantly less than each subsystem when the coupling strength β increases. The simulation results demonstrate V s of each subsystem is better than the case of the entire system, since the synchronization of one system benefits the synchronization of another system, but the synchronized direction of each subsystem is different from each other. In order to show the optimal synchronization V so as a function of system size n, we perform the the numerical simulations in Fig 8(d). Our simulation results suggest that, for each subsystem and the entire system, the optimal synchronization V so maintain a constant value as the system size n increases.

Conclusion
SPA systems usually interact or depend on each other, forming networks of SPA systems. In this work, we propose a model to mimic the synchronization of the network of SPA systems coupled by different relationships: symbiosis, predator-prey, and competition. The synchronization process in the networks of SPA systems could be influenced by both the properties of single SPA systems (the size of each subsystem n and the absolute velocity) and the relationship between different SPA systems (the coupling strength β and the radius R). Moreover, the networks coupled SPA systems coupled by different relationships show different behaviors: (1) the system coupled with symbiotic relationships shows a complete synchronization, as each subsystem reaching a complete synchronization, under high absolute velocity or high radius; (2) the system coupled by the predator-prey relationships shows an optimal but not complete synchronization, while each subsystem can arrive at a complete synchronization; (3) the system of coupled competitive SPA systems shows an optimal synchronization for small coupling strength between SPA systems, while no optimal synchronization for each subsystems. These interesting results can significantly improve our understanding of the synchronization principles of complex systems.