# Towards solving model bias in cosmic shear
forward modeling

Benjamin Remy, François Lanusse and Jean-Luc Starck
Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-
sur-Yvette, France
Corresponding author<EMAIL_ADDRESS>
###### Abstract

As the volume and quality of modern galaxy surveys increase, so does the
difficulty of measuring the cosmological signal imprinted in galaxy shapes.
Weak gravitational lensing sourced by the most massive structures in the
Universe generates a slight shearing of galaxy morphologies called cosmic
shear, key probe for cosmological models. Modern techniques of shear
estimation based on statistics of ellipticity measurements suffer from the
fact that the ellipticity is not a well-defined quantity for arbitrary galaxy
light profiles, biasing the shear estimation. We show that a hybrid physical
and deep learning Hierarchical Bayesian Model, where a generative model
captures the galaxy morphology, enables us to recover an unbiased estimate of
the shear on realistic galaxies, thus solving the model bias.

##  1 Introduction

Over the next couple of years, a new generation of wide-field optical galaxy
surveys will be coming online, starting with the ESA Euclid mission [1] and
the Vera C. Rubin Observatory Legacy Survey of Space and Time [2], soon
followed by the Roman Space Telescope [3]. At the core of the cosmology
science that these surveys will enable is the measurement of the weak
gravitational lensing effect, a very subtle shearing of the apparent images of
distant galaxies in the presence of massive structures along the line of
sight. On cosmological scales, this effect, also known as cosmic shear, is an
invaluable probe of the total matter content of the Universe, and will play a
key role in further constraining and testing our cosmological model over the
next decade.

Despite over 20 years of research in this field, unbiased measurements of
gravitational shear from observed galaxy images remains one of the most
challenging aspects of the scientific exploitation of these surveys. The
prevailing strategy for measuring shear from data has been to build estimators
of the ellipticity of galaxy light profiles, due to the fact that one can in
theory relate how a galaxy ellipticity $e$ changes under a gravitational shear
$\gamma$, as $e_{\text{obs}}\simeq e_{i}+\gamma$. Under the assumption that
galaxies should have no preferential intrinsic ellipticity, a simple estimator
for $\gamma$ can therefore be constructed from $\hat{\gamma}=\langle
e_{\text{obs}}\rangle=\langle e_{i}\rangle+\gamma$. However, to this day, no
shear measurement algorithm has proven to be unbiased.

In fact, state-of-the-art approaches for shear measurement [4, 5] no longer
attempt to build an unbiased shear estimator from images, but rather calibrate
the shape measurement methods estimating a multiplicative and an additive
bias. There are many factors contributing to these biases, but a fundamental
issue comes from the fact that the ellipticity $e$ is not a well defined
quantity for arbitrary galaxy light profiles.

In this paper we propose a new paradigm for cosmic shear estimation where we
no longer rely on ellipticities anymore. We propose to forward model the
galaxy images at the pixel level, similarly to [6], but using non-parametric
galaxy light profiles. The main idea is to train a deep generative model to
learn galaxy morophologies, which can then be used as part of a hierarchical
Bayesian model, a hybrid deep learning and physical model, simulating the
observations and used to infer the shear.

The purpose of this paper is first to demonstrate that using parametric light
profiles as part of the forward model induces a bias in the shear estimation,
which we call a model bias, and then to show that this bias can be removed by
using a deep generative model for the morphology. We show in this preliminary
study, that using a hybrid deep learning-physical Bayesian forward model can
not only solve the model bias, but also provide uncertainty quantification to
the shear estimation.

|
---|---
(a) Realistic galaxies fitted by a parametric profile | (b) Realistic galaxies fitted by a deep morphology model
Figure 1: The first row of images correspond to observed galaxies. The second
row corresponds to the light profiles Maximum a Posteriori fits, convolved
with their PSFs. The third row shows residual images between the observations
and the model fits. Column (a) is using Sérsic model as light profile and
column (b) is using deepgal model as light profile.

##  2 Forward modeling weak lensing observations at the pixel level

From a Bayesian perspective, making the joint inference of cosmic shear and
galaxy morphologies means estimating the joint posterior distribution
$p(\mathcal{G},\gamma|\mathcal{D})$, where $\mathcal{G}=\\{z_{i}\\}_{i=1\dots
N}$ is the set of morphology parameters and $\mathcal{D}=\\{d_{i}\\}_{i=1\dots
N}$ is the set of observed data. The Bayes identity tells us:

$p(\mathcal{G},\gamma|\mathcal{D})\propto
p(\mathcal{D}|\mathcal{G},\gamma)p(\mathcal{G},\gamma),$ (1)

where $p(\mathcal{D}|\mathcal{G},\gamma)$ is the data likelihood given by the
forward model described bellow. Since the intrinsic morphology of galaxies are
independent of the cosmic shear, then
$p(\mathcal{G},\gamma)=p(\mathcal{G})p(\gamma)$, where $p(\mathcal{G})$ and
$p(\gamma)$ are respectively the prior distributions of the morphological
parameters and of the cosmic shear.

We show the probabilistic graphical model in Figure 2. Nodes denote random
variables and arrows conditional relation between them. White nodes
corresponds to random variable on which we do the inference on, conditioned on
the shaded nodes corresponding to observed variables. Dot nodes correspond to
fixed variables, e.g. neural network weights $\theta$. In this paper we
compare two forward models: (1) using Sérsic law and (2) using the deepgal
model to generate the galaxy light profile. Once the galaxy light profiles are
drawn, we shear images with a constant value $\gamma$ shared by all the
galaxies. Each sheared galaxy $i$ is then convolved by its Point Spread
Function (PSF) $\Pi_{i}$. Assuming Gaussian noise corruption, the resulting
image is therefore evaluated by a Gaussian likelhood of variance
$\sigma_{\text{pix}}^{2}$.

Figure 2: Probabilistic Graphical Model.

Sérsic model  The light profile of the Sérsic model is given by the following
parametrization: $I(r)=F\times
I_{0}\exp(-b_{n}[(\frac{r}{r_{\text{hlr}}})^{\frac{1}{n}}-1]),$ where
$r_{\text{hlr}}$ is the half-light-radius, $I_{0}$ is a normalization factor
such that the integration of $I$ over all the space is $F$, which corresponds
to the flux of the galaxy. $n$ is called the Sérsic index and $b_{n}$ is a
dimensionless constant depending on $n$. The galaxy is then sheared with
parameters $e=(e_{1},e_{2})$ corresponding to its ellipticity. Finally, since
the real galaxy may not be perfectly centered in the postage stamp, we add two
shift parameters $s_{x}$ and $s_{y}$. Therefore, each galaxy $i$ has a set of
parameter $z_{i}=\\{n,r_{\text{hlr}},F,e_{1},e_{2},s_{x},s_{y}\\}_{i}$. Prior
distributions and parameters follow empirical distributions of the COSMOS
shape catalog and are described in subsection A.1.

