Small inter-event times govern epidemic spreading on temporal networks

Naoki Masuda 2 and Petter Holme Department of Mathematics University at Buffalo, State University of New York, Buffalo, NY 14260-2900, USA Computational and Data-Enabled Science and Engineering Program, University at Buffalo, State University of New York, Buffalo, NY 14260-5030, USA∗ Tokyo Tech World Research Hub Initiative (WRHI), Institute of Innovative Research, Tokyo Institute of Technology, Yokohama 226-8503, Japan


I. INTRODUCTION
Despite the continuous advances in medicine, infectious diseases is still a major burden to the global health. It is, however, an area where physics modeling can be of tangible support to society. Theoretical epidemiology has developed several core concepts that is guiding today's medical epidemiologists, such as: epidemic thresholds, herd immunity and the basic reproductive number [1]. Temporal network epidemiology [2]-studying how structures in the time and network of contacts between people affect the spread of infectious diseases-is an emerging area that can help improving epidemic forecasting and intervention.
One salient feature of many empirical contact data sets is that the times of contacts cannot be described as a Poisson process [3]. Specifically, the interevent times (i.e., the time between contacts), both of pairs of individuals and of the individuals themselves, follow rightskewed fat-tailed distributions, which are wider than the exponentially distributed ones of a Poisson process. Another (even more well-studied) similarly shaped distribution in contact networks is that of the degree (number of neighbors) [4]. The heterogeneities of these distribution affect the dynamics of epidemic outbreaks. As for degree-distributions, fat tails are believed to facilitate disease spreading [5]. This effect can be attributed to the highest-degree individuals. They shrink the distances in the networks, and have many more people who can infect them and whom they can infect. For interevent times, the effect of a fat-tailed distribution can depend on details of the model-some papers find heterogeneous distributions facilitates spreading [6][7][8][9][10], whereas others find the opposite effect [11][12][13][14][15][16].
Fat-tailed distributions can be approximated by a mixture of a small number of exponential distributions [17][18][19]. We use this observation to develop an analytical framework to calculate epidemic thresholds for the SIR * naokimas@buffalo.edu model when inter-event times obey a distribution whose variance can be tuned from Poissonian tails to infinitely fat ones.

II. MODEL
We start by analyzing the following susceptibleinfected-recovered (SIR) model with a general distribution of inter-event times attached on a static network with N nodes. Each node is either susceptible, infected or recovered at any point of continuous time. Each node i (1 ≤ i ≤ N ) carries an independent and identical point process whose probability density of inter-event times, denoted by τ , is given by ψ(τ ). When an event occurs at the ith node, the node is activated and selects one of its neighbors, denoted by j, uniformly randomly. The network is assumed to be undirected and unweighted. If either the ith or jth node is infected and the other is susceptible, the susceptible node becomes infected with probability β. Any infected node recovers according to a Poisson process with rate µ. Once a node has recovered, it will not be reinfected or infect other nodes. The present SIR model is node-centric in the sense that infection events are triggered by activation of nodes. The same [8][9][10] or similar [6] node-centric infection mechanisms have been used in the literature. If ψ(τ ) is exponential, each node process is Poissonian which reduces our model to the standard Markovian SIR model. In this case, the epidemic threshold in the well-mixed population is given by 2β/( τ µ) = 1, where the factor 2 results from the fact that each edge is selected for possible infection at a rate of 2/ τ owing to the activation of either of the two nodes incident to the edge.

