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. ). Then, in Sec.  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. , 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.
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 \(\lbrace  a_{ij}\rbrace \)  take value 1 \((0)\)  if the connection exists (does not exist), and \(\lbrace  d_{ij}\rbrace \)  take finite positive values corresponding to the spatial distance between nodes \(i\)  and \(j\)  [1]}. These two matrices fully determine the shortest paths matrix \(\mathcal {L}\) , whose elements \(\lbrace  l_{ij}\rbrace \)  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 stepsIt 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 [2]}. By computing the average over all pairs of nodes we obtain the global efficiency as proposed in [3]}, (despite an alternative definition has been given in [4]}):
\(E_{\text{glob}}= \frac{1}{N\left( N-1\right)} \sum _{i\ne j} \frac{d_{ij}}{l_{ij}} \,.\) 
\(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 [4]} 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. [3]} 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\ne m \in \Gamma _i} \frac{d_{jm}}{l_{jm/i}} \,,\) 
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 [7]}, [8]}, [9]}, [10]}, [11]}, [12]}, we selected a few models as benchmarks, thus encompassing a wide spectrum of possibilities:

the Minimum Spanning Tree (MST) [2]};

the Greedy Triangulation (GT), a maximally connected triangulation minimizing its total length [14]};

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 [15]};

the Gastner-Newman model (GN), a growth model where edges have an effective length which combines Euclidean and topological distances altogether [8]}.

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. REF  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. ), showed that the distribution of distances among the points is compatible with a random uniform distribution.
<FIGURE>In Fig. REF , 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. REF  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. REF  and REF a).
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}} \,.\) 
which is the Euclidean distance to the \((1,1)\)  coordinates in the two-dimensional space displayed in Fig. REF (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 [9]}, [18]}. 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. REF b.
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. (REF ) including a tunable multiplicative factor to account for asymmetric contributions.
<FIGURE>
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.
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. REF  a and REF  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\lbrace \frac{\tilde{E}_{\text{int}} - E_{\text{int}}}{\tilde{L}-L}\right\rbrace = \max _{\lbrace i,j\rbrace }\left\lbrace \frac{\tilde{E}_{\text{int}} - E_{\text{int}}}{d_{ij}}\right\rbrace  \,,\) 
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 discoveredTherefore, it is technically a quasi-optimal algorithm but we avoid that nomenclature for the sake of clarity..
In Fig. REF , 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. REF a. In particular, averages are computed over one hundred random distributions of \(N=100\)  nodes within the unit square. Panel (a) of Fig. REF  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>
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 \lbrace  200, 300, 400 \rbrace \)  nodes. Figure REF a 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. REF , 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 }} \,,\) 
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 }} \,.\) 
We have found that \(\alpha =1 / 3\) .
As we can observe in Fig. REF b, 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. REF a).
In order to avoid this difficulty, we restricted the range of admissible values of the normalized length to \(L^{\prime }\ge 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} \,.\) 
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. (REF ). Including shorter total link lengths (i.e., \(L^\prime \le 0.91\) ) would need a separate specific discussion which goes beyond the scope of the present work.
<FIGURE>
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 [19]}. 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 [20]}. 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) [21]} and Southern Etruria (SE) [15]}. 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 [23]}. 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 [24]}. 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 [25]}. 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. [26]}.

The main topological features of all these networks are reported in Tab. REF . Such table reveals the diversity of the networks under analysis.
Direct comparison of efficiencies
Efficiency values in Tab. REF  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 [9]}, [11]}, [29]}, [18]}. 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}}\) .
<TABLE>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}\) .

Systematic comparison considering total lenght
For each empirical network, we generate an optimized one using the algorithm presented in Sec.  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 \lbrace  loc, glob \rbrace \)  among the efficiencies of the empirical and optimal networks. In Fig. REF a 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. REF b). 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>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. REF , 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. (REF ) 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. REF  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. REF  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. (REF ). 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 [31]}.
<FIGURE>
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 [21]}, [15]}. 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.

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 [34]}, 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 [35]}: Latium Vetus, Southern Etruria, Catalonia railway and Hispania roads. The rest of the datasets are available online in the following references: Rome railway [25]} and Power grid [26]}. Regarding Cities datasets, please contact the author (see Acknowledgements section). The code used for running the algorithm can also be found in [35]}.
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.

The Gastner-Newman model
We present here a brief description of the Gastner-Newman (GN) model to generate spatial networks introduced in [8]}. 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)\)  [40]}, 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}\,.\) 
Where \(a_{ij}\)  is the element of the adjacency matrix, \(\cal {A}\) , of the graph [1]}. The total usage cost, \(Z_\lambda \) , instead, is:
\(Z_\lambda = \sum _{i=1}^N \sum _{j=i+1}^N \tilde{l}_{ij}(\lambda ) \,,\) 
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\)  [1]}. The path length depends on a parameter \(\lambda \)  such that:
\(\tilde{l}_{ij}(\lambda ) = \lambda \sqrt{N} \, d_{ij} + (1 - \lambda ) \,,\) 
with \(0 \le \lambda \le 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 \,,\) 
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 [43]}. 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. (REF ). Considering the same nodes layout of \(G^\star \) , we build a GN network \(G\)  through the following steps.

Create the complete graph \(G_0\) , and compute its cost \(C(0)\) .

At each step, \(t\) , perform with equal probability one of these two operations:

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).

Rewiring:
Choose an edge \(e_{ij}\)  at random. Then, choose a node \(k \ne {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}\) .


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

Accept the move with the following probability \(p\) :
\(\nonumber p ={\left\lbrace \begin{array}{ll}\exp \Bigl [ -\beta \bigl ( C(t) - C(t-1)\bigr ) \Bigr ] & \text{ if } C(t) > C(t-1)\\1 & \text{ otherwise }\end{array}\right.}\) 

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).


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


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 [44]}; 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 [45]}, [14]} we consider – instead – its greedy approximation known as Greedy Triangulation (GT).
<FIGURE>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 [47]}). Given the schematization displayed in Fig. REF , for each edge/segment, \(\overline{AB}\) , connecting nodes \(A,B \in G_0\)  we store the following information:
\(\nonumber \left\lbrace \text{ID}_A, x_A, y_A, x_M, y_M, \text{ID}_B, x_B, y_B, r, m \right\rbrace \,.\) 
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 \le x_B\)  and \(x_C \le x_D\) , – we compute the distance between their middle points \(d_{12} = d_{M_{\overline{AB}} \, M_{\overline{CD}}}\) . Then:

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

If \(d_{12} \le (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}\) 
Where \(x_\alpha , y_\alpha \; \alpha \in \lbrace 1,2\rbrace \)  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} \le (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. (REF ) 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\) .

Equitable Efficiency Model
In this section, we present the essential traits of the Equitable Efficiency Model (EEM) introduced by Prignano et al. in [15]}. 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 \rightarrow \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.
