Orientations and matrix function-based centralities in multiplex network analysis of urban public transport

We study urban public transport systems by means of multiplex networks in which stops are represented as nodes and each line is represented by a layer. We determine and visualize public transport network orientations and compare them with street network orientations of the 36 largest German as well as 18 selected major European cities. We find that German urban public transport networks are mainly oriented in a direction close to the cardinal east-west axis, which usually coincides with one of two orthogonal preferential directions of the corresponding street network. While this behavior is present in only a subset of the considered European cities it remains true that none but one considered public transport network has a distinct north-south-like preferential orientation. Furthermore, we study the applicability of the class of matrix function-based centrality measures, which has recently been generalized from single-layer networks to layer-coupled multiplex networks, to our more general urban multiplex framework. Numerical experiments based on highly efficient and scalable methods from numerical linear algebra show promising results, which are in line with previous studies. The centrality measures allow detailed insights into geometrical properties of urban systems such as the spatial distribution of major transport axes, which can not be inferred from orientation plots. We comment on advantages over existing methodology, elaborate on the comparison of different measures and weight models, and present detailed hyper-parameter studies. All results are illustrated by demonstrative graphical representations.

1. Introduction.By the year 2050 two-thirds of the growing world population are expected to live in urban areas [21].The sustainable planning and development of cities will greatly impact global economic, environmental, and social challenges, which demand interdisciplinary approaches to urban science [2].
In recent decades, the interdisciplinary field of complex network science has provided models and methods, which today impact everyday life [51,71,6,23,54].Network modeling approaches for urban systems can be traced back almost 300 years [33] and the abstraction of urban areas into geometrical models continues to form the basis of modern urban science [57,56,12,10,7,25,24,9,8,60].The combination of complex network models, urban science applications, and mathematical methodology from the field of numerical linear algebra is at the heart of this paper.
The main contribution of this paper is the development of a general methodology for two aspects of the spatial analysis of urban public transport networks.Exemplary results are given for major German and European cities. Figure 1 gives a schematic overview of the developed methods, which we describe further in the following paragraphs.This paper is accompanied by publicly available python implementations of all developed methods 1 .
The first aspect of geometrical city modeling studied in this paper is represented by the determination of orientations of complex urban networks, which reflect different urban organization patterns ranging from top-down planning to self-organization Fig. 1.Schematic overview of the methods introduced in this paper.Both aspects of geometrical city modeling rely on GTFS timetable data, which is described in the section Data.The methodology for the computation of public transport orientations using sequences of single-layer networks is described in the section Methodology.Details on the construction of the multiplex network representation of urban public transport systems are given in the section Multiplex network model.The methodology for the determination of the centrality of stops and lines of the multiplex networks using matrix function-based centrality measures is presented in the sections Definition and Numerical methods.
dynamics [12,10,9].Earlier works focused on the study of orientations of street networks [25,24,41,52,20].Our contribution to this aspect of geometrical urban analysis is to provide methodology and example results of orientations of public transport systems.We compare these to street network orientations obtained with existing methodology [19,20] and find interesting relations between the orientations of the two modes of transportation for the largest German and several major European cities.
The second aspect studied in this paper is the identification and ranking of the most central nodes of urban public transport networks.Complex network scientists have intensively studied a variety of different centrality measures in recent decades.Among the most prominent examples are degree, betweenness [34], closeness [35], and versions of eigenvector centrality [22,23,48,54].Urban science has been one of many disciplines to pose the natural question of identifying and ranking the most important entities of complex networks [27,26,57,58,66,53,3,42,44,28].
Our approach to the computation of centralities of urban public transport networks relies on two main ingredients: matrix function-based centrality measures and multiplex networks.The class of matrix function-based centrality measures [45,32,31,30,14,13] has attracted a lot of attention in recent years but has not been used in the urban science literature yet.These measures have the interesting property to interpolate between the concepts of local degree and global eigenvector centrality [15].Furthermore, these measures require less assumptions on the graph structure than classical eigenvector centrality, which makes them applicable to a wider range of problems.
For the network modeling we rely on a special class of multilayer networks.Re-search on multilayer networks exploded since 2014 when two survey papers unified terminology and gleaned the main aspects of a field, which had before been studied across various disciplines [47,18].Multilayer networks provide the means for building increasingly realistic models of complex systems as they allow entities to interact in different ways and on various levels.Urban science was one discipline to witness the rise of multilayer network modeling approaches in recent years [62,8,4,5,73,28].Some recent works provide generalizations of centrality measures well-studied on single-layer graphs to the case of different multilayer architectures [29,63,69,67,64,72,65,16]. Most importantly for this paper, the class of matrix functionbased centrality measures has very recently been generalized to the case of layercoupled multiplex networks.We adapt that methodology and study its applicability for the more general multiplex models of urban public transport networks.Highly efficient and scalable methods from numerical linear algebra enable the rapid and stable approximation of the centralities for small-to large-scale networks.We present various numerical results for multiplex networks with the number of physical nodes in the thousands and the number of layers in the hundrets and extensively study the influence of the different hyper-parameters involved in the model.

