{'fraction_non_alphanumeric': 0.03248332534915429, 'fraction_numerical': 0.011902568448981096, 'mean_word_length': 4.546188432454527, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Random walk informed community detection reveals heterogeneities in large networks

Random walks on networks are widely used to model stochastic processes such as search strategies, transportation problems or disease propagation. A prominent biological example of search by random walkers on a network is the guiding of naive T cells by the lymphatic conduits network in the lymph node. Motivated by this case study, we propose a general framework to find network heterogeneities, which we define as connectivity patterns that affect the random walk. We propose to characterize and measure these heterogeneities by i) ranking nodes and ii) detecting communities in a way that is interpretable in terms of random walk, moreover, we propose iii) an approximation to accurately and efficiently compute these quantities on large networks. The ranking parameter we propose is the probability of presence field, and the community detection method adapts previously defined diffusion coordinates. In addition, we propose an interactive data visualization platform to follow the dynamics of the random walks and their characteristics on our datasets, and a ready-touse pipeline for other datasets upon download. We first showcase the properties of our method on toy models. We highlight this way the efficiency of our methods at contrasting two highly similar networks (identical degree distribution, same number of nodes). Moreover, we show numerically that the ranking and communities defined in this way are not redundant with any other classical methods (centralities, global mean first passage, louvain, node2vec). We then use our methods to characterize the lymph node conduits network. We show that the lymph node conduits network appears homogeneous and therefore has a global structure that promotes a uniform exploration of space by T-cells.


I. INTRODUCTION
Random walks on networks are a widely used model to describe search strategies [1,2], transportation problems [3], transmission in epidemiology [4,5] or diffusion of information [6]. In this model, random walkers hop from node to node while choosing with an uniform probability the edges on which to travel. The structure of the underlying network, such as its degree distribution and connectivity pattern, will thus determine how the random walk evolves over time. Can this connectivity pattern favour the exploration of some nodes over others? We will refer to any feature of the connectivity pattern which can bias the exploration pattern of random walkers as 'local heterogeneities'. In this paper, we develop analytical and computational tools to define local heterogeneities for any network by ranking nodes and computing communities based only on its connectivity and the features of the random walk.
The lymph node conduits network (LNCN) offers a prominent example of a biological network whose structure can affect its function [7]. Upon an infection, dendritic cells bring the virus' antigen to the lymph node. The encounter between these dendritic cells and the small subset of the naive T cells able to react to the antigen (one out of 1 000 000) triggers the proliferation of these relevant T cell. This process can be modeled as a search problem [8]. The dendritic cells stay still at a certain location on the network, and the naive T cells scan the lymph node for dendritic cells to test their specificity. The lymph node is spanned by a conduits (pipes conveying lymph) network and the naive T cells were shown to use the conduits as substrates for their migration [9]. Thus, without additional hypothesis, the modality of T cells search behaviour is a random walk guided by the LNCN. Does the network connectivity optimize the search by making some regions more accessible? This case study is the initial motivation to undertake the analysis of random walk dynamics on networks.
Well studied statistical properties of random walks on network such as the stationary state, cover time, relaxation time, mean first passage time, recurrence time or exit probability are heavily dependent on the degree distribution of the nodes [10]. In the case of the LNCN however, most of the nodes are of degree 3 [7] [11]. Thus, classical methods might not be adapted to detect local heterogeneities in a network such as the LNCN. In addition, the LNCN is a very large network with approximately 200 000 nodes, thus random walks based analysis methods [12] [13] are not tractable without an approximation.
The aim of our study is to propose a way to analyze the dynamics of random walks on potentially large networks using only the connectivity information. We consider only unbiaised discrete time random walk on an unweighted undirected network. We propose i) a random-walk based ranking of nodes, ii) a random walk-based interpretable community detection method, iii) a way to accurately approximate the results of i) and ii) for large networks, and iv) an interactive visualization platform to explore random walk on actual datasets. More precisely, we define and use a probability of presence field based on the random walk as a node ranking parameter (i). Moreover, we attribute to each node diffusion coordinates as defined by Pons and Latapy [12] and define clusters in the diffusion space that enable us to find communities with similar random walk associated properties (ii). We first show on three toy models that the ranking method is not reducible to any of existing classical centrality measures. Two of the toy models, described later, are very similar, and most of their nodes are degree 3. We show that the community detection method we propose is able to distinguish between these two toy model networks based on the characteristics of their respective communities. In addition, these communities are different than the ones made by clustering in other embeddings spaces previously proposed such as the louvain algorithm or the Node2vec algorithm [14]. Finally, we apply our framework to the lymph node (LNCN). Since the dataset that we have is a large network with around 200 000 nodes, we use the approximation (iii) defined similarly to the diffusion maps approach [15,16]. Applying the pipeline to the LNCN shows communities with no distinguishable feature. Thus, the LNCN network appears homogeneous, and does not seem to lead or trap the T cells into specific regions. All our analyses are available on the interactive visualization platform (iv) available at https://randomwalknet.centuri-engineering.univ-amu.fr/ and the pipeline can be executed on any other datasets upon download.

II. MATERIAL AND METHODS
To study the influence of the connectivity of a network on the exploration pattern of a random walker, we extract the connectivity information and ignore the spatial coordinates and edges lengths. Therefore, we address a discrete-time random walk on an unweighted, undirected and unbiaised graph. The dynamics of a random walk on a network can be derived analytically [12,15]. In this section, we recall this analytical framework, and introduce the measures we propose for ranking and communities detection. Then, we show how they can be approximated to be used on a large network. Finally, we introduce the toy models on which the pipeline will be tested, and for which all results can be visualized with an online graphical tool.

