Generative models of simultaneously heavy-tailed distributions of inter-event times on nodes and edges

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 as contagion processes have been intensively studied. The literature supports that such heavy-tailed distributions are present for inter-event times associated with both individual nodes and individual edges in networks. However, the simultaneous presence of heavy-tailed distributions of inter-event times for nodes and edges is a non-trivial phenomenon, and its origin has been elusive. In the present study, we propose a generative model and its variants to explain this phenomenon. We assume that each node independently transits between a high-activity and low-activity state according to a continuous-time two-state Markov process and that, for the main model, events on an edge occur at a high rate if and only if both end nodes of the edge are in the high-activity state. In other words, two nodes interact frequently only when both nodes prefer to interact with others. The model produces distributions of inter-event times for both individual nodes and edges that resemble heavy-tailed distributions across some scales. It also produces positive correlation in consecutive inter-event times, which is another stylized observation for empirical data of human activity. We expect that our modeling framework provides a useful benchmark for investigating dynamics on temporal networks driven by non-Poissonian event sequences.

Heavy-tailed distributions of IETs are commonly found for single nodes [6,7,37,38,[41][42][43] and single edges [1,4,9,44,45].Such a distribution for a node implies that the sequence of event times for the node, regardless of the identity of the neighbor, obeys a non-Poissonian, heavy-tailed statistics.In fact, the distribution of IETs for both nodes and edges in a single data set are often heavy-tailed (see Section II for examples).However, the presence of heavy-tailed IET distributions for both edges and nodes in the same network is not trivial.Consider a node, denoted by v, that interacts with its k neighbors, and assume that the sequence of IETs on each edge incident to v independently obeys a heavy-tailed distribution.The superposition of the k event sequences yields the sequence of v's events.This situation is illustrated in Fig. 1.In Fig. 1(a), node v has k = 5 neighbors.We have produced the sequence of events on each of the five edges assuming a power-law distribution of IETs, as shown in Fig. 1(b).The sequence of events shown in the bottom of Fig. 1(b) is that for v, which one obtains by superposing the sequences of events on the k edges.In general, the distribution of IETs for v is less heavy-tailed than that for a single edge.This is because the superposition of independent point processes (more precisely, renewal processes) obeying heavy-tailed distributions roughly approaches, albeit not precisely, a Poisson process as one increases k [46,47]. of events generated by a power-law distribution of IET on each edge of the star network and the sequence of events on node 1 in (a).The CV of IETs, which we calculate from the first approximately 10 6 events on each edge, is shown in (b) for the five edges and node 1.We used the power-law distribution of IETs for edges given by p(τ ) = (α − 1)/(1 + τ ) α , where α = 3.5.According to an equilibrium renewal process [3,48], we draw the time to the first event on each edge from the distribution of waiting times given by p Recently, we proposed a model that generates heavytailed distributions of IETs for both edges and nodes [49].The model assumes that each node is activated at discrete times according to a renewal process that draws the inter-activation time from a power-law distribution.Then, at each time step, activated nodes are uniformly randomly selected to communicate with simultaneously activated neighbors.By construction, this model produces a heavy-tailed distribution of IETs for individual nodes.Although it is less trivial, IETs for edges also obey approximately heavy-tailed distributions.However, this model does not explain why we find heavy-tailed distributions of IETs for both nodes and edges in empirical data.Realistic mechanisms of the simultaneous presence of heavy-tailed IET distributions for nodes and edges are underexplored [4].
In the present study, we propose a model of timestamped event sequences on networks that is based on a latent state dynamics of nodes.In our model, each node switches between two states called the high-activity and low-activity states according to a Markov process in continuous time.We assume that if both nodes are in the high-activity state, events occur according to a Poisson process at a higher rate, otherwise at a lower rate.The rationale behind the model is that events between two individuals may happen more frequently when both individuals are motivated to interact with others than otherwise.We show that our model produces distributions of IETs with large dispersions for both single nodes and edges, resembling empirical data.

II. SIMULTANEOUSLY LARGE VARIABILITY OF INTER-EVENT TIMES ON NODES AND EDGES IS COMMON AND NON-TRIVIAL
Before presenting and analyzing our model, in this section we provide empirical evidence that heavy-tailed distributions of IETs are simultaneously present for edges and nodes in single data sets.We use data collected by the SocioPatterns collaboration [50].The survival function of the IETs (i.e., probability that the IET, τ , is larger than the specified value) for individual edges and nodes for social contact data in a primary school [50], which we refer to as PrimarySchool, is shown in Fig. 2(a), and 2(b), respectively.(See Appendix A for the results for the other data sets.)We only used edges with at least 100 events to calculate the distribution and the following statistics.The relatively slow decay in Fig. 2 suggests heavy-tailed distributions for both edges and nodes across some scales of τ .
We quantified the dispersion of IET distributions by the coefficient of variation (CV).The CV of a distribution is defined as the standard deviation divided by the mean.For an IET distribution, one can write where • represents the average over edges or nodes.A Poisson process produces an exponential IET distribution, which yields CV = 1.A periodic process yields CV = 0.A heavy-tailed distribution yields a large value of CV.Table I shows the mean and standard deviation of the CV of the IET distribution for edges and nodes for each data set.In this and the following analyses of the empirical data, we treated the data with two steps.First, we aggregated consecutive events with 20 s duration; the temporal resolution of the original data set is 20 s, i.e., the social contacts are measured every 20 s.We perform this first step because we are analyzing the IETs but not the duration of events.For instance, we aggregated a sequence of event times {20, 40, 60, 100, 200, 220} (in s) between two nodes into a sequence with three 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 = 100 s 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, we removed 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 seems to be common.We hypothesize that the correlation of the IET between different edges sharing a node contributes to large values of the CV for individual nodes.Therefore, for each data set, we uniformly randomly shuffled the IETs on each edge within each day of recording, preserving the time of the first and last events of each day on the edge.This shuffling method is equivalent to the timeline inter-event shuffling in an instant-event temporal network (formally named P[π L (∆τ ), t 1 ]) described in Ref. [56].Because this shuffling procedure preserves the distribution of IETs on each edge, the edge's CV is unchanged.However, it affects the IETs and hence the CV for individual nodes.The node's CV calculated from the shuffled data and the relative deviation defined by ∆ = 100% × (shuffled − original)/original, are shown in Table I.For all data sets, the CV for the nodes consistently decreases when one shuffles the IETs.Therfore, the shuffling destroys the heavy-tailed nature of the IET sequences for individual nodes.
To further support that the simultaneous presence of heavy-tailed distributions of IETs on edges and nodes is nontrivial, we assess a common approach to generate event sequences on each edge according to an independent renewal process with a power-law distribution, p(τ ), of IETs.We call this model the stochastic temporal network model [3].Consider the power-law distribution of IETs given by where α is a parameter.The first and second moments of p(τ ) are given by τ = 1/(α − 2), where α > 2, and , where α > 3, respectively.Therefore, the CV of IETs is given by where α > 3. We ran simulations with α = 3 and α = 3.5 for a node with k = 2, k = 5, and k = 10 neighbors, and calculated the CV for IETs on the node for each combination of α and k.Equation (3) predicts an edge's CV equal to 2.24 when α = 3.5 and its divergence for α ≤ 3, which is consistent with the numerical results shown in Table II.By contrast, the CV for the node is considerably smaller than that for edge and decreases towards 1 as k increases.Therefore, the stochastic temporal network model does not produce simultaneously heavy-tailed distributions of IETs for edges and nodes.

III. MODEL
We propose a model of node behavior that aims to simultaneously produce large CVs for IETs on individual edges and nodes.Consider a static network.The model is based on two main assumptions.First, we assume that each node stochastically switches between two states, called the high-activity state (denoted by h) and the low-activity state (denoted by ).The plausibility of this assumption is supported by various empirical and modeling studies [35-39, 41, 57-59].Second, we assume that 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's CV calculated after the timeline inter-event shuffling.The data originate from a primary school (PrimarySchool) [51,52], a scientific 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 the following 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.II.CV of IETs on edges and nodes obtained from the stochastic temporal network model.We assume an equilibrium renewal process, so we draw the time to the first event on each edge from the distribution of waiting times, as in Fig. 1.

PrimarySchool
We calculated the mean and standard deviation of the CVs on the basis of 100 realizations of the simulation.We stop each realization when all edges have obtained at least 10 6 events.
have a contact event much more likely when both nodes are in the high-activity state than otherwise.The intuition behind this assumption is that a pairwise human contact event may be much more likely to occur when both individuals are motivated to interact than otherwise.
Each node switches between h and according to a two-state continuous-time Markov process.In other words, the node switches to the opposite state according to a Poisson process whose rate depends on the current state.We denote by r h→ the rate at which a node in state h changes to state , and similar for r →h .We assume that different nodes share the same r h→ and r →h values but are associated with independent twostate Markovian processes.
If both adjacent nodes are in the h state, then events on the edge are assumed to occur according to a Poisson process at a higher rate denoted by λ h .Otherwise, the edge produces events according to a Poisson process at a lower rate denoted by λ (< λ h ).This process is schematically shown in Fig. 3.
The master equation for the probability that a node is in state h, denoted by p h (t), where t represents time, is given by Schematic illustration of the model.The events between two nodes occur at a higher rate λ h if and only if both are in the high-activity state (time windows shown in red) Otherwise, events occur at a lower rate λ (shown in gray).
Therefore, the stationary probability of finding a node in state h is given by The stationary probability of finding a node in state is IV. RESULTS

A. The model produces simultaneously large variability of IETs on edges and nodes
We numerically simulated the model to generate sequences of IETs and computed the survival function and the CV of IETs for individual edges and nodes.Apart from the structure of the static network, our model has four parameters, r h→ , r →h , λ h , and λ .Equation (5) yields r →h = r h→ p * h /(1 − p * h ).We write λ in terms of λ h as where 0 < γ < 1.
Simultaneously changing (r h→ , r →h , λ h , λ ) to (cr h→ , cr →h , cλ h , cλ ), where c > 0, is equivalent to not changing these four parameters and changing the time from t to ct.Therefore, without loss of generality, we set λ h = 1, unless we state otherwise.In the following simulations, for fixed values of r h→ , we varied γ and p * h .Using the Gillespie algorithm [60], we generated events on the edges until all edges had at least 10 6 events.
We used a star network composed of a node v and its k neighbors.Initially, each of the k + 1 nodes is independently in state h or with probability p * h or (1 − p * h ), respectively.Then, we run a continuous-time Markov process independently for each node.At each time, depending on the state of each edge, we generate events on the edge at rate λ h or λ .It should be noted that a next event on an edge may not be simply produced as a single Poisson process.For example, suppose that the two nodes connected by an edge are both in state h.Then, the time to the next event is drawn from the exponential distribution p(τ ) = λ h e −λ h τ .However, if either node switches to state before the next event occurs, then one has to discard the scheduled time to the next event, which was generated from p(τ ) = λ h e −λ h τ , and redraw the time to the next event from p(τ ) = λ e −λ τ .
We consider a star network composed of a node v and its two neighbors as an example.We set r h→ = 2×10 −5 , r →h = 4.7 × 10 −5 , λ h = 6 × 10 −3 , and λ = 3.5 × 10 −4 , which yields p * h ≈ 0.7 and γ ≈ 0.06.The survival function of IETs produced by the model is shown by the red lines for the two edges and node v in Fig. 2(a) and Fig. 2(b), respectively.The survival function for both the two edges and v decays more slowly than exponentially, roughly consistent with the non-Poissonian behavior observed in the empirical data (black lines in Fig. 2).For our model, the CV for the two edges is equal to 2.8 and 3.0, and that for v is equal to 2.6.These values of CV are 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→ (i.e., r h→ = 10 −4 , r h→ = 10 −3 , and r h→ = 10 −2 ) and three values of k (i.e., k = 2, k = 5, k = 10).We show the CV values in Fig. 4. The figure indicates that the model produces large CV values simultaneously for edges and nodes in a broad parameter region.In particular, the CV is large when γ (= λ /λ h ) is small and p * h is near 0.7, across the range of r h→ and k.

B. Analytical evaluation of the CV of inter-event times
In this section, we provide an analytical account for the CV values observed in section IV A. The state of the edge is specified by the states of two nodes forming the edge.We denote by h 2 when both nodes are in state h, and likewise for h and 2 .The edge state obeys a three-state Markov process in the state space S 1 = {h 2 , h , 2 }, where we do not distinguish between h and h, and represent both by h .In our model, the IET distribution conditioned on the edge's state is given by g(τ |h 2 ) = λ h e −λ h τ and g(τ | 2 ) = g(τ |h ) = λ e −λ τ , where g(τ |•) represents the distribution of IETs conditioned on the edge's state.
The Markov process of node activity is independent for different nodes.Therefore, two nodes are simultaneously in state h in the equilibrium with probability p * 2 h .The mean number of events produced when both nodes are in state h is given by λ h p * 2 h T , where T is the observation time.The mean number of events produced when either node is in state is given by λ (1 − p * 2 h )T .Therefore, an IET is produced at rate λ h and λ with probability . By combining these contributions and ignoring IETs during which the edge's state changes, we obtain the probability density function (PDF) of IETs for an edge as a mixture of two exponential distributions as The first two moments of this PDF are given by and (9) By substituting Eqs. ( 8) and (9) into Eq.(1), we obtain where CV edge is the CV for the edge's IETs.Now we consider a node v with k neighbors.From the viewpoint of node v, the sequence of events is a superposition of the events over its k edges.Because the k nodes are statistically the same for v, it is sufficient to consider a 2(k + 1)-state Markov process with state space where the first set 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 probability p * h in the equilibrium, the IET distribution for v depends on the states of v's neighbors.The probability of finding exactly k h neighbors in state h is given by where k k h is the binomial coefficient.In this situation, v experiences events produced by a Poisson process at rate k h λ h + (k − k h )λ because the superposition of the k independent Poisson processes on edges with rate λ h or λ is a Poisson process with the summed rate.Second, if v is in state , which happens with probability 1 − p * h , all edges produce events at rate λ .In this situation, v experiences events produced by a Poisson process at rate kλ .Therefore, an event By combining these contributions, we derive the PDF for IETs on a node with k neighbors as The first two moments of the PDF are given by and By substituting Eqs. ( 12) and ( 13) into Eq.( 1), we obtain the CV for node v as Equation ( 14) reduces to Eq. ( 10) when k = 1.
For any k, Eq. ( 14) yields lim γ→0 CV k → ∞ and lim γ→1 CV k = 1.In addition, the derivative of Eq.( 14) with respect to γ is negative for any 0 < γ < 1.Therefore, CV k monotonically decreases towards 1 as γ → 1, which is consistent with Fig. 4. Next, equating the derivative of the right-hand side of Eq. ( 10) with respect to p * h to zero yields p * h = 1/ √ 2 ≈ 0.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 ≈ 0.7 regardless of k and γ.Therefore, the CV k is large when γ small and p * h ≈ 0.7, which is consistent with Fig. 4.
To assess the accuracy of the theory, we calculated the relative error defined by [(theoretical CV) -(numerical CV)]/(numerical CV).Figures 5(a)-5(d) show the relative error when r h→ = 10 −4 .In this case, the relative error is small across the entirety of our parameter region.When r h→ = 10 −3 (see Figs. 5(e)-5(h)), and r h→ = 10 −4 (see Figs. 5(i)-5( )), the relative error is large when γ is small and p * h is large.The relative error increases as r h→ increases for the following reason.A large value of r h→ and moderate value of p * h (i.e., p * h that is not too close to 0 or 1) implies large r →h , because r →h = r h→ p * h /(1 − p * h ).When r →h is large, the mean IET for edges that are produced at event rate λ , which is equal to 1/λ , is longer than the typical duration of the low-activity state of the edge (i.e., either h or 2 ), which is proportional to 1/r →h .Note that, the low-activity state of an edge finishes only when both of the two nodes have transited from state to h, which occurs at rate r →h for each node.Therefore, an edge is populated with sufficiently many IETs produced in the low-activity state of the edge if and only if 1/r →h 1/λ , which leads to r →h λ .According to our parametrization, we obtain r →h = r h→ p * h /(1 − p * h ) and λ = γλ h = γ, because we set λ h = 1.Therefore, an edge is populated with sufficiently many IETs produced in its low-activity state if r h→ (1 − p * h )γ/p * h .For example, when γ = 0.01 and p * h = 0.9, we need r h→ (1 − p * h )γ/p * h ≈ 10 −3 .This condition is violated when r h→ = 10 −3 or 10 −2 .For these r h→ values, the probability that events occur in the low-activity state of the edge is small.However, our analytical derivation of the CV assumes that sufficiently many events and hence IETs occur in the typical duration of both lowactivity and high-activity states of the edge at respective rates, i.e., λ and λ h .This explains the discrepancy between the theoretical and numerical results observed in Figs.5(e)-( ).
However, the model is still capable of producing large CV values even if few IETs are produced in the lowactivity state of the edge.In this situation, the IET between the last event of a high-activity period of the edge and the first event of the next high-activity period would generate a long IET, contributing to a relatively heterogeneous distribution of IETs.This regime is not predicted by our analytical solution.
If r h→ ≈ λ h , then few IETs are produced during the h 2 state of the edge on average.Therefore, a large CV value requires r h→ λ h , which is satisfied in Figs. 4  and 5 because the largest value of r h→ that we use is 10 −2 and we have set λ h = 1.

C. Correlation between consecutive inter-event times
In human activities, IETs for both edges and nodes are often positively correlated, i.e., long IETs tend to be followed by long IETs and vice versa [1,4,41,61].To examine this property, we compute the memory coefficient, M , of a sequence of IETs [61], defined as where n is the number of IETs in the sequence, m 1 and σ 1 are the average and standard deviation of {τ 1 , τ 2 , . . ., τ n−1 }, respectively, and m 2 and σ 2 are the average and standard deviation of {τ 2 , τ 3 , . . ., τ n }, respectively.The memory coefficient measures the correlation coefficient of consecutive IETs, (τ i , τ i+1 ).
Figure 6 shows the memory coefficient for the empirical data.In all data sets, the edges show predominantly positive memory coefficients with values lying mostly between 0 and 0.1 (see Fig. 6(a)).The memory coefficient for nodes is also predominantly positive and tends to be larger than that for the edges (see Fig. 6(b)).These values are in accordance with previous results for various human activities [61].
The memory coefficient for sequences of IETs generated by our model is shown in Fig. 7.The figure indicates that the model produces positive M for both edges and nodes, which is qualitatively consistent with the empirical data.However, the memory coefficient values produced by the model are substantially larger than the empirical values.

D. Variants of the model
In our original model, events occur on the edge at a higher rate if and only if both nodes forming the edge are in state h.In this section, we study two variants of the model.In the first variant, we assume that event on an edge occur at the higher rate, λ h , if either node connected to the edge, not necessarily both nodes, is in state h.Events on the edge occur at the lower rate λ if and only if both nodes are in state .We call this variant the OR model.An interpretation of the OR model is that, if an individual wants to interact with a neighbor, he/she can do so at the higher event rate regardless of whether or not the neighbor wants to interact.
For the OR model, we calculated the CV for IETs on edges and nodes by scanning the same values of p * h , γ, r h→ , and k as those used in Fig. 4. The results are shown in Fig. 8.As in the original model, the OR model produces large CV values for both edges and nodes when γ is small.With respect to p * h , the OR model produces large CV values when p * h ≈ 0.3 for the edge and that 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 that the OR model produces large CV values for the node shrinks as k increases.
In the second variant of the model, we assume that the node's h and states independently contribute λ h and λ , respectively, to the event rate of the edge.In other words, events occur on the edge at rate 2λ h if both nodes are in the h state, 2λ if both nodes are in the state, and λ h + λ if one node is in the h state and the other node is in the state.An interpretation of this variant of the model, which we call the IND model (named after "independent"), is that the state of each individual independently contributes to the frequency of events between two individuals.
The CV values for the IND model are shown in Fig. 9. Similarly to the original and the OR models, the IND model produces large CV values when γ is small.For any given γ, the IND model produces large values of CV when 0.3 p * h 0.4 for the edge.For the node, the p * h value that maximizes the CV decreases as k increases.This behavior is consistent with the analytical result shown in Appendix C. Therefore, the IND model behaves similarly to the OR model, i.e., it produces large CV values when γ is small and p * h ≈ 0.3, and the parameter region in which the node's CV is large shrinks as k increases.
To quantitatively compare the three models, we calcu- ) pairs for which the CV value is larger than two; a CV value larger than two is consistent with the results for empirical data (see Section II).A large fraction value implies that a CV value larger than two is robustly produced for various parameter combinations.The result for this analysis is shown in Fig. 10(b), where we again set r h→ = 10 −4 and vary k.The figure indicates that the original model has a larger fraction of the parameter region with CV larger than two than the OR and IND models.The fraction remains roughly constant as k increases for the original model, whereas it rapidly decreases as k increases for the OR and the IND models.Figures 10(a) and 10(b) altogether suggest that the original model is more capable of producing large CVs of IETs on both edges and nodes than the OR and IND models, particularly when the node has a large degree.The results are qualitatively the same for larger r h→ values, i.e., r h→ = 10 −3 (Figs.10(c) and (d)) and r h→ = 10 −2 (Figs.10(e) and (f)).

