Correlated bursts in temporal networks slow down spreading

Spreading dynamics has been considered to take place in temporal networks, where temporal interaction patterns between nodes show non-Poissonian bursty nature. The effects of inhomogeneous interevent times (IETs) on the spreading have been extensively studied in recent years, yet little is known about the effects of correlations between IETs on the spreading. In order to investigate those effects, we study two-step deterministic susceptible-infected (SI) and probabilistic SI dynamics when the interaction patterns are modeled by inhomogeneous and correlated IETs, i.e., correlated bursts. By analyzing the transmission time statistics in a single-link setup and by simulating the spreading in Bethe lattices and random graphs, we conclude that the positive correlation between IETs slows down the spreading. We also argue that the shortest transmission time from one infected node to its susceptible neighbors can successfully explain our numerical results.

where μ 1 (μ 2 ) and σ 1 (σ 2 ) are the average and the standard deviation of the first (last) n−1 IETs, respectively. Positive M implies the tendency of large (small) IETs being followed by large (small) IETs. Negative M points towards the opposite, while M = 0 for uncorrelated IETs. In our setup, both P(τ) and M are inputs of the model, requiring us to consider M as a parameter rather than an estimator. Then by controlling the shape of P(τ) and the value of M for interaction patterns between nodes, one can study the effects of correlated bursts on the spreading in temporal networks. By analyzing the contagion dynamics on a single link, and then by simulating the spreading in regular and random temporal networks, we conclude that the positively-correlated inhomogeneous IETs slow down the spreading.

Models
In order to study the spreading dynamics, we consider one of the extensively studied epidemic processes, i.e., susceptible-infected (SI) dynamics 31 : A state of each node in a network is either susceptible (S) or infected (I), and an infected node can infect a susceptible node by the contact with it. Here we assume that the contact is instantaneous. One can study a probabilistic SI (PSI) dynamics, in which an infected node can infect a susceptible node with probability η (0 < η < 1) per contact, as depicted in Fig. 1(a), i.e.
Due to the stochastic nature of infection, multiple IETs can be involved in the contagion, hence the correlations between IETs in the contact patterns can influence the spreading behavior. The case with η = 1 corresponds to the deterministic version of SI dynamics, which we call one-step deterministic SI (1DSI) dynamics: A susceptible node is immediately infected after its first contact with an infected node, see Fig. 1(b). This dynamics can be described by Finally, for studying the effect of correlations between IETs on the spreading in a simpler setup, we introduce two-step deterministic SI (2DSI) dynamics as a variation of generalized epidemic processes [37][38][39][40] , see Fig. 1(c). Here a susceptible node first changes its state to an intermediate state (S′) upon its first contact with an infected node; it then becomes infected after the second contact with the same or another infected node. This can be written as S  I  S  I  I   ,  2 For modeling the interaction structure in a population, we focus on Bethe lattices as networks of infinite size, while regular random graphs and Erdös-Rényi random graphs will be later considered for finite network models. As for the temporal contact patterns, we assume that the contacts between a pair of nodes or on a link connecting these nodes are instantaneous and undirected. Moreover, the contact pattern on each link is assumed to be independent of the states of two end nodes as well as of contact patterns on other links. The contact pattern on each link is modeled by a statistically identical event sequence with inhomogeneous and correlated IETs. For this, the shape of IET distribution P(τ) and the value of memory coefficient M are inputs of our model. Firstly, we consider a power-law IET distribution with a lower bound τ min and an exponential cutoff τ c as follows:

