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

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.

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.
Modeling the stochastic dynamics that occur in many natural and technological systems has long depended on the Markovian assumption. In a Markov process, the probabilities of the occurrence of future events depend only on the present state of the system, being independent of the prior history. This memoryless property implies that such dynamics can be described by Poisson processes with fixed rates, which are characterized by an exponential distribution of the inter-event time between consecutive events [1]. The mathematical tractability of Markov processes enables great simplifications in problem formulation, leading to spectacular successes in the description of many dynamical processes unfolding on networks [2] and in other complex systems.
The dominance of the Markovian modeling framework has recently been challenged by the increasing availability of time-resolved data on different kind of interactions, ranging from human activity patterns, including communication and mobility [3][4][5], to natural phenomena [6,7], biological processes [8], and biochemical reactions [9]. These empirical observations have revealed correlated sequences of events with heavy tailed interevent time distributions [10], a clear signature that the homogeneous temporal process description is inadequate and that non-Markovian dynamics lie at the core of such interactions.
Meanwhile, the interest in non-Markovian dynamical processes within the complex systems community has blossomed, from the points of view of both mathematical modeling [11][12][13][14][15] and numerical simulation [16,17]. Particular attention has been devoted to epidemic spreading on complex networks, representing the diffusion of information or disease in a population [18]. Recently, it has been shown that a non-Markovian infection dynamics dramatically alter the Susceptible-Infected-Susceptible (SIS) spreading process [16,19]. non-Markovian effects are now known to give qualitatively new behavior in information spreading, e.g., on social networks, as revealed by measurements of inter-event times for email responses [20,21] and retweets on Twitter [22]. In the context of epidemiology, the non-Markovian assumption is particularly relevant, as empirical measurements of real diseases-smallpox, measles, ebola-indicate that the distribution of infectious periods is far from being exponential [23][24][25][26].
In this Letter, we consider the non-Markovian SIS epidemic model and show that its steady-state is equivalent to a Markovian one with an effective infection rate λ ef f , thus encoding all the non-Markovian effects into a single parameter. Interestingly, this result is independent of the underlying network topology. Our mathematical formalism demonstrates the existence of the effective rate λ ef f , allowing us to compute it by means of numerical simulations, and enables us to derive an approximate analytic expression λ app , in very good agreement with λ ef f . The approximate value λ app is expected to converge to λ ef f close to the epidemic threshold, and therefore the critical point and the set of critical exponents of the non-Markovian SIS dynamics, when expressed in terms of λ app , are the same of those of the Markovian case, as we show by means of a finite size scaling analysis. It is worth remarking that our formalism is not restricted to the SIS model and can be easily extended to any non-Markovian dynamics with a finite set of discrete states, allowing the determination of the extent to which such dynamics can be reduced to a Markovian equivalent (with redefined parameters) or whether the non-Markovian dynamics are fundamentally different.
Let us consider an undirected and unweighted network topology defined by an adjacency matrix a ij , with i, j = 1, · · · , N , and a general, non-Markovian SIS dynamics running on top. In this model, nodes exist in either of two states, susceptible or infected. Infected nodes decay spontaneously to the susceptible state after a random time t distributed as ψ R (t), that is, recovery from the illness does not confer any long lasting immunity, a characteristic present in some sexually transmitted diseases [27]. Susceptible nodes may become infected upon contact with infected neighbors, and we assume that each infectious (or active) link, connecting an infected node with a susceptible one, hosts statistically independent stochastic infection processes, each one controlled by the same interevent distribution ψ I (t). In an active link isolated from the rest of the system, the susceptible node becomes infected after a random time t has elapsed since the infection was initiated, with t distributed as ψ I (t). If a susceptible node is connected to more than one infected neighbor, infection processes take place independently along each infectious link.
Distributions ψ R (t) and ψ I (t) allow us to evaluate the (time-dependent) hazard rates, defined as the probability per unit of time that, given that the event did not take place by a time t since the process was initiated, it takes place in the time interval between t and t + dt [28]. The recovery and infection hazard rates are defined as δ(t) = ψ R (t)/Ψ R (t) and λ(t) = ψ I (t)/Ψ I (t), where Ψ R (t) and Ψ I (t) are the corresponding survival probabilities, that is, the probability that a given event takes a time longer than t. When temporal processes follow Poisson (Markovian) statistics, both distributions are exponential and the corresponding hazard rates are constants.
The SIS dynamics can be fully described by a set of binary stochastic processes {n i (t)}, i = 1, · · · , N ; defined as n i (t) = 1 if node i is infected at time t and zero if it is susceptible. The exact stochastic evolution of these processes can be written as In this equation, the first term in the sum of the right hand side is different from zero only when node i is infected and accounts for its recovery during the time interval (t, t + dt). To achieve this, the stochastic process ξ i (t, dt) is defined to be equal to zero with probability dtδ[t i (t)] and one otherwise, where t i (t) is the time elapsed, at time t, since node i became infected. Similarly, the second term in the sum of the right hand side of Eq. (1) accounts for the infection of susceptible node i by one of its infected neighbors during the time interval (t, t + dt). The the stochastic process η i (t, dt) is defined to be equal to one with probability dt j a ij n j (t)λ[τ ji (t)] and zero otherwise, where τ ji (t) is the time elapsed since the infection process of node j to node i started. Note that we implicitly assume that each infected neighbor defines a statistically independent random process so that the total infection hazard rate is simply the sum of the infection hazard rates of each individual process. Note also that in this formulation t i (t) and τ ji (t) are themselves stochastic processes.
The average of Eq. (1), first conditioned to the knowledge of the stochastic processes {n i , t i , τ ji } at time t, and then over the unconditional values, allows us to write the following differential equation for the probability of node i to be infected at time t, ρ i (t) ≡ n i (t) [29] (2) The first term in Eq. (2) can be rewritten as (see Supplementary Information, SI) In the limit t → ∞, the only information we have about t i , given that node i is infected, is that the recovery time of node i after infection is longer than t i . This implies that the probability density of t i is given by Ψ R (t i )/ t R , where t R is the average recovery time [28]. By combining this result with the form of the recovery hazard rate, we can write Similarly, the terms on the right hand side of Eq. (2) can be written as From this equation, we observe that the evolution of the density ρ i (t) depends on the evolution of two-point correlation functions, ρ ij (t) = n i (t)n j (t) , that appear 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 ρ i1···in (see SI): where we omit the dependence on t for brevity and we define the sets of nodes I ≡ {i 1 , i 2 , · · · i n } and I i ≡ I \ i. Eq. (6) can be written aṡ where ρ i1···j···in and ρ i1···inj are the n and (n + 1)-point correlation functions of the sets I i ∪ {j} and I ∪ {j}, respectively, and where we have also defined and  However, what makes our formulation interesting is the fact that all the non-Markovian effects of the dynamics are encoded in the termsδ i andλ ji . As we shall show later, under certain conditions these parameters take constant values independent of the nodes, that is,δ i =δ and λ ji =λ. In this case, the dynamics, even if strongly non-Markovian, can be described by a Markovian one on the same network, using effective parametersδ and λ. In this way, the considerable complexity of the non-Markovian effects is reduced to the evaluation of such effective parameters.
To proceed further, we need to define the details of the pairwise interaction that rules the infection process. We assume that the infection process between an infected node j and a susceptible node i depends on the state of node j alone, i.e. when a node j becomes infected, it starts an infection process independently to each of his neighbors, regardless of their state, according to a renewal process with inter-event time distribution ψ I (t). One can think of this process as a series of firing events, separated by random times t I , starting when node j becomes infected, so that when one such event takes place at a time that neighbor node i is susceptible, node i becomes infected (see Fig. 1). We also assume that the recovery process of an infected node depends on its state alone, i.e., when a node becomes infected, it starts a recovery process with random time t R , distributed as ψ R (t).
Within this framework, the average of λ[τ ji (t)] conditioned to the state of the system can be derived by noting that, at time t, the time elapsed since the infection process of node j to node i started, τ ji (t), is not truly independent of the state of i. That is, if the time elapsed since node i has recovered is τ R i , then it holds that τ ji > τ R i . Therefore, τ ji depends on the state of i but not on the state of any other neighbor and, consequently, Prob(τ ji ; t|n i = 0, n j = 1, {n k = 1, k ∈ I i }) = Prob(τ ji ; t|n i = 0, n j = 1) (see SI for a detailed proof). This implies that we can then define an effective infection rate λ ef f as where φ(τ ji ) ≡ lim t→∞ Prob(τ ji ; t|n j = 1, n i = 0) is the probability density of τ ji , where τ ji is the time elapsed since the start of the infection process from node j to node i, given that node i is susceptible and node j is infected and λ[τ ji (t)] is averaged over all active links i − j in the network.
Concerning the average of δ[t i (t)], we note that in general, in the long time limit, Eq. (8) does not reduce to Eq. (4), since Prob(t i ; t|{n k = 1, k ∈ I}) = Prob(t i ; t|n i = 1). This is due to the fact that, especially for low-degree nodes, the time t i (t) may depend on the status of node i's neighbors: if at time t node j is infected and connected only to node i (k j =1), node j must have been infected by node i, therefore t i has to be larger than the time elapsed since the infection process of node j started. Thus, if the recovery process follows a non-Markovian dynamics it is not possible to define an effective parameter δ ef f . Therefore, hereafter we consider only the case of Markovian recovery, which implies δ i =δ = t R −1 , so that the non-Markovian SIS dynamics can be reduced to a Markovian one with parametersδ and λ ef f .
Although the probability density φ(τ ji ) can be easily measured in numerical simulation, it is too cumbersome to be computed analytically, even in the simplest case of Markovian recovery (see SI). Therefore, we also evaluate an approximate effective infection rate λ app . To do so, we first notice that λ ef f can also be written as where ψ(τ ji ) ≡ lim t→∞ Prob(τ ji , n i = 0; t|n j = 1) is the joint probability that node i is susceptible and the time elapsed since the last infection attempt from j to i is equal to τ ji , given that node j is infected at a given observation time t → ∞, and where N is the normalization factor N = ∞ 0 ψ(τ ji )dτ ji . In the SI, we derive an approximate analytic expression for ψ(τ ji ) that, combined with Eq. (11), allows us to derive the following expression for the approximate effective rate (12) where ψ I (u) ≡ L{ψ I (t)} and Ψ I (u) ≡ L{Ψ I (t)} are the Laplace transforms of ψ I (t) and Ψ I (t), respectively, and k is the average degree of the network substrate. We check the validity of the effective infection rates, λ ef f and λ app , by means of extensive numerical simulations of the non-Markovian SIS dynamics, see SI. We consider a Poissonian (Markovian) recovery process with rateδ and an infection process with a Weibull inter-event time distribution, that is, (13) with parameter α I controlling the power-law start and tail of the infection inter-event time distribution. We choose b = t I [Γ(1 + 1/α I ))] −1 , so that t I is the average infection time. Hereafter, and without loss of generality, we set the time scale toδ = 1. Once the system has reached its steady state, we evaluate λ ef f by selecting random time instants along the process. For each time instant, we select all active links, measure the corresponding values of τ ji , and calculate λ ef f as the average of the hazard rates λ(τ ji ) = α I b τji b α I −1 . The approximate infection rate λ app is calculated from Eq. (12) by integrating numerically the Laplace transforms ψ I (δ) and Ψ I (δ) with δ = 1. Figure 2 shows the prevalence ρ st at the steady state as a function of λ ef f and λ app , for different values of α I and different network substrates: a two dimensional lattice with periodic boundary conditions, an Erdős Rényi (ER) graph with k = 8, a random degree regular (RDR) network with k = 8, and a scale-free (SF) network with exponent γ = 2.5. One can see that different curves of the prevalence, corresponding to different forms of the infection inter-event time distribution collapse onto one another when plotted as a function of λ ef f . This result is particularly noteworthy since two infection processes with the same average infection time t I but different forms of ψ I (t) are known to behave very differently [16], showing huge differences in the prevalence ρ st for the same average infection time. This is particularly true in the case of highly heterogeneous processes, such as the one controlled by α I = 0.25, with a very skewed form of the inter-event time distribution ψ I (t), and by α I = 10, which corresponds to an almost-periodic process.
The curves plotted as functions of the approximate infection rate λ app are also almost indistinguishable from the others, showing that λ app is a very accurate approximation of the exact effective rate, for every underlying network topology. In the SI we also show that λ app is considerably different from the mean-field approximation proposed in [16], and far more accurate in describing extreme cases such as α I = 0.25 and α I = 10, see Supplementary Fig. 4. Interestingly, as we show in the SI, Eq. (12) is expected to converge to λ ef f in the limit of low prevalence ρ st 1 and, thus, close to the epidemic threshold, λ c . This implies that the exact critical point λ c of the non-Markovian SIS dynamics can be evaluated by means of Eq. (12). Using the same argument, we also conclude that the set of critical exponents of the non-Markovian dynamics are the same as those of the Markovian one.
We check our hypothesis and evaluate the behavior of the non-Markovian SIS dynamics and its critical properties by performing a finite size scaling (FSS) analysis. We obtain the epidemic threshold λ c , evaluated by means of Eq. (12), and the set of critical exponents β, ν ⊥ and δ for a non-Markovian SIS dynamics with α I = 0.5 and α I = 2, on top of two-dimensional lattice and degree regular networks, by means of the lifespan method proposed in Ref. [30], see SI and Supplementary Fig. 5. Table I shows that the critical point and exponents for these cases are in very good agreement with corresponding ones known in literature for Markovian SIS dynamics.
In conclusion, we have demonstrated that non-Markovian SIS dynamics on arbitrary network topologies can be understood in terms of equivalent Markovian dy-namics on the same substrates. This simplification of the temporal nature of discrete-state processes promises to find application in the wide variety of areas where non-Markovian aspects are recognized as increasingly influential.

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 1 (t), n 2 (t), · · · , n N (t)), can be written as Let us consider a set of n nodes I ≡ {i 1 , i 2 , · · · i n }. The correlation function between these n nodes readṡ 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 ξ i and η i in Eq. (1) of the main text for different nodes. The first term in Eq (15) can be written, by means of Eq (14) as where {I k } is the set of all subsets of I containing k nodes. Because of the term dt k , however, the expansion to linear order in dt is a reduced sum over k < 2, and thus Therefore the n−point correlation function readṡ from which Eq (6) of the main text follows immediately.

