Diffusion on networked systems is a question of

Full text

Turn on search term navigation

ARTICLE

Received 4 Aug 2014 | Accepted 29 Apr 2015 | Published 9 Jun 2015

Network science investigates the architecture of complex systems to understand their functional and dynamical properties. Structural patterns such as communities shape diffusive processes on networks. However, these results hold under the strong assumption that networks are static entities where temporal aspects can be neglected. Here we propose a generalized formalism for linear dynamics on complex networks, able to incorporate statistical properties of the timings at which events occur. We show that the diffusion dynamics is affected by the network community structure and by the temporal properties of waiting times between events. We identify the main mechanismnetwork structure, burstiness or fat tails of waiting timesdetermining the relaxation times of stochastic processes on temporal networks, in the absence of temporalstructure correlations. We identify situations when ne-scale structure can be discarded from the description of the dynamics or, conversely, when a fully detailed model is required due to temporal heterogeneities.

DOI: 10.1038/ncomms8366

Diffusion on networked systems is a question of time or structure

Jean-Charles Delvenne1, Renaud Lambiotte2 & Luis E.C. Rocha2,3

1 ICTEAM and CORE, University of Louvain, 4 Avenue Lematre, Louvain-la-Neuve B-1348, Belgium. 2 Department of Mathematics and naXys, University of

Namur, 8 Rempart de la Vierge, Namur B-5000, Belgium. 3 Department of Public Health Sciences, Karolinska Institutet, 18A Tomtebodavagen, Stockholm

S-17177, Sweden. Correspondence and requests for materials should be addressed to J.-C.D. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366

The relation between network structure and dynamics has attracted the attention of researchers from different disciplines over the years15. These works are rooted in

the observation that in contexts as diverse as the Internet, society and biology, networks tend to possess complex patterns of connectivity, with a signicant level of heterogeneity4. In addition, a broad range of real-world dynamical processes, from information to virus spreading, is akin to diffusion. If the effect of structure, such as communities or degree heterogeneity, on diffusive processes is now well known6,7, the impact of the temporal properties of individual nodes is poorly understood1,8. Yet, empirical evidence indicates that real-world networked systems are often characterized by complex temporal patterns of activity, including a fat-tailed distribution of times916, correlations between events17,18 and non-stationarity1921. This is in remarkable contrast to a vast majority of mathematical models, which assume homogeneous interaction dynamics on networks1,2,22,23. In this work, we focus on the impact of structure and the waiting-time distribution on dynamics, and set aside any other temporal properties. This subject has triggered intense theoretical work in recent years, for example, in relation to anomalous diffusion, but predominantly on lattice-like, annealed24 and random structures20,23,2528. These limitations leave a fundamental question open, with important applications for the modelling of temporal networks: what are the effects of complex patterns, simultaneously in time and in structure, on the approach to equilibrium of diffusion processes?

To comprehend the interplay between temporal and structural patterns, we focus on a broad class of linear multi-agent systems describing N interacting nodes, dened by

Dx Lx 1 where xi, the ith component of x, represents the observed state of node i. The N N real matrix L encodes the mutual inuences in

the network, with non-zero entries indicating the presence of a link. D is either d/dt, the delay Dxi(t) xi(t 1), or any other

causal operator acting linearly on the trajectory of each entry xi(t). This equation couples network structure (represented by L)

and time evolution (represented by D) by describing a system where every node i has a state xi(t) Fui(t), where

ui(t) P

j Lijxj(t) represents the input, or inuence of the neighbouring states on node i. The operator F is the so-called transfer function29,30, dened as the inverse of D (see Methods).

A classic example is heat diffusion on networks, where every node has a temperature xi evolving according to the Fourier Law

m dx

dt Lx; 2 where m is the characteristic time of the dynamics and L is a

Laplacian of the network. The same set of equations can represent, possibly up to a change of variables, a basic model for the evolution of peoples opinions31, robots positions in the physical space30,32, approach to synchronization33,34 or the dynamics of a continuous-time random walker35our main example from now. In any case, the dynamical properties of the system are described by the spectral properties of the coupling matrix. The constraints imposed by the conservation of probability imply that the Laplacian dynamics is characterized by a stationary state, associated to the dominant eigenvector of L, which we will assume to be unique, as is the case in a large class of systems, for example, strongly connected networks. A key quantity is thus the second dominant eigenvalue, also called the spectral gap6, which provides us with the relaxation time to stationarity, usually called mixing time36 for stochastic diffusion processes. The spectral gap determines the speed of convergence to the stationary state and measures the effective size of the

system in terms of dynamics. The spectral gap is also related to important structural and dynamical properties of the system, such as the existence of bottlenecks and communities in the underlying network7,30.

In this study, we generalize the concept of spectral gap and of mixing time to random processes with general causal operator D, and focus in detail on operators with long-term memory, naturally emerging in diffusion with bursty dynamics. After showing connections between the theory of random walks and that of multi-agents systems, we identify the temporal and structural mechanisms driving the asymptotic dynamics of the system and provide examples when each mechanism prevails. By doing so, we show that the form of the temporal operator D may either slow down or accelerate mixing as compared with the differential operator equation (2). The results are further exploited to assess the possibility of coarse-graining the dynamics based on the network community structure and tested using numerical simulations on synthetic temporal networks calibrated with empirical data.

ResultsRandom walks with arbitrary waiting times. The generalized dynamics of the random walker, illustrated in Fig. 1a, is dened as follows. A walker arriving at a node i jumps towards a neigh-bouring node within a time interval [Dt,Dt d] with probability

r(Dt)d (for small d). In line with a standard discrete-time random walk process, the jump is directed towards a neighbour j with probability Pij. The probability density function r(Dt) is called the waiting- time distribution of the walker. At each hop, the waiting time Dt is reset to zero and consecutive waiting times are independent. The evolution of probability of the presence of a stochastic process in each state is ruled by the so-called master equation, or Kolmogorov forward equation, well known for this family of generalized walks24,37,38. We prefer to adopt here the equivalent, dual, viewpoint of Kolmogorov backward equation39, which belongs to the class of processes dened by equation (1) and can be analysed using the toolbox typical to multi-agent systems30,32, such as the transfer function formalism and eigenmode decomposition.

3 4

1 2 3 1

Time

Equilibrium

Homogeneous

Heterogeneous

4 4 1 1

12 3 1

4 41 1

Structure Timing

Figure 1 | Diffusion on temporal networks. (a) A random walk process can be illustrated by a letter (or banknote and so on) being randomly passed from neighbour to neighbour on a social network. Temporal patterns of waiting times between arrival and departure of the letter may be homogeneous (for instance, discrete or exponentially distributed times) or heterogeneous (for instance, bursts). (b) The relaxation time measures the characteristic time to reach equilibrium from any starting condition. (c) The competition between structure and temporal patterns regulates the relaxation time, or mixing time, of stochastic processes.

2 NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366 ARTICLE

We assign a xed real-valued observable xobs,i to every node i

and consider the value observed at time t by a walker starting initially from node j. This value is a random variable, with an expected value denoted by xj(t), taking initial value xj(0) xobs,j.

The random walk is ergodic and mixing on a strongly connected aperiodic network (aperiodicity is only required for the discrete time, when r(Dt) is a Dirac distribution). In this case, x(t)

describes a consensus dynamics, meaning that the individual values asymptotically converge to one another, xi(N) xj(N). If

we choose the observable xobs,j 0 at all nodes j, except l at node

k, then xi(t) is precisely the probability for the random walker starting at i to be found at k at time t. Therefore, the Kolmogorov backward equation, a consensus dynamics, embeds in particular the evolution of the probability of the presence on nodes.

The walker, assumed to have just hopped at time zero and nding itself at node i (another origin of time would not change the asymptotic decay times, of main interest in this work), hops again at time Dt with probability density r(Dt) and moves towards a neighbour j with probability Pij. The expected value observed by the walker at time t is xobs,i if it is still waiting for its