Two-step deterministic contagion. Single-link transmission.
For understanding the spreading behavior on temporal networks, we first focus on how long it takes for the infection to transmit across a single link, say, from a node u to its neighbor v. If u gets infected from its neighbor other than v in time t u , and later it infects v in time t v , the time interval between t u and t v defines the transmission time r ≡ t v − t u . Here we assume that v is not affected by any other neighbors than u, for the sake of simplicity. In order for the infected u to infect the susceptible v, u must wait at least for the next contact with v. This waiting or residual time is denoted by r 0 , see Fig. 1.
For the one-step deterministic SI dynamics, r = r 0 . Due to the independence of contact patterns of neighboring links, we can consider the infection of u to occur at random in time. In addition, since only the IETs larger than r contribute to the probability of having the transmission time r, we have the transmission time distribution as r 1D with μ denoting the mean IET, see Methods for the detailed derivation. The average transmission time is directly obtained as where σ 2 denotes the variance of IETs 6 . Note that a larger variance of IETs results in a larger average transmission time, expected to slow down the spreading. In general, r can be given as the sum of r 0 and subsequent IETs, as depicted in Fig. 1(a,c), for the generalized epidemic processes, including our two-step deterministic SI (2DSI) dynamics. In the case with 2DSI dynamics, the transmission process involves two consecutive IETs. If the infection of u occurs during the IET of τ i , then the transmission time is written as with τ i +1 denoting the IET following τ i . Information on the correlations between τ i and τ i+1 is carried by the joint distribution P(τ i , τ i+1 ) or the conditional distribution P(τ i +1 |τ i ). Using P(τ i +1 |τ i ) with τ i +1 = r − r 0 , Q 1D (r) in equation (4) can be extended to obtain the transmission time distribution for the 2DSI case as In order to relate this result to the memory coefficient in equation (2), we define a parameter as to finally obtain the analytical result of the average transmission time: In the case with M = 0 for uncorrelated IETs, one gets 〈r〉 2D = 〈r〉 1D + μ.
We remark that our result in equation (11) is valid for arbitrary functional forms of IET distributions as long as their mean and variance are finite. M is coupled with σ 2 /μ, implying that the impact of correlations between IETs becomes larger with broader IET distributions. More importantly, we find that a stronger positive correlation between consecutive IETs leads to a larger average transmission time. This can be understood in terms of the role of the variance of IETs in the average transmission time, as shown in the 1DSI case. That is, the variance of the sum of two consecutive IETs is amplified by the positive correlation between those IETs. Based on the result of the single-link analysis, the positive correlation between IETs is expected to slow down the spreading in a population.
Spreading in Bethe lattices. In order to investigate the effects of correlations between IETs on the spreading in a population, we numerically study spreading dynamics in a Bethe lattice, i.e., a regular tree of infinite size in which every node has exactly k neighbors. One can relate this dynamics to the early-stage dynamics in regular random graphs, in which cycles are rare as long as the network size is sufficiently large. As mentioned, the contact pattern on each link is modeled by an independent and identical point process with the same P(τ) and M. Precisely, for each link, we draw n random values from P(τ) to make an IET sequence T ≡ {τ i } i = 1, …,n , for sufficiently large n. Using equation (2), we measure the memory coefficient from T, denoted by M ∼ . Two IETs are randomly chosen in T and swapped only when this swapping makes ∼ M closer to M, i.e., the target value. By repeating the swapping, we obtain the IET sequence whose ∼ M is close enough to M, and from this IET sequence we get the sequence of contact timings for each link. Then the temporal network can be fully described by a set of contact timings for all links. Each simulation begins with one node infected at random in time, which we set as t = 0, while all other nodes are susceptible at this moment. For each simulation, we measure the number of infected nodes as a function of time, I(t). The average number of infected nodes 〈I(t)〉 is found to exponentially increase with time, e.g., as shown in Fig. 2(a): where a = a(k, α, M) denotes the exponential growth rate, known as the Malthusian parameter 44 . a turns out to be a decreasing function of M, indicating the slowdown of spreading due to the positive correlation between IETs, see Fig. 2 The slowdown can be more clearly presented in terms of the relative growth rate a/a 0 with a 0 ≡ a(M = 0) for all cases of k and α, as shown in Fig. 2(e-g). We summarize the main observations from the numerical simulations as follows: (iv) The deviation of a/a 0 from 1 tends to be larger for smaller α.
The observation (i) is expected from equation (11), so is (ii) as both μ and σ 2 /μ decrease with α. (iii) is trivial. (iv) implies that the effect of M becomes larger for smaller α, which can be roughly understood by a larger value of σ 2 /μ coupled to M in equation (11). We remark that equation (11) is the result for a single-link transmission, requiring us to study the transmission time in networks.
In order to comprehensively understand the above observations, in particular, the k-dependence of a, we need to study the effect of time-ordering between infections to different neighbors 45 . For this, we introduce the shortest transmission time as where r (j) for j = 1, …, k − 1 denotes the transmission time from an infected node to its jth neighbor. Here we focus on the average of r s , denoted by 〈r s 〉, which is a function of k, α, and M. In Fig. 3(a), we numerically find that a〈r s 〉 is independent of M, implying that the effect of the correlations between IETs on spreading can be fully understood by 〈r s 〉. Then we write a as follows: Here h(k, α) is generally expected to be a function of k and α, while only its k-dependence is clearly shown in Fig. 3(b), where h(k, α) increases with k. In Fig. 3(c) we observe that as k increases, 〈r s 〉 algebraically decays before converging to a constant, enabling us to assume that s where f, g, and δ are non-negative constants independent of k. By fitting the numerical results of 〈r s 〉 using equation (15), we find how these constants depend on α and M, as summarized in Fig. 3(d-f). Firstly, we find that f is overall independent of M. In the limit of k → ∞, 〈r s 〉 should asymptotically approach the smallest possible transmission time, denoted by r min , leading to f = r min . For the 2DSI dynamics and by our setup, f = τ min = 1 is expected, while the estimated values of f show systematic deviations from 1, possibly due to finite-size effects of k. Secondly and most importantly, g turns out to linearly increase with M such that with positive coefficients g 0 and g 1 , eventually leading to the linear dependence of 〈r s 〉 in equation (15) on M. Moreover, both g 0 and g 1 are found to decrease with α, shown in the inset of Fig. 3(e). These findings are comparable to the analytical result of average transmission time in equation (11). Finally, the estimated values of δ appear to slightly increase with M, while we consider δ to be constant of M in our argument. In sum, we rewrite 〈r s 〉 in equation (15) as s m in 0 1 Combining a in equation (14) and 〈r s 〉 in equation (17), we obtain the relative growth rate a/a 0 as 0 1 min i.e., a/a 0 linearly but slightly decreases with M, showing a good agreement with the numerical results for large α in Fig. 2(g).
Conclusively, it turns out that the statistical properties of the shortest transmission time can account for the numerical results, while more refined approach needs to be taken for better understanding the effect of network structure on spreading, e.g., k-dependence of a in the case of Bethe lattices.

