Dynamical properties of model communication networks

We study the dynamical properties of a collection of models for communication processes, characterized by a single parameter $\xi$ representing the relation between information load of the nodes and its ability to deliver this information. The critical transition to congestion reported so far occurs only for $\xi=1$. This case is well analyzed for different network topologies. We focus of the properties of the order parameter, the susceptibility and the time correlations when approaching the critical point. For $\xi<1$ no transition to congestion is observed but it remains a cross-over from a low-density to a high-density state. For $\xi>1$ the transition to congestion is discontinuous and congestion nuclei arise.


I. INTRODUCTION
The interaction between the elements of social, technological, biological, chemical and physical systems usually defines complex networks. The study of topological properties in such networks has recently generated a lot of interest among the scientific community [1,2,3,4,5]. Part of this interest comes from the attempt to understand the behavior of technology based communication networks such as the Internet [6], the World Wide Web [7,8], e-mail networks [9] or phone call networks [10]. The study of communication processes is also of interest in other fields, notably in the design of organizations [11,12]. It is estimated that more than one-half of the U.S. work force is dedicated to information processing, rather than to make or sell things in the narrow sense [11].
Tools taken from statistical mechanics are used to understand not only the topological properties of these communication networks, but also their dynamical properties. Particularly interesting is the phenomenon of congestion. In has been observed, both in real networks [13] and in model communication networks, [14,15,16,17] that the system only behaves efficiently when the amount of information handled is small enough. The network collapses above a certain threshold and some information is accumulated and remains undelivered over large time periods-or it is simply lost. This transition from a free to a congested regime is indeed a phase transition and could be related to the appearance of the 1/f noise observed in Internet flow data [18,19].
Understanding congestion is also important from an engineering point of view. For instance, in October of 1986, during the first Internet collapse reported in the literature, the speed of the connection between the Lawrence Berkeley Laboratory and the UC at Berkeleyseparated by 200 meters-dropped by a factor of 100 [13]. Although the problem of congestion, and in particular its prevention and control [13], has been studied because of its implications in digital communications, a deep understanding of the physics of congestion for general communication processes and beyond particular protocols is still lacking.
A general collection of models that captures the essential features of communication processes has been recently proposed [16]. In these models, agents-nodesare organized in a hierarchical network-a Cayley treeand interchange information packets. Each agent has a certain capability that decreases as the number of packets to deliver increases. When the capability is inversely proportional to the number of accumulated packets, a continuous phase transition was found between the free and the congested phases. This transition was characterized by means of an order parameter.
The aim of this paper is to study the collection of models mentioned above, since only hierarchical lattices and a very particular congestion behavior-capability inversely proportional to the number of accumulated packetswere considered. An extension to consider the existence of cost to maintain communication channels has already been done [20]. First, the network model is extended beyond purely hierarchical Cayley trees. In particular, 1-dimensional and 2-dimensional lattices are considered. It should be noted that simple models such as Cayley trees and regular lattices can capture the main characteristics of the dynamics of information flow in complex networks. Hierarchical trees have been considered to model the TCP/IP protocol in the Internet [18,19]. On the other hand, computer based communication networks have been described by placing routers and servers in square lattices [14,17].
Moreover it is shown that, independently of the topology, the collection of models can be split into three groups according to how the network collapses. In the first group, agents deliver more packets as they are more congested-although their capability, as defined in [16], decreases-, and the network never collapses. In the second group, agents deliver always the same number of packets independently of their load-number of packets to deliver-; this behavior leads to a continuous phase transition as reported for hierarchical lattices in [16]. Finally, when agents deliver fewer packets as their loads increase, the transition to the congested phase is discontinuous and the network collapses in an inhomogeneous way giving rise to congestion nuclei. To characterize these behaviors, the order parameter defined in [16] is used, but also the power spectrum of the fluctuations and a susceptibility-like function. Most of the current effort is focused on the quantitative study of the continuous transition and the critical behavior associated to such a transition. The paper is organized as follows. Section 2 presents the collection of models and describes some details concerning the extension beyond hierarchical lattices, including 1D and 2D lattices. In section 3, these models are studied numerically and analytically. The scaling of the critical point with the size of the system, the behavior of the order parameter, and the divergence of the characteristic time at the critical point are studied in detail. Finally, section 4 includes the discussion of the results and the conclusions.

