Effects of topological structure and destination selection strategies on agent dynamics in complex networks

A complex network comprises numerous nodes connected by edges. Elucidating the effect of agent decision-making on complex networks is of utmost importance. In the real world, complex networks exist on which agents travel, as exemplified by traffic transport systems [1], socio-physical networks [2] such as crowds, virtual networks on the internet, and biological brain networks. Here, agents represent objects such as parcels, self-driven objects such as humans or vehicles, or the units of data or signals. In the case of agents' self-intelligence and that of receiving instructions from a manager, agent destination selection strategies can significantly affect the overall network behavior. Particularly, avoiding congestion [3–5] can result in different agent mobilities. An example is provided in our previous study [6], which explored the effect of agents avoidance congestion on agent mobility and distribution among multiple nodes of a star topology. We demonstrated the existence of an optimized number of nodes that resulted in the most uniform agent distribution. Regarding the travel efficiency of the star topology, we found that agents conducting congestion avoidance linearizes the travel time for all agents visiting all nodes; however, the travel time increases exponentially with the number of nodes when they do not avoid congestion [6]. A star topology is one of the simplest graphs because it consists of a primary node and multiple secondary nodes. Complex networks are composed of numerous star topologies; the congestion avoidance behavior of agents may significantly impact agent dynamics on complex networks.

In contrast, the static structure of a network often affects the agent dynamics in the network, and such a mechanism may be observed in the human brain. Several studies [7, 8] have reported that the topological structure of brain networks forms interaction and signaling patterns among neurons and brain regions. The resulting communication mechanism is critical to the brain; however, their detailed mechanism remains unclear. Studies on synchronization and consensus problems in multi-agent communication have also been reported [9, 10]. Other scientists have reported the effects of network structures on information spreading [11, 12]. In addition, other related studies have utilized random walk models [13–18] to simulate the dynamics of crowds and neurotransmitters, such as dopamine, or to search for the topologies of complex networks. Clarifying the impact of destination selection strategies on agent dynamics in complex networks and the effect of static network structure on agent dynamics in complex networks is significant for various applications.

This study investigates the agent dynamics on three complex networks: Barabási–Albert (BA) [19–21], Erdos–Rényi (ER) [22], and Watts–Strogatz (WS) [23, 24] models under the following rules: agents (a) randomly select a destination among adjacent nodes, (b) exclude the most congested node among neighboring nodes as a destination and randomly select a destination among the remaining nodes, or (c) select the sparsest adjacent node as a destination. The BA, ER, and WS models are mathematical models designed in the network science field to reproduce the behavior and properties of real-world complex networks. First, the BA model is a network model in which nodes are preferentially added to the existing nodes with larger connections, leading to the formation of hubs and a power-law degree distribution. Examples of the BA models vary from the World Wide Web, air transport networks, and biological networks [25]. Second, the ER model is a random model that generates networks with a Poisson degree distribution, in which nodes are randomly connected with a fixed probability. The resulting network has fewer clusters and a homogeneous degree distribution. Finally, the WS model has a high clustering coefficient and short path lengths, reproducing the small-world nature of social networks [26]. A set of two nodes adjacent to a node is likely to be connected to each other ('small-world').

As evaluation criteria for agent dynamics in complex networks, we introduce the following two characteristics: the hunting rate, defined as the rate of change of the agent amounts in each node per unit time step, and the imbalance of the agent distribution among nodes. As explained previously, we measured these characteristics using Rules (a)–(c). These rules primarily assume decision-making situations that we frequently observe in the real world when agents move from current locations to destinations. Compared with the full random walk scenario in Rule (a), avoiding the most congested site in Rule (b) is a typical human behavior observed in congested crowds. Congestion-avoidance behavior itself can also be observed in social insects such as ants and bees [27, 28]. Selecting the sparsest adjacent site as a destination would be an evacuation behavior during an emergency. Hence, the primary application of this study is in collective dynamics. However, Rules (a)–(c) may also be applied to other control systems. In air-ground traffic systems, the air controller often instructs aircraft to avoid congested spots in taxiing [29–31]. In this way, the target scenarios under Rules (a)–(c) have great applicability to a wide range of network problems. However, our primary interest is in transportation networks and pedestrian networks, such as visitors to event venues, where the number of route branches would be at most about 20. Therefore, we focus on small complex networks with node degrees ranging from 0 to 20 to study agent behavior in transportation and pedestrian networks As a method for this study, we developed a network of cellular automata (CA) [32, 33], where each cell represents a node. Specifically, each agent on a node repeatedly moves to its adjacent nodes according to Rules (a)–(c), which is a straightforward but promising approach. Interestingly, in multivariate statistics, congestion-avoidance behavior is described as a mutual correlation between agents via congestion information [6], where agents interact with each other. Specifically, agents are prompted by the provided congestion information to congregate at less congested nodes while avoiding congested nodes. From an agent-agent interaction perspective, we see this situation as agents being attracted to other agents at less congested nodes while being repelled by agents at congested nodes; this is a kind of interaction among agents via the provided congestion information. Hence, the target system was a multi-agent system. However, computing the multiagent reinforcement learning (MARL) systems gets computationally extensive as the number of agents increases [34, 35]. In addition, it is still under discussion how to reproduce the functionality of congestion avoidance in a probability transition function. Therefore, we performed CA simulations to investigate complex networks. Section 2 details the proposed methodology.

The remainder of this study is organized as follows. Section 2 describes the design variables of the target problem. In addition, we represent two characteristic values, hunting and imbalance rates, that is, the criteria for agent mobility in complex networks. Next, we describe the numerical conditions of the benchmark tests. Section 3 reports the results of the agent simulations in the aforementioned network models: the BA, ER, and WS models for Rules (a)–(c). Section 4 highlights the four significant findings of our results. (1) The degree distribution, that is, the topological structure of a network, precisely determines the agent distribution among nodes in the steady state when agents perform a fully random walk. We demonstrated that the theory of a grand canonical ensemble in statistical physics describes this fact well. In addition, the agents' decisions to select their destination can alter the overall agent distribution. (2) Agents deciding to always select the sparsest adjacent node makes both the hunting and imbalance rates on a complex network significantly high compared with random-walk scenarios, regardless of the network type, particularly when the network has a high degree and activity rate. However, this is not necessarily true when the network has a significantly low degree. (3) The agents' decision to exclude the most congested node among neighboring nodes as a destination increases the hunting rate while decreasing the imbalance rate when the agent activity is low. Conversely, when agent activity is high, their decision increases both the hunting and imbalance rates, which is true for almost all cases of the ER and WS models owing to the relatively uniform degree distribution but not necessarily true for the BA model because of the localized distribution. (4) The hunting and imbalance rates show periodic undulations with respect to time. We discuss the application of these findings in real-world physical and biological networks. Finally, section 5 summarizes this study.