rst hop (toDt) and otherwise

Pj Pijxj(t Dt), by induction on the number of hops. Therefore, we obtain the vector equation

x t

xobs Z

which can be much higher than the average inter-contact time mcontact in case of bursts. This fact has been used in the literature to deduce that burst contact statistics slowdown diffusion in a complex network43. Nevertheless, it should not be confounded with the results presented in this paper, which will focus on the statistical properties of the waiting-time distribution and identify, among others, that its variance plays a signicant role on the asymptotic behaviour of the walker. In the above scenario, the variance of r(Dt) depends on the third moment of the inter-contact time distribution rcontact(t) and is associated to a

mechanism distinct from that of the bus paradox (equation (7)).

Eigenmodes for heterogeneous temporal operators. Equation (1) typically takes the form of an integro-differential equation. However, it simplies into the differential equation (2) in the case of a memoryless random walker. Memoryless refers to the case when the probability of hopping between Dt and Dt d, knowing

that the walker has waited at least Dt, is independent of Dt. This leads to an exponential waiting-time distribution45 r(Dt) e Dt/m/m, in other words an unconditional probability

e Dt/md/m of jumping between Dt and Dt d, for small d and for

mean waiting time m. In that case, we nd D(s) 1/r(s) 1 ms

in the Laplace domain, indeed recovering equation (2) in the time domain. The differential equation can then be analysed by changing the variables x to a linear combination of the eigenvectors vk of the Laplacian L, of eigenvalues

0 04l1,l2,...,lN 1Z 2 as follows30: x t

r Dt

1 dDt Z

Px t Dt

dDt; 3

with the discrete transition matrix P of entries Pij. The

convolution in time in the last term calls for a Laplace transform

x s

r Dt

x t

1 e stdt: 4

For simplicity, we use the same notations for functions in the time and in the Laplace domain, only distinguished by their variable, namely t or Dt for time and s for Laplace. This is justied as time and Laplace domain representations encode the same single physical object, for example, a probability or an observable. The same holds for operator D, thought of as acting in the time domain (for instance D d/dt) or the Laplace domain (D s)

according to the context.

The Laplace transform r(s) is the moment generating function of the waiting-time distribution r(Dt), as it encodes the moment in its Taylor series rs 1 ms m2 s2

s22 , where m is

the expected waiting time (rst moment), s2 is the variance and m2 s2 is the second moment. Using the fact that convolution

(respectively, integration from 0 to t) in the time domain corresponds to the usual product (respectively, division by s) in the Laplace domain40, equation (3) reduces to

x s

1 r s

s xobs r s

Px s

; 5

or equivalently

1 rs

x s

1s x t 0

Lx s

; 6

where we have made the dependence on the initial condition explicit by using the relation xobs x(t 0) and where L PI

denotes the (normalized) Laplacian of the network. This is an instance of equation (1), which shows that an inputoutput relationship is often best expressed in the Laplace domain rather than in the temporal domain. In this case, the transfer function F(s) is dened by the algebraic relation F 1(s) 1/r(s) 1 D(s), up to the initial condition term,

implicit in equation (1). See Methods for a derivation of equation (1) in a more general context.

From temporal networks to random walks. Diffusive processes often take place on temporal networks where individuals initiate from time-to-time short-lived contacts with their neighbours. A random walker can represent, for example, a letter or a banknote passing from hand to hand through rst contact initiated by the current node. The formalism described above focuses on the statistical properties of the waiting times of a walker on a node15,25, and not of the inter-contact times (the times between two subsequent contacts from a given node to another), as often considered in the literature20,27,28,4143. To illustrate this difference, let us consider an idealized scenario, where the network looks locally similar to a directed tree, to avoid indirect correlations due to cycles44, where the inter-contact time t between two contacts initiated by the same node is characterized by the same probability distribution rcontact(t), and where

activations on different edges are an independent random process. The corresponding waiting-time distribution r(Dt) for the random walker can be determined from rcontact45,46.

For example, the classic inspection paradox, or bus paradox, observes that the waiting time has a mean

m 1=2mcontact 1 s2contact

m2contact

zk t vk: 8

For simplicity we have supposed that the underlying network is undirected, connected and, in case of a discrete-time random walk, non-bipartite. Thus, the eigenvalues are real and the stationary state is uniquely dened6. Every zk(t), solution of

Dzk(t) lkzk(t), is called an eigenfunction of the operator D and

here takes the form of a decaying exponential. The problem is thus solved by decoupling structural and temporal variables, rst by identifying structural eigenvectors (vk), and then how they evolve in time (zk(t)). The resulting fundamental solutions zk(t)vk for equation (1) are called modes, or eigenmodes, of the system.

A similar analysis can be performed in the Laplace domain in the case of an arbitrary waiting-time distribution r(Dt), and thus when equation (1) has an arbitrary temporal operator D. In that case, the elementary solutions zk(t), associated to zk(s), need not

1 rs

NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366

be exponential functions and are obtained as solutions of equation (6), where the Laplacian is reduced to its eigenvalue lk

zks

1 1

lk 1=rs 1

zkt 0s ; 9

As before, any trajectory of the system can be expressed as a linear superposition of some or all the N modes.

Characteristic decay times of eigenmodes. Despite their non-exponential nature, a broad class of eigenfunctions dened by equation (9) still have a characteristic time tk describing the asymptotic decay of zk(t) to equilibrium as e t=tk (Fig. 1b). We nd that the decay time can be accurately estimated by performing a small s expansion in the Laplace domain (see Methods for the derivation and range of validity):

m lk

j j 1 b ;

where

2m2 11

is a measure of the burstiness of the temporal process based on the rst and second moments of its distribution. Burstiness b equals to zero for a Poisson process (memoryless waiting times, r(Dt) e Dt/m/m) and ranges from 1/2 (if r(Dt) is a Dirac

distribution) to arbitrarily large positive values (for highly bursty activity). This expression emerges naturally from the dynamical process and can be viewed as a measure of burstiness, equivalent to the commonly used burstiness measure47.

The estimate provided by equation (10) can be further tightened whenever the distribution r(Dt) of waiting times contains a fat tail, possibly softened by an exponential tail. The archetypical example is a power law with soft cutoff, rDt / Dt A ge Dt=ttail, a frequent model in human

dynamics supported by empirical evidence9,14,41. More generally, rDt / r0Dte Dt=ttail, where r0 is a fat tailed distribution, that is, decreasing subexponentially. The non-analytic point created in r(s) by the fat tail leads to an additional term in the characteristic time (see Methods)

tk max m lk

j j 1 b ;

ttail

which can be approximated as follows, if |lk| 1 and b have different orders of magnitude

max m lk

j j 1; mb; ttail

s2 m2

Mixing times and dominating mechanisms. Of particular importance is the slowest non-stationary mode (k 1), or mixing

mode, the characteristic decay time of which represents the worst-case relaxation time of any initial condition to stationarity. We call this time the mixing time of the process, in generalization of the classic memoryless case, where it is given by tmix mE 1,

determined by the so-called spectral gap E l140 (refs 6,36).

This quantity is related to the presence of bottlenecks (that is, weak connections between groups of highly connected nodes, a.k.a. network communities) in the network (via Cheegers inequality48,49). It is approximated as

tmix

max mE 1; mb; ttail

, comparing the exact mixing time (computed

with equation (18) in Methods), with what it would be with memoryless waiting times mE 1 , takes value in [0.6,1]. This

indicates a limited speed up of the mixing due to negative burstiness, bo0, whereas the structure plays a major role through

E traversing orders of magnitude.

