WITHDRAWN: Biased random walk in time-varying network

Accepted Manuscript Biased random walk in time-varying network Yan Wang, Xinjian Xu, Zongcai Jiang

PII: DOI: Reference:

S0375-9601(15)00075-4 10.1016/j.physleta.2015.01.014 PLA 23051

To appear in:

Physics Letters A

Received date: 29 September 2014 Revised date: 17 January 2015 Accepted date: 17 January 2015

Please cite this article in press as: Y. Wang et al., Biased random walk in time-varying network, Phys. Lett. A (2015), http://dx.doi.org/10.1016/j.physleta.2015.01.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights

• We consider the concurrent dynamics of the biased random walk and the time-varying networks at the same time scale. • Adjust the weight parameter θ to get the effective navigation. • In out-linking network, θ = 1.0 is the more effective navigation. • In in-linking network, θ = −1.0 is the more effective navigation.

Biased random walk in time-varying network Yan Wang1 , Xinjian Xu2∗, Zongcai Jiang3 1

Department of Computer Science, Guangdong Polytechnic Normal University, Guangzhou, 510665 2

3

Department of Mathematics, Shanghai University, Shanghai, 200444

Department of Mathematics and Information Science, Henan University of Economics and Law, Zhengzhou, 450046

Abstract We investigate the concurrent dynamics of biased random walks and the time-varying network (this network model is proposed in Sci. Rep. srep00469(2012)), where the preferential transition probability is in terms of the edge-weighting parameter. We also obtain the analytical expressions for stationary distribution, global mean first passage time (GMFPT) and coverage in directed and undirected networks, all of which depend on the weight parameter. Appropriately adjusting this parameter, more effective search strategy can be obtained compared with the unbiased random walk, whether in directed or undirected networks. Since network weights play a significant role in the diffusion process. PACS numbers: 89.75.Hc, 05.40.Fb

1

Introduction

Random walks are the natural framework to study diffusion, transport, navigation, and search processes on complex networks which have attracted increasing interest in the past few years[1, 2, 3, 4], because it has relevant applications in such contexts as spreading dynamics (i.e., virus or opinion spreading) and searching, animal and human mobility and so on. However, most previous work hold in two different limits. In the first case, the network is considered static. The time scale describing the evolution of its structure (τN ) is much slower than the time scale of the dynamical processes (τDP ), i.e., τN  τDP . In the second case, instead the networks are modeled as annealed[5, 6, 7]. It evolves on a much faster timescale, i.e., τDP  τN , and the dynamical processes effectively sense only time-averaged properties of the underlying network. While many real networks are dynamic structures in which connections appear, disappear, or are rewired on various time scales[8, 9, 10], thus they can have profound consequences for the dynamical processes taking place upon them at the same time scale. Nevertheless, these two limits are not appropriate for this phenomena, i.e., τN ∼ τDP . For this reason, it becomes very interesting to consider the concurrent dynamics of the random walk and the time-varying networks[11] thus evolve at the same time scale. Here the considered network is based on the activity potential ai of each node i, which is extracted from a distribution F (a) and describes its probability to establish links per unit time. Then we define the concurrent dynamics as follows: (a) At each discrete time step t, the instantaneous network G∇t starts with N disconnected nodes; (b) With probability ai ∇t, each vertex i becomes active and generates m links that are connected to m other randomly selected nodes. Nonactive nodes can still receive connections from other active vertices, thus the instantaneous network G∇t is generated. At the same time step t, the walker diffuses for a time ∇t; (c) At the next time step t + ∇t, all the edges in the network G∇t are deleted and the process starts over again to generate the network Gt+∇t , and the walker also diffuses for a time ∇t. All interactions have a constant duration τG = ∇t defining the time scale which describe the evolution of the network. Ref.[12] has dealt with the unbiased random walk upon this activity-driven model in the limit ∇t → 0. They derived the exact formulas of the stationary state of the random walk and the mean first passage time (MFPT) in both directed and undirected time-varying networks, it has also ∗ corresponding