Deepgal model  deepgal is a VAE-based model, proposed in [7], generating
complex galaxy light profiles. In this model, the observations are generated
following a random process mapping a latent space representation
$z_{i}\in\mathbb{R}^{16}$ to the observed image $d_{i}$ following the
probabilistic model illustrated in Figure 2. The neural network weights
$\theta$ of the VAE are trained maximizing the Evidence Lower Bound (ELBO) of
the data. Moreover, another generative model was trained on top of the VAE to
learn latent space distribution using Masked Auto Regressive Flows [8]. This
model was trained on the HST/ACS COSMOS survey postage stamps.

In the weak lensing regime, the shear is very small, therefore we can assume
that it does not change the morphological statistics of the training sample.
Thus, the model can be trained once, a priori, over a large catalog of
galaxies and be used to perform the joint inference of the latent variables
and the cosmic shear. For both light profile models, we denote the generated
image $\mathcal{I}_{0}$.

Cosmic shear  In this model every galaxy image is sheared by a constant value
$\gamma=\gamma_{1}+i\gamma_{2}$ shared by all the galaxies. The image shearing
transform, denoted $\mathcal{I}_{0}\otimes\gamma$, is given by the following
pixel coordinate change: $\big{(}\begin{smallmatrix}x^{\prime}\\\
y^{\prime}\end{smallmatrix}\big{)}=S(\gamma)\big{(}\begin{smallmatrix}x\\\
y\end{smallmatrix}\big{)}$, such that
$S(\gamma)=\frac{1}{\sqrt{1-|\gamma|^{2}}}\big{(}\begin{smallmatrix}1+\gamma_{1}&\gamma_{2}\\\
\gamma_{2}&1-\gamma_{1}\end{smallmatrix}\big{)}$.

Convolution with the PSF  The Point Spread Function (PSF) model is assumed to
be known and exact for each galaxy. We draw the PSF evaluated at the galaxy
$i$ position, denoted by $\Pi_{i}$, in Fourier space and then multiply it to
the Fourier image of the light profile before returning in the real space,
resulting in an image $\mathcal{I}$, such that
$\mathcal{\tilde{I}}=(\mathcal{\tilde{I}}_{0}\otimes\gamma)\cdot\tilde{\Pi}$.

Noise corruption We assume a known fixed gaussian noise identical for all the
galaxies. For simplicity, all of the galaxies in our observations are
corrupted by a constant Gaussian noise of variance $\sigma_{\text{pix}}^{2}$,
as well as for our forward model. In a more realistic scenario where we would
have information about the noise correlation properties for each postage
stamp, one could consider a noise model per galaxy. The likelihood is
therefore given by

$p(\mathcal{D}|\mathcal{G},\gamma)\propto\prod_{i=1,\dots,N}\mathcal{N}\left(d_{i}\Big{|}\mathcal{F}^{\dagger}\left[\left(\mathcal{\tilde{I}}_{0,i}\otimes\gamma\right)\cdot\tilde{\Pi}_{i}\right],\sigma^{2}\right),~{}~{}~{}\mathcal{F}^{\dagger}~{}\text{is
the inverse Fourier operator}.$ (2)

##  3 Illustration on realistic data

Figure 3: Shear marginal posterior contours against ground truth
$\gamma=(-0.03,0.03)$. The left panel shows the posterior using the Sérsic
light profile over the parametric observations. The middle panel shows the
posterior using the Sérsic light profile over the realistic galaxies. The
right panel shows the posterior using the deepgal model over the realistic
galaxies.

###  3.1 Simulating observations

We designed a toy model where galaxies are assumed to be already detected,
deblended and centered in postage stamps, so that we only need to simulate the
light profile, shearing transform and convolution with the PSF. In order to
demonstrate the model bias, we distinguish two sets of observations. (1)
generated with parametric galaxy light profiles using a Sérsic model and (2)
using realistic galaxies based on the galaxy postage stamps from the COSMOS
HST Advanced Camera for Surveys (ACS) field [9, 10, 11]. We used the postage
stamps from [12] and Sérsic fit parameters compiled in [13]. Because the
unconvolved galaxies are very noisy, we used the Auto-Encoder from [7] to get
denoised version of unconvolved galaxies, therefore this is not a
demonstration on real galaxies yet. However, the purpose of this paper is to
demonstrate the model bias, so we find that realistic galaxy morphologies that
cannot be captured by analytic profiles are sufficient to illustrate the
argument.

Each set of observation is a field of 256 galaxies, generated from the COSMOS
catalog, with the following cuts: $n>0.4$, $r_{\text{hlr}}>0.2$ and
$23.5<i^{+}<25.2$ so rather faint galaxies. That way the entire galaxy is in
the postage stamp, is not too peaky and not too resolved to be close to a
Sérsic model. Every postage stamps were corrupted by a Gaussian noise of
standard deviation $\sigma_{\text{pix}}=0.01$. In order to simulate shearing,
we removed any shear present in the selected sample from the COSMOS catalog,
by adopting the same strategy than in the GREAT3 Challenge [14], i.e. doubling
each galaxy, and rotating the doubled galaxies by 90°, as well as their PSF,
to cancel the overall sample shear and then apply the same constant shear to
every galaxy $\gamma=(-0.03,0.03)$.

###  3.2 Approximate inference with mean-field VI

In order to estimate the posterior distribution
$p(\mathcal{G},\gamma|\mathcal{D})$, we relied on approximate inference using
mean-field Variational Inference (VI), where we approximate the posterior
distribution using a parametric distribution
$q_{\phi}(\mathcal{G},\gamma)=q_{\phi_{\gamma}}(\gamma)\prod_{i}\prod_{j}q_{\phi_{i,j}}(z_{i,j})$,
parameterized by $\phi$. We used Normal distributions $q_{\phi_{i,j}}$ to
approximate the posteriror of parameters $i$ of the $j^{\text{th}}$ galaxy, so
the parameters $\phi_{i,j}=(\mu_{i,j},\sigma_{i,j})$ correspond to means and
standard deviations. Similarly, $q_{\phi_{\gamma}}$ is a 2d-diagonal
multivariate Gaussian. The parameters $\phi_{i,j}$ and $\phi_{\gamma}$ are
then optimized maximizing the ELBO:

$\phi^{*}=\arg\max_{\phi}~{}~{}\mathbb{E}_{(\mathcal{G},\gamma)\sim
q_{\phi}}\left[\log p(\mathcal{G},\gamma,\mathcal{D})-\log
q_{\phi}(\mathcal{G},\gamma)\right]$ (3)

This remains a very simplistic approach, which does not assume any correlation
between variables. A more accurate model would condition the morphological
parameters on the shear. In principle, the most accurate method would be to
run a Markov Chain Monte Carlo, but several experiments showed that sampling
based-methods are too computationally expensive for such high-dimensional
problem.

###  3.3 Identifying and solving the model bias

We ran three sets of experiments to demonstrate the model bias. First fitting
the fully parametric model to the parametric observations. Then fitting the
same parametric model to the realistic galaxies and then fitting the hybrid
model to the realistic galaxies. We used a Sérsic law as light profile as
proposed in [6]. One could argue that there exist more complex parametric
profiles such as the bulge+disk model which could reduce the aforementioned
bias, but we propose here a model which does not rely on ellipticities and in
principle can capture any complexity of galaxy morphologies.

