Equivalence between non-Markovian and Markovian dynamics in epidemic spreading processes

0.65 MB·

by  Michele Starnini

Equivalence between non-Markovian and Markovian dynamics
in epidemic spreading processes
Michele Starnini,1,2,∗ James P. Gleeson,3 and Mari´an Bogun˜´a1,2,†
1Departament de F´ısica de la Mat`eria Condensada,
Universitat de Barcelona, Mart´ı i Franqu`es 1, E-08028 Barcelona, Spain
2Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, Barcelona, Spain
3MACSI, Department of Mathematics and Statistics, University of Limerick, Ireland
A general formalism is introduced to allow the steady state of non-Markovian processes on
networks to be reduced to equivalent Markovian processes on the same substrates. The example
of an epidemic spreading process is considered in detail, where all the non-Markovian aspects are
shown to be captured within a single parameter, the effective infection rate. Remarkably, this
result is independent of the topology of the underlying network, as demonstrated by numerical
simulations on two-dimensional lattices and various types of random networks. Furthermore, an
analytic approximation for the effective infection rate is introduced, which enables the calculation of
the critical point and of the critical exponents for the non-Markovian dynamics.
7
1
0 Modeling the stochastic dynamics that occur in many the non-Markovian assumption is particularly relevant,
2 natural and technological systems has long depended on as empirical measurements of real diseases—smallpox,

n the Markovian assumption. In a Markov process, the measles, ebola—indicate that the distribution of infec-
a probabilities of the occurrence of future events depend tious periods is far from being exponential [23–26].
J
onlyonthepresentstateofthesystem,beingindependent In this Letter, we consider the non-Markovian SIS epi-