TIME τij OF ACTIVE LINK i − j DOES NOT DEPEND ON THE STATES OF OTHER NODES DIFFERENT FROM i AND j
A critical step in our approach is to prove that Prob(τ ji ; t|n i = 0, n j = 1, {n k = 1, k ∈ I i }) = Prob(τ ji ; t|n i = 0, n j = 1).
The probability in the left hand side of this equation can be written as where φ(τ R i , τ I j , {τ I k }; t) is the joint probability density, at time t, that given that node i is susceptible and nodes j and {k ∈ I i } are infected, the time elapsed since i recovered is τ R i and the times elapsed since j and {k ∈ I i } became infected are τ I j and {τ I k }, respectively. By Bayes' rule, φ(τ ji |τ R i , τ I j , {τ I k }) is the probability density of the time τ ji conditioned on the times τ R i , τ I j , {τ I k }. However, it is easy to see that since infection events take place in active links independently, once τ R i and τ I j are fixed, τ ji is totally independent of the elapsed times since nodes other than j became infected. Therefore, which directly gives the result in Eq. (20).

GENERAL FORMALISM FOR λ ef f
Using the result in Eq. (22), at the steady state the probability density φ(τ ji ) of the time elapsed since the infection 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 where φ(τ I j , τ R i ) is the joint probability that the time elapsed since j became infected is equal to τ I j and the time elapsed since i recovered is equal to τ R i . If we assume that the two process are uncorrelated, φ(τ I j , τ R i ) can be factorized into φ(τ I j , τ R i ) = φ I (τ I j )φ R (τ R i ), and Eq. (23) reduces to where Θ(t) is the Heaviside step function, φ I (τ I j ) is the probability that the time elapsed since j became infected is equal to τ I j and φ R (τ R i ) is the probability that the time elapsed since i recovered is equal to τ R i . The conditional where n is the number of infection attempts of node j to node i, P n (τ ) is the probability that the time elapsed since node j became infected and the moment of his n-th fire is equal to τ , and Ψ I (τ I j − τ R i − τ ) is the probability that the time elapsed between the n-th fire and the n + 1-th fire is greater than τ ji − τ R i . As computed in Eq. (4) of the main text, the probability that the time elapsed since j became infected is equal to τ I j is simply The survival probability Ψ R (τ I j ) of recovery events can be written as where ω(u) is the inverse Laplace transform of Ψ R (τ I j ). In Laplace space, the probability distribution P n (τ ) has a convenient form, P n (u) = ψ I (u) n , where ψ I (u) is the Laplace transform of ψ I (t). By inserting Eqs. (26) and (27) into Eq. (24) we obtain By inserting the form of φ(τ ji ) into Eq. (10) of the main text, we obtain an expression for the infection rate λ ef f where λ I is the infection hazard rate, Φ R (τ R i ) and Ψ I (t) are the survival probabilities of φ R (τ R i ) and ψ I (t), respectively and φ R * Ψ I is the convolution between φ R (τ R i ) and Ψ I (t). At this point, some ansatz regarding the form of φ R (τ R i ), the probability that the time elapsed since i recovered is equal to τ R i , is needed to continue. We note that if one does not consider the state of node i in the probability φ R (τ ji ), which corresponds to inserting φ R (τ R i ) = δ(τ R i ) into Eq. (28), one obtains In a, node j has attempted to infect i at least once at time t − τji. This implies that node i must necessarily be infected at that 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.
This effective infection rate λ mf , already found in Cator et al. [19] by using a mean field approximation, is now obtained within a more general formalism. A different possibility is to consider φ R (τ R i ) equal to an exponential distribution, 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 λ ef f . One can see that the exponential ansatz for the form of φ R (τ R i ) is correct for large values of 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 ψ(τ ji ) ≡ lim t→∞ p(τ ji , n i = 0; t|n j = 1), which is the join probability that, given that node j is infected at the observation time, node i is susceptible and the time elapsed since the last infection attempt from j to i is equal to τ ji . Our approximation consists of 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 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 − τ ji , node i either was already infected (in which case the infection attempt is 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 − τ ji and, thus, the probability that node i recovers before the observation time t is 1 − Ψ R (τ ji ). If 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 − Ψ R (τ ji ). If node j does not attempt to infect node i, we cannot know for certain the state of node i. However, given that node j became infected at time t − τ I j , one of his 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 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 j , and the probability that node i is then susceptible at the observation time t is 1 − Ψ R (τ I j ). With probability 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 − Ψ R (τ ji ). Instead, if node j does not attempt to infect node i, this probability reads (1 − Ψ R (τ I j ))/k + (k − 1)/k = (k − Ψ R (τ I j ))/k. In the limit of low prevalence, we expect this approximation to be exact. In the following, we also approximate the value of k by the average degree, k .
The probability density ψ(τ ji ) can be written as where again φ I (τ I j ) is the probability that the time elapsed since j became infected is equal to τ I j . The conditional probability ψ(τ ji |τ I j ) is 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 j , the probability ψ(τ ji ) reads If we restrict to the case of Markovian recovery, we can use its memoryless property, Ψ R (τ ji + τ ) = Ψ R (τ ji )Ψ R (τ ). By using the convenient Laplacian form of the probability distribution P n (τ ), one can obtain The normalization of ψ(τ ji ) reads therefore, by inserting Eq. (35) and Eq. (36) into Eq. (11) of the main text, one finally obtains the approximate infection rate λ app presented in Eq. (12) of the main text. 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 λ mf given by Eq. (30) for two extreme values of the exponent α I controlling the interevent time infection distribution, α I = 0.25 and α I = 10. One can see that the curves for λ mf do not collapse onto one another, especially for the lattice and SF network substrate.

NUMERICAL SIMULATIONS OF THE NON-MARKOVIAN SIS DYNAMICS
To check the validity of the effective infection rates λ ef f and λ app , we run extensive numerical simulations of the non-Markovian SIS dynamics. For each value of the average infection time t I and fixed average recovery time t R =δ −1 = 1, we simulate the non-Markovian SIS dynamics by implementing an algorithm based on a queue of infection and recovery events. At time t = 0, all nodes are in a susceptible state, and a set of f N 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 R from the distribution ψ R (t), and pushes his recovery event at time t + t R to the queue. He also starts k independent infection processes to his k neighbors. In each infection process to a neighbor j, an infection event from node i to node j is scheduled at time t + t 1 I , where t 1 I is randomly extracted from the distribution ψ I (t), only if t 1 I < t R , that is if node i is still infected at time t + t 1 I . A second infection event from node i to node j is scheduled at time t + t 1 I + t 2 I , where t 2 I is randomly extracted from the distribution ψ I (t), only if t 1 I + t 2 I < t R , and so on until n (with possibly n = 0) infection events are generated and pushed to the queue. 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 max , set conveniently. In order to measure the prevalence in the steady state ρ st and the effective infection rate λ ef f , we sample N s = 10 4 time instants uniformly chosen in [T min , T max ], with T min chosen such that the stationary state is reached long before it. For each time instant, we measure the prevalence and the values of τ ij for each active link between nodes i and j, so as to calculate λ ef f by 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 α I = 0.5 and α I = 2, and two different network substrates, 2D lattice and 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 T and average square lifetime T 2 over a number of runs N run , as a function of the average infection time t I (corresponding to an effective rate λ app , hereafter λ for brevity) close to the epidemic threshold, for different sizes N . We set N run = 10 5 for lattice, N run = 10 6 for RDR networks. For each value of N , λ, α I and network substrate we fit the curves of T and T 2 to find the peaks T p and T 2 p and their corresponding values of λ 1 p and λ 2 p . We set λ p as the average of λ 1 p and λ 2 p , provided that λ p falls within the λ 1 p and λ 2 p standard errors. The corresponding endemic probability at the peak P p is interpolated from the data.
The set of equations we used to evaluate the critical point λ c and critical exponents β, δ, ν ⊥ are We first evaluate the critical threshold λ c by plotting the endemic probability P (λ, N ) as a function of N , for several values of λ close to λ c , see the first row of Fig. 5. Through Eq (37), we estimate the value of λ c to be the one which produces the best fit of P (λ c , N ) as a power-law. We then plot the endemic probability at the peak P p as a function of the size N , see the second row of Fig. 5, and the difference |λ c − λ p | as a function of N , see the third row of Fig. 5. We estimate ν ⊥ by means of Eq (39). By means of Equations (37) and (38) we estimate β/ν ⊥ equal to the average of the fits of P (λ c , N ) and P p (N ), provided that β/ν ⊥ falls within the standard errors of the two fits, and so we calculate β, knowing the value of ν ⊥ . Finally, we plot the height of the peaks T p and T 2 p as a function of N , see the fourth row of Fig. 5. We estimate γ 2 and γ 1 for lattices and γ 2 for RDR networks by means of Eq (40), knowing the value of ν ⊥ . For lattices, we calculate δ by means of the equivalence γ n ∼ n − δ while for RDR networks we first check that T p diverges logarithmically as a function of N and then we calculate δ by means of the equivalence γ n = n − δ. The results are reported in Table 1 of the main text.