author:E-mail address: [email protected], [email protected]

shown strikingly different behavior occurred in contrast to that happening in quenched and annealed topologies[13]. Ref.[14] considered the behavior of random walks on the time-integrating activity-driven network and provided the numerical and analytical results of the stationary probability. In this work, without loss of generality, we set ∇t = 1 and present a biased random walks, where the weight ωij of the edge linking nodes i and j has strong correlations with their degrees ki and kj , and is assumed to be ωij = (ki kj )θ , with θ being a controllable parameter. Such a form of edge weights can be observed in various real weighted networks[15, 16]. Ref.[11] finds that degree distribution of the integrated time varying network takes the form which is similar to activity distribution, i.e., P (k) ∼ P (a), thus a k ∼ a linear correspondence should be apparent. Consequently, the weight connecting nodes i and j takes value ωij = (ki kj )θ ∝ (ai aj )θ when vertices i and j are connected, and zero otherwise. In order to quantify the searching efficiency of the network, we consider the coverage C(t)[17, 18, 19] which is defined as the number of different vertices that have been visited by the walker at time t, for the walker initially taken over the stationary distribution of the set of starting points. Furthermore, based on matrix spectral decomposition[20] and generation function theories[21], we will derive analytically the formulas of the stationary distribution, coverage, MFPT Tij from node j to node i, and the global mean first-passage time (GMFPT) to a target node i, i.e., < Ti >. Here the space average < Ti > is taken over the stationary distribution of the set of starting points, which is some distinct from the one observed on their time-aggregated counterparts. This paper is structured as follows. In Sec. II we investigate the biased random walk dynamics on complex network, recall some basic definitions, present an general expressions of stationary distribution, coverage C(t), MFPT, GMFPT. In Sec. III the previous approach is considered on undirect time varying network. In Sec. IV, all qualities can be extended to the case of direct networks. Using detailed balance condition, we change the asymmetric transition matrix into symmetric form and the exact formulas of MFPT, GMFPT can be obtained. Finally, we summarize our results and comment on some perspectives In Sec. V.

2

Overview of the random walk process

Here we will present some results of discrete-time random walk on general networks which is connected and non-bipartite, then there is a unique stationary probability. Let pij (t) be the probability that the walker is on the node i at time t from starting node j, and Πi→k be the transition probability rate that the walker moves from node i to node k per unit time. Then pij (t) obeys the master equation pij (t + 1) − pij (t) =

N  k=1,k=i

pkj (t)Πk→i −

N  k=1,k=i

pij (t)Πi→k

(2.1)

The first part of the right-hand side of Eq. (2.1) represents the transition probability from all other nodes to node i, and the second part represents the transition probability from i to other nodes. If we introduce the transition matrix ⎧ i = k; ⎨ Πk→i , N  (2.2) Mik = Πi→s , i = k. ⎩ 1− s=1,=i

→ − and probability vector P (t) = vol(p1j (t), p2j (t), · · · , pN j (t)), then Eq. (2.1) can be in a more concise form → − → − P (t + 1) = M P (t). (2.3) t Let M t express the tth power of matrix M . Its ijth entry denoted by Mij represents the probability of which a walker moves from node j to node i in t steps. Introduce Dirac notation in the following and using the method of the matrix spectral decomposition when eigenvectors of M span the whole space[20], we can solve

pij (t) =< i|M t |j >=

N  s=1

λts < i|ψs >< χs |j >,

(2.4)

where λi is the ith eigenvalue of the transition matrix M , < χi | and |ψi > are the normalized orthogonal left eigenvector and right eigenvector belonging to λi . These eigenvectors satisfy < χi |ψj >= δij , N  |ψi >< χi | = I,

(2.5)

i=1

where I denotes the identity matrix of order N .

2.1

Expressions for stationary distribution and coverage