On the other hand, competition between structure and time appears in scenarios of high temporal heterogeneity. For example, only strong communities are able to dominate power-law waiting times (Fig. 2e,f). The effect of structure on the mixing times is otherwise removed as burstiness (Fig. 2e) or tail (Fig. 2f) becomes the leading mechanism, scaling the mixing times up to 14-fold in the shown congurations. The transitions between the different mechanisms for a range of power-law congurations are presented in Supplementary Fig. 1 and Supplementary Note 1.

Model reduction. The use of coarse-graining through time-scale separation, which is the separate treatment of fast and slow dynamics that coexist inside a system, is crucial to reduce the complexity of systems made of a large number of interacting entities5052. This procedure is well known for differential equations such as equation (2). In this case, it consists in neglecting fast decaying modes, for example, with decay time less than a certain threshold ttreshan approximation invalid for early times but acceptable for times signicantly larger than ttresh. Only

the dominant modes, thus fewer variables, are left in the reduced model. Decreasing the threshold time, one produces a full hierarchy of increasingly more accurate models, but also with increasingly more variables. Reduced models have a clear interpretation from the structural point of view, as fast modes typically correspond to the fast convergence of the probabilities of the presence on nodes to a quasi-equilibrium within a network community. This process is followed by a slow equilibration of the population of random walkers trapped in each community to a global equilibrium5254. Each new reduced variable can therefore be interpreted as the slow-varying probability of the presence in each community, as if it had been collapsed into a single node. The hierarchy of increasingly more accurate and complex reduced models corresponds in this structural picture to a hierarchy of increasingly ner partitions into communities. Given that the decay times of the different modes in equation (2)

This expression shows that the asymptotic dynamics of each mode is determined by the competition between three factors: a structural factor mE 1 associated to the spectral properties of

the Laplacian of the underlying network and two temporal factors associated to the shape of the waiting-time distribution, namely its burstiness (mb) and its exponential tail (ttail). The slowest

(largest) of these factors dominates the asymptotic dynamics (Fig. 1c).

We emphasize that the burstiness and the fat tail effects are not necessarily related46. For instance, a power law r(Dt)p(Dt 1) 3 restricted to times DtrT (sharp cutoff), is

arbitrarily bursty but has no tail at all (ttail 0) for large but nite

T. On the other hand, the delayed power law r(Dt)p(Dt T) 4

restricted to times DtZT 1 has a fat tail (it decreases

subexponentially, ttail N) but low burstiness b for large T, as

the mean time m increases without bound and the variance remains constant. In the latter case, as for all pure power laws, r(Dt)BDt g for large Dt, the mixing time tmix, as ttail, is actually innite, reecting that mixing or relaxation to stationarity occurs only with subexponential convergence. In general, the properties of the tail of the distribution depend on the high-order moments of the distribution and not only on the rst two moments as captured in the coefcient of burstiness.

In Fig. 2, we study the mechanism dominating the mixing time in toy synthetic temporal networks with different waiting-time distributions. Not surprisingly, structure is the driving mechanism when waiting times are narrowly distributed around a mean value, as in the case of Dirac (discrete time, Fig. 2c) and Erlang distributions (resembling a discrete-time distribution with small uctuations, Fig. 2d). For those, the slowdown factor Y tmix= mE 1

4 NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366 ARTICLE

Strong communities Weak communities No communities

0 1

Increasing rewiring probability

Structure Tail

Burstiness

. Among those, any mode can only become negligible with respect to another after a very long time. For all practical purposes, the system exhibits an inherently complex dynamics with no timescale separation, as the ne-scale structure (the nest being at the level of a single node) and the large-scale network communities have equal or similar impact on the random walk dynamics. As a consequence, only two reduced models are available, one in which all N modes are considered and another described solely by the stationary state.

In general, only the Laplacian eigenvalues smaller in magnitude than m/ttail determine corresponding decay times. Just similar to that in classical memoryless case, they generate a hierarchy of timescales and reduced models corresponding to multiple levels of community structures. For instance, if ve eigenvalues are smaller than m/ttail, thus 0 l04l14?4l44 m/ttail, then

there is a reduced diffusion model capturing the movement of the random walker between ve aggregated nodes. The mixing, internal to each community and decaying as e t=ttail, is considered

instantaneous at large-enough timescales. Coarser aggregation, for example, based on two communities and two eigenvalues l0,

l1, may be relevant, although valid only at even larger times. However, aggregation based on any ner partitioning (other than into one-node communities) has the same accuracy as the ve-community model and thus little practical value as a reduced model. The degeneracy of characteristic decay times therefore limits the number of useful reduced models. This implies that a decomposition into communities is not necessarily associated to a timescale separation, or a reduced model, of the dynamics. For this reason, only models incorporating 15 or N-dominating modes are adequate for this particular example. This smaller choice of reduced models has ambivalent consequences. It limits the set of resolutions at which to describe the dynamics, but also provides a natural level for community structure (an open problem in multi-scale community detection5761), dened as the nest partitioning yielding a reduced model for the system.

Numerical analysis. Equation (14) shows that the mixing time depends in rst approximation only on the mean, variance and tail of the waiting-time distribution, whereas other properties of the distribution are irrelevant. With numerical simulations, we validate this approximation and provide some quantitative intuition on the competition between the three factors regulating the mixing. We construct synthetic temporal networks respecting our assumptions, for example, stationarity or the absence of correlations and calibrate them with the static structure and the inter-contact time distribution observed on a number of empirical data sets. The empirical networks used for calibration correspond to face-to-face interactions between visitors in a museum (SPM), between conference attendees (SPC)13 and between hospital staff (SPH)16; email communication within a university (EMA)10; sexual contacts between sex workers and their clients (SEX)14; and communication between members of a dating site (POK)12 (see Methods). On these networks, we observe the waiting-time distribution and the spectral gap, which allows to compute both the exact mixing time, with equation (18) in Methods, and its approximation in equation (14). The results, reported in Table 1, show a good agreement, except for SPC, as analysed later in this section.

Reference

0.6

Exponential

Discrete

Sharp power-law

Speedup

0.6

Speedup

0.6

Erlang law

Soft power-law

Slowdown

0 1

Increasing rewiring probability

Figure 2 | Network structure or waiting times regulate mixing times. We illustrate the ndings with a simple network composed of two communities, where the community structure is modulated by rewiring edges starting from two separated random graphs, as shown in a. Each panel (bf) shows (on the top bar) for each value of the rewiring probability the slowdown factor Y, ratio of the exact mixing time for a given waiting-time distribution

and for the exponential case (Poisson process). Each panel corresponds to a different waiting-time distribution and the colourbar scale on the right represents the values of Y for the respective colourmap. Y close to one

indicates that non-trivial waiting-time patterns can be neglected in the determination of the mixing time. On the other hand, a high Y indicates

that the details of the network structure are overshadowed by the temporal dynamics. (b) The reference case, that is, the exponential waiting times. We observe a speedup (Yo1) for (c) Dirac (discrete) and (d) Erlang

distributions, and a slowdown (Y41) for (e) a power law with sharp cutoff

and for (f) a power law with soft cutoff. From left to right, we show how this factor depends on the network structure by progressively removing the community structure (increasing E). The bottom bar of each panel shows,

according to equation (14), which is the dominant factor, either structure, burstiness or tail, regulating the relaxation time for each model of waiting times. See Methods for details on the synthetic networks and waiting-time distributions.

correspond to the Laplacian eigenvalues, the k-community partition is unsurprisingly found to be encoded into the k dominant Laplacian eigenvectors, in a way that is decoded by spectral clustering algorithms55,56.

This methodology can be extended to a more general temporal operators D, where the successive decay times now depend on the Laplacian eigenvalues and the temporal characteristics following equation (13). However, the effect of temporality may

have non-trivial consequences, as it may limit the number of reduced models. As an extreme scenario, let us consider a complex network in which the transient dynamics is dominated by the fat tail of the waiting-time distribution, such that all transient modes decay for large times as e t=ttail. Following

