# Comparing spatial networks: A ’one size fits all’ efficiency-driven approach

Ignacio Morer<EMAIL_ADDRESS>Departament de Fisica de la Matèria
Condensada, Universitat de Barcelona, Barcelona, Spain  Universitat de
Barcelona Institute of Complex Systems (UBICS) Universitat de Barcelona,
Barcelona, Spain  Alessio Cardillo<EMAIL_ADDRESS>Institut Català
de Paleoecologia Humana i Evolució Social (IPHES), E-43007 Tarragona, Spain
Department of Engineering Mathematics, University of Bristol, Bristol, BS8
1UB, United Kingdom  Department of Computer Science and Mathematics,
Universitat Rovira i Virgili, E-43007 Tarragona, Spain  GOTHAM Lab – Institute
for Biocomputation and Physics of Complex Systems (BIFI), Universidad de
Zaragoza, E-50018 Zaragoza, Spain  Albert Díaz-Guilera<EMAIL_ADDRESS>Departament de Fisica de la Matèria Condensada, Universitat de Barcelona,
Barcelona, Spain  Universitat de Barcelona Institute of Complex Systems
(UBICS) Universitat de Barcelona, Barcelona, Spain  Luce Prignano
<EMAIL_ADDRESS>Departament de Fisica de la Matèria Condensada,
Universitat de Barcelona, Barcelona, Spain  Universitat de Barcelona Institute
of Complex Systems (UBICS) Universitat de Barcelona, Barcelona, Spain  Sergi
Lozano<EMAIL_ADDRESS>Universitat de Barcelona Institute of Complex Systems
(UBICS) Universitat de Barcelona, Barcelona, Spain  Institut Català de
Paleoecologia Humana i Evolució Social (IPHES), E-43007 Tarragona, Spain  Àrea
de Prehistòria, Universitat Rovira i Virgili, Tarragona, Spain  Departament
d’Història Econòmica, Institucions, Política i Economia Mundial, Universitat
de Barcelona, Barcelona, Spain

###### Abstract

Spatial networks are a powerful framework for studying a large variety of
systems belonging to a broad diversity of contexts: from transportation to
biology, from epidemiology to communications, and migrations, to cite a few.
Spatial networks can be described in terms of their total cost (_i.e._ , the
total amount of resources needed for building or traveling their connections).
Here, we address the issue of how to gauge and compare the quality of spatial
network designs (_i.e._ , efficiency vs. total cost) by proposing a two-step
methodology. Firstly, we assess the network’s design by introducing a quality
function based on the concept of network’s efficiency. Second, we propose an
algorithm to estimate computationally the upper bound of our quality function
for a given network. Complementarily, we provide a universal expression to
obtain an approximated upper bound to any spatial network, regardless of its
size. Smaller differences between the upper bound and the empirical value
correspond to better designs. Finally, we test the applicability of this
analytic tool-set on spatial network data-sets of different nature.

##  I Introduction

A large variety of systems, both natural and artificial, are composed by
interconnected units embedded in space. All these systems can be mapped onto
spatial networks hagget-chorley-1969 ; barthelemy-phys_rep-2011 , a powerful
framework which provides the mathematical and conceptual tools to study them
formally. Such a framework, built on basic common features, applies to a broad
diversity of contexts: from transportation guimera-pnas-2005 ; gastner-
jstat-2006 ; porta-epb-2006 ; kaluza-jrosocint-2010 ; rinaldo-pnas-2014  to
biology lui-jchemphys-2013 ; sandhu-cell_rep-2012 , from neuroscience
bullmore-nat_rev_neu-2012  to animal behavior perna-jrosocint-2014 ;
bottinelli-jrosocint-2015 ; cook-etal-2014 , from epidemiology isella-
jtheobio-2010 ; balcan-pnas-2009 ; poletto-jtheobio-2013  to communications
blondel-epjds-2015  and migrations porat-epl-2016 , to cite a few.

Spatial networks are those whose nodes have associated spatial coordinates.
Consequently, links – that is, connections between nodes – are characterized
by the distance between the pair of nodes they connect. Such a distance
usually, translates into a _cost_ , standing for the amount of resources
needed for building or traveling (or both) a given connection. For example,
the cost to build a road connecting two cities increases as its extremities
are farther from each other.

Most of the literature on the topic (gastner-jstat-2006 ; gastner-epjb-2006 ;
louf-pnas-2013 ) considers the case of edges’ costs directly proportional to
their length (_i.e._ , distances between connected nodes). Despite other
options are possible (_e.g._ , cost proportional to a monotonically increasing
function of the connection’s length), this is generally considered a good
description of almost the totality of the systems representable as spatial
networks. Notice that construction or operative costs can also be assigned to
nodes in particular circumstances (_e.g._ , routers in communication
networks). In this article, we neglect such node costs. By doing this, spatial
networks can be characterized in terms of a _total cost_ defined as the sum of
their edges’ costs.

When comparing two spatial networks with the same node layout it is reasonable
to think that the one with a higher total link length can achieve better
efficiency. Likewise, the same total link length might not be able to
communicate in an equally efficient way two sets of nodes, that are different
in number and relative location. Our goal is to enable a meaningful comparison
between two spatial networks even in such cases. To reach that goal, we assume
the total cost to be an external constraint (_i.e._ , determined by external
factors fixing the amount of resources available for building connections),
and focus on assessing to what extent such resources have been employed
profitably. Specifically, we frame the quality of a network with regards to
the total cost and establish a reference of how better could it be under that
constraint, similarly to what was done in zamora-commphys-2019  for non-
spatial networks. It is also a very similar approach to that of Cardillo _et
al._ in which real networks (urban street patterns) are compared by rescaling
both their efficiencies and costs with values obtained from suitable null-
cases cardillo-pre-2006 .

Here, we propose a two-step methodology:  (1) to assess the design of
different networks by means of a quality function;  (2) to compare those
values with an upper bound estimated computationally and therefore calculate
how far each system lies from an optimal, equally constrained network.

The paper is organized as follows. After characterizing the behavior of some
reference models of spatial networks, we introduce the concept of _integrated
efficiency_ , $E_{\text{int}}$, and set it as our quality function (Sec. II).
Then, in Sec. III we devise an algorithm to estimate the maximum value that
such a metric can take (upper bound) for a given set of node positions (node
layout) and a given total link length (constrained total cost). We then
provide an approximate universal relation that allows to determine the upper
bound of the integrated efficiency as a function of the average distance
between nodes and the total link length, for any number of nodes. In Sec. IV,
we measure the quality of several empirical systems and put them in context
with their optimal counterparts, to enable a meaningful comparison regardless
of their individual characteristics. Finally, as an example, we present some
potential applications.

##  II Global, Local, and Integrated Efficiencies

Given a spatial network $G$ with $N$ nodes, its structure is completely
determined by the adjacency and distance matrices, $\mathcal{A}$ and
$\mathcal{D}$. Their corresponding elements $\\{a_{ij}\\}$ take value $1$
$(0)$ if the connection exists (does not exist), and $\\{d_{ij}\\}$ take
finite positive values corresponding to the spatial distance between nodes $i$
and $j$ latora-book-2017 . These two matrices fully determine the shortest
paths matrix $\mathcal{L}$, whose elements $\\{l_{ij}\\}$ stand for the length
of the shortest path between nodes $i$ and $j$. $l_{ij}$ is a straight sum of
the weights of the links in the path, no matter the number of steps111It may
happen that a path including more links is shorter than another with fewer but
longer steps..

There are several ways to assess the quality of a spatial network. We chose to
use the concept of _network efficiency_ which is, essentially, a comparison
between the spatial distance and the shortest path. The efficiency in the
communication between two nodes $i$ and $j$, $E_{ij}$ corresponds to the ratio
between the spatial distance and the shortest path length,
$E_{ij}=d_{ij}/l_{ij}$, commonly known as _detour index_ or _route factor_
barthelemy-phys_rep-2011 . By computing the average over all pairs of nodes we
obtain the _global efficiency_ as proposed in vragovic-pre-2005 , (despite an
alternative definition has been given in latora-marchiori-2001 ):

$E_{\text{glob}}=\frac{1}{N\left(N-1\right)}\sum_{i\neq
j}\frac{d_{ij}}{l_{ij}}\,.$ (1)