A. Analytical derivations
is the set of nodes and E = {(i, j)} is the set of edges. The connectivity information is encoded in the adjacency matrix A which is defined as A(i, j) = 1 if (i, j) ∈ E and A(i, j) = 0 otherwise. We consider a random walker that jumps from a node to an adjacent one at each time step.
1. Probability of presence departing from a given node At a given node i, the probability of choosing any of the adjacent node is 1/d i where d i is the degree of node i. Such a process is described with the transition matrix . The matrix T contains the transition probabilities of going from node i to node j in one time step: Figure  1A illustrates the approach.
The dynamics of the random walk on network is described by p(y, t|x), the probability to go from x to y in t time steps. This quantity is related to the transition matrix: for any time t = q ∈ N, we have p(y, q|x) = T q (x, y). The spectral decomposition of T offers another expression of T q .
We first note that T is similar [17] to a symmetric matrix T s defined as T s = D 1/2 T D −1/2 , therefore T and T s share the same eigenvalues. Since T s is symmetric, it is diagonalizable and has a set of n real eigenvalues {λ j } j∈[0,N −1] such that T s = V ΛV −1 . The left and right eigenvectors of T are defined as where Λ is the diagonal matrix whose components are {λ j } j∈[0,N −1] . The eigenvectors of the transition matrix T and the symmetrized transition matrix T s are related according to the relations Φ = D 1/2 V and Ψ = D −1/2 V . [15,16,18]. For any finite time t, p(y, t|x) can thus be decomposed into This decomposition is illustrated on Figure 1B.
2. Probability of presence integrating on all departure nodes, as a ranking parameter We can characterize the network by studying the probability of presence field at any time step, independently of the departure node. We introduce the probability of presence field p f , which is, for each node y, the sum of all the probabilities of reaching y from every departing node x in t time steps: The properties of this field as a ranking method will be developed in section III B.

Community detection
A community within a network is well defined as a set of nodes that share similar roles. Here we consider that sharing similar roles means that random walkers starting from these nodes follow similar trajectories. This is expressed by the diffusion distances, which can be computed for any time t. The diffusion distance between two nodes is small if random walkers starting at these two nodes are likely to be at the same location at time t.
We will use the diffusion distance defined in the following way and schematized in Figure 1C : where w(y) = 1/φ 0 (y) [16]. Following [15,16], we have that this distance can be computed as the euclidean distance in the diffusion space, i.e. in the representation space given by the eigenvectors of the transition matrix.
Methodological steps of our approach for the study of random walks on spatial networks at short, medium and large times scales A) The network adjacency matrix is turned into a transition matrix that models the dynamic on the graph. B) The transition matrix is then diagonalized, which enables to compute the transition probabilities in a parameterized way. This decomposition leads to an approximation, by taking only the first few eigenvalues and eigenvectors. C) Using only the left eigenvectors and the eigenvalues, we define the diffusion coordinates and the diffusion distance.

Approximation
When the network considered has a large number of nodes, for example n = 200, 000, the transition matrix has dimension n × n and requires a very large amount of storage (320 Go as a double array for n = 200, 000). Therefore, only sparse matrices of this size can be manipulated. A transition matrix is usually initially sparse but it has to be elevated at power t to encode the transition probabilities at time t (Figure 1.B), in the case of large networks, a few exponentiation are usually enough to loose the properties of sparsity. The spectral decomposition, provided the prior computation of the Ψ and Φ (Figure 1.B), enables to compute T t by only term-wise power elevation of the eigenvalues in Λ. However, the computation of Ψ and Φ requires the diagonalization of the symmetrized transition matrix T s , which cannot be completed fully for large n.
Therefore, we need approximated methods for large networks. In the literature addressing random walks on networks we could not find such an approximation. However, spectral decomposition of transition matrices are also used in the characterization of geometry of large datasets in high dimension space. The method called diffusion maps [16] is used to identify non linear coordinates of manifold, and include an approximation of the probability of presence p(y, t|x) that can be usefully adapted to our problem. The approximation consists in truncating the spectral decomposition to k dimensions which is given by:p and consequently, the transition matrix can be approximated asT with Ψ k = (ψ 0 , .., ψ k ), Φ k = (φ 0 , .., φ k ) and Λ k = diag (λ 0 , ..., λ k ). Because 1 = |λ 0 | > |λ 1 | ≥ ... ≥ |λ N −1 | ≥ 0, the larger t is, the more negligible the last terms (large j) are. To compare the actual T t and the approximatedT t k , we considered the spectral norm of T t −T t k 2,2 which can be either computed analytically as T t 2,2 = λ t k+1 or with the power iteration method [19], see supplementary material for more details. Numerical results are shown in section III D.
Similarly, the diffusion distance can be optimally approximated using the first k principal components [16]: B. Interactive visualization tool The method developed in this paper enables to monitor the probability density of a random walk from any starting point after any number of time steps p(y, t|x). As tracking the density of probability distributions over time is useful to explore visually the topology of a network, we developed an interactive visualization platform available at this address: https://randomwalknet.centuri-engineering.univ-amu.fr/. The visualization platform is divided into three parts. The first part contains 2 animations, one is the probability of presence field of a random walk starting from a given point, and the other is the probability field integrating on all departure nodes. The second part contains the analyses of the networks: betweenness centrality, closeness centrality, eigenvector centrality. The third part shows the results of the community detection algorithm.