First, as illustrated in Figure 3 (a), the inferred shear using the parametric
model on the parametric observations is unbiased. This was indeed expected
since the observations and the model are exactly the same. Secondly, Figure 3
(b) shows that the parametric model induces a bias in the shear estimation
when the galaxy morphologies are no longer elliptical. Besides, one can
observe the non-zero residuals in the Figure 1 fit to real galaxies,
illustrating the model discrepancy when applied to complex morphologies.
Finally, Figure 3 (c) shows that when we substitute the parametric light
profile model by the deepgal model, the shear estimation becomes unbiased. One
can notice that the uncertainties using the generative model are smaller than
for the Sérsic model, but these contours are not representative of the error
bars due to the mean-field approximation.

##  4 Discussion and conclusion

We have presented in this work the necessity of an accurate light profile
model for forward modeling of the cosmic shear. We have shown that a hybrid
physical and deep learning HBM has the potential to solve the model bias,
being able to capture the complex morphology of galaxies and simulating
physical processes such as shearing transform and convolution by the PSF. We
identify many directions of research towards a real cosmological analysis by
forward modeling. On the one hand, many assumptions could be relaxed by
building a more complex hierarchical model, e.g. parametrizing the PSF model
and jointly inferring it along with the shear. Concerning the detection and
deblending of galaxies, similar work have already been proposed, e.g. BLISS
[15]. On the other hand, accurate uncertainty quantification could be enabled
by scaling MCMCs for HBMs or using neural density estimator to capture the
shape of the joint posterior distribution.

##  5 Acknowledgments

BR acknowledges support by the Centre National d’Etudes Spatiales (CNES). This
work is supported by a public grant overseen by the French National Research
Agency (ANR) through the program UDOPIA, project funded by the
ANR-20-THIA-0013-01. This work was granted access to the HPC resources of
IDRIS under the allocation 2022-102038 made by GENCI.

##  6 Broader impact

The methodology described here finds many useful applications within and
outside of cosmology. Hierarchical Bayesian Models are very efficient tools
for statistical inference but can suffer form model misspecification.
Combining deep generative models and physical models seems to be a promising
approach to solve this problem. We believe this work does not entail any
negative consequences or ethical issues.

## References

  * [1] Rene Laureijs, J Amiaux, S Arduini, J-L Augueres, J Brinchmann, R Cole, M Cropper, C Dabin, L Duvet, A Ealet, et al.  Euclid definition study report.  arXiv preprint arXiv:1110.3193, 2011.
  * [2] Zeljko Ivezić et al.  LSST: From Science Drivers to Reference Design and Anticipated Data Products.  The Astrophysical Journal, 873(2):111, 3 2019.
  * [3] D Spergel, N Gehrels, C Baltay, D Bennett, J Breckinridge, M Donahue, A Dressler, BS Gaudi, T Greene, O Guyon, et al.  Wide-field infrarred survey telescope-astrophysics focused telescope assets wfirst-afta 2015 report.  arXiv preprint arXiv:1503.03757, 2015.
  * [4] Eric Huff and Rachel Mandelbaum.  Metacalibration: direct self-calibration of biases in shear measurement.  arXiv preprint arXiv:1702.02600, 2017.
  * [5] Erin S Sheldon and Eric M Huff.  Practical weak-lensing shear measurement with metacalibration.  The Astrophysical Journal, 841(1):24, 2017.
  * [6] Michael D Schneider, David W Hogg, Philip J Marshall, William A Dawson, Joshua Meyers, Deborah J Bard, and Dustin Lang.  Hierarchical probabilistic inference of cosmic shear.  The Astrophysical Journal, 807(1):87, 2015.
  * [7] François Lanusse, Rachel Mandelbaum, Siamak Ravanbakhsh, Chun-Liang Li, Peter Freeman, and Barnabás Póczos.  Deep generative models for galaxy image simulations.  Monthly Notices of the Royal Astronomical Society, 504(4):5543–5555, 2021.
  * [8] George Papamakarios, Theo Pavlakou, and Iain Murray.  Masked autoregressive flow for density estimation.  Advances in neural information processing systems, 30, 2017.
  * [9] Anton M Koekemoer, H Aussel, D Calzetti, P Capak, M Giavalisco, J-P Kneib, A Leauthaud, O Le Fevre, HJ McCracken, R Massey, et al.  The cosmos survey: Hubble space telescope advanced camera for surveys observations and data processing.  The Astrophysical Journal Supplement Series, 172(1):196, 2007.
  * [10] Nick Scoville, H Aussel, M Brusa, Peter Capak, C Marcella Carollo, M Elvis, M Giavalisco, L Guzzo, G Hasinger, C Impey, et al.  The cosmic evolution survey (cosmos): overview.  The Astrophysical Journal Supplement Series, 172(1):1, 2007.
  * [11] Nick Scoville, RG Abraham, H Aussel, JE Barnes, A Benson, AW Blain, D Calzetti, A Comastri, P Capak, C Carilli, et al.  Cosmos: Hubble space telescope observations.  The Astrophysical Journal Supplement Series, 172(1):38, 2007.
  * [12] Alexie Leauthaud, Richard Massey, Jean-Paul Kneib, Jason Rhodes, David E Johnston, Peter Capak, Catherine Heymans, Richard S Ellis, Anton M Koekemoer, Oliver Le Fevre, et al.  Weak gravitational lensing with cosmos: Galaxy selection and shape measurements.  The Astrophysical Journal Supplement Series, 172(1):219, 2007.
  * [13] Rachel Mandelbaum, Claire Lackner, Alexie Leauthaud, and Barnaby Rowe.  Cosmos real galaxy dataset.  Zenodo, January 2012.
  * [14] Rachel Mandelbaum, Barnaby Rowe, James Bosch, Chihway Chang, Frederic Courbin, Mandeep Gill, Mike Jarvis, Arun Kannawadi, Tomasz Kacprzak, Claire Lackner, et al.  The third gravitational lensing accuracy testing (great3) challenge handbook.  The Astrophysical Journal Supplement Series, 212(1):5, 2014.
  * [15] Derek Hansen, Ismael Mendoza, Runjing Liu, Ziteng Pang, Zhe Zhao, Camille Avestruz, and Jeffrey Regier.  Scalable bayesian inference for detection and deblending in astronomical images.  arXiv preprint arXiv:2207.05642, 2022.

##  Appendix A Appendix

###  A.1 Prior distributions and parameters for Sérsic light profile

The following distributions were fitted from the COSMOS shape catalog [12,
13].

  * •

$\log_{10}r_{\text{hlr}}\sim\mathcal{N}(-0.68,0.30)$

  * •

$\log_{10}n\sim\mathcal{N}(0.1,0.39)$

  * •

$\log_{10}F\sim\mathcal{N}(-1.97,0.53)$

  * •

$e_{1/2}\sim\mathcal{N}(0,0.28)$

  * •

$\gamma_{1/2}\sim\mathcal{N}(0,0.09)$

  * •

$s_{x/y}\sim\mathcal{N}(0,1)$


# Emergence of spatial transitions in urban congestion dynamics

