Generative models of simultaneously heavy-tailed distributions of inter-event times on nodes and edges
GGenerative models of simultaneously heavy-tailed distributionsof inter-event times on nodes and edges
Elohim Fonseca dos Reis, Aming Li,
2, 3 and Naoki Masuda
1, 4, ∗ Department of Mathematics, State University of New York at Buffalo, Buffalo, New York, USA Department of Zoology, University of Oxford, Oxford, UK Department of Biochemistry, University ofOxford, Oxford, UK Computational and Data-Enabled Science and Engineering Program,State University of New York at Buffalo, Buffalo, New York, USA (Dated: September 1, 2020)Intervals between discrete events representing human activities, as well as other types of events,often obey heavy-tailed distributions, and their impacts on collective dynamics on networks such ascontagion processes have been intensively studied. The literature supports that such heavy-taileddistributions are present for inter-event times associated with both individual nodes and individualedges in networks. However, the simultaneous presence of heavy-tailed distributions of inter-eventtimes for nodes and edges is a non-trivial phenomenon, and its origin has been elusive. In the presentstudy, we propose a generative model and its variants to explain this phenomenon. We assume thateach node independently transits between a high-activity and low-activity state according to acontinuous-time two-state Markov process and that, for the main model, events on an edge occurat a high rate if and only if both end nodes of the edge are in the high-activity state. In otherwords, two nodes interact frequently only when both nodes prefer to interact with others. Themodel produces distributions of inter-event times for both individual nodes and edges that resembleheavy-tailed distributions across some scales. It also produces positive correlation in consecutiveinter-event times, which is another stylized observation for empirical data of human activity. Weexpect that our modeling framework provides a useful benchmark for investigating dynamics ontemporal networks driven by non-Poissonian event sequences.
I. INTRODUCTION
Dynamics contacts as well as the static structure ofsocial contact networks govern how humans or animalsgather, communicate, and act. Many techniques fromtemporal networks have been proven useful for describ-ing and utilizing data of time-varying networks [1–5].The time between two consecutive contacts, called theinter-event time (IET), is a key quantity to character-ize temporal networks and dynamics on them. Myriadhuman activities, such as online chats, email correspon-dence, mobility, web browsing, and broker trading, haveheavy-tailed distributions of IETs [1, 4, 6, 7]. This obser-vation implies that sequences of discrete events that anindividual node or edge in a network experiences obeysnon-Poissonian statistics. By contrast, in most cases,stochastic processes on static networks implicitly assumethat events such as infection or broadcasting occur ac-cording to Poisson processes, with which IETs obey anexponential distribution. Therefore, the non-Poissoniannature of event sequences in empirical data inevitablyurges us to reconsider our understanding of stochastic dy-namical processes on networks. In fact, effects of heavy-tailed distributions of IETs on epidemic processes [8–15],opinion dynamics [16–19], evolutionary game dynamics[20], and cascade processes [21–24], random walks [25–28], to name a few, have been studied. There are also anumber of generative mechanisms and descriptive models ∗ naokimas@buffalo.edu for heavy-tailed distributions of IETs including priorityqueuing models [6, 7, 29–34], mixture of exponentials [35–37], those supplied by circadian and weekly rhythms [38],and other self-exciting processes [39, 40].Heavy-tailed distributions of IETs are commonly foundfor single nodes [6, 7, 37, 38, 41–43] and single edges[1, 4, 9, 44, 45]. Such a distribution for a node impliesthat the sequence of event times for the node, regardlessof the identity of the neighbor, obeys a non-Poissonian,heavy-tailed statistics. In fact, the distribution of IETsfor both nodes and edges in a single data set are oftenheavy-tailed (see Section II for examples). However, thepresence of heavy-tailed IET distributions for both edgesand nodes in the same network is not trivial. Considera node, denoted by v , that interacts with its k neigh-bors, and assume that the sequence of IETs on each edgeincident to v independently obeys a heavy-tailed distri-bution. The superposition of the k event sequences yieldsthe sequence of v ’s events. This situation is illustrated inFig. 1. In Fig. 1(a), node v has k = 5 neighbors. We haveproduced the sequence of events on each of the five edgesassuming a power-law distribution of IETs, as shown inFig. 1(b). The sequence of events shown in the bottomof Fig. 1(b) is that for v , which one obtains by superpos-ing the sequences of events on the k edges. In general,the distribution of IETs for v is less heavy-tailed thanthat for a single edge. This is because the superpositionof independent point processes (more precisely, renewalprocesses) obeying heavy-tailed distributions roughly ap-proaches, albeit not precisely, a Poisson process as oneincreases k [46, 47]. a r X i v : . [ phy s i c s . s o c - ph ] A ug FIG. 1. Schematic illustration of the superposition of event sequences on edges. (a) Star network of six nodes. (b) Sequenceof events generated by a power-law distribution of IET on each edge of the star network and the sequence of events on node1 in (a). The CV of IETs, which we calculate from the first approximately 10 events on each edge, is shown in (b) for thefive edges and node 1. We used the power-law distribution of IETs for edges given by p ( τ ) = ( α − / (1 + τ ) α , where α = 3 . p w ( t ) = ( α − / (1 + t ) α − . Recently, we proposed a model that generates heavy-tailed distributions of IETs for both edges and nodes [49].The model assumes that each node is activated at dis-crete times according to a renewal process that drawsthe inter-activation time from a power-law distribution.Then, at each time step, activated nodes are uniformlyrandomly selected to communicate with simultaneouslyactivated neighbors. By construction, this model pro-duces a heavy-tailed distribution of IETs for individualnodes. Although it is less trivial, IETs for edges also obeyapproximately heavy-tailed distributions. However, thismodel does not explain why we find heavy-tailed distri-butions of IETs for both nodes and edges in empiricaldata. Realistic mechanisms of the simultaneous presenceof heavy-tailed IET distributions for nodes and edges areunderexplored [4].In the present study, we propose a model of time-stamped event sequences on networks that is based ona latent state dynamics of nodes. In our model, eachnode switches between two states called the high-activityand low-activity states according to a Markov process incontinuous time. We assume that if both nodes are inthe high-activity state, events occur according to a Pois-son process at a higher rate, otherwise at a lower rate.The rationale behind the model is that events betweentwo individuals may happen more frequently when bothindividuals are motivated to interact with others thanotherwise. We show that our model produces distribu-tions of IETs with large dispersions for both single nodesand edges, resembling empirical data.
II. SIMULTANEOUSLY LARGE VARIABILITYOF INTER-EVENT TIMES ON NODES ANDEDGES IS COMMON AND NON-TRIVIAL
Before presenting and analyzing our model, in this sec-tion we provide empirical evidence that heavy-tailed dis-tributions of IETs are simultaneously present for edgesand nodes in single data sets. We use data collected bythe SocioPatterns collaboration [50]. The survival func-tion of the IETs (i.e., probability that the IET, τ , is largerthan the specified value) for individual edges and nodesfor social contact data in a primary school [50], whichwe refer to as PrimarySchool, is shown in Fig. 2(a), and2(b), respectively. (See Appendix A for the results forthe other data sets.) We only used edges with at least100 events to calculate the distribution and the follow-ing statistics. The relatively slow decay in Fig. 2 sug-gests heavy-tailed distributions for both edges and nodesacross some scales of τ .We quantified the dispersion of IET distributions bythe coefficient of variation (CV). The CV of a distributionis defined as the standard deviation divided by the mean.For an IET distribution, one can writeCV = (cid:115) (cid:104) τ (cid:105)(cid:104) τ (cid:105) − , (1)where (cid:104)·(cid:105) represents the average over edges or nodes. APoisson process produces an exponential IET distribu-tion, which yields CV = 1. A periodic process yieldsCV = 0. A heavy-tailed distribution yields a large valueof CV. Table I shows the mean and standard deviationof the CV of the IET distribution for edges and nodes foreach data set. Table I indicates that all data sets yieldCV values considerably larger than 1 for both edges andnodes. Therefore, the simultaneous presence of hetero-geneous distributions (i.e., with a larger dispersion than FIG. 2. Survival function, P ( τ ), of IETs on (a) edges and (b) nodes for the PrimarySchool data set (black) and simulationof our model (red). In this and the following analyses of the empirical data, we treated the data with two steps. First, weaggregated consecutive events with 20 s duration; the temporal resolution of the original data set is 20 s, i.e., the social contactsare measured every 20 s. We perform this first step because we are analyzing the IETs but not the duration of events. Forinstance, we aggregated a sequence of event times {
20, 40, 60, 100, 200, 220 } (in s) between two nodes into a sequence withthree contact events as {
20, 100, 200 } . In other words, the first event at t = 20 s lasts for 60 s, the second event at t = 100s lasts for 20 s, and the third event at t = 200 s lasts for 40 s. Second, to circumvent the effects of the circadian rhythm, weremoved IETs larger than eight hours. In both (a) and (b), we only considered edges that had at least 100 events. the Poisson case) of both edge’s and node’s IETs seemsto be common.We hypothesize that the correlation of the IET be-tween different edges sharing a node contributes to largevalues of the CV for individual nodes. Therefore, foreach data set, we uniformly randomly shuffled the IETson each edge within each day of recording, preservingthe time of the first and last events of each day on theedge. This shuffling method is equivalent to the time-line inter-event shuffling in an instant-event temporalnetwork (formally named P[ π L ( ∆ τ ) , t ]) described inRef. [56]. Because this shuffling procedure preserves thedistribution of IETs on each edge, the edge’s CV is un-changed. However, it affects the IETs and hence theCV for individual nodes. The node’s CV calculated fromthe shuffled data and the relative deviation defined by∆ = 100% × (shuffled − original) / original, are shown inTable I. For all data sets, the CV for the nodes consis-tently decreases when one shuffles the IETs. Therfore,the shuffling destroys the heavy-tailed nature of the IETsequences for individual nodes.To further support that the simultaneous presence ofheavy-tailed distributions of IETs on edges and nodesis nontrivial, we assess a common approach to generateevent sequences on each edge according to an indepen-dent renewal process with a power-law distribution, p ( τ ),of IETs. We call this model the stochastic temporal net-work model [3]. Consider the power-law distribution ofIETs given by p ( τ ) = α − τ ) α , (2) where α is a parameter. The first and second momentsof p ( τ ) are given by (cid:104) τ (cid:105) = 1 / ( α − α >
2, and (cid:104) τ (cid:105) = 2 / ( α − α − α >
3, respectively. There-fore, the CV of IETs is given byCV = (cid:112) ( α − α − , (3)where α >
3. We ran simulations with α = 3 and α = 3 . k = 2, k = 5, and k = 10 neighbors, andcalculated the CV for IETs on the node for each combina-tion of α and k . Equation (3) predicts an edge’s CV equalto 2.24 when α = 3 . α ≤
3, whichis consistent with the numerical results shown in Table II.By contrast, the CV for the node is considerably smallerthan that for edge and decreases towards 1 as k increases.Therefore, the stochastic temporal network model doesnot produce simultaneously heavy-tailed distributions ofIETs for edges and nodes. III. MODEL
We propose a model of node behavior that aims tosimultaneously produce large CVs for IETs on individ-ual edges and nodes. Consider a static network. Themodel is based on two main assumptions. First, we as-sume that each node stochastically switches between twostates, called the high-activity state (denoted by h ) andthe low-activity state (denoted by (cid:96) ). The plausibilityof this assumption is supported by various empirical andmodeling studies [35–39, 41, 57–59]. Second, we assumethat two nodes adjacent by an edge in the static network TABLE I. CV of IETs on edges and nodes for SocioPatterns data sets. “Node CV, shuffled” corresponds to the node’sCV calculated after the timeline inter-event shuffling. The data originate from a primary school (PrimarySchool) [51, 52], ascientific conference (SFHH) [50], a workplace (Office15) [50], a hospital (Hospital) [53], and a high school in two different years(HighSchool12 [54] and HighSchool13 [55]). The CV values shown are the average ± standard deviation. In this table and thefollowing figures and tables using the same data sets, we used the edges that have at least 100 events; we used nodes that have,among its incident edges, at least one edge with at least 100 events and all the other edges with at least 10 events.PrimarySchool SFHH Office15 Hospital HighSchool12 HighSchool13Edge CV, original 2 . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± . . ± .
5∆ (i.e., relative deviation) − − − − − − events. α k Edge CV Node CV2 3 . ± . . ± .
03 5 3 . ± . . ± .
010 3 . ± . . ± .
02 2 . ± . . ± . . ± . . ± .
010 2 . ± . . ± . have a contact event much more likely when both nodesare in the high-activity state than otherwise. The intu-ition behind this assumption is that a pairwise humancontact event may be much more likely to occur whenboth individuals are motivated to interact than other-wise.Each node switches between h and (cid:96) according toa two-state continuous-time Markov process. In otherwords, the node switches to the opposite state accordingto a Poisson process whose rate depends on the currentstate. We denote by r h → (cid:96) the rate at which a node instate h changes to state (cid:96) , and similar for r (cid:96) → h . Weassume that different nodes share the same r h → (cid:96) and r (cid:96) → h values but are associated with independent two-state Markovian processes.If both adjacent nodes are in the h state, then eventson the edge are assumed to occur according to a Pois-son process at a higher rate denoted by λ h . Otherwise,the edge produces events according to a Poisson processat a lower rate denoted by λ (cid:96) ( < λ h ). This process isschematically shown in Fig. 3.The master equation for the probability that a node isin state h , denoted by p h ( t ), where t represents time, isgiven by d p h ( t )d t = r (cid:96) → h [1 − p h ( t )] − r h → (cid:96) p h ( t ) . (4) FIG. 3. Schematic illustration of the model. The eventsbetween two nodes occur at a higher rate λ h if and only ifboth are in the high-activity state (time windows shown inred) Otherwise, events occur at a lower rate λ (cid:96) (shown ingray). Therefore, the stationary probability of finding a node instate h is given by p ∗ h = r (cid:96) → h r h → (cid:96) + r (cid:96) → h . (5)The stationary probability of finding a node in state (cid:96) is p ∗ (cid:96) = 1 − p ∗ h . IV. RESULTSA. The model produces simultaneously largevariability of IETs on edges and nodes
We numerically simulated the model to generate se-quences of IETs and computed the survival function andthe CV of IETs for individual edges and nodes. Apartfrom the structure of the static network, our model hasfour parameters, r h → (cid:96) , r (cid:96) → h , λ h , and λ (cid:96) . Equation (5)yields r (cid:96) → h = r h → (cid:96) p ∗ h / (1 − p ∗ h ). We write λ (cid:96) in terms of λ h as λ (cid:96) = γλ h , (6)where 0 < γ <
1. Simultaneously changing( r h → (cid:96) , r (cid:96) → h , λ h , λ (cid:96) ) to ( cr h → (cid:96) , cr (cid:96) → h , cλ h , cλ (cid:96) ), where c >
0, is equivalent to not changing these four parametersand changing the time from t to ct . Therefore, withoutloss of generality, we set λ h = 1, unless we state oth-erwise. In the following simulations, for fixed values of r h → (cid:96) , we varied γ and p ∗ h . Using the Gillespie algorithm[60], we generated events on the edges until all edges hadat least 10 events.We used a star network composed of a node v and its k neighbors. Initially, each of the k + 1 nodes is inde-pendently in state h or (cid:96) with probability p ∗ h or (1 − p ∗ h ),respectively. Then, we run a continuous-time Markovprocess independently for each node. At each time, de-pending on the state of each edge, we generate events onthe edge at rate λ h or λ (cid:96) . It should be noted that a nextevent on an edge may not be simply produced as a sin-gle Poisson process. For example, suppose that the twonodes connected by an edge are both in state h . Then,the time to the next event is drawn from the exponentialdistribution p ( τ ) = λ h e − λ h τ . However, if either nodeswitches to state (cid:96) before the next event occurs, thenone has to discard the scheduled time to the next event,which was generated from p ( τ ) = λ h e − λ h τ , and redrawthe time to the next event from p ( τ ) = λ (cid:96) e − λ (cid:96) τ .We consider a star network composed of a node v andits two neighbors as an example. We set r h → (cid:96) = 2 × − , r (cid:96) → h = 4 . × − , λ h = 6 × − , and λ (cid:96) = 3 . × − ,which yields p ∗ h ≈ . γ ≈ .
06. The survival func-tion of IETs produced by the model is shown by thered lines for the two edges and node v in Fig. 2(a) andFig. 2(b), respectively. The survival function for both thetwo edges and v decays more slowly than exponentially,roughly consistent with the non-Poissonian behavior ob-served in the empirical data (black lines in Fig. 2). Forour model, the CV for the two edges is equal to 2.8 and3.0, and that for v is equal to 2.6. These values of CVare statistically within the ranges of the CV for the Pri-marySchool data set (see Table I).To examine different parameter values of the model,we varied γ and p ∗ h for each of three values of r h → (cid:96) (i.e., r h → (cid:96) = 10 − , r h → (cid:96) = 10 − , and r h → (cid:96) = 10 − ) and threevalues of k (i.e., k = 2, k = 5, k = 10). We show theCV values in Fig. 4. The figure indicates that the modelproduces large CV values simultaneously for edges andnodes in a broad parameter region. In particular, theCV is large when γ (= λ (cid:96) /λ h ) is small and p ∗ h is near 0.7,across the range of r h → (cid:96) and k . B. Analytical evaluation of the CV of inter-eventtimes
In this section, we provide an analytical account forthe CV values observed in section IV A. The state ofthe edge is specified by the states of two nodes form-ing the edge. We denote by h when both nodes arein state h , and likewise for h(cid:96) and (cid:96) . The edge stateobeys a three-state Markov process in the state space S = { h , h(cid:96), (cid:96) } , where we do not distinguish between h(cid:96) and (cid:96)h , and represent both by h(cid:96) . In our model, theIET distribution conditioned on the edge’s state is givenby g ( τ | h ) = λ h e − λ h τ and g ( τ | (cid:96) ) = g ( τ | h(cid:96) ) = λ (cid:96) e − λ (cid:96) τ ,where g ( τ |· ) represents the distribution of IETs condi-tioned on the edge’s state.The Markov process of node activity is independent fordifferent nodes. Therefore, two nodes are simultaneouslyin state h in the equilibrium with probability p ∗ h . Themean number of events produced when both nodes arein state h is given by λ h p ∗ h T , where T is the observationtime. The mean number of events produced when eithernode is in state (cid:96) is given by λ (cid:96) (1 − p ∗ h ) T . Therefore,an IET is produced at rate λ h and λ (cid:96) with probability λ h p ∗ h / Ω and λ (cid:96) (1 − p ∗ h ) / Ω , respectively, where Ω = λ h p ∗ h + λ (cid:96) (1 − p ∗ h ). By combining these contributionsand ignoring IETs during which the edge’s state changes,we obtain the probability density function (PDF) of IETsfor an edge as a mixture of two exponential distributionsas f edge ( τ ) = λ h p ∗ h Ω λ h e − λ h τ + λ (cid:96) (1 − p ∗ h )Ω λ (cid:96) e − λ (cid:96) τ . (7)The first two moments of this PDF are given by (cid:104) τ (cid:105) edge ≡ (cid:90) ∞ τ f edge ( τ ) dτ = 1Ω . (8)and (cid:104) τ (cid:105) edge ≡ (cid:90) ∞ τ f edge ( τ ) dτ = 2Ω (cid:20) p ∗ h λ h + (1 − p ∗ h ) λ (cid:96) (cid:21) . (9)By substituting Eqs. (8) and (9) into Eq. (1), we obtainCV edge = (cid:115) p ∗ h (1 − p ∗ h )(1 − γ ) γ , (10)where CV edge is the CV for the edge’s IETs.Now we consider a node v with k neighbors. From theviewpoint of node v , the sequence of events is a super-position of the events over its k edges. Because the k nodes are statistically the same for v , it is sufficient toconsider a 2( k + 1)-state Markov process with state space S k = { h, (cid:96) } × { h k , h k − (cid:96), . . . , h(cid:96) k − , (cid:96) k } , where the firstset in the product of the two sets represents the state of v ,and the second set represents the states of v ’s neighbors.First, if v is in state h , which happens with proba-bility p ∗ h in the equilibrium, the IET distribution for v depends on the states of v ’s neighbors. The probabil-ity of finding exactly k h neighbors in state h is given by (cid:0) kk h (cid:1) p ∗ k h h (1 − p ∗ h ) k − k h , where (cid:0) kk h (cid:1) is the binomial coeffi-cient. In this situation, v experiences events producedby a Poisson process at rate k h λ h + ( k − k h ) λ (cid:96) becausethe superposition of the k independent Poisson processeson edges with rate λ h or λ (cid:96) is a Poisson process withthe summed rate. Second, if v is in state (cid:96) , which hap-pens with probability 1 − p ∗ h , all edges produce events atrate λ (cid:96) . In this situation, v experiences events producedby a Poisson process at rate kλ (cid:96) . Therefore, an event FIG. 4. CV values for IETs generated by our original model. We set r h → (cid:96) = 10 − in panels (a), (b), (c), and (d), r h → (cid:96) = 10 − in panels (e), (f), (g), and (h), and r h → (cid:96) = 10 − in panels (i), (j), (k), and ( (cid:96) ). We simulated an edge (panels (a), (e), and (i)),a node with k = 2 neighbors (panels (b), (f), and (j)), a node with k = 5 neighbors (panels (c), (g), and (k)), and a node with k = 10 neighbors (panels (d), (h), and ( (cid:96) )). For each set of parameter values, we generated IETs on the edges until all edgeshad at least 10 events. occurs when v is in state h and it has k h neighbors instate h with probability p ∗ h (cid:0) kk h (cid:1) p ∗ k h h (1 − p ∗ h ) k − k h [ k h λ h +( k − k h ) λ (cid:96) ] e − [ k h λ h +( k − k h ) λ (cid:96) ] τ / Ω k and when v is in state (cid:96) with probability (1 − p ∗ h )( kλ (cid:96) ) e − kλ (cid:96) τ / Ω k , where Ω k = kλ h p ∗ h + kλ (cid:96) (1 − p ∗ h ).By combining these contributions, we derive the PDFfor IETs on a node with k neighbors as f k ( τ ) = p ∗ h Ω k k (cid:88) k h =0 (cid:18) kk h (cid:19) p ∗ k h h (1 − p ∗ h ) k − k h [ k h λ h + ( k − k h ) λ (cid:96) ] e − [ k h λ h +( k − k h ) λ (cid:96) ] τ + 1 − p ∗ h Ω k ( kλ (cid:96) ) e − kλ (cid:96) τ = e − kλ (cid:96) τ Ω k (cid:40) p ∗ h (cid:104) − p ∗ h (1 − e − ( λ h − λ (cid:96) ) τ ) (cid:105) k (cid:34) k ( λ h − λ (cid:96) ) p ∗ h e − ( λ h − λ (cid:96) ) τ ( λ h − λ (cid:96) + 2)1 − p ∗ h (1 − e − ( λ h − λ (cid:96) ) τ ) ++ k ( k − λ h − λ (cid:96) ) p ∗ h e − λ h − λ (cid:96) ) τ (cid:2) − p ∗ h (1 − e − ( λ h − λ (cid:96) ) τ ) (cid:3) + ( kλ (cid:96) ) (cid:35) + (1 − p ∗ h )( kλ (cid:96) ) (cid:41) . (11)The first two moments of the PDF are given by (cid:104) τ (cid:105) k ≡ (cid:90) ∞ τ f k ( τ ) dτ = 1Ω k (12) and (cid:104) τ (cid:105) k ≡ (cid:90) ∞ τ f k ( τ ) dτ == 2 p ∗ h Ω k k (cid:88) k h =0 (cid:18) kk h (cid:19) p ∗ k h h (1 − p ∗ h ) k − k h k h λ h + ( k − k h ) λ (cid:96) + 2(1 − p ∗ h ) kλ (cid:96) Ω k . (13)By substituting Eqs. (12) and (13) into Eq. (1), we obtain the CV for node v asCV k = (cid:118)(cid:117)(cid:117)(cid:116) k [ p ∗ h + (1 − p ∗ h ) γ ] (cid:34) p ∗ h k (cid:88) k h =0 (cid:18) kk h (cid:19) p ∗ k h h (1 − p ∗ h ) k − k h k h + ( k − k h ) γ + 1 − p ∗ h kγ (cid:35) − . (14)Equation (14) reduces to Eq. (10) when k = 1.For any k , Eq. (14) yields lim γ → CV k → ∞ andlim γ → CV k = 1. In addition, the derivative of Eq.(14) withrespect to γ is negative for any 0 < γ <
1. Therefore,CV k monotonically decreases towards 1 as γ →
1, whichis consistent with Fig. 4. Next, equating the derivativeof the right-hand side of Eq. (10) with respect to p ∗ h tozero yields p ∗ h = 1 / √ ≈ .
71 for any γ . The value of p ∗ h at which the derivative of the right-hand side of Eq. (14)with respect to p ∗ h is equal to zero depends on k and γ ,but we numerically obtain p ∗ h ≈ . k and γ .Therefore, the CV k is large when γ small and p ∗ h ≈ . r h → (cid:96) = 10 − . In this case, the relativeerror is small across the entirety of our parameter re-gion. When r h → (cid:96) = 10 − (see Figs. 5(e)–5(h)), and r h → (cid:96) = 10 − (see Figs. 5(i)–5( (cid:96) )), the relative error islarge when γ is small and p ∗ h is large.The relative error increases as r h → (cid:96) increases for thefollowing reason. A large value of r h → (cid:96) and moderatevalue of p ∗ h (i.e., p ∗ h that is not too close to 0 or 1) im-plies large r (cid:96) → h , because r (cid:96) → h = r h → (cid:96) p ∗ h / (1 − p ∗ h ). When r (cid:96) → h is large, the mean IET for edges that are producedat event rate λ (cid:96) , which is equal to 1 /λ (cid:96) , is longer thanthe typical duration of the low-activity state of the edge(i.e., either h(cid:96) or (cid:96) ), which is proportional to 1 /r (cid:96) → h .Note that, the low-activity state of an edge finishes onlywhen both of the two nodes have transited from state (cid:96) to h , which occurs at rate r (cid:96) → h for each node. Therefore,an edge is populated with sufficiently many IETs pro-duced in the low-activity state of the edge if and only if1 /r (cid:96) → h (cid:29) /λ (cid:96) , which leads to r (cid:96) → h (cid:28) λ (cid:96) . According toour parametrization, we obtain r (cid:96) → h = r h → (cid:96) p ∗ h / (1 − p ∗ h )and λ (cid:96) = γλ h = γ , because we set λ h = 1. Therefore,an edge is populated with sufficiently many IETs pro-duced in its low-activity state if r h → (cid:96) (cid:28) (1 − p ∗ h ) γ/p ∗ h .For example, when γ = 0 .
01 and p ∗ h = 0 .
9, we need r h → (cid:96) (cid:28) (1 − p ∗ h ) γ/p ∗ h ≈ − . This condition is vio-lated when r h → (cid:96) = 10 − or 10 − . For these r h → (cid:96) val-ues, the probability that events occur in the low-activitystate of the edge is small. However, our analytical deriva-tion of the CV assumes that sufficiently many events andhence IETs occur in the typical duration of both low-activity and high-activity states of the edge at respective rates, i.e., λ (cid:96) and λ h . This explains the discrepancy be-tween the theoretical and numerical results observed inFigs. 5(e)–( (cid:96) ).However, the model is still capable of producing largeCV values even if few IETs are produced in the low-activity state of the edge. In this situation, the IETbetween the last event of a high-activity period of theedge and the first event of the next high-activity periodwould generate a long IET, contributing to a relativelyheterogeneous distribution of IETs. This regime is notpredicted by our analytical solution.If r h → (cid:96) ≈ λ h , then few IETs are produced during the h state of the edge on average. Therefore, a large CVvalue requires r h → (cid:96) (cid:28) λ h , which is satisfied in Figs. 4and 5 because the largest value of r h → (cid:96) that we use is10 − and we have set λ h = 1. C. Correlation between consecutive inter-eventtimes
In human activities, IETs for both edges and nodesare often positively correlated, i.e., long IETs tend to befollowed by long IETs and vice versa [1, 4, 41, 61]. To ex-amine this property, we compute the memory coefficient, M , of a sequence of IETs [61], defined as M ≡ n − n − (cid:88) i =1 ( τ i − m )( τ i +1 − m ) σ σ , (15)where n is the number of IETs in the sequence, m and σ are the average and standard deviation of { τ , τ , . . . , τ n − } , respectively, and m and σ are the av-erage and standard deviation of { τ , τ , . . . , τ n } , respec-tively. The memory coefficient measures the correlationcoefficient of consecutive IETs, ( τ i , τ i +1 ).Figure 6 shows the memory coefficient for the empiri-cal data. In all data sets, the edges show predominantlypositive memory coefficients with values lying mostly be-tween 0 and 0.1 (see Fig. 6(a)). The memory coefficientfor nodes is also predominantly positive and tends to belarger than that for the edges (see Fig. 6(b)). These val-ues are in accordance with previous results for varioushuman activities [61].The memory coefficient for sequences of IETs gener-ated by our model is shown in Fig. 7. The figure in-dicates that the model produces positive M for bothedges and nodes, which is qualitatively consistent with FIG. 5. Relative error between the analytically and numerically evaluated CV. We set r h → (cid:96) = 10 − in panels (a), (b), (c),and (d), r h → (cid:96) = 10 − in panels (e), (f), (g), and (h), and r h → (cid:96) = 10 − in panels (i), (j), (k), and ( (cid:96) ). The results are for anedge (panels (a), (e), and (i)), a node with k = 2 neighbors (panels (b), (f), and (j)), a node with k = 5 neighbors (panels (c),(g), and (k)), and a node with k = 10 neighbors (panels (d), (h), and ( (cid:96) )). the empirical data. However, the memory coefficient val-ues produced by the model are substantially larger thanthe empirical values. D. Variants of the model
In our original model, events occur on the edge at ahigher rate if and only if both nodes forming the edgeare in state h . In this section, we study two variants ofthe model. In the first variant, we assume that eventon an edge occur at the higher rate, λ h , if either nodeconnected to the edge, not necessarily both nodes, is instate h . Events on the edge occur at the lower rate λ (cid:96) ifand only if both nodes are in state (cid:96) . We call this variantthe OR model. An interpretation of the OR model isthat, if an individual wants to interact with a neighbor,he/she can do so at the higher event rate regardless ofwhether or not the neighbor wants to interact.For the OR model, we calculated the CV for IETs onedges and nodes by scanning the same values of p ∗ h , γ , r h → (cid:96) , and k as those used in Fig. 4. The results areshown in Fig. 8. As in the original model, the OR modelproduces large CV values for both edges and nodes when γ is small. With respect to p ∗ h , the OR model produceslarge CV values when p ∗ h ≈ . p ∗ h value that maximizes the CV decreases as k increases. This behavior is consistent with the analytical prediction(Appendix B). Unlike the original model, the region thatthe OR model produces large CV values for the nodeshrinks as k increases.In the second variant of the model, we assume thatthe node’s h and (cid:96) states independently contribute λ h and λ (cid:96) , respectively, to the event rate of the edge. Inother words, events occur on the edge at rate 2 λ h if bothnodes are in the h state, 2 λ (cid:96) if both nodes are in the (cid:96) state, and λ h + λ (cid:96) if one node is in the h state andthe other node is in the (cid:96) state. An interpretation ofthis variant of the model, which we call the IND model(named after “independent”), is that the state of eachindividual independently contributes to the frequency ofevents between two individuals.The CV values for the IND model are shown in Fig. 9.Similarly to the original and the OR models, the INDmodel produces large CV values when γ is small. For anygiven γ , the IND model produces large values of CV when0 . (cid:47) p ∗ h (cid:47) . p ∗ h valuethat maximizes the CV decreases as k increases. Thisbehavior is consistent with the analytical result shown inAppendix C. Therefore, the IND model behaves similarlyto the OR model, i.e., it produces large CV values when γ is small and p ∗ h ≈ .
3, and the parameter region inwhich the node’s CV is large shrinks as k increases.To quantitatively compare the three models, we calcu- FIG. 6. Box plots of the memory coefficient for (a) edgesand (b) nodes for the empirical data sets. The box shows themedian (orange line), the first quartile ( Q ), and the thirdquartile ( Q ); the whiskers show the minimum ( Q − . × IQR) and the maximum ( Q + 1 . × IQR) values excludingoutliers, where IQR = Q − Q . The open circles are theoutliers. The green triangles are the sample means. lated two quantities. First, we compute the largest CVvalue produced by each model when we vary γ and p ∗ h in the parameter region used in Figs. 4, 8, and 9. Thelargest CV value is compared among the three modelsin Fig. 10(a), where we set r h → (cid:96) = 10 − and vary thenode’s degree k . The figure indicates that the originalmodel consistently produces the largest CV values as k increases, although the OR model produces comparablylarge CV values up to k = 3. Unlike the original model,the largest CV value produced by the OR and IND mod-els visibly decreases as k increases. Second, we computethe fraction of the ( γ , p ∗ h ) pairs for which the CV valueis larger than two; a CV value larger than two is consis-tent with the results for empirical data (see Section II). Alarge fraction value implies that a CV value larger thantwo is robustly produced for various parameter combina-tions. The result for this analysis is shown in Fig. 10(b),where we again set r h → (cid:96) = 10 − and vary k . The figureindicates that the original model has a larger fractionof the parameter region with CV larger than two than the OR and IND models. The fraction remains roughlyconstant as k increases for the original model, whereas itrapidly decreases as k increases for the OR and the INDmodels. Figures 10(a) and 10(b) altogether suggest thatthe original model is more capable of producing largeCVs of IETs on both edges and nodes than the OR andIND models, particularly when the node has a large de-gree. The results are qualitatively the same for larger r h → (cid:96) values, i.e., r h → (cid:96) = 10 − (Figs. 10(c) and (d)) and r h → (cid:96) = 10 − (Figs. 10(e) and (f)). V. DISCUSSION
We have started from the observation that heavy-taileddistributions of IETs are simultaneously present for bothindividual nodes and edges in the same empirical data.We have proposed a continuous-time model and its vari-ants for generating discrete events on edges that replicatethis behavior to different extents. Our main model cru-cially assumes that each node alternates between high-activity and low-activity states in a Markovian manner.We showed that the original model, which requires thatboth nodes are in the high-activity state for the edge tohave frequent events, is capable of producing large CVvalues for both individual nodes and edges in a broadparameter region. The other two variants of the modelare also capable of producing reasonably large CV valuesto some extent. The proposed models allow interpreta-tions. For example, in the original model, two nodes arelikely to interact if and only if both of them feel like in-teracting with others.We have derived analytical solutions for our models bydiscarding the effects of IETs that contain state transi-tions of the edge. The analytical solution was accuratewhen two conditions were met. First, the ratio of thehigh- to low-activity event rates (i.e., 1 /γ ) should not beextremely large, such as 100. This condition is probablynot unrealistic. Second, the state transition rate of thenode should be sufficiently small compared to the eventrates on edges. Under this condition, an epoch of thehigh- or low-activity state of an edge is typically longenough to host sufficiently many events at the constantrate, which is either λ h or λ (cid:96) . This is the situation thatthe theory in Section IV B assumes. Otherwise, a largefraction of IETs contains transitions of the edge’s state,which cause a systematic discrepancy of the analyticalexpression from the numerical results.In our models, each node switches between two states.A possible extension of this assumption is to the case ofmore than two states for each node. Then, dependingon how such a model translates the nodes’ states intothe event rate, the distribution of IETs on edges may beapproximately a mixture of more than two exponentialdistributions, which may resemble or actually produceheavy-tailed distributions [36, 57, 62, 63]. In fact, a mix-ture of a small number of exponential distributions, in-cluding the case of just two exponential distributions, is0 FIG. 7. Memory coefficient for IETs generated by our original model. We set r h → (cid:96) = 10 − in panels (a), (b), (c), and (d), r h → (cid:96) = 10 − in panels (e), (f), (g), and (h), and r h → (cid:96) = 10 − in panels (i), (j), (k), and ( (cid:96) ). The results are for an edge (panels(a), (e), and (i)), a node with k = 2 neighbors (panels (b), (f), and (j)), a node with k = 5 neighbors (panels (c), (g), and (k)),and a node with k = 10 neighbors (panels (d), (h), and ( (cid:96) )). For each set of parameter values, we generated IETs on the edgesuntil all edges had at least 10 events. often sufficient for approximating many empirical heavy-tailed distributions of IETs [36, 37, 57, 62]. Therefore,one should carefully assess trade-offs between the com-plexity of extended models and the explanatory powerof the model that one gains by assuming more states fornodes.The distribution of IETs affects how disease and infor-mation spread across contact networks [8–15]. Becausemany time-stamped event data probably have heavy-tailed distributions of IETs for both individual nodesand edges, the temporal network models proposed in thepresent study are expected to be useful for modeling dy-namical processes on temporal networks including con-tagion processes. It seems that model-based studies ofcontagion processes on temporal networks have not paidmuch attention to the simultaneous presence of heavy-tailed distributions of IETs on nodes and edges [4]. Howthis property affects key indicators of contagion processessuch as the epidemic threshold, the final epidemic size,and equilibrium fraction of infected nodes, as well as in-dicators of other dynamical processes, warrants futurework. ACKNOWLEDGMENTS
Appendix A: Distributions of inter-event times forthe other data sets
Figure 11 shows the survival function of IETs for thedifferent data sets.
Appendix B: Analytical evaluation of the CV ofinter-event times for the OR model
In this section, we analytically examine the CV of IETsfor the OR model. In the OR model, events occur at the1
FIG. 8. CV values for IETs generated by the OR model. We set r h → (cid:96) = 10 − in panels (a), (b), (c), and (d), r h → (cid:96) = 10 − inpanels (e), (f), (g), and (h), and r h → (cid:96) = 10 − in panels (i), (j), (k), and ( (cid:96) ). The results are for an edge (panels (a), (e), and(i)), a node with k = 2 neighbors (panels (b), (f), and (j)), a node with k = 5 neighbors (panels (c), (g), and (k)), and a nodewith k = 10 neighbors (panels (d), (h), and ( (cid:96) )). For each set of parameter values, we generated IETs on the edges until alledges had at least 10 events. lower rate λ (cid:96) if and only if the two nodes forming anedge are in state (cid:96) . Therefore, using the same reasoningas that for the original model in Sec. IV B, one obtainsthe PDF of IETs for an edge as follows: f OR edge ( τ ) = λ (cid:96) p ∗ (cid:96) Ω OR λ (cid:96) e − λ (cid:96) τ + λ h (1 − p ∗ (cid:96) )Ω OR λ h e − λ h τ , (B1)where Ω OR = λ (cid:96) p ∗ (cid:96) + λ h (1 − p ∗ (cid:96) ) and, for convenience, wehave used p ∗ (cid:96) = 1 − p ∗ h .The first two moments of this PDF are given by (cid:104) τ (cid:105) edge ≡ (cid:90) ∞ τ f OR edge ( τ ) dτ = 1Ω OR (B2) and (cid:104) τ (cid:105) edge ≡ (cid:90) ∞ τ f OR edge ( τ ) dτ = 2Ω OR (cid:20) p ∗ (cid:96) λ (cid:96) + (1 − p ∗ (cid:96) ) λ h (cid:21) . (B3)By substituting Eqs. (B2) and (B3) into Eq. (1), weobtain CV OR edge = (cid:115) p ∗ (cid:96) (1 − p ∗ (cid:96) )(1 − γ ) γ , (B4)where CV OR edge is the CV for the edge’s IETs for the ORmodel.We proceed with the same steps as those in Sec. IV B tocalculate the CV for a node with k neighbors as follows.The PDF for a node with k neighbors is given by f OR k ( τ ) = p ∗ (cid:96) Ω OR k k (cid:88) k (cid:96) (cid:18) kk (cid:96) (cid:19) p ∗ k (cid:96) (cid:96) (1 − p ∗ (cid:96) ) k − k (cid:96) [ λ (cid:96) k (cid:96) + λ h ( k − k (cid:96) )] e − [ λ (cid:96) k (cid:96) + λ h ( k − k (cid:96) )] τ + 1 − p ∗ (cid:96) Ω OR k ( kλ h ) e − kλ h τ , (B5)where k (cid:96) is the number of neighbors in the (cid:96) state, andΩ OR k = kλ (cid:96) p ∗ (cid:96) + kλ h (1 − p ∗ (cid:96) ). The first two moments of f OR k ( τ ) are given by (cid:104) τ (cid:105) k ≡ (cid:90) ∞ τ f k ( τ ) dτ = 1Ω OR k (B6)2 FIG. 9. CV values for IETs generated by the IND model. We set r h → (cid:96) = 10 − in panels (a), (b), (c), and (d), r h → (cid:96) = 10 − in panels (e), (f), (g), and (h), and r h → (cid:96) = 10 − in panels (i), (j), (k), and ( (cid:96) ). The results are for an edge (panels (a), (e), and(i)), a node with k = 2 neighbors (panels (b), (f), and (j)), a node with k = 5 neighbors (panels (c), (g), and (k)), and a nodewith k = 10 neighbors (panels (d), (h), and ( (cid:96) )). For each set of parameter values, we calculated the CV as an average over100 realizations of the simulation, each simulation with 10 events. and (cid:104) τ (cid:105) k ≡ (cid:90) ∞ τ f k ( τ ) dτ == 2 p ∗ (cid:96) Ω OR k k (cid:88) k (cid:96) =0 (cid:18) kk (cid:96) (cid:19) p ∗ k (cid:96) (cid:96) (1 − p ∗ (cid:96) ) k − k (cid:96) k (cid:96) λ (cid:96) + ( k − k (cid:96) ) λ h + 2(1 − p ∗ (cid:96) ) kλ h Ω OR k . (B7)By substituting Eqs. (B6) and (B7) into Eq. (1), we ob-tain the CV for node v asCV OR k = (cid:118)(cid:117)(cid:117)(cid:116) k [ p ∗ (cid:96) + (1 − p ∗ (cid:96) ) γ ] (cid:34) p ∗ (cid:96) k (cid:88) k (cid:96) =0 (cid:18) kk (cid:96) (cid:19) p ∗ k (cid:96) (cid:96) (1 − p ∗ (cid:96) ) k − k (cid:96) k (cid:96) + ( k − k (cid:96) ) γ + 1 − p ∗ (cid:96) kγ (cid:35) − , (B8)which generalizes Eq. (B4).The behavior of Eq. (B8) with respect to γ isqualitatively the same as that of the original model(i.e., Eq. (14)). In other words, lim γ → CV OR k → ∞ ,lim γ → CV OR k = 1, and CV OR k monotonically decreasesas γ increases, which is consistent with Fig. 8. Theextremum of Eq. (B4) with respect to p ∗ (cid:96) occurs at p ∗ (cid:96) = 1 / √ ≈ .
71, hence, p ∗ h ≈ .
29, for any γ . There-fore, the CV for the edge is large when γ is small and p ∗ h ≈ .
3. By contrast, the value of p ∗ h that maximizesEq. (B8) for a given γ value strongly depends on k . InFig. 12(a), we numerically inspect this dependence. Thevalue of p ∗ h that maximizes Eq. (B8) decreases as k in-creases, and, for a given value of k , it monotonically in-3 FIG. 10. Comparison of CV values obtained from the origi-nal, OR, and IND models. Panels (a), (c), and (e) show thelargest CV value in the ( γ , p ∗ h ) parameter region explored inFigs. 4–9. Panels (b), (d), and (f) show the fraction of the ( γ , p ∗ h ) pairs for which the CV values is larger than 2. Because k represents the degree of the node, k = 1 corresponds to thecase of the single edge. We set r h → (cid:96) = 10 − in (a) and (b), r h → (cid:96) = 10 − in (c) and (d), and r h → (cid:96) = 10 − in (e) and (f). creases as γ increases. These results are consistent withFig. 8. Appendix C: Analytical evaluation of the CV ofinter-event times for the IND model
In this section, we analytically examine the CV of IETsfor the IND model. In the IND model, events on an edgeoccur at rate 2 λ h if both nodes are in state h , at rate λ h + λ (cid:96) if one node is the h and the other node is in the (cid:96) state, and at rate 2 λ (cid:96) if both nodes are in the (cid:96) state.Therefore, the PDF of IETs for an edge is given by f IND edge ( τ ) = 2 λ h p ∗ h Ω IND λ h e − λ h τ ++ 2( λ h + λ (cid:96) ) p ∗ h (1 − p ∗ h )Ω IND ( λ h + λ (cid:96) ) e − ( λ h + λ (cid:96) ) τ ++ 2 λ (cid:96) (1 − p ∗ h ) Ω IND λ (cid:96) e − λ (cid:96) τ , (C1)where Ω IND = 2 λ h p ∗ h + 2( λ h + λ (cid:96) ) p ∗ h (1 − p ∗ h ) + 2 λ (cid:96) p ∗ h . FIG. 11. Survival function, P ( τ ), of IETs on edges andnodes for the different data sets. The first two moments of f IND edge ( τ ) are give by (cid:104) τ (cid:105) edge ≡ (cid:90) ∞ τ f IND edge dτ = 1Ω IND (C2)and (cid:104) τ (cid:105) edge ≡ (cid:90) ∞ τ f ORedge ( τ ) dτ = 2Ω IND (cid:20) p ∗ h λ h + 2 p ∗ h (1 − p ∗ h ) λ h + λ (cid:96) + (1 − p ∗ h ) λ (cid:96) (cid:21) . (C3)4 FIG. 12. Values of p ∗ h that yield the maximum of (a) CV OR k and (b) CV IND k for different values of γ and k . By substituting Eqs. (C2) and (C3) into Eq. (1),we obtainCV
IND edge = (cid:115) p ∗ h [1 − (2 − γ ) p ∗ h + (1 − γ ) p ∗ h ] γ (1 + γ ) , (C4)where CV IND edge is the CV for the edge’s IETs for the INDmodel.The PDF for IETs on a node with k neighbors is givenby f IND k ( τ ) = k (cid:88) k h =0 (cid:18) kk h (cid:19) p ∗ k h h (1 − p ∗ h ) k − k h Ω IND k (cid:26) p ∗ h [ λ h ( k + k h ) + λ (cid:96) ( k − k h )] e − [ λ h ( k + k h )+ λ (cid:96) ( k − k h )] τ ++ (1 − p ∗ h ) [ λ h k h + λ (cid:96) (2 k − k h )] e − [ λ h k h + λ (cid:96) (2 k − k h )] τ (cid:27) , (C5)where Ω IND k = 2 k [ λ h p ∗ h + λ (cid:96) (1 − p ∗ h )].The first two moments of f IND k ( τ ) are given by (cid:104) τ (cid:105) k ≡ (cid:90) ∞ τ f k ( τ ) dτ = 1Ω IND k (C6) and (cid:104) τ (cid:105) k ≡ (cid:90) ∞ τ f k ( τ ) dτ == k (cid:88) k h =0 (cid:18) kk h (cid:19) p ∗ k h h (1 − p ∗ h ) k − k h Ω IND k (cid:40) p ∗ h λ h ( k + k h ) + λ (cid:96) ( k − k h ) ++ 2(1 − p ∗ h ) λ h k h + λ (cid:96) (2 k − k h ) (cid:41) . (C7)By substituting Eqs. (C6) and (C7) into Eq. (1), weobtain the CV for node v as5CV IND k = (cid:118)(cid:117)(cid:117)(cid:116) k [ p ∗ h + γ (1 − p ∗ h )] k (cid:88) k h =0 (cid:18) kk h (cid:19) p ∗ k h h (1 − p ∗ h ) k − k h (cid:34) p ∗ h k + k h + γ ( k − k h ) + 1 − p ∗ h k h + γ (2 k − k h ) (cid:35) − , (C8)which generalizes Eq. (C4).The behavior of Eq. (C8) with respect to γ is qualita-tively the same as that of the original model. In otherwords, lim γ → CV IND k → ∞ , lim γ → CV IND k = 1, and the CV IND k monotonically decreases as γ increases, which are consis-tent with Fig. 9. The p ∗ h value that maximizes Eq. (C4)is given by p ∗ h = 2 − γ − (cid:112) − γ + γ − γ ) , (C9) which is plotted as the blue line in Fig. 12(b). The nu-merically obtained p ∗ h value that maximizes Eq. (C8) de-creases as k increases, and increases as γ increases (seeFig. 12(b)). These results are consistent with Fig. 9. [1] P. Holme and J. Saram¨aki, Phys Rep , 97 (2012).[2] P. Holme, Eur. Phys. J. B , 234 (2015).[3] N. Masuda and R. Lambiotte, A Guide to Temporal Net-works (World Scientific, Singapore, 2016).[4] M. Karsai, H.-H. Jo, K. Kaski, et al. , Bursty human dy-namics (Springer, Berlin, 2018).[5] P. Holme and J. Saram¨aki,
Temporal Network Theory (Springer, Cham, 2019).[6] A.-L. Barab´asi, Nature , 207 (2005).[7] A. V´azquez, J. G. Oliveira, Z. Dezs¨o, K.-I. Goh, I. Kon-dor, and A.-L. Barab´asi, Phys. Rev. E , 036127 (2006).[8] B. Min, K.-I. Goh, and A. V´azquez, Phys. Rev. E ,036102 (2011).[9] M. Karsai, M. Kivel¨a, R. K. Pan, K. Kaski, J. Kert´esz,A.-L. Barab´asi, and J. Saram¨aki, Phys. Rev. E ,025102(R) (2011).[10] L. E. C. Rocha, F. Liljeros, and P. Holme, PLoS Comput.Biol. , e1001109 (2011).[11] G. Miritello, E. Moro, and R. Lara, Phys. Rev. E ,045102(R) (2011).[12] N. Masuda and P. Holme, F1000prime Reports , 6(2013).[13] H.-H. Jo, J. I. Perotti, K. Kaski, and J. Kert´esz, Phys.Rev. X , 011041 (2014).[14] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, andA. Vespignani, Rev. Mod. Phys. , 925 (2015).[15] N. Masuda and P. Holme, Temporal Network Epidemiol-ogy (Springer, Berlin, 2017).[16] Y. Wu, C. Zhou, M. Chen, J. Xiao, and J. Kurths, Phys-ica A , 5832 (2010).[17] T. Takaguchi and N. Masuda, Phys. Rev. E , 036115(2011).[18] J. Fern´andez-Gracia, V. M. Egu´ıluz, and M. San Miguel,Phys. Rev. E , 015103(R) (2011).[19] R. Nishi and N. Masuda, EPL , 48003 (2014).[20] A. Li, L. Zhou, Q. Su, S. P. Cornelius, Y.-Y. Liu,L. Wang, and S. A. Levin, Nat. Commun. , 1 (2020).[21] F. Karimi and P. Holme, Physica A , 3476 (2013).[22] T. Takaguchi, N. Masuda, and P. Holme, PLoS ONE ,e68629 (2013).[23] V.-P. Backlund, J. Saram¨aki, and R. K. Pan, Phys. Rev.E , 062815 (2014). [24] S. Unicomb, G. I˜niguez, J. P. Gleeson, and M. Karsai,Preprint arXiv:2007.06223 (2020).[25] T. Hoffmann, M. A. Porter, and R. Lambiotte, Phys.Rev. E , 046102 (2012).[26] M. Starnini, A. Baronchelli, A. Barrat, and R. Pastor-Satorras, Phys. Rev. E , 056115 (2012).[27] L. Speidel, R. Lambiotte, K. Aihara, and N. Masuda,Phys. Rev. E , 012806 (2015).[28] N. Masuda, M. A. Porter, and R. Lambiotte, Phys. Rep. , 1 (2017).[29] A. V´azquez, Phys. Rev. Lett. , 248701 (2005).[30] G. Grinstein and R. Linsker, Phys. Rev. Lett. , 130201(2006).[31] G. Grinstein and R. Linsker, Phys. Rev. E , 012101(2008).[32] N. Masuda, J. S. Kim, and B. Kahng, Phys. Rev. E ,036106 (2009).[33] J. G. Oliveira and A. V´azquez, Physica A , 187(2009).[34] H.-H. Jo, R. K. Pan, and K. Kaski, Phys. Rev. E ,066101 (2012).[35] N. Masuda and P. Holme, Phys. Rev. Res. , 023163(2020).[36] M. Okada, K. Yamanishi, and N. Masuda, R. Soc. OpenSci. , 191643 (2020).[37] Z.-Q. Jiang, W.-J. Xie, M.-X. Li, W.-X. Zhou, andD. Sornette, J. Stat. Mech. , 073210 (2016).[38] R. D. Malmgren, D. B. Stouffer, A. E. Motter, andL. A. N. Amaral, Proc. Natl. Acad. Sci. U.S.A. ,18153 (2008).[39] R. D. Malmgren, D. B. Stouffer, A. S. L. O. Campanharo,and L. A. N. Amaral, Science , 1696 (2009).[40] N. Masuda, T. Takaguchi, N. Sato, and K. Yano, in Temporal Networks , edited by P. Holme and J. Saram¨aki(Springer, Berlin, 2013) pp. 245–264.[41] M. Karsai, K. Kaski, A.-L. Barab´asi, and J. Kert´esz, Sci.Rep. , 397 (2012).[42] T. Takaguchi, M. Nakamura, N. Sato, K. Yano, andN. Masuda, Phys. Rev. X , 011008 (2011).[43] J.-P. Eckmann, E. Moses, and D. Sergi, Proc. Natl. Acad.Sci. U.S.A. , 14333 (2004).[44] M. Karsai, K. Kaski, and J. Kert´esz, PLoS ONE (2012). [45] J. Saram¨aki and E. Moro, Eur. Phys. J. B , 164 (2015).[46] B. Lindner, Phys. Rev. E , 022901 (2006).[47] H. Cˆateau and A. D. Reyes, Phys. Rev. Lett. , 058101(2006).[48] D. R. Cox, Renewal theory (Methuen, London, 1962).[49] T. Hiraoka, N. Masuda, A. Li, and H.-H. Jo, Phys. Rev.Res. , 023073 (2020).[50] M. G´enois and A. Barrat, EPJ Data Sci. , 11 (2018).[51] V. Gemmetto, A. Barrat, and C. Cattuto, BMC Infect.Dis. , 695 (2014).[52] J. Stehl et al. , PLoS ONE , e23176 (2011).[53] P. Vanhems et al. , PLoS ONE , e73970 (2013).[54] J. Fournet and A. Barrat, PLoS ONE , e107878 (2014).[55] R. Mastrandrea, J. Fournet, and A. Barrat, PLoS ONE , e0136497 (2015).[56] L. Gauvin, M. G´enois, M. Karsai, M. Kivel¨a, T. Tak- aguchi, E. Valdano, and C. L. Vestergaard, PreprintarXiv:1806.04032 (2018).[57] V. Raghavan, G. Ver Steeg, A. Galstyan, and A. G.Tartakovsky, IEEE Trans. Comput. Social. Syst. , 89(2014).[58] A. E. Clementi, C. Macci, A. Monti, F. Pasquale,and R. Silvestri, in Proc. Twenty-seventh ACM Sympo-sium on Principles of Distributed Computing , edited byR. Bazzi and B. Patt-Shamir (Association for ComputingMachinery, New York, 2008) p. 213.[59] A. E. Clementi, C. Macci, A. Monti, F. Pasquale, andR. Silvestri, SIAM J. Disc. Math. , 1694 (2010).[60] D. T. Gillespie, J. Phys. Chem. , 2340 (1977).[61] K.-I. Goh and A.-L. Barab´asi, EPL , 48002 (2008).[62] A. Feldmann and W. Whitt, Perform. Evaluation , 245(1998).[63] N. Masuda and L. E. C. Rocha, SIAM Rev.60