C. Datasets
Our analysis tools to define ranking and communities are first showcased on a few toy models. The aim is to understand their features on smaller and well understood networks before applying it to the LNCN. In this section we describe the features of all the networks mentioned in our study, the toy models and the LNCN. Two networks were particularly helpful in understanding and testing the ranking and community detection methods described in the next sections. These two networks were obtained as outputs of a numerical minimal model for city morphogenesis [20] (details on the model can be found in the supplementary material), nodes being crossings and edges being streets portions. They both have the same number of nodes N=2050, number of edges E=3073 and degree distribution (degree 3 for all nodes except from the 4 nodes on the corners). However, the model parameter favours homogeneous density of nodes in one case (Homogeneous city network HCN) and in the existence of regions of high density of nodes in the other case (Polar city network PCN). Note that density here is understood according spatial distance (number of nodes per unit surface). There is no trivial intuitive relationship between density and connectivity as shown in Figure 2.B and C. On their upper panels the nodes are represented with their positions in physical space as generated by the model. For the HCN, the density appears uniform whereas for the PCN, there are two high density regions, near the left and right borders. On the bottom panels, the same networks are represented using only the connectivity information : which nodes are connected to which and which ones are not. The layout is obtained using a force-directed drawing algorithm [21], which minimizes the occurrence of edges crossing and edges as equally long as possible. We note that on the force-directed representation of the PCN the high density zones do not appear as distinguishable regions.

Real-World Transportation Network (RWTN)
This network serves as a an example of network that, contrary to the HCN and PCN has nodes that have higher degree than the others. This network represents the connectivity of the London transportation network (see Figure 2 A.) from [3], with 369 nodes that correspond to the train stations and 441 edges that correspond to the connection between them.

The lymph node conduit network (LNCN)
The network's connectivity, published in [7], was extracted from a segmented 3D microscopy acquisition of a whole mouse popliteal lymph node in which the conduits are made fluorescent by injection of labelled molecular tracer into the lymphatic vessels. The conduit network is restricted to the T zone in which T-cells are present. This is a 3D network made of 192,386 nodes and 274,906 edges. Most of its nodes are degree 3 (72%) , and 99% of nodes having degree between 1 and 4 (see degree distribution in supplementary material). There is no apparent distinctive high concentration of high degree centrality in any specific region. The characterization of this dataset is performed in the section III D.

A. Characteristics of the dynamics of random walks on the toy networks
When considering a random walk on a network, one of the first question is to understand its time scale. After a certain time, the probability of presence reaches the stationary probability distribution, which does not depend on the departure nodes. Indeed, it follows from equation 2 that the larger t is, the smaller the terms j > 0 will be, until they are negligible and the sum converges to p(y, t = ∞|x) = φ 0 (y). The speed at which this convergence is reached is a property of the considered network.

Time scales
Relaxation time It is the time at which the difference between the probability field and the stationary field is reduced significantly by a constant factor. This time is governed by the ratio between the magnitude of second largest eigenvalue of the transition matrix |λ 1 | and λ 0 = 1, τ = 1 1−λ1 [22] (See supplementary for the derivation). Convergence time Since the stationary probability field is reached when the probability field does not depend anymore on the departure nodes, we propose to measure the dependence of the probability of presence field with respect to departure nodes across time. This dependence should decrease over time and reach a plateau. We define the convergence time as the time step where the plateau starts. The dependence, or deviation D(t) is defined as the mean over all arrival nodes y of the standard deviation of p(y, x|t) over x (details in Supplementary). Table I summarizes the relaxation times and convergence times of the considered toy networks. Interestingly, the homogeneous city model (HCN) converges faster to its stationary distribution than the polar network (PCN) even though they have the same number of nodes and same uniform degree distributions. The relaxation time is linked to the "conductance" of the network which measures how easy it is to cross between the two most separate communities in the network [10]. This gives a first confirmation that there is a difference in higher order connectivity patterns between the two networks, even with an identical degree distribution and number of nodes and edges.

Probability of presence fields for the toy datasets
We consider now the probability of presence fields. To represent it in a informative manner, we chose t significantly smaller than the relaxation time τ as a relevant time step to show the probability of presence fields for each network.  interpret that it takes a very long time to penetrate to the center of this zone for a random walker. On the contrary, the second high density zone, on the left, has a high probability of presence.

B. Node ranking based on the probability of presence
The importance of nodes can be defined using a variety of criteria, which are measured by different indices. These measures assess the influence of each node either with respect to the shortest paths between nodes, such as the classical closeness and betweenness centralities, either with respect to random walk, such as the random walk betweenness centrality [23] [24], PageRank centrality [25], or 'Global mean first passage time' (GMFPT) [26].
Classical measures relying on random walk (PageRank and GMFPT), rely on the calculation of the stationary probability of presence. Thus, these measures are heavily dependent on the degree of the nodes. Here we propose to use p f (y, t), the probability of presence integrated on all departure nodes, as a ranking method. This ranking depends on time, allowing to access short time scales, before the effect of the degree of the nodes becomes dominant.
We compute the probability of presence on the toy networks HCN and PCN and on the RWTN. We show that the probability of presence field offers a ranking method that is not reducible to a classical centrality measure relying on shortest paths, and is also different from the GMFPT.