A.L.Aniello Lampo  J.B.-H.Javier Borge-Holthoefer  S.G.Sergio Gómez
A.S.-R.Albert Solé-Ribalta  Internet Interdisciplinary Institute (IN3),
Universitat Oberta de Catalunya, Barcelona, Catalonia, Spain  Departament
d’Enginyeria Informàtica i Matemàtiques, Universitat Rovira i Virgili,
Tarragona, Catalonia, Spain  URPP Social Networks, University of Zurich,
Zurich, Switzerland

###### Abstract

The quantitative study of traffic dynamics is crucial to ensure the efficiency
of urban transportation networks. The current work investigates the spatial
properties of congestion, that is, we aim to characterize the city areas where
traffic bottlenecks occur. The analysis of a large amount of real road
networks in previous works showed that congestion points experience spatial
abrupt transitions, namely they shift away from the city center as larger
urban areas are incorporated. The fundamental ingredient behind this effect is
the entanglement of central and arterial roads, embedded in separated
geographical regions.

In this paper we extend the analysis of the conditions yielding abrupt
transitions of congestion location. First, we look into the more realistic
situation in which arterial and central roads, rather than lying on sharply
separated regions, present spatial overlap. It results that this affects the
position of bottlenecks and introduces new possible congestion areas.
Secondly, we pay particular attention to the role played by the edge
distribution, proving that it allows to smooth the transitions profile, and so
to control the congestion displacement. Finally, we show that the
aforementioned phenomenology may be recovered also as a consequence of a
discontinuity in the node’s density, in a domain with uniform connectivity.
Our results provide useful insights for the design and optimization of urban
road networks, and the management of the daily traffic.

Urban systems,

Road networks,

Congestion,

Phase transitions,

###### keywords:

Research

##  1 Introduction

The sustainability of urban life represents one of the greatest challenges of
our time. This has sparked a lot of interest among researchers coming from
several areas, and has led to the birth of the so-called science of cities
[1]. In this context, the engineering and optimization of road transportation
networks stands as a crucial task to ensure the efficiency of traffic dynamics
and to control the related congestion emergence.

In the literature, the study of road transportation networks has been
addressed following different methods, according to the purpose and the system
characteristics. For instance, the phenomenology emerging on arterial roads
(also named inter-urban, peripheral or high-capacity roads), characterized by
long segments and limited inter-connections, has been treated by means of
fluid models [2], or the fundamental diagram of traffic flow [3, 4]. The case
of central streets (also called intra-urban), instead, has been traditionally
analyzed by means of network science [5, 6, 7].

These two types of road networks, which are usually studied independently, are
increasingly intertwined as cities sprawl over suburban areas. Therefore, an
entanglement perspective is of great importance to properly deal with the
dynamics of vehicles flux, as well as other related phenomena. In [8], we
pursued this path and focused on the spatial behavior of congestion onset. We
found that the structural discontinuity between the central and arterial roads
yields an abrupt transition of the location where congestion occurs. In
particular, congestion points shift away from the city center as larger areas
are included in the urban spatial domain. Remarkably, such property is also
detected in almost one hundred real cities road networks worldwide.

This result represents an important resource to improve the efficiency of road
networks. The displacement of the traffic bottlenecks towards the peripheral
zones, indeed, is crucial to reduce the pressure on the city center, and to
avoid the degradation of the transportation system. Therefore, the possibility
to control and manipulate such effect constitutes a notable task which has
never been performed before. The current paper goes in this direction.
Specifically, we present a series of extensions of the theoretical setup in
[8], aimed to approach more realistic situations, and shed light on the
fundamental ingredients ruling the observed phenomenon.

In [8], the spatial separation between urban center and the surrounding
periphery was recognized to be a necessary requirement for the emergence of
congestion transitions. There, the road transportation networks are treated as
the inextricable blending of two different structures, that are however
sharply separated in space. Here, instead, we look into the consequences of a
spatial overlap between center and periphery, namely we introduce a hybrid
region where the dense and highly regular frame of central roads coexists with
that of the arterial ones. This is a crucial point to match real contexts,
since the transition between center and periphery may not be sharp, but it is
expected to occur in a wider spatial range.

Along this line, we investigate the role of the edge distribution,
implementing the most general case of a non-uniform spatial-dependent density.
This turns to be a key object to manipulate the location of congestion onset,
allowing to smooth its shift from the center. The edge distribution, i.e., the
number of edges over the unit area, may also be manipulated by modifying the
number of nodes. Thus, we pay great attention to the situation in which
connectivity is assumed to be uniform, while urban center and periphery are
constructed as regions with a high and low concentration of nodes,
respectively. In this way, we recover the abrupt transitions of congestion,
and highlight the importance of both node and edge distributions to the
control of traffic bottlenecks.

The manuscript is organized as follows. In Sec. 2 we present the DT-MST random
planar model developed in [8], constituting a practical benchmark to deal with
the interplay of central and arterial roads. The former are described using a
Delaunay triangulation (DT), while the latter are treated using a minimum
spanning tree (MST), both embedded on separated spatial domains. In Sec. 3 we
employ this model to study the behavior of the congestion radius as a function
of the urban size, showing that it is possible to detect two regimes, which
correspond to bottlenecks located in the center, and in its connection with
the periphery, separated by an abrupt transition. In Sec. 4 we consider the
more realistic situation where MST and a DT are not sharply separated and can
share a finite portion of space. We find that this yields a second transition,
as a consequence of the further structural discontinuity introduced by the new
hybrid region. In Sec. 5 we extend the previous cases by including the
possibility of a non-uniform edge distribution, and show that it leads to
smoother transitions. Finally, in Sec. 6 we investigate the inverse situation
in which center and periphery emerge as areas with different node densities,
but maintain the same connectivity. This situation enables the recovery of the
abrupt transition and highlight the role of nodes distribution. The importance
of these results for urban planning, as well as possible future extensions,
are discussed in Sec. 7.

##  2 A model for center-periphery entanglement

Figure 1:  Graphical representation of a DT (left), MST (middle), and a DT-MST
random planar model (right). Networks are built over a uniform distribution of
$N=4000$ in a square with $L=45$ side. The DT-MST random planar model is
obtained by keeping the DT in a circle with radius $R^{DT}=24$, and the MST
outside.

Road transportation networks consist of networks whose edges represent roads,
while the points where they intersect constitute the nodes. Road networks are
a paradigmatic example of planar graphs, namely graphs which can be embedded
on a surface in such a way that no edges cross each other; of course, real
road networks may have crossing edges, but they always represent a very small
fraction of the whole system. In this section we present a general class of
planar graphs that, while being as free as possible of any microscopic
constraints, keeps the distinctive feature of an entanglement between central
and arterial roads.

The most natural choice to model urban center, characterized by a dense web of
streets resembling short segments, is provided by a Delaunay triangulation
(DT). For a given set of points in a plane, DT is a triangulation over the
plane, i.e., any triad of neighboring points defines a triangle, such that no
point lies within the circumcircle of any triangle. An example of DT is
depicted in Fig. 1 for a random uniform distribution of points. Roughly
speaking, DT yields a network structure characterized by a high number of
shorts edges, uniformly spread over all the nodes (regularity), with a great
redundancy of paths. Remarkably, there are no edge overlapping, i.e., DT
preserves planarity. Similar properties could be recovered by considering
other types of regular graphs, such as grids, but these would impose a spatial
constraint on the nodes positions, making the model contrived. However, as
widely discussed in [8], this kind of structure produces equivalents results
for congestions.

