# Emergence of complex network topologies from  flow-weighted optimization of
network efficiency

Sebastiano Bontorin  Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo
(TN), Italy  Department of Physics, University of Trento, Via Sommarive 14,
38123 Povo (TN), Italy  Giulia Cencetti  Fondazione Bruno Kessler, Via
Sommarive 18, 38123 Povo (TN), Italy  Riccardo Gallotti  Fondazione Bruno
Kessler, Via Sommarive 18, 38123 Povo (TN), Italy  Bruno Lepri  Fondazione
Bruno Kessler, Via Sommarive 18, 38123 Povo (TN), Italy  Manlio De Domenico
<EMAIL_ADDRESS>University of Padua, Via Francesco Marzolo 8,
35131, Padua, Italy  Padua Center for Network Medicine, University of Padua

###### Abstract

Transportation and distribution networks are a class of spatial networks that
have been of interest in recent years. These networks are often characterized
by the presence of complex structures such as  central loops paired with
peripheral branches, which can appear both in natural and man-made systems,
such as subway and railway networks. In this study, we investigate the
conditions for the emergence of these non-trivial topological structures in
the context of human transportation in cities. We propose a minimal model for
spatial networks generation, where a network lattice acts as a spatial
substrate and edge velocities and distances define an effective temporal
distance which quantifies the efficiency in exploring the  urban space.
Complex network topologies  can be recovered from the optimization of joint
network  paths and we study how the interplay between a  flow probability
between two nodes in space and the associated travel cost influences the
resulting optimal network. In the perspective of urban transportation we
simulate these flows by means of human mobility models to obtain Origin-
Destination matrices. We find that  when using simple  lattices, the  obtained
optimal topologies transition from tree-like structures to more regular
networks, depending on the spatial range of  flows.  Remarkably, we find that
branches paired to large loops structures appear as optimal structures when
the network is optimized for an interplay between heterogeneous mobility
patterns of small range travels and longer range ones typical of commuting.
Finally, we  show that our framework is able to recover the statistical
spatial properties of the  Greater London Area subway network.

##  I Introduction

Cities represent one of the most fascinating man-made complex systems,
exhibiting complex features ranging on different scales: from their structure
and dynamical behavior, up to the scaling of socio-economic factors with their
size Batty (2008); Barthelemy (2019); Bettencourt and West (2010); Pan  _et
al._ (2013); Arcaute  _et al._ (2015). These features represent a strong hint
towards the existence of universal underlying mechanics behind apparently very
different cities Bettencourt (2013); Bassolas  _et al._ (2019); Bettencourt
(2020). Out of these structural properties, one of the most relevant, as it
plays a fundamental role mediating the complex interplay  between human
dynamics Schläpfer  _et al._ (2021); Barbosa  _et al._ (2018) and mobility in
urban context, are transportation networks Alessandretti  _et al._ (2022); Lee
_et al._ (2017); Gallotti  _et al._ (2021); Barthélemy (2011); Morris and
Barthelemy (2012). These networks are a class of spatial networks whose
properties have been investigated in the literature during the  last two
decades Barthélemy (2011); Gastner and Newman (2006). In particular, they have
been studied under the lens of optimality conditions and minimization of cost-
based functionals Gastner and Newman (2006), in order to identify specific
features behind efficient networks. The concept of optimal networks Barthelemy
(2019) and energy-like minimization Barthélemy and Flammini (2006) has its
natural understanding in the physics language. States of a system which
minimize a functional defining trade-offs between system’s observables (e.g.,
free energy) represent the most likely to be observed states of many real
world systems. While in some complex systems, such as cities, these physical
variables can not be derived from first principles, these analogies and
concepts can still offer a valid perspective and provide an embedding of these
systems in a space where the interplay between their structure and dynamics
can be unfolded and better understood. Simple laws have been studied Louf  _et
al._ (2013); Gastner and Newman (2006) to better understand the emergence of
hierarchy and the role of traffic in the network state. Moreover, global and
local optimization criteria lie in the evolution of man-made systems where
policy makers and planners can adopt some of these criteria in their plans
Barthélemy (2011).

Transportation networks are often characterized by the presence of complex
structures Banavar  _et al._ (1999); Pei  _et al._ (2022); Gallotti and
Barthelemy (2015) such as loops paired with branches Roth  _et al._ (2012),
which can appear both in natural Tero  _et al._ (2010) and man-made systems
Barthélemy (2011),  like railway and subway networks. These structures
represent the key topological elements behind efficient public transportation
systems Pei  _et al._ (2022). In this study we investigate the conditions for
the emergence of these non-trivial structures Louf and Barthelemy (2014); Louf
_et al._ (2013) in the context of human transportation in cities. We aim to
reconstruct these topologies by means of an optimal configuration Ibrahim  _et
al._ (2021) of the network state. Under the assumption of a fixed total cost
and a limited set of high-capacity connections (e.g., a constraint in the
expenditure available on infrastructure), the optimal configuration is the
assignation of connections’ velocities, or edges’ weights, such that the joint
amount of time required to travel between two nodes is minimized for all pairs
of nodes. Moreover, as these networks represent the arteries in urban
exploration/navigation via public transportation, we study the role of traffic
Zhang  _et al._ (2015); Louf  _et al._ (2013) in weighting the importance of
specific set of connected edges (paths). We model the urban morphological
structures which generate heterogeneous distributions of human mobility in
space, biasing these optimal networks to converge towards specific topological
features. We aim to explore the minimum requirements and the conditions for
these optimization processes to reproduce the empirical structures
aforementioned.

At variance with the recent works on network efficiency, we adopt some
fundamentally different modeling choices. We evaluate the efficiency in terms
of time necessary to explore the network, where edges’ weights act as  travel
speeds. We also weight path efficiency by the traffic probability between
nodes. The underlying network lattice (as represented in its simplest form by
the triangular lattice in the next sections) acts as a substrate that allows
the network to evolve Szell  _et al._ (2022) and possibly exhibit the network
topological features typical to real world systems. On this framework, we show
how  introducing simple probabilities biasing the optimal efficiency between
points in space force a transition between a tree-like topology and a network
resembling a simple lattice. We show also that the modeling of traffic-like
flows forces the emergence of preferential shared paths in space. The optimal
configuration of these shared paths leads to complex topologies, which
ultimately shows features seen in real systems. Features such as a bi-modality
in the edges’ velocity distribution, characteristic of multi-layered
transportation, and the central core with loops paired with branches typical
of subway systems Roth  _et al._ (2012); Pei  _et al._ (2022) are recovered.
We finally show an application of the model within the  Greater London Area,
finding similarities of the optimal model with its  London Underground
network.

Figure 1:  Spatial network model for urban morphology: A) Mapping population
distribution and urban transportation network to a minimal spatial network
where nodes encode urban features. Example with hexagonal tiling mapped to the
triangular lattice. B) Network-based distance $c_{\Pi_{ij}}$ versus euclidean
distance $d_{ij}$; edges weights/velocity $w_{e}$ are depicted as widths. C)
Edges weights of the lattice substrate are optimized via simulated annealing
to unravel spatial features of the optimal transportation network.

##  II Framework for Urban Spatial Morphology

We introduce here a general framework for spatial networks with the aim of
recovering a minimal model for urban morphologies that encode both
transportation properties and urban features such as population and density of
points of interest (POIs). To this aim, we begin from the definition of a
network substrate which acts as an effective discretization of the spatial
dimension. Its simplest form can be found in an hexagonal 2-dimensional tiling
Birch  _et al._ (2007) and its planar dual, the triangular lattice. More
formally, in this network substrate each tile in space is represented by a
node, connected to its set of neighboring nodes (see Fig. 1). The existence of
a physical edge between nodes/tiles $i$ and $j$ is encoded in the adjacency
matrix $A$ where $A_{ij}=1$ if the regions are neighbors in the lattice.
Distances and metrics are therefore computed on top of this network substrate
and are biased by nodes’ features.

Figure 2:  Optimization of synthetic networks: The role of $\beta$ is studied
for two  network models: the triangular lattice (HEX) and the non-spatial (ER)
network. A) Heatmap of target nodes probabilities $p_{ij}$ from source node
(yellow) under two different $\beta$ values: as the penalty parameter grows,
farther nodes are more penalized and  flows tend to stay close to the source.
B) Samples of the associated optimized network states: when  flows are not
affected by distances ($\beta$ = 0.1) source nodes target all the other nodes
in the network with approximately equal probability, the optimal network
converges to a tree-like structure. With larger $\beta$ ($\beta$ = 5.0), trip
probabilities are more localized and the presence of loops appear in the
optimal structure.

