Burstiness in activity-driven networks and the epidemic threshold
Marco Mancastroppa Alessandro Vezzani Miguel A. Muñoz Raffaella Burioni





^1 Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università degli Studi di Parma, Parco Area delle Scienze, 7/A 43124 Parma, Italy
^2 INFN, Gruppo Collegato di Parma, Parco Area delle Scienze 7/A, 43124 Parma, Italy
^3 IMEM-CNR, Parco Area delle Scienze, 37/A, 43124 Parma, Italy
^4 Departamento de Electromagnetismo y Física de la Materia e Instituto Carlos I
  de Física Teórica y Computacional. Universidad de Granada.
  E-18071, Granada, Spain
marco.mancastroppa@unipr.it, alessandro.vezzani@unipr.it,  mamunoz@onsager.ugr.es and raffaella.burioni@unipr.it



We study the effect of heterogeneous temporal activations on epidemic spreading in temporal networks. We focus on the susceptible-infected-susceptible (SIS) model on activity-driven networks with burstiness.
By using an activity-based mean-field approach, we derive a closed analytical form for the epidemic threshold for arbitrary activity and inter-event time distributions. We show that, as expected, burstiness lowers the epidemic threshold while its effect on prevalence is twofold.  In low-infective systems burstiness raises the average infection probability, while it weakens epidemic spreading for high infectivity. Our results can help clarify the conflicting effects of burstiness reported in the literature. We also discuss the scaling properties at the transition, showing that they are not affected by burstiness. 




Keywords: Network dynamics, Epidemic modelling, Phase transitions into absorbing states, Critical phenomena of socio-economic systems.





§ INTRODUCTION

The heterogeneous distribution of times between two consecutive actions – often denoted as burstiness – is a typical signature of time-resolved records of human activities <cit.>.  Mobile phone calls, social contacts and many other complex real-world dynamics exhibit patterns of enhanced activity within short time periods, followed by long lags of inactivity.  Actually, the frequency of events shows large variability and strong temporal fluctuations <cit.> and the event dynamics strongly differs from a Poisson process, where the rate at which events occur is constant. In “bursty” processes, the probability distribution for the inter-event time between two successive actions is not an exponential as in Poisson processes but, typically, it features a fat tail, easily measured from large datasets <cit.>.

Recent studies have addressed the problem of the origin of this highly heterogeneous behavior, arising from the decision-based queuing process that humans apply in performing tasks and allocating their priority <cit.>. Besides, a large amount of work has been devoted to clarify the effects of intermittent patterns on the temporal structures of interactions, as described by temporal networks <cit.>. Bursty dynamics can indeed influence the structure of links on the local and on the global scale <cit.>. More importantly, heterogenous temporal patterns in the evolution of time-varying networks can affect in a non-trivial way dynamical processes. Analytical arguments <cit.> and numerical approaches <cit.> have indeed shown that burstiness can significantly modify processes mediated by interactions, such as random walks, epidemics, information diffusion, consensus formation, percolation. Among these, epidemic spreading is one of the most representative and widely applicable example <cit.>.  Analytical results on epidemics in temporal network have focused on the epidemic threshold <cit.>, in general pushed to lower values by bursty effects. In other works, burstiness is introduced on a static network as a non-Markovian effect in infection (or recovery) processes <cit.>.

Interestingly, when comparing bursty with Poisson dynamics in epidemics, 
conflicting observations have been reported. On the one hand,
numerical results <cit.> and modeling techniques <cit.> provided strong evidence of a slowing down for the late time spreading in the presence of burstiness, while opposite  effects are reported in the early-time dynamics <cit.>. In <cit.> and <cit.>, burstiness is found to slow down spreading in early times, while other works observe the opposite effect <cit.>.

To understand such conflicting observations and keep track at the same time of the temporal evolution of interactions, in this paper we focus on the Susceptible-Infected-Susceptible (SIS) process in the presence of burstiness on activity-driven networks <cit.>. In activity-driven networks,  the propensity to engage
an interaction is moved from the links to nodes (N), by assigning to each node i its activity potential a_i measuring the typical number of activations (link formations) per time performed by agent i. We then model the bursty behavior by a fat tailed distribution of inter-event times between two successive activations of the same node. Using an activity based mean-field approach, which is here exact in the large N limit, we find a closed analytical form for the epidemic threshold for general activity and inter-event time distributions including, as a particular case, bursty dynamics. The analytical results are in excellent agreement with extensive numerical simulations. We show that, as expected, burstiness lowers the epidemic threshold, while its effect on prevalence is twofold: burstiness tends to raise the average stationary infection probability  in low infective systems, therefore strengthening the epidemic, while it weakens the epidemic in high infective systems, lowering the prevalence. This result can help to clarify the conflicting effects of burstiness reported in the literature. We also analyze in details the scaling behavior at the transition, showing that the critical exponents are mean-field exponents of the directed percolation (DP) universality class and that they are unchanged in the presence of burstiness.

The paper is organized as follows. In Section <ref> we describe the activity-driven network with bursty dynamics and we define the epidemic process unfolding on top of it. In Section <ref> we present the analytical approach and we derive an exact expression for the epidemic threshold, comparing the result with extensive numerical simulations. In Section <ref> we analyze the effects of burstiness on disease spreading, focusing on the epidemic threshold as well as on prevalence, and we discuss the scaling properties. Finally, Section <ref> contains a summary of our results and some perspectives of our work.



§ EPIDEMICS ON ACTIVITY-DRIVEN NETWORKS



 §.§ Activity-driven temporal networks


We introduce the activity-driven network by associating to each node i=1,…,N an activity a_i drawn from a distribution ρ(a_i). The activation rate of node i is defined by its inter-event time distribution Ψ_a_i(τ) <cit.>: this sets the statistic of  time intervals between consecutive activations of node i, so that the average inter-event time equals the inverse of the activity of node i, i.e. ⟨τ⟩_i=a_i^-1. The network evolves  according to a Gillespie-like algorithm <cit.>. We initialize the system at t=0, drawing from Ψ_a_i(t) the first activation time, t_i, of each node. Then the agent i with the lowest t_i becomes active and create m links connecting to 
m randomly chosen nodes. We fix the next activation time t_i for node i as t_i + τ
with τ drawn from Ψ_a_i(τ). Finally all  links are removed and the process is iterated, activating the node with smallest activation time t_i.

Empirical measures on real datasets provide a strong evidence of heterogeneous activity patterns, with ρ(a_i) being a broad, heterogeneous function with large fluctuations <cit.>. We model such a heterogeneity in the activity distribution as a power law: 

    ρ(a_i) = ν  a_m^ν a_i^-(ν+1)

for a_i>a_m and ρ(a_i)=0 for a_i<a_m, where a_m is a lower cut-off.  However, as we show in the next sections, our results hold for any ρ(a_i) and therefore for any activity-driven network.