Shortest-paths based centrality measures for each network
We computed three different centrality measures based on the shortest paths, to later compare them with the probability of presence fields. Two of them are classical measures : closeness centrality and betweenness centrality [27]. We added another measure which we call 'Maximal remoteness', its definition is below.
We show that for a given network, the different centrality measures (detailed calculations in supplementary) give very different rankings from each other as shown in Figure 3. Furthermore, for a given centrality measure, we see major differences between the (HCN) and the (PCN). We will show in section III B 3 that p f (y, t) provides a ranking among nodes that is not redundant with any of the classical centrality measure.
The results obtained for each centrality measure on the considered networks are recapitulated below and shown in Figure 3: • The betweenness centrality (calculated with the NetworkX Python library [28] networkx.algorithms.centrality) is the number of shortest paths between two nodes that pass through that given node. The high betweenness centrality nodes are scattered within each of the networks. They don't seem to be correlated with particular spatial structures.
• The closeness centrality (calculated with the NetworkX Python library [28] networkx.algorithms.centrality) is the inverse of the sum of the length of the shortest path to all the nodes in the network. The high closeness centrality nodes correlates with the center of the HCN, and seem to correlate with one of the poles (right) in the PCN and anti-correlate with the other one (left).
• The maximal remoteness is a measure we propose to define as the length of the shortest path between u and v. In the HCN it has low values in the center of the homogeneous network and high values on the border. Interesingly, in the PCN, it highlights with high values both of the high density zones. It is the only centrality measure that seem to distinguish both high density zones with high values.
Interestingly, except for the maximal remoteness centrality, which is not a classical measure, the centrality measures fail at highlighting the two high density zones simultaneously (left one and right one) in the PCN. Indeed, the two high density zones appear with different extreme value ranges: one has very low centrality value and the other very high.

Random-walk based centrality measures
The eigenvector centrality [27], similar to PageRank ranking [25], was computed as the eigenvector of the adjacency matrix with the highest eigenvalue Ax = λx . It is interpreted as a centrality measure in which a node can have a high centrality value if it is connected to few other nodes but those nodes have a high centrality value. The results are shown on Figure 3.D. For the HCN, the eigenvector centrality results in a very symmetrical pattern. For the PCN, only one of the pole correlates well with high values of the centrality.
Another way to rank the nodes is to compute how much time in average it takes for a random walker to reach it, departing from any other node of the network. This is known as the 'Global mean first passage time' (GMFPT) [26]. It is also the inverse of what was previously called random walk centrality [13]. We show that the GMFPT gives a different ranking than the probability field. As shown in Figure 3.F, the high values are on the borders of the networks (Calculation in Supplementary). In the HCN, the four corners would be accessed last in average by a random walker. In the PCN, the four corners also have the highest values but there is also the left and bottom borders of the square. The polar region on the right has low values except in the very center.
Contrary to the other measures, high value of GMFPT means low influence.
Finally, we computed the Random walk betweenness. For a given node, this is the number of random walk paths between any two nodes that pass through the given node [23]. The corners appear with low values, and the rest appear with higher values that are scattered, so that no compact region is distinguishable.
3. The probability of presence field is not reducible to any of the previous centrality measures For each time point t the probability of presence field p f (y, t) is computed as in equation (3) and is compared to the normalized centrality values. The correlation coefficient (R) and normalized mutual information (NMI) are computed for each time point and shown in Figure  4 (See supplementary for the calculation of NMI). We exclude from the computation of R and NMI the four corners of the HCN and PCN which are degree 2, to focus on the ranking of nodes which are all the same degree (degree 3). Figure 4 shows all relations (R and NMI) between each centrality measure (6 measures, either shortest paths or random walk based) and p f (y, t), for each of the three toy networks, with respect to time t. Even though some of the centrality measurements have high correlation with p f (y, t) in certain networks, none of the centralities maintains this high correlation in all networks. For instance, the centrality that relates most closely to both p f (y, t) for HCN and PCN is the eigenvector centrality. But eigenvector centrality is poorly correlated to p f (y, t) for RWTN. On the contrary, betweenness centrality is well correlated to p f (y, t) on RWTN, but not on HCN and PCN. Thus, we demonstrate that p f (y, t) is not reducible to any of the centrality measurements.
C. Detecting, characterizing and interpreting communities in the network using diffusion distances Are there communities of nodes that show distinct random walk guiding features due to their connectivity? We presented in section II A 3 a method to define communities based on their random walk based properties [12] [13]. We use the coordinate of each node in the diffusion space to define clusters. It is particularly interesting for the two toy models HCN and PCN, because they have exactly the same number of nodes, edges, and the same degree distribution (degree 3 at almost all nodes). The proposed community detection method will appear valuable if it is able to distinguish between the two instances, one having high density zones and the other not. More precisely, we expect to find identical communities in the HCN and more diverse communities in the PCN. Then, we compare this clustering to the classical Louvain algorithm that optimizes modularity [29] and the Node2vec   algoritm [14] and show that they yields different communities.
1. How does diffusion distance relate to topological distance?
We base our community detection method on the diffusion distance. How does this distance differ from the topological distance (length of the shortest path between the pair of nodes)? In both HCN and PCN, we find a good correlation between both distances before convergence, as shown in the right panel of Figure 5.A. However, this correlation is better in the HCN than in PCN, suggesting that the diffusion distances carry information that can distinguish between these two models. This makes it a promising tool to address the question : How to find regions with distinct connectivity within the network when the network has the same degree at all nodes?

