First-passage times to quantify and compare structural correlations and heterogeneity in complex systems

Virtually all the emergent properties of a complex system are rooted in the non-homogeneous nature of the behaviours of its elements and of the interactions among them. However, the fact that heterogeneity and correlations can appear simultaneously at local, mesoscopic, and global scales, is a concrete challenge for any systematic approach to quantify them in systems of different types. We develop here a scalable and non-parametric framework to characterise the presence of heterogeneity and correlations in a complex system, based on the statistics of random walks over the underlying network of interactions among its units. In particular, we focus on normalised mean first passage times between meaningful pre-assigned classes of nodes, and we showcase a variety of their potential applications. We found that the proposed framework is able to characterise polarisation in voting systems, including the UK Brexit referendum and the roll-call votes in the US Congress. Moreover, the distributions of class mean first passage times can help identifying the key players responsible for the spread of a disease in a social system, and comparing the spatial segregation of US cities, revealing the central role of urban mobility in shaping the incidence of socio-economic inequalities.

The elements of a variety of complex systems can be naturally associated to one of a small number of classes or categories.Typical examples include the organisation to which an individual belongs [1], the political party of a voter [2], the income level of a household [3], or the functional group of a neuron [4].Quite often, the co-existence of nodes belonging to different classes and the interactions among those classes play a fundamental role in the functioning of a system.For instance, economic and ethnic segregation in cities is known to be associated to the emergence of social inequalities [5,6].Similarly, the organisation of neural cells in the brain and the intricate patterns of relations among different functional areas are known to be responsible for the large variety of cognitive tasks that we as humans are able to perform [7,8].However, obtaining a robust and non-parametric quantification of the heterogeneity of class distributions, especially in systems consisting of a large number of interconnected discrete components, is still an outstanding problem.
We propose here a principled methodology to quantify the presence of correlations and heterogeneity in the distribution of classes or categories in a complex system.The method is based on the statistics of passage times of a uniform random walk on the graph of interconnections among the units of the system.For instance, in the case of a social system we can construct a graph among individuals based on the observed relations or contacts among them.Similarly, in a urban system we can consider the network of adjacency between census tracts or the connections among census tracts due to human mobility flows.The method we propose moves from the classical research on mean first passage times between pairs of nodes in a graph [9][10][11][12], and focuses instead on the distribution of Class Mean First Passage Times (CMFPT) i.e., the expected number of steps needed to a random walker to visit for the first time a node of a certain class β when it starts from a node of class α.By normalising class mean first passage times with respect to a null-model where classes are reassigned to nodes uniformly at random, we can effectively quantify and compare the heterogeneity of class distributions in systems of different nature, size, and shape.
We first test the framework on a variety of systems with simple geometries and ad-hoc class assignments, and then we use it in three real-world scenarios, namely the quantification of polarisation in the Brexit referendum and in the US Congress since 1926, the role of face-to-face interactions among individuals in the spread of an epidemics, and the relation between economic segregation and prevalence of crime in the 53 US cities with more than one million citizens.

I. MODEL
Let us consider a graph G(V, E) with |V | = K edges on |V | = N nodes with adjacency matrix A = {a ij }, and a given colouring function f : V → C, C = {1, 2, . . ., C}, that associates each node i of G to a discrete label f (i) = c i .Let us also consider a random walk on G, defined by the transition matrix Π = {π ij } where π ij is the probability that a walker at node i jumps from node i to node j in one step (see Fig. 1(a)).In general, Π could be any row-stochastic transition matrix, but in the following we will consider only uniform random walks.We are interested here in the statistical properties of the symbolic dynamics W i = {c i0 , c i1 , . ..} of node labels or classes visited by the walk W when starting from i.This dynamics obviously contains information about the existence of correlations and heterogeneity in the distribution of colours.
In the regular square lattice shown in Fig. 1(b), for example, a function f associates 5 distinct colours to nodes of the graph uniformly at random, while in Fig. 1(c), instead, the colouring function induces a small number of homogeneous clusters.We expect that, for long-enough times, all the trajectories of random walks starting from each of the N nodes of the graph in panel (b) will be associated with the same symbolic dynamics, thus becoming indistinguishable.Indeed, that system has neither structural heterogeneity, since all the nodes have the same degree (except for the nodes at the border of the grid), nor inhomogeneity in class assignments, since the probability for a node to belong to a certain class does not depend on its position in the graph or on the classes of its neighbours.In particular, a random walker starting at any node belonging to, say, the light-blue class (see nodes 1 and 2 in panel (b)) will require on average a small number of steps to hit a node in the orange cluster, since orange nodes can be found in the vicinity of any blue node.
If the colouring function induces compact clusters of nodes belonging to the same class, as in the lattice shown in Fig. 1(c), then the statistical properties of the symbolic dynamics W i will heavily depend on the starting node i, despite the fact that almost all the nodes have identical degree.In particular, a random walker starting in the blue cluster at the top-left corner of the graph (node 3) will in general require a very large number of steps before hitting for the first time an orange node.Conversely, a random walker starting at node 4 will hit a node in the orange cluster in a much smaller number of steps, just because one of its immediate neighbours is indeed orange.