0 of the prior history. This memoryless property implies demic model and show that its steady-state is equivalent
1
that such dynamics can be described by Poisson pro- to a Markovian one with an effective infection rate λ ,
eff
] cesses with fixed rates, which are characterized by an thus encoding all the non-Markovian effects into a sin-
h exponential distribution of the inter-event time between gle parameter. Interestingly, this result is independent
p
consecutive events [1]. The mathematical tractability of of the underlying network topology. Our mathematical
-
c Markov processes enables great simplifications in prob- formalism demonstrates the existence of the effective rate
o lem formulation, leading to spectacular successes in the λ , allowing us to compute it by means of numerical
eff
s
description of many dynamical processes unfolding on simulations, and enables us to derive an approximate an-
.
s
networks [2] and in other complex systems. alytic expression λ , in very good agreement with λ .
c app eff
i The dominance of the Markovian modeling framework The approximate value λ is expected to converge to
s app
y hasrecentlybeenchallengedbytheincreasingavailability λeff close to the epidemic threshold, and therefore the
h of time-resolved data on different kind of interactions, critical point and the set of critical exponents of the non-
p ranging from human activity patterns, including commu- Markovian SIS dynamics, when expressed in terms of
[
nication and mobility [3–5], to natural phenomena [6, 7], λapp, are the same of those of the Markovian case, as we

1 biological processes [8], and biochemical reactions [9]. showbymeansofafinitesizescalinganalysis. Itisworth
v These empirical observations have revealed correlated remarking that our formalism is not restricted to the SIS
5 sequences of events with heavy tailed interevent time dis- model and can be easily extended to any non-Markovian
0
tributions [10], a clear signature that the homogeneous dynamics with a finite set of discrete states, allowing the
8
temporal process description is inadequate and that non- determination of the extent to which such dynamics can
2
0 Markovian dynamics lie at the core of such interactions. be reduced to a Markovian equivalent (with redefined
. Meanwhile, the interest in non-Markovian dynamical parameters) or whether the non-Markovian dynamics are
1
0 processes within the complex systems community has fundamentally different.
7 blossomed, from the points of view of both mathematical Let us consider an undirected and unweighted net-
1 modeling [11–15] and numerical simulation [16, 17]. Par- work topology defined by an adjacency matrix a , with
ij
:
v ticular attention has been devoted to epidemic spreading i,j = 1,··· ,N, and a general, non-Markovian SIS dy-
i on complex networks, representing the diffusion of infor- namics running on top. In this model, nodes exist in
X
mation or disease in a population [18]. Recently, it has either of two states, susceptible or infected. Infected
r
a beenshownthatanon-Markovianinfectiondynamicsdra- nodesdecayspontaneouslytothesusceptiblestateaftera
matically alter the Susceptible-Infected-Susceptible (SIS) randomtimetdistributedasψ (t),thatis,recoveryfrom
R
spreadingprocess[16, 19]. non-Markovianeffects are now the illness does not confer any long lasting immunity, a
known to give qualitatively new behavior in information characteristic present in some sexually transmitted dis-
spreading,e.g.,onsocialnetworks,asrevealedbymeasure- eases [27]. Susceptible nodes may become infected upon
mentsofinter-eventtimesforemailresponses[20,21]and contact with infected neighbors, and we assume that each
retweets on Twitter [22]. In the context of epidemiology, infectious (or active) link, connecting an infected node
2
with a susceptible one, hosts statistically independent i to be infected at time t, ρ (t)≡(cid:104)n (t)(cid:105) [29]
i i
stochastic infection processes, each one controlled by the
N
same interevent distribution ψ (t). In an active link iso- (cid:88)
I ρ˙ (t)=−(cid:104)n (t)δ[t (t)](cid:105)+ a (cid:104)[1−n (t)]n (t)λ[τ (t)](cid:105).
lated from the rest of the system, the susceptible node i i i ij i j ji
j=1
becomes infected after a random time t has elapsed since
(2)
the infection was initiated, with t distributed as ψ (t). If
I The first term in Eq. (2) can be rewritten as (see Supple-
a susceptible node is connected to more than one infected
mentary Information, SI)
neighbor, infection processes take place independently
along each infectious link. (cid:104)n (t)δ(t (t))(cid:105)=ρ (t)(cid:104)δ[t (t)]|n =1(cid:105). (3)
i i i i i
Distributions ψR(t) and ψI(t) allow us to evaluate the Inthelimitt→∞,theonlyinformationwehaveaboutti,
(time-dependent) hazard rates, defined as the probabil- given that node i is infected, is that the recovery time of
ity per unit of time that, given that the event did not node i after infection is longer than t . This implies that
i
take place by a time t since the process was initiated, it theprobabilitydensityoft isgivenbyΨ (t )/(cid:104)t (cid:105),where
i R i R
takes place in the time interval between t and t+dt [28]. (cid:104)t (cid:105) is the average recovery time [28]. By combining this
R
The recovery and infection hazard rates are defined as result with the form of the recovery hazard rate, we can
δ(t)=ψR(t)/ΨR(t) and λ(t)=ψI(t)/ΨI(t), where ΨR(t) write
and Ψ (t) are the corresponding survival probabilities,
I (cid:90) ∞ Ψ (t ) ρst
that is, the probability that a given event takes a time lim(cid:104)n (t)δ(t (t))(cid:105)=ρst R i δ(t )dt = i ,
longer than t. When temporal processes follow Poisson t→∞ i i i 0 (cid:104)tR(cid:105) i i (cid:104)tR(cid:105)
(4)
(Markovian) statistics, both distributions are exponential
where we have defined ρst =lim ρ (t). Similarly, the
and the corresponding hazard rates are constants. i t→∞ i
terms on the right hand side of Eq. (2) can be written as
The SIS dynamics can be fully described by a set of
binary stochastic processes {ni(t)}, i=1,··· ,N; defined (cid:104)[1−ni(t)]nj(t)λ[τji(t)](cid:105)=
(5)
as n (t)=1 if node i is infected at time t and zero if it
i
=(cid:104)[1−n (t)]n (t)(cid:105)(cid:104)λ[τ (t)]|n =0,n =1(cid:105).
is susceptible. The exact stochastic evolution of these i j ji i j
processes can be written as From this equation, we observe that the evolution of
the density ρ (t) depends on the evolution of two-point
i
correlation functions, ρ (t) = (cid:104)n (t)n (t)(cid:105), that appear
ij i j
n (t+dt)=n (t)ξ (t,dt)+[1−n (t)]η (t,dt). (1)
i i i i i in the second term of Eq. (2). Using similar arguments
to those used to derive Eq. (2), we can write an exact
differential equation for the n-point correlation function
Inthisequation,thefirstterminthesumoftherighthand
ρ (see SI):
side is different from zero only when node i is infected i1···in
and accounts for its recovery during the time interval  
N
(its,dte+findte)d. tToobaecheiqeuvaeltthoisz,etrhoewstiothchparsotibcapbriolicteyssdtξδi[(tt,(dt)t]) ρ˙i1···in =(cid:88)(cid:104)−δ(ti)ni+(1−ni)(cid:88)aijnjλ(τji) (cid:89) nk(cid:105),
i i∈I j=1 k∈Ii
and one otherwise, where ti(t) is the time elapsed, at (6)
time t, since node i became infected. Similarly, the sec- where we omit the dependence on t for brevity and we
ond term in the sum of the right hand side of Eq. (1) define the sets of nodes I ≡{i ,i ,···i } and I ≡I\i.
1 2 n i
accounts for the infection of susceptible node i by one of Eq. (6) can be written as
its infected neighbors during the time interval (t,t+dt).
The the stochastic process η (t,dt) is defined to be equal ρ˙ = −ρ (cid:80) δ˜+
(cid:80)i i1···in i1···in i∈I i
to one with probability dt a n (t)λ[τ (t)] and zero (7)
j ij j ji
otherwise, where τ (t) is the time elapsed since the infec- + (cid:80) (cid:80)N a λ˜ [ρ −ρ ]
ji i∈I j=1 ij ji i1···j···in i1···inj
tion process of node j to node i started. Note that we
implicitly assume that each infected neighbor defines a where ρi1···j···in and ρi1···inj are the n and (n+1)-point
statisticallyindependentrandomprocesssothatthetotal correlation functions of the sets Ii ∪{j} and I ∪{j},
infection hazard rate is simply the sum of the infection respectively, and where we have also defined
hazard rates of each individual process. Note also that in δ˜(t)≡(cid:104)δ[t (t)]|{n =1,j ∈I}(cid:105) (8)
this formulation t (t) and τ (t) are themselves stochastic i i j
i ji
processes. and
The average of Eq. (1), first conditioned to the knowl- λ˜ (t)≡(cid:104)λ[τ (t)]|n =0,n =1,{n =1,k ∈I }(cid:105). (9)
ji ji i j k i
edge of the stochastic processes {n ,t ,τ } at time t, and
i i ji
then over the unconditional values, allows us to write the Equations (7), (8), and (9) are the central result of our
following differential equation for the probability of node paperastheyfullydescribethedynamicsoftheepidemic.
3
e j tR then define an effective infection rate λeff as
of nod tI λeff ≡(cid:90)tl→im∞∞(cid:104)λ[τji(t)]|ni =0,nj =1(cid:105)
ate = φ(τji)λ(τji)dτji, (10)
st 0
time where φ(τji) ≡ limt→∞Prob(τji;t|nj = 1,ni = 0) is the
de i probability density of τji, where τji is the time elapsed
o since the start of the infection process from node j to
n
of node i, given that node i is susceptible and node j is
ate infectedandλ[τji(t)]isaveragedoverallactivelinksi−j
st in the network.
time Concerning the average of δ[t (t)], we note that in
i
Figure 1. Sketch of the infection mechanism from node j to general, in the long time limit, Eq. (8) does not re-
node i. Infection events triggered by node j (represented by duce to Eq. (4), since Prob(ti;t|{nk = 1,k ∈ I}) (cid:54)=
stars in the figure) are ineffective if node i is already infected. Prob(t ;t|n = 1). This is due to the fact that, espe-
i i
cially for low-degree nodes, the time t (t) may depend
i
on the status of node i’s neighbors: if at time t node j
is infected and connected only to node i (k =1), node j
j
However, what makes our formulation interesting is the
must have been infected by node i, therefore t has to be
i
factthatallthenon-Markovianeffectsofthedynamicsare
larger than the time elapsed since the infection process
encoded in the terms δ˜ and λ˜ . As we shall show later,
i ji of node j started. Thus, if the recovery process follows
under certain conditions these parameters take constant
a non-Markovian dynamics it is not possible to define
values independent of the nodes, that is, δ˜ = δ˜ and
i an effective parameter δ . Therefore, hereafter we con-
λ˜ = λ˜. In this case, the dynamics, even if strongly eff
ji sider only the case of Markovian recovery, which implies
non-Markovian, can be described by a Markovian one δ˜ =δ˜=(cid:104)t (cid:105)−1,sothatthenon-MarkovianSISdynamics
on the same network, using effective parameters δ˜ and cian be redRuced to a Markovian one with parameters δ˜
λ˜. In this way, the considerable complexity of the non-
and λ .
eff
Markovian effects is reduced to the evaluation of such
Although the probability density φ(τ ) can be easily
ji
effective parameters.
measured in numerical simulation, it is too cumbersome
To proceed further, we need to define the details of the to be computed analytically, even in the simplest case of
pairwise interaction that rules the infection process. We Markovian recovery (see SI). Therefore, we also evaluate
assume that the infection process between an infected anapproximateeffectiveinfectionrateλapp. Todoso, we
node j and a susceptible node i depends on the state of first notice that λeff can also be written as
nodejalone,i.e. whenanodejbecomesinfected,itstarts (cid:90) ∞
aninfectionprocessindependentlytoeachofhisneighbors, λeff =N−1 ψ(τji)λ(τji)dτji, (11)
regardless of their state, according to a renewal process 0
with inter-event time distribution ψI(t). One can think where ψ(τji) ≡ limt→∞Prob(τji,ni = 0;t|nj = 1) is the
of this process as a series of firing events, separated by joint probability that node i is susceptible and the time
random times t , starting when node j becomes infected, elapsed since the last infection attempt from j to i is
I
so that when one such event takes place at a time that equal to τji, given that node j is infected at a given ob-
neighbor node i is susceptible, node i becomes infected servation time t→∞, and where N is the normalization
(cid:82)∞
(see Fig. 1). We also assume that the recovery process of factor N = 0 ψ(τji)dτji. In the SI, we derive an ap-
an infected node depends on its state alone, i.e., when a proximate analytic expression for ψ(τji) that, combined
node becomes infected, it starts a recovery process with withEq.(11), allowsustoderivethefollowingexpression
random time t , distributed as ψ (t). for the approximate effective rate
R R
Within this framework, the average of λ[τji(t)] condi- ψ(cid:98)I(2δ˜)+(cid:104)k(cid:105)ψ(cid:98)I(δ˜)(cid:104)1−ψ(cid:98)I(2δ˜)(cid:105)(cid:104)ψ(cid:98)I(δ˜)−1(cid:105)−1
tionedtothestateofthesystemcanbederivedbynoting λ =
that,attimet,thetimeelapsedsincetheinfectionprocess app [(cid:104)k(cid:105)−1]Ψ(cid:98)I(δ˜)
ofnodej tonodeistarted,τ (t),isnottrulyindependent (12)
ji
of the state of i. That is, if the time elapsed since node i where ψ(cid:98)I(u)≡L{ψI(t)} and Ψ(cid:98)I(u)≡L{ΨI(t)} are the
has recovered is τR, then it holds that τ >τR. There- Laplace transforms of ψ (t) and Ψ (t), respectively, and
i ji i I I
fore, τ depends on the state of i but not on the state (cid:104)k(cid:105) is the average degree of the network substrate.
ji
of any other neighbor and, consequently, Prob(τ ;t|n = We check the validity of the effective infection rates,
ji i
0,n = 1,{n = 1,k ∈ I }) = Prob(τ ;t|n = 0,n = 1) λ and λ , by means of extensive numerical simu-
j k i ji i j eff app
(see SI for a detailed proof). This implies that we can lations of the non-Markovian SIS dynamics, see SI. We
4
1 1 Topology α λ β ν δ
LATT RDR I c ⊥
0.8 0.8 0.5 0.4194 0.596 0.739 0.453
Lattice 2.0 0.4200 0.594 0.727 0.445
0.6 0.6
st st 1.0 0.4122 0.583 0.733 0.4505
r 0.4 r 0.4 0.5 0.3438 1.01 2.06 1.04
RDR 2.0 0.3491 1.01 2.06 1.04
0.2 0.2
1.0 0.3452 1 2 1
0 0
0 0.250.50.75 1 1.251.51.75 2 0 0.250.50.75 1 1.251.51.75 2
l , l l , l
eff app eff app
1 1
Table I. Comparison between the critical point λ and critical
ER SF c
0.8 0.8 exponentsβ,ν⊥ andδ fornon-Markovian(NM)SISdynamics
0.6 a =0.25 0.6 withdifferentexponentαI,andtheMarkoviancase(forwhich
st a I=0.5 st αI =1),ondifferentunderlyingnetworktopologies,2Dlattice
r 0.4 a I=1.5 r 0.4 and RDR network. The critical point λc in NM SIS dynamics
0.2 a I=10 0.2 is evaluated by means of Eq. (12).
I
0 0
0 0.250.50.75 1 1.251.51.75 2 0 0.250.50.75 1 1.251.51.75 2
l , l l , l
eff app eff app
with the same average infection time (cid:104)t (cid:105) but different
I
forms of ψ (t) are known to behave very differently [16],
I
Figure 2. Steady-state prevalence ρst as a function of the showing huge differences in the prevalence ρst for the
effective infection rate, for different values of the exponent same average infection time. This is particularly true in
αI controlling the interevent time infection distribution and the case of highly heterogeneous processes, such as the
different network substrates. Symbols represent the effective
one controlled by α =0.25, with a very skewed form of
infection rate λ , extracted by numerical simulations, con- I
eff the inter-event time distribution ψ (t), and by α =10,
tinuous lines represent the approximate rate λ . I I
app
which corresponds to an almost-periodic process.
The curves plotted as functions of the approximate
consider a Poissonian (Markovian) recovery process with infection rate λapp are also almost indistinguishable from
rate δ˜and an infection process with a Weibull inter-event the others, showing that λapp is a very accurate approx-
time distribution, that is, imation of the exact effective rate, for every underlying
network topology. In the SI we also show that λ is
app
α (cid:18)t(cid:19)αI−1 considerably different from the mean-field approximation
ψ (t)= I e−(t/b)αI, ψ (t)=δ˜e−δ˜t (13)
I b b R proposed in [16], and far more accurate in describing
extreme cases such as α = 0.25 and α = 10, see Sup-
I I
withparameterα controllingthepower-lawstartandtail plementary Fig. 4. Interestingly, as we show in the SI,
I
of the infection inter-event time distribution. We choose Eq. (12) is expected to converge to λ in the limit of
eff
b = (cid:104)t (cid:105)[Γ(1+1/α ))]−1, so that (cid:104)t (cid:105) is the average in- low prevalence ρst (cid:28)1 and, thus, close to the epidemic
I I I
fection time. Hereafter, and without loss of generality, we threshold, λ . This implies that the exact critical point
c
set the time scale to δ˜=1. Once the system has reached λ of the non-Markovian SIS dynamics can be evaluated
c
its steady state, we evaluate λ by selecting random by means of Eq. (12). Using the same argument, we
eff
timeinstantsalongtheprocess. Foreachtimeinstant,we also conclude that the set of critical exponents of the
select all active links, measure the corresponding values non-Markovian dynamics are the same as those of the
of τ , and calculate λ as the average of the hazard Markovian one.
ji eff
rates λ(τji) = αbI (cid:0)τbji(cid:1)αI−1. The approximate infection We check our hypothesis and evaluate the behavior of
rate λ is calculated from Eq. (12) by integrating nu- the non-Markovian SIS dynamics and its critical prop-
app
merically the Laplace transforms ψ(cid:98)I(δ˜) and Ψ(cid:98)I(δ˜) with erties by performing a finite size scaling (FSS) analysis.
δ˜=1. We obtain the epidemic threshold λc, evaluated by means
Figure 2 shows the prevalence ρst at the steady state of Eq. (12), and the set of critical exponents β, ν⊥ and
as a function of λeff and λapp, for different values of δ for a non-Markovian SIS dynamics with αI =0.5 and
αI and different network substrates: a two dimensional αI =2,ontopoftwo-dimensionallatticeanddegreeregu-
latticewithperiodicboundaryconditions,anErdo˝sR´enyi larnetworks,bymeansofthelifespanmethodproposedin
(ER) graph with (cid:104)k(cid:105)=8, a random degree regular (RDR) Ref.[30], seeSIandSupplementaryFig. 5. TableIshows
network with (cid:104)k(cid:105)=8, and a scale-free (SF) network with that the critical point and exponents for these cases are
exponent γ = 2.5. One can see that different curves of in very good agreement with corresponding ones known
the prevalence, corresponding to different forms of the in literature for Markovian SIS dynamics.
infection inter-event time distribution collapse onto one In conclusion, we have demonstrated that non-
another when plotted as a function of λ . This result Markovian SIS dynamics on arbitrary network topologies
eff
is particularly noteworthy since two infection processes can be understood in terms of equivalent Markovian dy-
5
namics on the same substrates. This simplification of [22] K. Lerman, Future Internet 8, 21 (2016).
the temporal nature of discrete-state processes promises [23] N. T. J. BAILEY, Biometrika 43, 15 (1956),
to find application in the wide variety of areas where http://biomet.oxfordjournals.org/content/43/1-
2/15.full.pdf+html.
non-Markovian aspects are recognized as increasingly in-
[24] M. Eichner and K. Dietz, American Jour-
fluential.
nal of Epidemiology 158, 110 (2003),
We acknowledge support from the James S. McDon- http://aje.oxfordjournals.org/content/158/2/110.full.pdf+html.
nell Foundation; the ICREA Academia prize, funded by [25] H. NISHIURA and M. EICHNER, Epidemiology and
the Generalitat de Catalunya; the MINECO projects Infection 135, 1145 (2006).
no. FIS2013-47282-C2-1-P and FIS2016-76830-C2-2-P; [26] G.ChowellandH.Nishiura,BMCMedicine12,1(2014).
[27] M.KeelingandP.Rohani,ModelingInfectiousDiseasesin
Generalitat de Catalunya grant no. 2014SGR608; and
Humans and Animals (Princeton University Press, 2007).
Science Foundation Ireland grant no. 11/PI/1026.
[28] D. R. Cox, Renewal Theory (Methuen, London, 1967).
[29] Notethatwiththisdefinition,theprevalenceofthedisease
at time t is simply given by ρ(t)=N−1(cid:80)N ρ (t).
i=1 i
[30] M.Bogun˜a´,C.Castellano, andR.Pastor-Satorras,Phys.
Rev. Lett. 111, 068701 (2013).
∗ Corresponding author: [email protected]
† Corresponding author: [email protected]
[1] F. A. Haight, “Handbook of the poisson distribution,”
(1967).
[2] A.Barrat,M.Barth´elemy, andA.Vespignani,Dynamical
Processes on Complex Networks (Cambridge University
Press, Cambridge, 2008).
[3] A.-L. Barabasi, Nature 435, 207 (2005).
[4] J. G. Oliveira and A.-L. Barabasi, Nature 437, 1251
(2005).
[5] M. C. Gonzalez, C. A. Hidalgo, and A.-L. Barabasi,
Nature 453, 779 (2008).
[6] A. Corral, Phys. Rev. Lett. 92, 108501 (2004).
[7] M. S. Wheatland, P. A. Sturrock, and J. M. McTiernan,
The Astrophysical Journal 509, 448+ (1998).
[8] T. Kemuriyama, H. Ohta, Y. Sato, S. Maruyama,
M.Tandai-Hiruma,K.Kato, andY.Nishida,Biosystems
101, 144 (2010).
[9] D. Bratsun, D. Volfson, L. S. Tsimring, and J. Hasty,
Proceedings of the National Academy of Sciences of
the United States of America 102, 14593 (2005),
http://www.pnas.org/content/102/41/14593.full.pdf.
[10] M.Karsai, K.Kaski,A.-L.Baraba´si, andJ.Kert´esz,Sci.
Rep. 2 (2012).
[11] A. Moinet, M. Starnini, and R. Pastor-Satorras, Phys.
Rev. Lett. 114, 108701 (2015).
[12] M. Karsai, N. Perra, and A. Vespignani, Sci. Rep. 4,
srep04001 (2014).
[13] G. Garc´ıa-P´erez, M. Bogun˜a´, and M. A´. Serrano, Scien-
tific Reports 5, 9714 EP (2015).
[14] H.-H.Jo,J.I.Perotti,K.Kaski, andJ.Kert´esz,Physical
Review X 4, 011041 (2014).
[15] I. Z. Kiss, G. R¨ost, and Z. Vizi, Phys. Rev. Lett. 115,
078701 (2015).
[16] P. Van Mieghem and R. van de Bovenkamp,
Physical Review Letters 110 (2013), 10.1103/Phys-
RevLett.110.108701.
[17] M.Bogun˜a´,L.F.Lafuerza,R.Toral, andM.A.Serrano,
Phys. Rev. E 90, 042108 (2014).
[18] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and
A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
[19] E. Cator, R. van de Bovenkamp, and P. Van Mieghem,
Physical Review E 87, 062816 (2013).
[20] J. L. Iribarren and E. Moro, Physical review letters 103,
038702 (2009).
[21] J.L.IribarrenandE.Moro,PhysicalReviewE84,046116
(2011).
6
SUPPLEMENTARY INFORMATION
DERIVATION OF THE DIFFERENTIAL EQUATION FOR THE n−POINT CORRELATION FUNCTION
The average of Eq (1) of the main text given the state of the system at time t, n(t)≡(n (t),n (t),··· ,n (t)), can
1 2 N
be written as
(cid:88)
(cid:104)n (t+dt)|n(t)(cid:105)=n (t)+dtA (t), A (t)=−n (t)δ[t (t)]+(1−n (t)) a λ[τ (t)]n (t). (14)
i i i i i i i ij ji j
j
Let us consider a set of n nodes I ≡{i ,i ,···i }. The correlation function between these n nodes reads
1 2 n
(cid:42)(cid:34) (cid:35)(cid:43)
1 (cid:89) (cid:89)
ρ˙ (t)= (cid:104)n (t+dt)|n(t)(cid:105)− n (t) , (15)
i1···in dt i i
i∈I i∈I
where the outer average is over the state of the system at time t. Notice also that the factorization in the first term of
this equation is a direct consequence of the independence of the random variables ξ and η in Eq. (1) of the main text
i i
for different nodes. The first term in Eq (15) can be written, by means of Eq (14) as
n
(cid:89)(cid:104)n (t+dt)|n(t)(cid:105)=(cid:88) (cid:88)(dt)k (cid:89) n (t) (cid:89) A (t), (16)
i i l
i∈I k=0{Ik} i∈Ik l∈Ik
where {I } is the set of all subsets of I containing k nodes. Because of the term dtk, however, the expansion to linear
k
order in dt is a reduced sum over k <2, and thus
(cid:89) (cid:89) (cid:88) (cid:89)
(cid:104)n (t+dt)|n(t)(cid:105)= n (t)+dt A (t) n (t). (17)
i i i k
i∈I i∈I i∈I k∈I\i
Therefore the n−point correlation function reads
(cid:88) (cid:88) (cid:89) (cid:88) (cid:89)
ρ˙ (t)=(cid:104) (1−n (t)) a λ[τ (t)]n (t) n (t)− n (t)δ(t (t)) n (t)(cid:105) (18)
i1···in i ij ji j k i i k
i∈I j k∈I\i i∈I k∈I\i
(cid:42)  (cid:43)
(cid:88) (cid:88) (cid:89)
= (1−ni(t)) aijλ[τji(t)]nj(t)−ni(t)δ(ti(t)) nk(t) , (19)
i∈I j k∈I\i
from which Eq (6) of the main text follows immediately.
TIME τ OF ACTIVE LINK i−j DOES NOT DEPEND ON THE STATES OF OTHER NODES
ij
DIFFERENT FROM i AND j
A critical step in our approach is to prove that
Prob(τ ;t|n =0,n =1,{n =1,k ∈I })=Prob(τ ;t|n =0,n =1). (20)
ji i j k i ji i j
The probability in the left hand side of this equation can be written as
(cid:90) (cid:90)
(cid:89)
Prob(τ ;t|n =0,n =1,{n =1,k ∈I })= ··· φ(τR,τI,{τI};t)φ(τ |τR,τI,{τI})dτRdτI dτI, (21)
ji i j k i i j k ji i j k i j k
k∈Ii
where φ(τR,τI,{τI};t) is the joint probability density, at time t, that given that node i is susceptible and nodes j
i j k
and {k ∈I } are infected, the time elapsed since i recovered is τR and the times elapsed since j and {k ∈I } became
i i i
infected are τI and {τI}, respectively. By Bayes’ rule, φ(τ |τR,τI,{τI}) is the probability density of the time τ
j k ji i j k ji
conditioned on the times τR,τI,{τI}. However, it is easy to see that since infection events take place in active links
i j k
independently, once τR and τI are fixed, τ is totally independent of the elapsed times since nodes other than j
i j ji
became infected. Therefore,
φ(τ |τR,τI,{τI})=φ(τ |τR,τI), (22)
ji i j k ji i j
which directly gives the result in Eq. (20).
7
GENERAL FORMALISM FOR λ
eff
Using the result in Eq. (22), at the steady state the probability density φ(τ ) of the time elapsed since the infection
ji
process of node j to node i started, given that node i is susceptible and node j is infected, can be written in general as
(cid:90) (cid:90)
φ(τ )= φ(τ |τI,τR)φ(τI,τR)dτIdτR, (23)
ji ji j i j i j i
where φ(τI,τR) is the joint probability that the time elapsed since j became infected is equal to τI and the time
j i j
elapsedsinceirecoveredisequaltoτR. Ifweassumethatthetwoprocessareuncorrelated, φ(τI,τR)canbefactorized
i j i
into φ(τI,τR)=φ (τI)φ (τR), and Eq. (23) reduces to
j i I j R i
(cid:90) ∞ (cid:90) ∞ (cid:110) (cid:111)
φ(τ )= dτIφ (τI) dτRφ (τR) Θ(τI −τR)φ(τ |τR ≤τI)+Θ(τR−τI)φ(τ |τR >τI) , (24)
ji j I j i R i j i ji i j i j ji i j
0 0
where Θ(t) is the Heaviside step function, φ (τI) is the probability that the time elapsed since j became infected is
I j
equal to τI and φ (τR) is the probability that the time elapsed since i recovered is equal to τR. The conditional
j R i i
probability φ(τ |τR >τI) is simply φ(τ |τR >τI)=δ(τ −τI), and
ji i j ji i j ji j
(cid:90) ∞ ∞
(cid:88)
φ(τ |τR ≤τI)= Θ(τI −τR−τ)δ(τ −(τI −τ))Ψ (τI −τR−τ) P (τ)dτ, (25)
ji i j j i ji j I j i n
0 n=0
where n is the number of infection attempts of node j to node i, P (τ) is the probability that the time elapsed since
n
node j became infected and the moment of his n-th fire is equal to τ, and Ψ (τI −τR−τ) is the probability that the
I j i
time elapsed between the n-th fire and the n+1-th fire is greater than τ −τR. As computed in Eq. (4) of the main
ji i
text, the probability that the time elapsed since j became infected is equal to τI is simply
j
φ (τI)=δ˜Ψ (τI). (26)
I j R j
The survival probability Ψ (τI) of recovery events can be written as
R j
(cid:90) ∞
ΨR(τjI)= ω(u)e−uτjIdu=ω(cid:98)(τjI), (27)
0
where ω(u) is the inverse Laplace transform of Ψ (τI). In Laplace space, the probability distribution P (τ) has a
R j n
(cid:104) (cid:105)n
convenient form, P(cid:98)n(u)= ψ(cid:98)I(u) , where ψ(cid:98)I(u) is the Laplace transform of ψI(t). By inserting Eqs. (26) and (27)
into Eq. (24) we obtain
(cid:40) (cid:41)
(cid:90) ∞ (cid:90) ∞ 1
φ(τ )=δ˜ dτRφ (τR) due−uτjiω(u) θ(τR−τ )+θ(τ −τR)Ψ (τ −τR) . (28)
ji i R i i ji ji i I ji i
0 0 1−ψ(cid:98)I(u)
By inserting the form of φ(τ ) into Eq. (10) of the main text, we obtain an expression for the infection rate λ
ji eff
 