The inter-event time distribution sets the dynamics of activation for each site. For example,
an exponential Ψ_a_i(τ) =a_i e^-a_i τ leads to a Markovian Poisson process as in standard continuous time activation dynamics.
When the inter-event time distribution of a single node is heterogeneous, with trains of activation events with short time separation, followed by long inactivity, 
we take into account the bursty behavior choosing a power-law (Pareto) distribution of inter-event times with a cut-off ξ_i: 

    Ψ_a_i(τ) = αξ_i^ατ^-(α+1)

for τ> ξ_i and Ψ_a_i(τ) =0 for τ<ξ_i.  More specifically; for α>1, we define ξ_i = α-1/α a_i so that ⟨τ⟩_i=a_i^-1, and the distribution of ξ_i is Φ(ξ_i) ∼ρ(α-1/αξ_i) ξ_i^-2. Doing so, we can effectively compare systems with the same average inter-event time, fixed by the inverse of the node activity, but with different burstiness distributions Ψ_a_i(τ).  In the case 0<α<1, ⟨τ⟩_i diverges and the activity is not properly defined.  We can nevertheless introduce the inter-event time distribution according to Eq. (<ref>). As in <cit.> we set ξ_i=a_i^-1, with the cut off ξ_i distributed among the nodes according to Φ(ξ_i) = ξ_i^-2ρ(ξ_i^-1). However, in this case a_i does not correspond to the activity of the node, i.e. the inverse of average activation time, and non-trivial aging effects may be present <cit.>. In particular, for the SIS model, since node activation becomes always more unlikely as time increases, we will show that epidemics do not spread in the system.






 §.§ The SIS epidemic model on activity-driven networks


We consider the epidemic Susceptible-Infected-Susceptible (SIS) model. Each node can be infected (I) or susceptible (S-healthy) and
the dynamics is based on two elementary processes: contagion, when a susceptible individual, coming in contact with an infected person, has a probability λ to contract the infection S+I  2I; spontaneous recovering, when an infected turns susceptible with a rate μ: I  S.  The SIS model on activity-driven networks is defined by a Gillespie-like algorithm as follows <cit.>):  

  * Up to a relaxation time, t_0, the network evolves in the absorbing state where all nodes are susceptible, so that the node activation dynamics relaxes to its equilibrium.  

  * At time t = t_0 we start the epidemics setting the initial condition, for which we divide the population into a configuration of susceptible (S) and infected (I).    

  *  Agent i with the lowest activation time t_i activates. Infected nodes at time t are recovered at t_i with probability 1-e^-μ (t_i-t), i.e. the recovery is a Poisson process. Then we set t=t_i.  

  * The active node i generates m links with other randomly-chosen nodes. If the nodes involved in the contacts (links) are both infected (I) or susceptible (S) nothing happens, otherwise a contagion can occur with probability λ: S+I  2I. Here we consider m = 1, with no loss of generality. 

  * A new inter-event time for the active node i is drawn from Ψ_a_i(τ), hence the new activation time t_i would be t_i + τ. Then all its links are removed and the process is iterated from point <ref>.  


As well known, the SIS model features a phase transition in the class of directed percolation (DP) between an absorbing phase, where the system evolves towards a state with the whole population susceptible (healthy), and an endemic phase, with  the epidemics spreading persistently in the population <cit.>. We consider a control parameter r  (e.g. r=λ/μ on static networks, r=λ⟨ a ⟩/μ on activity-driven networks with Poisson dynamics <cit.>), so that the phase transition occurs at a critical value  r_C, the epidemic threshold: if r<r_C the system is in the absorbing phase, while if r> r_C the system is in the endemic phase. This is the quantity that we will  investigate in the next section.



§ THE EPIDEMIC THRESHOLD



 §.§ Analytical results

We consider the SIS model on activity-driven networks with arbitrary ρ(a) and Ψ_a(τ).  If the average inter-event time ⟨τ⟩ is finite (i.e. α>1 considering the power-law distribution of Eq. (<ref>)), nodes can be identified according to their activity, so we divide the population into classes with the same a_i=a and study the probability that a node belonging to a class a is infected at time t. This corresponds to an activity-based mean-field approach <cit.>, which turns out to be exact for the present system given that connections are continuously reshuffled at each time step, destroying local correlations.


At each contact event occurring at time t, we consider the following three quantities to fully characterize the system state:

- Q_a(t): the probability  for a node in the class a, that activates at time t, to be infected right before t;

-  P_a(t): the probability for a non-active node, with activity a, to be infected at t;