$E_{\text{glob}}$ quantifies the ability of the system as a whole to
communicate efficiently among its elements. Additionally, it is also relevant
to assess the fault tolerance of the system’s communicability at local level.
To this aim, Latora and Marchiori latora-marchiori-2001  introduced the so-
called _local efficiency_ , $E_{\text{loc}}$. Such an indicator measures how
efficient is the communication in the local neighborhood of a node $i$ after
its removal, thus accounting for the robustness of the connectivity against
local damage.. In this paper, we adopt a modified $E_{\text{loc}}$ proposed by
Vragovic _et al._ vragovic-pre-2005  that measures the efficiency of the
communication between any two neighbors $j$ and $m$ of node $i$ – no matter
how far they are – considering all possible paths connecting them:

$E_{\text{loc}}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{k_{i}\left(k_{i}-1\right)}\sum_{j\neq
m\in\Gamma_{i}}\frac{d_{jm}}{l_{jm/i}}\,,$ (2)

where $\Gamma_{i}$ represents the local sub-graph of neighbors of node $i$ and
$l_{jm/i}$ is the length of the shortest path joining nodes $j$ and $m$ in
absence of $i$, and $l_{jm/i}=\infty$ when $j$ and $m$ belong to disconnected
components. Finally, $k_{i}$ is the degree (_i.e._ , number of connections) of
node $i$.

Given a certain layout of the nodes in space, there are multiple connectivity
patterns presenting (approximately) the same _total length_ ,
$L_{\text{tot}}=\sum_{i,j}a_{ij}\,d_{ij}$. Among the plethora of spatial
network models available in the literature gastner-jstat-2006 ; gastner-
epjb-2006 ; banavar-prl-2000 ; kaiser-pre-2004 ; barthelemy-jstat-2006 ;
barthelemy-prl-2008 , we selected a few models as benchmarks, thus
encompassing a wide spectrum of possibilities:

  * •

the Minimum Spanning Tree (MST) barthelemy-phys_rep-2011 ;

  * •

the Greedy Triangulation (GT), a maximally connected triangulation minimizing
its total length cardillo-pre-2006 ;

  * •

the Equitable Efficiency Model (EEM), a growth model adding one link at a time
ensuring that such move constitutes the best improvement in the communication
among any pair of nodes prignano-jas-2019 ;

  * •

the Gastner-Newman model (GN), a growth model where edges have an _effective
length_ which combines Euclidean and topological distances altogether gastner-
epjb-2006 .

A detailed explanation of all the models (except for the well-known Minimum
Spanning Tree) can be found in the Appendix.

Starting from a random distribution of nodes in a unit square, we built the
corresponding MST, GT, EEM, and GN networks and computed their
$E_{\text{glob}}$ and $E_{\text{loc}}$. The results are shown in Fig. 1 and
denote remarkable differences amidst the benchmarks. Preliminary work
confirmed that differences persist when one considers different spatial
distributions of points, that is, they do not depend critically on the details
of the node positions. Additionally, the analysis of spatial layouts of the
empirical networks considered in this study (see Sec. IV), showed that the
distribution of distances among the points is compatible with a random uniform
distribution.

Figure 1: Average global efficiency, $\langle E_{\text{glob}}\rangle$, (panel
a), and average local efficiency, $\langle E_{\text{loc}}\rangle$, (panel b)
as a function of the total length of the system, $L_{\text{tot}}$, for
different models: Greedy Triangulation (GT), Minimum Spanning Tree (MST),
Equitable Efficiency Model (EEM), and Gastner-Newman (GN) with $\lambda=0.85$
and $\gamma=0.06$. Dots refer to the properties of the final networks, while
the solid lines account for the evolution – whenever available – of those
properties throughout the growth process. The graphs have been built by
averaging uniform random distributions of $N=100$ points considering
$N_{real}=50$ different realizations.

In Fig. 1, points correspond to the networks grown until their total length is
almost the same as the GT one. Lines, instead, account for the intermediate
stages corresponding to the growth phase – _i.e._ , progressive addition of
edges – of the model (if available). A first glance at the panels reveals some
interesting features. The first one is the opposite behavior of EEM and GT,
with EEM performing better than GT in terms of $E_{\text{glob}}$, and the
other way around for $E_{\text{loc}}$. Another feature is the fact that MST
has a nonzero value of $E_{\text{glob}}$ but $E_{\text{loc}}=0$, although this
is expected since the removal of a single node in a tree implies the
impossibility of communicating between its neighbors. Finally, the position of
GN networks indicates higher values of both efficiencies but at a higher total
cost (_i.e._ , total link length). This is due to the fact that GN networks
are built using a cost function different from the mere spatial one.

The analysis of Fig. 1 highlights how differently the benchmark models behave
with respect to the two efficiencies. Such differences can be leveraged and
used to characterize each network using both $E_{\text{glob}}$ and
$E_{\text{loc}}$. Thus, we can represent each network with the pair of values
$(E_{\text{loc}},E_{\text{glob}})$, which corresponds to a point in a two
dimensional $[0,1]\times[0,1]$ diagram.