A network is represented as a graph consisting of nodes connected by edges. In this study, we focused on undirected graphs. We introduce three measures for complex networks. First, the degree represents the number of edges per node. The average degree (k) refers to the average degree value for all nodes. The path length indicates the number of minimum edges traced from one node to another. The average path length (l) represents the average path length of all nodes. In addition, we refer to a group of three nodes that are mutually connected to form a triangle, that is, the closed path with a length of two, as a cluster. We defined the clustering coefficient (C) as the ratio of the number of clusters to the total number of paths with a length of two. k, l, and C are the key characteristics representing the spatial structure of complex networks.

In this study, the target system is described as a set $(N_n, N_a, N_t, A_r, G(X))$, where N_n represents the number of nodes, N_a indicates the number of agents, and N_t represents the number of timesteps in each simulation. A_r indicates the probability of each agent hopping from one node to another connected node, and A_r represents the activity rate of the agents. G is a type of graph: a BA, ER, or WS model. The X of G symbolically suggests a set of design variables for each graph. The BA and ER models have a single design variable in addition to variables $(N_n, N_a, N_t)$. Specifically, the BA model has a design variable N_J that indicates the number of edges attached from a newly added node to the existing nodes in the creation process. The preference of an existing node to be connected to a newly-added node is given by their degrees; that is, the probability for the ith node to be selected as one side of N_J new edges is provided by the ratio of the ith node's degree k_i to the sum of the degrees of all the nodes: $k_{i}/\sum k_{i}$, which suggests that the node with higher degrees can be preferentially selected. This is the reason why the BA model has a so-called 'hub,' a gathering of densely connected nodes. The resulting degree distribution is known to have the scale of k⁻³, a power-law distribution; therefore, the BA network is called 'scale-free.' The average path length l and clustering coefficient C have scales of $\textrm{ln}(N_{n})/\textrm{ln ln} (N_{n})$ and $N_{n}^{-0.75}$, respectively, in the BA model [36]. In contrast, the ER model generates an edge between each pair of nodes with probability P_r . That is, each of the $N_{n}(N_{n}-1)/2$ combinations of N_n nodes creates an edge with probability P_r . When $P_{r} = 1$, each node has $N_{n}-1$ connections; conversely, if P_r is sufficiently small, isolated nodes with no connection to other nodes emerge. Because each pair of nodes has an edge with probability P_r and does not have it with $1-P_{r}$, the probability of the existing M edges among N_n nodes obeys the binominal distribution. Additionally, the average path length l and clustering coefficient C of the ER model have scales of $\textrm{ln}(N_{n})/\textrm{ln} (k)$ and $k/N_{n}$, respectively, [36]. The generation of complex networks with isolated nodes should be avoided because this study focuses on the agent dynamics in complex networks. We examined the connectivity of the generated static networks for different values of P_r between zero and one before performing agent simulations. We determined the suitable parameter range of P_r to be between 0.4 and 1.0 for the simulations described in section 3.

Unlike the BA and ER models, the WS model has two design variables: the number of nodes connected to each node in the initial state N_J and the probability of rewiring edges P_r . During the creation process, all nodes are arranged in a ring at regular intervals in the initial state. Each node has an edge to its neighboring nodes. Here, the range of neighbors depends on the size of N_J ; when $N_J = 2$, each node has edges only to the closest nodes in addition to the adjacent nodes on the ring. In contrast, it has four edges to the closest and second-closest nodes when $N_J = 4$. As the value of N_J increases, neighboring nodes become mutually connected. Thus, the clustering coefficient C increases with N_J . By rewiring the edges among the nodes with probability P_r , the clusters collapse; as P_r increases, C decreases. The average path length l decreases rapidly owing to the rewiring process, which connects a node with the non-neighboring nodes. This rapidly decreases the average path length, and eventually, the network approaches a random network structure. When a network has a small l while maintaining a high clustering coefficient C, the network is called a 'small-world' network because the neighboring nodes of a node are likely to be connected. That is, a person's friends are likely to be friends. As for the remaining key characteristics, the average path length l and clustering coefficient C of the WS model have scales of $\frac{N_{n}}{k}\cdot f(P_{r} k N_{n})$ and $\frac{3(k-2)}{4(k-1)}(1-P_r)^{3}$, respectively [36], where f(x) represents a function that becomes constant when $x \ll 1$ and becomes $\textrm{ln}(x)/x$ when $x \gg 1$. In summary, the independent variables of the BA model are $(N_n, N_a, N_t, A_r, N_J)$, those of the ER model are $(N_n, N_a, N_t, A_r, P_r)$, and those of the WS model are $(N_n, N_a, N_t, A_r, N_J, P_r)$. The physical significances of N_J and P_r differ according to the model. As described previously, for the ER model, P_r is confined between 0.4 and 1.0; otherwise, 0 and 1. We investigated target systems with different values of N_J between 1 and 10, assuming the number of junctions or branches of real-world traffic networks on which the agents travel. Activity rate A_r was set from 0.1 to 0.9 at the initial state. The method for determining the remaining parameters $(N_{n}, N_{a}, N_{t})$ is described below.