2.
Data.The focus of this paper lies on the analysis of urban public transport networks.For their realistic modeling we rely on the General Transit Feed Specification (GTFS) 2 , which defines a standardized format for public transport timetables used around the globe.For comparison and visualization purposes we also consider street networks.Here, we rely on the python package OSMnx3 [19] built on top of the open source platform OpenStreetMap 4 .OSMnx can, amongst other things, be used to generate graph structures representing street networks according to userspecified queries and provides various routines for their analysis.Note, that while OpenStreetMap also contains a public transport tagging feature, timetable information can not be expected to be as accurate as high-quality GTFS data [1].
We consider the 36 largest German cities by population as well as 18 selected major European cities for which we could obtain complete GTFS data.The three German cities Bochum, Gelsenkirchen, and Magdeburg had to be excluded from our studies as no complete GTFS data was available.
For Germany, several high quality GTFS feeds are publicly available5 on a daily basis and this paper builds on the local public transport data set 67 for April 22nd, 2021 (feed version "light-2021-04-22").This data set covers over 20 000 public transport lines serving over 450 000 stops across Germany.However, only required and conditionally required files (no optional ones) are provided8 , e.g., no explicit information about route frequencies or transfer options between stops is provided.
We also obtained several individual GTFS data sets for 18 selected major European cities 9 .However, as data integrity and availability varies among the different feeds, no common feed date could be guaranteed for all data sets.Instead, we downloaded the latest complete data set available as of May 21st, 2021 in each case.
For all considered cities, we use polygons representing the administrative bound-aries of the cities proper (excluding suburban areas) obtained via appropriate OSMnx queries to filter for all stops within the city limits.In GTFS terminology, all "routes" with at least one associated "trip" connecting at least two of the filtered "stops" are considered a valid line in the public transport network.In addition to the stops' spatial coordinates, we process the fields "arrival time" and "departure time" of each "trip" from the file "stop times" to determine travel times between connected pairs of stops.
The polygons used to filter public transport stops from the GTFS data sets are also deployed in all OSMnx routines used in the context of this paper.For our analysis, we require filtering for streets within each city's limits for two purposes: firstly, to produce street network orientation plots and, secondly, to create and plot street networks for visualization purposes later in this paper.
The following two sections provide details on how the data is transformed into different graph structures.While a sequence of single-layered graphs is sufficient for the determination of public transport network orientations we introduce a more sophisticated modeling approach based on multiplex networks for the determination of the most central stops and lines of the networks.
3. Public transport network orientations.Urban areas in large parts of the world have undergone rapid expansion in past decades and are expected to grow further [21].Global economic, environmental, and social challenges will be crucially impacted by future urban planning [2].One of many important aspects of the growth of urban systems is their geometrical expansion [57,56,12,10,7,25,24,9,8,60].Previous works have identified very different mechanisms for the spatial evolution of cities ranging from top-down planning to self-organization dynamics [12,10,9].A recent work beautifully illustrates these different mechanisms by means of orientation plots of urban street networks on a global scale [20].
In this section, we consider orientations of both street and public transport networks of the 36 largest German cities by population as well as 18 selected major European cities specified in the previous section.We describe the developed methodology in a first, and illustrate and discuss the obtained results in a second step.
3.1.Methodology.For the determination of street network orientations we use readily implemented routines 10 from the python package OSMnx [19] in combination with appropriate search queries, which follow the methodology described in [20].This methodology relies on the creation of a primal [57] undirected single-layer graph in which intersections are represented by nodes and streets by straight edges between them (ignoring curvature).The compass bearings of both directions (always including the reciprocal of any street bearing) of each street segment are recorded and added in an unweighted manner, cf.[52] for a discussion on weighting techniques.Afterwards, they are divided into 36 equal-sized bins with the cardinal directions located in the middle of the associated bins and plotted as a rose diagram.
Our first contribution to the geometrical analysis of urban systems is the investigation of orientations in public transport networks.To this end, we use GTFS data describing public transport timetables as specified in the section Data and we develop a methodology for the computation of public transport network orientations similar to that of street network orientations.
Our methodology also relies on the creation of single-layered graphs and the computation of bearings between pairs of nodes.In our case, nodes represent stops within the city limits and edges represent connections of public transport lines between stops, which we also assume to be straight and which are also not weighted (by, e.g., segment lengths or travel times).Deviating from the street network case, we consider these edges to be directed.The edges' directions are determined by the sequence of stops of, in GTFS terminology, each "trip" of each valid "route" of the urban public transport network.We realize this by successively processing all "trips" belonging to all valid "routes" as specified in the section Data.This procedure has the effect of assigning an implicit weight to all lines proportional to their operation frequency.An alternative approach would be to only process unique trips of all lines, which would, in the simplest case, lead to two copies of the same sequence of stops in reverse order.Interestingly, we empirically found only small deviations in the results of the two approaches, which we will comment on in the subsection Results and discussion.In both cases, most of the created directed graphs are chain graphs connecting only a small subset of the stops within the city limits.We then store the directed bearings (not including reciprocals) between all pairs of connected nodes of all graphs, divide them into 36 equal-sized bins, and visualize the results in a rose diagram as in the case of street network orientations.Note that the computations of orientations of different public transport lines are entirely independent, which offers great potential for parallelization.