equation (8), the approach to stationarity is described by the superposition of modes, all of the form fkte t=ttailvk, for various

functions fk(t) that decay more slowly than exponentially, as for example fkt / t g

NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366

Except for SEX42 and POK data sets, in the other cases the temporal heterogeneity substantially increases the mixing times (see slowdown factor Y in Table 1). Table 2 shows that in these networks the dominant factor regulating the mixing time depends on the characteristics of the system. Fat tails of the waiting-time distribution drive the relaxation for the cases of face-to-face contacts (SPM, SPC and SPH). However, the structure is the leading mechanism behind the networks corresponding to other situations of human communication (EMA, SEX and POK).

If communities are completely removed by randomizing the network structure, the link sparsity of SEX and POK networks guarantees that structure remains dominating, as the sparsity results in the inevitable creation of bottlenecks for diffusion even in a random network. On the other hand, in the EMA network,

which has a relatively dense connected structure, the absence of communities leads to temporal dominance (Table 2). Finally, when the contact times are uniformly distributed, we recover the well-known result that network structure, with or without communities (Table 2), is the main factor regulating the convergence to stationarity.

As the raw empirical temporal networks do not necessarily abide by our simplifying assumptions, one cannot validate our theoretical derivations, for example, by estimating the mixing time directly from simulations on empirical data. Indeed the periodic rythms or correlations between successive jumps or other temporal patterns may induce effects not captured by our formula, as discussed in Supplementary Table 1 and Supplementary Note 3.

Empirical data are collected during a nite time span T, leaving the choice between a model for r(Dt) that is limited to T (sharp cutoff) or extrapolated to innite times with a tail decay (soft cutoff). The latter may occur, for instance, if we have theoretical or practical reasons to prefer a fat-tail-based model (for instance power law) with exponential cutoff for r(Dt). This choice is of little impact on the results, provided a sufcient observation time T (Supplementary Note 2). When temporal patterns trump structure, we have a self consistency condition tmix m E 1 b

ttail (Table 2) for SPM and SPH,

expressing that both models lead to similar mixing times, albeit through the different mechanisms of either burstiness (sharp) or fat tail (soft). For SPC, self-consistency is not attained. This happens possibly because the observation time of B2 days is not sufciently large to dilute the irregularity on the observed inter-contact time distribution (thus, in the simulated waiting times) induced by the inactivity during night periods. Consequently, this makes the tail difcult to measure, as it has not clearly emerged from a still transient behaviour, and a signicant mismatch is observed between exact and approximate mixing (Table 2).

We evaluate the possible sizes of reduced models for the data as follows. As the eigenvalues lk of the Laplacian range between 0 and 2, the number of modes with structurally determined

decay times ( m/ttailolkr0) can be roughly evaluated to

Nm/2ttail on an N-node network, if eigenvalues are evenly spaced. A similar reasoning holds whenever burstiness dominates the tail effect. This analysis reveals three different scenarios in the structure-dominated data sets considered above: (i) SEX benets from a full hierarchy of reduced models, as m/2ttail41. (ii) The

dynamics on EMA can be approximated with just three modes, associated to three network communities, as 0 l04l14l24

m/ttail4l34y, while a more detailed description, unless at the

node level, would not gain any signicant accuracy as further modes are all degenerate with the same decay time m/ttail.

(iii) POK exhibits a practically full hierarchy of reduced models, as around one third of its modes are determined by structure, m/2ttail 0.33; therefore, the nest reduced model is at a

few-node community level. This spectral analysis is comforted by applying a multiscale community detection algorithm to the empirical networks, which nds a partition into three communities dominated by fat-tail effects for EMA, but not for POK or SEX (Fig. 3).

DiscussionWe have presented a unied mathematical framework to calculate the relaxation time to equilibrium in a wide variety of stochastic processes on networked temporal systems. Our formalism is able to refer to arbitrary linear multi-agent complex systems, including linearizations of non-linear dynamical models such as Kuramoto oscillators or non-Laplacian diffusion dynamics such as Susceptible-Infected-Recovered epidemics, as detailed in Methods. It is

Table 1 | Exact versus approximate mixing times.

SPM SPC SPH EMA SEX POK Y 1.77 1.81 1.54 1.19 1.02 1.02

Exact mixing 54 205 671 718 7,887 29,347

Approximate Mixing 61 313 664 603 7,741 28,685

Relative Difference (%) 13.0 52.7 1.0 16 1.9 2.3

EMA, email communication within a university; SEX, sexual contacts between sex workers and their clients; POK, communication between members of a dating site; SPC, face-to-face interactions between conference attendees; SPH, face-to-face interactions between hospital staff; SPM, face-to-face interactions between visitors in a museum.

The table shows that for six congurations of temporal networks, the slowdown factor

Y tmix= mE

1 that is the ratio of the exact mixing time, calculated from equation (18) using the simulated waiting-time distribution, to the mixing time we would have for the exponential case (Poisson process, with same mean). The results indicate strong (SPM, SPC and SPH), medium (EMA) and weak (SEX and POK) slowdown (Y41). The approximate mixing time, computed with equation (14), shows a good agreement with the exact value.

Table 2 | Dominating mechanisms on temporal networks.

Original

structure

SPM 31

Randomized

structure

[afii9830]1 [afii9848]tail

[afii9830]1

61 24 23 43

[afii9848]tail

a b

SPC SPH EMA

SEX POK

SPM

SPC SPH EMA

SEX POK

105 377 203 1523 12524

3.6

Original

times

86 293

84 32 505

265 888 226

82 2124

c d

114 437 603 7741 28685

4.714.011.8 163 3447 3204

86 293

84 32 505

0.5 0.5 0.5 0.5 0.5 0.5

313 664 207

98 1755

1.89.37.316.140.8 133

Randomized

times

1310.254.9 678 1399

0.5 0.5 0.5 0.5 0.5 0.5

1.89.37.316.140.8 133

Here, a shows the regulating potential of each mechanism according to equation (14), either structure mE 1

burstiness (mb), or tail (t ), for the synthetic networks based on the empirical network structures and on the waiting times obtained by simulating a random walk on these structures. Either structure (brown) or tail (yellow) dominates the mixing time, also when structure or contact times are randomized independently (bd).

6 NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366 ARTICLE

Original network Community structure

Partition

EMA

100%

Structure

80%

60%

Burstiness

40%

Tail

20%

2 3

4 4

6 6 7 7 8 9 12 13

No. communities in partition

SEX

POK

100%

80%

F mechanism

60%

40%

20%

2 3

1 32

7 9 10 14 21 23 28 30

4 6

2 3

4 4

5 6 6 7 7 9 10 10

No. communities in partition

F mechanismF mechanism

Figure 3 | Global mixing time versus mixing times in communities.(a) Although the mixing time of diffusion on a network may be dened by the community structure, diffusion in the communities themselves may be regulated by structural or temporal patterns. (b) The original network constructed from an empirical data set is divided into communities with a multi-resolution community detection algorithm, for different values of the resolution parameter (see Methods). For each partition, we identify the fraction of communities with more than 10 nodes (y axis), where either structure or tail dominates (according to equation (14)), and plot against the number of communities detected at each partition (x axis). This conrms our spectral analysis, which predicts a compact description of the diffusion dynamics on EMA into three modes describing a probability ow between three communities, each dominated by fat-tailed waiting times. Also in agreement with the spectral viewpoint, this is not the case for SEX or POK, unless at the level of a few-node communities.

in other terms xit R

t0 uit0 xi;0dtt0dt0 (for the Dirac impulse d(t)), then in

the Laplace domain we nd xi(s) 1/s(ui(s) xi,0); in this case, the Laplace domain