A. Class Mean First Passage Time
A quantity of interest in the study of a symbolic dynamics over graphs is the expected time needed to hit a node of a given colour α for the first time, usually know in the literature as the hitting time or mean first passage time (MFPT) to class α [13].We denote as T i,α the average mean first passage time from node i to any node of class α, i.e., the expected number of steps needed to a walk starting on i to visit for the first time any node j such that f (j) = α.Following the formalism to derive the mean first passage times between nodes in a graph, we can write a set of self-consistent equations for T i,α [14]: We can derive an analytic solution for Eq. ( 1) which depends only on the structure of the graph and on the colouring function f .Let us denote as T α the column vector of hitting times from nodes of class β = α to nodes of class α.By convention we set The self-consistent equation for hitting times can be written as: where Π α is the transition matrix of the walk where all the rows and columns corresponding to nodes of class α are set to zero.We denote by D α the indicator vector of nodes belonging to class α, and by D α = 1 N − D α the indicator vector of nodes not belonging to class α, i.e., {D α } i = 1 − δ f (i),α .This leads to the solution: The average Class MFPT τ αβ from class α to class β is computed as: where N α is the number of nodes belonging to class β.
The return time to class α, that is the expected number of steps needed to a walker starting on a node of class α to hit a node of class α (including its starting point) can be computed in a similar way.The forward equation for the hitting time to class α from a node of class α reads: where the first contribution accounts for the neighbours of node i which actually belong to class α, while the second contribution corresponds to walks passing through immediate neighbours of i not belonging to class α.The equation can be written in a compact form as: where R α is the vector of return times to class α, such that {R α } i = 0 if f (i) = α, and T α is the vector of MFPT to class α from nodes that do not belong to class α, as above.Here we denote by Π αα is the transition matrix of the walk restricted to nodes of class α, i.e., whose generic element π ij is set to 0 if either i or j do not belong to class α.Similarly, Π αα is the transition matrix restricted to links from nodes of class α to nodes not in class α.By solving Eq. ( 2) and Eq. ( 5) for the grid lattice shown in Fig. 1(b) we obtain the distribution of CMFPTs: Notice that, as expected, the CMFPTs in the lattice with compact clusters are in general much higher than those observed in the same lattice with uniformly random class assignments.Moreover, both cases we have in general τ αβ = τ βα , since CMFPTs depend primarily on the shape and size of clusters and on the actual fine-grain arrangement of colours.

B. Mean-Field approximation
The expressions for class mean first passage time and class return time provided in Eq. ( 2) and Eq. ( 5) are exact, but they have the drawback of being computationally intensive for graphs with a large number of nodes.It is possible to construct C-class mean-field approximations of these expressions, by representing the behaviour of all the nodes of a class with a single node, and looking at the graph of node classes.If we denote by π αβ the probability for a walker to jump in one step from any node of class α to any node of class β, e.g.expressed as the total fraction of edges from nodes of class α to nodes of class β, the general Mean-Field equation for CMFPT reads: which can be written in compact form as: where by definition {T MF α } α = 0 and Π α is the transition matrix where the row and column corresponding to class α are set equal to zero as above.The solution to this equation is: Notice that this equation is formally identical to Eq. ( 2), with the only difference that it deals with mean first passage times from all other classes to class α, while Eq. ( 2) provides the mean first passage times from all the nodes in other classes to class α.However, the mean-field approximation only retains information about the average probability of jumping between two classes in one time-step, which indeed discards all the (possibly relevant) structural information contained in the underlying graph.For instance, if we use the mean-field approximation for the lattice with random colour assignments shown in Fig. 1(b), we get the solutions: where the values of CMFPT between different classes are substantially different from those computed above using Eq.(2).

C. Normalisation of CMFPT distributions
As we will see in the following sections, the statistics of CMFPT depend substantially on the shape, size, and or- ganisation of class assignments.To allow a fair comparison of CMFPT in systems with different sizes and shapes, we compute the normalised inter-Class Mean First Passage Time τ αβ = τ αβ /τ null αβ .This quantity is the ratio between the expected number of time steps needed to a walk that starts on a node of class α to reach a node of class β for the first time, and the corresponding value in a null-model where classes are assigned to nodes uniformly at random, by preserving their relative abundance.Notice that if classes are distributed uniformly across the system, then τ αβ ≈ 1, ∀α, β, while in the presence of spatial correlations τ αβ will deviate from 1.By using τ αβ we properly take into account several common confounding factors, including an uneven abundance of colours, differences in size and shape, and the effect of borders, thus making it possible to compare different systems on common grounds.Depending on the size of the system, the computation of τ null αβ is based on the average over 10 3 −10 5 independent assignments of classes to the nodes of the original graph.

D. Computation of CMFPT distributions in real-world networks
In the following we will compute and use the distribution of CMFPT in a variety of synthetic and real-world networks.If a system is naturally represented by an unweighted network, the elements of the transition matrix of the random walker will be set to π ij = aij ki , where a ij = 1 is there is an edge from node i to node j, and k i = j a ij is the (out-) degree of node i.If the network is weighted, instead, the elements of the transition matrix will be π ij = wij si , where w ij is the weight of the edge connecting node i to node j, and s i = j w ij is the (out-) strength of node i.For networks with up to N ∼ 10 4 nodes, we solved Eq.(2) and Eq.( 5) exactly, using standard linear algebra packages.For larger graphs we reverted instead to Monte-Carlo simulations, where we estimated the value of T i,α for each node i of the graph as the average over 10 4 ∼ 10 6 random walks originating at that node, and then obtained τ αβ using Eq.(3).