Results and discussion.
The methodology for the computation of street network orientations described in the beginning of the section Methodology leads to the results displayed in Figure 2 for the 36 German cities under investigation.The results show that most large German cities' street networks have two orthogonal preferential directions, which is consistent with previous findings [24].How strongly pronounced these directions are, however, differs significantly: while, e.g., Halle (Saale) or Krefeld exhibit quite distinct cross-shaped patterns throughout the city area other cities like, e.g., Bielefeld or Mönchengladbach show almost equally distributed orientations.These differences likely reflect diverse historic city developments and urban planning approaches.Furthermore, the preferential directions seldomly coincide with the cardinal directions, but can often be linked to geographical constraints such as rivers or mountains.Some cities like, e.g., Lübeck or Rostock appear to comprise of two sets of orthogonal preferential directions prevailing in different regions of the city area.
At first glance, one might expect public transport network orientations to strongly correlate with street network orientations as buses and usually trams, which dominate the public transport system in most German cities, cf. the layer column in Table 1, operate on or alongside streets.However, as public transport stops are usually not located at street intersections but on street segments one could interpret public transport networks as dual graphs [56] of street networks in which nodes represent streets and edges represent intersections and in which only a subset of street segments is equipped with public transport stops.In the case of a line following the same street between two stops the orientations of both networks will coincide, but as soon as a public transport line takes a turn at a street intersection, the public transport bearing will point into a direction in between the two (often orthogonal) street orientations.
Figure 3 shows the public transport network orientations for the 36 considered German cities.In contrast to the street network orientations of the same cities in Figure 2 the public transport network orientations lack clear orthogonal preferential directions for most cities.Instead, one often observes only one blurred preferential  2. Street network orientations of the 36 largest German cities with valid GTFS data as described in the section Data.The plots were created using the OSMnx python package [19].
direction consisting of several bins of the rose diagram with no or only a weakly pronounced second orthogonal direction.Interestingly, this one blurred preferential direction often approximately coincides with the one preferential direction of the corresponding street network, which is closer to the east-west axis.None of the considered German cities has its most pronounced public transport direction closer than 45 degrees to the north-south axis.
Moreover, the blurred tilted east-west-like preferential directions of the public transport networks are often accompanied by a visually rather discontinuous distribution of bearings in the remaining directions.As there is a tendency for more continuous bearing distributions for cities with many public transport lines we conjecture that part of the recorded orientations can be considered random, e.g., due to the effect induced by taking turns described above, which decouples public transport orientations from street network orientations.Cities with more lines and consequently a larger sample size of different bearings seem to get closer to the actual, relatively smooth bearing distribution.For instance, the bearing distribution of Mainz with only 39 lines appears much more discontinuous than that of Hamburg with 253 lines.
By the choice of considering directed edges in public transport networks the corresponding orientation diagrams are nonsymmetric.However, the degree of nonsymmetry in Figure 3 is rather small indicating that most lines serve the same sequence of stops roughly equally in both travel directions.
Comparing the two approaches of considering all "trips" per "route" against considering only unique "trips" the orientation plots using only unique trips are visually almost identical to those displayed in Figure 3 for almost all 36 cities. Very few cities like Freiburg and Krefeld, which both comprise of few lines contain somewhat less pronounced tilted east-west axes, implying that the respective bearings belong to highly frequented routes.
For comparison, we also apply the methods described in section Methodology to 18 selected major European cities. Figure 4 shows their street orientation diagrams while Figure 5 illustrates the corresponding public transport orientations.The key Table 1 Multiplex network size of the 36 largest German cities with valid GTFS data as described in the section Data.Nodes correspond to stops and layers correspond to public transport lines, i.e. e.g., the public transport network of Augsburg consists of 5 tram lines, 63 bus lines, and a total of 68 public transport lines.Each line only serves a subset of all stops of the city, which is typically small compared to the full network.This entails that most node-layer pairs in the individual layers and consequently in the full networks are isolated, i.e, possess no incident edge.Hence, the number of non-isolated node-layer pairs and the number of intra-layer edges is similar for all considered networks.
observations made for the German cities hold true for roughly half of the considered European cities: street network orientations contain two orthogonal preferential orientations, which can sometimes be linked to geographical constraints, and main public transport network orientations tend to run along the east-west axis.Deviating from the latter observation, merely Nice clearly possesses a more pronounced north-south orientation.Conversely, street orientations of, e.g., Athens, Rome, and Stockholm are quite equally distributed, which could suggest that self-organizing dynamics dominated the spatial expansion of these cities [12,10,9].The fact that Athens and Rome both have a history of being important European cultural centers for millennia might support that hypothesis.Interestingly, some cities like Luxembourg City, Oslo, and Stockholm, whose street network orientations are quite equally distributed, still show a clearly pronounced (tilted) east-west axis in the public transport network orientation plots.Finally, we comment on some European cities with unique properties in at least one of the two orientation plots.The street network orientations of Brussels differ from other considered cities in the fact that the preferential directions are not orthogonal, but rather form an angle of 45 degrees.The public transport network orientations of Belgrade are remarkably nonsymmetric suggesting that a transport concept different from lines going back and forth along the same route with the same frequencies is in place.Lastly, Barcelona is the only city considered in this paper with two equally pronounced orthogonal preferential directions in both orientation plots.The latter observation is the result of a remarkable city planning effort by Cerdà in the year 1860 to connect the medieval core of Barcelona with surrounding villages by an extensive road network in the form of an orthogonal grid [55], which is certainly only one of many stories behind the graphics in Figures 2 to 5.
All results presented in this section are reproducible with publicly available python codes 11 .The code is designed according to the General Transit Feed Specification (GTFS) and should work on any valid GTFS data set making it easily utilizable for data of all regions around the world.
4. Multiplex network model.We now move to the second aspect of geometrical city modeling studied in this paper: the application of matrix function-based centrality measures to a multiplex network representation of the urban public transport systems.To this end, we introduce the multiplex framework employed to formalize these systems.We discuss the relation between network centralities and orientations in the section Discussion.
Network models have a long history in urban science: going back almost 300 years to Euler's Königsberg bridge problem [33] it could be argued that graph theory was motivated by an urban science problem.Especially the abstraction of street networks to (almost) planar single-layer graphs has been the basis for many mathematical models of cities ever since [7,8].More recently, following the advent of multilayer graphs in various scientific disciplines [47,18], urban science profited from the flexibility provided by multi-layered network structures to construct more realistic models of complex interactions [5,28,4,73,62].Some of these works consider multi-layered public transport networks in which layers correspond to either lines or modes of transportation.
In this section, we introduce a primal [57] multiplex representation of public transport networks in which layers represent lines.We assume the layers to be nodealigned so that each layer contains an instance of each physical node representing a stop within the city limits.We denote the number of physical nodes by n and the number of layers by L. Each copy of a physical node in any of the layers is called a node-layer pair.
We distinguish between two types of edges that can connect pairs of node-layer pairs: intra-and inter-layer edges.Intra-layer edges connect node-layer pairs belonging to the same layer whenever the corresponding line directly connects the two stops, i.e., if a line approaches stops A-B-C, then A and B as well as B and C are connected by an edge, but A and C are not.Inter-layer edges, on the other hand, connect instances of the same physical nodes belonging to different layers whenever both lines serve the corresponding stop, i.e., if the stop can be used to change between the lines.We discuss modeling approaches to assign weights to both intra-and inter-layer connections later in this section.
In our multiplex networks we restrict ourselves to undirected edges although directed models taking travel directions into account would appear suitable for the application, cf.our approach to network orientations in the section Public transport network orientations.Our choice accounts for the structure of the GTFS data, which contains multiple different stop-ID's, i.e., physical nodes per stop.This structure entails that, e.g., a bus stop with one stop on each side of a street but the same name is not represented by one but two physical nodes.Employing a strategy to aggregate stops with the same name is error-prone as it provokes unexpected behavior when different data providers use inconsistent stop naming logics, e.g., in-and excluding track numbers in stop names.In order to prevent in-and outbound traffic of a stop to be unequally distributed across two different physical nodes we use undirected edges, which equally represent in-and outbound traffic at the involved node-layer pair regardless of the current travel direction.In the example scenario described above this leads to both stop-ID's carrying the same in-and outbound information.This way of "counting each connection twice" appears preferable over each node-layer pair of the network carrying only half of the information available for the respective stop.
More formally, the multiplex networks described so far can be defined as the multilayer graph G = ( Ṽ, E (1) , . . ., E (L) , Ẽ) consisting of a common vertex set Ṽ for all layers, intra-layer edge sets E (l) for all layers l = 1, 2, . . ., L, and an inter-layer edge set Ẽ. Note that similar networks have been considered before, e.g., in [63,64,65,16] but in this paper we employ a different notion of inter-layer edges, which is determined by the data.
We choose a supra-adjacency matrix representation as the linear algebraic formulation of these multiplex networks [47,63,64,65,16].The two different types of edges are represented by two separate matrices.The multilayer intra-layer adjacency matrix A intra ∈ R nL×nL contains the edges representing connections by public transport lines and is defined as the block-diagonal matrix (4.1) , with 0 ∈ R n×n the matrix of all zeros.Each block-diagonal entry A (l) ∈ R n×n corresponds to the adjacency matrix of layer l.It contains the weight of the edge between the physical nodes i and j in layer l in the entry [A (l) ] ij and zero if no edge is present for i, j ∈ {1, 2, . . ., n} and l ∈ {1, 2, . . ., L}.
, where [O (lk) ] ij = 1, if i = j and physical node x i is connected between layers l and k, 0, otherwise, with l, k ∈ {1, 2, . . ., L} such that O (lk) ∈ R n×n contains ones only on a subset of its diagonal entries.Thus, the presence of inter-layer edges depends on whether the lines represented by layers l and k both stop at a given physical node.We set the blockdiagonal of A inter to zero matrices 0 ∈ R n×n as these entries would not represent a change of lines.
On the one hand, the above definition of A inter is more general than those in [63,64,65,16], where each O (lk) is a full identity matrix, as it includes all combinations of up to n zero diagonal entries.On the other hand, the construction in [63,64,65,16] allows each layer-layer pair to be coupled with a different weight, while A inter defined in Equation (4.2) is unweighted.We will introduce the parameter ω in Equation (4.3) to scale the inter-layer adjacency matrix, which can be interpreted as modeling a constant transfer time between all lines at all stops of the network.Note that all numerical methods would be equally applicable if each inter-layer edge was weighted individually.
The definitions in Equations (4.1) and (4.2) now allow us to define the supraadjacency matrix A ∈ R nL×nL of the multiplex network as (4.3) where ω ∈ R ≥0 is a scalar coupling parameter that trades off the relative importance of intra-and inter-layer weights.Note that by our choice of undirected edges discussed earlier in this section we have A = A T throughout this paper.We visualize the above definitions by a small example multiplex network in Figure 6.The left plot shows a multiplex network that consists of three layers (two tram lines and one bus line) taken from the full multiplex network of the German city Freiburg.Only the subset of physical nodes served by these lines is included in the figure.The right plot shows the sparsity structure of the corresponding supraadjacency matrix.Both illustrations distinguish the two types of edges by different colors.
We remark that the supra-adjacency matrix A of our public transport networks is characterized by a high degree of sparsity.This is true by the definitions in Equations (4.1) and (4.2), but our application is characterized by additional sparsity due to the facts that each line typically only serves a small subset of stops in the city and that most layers are approximately represented by chain graphs in which the number of edges is of the order of the number of non-isolated node-layer pairs.Table 1 illustrates the most important multiplex network properties of the 36 German cities considered in the section Public transport network orientations.It confirms the statements above: around 95 − 98% of the node-layer pairs in the networks are isolated, i.e., possess no incident edge and the number of intra-layer edges is close to the number of non-isolated node-layer pairs for all networks.
Finally, we comment on our approach to assigning weights to the edges of the multiplex networks.Many weighted urban network models rely on choosing either the geographical distance [27,26,57] or travel times [5,11] as weights, which is a sensible and straightforward choice for many problems like, e.g., routing algorithms.We, however, aim at applying matrix function-based centrality measures to the multiplex networks in which a node-layer pair is considered central if it is connected to many other node-layer pairs by short walks along edges with large weights.To decide what it means to be "close" to many node-layer pairs in a public transport network, we rely on the combination of two concepts: the frequency of the public transport line [5] and a Gaussian kernel applied to travel times.Gaussian kernels have become a well-established choice of similarity measures in many data-driven applications [59,68,61,17].
We define the frequency φ (l) ij ∈ N as the number of connections that line l offers between stops i and j per day.Furthermore, we denote the travel time between nodelayer pairs (i, l) and (j, l) with i, j ∈ {1, 2, . . ., n} and l ∈ {1, 2, . . ., L} by (∆t) ij leads to most weights being close to zero and σ ij leads to most weights being close to one.As this parameter scales the travel times of public transport lines in the urban network it can be interpreted as a "normalizing travel time".Numerical experiments on the role of the parameter σ are presented in the section Influence of the normalizing travel time.As described at the beginning of this section, edge weights are set to zero if the corresponding layer offers no direct connection.Note, that information about both line frequencies and travel times can be extracted from GTFS data.
We allow different combinations of both weighting concepts by defining the intralayer weights as ij , if (i, l) and (j, l) are connected, 0, ij equal to one of the following expressions: • 1, using no travel times and no frequencies, using only frequencies, , using only travel times, , using both travel times and frequencies.
With the fourth choice of weights we combine travel times and frequencies by multiplication of the two individual weights.This leads to a (linearly) proportional dependence of the weight on the frequency and a (non-linearly) anti-proportional dependence of the weight on the travel time.With this weighting approach, stops connected to many other stops via highly frequented routes with short travel times will be recognized as most central.We compare the different weighting approaches in the section Comparison of different weight models.As described in the section Data, files containing information on transfer options and transfer times between different public transport lines are not required but only optional in GTFS data.Hence, we do not possess knowledge of realistic transfer times, which are represented by inter-layer edge weights in our multiplex network formulation.We thus propose to utilize the unweighted inter-layer adjacency matrix defined in Equation (4.2) and to include the cost of changing lines into the coupling parameter ω from Equation (4.3) by defining a transfer time ∆t transfer ∈ R ≥0 , which we assume to be constant across all pairs of lines and all stops.For consistency, we also apply the same Gaussian kernel to ∆t transfer whenever we include travel times in the weights of intra-layer edges.Thus, the inter-layer weights are given by (4.5) Future public transport models should include realistic individual transfer times for all stops and all combinations of lines serving these stops.From an implementation point of view these transfer times would be easy to include in our multiplex model without causing any issues for our algorithms.We envision that more detailed knowledge of transfer times (including distances and walking times between different lines at the same stop) of practitioners can lead to more advanced city-specific multiplex public transport models in the future.
In this section, we consider matrix function-based centrality measures [13,31,32,14,45,16], which interpolate between local degree centrality and global eigenvector centrality [15].An advantage of this class of centrality measures over eigenvector centrality is its more general applicability. 12.Recently, the class of matrix functionbased centrality measures has been generalized to layer-coupled multiplex networks [16].We illustrate that the same methods are also applicable to the more general multiplex network framework introduced in the section Multiplex network model.In the remainder of this section, we briefly motivate and introduce matrix function-based centrality measures and present efficient numerical methods for their computation as well as numerical results including a discussion of the impact of certain hyperparameters.
5.1.Definition.Matrix function-based centrality measures are based on the application of the matrix exponential or the matrix resolvent function to the adjacency 12 Eigenvector centrality is defined by the eigenvector belonging to the largest eigenvalue of a suitable matrix like, e.g., the graph's supra-adjacency matrix, cf.[63,64,65].The unique existence of this eigenvector can, however, only be guaranteed if the assumptions of the Perron-Frobenius theorem are satisfied [65,Thm. 3.7].This restriction does not apply to matrix function-based centrality measures, cf.e.g.[16,Sec. 6.3] for an example.Note, however, that variants of eigenvector centrality circumventing this shortcoming exist, cf.e.g.[67].matrix of a network.We briefly motivate the definitions by considering walks on single-layer networks and refer the reader to [16,Sec. 3] for more details.At the end of this subsection, we comment on the generalization of the definitions to the case of the multiplex networks introduced in the section Multiplex network model.
In the case of an undirected and unweighted single-layer graph with n nodes, the entry [A] ij of the adjacency matrix A ∈ R n×n is 1 if an edge is present between nodes i and j and 0 otherwise.Note that in the formalism introduced in the section Multiplex network model this corresponds to L = 1 and A = A intra = A (1) .It is well-known from graph theory that an entry [A p ] ij of the pth power of such an adjacency matrix contains the number of walks of length p that exist between nodes i and j [30].While (local) degree centrality can be formulated in terms of the (first power of the) adjacency matrix and (global) eigenvector centrality in terms of the limit of adjacency matrix powers, matrix function-based centralities consider walks of all lengths.This approach can be represented by considering the adjacency matrix power series ∞ p=0 A p .Furthermore, one typically introduces a damping factor in this power series to control the magnitude of the entries of the matrix powers, which typically grows rapidly with p.The following two choices of damping factors lead to the power series of the matrix exponential (left) and the matrix resolvent function (right) where I ∈ R n×n denotes the identity matrix and β ∈ R >0 and 0 < α < 1/λ max 13 denote scalar hyper-parameters.These parameters control the trade-off of locality and globality in the considered walks on the network by assigning more or less weight to longer walks.The centrality of any given node i in the network is then revealed by the following two expressions of these matrix functions: e T i f (A)e i denotes the diagonal element of f (A), which contains the weighted sum of walks of all lengths, which start and end at node i; e T i f (A)1 denotes the sum of the ith row of f (A), which contains the weighted sum of all walks starting at node i (regardless of where they end).Here, e i = [0, . . ., 0, 1, 0, . . ., 0] T ∈ R n denotes the ith unit vector and 1 = [1, 1, . . ., 1] T ∈ R n the vector of all ones.This leads to the definitions of subgraph centrality [32] (left) and resolvent-based subgraph centrality [13] as well as total communicability [14] (left) and Katz centrality [45] (right) The same definitions hold true for weighted graphs and although extensions to directed networks exist [16], we restrict ourselves to undirected networks in this paper, which leads to symmetric adjacency matrices A T = A.
It has recently been proposed to generalize the above definitions to the case of layer-coupled multiplex networks by replacing the graph adjacency matrix by the supra-adjacency matrix of the corresponding multiplex networks [16,Sec. 4].We will demonstrate that all methods from [16] for the symmetric case A = A T still apply for the more general multiplex networks introduced in the section Multiplex network model.
Note that by the construction of the supra-adjacency matrix defined in Equation (4.3) the above quantities yield centrality values for all node-layer pairs, which allows to identify and rank the most central node-layer pairs of the network.Following [63], we call these values joint centralities and denote the joint centrality of physical node i in layer l by JC(i, l).However, we remark that the numerical methods introduced in the following subsection are unable to compute subgraph and resolventbased subgraph centralities of isolated node-layer pairs, i.e., node-layer pairs without any adjacent edge.We set the centrality values of these node-layer pairs to 1, which is consistent with Equation (5.1), in which the respective rows and columns of the supra-adjacency matrix power series are given by unit vectors.
Finally, we adapt the concept of marginal node (M N C) and marginal layer centralities (M LC) [63].These quantities are defined as and can be used to assess the importance of all physical nodes and layers of the multiplex network, which will play a key role in our interpretation of the results.