Nodes of this network  can encode spatial features at the urban scale, such as
population or amenities’ density in a given node. We therefore have a minimal
representation of a urban morphological structure (see Fig. 1), and a network
substrate that acts as a transportation system and can be optimized to
generate optimal transportation networks Szell  _et al._ (2022). The path-
based temporal distance on top of the transportation network acts as the
fundamental metric we aim to minimize. The rationale behind a network-based
distance is grounded on the assumption that in the context of public
transportation, urban systems are not navigated by considering geographical
distance but rather by  evaluating the travel-time between  departure and
arrival. More specifically, multi-layer transportation networks Gallotti and
Barthelemy (2015); Alessandretti  _et al._ (2022) are characterized by layers
having a hierarchical organization with different characteristic speeds
Gallotti  _et al._ (2016). Thus, an effective temporal distance becomes
fundamental in determining accessibility and efficiency in urban space
exploration.

Figure 3:  Loop Dimension vs $\beta$: Minimum cycle basis is used as a
network’s observable to study the appearance of loops. For each point the
median and its absolute deviation are shown. A) The average size (number of
edges) of the loops that constitute the cycle basis. B) Cycle basis dimension
as number of loops. The optimal network ranges from a tree structure to a
lattice-like with small loops, as the probability of long range movement
decreases (large $\beta$). A transition in the cycle basis property is
observed at  $\beta\sim 1.0$ for the triangular lattice under study, where the
optimal network results in an intermediate state with large loops.

In this model, we denote $e$ as an $edge$ in the network, and $w_{e}$ as the
associated edge’s weight which can be seen as a  velocity of the edge in the
transportation network. $d_{e}$ is the euclidean distance of edge $e$ between
the nodes it is connecting; here edge weights are visually mapped as widths of
the links. Information about edges distance in this framework can be relevant
when generalizing to the case of random spatial networks where edges have
different lengths. In the case of a general non-spatial network, where there
is no notion of spatial distances, the model can be adapted by fixing
$d_{e}=1$. Finally, we define $\Omega_{\Gamma_{ij}}$ as the set of paths
connecting the two nodes. We then maximize the efficiency of this underlying
substrate. The transportation efficiency between two nodes $i-j$ is computed
as a cost in terms of time Ren  _et al._ (2014),  and we do not account for
congestion, which can be introduced in a possible extension of this framework.
We find the path (a set of connected edges starting from source node $i$ and
ending in destination node $j$) with the smallest cumulative time, where the
time delay introduced by choosing en edge $e$ is measured as a ratio between
its euclidean distance and its weight, representing a proxy of speed. Refer to
Fig. 1 for a graphic depiction. Here $G(\\{w_{e}\\})$ is used to indicate the
network configuration with the associated set of edges weights $\\{w_{e}\\}$.
We therefore aim to find the assignment of weights $\\{w_{e}\\}$ as a trade-
off between network efficiency in transportation, by minimizing the set of
costs $\\{c_{ij}\\}$ of travelling between pair of nodes $i-j$ where each
element $c_{ij}$ is defined as:

$c_{ij}(\\{w_{e}\\})=\min_{\Pi\in\Omega_{\Gamma_{ij}}}\left[\sum_{e\in\Pi_{ij}}\frac{d_{e}}{w_{e}}\right],$
(1)

and in absence of further information, the optimization procedure is the
equivalent of minimizing the travel costs $\sum_{ij}c_{ij}$. Here we add a
novel ingredient, in which we couple the optimization of the network temporal
distances with a traffic  flow or probability between pairs of nodes.
Operationally, when dealing with real world Origin-Destination (OD) matrices,
this probability can be then mapped to a traffic $T_{ij}$ between two points.
$T_{ij}$ represents the probability of a person from node $i$ to travel to
node $j$, and a traffic can be recovered when information about populations in
source and target nodes is added. $T_{ij}$ effectively acts as a rank in the
importance of a specific path in the network. As paths connecting different
pairs of nodes may share common edges of the network substrate, complex
topologies emerge from the shared paths jointly optimizing the network
efficiency. The  flow-weighted transportation efficiency therefore becomes:

$E(G(\\{w_{e}\\}))=\frac{1}{N(N-1)}\sum_{i}^{N}\sum_{j\neq i}^{N}T_{ij}\cdot
c_{ij}(\\{w_{e}\\})$ (2)

We also require that the total network  infrastructure cost, defined as the
cumulative sum of edges weights per unit length, multiplied by edge distance
$C_{G}=\sum_{e\in G}d_{e}w_{e}$ is conserved. This is a generalization of a
standard optimization process,  in the sense that when $T_{ij}=1$,
$\forall(i,j)$, the efficiency is optimized for all possible trip pairs
$(i,j)$ with equal importance, where the Minimum Spanning Tree often
represents the optimal solution Barthélemy (2011).

Before tackling the problem of traffic-like (OD) flows, we study a simpler
definition of $T_{ij}$, which allows to understand the role of  distance in
the optimization process, in absence of other nodes’ features:

$T_{ij}\propto e^{-\beta d_{ij}}.$ (3)

The coefficient $\beta$ appearing in Eq. 3 is introduced as a penalizing
parameter and determines how relevant is the pair-wise distance $d_{ij}$ when
computing probabilities. We can understand it as the inverse of a
characteristic traveling distance for an agent on the network
$\beta\sim\frac{1}{d_{0}}$. While several alternatives on the integration of
distance in spatial-dependent probabilities (such as power-laws $T_{ij}\sim
d_{ij}^{-\gamma}$) can be employed, we focus on the exponential dependence as
it represents the foundational result from the maximum entropy derivation of
gravity flows Wilson (1975). The introduction of gravity-like flows will be
discussed in Section IV.

In the following, we introduce the application of the model on simple
substrates to explore  the role of $\beta$ in absence of spatial urban
features.

##  III Optimization of simple network substrates

In order to asses the role of the characteristic distance parameter $\beta$ in
the emergence of specific topologies, we compute networks statistics on a set
of generative models for both spatial and non spatial networks. As hexagonal
tiling of space is preferable when an analysis includes aspects of
connectivity Birch  _et al._ (2007), the first model we study is a triangular
lattice. The  reason behind this choice is that it represents the planar dual
Barthélemy (2011); Viana  _et al._ (2013) of the hexagonal lattice. Therefore,
as space is discretized in hexagonal tiles, the spatial network connecting its
centers is the triangular lattice,  which is isotropic and presents less
equivalent degenerate paths of a rectangular lattice. As a direct reference to
hexagonal tiling, we refer to this network as HEX (see Fig. 2). We also extend
the analysis also to the case of a random network model where nodes are not
embedded in a metric space. Specifically, we study an Erdős-Rényi (ER) network
topology, where the definition of distance between nodes $L_{ij}$ can be
defined in terms of topological shortest path distance Newman (2010).

Figure 4:  Minimal models of urban morphology and attractiveness distributions
under study (3-Points and Exponential decay): Morphology of POIs, where
attractiveness $W_{j}$ is mapped with color intensity (yellow being higher).
Optimized edges weights distributions $P(\\{w_{e}\\})$ are characterized by
the bi-modal nature that reveals the multi-layered structure of the optimal
transportation networks when close-range flows are paired with long-range
traffic typical of commuting towards city centers (Insets $P(T_{ij})$ with
peaks on large-flows due to POIs). Gaussian KDE is shown in orange as a visual
aid. A) 3-points polycentric distribution of POIs,  resembling the euclidean
Steiner Tree problem Dreyfus and Wagner (1971); Brazil  _et al._ (2013) for
three points. The network is optimized with $\beta=0.1$ and $\beta=4.0$,  and
shows the appearance of branches connecting the POIs paired to large loops in
the periphery. B) Optimal state and distribution of speeds with exponential
decay of $w_{j}$ from the center and an exemplifying result with $\beta=3.0$.
The optimal topology is characterized by a central loop paired with branches.

As a first benchmark we  simplify flows as a spatial probability
$T_{ij}=p_{ij}$ that decays exponentially with distance and does not consider
nodes features, the resulting equation for  $p_{ij}$ is:

$p_{ij}=\frac{e^{-\beta d_{ij}/\langle d\rangle}}{\sum_{k\neq
i}^{N-1}e^{-\beta d_{ik}/\langle d\rangle}}$ (4)