We used NetworkX [37], a most well-established Python package for complex networks. Specifically, we implemented the BA network using barabasi_albert_graph, the ER network using erdos_renyi_graph, and the WS network using watts_strogatz _graph. Parameter N_J of the BA model corresponds to the input parameter 'm' of barabasi_albert_graph. Because they determine the number of edges attached from a newly added node to the existing nodes, the BA model has a high degree as N_J increases. Similarly, parameter P_r of the ER model corresponds to the input parameter 'p' of erdos_renyi_graph that defines the probability of creating an edge between each pair of nodes; as P_r increases, the ER model has a high degree. Additionally, parameters N_J and P_r correspond to the 'k' and 'p' of watts_strogatz_graph, where k represents the number of nodes connected to each node at the initial state, and 'p' indicates the probability of rewiring edges. Thus, as N_J increased, the WS model exhibits a higher degree. Figure 1 shows the dependence of the average degree k on the design variables N_J for BA, P_r for ER, and N_J for WS. We can confirm that N_J increases for the BA or WS models and P_r increases for the ER model, and that each model has a higher degree. Reasonably, the WS model shows a step function because it connects $N_{J}-1$ neighbors if N_J is odd [23, 37]. In addition, we confirm that the average degree k is independent of P_r in the WS model. To be clear, the fact that the average degree k increases as the parameter N_J or P_r increases does not mean that the clustering coefficient C also increases as N_J or P_r increases. In fact, related studies suggest that the magnitude of the clustering coefficients asymptotically converges to the average clustering coefficient in certain cases [38]. In any case, we have confirmed the average degree k increases as N_J or P_r increases within the focused parameter range.

Zoom In Zoom Out Reset image size

As mentioned in the Introduction, we investigate the agent dynamics on three complex networks under the following rules: agents (a) randomly select a destination among adjacent nodes; (b) exclude the most congested node among neighboring nodes as a destination and randomly select a destination among the remaining nodes, or (c) select the sparsest adjacent node as a destination. The mathematical expressions for Rules (a)–(c) on the ith node are as follows:

$$\begin{eqnarray} D^{\left(i\right)}_{s}\left(X\right),~~X~: = ~ \begin{cases} ~\left\{x \in N_{d}^{\left(i\right)}~|~x = ~{\mathrm{rand}}~\left(n_{a}^{1},n_{a}^{2},\cdots,n_{a}^{k_{i}}\right)~\right\} & {\mathrm{if~~Rule~(a)}},\\[4pt] ~\left\{x \in N_{d}^{\left(i\right)}~|~x = ~\overline{{\mathrm{max}}}~\left(n_{a}^{1},n_{a}^{2},\cdots,n_{a}^{k_{i}}\right)~\right\}& {\mathrm{if~~Rule~(b)}},\\[4pt] ~\left\{x \in N_{d}^{\left(i\right)}~|~x = ~{\mathrm{min}}~\left(n_{a}^{1},n_{a}^{2},\cdots,n_{a}^{k_{i}}\right)~\right\} & {\mathrm{if~~Rule~(c)}},\\ &&\end{cases} \end{eqnarray} \tag{ 1 }$$

where $N^{(i)}_{d}$ represents the set of all nodes indexed from the viewpoint of the ith node, and k_i represents the ith node degree, that is, the number of nodes connected to the ith node. $n^{j}_{a}$ indicates the number of agents in the jth neighboring node. The function $\mathrm{\mathsf rand(\cdots})$ returns one of $\{n^{1}_a, n^{2}_a, \ldots, n^{k_i}_{a}\}$ at random; X has only one selected element in this case. The function $\mathrm{\mathsf{\overline{max}}}(\cdots)$ returns the complement of the function $\mathrm{\mathsf max}(\cdots)$ that returns the maximum value among $\{n^{1}_a, n^{2}_a, \ldots, n^{k_i}_{a}\}$; in this case, X usually has multiple elements. The function $\mathrm{\mathsf min}(\cdots)$ returns the minimum value among $\{n^{1}_a, n^{2}_a, \ldots, n^{k_i}_{a}\}$; X has multiple elements when two or more of $\{n^{1}_a, n^{2}_a, \ldots, n^{k_i}_{a}\}$ have the same minimum value. $D^{(i)}_{s}(\cdots)$ is a function that converts $n^{l}_a$ to parameter l, a local index of the neighboring node that will be selected as a destination, after randomly selecting one element in X when X has multiple elements.

We introduce the following two physical quantities as indicators of agent dynamics in a complex network: the hunting rate H, which is the change rate of agent counts in each node per unit time step, and imbalance rate L, which is the non-uniformity of the agent distribution among nodes. Using parameters N_n and N_a , the hunting rate H and imbalance rate L are expressed in the following forms:

$$\begin{eqnarray} H &: = \frac{1}{N_{a}} \sum\limits_{i = 1}^{N_{n}} \biggl | {N_{n}}_{\left(m\right)}^{\left(i\right)} - {N_{n}}_{\left(m-1\right)}^{\left(i\right)} \biggr |, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} L &: = \frac{1}{N_{n}} \sum\limits_{i = 1}^{N_{n}} \biggl | {N_{n}}_{\left(m\right)}^{\left(i\right)} - {N_{n}}^{avr} \biggr |, \end{eqnarray} \tag{ 3 }$$

where $N_{n}^{\mathrm {avr}} = N_{n}/N_{a}$ is the average number of agents per node. ${N_{n}}_{(m)}^{(i)}$ represents the number of agents on the ith node at the mth time step. Equation (2) shows that parameter H increases as the change rate of the number of agents per unit of time increases, suggesting that H is an indicator of the risk that agents thrust into a node in a short time. We can say that agent mobility in a complex network is moderate when the value of H is low in the stationary state. Conversely, when the system exhibits high values of H, the agents become rapidly mobile. We introduced $1/N_{a}$ into equation (2) for normalization based on the problem size. For equation (3), parameter L increases as the deviation of the number of agents from the average number of agents per node increases, which indicates that the agent distribution among nodes becomes uniform as L decreases, and the agent distribution becomes non-uniform as L increases. In addition, we adopted the parallel update scheme for the time integration in which the states of all nodes are simultaneously updated at every time step.