Peripheral arterial roads are high-capacity streets whose primary function is
to deliver traffic between different urban centers. As such, intersections are
often reduced to improve traffic flow, and the resulting peripheral street
frame resemble a sparse web of long segments. Still, we seek for a network
structure where these properties emerge automatically. The low connectivity
suggests to look into a tree structure, whose general implementation is
provided by a minimum spanning tree (MST).

The MST is defined as the subset of the edges (of a connected, edge-weighted
undirected graph) connecting the vertices in such a way that there are no
redundant paths, and minimizing the total edge weight. In our model, the
construction of the MST is obtained from a DT, in which the edge weights are
the reciprocals of their corresponding edge-betweenness, rather than the
Euclidean distance between their adjacent nodes. This is a crucial aspect to
reproduce the main routing backbones of arterial periphery webs, as shown in
Fig. 1 (middle). A distance-based MST, like that employed in [9], does not
consider routing flows, which are essential in transportation networks.

The model we propose to treat central and arterial peripheral roads in a
unified entangled framework is constituted by a DT core surrounded by a MST.
Summarizing, the resulting synthetic network is built according to the
following steps:

  * •

We consider a uniform random distribution of $N$ points (the nodes) in the
spatial domain $\mathcal{I}=\left[L,-L\right]\times\left[L,-L\right]$;

  * •

This set of points induces a DT over $\mathcal{I}$. We hold it within the
circle of radius $R^{(DT)}$, i.e., for points with radius $r\leq R^{(DT)}$;

  * •

Out of such a circle, for $r>R^{(DT)}$, we look for the MST resulting from the
DT, obtained by minimizing with respect to the reciprocals of the edge
betweenness of the links.

Our model, coined as DT-MST random planar model, is so defined by only three
parameters: the square spatial side length $L$; the number of nodes $N$; and
the DT radius, $R^{DT}$. In Fig. 1 (right) we present a realization of the
model for a given choice of these parameters.

The idea to model central and peripheral roads as a DT and a MST,
respectively, may be better understood in terms of betweenness distribution.
This has been carefully studied in [9] for a significant number of cities
worldwide, proving that it exhibits an invariant bimodal structure. Such
bimodal character is a first indicator of the existence of a structural
discontinuity between different topologies as we depart from the city center.
Specifically, as discussed in Sec. III of [8], the two branches of the
betweenness distribution may be fitted with those of a tree and a DT
respectively, providing an endorsement of the DT-MST model.

##  3 Spatial abrupt transitions of congestion

Figure 2:  Spatial behavior of congestion for a DT (top) and a MST (bottom)
networks, built on a uniform distribution of $N=6000$ nodes in a square with
side length $L=60$. Left column contains the graphic representations of the
networks. Hexagonal bins provide information about the occurrence of
congestion nodes in space, and their colors show the related frequency in
proportion to the intensity. In the right column we present the dependence of
the congestion radius $R_{c}$ on the patch radial size $R_{p}$, i.e., the
radius of the circle centered in the middle of the network, that is used to
define the considered subgraph within this radius. The congestion radius range
is divided in bins, and each point is located at the statistical mode
$R^{*}_{c}$ obtained with the distribution of $R_{c}$ after 100 realizations.
The color of the circular markers shows the probability of that value over the
experimental $R_{c}$ distribution. Shadow areas represent the variance of the
different realizations of $R_{c}$ values with respect to the bin average
value. Importantly, congestion mostly occurs in the center and no transitions
arise, suggesting that these are a consequence of the intertwining of
different graph structures.

We study first the behavior of congestion in the center-periphery intertwined
framework introduced above. The quantitative assessment of congestion may be
performed by recalling betweenness centrality. There exists a large set of
studies about traffic [5, 10, 11, 12, 13, 14, 15, 16] where the critical
injection rate of vehicles $\gamma_{c}$, i.e., the maximum rate at which
vehicles can enter the system without congesting it, can be evaluated as

$\gamma_{c}=\frac{N-1}{B_{n}^{\ast}}\,.$ (1)

where $N$ is the nodes number of the street network and $B_{n}^{\ast}$
represents the maximum node betweenness.

The analysis of congestion in terms of betweenness permits to evaluate how
congestion is affected by the underlying network structure. It is worth
highlighting that the definition of betweenness embodies a routing protocol
which implicitly assumes a traffic dynamic, relying on shortest-paths, and an
Origin-Destination matrix which in our case is assumed to be uniform and all-
to-all. Thus, combining these two ingredients the traffic model can be defined
with precision. Previous works [17, 18, 19, 20] show the exact correspondence
between the betweenness centrality as a measure to estimate node and edge
traffic load, and the results obtained through the use of agent-based models
(the actual traffic model).

Equation (1) suggests that the maximum betweenness node marks the congestion
onset. We are interested in its geographic location, that can be quantified by
the distance to the center (radius) of the maximum betweenness node, hereafter
called the congestion radius, $R_{c}$. In order to explore how the congestion
location is affected by the intertwining of arterial and central roads, we
calculate $R_{c}$ for subgraphs of the DT-MST synthetic networks, contained in
concentric circles with radius $R_{p}$. Low values for this patch radius
$R_{p}$, e.g., $R_{p}<R^{(DT)}$ refer to situations where the arterial
periphery is suppressed, while large values, i.e., $R^{(DT)}<R_{p}\lesssim L$
describe an urban context where arterial periphery is relevant.

Even though there is not a direct mapping, there are at least two different
urban processes that may be simulated by tuning the patch radius $R_{p}$: city
sprawl, and changes in daily traffic patterns. In this latter case, the patch
radius $R_{p}$ defines the city geographical domain where traffic dynamics is
concentrated, and out of which it is supposed to be negligible. Large values
for this quantity describe situations in which peripheral traffic flow is
significant.

The dependence of the congestion radius $R_{c}$ on the patch radius $R_{p}$ is
presented in Fig. 2 (right) for the DT and MST null models. In both
situations, the congestion radius is constant and exhibits a very low value
compared with the size of the system $L$. This means that congestion mainly
occurs in the city center, regardless of the selected urban area. Such
behavior is also reflected in the left panel, where hexagonal bins provide an
information about the location of maximum betweenness nodes.

Figure 3:  Spatial behavior of congestion nodes for the random planar models
introduced in Sec. 2 with $N=6000$, $L=60$, and $R^{(DT)}=24$. In the left
panel, the resulting network configuration, with MST and DT edges painted in
red and black, respectively. Hexagonal bins provide information about the
occurrence of congestion nodes in space, and their colors show the related
frequency in proportion to the intensity. In the right panel, we present the
dependence of the congestion radius $R_{c}$ on the patch radial size $R_{p}$,
i.e., the radius of the circle centered in the middle of the network, that is
used to define the considered subgraph within this radius. The congestion
radius range is divided in bins, and each point is located at the statistical
mode $R^{*}_{c}$ obtained with the distribution of $R_{c}$ after 100
realizations. The color of the circular markers shows the probability of that
value over the experimental $R_{c}$ distribution. Shadow areas represent the
variance of the different realizations of $R_{c}$ values with respect to the
bin average value.