Where $\langle d\rangle=\frac{1}{N(N-1)}\sum_{i\neq j}d_{ij}$ is the average
distance of points in the network and acts as a normalization factor
(euclidean distance $\langle r\rangle$ in case of a spatial network or
topological $\langle L\rangle$ for the ER network).

Therefore $p_{ij}$ encodes how much of the nearby space is explored by a
single source node. An illustration of the spatial dependence of target
probabilities and samples of the resulting optimal topologies are presented in
Fig. 2.

For a range of $\beta$ values the optimization process is performed on an
ensemble of these models. To assess the emergence of complex structures we
observe the number of loops that emerge in the optimal state. This measure is
relevant in the context of spatial networks, where loops break the symmetry
introduced by optimal structures such as trees. We compute the minimum cycle
basis set as a metric to observe the emergence of loops Kavitha  _et al._
(2007): i.e. the minimum set of loops ( where a single loop is encoded in a
set of edges that defines a  closed path in the graph) such that any other
closed path in the network can be reconstructed via combination of this cycle
basis Kavitha  _et al._ (2007). Specifically, we investigate the cycle basis
dimension (the number of loops that constitute this set) and the average loop
size, against a range of $\beta$ values.  This metric allows to quantify the
emergence of spatial topological features that differentiate the optimal state
from a tree structure. Results for these synthetic systems are presented in
Fig. 3. Additional boxplots are shown in  SM Figures 1-2. A tree-like topology
is recovered when the  flows probabilities are distributed uniformly across
all nodes in space (when $\beta\rightarrow 0$ and distance is  therefore not a
penalizing variable in Eq. 4), while loops emerge when farther targets become
less likely to be explored and the network is globally optimized for close-
range trips.  Notably, in Fig. 3 around $\beta\approx 1.0$, we observe a sharp
transition in the average loop size  in the HEX lattice under analysis:
connections appear between neighboring nodes which are far from the tree-root
as it becomes more efficient to have a direct link.  In this $\beta$ regime
the tree topology does not guarantee the most efficient configuration for
peripheral nodes, which have their high probability targets in their direct
neighborhood (see Eq. 4). Thus in the optimization process edges appear
between leaves nodes which are in separated branches: this ultimately breaks
the tree structure and leads to the emergence of large-scale loops. Eventually
the optimal network converges to a simpler structure with small loops as the
network is optimized for nodes to target only direct neighbors in the lattice.
Finally, in  SM Section 2 we show an application on the case of a single
target node in the perimeter of the lattice, where the model reproduces leaves
venation patterns Barthélemy (2011); Katifori  _et al._ (2010).

##  IV Spatial attractiveness and traffic-like flows

In the context of urban systems, optimal transportation networks need to be
devised to accommodate traffic  flows Zhang  _et al._ (2015) towards specific
areas of interest, e.g. due to high commercial and business land use density.
Hence we extend the efficiency optimization framework in the case where we
have more realistic traffic on top of the urban networks, as the presence of
nodes with high attractiveness (POIs) biases the  flows towards them. In urban
scenarios we adopt spatial-interaction models to mimic more traffic-like
flows.  In these models, flows are obtained via a gravity-like equation:
$T_{ij}\propto p_{i}p_{j}\exp{(-\beta d_{ij})}$ Barbosa  _et al._ (2018) which
can be derived from first principles via entropy maximization, thus
representing the most likely set of  flows to be observed. In the context of
urban exploration, the gravity equation can be mapped to a model for spatial
interaction Piovani  _et al._ (2017); Wilson (1975) where nodes with a given
attractiveness $W_{j}$ compete as possible targets for traffic:

$T_{ij}\propto\frac{1}{Z}P_{i}W_{j}\exp\left(-\beta d_{ij}\right)$ (5)

Normalization $Z$ accounts for all possible trip alternatives
${\sum_{k}W_{k}\exp\left(-\beta d_{ik}\right)}$. $P_{i}$ is the population
density in node $i$ and $W_{j}$ encodes a suitable definition of
benefit/attractiveness of node $j$ as a possible target Piovani  _et al._
(2017). $T_{ij}$ is therefore the fraction of population in node $i$
commuting/travelling on average to node $j$. To better understand the role of
nodes’ attractiveness, we start with the simplest assumption of equal
population distribution on all nodes: $P_{i}=1.0~{}\forall i$; we will
introduce more realistic population distribution in the next section with the
London case study.

Figure 5:  Optimal network model for  Greater London Area subway system:
Application of the efficiency optimization with realistic  flows on the urban
structure of the  Greater London Area. A) Urban morphology data is recovered
from Census and OSM and population and POI densities are mapped to the H3
tiling. B) Data is mapped to the triangular lattice, with nodes having
features which allow the calculation of traffic-like  flows, a sample OD
matrix is shown where $T_{ij}$ are computed with $\beta=0.35$. C) Optimal
network state for the London model, where only edges and nodes corresponding
to the second mode are shown ( see SM Section 3). Central core structure with
loops paired with peripheral branches can be visually seen.

We apply these models on the triangular lattice to unravel the optimal
topologies that emerge  when traffic probabilities are biased towards some
nodes having high attractiveness (simulating POIs) and we study two spatial
configurations for nodes’  $W_{j}$. In the first configuration high $W_{j}$ is
assigned to three nodes (POIs) placed at the vertices of an equilateral
triangle.  We study the 3-points distribution as it mimics a prototypical
polycentric distribution of city-centers, and it can be linked to the solution
of the euclidean Steiner Tree problem Dreyfus and Wagner (1971); Brazil  _et
al._ (2013). The Steiner Tree is a class of problems where given a set of N
points in a plane the goal is to find the set of lines connecting the points
with minimum cumulative length. In our case, the solution would lie in the
central node of the lattice being the Fermat point Brazil  _et al._ (2013) and
the Steiner node, which connects the three vertices of the high $W_{j}$
triangle, as illustrated in Fig. 4 panel A. The second case is a distribution
of $W_{j}$ that decays exponentially from the center, mimicking a more
realistic morphology for a urban monocentric structure. The two morphologies
are depicted in Fig. 4.

We find that due to nodes in the network biasing the traffic  flows, as it can
be seen in the insets  of Fig. 4 A, the traffic  flows get divided into two
types: a close range paired to a long range  set of flows, due to POI
polarization. We show  in Fig. 4 optimal solutions for  values of $\beta=0.1$,
$4.0$.  Interestingly, optimal solutions are characterized by three central
lines branching from the center (which therefore acts as Steiner node) and
connecting the three nodes with high attractiveness, therefore resembling the
solution of the Steiner tree problem. Moreover, in the case of more localized
flows ($\beta=4.0$) these  lines are also paired with large scale loops
connecting farther nodes. We also find that the heterogeneity of traffic
flows forces the appearance of a second mode in the distribution of speeds
$w_{e}$ (see Fig. 4). The two peaks in the optimal $P(w_{e})$ can be
interpreted as two different levels of speed, which suggests that the entire
process can be decomposed in two distinct mechanisms which can be mapped as a
bi-layer network: one layer at high capacity with long-range/commuting
trajectories and the other one at low velocity with short-range paths. These
two layers can be ideally separated,  hinting towards a possible extension of
the model to multilayer networks.

##  V Greater London Area: generative model for the subway system

We extend in this section the application of the model by integrating data
from a real urban structure. Specifically, we model the morphology of  Greater
London Area (GLA) on top of our framework and apply the efficiency
optimization process with the aim of understanding if the temporal efficiency
optimization of the spatial substrate paired with realistic flows is
sufficient to yield a transportation network with similar topological features
(such as a central core paired with peripheral branches Roth  _et al._ (2012))
as the London subway system. To extend the model to real urban scenarios, we
first obtain the distribution of amenities Hidalgo  _et al._ (2020) from
OpenStreetMap OpenStreetMap contributors (2017) and we use this density of
points in space as a proxy to estimate the attractiveness $W_{j}$ of a tile.
Census data for Greater London Area yards from 2014 is used to recover
population density $P_{i}$. These densities are then mapped to Uber’s H3
tiling to recover the spatial discretization in hexagonal tiles, such that we
can have direct mapping to the nodes on a triangular lattice, as in the
examples discussed in previous sections. We thus have the ingredients to
finally simulate the spatial interaction  flows $T_{ij}$ in Eq. 5. In  Fig. 5
the integration of urban data  describing the London’s morphology  in the
model is explained and we provide a depiction of the OD  matrix that arises
from the spatial interaction model. With the aim of reproducing real features,
we impose an upper limit on edge weight, so that the distribution of weights
is bounded during the optimization process: $w_{e}\in(0,w^{*})$. This better
simulates the upper bound in speed of real multilayer systems. Further
explanation of data recovery and integration in the model is provided in  SM
Section 3. We find ( see SM Figure 6) that $\\{w_{e}\\}$ distribution
displays a bi-modal shape, and this allows the analysis of the generated
network in a sub-graph defined by the set of high speed edges. In Fig. 5,
panel C, we show a sample result for $\beta=0.35$ of this sub-graph. The
characterization of the network into a central core paired with peripheral
branches as the optimal state can be visually observed.