operator is the multiplication by 1/s and vi is a Dirac impulse in the time domain, a constant in the Laplace domain that encodes the initial condition. One can as well invert the transfer operator, Di F 1i, and write dxi/dt ui(t) xi,0d(t) or

sxi(s) ui(s) xi,0; here, the Dirac impulse is seen as arising from differentiating the

discontinuous step of the state from rest to initial condition xi,0.

Connecting agents such that the input ui(t) of agent i is a weighted sum of other agents states,

Pj Lijxj(t) leads to a global dynamics dx/dt Lx x0d(t), or equivalently sx(s) Lx(s) x0, where L is the matrix of entries Lij, x is the vector of

states xi and x0 is the vector of initial conditions xi,0. When L is a Laplacian (rows summing to zero, non-negative off-diagonal weights), this is a simple consensus system where agents (for instance, robots or individuals) change their state (for instance, position, opinion modelled as a real number) towards an average of their neighbours states. With arbitrary interconnection weights and arbitrary transfer functions, one similarly obtains

Dx Lx v 15 where D is the diagonal matrix of transfer functions F 1i (D reduces to a scalar in case of identical agents Fi Fj) and v the vector of initial condition input vi. The

case when L remains Laplacian, but D is arbitrary, describes general, higher-order consensus systems where the agents converge to equal values on a strongly connected network through arbitrarily complicated internal dynamics, for example, modelling realistic robots or vehicles. For example, vehicles of mass m driven by a force and friction will obey mx a_

x Lx. An abuse of notation allows to drop

implicitly the initial condition and write simply dx/dt Lx and more generally

equation (1), Dx Lx, instead of equation (15) above. Discrete-time systems are

recovered with a discrete-derivative operator Dx(t) x(t 1) x(t).

Exact and approximate characteristic times for random walkers. We derive the approximate characteristic times (equation (12)) of decaying modes associated to a general random walk. The dominant mode associated to l0 0 is the stationary

also possible to accommodate non-uniformity of parameters as nodes are not identical in real systems. Our results are particularly relevant to improve the understanding of temporal networks, by highlighting the important interplay between structure and the temporal statistics of the network. Our formalism is different from previous studies of random walks on temporal networks that focus on homogeneous temporal patterns23 or do not account for the competition between structure and time15. We emphasize that questions related to ordering, not timing of events, such as the number of hops before stationarity or the succession of nodes most probably traversed by the walker, depend on the structure alone and not of hops timing statistics. Moreover, key mechanisms have been left aside in our modelling approach, such as the non-stationarity or periodicity of most empirical networks1921 and the existence of correlations between edge activations, and therefore preferred pathways of diffusion17,18,41. Important future work includes their incorporation in our mathematical framework and the identication of dominant mechanisms in empirical data.

We have shown that in the absence of some temporal correlations, the characteristic times of the dynamics are dominated either by temporal or by structural heterogeneities, as those observed in real-life systems. The competing factors are not only observed in different classes of networks modelled from empirical data but also at different hierarchical levels of the same network represented by its different scales of community

structure. We have also identied two contrasting metrics of the statistics of waiting times, burstiness and fat tails. We have shown that they regulate the dynamics on the network in different ways. In systems where temporal patterns are the dominating factor, the reduced models obtained by aggregation of communities as commonly used in practice are not necessarily relevant, as small-scale details are impenetrably intertwined with large-scale structure to form a complex global dynamics. In general, the temporal characteristics impose the natural description levels of the dynamics. Altogether, our results suggest the need for a critical assessment of a complexity/accuracy trade-off when modelling network dynamics. In some classes of real-world systems, the burden of increased model complexity may not compensate the incremental gain in knowledge, whereas other systems require the ne network structure as a key ingredient to any realistic modelling.

Methods

Overview. We provide a description of multi-agent linear dynamics, detail the derivation of approximation (12) and its range of validity, illustrate the generality of the framework beyond random walks and describe the empirical data and numerical implementations. In the Supplementary Information, we discuss on power laws and cut offs (Supplementary Note 1 and Supplementary Fig. 1), the evaluation of waiting-time distributions in empirical data (Supplementary Note 2) and the adequacy of our theoretical framework in empirical networks with correlations and non-stationarity (Supplementary Note 3).

Linear multi-agent systems. We derive equation (1) for general linear-networked systems, or multi-agent systems30, with an illustration on consensus dynamics. Every node i is modelled by a linear agent whose internal state is initially zero, and that converts an input signal ui with an operator Fithe so-called transfer functioninto a state signal xi Fiui. Here, ui represents a time trajectory ui(t) or

its Laplace transform ui(s), and similarly for the state trajectory xi. The transfer function Fi is an operator mapping input trajectories ui to state trajectories xi, which is requested to be linear, causal (if ui(t) 0 for all toT, then the same holds

for xi(t)) and time invariant (shifting ui(t) in time results in shifting xi(t)). Under those conditions, this operator takes the form xi(s) Fi(s)ui(s), a simple product of

functions, in the Laplace domain.

When using Laplace transforms, it is customary to use time domain trajectories that are zero for negative times. We account for the state initial condition xi(0) xi,0

with another input term vi that can be, for example, an impulse exciting the rest state to any desired initial condition at time zero. The agent dynamics therefore writes xi(s) Fi(s)(ui(s) vi(s)). For instance if Fi is the integration operator in

the time domain, xit xi;0 R

t0 uit0dt0 from the initial value at time 0, or

NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366

distribution, unique for an ergodic random walk. The relaxation time to the stationary distribution, generalizing the well-known mixing time of Markov chain theory36, is the characteristic time of the slowest decaying mode after the stationary mode, usually associated to l1, as we suppose here. An example of a case when l1 is not associated to the mixing mode, but rather lN, is the discrete random walk, where the distribution is a Dirac distribution at m, with zero variance, and the network is bipartite or close to bipartite.

The characteristic decay time of a time-domain function f t e t=t can

be found in the Laplace domain, as f(s) is dened and analytic over all s of real part larger than 1/t

[afii9845]approx

[afii9845](s)

(1+[afii9838]4)1

(1+[afii9838]3)1

(1+[afii9838]2)1

(1+[afii9838]1)1= (1+)1

(1+[afii9838]N1)1 (1+[afii9838]N2)1

decay, but undened or non-analytic in at least one point s of real part 1/t

decay. For example, the Laplace transform of e t=t is 1/(1 t

decay s),

decay, and the Laplace transform of t 1 ge t=t for

g41 is dened, but is non-analytic at s 1/t

decay.

with pole at s 1/t

[afii9848]tail

The eigenfunction zk(t) associated to lk is the solution ofD(s)zk(s) lkzk(s) D(s)zk(t 0)/s, derived in the text as equation (6), with

D(s) 1/r(s) 1 and L replaced with lk. The decay time of zk(t) is found if we nd

the right-most non-analytic point of

zks

1 Ds lk

[afii9848]

1 1 [afii9848]mix1

s zkt 0; 16

derived in the main text as equation (9). Note that for lka0, the singularity in s 0

is only apparent as D(s) ms (s2 m2)s2/2 y This expression is non-analytic

in s in exactly two circumstances: (i) D(s) is non-meromorphic in s, meaning that it cannot be dened analytically in a punctured neighbourhood of s (a neighbourhood of s s) or (ii) D(s) lk.

Case (i) happens when r(s) is non-meromorphic in s. A typical example is when rDt / r0Dte Dt=t , where r0(t) is a fat-tailed distribution, in other words

decreases more slowly than any exponential, leading to a non-meromorphic transfer functiona rare occurrence in standard systems theory. The non-analytic point is precisely s 1/ttail.

Case (ii) commands to solve the equation

rs 1= 1 lk

; 17 whose solution sk enters the expression for the exact characteristic decay time, combining the two cases:

tk max 1=sk; ttail