Probabilistic contagion. Single-link transmission.
In a more realistic scenario than the two-step deterministic contagion dynamics, the infection can be described by a stochastic process, i.e., probabilistic SI (PSI) dynamics: An infected node infects a susceptible node with probability η upon contact. Similarly to the deterministic cases, we begin with the analysis for a single-link transmission. If the infection of u occurs during the IET of τ i as depicted in Fig. 1(a), the transmission time for a successful infection of v after l failed attempts for l ≥ 0 is Information on the correlations between τ i , …, τ i+l is carried by the joint distribution P(τ i , …, τ i+l ) for l + 1 consecutive IETs. The distribution of r, denoted by Q P (r), can be written as the weighted sum of transmission time distributions for multi-step deterministic dynamics, similarly done in the previous work 35 : where Q l (r) denotes the distribution of transmission time after l failed attempts. Note that Q 0 (r) = Q 1D (r) in equation (4) and Q 1 (r) = Q 2D (r) in equation (7). Q l (r) for general l ≥ 1 is obtained using the joint distribution P(τ i , …, τ i +l ) as  where it is obvious from equation (21) that τ i ≥ r 0 and 0 ≤ r 0 ≤ r, and δ denotes a Dirac delta function. Then one gets the average transmission time as follows: For the details of the derivation, see Methods. We define the generalized memory coefficient between two IETs separated by j − 1 IETs 36 as We note that this result is valid for arbitrary functional forms of IET distributions as long as their mean and variance are finite. Similarly to the deterministic case in equation (11), the average transmission time for the PSI case turns out to be a linearly increasing function of the memory coefficient M.
Spreading in Bethe lattices. Next, we numerically examine the spreading behavior for the PSI dynamics with η < 1 in Bethe lattices. Similarly to the two-step deterministic case, we observe an exponential growth in the average number of infected nodes as well as the slowdown of spreading when the memory coefficient is positive. For example, the case with η = 0.1 is depicted in Fig. 4. As η approaches 1, the slowdown effect due to the correlated IETs becomes weak, as expected (not shown). Based on the results in Fig. 5, we make overall the same conclusions as in the 2DSI case: a is a decreasing (increasing) function of M (both α and k), and the deviation of a/a 0 from 1 tends to be larger for smaller α. The above observations in the PSI case can be understood by the same argument made for the 2DSI case, namely, the functional form of a in equation (14) with 〈r s 〉 in equation (17), but with some important differences: Firstly, the shortest possible transmission time is r min = 0, although the estimated values of f show systematic deviations from 0 in Fig. 5(d). This deviation is denoted by a small positive value ε, leading to f = ε. Secondly, δ appears to be an increasing function of M rather than a constant in Fig. 5(f), which we assume to be δ = δ 0 + δ 1 M with positive coefficients δ 0 and δ 1 . We therefore modify 〈r s 〉 in equation (17) as follows: We note that due to the positive δ 1 , the above 〈r s 〉 may decrease with M but only for sufficiently large k and M. However, we find no evidence for the decreasing behavior in the ranges of k and M studied in our paper. Using equation (30), the relative growth rate is obtained as Since ε is a small number, we consider only the case of g k   Fig. 4(e). In Fig. 4(e-g), we observe that the difference between curves of a/a 0 for different ks increases and then decreases as α increases from  Spreading in finite networks. So far we have focused on the spreading in Bethe lattices, i.e., regular networks of infinite size, which can also approximate the early-stage dynamics of spreading in finite networks as long as the cycles are rare. In addition to the early stage, the late-stage dynamics of spreading in finite networks has also been of interest 30,46 . For this, we employ two network models of size N: random regular graphs, in which every node has exactly k neighbors, and Erdös-Rényi random graphs, in which every possible pair of nodes is connected with a probability p, hence the average degree is p (N−1). On each of these graphs, both 2DSI and PSI dynamics are tested by the numerical simulations to measure the average numbers of infected nodes as a function of time, 〈I(t)〉. In all cases, we use networks of size N = 10 3 , and the results are averaged over 10 3 simulation runs with different initial conditions for each parameter set. In Fig. 6 we find that the positive correlation between IETs lowers the average number of infected nodes for the almost entire range of time, except for the final-stage dynamics. This final-stage dynamics can be quantified by the time it takes to fully infect the population, denoted by T full . As depicted in the bottom panels of Fig. 6, the average value of T full is clearly increasing with M when simulated in the regular random graphs, consistent with our previous results, while it appears to fluctuate around a constant in the case with ER graphs. Such fluctuating behavior of T full in the case with ER graphs can be understood in terms of the variation of degrees of nodes, in contrast to the regular random graphs with a fixed degree for all nodes: At the final stage, it can take a longer time for some lower-degree susceptible nodes to get infected from their infected neighbors. This large timescale for the infection of low-degree nodes and its large fluctuations can eventually dominate over the early-and intermediate-stage dynamics that show clear dependence on M. In another work using the apparently power-law IET distributions without exponential cutoff 33 , the positive correlation between IETs was reported to reduce T full , calling for more systematic approaches.