The model’s subgraph of high speed edges is compared to the real tube network
in the  Greater London Area Gallotti and Barthelemy (2015) to assess the
similarities between the optimal structure and the real subway system. We
quantify this similarity by means of spatial scaling laws Roth  _et al._
(2012), these are convenient to highlight the recovery of the central core
structure characterized by loops, paired with quasi mono-dimensional lines
branching from the core. We investigate the distribution of nodes stations
using  the profile function $N(r)$ that quantifies the total number of
stations at a distance $r$ from the network barycenter,  computed as the
average location of all station nodes Roth  _et al._ (2012). Results of this
scaling  analysis for the real and simulated networks are presented in Fig. 6.
The two scaling regimes indicate the separation of core and branches: the
scaling of $r^{2}$ in the core center and a second trend due to mono-
dimensional branches for $r>r_{C}$, where $r_{C}$ is the radius of the core
structure.  The second  trend can be computed analytically via an integral
curve for $N(r>r_{C})$ which can be approximated by a power law $r^{\gamma}$
($\gamma=1.25\pm 0.02$, see SM Section 4), as in Ref. Roth  _et al._ (2012).
The curve of $N(r)$ is consistent with the real network and confirms scaling
laws prediction from Roth  _et al._ (2012).

Figure 6:  Scaling properties of GLA Tube stations: Profile of the number of
stations (nodes in the optimal discretized network (see  SM Section 3)
reproducing GLA underground) versus the distance from the barycenter. The
scaling of $N(r)$ profile of the model is compared with the real network
system. Scaling properties predicted in Roth  _et al._ (2012) are verified,
finding the two different scaling regimes separated at $\frac{r}{r_{C}}\sim 1$
for core paired with branches systems, where $r_{C}$ is the core radius
(characterized by $r^{2}$ scaling), and $N_{C}$ is the number of stations in
the core.  The scaling exponent $\gamma=1.25\pm 0.02$ is obtained as a linear
fit of the integral curve Roth  _et al._ (2012) for $r>r_{C}$ (see SM Section
4 for more details).

##  VI Discussion

Starting from simple conditions on temporal efficiency on a spatial network
substrate, we show that network optimization paired with traffic-like  flows
weighting the importance of specific connections in space can reproduce
complex networks features from man-made transportation networks. Specifically,
we devise a framework for spatial networks where nodes can encode features of
urban systems and can ultimately lead to the study of optimal topologies in
real scenarios. A key novelty lies in the optimization process happening on a
spatial substrate, such that edges of the resulting optimal network are
optimized to improve the efficiency of the shared space by all nodes in the
network. We show how the probabilities of moving from one point to another in
space force a transition between a tree-like and a lattice-like topology in
the optimal network. Fixing certain target points in space with a higher
attractiveness for  flows can reproduce theoretical results such as the
Steiner tree solution or leaves venation patterns. We also show that extending
these probabilities using urban spatial information and traffic-like flows
modeling forces the emergence of shared preferential paths that are organized
as complex topologies, resulting from traffic weighted optimization of network
time efficiency, which ultimately exhibits the characteristics seen in real
systems. We recover features such as a bi-modality in the speed distribution
of the edges of the network, characteristic of multilayer transportation. Or
the appearance of a central core with loops coupled to branches typical of
underground systems, as in the case of the London underground system. We find
that branches paired to large loops structures appear as optimal structures
when the network is optimized for an interplay of traffic  flows mixed between
small range travels and longer range ones typical of commuting. This novel
framework for the optimization of spatial networks in urban contexts may show
further improvements and extensions to better accommodate the concepts of
multi-layer and shared space. It could be extended also to the case of inter-
cities transportation, where  specific nodes  in the network substrate
represent cities. To conclude, in this work the problem is addressed in a
theoretical way with the aim of reproducing and understanding some features
observed in real spatial networks, but future works can exploit this framework
as a basis to understand how to generate optimal transportation networks in a
urban planning scenario.

### Competing Financial Interests

The authors declare no competing financial interests

### Data Availability

The data used in this work are publicly available from the original references

### Code Availability

The code to perform the analysis will be available upon request.

## References

  * Batty (2008) M. Batty, Science 319, 769 (2008).
  * Barthelemy (2019) M. Barthelemy, Nature Reviews Physics 1, 406 (2019).
  * Bettencourt and West (2010) L. Bettencourt and G. West, Nature 467, 912 (2010).
  * Pan  _et al._ (2013) W. Pan, G. Ghoshal, C. Krumme, M. Cebrian, and A. Pentland, Nature Communications 4 (2013), 10.1038/ncomms2961.
  * Arcaute  _et al._ (2015) E. Arcaute, E. Hatna, P. Ferguson, H. Youn, A. Johansson, and M. Batty, Journal of The Royal Society Interface 12, 20140745 (2015).
  * Bettencourt (2013) L. M. A. Bettencourt, Science 340, 1438 (2013).
  * Bassolas  _et al._ (2019) A. Bassolas, H. Barbosa-Filho, B. Dickinson, X. Dotiwalla, P. Eastham, R. Gallotti, G. Ghoshal, B. Gipson, S. A. Hazarie, H. Kautz, O. Kucuktunc, A. Lieber, A. Sadilek, and J. J. Ramasco, Nature Communications 10 (2019), 10.1038/s41467-019-12809-y.
  * Bettencourt (2020) L. M. A. Bettencourt, Science Advances 6 (2020), 10.1126/sciadv.aat8812.
  * Schläpfer  _et al._ (2021) M. Schläpfer, L. Dong, K. O’Keeffe, P. Santi, M. Szell, H. Salat, S. Anklesaria, M. Vazifeh, C. Ratti, and G. B. West, Nature 593, 522 (2021).
  * Barbosa  _et al._ (2018) H. Barbosa, M. Barthelemy, G. Ghoshal, C. R. James, M. Lenormand, T. Louail, R. Menezes, J. J. Ramasco, F. Simini, and M. Tomasini, Physics Reports 734, 1 (2018).
  * Alessandretti  _et al._ (2022) L. Alessandretti, L. G. N. Orozco, M. Saberi, M. Szell, and F. Battiston, Environment and Planning B: Urban Analytics and City Science , 239980832211081 (2022).
  * Lee  _et al._ (2017) M. Lee, H. Barbosa, H. Youn, P. Holme, and G. Ghoshal, Nature Communications 8 (2017), 10.1038/s41467-017-02374-7.
  * Gallotti  _et al._ (2021) R. Gallotti, P. Sacco, and M. D. Domenico, Complexity 2021, 1 (2021).
  * Barthélemy (2011) M. Barthélemy, Physics Reports 499, 1 (2011).
  * Morris and Barthelemy (2012) R. G. Morris and M. Barthelemy, Physical Review Letters 109 (2012), 10.1103/physrevlett.109.128703.
  * Gastner and Newman (2006) M. T. Gastner and M. E. J. Newman, Physical Review E 74 (2006), 10.1103/physreve.74.016117.
  * Barthélemy and Flammini (2006) M. Barthélemy and A. Flammini, Journal of Statistical Mechanics: Theory and Experiment 2006, L07002 (2006).
  * Louf  _et al._ (2013) R. Louf, P. Jensen, and M. Barthelemy, Proceedings of the National Academy of Sciences 110, 8824 (2013).
  * Banavar  _et al._ (1999) J. R. Banavar, A. Maritan, and A. Rinaldo, Nature 399, 130 (1999).
  * Pei  _et al._ (2022) A. Pei, F. Xiao, S. Yu, and L. Li, Scientific Reports 12 (2022), 10.1038/s41598-022-12053-3.
  * Gallotti and Barthelemy (2015) R. Gallotti and M. Barthelemy, Scientific Data 2 (2015), 10.1038/sdata.2014.56.
  * Roth  _et al._ (2012) C. Roth, S. M. Kang, M. Batty, and M. Barthelemy, Journal of The Royal Society Interface 9, 2540 (2012).
  * Tero  _et al._ (2010) A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki, Science 327, 439 (2010).
  * Louf and Barthelemy (2014) R. Louf and M. Barthelemy, Scientific Reports 4 (2014), 10.1038/srep05561.
  * Ibrahim  _et al._ (2021) A. A. Ibrahim, A. Lonardi, and C. D. Bacco, Algorithms 14, 189 (2021).
  * Zhang  _et al._ (2015) X. Zhang, A. Adamatzky, F. T. Chan, Y. Deng, H. Yang, X.-S. Yang, M.-A. I. Tsompanas, G. C. Sirakoulis, and S. Mahadevan, Scientific Reports 5 (2015), 10.1038/srep10794.
  * Szell  _et al._ (2022) M. Szell, S. Mimar, T. Perlman, G. Ghoshal, and R. Sinatra, Scientific Reports 12 (2022), 10.1038/s41598-022-10783-y.
  * Birch  _et al._ (2007) C. P. Birch, S. P. Oom, and J. A. Beecham, Ecological Modelling 206, 347 (2007).
  * Gallotti  _et al._ (2016) R. Gallotti, A. Bazzani, S. Rambaldi, and M. Barthelemy, Nature Communications 7 (2016), 10.1038/ncomms12600.
  * Ren  _et al._ (2014) Y. Ren, M. Ercsey-Ravasz, P. Wang, M. C. González, and Z. Toroczkai, Nature Communications 5 (2014), 10.1038/ncomms6347.
  * Wilson (1975) A. Wilson, Transportation Research 9, 167 (1975).
  * Viana  _et al._ (2013) M. P. Viana, E. Strano, P. Bordin, and M. Barthelemy, Scientific Reports 3 (2013), 10.1038/srep03495.
  * Newman (2010) M. Newman,  _Networks_ (Oxford University Press, 2010).
  * Dreyfus and Wagner (1971) S. E. Dreyfus and R. A. Wagner, Networks 1, 195 (1971).
  * Brazil  _et al._ (2013) M. Brazil, R. L. Graham, D. A. Thomas, and M. Zachariasen, Archive for History of Exact Sciences 68, 327 (2013).
  * Kavitha  _et al._ (2007) T. Kavitha, K. Mehlhorn, D. Michail, and K. E. Paluch, Algorithmica 52, 333 (2007).
  * Katifori  _et al._ (2010) E. Katifori, G. J. Szöllősi, and M. O. Magnasco, Physical Review Letters 104 (2010), 10.1103/physrevlett.104.048704.
  * Piovani  _et al._ (2017) D. Piovani, C. Molinero, and A. Wilson, PLOS ONE 12, e0185787 (2017).
  * Hidalgo  _et al._ (2020) C. A. Hidalgo, E. Castañer, and A. Sevtsuk, Habitat International 106, 102205 (2020).
  * OpenStreetMap contributors (2017) OpenStreetMap contributors, “Planet dump retrieved from https://planet.osm.org ,” https://www.openstreetmap.org (2017).