; 18 which can be approximated from an expansion of r(s). It is noteworthy that r(s) 1 ms (m2 s2)s2/2 ? is also called the moment-generating function,

whose successive derivatives around s 0 are the moments of the distribution, up

to sign. We consider the Pad (that is, rational) approximation:

2m s2 m2

2m s

Figure 4 | The exact and approximate characteristic times of a waiting-time probability distribution. Computation of the characteristic time of the distribution rs / r0te t=t (blue curve) with a fat-tailed r0

and a non-analytic point at s 1/t

tail. The exact characteristic decay times tk are found from the eigenvalues and the blue curve, following

equation (18). In this example, only three modes, including the stationary, are inuenced by the structure of the network, whereas the other modes collapse to a single decay time ttail. The dynamics can thus be

approximated by three modes, typically describing aggregate probability ows between three network communities. The approximate solution equation (12) follows the red curve, constructed from the Pad approximant rapprox, see equation (19). As eigenvalues lk change from 0 to 2, the

values 1/(1 lk) can take any real value 41 or o 1.

2 2 ; 19

equivalent for small s to a second-order Taylor approximation in terms of accuracy. This approximates the transfer function D 1(s) F(s) r(s)/(1 r(s)) as

1ms b; 20

where the adimensional term b s m2m takes its minimum value at 1/2 for the

discrete-time random walk, zero for a memoryless process and arbitrary large values for highly heterogeneous distributions. It is therefore a suitable measure of the burstiness of the process41,47. Such a transfer function is called Proportional-Integral, a common class in systems theory29. The equation D(s) lk is

approximately solved by 1/skEm(|lk| 1 b). A possible non-analytic point at

s 1/ttail caps the characteristic decay time of eigenfunction k to tk max m lk

j j 1 b ;

ttail

_ o X

Lijsin yj yi

;

If lk 1 and b are of different orders of magnitude, one may further approximate as

tk max m lk

j j 1 mb ;

ttail

in which the inuence of the fat tail, the structure and the burstiness are decoupled. The mixing time is associated with l1. Figure 4 shows that typically half of the modes have structurally determined decay times, namely those with positive 1 lk,

if eigenvalues are uniformly spaced between 0 and 2, in apparent contradiction

with the back of the envelope calculation Nm/2ttail proposed in the main text, which can reach N. This is a consequence of the progressive degradation of the Pad approximation for larger s (larger eigenvalues).

Even for small eigenvalues, for example, for l1 E, the approximation

(equation (12)), thus also equation (14), has a limited range of validity. The underlying Pad approximation (equation (19)) is valid for small-enough s, while for large-enough s the exact solution is ttail. The small s behaviour is captured by

the rst and second moments, and the large s behaviour, in other terms the tail ttail,

determines the growth of high-order moments. This double approximation based on high and low moments may explain why the approximation is successful on diverse data sets (see Fig. 3a). However, it is likely to fail precisely when the intermediate moments (third, fourth and so on), not covered by the approximation, dominate the behaviour of the moment generating function r(s).

This occurs for instance in the case of a power law of exponent 43 with sharp cutoff at a large time: although the rst and second moment remain bounded and

the tail is absent, the intermediate moments grow without bound with the cutoff time. However, such distributions are rarely used as models for real-life data. On the other hand, numerical experiments showed excellent accuracy for a wide range of power laws with soft cutoff.

Beyond random walks. Random walks are many times formally identical to the asymptotic behaviour of non-linear dynamical models and also serve as a prototype for a wider class of dynamics that includes epidemic spreading. For example, power networks may be modelled as a network of second-order Kuramoto oscillators33 whose state in every node is a phase (an angle) yi inuenced by its neighbours

Iyi a

where I is the inertia, a the friction and o describes the natural frequency of the oscillation. In this construction, we assume identical parameters for each node. After the change of variables fi yi ot/a, the linearization around the

synchronized equilibrium (fi fj is constant for all times) is writtenID

fi aD

fi X

LijDfj; 22

which is formally identical to a consensus equation structured by the Laplacian L33. Hence, it supports the same formalism as random walks, after the operator D is changed accordingly. The linear dynamics associated to the asymptotic behaviour is often a crucial rst step in understanding the non-linear dynamics of the network at some time scale.

Another class of dynamical systems where our results apply are linear dynamics given by equation (1), but where the interaction matrix L does not have a Laplacian structure and where the dominant eigenmode is not necessarily stationary. This is the case, for instance, for epidemic processes where infected individuals diffuse on a meta-population network, where nodes are large populations (cities and so on), and the individuals may either reproduce (by contamination of a healthy individual) or die (or recover). One classic model for this process is a multi-type branching process. This is akin to a random walk whose number of walkers is not preserved in time, as random events are not limited to hops but also death or split. This kind of model also emerges as the linearization of classic compartmental epidemic models such as Susceptible-Infected or Susceptible-Infected-Recovered62. The waiting time between two consecutive events may also be described by an

8 NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366 ARTICLE

arbitrary distribution r(Dt). The expected number of walkers in time may again be described by an equation Dx Lx, where D depends on r in the same way as

before, but now L is no more a Laplacian and in particular may have negative and positive eigenvalues. If the epidemics is supercritical, then the system is unstable. This implies the existence of unstable eigenmodes (lk40) where the number of walkers grow exponentially as et=t , with tk now approximated by m l 1k b

In this case, the fat tail of the waiting-time distribution does not play any role for the unstable eigenmodes. From this formula, we see that although the decay of stable eigenmodes tends to slow down due to burstiness and possibly fat tail, the unstable eigenmodes are boosted by burstiness. As a consequence, it leads in all eigenmodes to a long-term increase of the infected population due to the temporal heterogeneity.

Numerical analysis. We use six data sets of empirical networks: SPM, SPC13 and SPH16; EMA10; SEX14; and POK12 (see Table 3).

The synthetic networks used in Fig. 2 are formed by 1,000 nodes and 10,000 links equally divided into two groups, initially disconnected (rewiring probability equal to zero). Within each group, pairs of nodes are formed between uniformly chosen nodes. A fraction of links are rewired to weaken communities. Rewiring consists of uniformly choosing two pairs of nodes and swapping the pairs. Rewiring probability equal to 1 removes any signicant community structure from the network. Communities refer to groups of nodes with more connections between themselves than with nodes at other groups. The spectral gap is calculated for each network conguration. The exemplied waiting-time distribution arer(Dt) exp( Dt/m)/m with m 1 (exponential); 1.0 (discrete); pDtk 1

exp( kDt/m) with k 3 and m 3 (Erlang law); pDt g with g 2 and cutoff at

T 1,000 (power law with sharp cutoff); pDt gexp( Dt/ttail) with g 3 and

ttail

Table 3 | Summary statistics of the empirical networks.

SPM SPC SPH EMA SEX POK N 72 113 75 3,186 11,416 28,295

L 6,980 20,818 32,424 309,120 33,645 528,869 dt min min min h day h

T (days) 1 2.5 4 82 1,000 512 m 15.4 80.5 288 61.5 92.9 1,205

s 30.8 143 502 119 121 1,634

E 0.505 0.709 0.66 0.102 0.012 0.042

ttail 61 313 664 207 98 1,755

For the largest connected component of the empirical networks used in this study, we display the number of vertices (N), number of links (L), temporal resolution used in the simulations (dt), total duration of the network data (T), average (m) and s.d. (s) of waiting times for the simulated random walk, in the same unit as dt, and the spectral gap of the original unweighted network E.

20 (power law with soft cutoff).