Since M is a left stochastic matrix with each column summing to 1, then λ1 = 1 and |λs | < 1(s ≥ 2) hold under non-bipartite graph, let left eigenvector belonging to λ1 be < χ1 |i >= 1, ∀i = 1, 2, · · · N . Taking the limit t → ∞ on both hands of Eq. (2.4), we get lim pij (t) =< i|ψ1 >, ∀j ∈ Z + (j = i),

t→∞

(2.6)

thus |ψ1 > is the stationary distribution. On the other hand, the stationary distribution satisfies balance equations < i|ψ1 >= < i|ψ1 >= N  j=1,j=i

N 

Mij < j|ψ1 >,

j=1 N 

j=1,j=i

Πj→i < j|ψ1 > +(1 −

Πj→i < j|ψ1 >=

N  j=1,j=i

N  k=1,k=i

Πi→k ) < i|ψ1 >,

(2.7)

Πi→j < i|ψ1 >,

that is N 

Πj→i < j|ψ1 >=

j=1

N 

Πi→j < i|ψ1 > .

(2.8)

j=1

The left-hand side represents the total flow from all nodes j into node i, while the right-hand side represents the total flow out of node i into all other states j. Thus (2.8) equates the flow into and out of nodes. The coverage C(t) is defined as the number of different vertices that have been visited by the walker at time t, it can be interpreted as the searching efficiency of the network, measuring the number of different individuals that can be reached in a given number of time steps, averaged for the walker initially taken over the stationary distribution, that is N C(t) 1  = [1 − (1− < i|ψ1 >)t ], N N i=1

(2.9)

where < i|ψ1 > is ith component of stationary distribution that the walker occupies the node i, or in other words, it means the walker will land on node i after a jump, because < i|ψ1 >=

N 

Mij < j|ψ1 > .

(2.10)

j=1

2.2

MFPT between two nodes and GMFPT

How fast is the random walk motion? To answer to this question, we study the MFPT and GMFPT. The first-passage probability qij (t) from j to i after t steps satisfies the relation pij (t) =

t  s=0

qij (t)pii (t − s).

(2.11)

Introducing the generating functions ∞ 

 ij (x) = Q

qij (t)xt ,

(2.12)

pij (t)xt ,

(2.13)

t=0

and Pij (x) =

∞  t=0

we obtain  ij (x) = Pij (x) , Q ii (x) P

(2.14)

so the MFPT from node j to node i is given by ∞

Tij =

  ij (x) dQ |x=1 = tqij (t). dx t=0

(2.15)

Substituting (2.4) into (2.13) and (2.14), after some detailed calculation, one obtains N

Tij =

 1 1 [< i|ψs >< χs |i > − < i|ψs >< χs |j >] , < i|ψ1 > s=2 1 − λs

(2.16)

and GMFPT to a target node i, i.e., < Ti >, for the walker initially taken over the stationary distribution of the set of starting points is < Ti >= =

3

1 1−

N  j=1,j=i

1 (1−)

< j|ψ1 > Tij N  s=2

(2.17) < i|ψs >< χs |i >

1 1−λs .

Results of undirect time varying networks