Originally, 'hunting' was used as an electrical engineering term referring to the phenomenon of unexpected oscillations in the time direction of a characteristic value in an electrical system, such as a power system [39], of which a motor voltage disturbance is a typical example. The term is now being used in other fields and is becoming popular to refer to a phenomenon in which a certain characteristic of a complex system causes unexpected oscillations in the time direction. We refer to this as the intensity of variation of the agents at each node per unit of time in complex networks. In a real pedestrian network, the greater the number of agents moving along a given path in the network in the amount per unit of time, the more they can collide with each other, leading to an increase in accidents such as crowd avalanches in the most congested case. In this sense, parameter H in equation (2) can be said to be an objective metric that corresponds to the risk measure for crowd accidents in such real networks. On the other hand, the imbalance rate L in equation (3) is a measure of how much an agent distribution per unit time deviates from the case where agents are evenly distributed at each node. We are interested in how much L changes when the rule is changed rather than the value itself. For example, changing Rule (a) to Rule (b) causes agents to avoid the most congested nodes. Each agent is a free agent under Rule (a), but under Rule (b), it becomes a decision entity that wants to avoid the most congested node. If we approximately consider the behavior of agents visiting an event site consisting of multiple booths as random walks, the variation of L when we change Rule (a) to Rule (b) may serve as an indicator of how people's collective behavior changes when they are provided with congestion information and change from randomly wandering free agents to decision-making entities. Meanwhile, Rule (c), where agents always choose the most sparse node, reflects the tendency of pedestrians to avoid congestion during evacuations. We also know from our experience that we often avoid congested places when visiting amusement parks, shopping malls, or attending multiple events. Compared to Rule (b), the will to avoid congestion in Rule (c) is stronger because the agent in Rule (b) only avoids the most congested node and has the right to choose one of the other neighbors, whereas the agent in Rule (c) always makes a coercive and exclusive decision to choose the most vacant node. By comparing Rule (b) and Rule (c), we can see how the different degree of congestion avoidance decision manifests itself in the collective dynamics of the agents. In this way, Rules (a)–(c) are all associated with real-world scenarios.

Many network properties discussed in previous studies are mainly related to spatial structures such as clustering coefficients and fractal dimension [40–42]. For multi-agent dynamics where agents perform mutual interactions such as attraction or repulsion, there is a study that presents statistical physical metrics, including self-organization and complexity [43]. In contrast, our metrics in equations (2) and (3) are oriented toward estimating agent mobility and distribution in a social system where people with self-intelligence of select destinations, such as visitors to an event site.

First, we measured the values of H and L averaged over the elapsed time step for several cases of $(N_n, N_a)$ as a preliminary test. This is to identify the effect of interactions between agents due to congestion avoidance behavior on the long-term evolution of these two metrics. Consequently, we found that H and L reached a steady state shortly after the simulation started in all cases of $(N_n, N_a)$, except that the time step required for the system to reach a steady state minorly differs in each case. Figure 4(A) shows the dependencies of H on the time steps for $(N_n, N_a) = (30, 2000)$ for the BA model with different values of N_J , as shown in figure 5(A) for the ER model for different values of P_r , and figures 6(A), 7(A) and 8(A) for the WS models with different values of N_J and probabilities $P_{r} = 0.1$, $P_{r} = 0.3$, and $P_{r} = 0.9$, respectively. Each column (A) of figures 4–8 displays the cases of Rules (a)–(c) from the left. The same observations were made for L, as displayed in each column (C) of the corresponding figures for each model. These results confirm that the value of N_t required to reach a steady state can be specified as $(N_n, N_a)$, which is approximately 600 for $(N_n, N_a) = (30, 2000)$. We also confirmed that the coefficients of variation for the time series data for the variation of the averaged values of H and L in a steady state area shown in (A) and (C) of figures 4–8 were less than 0.5 percent. Accordingly, we start our analysis with the set $(N_{n}, N_{a}, N_{t}) = (30, 2000, 600)$ and then discuss the extendibility of the obtained results to other cases of $(N_{n}, N_{a})$.

We equally distributed the agents among nodes at the initial state. We then conducted agent simulations wherein agents walked around among nodes according to Rules (a)–(c) under the condition of $(N_n, N_a, N_t) = (30, 2000, 600)$ for the BA, ER, and WS models for different values of N_J and P_r . Consequently, the agent and degree distributions were observed to agree exactly irrespective of the network type and design variables N_J and P_r . Figure 2 shows example results of the agent simulations when setting N_J to 0.5, P_r to 0.5, and $(N_J, P_r) = (0.5,0.5)$ for the WS models. In figure 2(F), the top row illustrates the generated network graphs, and the bottom row presents the corresponding degree distributions. In addition, figure 2(A) shows the agent distribution in the steady state for Rule (a). Compared with (A) and (F), the agent and degree distributions were confirmed to match in all three network model cases. The blue bars in figure 2(E) show the values of R², which is the coefficient of determination between the agent and degree distributions for Rule (a); the R² value becomes one with significant accuracy in all three models.

Zoom In Zoom Out Reset image size

Let us elaborate on the agreement between the agent and degree distributions using statistical mechanics. We assume that the agent distribution reaches an equilibrium state, implying that the flows of agents per node are equal. Next, we consider the ith node, which has a degree of k_i , that is, k_i connected nodes. Upon reaching the equilibrium state, we can consider the ith node a thermodynamic system connected to multiple reservoirs with agents as particles flowing into/out of the system; the system is a grand canonical ensemble [44]. Thus, we consider the system an absorption problem. Specifically, we regard the ith node with k_i neighboring nodes as a node with k_i adsorption sites. When a single agent (particle) flows into the ith node, one of the k_i adsorption sites is filled with the particle. Assuming that the energy state of a node is proportional to the number of agents when the system has only agents, the energy state of the system can be represented as E₀. This assumption is valid because the state of a node is designated only by the number of agents in the target system. However, the total number of ways in which k_m sites among the k_i sites are filled with agents corresponds to the number of combinations of selecting k_m sites from k_i nodes, expressed as $C_{i} = {k_{i}!}/{{k_{m}!} {(k_{i}-k_{m})!} }$. Because each state has the same energy E₀ multiplied by k_i as $E_{0} k_{i}$, the resulting grand partition function is expressed as $\Xi = [1 + \textrm{e}^{\beta(E_{0}+\mu)}]^{k_{i}}$, where β and µ denote the inverse temperature and chemical potential of the node. Using a series of formulas for the grand partition function, the expected number of sites filled with agents, that is, the expected number of agents in the ith node, is obtained by the partial derivative of the natural logarithm of Ξ scaled by β⁻¹ with respect to the chemical potential µ as follows:

$$\begin{eqnarray} &&\langle N_i\rangle~: = ~\frac{1}{\beta} \frac{\partial}{\partial \mu} {\mathrm{ln}} \Xi~ = ~k_i \frac{{\mathrm{e}}^{\beta\left(E_0+\mu\right)}}{1 + {\mathrm{e}}^{\beta\left(E_0+\mu\right)}}. \end{eqnarray} \tag{ 4 }$$

Here, all nodes have the same parameters β and µ because the activity rates A_r of all nodes are the same. Hence, the expected distribution of the number of agents is expressed as follows:

$$\begin{eqnarray} &&\left\{N_{1}, N_{2}, \ldots, N_{N_n}\right\}~ = ~D \left\{k_{1}, k_{2}, \ldots, k_{N_n}\right\},~~~~D~: = ~\frac{{\mathrm{e}}^{\beta\left(E_0+\mu\right)}}{1 + {\mathrm{e}}^{\beta\left(E_0+\mu\right)}}~ = ~{\mathrm{const}}. \end{eqnarray} \tag{ 5 }$$

Therefore, the agent distribution corresponds to the degree distribution scaled by a constant coefficient D; they agree with each other, and the R² value becomes one after normalization. As described later, characteristic quantities, such as the hunting and imbalance rates, show slight periodic undulations in the time direction; however, these undulations are typically less than one percent and are negligible by assuming an equilibrium state.

The advantage of proving the agreement between the agent and degree distributions using the grand canonical ensemble in statistical physics is that equation (5) can be applied in infinite networks because our description requires only the assumption of equilibrium states for each node and an identical activity rate for all nodes. In contrast, the previous explanation assumes a finite network with a stationary density distribution to define the transition probability matrix [17]. Infinite networks differ significantly from finite networks in the random walk behaviors of the agents, such as that agent recursion is not guaranteed. Our description shows that in real-world systems, even networks that can be considered infinite networks, such as neural networks, supercrystals [45–47], and biomolecular networks [48], are guaranteed to match these two distributions as long as the above conditions are satisfied. Furthermore, because the activity rates among the nodes are the same, this model is considered suitable for homogeneous networks. This is the first study to prove that these two distributions coincide in infinite networks using the grand canonical ensemble of statistical physics. Accordingly, we demonstrated that the degree distribution precisely determines the agent distribution when the agents perform a fully random walk.

Two significant findings were obtained. First, the accuracy of the agreement between the agent and degree distributions depends on the ratio of the number of agents to the number of nodes $N_{a}/N_{n}$. Figure 2(D) shows the dependence of R² on $N_{a}/N_{n}$ for the two cases, $N_n = 30$ and $N_n = 300$. R² values are confirmed to converge to one when $N_{a}/N_{n}$ is approximately greater than 60 in both cases. That is, the agent and degree distributions become congruent although the number of nodes increases while the number of agents is sufficient, such as $N_{a}/N_{n} \gt 60$. Second, the agents' congestion-avoiding behavior alters the agent distribution from one with a fully random walk in the steady state. The red bars in Rule (b) of figure 2(E) show that the R² value decreases compared with Rule (a). Although Rule (b) differs only from Rule (a) in that each agent excludes the most congested adjacent node as a destination before it starts walking randomly, this difference alters the overall agent distribution. A similar observation is true and more distinctive in Rule (c), where agents always select the sparsest adjacent node. The R² value then further decreases compared with Rule (b), as shown by the white bars in figure 2(E).

Next, we measured the hunting rate H and imbalance rate L averaged for the elapsed time step by varying the N_J values between 1 and 9 for the BA model, P_r values ranging from 0.4 to 1.0 for the ER model, and N_J values between 2 and 9, with P_r set to 0.1, 0.5, or 0.9 for the WS model. We collected 200 simulation results for each rule and compared Rules (a)–(c) using 600 different datasets. We report the effect of rule differences on the agent dynamics in complex networks extracted from these datasets. First, our findings show that for a network with a high degree caused by the high values of N_J or P_r , which are approximately 8–9 for N_J and 0.8–1.0 for P_r , Rule (c), with which agents always select the sparsest adjacent node, drastically increases the hunting rate H and imbalance rate L compared with the random-walk scenarios, as in (a) and (b). We can observe such surges in H and L in Rule (c) for all three network types in the high degree ranges, as shown in all subfigures of figures 4–8. The exceptional cases are the imbalance rate L of the BA and WS models at $A_{r} = 0.1$ in figures 4 and 6–8, for which the activity rates were sufficiently low. Second, the aforementioned surges in H and L with Rule (c) were most distinctive with the ER model, as shown in figure 5, for all measured cases. In contrast, they were not necessarily evident in the BA and WS models when the network has a low degree owing to small values of N_J and P_r such as $N_{J} = 2$ and $P_{r} = 0.4$. From a collective dynamics viewpoint, the results can be interpreted as follows: when each individual agent intentionally selects the sparsest adjacent node as their destination, both the hunting and imbalance rates drastically increase regardless of the network type, as long as the network has a high degree and activity rate, compared with the cases in which agents walk randomly. However, this phenomenon is not necessarily observed when the network has a low degree or activity rate.

Moreover, compared with random walk-based Rules (a) and (b), we found that Rule (b), with which agents exclude the most congested adjacent nodes as a destination, increases the hunting rate H and concurrently decreases the imbalance rate L compared with Rule (a) when the activity rate A_r of the agents is small, such as $A_{r} = 0.1$. This trade-off relationship is valid for all network types when $A_r = 0.1$ for almost all N_J and P_r parameter ranges. For example, we compare the cases at $A_r = 0.1$ for (B) and (D) in figure 5. We can confirm that the red bar for Rule (b) is greater than the blue bar for Rule (a) for the hunting rates in (B), whereas the relations are inversed for the imbalance rates of the corresponding cases in (D). Similar observations are made for the other network models when $A_r = 0.1$ and $A_r = 0.3$, as shown in figures 4–8. Conversely, Rule (b) increases both H and L compared with Rule (a) when A_r is sufficiently large, such as $A_{r} = 0.9$. From a collective dynamics perspective, these observations can be explained as follows. In principle, avoiding the most congested nodes effectively mitigates the imbalance in agent distribution as the agents are distributed. However, agents move to less-congested nodes simultaneously. Thus, the congestion avoidance behavior in Rule (b) reduces the imbalance rate and increases the hunting rate. This is why the hunting and imbalance rates exhibit a trade-off relationship when the activity is low. However, when group activity is too high, the hunting rate increases, and the magnitude of the dispersion effect from avoiding congestion remains at the same level. In contrast, most agents at the currently congested node move to another node the next time owing to high activity; however, the agents do not select the currently congested nodes as their destinations because they determine them based on the congestion information at the current time. Consequently, the congested node becomes one of the sparsest adjacent nodes the next time. This mismatch between destinations increases the imbalance rate. Accordingly, congestion avoidance behavior can significantly increase the hunting and imbalance rates. These observations are distinctive for the ER and WS models because they have a relatively uniform degree distribution, which are not necessarily true for the BA model because of its localized distribution. In the BA model, the agent distribution is localized owing to the localized degree distribution. Nodes with higher degrees tend to have a large number of agents, even with Rule (b), as shown in (B) in figure 2. Therefore, nodes with higher degrees can generate more significant destination mismatches, leading to high imbalance rates. Because the agent distribution depends on the degree distribution, the obtained imbalance rate differs depending on the parameter N_J , as shown in the subfigure for $A_{r} = 0.9$ in (D) of figure 4.