(cid:90) ∞  (cid:92) (cid:86) 1 
λ = duω(u) [λ Φ ](u)+[λ [φ ∗Ψ ]](u) , (29)
eff I R I R I
0  1−ψ(cid:98)I(u)
whereλ istheinfectionhazardrate, Φ (τR)andΨ (t)arethesurvivalprobabilitiesofφ (τR)andψ (t), respectively
I R i I R i I
and φ ∗Ψ is the convolution between φ (τR) and Ψ (t). At this point, some ansatz regarding the form of φ (τR),
R I R i I R i
the probability that the time elapsed since i recovered is equal to τR, is needed to continue. We note that if one does
i
not consider the state of node i in the probability φ (τ ), which corresponds to inserting φ (τR)=δ(τR) into Eq.
R ji R i i
(28), one obtains
λ =δ˜(cid:90) ∞ω(u) ψ(cid:98)I(u) du. (30)
mf
0 1−ψ(cid:98)I(u)
8
a b
j j

e t e t =tI
d ji d ji j
o o
n n

f f
o o

e e
t t
a a
t t
s s
time time
i i

e e
d observation time t d observation time t
o o
n n

f f
o o

e e
t t
a a
t t
s s
time time
Figure 3. Sketch of two possible ways to obtain the time τ and, simultaneously, node i infected at the observation time. In a,
ji
node j has attempted to infect i at least once at time t−τ . This implies that node i must necessarily be infected at that
ji
moment and, thus, it has to recover before the observation time. In b, node j has not attempted to infect i since it became
infected by node i. In this case, we know that i was infected when j became infected and, again, it has to recover before the
observation time.
Thiseffectiveinfectionrateλ ,alreadyfoundinCatoretal.[19]byusingameanfieldapproximation,isnowobtained
mf
within a more general formalism. A different possibility is to consider φ (τR) equal to an exponential distribution,
R i
φR(τiR)=we−wτiR, (31)
with rate w. The rate w can be written as a simple function of the prevalence ρ, w =δρ/(1−ρ). However, even if it
were not possible to find a closed analytic form for the effective infection rate, one can resort to numerical simulation
in order to compute λ . One can see that the exponential ansatz for the form of φ (τR) is correct for large values of
eff R i
the prevalence ρ, but it fails for low prevalence, thus close to the epidemic threshold.
APPROXIMATIONS TO λ
app
To find an infection rate which is accurate and analytically treatable close to the epidemic threshold, we follow a
different approach. As stated in the main text, we consider here the probability density ψ(τ )≡lim p(τ ,n =
ji t→∞ ji i
0;t|n =1),whichisthejoinprobabilitythat,giventhatnodej isinfectedattheobservationtime,nodeiissusceptible
j
and the time elapsed since the last infection attempt from j to i is equal to τ . Our approximation consists of
ji
estimating the probability that node i is susceptible at the observation time t, which depends on the time instant at
which node i became infected, this time instant being unknown in principle.
We first consider the case of low prevalence, ρst (cid:28)1. Then, let us consider separately the cases in which node j
attempts at least once to infect node i, n>0, and the case of no attempts, n=0. In the first case, at the time of
the last infection event from j to i, t−τ , node i either was already infected (in which case the infection attempt is
ji
ineffective) or it became infected by this event (see Fig. 3a). In both cases, we are certain that node i is in an infected
state at time t−τ and, thus, the probability that node i recovers before the observation time t is 1−Ψ (τ ). If
ji R ji
the prevalence is low, the probability that node i is subsequently infected by one of its neighbor (other than j) and
then recovers before the observation time is also very low, and we assume it to be zero. With these assumption, the
probability that node i is susceptible at time t is simply 1−Ψ (τ ). If node j does not attempt to infect node i, we
R ji
cannot know for certain the state of node i. However, given that node j became infected at time t−τI, one of his
j
neighbors must have infected him. Let us consider that node j has degree k. If the prevalence is low, it is very unlikely
to find more than one neighbor of node j infected simultaneously and we assume that only one of his neighbors was
9
infected and infected him. With probability 1/k, such infected node is node i (See Fig. 3b), so that i is infected at time
t−τI, and the probability that node i is then susceptible at the observation time t is 1−Ψ (τI). With probability
j R j
1−1/k, the infected node is a neighbor other than i and, thus, we assume that node i was susceptible at time t−τI
j
and it will remain in this state until time t with probability equal to one. Summing up, if node j attempts at least
once to infect node i, then the probability that it is susceptible at the observation time is 1−Ψ (τ ). Instead, if node
R ji
j does not attempt to infect node i, this probability reads (1−Ψ (τI))/k+(k−1)/k =(k−Ψ (τI))/k. In the limit
R j R j
of low prevalence, we expect this approximation to be exact. In the following, we also approximate the value of k by
the average degree, (cid:104)k(cid:105).
The probability density ψ(τ ) can be written as
ji
(cid:90) ∞
ψ(τ )= ψ(τ |τI)φ (τI)dτI, (32)
ji ji j I j j
0
where again φ (τI) is the probability that the time elapsed since j became infected is equal to τI. The conditional
I j j
probability ψ(τ |τI) is
ji j
ψ(τ |τI)=δ(τ −τI)Ψ (τI)(cid:34)(cid:104)k(cid:105)−ΨR(τjI)(cid:35)+(cid:90) τjI δ(τ −(τI −τ))[1−Ψ (τ )]Ψ (τI −τ)(cid:88)∞ P (τ)dτ (33)
ji j ji j I j (cid:104)k(cid:105) ji j R ji I j n
0 n=1
where the first term accounts for the case in which there are no infection attempts from j to i, n=0, while the second
term accounts for the case n>0. By inserting Eq. (26) into Eq. (32) and integrating over τI, the probability ψ(τ )
j ji
reads
ψ(τ )=δ˜Ψ (τ )(cid:40)Ψ (τ )(cid:20)(cid:104)k(cid:105)−ΨR(τji)(cid:21)+[1−Ψ (τ )](cid:90) ∞Ψ (τ +τ)(cid:88)∞ P (τ)dτ(cid:41). (34)
ji I ji R ji (cid:104)k(cid:105) R ji R ji n
0 n=1
If we restrict to the case of Markovian recovery, we can use its memoryless property, Ψ (τ +τ)=Ψ (τ )Ψ (τ). By
R ji R ji R
using the convenient Laplacian form of the probability distribution P (τ), one can obtain
n
Ψ (τ )Ψ (τ )(cid:110) (cid:104) (cid:16) (cid:17) (cid:105)(cid:111)
ψ(τji)=δ˜ I1−jiψ(cid:98)IR(δ˜)ji 1−ΨR(τji) (cid:104)k(cid:105)−1 1−ψ(cid:98)I(δ˜) +ψ(cid:98)I(δ˜) . (35)
The normalization of ψ(τ ) reads
ji
(cid:90) ∞ δ˜ (cid:110) (cid:104) (cid:16) (cid:17) (cid:105)(cid:111)
N = 0 ψ(τji)dτji = 1−ψ(cid:98)I(δ˜) Ψ(cid:98)I(δ˜)−Ψ(cid:98)I(2δ˜) (cid:104)k(cid:105)−1 1−ψ(cid:98)I(δ˜) +ψ(cid:98)I(δ˜) , (36)
therefore, by inserting Eq. (35) and Eq. (36) into Eq. (11) of the main text, one finally obtains the approximate
infection rate λ presented in Eq. (12) of the main text.
app
In Fig. 4, we show the steady-state prevalence ρst as a function of the approximate effective infection rate λ ,
app
given by Eq. (12) of the main text, and the mean field effective rate λ given by Eq. (30) for two extreme values of
mf
the exponent α controlling the interevent time infection distribution, α =0.25 and α =10. One can see that the
I I I
curves for λ do not collapse onto one another, especially for the lattice and SF network substrate.
mf
NUMERICAL SIMULATIONS OF THE NON-MARKOVIAN SIS DYNAMICS
To check the validity of the effective infection rates λ and λ , we run extensive numerical simulations of the
eff app
non-Markovian SIS dynamics. For each value of the average infection time (cid:104)t (cid:105) and fixed average recovery time
I
(cid:104)t (cid:105) = δ˜−1 = 1, we simulate the non-Markovian SIS dynamics by implementing an algorithm based on a queue of
R
infection and recovery events.
At time t=0, all nodes are in a susceptible state, and a set of fN randomly chosen nodes, with f =0.5, change
their state to the infected one. In the algorithm, whenever a node i changes his state from susceptible to infected at
time t, he first randomly extracts his recovery time t from the distribution ψ (t), and pushes his recovery event at
R R
time t+t to the queue. He also starts k independent infection processes to his k neighbors. In each infection process
R
to a neighbor j, an infection event from node i to node j is scheduled at time t+t1, where t1 is randomly extracted
I I
10
1
LATT RDR
0.8
0.8
0.7
0.6 0.6
st st
r 0.5 r 0.4
0.4
0.2
0.3
0.2 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
l , l l , l
app mf app mf
1
ER 0.8 SF
0.8
0.6
0.6 a =0.25
I
st a =10 st 0.4
r 0.4 I r
0.2 0.2
0 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
l , l l , l
app mf app mf
Figure 4. Prevalence ρ as a function of the approximate effective infection rate λ (points) and the mean field effective rate
app
λ (continuous line), for different values of the exponent α and different network substrate.
mf I
from the distribution ψ (t), only if t1 <t , that is if node i is still infected at time t+t1. A second infection event
I I R I
from node i to node j is scheduled at time t+t1+t2, where t2 is randomly extracted from the distribution ψ (t), only
I I I I
if t1+t2 <t , and so on until n (with possibly n=0) infection events are generated and pushed to the queue.
I I R
The queue is pulled by following the time order of the events. If the pulled event is the recovery of node i, i changes
his state from infected to susceptible. If the pulled event is an infection event from node i to node j, and j is already
in a infected state, nothing happens, otherwise node j changes his state from susceptible to infected and schedules his
recovery and infection events, pushing them to the queue. The queue is pulled until either no more events are left (and
so all nodes are susceptible) or the time reaches a time T , set conveniently. In order to measure the prevalence
max
in the steady state ρst and the effective infection rate λ , we sample N =104 time instants uniformly chosen in
eff s
[T ,T ], with T chosen such that the stationary state is reached long before it. For each time instant, we
min max min
measure the prevalence and the values of τ for each active link between nodes i and j, so as to calculate λ by
ij eff
means of Eq. (10) of the main text.
We have double checked our event queue algorithm by simulating the non-Markovian SIS dynamics with a non-
Markovian Gillespie algorithm [17], which is much slower, and we obtained identical results for the prevalence and the
effective infection rate.
EPIDEMIC THRESHOLD AND CRITICAL EXPONENTS
We run extensive numerical simulations of the non-Markovian SIS dynamics in order to evaluate its behavior close
to the epidemic threshold. We consider α =0.5 and α =2, and two different network substrates, 2D lattice and
I I
RDR network. We address the critical properties by means of the lifespan method [30], in which the infection starts
with a single infected node. In the lifespan method, each realization is characterized by its lifetime, T, and its coverage,
C, defined as the number of distinct nodes that have become infected at least once. We let each realization run until
either the coverage C reaches a certain threshold C∗ (and we consider it endemic), or the realization dies out, and
we measure its lifetime T. We set C∗ =ΘN, with Θ=0.9. We then measure the probability of having an endemic
realization P, the average lifetime (cid:104)T(cid:105) and average square lifetime (cid:104)T2(cid:105) over a number of runs N , as a function of
run
the average infection time (cid:104)t (cid:105) (corresponding to an effective rate λ , hereafter λ for brevity) close to the epidemic
I app