In the following, we define discrete random walk process on time varying networks as follows: at each time step t the instantaneous network Gt is generated, and the walker diffuses for per unit time. After diffusion, at time t + 1, a new network Gt+1 is generated. The concurrent dynamics of the random walker and the network thus take place with the same time scale, namely, the walker can get trapped in temporarily isolated nodes i until a new link involving i occurs at the successive time step, but it never occur in static networks. Ref. [12] has analyzed unbiased random walks occurring on time varying network, here we will investigate the case of biased random walk, introducing the special weight-assigning form ωij = (ai aj )θ if nodes i and j are connected. Note that at any time t, ij−edge is generated by two cases: (I) node i becomes active and chooses to connect to node j (with probability ai m N ); (II) node ). In these two cases, the instantaneous j becomes active and connects to node i (with probability aj m N average degree of node i are kiI = m + kiII = 1 +

N  j=1,=i N  j=1,=i

ai aj m N = m + m < a > −m N ≈ m + m < a >,

(3.18) ai aj m N = 1 + m < a > −m N ≈ 1 + m < a >,

 ωis = respectively. Moreover, the instantaneous average strength si of node i is defined by si = s∈Ω(i)  θ  θ as = aθi [ as p(as |ai )]ki , where the set Ω(i) is the neighboring nodes of i, Ω is the set of all aθi s∈Ω(i)

as ∈Ω

nodes. p(as |ai ) is the conditional probability that a node with activity ai is connected to a node with activity as . In uncorrelated networks, F (as )N kas F (as )N m = F (as ), mN = mN F (as )N kas F (as )N mas as F (as ) = = mN mN ,

pI (as |ai ) = pII (as |ai ) =

noting Eqs. (3.18), (3.19) we have  θ I as p (as |ai )]kiI = aθi < aθ > m(1+ < a >) sIi = aθi [ as ∈Ω  θ II θ+1 θ as p (as |ai )]kiII = aθi > (1 + m < a >). sII i = ai [

(3.19)

(3.20)

as ∈Ω

Thus transition probability rate from i to j can be obtained Πi→j =

maθ+1 ai aθj mai ωij maj ωij j + , + = I II θ θ+1 N si N si N < a > (1+ < a >) N < a > (1 + m < a >)

(3.21)

in a similar way, we get Πj→i =

N<

aθ

maθ+1 aj aθi i + . θ+1 > (1+ < a >) N < a > (1 + m < a >)

Plugging (3.21) (3.22) into (2.8) and using

N  i=1

< i|ψ1 >=

(3.22)

< i|ψ1 >= 1 yields

maθ+1 aθi φ + N (1+m) N (1+) ai m 1+ + 1+m

,

(3.23)

m 1+m ),

(3.24)

with φ=

1−α1 α2 , N 

α1 = α2 =

i=1 N  i=1