In summary, we made the following two observations regarding decision-making strategies. First, the agents' decision to always select the sparsest adjacent node among neighboring nodes makes both the hunting and imbalance rates significantly higher compared with the random-walk scenarios, particularly when the network has a high degree and activity rate; this is observed in all three network models. Second, the agents' decision to exclude the most congested adjacent node as a destination increases the hunting rate while decreasing the imbalance rate when agent activity is low. Conversely, when agent activity is high, their decisions increase both the hunting and imbalance rates. The second observation is distinctive for the ER and WS models because they have relatively uniform degree distribution, but this is not necessarily valid for the BA model because of the localized agent distribution. Figure 3 summarizes these observations for the ER model. Here, parameters H and L are respectively scaled by the larger of the measured values in each comparison. It was confirmed from figures 3–8 that relationship (1) in figure 3 was true for all three network models. Moreover, relationships (2) and (3) in figure 3 were true for the ER and WS models; however, they were not necessarily observed in the BA model.

Zoom In Zoom Out Reset image size

As a further step in our research, we performed a Fourier analysis of the time-dependent variation in the hunting and imbalance rates. In the previous analysis, the average values of equations (2) and (3), over time. However, in real-world systems, we are interested in the temporal variation of these characteristic quantities, in addition to the mean, because short-term fluctuations in these physical quantities can cause state transitions, such as crowd avalanches and agent collisions. We obtained the autocorrelation functions of the time series data in equations (2) and (3) by multiplying the data using the Fourier transform and inverse Fourier transform of the power spectrum. Figure 9 displays all the results of the obtained autocorrelation functions for 200-time series data, where S_h and S_l denote the values of the autocorrections of H and L, respectively. The vertical axes of the respective subfigures in the left columns show the results of S_h , and those of the subfigures in the right columns display those of S_l . The horizontal axis in each subfigure represents the scaled time T_s , which is obtained by dividing the number of time steps by the total number of time steps N_t . Each data is shifted by $N_t/2$ to center the time base point. Each subfigure is enlarged such that the midpoint of the maximum and minimum values is located at the center of the subfigure in the vertical direction, and the maximum and minimum values fit into the figure. The common denominator for all the data is that the autocorrelation function shows periodic undulations over multiple time steps for several time periods; however, these undulations are typically less than one percent. The patterns of the periodic undulation were case-by-case. Some results are expanded and highlighted, and other results are shown as thumbnails for reference. Accordingly, we confirmed that the hunting and imbalance rates exhibit periodic undulations over time.

First, we summarize our findings. Our simulation study supports the following findings on the effect of the topological network structure and the three different effects of agents' destination selection on the collective dynamics of agents in three complex networks: the BA, ER, and WS models.

• (A)  
  The degree distribution, that is, the topological structure of a network, precisely determines the shape of the agent distribution among nodes in the steady state when agents perform a full random walk. We demonstrated that the theory of a grand canonical ensemble in statistical physics describes this well. In addition, their decisions on selecting destinations alter the overall agent distribution.
• (B)  
  The agents' decision to always select the sparsest adjacent node among neighboring nodes makes both the hunting and imbalance rates significantly high compared with the random-walk scenarios, particularly when the network has a high degree and activity rate; this is observed in all three network models.
• (C)  
  The agents' decision to exclude the most congested adjacent node as a destination increases the hunting rate while decreasing the imbalance rate when agent activity is low. Conversely, when agent activity is high, their decisions increase both hunting and imbalance rates. This is observed for the ER and WS models because they have relatively uniform degree distribution, and it is not necessarily valid for the BA model because of the localized agent distribution.
• (D)  
  The hunting and imbalance rates show slight periodic undulations with respect to time.

Regarding (A), many studies would have empirically acknowledged in their respective fields that the degree distribution of a network can contribute to the formation of an agent distribution, even if the network is an infinite network such as a neural, supercrystal, or biomolecular network. However, mathematically proving the agreement between agent and degree distributions is difficult because we cannot define a stationary probability density distribution for infinite networks owing to its openness; it follows that a transition probability matrix to the network cannot be defined. This study is the first to quantitatively ensure that the degree distribution determines the shape of the agent distribution for general cases, including infinite networks, using a theoretical model of a grand canonical ensemble in statistical physics. In addition, equation (5) assumes a grand canonical ensemble in each node and the condition that the activity rate is the same for all nodes, and the latter allows us to derive equation (5) from equation (4). Here, we can rephrase these conditions as each node reaching an equilibrium state where the amount of inflow/outflow is at least quasi-statically equal in each node and its magnitude is equal everywhere. They can be satisfied as long as there are fully random walks in complex networks that have no isolated nodes. A good example is the PageRank [49, 50] distribution over web pages, which is known to be a fully random walk problem. However, this proof does not hold if there are differences in the scale of the activity rate in the respective nodes or if the independence of the agents is broken by behaviors such as congestion avoidance. Nevertheless, the statistical mechanical explanation of the correspondence between agent and degree distributions discussed here can be applied to real cases, such as epidemic metapopulation models, which in certain cases assume a fully random walk in a network connecting subpopulations [51].