III. RESULTS
To analytically examine effects of a fat-tailed ψ(τ ) on the epidemic threshold, we assume an infinite well-mixed population and represent fat-tailed distributions by a mixture of two exponential distributions with rates λ and λ h (≥ λ ). In other words, we set where 0 < a < 1 is the mixture weight. We refer to Eq. (1) as a two-component exponential mixture model (EMM2). A large τ tends to be produced when the exponential distribution with the lower rate, λ , is selected with probability a, and vice versa. Across a scale of τ , EMM2 mimics a fat-tailed ψ(τ ) if a is small and λ λ h . Although EMM2 has, strictly speaking, a short tail, it is often practically good at approximating empirical distributions of inter-event times [17][18][19]21].
We analyze the SIR dynamics combined with EMM2 inter-event times by assuming every node to take either of the two states, low or high, depending on which of λ or λ h is currently employed for generating the next activation time. Each node is assumed to transit between the two states. Upon its activation, each node interacts with a randomly selected different node, possibly transmitting the infection between them. Then, the activated node draws its new rate (i.e., either λ or λ h ), hence the new state, according to which the time to the next activation of the node is generated. Because inter-event times are independently drawn from Eq. (1) each time, the node upon the activation visits the low and high states with probability a and 1 − a, respectively, regardless of the current state of the node.
We denote the fraction of the susceptible nodes in the low and high states among the N nodes by S and S h , respectively, and the fraction of the infected nodes in the low and high states by ρ and ρ h , respectively. The fraction of the recovered nodes is equal to 1−S −S h −ρ − ρ h . The master equation representing the SIR dynamics is given by where t denotes the time. The first term on the righthand side of Eq. (2a) corresponds to infection; λ S (ρ + ρ h )β is the product of the rate λ S at which susceptible nodes in the low state is activated, the probability ρ + ρ h that an infected node is selected as partner to interact with, and the probability β that infection happens. The rate (λ ρ + λ h ρ h )S β is the product of the rate λ ρ + λ h ρ h at which infected nodes are activated, the probability S that a susceptible node is selected as partner, and the probability β that infection happens. The second term on the right-hand side of Eq. (2a) represents the transition of a susceptible node from the high to the low state. The third term represents the transition of a susceptible node from the low to the high state. Equations (2b), (2c), and (2d) were similarly derived.
To derive the epidemic threshold, we consider the initial condition in which most nodes are susceptible and a small fraction of nodes is infected. The zerothorder continuous-time Markov process representing the time course of each node's state has the transition rate from the low to high state of λ (1 − a) and that from the high to low state of λ h a. Therefore, the stationary state of the two states before an infinitesimal fraction of infection is introduced to the population is given by We perturb this stationary state by setting S = S * + ∆S , S h = S * h + ∆S h , ρ = ∆ρ and ρ h = ∆ρ h , where ∆S , ∆S h , ∆ρ , and ∆ρ h are small and satisfy ∆S + ∆S h + ∆ρ + ∆ρ h = 0 due to the conservation of the mass at t = 0. Denote by J the Jacobian matrix of the linearized dynamics around the disease-free initial condition. If any of the eigenvalues of J is positive, a macroscopic number of infected (and recovered) nodes will result, which is the condition that determines the epidemic threshold. As shown in the Supplemental Information, the condition under which a large-scale epi-  [20]. For PL1, we set α = 1.5 or 2.5. Assuming the equilibrium, we drew the distribution of waiting times, i.e. the time to the first activation from ψ w (τ ) = (α − 1)/(1 + τ ) α , independently for each node. For the Pareto distribution in (a), we set α = 1.3 and α = 1 + 6/5 and drew the time to the first activation for each node; we also set τ0 = 1 without loss of generality, because changing τ0 corresponds to rescaling the time. For EMM2, we set a = 0.3, CV = √ 5, and drew the time to the first activation from For each of the six distributions of inter-event times, we scanned values of µ to span a range of λ eff . We measured Ω as the average over n = 5, 000 simulations. Error bars are omitted for clarity. They would be 3-5 times thicker than the lines. demic spreading can occur is given by where µ c is the epidemic threshold in terms of the recovery rate, is the mean inter-event time for each node and Equation (3a) indicates that, when λ = λ h such that ψ(τ ) is an exponential distribution, one obtains b = 0 such that the epidemic threshold in terms of the effective infection rate [22], i.e., λ eff ≡ 2β/ τ µ, is equal to 1. Crucially, Eq. (3a) also implies that µ c > 2β/ τ , or equivalently, λ eff < 1, if and only if b > 0. This is because the two activation rates are expressed as and λ h > 0 guarantees γ > 0. Because b represents the extent to which τ is heterogeneously distributed (related to the coefficient of variation, CV, via b = (CV 2 − 1)/2), these results mean that infection spread more easily for heterogeneously distributed than exponentially distributed inter-event times, within the framework of EMMs.
To validate the theory, we numerically integrated the master equations (Eqs. (2a)-(2d)) with the initial conditions where ρ 0 = 10 −3 is the initial fraction of infected nodes.
For the subsequent analysis, we set a = 0.1 and β = 1.
The final fraction of recovered nodes, also known as the final outbreak size, denoted by Ω, is estimated as the fraction of recovered nodes at a sufficiently large time, t = 10 3 . In Fig. 1, we show Ω as a function of the effective infection rate 2( τ µ) −1 and b. Note that, for EMM2 it holds that b ≥ 0 [23] and b < (1−a)/a = 9. Figure 1 confirms the theoretically derived epidemic threshold (solid line) and indicates that Ω increases with b. This result is consistent with the previous numerical results obtained with the same [9] or similar [6] node-centric SIR models, the theoretical results derived from branching processes under the assumption of tree-like networks [8], and the theoretical results derived for the SIS model for a variant of the activity-driven network model of temporal networks [10]. Equation (3a) suggests that µ c diverges (i.e. λ eff → 0) as γ → ∞. For such a distribution of inter-event times, a large-scale epidemic spreading can occur even when λ eff is tiny. This phenomenon is reminiscent of the vanishing epidemic threshold in the Poissonian SIR and SIR dynamics in scale-free networks [5]. However, the mechanism is different. In the present model, the epidemic threshold approaches zero when b → (1 − a)/a, λ → a/ τ , and λ h → ∞. Therefore, the epidemic threshold vanishes because infinitesimally short inter-event times are produced with a finite probability, 1 − a, not because a fat-tailed distribution occasionally produces large values as in the case of epidemic spreading in static scale-free networks [5]. A similar result was found in Ref. [8].
To further establish the importance of the short-τ end of the distribution, we carried out simulations of the SIR dynamics on the complete graph with N = 10 3 nodes with different power-law ψ(τ ) having different lower bounds. We started each simulation from the initial condition in which one randomly selected node is infected and all the other nodes are susceptible. We drew the time to the first activation of each node from the waiting-time distribution corresponding to ψ(τ ) [8,24].
Consider a power-law distribution of inter-event times, ψ(τ ) = α/(1 + τ ) α+1 , which we denote by PL1. PL1 has support τ ≥ 0 and therefore has some probability mass near τ = 0, being able to produce short inter-event times with some probability. With α = 2.5, PL1 has b = 1/(α − 2) = 2/3 and CV = α/(α − 2) = √ 5. The Pareto distribution is another popular power-law distribution, is given by ψ(τ ) = α τ α 0 /τ α +1 for τ ≥ τ 0 and ψ(τ ) = 0 otherwise, and is not capable of producing inter-event times shorter than τ 0 . For this distribution, one obtains b = (α − 1) 2 / [2α (α − 2)] and CV = [α (α − 2)] −1/2 . The b value of 2/3 is produced by α = 1 + 6/5 ≈ 2.095. Thus, we can compare two power-law distributions sharing the b (and CV) value but are considerably different at both small and large τ . In Fig. 2, we show Ω for a range of λ eff for the two power-law ψ(τ ) and the exponential ψ(τ ) with equal mean τ . The figure shows that that PL1 enhances epidemics as compared to the Poisson case, at least near the epidemic threshold, whereas the Pareto distribution suppresses epidemics. We believe that this is because the Pareto distribution is not capable of generating short inter-event times, whereas PL1 is. Note that the present Pareto distribution has a longer tail than PL1, i.e., α < α. To further support this claim, we also simulated the case in which ψ(τ ) is EMM2 with a = 0.3 and CV = √ 5, which automatically determines λ and λ h via Eqs. (4a) and (4b). The outbreak size with this ψ(τ ) is shown in Fig. 2. We find that Ω is much larger than in the case of power-law distributions with the same mean across the entire range of λ eff . We emphasize that EMM2 does not have a fat tail in a rigorous sense. Qualitatively, the same results hold true even when τ 2 diverges. Results for PL1 with α = 1.5 and those for the Pareto distribution with α = 1.3 are shown by the dotted lines in Fig. 2. To summarize these results, although we expressed our theoretical results in terms of b, a larger value of b or a fat tail of ψ(τ ) is not a key factor impacting the epidemic threshold or the outbreak size. It is the size and the probability mass of τ at small τ values that play the key role. We also confirmed qualitatively similar results for networks with a fat-tail degree distribution generated by the Barabási-Albert (BA) model [25] and the largest connected component of an email social network [20] (Figs. 2(b) and 2(c)).
Many previous studies of the effect of inter-event times on epidemics found that fat-tailed distributions slow down the dynamics [11-13, 15, 16]. These studies all used edge-centric models, where the next contact is generated for individual edges rather than nodes. In edgecentric models, delay or reduction in epidemic spreading is attributed to the waiting-time paradox, with which a newly infected node has to wait for long time (i.e., averaging waiting time larger than the mean inter-event time, τ ) before infecting other nodes. In the Supplemental Information we show that also in this situation it is the short-time end of the distribution of inter-event times that controls the dynamics. With extreme distributions having a large probability mass near 0, the outbreak size for a standard edge-centric SIR model can be even magnified compared to the Poissonian case with the same mean τ .