Numerical methods.
We now discuss strategies for the numerical evaluation of the matrix function-based centrality measures defined in Equations (5.2) and (5.3).For small network sizes, many software packages provide accurate algorithms for the explicit evaluation of the matrix exponential and the solution of a regular linear system (an efficient numerical implementation of, e.g., (I − αA) −1 1 would solve the linear system (I − αA)x = 1 for the vector x).These methods, however, quickly become infeasible for medium to large network sizes and we briefly present efficient and scalable approximations based on Krylov subspace methods.More details can be found in [13] for single-layer networks and in [16,Sec. 5] for multiplex networks.
As we only consider undirected graphs, i.e., symmetric supra-adjacency matrices in this paper we rely on the symmetric methods introduced in [16, Sec.5], which are based on the symmetric Lanczos method [49,39].In its kth iteration, this method constructs the kth column of an orthogonal basis Q k ∈ R nL×k of the Krylov subspace (5.5) has tridiagonal form.This approximation typically achieves a high accuracy for k nL, which makes the eigendecomposition T k = S k Θ k S T k easy to compute with standard methods.These two matrix factorizations can then be combined to compute total communicability and Katz centrality by evaluating the quantity where f is applied elementwise to the eigenvalues of T k [43].The ith entry of the resulting vector then corresponds to the centrality value of the ith node-layer pair.
For subgraph and resolvent-based subgraph centrality we rely on an elegant relation between the symmetric Lanczos method, orthogonal polynomials, and Gauss quadrature discussed by Golub and Meurant [38,40,36,37], which can be used to compute lower and upper bounds on the sought quantities.The final result of this approach yields a lower Gauss quadrature bound, which reads (5.7) where the computation of T k = S k Θ k S T k relies on the basis of the Krylov subspace from Equation (5.5) with v = e i .We refer to [38,40,36,37] and [16,Sec. 5.2.1] for theoretical background on this approximation as well as details on the construction of Gauss-Radau and Gauss-Lobatto rules, which yield an additional lower as well as two upper bounds on e T i f (A)e i .Note that due to the high degree of sparsity in the supra-adjacency matrices of our multiplex public transport networks, Gauss quadrature rules can only be applied to the small subset of non-isolated node-layer pairs.As described in the previous subsection, the centrality value of the remaining node-layer pairs is set to 1.
All numerical methods introduced in this subsection scale linearly in the number of node-layer pairs under the assumption of sparsity in the supra-adjacency matrix, cf.[16,Sec. 5], and usually obtain highly accurate approximations in only around 10 Lanczos iterations.However, the quantities e T i f (A)e i require the approximation of a separate quantity for each non-isolated node-layer pair, which makes them computationally more demanding than the quantities f (A)1, which only need to be computed once.As discussed in [16,Sec. 6.1], the convergence of Equation (5.7) is usually faster than that of Equation (5.6).A python implementation of the above methods is available at https://github.com/KBergermann/Urban-multiplex-networks.