II. THE MODEL
The model considers three basic components: (i) the physical support for the communication process-agents and communication channels-, (ii) the discrete information packets that are interchanged, and (iii) the limited capability of the agents to handle such packets.
The communication network is mapped onto a graph where nodes mimic the communicating agents (for instance, employees in a company, routers and servers in a computer network, etc.) and the links between them represent communication lines. In particular, the three different topologies depicted in Fig. 1 are considered: 1D lattices of length L, 2D square lattices of side length L and size S = L 2 , and hierarchical Cayley trees with branching z and a number of generations or levels m, hereafter denoted (z, m).
The dynamics of the model is as follows. At each time step t, an information packet is created at every node with probability p. Therefore p is the control parameter: small values of p correspond to low density of pack-ets and high values of p correspond to high density of packets. When a new packet is created, a destination node, different from the origin node, is chosen randomly in the network. Thus, during the following time steps t + 1, t + 2, . . . , t + T , the packet travels toward its destination. Once the packet reaches the destination node, it is delivered and disappears from the network.
The time that a packet remains in the network is related not only to the distance between the source and the target nodes, but also to the amount of packets in its path. In particular, at each time step, all the packets move from their current position, i, to the next node in their path, j, with a probability q ij . This probability q ij is called the quality of the channel between i and j, and was defined in [16] as where k i represents the capability of agent i at each time step. The quality of a channel is, thus, the geometric average of the capabilities of the two nodes involved, so that when one of the agents has capability 0, the channel is disabled. High qualities (q ij ≈ 1) imply that packets move easily while low qualities (q ij ≈ 0) imply that it takes a long time for a packet to jump from one node to the next. The general equation proposed for k i is: where n i is the total number of packets currently at node i. The function f (n) determines how the capability evolves when the number of packets at a given node changes, with ξ ≥ 0. Equation (3) defines a complete collection of models with agents that behave differently depending on the exponent ξ. The average number of packets delivered during one time step by a node i to another node j is proportional to n i /(n ξ/2 i n ξ/2 j ). Assuming that n i ∝ n j the former expression reads n 1−ξ i . The proportionality is exact both in the hierarchical lattice and in large enough 1D and 2D lattices, where adjacent nodes are statistically equivalent.
Therefore, for ξ < 1 the number of delivered packets increases with the number of accumulated packets. For ξ > 1 the number of delivered packets decreases as the number of accumulated packets increases. Finally, for the particular case ξ = 1, the number of delivered packets is independent of the number of accumulated packets. Note that this particular case is consistent with simple models of queues [14].
The last point that needs to be explained to completely define the current model is the routing algorithm or, in other words, the set of rules that the nodes follow to select where to send a certain packet. In all the topologies considered (1D, 2D and Cayley) the packets follow paths of minimum length from their origin to their destination-open boundary conditions are set in both 1D and 2D networks. In 1D lattices and Cayley trees it is trivial to follow a minimum path because there is only one minimum path between two arbitrary nodes. In 2D lattices, however, there are many minimum paths. If a node can choose between two neighbors when sending a given packet, and both neighbors belong to a minimum path between the origin and the destination of the packet, one of them is chosen randomly with equal probability. This algorithm is indeed the simplest one and the interpretation of the results is clearer than for more complex routing algorithms. However, congested nodes will still have lower probability of receiving packets because of the definition of the quality of a channel. Therefore packets will avoid those nodes to some extent, as would happen in more sophisticated routing algorithms.