In Fig. 3 we repeat the same experiment but for the DT-MST random planar model
introduced in Sec. 2. The relationship between congestion radius and patch one
is not constant anymore, but depicts a step pattern configuration, namely an
abrupt transition between two regimes associated to congestion at the DT
center and at its connection with the MST. The former emerges for small
patches size, while the latter arises at large ones, i.e., when the arterial
roads contribute significantly to the network traffic. The comparison with
Fig. 2 suggests that this unexpected behavior is due to the interplay of the
composing structures, and in particular to the spatial discontinuity produced
in the network. Note that one could also implement the DT-MST entanglement by
homogeneously overlapping the edges of the former an arboreal skeleton,
without any spatial separation, and, as discussed in [8] (Fig. S10), no
transition is detected. This is a crucial aspect that will be explored in the
rest of the paper.

##  4 Center-periphery entanglement with spatial overlap

The spatial discontinuity between the two topologies, resulting from their
embedding on separated regions, is a fundamental ingredient to obtain the
abrupt transition. So far we have considered a sharp separation between DT and
MST. However, this is an extremely idealized situation, since the transition
between a dense center to a sparse periphery usually occurs in a wider spatial
range.

In this section we maintain the general separation between DT and MST but
include the possibility of a spatial overlap in a finite interval. This is
obtained by introducing a hybrid region with intermediate connectivity
properties. The new region is ring-shaped, with small and large radius
respectively given by $R^{(DT)}$ and $R^{(DT)}+\Delta R$, where $\Delta R>0$
denotes the width of the ring. In this region, we cover the MST skeleton with
a density $\pi\in\left[0,1\right]$ of DT edges, sampled from a uniform
distribution. The extreme bound $\pi=0$ refers to the situation in which no
edge has been added, while $\pi=1$ corresponds to a completely saturated
region.

This class of networks is shown in the left column of Fig. 4, for different
values of the parameter $\pi$. The corresponding relationship between the
congestion radius and the patch one is presented in the right column. The
cases with extreme values $\pi=0$ and $\pi=1$ fall into the category treated
in the previous section, with DT radius equal to $R^{(DT)}$ and
$R^{(DT)}+\Delta R$, respectively. These exhibit one abrupt transition between
two regimes related to congestion in the center and in its frontier with the
periphery.

Interesting phenomenology arises for intermediate values of $\pi$. In these
cases, the hybrid region is not saturated and differs from both the DT and
MST. This induces a further structural discontinuity in the network, which
affects the spatial behavior of congestion. Particularly, a second abrupt
transition stands out: as $R_{p}$ increases, congestion points abruptly shift
from the center to the DT frontier ($R_{c}\approx R^{(DT)}$), and then to that
of the hybrid ring ($R_{c}\approx R^{(DT)}+\Delta R$). The spatial overlap,
thus, induces a further possible congestion location.

This experiment sheds light on the mechanism producing the abrupt transition.
The structural discontinuity required to obtain such a phenomenon consists of
a sharp change of the spatial density of edges, i.e., number of edges over the
unit area. Then, one may suppose that such quantity permits to manipulate the
abruptness of the congestion transition; this is the aim of Sec. 5.
Additionally, the spatial edge density can also be modified by acting on the
nodes distribution, rather than network connectivity; this will be explored in
Sec. 6.

Figure 4:  Spatial behavior of the congestion nodes for the setup introduced
in Sec. 4. The figure is organized similarly to Figs. 2 and 3: left panels
show the network configurations, where blue edges refer to the $\pi$ density
of further $DT$ links we add on top of the MST skeleton in the $[R,R+\Delta
R]$ hybrid region. Still, hexagonal bins describe the occurrence of congestion
in space. Right panels presents the dependence of congestion radius on the
patch one. Here, circles are located at the statistical mode over an ensemble
of 100 realizations of different random distributions of $N=3000$ nodes in a
square of side equal to $L=20$, with $R^{(DT)}=5$, and $\Delta R=5$. The
figure highlights the emergence of a further intermediate congestion regime,
as a consequence of the hybrid region marked in blue.

The intermediate congestion regime emerged in Fig. 4 portrays a similar
pattern to that detected in [8], although in a different framework: GT-model.
This consists of square grid, representing the city center, connected to four
regular trees, one at each side, representing the urban periphery. In this
context, three congestion regimes arise, localized respectively in the grid
center, grid side and the tree root. This situation is equivalent to that
discussed in the current section. the GT-model, indeed, defines three regions
with different edges density: the areas where are embedded the grid, the trees
and the grid-tree connections. The discontinuity between these regions lead to
a multiple congestion pattern.

##  5 Center-periphery entanglement with non-uniform edges distribution

In the previous section we recognized the relationship between the observed
phenomenology and the connectivity properties of a road network. Particularly,
abrupt transitions follow from sharp discontinuities in the edge density. This
suggests that bottlenecks location may be controlled by properly tuning such
quantity. This is the purpose of the current section.

So far, the spatial discontinuity in the edge density emerged spontaneously as
a consequence of the intertwining of different topologies. However, within a
given region the connectivity has been considered constant in space.
Particularly, the hybrid region introduced in Sec. 4 was constructed by
covering the MST skeleton with DT edges sampled from a uniform spatial
distribution. This constitutes a strong idealization because road density is
in general not uniform over the city. Specifically, the transition from a
dense center to a sparse arterial periphery depicts a non-trivial profile of
the street distribution, interpolating between that of the DT and the MST one.

We propose an extension of the setup in Sec. 4, where the DT edges in the
hybrid region are sampled from a non-uniform distribution. We assume that an
edge with center at radius $r$ has a probability to be added given by

$P(r)=\alpha\exp(-kr)+\beta,\quad r\in\left[R^{(DT)},R^{(DT)}+\Delta
R\right],$ (2)

where $k\neq 0$ is a real constant ruling the damping in space of the
exponential. The $\alpha$ and $\beta$ coefficients are fixed by the boundary
conditions

$P(R^{(DT)})=1,\quad P(R^{(DT)}+\Delta R)=0\,,$ (3)

which lead to

$\alpha=\frac{1}{e^{-kR^{(DT)}}-e^{-k(R^{(DT)}+\Delta
R)}},\quad\beta=\frac{1}{2}\left[1-\frac{e^{-kR^{(DT)}}+e^{-k(R^{(DT)}+\Delta
R)}}{e^{-kR^{(DT)}}-e^{-k(R^{(DT)}+\Delta R)}}\right].$ (4)

Additionally, for generalization purposes and consistency with the method in
Sec. 4, we ensure that the final network includes a fraction $\pi$ of all the
candidate edges.

