On the role of spatial dynamics and topology on network flows

Flow networks are inherently liable to congestion. The ability of these networks to handle demand at reasonable levels is crucial as otherwise a congested phase of transport affects the performance across the entire network. Therefore, it is of prime interest to analyze how and where networks begin to undergo a transition to a congested state and the dynamics of its response. Dissecting flow patterns is essential to address this problem. In this context, flow of data packets in the internet is well understood, as analyses of their traffic dynamics and phase transitions are abundant [1–6]. The transition point to congestion is in this case well established through analytical solutions and simulations [7–11]. To interpret the role of a network backbone for managing flows, optimal paths and minimum spanning trees have been studied [12–14]. The most relevant metric to determine the vulnerability of internet flows is the betweenness centrality for shortest paths, because it determines the critical element generating congestion in the network [15, 16].

All of these previous studies have overlooked the role of space and time in particle flows in the networks; a reasonable assumption in internet applications. In this case the travel time between nodes is negligible and queues accumulate in the nodes which have limited capacity to process the packets and congestion occurs. In contrast, in transportation networks travel time of vehicles or individuals is crucial. Queues are formed on the links which have limited spatial capacity, and in turn their velocity depends on the density of the travelers on each link. In this context, transportation research focuses on capturing traffic flow by making use of the fundamental diagrams that empirically relate flow, density and speed in road segments and utilizing them in macroscopic link models [17, 18]. The cell transmission model [19, 20] and the simple point-queue models (PQMs) [21–23] are well established among such traffic flow models. Alternatively, cellular automata models for vehicular traffic [24, 25] have also been used to mimic traffic flow behavior along with many other discrete stochastic models [26, 27]. In none of these cases, the interaction of the spatial dynamics with the network topology has been addressed. Our goal here is to build on the simple model of the internet by adding temporal and spatial dimensions of particle flow and hence construct a framework to analyze congestion in spatial networks from a network science perspective. The proposed framework can be extended to flows in other kinds of spatial networks [28–30], and more importantly could open new avenues of the research on urban road networks that would go beyond modeling and topological analysis for the statistical physics community [31–33].

We begin by recalling the scheme in [1–5, 7–11] and refer to it as the internet model (IM): the network is loaded with R identical particles at each timestep t with randomly assigned origins and destinations. A fixed shortest path routing table guides the particles toward their destination. The nodes can transmit as many particles per timestep as their outflow capacity, C, and travel between two nodes takes a unit timestep. Queues of particles form at the nodes, and they can grow infinitely large. The particles are exempt from joining the queue at their destination and are removed from the system upon arrival. The network response is measured by the order parameter H [9]

$$\begin{equation} H\left(R\right) = \lim \limits_{t \to T} \frac{\left\langle \Delta W \right\rangle}{R \Delta t}, \end{equation} \tag{ 1 }$$

where W denotes the number of particles in the system, $\left \langle \Delta W \right \rangle$ is the average change in the number of particles still in the system after a timestep, Δt is the unit timestep and T is the length of the simulation. Figures 1(a) and (b) depict the IM. For low values of R, the network reaches a rate of particle arrival equal to the loading rate. W remains constant and consequently H = 0. Conversely, if R exceeds a certain threshold R_c, a linear increase in W with a slope of H is observed due to excessive queueing. This behavior maps a second-order phase transition to congestion.

Zoom In Zoom Out Reset image size

Spatial networks, being embedded in two-dimensional space, give rise to three heterogeneities that need to be captured. The first one is non-uniform travel times. This differs from the internet where data packets hop from one node to the next in a single timestep. Second, these networks carry flows along the links. A particle in a spatial network has a specific position on the link it is traveling on. The third source of heterogeneity is a consequence of this: particles occupy physical space and gradually fill the segment. Among the various models of traffic flow in the transportation literature aimed to address these issues, the PQM used in [21–23, 34] is an adaptation of the IM that shifts flow from the nodes to the links, incorporating the non-uniform travel time distribution and thereby making the flow analysis very similar to that of the internet. The particles traverse the link freely by hopping through the τ slots of unit travel time to join a queue at the end of the link from which they will be discharged at the outflow capacity. The total travel time consists of the free travel time and the delay, namely, the timespan between the particle entering the queue and exiting it. The spatial PQM (SPQM) incorporates a single additional constraint to the PQM: every link has an upper limit for the number of particles it can hold at once. We will refer to this value as the volume capacity of the link, V . The links cannot accept any new particles when they reach their volume capacity, as illustrated in figure 1(c). This additional constraint has a crucial effect on the nature of the network response, as links at volume capacity clog the upstream links and cause them to succumb to congestion as well. This spreading of the congestion occurs at rates that depend on the loading rate and the network topology along with specific link properties. The rate of particles unable to travel determines the speed with which the congestion spills, which makes the spreading process non-binary unlike traditional spreading models in the literature. Although the SPQM share some aspects with several directed percolation models, the movement of non-identical particles along predetermined spatial shortest paths with non-binary spreading is uncommon and therefore relatively unstudied [35, 36].