The application of findings (A)–(C) to real-world collective dynamics can be described as follows. First, (A) allows us to estimate the agent distribution from the degree distribution, which is the topological structure of the network, in situation in which agents can be regarded as performing fully random walks, such as in crowds, metapopulations, and the dynamics of neurotransmitters in the brain. Once we identify the hubs of the nodes from the degree distribution, we can predict the congested areas of crowds because of the agreement between the agent and degree distribution, which may help prevent crowd incidents. (A) also helps us construct facilities where people congregate. By referring to (A), we can design these facilities to avoid the formation of dense crowds by eliminating the hubs in the network. (A) is also applicable in systems that guide vehicles, that is, for agents in highway parking and air and ground transportation systems with multiple lanes for arriving and departing aircraft. If the degree of the distribution of parking and ground traffic conduits is known, we can predict the locations where congestion may occur.

(B) and (C) significantly contribute to accident prevention. Our results suggest that agents avoiding the most congested adjacent nodes concurrently increases the hunting rate and decreases the imbalance rate when the crowd is inactive. Therefore, if crowd managers need to prioritize dispersing people in situations where people are not very active or stuck, instructing agents to avoid the most congested sites is effective in distributing people while admitting an increase in the movement of agents per unit time. However, we revealed that intentionally conducting avoidance behaviors increases both the hunting and imbalance rates when the crowding activities are high; this instruction should be averted in such cases. In addition, we emphasize that the agents' decision to select the sparsest adjacent node should be avoided, except for in some sufficiently low-degree networks, based on (B). This selection makes both the hunting and imbalance rates significantly high on a complex network, regardless of the network type. Findings (A)–(C) can be applied to problems involving controlled agents in complex networks. For several ground traffic controls at airports, the airport's geometric structure may create areas where taxis are likely to crowd. According to (A), airport controllers can consciously instruct pilots in advance to avoid concentrating on taxing aircraft. At this time, directing aircraft to the sparsest adjacent spots may significantly further imbalance the airport spot occupancy and number of aircraft movements per unit time owing to (B). Furthermore, we can understand from (C) that weighing the advantages of directing aircraft to move to other locations would be beneficial in reducing the imbalance rate against the disadvantages of fuel consumption and additional costs involved in moving aircraft in a short period.

Nevertheless, we must discuss the difference between the real cases and the cases with Rules (a)–(c) when applying our results to real cases. Recall an important result. In Rule (b), avoiding the most congested nodes effectively mitigates the imbalance in agent distribution as agents are distributed among remaining nodes. However, the number of agents in less congested nodes simultaneously increases. Therefore, the congestion avoidance behavior in Rule (b) reduces the imbalance rate and increases the hunting rate when the activity rate is low; in other words, we can see a trade-off relationship between Rules (a) and (b), as shown in figure 3, when the activity is low. However, as the group activity becomes high, the hunting rate increases and the magnitude of the dispersal effect of avoiding congestion remains at the same degree of avoiding the most congested node. Consequently, both the hunting rate and the imbalance rate increase under Rule (b) compared to Rule (a) when the activity rate is high. In real cases, however, agents may still choose the most congested node as their destination if the node is still attractive, even considering its congestion. This is exemplified by the situation where agents can get some incentive (e.g. food or services limited at the node). Thus, if practical conditions exist such that the node has a priority, an external influence that forces agents to the node, or a high popularity among agents, the results will differ from the case of Rule (b). We should consider these points when extending Rule (b) to real cases, where the complexity of decision-making holds.

Recall another important finding about agent dynamics with Rule (b) when the activity rate is high. In this case, agents at the currently most congested node move to another node at the next time step. However, agents at the other nodes do not choose the currently most congested node as their destination to avoid congestion; this decision-making is based on the congestion information at the current time. As a result, the most congested node becomes one of the most sparse neighboring nodes the next time. In short, destination mismatches increase the congestion rate. Accordingly, Rule (b) significantly increases hunting and imbalance rates compared to Rule (a). These observations were particularly evident for the ER and WS models, which have relatively uniform degree distributions. Importantly, the aforementioned destination mismatches occur because agents simultaneously decide their destinations based on the current congestion information. However, in real cases, agents can make more comprehensive decisions based on current and past time congestion information, namely, predicting whether the congestion at the most congested node is transient or persistent based on experience. In short, agents can choose the currently most congested node as their destination when they know that the node will become emptier at the next time. This strategic decision with a higher level of thinking can prevent the simultaneous migration of large agents, which mitigates the hunting rate and the imbalance rate. Therefore, we must consider this point when extending the obtained fact that destination mismatches increase the hunting and imbalance rates to real cases.

Finally, (D) may contribute to the field of neuroscience more than engineering. The quasi-periodic patterns of neural activity in the brain contribute to typical brain functions [52, 53]. In addition, this function has been reported to be reduced in patients with ADHD [54]. Here, figure 9 reported that the phase of the periodic undulation of network characteristics differs depending on the network type and moving rules of the agents. If we consider the brain as a complex network and the agents as signals produced by neuronal interactions, the periodic undulation of network properties reported in (D) may represent the simplest form of periodic patterns of neural activity. Therefore, the correspondence between our results and those of real-world systems must be further investigated in detail. Nonetheless, in this manner, we can say that our findings are helpful for the agent mobility problem in complex networks. As shown in (A) and (C) in figures 4–8, the averaged imbalance rate and the hunting rate reach a steady state within the time and space scale. However, there was a possibility that we missed the slight periodic undulations in the original time series, which is negligible after the averaging operation. In terms of risk management in the crowd, it is more important to detect the variation of H and L in local time and space than the variation of their averages. This is because the abrupt increase in H in a small period of time indicates the occurrence of rapid mass migration, which could cause accidents such as collisions among agents in real cases. If we find the existence of such a slight periodic undulation and the correlation between its periodic behavior and the abrupt variation of H or L in a short time, we can use it for crowd control by predicting the occurrence of accidents in advance. However, to the best of our knowledge, no study has investigated the periodic undulations of the feature value of multiple agents in complex networks using CA simulations. Therefore, as a preliminary step for future studies, we focused on revealing the existence of such periodic undulations in the respective rules. As shown in figure 9, this study confirmed the existence of slight periodic undulations to a certain extent. Recall that we focused on the Markovian process (i.e. the state at n + 1 time step determined only from n time steps) for all analyses in this research. In this respect, a similar discussion on non-Markovian agent networks [55, 56] may allow us to discuss the system closer to real cases; this should be done in future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a prelude to future studies, we discuss the resilience of the system under our scenarios. We performed agent simulations in which agents follow Rules (a)–(c) on BA, ER, and WS network models with a functionality of disconnecting nodes in a certain time period during the simulation; the results are shown in figure 10. Note that we set the parameters $(N_n, N_a, N_t)$ to $(30, 2000, 600)$, similar to the other cases. In the simulations, we randomly selected nodes with a probability of 50$\%$ at the time of 300 time steps, then isolated the selected nodes by forbidding agents to travel to them during the time from 300 to 400 time steps. We then reconnected these isolated nodes to networks by allowing agents to access them when the simulation reached 400 time steps. We measured the time variation of the mean value of the hunting and imbalance rates in all measurable cases for BA and ER, and P_r = 0.5 for WS, similar to figures 4, 5 and 7; therefore, the measured cases of H and L for each network model correspond to (A) and (C) in figures 4, 5 and 7, respectively.