Results.
In this subsection, we present numerical results of the different matrix function-based centrality measures applied to the multiplex network representation of several German urban public transport networks.For the interpretation of the resulting rankings of node-layer pairs (representing stop-line pairs), we mainly rely on marginal node centralities defined in Equation (5.4).This approach corresponds to the identification and ranking of the most central locations in urban networks, which has been the subject of many earlier studies with different centrality measures [27,26,57,58,66,53,3,28].We compare different weight models and matrix function-based centrality measures, which were introduced in the preceding sections.Furthermore, we study the impact of all involved hyper-parameters and close with a discussion of the obtained results including their relation to public transport orientations.
Previous results from the literature on the application of centrality measures to urban transport networks indicate qualitative differences in the distribution of the most central nodes between the cases of public transport [58,66,28] and street networks [27,26,57,53,3].As street networks are usually modeled as (almost) planar graphs in which the "closeness" of nodes is determined by their geographical distance the distribution of central nodes is usually characterized by a relatively smooth transition from more to less central nodes in terms of this geographical distance.For instance, one often observes approximately circular shapes of descending centrality around the most central node of such a network [27,26,57,53,3].Conversely, in public transport networks the "closeness" of nodes is typically tied to different indicators such as travel times or line frequencies.This results in a less smooth distribution of node centralities with respect to the geographical position of the nodes.For instance, stops, which are geographically very close to the most central stop of a public transport network can be classified as non-central if the two stops are not directly con- nected by any line [58,66,28].This general qualitative behavior is confirmed by the multiplex matrix function-based centrality measures and can be observed throughout the subsequent results.Note that a combination of the different modeling approaches by, e.g., adding car or walking layers to public transport networks could mitigate this effect and constitute an interesting future research direction.
To illustrate this general behavior for our modeling approach we start by considering marginal layer and marginal node subgraph centralities of the German city Halle (Saale) as a first example.Figure 7a illustrates the ranking of lines while Figure 7b illustrates the ranking of stops in the corresponding multiplex public transport network.All marginal centrality plots in this paper employ a heat map color scheme consisting of six categories representing centrality value sub-intervals of equal length.While there is a tendency of geographically central stops to be classified as central public transport stops there are various exceptions in the form of dark blue (least central) stops in or around the city's geographical center, cf. Figure 7b.Instead, the comparison with Figure 7a shows that there is a tendency of central stops to line up along important transport axes of the urban public transport system.It is also interesting to note that the orientations of the most central lines in Figure 7a coincide with the preferential directions of the public transport network identified in Figure 3.