V. DISCUSSION
We have started from the observation that heavy-tailed distributions of IETs are simultaneously present for both individual nodes and edges in the same empirical data.We have proposed a continuous-time model and its variants for generating discrete events on edges that replicate this behavior to different extents.Our main model crucially assumes that each node alternates between highactivity and low-activity states in a Markovian manner.We showed that the original model, which requires that both nodes are in the high-activity state for the edge to have frequent events, is capable of producing large CV values for both individual nodes and edges in a broad parameter region.The other two variants of the model are also capable of producing reasonably large CV values to some extent.The proposed models allow interpretations.For example, in the original model, two nodes are likely to interact if and only if both of them feel like interacting with others.
We have derived analytical solutions for our models by discarding the effects of IETs that contain state transitions of the edge.The analytical solution was accurate when two conditions were met.First, the ratio of the high-to low-activity event rates (i.e., 1/γ) should not be extremely large, such as 100.This condition is probably not unrealistic.Second, the state transition rate of the node should be sufficiently small compared to the event rates on edges.Under this condition, an epoch of the high-or low-activity state of an edge is typically long enough to host sufficiently many events at the constant rate, which is either λ h or λ .This is the situation that the theory in Section IV B assumes.Otherwise, a large fraction of IETs contains transitions of the edge's state, which cause a systematic discrepancy of the analytical expression from the numerical results.
In our models, each node switches between two states.A possible extension of this assumption is to the case of more than two states for each node.Then, depending on how such a model translates the nodes' states into the event rate, the distribution of IETs on edges may be approximately a mixture of more than two exponential distributions, which may resemble or actually produce heavy-tailed distributions [36,57,62,63].In fact, a mixture of a small number of exponential distributions, including the case of just two exponential distributions, is often sufficient for approximating many empirical heavytailed distributions of IETs [36,37,57,62].Therefore, one should carefully assess trade-offs between the complexity of extended models and the explanatory power of the model that one gains by assuming more states for nodes.
The distribution of IETs affects how disease and information spread across contact networks [8][9][10][11][12][13][14][15].Because many time-stamped event data probably have heavytailed distributions of IETs for both individual nodes and edges, the temporal network models proposed in the present study are expected to be useful for modeling dynamical processes on temporal networks including contagion processes.It seems that model-based studies of contagion processes on temporal networks have not paid much attention to the simultaneous presence of heavytailed distributions of IETs on nodes and edges [4].How this property affects key indicators of contagion processes such as the epidemic threshold, the final epidemic size, and equilibrium fraction of infected nodes, as well as indicators of other dynamical processes, warrants future work.lower rate λ if and only if the two nodes forming an edge are in state .Therefore, using the same reasoning as that for the original model in Sec.IV B, one obtains the PDF of IETs for an edge as follows: where ) and, for convenience, we have used p * = 1 − p * h .The first two moments of this PDF are given by (B3) By substituting Eqs.(B2) and (B3) into Eq.(1), we obtain where CV OR edge is the CV for the edge's IETs for the OR model.
We proceed with the same steps as those in Sec.IV B to calculate the CV for a node with k neighbors as follows.The PDF for a node with k neighbors is given by where k is the number of neighbors in the state, and The first two moments of f OR k (τ ) are given by and By substituting Eqs.(B6) and (B7) into Eq.( 1), we obtain the CV for node v as which generalizes Eq. (B4).
The behavior of Eq. (B8) with respect to γ is qualitatively the same as that of the original model (i.e., Eq. ( 14)).In other words, lim Therefore, the CV for the edge is large when γ is small and p * h ≈ 0.3.By contrast, the value of p * h that maximizes Eq. (B8) for a given γ value strongly depends on k.In Fig. 12(a), we numerically inspect this dependence.The value of p * h that maximizes Eq. (B8) decreases as k increases, and, for a given value of k, it monotonically in- The PDF for IETs on a node with k neighbors is given by where The first two moments of f IND k (τ ) are given by and By substituting Eqs.(C6) and (C7) into Eq.( 1), we obtain the CV for node v as which generalizes Eq. (C4).The behavior of Eq. (C8) with respect to γ is qualitatively the same as that of the original model.In other words, lim which is plotted as the blue line in Fig. 12(b).The numerically obtained p * h value that maximizes Eq. (C8) decreases as k increases, and increases as γ increases (see Fig. 12(b)).These results are consistent with Fig. 9.