Characterization of random walk-based clusters in the city streets networks (HCN&PCN)
Clustering in the diffusion space and quantification As described in section II A 3, we define communities in the HCN and PCN by k-means clustering in the diffusion space. The time t chosen to compute the diffusion coordinates is t = 500, the same as for the computation of p f (y, t). The number of clusters was chosen to be 4 for both models, HCN and PCN, with Dunn indexes being respectively 2.8 and 2.2 (highest Dunn index in the HCN and second highest Dunn index in the PCN) and satisfactory cluster sizes.
These clusters are shown in Figure 5B-C for the HCN and Figure 5G-H for the PCN. We show both the clusters in the initial space with the coordinates generated by the Pousse et al. model [20] (B and G), and the clusters in the umap representation of the N-d diffusion space (C and H).
In the PCN, two clusters (red and yellow in Figure  5) that were defined in the diffusion space correspond visually to the high density zones, whereas in the HCN, the four clusters simply divide the space equally in quadrants.
Clusters characteristics • Intra-cluster diffusion distances variability: For each cluster, for all the pairs of nodes within the cluster (x 1 , x 2 ) we compute the diffusion distance ||Ψ t,k (x 1 ) − Ψ t,k (x 2 )|| 2 . In the PCN, the histograms of the pairwise diffusion distances ( Figure  5.K) show that the maximal intra-cluster diffusion distances are smaller in the poles clusters than the non poles clusters, whereas they are identical in the HCN (Figure 5.F).
• Inter-cluster diffusion characteristics variability: For each cluster, the mean entry and exit probabilities were computed respectively as p in (C) = Overall, we show that the characteristics of the clusters obtained in the diffusion space are successful at differentiating the HCN from the PCN and at distinguishing the two high density zones, the poles, in the PCN as two clusters with different properties than the non-poles clusters.

Diffusion coordinates-based communities are different from Louvain and node2vec communities
We compared the communities defined in the previous section with two other community detection methods, the Louvain algorithm [29] representative of the modularitybased family and the node2vec algorithm [14] which, like our method, relies on diffusion but does not use the diffusion space.
Louvain algorithm This algorithm for community detection based on modularity is widely used. We applied the Louvain algorithm to our toy models network, HCN and PCN. The best modularity score is reached for a partition of nodes into 7 clusters, both for HCN and HCN and PCN. In order to compare the Louvain algorithm and the clustering in the diffusion method, we performed 7-means clustering in the diffusion space for both networks at time=500. After performing an optimal matching of clusters between the ones made in the diffusion space and the ones made by Louvain algotithm that minimizes the Jaccard distances, we show that both clustering give different results, as shown in Figure 6 with mean Jaccard distances being J dist = 0.61 and J dist = 0.49 respectively for the HCN and PCN. Node2vec algorithm The Node2vec algorithm [14] provides an embedding of the nodes of the network that preserves the neighbourhoods of the nodes, the neighborhoods being sampled by simulated biased random walks. The bias is defined by setting p and q, respectively the return parameter and the in-out parameter. At each step, the probability of going back to the previous node is 1 p and the probability of going to any other neighbour-ing node is 1 q provided that they are not also immediate neighbours of the previous node. Random walks are simulated and samples of k successively occupied nodes are recorded. Theses sampling define the neighbourhoods of nodes. A stochastic descent is run to embed the nodes in a d dimensions space while preserving as best as possible the neighbourhood relationships. We computed the embeddings using the PecanPy Python implementation [30]. We explored the parameters to optimize the embeddings so to get clusters as similar as possible to ours with the procedure described in Supplementary. The 10 combinations with lowest scores all present rather compact clusters in space, but they all show a band-like structure. The clusters obtained with the set of parameters {d = 10, q = 0.1, p = 5, k = 50}, optimized for the PCN, are shown on Figure 6.E-F. Thus, the embeddings by Node2vec lead to clusterings than the ones obtained with our method.

D. Application to a large biological network: the Lymph Node Conduit Network LNCN
The intial motivation of our study was to assess the potential influence of the connectivity of the LNCN on the search process of the naive T cells for the antigen carrying dendritic cells. The encounter between these dendritic cells and the subset of relevant T cells triggers the proliferation of those T cells. Thus the modality of this search is crucial to determine how fast the adaptive response can be triggered [8]. The conduits are narrow lymphatic tubes that convey molecules but no cells. They are ensheated and constructed by fibroblastic reticular cells (FRCs), which control T-cell life in many ways. The naive T-cells migrate in close association with the FRCs, therefore using the conduits as substrate for their migration, like on 3D routes network [9].
The relation between the network and its functions have been studied before. The FRC network was shown to have small-world properties and shows high robustness to perturbation on a portion of a mouse lymph node [31]. Besides, a reconstruction of complete conduits network showed a higher density of conduits in the superficial regions than in the deep regions [7]. The authors ran an agent-based model on parts of the reconstructed network to understand the relation between spatial density of nodes and the motility coefficient of naive T cells, which measures how fast the agent moves away from its departure point (in analogy with diffusion coefficient). They showed that there is no significant difference in dense and sparse regions in terms of motility coefficient. The question we address is complementary to their study. They have determined the role of spatial density of nodes, we address the question of how does the connectivity of the network influence the exploration behaviour of naive T-cells?