II. SIMPLE GEOMETRIES
We compute here the distribution of CMFPT for some simple geometries with simple class assignments.These examples aim at showing that CMFPT depends heavily on class assignment and on the way nodes belonging to the same classes are arranged.In all these cases, the symmetric nature of colour assignments will allow us to perform the computations on a minimal weighted graph whose distribution of CMFPT is identical to that of the original system.Notice that those minimal weighted graphs provide exact solutions for the original colour assignment they represent, and should not be confused with mean-field approximations.
The first example is that of a 2-dimensional square lattice with periodic boundary conditions (a torus), whose nodes are organised in alternate stripes of black and white nodes, as shown in Fig. 2(a).Thanks to the symmetric nature of this specific colour assignment, a uniform random walk on that graph is effectively equivalent to a random walk on the weighted minimal two-node graph shown in the lower panel of Fig. 2(a), with the transition matrix: The hitting time to the black node in the minimal equivalent graph when starting from the white node can be written as: which gives T •,• = 2 and, by symmetry, also T •,• = 2.
As a second example we consider a finite chain of white nodes surrounded by black nodes, as shown in the top panel of Fig 2(b).We are interested here in showing how the length L of a linear cluster of a given colour influences the distribution of CMFPT across the cluster, so we will focus on the CMFPT from nodes of class • to nodes of class •.In this case the system has a mirror symmetry, which effectively allows us to focus on the L 2 on either side of the chain.The only caveat is that weighted minimal graphs associated to chains with even length L are slightly different from those associated to chains with odd length, as shown in Fig. 2(b).For L ≥ 6, the closed expression for the distribution of CMFPT from each node in the chain is: , and A k and B k are two rational sequences whose form depends on whether the length of the chain is even or odd.More details about the derivation are provided in Appendix A 1, where we also compute the distributions for 1 ≤ L < 6.In Fig. 2(d) we report the distribution of T k,• across the chain, while in Fig. 2(e) we show the average CMFPT to nodes of class •.Notice that when L 1, T • converges towards 2, which is the same value obtained in the case of a lattice with rows of alternate colours seen above, as expected.
As a final example we consider the infinite cylinder shown in the top panel of Fig. 2(c), where the nodes in the upper M + 1 rows are of class • and those in the bottom row are of class •.The aim of this example is to show the behaviour of CMFPT as a cluster becomes deeper, i.e., as the nodes in the cluster of a certain class are placed farther away from the frontier with the other cluster.For our purposes, this geometry is effectively equivalent to the linear chain of nodes shown in the bottom half of Fig 2(c), where each row of the original graph is represented by a single node in the minimal graph, with an appropriately-weighted self-loop.We can write a set of self-consistent equations for the hitting time to class • for a walk started on each of the rows k = 0, 1, . . ., M : whose solution is: and The full derivation is reported in Appendix A 2. In Fig. 2(f) we show the scaling of average mean first passage time to class •, defined as: 8M 2 + 33M (11) where the approximation is accurate for M 1.This means that the average MFPT to class • scales as the square of the height of the cylinder.In the same panel we also show the scaling of the second moment of T k,• : ) which indicates that the standard deviation of T k,• is a function that grows as O(M 2 ) as well.In other words, the deeper a cluster the much higher the CMFPT to its border, and the much wider the distribution of the MFPT from any single node of the cluster to the border.