The synthetic temporal networks used in Table 2 and Fig. 3 are constructed by using the unweighted version of the empirical networks as underlying xed structures. The randomized networks used in Table 2 are obtained by uniformly selecting two pairs of nodes of the original network and then swapping the respective contacts. In the dynamic network, a node and its links alternate between inactive and active states. To t this synthetic temporal network to our theoretical assumptions, we assume that all links of a node are activated simultaneously, the system is time invariant (no daily or weekly cycles), exhibits no temporal correlations between successive jumps and node-activation times are sampled from the same distribution of interactivation times (node homogeneity). The active state of a node and its links lasts for one time step dt, and then the node and links return to the inactive state. The distribution of interactivation times (a.k.a. inter-contact times) corresponding to the original times is obtained by pooling the times between two subsequent activations of the same node, in other words two consecutive contacts established between the node and its neighbourhood, as observed in the empirical networks20,42. The randomized times in Table 2 correspond to activation times sampled from exponential distributions with the same mean as the corresponding empirical cases. To obtain the waiting-time estimates presented in Table 2, we simulate a random walk in these synthetic networks. A random walker starts in a randomly chosen node. As time goes by, a walker remains in the node until its new activation and then hops to a uniformly chosen neighbour. The waiting time is thus the time elapsed between the arrival and departure of the walker in a node. We simulate a single random walker and let it hop 300,000 times, starting at 10 uniformly chosen nodes; hence, we collect a total of 3 million points for the statistics of waiting times. Besides that, for each starting conguration we discard the initial 5,000 hops to remove the stochastic transient. The waiting-time distribution observed in the simulations provides us with estimates of tmix, m, s2

and ttail. The mixing time tmix is estimated from the distribution expressedin the Laplace domain through equation (18) (for k 1, the mixing mode). This

equation is exact, as the synthetic data has been constructed so as to satisfy the conditions of validity under which it is proved. One could also estimate the mixing time from the convergence of random trajectories starting in a denite node towards stationarity; however, this exercise is computationally demanding. The exponential tail decay time ttail is estimated using least-square tting by a line in a lin-log plot of the second half of the waiting-time distribution data, thus between Tmax/2 and Tmax, where Tmax is the largest observed waiting time. Indeed, an exponential decay translates into a straight line in a lin-log plot, the slope of which determines the decay time. The choice of the interval [Tmax/2, Tmax] to perform a

linear t is steered by the want for a simple criterion uniform across data sets. This may be inappropriate for data sets where the observation time forces a cutoff on the data sets before the tail decay has had time to emerge, for instance in SPC, where a more careful denition of tail, for example, on a smaller interval, may be more adequate.

Let us interpret some elementary observations that can be made on Table 2. From (a) to (b), ttail is only approximately preserved, because it is recomputed from the simulated waiting-time distribution on a different network, thus a slightly different distribution, even though the inter-contact times are sampled from the same distribution in both panels. In (c) and (d), the inter-contact time distribution is memoryless (exponential), thus identical to the waiting-time distribution regardless of the network, and ttail coincides with m. From (a) to (c), the quantity mE 1 drops sharply, although E is clearly preserved as the network is unaltered. The reason is that although the mean inter-contact time is preserved by construction,

the mean waiting time also depends on the variance of the inter-contact times, in virtue of the bus paradox, and this variance is clearly strongly affected by the time randomization. The drop in mE 1 from (b) to (d) is associated to the same mechanism.

To identify the community structure of the empirical networks in Fig. 3, we use the partition stability method54 using a freely available implementation63. The method uses a simple diffusion process such that different communities are detected at different time scales according to the potential of the communities to trap the diffusion at the given time scale. Optimal communities have been derived for values of the resolution parameter decreasing from 102,101.95,101.9,101.85,y, with increasingly ne partitions. After identifying the relevant communities, we discard the inter-community links and calculate the spectral gap of each individual community of at least ten nodes.

References

1. Bansal, S., Read, J., Pourbohloul, B. & Meyers, L. A. The dynamic nature of contact networks in infectious disease epidemiology. J. Biol. Dyn. 4, 478489 (2010).2. Barrat, A., Barthlemy, M. & Vespignani, A. Dynamical Processes on Complex Networks (Cambridge University Press, 2012).

3. Moody, J. The Oxford Handbook of Analytical Sociology. pages 447474 (Oxford University Press, 2009).

4. Newman, M. E. J. Networks: An Introduction (Oxford University Press, 2010).5. Vespignani, A. Modelling dynamical processes in complex socio-technical systems. Nat. Phys. 8, 3239 (2012).

6. Chung, F. R. K. Spectral Graph Theory (American Mathematical Society, 1996).7. Lovsz, L. Random walks on graphs: a survey. Boy. Soc. Math. Stud. 2, 146 (1993).

8. Holme, P. & Saramaki, J. Temporal networks. Phys. Rep. 519, 97125 (2012).9. Barabsi, A. -L. The origin of bursts and heavy tails in human dynamics. Nature 435, 207211 (2005).

10. Eckmann, J. -P., Moses, E. & Sergi, D. Entropy of dialogues creates coherent structures in e-mail trafc. Proc. Natl Acad. Sci. 101, 1433314337 (2004).11. Haerter, J. O., Jamtveit, B. & Mathiesen, J. Communication dynamics in nite capacity social networks. Phys. Rev. Lett. 109, 168701 (2012).

12. Holme, P., Edling, C. R. & Liljeros, F. Structure and time-evolution of an Internet dating community. Soc. Net. 26, 155174 (2004).

13. Isella, L. et al. Whats in a crowd? Analysis of face-to-face behavioral networks.J. Theor. Biol 271, 166180 (2011).14. Rocha, L. E. C., Liljeros, F. & Holme, P. Information dynamics shape the sexual networks of Internet-mediated prostitution. Proc. Natl Acad. Sci. 107, 57065711 (2010).

15. Starnini, M., Baronchelli, A., Barrat, A. & Pastor-Satorras, R. Random walks on temporal networks. Phys. Rev. E 85, 056115 (2012).

16. Vanhems, P. et al. Estimating potential infection transmission routes in hospital wards using wearable proximity sensors. PLOS One 8, e73970 (2013).

17. Rosvall, M., Esquivel, A. V., Lancichinetti, A., West, J. D. & Lambiotte, R. Memory in network ows and its effects on spreading dynamics and community detection. Nat. Commun. 5, 4630 (2014).

18. Scholtes, I. et al. Causality-driven slow-down and speed-up of diffusion in non-Markovian temporal networks. Nat. Commun. 5, 5024 (2014).

19. Holme, P. & Liljeros, F. Birth and death of links control disease spreading in empirical contact networks. Sci. Rep. 4, 4999 (2014).

20. Rocha, L. E. C. & Blondel, V. D. Bursts of vertex activation and epidemics in evolving networks. PLoS Comput. Biol. 9, e1002974 (2013).

NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8366

21. Horvth, D. X. & Kertsz, J. Spreading dynamics on networks: the role of burstiness, topology and non-stationarity. N. J. Phys. 16, 073037 (2014).

22. Malmgren, R. D., Stouffer, D. B., Motter, A. E. & Amaral, L. A. N. A Poissonian explanation for heavy tails in e-mail communication. Proc. Natl Acad. Sci. 105, 1815318158 (2008).

23. Perra, N. et al. Random walks and search in time-varying networks. Phys. Rev. Lett. 109, 238701 (2012).

24. Klafter, J. & Sokolov, I. M. First Steps in Random Walks: From Tools to Applications (Oxford University Press, 2011).

25. Iribarren, J. L. & Moro, E. Impact of human activity patterns on the dynamics of information diffusion. Phys. Rev. Lett. 103, 038702 (2009).

26. Jo, H. -H., Perotti, J. I., Kaski, K. & Kertsz, J. Analytically solvablemodel of spreading dynamics with non-Poissonian. Phys. Rev. X 4, 011041 (2014).

27. Min, B., Goh, K. -I. & Vazquez, A. Spreading dynamics following bursty human activity patterns. Phys. Rev. E 83, 036102 (2011).