III. RESULTS
For certain parameters of the model (in particular ξ = 1 [16]) and for similar models of traffic [14], a transition from a free regime, where all the packets reach their destination, to a congested, regime where some packets are accumulated in the network, has been found. This transition is properly described by means of the following order parameter η(p) introduced previously [16]: In this equation ∆N = N (t + ∆t) − N (t) and . . . indicates average over time windows of width ∆t. These averages can be over one or many realizations, yielding the same result. Essentially, the order parameter represents the ratio between undelivered and generated packets calculated at long enough times such that ∆N ∝ ∆t. Thus, η is only a function of the probability of packet generation per node and time step, p. The power spectrum S(f ) of the total number of packets in the network N (t) = i n i is also used here to further understand the phase transition. By means of the power spectrum the behavior of the time correlations of the system can be studied.
Let us study separately the critical case ξ = 1 and the non critical cases ξ < 1 and ξ > 1.
A. The critical case ξ = 1 A continuous transition between the free regime and the congested regime occurs in hierarchical networks for ξ = 1, as reported in [16]. For small values of p, all the packets reach their destination and the total number of packets N (t) fluctuates around a finite value. In this case the order parameter is η = 0. However, as p increases, a critical point p c is reached, where the fluctuations in N (t) become very large and the characteristic time of the system diverges (critical slowing down). Beyond this point, some packets remain undelivered and N (t) grows linearly with t. The same qualitative behavior is observed here for 1D and 2D lattices, although there are quantitative differences.
From an engineering point of view it is interesting to study first the behavior of p c as the number of nodes in the network S changes, because it will provide valuable information about the scalability of the network. Note that N c = p c S is the maximum number of packets that the network can handle per time step and, thus, it is a measure of the capacity of the network.
For the hierarchical Cayley tree, a mean field calculation of p c was obtained in [16]: where the approximation holds for z m−1 ≫ 1.
It is also possible to derive a mean field expression of p c for the 1D lattice. Since the most congested node is, from symmetry arguments, the central one-the node at ℓ = L/2-, the network will collapse when the amount of packets received by this central node is higher than the maximum number of packets that it is able to deliver. Since in a large enough network it is safe to assume that the central node will be congested similarly to its neighbors, n ℓ−1 = n ℓ = n ℓ+1 , the maximum number of delivered packets should be 1. On the other hand, the number of packets arriving to the central node at each time step is the number of packets that are generated at each time step at the left half of the network and have their destination at the right half and conversely, this is pL/2. Then the critical condition is given by These mean field expressions (5) and (6) can be compared with simulations. Nevertheless, the fluctuations of N (t) become very large near p c and it is difficult to calculate the value of the order parameter. Instead, a susceptibility-like function χ(p) can be defined by analogy with equilibrium thermal critical phenomena, and used to estimate more accurately the value of the critical probability of packet generation, p c . The susceptibility χ is related to the fluctuations of the order parameter by where T is the width of a time window, and σ η (T ) is the standard deviation of the order parameter estimated from the analysis of many different time windows of width T . Thus a calculation implies a long realization of N (t), its division into windows of width T , calculation of the average value of the order parameter in each window and finally the determination of the standard deviation of these values. The susceptibility shows (Fig. 2) a singularity at p c as T grows, as expected.  (6) and (5) respectively.) Fig. 3 shows that there is good agreement between expressions (5) and (6) and the values of p c obtained numerically by means of susceptibility measures for different network sizes.
For the 2D lattice it is more difficult to obtain even a mean field expression for p c . However, since for 1D lattices and Cayley trees the scaling relation p c ∝ S −1 holds-where S is the size of the network-, one may expect the same behavior for the 2D lattice. Using the susceptibility to numerically determine p c from simulations, one finds that this turns out to be incorrect. Although it is difficult to obtain a precise value of the exponent, Fig 4 shows that it is close to 0.6 instead of 1.0: This result suggests that the existence of multiple paths to get from the origin to the destination has important consequences, not only in shifting the value of p c , but actually changing its critical scaling behavior. The behavior of the order parameter is studied next. It is possible to derive an analytical expression for the simplest 1D case where there are only two nodes that exchange packets. Since from symmetry considerations n 1 = n 2 , the average number of packets eliminated in one time step is 2, while the number of generated packets is 2p. Thus p c = 1 and with the present formulation of the model it is not possible to reach the super-critical congested regime. However, p can be extended to be the average number of generated packets per node at each step (instead of a probability) and in this case it can actually be as large as needed. As a result, the order parameter for the super-critical phase is η = (p − 1)/p. As observed in Fig. 5, the general form fits very accurately the behavior of the order parameter not only for this simple 1D lattice with L = 2 or the Cayley tree [16], but also for any 1D lattice. Two-dimensional lattices again deviate from this behavior, although the deviation is small. In particular, near p c , equation (9) implies and thus the critical exponent for the order parameter is equal to 1 at least for the 1D lattice and the Cayley tree.
The analysis of the power spectrum shows that in the sub-critical regime, i.e. in the free phase, the spectrum is well fitted by a Lorentzian characterized by a frequency, f c , corresponding to exponentially decaying correlations with a characteristic time τ , as depicted in Fig. 6. As p approaches p c , it is observed that f c → 0 and the power spectrum becomes 1/f 2 for the whole range of frequencies. Alternatively the characteristic time diverges as τ → ∞ (critical slowing down). This qualitative behavior is common to all the topologies of the underlying network, as shown in Fig. 6.  It is interesting to study how the characteristic frequency drops to 0 for each network topology. Near the critical point, one expects the scaling behavior The value of the critical exponent γ can be estimated by fitting Eq. (11) to values of f c (p) close enough to the critical point, as shown in Fig. 6. Note that we fit p c and γ simultaneously. This procedure yields very accurate values of p c but the values of γ are subject to large errors. Fig. 6 yields γ ≈ 0.9 for a 1D network with L = 100, γ ≈ 2.5 for a 2D network with L = 7, and γ ≈ 2 for a (7,5) Cayley tree. The determination of γ is interesting not only from an academic point of view, but also from an engineering perspective. Indeed, this exponent is related to divergences of other relevant quantities near the critical point. Any characteristic time τ -the average time to deliver a packet, for instance-will diverge as and similarly the total number of packets since from Little's law [21]. This law states that, in steady state, the number of delivered packets and the number of generated packets are equal.
The estimation of γ is particularly interesting in electronic communication protocols. Indeed, Eq. (12) is used to determine the waiting time before a packet is considered lost in the network and therefore sent again [13]. In practice, the exponent γ = 1 predicted by classical queue theory [21] is assumed, while our current results suggest that more complex settings can lead to exponents even larger than 2.
B. Non critical cases ξ < 1 and ξ > 1 We have shown that the number of packets delivered by node i is n 1−ξ i and thus, when ξ < 1, it increases with the number of packets that this node accumulates. It is difficult to imagine a real scenario with this characteristic. However, this case has been included to understand the critical behavior when ξ = 1, i.e., to show the relationship between criticality and the amount of packets that can be delivered when load increases. As a consequence of the increase of the deliver capability with the load, the transition to collapse will never occur because, at some point in time, the number of accumulated packets will be large enough and the number of delivered and created packets will balance each other. Thus, the order parameter will be zero for any value of the control parameter p, and the correlations will decay exponentially. As shown in Fig. 7, the characteristic frequency tends asymptotically to f * c as p increases. This asymptotic value depends on the size of the system.
For a 1D lattice with a high density of packets (p → 1), the number of packets that are delivered by a node is n 1−ξ i while the number of packets that are being delivered to this node is proportional to L (for instance, for the central node, this number is simply pL/2). Therefore, n i ∝ L 1/(1−ξ) . The total number of packets is N = i n i ∼ L 1+1/(1−ξ) and according to Little's Law On the other hand, for p → 0 the scaling of the characteristic frequency is given by since the packets success to jump from one node to the next at each time step, and therefore the characteristic time is directly the average path length between nodes. Therefore, although there is no phase transition in this case ξ < 1, there is a cross-over from a low density behavior to a high density behavior, as shown in Fig. 7. This crossover is also observed in 2D lattices and Cayley trees.
The phase transition observed for ξ = 1 is recovered when ξ > 1. Above a certain threshold, some packets are accumulated in the network and the order parameter differs from 0. However, the number of packets delivered by a node i in this case ξ > 1 decreases with the number of packets accumulated at that node. Therefore, when some packets are accumulated, n i grows and finally no packets are delivered at all. Thus suddenly above the transition, which is discontinuous, the order parameter becomes 1.
The change in the order of the phase transition affects the spreading of the collapse over the network. In the critical case ξ = 1, the collapse starts at the most central node and then spreads from this point to the rest of the network. In this case ξ > 1, the reinforcement effectthe fact that the more collapsed a node is, the more collapsed will get in the future-leads to the formation of many congestion nuclei generated by fluctuations, that spread over the whole network. Figure 8 illustrates the formation of these congestion nuclei for 2D lattices with ξ = 5 and p = 0.001, and ξ = 2 and p = 0.01