Comparison of different measures.
We continue our discussion with the comparison of the four matrix function-based centrality measures for the multiplex public transport networks defined in Equations (5.2) and (5.3) at the example of marginal node centralities of Cologne.
Figure 8 reflects the typical behavior of matrix function-based centrality measures that rankings produced by the four measures are similar but not equal [13,16].For the multiplex public transport networks considered in this paper one typically observes very similar results for SC and SC res as well as for T C and KC.However, SC and SC res typically show a centrality value distribution that is more uniformly distributed across the six centrality value categories, which results in more stops being classified into the top three categories (marked in red and orange) than it is the case for T C and KC.The difference between the two groups of measures is that SC and SC res are defined by the matrix functions' diagonal entries, while off-diagonal entries of the matrix functions are additionally included in T C and KC.This corresponds to SC and SC res only considering closed walks on the networks whereas T C and KC consider all walks.It thus follows that the inclusion of the off-diagonals amplifies the quantitative ranking of the stops due to a steeper distribution of the off-diagonals in comparison to the diagonal entries of the matrix functions.5.3.2.Comparison of different weight models.We also compare the different weight models discussed in the section Multiplex network model.The modeling approaches include all combinations of in-and excluding travel times and line frequencies in the intra-layer weights w (l) ij as specified in Equation (4.4).Whenever intra-layer travel times are included in w (l) ij we also include the transfer time between lines in the coupling parameter ω as defined in Equation (4.5).
Figure 9 presents marginal node resolvent-based subgraph centralities of Stuttgart in the different scenarios.It shows that the inclusion of line frequencies has a larger impact on the centralities than travel times do.Interestingly, the exclusion of line frequencies has the effect of advantaging stops with large inter-layer degrees, i.e., stops with many transfer options.In the case of Figures 9a and 9c this leads to "Stuttgart Universität" (Stuttgart university) being classified as the most central stop of the city albeit being geographically located near the city limit.Later in the section Impact of transfer times we discuss an interesting property of weight models without frequencies regarding the influence of transfer times on marginal node centralities.
Figures 9b and 9d, which include line frequencies, show what is likely to be a more expectable stop ranking for Stuttgart.Here, the increased frequencies of lines serving the geographical center of the city shift the highest ranked stops to this central region.
Overall, Figure 9 suggests that travel times do not have a large impact on the obtained stop rankings.However, their inclusion adds an important ingredient for a realistic modeling of urban public transport networks as it allows the inclusion of transfer times between lines.
In the following three subsections, we discuss the influence of all involved hyperparameters on the obtained stop rankings.These parameters include the normalizing travel time σ, the matrix function parameters α and β, and the transfer time ∆t transfer .

Influence of the normalizing travel time.
Numerical experiments suggest that the normalizing travel time σ only has a small impact on the obtained rankings in a small parameter range: while increasing σ from 0.2 to 2 (corresponding to 12 seconds and 2 minutes normalizing travel time, respectively) slightly increases the number of highly ranked stops in the city center, a further increase of σ to 20 or 200 minutes has no noticeable effect.