Conclusions
Spreading dynamics in temporal networks has been extensively studied for tackling the important issue of what features of temporal networks are most relevant to the speed of spreading taking place in such networks. One of the widely studied features is the inhomogeneity of interevent times (IETs), typically represented by heavy-tailed IET distributions, in the temporal interaction patterns between nodes. Although the impact of the inhomogeneous IETs on the spreading has been largely explored, yet little is known about the effects of correlations between IETs on the spreading. It is partly because many previous works have considered the contagion dynamics with the assumption of the immediate infection upon the first contact between susceptible and infected nodes, hence without the need to consider the correlated IETs. However, since temporal correlations in the interaction patterns can be fully understood both by IET distributions and by correlations between IETs 18 , the effects of inhomogeneous and correlated IETs on the spreading need to be systematically studied for better understanding the dynamical processes in complex systems. For this, we consider two contagion dynamics, i.e., two-step deterministic SI and probabilistic SI dynamics, naturally involving multiple consecutive IETs. For both dynamics, we derive analytical expressions of average transmission times 〈r〉 for a single-link setup, which turn out to linearly increase with the memory coefficient M as shown in equations (11) and (29). Therefore, the positive correlation between IETs is expected to slow down the spreading. By performing numerical simulations of the contagion dynamics in regular networks of infinite size and random graphs of finite size, we conclude that the positive correlation between IETs indeed slows down the spreading, compared to the case of uncorrelated IETs but from the same IET distributions.
The numerically obtained spreading speed, e.g., in Bethe lattices of degree k, could be successfully explained by means of the statistics of the shortest transmission time among k − 1 transmission times from one infected node to its k − 1 susceptible neighbors. In the case when IETs in the interaction patterns are largely homogeneous, the average transmission time 〈r〉 will serve as a representative timescale that determines the spreading speed. However, in the other case with inhomogeneous IETs or heavy-tailed IET distributions, the transmission time to each of k − 1 neighbors will be heterogeneously distributed, implying that neighbors infected earlier tend to spread the disease or information more quickly, hence more broadly, than those infected later. In this sense, the majority of the infected nodes can be largely explained by the descendants of early-infected neighbors, and the characteristic timescale of spreading speed can also be set by the average shortest transmission time 〈r s 〉, rather than 〈r〉. Unfortunately, as the analysis of 〈r s 〉 appears not to be straightforward, more detailed and rigorous understanding of the behavior of 〈r s 〉 is left for future works.
So far we have focused on the positive memory coefficient based on the empirical observations, while the effects of the negative memory coefficient on spreading could also be studied. Our analytic results of average transmission times in equations (11) and (29) and numerical results on temporal networks enable us to predict the accelerated spreading due to the negativity of the memory coefficient. We also note that the shape of IET distributions may restrict the range of memory coefficient 42 , which should be properly considered when designing models.
Finally, we remark that in addition to the memory coefficient, the correlations between IETs have also been identified by other methods, e.g., in terms of bursty trains, which can detect long-range memory effects between IETs 5 : The number of events in each bursty train, i.e., the burst size, has been described by heavy-tailed distributions. Regarding this, the relation between memory coefficient and burst size distributions was recently studied 19 . Our approach can be extended by incorporating such heavy-tailed burst size distributions. We also note that more realistic network structures can be adopted for modeling temporal networks, such as networks with heterogeneous degrees 47 and community structure 48 among other network properties, e.g., stylized facts in social networks 49 .

Derivation of transmission time distributions.
We first derive the transmission time distribution for the single-link transmission in the case with 1DSI. The transmission time r is equivalent to the time interval between a random arrival and the next event or contact. The probability that a random arrival falls within an IET of τ is τP(τ)/μ, where the mean IET μ is for the normalization. Given the arrival within the IET of τ, the probability of the transmission time being r is 1/τ. Then integrating such probability over the possible range of τ ≥ r gives the transmission time distribution as ∞ ∞ The transmission time distribution for the 2DSI case can be derived similarly. Given the random arrival within an IET of τ i and the residual time r 0 to the next event, the transmission time r is r 0 + τ i+1 . The probability of having τ i+1 is given by the conditional distribution P(τ i+1 |τ i ), i.e., P(r − r 0 |τ i ). Together with the probabilities we already have for the 1DSI case, we obtain the transmission time distribution for the 2DSI case by integrating the probability over τ i ≥ r 0 and 0 ≤ r 0 ≤ r as  The last expression is rewritten by introducing the integration with respect to τ i+1 :