IV. DISCUSSION AND CONCLUSIONS
A collection of models recently proposed for hierarchical networks [16] has been examined in detail for several network topologies including 1D and 2D lattices, to characterize the phase transition to collapse. New congestion scenarios have been analyzed.
First, it has been shown that the congestion behavior is governed by the ability of agents to deliver information packets when their load increases. When agents deliver packets at a constant ratio independently of their loadnumber of packets to deliver-, a continuous transition as the one reported in [16] for Cayley trees is also observed for 1D and 2D lattices, the order parameter being the fraction of accumulated packets. When the number of delivered packets increases with the load, no phase transition is observed. Conversely, when the number of delivered packets decreases with the load, the order parameter jumps from zero to one and the transition becomes discontinuous. These different behaviors are tuned by a single parameter ξ, which is the exponent that determines how the capability of nodes evolves with the number of accumulated packets. When ξ < 1 there is no transition to the congested regime, for ξ = 1 the transition is continuous and for ξ > 1 the transition is discontinuous. Note that the continuous transition (reported in [16] and in some models of queues [15,17]) is only a particular case between a no-congestion behavior and a discontinuous transition behavior. Thus, the critical behavior is intimately related to the independence between load and deliver capability of nodes.
These different behaviors have been analyzed separately. The critical case ξ = 1 presents the most interesting behavior. To properly understand the physics of the collapse process it is necessary to know how the system approaches the critical point. The order parameter becomes useless in this region due to the emergence of large fluctuations, and then we have defined a susceptibilitylike function that is related, by analogy to equilibrium critical phenomena, to fluctuations of the order parameter. This quantity shows a peak at the critical point and allows to study the dependence on the size and topology of the system. For 1D lattices and Cayley trees, it has been shown (by means of simulation and also with mean field calculations) that p c scales as the inverse of the size of the system, p c ∝ S −1 . However, for 2D lattices a dif-ferent scaling p c ∝ S −0.6 is obtained. It is suggested that the existence of multiple paths from the origin to the destination of the packets is the responsible for this change in the scaling behavior. On the other hand, the study of the power spectrum of the number of the packets in the network as a function of time has provided valuable information about the variation of temporal correlations in the system. In particular, below p c correlations decay exponentially, with a characteristic time τ that diverges at p c . The exponent of this divergence is significantly larger than what is expected from queue theory. This can be relevant to accurately forecast the temporal behavior in complex communication settings.
For ξ < 1, no phase transition is observed. Instead, a cross-over from a low-density to a high-density regime occurs. In the low density regime, the characteristic frequency is determined by the average distance between nodes, while in the high density regime the characteristic frequency is determined by the capability of agents to deliver information packets. In 1D lattices, the cross-over changes the scaling from f c ∝ L −1 at low densities to f c ∝ L −1/(1−ξ) at high densities.
Finally, when ξ > 1 the transition to the congested regime is discontinuous, since the number of delivered packets decreases with the load and for long times no packets are delivered at all. Thus the order parameter jumps from 0 to 1. Moreover, due to the reinforcement effect-the fact that the more collapsed a node is, the more collapsed it will get in the future-leads to the formation of many congestion nuclei generated by fluctuations, that spread over the whole network. The existence and distribution of such congestion nuclei can also be of interest from an engineering point of view.