This procedure introduces a bias in the construction of the hybrid region,
whose edge connectivity interpolates now between that of the DT and the MST.
In Fig. 5 (left) the grey solid line refers to Eq. (2) for a given choice of
the parameter $k$. The blue histogram portrays the final distribution of a
fraction $\pi=0.5$ of the DT edges: their density approximately reproduces the
profile in Eq. (2) and is not constant in space but decreases exponentially as
a function of the radius. Of course, when $\pi=1$ the hybrid region saturates
and gets indistinguishable from the DT core, which extends until
$r=R^{(DT)}+\Delta R$. In this case the system reduces to the DT-MST model in
Sec. 2. In the right panel we present the dependence of the congestion radius
$R_{c}$ as a function of the patch size $R_{p}$. Black (purple) dashed lines
refer to $\pi=1$ ($\pi=0.05$), corresponding to a hybrid region in which all
the edges (approximately no one) have been added, and reproduces the abrupt
transition pattern discussed in Sec. 2. A new behavior arises at $\pi=0.5$:
the transition between different congestion regimes is not abrupt anymore.
Congestion radius, approximately equal to zero at small parches, smoothly
grows until $r=R^{(DT)}+\Delta R$, as a consequence of the continuous edge
density degrowth shown in the left panel.

Figure 5:  Congestion radius pattern for the setup introduced in Sec. 5, with
$N=3000$, $L=20$, $R^{(DT)}=5$, and $\Delta=10$. Hybrid region is generated by
overlapping a fraction $\pi=0.5$ of DT edges sampled from the distribution in
Eq. (2), with $k=0.1$. In the left panel, the blue histogram describes the
final distribution of these edges. The right panel reports the congestion
radius $R_{c}$ as a function of the parch radius $R_{p}$. The black and purple
dashed line report the result related to $\pi=1$ and $\pi=0.05$ respectively,
where all the edges and approximately no one have been added, and match the
situation described in Sec. 4. The dotted line refers to the edge distribution
shown in the right panel, where circles are located at the statistical mode
over an ensemble of 100 realizations. In this case transition exhibits a
smooth, rather than abrupt, profile as a consequence of the non-tight
discontinuity in the edge density.

A similar behavior is detected in Fig. 6, for a different choice of
parameters. Both make use of $k=0.1$, and this deserves to be commented. High
values for $|k|$ make the function in Eq. (2) elbow-shaped, leading to a
quasi-sharp second jump in the edge density, which would reproduce the
situation analyzed in Sec. 4. In general, the density profile in Eq. (2), with
the coefficients in Eq. (4), approaches a constant function with value 1 for
$k\rightarrow-\infty$, and with value 0 for $k\rightarrow+\infty$. These two
cases fall into the situation described in Sec. 4, with a DT radius equal to
$R^{(DT)}+\Delta R$ and $R^{(DT)}$, respectively.

Figure 6:  Congestion radius pattern for the setup introduced in Sec. 5, with
$N=6000$, $L=20$, $R^{(DT)}=5$, and $\Delta=10$. Hybrid region is generated by
overlapping a fraction $\pi=0.4$ of DT edges sampled from the distribution in
Eq. (2), with $k=0.1$. In the left panel, the blue histogram describes the
final distribution of these edges. The right panel reports the congestion
radius $R_{c}$ as a function of the parch radius $R_{p}$. The black and purple
dashed line report the result related to $\pi=1$ and $\pi=0.05$ respectively,
where all the edges and approximately no one have been added, and match the
situation described in Sec. 4. The dotted line refers to the edge distribution
shown in the right panel, where circles are located at the statistical mode
over an ensemble of 100 realizations. In this case transition exhibits a
smooth, rather than abrupt, profile as a consequence of the non-tight
discontinuity in the edge density.

##  6 Center-periphery entanglement with non-uniform nodes distribution

The analysis of the previous sections suggests that the existence of multiple
congestion regimes follows from different connectivity regions in the road
networks. In the most simple case of Sec. 2, this has been implemented by
entangling a DT core and a MST periphery, yielding a discontinuity in the edge
density and so to different degrees of paths redundancy. We argue that there
exist alternative implementations to reproduce this situation.

Figure 7:  Spatial behavior of congestions for the network model introduced in
Sec. 6. Left column shows the network configuration, where center and
peripheral edges are painted in black and red, respectively. Here, hexagonal
bins provide an information about the location of maximum betwenness nodes.
Right column presents the dependence of the congestion radius on the patch
one. Here, circles are located at the statistical mode over an ensemble of 100
realizations of different random distributions of nodes in a square of side
equal to $L=100$, with $R^{(c)}=5$, $\rho_{c}=\frac{N}{\pi(R^{(C)})^{2}}$,
($N=10000$). Rows refer to different density ratio. The figure depicts a clear
transition pattern as that detected in Fig. 3

.

So far, the network model has been constructed over an uniform distribution of
nodes, i.e., a constant nodes density. This leads to considering different
topologies in order to produce the wished change in the spatial density of
edges. One could also look into the inverse situation, in which center and
periphery result from different nodes densities, high in the center and low in
the periphery, but exhibiting the same connectivity properties (i.e., the same
average degree).

The network we consider in this section has a DT structure, and a different
nodes density between center and periphery. The former, defined as a circle
with radius $R^{(C)}$, has density $\rho_{c}$, while that of the latter is
$\rho_{p}<\rho_{c}$. In the limit in which $\rho_{c}=\rho_{p}$ the network
reduces to a uniform DT as that presented in Fig. 2. We assume that both
$\rho_{c}$ and $\rho_{p}$ are constant, i.e., they are uniform within their
own regions.

In Fig. 7 we present the spatial behavior of congestion for different center-
periphery density ratios $\rho_{c}/\rho_{p}$. It is possible to note a clear
abrupt transition between a regime describing congestion in the center, and
its connection with the periphery. We remark once more that the periphery is
not modeled as a MST, but the whole network results from a DT over the nodes
distribution. This confirms the possibility to produce an abrupt displacement
of traffic bottlenecks as a consequence of a discontinuity in the nodes
density, without any change in the topology.

An important aspect of the behavior in Fig. 7 concerns the pattern at
different density ratios $\rho_{c}/\rho_{p}$. The minimum patch radius value
allowing to observe congestion in the center-periphery connection, grows as
the density ratio increases. This is due to the fact that, in addition to a
structural discontinuity, the transition requires a large peripheral flux
(i.e., a sufficiently large number of nodes) in order to emerge. For a given
choice of the system parameters, such flux reduces as the density ratio
decreases. We expect, indeed, that when the peripheral number of nodes is low
enough, $\rho_{c}/\rho_{p}\rightarrow\infty$, no transition arises, since we
approach a network compound by a DT core with a few surrounding points.
Therefore, the experiment discussed in this section is constrained to a
certain range values of $\rho_{c}/\rho_{p}$: large enough to produce a
structural discontinuity, low enough to ensure a relevant peripheral flux.

##  7 Conclusions and Perspectives

In [8] we showed that traffic bottlenecks shift away from the city center as
larger concentric areas are included in the urban spatial domain, giving rise
to abrupt congestion transitions. The current work sheds lights on the main
mechanisms which permit to control the phenomena, and provides useful insights
in urban planning and the management of traffic.