Probability of presence fields in the LNCN
With the online tool introduced in section II B, it is possible to visualize the probability fields departing from a given node at different time steps, as shown in Figure  7.A. The relaxation time is τ = 17797. With an estimation that the naive T-cells move with an average speed of 13µm.min −1 [32] and that the average length of an edge is 10µm [31], we can make the rough estimation that one time step represents 0.77 min. Thus, the relaxation time for the LNCN, is about 208 hours (8.5 days).
We approximated the probability of presence field (equation (8)) and diffusion coordinates (equation (10)) by using a partial spectral decomposition with only k = 2000 largest eigenvalues and corresponding eigenvectors. Figure 7.B shows for k = 2000, the normalized spectral norm between the exact transition matrix T t and the approximated T t k for different t (see Supplementary 3.2 for the complete calculation). This distance decreases exponentially. For k=2000, the distance is smaller than 1% from time steps t ≥ 250. As a side note, it is possible to accurate approximation (less than 1% error) for time steps as small as t = 100 for k=5000.
In the following we will define communities by using diffusion coordinates at t = 5000 (about 64 hours), which is significantly smaller than the relaxation time, and for which the normalized spectral norm is the order of 10 −9 .
The probability of presence field at t = 5000, shown in Supplementary, did not show any difference between regions of the network.

The lymph node network communities are all similar
Intra-clusters diffusion distances are similar for all clusters We computed 500 clusters with k-means clustering in the diffusion space at time t = 5000. The cluster size varies among the clusters from 7 to 1207 nodes and has an average value of 385 nodes. The 7 smallest clusters which represent 0.05% of the nodes, are also the ones with highest mean intra-cluster pairwise diffusion distances. The intra-cluster pairwise diffusion distances weakly anti-correlate with the cluster size with R=-0.52 and NMI=0.17. Apart from the 7 smallest cluster, there are little differences among the remaining clusters in the pairwise diffusion distances within each cluster, most the values being spread continuously between 0.05 and 0.2. It does not appear that there are distinct populations of clusters with well separated values.
Inter-clusters diffusion characteristics are similar for all clusters For each cluster, like previously for HCN and PCN, we computed the mean probabilities to enter and to exit from the cluster. The mean entry and exit probabilities are roughly equal to each other for each cluster. Setting aside the clusters with extreme small sizes which show high standard deviations (first 30 and last 30 clusters which represent together 1.2% of the total number of nodes), the entry and exit probabilities vary little: almost all clusters are centered around 0.5 × 10 −5 with standard deviation about 10 −5 .
We defined a mixing index which describe how the clusters obtained in the diffusion space are separated in the physical space. This index is calculated for each cluster as the number of edges between a node of this cluster to a node of another cluster (exogeneous) divided by the number of edges connecting nodes inside the cluster of interest(endogeneous). The mixing index is expected to have a scaling relationship to the size of the cluster. Indeed, the number of endogeneous edges is expected to vary linearly with the cluster size n c (analogous to the volume of the cluster) whereas the number of exogeneous edges is expected to vary with n Neighbourhood relationhips are the same in the physical space and in the diffusion space The mixing index is rather small (around 10%) for most clusters except for the smallest ones (outliers for clusters size). This low value of mixing index shows that the nodes which are neighbours in space are also close in the diffusion space. Since the diffusion distance describes similarity in trajectories of random walks, we can interpret it as the fact that random walkers departing from neighboring nodes will follow similar trajectories.
Overall our results on the LNCN indicate that the structure of the lymph node is highly homogeneous with no region-specific behaviour of a random walker. Therefore, to the initial question: does the structure of the conduit network influence its guidance function and affect the ability of T-cells to interact with antigens? We can answer that the relative uniformity of the clusters characteristics presented on Figure 7 establishes that the structure of the network is promoting a non biased search of the T-cells over the entire lymph node.

IV. CONCLUSIONS AND DISCUSSION
We proposed in this article a framework to find heterogeneities in networks connectivity that can bias a random walker exploration pattern. In particular, we show the power of this framework to define interpretable heterogeneities even in networks with uniform degree distributions. i)We derived the probability of presence field and used it as a ranking parameter for nodes. We showed that this ranking is not redundant with any other classical centrality measure. ii) We showed that the clusters obtained in the diffusion space provide interpretable communities of nodes with comparable attributes between communities. (iii) We applied this pipeline to the LNCN which guides naive T cells migration. We showed that it has a connectivity that promotes a uniform behavior in all the areas of the organ. This finding suggests that the structure of the network is designed in a way that wouldn't FIG. 7. A) probability of presence departing from a given node B) Estimation of the spectral norm relative error of the approximation k = 2000 C)-G) Clusters analysis in the LNCN dataset (t = 5000). C) Mean and standard deviation for the pairwise diffusion distances for each cluster, ordered in increasing mean value. D) cluster sizes ordered like in C E) Entry and exit probability for each cluster F) cluster sizes, ordered like in E G)Mixing rate for each cluster with respect to the cluster size and the fitted estimated values from a scaling model calculated as mixing = 0.7195 × 1 n 1/3 c bias random exploration of naive T-cells in the whole space. This does not exclude that the structure is optimized for the T cells search for antigen-carrying dendritic cells. Furthermore, the LNCN has other functions than guiding the T cells migration, which can also constrain its connectivity. They convey lymphatic fluid containing crucial molecules of the immune response, and they must be remodelled in a short time in case of inflammation [33].
In this work, we have restricted our attention to random walk with no memory on unweighted and undirected network that could be embedded in 2D or 3D. Going further, it would be interesting to include more complex behaviors, such as a dependency on the path already ex-plored, which is translated as non-markovian behaviors [34] and can be well adapted to describe animal behavior [35]. Another specificity of our study is that we ignore the edge lengths. Including them in the framework would require to address continuous time random walk to account for the fact that some edges take more time than others to be crossed, which can also be formulated as: the waiting time between two jumps will be long if the second jump involves crossing a long edge. Another possibility would be, instead of a continuous time random walk, extend the discrete-time random walk framework by adding intermediary nodes that divide long edges into several fixed lengths edges, which all take one time step to cross. However, the process won't be Markovian anymore since no return should be allowed once engaged into a series of adjacent intermediary nodes that belong to the same initial edge.
Finally, we have only considered spatial networks (which can be embedded in 2D and 3D), but this framework is readily amenable for use and prediction in higher dimensional space. It could benefit various field where network based inference of dynamical systems is important. One interesting application concerns the study of trajectories in the space of gene expression, for example when studying differentiation trajectories in developmental biology. In this setting, individual cells can be embedded in a high dimensional space, where the coordinate represent the expression of individual genes. Cells obtained from a developing organism often have a lowdimensional structure, a manifold, which can be represented as a network. The dynamics on this structure are then well described as random walks on this network [36].
T is non symmetric so its spectral decomposition is not guaranteed. Let T s = D 1/2 T D −1/2 the symmetrized transition matrix. T s can be composed as T s = V ΛV T , where Λ is a diagonal matrix with the eigenvalues of T s , which as the same as the eigenvalues of T and V is an orthogonal matrix (V V T = Id) which columns are the eigenvectors of T s . Then, gives t, we calculate it as the standard deviation of p(y, x|t) over x (dependence on x), averaged over y.
D decreases as time increases down to a plateau starting at an elbow on the log scale representation. We identify this elbow as the convergence time when the decrease is no longer detectable numerically.