Zoom In Zoom Out Reset image size

The following results were obtained. First, it was found that the hunting rate H tends to decrease and the imbalance rate L tends to increase as soon as the isolation of a node starts, regardless of the network type. When the network is restored, both the hunting rate H and the imbalance rate L are found to recover gradually (in some cases not completely within the measurement time). The reason why H decreases and L increases with node isolation can be explained as follows. When a node is isolated, there are fewer candidates as destinations for agents. The imbalance rate represents the average deviation between the number of agents per node when agents are evenly distributed and the actual number of agents in each node. However, as the number of agents at the isolated node is forced to always be zero, these deviations increase, and the imbalance rate increases compared to before the isolation operation. Meanwhile, the hunting rate represents the average change from the previous time in the number of agents for each node. Since node isolation reduces the number of candidate nodes as destinations for agents to visit, agents are forced to travel to remaining nodes other than the isolated nodes. Because the same number of agents is spread over fewer nodes, there will always be a certain number of agents per node, and the frequency of situations where there are no agents at all at a node will decrease. In addition, since the number of agents in the isolated nodes is zero, their hunting rate is always zero. Thus, the overall hunting rate is reduced by isolation. As mentioned above, the hunting rate is a measure of the average amount of migration per unit of time, so it is desirable for the hunting rate to be small in order to reduce the risk of pedestrian collisions and other occurrences in real systems. In a particular transport network, if the hunting rate is large, it may be possible to reduce it by daring to shut down some nodes. Of course, this would need to be verified in more realistic scenarios and networks.

We examined the agent dynamics in three complex networks: the BA, ER, and WS models under the following rules: agents (a) randomly select a destination among adjacent nodes, (b) exclude the most congested adjacent node as a potential destination and randomly select a destination among the remaining nodes, or (c) select the sparsest adjacent node as a destination. We measured two significant characteristics in the wide parameter ranges of the network characteristic parameters: the hunting rate, which is the change rate of agent amounts in each node per unit time step, and imbalance rate, which is the non-uniformity of agent distribution among nodes.

Our analysis suggests that the degree distribution, which represents the topological structure of a network, precisely determines the shape of the agent distribution among nodes in steady state when agents perform a full random walk. The model of a grand canonical ensemble in statistical physics was determined to describe this well. However, the decision to select a destination altered the overall agent distribution. Particularly, the agents' decision to always select the sparsest adjacent node among neighboring nodes makes both the hunting and imbalance rates significantly higher compared with random-walk scenarios, particularly when the network has a high degree and activity rate, which was observed in all three network models. In contrast, the agents' decision to exclude the most congested adjacent node as a destination increases the hunting rate while decreasing the imbalance rate when agent activity is low. Conversely, when agent activity is high, their decisions increase both the hunting and imbalance rates. This was observed for the ER and WS models because they had a relatively uniform degree distribution, and was not necessarily valid for the BA model. In addition, the hunting and imbalance rates show periodic undulations with respect to time.

In this paper, we have studied the periodic undulation of two characteristic values of agent behavior, as summarized in figure 9. Although we have focused on three representative complex networks (BA, ER, and WS models), a similar discussion on the periodic undulation of agent dynamics in differently structured networks, such as multifractal networks [41], would be interesting. This is because such a regular and hierarchical structure could resonate with the periodic undulation of agents and produce unobserved agent behavior in the whole network. Thus, we can say that studying such multifractal networks is a direction for future research. Furthermore, if we are targeting applied real-world cases, it would be better to consider a model where each node or edge has different weights [57] rather than focus on the standard model of three complex networks. How different results are obtained compared to the unweighted case treated in this study is also of interest for future studies. To create realistic biological or social networks for our multi-agent analysis, a methodological model that can infer unobserved parts of the network structure would be useful [58]. The study of agent dynamics with congestion-avoidance behavior in such a real network should be investigated in the future. In this way, our results have great potential and applicability to various problems involving agent mobility in the complex networks in related fields.

This study was supported by the JST-Mirai Program Grant Number JPMJMI20D1, Japan, as well as JSPS KAKENHI Grant Numbers JP21H01570 and 21H01352, and partly supported by MEXT as 'Program for Promoting Research on the Supercomputer Fugaku: Exploring Next Generation Aerospace Mobility and its Extension to Social System via Supercomputer Fugaku (Project ID: hp230198 and JPMXP1020230215)'. The authors thank Editage (www.editage.jp) for their English language editing services. The authors express special thanks to the administrative staff at Itoh Laboratory and RCAST.

All data generated and analyzed during this study are included in this published article. The data cannot be made publicly available upon publication because the cost of preparing, depositing and hosting the data would be prohibitive within the terms of this research project. The data that support the findings of this study are available upon reasonable request from the authors.

S T, D Y, E I, and K N designed the study; S T conducted the modeling and numerical experiments and data analysis; E I, D Y, and K N provided funding, project administration, and resources; and S T drafted the manuscript. All authors contributed to the feedback and editing of the manuscript.

The authors declare no conflicts of interest.