Modeling temporal networks with bursty activity patterns of nodes and links

,


I. INTRODUCTION
Temporal networks have become an important framework to understand the dynamics of complex systems over the last decade [1][2][3].By integrating the topological knowledge of a system, described by a graph, with the information about the temporal nature of the interaction between its components, represented by time series, we can precisely track who interacts with whom and when.The interaction dynamics can be captured at several different levels.First, the interaction between each pair of nodes specifies the dynamics of the link.Second, by aggregating the interaction between a node and all of its neighbors, one obtains the dynamics of the node, which shows how the node interacts with others.Lastly, by collecting all the interaction between every pair of nodes, one can tell about the dynamics of the entire system.For instance, in communications systems where people interact by sending messages, the link dynamics corresponds to the message correspondence between a pair of individuals, while the node dynamics corresponds to the inbox of messages sent or received by an individual.
Human communication is known to exhibit non-Markovian, inhomogeneous temporal patterns which are commonly referred to as being bursty [4,5].When each communication event is instantaneous or lasts for a short period so that its duration can be neglected compared to other time scales, one can regard the communication sequence as a realization of a point process.The burstiness * hang-hyun.jo@apctp.org of a point process is mainly characterized by a heavytailed distribution of time intervals between consecutive events, or inter-event times (IETs), in contrast to Poisson processes for which the IET distributions are exponential.
Interestingly, empirical data suggest that in communications systems, the communication sequences of nodes and of links are both characterized by power-law distributions with a similar scaling exponent [6,7].This cannot be taken for granted for the following reason.As mentioned above, the communication sequence of a node is the superposition of the communication events on all the links between the node and its neighbors.However, a superposition of independent renewal processes does not retain the statistics of the original processes in general.In fact, the IET distribution for the superposed process tends to an exponential distribution in the limiting case where the number of independent source processes is large [8][9][10].Therefore, the observation that both node dynamics and link dynamics are bursty suggests the presence of correlations across communication processes on different links.Such link-link correlations can have a significant impact on the dynamical processes taking place in the network [7,[11][12][13], but their origin has yet to be understood.
Here, we study the mechanisms behind the burstiness in node and link activity patterns by considering a model in which the nodes are activated randomly in time with non-Poissonian statistics and two nodes may communicate if and only if they are simultaneously activated.In Sec.II, we introduce two variants of the model with different communication rules.In Sec.III, we report results of numerical simulations performed on networks with var-ious topologies.We show that, for both models and for all the networks, the communication patterns are characterized by heavy-tailed IET distributions for both nodes and links.Section IV is devoted to explaining the origin of the burstiness in node and link activity patterns for each model.We describe the behavior of the model in a system of two nodes by relating it to the statistics of the sum of a random number of random variables.We use the same approach to derive the activity patterns of nodes and links in larger networks.Finally, we conclude our work in Sec.V.

II. MODEL
We consider a network of size N with a given structure, in which each node is activated randomly at discrete times and its activation pattern is modeled by a renewal process with a given inter-activation time (IAT) distribution, denoted by P (r).In order to start the activation process at equilibrium, the first activation time t 0 ≥ 0 of each node is assigned according to the residual time distribution [14] where r denotes the average IAT.The node is then activated at times t l = t 0 + l l =1 r l for l = 1, 2, • • • , where each IAT, denoted by r, is independently drawn from P (r) [15].In our work, we adopt a power-law IAT distribution, where ζ(α) ≡ ∞ x=1 1/x α is the Riemann zeta function.We choose α > 2 to make Eq.(1) converge.
At each time step t (0 ≤ t ≤ T ), pairs of activated nodes communicate with each other.As depicted in Fig. 1, here we consider two variants of the model.The first variant, which we call the polyvalent model, assumes that an activated node communicates with all the activated neighbors.The case where an activation does not lead to communication only occurs when the node does not have any simultaneously activated neighbors.In the second variant, which we refer to as the monovalent model, an activated node is randomly paired with one of its activated neighbors to have at most one communication partner at the same time, as is the case for one-to-one phone calls.An activated node cannot communicate with others if none of the neighbors are simultaneously activated, or if all the simultaneously activated neighbors are already paired with other nodes.
In both models, the communication events between a pair of adjacent nodes can be associated not only with the link but also with the nodes.In other words, we can define a communication event sequence for each link as well as for each node, the latter being the union of the communication events on all the links attached to the node.Hereafter, we refrain from the wording "inter-event time" to avoid confusion and instead use "inter-communication time" (ICT) to represent the time interval between consecutive communication events on the sequence affiliated with a node or with a link.Note that the communication sequence of a node does not agree with the activation sequence of the node in general, because an activated node may not communicate with anyone, as shown in Fig. 1.We denote the ICTs by τ .Whenever necessary, subscripts i and ij will distinguish between node i's properties and link ij's properties; sub-or superscripts p and m will indicate variables and functions related to the polyvalent and monovalent models, respectively.In the following sections, we discuss the statistics of node and link ICTs.

III. NUMERICAL RESULTS
We carry out numerical simulations for synthetic networks with different topologies such as complete graphs, random regular graphs, and scale-free graphs, as well as Zachary's karate club network [16] as an example of a real-world network.Figure 2 summarizes the node and link ICT distributions ψ(τ ).Here we set α = 2.5.
The polyvalent model yields the node and link ICTs both of which are distributed almost indistinguishably from the power-law IAT distributions for the various network topologies.In contrast, the monovalent model results in different communication patterns depending on the network structure.For homogeneous graphs such as complete and random regular graphs, the node ICT distributions are almost identical to those of IATs, while the link ICT distributions show a hump at short time scales that is not in the power-law IAT distributions.As  we make the network sparser by reducing the degree of random regular graphs, the hump becomes smaller and the range of τ in which the distribution approximately follows a power law becomes wider, implying that the sparseness of networks plays an important role in realizing bursty communication patterns on links.
For scale-free graphs, both the node and link ICTs differ from the IATs in terms of the distribution.To examine the effect of structural heterogeneity, we group the nodes by degree and consider the node ICT distribution for each group.Figure 3(a) shows that the deviation between the ICT and IAT distributions is larger for nodes with smaller degrees.This can be understood intuitively as follows: When activated, a node with more neighbors is more likely to find a communication partner, i.e., a neighbor who is simultaneously activated and available.In the extreme case where the degree of a node is large enough so that the node almost always finds a partner whenever activated, the node ICT distribution would coincide with the IAT distribution.On the other hand, the ICT distribution would deviate more from the IAT distribution for nodes with less neighbors because of the difficulty of finding a partner.
Similarly, we group the links by the product of the degrees of the end nodes of the link, k i and k j , and measure the distribution of τ ij conditioned on each value of k i k j .As shown in Fig. 3(b), the ICT distributions for links with larger values of k i k j show larger deviations from the IAT distribution.Suppose that two adjacent nodes are simultaneously activated.If the degrees of the nodes are larger, then the nodes are likely to have larger numbers of other simultaneously activated neighbors that they potentially communicate with.Therefore, the probability that the two nodes communicate with each other decreases.In contrast, if the degrees are small, then the two nodes are more likely to communicate with each other.We discuss these intuitions more quantitatively in the following Section.
Finally, for Zachary's karate club network, the results are similar to those for the homogeneous graphs.The node ICT distribution shows a small deviation from the IAT distribution, which is due to the degree heterogeneity.

IV. ANALYTICAL RESULTS
In this section, we provide an analytical examination of the behavior of the models.For that, we start by considering a minimal system that consists of a pair of adjacent nodes, which we call a dimer.We show that the number of IATs that compose each ICT is a random variable that follows a power-law distribution.Thus, we can describe an ICT as the sum of a random number of random variables, and we show that the sum is also distributed as a power law.Then, we derive the statistics of the link and node ICTs for the polyvalent model directly from the results for a dimer.Finally, we describe the monovalent communication as a result of random success and failure of the polyvalent communication.This leads to expressions of link and node ICTs for the monovalent model as a geometric sum of polyvalent link and node ICTs, respectively.We obtain power-law statistics in this case as well.

A. Case of dimers
A dimer is a pair of nodes connected only to each other.Under this setup, the monovalent model is equivalent to the polyvalent model because the nodes have no other nodes to communicate with except for each other.Moreover, the communication sequences of the two nodes are identical to each other as well as to that of the link connecting them.
In general, each of the two nodes of the dimer, denoted by i and j, can be activated more than once between two consecutive communication events, as sketched in Fig. 4(a).This leads to expressions of an ICT, denoted by τ , as the sum of successive IATs of each node: where r ω,n denotes the n th IAT of node ω (ω = i, j) within the ICT and n ω denotes the number of times that node ω is activated between the two communication events.We call n ω an inter-communication activation number (ICAN) of node ω.The random variables n i and n j will have the same statistics by symmetry, which allows us to focus on node i's point of view from now on.Keeping Eq. ( 3) in mind, our goal is (i) to derive the statistics of n i and (ii) to compute the statistics of τ as the sum of the n i independent random variables r i .
Let us consider the activation processes of the two nodes between two consecutive communication events as follows (see the right panel of Fig. 4(a)).Node j is activated and communicates with node i at time t j,0 for the first time, and then activated at t j,1 , t j,2 , . . ., t j,nj −1 until it communicates with i again at t j,nj for the second time.The number of activations of node i between the two consecutive activations of node j at time t j,n −1 and t j,n is denoted by ñi,n , where 1 ≤ n ≤ n j .The ICAN is then written as the sum of the numbers of activations in each segment indexed by n , In order to derive the distribution of ICANs, we map the renewal process of the activations of node j to an inhomogeneous Bernoulli process, in which the activation probability is a time-dependent parameter.Particularly, we adopt the framework of mapping a continuous renewal process into an event-modulated Poisson process [17].An event modulated Poisson process is a process where the event rate λ is independently redrawn from a distribution F (λ) after every event and remains constant until the next event occurs with that rate.The cumulative IET distribution is then shown to be the Laplace transform of F (λ) [17].
In our case, we should instead consider an activationmodulated Bernoulli process since time is discrete.In this framework, node j is activated at each time step with a probability λ j , which is independently redrawn from a distribution F (λ j ) upon every activation.Then, from node i's point of view, each of its ñi,n activations between two consecutive activations of node j can be considered as an independent Bernoulli trial with the success probability λ j,n that node j is activated at the same time (see the right panel of Fig. 4(a)).Now, we hypothesize that a large n i is likely to be dominated by the numbers of activations that occur within a few long IATs of node j governed by small activation rates.This leads to a simplification of the argument: instead of calculating n i by a combination of processes with different activation rates, we regard the process between two communication events as almost entirely homogeneous and approximate the distribution of large n i by a geometric distribution with parameter In other words, we replace the activation-modulated process by a communication-modulated process.The distribution of ICANs is then given as where G(λ j ) is the distribution of the activation probability for each ICT.The approximation in the second line is derived from the fact that the tail behavior of φ(n i ) will be dominated by the contributions from small λ j .We put a small finite cutoff and perform the integration up to this value.In order to estimate the functional form of G(λ), we go back to the original activation-modulated picture.For a given activation rate λ, the IATs are distributed as p(r|λ) = λe −λr . ( t j,0 t j,1 t j,2 t j,3 t j,4 λ j,3 λ j,4 4. Schematic illustrations of communication sequences between a pair of nodes i and j.The red and white rectangles represent activations with and without communication, respectively.The solid vertical lines identify communication events between the pair.In (b) and (c), nodes are activated with the same temporal pattern as in the left panel of (a).(a) The case where the two nodes form a dimer.Variable ni denotes the number of activations (represented by rectangles) of node i between two communication events (red rectangles).The right panel in (a) shows how the activation pattern of node j can be mapped to an activation-modulated Poisson process where each IAT is associated with an activation probability λj.Variable ñi denotes the number of activations of i in each period segmented by activations of j.(b) The polyvalent model.The communication pattern on the link is the same as that of the dimer, while each node also communicates with other neighbors and experiences more frequent communication with smaller ICANs and ICTs, which are denoted by n p i and τ p i , respectively.(c) The monovalent model.Simultaneous activation of adjacent nodes may fail to trigger communication, as indicated by light blue rectangles and vertical dotted lines.Variable mij denotes the number of events that nodes i and j are simultaneously activated (represented by vertical lines) between two consecutive communication events on link ij (solid vertical lines), while mi denotes the number of activations of i that concurred with activations of any of its neighbors (colored rectangles) between two consecutive communication events of i with one of its neighbors (red rectangles).See main text for details.
Conversely, the distribution of λ conditional on r is given by the Bayes' theorem, As we do not assume anything about the prior p(λ) except λ > 0, we adopt the non-informative density, which is uniform throughout its domain.Eq. ( 8) then reads with the normalization factor r 2 .When the IAT distribution scales as P (r) ∼ r −α at the tail, the following scaling holds for small values of λ: This is consistent with the fact that for the eventmodulated Poisson processes, P (r) will have a power-law tail with exponent α when F (λ) is a gamma distribution with shape parameter α − 1 [17], which scales as F (λ) ∼ λ α−2 for small λ.We assume that G(λ) F (λ) for 0 < λ ≤ and plug the scaling into Eq.( 6) to obtain This derivation tells us that the statistics of the ICANs of node i is determined by the activation process of node j.If the IAT distributions for node i and j are characterized by different exponents α i and α j , respectively, then φ(n i ) ∼ n

−αj i
. We now turn to our second question regarding the statistics of ICTs as the sum of an ICAN of IATs.Since the ICAN and IAT are independent random variables, we exploit the analytical results in Ref. [18]: We consider the following sum where the summands r and the number of summands n are independent random variables and they both follow power-law distributions, P (r) ∼ r −α and φ(n) ∼ n −β .Then, τ also asymptotically obeys a power-law distribution ψ(τ ) ∼ τ −α where α = min{(α − 1)(β − 1) + 1, α, β}.
In our case, since the IAT r i and ICAN n i in Eq. ( 3) are shown to have the same scaling exponent as β = α, the ICT distribution also follows a power law with the same exponent α = α, that is,

B. Polyvalent model
The polyvalent model assumes that communication occurs on a link every time when the two end nodes are activated at the same time, irrespective of the states of other nodes in the system.Therefore, the dimer case discussed in the previous subsection directly translates into the communication patterns on links (compare the left panel of (a) to (b) in Fig. 4).Indeed, Fig. 5(a) shows that the link ICAN distributions follow power laws, and that the scaling given by Eq. ( 11) agrees well with the numerical results.The distribution of polyvalent link ICTs is the same as that of ICTs for the dimer case and given by In order to investigate the node communication processes, we introduce a time frame defined by counting the number of activations of a node, formally expressed as where t is the wall-clock time, t i are the times that node i is activated, and θ(•) is the Heaviside step function.
Simply put, the activation-based time ν i is measured by a clock that ticks one unit forward upon every activation of node i.This time transformation t → ν i (t) rescales an ICT τ = t − t into an ICAN n i = ν i (t ) − ν i (t ), meaning that an ICAN is an "inter-communication time" for the processes in the time frame ν i .We note that similar concepts of time frame transformation, named "relative clock" and "activity clock," are used in recent studies [19][20][21].
In the wall-clock time frame, the communication processes on adjacent links ij and ij are correlated because of the underlying activation process of node i.However, in the activation-based time frame ν i , in which the activations of node i are regularized, the two communication processes are independent because when i is activated, communication between nodes i and j depends only on whether j is simultaneously activated and is not affected by the behavior of node j .Since the communication process on each link is characterized by power-law distributed ICANs as in Eq. ( 11), the node communication process as the superposition of independent link processes has a thinner-tailed ICAN distribution, that is, with β > α (see Fig. 5(a)).
By taking the same approach as in the previous subsection, we write a node ICT as follows: where P (r) ∼ r −α .Then, by following the scaling relation given by Eq. ( 13), we find that the polyvalent node ICTs are distributed as

C. Monovalent model
In contrast to the polyvalent model, simultaneous activation of two adjacent nodes do not necessarily trigger communication between them in the monovalent model.In order to study the statistics of link ICTs, we first need to account for the random pattern of successful communication when a pair of adjacent nodes are simultaneously activated.
Suppose that node i with degree k i is activated at a time step along with κ i + 1 activated neighbors including node j.Because activation processes of different nodes are independent of each other, variable κ i is binomially distributed as where ρ = 1/ r denotes the probability that each neighbor of node i is activated when node i is activated.If κ i > 0, the communication between nodes i and j occurs only if i selects j as the counterpart as a result of random matching.Although the probability of selecting each of the activated neighbors is not uniform in general, we assume the uniformity for simplicity so that the probability that node j is selected is equal to 1/(κ i + 1).The same goes for node j for its κ j + 1 activated neighbors including node i.Then, the probability that simultaneous activation of nodes i and j leads to communication between them is approximated by This form reduces to q ij = 1 for the dimer case of k i = k j = 1, in which they communicate with each other every time they are simultaneously activated.In the limit where k i , k j 1, we have q ij 1/ρ 2 k i k j .Equation ( 21) is, on the whole, numerically supported as shown in Fig. 5(c), although deviations and fluctuations are notable.We think these deviations originate from perturbation by higher-order effects involving more than two nodes, which violates the uniformity assumption.
Let m ij be the number of times that adjacent nodes i and j are simultaneously activated between two consecutive communication events, including their activation that triggered the latter of the two communication events but excluding the one that triggered the former (see Fig. 4(c)).We call m ij an inter-communication coactivation number (ICCAN) of link ij.Because random pairing is done independently at each time step, m ij is geometrically distributed (see Fig. 5(d)) as As depicted in Fig. 4(c), a monovalent link ICT, denoted by τ m ij , is equal to the sum of m ij successive polyvalent link ICTs, The distribution of τ m ij is written as where is the probability that a monovalent link ICT is equal to τ and it is composed of m polyvalent link ICTs.Here δ(•) denotes the Dirac delta function.
An analytical evaluation of the discrete power-law distribution ψ p (τ ) is not straightforward.Instead, we consider a continuous counterpart given by where α > 1.The Laplace transform of Eq. ( 26) is given as where Γ(•, •) denotes the upper incomplete gamma function.In the asymptotic limit of s → 0, one gets , where Γ(•) is the gamma function, and c ≡ (α−1)/(2−α).By only keeping the leading terms of expansion with respect to s, we have The inverse Laplace transform of Eq. ( 29) in the limit of τ → ∞ yields Equation ( 30) is valid for any values of α because considering higher-order terms in the expansion of Eq. ( 29) does not affect the asymptotic form given by Eq. ( 30).This result indicates that the link ICT distribution for the monovalent model has a power-law tail with the same exponent as the link ICT distribution for the polyvalent model, which is consistent with the numerical results presented in Fig. 5(b).At the same time, the geometric distribution of ICCANs contributes to the hump at the bulk part of the monovalent link ICT distribution.For a dense network where the degrees are generally large, q ij is small and the geometric decay of χ in Eq. ( 22) is slow; this is why the hump is larger in denser networks as shown in Fig. 2. Lastly, we discuss the node ICT distribution for the monovalent model.Unlike the polyvalent case, the monovalent communication events on adjacent links (i.e., links sharing a node) are not independent of each other even in the activation-based time frame because communication between a pair of nodes forbids the nodes to communicate with other nodes at the same time.Nevertheless, Fig. 5(e) shows that node ICCANs m i , i.e., the numbers of times that node i is activated simultaneously with any of its neighbors until it communicates with one of them, are geometrically distributed.This observation suggests that the probability that a node succeeds to communicate with another node is constant every time it is simultaneously activated with at least one of its neighbors.
A monovalent node ICT, denoted by τ m i , can be written as the sum of m i successive polyvalent node ICTs as follows: Equation ( 31) is analogous to the relation between the monovalent and polyvalent link ICTs given by Eq. ( 23).Again using the scaling relation given by Eq. ( 30), one obtains the monovalent node ICT distribution with a power law at its tail as follows: This result is in good agreement with the numerically obtained scaling relations shown in Fig. 5(b).

V. CONCLUSION
In order to explain the origin of the bursty activity patterns of nodes and links observed in empirical communication systems, we have proposed a temporal network model where the nodes communicate with each other according to their non-Poissonian random activation.The two variants of the model that we discussed are both able to reproduce heavy tails in the inter-communication time (ICT) distributions for nodes and links for various network topologies.We have shown that the polyvalent ICTs are power-law distributed because each of them is a sum of inter-activation times (IATs) where both the summands and the number of summands are power-law distributed random numbers, which stem from the node activation processes.A monovalent ICT is described as a sum of polyvalent ICTs, where the number of summands is geometrically distributed because the activated nodes are paired randomly and independently at each time step.As a result, an exponential hump appears in the bulk part of the ICT distributions, especially prominently for small-degree nodes and for links between largedegree nodes; nevertheless, the tail part of the distributions follows a power law with the same scaling exponent as the power law in the IAT distribution.
The superposition of independent communication sequences with power-law distributed ICTs does not yield a sequence with the ICTs distributed as a power law with the same scaling exponent.Therefore, the assumption of independence between links is unable to account for the real-world observations.Our results suggest a possible mechanism behind the reconciliation between bursty dynamics of nodes and of links: Link-link correlations emerge as a result of underlying node activations, each of which may or may not realize actual communication.
Further steps can be taken in this line of research.In this work, we have considered a homogeneous population of nodes that shares the same activation statistics.In reality, the activity levels of nodes and the weights of links, i.e., the frequency of communication between pairs of individuals, are heterogeneous [22].It would be straightforward to include such heterogeneity into our modeling framework if we consider nodes endowed with different scaling exponents of the IAT distribution.We have also assumed that a node behaves in an equal manner toward every neighbor.However, empirical data show that individuals allocate their efforts to communicate with others unevenly among alters [23].This effect can be taken into account in the monovalent model by setting biases toward certain links when pairing communication partners.Another possible extension is to implement communication among a group of nodes, which corresponds to "conference calls" or "group chats," in a direction similar to Ref. [24].One can also tailor the temporal structure of node activation patterns to account for the empirical observation of long-range correlated node ICTs in human communication [25,26].Future work also includes how the presence of link-link correlations affects dynamical processes taking place in temporal networks and the associated network control problems [27].

FIG. 1 .
FIG. 1. Schematically illustrated snapshots of communication according to the two models given the same set of activated nodes, which are enclosed by thick solid lines.The red filled circles and red solid lines represent the nodes and links with a communication event, respectively.(a) The polyvalent model assumes that the activated nodes communicate with all the neighbors that are simultaneously activated.(b) In the monovalent model, each activated node communicates with at most one neighbor.

FIG. 2 .
FIG.2.The inter-communication time (ICT) distributions ψ(τ ) of nodes and links for the polyvalent model (top panels) and for the monovalent model (bottom panels).The network structures are, from left to right, a complete graph (a, f), random regular graphs with degree k = 15 (b, g) and with k = 6 (c, h), a scale-free graph with degree distribution ∝ k −2.1 (d, i), and Zachary's karate club network (e, j).The parameters used are N = 100 and T = 10 7 for the complete and random regular graphs, N = 1000 and T = 10 6 for the scale-free graph, and N = 34 and T = 10 7 for Zachary's karate club network.The inter-activation time (IAT) distribution, which follows Eq. (2) with α = 2.5, is represented by the dashed line.

FIG. 3 .
FIG. 3. (a) Node ICT distributions grouped by degree ki.(b) Link ICT distributions grouped by the product of the degrees of the end nodes, kikj.Both figure panels are based on the simulation results for the monovalent model on the scale-free graph shown in Fig. 2(i).

FIG. 5 .
FIG. 5. (a, b)The scaling exponents of ICT and ICAN distributions, each for nodes and links, as a function of the scaling exponent α of IAT distributions in the cases of (a) the polyvalent and (b) monovalent models, respectively.The results are obtained from numerical simulations run for T = 4 × 10 7 time steps on a random regular graph with N = 100 and k = 6.The dotted lines represent identity.The insets show the ICAN distributions when α = 3.2.(c) The probability qij of communication between pairs of simultaneously activated adjacent nodes of degree ki and kj in a scale-free graph.The parameters are the same as in Fig.2(i).The continuous lines, for each ki, represent the theoretical curves given by Eq. (21) with a numerically obtained value of ρ. (d, e) The distributions of (d) link ICCANs and (e) node ICCANs, both obtained from simulations of the monovalent model with α = 3.2.The network topology is the same as in (b).