## Supplementary Material for “Emergence of complex networks topologies from
flow-weighted optimization of network efficiency”

##  S1  Application On Simple Network Substrates

In this Section we report more detailed results for the optimization of simple
network substrates that was discussed in the Main text, specifically for the
triangular (HEX) lattice (planar dual of the hexagonal tiling of space) and
Erdős-Rényi (ER) network. The aim is to explore the resulting topologies that
emerge both in spatial and non spatial networks when simple probabilities are
taken into consideration (see Main text).

###  S1.1 HEX Lattice and ER Network

Results of model application on the triangular (HEX) lattice and results of
model application on an Erdős-Rényi (ER) non-spatial network. For the non-
spatial generative model, an ensemble of 20 networks with N = 30 nodes and
edge probability $\rho=0.2$ is generated. For the spatial case, the
optimization process is repeated 20 times for each value of $\beta$ on the
triangular lattice with N = 37 nodes.

SM Fig. 1 and SM Fig. 2 show boxplots for the distributions of metrics
computed on the cycle basis in panels A on both figures. In panels B, samples
of the optimized network states are shown for different values of $\beta$.

Supplementary Figure S1:  Cycle basis properties and samples - HEX lattice.
Boxplot statistics for the cycle basis dimension and optimal network samples
for triangular lattice across different $\beta$ values.  Supplementary Figure
S2:  Cycle basis properties and samples - ER network. Boxplot statistics for
the cycle basis dimension and optimal network samples for an Erdos-Renyi (ER)
network across different $\beta$ values.

##  S2 Leaves patterns

Here we show an application of the model to reproduce leaf venation patterns
Katifori  _et al._ (2010). A single attracting node (single target or sink) is
considered at one of the perimeter nodes of the lattice, and the substrate is
optimized using all the other nodes as sources.

Supplementary Figure S3:  Optimal networks resembling leaves’ veins patterns.
Optimal state when a single target (orange node) in a single spatial extremity
is considered. Efficiency is optimized for all nodes in space to reach the
target. The resulting optimal state resembles tree-like patterns found in
leaves (A), while the distribution of edges weights is shown in (B).

##  S3 Application on London Tube network

In this section we discuss more in detail the application of the model on the
Greater London Area urban morphology. We show that the London subway network
spatial properties can be recovered by means of the optimal network state.
In particular, to simulate a real transportation system, an upper bound on
edges travel velocity is imposed.

###  S3.1 Data Integration: Census and OpenStreetMap Data

To extend the model to real urban scenarios, we gather data regarding the
urban morphological structure from OpenStreetMap (OSM) OpenStreetMap
contributors (2017) and Census.

####  S3.1.1 OSM Data

To model the attractiveness of nodes in the lattice, we use the density of
Points Of Interest in the urban space. We use amenities Hidalgo  _et al._
(2020) points in OSM as a proxy for Points Of Interest (POIs), and a node $j$
attractiveness ($W_{j}$) therefore encodes the density of amenities in space.

The bounding box for Greater London Area (GLA) is obtained via OSM
(https://wiki.openstreetmap.org/wiki/Bounding_Box) and POI densities are
recovered inside this box. Specifically we use the following amenities sub-
categories:

’cafe’,’college’,’library’,’school’,’university’,’kindergarten’,’restaurant’,’pub’,
’fast
food’,’bar’,’bank’,’dentist’,’pharmacy’,’hospital’,’clinic’,’doctor’,’arts
centre’,’cinema’,’community centre’,’police’,’post office’,’marketplace’

In Fig. S4 the amenities points recovered from OSM are plotted in the bounding
box.

Supplementary Figure S4:  OSM amenities query for GLA bounding box. Amenities
retrieved from OSM (https://openstreetmap.org/copyright) for the bounding box
defined by the following longitude $\left[-0.6,0.35\right]$ and latitude
interval $\left[51.25,51.75\right]$.

####  S3.1.2 Census Data

Supplementary Figure S5:  Census data retrieval. A) London’s wards data from
Census 2014 (https://data.london.gov.uk/dataset/ward-profiles-and-atlas). B)
Points are generated in space following wards polygons with a density
proportional to Census data. C) Points are mapped to HEX3 tiles with
associated densities of points. Finally, a restriction to a disc is used to
enforce symmetry in the distribution of nodes and to ease the computational
load.

Points are mapped to density in space by using HEX3 tiling
(https://eng.uber.com/h3/) with spatial resolution $RES$ = 8. This spatial
discretization process via tiling is of particular relevance to map this
information to the HEX triangular lattice model which was described in the
Main text. Both Fig. S5 and Fig. S4 show how model information $P_{i}$ and
$W_{j}$ are recovered from data and then mapped to tiles as in shown in panel
C in Fig. S5. Tiles covering the urban area are then mapped to the spatial
network model as nodes in the triangular lattice (see Fig. 5 in Main text).

For computational simplicity and to preserve the isotropy of the lattice
substrate from its central point, we restrict our analysis to a disc centered
in the London’s region with highest attractiveness, that lies approximately in
the City Of London district. Results are robust against integration of the
remaining GLA region. As the optimal edges weights are influenced by the
traffic on the substrate, the discarded regions in the Greater London Area do
not add relevant contributions when compared to more central regions with
higher population density.

###  S3.2 Results with distributions for traffic and bounded weights
distribution

To simulate a real system we limit upper edges weights to a fixed value
$w^{*}=7*w_{init}$ where $w_{init}$ is the initial edges weight assigned to
the network state ($w_{init}=1.0$) before optimization process, such that
$\sum_{e\in G}d_{e}w_{e,init}=C$. In this simulation we work with
$\beta=0.35$.

Supplementary Figure S6:  Mapping the data to Hex lattice + GLA optimized
network. A) Urban morphology data mapped to the HEX tiling and then mapped to
the triangular (HEX) lattice. At this step, OD matrix is generated and the
network is then optimized. B) Resulting optimal network state with its
distribution of edges weights showing a bi-modal shape. Only fast edges having
weight larger than a threshold ($w_{e}>5$) are kept to isolate the sub-graph
constituted by a high velocity set of edges, such as a subway system.
Supplementary Figure S7:  Pruned network vs real tube. A) London subway
network, restricted to the disc under study. B) GLA model output, limited to
the nodes and edges which have a weight $w_{e}>5.0$, where this value was
chosen from the previously plotted $\\{w_{e}\\}$ distribution as a threshold
to separate the second “mode” with fast edges. This “discretized” network is
the model on which statistical measures are performed.