From local to global.
As proven for single-layer networks in [15] and as confirmed empirically for layer-coupled multiplex networks in [16,Sec. 6], all four matrix function-based centrality measures defined in Equations (5.2) and (5.3) contain degree and eigenvector centrality as limit cases of the parameters α or β.In this subsection, we empirically confirm this behavior for marginal node centralities of our more general multiplex network models.
We illustrate this at the example of marginal node Katz centralities of Düsseldorf in Figure 10.As discussed before, the admissible parameter range for α for this measure is given by (0, 1/λ max ), where λ max is the largest eigenvalue of the supraadjacency matrix.In Figure 10a, α = 0.01/λ max is chosen to be close to the lower end of that interval, which corresponds to being close to degree centrality in which only direct neighbors of each node-layer pair are considered.The other extreme is represented by Figure 10c in which α = 0.99/λ max is chosen to be close to the upper end of the admissible interval, which corresponds to being close to eigenvector  (l) ] ij in the case that nodes i and j are connected in layer l.Otherwise, [A (l) ] ij is set to 0 in all cases.The plots in the first column are obtained without frequencies, whereas the plots in the second column are obtained with frequencies.The plots in the first row are obtained without a Gaussian kernel applied to travel times, whereas the plots in the second row are obtained with a Gaussian kernel applied to travel times.Furthermore, the same Gaussian kernel is also applied to inter-layer transfer times in the second row.The parameters are chosen as α = 0.5/λmax, σ = 5, and ∆t transfer = 5.The street network in the background of the plots is created with the OSMnx python package [19].
centrality in which nodes are considered central if their closest neighbors are also central.
It can be seen in Figure 10a that the local case is characterized by a relatively uniform distribution of marginal node centralities into the six centrality categories.Furthermore, relatively central stops are geographically distributed across the whole city.In this situation, each local neighborhood can possess its own central stops, regardless of their relative importance for the whole city.Conversely, eigenvector centrality in Figure 10c considers the stationary distribution of any initial distribution of walkers on the node-layer pairs of the network.This often leads to localization effects [50], which entail a centrality distribution that is almost uniform for all but a few very central stops.The strength of matrix function-based centrality measures now lies in the fact that the choice of α allows to continuously interpolate between these two established concepts of centrality measures.In Figure 10b, the choice of α = 0.75/λ max illustrates this property by being "visually in between" the two extreme cases.5.3.5.Impact of transfer times.In Equation (4.5) we specified our approach to modeling a constant transfer time across all lines and stops of the network by means of the coupling parameter ω, which is defined by a Gaussian kernel applied to the transfer time ∆t transfer .We start by presenting an example of a multiplex network of Chemnitz in which intra-layer edges are weighted with travel times and line frequencies.However, an interesting property more frequently emerges when line frequencies are excluded from the intra-layer weight model.We conclude the results section with a parameter study of the coupling parameter ω in the situation without frequencies.
Figure 11 illustrates marginal node total communicabilities of Chemnitz with the constant transfer time varying between 0 and 15 minutes.Assigning no cost for changing lines in Figure 11a leads to Chemnitz's stop "Zentralhaltestelle" to be ranked as the most central stop of the city by a large margin.The role of this stop is encoded in its name, which literally translates to "central stop".The public transport system of Chemnitz is organized such that most of the city's lines stop at this geographically central stop, i.e., it offers by far the most opportunities to change lines.However, this characteristic forfeits its significance as a cost for changing lines is introduced.Figure 11b shows that many other stops in-and outside of the geographical city center are classified as central when the transfer time is increased to 5 minutes.Figure 11c illustrates that a further increase from 5 to 15 minutes (and above) has no significant additional impact.This behavior, however, emerges only for a relatively small normalizing travel time of σ = 1.For larger parameters σ or other cities with a more uniform distribution of inter-layer degrees this behavior is less pronounced when line frequencies are included in the weight model.
Excluding line frequencies from the weight model leads to a stronger localization of marginal node centralities in the situation ∆t transfer → 0 across different networks and a larger range of σ.The effect of the variation of ∆t transfer and subsequently ω on marginal layer and marginal node centralities is illustrated in Figure 12.Here, we observe an interesting clustering behavior of marginal node centralities in the strong coupling limit in which ω approaches its maximum value of 1 corresponding to ∆t transfer = 0.These clusters are determined by the stops' inter-layer degrees, i.e., by the number of lines stopping at the corresponding stop.This behavior reflects the accessibility of a larger part of the network when a larger number of lines can be used within a constant total travel time (defined as the sum of intra-and inter-layer travel

times).
Another interesting observation in Figure 12 concerns the peak in most marginal layer and marginal node centralities between ω = 0.01 and ω = 0.1.This phenomenon has been similarly encountered in other applications before [16,Sec. 6.3] and is an interesting question for future research.5.3.6.Discussion.The preceding subsections gave a detailed account of the differences between the introduced centrality measures and weight models as well as the influence of all involved hyper-parameters.Earlier in this paper, we commented on conclusions obtained from public transport network orientations and their relation to the corresponding street network orientations.We now discuss the relation between these two aspects of urban public transport networks.
Figure 13 provides an overview of marginal node Katz centralities for all 36 German cities from Figure 3 using a common set of hyper-parameters, which allows a direct comparison of the results.While we identified a clear pattern in the public transport network orientations across all German cities the centrality plots in Figure 13 give a mixed picture.Here, we observe various different distributions of central stops both in terms of the relative number of central stops per city as well as their geographical arrangement.
In Figure 14 we focus on four example cities with rather different centrality distributions, which are not necessarily reflected in the corresponding public transport orientation plots.Figure 14a shows that Bielefeld has only a small number of central stops located at the geographical city center.This center is the origin of various sequences of light blue stops (indicating the second to last central category), which seem to equally spread out into all directions and not only into the preferential directions of the corresponding orientation plot.Figure 14b reveals that several central stops are spread out over the city area of Munich.However, the connections between these central stops can not clearly be linked to the preferential direction of the corresponding orientation plot.In Figures 14c and 14d and several other examples we observe highly central sequences of stops lining up along major public transport axes.In some cases, the directions of these axes can be linked to the most pronounced bearings of the corresponding orientation plots, cf.Karlsruhe in Figure 14d but also Braunschweig, Erfurt, Halle (Saale), Hanover, Münster, and Rostock.This pattern, however, can not be observed in other examples like Duisburg in Figure 14c, Cologne, or Stuttgart.
We conclude that orientation plots of street and public transport networks are capable of revealing interesting high-level properties of urban systems.For a more detailed analysis of the geometrical properties of cities like, e.g., the spatial distribution of major transport axes additional measures must be taken into account.This manuscript proposes multiplex matrix function-based centralities as one such measure.