The critical loading rate R^IM_c has been shown [7] to be equal to R^IM_c = N(N − 1)(C_max/B^N_max) where N is the network size measured by the number of nodes, B^N_max is the maximum node betweenness and C_max is the outflow capacity of this node. This relation arises from the fact that the inflow to a node is proportional to its betweenness centrality [37]. At R_c, the inflow is equal to the outflow at the node with the minimum C_max/B^N_max value. For the PQM, we adjust this equation by replacing the node betweenness by an edge betweenness value, R^PQM_c = N(N − 1)(C_max/B^E*_max), where B^E*_max refers to the maximum modified betweenness. This modification is necessary since the particles do not join the queues in the final links, which should be omitted from the edge betweenness calculations.

Next, we introduce a framework to analytically calculate the entire transition curve to the congestion. For R > R_c, particle inflow at certain elements will be larger than the outflow. We define Rⁱ_c as the critical loading rate specific to element i. For large R, the outflow of the congested links is maximized to capacity, which in consequence affects the inflow to the links downstream. To account for this we define the delay factor, D_i(R), referring to the fraction of paths through i that are not suffering from delay as

$$\begin{equation} D_{i}(R) = \mathcal{H} ( R_{\rm c}^{i}-R)+\frac{C_{i}}{I_{i}(R)} \mathcal{H} ( R-R_{\rm c}^{i}), \end{equation} \tag{ 2 }$$

where $\mathcal {H}(x)$ is the Heaviside step function and C_i and I_i(R) are the outflow capacity and the inflow of element i for loading rate R. D_i(R) = 1 suggests no congestion for element i, whereas lower values indicate the levels of congestion. By using this definition, the inflow I_i(R) at a specific loading rate R can be quantified as

$$\begin{equation} I_{i}(R) = \frac{\sum_{k \in \Gamma(i)} \prod_{j \in k} R D_{j}(R)}{N(N-1)}, \end{equation} \tag{ 3 }$$

$$\begin{equation} H(R) = \sum_{i \in \mathcal{N}} (I_{i}(R)-C_{i}) \mathcal{H}(I_{i}(R)-C_{i}), \end{equation} \tag{ 4 }$$

where Γ(i) is the set of paths passing through element i. Equation (4) accounts for all of the delay factors of the elements upstream of element i by going through the shortest paths. Equations (2) and (3) form a set of coupled equations that can be solved to obtain the inflows for every element. H(R) is obtained by summing all of the positive values of I_i(R) − C_i. Figure 1(b) also reveals the exact solution for the simple network around the critical point with the dashed curve H = (7/42)(1 − R_c/R). In order to test different network topologies and examine the effect of space on criticality, we use a non-periodic lattice as a substrate and rewire each edge (i,j) with probability ρ to a new destination j* chosen with the probability proportional to d(i,j*)^α where d(i,j) denotes the Euclidean distance [38, 39]. Figure 2 reveals that the simulations and the analytical results perfectly coincide for transitions in both a two-dimensional non-periodic lattice of N = 1225 and C = 4, and its completely rewired instance. Critical rates decrease for smaller ρ as newly introduced shortcuts have higher betweenness values. As ρ increases, the value of α and its effect on network topology becomes more pronounced: rewired links in networks with lower α values are more localized and therefore maximum betweenness values are higher in these networks. Consequently the increase in H is sharper for lower α. Figures 2(b) and (d) show sharper transitions and lower R_c for decreasing α values for network instances with ρ = 1.

Zoom In Zoom Out Reset image size