28. Vazquez, A., Rcz, B., Lukcs, A. & Barabsi, A. -L. Impact of non-Poissonian activity patterns on spreading processes. Phys. Rev. Lett. 98, 158702 (2007).

29. Astrm, K. J. & Murray, R. M. Feedback Systems: An Introduction for Scientists and Engineers (Princeton University Press, 2010).

30. Fax, J. A. & Murray, R. M. Information ow and cooperative control of vehicle formations. Autom. Cont. 49, 14651476 (2004).

31. Blondel, V. D., Hendrickx, J. M., Olshevsky, A. & Tsitsiklis, J. N. Convergence in multiagent coordination, consensus, and ocking. Proc. 44th IEEE Conf. Decision Control, 2996-3000 (2005).

32. Jadbabaie, A., Lin, J. & Morse, A. S. Coordination of groups of mobile autonomous agents using nearest neighbor rules. Autom. Control IEEE Trans. 48, 9881001 (2003).

33. Drer, F. & Bullo, F. Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators. SIAM J. Cont. Opt. 50, 16161642 (2012).

34. Strogatz, S. H. From Kuramoto to Crawford: exploring the onset of synchronization of in populations of coupled oscillators. Phys. D 143, 120 (2000).

35. Yin, G. G. & Zhang, Q. Continuous-Time Markov Chains and Applications: A Two-time-scale Approach volume 37 (Springer, 2012).

36. Levin, D. A., Peres, Y. & Wilmer, E. L. Markov Chains and Mixing Times (American Mathematical Society, 2008).

37. Hoffmann, T., Porter, M. A. & Lambiotte, R. Generalized master equations for non-Poisson dynamics on networks. Phys. Rev. E 86, 046102 (2012).

38. Montroll, E. W. & Weiss, G. H. Random walks on lattices ii. J. Math. Phys. 6, 167181 (1965).

39. Kolmogoroff, A. Uber die analytischen Methoden in der Wahrscheinlichkeitsrechnung. Math. Ann. 104, 415458 (1931).

40. Dyke, P. An Introduction to Laplace Transforms and Fourier Series (Springer, 2014).

41. Kivela, M. et al. Multiscale analysis of spreading in a large communication network. J. Stat. Mech. 2012, P03005 (2012).

42. Rocha, L. E. C., Liljeros, F. & Holme, P. Simulated epidemics in an empirical spatiotemporal network of 50,185 sexual contacts. PLoS Comput. Biol. 7, e1001109 (2011).

43. Karsai, M. et al. Small but slow world: how network topology and burstiness slow down spreading. Phys. Rev. E 83, 025102 (2011).

44. Speidel, L., Lambiotte, R., Aihara, K. & Masuda, N. Steady state and mean recurrence time for random walks on stochastic temporal networks. Phys. Rev. E 91, 012806 (2015).

45. Kleinrock, L. Queueing Systems: Theory. I (Wiley-Interscience, 1975).46. Lambiotte, R., Tabourier, L. & Delvenne, J. -C. Burstiness and spreading on temporal networks. Eur. Phys. J. B 86, 7 (2013).

47. Goh, K. I. & Barabsi, A. -L. Burstiness and memory in complex systems. Europhys. Lett. 81, 4 (2008).

48. Cheeger, J. Problems in Analysis. pages 195199 (Princeton University Press, 1970).

49. Diaconis, P. & Stroock, D. Geometric bounds for eigenvalues of Markov chains. Ann. Appl. Probab. 1, 3661 (1991).

50. Gfeller, D. & De Los Rios, P. Spectral coarse graining and synchronization in oscillator networks. Phys. Rev. Lett. 100, 174104 (2008).

51. Kokotovic, P., Khalil, H. K. & Oreilly, J. Singular Perturbation Methods in Control: Analysis and Design volume 25 (Society for Industrial and Applied Mathematics, 1987).

52. Simon, H. A. & Ando, A. Aggregation of variables in dynamic systems. Econometrica 29, 111138 (1961).

53. Simonsen, I. Diffusion and networks: a powerful combination! Phys. A 357, 317330 (2005).

54. Delvenne, J. -C., Schaub, M. T., Yaliraki, S. N. & Barahona, M. Dynamics On and Of Complex Networks volume 2 pages 221242 (Springer, 2013).

55. Shen, H. -W. & Cheng, X. -Q. Spectral methods for the detection of network community structure: A comparative analysis. J. Stat. Mech. 2010, P10020 (2010).

56. Von Luxburg, U. A tutorial on spectral clustering. Stat. Comput. 17, 395416 (2007).

57. Reichardt, J. & Bornholdt, S. Statistical mechanics of community detection. Phys. Rev. E 74, 016110 (2006).

58. Ronhovde, P. & Nussinov, Z. Multiresolution community detection for megascale networks by information-based replica correlations. Phys. Rev. E 80, 016109 (2009).

59. Fortunato, S. Community detection in graphs. Phys. Rep. 486, 75174 (2010).60. Karrer, B., Levina, E. & Newman, M. E. J. Robustness of community structure in networks. Phys. Rev. E 77, 046119 (2008).

61. Delmotte, A., Tate, E. W., Yaliraki, S. N. & Barahona, M. Protein multi-scale organization through graph partitioning and robustness analysis: Application to the myosin light chain interaction. Phys. Biol. 8, 055010 (2011).

62. Keeling, M. J. & Ross, J. V. On methods for studying stochastic disease dynamics. J. R. Soc. Interface 5, 171181 (2008).

63. Delmotte, A., Schaub, M. T., Yaliraki, S. N. & Barahona, M. Community Detection using the stability of a graph partition https://github.com/michaelschaub/PartitionStability

Web End =https://github.com/ https://github.com/michaelschaub/PartitionStability

Web End =michaelschaub/PartitionStability (2012).

Acknowledgements

We thank Michael Schaub for carefully reading the manuscript. J.C.D. and R.L. thank nancial support of IAP DYSCO and ARC Mining and Optimization of Big Data Models. R.L. acknowledges support from FNRS and from FP7 project Optimizr.L.E.C.R. is an FNRS charg de recherches and thanks VR for nancial support.

Author contributions

J.C.D., L.E.C.R. and R.L. conceived the project and wrote the manuscript. J.C.D. derived the analytical results. L.E.C.R. performed the numerical simulations and analysed the data.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Delvenne, J.-C. et al. Diffusion on networked systems is a question of time or structure. Nat. Commun. 6:7366 doi: 10.1038/ncomms8366 (2015).

10 NATURE COMMUNICATIONS | 6:7366 | DOI: 10.1038/ncomms8366 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

Word count: 11398

Show less

Abstract

Translate

Network science investigates the architecture of complex systems to understand their functional and dynamical properties. Structural patterns such as communities shape diffusive processes on networks. However, these results hold under the strong assumption that networks are static entities where temporal aspects can be neglected. Here we propose a generalized formalism for linear dynamics on complex networks, able to incorporate statistical properties of the timings at which events occur. We show that the diffusion dynamics is affected by the network community structure and by the temporal properties of waiting times between events. We identify the main mechanism--network structure, burstiness or fat tails of waiting times--determining the relaxation times of stochastic processes on temporal networks, in the absence of temporal-structure correlations. We identify situations when fine-scale structure can be discarded from the description of the dynamics or, conversely, when a fully detailed model is required due to temporal heterogeneities.

Details

Title

Diffusion on networked systems is a question of time or structure

Author

Delvenne, Jean-charles; Lambiotte, Renaud; Rocha, Luis E C

Pages

7366

Publication year

2015

Publication date

Jun 2015

Publisher

Nature Publishing Group

e-ISSN

20411723

Source type

Scholarly Journal

Language of publication

English

DOI

https://doi.org/10.1038/ncomms8366

ProQuest document ID

1686754602

Diffusion on networked systems is a question of time or structure

Jump to:

Full text

Abstract

Details

Suggested sources