London transportation network -RWTN
We used the data assembled by CoMuNe Lab [3]. The nodes represent train stations and edges route between stations considering Underground, Overground and DLR (Docklands Light Railway). The data was initially a weighted network. We simplified it into an unweighted network (two stations are connected by an edge if they are only one station apart in any of the train lines). The network is comprised of 369 nodes and 441 edges. The degree distribution is mainly centered on 2 as shown in Figure 8 because of long underground lines, with a few number of stations having degree 5,6,7 corresponding to transportation hubs.

City street network generative model -HCN & PCN
HCN and PCN are two instances of generative model for city streets morphogenesis [20]. The networks are constructed by iterative subdivisions of an initial rectangle.
The streets segments are the rectangle boundaries. At each generation, the rectangle with the highest score U rec is divided. U rec reads: where d(C rec , i) is the distance from the center of mass of the cell of interest to the node i and L rec is the length of the cell.
The first term favours further subdivision of the high density zones while the second term favours the divisions in the low density zones. Their respective weights are tuned by the choice of α and β. When α > β the morphogenetic process is leading to cities with dense poles. On the other hand when α < β the morphogenetic process is leading to homogeneous streets density.
The HCN (Homogeneous City Network) and PCN (Polar City Network) were generated respectively with α < β and α > β.
Both networks contain 2040 nodes and 3073 edges. The two networks have exactly the same degree distribution with all nodes having degree 3 except from the 4 corners which have degree 2 as shown in Figure8.

Lymph node conduit network -LNCN
We used a complete reconstruction of an entire mouse popliteal lymph node conduits network (850 x 750 x 900 µm) obtained from microscopy data [7]. The reconstructed network contains 192,386 nodes and 274,906 edges, its degree distribution is shown in Figure8. The probability of presence field, p f (y, t), was estimated on a sampling departure nodes as 10% of the total number of nodes. Figure 9A-F show the evolution in time of the distributions of p f (y, t), converging to the degrees distribution. The p f (y, t) field in space is shown for times t=100 and t=5000 on Figure9G-H. For t=5000, the time of reference chosen, there is no distinct regions. The variation appear very local. The shortest-path based centralities could not be calculated on the LNCN because of the computational complexity (the complexity of finding the shortest path is at best O(V.E)).
C. Numerical approximation of the spectral decomposition for the LNCN

Numerical approximation of p(y, t|x)
The computation of the full spectral decomposition of a symmetric matrix A of size N × N requires O(N 3 ) floating point operations and O(N 2 ) coefficients to store, which becomes intractable even for medium-sized N . In particular, this algorithm fails for the LNCN dataset where N = 200, 000.
However approximating T t by T t k = V k Λ t k V T k only requires the computation of the truncated decomposition V k Λ k V T k where V k is of size N × k and contains the eigenvectors associated to the k largest eigenvalues and Λ k = diag(λ 0 , . . . , λ k ) is a diagonal matrix of size k × k. This observation can be leverage to implement spectral decomposition algorithm that are fast provided that k ≪ N . We used the eigsh function from the scipy library which implements a Lanczos method [19,37]. This method is iterative and mainly needs the computation of k matrix-vector products and an orthogonalization step of complexity O(k 2 N ). For general matrices, the cost of k matrix-vector products is O(kN 2 ), but for sparse matrices, this complexity reduces to O(kLN ) where L is the average number of non-zero coefficients in each row. The total complexity of the Lanczos method for sparse matrices is therefore bounded by O(k 2 N + kLN ). For the LNCN dataset, the matrix T s contains roughly L = 3 non-zeros per row, leading to a very efficient diagonalization procedure.