III. SYNTHETIC COLOURINGS IN TWO-DIMENSIONAL LATTICES
We study here the distribution of CMFPTs in a twodimensional square lattice with N = L × L nodes.Here each node is associated to one of two possible classes, namely • or •, depending on their position in the lattice.Without lack of generality, we set the relative abundance of • nodes r = N• N .Then we assign N • nodes to class • sampling their coordinates (x, y) from a symmetric twodimensional Gaussian distribution centred in the middle of the lattice, with standard deviation equal to σ: By tuning r and σ this model allows for a continuous transition between homogeneous distributions of colours and patterns with strong class segregation.In Fig. 3(a)-(f) we show sample sketches of how the spatial distribution of colours looks like as a function of σ and r.For very large values of σ, the picture approaches a homogeneous distribution, regardless of the relative abundance of the two colours.As σ decreases, instead, the nodes of class • will become more strongly clustered around the centre.The role of the relative abundance of colours is evident from the comparison of panel (b) and panel (c), which are two configurations with the same σ respectively for r = 0.9 and r = 0.5.In particular, we note that the relative abundance of the two colours also has a non-trivial role in determining the degree of mixing between the two classes, since more • nodes can be found inside the • cluster for r = 0.5 than for r = 0.9.
A summary of the interplay between these two parameters is reported in Fig. 3(g), where we show the quantities τ •→• and τ •→• / τ •→• for a variety of values of r and σ (the points in panel (g) corresponding to the configurations in panels (a)-(f) are labelled accordingly).For large values of σ, we expect a more homogeneous distributions of the two colours, (see Fig. 3 It is worth noting that the ( τ ) phase space provides a very intuitive interpretation and a powerful visualisation of the heterogeneity of distributions of two classes, and it would thus be quite useful to characterise the spatial heterogeneity of generic binary feature distributions and point statistics [15,16].

IV. POLARISATION AND SEGREGATION IN VOTING DYNAMICS
Among the variety of social dynamics that exhibit complex behaviours, voting is possibly one of the most interesting.And not just because anything concerning politics can spur endless and ferocious discussions, but also because voting patterns are the result of the interplay of a variety of factors that are generally difficult to model in an accurate way, including socio-economic and cultural background and spatial and temporal correlations [2,17].We focus here on two examples where voting dynamics result in the emergence of heterogeneity and correlations, namely the spatial clustering of Leave/Remain voters in the Brexit referendum and the polarisation of opinions of roll-call votes in the US Congress.

A. Spatial heterogeneity of Brexit vote
Here we show how inter-class mean first-passage times can be used to quantify the spatial segregation of voting patterns.We consider the results of the so-called "Brexit" referendum, held in the United Kingdom in 2016 to decide whether to leave the European Union.The referendum had a turnout of 72.2%, and 51.9% of the voters expressed the preference to leave the EU.The results of the referendum have been analysed in several works [18] which have outlined interesting correlations of voting preference with a variety of socio-economic indicators, including age, income level, unemployment, and level of education [19].One of the most intriguing aspects of the results was that the vote was highly segregated.Indeed, highly urbanised areas, as well as the majority of constituencies in Scotland, voted preferentially for Remain, while the rest of the country expressed a preference to Leave.
We constructed the planar graph of constituencies in Great Britain (i.e., all the mainland constituencies in England, Wales, and Scotland, leaving out the few constituencies in Northern Ireland), where each node is associated to a constituency and a link between two nodes exists if the corresponding constituencies border each other.We assigned each node to either Leave or Remain according to which party won the majority of votes in the corresponding constituency.We then computed the normalised average CMFPT pattern for Remain (R) and Leave (L), we obtained the values shown in Table I It is worth noting that while the normalised return time to each class is close to 1 in both cases, i.e., it is consistent with the corresponding null model where classes are reassigned to nodes uniformly at random, the proper inter-class MFPTs exhibit a pronounced disparity between the two classes.In particular, the normalised CMFPT from Leave to Remain τ LR is much smaller (2.827) than its counterpart τ RL (12.304).This means that, on average, in this graph is much easier for a random walker starting at a node whose citizens expressed  A possible interpretation of this result is the presence of a structural reinforcement of segregation.Indeed, if we assume that voters could influence each other by discussing the matter of the referendum with other voters holding opposite opinions, then a person determined to vote for Remain would have had a much harder time finding a Leave supporter to convince.On the contrary, Leave supporters would have been able to find Remain supporters much more easily, as Leave constituencies are on average closer (in terms of MFPT) to other Remain constituencies than the other way around.
However, this picture is completely reversed if we take into account the actual percentage of votes for Leave and Remain in each constituency, instead of noting only which party got a majority.We considered an ensemble of colour assignments obtained by assigning colour L to node i with probability p L (i) equal to the percentage of Leave voters in constituency i.For instance, if in a given constituency i we had 45% of votes for Leave, then node i will be assigned to class L with probability 0.45, and to class R with probability 0.55.We computed the inter-class Mean First Passage Times in this ensemble of colourings, and normalised them by the corresponding values in a null-model where we re-assigned vote proportions among nodes uniformly at random.The results are reported in Table I It is evident that, by taking into account the actual distribution of voters in each constituency, the pattern of inter-class MFPT becomes practically indistinguishable from the one we would observe in the corresponding null-model.This means that there was indeed no significative spatial segregation effect in the Brexit vote, and indicates that indeed the reasons in support for Leave or Remain were most probably linked to socio-economic characteristics, rather than to geographical ones.

B. Polarisation of roll-call votes
In this section we show how τ αβ can be used to quantify the level of polarisation in roll-call votes in the US congress [20][21][22], and to keep track of its evolution over time.We considered the full data set of affiliation and single roll-call votes of the members of the US Congress, and we built a weighted graph among members of each chamber in each term between 1929-2016, where the weight of the edge connecting two nodes is equal to the number of times the votes of the corresponding members in a roll-call coincided.We assigned each node to either Republicans, Democrats, or Others, according to the party to which they belonged, but we focused exclusively on the CMFPT between Republicans and Democrats, since members of Other parties are normally a rather small minority, if present at all.
It is worth noting that, at difference with the synthetic networks we have studied so far, these networks do not admit a natural embedding in a metric space.Intuitively, we expect that members of both the Senate and the House of Representatives would, in general, be more likely to vote as other members of their party, giving rise to somehow definite clusters.However, the situation is not always that clear.In Fig. 4(a)-(c) we show the networks of Senate members observed in three different terms from the last 30 years.Indeed, it is evident that stronger connections between members of the same party appear as time passes.Moreover, when comparing the 112th and 115th congresses, we can see that the party with majority tends to be more heavily connected than the other one.This is also reflected in the relative size of nodes, which is proportional to the corresponding CMFPT to the other party.Just by looking at these three graphs we would be inclined to think that the level of polarisation in the US Congress seems to increase over time.
For each graph, we also show the distribution of normalised mean first passage times T i,α from each node i to each of the parties, respectively for Democrats and Republicans.By looking at these distributions on the same scale, it is evident that the polarisation, intended as the relative distance between nodes of different parties as measured by inter-class MFPT, has increased substantially in the last thirty years.In particular, the separation between the intra-and inter-class MFPT distribu-tions has increased dramatically, to the point that in the 115th Congress the distributions of intra-class and interclass MFPT are clearly separated.When comparing the 112th and 115th legislatures we also observe that the changes in the shape and position of the MFPT distributions seem to depend on which party holds the majority of the seats.In particular, when Democrats have the majority, the CMFPT to Republicans across nodes is higher than from Republicans to Democrats.Conversely, when Republicans are ruling the situations is inverted, pointing out that the party that holds the majority seems to be the main driver of polarisation.
In Fig. 4(d)-(e), we show the evolution of the average τ αβ between the two main US parties in both the Senate and the House of Representatives.The dashed lines indicate the intra-class MFPT, and indeed support the intuition that polarisation has increased dramatically in recent years.The plots of return times τ αα instead (solid lines), are far more stable, and the small oscillations we observe depend only on which is the ruling party, since that one will normally yield lower values of τ αα .
We quantify the overall polarisation in the Senate and the US House of Representatives by computing the average between the inter-class passage times ( τ D→R + τ R→D )/2, which provides an estimate of how close are the voting behaviours of the two main parties in each Congress.The temporal evolution of ( τ D→R + τ R→D )/2 is shown in Fig. 4(f), where each point is coloured according to the party holding the majority of seats in that term.Again, we observe a clear increase of polarisation after the 1970s in both branches of the Congress, which is in very good agreement with prior works [20][21][22].Notice that there are some interesting patterns there, like for instance the term starting in January 2001, which displays the lowest polarisation in the last 25 years, most likely due to the 9/11 terrorist attacks later that year.
We inspect next if, as hypothesised, the increasing polarisation observed in Fig. 4(f) is due to the party holding the majority of seats in each term.We computed the ratio of inter-class passage times τ M aj→M in / τ M in→M aj , where τ M aj→M in is the normalised inter-class MFPT from the party with the majority to the party with the minority and τ M in→M aj is the normalised inter-class MFPT from the minority to the majority.The values of τ M aj→M in / τ M in→M aj (Fig. 4(g)) are indeed relatively stable over time, with the vast majority of points lying along or above the solid grid line corresponding to absence of polarisation.This indicates that in most of the terms over the last 80 years the party holding the majority has been the main driver of polarisation.These results demonstrate that τ αβ is indeed a robust measure of polarisation.We argue that the same framework could be easily used in other contexts, including the polarisation of discussions in (online) social networks [23,24] or the flip of candidates between different political parties [25].In particular, it could be also possible to identify those agents or individuals who contribute more to polarisation by looking at the ranking of nodes by their values of CMFPT to other classes.

V. CONTACT ASSORTATIVITY AND RELATION WITH EPIDEMICS SPREADING
Understanding the mixing between groups of individuals in a network can provide a lot of information about the properties of social dynamics, including the role played by different individuals in the transmission of diseases [1,[26][27][28].As a second case study, we use CMF-PTs to identify how different groups of people interact in three face-to-face contact networks, namely the contacts in a hospital, in a school, and in an enterprise, obtained from the SocioPatterns project data set [1,[29][30][31].
The definition of groups or classes is specific of each system, e.g., the role of a person in the case of hospitals, the class in the case of schools, and the department in which a person works for the network of contacts in an enterprise environment.The weight of the undirected edge connecting two nodes in each graph is equal to the number of contacts between the corresponding individuals.By looking at the dynamics of two simple epidemic models, namely a Susceptible-Infected-Susceptible (SIS) and a Susceptible-Infected-Recovered (SIR), we show here that the distribution of CMFPT in each graph provides relevant information about the dynamics of disease spreading in the system.For both epidemic models, if an individual i is infected it selects one of its neighbours j with probability w ij / j w ij and infects it with a probability β.Afterwards, each of the infected individuals will either recover with probability µ in the case of the SIR model, or become susceptible again in the case of the SIS model.
In Fig. 5(a)-(c) we report the matrices of τ αβ respectively for the school, the enterprise, and the hospital.In the school network, we observe a consistent pattern of lower values of normalised intra-CMFPT τ αα , which is most probably due to the much higher number of contacts among individuals in the same classroom compared to individuals in other classrooms.Notwithstanding this general pattern, some of the classes are more tightly connected than other close-by classes, as in the case of ce1b and ce2b.Also in the case of the enterprise network we distinguish a clear pattern where for each class α the value of τ αα is normally much smaller than τ αβ for α = β.Yet, there are some interesting deviations, as in the case of SFLE.This department plays a similar role to that of teachers in the case of the school network (i.e., contain people that tend to interact with more than one class), and all the other departments exhibit similar values of inter-CMFPT to that class.Finally, in the hospital contact network patients seem to be the most isolated, while the paramedical staff and display lower CMFPT to all the other classes.
As we show in the following, we found a quite interesting relation between the pattern of CMFPT in each network and the spreading dynamics on the same graph.We run a large number of simulations of the SIR and SIS models, seeding the disease in each of the nodes of a network, and calculating the number of time-steps needed to the spread to reach the peak in each of the classes, as a function of the group to which the seed node belongs.In Fig. 5(d)-(f) we report, as a function of τ αβ , the average number of steps until the peak of the epidemic t peak in class β is reached in a SIR model where the seed is a node in class α.Interestingly, we found that the time to the peak in class β is an increasing function of τ αβ , and the rank correlation between the two variables is always pretty large and significant in the three systems.
The results on the SIS model somehow complement the picture observed in the case of SIR.In Fig. 5(g)-(i) we show that the class mean return times τ αα are strongly correlated with the fraction ρ α of infected individuals of that class in the stationary state.In particular, the lower the value of τ αα , the larger the fraction of infected individuals in class α in the endemic state, suggesting that the steady-state dynamics is indeed predominantly driven by interactions among individuals in the same class.In Table II we report these correlations for a wider range of β and µ.As shown in detail in Fig. 7 and Fig. 8, the observed correlations between ρ α and τ αβ are consistently higher than the correlation with either the total number of edges between two classes or the total fraction of edges from class α to class β.Those results as well as the correlations with a wider range of parameters and contact networks can be found in Appendix B.

VI. RESIDENTIAL AND DYNAMICAL URBAN ECONOMIC SEGREGATION
Socio-economic segregation has an enormous impact on city livability, and many different measures to quantify it have been proposed in the past years.However, most of those measures, with a few notable exceptions [32], focus on first-neighbour information, and disregard the role of the mobility of citizens [33][34][35][36].We test here the potential of CMFPTs to quantify urban economic segregation in large metropolitan areas, taking into account the daily mobility patterns of individuals.
We considered the 53 US cities with more than one million inhabitants, and census information about the number of households in each of the 16 income categories defined by the US Census Borough (see Table III and the 2017 American Community survey [37]), where class 1 is lowest income and class 16 is the highest one.We constructed two different graphs among census tracts, namely the graph of tract adjacency and the graph of daily workplace commuting [38].The former graph is undirected and unweighted, and is better suited to measure the so-called residential segregation, i.e., the extent to which people with similar levels of income tend to live in close-by areas.The commuting graph, instead, is directed and weighted so that the weight of a link going from i to j is given by the sum of the residents in i working in j and the residents of j working in i in order to mimick the daily mobility of citizens.
The fact that households of more than one category are present in each block group prevents us from assigning a single class to each unit, so we computed the distri-bution of CMFPT by averaging over a large number of realisations of class assignments.In each realisation, the class of each node in the graph is sampled from the distribution of household in the corresponding census tract, so that the probability that node i is assigned to class α in a given realisation is equal to m i,α / ∀β m i,β , where m i,α is the number of people of category α living in node i.We took a slightly different approach for the commuting graph as described in [39], where the population attributed to node i is a combination of the resident population at i and the number of commuters working at i: where ω ji is the weight from node j to node i in the commuting graph, which is equal to the total number of people who live in j and commute to their workplace in i.By doing so we take into account the fact that commuters effectively contribute to the diversity of an area, as they spend a considerable amount of time in there and actually interact with other commuters coming from different areas as well.For each realisation of a class assignment, we run 10 4 − 10 6 random walkers from each of the nodes.Then, we obtain the CMFPT τ αβ for each ordered pairs of classes (α, β) by averaging over all nodes and realisations of class assignments, and we analyse the normalised CMFPT τ αβ = τ αβ /τ null αβ .In Fig. 6(a)-(b) we show the profiles of τ αβ from each of the 16 classes computed over the adjacency graph, respectively for Detroit and Boston.It is worth noting that the categories at the two extremes (i.e., the poorest and the wealthiest ones) exhibit a quite similar pattern in the two cities.They are both characterised by larger values of CMFPT from any of the other classes, meaning that those two classes are in general more isolated from the rest of the population, with most of the highincome classes appearing slightly more isolated in Detroit than in Boston.Moreover, in both cities all classes have a virtually identical value of CMFPT to class 9, and very similar values to class 8 and 10, which indicates that these middle-income classes play a pivotal role in the spatial distribution of income.However, despite the fact that the qualitative behaviour is similar in the two cities, there are some noticeable quantitative differences.First of all, the values of normalised CMFPT are significantly larger in Detroit than in Boston.Second, in Detroit we observe a strong dependence of CMFPT on the class of the starting node, while in Boston the average number of steps required to reach a class below 12 is almost constant regardless of the category of the origin node.These important quantitative differences would suggest that the spatial distribution of income in Detroit is more heterogeneous than in Boston, a conclusion which is in line with the classical literature about economic segregation in the US [3,[40][41][42].
However, in large metropolitan areas most of the daily activities of individuals happen far away from their home, due to urbanisation pressure and to decentralisation of productive sectors, so that residential segregation can hardly tell the whole story.Indeed, in the last years there has been an increasing interest in the quantification of social and economic segregation by taking into account mobility patterns [43][44][45].This is easily doable within our framework by letting the walkers move on the mobility graph instead of the adjacency network between census tracts.We have computed τ αβ upon the mobility graph of each US city in the data set, as obtained from workplace commuting information.The results are shown in Fig. 6(c)-(d), and provide an interesting picture of the differences between residential and mobility-focused segregation.First of all, in both cities τ αβ does not depend too much on the origin class α but instead on the destination class β.This is most likely due to the fact that people of different backgrounds commute to similar areas -i.e. the city centre and industrial sites-.Still, both cities display a distinct organisation of CMFPTs.In Detroit low income classes are much more isolated (higher τ αβ ) compared to high income classes, which exhibit systematically lower values of τ αβ .In the case of Boston, instead, τ αβ is almost flat with no important dependence on either the source or the destination class.
The profiles of τ αβ that we show in Fig. 6(a)-(d) provide an overall clear picture of the distribution of CMFPT in a metropolitan area, but do not allow to easily compare two cities in a systematic manner.Hence, we devised two synthetic indices ξ out and ξ in that summarise the information on spatial income heterogeneity in a single number.The idea behind these quantities is that an income class α is more heterogeneously distributed if there is a large difference between the CMFPT from α to the income classes immediately adjacent to α and the median CMFPT from α to any other class.We consider the discrepancy between the local and global median of CMFPT from class α: where τ out α is the median of τ αβ when α = β and τ out αnn is the median of τ αβ to its nearest neighbours β ∈ {α − 1, α, α + 1} [46] Similarly for the discrepancy between local and global median of CMFPT to class α: where τ in α and τ in αnn are now the median of τ βα when α = β and the median of τ βα from its nearest neighbours β ∈ {α − 1, α, α + 1}, respectively.
In Fig. 6(e) we show the ranking of US cities induced by ξ = (ξ in + ξ out )/2, that is the average discrepancy between local and global CMFPT from/to each income class.In general, larger values of ξ indicate more pronounced levels of segregation.On the left-hand side of the panel the cities are ranked according to ξ in the adjacency graph of census tracts, while on the righthand side the ranking is based on ξ in the commuting graph.Interestingly, Detroit is the first US city by residential segregation, with other cities traditionally known for their high levels of segregation like Milwaukee and Cleveland following closely.Conversely, Boston is at the bottom of the ranking.However, the ranking changes substantially if we consider instead the mobility graph, and what we call dynamic segregation, as reported in the right-hand side of Fig. 6(e).For instance, Baltimore (which is ranked pretty high for residential segregation) gets relegated to a mid-rank position, while cities like Buffalo or Indianapolis, where residential segregation is not that high, get to the top of the ranking of dynamic segregation.
The most interesting aspect of ξ is that it captures some of the most undesirable consequences of income segregation, namely the incidence of different types of crimes obtained from [47].In Fig. 6(f)-(h) we report the correlation between incidence per-capita of assaults, violent crimes, and robberies with the levels of residential segregation measured on the adjacency and on the commuting graph of the cities in our data set.Interestingly, both indices display a significant correlation with all three types of crime.Moreover, the indices computed over the commuting network display a stronger correlation in all three cases, reinforcing the idea that quantifying segregation by disregarding mobility can indeed lead to distorted conclusions.For instance, the city that appears on the top of the dynamical segregation ranking (Memphis) is the one displaying the highest incidence of violent crimes percapita among the largest US cities, although it is placed in the second quartile of the ranking by residential seg-regation.As a comparison, we report in Appendix C the correlations of traditional metrics of segregation (i.e., the Spatial Gini coefficient and the Moran's I index) with the same crime indicators, showing that they are not able to attain values as high as those obtained with ξ .Besides the lower correlations, it is also important to note that both the Spatial Gini Coefficient and the Moran's I index are symmetric quantities, and as such they fail to capture the intrinsic asymmetry between income classes revealed by Fig. 6(a)-(d).
Overall, our methodology based on the diffusion of random walks is not only a natural extension of the latter multi-scalar approaches introduced to characterise residential segregation [48], but it also allows us to define a dynamical segregation that includes mobility into the analysis as it has been recently discussed for instance in Refs.[45,49].

VII. CONCLUSIONS
Although being a relatively simple and intuitive tool, random walks have been extensively and successfully used to model and characterise transportation [50,51], biological [52,53], and financial systems [54].They have proven useful to detect meaningful structural properties in complex networks [14,55], including communities [56][57][58], node roles [59], navigability [60], and temporal variability [61].Moreover, interesting insights about the relation between the structure and dynamics on a network have come from the analysis of transient and long trajectories of random walks on graphs, including the statistics of first passage times and coverage times [12,13,[62][63][64], or the systematic study of their fluctuations [65].However, the potential usefulness of random walks to quantify the heterogeneity of class distributions on networks, either in spatially-embedded systems or in high-dimensional networks, has only recently been hinted to [39,65,66].
We have shown here that the information captured by the distribution of inter-class mean first passage times can be used not just as a way to detect the presence of anisotropy and correlations in the properties of nodes, but also as a reliable proxy for the dynamics and emergent behaviours of a complex system.One of the most interesting aspects of the measures of heterogeneity, polarisation, and segregation that we have introduced in this work is that they take into account microscopic, meso-scopic and global relations among classes, due to the fact that in principle random walks integrate information about paths of all possible lengths.Another relevant property of the measures of segregation based on CMFPT is that they are non-parametric and correctly normalised with respect to a meaningful null-model, hence allowing us to compare on equal grounds the heterogeneity of class distributions in systems of different sizes, which is where most of the classical indices of segregation fail [34,48].Even more importantly, the profiles of inter-class mean-first passage times are not symmetric with respect to classes, and provide fine-grained information about which classes are most responsible for the emergence of polarisation and heterogeneity.In this respect, it would be worthy exploring how the simple measure of polarisation that we proposed can be extended to the case of more than two classes.
The fact that measures of class heterogeneity based on random walk statistics correlate quite well with some of intrinsic dynamics happening in social network (i.e., the spread of an epidemic) and with some other exogenous processes mediated by the underlying graph (i.e., the incidence of crime in a city) confirm that the profiles of CMFPT are a useful toolbox for targeted mitigation of the undesired effects of these dynamics.For instance, the groups of a social network that are more central according to τ αβ might be the best candidates for early vaccinations to slow-down an epidemic.At the same time, the definition of dynamical segregation based on walks on the mobility graphs, and the fact that it correlates quite substantially with crimes, potentially paves the way for a re-definition of the traditional role attributed to residential segregation, in favour of a more balanced view that takes into account the activity patterns of citizens together with the spatial distribution of their dwellings.
The generality of the methodology proposed in this paper, and its applicability to different classical problems in complexity science, establish a concrete link between classical statistical physics and modern complexity science, and have the potential to provide new interesting insights about the relation between structure and dynamics of complex systems.
Notice that the forward equation for T 0,• is somehow simpler, and reads: which can be rewritten as: If we now consider Eq. (A25) for k = M − 1: and we plug it into Eq.(A26), we get: which can be solved for T 0,• , yielding: By substituting T 0,• back in Eq. (A27), and recursively using the results to obtain T 2,• , T 3,• , . .., it is easy to find that: and .
Appendix B: Supplementary results for the epidemic spreading in contact networks We report here the correlation between some simple structural measures and the same epidemic variables considered in Fig. 5.In particular, we have considered two main quantities, namely, the total number of contacts between individuals of group α and group β and the normalised contacts of the members of group α with those of group β.The first quantity is symmetric and is just the sum of the total number of contacts between the members: while the second one: is normalised by the total number of contacts of the group of origin and, therefore, is in general not symmetric.In Fig. 7 we show the correlation between the total number of steps until the peak and the total number of contacts between each pair of groups, as well as the correlation between the fraction of population in each group infected in the endemic state (SIS model) and the total contacts between members of the same group.Despite correlations are relatively high and significant in some cases, none of them is higher than those obtained for τ αβ .
Similarly, we show in Fig. 8 the correlation of both the SIR and the SIS models with the normalised number of contacts.Correlations look slightly higher than in the case of the total number of contacts, yet they are still considerablu lower than the values obtained for τ αα .Particularly in the case of the SIS model (stationary dynamics), we do not observe the clear alignment that appears in Fig. 5.
In Fig. 5 we reported the correlation coefficients for only 3 of the 7 contact networks we have analysed.Table II includes additional results of Pearson's r p and Spearman's r s for many combinations of the parameters and a for wider variety of contact networks.Interestingly, the results are compatible with those reported in Fig. 5, and confirm that correlations between CMFPT and critical epidemiological variable are large and significative.given by which is generally defined for any type of network.However.there is a slight difference between the adjacency and the commuting graph.In the case of the adjacency graph we have used a row normalised approach for the weights so that each cell counts the same in the global   average, while for the commuting graph no normalisation was applied so that a cell with more (in)outgoing links will count more.In Fig. 9, we provide the values of both indices for each of the income categories analysed throughout this paper.As can been seen, the behaviour observed is quite comparable to the one in Fig. 6, with the low and high-income classes displaying significantly higher segregation.Moreover, segregation seems to decrease when computed over the commuting graph.
To extract a single metric for each of the cities studied, we have taken the median over all income categories and we have compared the correlations with those observed for the index extracted from CMFPT.We report in Fig. 10 and Fig. 11 the correlations for both quantities, the Spatial Gini and the Moran's I, respectively.While a significant correlation appears in some cases, the highest one observed for the Moran's I calculated on the adjacency graph is still lower than the lowest correlations observed with our index ξ .This confirms the higher explicative power of a metric based on the diffusion on random walks and the crucial role of including mobility  into the analysis.

FIG. 2 .
FIG. 2. (a)-(c): Three simple 2-colour geometries on infinite 2D lattices (top), and the equivalent minimal graphs used to compute inter-class mean first passage times (bottom).When nodes are organised in rows of alternate colours (a) the equivalent graph contains only two nodes with symmetric connections.The equivalent graph of a chain of • nodes surrounded by • nodes (b) is equivalent to a comb graph, while an infinite cylinder of with M + 1 rows of • nodes (c) reduces to the chain with M + 1 • nodes.In the case of the chain graph in (b), the CMFPT of a node as a function of its distance k from the endpoint converges to 2 as k increases (d).Similarly, the average CMFPT to • converges to 2 when the length L of the chain increases (e).In panel (f) we show the first two moments of T k,• in an infinite cylinder, as a function of the height of the cylinder.

FIG. 3 .
FIG. 3. (a)-(f) Snapshots of the model for different values of the parameters σ and ratio.(a) σ = 28 and ratio=0.5(b) σ = 14 and ratio=0.9(c) σ = 14 and ratio=0.5(d) σ = 9 and ratio=0.7 (e) σ = 8 and ratio=0.1 (f) σ = 7 and ratio=0.3.(g) Transition in the coordinates τ•→•/ τ•→• as a function of τ•→• for different values of σ and the ratio between colours.The colour of each point corresponds to the value of σ in the Gaussian first and type of marker to the ratio between colours.Simulations were performed in a regular lattice of size 60 × 60.
(a)), and indeed we have τ •→• 1, meaning that the relative distribution of CMFPT from • nodes is compatible with the one observed in the corresponding null model.At the same time, the ratio τ •→• / τ •→• is close to 1 as well, meaning that the normalised CMFPTs of the two classes are indistinguishable.As σ decreases, the • cluster becomes more prominent, but the actual relation between τ •→• and τ •→• / τ •→• will depend on the value of r.In particular, if r > 0.5, i.e., nodes of class • are the majority (see Fig. 3(b) and Fig. 3(d)), the ratio τ •→• / τ •→• normally remains larger than 1.This is again expected, since a deeper cluster of • nodes causes an increase in the CMFPT from • to • nodes, along the same lines of the increase in CMFPT observed in the simple geometry shown in Fig. 2(c).Conversely, if r < 0.5 then • nodes are interspersed within a large cluster of • nodes (see Fig. 3(e)), making it harder for a walker started at a • node to find a • node.This results in values of τ •→• larger than 1, i.e., much longer than in the corresponding null-model, and in τ •→• / τ •→• < 1.There are some particular situations (see Fig.3(c),(f)) in which despite the increase in τ •→• , the ratio τ •→• / τ •→• remains close to 1.This is due to the fact that in these cases the two classes are distributed in a similar fashion, and the relative depths of the two clusters are indeed comparable.

FIG. 4 .
FIG. 4. Roll-call polarisation in the Senate and the U.S. House of Representatives.(a)-(c) Three networks connecting the members of the senate where edge weights correspond to the number of calls in which the corresponding two members voted the same on a roll call, and the corresponding distribution of CMFPT from the nodes (members) to each of the parties.(a) 101st Congress (1989-1990), (b) 112th Congress (2011-2012) and (c) 115th Congress (2017-2018).(d)-(e) Temporal evolution of τ αβ in (d) the Senate and (e) the U.S. House of Representatives.(f) Temporal evolution of polarisation in the the Senate and the U.S. House of Representatives, as measured by ( τD→R + τR→D)/2.(f) Temporal evolution of τMaj→Min/ τMin→Maj.The fact that this quantity is larger than one in most of the terms indicates that the party with a majority is the one responsible for the largest contribution to polarisation in the corresponding term.

5 FIG. 5 .
FIG. 5. Class mean first passage times and the spread of epidemics in contact networks.(a)-(c) Class mean first passage times in contact networks corresponding to a school, a hospital and an enterprise.(d)-(f) Connection between τ αβ and the time steps until the peak of the epidemic t peak in a SIR model.(g)-(i) Connection between the return times ταα and the fraction of infected individuals in each class in the stationary state.Simulations were performed using the parameters β = 0.7 and µ = 0.1.

FIG. 6 .
FIG. 6. Class mean first passage times and urban income inequality.(a)-(d) Class mean first passage times τ αβ between each of the income categories in Detroit and Boston when computed over the (a,b) adjacency (residential segregation) and the (c,d) commuting graphs (dynamical segregation.(e) Changes in the ranking of values between the residential and dynamical segregation.(f)-(h) Correlation between the index of segregation ξ and the criminality in US cities. (f) Violent crimes, (g) robbery crimes and (h) assault crimes.

FIG. 7 . 60 FIG. 8 .
FIG. 7. Epidemic spreading in contact networks and inter-group contacts.(a)-(c) Connection between the total number of contacts between classes and the time steps until the peak of the epidemic in a SIR model.(d)-(f) Connection between the total number of contacts within each class and the fraction of infected individuals in each class in the stationary state.

FIG. 9 .
FIG. 9. Segregation of each income category with traditional segregation indices when calculated on (a)-(b) the adjacency graph and (c)-(d) the commuting graph

FIG. 10 .
FIG. 10.Economic segregation and criminality.(a)-(c) Correlation between the Spatial Gini and the criminality in US cities.(a) Violent crimes, (b) robbery crimes and (c) assault crimes.

FIG. 11 .
FIG. 11.Economic segregation and criminality.(a)-(c) Correlation between the Moran's I index and the criminality in US cities.(a) Violent crimes, (b) robbery crimes and (c) assault crimes.

TABLE III .
Corresponding income branch of each class or category in US dollars.