Conclusion.
We studied two geometrical aspects of urban public transport networks: orientations and centralities.We determined orientations of directed public transport networks of the 36 largest German as well as 18 major European cities and compared them to orientations of undirected street networks.All considered German cities revealed two more or less pronounced orthogonal preferential street network directions, which can often be linked to geographical constraints.We found that most considered German public transport network orientations concentrate around the one of the two preferential street network directions, which is closer to the cardinal east-west axis.The same qualitative behavior could only be observed for a small subset of the considered European cities.However, north-south-like preferential public transport directions remained rare.
Furthermore, we formally introduced urban public transport multiplex networks in which nodes correspond to stops and layers to lines and applied multiplex matrix function-based centrality measures in order to identify and rank the most central lines and stops of the considered German cities.These measures generate rankings, which are consistent with previous investigations.In addition, they offer the benefit of being able to flexibly choose the desired degree of locality, possessing efficient and scalable numerical implementations, and being applicable to a wide range of problems.The influence of different hyper-parameters, which all have meaningful interpretations in terms of the urban science application, was thoroughly studied.Our study showed that matrix function-based centralities are capable of revealing insights into geometrical aspects of urban systems on a more granular level than orientation plots are.
We believe that multiplex models of cities, ideally combining various aspects like, e.g., various modes of transportation or additional aspects of urban life can contribute to a better understanding of urban systems.We hope that the presented methodology can add to urban scientists' toolkits and that city-specific modeling as well as parameter-tuning of domain experts yield further contributions to the economic, environmental, and social challenges lying ahead.
Fig.2.Street network orientations of the 36 largest German cities with valid GTFS data as described in the section Data.The plots were created using the OSMnx python package[19].

Fig. 3 .
Fig.3.Public transport network orientations of the 36 largest German cities with valid GTFS data as described in the section Data.

Fig. 6 .
Fig.6.Example public transport multiplex network with L = 3 selected layers of Freiburg (corresponding to two tram layers and one bus layer) and n = 126 stops.The graphics were created with the python package pymnet[46].Left: Visualization of the multiplex network.Intra-layer edges are marked dark blue and inter-layer edges are marked red.Right: Sparsity structure of the corresponding supra-adjacency matrix.Dark blue entries represent intra-layer edges and red entries represent inter-layer edges.

R
≥0 and use a Gaussian kernel to define the similarity of them as exp − (∆t) 1].The scalar parameter σ ∈ R >0 determines the distribution of these similarities, e.g., σ

Fig. 7 .
Fig. 7. Example of marginal layer and marginal node subgraph centralities of Halle (Saale).Central lines and stops are marked red and non-central lines and stops are marked blue.The entries in the supra-adjacency matrix are weighted with travel times and frequencies.The parameters are chosen β = 0.5/λmax, σ = 5, and ∆t transfer = 5.The street network in the background of the plots is created with the OSMnx python package [19].Left: Marginal layer centralities (corresponding to a ranking of the lines).Right: Marginal node centralities (corresponding to a ranking of the stops).

Fig. 8 .
Fig. 8.Comparison of the four different matrix function-based centralities at the example of marginal node centralities of Cologne.Central stops are marked red and non-central stops are marked blue.The entries in the supra-adjacency matrix are weighted with travel times and frequencies.The parameters are set to α = β = 0.5/λmax, σ = 5, and ∆t transfer = 5.The street network in the background of the plots is created with the OSMnx python package [19].

Fig. 9 .
Fig. 9. Comparison of different weight models at the example of marginal node resolvent-based subgraph centralities of Stuttgart.Central stops are marked red and non-central stops are marked blue.The subcaptions denote the respective intra-layer weights [A(l) ] ij in the case that nodes i and j are connected in layer l.Otherwise, [A(l) ] ij is set to 0 in all cases.The plots in the first column are obtained without frequencies, whereas the plots in the second column are obtained with frequencies.The plots in the first row are obtained without a Gaussian kernel applied to travel times, whereas the plots in the second row are obtained with a Gaussian kernel applied to travel times.Furthermore, the same Gaussian kernel is also applied to inter-layer transfer times in the second row.The parameters are chosen as α = 0.5/λmax, σ = 5, and ∆t transfer = 5.The street network in the background of the plots is created with the OSMnx python package[19].

Fig. 10 .
Fig. 10.Marginal node Katz centrality plots of Düsseldorf with varying parameter α.The left figure corresponds to a local measure close to degree centrality, the right figure corresponds to a global measure close to eigenvector centrality, and the center figure corresponds to an interpolated intermediate situation.Central stops are marked red and non-central stops are marked blue.The entries in the supra-adjacency matrix are weighted with travel times and frequencies.The remaining hyper-parameters are chosen as σ = 5 and ∆t transfer = 15.The street network in the background of the plots is created with the OSMnx python package [19].

15 Fig. 11 .Fig. 12 .
Fig. 11.Marginal node total communicability plots of Chemnitz with varying transfer time ∆t transfer .Central stops are marked red and non-central stops are marked blue.The entries in the supra-adjacency matrix are weighted with travel times and frequencies.The remaining hyperparameters are chosen as β = 0.5/λmax and σ = 1.The street network in the background of the plots is created with the OSMnx python package [19].

Fig. 13 .
Fig. 13.Marginal node Katz centrality plots of the 36 largest German cities with valid GTFS data.Central stops are marked red and non-central stops are marked blue.The entries in all supraadjacency matrices are weighted with travel times and frequencies.The hyper-parameters are chosen as α = 0.5/λmax, σ = 5, and ∆t transfer = 5 for all cities.The street networks in the background of the plots are created with the OSMnx python package [19].

Fig. 14 .
Fig. 14.Comparison of selected German cities' margial node Katz centralities and public transport orientations.The entries in all supra-adjacency matrices are weighted with travel times and frequencies.The hyper-parameters are chosen as α = 0.5/λmax, σ = 5, and ∆t transfer = 5 for the centrality plots of all cities.The street networks in the background of the centrality plots are created with the OSMnx python package [19].