At R_c, the link that triggers congestion, also referred to as the critical element, is expected to fluctuate between free flow and congested phases. Figure 3(a) exhibits the frequency distribution of the timespans at which this most critical element operates at its outflow capacity, as an indicator of the temporality of the phase transition. The results show that these timespans follow a power law with exponents of −0.58 ± 0.04 for the IM and −0.48 ± 0.04 for the PQM, independent of the network topology. SPQM exhibits a different behavior. For low volume capacities the links tend to fill up, causing the links upstream to fill as well. The fluctuations cause a gridlock, a condition where all of the elements of a cycle are completely filled and hence flows stop. In case of a gridlock, the order parameter increases very sharply. In figure 3(b), for varying volume capacities we measure the average number of timesteps it takes for a gridlock to occur, t^g, normalized by the length of the simulation T. It can be observed that for lower V , the gridlocks are observed relatively quickly. For large volume capacities, PQM and SPQM have the same R_c, which we will refer to as the PQM-limit. In either model, at steady state, the queues are not expected to be necessarily empty but rather steady in their size. If the volume capacity of a link is smaller than this steady state queue size, the particles will be blocked in the upstream link which consequently may suffer from a decrease in its outflow due to this clogging effect. This suggests a lower critical point for a network with active volume capacity constraints, correspondingly R^PQM_c ⩾ R^SPQM_c. Therefore, to realize the congestion-free transport to the fullest, the PQM-limit should be aimed at. Figure 3(c) illustrates the response of a non-periodic lattice with N = 625 for decreasing levels of volume capacity and reveals that both the critical loading rate and the nature of the transition are affected by V . Figure 3(d) shows the effect of the volume levels on R_c. In PQMs, the number of particles on a link is the sum of those in the queue, which was shown to be proportional to the modified betweenness, and those that are traveling. The expected inflow to a single slot of travel on a link is proportional to the actual betweenness of that link, hence the expected number of particles traveling is τB^E. Therefore for small V and strictly deterministic inflows, a link reaches its maximum volume when $V=R_{\mathrm { c}}\left (\tau B^{\mathrm { E}}+B^{{\mathrm { E}}*}\right )_{{\mathrm { max}}}/N(N-1)$. The inflow to the link is proportional to the actual routing betweenness hence the critical point can be expressed as

$$\begin{equation} R_{\rm c} = N(N-1)\min \left( \frac{C}{B^{{\rm E}*}_{{\rm max}}}, \frac{V}{(\tau B^{\rm E}+B^{{\rm E}*})_{{\rm max}}} \right), \end{equation} \tag{ 5 }$$

where B^E is the edge betweenness. Consequently a linear increase in R_c is observed up to the point where V =C(τB^E + B^E*)_max/B^E*_max. However, the stochastic nature of the model, along with the first-order nature of its transition, causes fluctuations that force R_c/C below this analytical bound as it converges to the point-queue limit.

Zoom In Zoom Out Reset image size

Figure 4 depicts the transitions for the PQM and the SPQM for the San Francisco road network with N = 1152 and an average degree of 3.2. The network is discretized by unit travel times of 10 s. The outflow capacity of a road segment is obtained by using the speed limit and the number of lanes. Volume capacities are estimated for every road segment assuming that the volume capacity is reached when the speed drops to half of the speed limit. Under these assumptions, the PQM-limit is not reached as R^PQM_c = 40 (14 400 vehicles h⁻¹) and R^SPQM_c = 30 (10 800 vehicles h⁻¹). To capture the network response in the SPQM, segment volumes are recorded at different timesteps of the simulation t = 720 (2 h), t = 1440 (4 h) and t = 2160 (6 h) for R = 36, slightly above the critical load. Figure 4 illustrates the network response by mapping road segments that have reached 80% of their volume capacities in the given time periods by color. Congestion originates from an artery leading to the downtown area and anisotropically spreads to other regions. After 6 h, most network elements are suffering from congestion.

Zoom In Zoom Out Reset image size

Our findings suggest that in the SPQM, the critical road segment cannot accept incoming vehicles as it saturates to its volume capacity before it reaches its outflow capacity. This outcome can be traced back to the inherence of congestion in downtown areas: cities with high population densities have concentrated spatial demand distributions, which result in the inadequacy of urban space to accommodate such concentrated flows. This work is a step further on a systems analysis applied to congestion in roads. In further studies, population and facility distributions can be modeled. In the applications domain an open question is to know how the studied transition to congestion is influenced by introducing empirical origin–destination matrices that represent the real population's travel demand [40, 41] and how network topology contributes toward reducing- or maybe aggravating-congestion. While in the domain of phase transitions, analyzing the set of critical exponents using directed percolation as a benchmark case remains an intriguing and open question [35, 36].

This work was funded by New England UTC Year 24 grant, awards from NEC Corporation Fund, the Solomon Buchsbaum Research Fund, BMW, UPS and the National Natural Science Foundation of China (no. 51208520). The authors would like to thank Zoltán Toroczkai and Ruben Juanes for useful comments and technical support.