##  S4 Scaling of network stations

The spatial organization of the transportation network can be inspected by
taking into consideration the number of nodes’ stations $N(r)$ up to a
distance $r$ from the barycenter of stations. In the Main text we analyze the
scaling regimes for the simulated and real London Tube network Gallotti and
Barthelemy (2015) following the analysis employed in Roth  _et al._ (2012).
Specifically, they obtain functional forms for the scaling properties of
$N(r)$ for different distances from the barycenter. We report here the
scalings obtained in Roth  _et al._ (2012):

$N(r)\sim\begin{cases}\rho_{C}\pi r^{2}&\text{ for }r<r_{C}\\\ \rho_{C}\pi
r_{C}^{2}+\mathcal{N}_{B}\int_{r_{C}}^{r}\frac{dr}{\Delta(r)}&\text{ for
}\quad r_{C}<r<r_{\max}\\\ N&\text{ for }r>r_{\max}\end{cases}$ (S1)

Specifically, they show that in the large distance regime ($r>r_{C}$ and
$r<r_{\max}$ ) the number of stations can be approximated by adding the
integral curve $\mathcal{N}_{B}\int_{r_{C}}^{r}\frac{dr}{\Delta(r)}$ to a
constant term. In Roth  _et al._ (2012) it is also reported that the large
distance behavior can be also, in general, approximated by a scaling law.
Therefore we plot here the computation of the integral curve against rescaled
values of $r$, and show that in that regime it can be approximated by a power
law, and its exponent can be obtained via a linear fit. The $N(r)$ curve in
the Main text was computed on the real London Tube restricted to the disc area
as presented in SM Fig. 7. The associated exponent $\gamma=1.25\pm 0.02$ was
used in the Main text to highlight the secondary scaling for $r>r_{C}$. As
discussed in Ref. Roth  _et al._ (2012), due to $\Delta(r)$ in SM Eq. S1 being
often noisy, this scaling property is not often well reproduced in empirical
networks, and it is often restricted to a small region of $r$ values. In Fig.
6 of the Main text, we see that the secondary scaling with the exponent
$\gamma$ approximates $N(r)$ in a limited interval for $r>r_{C}$. The value of
the exponent in the branches region is expected to be $\gamma<2.0$ Roth  _et
al._ (2012) and our result is consistent with this.

Supplementary Figure S8:  Scaling of the integral curve. Log-Log plot of the
integral curve and its linear fit. With the associated value of the $\gamma$
exponent approximated, to highlight a secondary scaling region as mentioned in
Roth  _et al._ (2012).

## References

  * Batty (2008) M. Batty, Science 319, 769 (2008).
  * Barthelemy (2019) M. Barthelemy, Nature Reviews Physics 1, 406 (2019).
  * Bettencourt and West (2010) L. Bettencourt and G. West, Nature 467, 912 (2010).
  * Pan  _et al._ (2013) W. Pan, G. Ghoshal, C. Krumme, M. Cebrian, and A. Pentland, Nature Communications 4 (2013), 10.1038/ncomms2961.
  * Arcaute  _et al._ (2015) E. Arcaute, E. Hatna, P. Ferguson, H. Youn, A. Johansson, and M. Batty, Journal of The Royal Society Interface 12, 20140745 (2015).
  * Bettencourt (2013) L. M. A. Bettencourt, Science 340, 1438 (2013).
  * Bassolas  _et al._ (2019) A. Bassolas, H. Barbosa-Filho, B. Dickinson, X. Dotiwalla, P. Eastham, R. Gallotti, G. Ghoshal, B. Gipson, S. A. Hazarie, H. Kautz, O. Kucuktunc, A. Lieber, A. Sadilek, and J. J. Ramasco, Nature Communications 10 (2019), 10.1038/s41467-019-12809-y.
  * Bettencourt (2020) L. M. A. Bettencourt, Science Advances 6 (2020), 10.1126/sciadv.aat8812.
  * Schläpfer  _et al._ (2021) M. Schläpfer, L. Dong, K. O’Keeffe, P. Santi, M. Szell, H. Salat, S. Anklesaria, M. Vazifeh, C. Ratti, and G. B. West, Nature 593, 522 (2021).
  * Barbosa  _et al._ (2018) H. Barbosa, M. Barthelemy, G. Ghoshal, C. R. James, M. Lenormand, T. Louail, R. Menezes, J. J. Ramasco, F. Simini, and M. Tomasini, Physics Reports 734, 1 (2018).
  * Alessandretti  _et al._ (2022) L. Alessandretti, L. G. N. Orozco, M. Saberi, M. Szell, and F. Battiston, Environment and Planning B: Urban Analytics and City Science , 239980832211081 (2022).
  * Lee  _et al._ (2017) M. Lee, H. Barbosa, H. Youn, P. Holme, and G. Ghoshal, Nature Communications 8 (2017), 10.1038/s41467-017-02374-7.
  * Gallotti  _et al._ (2021) R. Gallotti, P. Sacco, and M. D. Domenico, Complexity 2021, 1 (2021).
  * Barthélemy (2011) M. Barthélemy, Physics Reports 499, 1 (2011).
  * Morris and Barthelemy (2012) R. G. Morris and M. Barthelemy, Physical Review Letters 109 (2012), 10.1103/physrevlett.109.128703.
  * Gastner and Newman (2006) M. T. Gastner and M. E. J. Newman, Physical Review E 74 (2006), 10.1103/physreve.74.016117.
  * Barthélemy and Flammini (2006) M. Barthélemy and A. Flammini, Journal of Statistical Mechanics: Theory and Experiment 2006, L07002 (2006).
  * Louf  _et al._ (2013) R. Louf, P. Jensen, and M. Barthelemy, Proceedings of the National Academy of Sciences 110, 8824 (2013).
  * Banavar  _et al._ (1999) J. R. Banavar, A. Maritan, and A. Rinaldo, Nature 399, 130 (1999).
  * Pei  _et al._ (2022) A. Pei, F. Xiao, S. Yu, and L. Li, Scientific Reports 12 (2022), 10.1038/s41598-022-12053-3.
  * Gallotti and Barthelemy (2015) R. Gallotti and M. Barthelemy, Scientific Data 2 (2015), 10.1038/sdata.2014.56.
  * Roth  _et al._ (2012) C. Roth, S. M. Kang, M. Batty, and M. Barthelemy, Journal of The Royal Society Interface 9, 2540 (2012).
  * Tero  _et al._ (2010) A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki, Science 327, 439 (2010).
  * Louf and Barthelemy (2014) R. Louf and M. Barthelemy, Scientific Reports 4 (2014), 10.1038/srep05561.
  * Ibrahim  _et al._ (2021) A. A. Ibrahim, A. Lonardi, and C. D. Bacco, Algorithms 14, 189 (2021).
  * Zhang  _et al._ (2015) X. Zhang, A. Adamatzky, F. T. Chan, Y. Deng, H. Yang, X.-S. Yang, M.-A. I. Tsompanas, G. C. Sirakoulis, and S. Mahadevan, Scientific Reports 5 (2015), 10.1038/srep10794.
  * Szell  _et al._ (2022) M. Szell, S. Mimar, T. Perlman, G. Ghoshal, and R. Sinatra, Scientific Reports 12 (2022), 10.1038/s41598-022-10783-y.
  * Birch  _et al._ (2007) C. P. Birch, S. P. Oom, and J. A. Beecham, Ecological Modelling 206, 347 (2007).
  * Gallotti  _et al._ (2016) R. Gallotti, A. Bazzani, S. Rambaldi, and M. Barthelemy, Nature Communications 7 (2016), 10.1038/ncomms12600.
  * Ren  _et al._ (2014) Y. Ren, M. Ercsey-Ravasz, P. Wang, M. C. González, and Z. Toroczkai, Nature Communications 5 (2014), 10.1038/ncomms6347.
  * Wilson (1975) A. Wilson, Transportation Research 9, 167 (1975).
  * Viana  _et al._ (2013) M. P. Viana, E. Strano, P. Bordin, and M. Barthelemy, Scientific Reports 3 (2013), 10.1038/srep03495.
  * Newman (2010) M. Newman,  _Networks_ (Oxford University Press, 2010).
  * Dreyfus and Wagner (1971) S. E. Dreyfus and R. A. Wagner, Networks 1, 195 (1971).
  * Brazil  _et al._ (2013) M. Brazil, R. L. Graham, D. A. Thomas, and M. Zachariasen, Archive for History of Exact Sciences 68, 327 (2013).
  * Kavitha  _et al._ (2007) T. Kavitha, K. Mehlhorn, D. Michail, and K. E. Paluch, Algorithmica 52, 333 (2007).
  * Katifori  _et al._ (2010) E. Katifori, G. J. Szöllősi, and M. O. Magnasco, Physical Review Letters 104 (2010), 10.1103/physrevlett.104.048704.
  * Piovani  _et al._ (2017) D. Piovani, C. Molinero, and A. Wilson, PLOS ONE 12, e0185787 (2017).
  * Hidalgo  _et al._ (2020) C. A. Hidalgo, E. Castañer, and A. Sevtsuk, Habitat International 106, 102205 (2020).
  * OpenStreetMap contributors (2017) OpenStreetMap contributors, “Planet dump retrieved from https://planet.osm.org ,” https://www.openstreetmap.org (2017).