IV. CONCLUSIONS
The literature about how fat-tailed interevent time distributions affect epidemic spreading on networks has been inconclusive. In the present study, we have resolved this issue by attributing the ease of spreading to the balance between the waiting-time paradox and the abundance of short inter-event times. Our results account for both cases where infection events are triggered by node activation (i.e., node-centric) or edge activation (i.e., edge-centric). In particular, in the node-centric SIR model (that is only weakly influenced by the waiting-time paradox), we have found that it is the frequency of short inter-event times, not a fat tail of ψ(τ ), that determines the epidemic threshold and Ω. We have reached these conclusions via a novel technique using a mixture of exponential distributions, capable of mimicking fat-tailed distributions. We expect that the same technique is applicable for analyzing other dynamical processes such as stochastic opinion dynamics, complex contagions, and coevolutionary dynamics in empirical and model temporal networks. This approach also lends itself well to simulation studies as the dynamics for well-mixed populations is exactly solvable by integrating the master equations. One can achieve even more realism with the use of exponential mixture models with more than two components [17,19].  For the UMM2, we set a = 1, which yields ≈ 1.389.

Appendix A: Derivation of the epidemic threshold
The SIR dynamics linearized around the disease-free steady state is given by where J is the Jacobian matrix, and The leading principal minor of J of order 2 gives two eigenvalues of J, which are equal to 0 and − [λ (1 − a) + λ h a] < 0. The zero eigenvalue reflects the constraint ∆S + ∆S h + ∆ρ + ∆ρ h = 0.