FIG. 1 .
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 node 1 in (a).The CV of IETs, which we calculate from the first approximately 10 6 events on each edge, is shown in (b) for the five edges and node 1.We used the power-law distribution of IETs for edges given by p(τ ) = (α − 1)/(1 + τ ) α , where α = 3.5.According to an equilibrium renewal process[3,48], we draw the time to the first event on each edge from the distribution of waiting times given by p w (t) = (α − 2)/(1 + t) α−1 .

FIG. 2 .
FIG.2.Survival function, P (τ ), of IETs on (a) edges and (b) nodes for the PrimarySchool data set (black) and simulation of our model (red).In this and the following analyses of the empirical data, we treated the data with two steps.First, we aggregated consecutive events with 20 s duration; the temporal resolution of the original data set is 20 s, i.e., the social contacts are measured every 20 s.We perform this first step because we are analyzing the IETs but not the duration of events.For instance, we aggregated a sequence of event times {20, 40, 60, 100, 200, 220} (in s) between two nodes into a sequence with three 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 = 100 s 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, we removed IETs larger than eight hours.In both (a) and (b), we only considered edges that had at least 100 events.

FIG. 6 .
FIG. 6. Box plots of the memory coefficient for (a) edges and (b) nodes for the empirical data sets.The box shows the median (orange line), the first quartile (Q1), and the third quartile (Q3); the whiskers show the minimum (Q1 − 1.5 × IQR) and the maximum (Q3 + 1.5 × IQR) values excluding outliers, where IQR = Q3 − Q1.The open circles are the outliers.The green triangles are the sample means.

FIG. 10 . 8 .IND 1 2λ h e −2λ h τ + + 2 ( 1 (IND 1 =FIG. 11 .
FIG. 10.Comparison of CV values obtained from the original, OR, and IND models.Panels (a), (c), and (e) show the largest CV value in the (γ, p * h ) parameter region explored in Figs.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 the case of the single edge.We set r h→ = 10 −4 in (a) and (b), r h→ = 10 −3 in (c) and (d), and r h→ = 10 −2 in (e) and (f).

= 1 ,
and the CV IND k monotonically decreases as γ increases, which are consistent with Fig. 9.The p * h value that maximizes Eq. (C4) is given by