Estimation of error
In this section we describe the computation of the approximation error between T t to T t k as defined in equation (9) in the main part of the paper.
Let A be a matrix of size M × N , its spectral norm A 2,2 is given by Importantly, the spectral norm of a matrix is equal to its largest singular value i.e. A 2,2 = σ 1 (A). This norm induces the spectral distance A − B 2,2 between two matrices of size M × N and measures the maximal error that can be made between the two vectors Ax and Bx with x ∈ R N .
Large graphs prevents the matrix ∆ T ∆ to be stored in memory, we therefore need an indirect way of computing its largest eigenvalue. The power iteration method [19, Section 8.2.1 p406] computes the largest eigenvalue of a given matrix using only its matrix-vector products, which makes it suitable for large scale graphs.
We implemented the power method on ∆ T ∆ through the two the matrix-vector products with ∆ T and ∆. This products are made of: • t matrix-vector products with the sparse matrix T and T T . For a matrix of LN non-zero coefficients, this step costs O(tLN ) operations.
• matrix-vector products with T t k = ΨΛ t k Φ T and ( T t k ) T = ΦΛ t k Ψ T which both cost O(kN ) operations.
For the LNCN dataset, the matrix T contains roughly L = 3 non-zeros per row, leading to a fast computation of ∆ 2,2 .
We define the relative error as ∆ 2,2 T t 2,2 , with T t 2,2 also computed with power iteration method.

Betweenness centrality
The betweenness centrality value of a node u is the number of shortest paths between two nodes that passes through that given node.
where σ(s, t) denotes the set of paths connecting s and t and σ(s, t|v) is the set of paths connecting s and t and going through node u. The betweenness centrality was calculated using the NetworkX Python library [28] networkx.algorithms.centrality.betweenness centrality implemented from [38] 2. Closeness centrality The closeness centrality of a node is computed as the inverse of the sum of the length of the shortest path to all the nodes in the network.
where d(u, v) is the length of the shortest path between nodes u and v, and n the number of nodes that can reach u. The closeness centrality was calculated using the NetworkX Python library [28] networkx.algorithms.centrality.closeness centrality implemented from [39].

Maximal remoteness centrality
We define the maximal remoteness centrality as C maxrm (u) = max v d(u, v) with d(u, v) the length of the shortest path between u and v. This is not a classical measure but it captured well the two high density zones in the polar model. E. Classical random-walk based centralities

Eigenvector centrality
The eigenvector centrality [40] is the score obtained by computing the eigenvector of the adjacency matrix with the highest eigenvalue Ax = λx using eigs function in Matlab. It is interpreted as a centrality measure in which a node can have a high centrality value if it is connected to few other nodes but those nodes have a high centrality value.

Random walk betweennness
The random walk betweenness [23], also known as current-flow betweenness centrality [24]. For a given node, it measures the number of random walk paths between two nodes that pass through that node. It was computed using NetworkX library, using "Current-Flow Betweenness". where H(X) and H(Y ) are the entropies of X and Y and H(X, Y ) is the joint entropy. The entropy is H(X) = x∈X p(x) log p(x).
The probabilities p(x) were discretized by dividing the values x of X into histogram bins using 'histcounts' automated binning algorithm from Matlab so that H(X) = N bins i ni ntot log ni ntot , n i being the number of hits in bin i and n tot the total number of values in X. The joint probability p(x, y) was discretized using 2-d histogram with number of bins equal to the minimum between the automated bins number for X and for Y N bins = min(N bins(X) , N bins(Y ) ) H(X, Y ) = N bins i N bins j nij ntot log nij ntot . When N bins = 1, the NMI value is chosen to be 0.

G. Interactive visualization
The code is developed in Python. There are 2 files : the first script "graph object.py" define a class object for a network. All the fonction needed to either load a graph and compute analysis are defined in. The main script aims to both compute the network analysis and create figures. At the end of the script, the interactive application is deployed to visualize all the analysis. We used the module plotly for the figures and dash for the app.

H. Node2vec parameters optimization
The optimization was done with a combination of two criteria that evaluate the 4-clusters done on the ddimensions embedding. The first one is the mixing index < m({C}) >=< exogeneous endogeneous > C and the second one is the mean Jaccard distances between their clusters and ours, after optimizing the pairing of clusters. We tested all the combinations of parameters in the sampling p ∈ {0.05, 0.1, 0.5, 1, 5, 10, 20}, q ∈ {0.05,0.1,0.5,1,5,10,20}, k ∈ {3, 5, 10, 20, 50} and d ∈ {2, 3, 5, 10, 50, 100, 200}. We compared the two populations p < min(q, 1) and p > max(q, 1), with the other parameters varying. The scores s({C}) = <m({C})>+<J({C})> 2 ( Figure 10.A) were smaller for the population p > max(q, 1), which corresponds to the no-return configurations. We kept only the combination of parameters with p > max(q, 1). Among these combinations, we computed the correlations between the combined score s({C}) and each of the other parameters p,q,k and d. There was a weak anticorrelation with k (R=-0.55) so we set k = max({k} = 50). Then, we ranked the remaining combinations by ascending scores and we plotted the corresponding scores. The curve shows a steep slope then a plateau ( Figure  10.B). With return:p < min(q, 1). The score S({C}) to be minimized combines the mixing index and the Jaccard distance to the clusters make on the diffusion space at t=500 shown previously. C)4-means clustering done on Node2vec embedding with parameters {d = 10, q = 0.1, p = 5, k = 50}, which gives the lowest score.