Appendix B: Edge-centric SIR model
In this Letter we conclude that for the node-centric SIR model a fat-tailed ψ(τ ) can enhance spreading, and this is caused by the short-time end of the distribution. In contrast, most of the literature support that epidemic spreading is suppressed by a fat-tailed distribution of inter-event times [11-13, 15, 16]. In fact, the SIR or the susceptible-infected (SI) models considered in those studies are substantially different from the present SIR model in that events are generated on edges as opposed to nodes. Therefore, we call them edge-centric epidemic processes, akin to classification of random walks on temporal networks [26,27]. In edge-centric epidemic processes, fat-tailed ψ(τ ) is known to suppress epidemic spreading due to the so-called waiting-time paradox [13]. In these processes, if a node v i infects its neighbor v j , the time t to the next infection event on edge (v j , v ), where = i, obeys the waiting-time distribution given by rather than ψ(τ ) [26,28,29]. The mean waiting time is given by τ 2 /(2 τ ), where · is the expectation under ψ(τ ), and it is always larger than the half of the mean inter-event time, τ /2, unless ψ(τ ) is the delta distribution.
In particular, when ψ(τ ) is fat-tailed, the mean waiting time is much larger than τ /2 because τ 2 τ 2 . Then, typical attempts to relay infection to a new node would take much longer time than τ . This is the reason why fat-tailed ψ(τ ) suppresses epidemic spreading in edge-centric epidemic processes.
In contrast, the node-centric SIR model, which we mostly focus on in the main text, circumvents the waiting-time paradox. If a susceptible node v i is activated and infected from v j , then v i draws a time to the next activation from distribution ψ(τ ), not from ψ w (τ ), potentially infecting v i 's other neighbors. In this manner, infection can spread rapidly in arbitrary networks.
We compare the final size for the edge-centric SIR model when ψ(τ ) is exponential or PL1 with α = 2.5, sharing the mean, τ = 1, in Fig. 3(a). Confirming the results of the previous studies, the final size is substantially smaller in the case of PL1 than the exponential distribution. Given our insights we have obtained for the node-centric SIR model, is it possible to enhance epidemic spreading in the edge-centric SIR model by devising ψ(τ )?
A waiting-time distribution ψ w (τ ) may have a certain amount of probability mass at small τ depending on ψ(τ ). We consider that in conventional edge-centric SIR models, the effect of the waiting-time paradox outweighs that of the probability mass of waiting times at small τ such that epidemic spreading would be suppressed by a long-tailed ψ(τ ). Therefore, we considered a ψ(τ ) which is little influenced by the waiting-time paradox while it has a relatively high probability mass at small waiting times. We created such a distribution as a mixture of two uniform densities, which we call the UMM2 (UMM for uniform mixture model). We define the UMM2 by otherwise. (B2) Here we use UMM2 rather than EMM2 to suppress the effect of the waiting-time paradox; unless is large, τ 2 is smaller for a UMM2 than an EMM2 given the same mean τ . Then, we imposed τ 2 = 2, such that b = 0 and CV = 1 as in the exponential distribution, which gives = 2 3(1 − a) − 3a(2 − a) /(3 − 2a). In this way, one can compare the UMM2 and the exponential distribution that share the mean inter-event time τ = 1 and the mean waiting time τ 2 /(2 τ ) = 1, while the two distributions have different probability masses at small waiting times. For the exponential distribution, one obtains ψ(τ ) = ψ w (τ ) = exp(−τ ). When ψ(τ ) is the UMM2, one obtains Because ψ w (0) = 1 for both exponential distribution and UMM2, by letting large, one can create a relatively large probability mass of the waiting-time for the UMM2 at small τ values as compared to the exponential distribution ( Fig. 3(b)).
We set a = 0.1, which yields ≈ 1.389. The final size for this UMM2 is shown by the red line in Fig. 3(a). The figure suggests that the final size is larger for this distribution as compared to the exponential distribution. We note that the two distributions have the same τ and τ 2 . variant we consider, an infection happens with probability β if and only if v i is susceptible and v j is infected. In the second variant, infection happens with probability β if and only if v i is infected and v j is susceptible.
These node-centric SIR models are partially affected by the waiting-time paradox in different manners. In the first variant, if v i is infected by v j , the time until v i contacts its neighbor v for possibly infecting v (if v is susceptible) obeys the waiting-time distribution. Therefore, possible infection of v by v i tends to be delayed, such that epidemic spreading may be suppressed by a fat-tailed ψ(τ ). In contrast, epidemic spreading may be enhanced in this model because, if v i does not get infected by v j , then v i may draw a short inter-event time according to ψ(τ ) to be infected by its neighbor. In the second variant, if v i infects v j , then the time to the next activation of v j , which may let v j to infect another node, obeys the waiting-time distribution. Therefore, this secondary infection process may happen late. In contrast, epidemic spreading in the same model may be enhanced because, after infecting v j , node v i draws the time to its next activation according to ψ(τ ) to possibly infect its different neighbor.
Consider the first variant of the node-centric SIR model. When the population is infinite and well-mixed and ψ(τ ) is an EMM2, the stochastic dynamics of epidemic spreading are given by Similarly, the stochastic dynamics of the second variant of the node-centric SIR model are given by The eigenvalues of the Jacobians of these systems evaluated at the disease-free equilibrium determine the epidemic threshold. For both SIR models, the eigenequation to determine the two eigenvalues resulting from the leading principal minor of size two is the same that of the original SIR model. In addition, for both variants, the other two eigenvalues of the Jacobian are given by −µ + Λ 1 and −µ + Λ 2 , where Λ 1 and Λ 2 are the solutions of The difference in the multiplicative factor 2, which is present in the two terms proportional to β in Eq. (A3) and absent in Eq. (C3), originates from the fact that the original model, but not the two variants, allows bidirectional infection between an activated node and its neighbor. Therefore, this difference is not essential. The crucial difference between the original model and the two variants are the last term in Eq. (A3). This term originates from the simultaneous presence of the two infection processes on the same edge, where an activated infected node can infect a susceptible neighbor and an activated susceptible node can be infected by an infected neighbor. By substituting Eqs. (3b), (3d), (4a), and (4b) into Eq. (C3), one obtains Λ = β/ τ , −1/ τ . Therefore, one obtains µ c = β/ τ , or λ eff ≡ β/ τ µ c = 1, which implies that the distribution of inter-event times does not affect the epidemic threshold.
In fact, numerically obtained final sizes are smaller for PL1 or EMM2 than for the exponential ψ(τ ) with the same mean on the complete graph and different networks (Figs. 4 and 5).