Any topology lies inside this square, its position depending on its specific
features. For instance, a set of isolated nodes (networks with no links) will
lie at the lower left corner $(0,0)$, while the upper right one $(1,1)$
corresponds to the complete graph which, by definition, has the maximum
possible efficiency. Topologies having $(E_{\text{loc}},0)$ or
$(E_{\text{loc}},1)\;\forall E_{\text{loc}}\in\;]0,1[$ are not allowed since
$E_{\text{glob}}=0$ and $E_{\text{glob}}=1$ can be obtained exclusively by a
set of isolated nodes or a complete graph, respectively. Topologies falling on
the $E_{\text{loc}}=0$ line correspond to tree-like (acyclic) graphs, while
those falling on the $E_{\text{loc}}=1$ line are ensembles of disconnected
complete subgraphs.

Every model of link growth – _i.e._ , a model that builds networks by
progressively adding connections to a set of initially isolated nodes – draws
a trajectory in the diagram starting from $(0,0)$ and, if not bounded to stop
earlier, reaching $(1,1)$ (see Figs. 1 and 2a). In this sense, we can regard
each real network as an intermediate stage of an unknown growth model, ideally
connecting the point $(0,0)$ to $(1,1)$. Such a framework provides us with a
metric to assess directly how efficient a given topology is from an overall
viewpoint: the normalized distance between the point representing the
considered topology and the upper right corner of the diagram (_i.e._ , the
final target of any network growth model).

Therefore, we adopt this metric, which we call _integrated efficiency_ , as a
comprehensive measure of the efficiency of spatial networks:

$E_{\text{int}}=1-\sqrt{\frac{(1-E_{\text{glob}})^{2}+(1-E_{\text{loc}})^{2}}{2}}\,.$
(3)

which is the Euclidean distance to the $(1,1)$ coordinates in the two-
dimensional space displayed in Fig. 2(a). The above definition satisfies a
crucial general consideration about the efficiency of real spatial networks:
They perform reasonably well at both local and global scale banavar-prl-2000 ;
katifori-prl-2010 . Indeed, this specific formulation encapsulates equally
both scales by rewarding the balance between the two efficiencies. Consider
the alternative, much simpler, measure
$E_{\text{int}}^{\prime}=(E_{\text{glob}}+E_{\text{loc}})/2$. Three
hypothetical “topologies” located at coordinates $(0,1)$, $(1,0)$, and
$(0.5,0.5)$, respectively, would score the same in terms of
$E_{\text{int}}^{\prime}$. On the contrary, the proposed measure
$E_{\text{int}}$ takes a higher value in the third case, enhancing the balance
between $E_{\text{glob}}$ and $E_{\text{loc}}$. The behavior of
$E_{\text{int}}$ for benchmark models as a function of $L_{\text{tot}}$ is
displayed in Fig. 2b.

For simplicity, we consider a formulation of $E_{\text{int}}$ where the
contributions of $E_{\text{loc}}$ and $E_{\text{glob}}$ are equal. However, it
is possible to alter Eq. (3) including a tunable multiplicative factor to
account for asymmetric contributions.

Figure 2: (panel a) Comparison between the averages local, $\langle
E_{\text{loc}}\rangle$, and global, $\langle E_{\text{glob}}\rangle$,
efficiencies for the GT, MST, EEM, and GN topologies. The values of the other
parameters are the same of those used in Fig. 1. Continuous black lines stand
for $E_{\text{int}}$ equilines. (panel b) Average integrated efficiency,
$\langle E_{\text{int}}\rangle$, as a function of the total length of the
system, $L_{\text{tot}}$.

##  III Comparing network designs

The measure introduced in the previous section informs whether a certain
topology is more, or less, efficient than another. Our final goal, however, is
to compare the design of spatial networks in terms of resource allocation,
something that is conceptually different. When we compute a network’s
integrated efficiency, we are calculating – by definition – how far it lies
from the complete graph, which is the only graph reaching maximum values of
$E_{\text{loc}}$ and $E_{\text{glob}}$ simultaneously. Nonetheless, there is a
limit to how close a system can get to such an extreme, which is conditioned
by the total cost $L_{\text{tot}}$. A simple solution to this issue would be
to divide the integrated efficiency, $E_{\text{int}}$, by
$L_{\text{tot}}/L_{\text{cg}}$, where $L_{\text{cg}}$ is the cost of the
complete graph with the same node layout. However, such a rescaling procedure
implies assuming a linear dependence of the efficiency on the total link
length which in general is not true. For the purpose of a fair comparison, it
is essential to avoid arbitrary assumptions that may introduce biases.

###  III.1 Quasi-exact comparisons. A numerical recipe.

Our proposal is based on a very simple idea: to build the best possible
network for a given amount of resources (_i.e._ , a given total cost
$L_{\text{tot}}$) and fixed node layout. In this manner, for every network
under study, we can construct its optimal counterpart and hence compute
$\Delta E_{\text{int}}=E_{\text{int}}^{\text{opt}}-E_{\text{int}}$, where
$E_{\text{int}}$ is the integrated efficiency of the original network and
$E_{\text{int}}^{\text{opt}}$ that of its optimal counterpart. This value
quantifies the room for improvement in a design with the same amount of
resources. More importantly, if we consider two systems with different total
cost and spatial scale, $\Delta E_{\text{int}}$ enable to perform an indirect
but fair comparison between them.

Our goal is not to build an improved version of real networks – which should
display all their ideal characteristics – but, rather, to generate appropriate
benchmarks by computing the maximum value that a given object function could
take under certain constraints.

As a consequence, we have designed a growth algorithm that works as follows:
we start considering an empty graph $G$ with $N$ nodes. Then, we add edges
iteratively until the total length of the graph reaches the desired one.
Adding edges means adding consecutive line segments in the
$(L,E_{\text{int}})$ plane (see Figs. 3 a and 4 a), connecting the origin with
the $(L_{\text{tot}},E_{\text{int}}(L_{\text{tot}}))$ point. Since we want to
maximize the value of $E_{\text{int}}(L)$ for $L=L_{\text{tot}}$, at each
iteration we have to look for the segment with the maximum vertical slope,
_i.e._ , the edge maximizing the ratio between the variation of
$E_{\text{int}}$ and the increase in total cost:

$\max\left\\{\frac{\tilde{E}_{\text{int}}-E_{\text{int}}}{\tilde{L}-L}\right\\}=\max_{\\{i,j\\}}\left\\{\frac{\tilde{E}_{\text{int}}-E_{\text{int}}}{d_{ij}}\right\\}\,,$
(4)

where $L$ and $E_{\text{int}}$ are the total weight of links and the
integrated efficiency at the current step, respectively, and $\tilde{L}$ and
$\tilde{E}_{\text{int}}$ stand for the same quantities after adding the link
$(i,j)$. Since the identification of the edge to add involves the evaluation
of the contribution of all the possible candidates, the overall procedure is
completely deterministic.

Even though it is not possible to ensure that the topologies produced by such
an algorithm reach the maximum possible value of the integrated efficiency,
there are strong hints that they are very close to it. The optimization of
$E_{\text{int}}$ is a non-Markovian process and there exists the chance that
different choices, locally not optimal, could lead to a better final result.
To address this issue, we have explored the possibility to use alternative
search methods based on simulated annealing. Our conclusion was that slightly
higher values of $E_{\text{int}}$ may be reached by a limited number of
alternative topologies – perhaps with a different balance between the local
and global efficiencies –, requiring a considerably higher computational cost
(_i.e._ , exploring the space of the configurations more exhaustively) in
order to be discovered222Therefore, it is technically a _quasi-optimal_
algorithm but we avoid that nomenclature for the sake of clarity..

In Fig. 3, we display the average behavior of $E_{\text{glob}}$,
$E_{\text{loc}}$, and $E_{\text{int}}$ against $L_{\text{tot}}$ for the
networks generated using our optimal algorithm. We also report the evolution
of the global and local efficiencies across the growth process according to
the bi-dimensional representation adopted in Fig. 2a. In particular, averages
are computed over one hundred random distributions of $N=100$ nodes within the
unit square. Panel (a) of Fig. 3 tells us that the algorithm first favors
increases in $E_{\text{loc}}$ up to a point where it is not possible to
increase $E_{\text{int}}$ further at the expense of $E_{\text{loc}}$. Unlike
$E_{\text{glob}}$, this metric has a non-monotonous behavior. For example, a
system made of separated cliques has a $E_{\text{loc}}=1$, and adding any
other link will reduce $E_{\text{loc}}$. This is the evolution’s least
interesting phase since, at this stage, systems are mainly composed by many
connected components while the overall connectivity is extremely low. After
the peak, the algorithm begins to link these isolated components at the
expense of $E_{\text{loc}}$ increasing $E_{\text{glob}}$. When the total cost
is roughly equal to that of the MST, the curves of $E_{\text{loc}}$,
$E_{\text{glob}}$, and $E_{\text{int}}$ merge together and start to behave in
the same way.

Figure 3: (panel a) Evolution of the local ($E_{\text{loc}}^{\text{opt}}$),
global ($E_{\text{glob}}^{\text{opt}}$), and integrated
($E_{\text{int}}^{\text{opt}}$) efficiencies with respect to the total length
of the system, $L_{\text{tot}}$, for the networks generated using our optimal
model. (panel b) Evolution of the average local ($\langle
E_{\text{loc}}\rangle$), and global ($\langle E_{\text{glob}}\rangle$)
efficiencies for the same networks. The green star denotes the values at which
90% of the nodes belong to the giant component. Results are averages over
$N_{real}=100$ different layouts of $N=100$ nodes uniformly distributed on the
unit square. Continuous black lines stand for $E_{\text{int}}$ equilines.

###  III.2 Approximated comparisons. An estimated universal upper bound.

The above algorithm provides an optimal network counterpart usable as a
reference. However, in practice, the increase of the system size $N$ severely
affects the runtime of the algorithm – in particular of the calculation of
$E_{\text{loc}}$, – jeopardizing the applicability of our methodology to large
size systems. To overcome this limitation, we determined the expected value of
$E_{\text{int}}$ for any value of $L_{\text{tot}}$ and any $N$, in the case of
layouts of nodes randomly distributed in a square. First, we studied the
behavior of the integrated efficiency against the overall cost, for several
layouts of $N\in\\{200,300,400\\}$ nodes. Figure 4a displays the behavior of
the $\langle E_{\text{int}}\rangle$ curves, confirming that the results are
robust across sizes. Specifically, the increase in size translates into a
shift of the $E_{\text{int}}$ curve towards higher values of $L_{\text{tot}}$.

By rescaling the $x$ coordinate of plots in Fig. 4, we collapse them into a
single, “universal” one which is the same for any value of $N$. This leads to
the definition of a _normalized total cost_ of a network $G$, $L^{\prime}$,
which reads:

$L^{\prime}=\frac{L_{\text{tot}}/\langle d\rangle}{\left(L_{\text{cg}}/\langle
d\rangle\right)^{\alpha}}\,,$ (5)

where $\langle d\rangle$ stands for the average spatial distance among nodes
and $L_{\text{cg}}$ is the length of the complete graph having the same node
layout as $G$. Such a normalized length can be rewritten as a combination of
two variables: $\langle d\rangle$ and $N$, since the length of a complete
graph is $L_{\text{cg}}=\langle d\rangle\cdot\frac{N(N-1)}{2}\approx\langle
d\rangle\frac{N^{2}}{2}$. Hence, we obtain:

$L^{\prime}\simeq\frac{L_{\text{tot}}}{\langle d\rangle\cdot N^{2\alpha}}\,.$
(6)

We have found that $\alpha=1/3$. As we can observe in Fig. 4b, the rescaling
of $L_{\text{tot}}$ returns perfectly overlapped curves. Such a new,
universal, curve allows us to compute the expected maximum value of
$E_{\text{int}}$ for a system with a given $L^{\prime}$, providing an
approximate upped bound to perform an indirect comparison between different
systems. Specifically, given an empirical network with a certain $L^{\prime}$,
we can use the difference between its actual value of integrated efficiency
and its expected maximum value, as a proxy of the system’s performance.

To improve the usability of this upper bound in real-world applications, it
would be extremely useful to provide an analytical expression for the
dependence of the expected maximum value of $E_{\text{int}}$ on $L^{\prime}$.
However, for $L^{\prime}<1$, networks are usually disconnected and the
behavior of the optimal integrated efficiency is very noisy (see also Fig.
3a). In order to avoid this difficulty, we restricted the range of admissible
values of the normalized length to $L^{\prime}\geq 0.91$ (thus discarding the
region where the slope of the curve is almost zero) and fitted our numerical
data to the relation:

$\bar{E}_{\mathrm{int}}^{\mathrm{opt}}(L^{\prime})=1-c_{1}(L^{\prime}-c_{2})^{-c3}\,.$
(7)

Using non-linear least squares, we found that $c_{1}=0.187$, $c_{2}=0.406$,
and $c_{3}=1.211$. In this way, given any real network, it is possible to
calculate the expected upper bound for its integrated efficiency without
generating any artificial network, but simply from the average distance
between its nodes and its total cost, through $L^{\prime}$ and Eq. (7).
Including shorter total link lengths (_i.e._ , $L^{\prime}\leq 0.91$) would
need a separate specific discussion which goes beyond the scope of the present
work.

Figure 4: (panel a) Evolution of the average integrated efficiency, $\langle
E_{\text{int}}^{\text{opt}}\rangle$, with respect to the total length
$L_{\text{tot}}$ for systems with different number of nodes $N$. Solid lines
account for the evolving systems, while symbols denote the final networks.
(panel b) Same quantity of panel a, but displayed as a function of the
rescaled length $L^{\prime}$. The curves corresponding to different network
sizes overlap perfectly.

##  IV Applications

To illustrate the applicability of the proposed analytical tools, here we use
them to compare some real spatial networks based on their integrated
efficiency. Specifically, we consider seven different collections of networks,
a few of them (UK-Flights, Cities, and Latium Vetus/Southern Etruria)
correspond to successive snapshots of the same evolving system. The
collections are:

UK Flights

Time-varying network of domestic flights in the United Kingdom between years
1990 and 2003 uk-airlines . Nodes correspond to airports, while an edge
between two airports accounts for the distance among them. For each
year/graph, we keep only those routes with, at least, 5000 carried passengers
across the year.

Cities

The evolution of urban street patterns of a small region in northern Italy,
captured in four snapshots between 1955 and 2007 strano-scirep-2012 . Nodes
correspond to the intersection between two streets or dead ends, while the
weight of an edge corresponds to the length of the street connecting two
nodes. For computational reasons, we consider only a small – rectangular –
sample of the whole dataset centered around a single village.

Latium Vetus/Southern Etruria

The networks of trails among settlements between 950 and 509 BC (Iron Age) in
two regions of Italy, namely: Latium Vetus (LV) fulminante-front-2017  and
Southern Etruria (SE) prignano-jas-2019 . Nodes represent settlements, while
an edge denotes a direct route connecting them. From older to earlier, we have
five snapshots for both regions: Early Iron Age 1 Early , Early Iron Age 1
Late, Early Iron Age 2, Orientalizing Age and Archaic Age. We label these
snapshots chronologically as (LV$x$/SE$x$) with $x\in[1,5]$. The snapshots do
not have the same duration, but the properties of the system are more or less
“stable” within each snapshot.

Catalonia railway

This network describes the current regional railway network in Catalonia in-
preparation1 . Nodes correspond to aggregated groups of contiguous towns,
while edges denote the length of the railway line connecting them.

Hispania roads

The networks of trails among cities and towns in Hispania (Iberian peninsula)
during the Roman Empire in-preparation2 . As for the Latium Vetus and Etruria
collections, nodes represent settlements, while an edge denotes a direct route
connecting them.

Rome railway

The network of rail connections in Rome, where nodes represent stops/stations
and link constitutes a direct connection between two nodes kujala-scidata-2018
. Weights correspond to the geodesic distance between both ends of the link.
The original data-set splits many stops in two, each one corresponding to the
two ways of the line. We simplify the network by merging stops having the same
name into a single node.

Power grid

A simplified model of the power grid network of Italy, where transmission
lines are assumed bidirectional and identical. It is a subset of a model of
the continental European electricity system, comprising 1494 buses and 2156
lines. power_dataset .

The main topological features of all these networks are reported in Tab. 1.
Such table reveals the diversity of the networks under analysis.

###  IV.1 Direct comparison of efficiencies

Efficiency values in Tab. 1 shows interesting features among the variety of
datasets. First, we see that $E_{\text{loc}}$ and $E_{\text{glob}}$ in the UK-
Flights fall within a narrow range of values centered around $0.6$, despite
the edge density $\rho$ is fairly high. On the other side, we noticed that all
the rest of networks are all _terrestrial_ , and tend to be more efficient at
a global rather than local scale. This is in line with the principles behind
the design and growth of such kind of networks, which tend to privilege tree-
like structures spanning the whole system at the expenses of resilience
banavar-prl-2000 ; barthelemy-jstat-2006 ; bohn-prl-2007 ; katifori-prl-2010 .
In this sense, terrestrial infrastructure networks are likely to show fairly
high $E_{\text{glob}}$, since they provide paths among all nodes with little
chance to large route factors. Nevertheless, exceptions are found. For
example, the rail network of the city of Rome shows two connected components,
dragging down the value of $E_{\text{glob}}$ compared to other similar
systems. On the other end, terrestrial networks are more vulnerable at the
local level, as denoted by their values of $E_{\text{loc}}$.

| $N$ | $K$ | $\rho$ | $L_{\text{tot}}$ | $L_{\text{cg}}$ | $\langle d\rangle$ | $L^{\prime}$ | $E_{\text{loc}}$ | $E_{\text{glob}}$ | $E_{\text{int}}$ | Ref.
---|---|---|---|---|---|---|---|---|---|---|---
|  |  | $(\%)$ |  $(Km)$ or $(m)$ |  |  |
UK Flights |  | uk-airlines
1992 | 41 | 130 | $15.85$ | $44833.310$ | $328656.403$ | $400.800$ | $12.149$ | $0.663$ | $0.636$ | $0.649$ |
1995 | 39 | 141 | $19.03$ | $50215.194$ | $304749.131$ | $411.267$ | $13.729$ | $0.688$ | $0.680$ | $0.684$ |
1998 | 41 | 146 | $17.80$ | $54019.694$ | $333391.976$ | $406.576$ | $14.431$ | $0.664$ | $0.669$ | $0.666$ |
2001 | 39 | 141 | $19.03$ | $52949.882$ | $305950.002$ | $412.888$ | $14.419$ | $0.553$ | $0.614$ | $0.582$ |
Cities |  | strano-scirep-2012
1955 | 29 | 36 | $8.87$ | $4474.270$ | $148386.739$ | $365.485$ | $1.692$ | $0.343$ | $0.819$ | $0.518$ |
1980 | 71 | 98 | $3.94$ | $12383.492$ | $1417732.535$ | $570.516$ | $1.618$ | $0.407$ | $0.824$ | $0.563$ |
1994 | 80 | 110 | $3.48$ | $13111.289$ | $1763976.734$ | $572.534$ | $1.587$ | $0.432$ | $0.818$ | $0.579$ |
2007 | 90 | 124 | $3.10$ | $14152.991$ | $2385978.985$ | $609.290$ | $1.485$ | $0.423$ | $0.819$ | $0.572$ |
Latium
Vetus |  | fulminante-front-2017
LV1 | 93 | 198 | $4.63$ | $1359.614$ | $111083.691$ | $25.966$ | $3.249$ | $0.678$ | $0.890$ | $0.760$ |
LV2 | 93 | 198 | $4.63$ | $1318.020$ | $108710.242$ | $25.411$ | $3.218$ | $0.648$ | $0.875$ | $0.736$ |
LV3 | 107 | 239 | $4.21$ | $1625.535$ | $143007.035$ | $25.217$ | $3.637$ | $0.679$ | $0.871$ | $0.755$ |
LV4 | 93 | 220 | $5.14$ | $1637.569$ | $102742.846$ | $24.017$ | $4.231$ | $0.687$ | $0.876$ | $0.762$ |
LV5 | 78 | 189 | $6.29$ | $1471.613$ | $75232.370$ | $25.052$ | $4.107$ | $0.731$ | $0.917$ | $0.801$ |
Southern
Etruria |  | prignano-jas-2019
SE1 | 116 | 199 | $2.98$ | $1586.445$ | $271499.897$ | $40.705$ | $2.082$ | $0.507$ | $0.875$ | $0.641$ |
SE2 | 115 | 207 | $3.16$ | $1660.375$ | $265460.422$ | $40.497$ | $2.204$ | $0.571$ | $0.887$ | $0.686$ |
SE3 | 130 | 235 | $2.80$ | $1811.973$ | $344291.037$ | $41.060$ | $2.183$ | $0.606$ | $0.869$ | $0.706$ |
SE4 | 168 | 311 | $2.22$ | $2062.108$ | $605568.932$ | $43.169$ | $1.989$ | $0.851$ | $0.649$ | $0.730$ |
SE5 | 179 | 366 | $2.30$ | $2300.959$ | $679017.117$ | $42.622$ | $2.153$ | $0.868$ | $0.714$ | $0.777$ |
Catalonia |  | in-preparation1
Railway | 34 | 37 | $6.60$ | $1040.580$ | $54481.616$ | $97.115$ | $1.299$ | $0.233$ | $0.637$ | $0.400$ |
Hispania |  | in-preparation2
Roads | 89 | 127 | $3.24$ | $10115.003$ | $1742076.283$ | $444.861$ | $1.453$ | $0.460$ | $0.812$ | $0.595$ |
Rome |  | kujala-scidata-2018
Railway | 80 | 103 | $3.26$ | $377.297$ | $39333.959$ | $12.447$ | $2.066$ | $0.246$ | $0.538$ | $0.375$ |
Power grid |  | power_dataset
Italy | 139 | 207 | $2.16$ | $11914.776$ | $4087233.378$ | $426.153$ | $1.316$ | $0.467$ | $0.794$ | $0.596$ |
Table 1: Summary of the topological indicators for all the empirical datasets.
For each network, we report its number of nodes, $N$, of edges, $K$, the edge
density, $\rho$, the total length for the empirical, $L_{\text{tot}}$, and
complete graph, $L_{\text{cg}}$, as well as the average spatial distance among
the nodes, $\langle d\rangle$, the rescaled length of the system,
$L^{\prime}$, and their local $E_{\text{loc}}$, global $E_{\text{glob}}$, and
integrated $E_{\text{int}}$ efficiencies. Finally, for each network dataset,
we report the bibliographic source of the data.

Beyond comparing the efficiencies of the different spatial networks in a
descriptive way, our methodology allows for a systematic comparison of diverse
spatial networks considering both their $E_{\text{int}}$ and $L_{tot}$.

###  IV.2 Systematic comparison considering total lenght

For each empirical network, we generate an optimized one using the algorithm
presented in Sec. III while preserving the node layout (_i.e._ , their
position) and the total length, $L_{\text{tot}}$.

First, we analyse the differences on both local and global efficiencies
separately. We compute the differences $\Delta
E_{X}=E^{\text{opt}}_{X}-E_{X}\,,\;X\in\\{loc,glob\\}$ among the efficiencies
of the empirical and optimal networks. In Fig. 5a we report the values of
$\Delta E_{\text{loc}}$ and $\Delta E_{\text{glob}}$ for all the cases under
scrutiny. It is worth noting that $\Delta E_{\text{loc}}$ is always negative,
while this is not the case of $\Delta E_{\text{glob}}$. With the exception of
UK-Flights and Rome Railway, all the other collections tend to have values of
$E_{\text{glob}}$ close to the optimal counterpart one (_i.e._ , $\Delta
E_{\text{glob}}=0$). A closer inspection to $\Delta E_{\text{glob}}$
highlights interesting features. One is that for Cities, Hispania, and
Catalonia Railway networks $\Delta E_{\text{glob}}<0$ (_i.e._ , they are more
efficient than their optimized counterparts). Latium Vetus, Etruria and the
Italian Power Grid, instead, fall very close to the optimum. Another
interesting feature is that UK-Flights and Rome Railway networks are sub
optimal both locally and globally, confirming our guess about the existence of
criteria behind their design other than merely spatial ones.

We sort all the datasets according to their $\Delta
E_{\text{int}}=E_{\text{int}}^{\text{opt}}-E_{\text{int}}$ and check how far
the real networks lie from the upper bound (Fig. 5b). The extent of
efficiency’s difference tells us how much better the systems could have
performed consuming the same amount of resources. The differences range from
$\Delta E_{\text{int}}\approx 15\%$ (City 2007, LV5) to above $50\%$ (Rome
Railway), while the majority of the real networks are about $20\%$ less
efficient than their optimal counterparts.

Figure 5: (panel a) Global efficiency difference, $\Delta
E_{\text{glob}}=E^{\text{opt}}_{\text{glob}}-E_{\text{glob}}$ , with respect
to the same quantity calculated for the local efficiency, $\Delta
E_{\text{loc}}$ for all the dataset reported in Tab. 1. The hue of the color
is used to order the dataset in a chronological way. (panel b) Ranking of the
empirical systems according to the integrated efficiency difference $\Delta
E_{\text{int}}=E_{\text{int}}^{\text{opt}}-E_{\text{int}}$. Points for data
SE4 and SE5 are not shown since, as mentioned in Sec. III.2, we were unable to
generate the optimal counterparts.

We have thus proven that it is relevant to consider the upper bound of the
integrated efficiency of each real network since it provides novel,
complementary information with respect to the mere value of $E_{\text{int}}$.
However, as discussed in Sec. III.2, the computational cost of determining
such upper bound increases rapidly with the size and total cost of the system
under scrutiny. Therefore, it is interesting to assess whether replacing the
value of $E_{\text{int}}^{\text{opt}}$ with its expected value provided by Eq.
(7) leads to similar results. By doing so, we are disregarding the details of
the node layouts, while still considering their overall characteristics
through the average node distance. In Fig. 6 we report the values of the
empirical $E_{\text{int}}$ as a function of the rescaled length $L^{\prime}$.
As expected, all the values of $E_{\text{int}}$ in empirical systems lie way
below the optimal curve. Generally speaking, networks with higher connectivity
are more efficient than those with less links, albeit being more costly.
However, this is not always true. By looking at Fig. 6 we notice that UK
Flights networks attain approximately the same efficiency of Southern Etruria
ones, but with a much higher cost.
We rank networks according to $\Delta
E_{\text{int}}^{\prime}=\bar{E}_{\mathrm{int}}^{\mathrm{opt}}-E_{\text{int}}$,
the difference between their $E_{\text{int}}$ and the corresponding value
computed through Eq. (7). We find that rankings based on $\Delta
E_{\text{int}}^{\prime}$ and on $\Delta E_{\text{int}}$ provide a Spearman
brank correlation of $r_{s}=0.618$ ($\text{p-val}=0.0037$). This indicates
that the curve is a valid alternative to ranking networks according to the
actual difference between optimal and empirical integrated efficiencies.
Nevertheless, the correlation’s value highlights non negligible discrepancies
between the rankings. The reason for such differences is that real spatial
layouts differ from layouts in our simulations. Random distribution of nodes
in a square show very little fluctuations as indicated by the shadow of the
curve. On the contrary, a real network’s layout can be far from this
distribution, thus affecting the output of a spatial network model
pablo_marti-scirep-2017 .

Figure 6: Integrated efficiency, $E_{\text{int}}$ as a function of the
rescaled length $L^{\prime}$ for the same datasets. The filled area for
$L^{\prime}\leq 0.91$ denotes the region for which the system has not fully
percolated yet. The hue of the color is used to order the dataset in a
chronological way. The continuous black line stand for the optimal curve
$E_{\text{int}}^{\text{opt}}(L^{\prime})$.

###  IV.3 Application to networks’ comparison in meaningful scenarios

We have shown that our methodology can be used with diverse spatial networks
in a plethora of application domains. Now, we aim at presenting how such
comparison can be used to address relevant research questions. In the
following, we use some of the empirical networks presented above to highlight
two meaningful ways of applying our methodology. The first example is an
archaeological case of study based on the Latium Vetus and Southern Etruria
datasets fulminante-front-2017 ; prignano-jas-2019 . These two data-sets
describe the pathway networks connecting settlements in two neighboring
regions with similar characteristics through time. Since there are obvious
differences in terms of covered distance, system size and total cost (length)
among the ten datasets, the networks corresponding to the two regions cannot
be compared directly. On the contrary, our methodology allows assessing which
region hosted better designed infrastructures by looking at the efficiency
difference to their optimal counterparts. A second example would involve
analysing the evolution of roads in the Cities’ dataset. There, we can assess
the efficiency in the area and the evolution of a system that increases both
in size and total cost, but for which we initially cannot know how better
could have performed under those constraints. Despite the strong changes in
structural terms (_e.g._ , size increase and density decrease), our
methodology would allow for a longitudinal comparison throughout the temporal
snapshots in the data-set. Even when the computation of the optimal is
jeopardized by the system size (like in 1994 and 2007) we could use the upper
bound curve to obtain their estimation.

##  V Conclusions

This manuscript provides tools to compare the design of different spatial
networks in terms of the duality _resources employed_ and _quality reached_.
Following an approach similar to zamora-commphys-2019 , such a comparison is
performed indirectly by contrasting the quality improvement that each one of
the networks’ designs could reach.

The paper presents the different components of our tool-set progressively.
First, we have introduced the notion of _integrated efficiency_ ,
$E_{\text{int}}$, as a metric to quantify spatial networks performance at
global and local scale simultaneously, while rewarding the balance between the
two. Second, we propose an algorithm to computationally estimate the upper
bound of our quality function for a given specific network: we have devised an
algorithm to generate networks with maximal $E_{\text{int}}$
($E_{\text{int}}^{\mathrm{opt}}$) with the same node layout and total cost as
in the original network. The smaller is the difference between such an upper
bound and the empirical value, the higher we consider the design quality of
the network under analysis to be.

Since the high computational cost of the optimal network algorithm may hinder
its applicability on large networks, we provide a universal expression for
approximated upper bound to any network. Considering a setting of $N$ nodes
randomly distributed in a unit square, we computed the expected maximal value
of $E_{\text{int}}$ ($\bar{E}_{\mathrm{int}}^{\mathrm{opt}}$) as a function of
the total cost $L_{\text{tot}}$. Then, by defining a rescaled total link
length, $L^{\prime}$, we successfully collapsed the
$\bar{E}_{\mathrm{int}}^{\mathrm{opt}}$ versus $L_{\text{tot}}$ curves for
different sizes onto a single one. In this way, we have been able to express
$\bar{E}_{\mathrm{int}}^{\mathrm{opt}}$ as a function of the number of nodes
$N$, the average distance between nodes $\langle d\rangle$, and the total cost
$L_{\text{tot}}$ of the network under study.

Finally, to test the applicability of our method, we have analyzed the
performance of a heterogeneous set of spatial networked systems. We have
checked that our approach provides new information beyond the mere comparison
between two networks’ efficiency. We provide also several examples of
application among the included datasets, showing how it helps to compare
networks with different size and total length in a meaningful way.

In conclusion, we have shown that a meaningful comparison of spatial networks
cannot be exempt from the definition of proper upper bounds with specific cost
constraints. This can be done (almost) exactly, by running our maximal
efficiency algorithm, or approximately, thanks to the universal curve. The
latter constitutes a good approximation for systems whose size does not make
the computation of $E_{\text{int}}^{\text{opt}}$ feasible. However, the
particularities of the layout (especially for low $L^{\prime}$) may affect the
precision of the method. In the future, it will be worth exploring how the
specificities of a layout affects a systems’ $E_{\text{int}}^{\text{opt}}$
with respect to the value provided by the curve. Another direction to pursue
is to explore the effect of weighting differently the contributions of the
local and global efficiencies in the computation of $E_{\text{int}}$.

## Data availability

The following datasets are available at data_repo : Latium Vetus, Southern
Etruria, Catalonia railway and Hispania roads. The rest of the datasets are
available online in the following references: Rome railway kujala-scidata-2018
and Power grid power_dataset . Regarding Cities datasets, please contact the
author (see Acknowledgements section). The code used for running the algorithm
can also be found in data_repo .

###### Acknowledgements.

The authors thank Simon Carrignon for his help with computational issues, and
Michael Gastner for helping with the implementation of the GN model. The
authors also thank Emanuele Strano, Francesca Fulminante, Lluc Font-Pomarol,
Pau de Soto, Tom Brughmans and Sergi Valverde for providing the data on
evolving cities, archeological networks, and the power grid.
AC acknowledges the support of the Spanish _Ministerio de Ciencia e
Innovacion_ (MICINN) through Grant IJCI-2017-34300. AC and SL acknowledge the
financial support of _Ministerio de Economia y Competitividad_ (MINECO)
through grant RYC-2012-01043. IM, LP and AD-G acknowledge the financial
support of the _European Research Council_ within the Advanced Grant EPNet
(340828). AD-G acknowledges financial support from FIS2015-71582-C2-2 -P
(MINECO/FEDER), and the Generalitat de Catalunya (Project No.2014SGR-608).

IM, AC, AD-G, LP and SL designed research and wrote the article; IM and AC
wrote code and performed simulations; IM analyzed numerical results.

##  Appendix A The Gastner-Newman model

We present here a brief description of the Gastner-Newman (GN) model to
generate spatial networks introduced in gastner-epjb-2006 . The main idea
behind the model is that there are two types of “costs” associated with a
given network: one related with the construction of the infrastructure, and
another related with its usage. Given a graph with $N$ nodes and $K$ edges,
$G(N,K)$ boccaletti-phys_rep-2006 , we consider its embedding in the
bidimensional space, $\mathbb{R}^{2}$. We denote with $d_{ij}$ the Euclidean
distance between nodes $i$ and $j$, respectively. Therefore, the total _cost
of construction_ of graph $G$, $T$, reads:

$T=\sum_{i=1}^{N}\sum_{j=i+1}^{N}a_{ij}\,d_{ij}\,.$ (8)

Where $a_{ij}$ is the element of the adjacency matrix, $\cal{A}$, of the graph
latora-book-2017 . The total _usage cost_ , $Z_{\lambda}$, instead, is:

$Z_{\lambda}=\sum_{i=1}^{N}\sum_{j=i+1}^{N}\tilde{l}_{ij}(\lambda)\,,$ (9)

with $\tilde{l}_{ij}$ being the _shortest path length_ between nodes $i$ and
$j$ which, in turn, is the sum of the lengths of the edges forming the _path_
between $i$ and $j$ latora-book-2017 . The path length depends on a parameter
$\lambda$ such that:

$\tilde{l}_{ij}(\lambda)=\lambda\sqrt{N}\,d_{ij}+(1-\lambda)\,,$ (10)

with $0\leq\lambda\leq 1$. Parameter $\lambda$ accounts for the users’
perception of distances. For $\lambda=0$, users give more importance to paths
made of few hops, while for $\lambda=1$ they pay more attention to shorter
paths (in terms of distance). Finally, the total cost, $C$, of the whole graph
is:

$C(G)=T+\gamma Z_{\lambda}\,,$ (11)

where parameter $\gamma\in[0,1]$ determines the relative weight between
construction and usage cost. The GN algorithm generates _optimal_ spatial
networks minimizing cost $C$. The cost optimization can be implemented using
either greedy or simulated annealing techniques newman-book-2012 . We decided
to implement the latter, since it ensures higher probabilities of finding the
optimal network. In our case, we were interested in building GN networks with
a given total cost. Hence, given a spatial network $G^{\star}$, we compute its
cost $C^{\star}$ according to Eq. (11). Considering the same nodes layout of
$G^{\star}$, we build a GN network $G$ through the following steps.

  1. 1.

Create the complete graph $G_{0}$, and compute its cost $C(0)$.

  2. 2.

At each step, $t$, perform with equal probability one of these two operations:

    1. (a)

Add/Remove an edge:
Choose a random pair of nodes $i$ and $j$, and if they are connected (_i.e._ ,
$\exists\,e_{ij}$) we remove the corresponding edge. Otherwise, we add the
edge. The removal can take place unless one of the two nodes has degree one
(_i.e._ , otherwise the node will get disconnected).

    2. (b)

Rewiring:
Choose an edge $e_{ij}$ at random. Then, choose a node $k\neq{i,j}$ at random
and create the edge $(i,k)$ or $(j,k)$ – if it does not exist already.
Finally, remove the edge $e_{ij}$.

  3. 3.

Compute the cost of the resulting graph, $C(t)$, and then:

    1. (a)

Accept the move with the following probability $p$:

$p=\begin{cases}\exp\Bigl{[}-\beta\bigl{(}C(t)-C(t-1)\bigr{)}\Bigr{]}&\text{
if }C(t)>C(t-1)\\\ 1&\text{ otherwise }\end{cases}$

    2. (b)

increase the value of $\beta$ of a quantity $\Delta\beta=1+3\cdot 10^{-6}$
(with $\beta(0)=\tfrac{0.1}{C_{MST}}$, and being $C_{MST}$ the cost of the
minimum spanning tree of the nodes layout).

  4. 4.

Repeat stages 2 and 3 either $t_{\infty}=1.5\times 10^{6}$ times, or until
$C\simeq C^{\star}$.

##  Appendix B The Greedy Triangulation

Here we provide a brief description on how to compute the _Greedy
Triangulation_ graph of a given layout of nodes embedded into a bidimensional
metric space $\mathbb{R}^{2}$. Given a set of $N$ nodes embedded into a two
dimensional space, the most connected (planar) triangulation graph $G^{P}$ has
– at most – $K=3N-6$ edges bondy-book-2008 ; and the maximally connected
triangulation minimizing its total length, $L$, is called the _minimum
weighted triangulation_ (MWT). Since no polynomial time algorithm is known to
compute the MWT, following buhl-epjb-2004 ; cardillo-pre-2006  we consider –
instead – its _greedy_ approximation known as Greedy Triangulation (GT).

Figure 7: Schematic representation of a segment $\overline{AB}$.

Taking $N$ points embedded into a bidimensional Euclidean space, we compute
first the corresponding complete weighted graph $G_{0}$ with
$K=\tfrac{N(N-1)}{2}$ edges. Then, we prepare a list of $G_{0}$’s edges,
sorted in ascending order of length (_e.g._ , using quicksort cormen-book-2001
). Given the schematization displayed in Fig. 7, for each edge/segment,
$\overline{AB}$, connecting nodes $A,B\in G_{0}$ we store the following
information:

$\left\\{\text{ID}_{A},x_{A},y_{A},x_{M},y_{M},\text{ID}_{B},x_{B},y_{B},r,m\right\\}\,.$

Where $\text{ID}_{i},x_{i},y_{i}$ are the ID, and $x,y$ coordinates of
node/point $i$, $M$ is the middle point of $\overline{AB}$,
$r=\tfrac{d_{AB}}{2}$ is its semi-length, and $m$ is its angular coefficient
(_i.e._ , the tangent of the angle between $\overline{AB}$ and the $x$-axis).
In its essence, the algorithm to compute the GT graph, $G_{GT}$, parses the
sorted edgelist of $G_{0}$ and checks whether each candidate edge
$e^{\star}\in G_{0}$ belongs to the GT or not. After adding the shortest
edge/segment of $G_{0}$ to the empty set of edges of $G_{GT}$, we check if any
other edge of $G_{0}$’s edge list intersects (and eventually how) with those
of $G_{GT}$. To check if two segments $\overline{AB}\in G_{GT}$ and
$\overline{CD}\in G_{0}$ intersect, – and assuming that $x_{A}\leq x_{B}$ and
$x_{C}\leq x_{D}$, – we compute the distance between their middle points
$d_{12}=d_{M_{\overline{AB}}\,M_{\overline{CD}}}$. Then:

  1. 1.

If $d_{12}>(r_{AB}+r_{CD})$ the two segments do not intersect, and thus
$\overline{CD}$ potentially belongs to $G_{GT}$.

  2. 2.

If $d_{12}\leq(r_{AB}+r_{CD})$ the two segments may intersect, and we need to
perform further checks to include/exclude the edge from $G_{GT}$.

To perform such checks, we have to compute the coordinates of the intersection
point, $(X,Y)$, which are:

$\begin{split}X&=\dfrac{y_{1}-m_{AB}\;x_{1}-\left(y_{2}-m_{CD}\;x_{2}\right)}{m_{CD}-m_{AB}}\,,\\\
Y&=\dfrac{m_{CD}\;y_{1}-m_{AB}\;y_{2}+m_{AB}\;m_{CD}\left(x_{2}-x_{1}\right)}{m_{CD}-m_{AB}}\,.\end{split}$
(12)

Where $x_{\alpha},y_{\alpha}\;\alpha\in\\{1,2\\}$ are the coordinates of the
middle point of $\overline{AB}$ if $\alpha=1$, or $\overline{CD}$ if
$\alpha=2$. If $X\in[\min(x_{A},x_{C}),\max(x_{B},x_{D})]$, then the two
segments will cross for sure (since $d_{12}\leq(r_{AB}+r_{CD})$) and
$\overline{CD}$ could be discarded (a similar criterion could be established
for $Y$). For each candidate edge $\overline{CD}$, we repeat the procedure
described above for all the edges $\overline{AB}$ already present in $G_{GT}$.
However, there are some exceptions to the “intersection” rule. In particular,
segments sharing one vertex do “technically” intersect, but without breaking
$G_{GT}$’s planarity and, hence, might belong to $G_{GT}$. Another case
requiring special attention is that of segments either parallel to one axis or
perpendicular to them. In such case the check on either $X$ or $Y$ alone is
not enough, and we must ensure that both $X$ and $Y$ fall outside their
respective intervals, instead. Lastly, for parallel segments (_i.e._ ,
$m_{CD}=m_{AB}$), the relations in Eqs. (12) have a singularity and, thus,
cannot be used to compute the coordinates of the intersection point. If
segment $\overline{CD}$ does not cross any of those of $G_{GT}$, it can be
added to $G_{GT}$ and we proceed to check the next candidate of $G_{0}$’s
edgelist. The algorithm stops either when $3N-6$ edges have been added to
$G_{GT}$, or if no more candidates to check are available. In the latter case,
the number of edges of the GT will be lesser than $3N-6$. This is due to the
fact that some edges between nodes laying at the outskirt of the node layout
might be added without breaking the planarity. However, such edges cannot be
represented as straight lines, and thus their intersection cannot be computed
using the above mentioned method. The amount of missing edges is approximately
in the order of $\sqrt{N}\ll 3N-6$.

##  Appendix C Equitable Efficiency Model

In this section, we present the essential traits of the Equitable Efficiency
Model (EEM) introduced by Prignano _et al._ in prignano-jas-2019 . EEM is a
growth model which builds networks from empty and static spatial node layouts.
In its essence, the model adds one link at a time ensuring that such addition
constitutes the best improvement in the efficiency of communication among any
pair of nodes.

Given a node layout embedded in a two dimensional space, $\mathbb{R}^{2}$, at
each step we calculate the route factor, $E_{ij}$, (_i.e._ , the ratio between
the spatial distance, $d_{ij}$, and the shortest path length, $l_{ij}$)
between all the pair of nodes $i$ and $j$. According to its definition,
$E_{ij}\in[0,1]\,\forall\,i,j$; where $E_{ij}=0$ when $i$ and $j$ belong to
different components of the system (_i.e._ , $l_{ij}=\infty$), and $E_{ij}=1$
when they are directly connected, instead. After computing all the values of
$E_{ij}$, we sort them in ascending order. The connection having the smallest
$E_{ij}$ is added to the network, and the above procedure is repeated
iteratively until the graph has a total length, $L_{\text{tot}}$, equal to the
desired one. However it is worth noting that, according to the definition of
route factor, $E_{ij}=0$ for all nodes belonging to different components,
regardless of their distance. To ensure a parsimonious usage of resources, and
avoid an arbitrary selection of one of the unconnected pairs, we ideally
replace $l_{ij}=\lim_{\Lambda\to\infty}\Lambda$ with $l_{ij}=\Lambda$ where
$\Lambda\gg L_{\text{cg}}$ is a large, but finite, length. This replacement
implies that the route factor between pairs of nodes belonging to different
components will be ranked according to their spatial distance. Therefore,
until the graph has one single component, the algorithm will select
connections between unreachable nodes, starting from those that are physically
closer to each other. This means that the set of links connecting the nodes
into a single component is nothing else than the Minimum Spanning Tree (MST)
of the layout under consideration.

## References

  * (1) P. Haggett, and R. J. Chorley, _Network Analysis in Geography_ (Edward Arnold, London, 1969). ISBN:0713154594
  * (2) M. Barthélemy, Spatial networks. Physics Reports, 499(1-3), 1–101 (2011). https://doi.org/10.1016/j.physrep.2010.11.002
  * (3) R. Guimerá, S. Mossa, A. Turtschi, and L. A. N. Amaral, Proc. Natl. Acad. Sci. USA, 102, 7794–9 (2005). https://doi.org/10.1073/pnas.0407994102
  * (4) M. T. Gastner and M. E. J. Newman Journal of Statistical Mechanics: Theory and Experiment, P01015 (2006). https://doi.org/10.1088/1742-5468/2006/01/P01015
  * (5) S. Porta, P. Crucitti, and V. Latora, Environment and Planning B: Planning and Design, 33, 705–725 (2006). https://doi.org/10.1068/b32045
  * (6) P. Kaluza, A. Kölzsch, M. T. Gastner, and B. Blasius, J. R. Soc. Interface, 7, 1093–1103 (2010). https://doi.org/10.1098/rsif.2009.0495
  * (7) A. Rinaldo, R. Rigon, J. R. Banavar, A. Maritan, and I. Rodriguez-Iturbe, Proc. Natl. Acad. Sci. USA 111, 2417–24 (2014). https://doi.org/10.1073/pnas.1322700111
  * (8) S. Lui, and G. Tiana, The Journal of Chemical Physics, 139, 155103 (2013). https://doi.org/10.1063/1.4826096
  * (9) K. S. Sandhu, G. Li, H. M. Poh, Y. L. K. Quek, Y. Y. Sia, S. Q. Peh, F. H. Mulawadi, J. Lim, M. Sikic, F. Menghi, A. Thalamuthu, W. K. Sung, X. Ruan, M. J. Fullwood, E. Liu, P. Csermely, and Y. Ruan, Cell Reports, 2, 1207–1219 (2012). https://doi.org/10.1016/j.celrep.2012.09.022
  * (10) E. T. Bullmore, and O. Sporns, Nature Rev. Neurosci., 13, 336–49 (2012). https://doi.org/10.1038/nrn3214
  * (11) A. Perna, and T. Latty, J. R. Soc. Interface 11, 20140334 (2014). https://doi.org/10.1098/rsif.2014.0334
  * (12) A. Bottinelli, E. van Wilgenburg, D. J. T. Sumpter, and T. Latty, J. R. Soc. Interface 12, 20150780 (2015). https://doi.org/10.1098/rsif.2015.0780
  * (13) Cook Z, Franks D W, Robinson E J H 2014 Behavioral Ecology and Sociobiology 68 3 509–517 https://doi.org/10.1007/s00265-013-1665-8
  * (14) L. Isella, J. Stehlé, A. Barrat, C. Cattuto, J.-F. Pinton, and W. Van den Broeck, Jour. Theo. Bio., 271, 166–180 (2010). https://doi.org/10.1016/j.jtbi.2010.11.033
  * (15) D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco, and A. Vespignani, Proc. Natl. Acad. Sci. USA, 106, 21484–9 (2009). https://doi.org/10.1073/pnas.0906910106
  * (16) C. Poletto, M. Tizzoni, and V. Colizza, Jour. Theo. Bio., 338, 41–58 (2013). https://doi.org/10.1016/j.jtbi.2013.08.032
  * (17) V. D. Blondel, A. Decuyper, and G. Krings, Eur. Phys. Jour. Data Science, 4, 10 (2015). https://doi.org/10.1140/epjds/s13688-015-0046-0
  * (18) I. Porat, and L. Benguigui, Europhys. Lett., 115, 18002 (2016). https://doi.org/10.1209/0295-5075/115/18002
  * (19) M. T. Gastner, and M. E. J. Newman, Eur. Phys. Jour. B, 49, 247–252 (2006). https://doi.org/10.1140/epjb/e2006-00046-8
  * (20) R. Louf, P. Jensen, and M. Barthélemy, Proc. Natl. Acad. Sci. USA, 110, 8824–8829 (2013). https://doi.org/10.1073/pnas.1222441110
  * (21) G. Zamora-López, and R. Brasselet, Communication Physics, 2, 144 (2019). https://doi.org/10.1038/s42005-019-0239-0
  * (22) A. Cardillo, S. Scellato, V. Latora, and S. Porta, Phys. Rev. E, 73, 66107 (2006). https://doi.org/10.1103/PhysRevE.73.066107
  * (23) V. Latora, V. Nicosia, and G. Russo, _Complex Networks: Principles, Methods and Applications_ (Cambridge University Press, 2017). https://doi.org/10.1017/9781316216002
  * (24) I. Vragović, E. Louis, and A. Díaz-Guilera, Phys. Rev. E,71, 036122 (2005). https://doi.org/10.1103/PhysRevE.71.036122
  * (25) V. Latora, and M. Marchiori, Phys. Rev. Lett., 87, 198701 (2001). https://doi.org/10.1103/PhysRevLett.87.198701
  * (26) J. R. Banavar, F. Colaiori, A. Flammini, A. Maritan, and A. Rinaldo, Phys. Rev. Lett., 84, 4745 (2000). https://doi.org/10.1103/PhysRevLett.84.4745
  * (27) M. Kaiser, and C. C. Hilgetag, Phys. Rev. E, 68, 036103 (2006). https://doi.org/10.1103/PhysRevE.69.036103
  * (28) M. Barthélemy, and A. Flammini, Journal of Statistical Mechanics: Theory and Experiment, L07002 (2006). https://doi.org/10.1088/1742-5468/2006/07/L07002
  * (29) M. Barthélemy, and A. Flammini, Phys. Rev. Lett., 100, 138702 (2008). https://doi.org/10.1103/PhysRevLett.100.138702
  * (30) L. Prignano, I. Morer, F. Fulminante, and S. Lozano, Journal of Archaeological Science, 105, 46–58 (2019). https://doi.org/10.1016/j.jas.2019.02.007
  * (31) E. Katifori, G. J. Szöllősi, and M. O. Magnasco, Phys. Rev. Lett., 104, 048704 (2010). https://doi.org/10.1103/PhysRevLett.104.048704
  * (32) United Kingdom Civil Aviation Authority: Airline data historic. https://www.caa.co.uk/Airline_data_historic/. Data available at: http://www.bifi.es/~cardillo/data.html (accessed 01/02/2018).
  * (33) E. Strano, V. Nicosia, V. Latora, S. Porta, and M. Barthélemy, Sci. Rep., 2, 296 (2012). https://doi.org/10.1038/srep00296
  * (34) F. Fulminante, L. Prignano, I. Morer, and S. Lozano, Front. Dig. Hum., 4, 1–12 (2017). https://doi.org/10.3389/fdigh.2017.00004
  * (35) P. De Soto, T. Brughmans, I. Morer, and L. Prignano, _In preparation_. (2018).
  * (36) I. Morer, L. Font-Pomarol, S. Lozano, and L. Prignano, _In preparation_. (2018).
  * (37) R. Kujala, C. Weckström, R.K. Darst, M.N. Mladenović, and J. Saramäki, Sci. Data, 5, 180089 (2018). https://doi.org/10.1038/sdata.2018.89
  * (38) T.V. Jensen, H. de Sevin, M. Greiner, and P. Pinson, The RE-Europe data set, Zenodo, (2015). Available at: http://doi.org/10.5281/zenodo.35177
  * (39) S. Bohn, and M. O. Magnasco, Phys. Rev. Lett., 98, 088702 (2007). https://doi.org/10.1103/PhysRevLett.98.088702
  * (40) F. Pablo-Martí, and A. Sánchez, Sci. Rep., 7, 4498 (2017). https://doi.org/10.1038/s41598-017-04892-2
  * (41) Ancillary files available at: https://arxiv.org/abs/1807.00565
  * (42) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. Hwang, Phys. Rep., 424, 175–308 (2006). https://doi.org/10.1016/j.physrep.2005.10.009
  * (43) M.E.J. Newman, _Computational Physics_. (Createspace Independent Publishing Platform, Scotts Valley, 2012). ISBN-13: 978-1480145511
  * (44) J.A. Bondy, U.S.R. Murty, _Graph Theory_. (Springer, 2008). ISBN:1846289696
  * (45) J. Buhl, J. Gautrais, R.V. Solé, P. Kuntz, S. Valverde, J.L. Deneubourg, and G. Theraulaz, Eur. Phys. Jour. B, 42, 123–129 (2004). https://doi.org/10.1140/epjb/e2004-00364-9
  * (46) T.H. Cormen, C.E. Leiserson, R.L. Rivest, and C. Stein, _Introduction to Algorithms_. (MIT Press, 2001). ISBN:0262032937