# Time Series Analysis of Urban Liveability

Alex Levering
Wageningen University
Wageningen, the Netherlands
<EMAIL_ADDRESS>
Diego Marcos
INRIA
Montpellier, France
Devis Tuia
EPFL
Sion, Switzerland

###### Abstract

In this paper we explore deep learning models to monitor longitudinal
liveability changes in Dutch cities at the neighbourhood level. Our
liveability reference data is defined by a country-wise yearly survey based on
a set of indicators combined into a liveability score, the Leefbaarometer. We
pair this reference data with yearly-available high-resolution aerial images,
which creates yearly timesteps at which liveability can be monitored. We
deploy a convolutional neural network trained on an aerial image from 2016 and
the Leefbaarometer score to predict liveability at new timesteps 2012 and
2020. The results in a city used for training (Amsterdam) and one never seen
during training (Eindhoven) show some trends which are difficult to interpret,
especially in light of the differences in image acquisitions at the different
time steps. This demonstrates the complexity of liveability monitoring across
time periods and the necessity for more sophisticated methods compensating for
changes unrelated to liveability dynamics.

##  1 Introduction

Cities are dynamic and complex entities with many different and competing
functions. One such function is the provision of living space to its
residents, which is commonly referred to as the liveability of the city. It is
commonly defined as "the degree to which its provisions and requirements fit
with the needs and capacities of its members" [1]. It therefore promotes the
human-centric design of cities such that they are made into places where
people want to live. A failure to provide liveable environments can result in
a variety of detrimental health effects, such as higher mortality rates [2]
and higher morbidity rates [3]. As such, there is a vested interest from
policymakers to provide liveable environments for residents.

Measuring liveability is generally a costly and time-consuming process as it
is done through the surveying of residents. Recent research has attempted to
predict liveability at a larger scale by means of remote sensing imagery [4,
5]. Such models are now delivering promising results which potentially allow
for monitoring without requiring the acquisition of expensive reference data.
The next step in operationalizing liveability monitoring is to explore the
potential for such models to predict over multiple timesteps. If such models
can reliably forecast liveability, then data-constrained countries may
consider the adoption of remote sensing for monitoring the liveability of the
territory.

In this paper we discuss the potential for deep learning on aerial images to
monitor the liveability of cities across multiple time steps. We discuss the
creation of a liveability monitoring dataset over The Netherlands and test the
applicability of pre-trained models for liveability estimation to demonstrate
the need for more sophisticated time-series modelling approaches accounting
for domain shifts related to different acquisition conditions, camera, etc. We
study multi-temporal liveability estimation in two Dutch cities and use the
Leefbaarometer, a Dutch project providing longitudinal studies for liveability
in the Netherlands, as reference data to evaluate our strategy.

##  2 Liveability Monitoring data

Assessing the liveability of cities has been an actively studied topic over
the years, with more recent literature focusing on predictive modeling. Recent
work has made use of advances in deep learning models to predict liveability
factors directly from overhead imagery, proving that liveability factors can
be accurately predicted with remote sensing [4, 6]. However, liveability
monitoring requires multiple timesteps of both the reference data for
neighbourhoods and the aerial images, and for this reason there have thus far
been limited research efforts which attempt to monitor liveability dynamics
(i.e. longitudinally over time) from overhead images. In the following, we
detail our two sources of data, which we used for monitoring at three time
steps: 2012, 2016 (data used to train the model) and 2020.

###  2.1 Reference data: the Leefbaarometer.

The Netherlands is a notable exception to the lack of reference data. The
Dutch government has monitored the liveability of residents for many years
through a project known as the _Leefbaarometer_ (LBM). In the LBM project,
researchers attempt to model the opinion of residents on their liveability as
a function of various publicly available datasources which can be grouped as
_liveability domains_. These domains encompass factors which affect
liveability, such as socio-economic variables, housing quality, and access to
amenities. The LBM project therefore measures both how liveable a
neighbourhood is, as well as how various domains contribute to that
liveability. The 3.0 version of the LBM project released in 2022 reaches an
$R^{2}$ of 0.75 in the reproduction of the residents responses with the five
selected domains (center column of Fig. 1). Furthermore, the LBM is verified
to be accurate by policy makers. In this version, the LBM model is calculated
yearly from 2012 until 2020, and is due to be updated for the recent and
coming years. In this work we use the liveability scores provided by LBM for
each 100 m $\times$ 100 m inhabited cell across the Netherlands.

Figure 1: Data construction process for the liveability time series monitoring
task. For each year there is a grid with liveability reference data available,
which can be used to match the yearly aerial images to the neighbourhood
liveability data. Each grid cell also has liveability domain score information
available, which is the contribution of each domain to the liveability of a
neighbourhood. As it is zero-centered, the overall score of 0.23 indicates
that this patch is above averagely liveable.

###  2.2 Aerial data.

We use very high resolution aerial imagery to predict the liveability scores.
Since 2016, the government of Netherlands releases a yearly cloud-free winter
aerial image with a uniform sensor and 25cm pixel resolution across the entire
country. Further yearly images are available before this timestep as well,
with a lower resolution and fewer timesteps. The continued availability of
both data sources means that the amount of data for this monitoring task will
continue to increase over the years. Fig. 1 illustrates a set of image
patches, together with the Liveability score and their intermediate domain
scores. To encourage the use of the dataset, we provide scripts and links to
the data to prepare it in the following Github repository:
https://github.com/ahlevering/LBMOverheadMonitor.

##  3 Monitoring using Pre-trained Models

We developed a machine learning model to predict liveability, as well as the
domain scores, in [5] and we use that model as a starting point to study the
complexity of monitoring multiple timesteps. For this purpose we use a pre-
trained model to predict the liveability at new timesteps. Our model is
trained on an aerial image from 2016 (details in [5]) and we apply it without
modifications to aerial images from 2012 and 2020. For both timesteps, we
compare the overall liveability score with the ground truth provided by the
LBM 3.0 project. While our model is trained on the 2.0 version of the LBM
project, there is compatibility with the 3.0 version as they both measure
liveability as the deviation to the national average. However, to maintain
compatibility, we resize the image to 500x500 pixels, rather than 600x600. For
this paper we do not consider the intermediate sub-scores as the definition of
the sub-scores changed between versions 2.0 and 3.0. We assess neighbourhood
shifts relative to one another by analyzing their relative ranking within the
city. If one neighbourhood improves then, by comparison, another neighbourhood
may seem less attractive. This also makes the prediction case less sensitive
to mean shifts when predicting on new images. We use the reference data from
the LBM 3.0 to determine if the model picks up on measured shifts in the
liveability of neighbourhoods.