maθ+1 ai i /( 1+ N (1+m) aθi

N (1+)

ai /( 1+ +

+

m 1+m ).

To support the results of the analytical treatment, we have performed the numerical simulation with N = 1000 nodes, m = 6, and activity distribution F (a) ∼ a−2 . All presented data in the following are obtained with double precision Fortran. In particular, we restricted the activity in the interval a ∈ [0.1, 1] to get great accuracy simulation. From Eq. (3.23), we see that the stationary probability is dominated by the weight parameter θ, differing substantially from the ones in static weighted network[22, 23], where it is a strictly linear increasing function of the node degree in log-log axis. The difference is due to the effects of the dynamical connectivity patterns that play a crucial role in the characterization of the stationary state. Left panel of Figure 1 shows the stationary probability < i|ψ1 > is a function of the activity for some values of θ. In the case of θ = 0.0, or 3.0, it increases with ai which means that it is easier to find the walker after a long time at a large-activity node than at a small-activity node. On the contrary, in the case of θ = −3.0, < i|ψ1 > is a decreasing function of ai , thus the final occupation probability for a node with a small activity is higher than another node with a large activity. Using (2.9) and (3.23) we get the coverage C(t) of time varying network, which is a function of time t (in Figure 1), averaged for different walks starting from different sources. In the cases of θ = −1.0, θ = 0.0, and θ = 1.0, it increases faster for negative θ than the case of nonnegative θ, which means weight parameter θ can adjust the searching efficiency in some extent. However, for large absolute values of θ, the searching efficiency is much poorer than the case of θ = 0.0 whether time is small or large. The reason is easy to understand. For small negative values of θ, i.e., θ = −1.0, aθ is a decreasing function of activity a (0 < a < 1), it reinforces the effect of the nodes with small activity which form a large group and weakens the effect of nodes with large activity. However, aθ increases with activity a increasing for positive values of θ, i.e., θ = 1.0, it enhances the accessibility between nodes with large activity by trapping the walker in a small group of adjacent high activity nodes, thus slows down the exploration of the network.

1.0

θ=-3.0 θ=0.0 θ=3.0

0.01

coverage

1E-3

0.6 0.4 0.2

(a) 0.0 1.0 0 0.8 1E-4

coverage

stationary probability

θ=0.0 θ=-1.0 θ=1.0

0.8

theory curve θ=-3.0 theory curve θ=0.0 theory curve θ=3.0

θ1000 =0.0 θ=-1.0 θ=1.0

2000

3000

4000

2000

3000

4000

0.6 0.4 0.2

1E-5 0.1

0.2

0.3

0.4

0.5

5000

(b)

0.6 0.7 0.8 0.9 1

0.0

activity

0

1000

t

5000

Figure 1: Plots of the stationary probability and the average coverage for undirect-link case. Left panel: the stationary probability as a function of activity for different values of θ, the solid and dotted lines represent the analytical prediction of Eqs. (3.23) and (3.24); Right panel: the average coverage is the function of activity for different values of θ, (a) is the numerical results of 1000 walkers in un-direct network obtained by integrating the activity-driven model for over 5000 time steps, and (b) is the analytical prediction of Eq. (2.9).

4

Results of direct time varying networks

There are two kinds of directed networks, that is, the active node i generates m-links which could be outgoing edges (Type I) or ingoing edges (Type II). If a walker can move to node j starting from node i just when i is active and creates outgoing link pointing to j (Type I) or j is active and creates ingoing link pointing to i (Type II). We only consider the results on uncorrelated networks, since F (a )N k

s as s )N m pout (as |ai ) = mN = F (amN = F (as ), kiout = m, a F (a )N k a F (as ) s s as s F (as )N m in in , ki = m < a >, p (as |ai ) = mN = mN = as

(4.25)

then transition probability rates are Πout i→j = Πin i→j =

mai ωij N sout i maj ωij N sin i

= =

ai aθj , N aθ+1 j m . N 1+m

Substituting (4.26) into (2.8) and using the normalization condition

(4.26) N  i=1

stationary probabilities < i|ψ1out >= < i|ψ1in >=

aθ−1 i , N θ+1 ai . N

< i|ψ1 >= 1, we are led to the

(4.27)

It reduces to the case of unbiased random walks if θ = 0, accordingly, Eq. (4.27) coincides with Eq. (7) in Ref.[12]. Furthermore, Eq. (4.27) also shows that the stationary distribution is dominated by the weight parameter θ. In type I network active nodes create outgoing links. Walkers are thus more likely to diffuse out of these nodes, meaning that the higher the activity, the smaller the stationary probability in that node. Through reinforcing the transition probability to a small group of large-activity nodes, i.e., setting θ = 1.0

0.01

0.01

stationary probability

stationary probability

θ=3.0 θ=1.0 θ=−1.0 θ=0.0

1E-3

(a) 1E-4 0.1

0.2

0.3

0.4

θ=3.0 θ=1.0 θ=−1.0 θ=0.0

1E-3

1E-4 0.1

0.5 0.6 0.7 0.80.9 1

0.3

0.4

0.5 0.6 0.7 0.80.9 1

0.4

0.5 0.6 0.7 0.80.9 1

activity

activity 0.01

0.01

θ=-3.0 θ=-1.0 θ=1.0

θ=-3.0 θ=-1.0 θ=1.0 staionary probability

staionary probability

(b) 0.2

1E-3

1E-3

1E-4

(c) 1E-4 0.1

0.2

0.3

activity

0.4

0.5 0.6 0.7 0.80.9 1

1E-5 0.1

(d) 0.2

0.3

activity

Figure 2: Plots of the stationary probability for two kinds of edge-linking. (a): Numerical simulations of 1000 walkers in out-direct network obtained by integrating the activity-driven model for over 200 time steps and averaged over 50 distinct realizations; (b): Analytical prediction of Eq. (4.27) for out-linking case; (c): Numerical simulations of 1000 walkers in in-direct network obtained by integrating the activitydriven model for over 100 time steps and averaged over 50 distinct realizations; (d): Analytical prediction of Eq. (4.27) for in-linking case.

1.0 θ=-3.0 θ=0.0 θ=1.0 θ=3.0

0.4

0.2

(a)

0.4

1000

2000

3000

4000

5000

θ=-3.0 θ=0.0 θ=1.0 θ=3.0

(b) 0.0 1000

0.4 0.2

2000

t

3000

4000

5000

(a) 1000

2000

3000

4000

3000

4000

5000

θ=−1.0 θ=0.0 θ=3.0 θ=−3.0

0.8

0.2

0

0.6

0.0 1.0 0

coverage

coverage

0.0 0.6 0

θ=3.0 θ=0.0 θ=−3.0 θ=−1.0

0.8

coverage

coverage

0.6

0.6 0.4 0.2

(b)

0.0 0

1000

2000

t

5000

Figure 3: Plots of the coverage for two kinds of edge-linking. Left panel: (a) Numerical simulations of 1000 walkers in out-direct network obtained by integrating the activity-driven model for over 5000 time steps; (b) Analytical prediction of Eq. (2.9) for out-linking case; Right panel: (a) Numerical simulations of 1000 walkers in in-direct network obtained by integrating the activity-driven model for over 5000 time steps; (b) Analytical prediction of Eq. (2.9) for in-linking case. makes the stationary distribution be uniform with < i|ψ1out >= N1 for node i. If θ > 1.0, < i|ψ1out > increases with ai , which means that it is easier to find the walker after a long time at a large-activity node than at a small-activity one. On the contrary, if θ < 1.0, < i|ψ1out > is a decreasing function of ai , thus the final occupation probability for a node with a small activity is higher than another node with a large activity (see Fig. 2(a) and Fig. 2 (b)). While in type II network, active nodes create ingoing links. Walkers are thus more likely to diffuse into high activity nodes, meaning that the higher the activity, the larger the stationary probability in that nodes. Through reinforcing the transition probability to a large group of low-activity nodes, i.e., setting θ = −1.0 makes the stationary distribution be uniform. < i|ψ1out > increases with ai for θ > −1.0 and decreases with ai for θ < −1.0 (see Fig. 2(c) and Fig. 2 (d)). Fig. 3 shows the coverage C(t) as a function of time t, averaged for different walks starting from different sources. θ = 1.0 and θ = −1.0 are the optimal search methods for type I and type II networks respectively, because C(t) can reach the equivalent values more earlier than that other cases of θ, From both panels of Fig. 3(a) and Fig. 3 (b), the numerical results are in good agree with the theoretical predictions. In terms of searching network efficiency, it is better for small |θ| compared to that for large |θ| in a long run. In order to study the transport dynamics in such networks, we mainly consider the average number of steps needed to reach or find a specific target, i.e., GMFPT. However, in type II network the eigenvectors of M can not span the whole space, but detailed balance condition < i|ψ1in > Mji =< j|ψ1in > Mij , ∀i, j = 1, 2, · · · N

(4.28)

is satisfied, thus the following matrix V can be easily verified to be symmetric in virtue of Eq. (4.28),  −1 1 < j|ψ1in > < i|ψ1in > (4.29) Vij =  Mij =  Mji = Vji , i.e., V = D 2 M D 2 , in in in < i|ψ1 > < i|ψ1 >< j|ψ1 >  1 where D 2 is a diagonal matrix with its ith diagonal entry equal to < i|ψ1in >. V is similar to M and thus has the same set of eigenvalues {λs } as M . Assume {|ϕs >} be the normalized and the mutually orthogonal eigenvectors of V , from Eq. (4.29) it is easy to verify that < i|ϕ1 >= < i|ψ1in >. (4.30)

Accordingly, the walker is at node i from node j at time t is  N < i|ψ1in > 

pij (t) = 

< j|ψ1in >

< i|ϕs >< j|ϕs > λts .

(4.31)

s=1

Using the previous approach we get GMFPT < Ti >= =

5

1 (1−)

1 (1−)

N  k=1,k=i N  k=1,k=i

< k|ψ1in > Tik < k|ψ1in >

N  s=2

√ 1 [< i|ϕs >2 − √ < k|ϕs >< i|ϕs >] 1−λ . in s

(4.32)



conclusions

In this paper, we have shown that dynamical topologies of network affect the behavior of a biased random walk performed on top of them at the same time scale, it is largely because of walker being trapped on the temporary isolated node until a new link involving this node occurs at the successive time step, thus it induced a slower dynamics, compared to quenched and annealed network. In the biased random walk, We consider the weight intensity to be ωij = (ki kj )θ ∝ (ai aj )θ when nodes i and j are connected, here θ is a control parameter of weights. For the special case of uncorrelated networks, we have calculated analytically the expressions of the stationary distribution, the coverage, GMFPT to a given target node in undirected and directed time varying network. The resultant formulas for all these interesting quantities are sensitive to the change of θ. These results are different from those in the case of the static networks. From the relation between C(t) and θ, we find that it is not the optimal search strategy for the case of θ = 0.0. Furthermore, it is also interesting for us to investigate the diffusion dynamics on other temporal network.

6

Acknowledgments

This work was jointly supported by NSFC/11331009, NSFC/1142600737, SHMEC/13YZ007 and STCSM/13ZR1416800, Natural Science research project of Henan Education department(Grant No.2011B110024).

References [1] J. D. Noh, H. Rieger, Phys. Rev. Lett. 92, 118701 (2004). [2] Z. G. Huang, X. J. Xu, Eur. Phys. J. B 51, 549 (2006). [3] V. Tejedor, O. B´ enichou, et al., Phys. Rev. E 80, R065104 (2009). [4] V. Tejedor, O. B´ enichou, et al.,Phys. Rev. E 83, 066102 (2011). [5] S. N. Dorogovtsev, A. V. Goltsev, et al., Rev. Mod. Phys. 80, 1275 (2008). [6] C. Castellano, R. Pastor-Satorras, Phys. Rev. Lett. 100, 148701 (2008). [7] M. Bogu˜ na ´, C. Castellano, et al, Phys. Rev. E 79, 036110 (2009). [8] P. Holme, J. Saram¨ aki, Phys. Rep. 519, 97 (2012). [9] A. L. Barab´ asi, Nature 435, 207 (2005). [10] J. G. Oliveira, A. L. Barab´ asi, Nature 437, 1251 (2005).

[11] N. Perra, B. Goncalves, et al., Sci. Rep. 2, 469 (2012). [12] N. Perra, A. Baronchlli, et al., Phys. Rev. Lett. 109, 238701 (2012). [13] A. Baronchelli, R. Pastor-Satorras, Phys. Rev. E 82, 011111 (2010). [14] B. Ribeiro, N. Perra, et al., Sci. Rep. 3, 3006 (2013). [15] A. Barrat, M. Barth´ elemy, et al., Proc. Natl. Acad. Sci. U. S. A. 101, 3747 (2004). [16] P. J. Macdonald, E. Almaas, et al, Euro. phys. Lett. 72, 308 (2005). [17] M. Starnini, A. Baronchelli, Phys. Rev. E 85, 056115 (2012). [18] E. W. Montroll, G. H. Weiss, J. Math. Phys. 6, 167 (1965). [19] A. Baronchelli, M. Catanzaro, et al., Phys. Rev. E 78, 011114 (2008). [20] L. E. Reichl, A modern course in statistical physics, John Wiley, New York, (1998). [21] H. S. Wilf, Generatingfunctionology, Academic Press, London (2nd edn.), (1994). [22] A. Fronczak, P. Fronczak, Phys. Rev. E 80, 016107 (2009). [23] Z. Z. Zhang, T. Shan, et al, Phys. Rev. E 87, 012112 (2013).