The analysis in [8] points to the structural discontinuity between central and
arterial roads as the fundamental ingredient behind the abrupt transition of
congestion points. The present paper clarifies the formal meaning of such
structural discontinuity and relates it to the spatial density of edges. We
show that, acting on this quantity, it is possible to modify the transition
pattern, and so to control the geographic occurrence of congestion. Briefly,
our results may be summarized as follows. First, the introduction of
additional jumps in the edge density is associated to the emergence of new
abrupt transitions, (Sec. 4). Second, these transitions get smoother if the
space-density of DT edges on the MST skeleton exhibits a proper continuous
profile (Sec. 5). Third, abrupt transitions may be obtained by means of
discontinuities in the nodes density, without any change in the network
topology.

These results extend the applicability of the DT-MST model and allows to
better interpret the empirical validation presented in [8] (Sec. VII). There,
the congestion radius behavior as a function of the parch size was studied for
street networks of 97 cities worldwide. Overall results showed that 52 cities
present clear, detectable at naked eye, regimes with abrupt transitions
between them. The existence of multiple transitions, rather than a single one,
may be understood as the consequence of various edge density discontinuity
points in the road networks. Still, 45 cities present detectable regimes with
smoother transitions which can be now understood as the results of smooth
change in the edge density.

The possibility to control the spatial distribution of the congestion points,
and in particular to remove them from the center, constitutes an interesting
tool for urban planning. The displacement of the traffic bottlenecks towards
peripheral areas is a positive effect to reduce pressure on the center, which
is usually not prepared to support high vehicles loads. Within this context,
excessively abrupt shift of bottlenecks indicates areas where the structural
transition is excessively sharp, and this must have implications, yet to be
studied, with respect to the efficiency (in whatever terms) of the
transportation system.

The setup we have developed constitutes a complete benchmark to address
several problems related with the spatial localization of bottlenecks. First,
the role of edge weights deserves to be better investigated. In this work we
consider these as the Euclidean distances between the endpoints, but it would
be interesting to wonder if the congestion location may be manipulated by
means of the weights distribution. Secondly, concerning the reproduction of
the aforementioned empirical results, much effort is needed to understand
which actually is the effective origin of transition beyond edge density:
discontinuity in nodes distribution, or in connectivity? Our framework permits
to construct synthetic networks with nodes and edges properties sampled from
the corresponding empirical distributions, thus constituting a simple
playground to scrutinize the relationships between structure and function in
urban road networks.

## Declarations

## Availability of data and materials

The codes used during the current study are available from the corresponding
author on reasonable request.

## Competing interests

The authors declare that they have no competing interests.

## Funding

A.L., J.B-H. and A.S-R. acknowledge the support of the Spanish MICINN project
(grant PGC2018-096999-A-I00). S.G. acknowledges financial support from Spanish
MINECO (grant PGC2018-094754-B-C21), Generalitat de Catalunya (grant
2017SGR-896), and Universitat Rovira i Virgili (grant 2019PFR-URV-B2-41).

## Authors’ contributions

All the authors contributed equally.

## Acknowledgements

Not applicable

## Abbreviations

  * •

DT: Delaunay Triangulation

  * •

MST: Minimum Spanning Tree

## References

  * [1] Batty, M.: Building a science of cities. Cities 29, 9–16 (2012). doi:10.1016/j.cities.2011.11.008. Current Research on Cities
  * [2] Helbing, D.: Traffic theory from first principles. Available at SSRN 2960224 (2015)
  * [3] Godfrey, J.: The mechanism of a road network. Traffic Engineering & Control 8(8) (1969)
  * [4] Daganzo, C.F., Geroliminis, N.: An analytical approximation for the macroscopic fundamental diagram of urban traffic. Transportation Research Part B: Methodological 42(9), 771–781 (2008). doi:10.1016/j.trb.2008.06.008
  * [5] Guimerà, R., Díaz-Guilera, A., Vega-Redondo, F., Cabrales, A., Arenas, A.: Optimal network topologies for local search with congestion. Phys. Rev. Lett. 89, 248701 (2002). doi:10.1103/PhysRevLett.89.248701
  * [6] Solé-Ribalta, A., Gómez, S., Arenas, A.: Decongestion of urban areas with hotspot pricing. Networks and Spatial Economics 18(1), 33–50 (2018). doi:10.1007/s11067-017-9349-y
  * [7] Chen, J., Hu, M.-B., Li, M.: Traffic-driven epidemic spreading in multiplex networks. Physical Review E 101(1), 012301 (2020)
  * [8] Lampo, A., Borge-Holthoefer, J., Gómez, S., Solé-Ribalta, A.: Multiple abrupt phase transitions in urban transport congestion. Phys. Rev. Res. (in press), (2021)
  * [9] Kirkley, A., Barbosa, H., Barthelemy, M., Ghoshal, G.: From the betweenness centrality in street networks to structural invariants in random planar graphs. Nature Communications 9, 2501 (2018). doi:10.1038/s41467-018-04978-z
  * [10] Zhao, L., Lai, Y.-C., Park, K., Ye, N.: Onset of traffic congestion in complex networks. Physical Review E 71(2), 026125 (2005)
  * [11] Echenique, P., Gómez-Gardeñes, J., Moreno, Y.: Dynamics of jamming transitions in complex networks. Europhysics Letters (EPL) 71(2), 325–331 (2005). doi:10.1209/epl/i2005-10080-8
  * [12] Yan, G., Zhou, T., Hu, B., Fu, Z.-Q., Wang, B.-H.: Efficient routing on complex networks. Physical Review E 73(4), 046108 (2006)
  * [13] Liu, Z., Hu, M.-B., Jiang, R., Wang, W.-X., Wu, Q.-S.: Method to enhance traffic capacity for scale-free networks. Physical Review E 76(3), 037101 (2007)
  * [14] Dong, J.-Q., Huang, Z.-G., Zhou, Z., Huang, L., Wu, Z.-X., Do, Y., Wang, Y.-H.: Enhancing transport efficiency by hybrid routing strategy. EPL (Europhysics Letters) 99(2), 20007 (2012)
  * [15] Tan, F., Wu, J., Xia, Y., Chi, K.T.: Traffic congestion in interconnected complex networks. Physical Review E 89(6), 062813 (2014)
  * [16] Manfredi, S., Di Tucci, E., Latora, V.: Mobility and congestion in dynamical multilayer networks with finite storage capacity. Physical review letters 120(6), 068301 (2018)
  * [17] Solé-Ribalta, A., Gómez, S., Arenas, A.: Decongestion of urban areas with hotspot pricing. Networks and Spatial Economics 18(1), 33–50 (2018)
  * [18] Solé-Ribalta, A., Gómez, S., Arenas, A.: Congestion induced by the structure of multiplex networks. Phys. Rev. Lett. 116, 108701 (2016). doi:10.1103/PhysRevLett.116.108701
  * [19] Solé-Ribalta, A., Gómez, S., Arenas, A.: A model to identify urban traffic congestion hotspots in complex networks. Royal Society open science 3(10), 160098 (2016)
  * [20] Solé-Ribalta, A., Arenas, A., Gómez, S.: Effect of shortest path multiplicity on congestion of multiplex networks. New Journal of Physics 21(3), 035003 (2019)