We consider two cities for our experiment chosen for their contrasting
qualities, namely Amsterdam and Eindhoven. Amsterdam is the biggest
conurbation of The Netherlands. It has a long history and it is unique in its
design compared to other cities in the country. Eindhoven is an industrial
city which urbanized only since the industrial revolution. Further, it
important to note that Amsterdam was included in the training cities in used
to build the model of [5], while Eindhoven was included in its test set and
was never seen during model training.

##  4 Results and discussion

###  4.1 Liveability rankings

In this section we discuss how our model predicted the liveability rankings in
Amsterdam and Eindhoven for two new timesteps unseen during training. Firstly
we present Kendall’s $\tau$ for both 2012 and 2020 in table 1. $\tau$ is a
coefficient calculating correlations between in two lists using ranking of
indices rather than the actual values of such indices. The sharp decrease in
performance on the new timesteps in Amsterdam suggests that model has
overfitted to this region during training. However, the minor decrease in
accuracy over Eindhoven suggests that liveability prediction models can remain
stable when predicting new timesteps, even in geographical regions not seen
during training. In Fig. 2 we show the predicted scores, the LBM 3.0 reference
values, and the changes between 2012 and 2020. We discuss the patterns for
each city separately.

Table 1: Kendall’s $\tau$ of the predictions compared to the LBM 3.0 reference scores. We also include the comparison to 2016, which is the year on which the model was trained. | 2012 $\tau$ | 2016 $\tau$ | 2020 $\tau$
---|---|---|---
Amsterdam | 0.316 | 0.656 | 0.337
Eindhoven | 0.487 | 0.537 | 0.478

Amsterdam

The predictions at different timesteps for Amsterdam show some correct change
patterns, but a closer inspection yields various issues. At timestep 2012, the
model predicts the overall trends correctly, but with some notable exceptions:
it especially mispredicts the liveability of the peripheral areas. Most
notably, the model ranks _Amstelveen_ (south-western peripheral city) among
the worst areas to live. In turn, it considers areas in the north of Amsterdam
to be more liveable than the reference data. The model’s predictions for the
2020 reference data are even more unexpected. Firstly, the model is better
able to predict the liveability of the periphery where it follows the overall
trend better. However, in turn it severely underpredicts the liveability of
the city center of Amsterdam. As a whole, the 2020 prediction shows swings in
liveability which are unrealistic given the time difference to the year in
which the original model was trained (2016). This result reveals that the
predictive qualities of liveability prediction models can break down when
introduced with new time series data. Our results corroborate previous
research which suggest that liveability factors of culturally distinct cities
are more difficult to generalize [6]. Moreover, improving predictions for the
center of Amsterdam through normalization experiments resulted in the
degradation of predictions in other cities, which indicates that there is a
domain shift problem which needs to be studied in more detail.

Figure 2: Liveability rankings as predicted by a pre-trained liveability
prediction model compared to the LBM 3.0 reference values. Neighbourhoods are
ranked from lowest (red) to highest (blue) ranking. We also show the changes
between 2012 and 2020 from purple (negative change) to green (positive change)
for both cities.

Eindhoven

For the city of Eindhoven the results are less erratic over the years.
Eindhoven is part of the test set in our original study [5], so it is
surprising that the model has more stable ranking predictions here compared to
Amsterdam. For both timesteps, the model predicts fewer strong swings in
liveability, and generally predicts changes in the same areas as the reference
data. However, the changes it predicts in some places trend in the wrong
direction. For instance, near _Veldhoven_ , which is the bottom left block of
grid cells, it correctly predicts that the north of the village has a positive
development trend. However, the model does not predict an overall improvement
in the entire area, as it considers the south of the village to be less
liveable in 2020. It is an encouraging result that the model is able to
provide stable ranking predictions which fluctuate in the same areas as the
reference data, but more research is needed to provide models which can more
accurately correctly predict the direction of these changing trends.

###  4.2 Liveability monitoring

With the release of the third version of the LBM and the availability of
aerial images both at yearly resolution (2012 onwards), it is now possible to
approach this problem as a time series prediction problem. The results above
show some pointers to inspire future research efforts in liveability modeling
with deep learning. First, our results suggest that such a task is not
trivial, and will likely require a more sophisticated approach. Approaches
rooted in deep domain adaptation [7] seem to be necessary, for example to
reduce the gap in the image space between the different years acquisitions for
example proposing image normalisation strategies that would align the years:
fine-tuning on the specific year data, projective methods [8, 9] or
adversarial approaches [10] are potential options to be explored, to be sure
that the discrepancies between monitoring years are due to changes in
liveability, rather than differences in acquisition conditions of the images.

##  5 Conclusions

In this paper we have studied the task of liveability monitoring from aerial
imagery over the entirety of The Netherlands and focused on the examples of
the cities of Amsterdam and Eindhoven. We created and made available a dataset
from the yearly liveability reference data and aerial images, which opens up
the potential for long-term monitoring of liveability. We have also tested the
applicability for pre-trained liveability prediction models to predict over
new timesteps, which has shown promising results for the general trends, but
also highlighted the necessity for more advanced methods for this task able to
account for domain differences and provide accurate liveability maps out of
the training domain.

## References

  * [1] Ruut Veenhoven, Joop Ehrhardt, Monica Sie Dhian Ho, and Astrid de Vries.  Happiness in nations: Subjective appreciation of life in 56 nations 1946–1992.  Happiness in nations: Subjective appreciation of life in 56 nations 1946–1992. Erasmus University Rotterdam, Rotterdam, Netherlands, 1993.  Pages: 365.
  * [2] M. Haan, G. A. Kaplan, and T. Camacho.  Poverty and health. Prospective evidence from the Alameda County Study.  American Journal of Epidemiology, 125(6):989–998, June 1987.
  * [3] Sharrelle Barber, DeMarc A. Hickson, Ichiro Kawachi, S. V. Subramanian, and Felton Earls.  Neighborhood Disadvantage and Cumulative Biological Risk Among a Socioeconomically Diverse Sample of African American Adults: An Examination in the Jackson Heart Study.  Journal of Racial and Ethnic Health Disparities, 3(3):444–456, 2016\.
  * [4] Esra Suel, Samir Bhatt, Michael Brauer, Seth Flaxman, and Majid Ezzati.  Multimodal deep learning from satellite and street-level imagery for measuring income, overcrowding, and environmental deprivation in urban areas.  Remote Sensing of Environment, 257:112339, May 2021.
  * [5] Alex Levering, Diego Marcos, Jasper van Vliet, and Devis Tuia.  Predicting the liveability of Dutch cities with aerial images and semantic intermediate concepts.  Remote Sensing of Environment, 287:113454, March 2023.
  * [6] Sanja Scepanovic, Sagar Joglekar, Stephen Law, and Daniele Quercia.  Jane Jacobs in the Sky: Predicting Urban Vitality with Open Satellite Data.  ACM Human-Computer Interaction, 5:1–25, April 2021.
  * [7] Mei Wang and Weihong Deng.  Deep visual domain adaptation: A survey.  Neurocomputing, 312:135–153, October 2018.
  * [8] Wolfgang Gross, Devis Tuia, Uwe Soergel, and Wolfgang Middelmann.  Nonlinear Feature Normalization for Hyperspectral Domain Adaptation and Mitigation of Nonlinear Effects.  IEEE Transactions on Geoscience and Remote Sensing, 57(8):5975–5990, August 2019.  Conference Name: IEEE Transactions on Geoscience and Remote Sensing.
  * [9] Onur Tasar, S L Happy, Yuliya Tarabalka, and Pierre Alliez.  ColorMapGAN: Unsupervised Domain Adaptation for Semantic Segmentation Using Color Mapping Generative Adversarial Networks.  IEEE Transactions on Geoscience and Remote Sensing, PP:1–16, March 2020.
  * [10] Benjamin Kellenberger, Onur Tasar, Bharath Bhushan Damodaran, Nicolas Courty, and Devis Tuia.  Deep Domain Adaptation in Earth Observation.  In Deep Learning for the Earth Sciences, pages 90–104. John Wiley & Sons, Ltd, 2021.  Section: 7 _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119646181.ch7.