- Z_a(t,t'): the probability for a node with activity a to be infected at time t, knowing that its last activation occurred at time t'<t.

The Master Equation for the evolution of Z_a(t,t') can be written as:

    ∂_t Z_a(t,t')=-μ Z_a(t,t') + λ (1-Z_a(t,t')) ∫_0^∞ da' ρ(a') a' Q_a'(t)


where first term accounts for the recovery process with rate μ, while the second term describes (for α>1) the contagion process; actually, this is the product of the probability to be susceptible at time t, (1-Z_a(t,t')), times the probability to be infected by an active node, which, for α>1, can be expressed as
λ∫_0^∞ da' ρ(a') a' Q_a'(t) [ The probability for node i to infect a node in the interval  [t,t+dt] is λ n_i(t)Q_i(t)/N, where n_i(t) is the number of activations of node i in dt. For α>1  we have n_i(t)= dt/⟨τ⟩_i =a_i dt, where ⟨τ⟩_i is the mean time between two activations.  Thus, the probability to be contacted by any infected node in dt is N^-1∑_i=1^N a_i   dt   Q_i(t) = dt ∫ da' ρ(a') a' Q_a'(t).]  Notice that this expression –together with the whole analytic approach– does not apply for 0<α<1, since the average inter-event time diverges and a_i dt is not the probability of activation of node i in the interval dt. We observe that in Eq. (<ref>) we only include passive contagion events because, by definition of Z_a(t,t'), the node does not activate between times t' and t.  

Eq. (<ref>) needs to be solved using the initial condition at time t':

    Z_a(t',t')=Q_a(t')+λ [1-Q_a(t')] ∫_0^∞ da' ρ(a') P_a'(t')

where, in analogy with the equation above, the first term corresponds to the case where the active node with  activity a is already infected at the contact time t'; the second term is the probability (1-Q_a(t')) that the node is susceptible before t', then it connects to an infected node with probability  ∫_0^∞ da' ρ(a') P_a'(t'), and it is infected with probability λ.


Eq. (<ref>) with initial condition (<ref>) can be solved in the stationary condition, in the hypothesis that the system reaches an asymptotic steady state (P_a^0,Q_a^0)  with 
(P_a(t),Q_a(t))  (P_a^0,Q_a^0). In this regime, also the node activation dynamics needs to reach the steady state to avoid aging effect in the network evolution. For α> 1 we expect this condition to be fulfilled if the relaxation time is much larger than the average activation time  (t_0 ≫⟨ a ⟩^-1).

In the stationary regime, one can consider the averages on the activities as independent of time
∫_0^∞ da ρ(a) a Q_a(t) ∼∫_0^∞ da ρ(a) a Q_a^0 = aQ and
∫_0^∞  da  ρ(a)   P_a(t)∼∫_0^∞ da ρ(a) P_a^0 = P
to obtain:

    Z_a(t,t')= [ Q_a(t') + λP (1-Q_a(t')) - λaQ/μ+λaQ] e^-(μ+ λaQ)(t-t') + λaQ/μ + λaQ


It is now possible to write explicitly Q_a(t) and P_a(t) as functions of Z_a(t, t')
using the relations: 
    Q_a(t)    =   ∫_0^∞ Z_a(t,t-τ') Ψ_a(τ') dτ' 
     P_a(t)    =   1/∫_0^∞ dτ' ∫_τ'^∞Ψ_a(τ”) dτ”∫_0^∞ dτ' Z_a(t,t-τ') ∫_τ'^∞ dτ”Ψ_a(τ”),

and plugging the solution (<ref>) into Eqs. (<ref>), one obtains:

    P_a(t)   =   λaQ/μ + λaQ + a [ λP - λaQ/μ + λaQ] 1-L_a(μ+ λaQ)/μ + λaQ
       +    a [1-λP] ∫_0^∞ dτ' Q_a(t-τ') e^-(μ+λaQ)τ'∫_τ'^∞ dτ”Ψ_a(τ”) 
    
    Q_a(t)   =   λaQ/μ + λaQ+ [ λP - λaQ/μ + λaQ] L_a(μ+ λaQ) 
       +    [1 - λP] ∫_0^∞ dτ'  Q_a(t-τ') e^-(μ+λaQ)τ'Ψ_a(τ')

where L_a(μ)=∫_0^∞ dτ  e^-μτΨ_a(τ) is the Laplace transform of Ψ_a(τ). In order to determine the epidemic threshold, one needs to study the stability of the absorbing state solution (P_a(t),Q_a(t))=(0,0)  ∀  a of the Eqs. (<ref>).


Solving Eqs. (<ref>) for t →∞, two self-consistency equations for P, aQ are obtained:

    P   = h(P,aQ)= λaQ/μ + λaQ + ∫_0^∞ da ρ(a)   a 1-L_a(μ+ λaQ)/(μ+λaQ)^2λμP/1-(1-λP)L_a(μ+ λaQ)
    aQ   = f(P,aQ)= ∫_0^∞ da  ρ(a)   a [ λaQ/μ + λaQ + ( λP - λaQ/μ + λaQ) L_a(μ+ λaQ) ]/1-(1-λP)L_a(μ+ λaQ)

 (see <ref> for a detailed derivation).
Eqs. (<ref>) can be considered as a two dimensional map: (P_n+1, aQ_n+1)= M(P_n, aQ_n). The fixed point  (P_n,aQ_n)=(0,0) corresponds to the solution (P_a(t),Q_a(t))=(0,0)  ∀  a in Eq. (<ref>). Then, the stability of the absorbing state can be studied by linearizing  M(P_n, aQ_n) in (P_n,aQ_n)=(0,0):

    [  P_n+1; aQ_n+1 ]∼
    J
    [  P_n; aQ_n ]
    =
    [                         λ⟨ a ⟩/μ                              λ/μ; λ∫_0^∞ da ρ(a) a L_a(μ)/1-L_a(μ)                         λ⟨ a ⟩/μ ][  P_n; aQ_n ]

The absorbing state is stable if all the eigenvalues of J are smaller than unity, i.e.:

    λ⟨ a ⟩/μ + λ√(1/μ∫_0^∞ da ρ(a) a L_a(μ)/1-L_a(μ))<1.



Eq. (<ref>) provides an analytical form of the epidemic threshold for the SIS model on an activity-driven network (m = 1) with arbitrary ρ(a) and Ψ_a(τ). As we will show, numerical simulations agree with the analytical estimate, suggesting that the activity based mean-field approach and the stationary hypothesis are exact for this model.

In particular, for an exponential distribution of inter-event times Ψ_a(τ) = a e^-a τ and an arbitrary activity distribution ρ(a), one can get the epidemic threshold from Eq.(<ref>): 

    λ⟨ a ⟩/μ = ⟨ a ⟩/⟨ a ⟩ + √(⟨ a^2 ⟩)

which is the known result for the SIS model on activity-driven networks with standard non-bursty dynamics <cit.>.
 
In the general case, the parameters λ, μ and ⟨ a ⟩ are coupled in a non-trivial way, since  Eq. (<ref>)
depends on the whole functional form of the activity and inter-event time distributions, so that the control parameter of the phase transition is not easily identified.
Hereafter, we will use the adimensional control parameter r= ⟨ a ⟩/μ. Notice that, unlike the exponential case, the critical value depends in a non-trivial way on λ.

For a bursty distribution of inter-event times Ψ_a(τ) ∼τ^-(α+1)  and the power law activity distribution ρ(a) ∼ a^-(ν+1) defined by Eqs. (<ref>) and (<ref>) respectively, Eq. (<ref>) becomes:

    λ r + λ√(ν ( r(ν-1)/ν)^ν∫_r (ν-1)/ν^∞ dx 1/x^ν α( α-1/α x)^αΓ(-α, α -1/α x)/1-α( α-1/α x)^αΓ( -α, α -1/α x)) <1

where Γ(s,x)=∫_x^∞ t^s-1e^-t dt is the upper incomplete Gamma function. In this case, the epidemic threshold, r_C, can be obtained by solving numerically Eq.(<ref>) for r.



 §.§ Numerical simulations

We verify with numerical simulations the analytic result (<ref>), checking the validity of our hypotheses. Henceforth, we fix λ=1, with no loss of generality.  We first let the network evolve for a time t_0 so that network dynamics is relaxed to the equilibrium of the activation dynamics and it has no memory of the initial condition at t=0 (when all sites are considered to be active). For α>1 this means that t_0 needs to be much larger than the average activation time t_0 ≫⟨ a ⟩^-1.

The epidemic model evolves following the steps described in Section <ref>. The epidemic threshold is numerically computed with the lifetime (or lifespan) method <cit.>. The lifetime τ_life is the average time that the system takes to reach the absorbing state when infecting at t_0 a single node of the system, discarding in the average  all the endemic realizations which remain infected for an infinite time. In simulations, we consider as endemic a realization whose coverage (number of distinct, ever infected nodes) reaches N/2. Observe that fluctuations can cause finite-time recovering even in the endemic phase so that the lifetime as a function of r displays a peak at the critical point r=r_C behaving like an effective susceptibility as discussed in <cit.>. Also, it has been shown that the most effective initial condition to obtain a good estimate of the threshold on an inhomogeneous network is to infect at t_0 only the most connected site<cit.>; this is the strategy we adopt here, infecting the most active one.

Figure <ref> shows the lifetime as a function of the control parameter for different network sizes N, setting all nodes with same activity (i.e. ρ (a) = δ (a-a_m)).
Left and right panels consider an exponential Ψ_a(τ) ∼ e^-a τ and a power-law Ψ_a(τ) ∼τ^-(α+1) distribution, respectively. In both cases, for sufficiently large N, the numerical data agree very well with the analytical prediction of the threshold.

Figure <ref> shows the lifetime as a function of r for a network where ρ(a) ∼ a^-(ν+1) and Ψ_a(τ) ∼τ^-(α+1) are power-law distributed according to Eqs. (<ref>) and (<ref>) respectively. Several values of the exponent α have been considered. When 0< α<1, as we discussed above, analytic results do not apply and we can introduce the inter-event time distribution Ψ_ξ_i(τ) ∼ξ_i^ατ^-(α+1) with a broadly distributed  cut-off ξ_i (as stated in Section <ref>). Simulations reveal that the lifetime does not feature a peak since the epidemics always reaches the absorbing state. In particular, the average recovering time μ^-1 is finite while node activation becomes always more unlikely as time increases, implying that, eventually, all nodes recover: an endemic or active phase does not exist in this case.

Finally, in Figure <ref> we compare numerically and analytically obtained epidemic thresholds in the case of bursty dynamics for two different activity distributions ρ(a), illustrating the excellent agreement between theory and numerics. 




§ BURSTINESS AND EPIDEMIC SPREADING



 §.§ Epidemic thresholds for various distributions


As explained above, for a generic distribution of inter-event times and activities, we use r= ⟨ a ⟩/μ as a control parameter for the phase transition, setting λ=1.  In particular, for a system with Ψ_a(τ) and ρ(a) given by Eqs. (<ref>) and (<ref>) respectively, the critical threshold r_C is obtained solving Eq. (<ref>); results are shown in Figure <ref>.  The behavior of r_C is qualitatively the same for all ν: it vanishes for α→ 1^+ and it saturates to a constant depending on ν for α→∞.  For α<1, as explained above, the epidemics always ends up in the absorbing state.


As expected, a broader activity distribution leads to a lower epidemic threshold and favors the propagation of the epidemics. In particular, for large ν the activity distribution tends to a δ-function while in the limit ν→ 2, i.e. when the activity fluctuations diverge, the epidemic threshold vanishes.

Let us now show that this is a general property for any inter-event distribution. Since ∫_0^∞τΨ_a(τ) dτ=a^-1, Ψ_a(τ) tends to a δ-function peaked in τ=0 when a→∞. Hence, it is possible  to expand e^-μτ around τ=0 when evaluating L_a(μ)  for a →∞, in this way: L_a(μ)  1
and 1-L_a(μ)=∫_0^∞ (1-e^-μτ)Ψ_a(τ) dτ∫_0^∞μτΨ_a(τ) dτ= μ/a. Plugging these expression into Eq. (<ref>) we obtain that the integral over the activities diverges for ν<2 due to the slow decay at large a.  Hence, Eq. (<ref>) cannot be satisfied 
for ν<2 and the absorbing state is never stable. 
The behavior for ν→ 2 and α→ 1  confirms that, keeping the average activity fixed, fluctuations and burstiness lower the epidemic threshold. 

In Figure <ref>, we show the epidemic threshold for a power-law activity distribution Eq. (<ref>) with ν=2.5 and compare several inter-event time distributions Ψ_a(τ). In the left panel the continuous (blue) line is the epidemic threshold r_C for inter-event times distributed according to the power law distribution with a cut-off in Eq. (<ref>), Pareto type I distribution. The plot shows that for large α the threshold becomes even higher than in the exponential case, i.e. the dotted (orange) horizontal line. Indeed, for α→∞, the power law distribution Ψ_a(τ) tends to a δ-function and the node dynamics is periodic <cit.>, hence more regular than the Poisson dynamics (see the green dash-dotted line for comparison in Figure <ref> - left panel). The dashed (red) curve in the left panel represents the epidemic threshold r_C for the Lomax (or Pareto type II) distribution Ψ_a(τ)= αξ_a^α (τ+ξ_a)^-(α+1) with τ>0: here ξ_a=(α-1)a^-1 so that ⟨τ⟩_a=a^-1, i.e. we introduce a different lower cut-off.  In this case, for α→∞ we have that Ψ_a(τ)→ a e^-aτ, so in this limit the epidemic threshold tends to the value of the Poisson process (dashed orange line). The Lomax distribution is more heterogeneous than an exponential distribution, therefore the epidemic threshold is smaller for all α. 

Another interesting choice for the inter-event time –usually
considered in epidemic spreading–
 is the Weibull distribution Ψ_a(τ)= βδ_a^-βτ^β-1 e^-(τ/δ_a)^β with τ>0 <cit.>: here δ_a = (a Γ(1+1/β))^-1 so that ⟨τ⟩_a = a^-1. The Weibull distribution  interpolates between a broad distribution for β→ 0^+ (it tends to the Zipf distribution <cit.>) and a δ-function for β→∞. Moreover, for β=1 it reduces to an exponential distribution. Notice that for any β>0 for large τ it decays faster than any power-law.
In the right panel of Figure <ref> we plot r_C for the Weibull distribution with a solid (red) line.
The threshold vanishes when β→ 0, it increases with β
and for β=1, r_C coincides with that of the exponential case (dotted orange line).
For β>1 the threshold increases and for β→∞ it converges to the threshold of a periodic dynamics (dash-dotted green line).

 In general, from Figure <ref>, we see that the value of r_C depends on the whole shape of the inter-event time distribution, displaying, however, some general features.  For broad distributions, i.e. small α and β, the threshold is small and it vanishes for α→ 1 and β→ 0.  Then, the threshold increases with α and β and it reaches a maximum value for a periodic dynamics i.e. if Ψ_a(τ) tends to a δ-function.  From here on, we consider a power-law distribution of inter-event times Ψ_a(τ) as defined in Eq. (<ref>).






 §.§ Effects of burstiness on prevalence
 
Within our analytic approach, beside the epidemic threshold value, a numerical solution of Eq. (<ref>) provides us with the value of P in the stationary state.
 This probability, called prevalence, is the order parameter of the phase transition: it is zero below the epidemic threshold (r<r_C) and is different from zero above it (r>r_C).
In Figure <ref> we plot the prevalence as a function of the control parameter r = ⟨ a ⟩/μ  for different values of the burstiness exponent α, when ρ(a) ∼δ(a-a_m) (left panel) and when ρ(a) ∼ a^-(ν+1) as in Eq. (<ref>) (right panel).  


Plots confirm that burstiness lowers the epidemic threshold as discussed above. However, they also illustrate that the global effect of  burstiness (α→ 1) is twofold: it raises the prevalence in low infective systems (r <1) strengthening the epidemic, while it weakens the epidemic in high infective systems (r > 1), lowering the prevalence. This behavior helps clarifying the conflicting effects of burstiness reported in the literature.  At λ=1 (see Fig. <ref>), the switch from a regime where the burstiness strengthens the epidemics to the regime where epidemics is weakened occurs at r=1 and P=1/2: in <ref> we show that, regardless the shape of Ψ_a(τ) and ρ(a), this property holds when r=1/λ. One can argue that the twofold effect of burstiness on prevalence is purely due to fluctuations. In this perspective, in <ref> we show that the same behavior can be induced by activity fluctuation in systems without burstiness.



 §.§ Critical scaling behavior


The study of so-called spreading experiments – in which a single seed of infection is placed into an otherwise quiescent network – is a standard tool to locate the phase transition between the active and the absorbing state. In such experiments one usually computes,  after the initial seed is placed, the probability P_d(t) that a run reaches the absorbing state in a time interval [t, t+dt], and the total number of infected nodes N(t) at a given time t. These quantities scale as power laws at criticality, with exponents δ and η, respectively, defined by P_d (t) ∼ t^-δ-1 and N(t) ∼ t^η <cit.>. In particular, in the DP universality class the mean-field exponents take the values δ=1 and η=0 <cit.>.  These are typically observed in high dimensional lattices, dense networks and fully connected systems.  In Figure <ref>, for ρ(a) ∼δ(a-a_m), we verify that the exponent δ in our activity driven model model is consistent with the mean field value δ≃ 1 both for exponential and power-law inter-event time distributions. A similar result can also be obtained for the exponent η≃ 0. 




In finite size systems, the probability of reaching the absorbing state is:

    P_d(t,L) ∼ t^-δ-1F(t L^-z)

where F(t L^-z) is an unspecified scaling function, controlling the finite-size effects,
 z is the dynamical exponent and L is –in general– a measure of the linear size of the system <cit.>.  The case of a lattice of size N we get N∼ L^d, where d is the dimension. Thus, in general, the average lifetime is τ_life=⟨ t ⟩∼∫ dt   t^-δ -1  t  F(t N^-z/d) ∼N^z(1-δ)/d.  Since in the mean-field case δ=1,  the integral results into a logarithm and one expects τ_life to grow logarithmically with network size τ_life∼log(N). 
Similarly, as τ_life is peaked around a central 
value, the position of the peak, τ_life^max, scales in the same way.

As shown in Figure <ref> this is, indeed, the case in numerical simulations for either exponential (Poisson) or power-law distributed inter-event time distributions keeping ρ(a) ∼δ(a-a_m) (see figure caption).


In the light of the above, let us adopt a finite-size scaling approach – similar to that proposed in <cit.> – to analyze lifespans (such as those shown in Figure <ref>) as a function of N and the distance Δ to the critical point: Δ≡(r-r_c)/r_c.  Since right at the critical point (Δ=0) lifespans grows as ln(N), slightly away from the epidemic threshold we get the finite-size scaling relation:

    τ_life(Δ,N) ∼ln(N) G(Δ N^γ)
  
where γ is expected to take its mean-field value 1/2 [Actually, the characteristic length-scale scales in mean field as L ∼Δ^-1/2 and N ∼ L^d, where in mean-field d=d_C=4 is the critical dimension, so that N∼Δ^-2 and, thus γ=1/2 <cit.>.]  Finally, in Figure <ref> we verify the scaling hypothesis of Eq. (<ref>) both in the exponential and in the power law case, showing a rather good curve collapse when rescaling (dividing) τ_life with ln(N) and Δ with a factor N^γ=1/2.  



We also measure the  critical exponent β, that describes the vanishing of the prevalence P at the epidemic threshold:

    P∼ (r-r_C)^β

for r>r_C. In <ref> we show that β = 1 for the Poisson process, regardless of ρ(a). Figure <ref> shows that also in the bursty case β≃ 1 for every α and activity distribution ρ(a).


Summarizing, despite the non-trivial temporal dynamics, critical exponents suggest that burstiness does not change the universality class of the absorbing-active phase transition for the SIS model on activity-driven network, when the average inter-event time is finite (α > 1). The system presents the same mean-field like scaling properties for non-Poissonian (bursty) and Poissonian dynamics, regardless of the exponent α, and shows the mean-field exponents of DP universality class. A similar equivalence is obtained in <cit.>, where burstiness is introduced on static networks as a non-Poissonian effect in infection process.



§ CONCLUSIONS


Our results provide a coherent understanding of the interplay between the bursty dynamics of an activity-driven network and the epidemic process unfolding on it.
Focusing on the SIS model, by means of an activity-based mean-field approach we analytically derive a closed form for the epidemic threshold on an activity-driven network: the threshold depends on the whole functional form of the activity and inter-event time distributions. The mean-field approach turns out to be exact in this system, as connections are continuously reshuffled destroying local correlations. We show that the analytical results are in excellent agreement with extensive numerical simulations.

The model captures two key aspects driving the evolution of social
networks: burstiness and heterogeneity in human activities. As is well known, heterogeneity produces a decrease in the epidemic threshold.
We show that burstiness also reduces the epidemic threshold, while the effect on prevalence is more subtle: in low-infective systems burstiness raises the prevalence, while it weakens the epidemics in high infective systems. We showed that this behavior can be attributed to fluctuations effects. Finally, despite the non-trivial temporal dynamics,  the values of the scaling exponents suggest that burstiness does not change the universality class of the absorbing transition.

Our results can help to clarify the effect of burstiness on dynamical process in temporal network and, in turn, to design efficient immunization strategies <cit.> and control protocols, based on temporal structure of the human interactions. The model can be extended to different epidemic processes and  the analytic results for the thresholds should hold also in the case of the SIR model,  while a different behavior is expected in the active phase <cit.>. We can also include in the model further features of realistic social networks, i.e. the presence of communities <cit.>, adaptive behavior <cit.> and memory <cit.>, these combined with burstiness can produce non-trivial effects <cit.>. 





§ EVALUATION OF THE EPIDEMIC THRESHOLD

Let us consider the system of equations (<ref>).
We set an initial condition, adding ε_a δ(t) to Q_a(t): i.e. we impose that at t=0 all nodes are active and a node with activity a is infected with probability ε_a.  For  Q_a(s) and P_a(s), the Laplace transforms of Q_a(t) and P_a(t), we obtain:


    Q_a(s) = ∫_0^∞ dt   e^-st  Q_a(t) =   1/s[ λaQ/μ + λaQ + ( λP - λaQ/μ + λaQ) L_a(μ+λaQ) ]
       + (1-λP) Q_a(s) ∫_0^∞ dτ' e^-(μ+λaQ+s)τ'Ψ_a(τ') + ε_a


    P_a(s) = ∫_0^∞ dt   e^-st  P_a(t) =   1/s[ λaQ/μ + λaQ + a ( λP- λaQ/μ + λaQ) 1-L_a(μ+ λaQ)/μ + λaQ]
       + a (1-λP) Q_a(s) ∫_0^∞ dτ' e^-(μ + λaQ +s )τ'∫_τ'^∞ dτ”Ψ_a(τ”)

We substitute Q_a(s) in P_a(s) and expand both for s ∼ 0, since for the calculation of the epidemic threshold we only need the asymptotic behavior (t →∞). We obtain at first order in 𝒪(1/s):


    Q_a(s) = C_1/s + 𝒪(1/s)


    P_a(s) = C_3/s + C_1/s a (1-λP)  1-L_a(μ+ λaQ)/μ + λaQ + 𝒪(1/s)


where

    C_1= [ λaQ/μ + λaQ + ( λP - λaQ/μ + λaQ)L_a(μ+ λaQ) ]/1-(1-λP)L_a(μ+ λaQ)    ,      L_a(s)= ∫_0^∞ dτ  e^-s τ Ψ_a(τ)



    C_3 = λaQ/μ+λaQ + a ( λP - λaQ/μ+λaQ) 1-L_a(μ+λaQ)/μ + λaQ

Going back to the temporal domain, in the asymptotic limit we obtain:

    Q_a(t)  C_1=Q_a^0     ,     
    P_a(t)  C_4=P_a^0

where


    C_4=C_3 + C_1 a(1-λP) 1-L_a(μ+ λaQ)/μ+λaQ = C_3 + a(1-λP) 1-L_a(μ+ λaQ)/μ+λaQ Q_a^0

In this way we obtain two self-consistency equations for Q_a^0 and P_a^0:

    P_a^0=λaQ/μ + λaQ + a 1-L_a(μ+ λaQ)/(μ + λaQ)^2λμP/1-(1-λP)L_a(μ+ λaQ)


    Q_a^0=[ λaQ/μ + λaQ + ( λP - λaQ/μ + λaQ)L_a(μ+ λaQ) ]/1-(1-λP)L_a(μ+ λaQ)

since P = ∫_0^∞ da  ρ(a)   P_a^0 and aQ= ∫_0^∞ da   a  ρ(a)   Q_a^0. Integrating the above expression over a we get the self consistence equations described in Eq. (<ref>).



§ PREVALENCE AT R=1/Λ


We can write Eqs. (<ref>) and (<ref>) for P_a^0 and Q_a^0 as:

    P_a^0(μ+λaQ)^2 = λaQ(μ+λaQ) + λμPa (1-L_a(μ+ λaQ))/1-(1-λP)L_a(μ+ λaQ)


    Q_a^0 (μ+λaQ) =  λaQ +  μλP L_a(μ+ λaQ)/1-(1-λP)L_a(μ+ λaQ)

We multiply Eq. (<ref>) by ρ(a) and Eq. (<ref>) by ρ(a)a/μ. Then we integrate over the activity obtaining for aQ and P:

    -μP+λ B P - λ^2/μaQ^2 (P-1) + λaQ-2λaQ P=0


    λ/μaQ^2 + ( 1-λ⟨ a ⟩/μ) aQ-λ H P=0

where B= ∫_0^∞ da  ρ(a) a(1-L_a(μ+ λaQ))/1-(1-λP)L_a(μ+ λaQ), H= ∫_0^∞ da  ρ(a) aL_a(μ+ λaQ)/1-(1-λP)L_a(μ+ λaQ).
Setting ⟨ a ⟩ / μ =r = 1/λ, Eq. (<ref>) gives:

    aQ^2 / μ=  H P

Plugging Eq. (<ref>) in Eq. (<ref>) and setting μ=λ⟨ a ⟩, we obtain:

    P(B-⟨ a ⟩ - λ H(P-1))+aQ(1-2P)=0

Since

    B-⟨ a ⟩ - λ H (P-1)    =   ∫_0^∞ da ρ(a) [ a(1-L_a(μ+ λaQ))/1-(1-λP)L_a(μ+ λaQ). 
          - a- . λPaL_a(μ+ λaQ)/1-(1-λP)L_a(μ+ λaQ)] + λ H 
       =   λ H (1-2P)

then Eq. (<ref>) becomes:

    (1-2P)(aQ+λ H P)=0

Since H>0, aQ≥0 and P≥0, Eq.  (<ref>) has only two solutions: the absorbing state (P,aQ)=(0,0) which is unstable because from Eq. (<ref>) we get r_C< 1/λ=r; and P=1/2 which represents an endemic state, stable at r= 1/λ.

In this way we proved that P=1/2 if r=1/λ regardless of the inter-event time distribution Ψ_a(τ) and of the activity distribution ρ(a).
Therefore, as observed in Fig. <ref>, the prevalences P calculated for different α cross at r=1/λ and so there is a switch at r=1/λ in the effect of burstiness on epidemic spreading. 



§ EXPONENTIAL INTER-EVENT TIME DISTRIBUTIONS


In this section we discuss two properties of the model with exponential inter-event time distribution.  First we calculate exactly the critical exponent β=1. 
In this case the probability P_a(t) is governed by the standard activity based mean-field differential equation:

    ∂_t P_a(t) = -μ P_a(t) + λ (1-P_a(t)) [a P(t)+aP(t)]

Imposing the steady state condition ∂_t P_a(t)=0 we get the stationary value P_a^0 of the P_a(t):

    0 = -μ P_a^0 + λ (1-P_a^0)(a P+ aP)

Averaging Eq. (<ref>)  over a we get aP as a function of P:

    aP=μ (λ r -1)/λ (2P-1)P

Plugging Eq. (<ref>) in Eq. (<ref>) we obtain for P_a^0:

    P_a^0 = λ a P (2 P-1)+ μ (λ r -1) P/λ a P (2 P-1)+ μ (λ r -1) P + μ (2P-1)

Multiplying both members for ρ(a) and integrating over the activity we obtain a self-consistency equation for P:

    P= P∫_0^∞ da ρ(a) 2 λ a/μP + λ r -1 - λ a/μ/2 λ a/μP^2 +  (λ r +1 - λ a/μ)P -1

Eq. (<ref>) admits the solution P=0 representing the absorbing state. Other possible solutions can be obtained by solving for P  the equation:

    1 = ∫_0^∞ da ρ(a) 2 λ a/μP + λ r -1 - λ a/μ/2 λ a/μP^2 +  (λ r +1 - λ a/μ)P -1

To evaluate β we assume r ∼ r_C and we can expand Eq. (<ref>) near the critical point,  for P∼ 0. We get:

    1 = ∫_0^∞       da ρ(a) [ 2 λ a/μP + λ r -1 - λ a/μ] [ -1+ P( λ a /μ - λ r -1 ) . 
          .
    + P^2 ( 2 λ^2 a /μ r - λ^2 r^2 -1 -λ^2 a^2/μ^2-2λ r ) + 𝒪(P^2) ]


Integrating over the activity and keeping only the leading orders in P we get:

    P[ P C(r) -  (r-r_C)D(r) ]=0

where 

    C(r)=1+ λ^3 ⟨ a^3 ⟩/μ^3 + 2 λ^3 r^3 - 3 λ^2 r^2 + 3 λ^2 ⟨ a^2 ⟩/μ^2 (1-λ r)

 and 

    D(r)=λ^2( ⟨ a^2 ⟩/λ^2 ⟨ a ⟩^2 -1 )( ⟨ a ⟩/λ (√(⟨ a^2 ⟩)-⟨ a ⟩ ) +r )

and for the exponential case

    r_C=  ⟨ a ⟩/λ (√(⟨ a^2 ⟩)+⟨ a ⟩ )


Since ⟨ a^2⟩>⟨ a ⟩^2 we have that D(r)>0. Moreover close to the critical point:

    C(r_C) = ⟨ a^3 ⟩ + 4 ⟨ a^2 ⟩ ^3/2 + 3 ⟨ a ⟩⟨ a^2 ⟩/(⟨ a ⟩ + √(⟨ a^2 ⟩))^3>0




Hence for r∼ r_C Eq. (<ref>) has two solutions: P=0 i.e. the absorbing state, and for r>r_C 

    P(r)=(r-r_C)D(r_C)/C(r_C)

which means that β=1.

We now consider the behavior of the prevalence P(r) in the active phase as a function of the activity fluctuations, i.e. we vary the exponent ν characterizing the distribution ρ(a) defined by Eq. (<ref>).  We notice that, according to Eq. (<ref>),  the threshold r_C increases by increasing ν, hence at small P close to the epidemic thresholds the prevalence grows with ν. Then the results of <ref> show that at r=1/λ> r_C, for any values of ν, P=1/2,  i.e. the prevalence is independent of activity fluctuations. Finally, for r>1/λ the opposite behavior is expected i.e.  the system presents a larger prevalence at large ν. The overall behavior is illustrated in Figure <ref>: large activity fluctuations and large burstiness give rise to a similar effect on the prevalence i.e. an increase at small r and a decrease at large r.



§ ACKNOWLEDGEMENTS

R.B. acknowledges funding from the INFN BIOPHYS project. M.A.M. acknowledges the Spanish Ministry of Science as well as the Agencia Española de Investigación (AEI) for financial support under grant FIS2017-84256-P (FEDER funds). We warmly thank V. Buendia for very useful comments.



§ REFERENCES



10
url<#>1#1urlprefixurl 



barabasi2005origin
Barabási A L 2005 The origin of bursts and heavy tails in human dynamics
  Nature http://dx.doi.org/10.1038/nature03459435 207

karsai2018bursty
Karsai M, Jo H H and Kaski K 2018 Bursty human dynamics (Springer)

goh2008burstiness
Goh K I and Barabàsi A L 2008 Burstiness and memory in complex systems Europhys. Lett. http://dx.doi.org/10.1209/0295-5075/81/4800281 48002

vazquez2006modeling
Vázquez A, Oliveira J a G, Dezsö Z, Goh K I, Kondor I and Barabási A L
  2006 Modeling bursts and heavy tails in human dynamics Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.73.03612773 036127

karsai2012universal
Karsai M, Kaski K, Barabási A L and Kertész J 2012 Universal features
  of correlated bursty behaviour Sci. Rep.
  http://dx.doi.org/10.1038/srep003972 397

malmgren2008poissonian
Malmgren R D, Stouffer D B, Motter A E and Amaral L A 2008 A poissonian
  explanation for heavy tails in e-mail communication Proc. Natl. Acad.
  Sci. U.S.A http://dx.doi.org/10.1073/pnas.0800332105105
  18153–18158

holme2013temporal
Holme P and Saramäki J 2013 Temporal networks (Springer)

ubaldi2017burstiness
Ubaldi E, Vezzani A, Karsai M, Perra N and Burioni R 2017 Burstiness and tie
  activation strategies in time-varying social networks Sci. Rep.
  http://dx.doi.org/10.1038/srep462257 46225

burioni2017asymptotic
Burioni R, Ubaldi E and Vezzani A 2017 Asymptotic theory of time varying
  networks with burstiness and heterogeneous activation patterns J. Stat.
  Mech. Theory Exp. http://dx.doi.org/10.1088/1742-5468/aa6ce72017 054001

moinet2015burstiness
Moinet A, Starnini M and Pastor-Satorras R 2015 Burstiness and aging in social
  temporal networks Phys. Rev. Lett.
  http://dx.doi.org/10.1103/PhysRevLett.114.108701114 108701

moinet2016aging
Moinet A, Starnini M and Pastor-Satorras R 2016 Aging and percolation dynamics
  in a non-poissonian temporal network model Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.94.02231694 022316

colman2015memory
Colman E R and Vukadinović ć Greetham D 2015
  Memory and burstiness in dynamic networks Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.92.01281792 012817

akbarpour2018diffusion
Akbarpour M and Jackson M O 2018 Diffusion in networks and the virtue of
  burstiness Proc. Natl. Acad. Sci. U.S.A.
  http://dx.doi.org/10.1073/pnas.1722089115115 E6996–E7004

vazquez2007impact
Vazquez A, Rácz B, Lukács A and Barabási A L 2007 Impact of
  non-poissonian activity patterns on spreading processes Phys. Rev.
  Lett. http://dx.doi.org/10.1103/PhysRevLett.98.15870298
  158702

iribarren2009impact
Iribarren J L and Moro E 2009 Impact of human activity patterns on the dynamics
  of information diffusion Phys. Rev. Lett.
  http://dx.doi.org/10.1103/PhysRevLett.103.038702103 038702

karsai2011small
Karsai M, Kivelä M, Pan R K, Kaski K, Kertèsz J, Barabàsi A L and
  Saramäki J 2011 Small but slow world: How network topology and burstiness
  slow down spreading Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.83.02510283 025102

rocha2011simulated
Rocha L E, Liljeros F and Holme P 2011 Simulated epidemics in an empirical
  spatiotemporal network of 50,185 sexual contacts PLOS Comput. Biol.
  http://dx.doi.org/10.1371/journal.pcbi.10011097 1–9

pastor2015epidemic
Pastor-Satorras R, Castellano C, Van Mieghem P and Vespignani A 2015 Epidemic
  processes in complex networks Rev. Mod. Phys.
  http://dx.doi.org/10.1103/RevModPhys.87.92587 925–979

valdano2015analytical
Valdano E, Ferreri L, Poletto C and Colizza V 2015 Analytical computation of
  the epidemic threshold on temporal networks Phys. Rev. X
  http://dx.doi.org/10.1103/PhysRevX.5.0210055 021005

valdano2018epidemic
Valdano E, Fiorentin M R, Poletto C and Colizza V 2018 Epidemic threshold in
  continuous-time evolving networks Phys. Rev. Lett.
  http://dx.doi.org/10.1103/PhysRevLett.120.068302120 068302

starnini2017equivalence
Starnini M, Gleeson J P and Boguñá M 2017 Equivalence between non-markovian
  and markovian dynamics in epidemic spreading processes Phys. Rev.
  Lett. http://dx.doi.org/10.1103/PhysRevLett.118.128301118
  128301

boguna2014simulating
Boguñá M, Lafuerza L F, Toral R and Serrano M A 2014 Simulating
  non-markovian stochastic processes Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.90.04210890 042108

liu2018burst
Liu Q and Van Mieghem P 2018 Burst of virus infection and a possibly largest
  epidemic threshold of non-markovian susceptible-infected-susceptible
  processes on networks Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.97.02230997 022309

van2013non
Van Mieghem P and Van de Bovenkamp R 2013 Non-markovian infection spread
  dramatically alters the susceptible-infected-susceptible epidemic threshold
  in networks Phys. Rev. Lett.
  http://dx.doi.org/10.1103/PhysRevLett.110.108701110 108701

cator2013susceptible
Cator E, Van de Bovenkamp R and Van Mieghem P 2013
  Susceptible-infected-susceptible epidemics on networks with general infection
  and cure times Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.87.01281187 012811

rocha2013bursts
Rocha L E and Blondel V D 2013 Bursts of vertex activation and epidemics in
  evolving networks PLOS Comput. Biol.
  http://dx.doi.org/10.1371/journal.pcbi.10029749 1–9

jo2014analytically
Jo H H, Perotti J I, Kaski K and Kertész J 2014 Analytically solvable model
  of spreading dynamics with non-poissonian processes Phys. Rev. X
  http://dx.doi.org/10.1103/PhysRevX.4.0110414 011041

masuda2013predicting
Masuda N and Holme P 2013 Predicting and controlling infectious disease
  epidemics using temporal networks F1000Prime Rep.
  http://dx.doi.org/10.12703/P5-65:6

perra2012activity
Perra N, Gonçalves B, Pastor-Satorras R and Vespignani A 2012 Activity
  driven modeling of time varying networks Sci. Rep.
  http://dx.doi.org/10.1038/srep004692 469

starnini2013topological
Starnini M and Pastor-Satorras R 2013 Topological properties of a
  time-integrated activity-driven network Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.87.06280787 062807

liu2014controlling
Liu S, Perra N, Karsai M and Vespignani A 2014 Controlling contagion processes
  in activity driven networks Phys. Rev. Lett.
  http://dx.doi.org/10.1103/PhysRevLett.112.118702112 118702

gillespie1976general
Gillespie D T 1976 A general method for numerically simulating the stochastic
  time evolution of coupled chemical reactions J. Comput. Phys.
  http://dx.doi.org/10.1016/0021-9991(76)90041-322 403 – 434

masuda2018gillespie
Masuda N and Rocha L E 2018 A gillespie algorithm for non-markovian stochastic
  processes SIAM Rev. http://dx.doi.org/10.1137/16M105587660 95–115

ubaldi2016asymptotic
Ubaldi E, Perra N, Karsai M, Vezzani A, Burioni R and Vespignani A 2016
  Asymptotic theory of time-varying social networks with heterogeneous activity
  and tie allocation Sci. Rep.
  http://dx.doi.org/10.1038/srep357246 35724

karsai2014time
Karsai M, Perra N and Vespignani A 2014 Time varying networks and the weakness
  of strong ties Sci. Rep.
  http://dx.doi.org/10.1038/srep040014 4001

henkel2008non
Henkel M, Hinrichsen H and Lübeck S 2009 Non-Equilibrium Phase
  Transitions vol 1 Absorbing Phase Transitions (Springer Science &
  Business Media)

tizzani2018epidemic
Tizzani M, Lenti S, Ubaldi E, Vezzani A, Castellano C and Burioni R 2018
  Epidemic spreading and aging in temporal networks with memory Phys. Rev.
  E http://dx.doi.org/10.1103/PhysRevE.98.06231598 062315

boguna2013nature
Boguñá M, Castellano C and Pastor-Satorras R 2013 Nature of the epidemic
  threshold for the susceptible-infected-susceptible dynamics in networks Phys. Rev. Lett.
  http://dx.doi.org/10.1103/PhysRevLett.111.068701111 068701

mata2015lifespan
Mata A S, Boguñá M, Castellano C and Pastor-Satorras R 2015 Lifespan
  method as a tool to study criticality in absorbing-state phase transitions
  Phys. Rev. E http://dx.doi.org/10.1103/PhysRevE.91.05211791 052117

munoz1999avalanches
Muñoz M A, Dickman R, Vespignani A and Zapperi S 1999 Avalanche and spreading
  exponents in systems with absorbing states Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.59.617559 6175–6179

lee2012exploiting
Lee S, Rocha L E, Liljeros F and Holme P 2012 Exploiting temporal network
  structures of human interaction to effectively immunize populations PLoS
  One http://dx.doi.org/10.1371/journal.pone.00364397 1–10

nadini2018epidemic
Nadini M, Sun K, Ubaldi E, Starnini M, Rizzo A and Perra N 2018 Epidemic
  spreading in modular time-varying networks Sci. Rep.
  http://dx.doi.org/10.1038/s41598-018-20908-x8 2352

funk2009spread
Funk S, Gilad E, Watkins C and Jansen V A 2009 The spread of awareness and its
  impact on epidemic outbreaks Proc. Natl. Acad. Sci. U.S.A.
  http://dx.doi.org/10.1073/pnas.0810762106106 6872–6877

rizzo2014effect
Rizzo A, Frasca M and Porfiri M 2014 Effect of individual behavior on epidemic
  spreading in activity-driven networks Phys. Rev. E
  http://dx.doi.org/10.1103/PhysRevE.90.04280190 042801




