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.
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 Haight (1967). 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 Barrat et al. (2008) 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 Barabasi (2005); Oliveira and Barabasi (2005); Gonzalez et al. (2008), to natural phenomena Corral (2004); Wheatland et al. (1998), biological processes Kemuriyama et al. (2010), and biochemical reactions Bratsun et al. (2005). These empirical observations have revealed correlated sequences of events with heavy tailed interevent time distributions Karsai et al. (2012), 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 Moinet et al. (2015); Karsai et al. (2014); García-Pérez et al. (2015); Jo et al. (2014); Kiss et al. (2015) and numerical simulation Van Mieghem and van de Bovenkamp (2013); Boguñá et al. (2014). Particular attention has been devoted to epidemic spreading on complex networks, representing the diffusion of information or disease in a population Pastor-Satorras et al. (2015). Recently, it has been shown that a non-Markovian infection dynamics dramatically alter the Susceptible-Infected-Susceptible (SIS) spreading process Van Mieghem and van de Bovenkamp (2013); Cator et al. (2013). 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 Iribarren and Moro (2009, 2011) and retweets on Twitter Lerman (2016). 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 BAILEY (1956); Eichner and Dietz (2003); NISHIURA and EICHNER (2006); Chowell and Nishiura (2014).
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 , 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 , allowing us to compute it by means of numerical simulations, and enables us to derive an approximate analytic expression , in very good agreement with . The approximate value is expected to converge to 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 , 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 , with , 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 distributed as , that is, recovery from the illness does not confer any long lasting immunity, a characteristic present in some sexually transmitted diseases Keeling and Rohani (2007). 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 . In an active link isolated from the rest of the system, the susceptible node becomes infected after a random time has elapsed since the infection was initiated, with distributed as . If a susceptible node is connected to more than one infected neighbor, infection processes take place independently along each infectious link.
Distributions and 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 since the process was initiated, it takes place in the time interval between and Cox (1967). The recovery and infection hazard rates are defined as and , where and are the corresponding survival probabilities, that is, the probability that a given event takes a time longer than . 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 , ; defined as if node is infected at time 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 is infected and accounts for its recovery during the time interval . To achieve this, the stochastic process is defined to be equal to zero with probability and one otherwise, where is the time elapsed, at time , since node became infected. Similarly, the second term in the sum of the right hand side of Eq. (1) accounts for the infection of susceptible node by one of its infected neighbors during the time interval . The the stochastic process is defined to be equal to one with probability and zero otherwise, where is the time elapsed since the infection process of node to node 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 and are themselves stochastic processes.
The average of Eq. (1), first conditioned to the knowledge of the stochastic processes at time , and then over the unconditional values, allows us to write the following differential equation for the probability of node to be infected at time , 111Note that with this definition, the prevalence of the disease at time is simply given by .
The first term in Eq. (2) can be rewritten as (see Supplementary Information, SI)
In the limit , the only information we have about , given that node is infected, is that the recovery time of node after infection is longer than . This implies that the probability density of is given by , where is the average recovery time Cox (1967). By combining this result with the form of the recovery hazard rate, we can write
where we have defined . 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 depends on the evolution of two-point correlation functions, , 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 -point correlation function (see SI):
where we omit the dependence on for brevity and we define the sets of nodes and . Eq. (6) can be written as
where and are the and -point correlation functions of the sets and , respectively, and where we have also defined
Equations (7), (8), and (9) are the central result of our paper as they fully describe the dynamics of the epidemic. However, what makes our formulation interesting is the fact that all the non-Markovian effects of the dynamics are encoded in the terms and . As we shall show later, under certain conditions these parameters take constant values independent of the nodes, that is, and . 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 and a susceptible node depends on the state of node alone, i.e. when a node 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 . One can think of this process as a series of firing events, separated by random times , starting when node becomes infected, so that when one such event takes place at a time that neighbor node is susceptible, node 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 , distributed as .
Within this framework, the average of conditioned to the state of the system can be derived by noting that, at time , the time elapsed since the infection process of node to node started, , is not truly independent of the state of . That is, if the time elapsed since node has recovered is , then it holds that . Therefore, depends on the state of but not on the state of any other neighbor and, consequently, (see SI for a detailed proof). This implies that we can then define an effective infection rate as
where is the probability density of , where is the time elapsed since the start of the infection process from node to node , given that node is susceptible and node is infected and is averaged over all active links in the network.
Concerning the average of , we note that in general, in the long time limit, Eq. (8) does not reduce to Eq. (4), since . This is due to the fact that, especially for low-degree nodes, the time may depend on the status of node ’s neighbors: if at time node is infected and connected only to node (=1), node must have been infected by node , therefore has to be larger than the time elapsed since the infection process of node started. Thus, if the recovery process follows a non-Markovian dynamics it is not possible to define an effective parameter . Therefore, hereafter we consider only the case of Markovian recovery, which implies , so that the non-Markovian SIS dynamics can be reduced to a Markovian one with parameters and .
Although the probability density 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 . To do so, we first notice that can also be written as
where is the joint probability that node is susceptible and the time elapsed since the last infection attempt from to is equal to , given that node is infected at a given observation time , and where is the normalization factor . In the SI, we derive an approximate analytic expression for that, combined with Eq. (11), allows us to derive the following expression for the approximate effective rate
where and are the Laplace transforms of and , respectively, and is the average degree of the network substrate.
We check the validity of the effective infection rates, and , 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,
with parameter controlling the power-law start and tail of the infection inter-event time distribution. We choose , so that is the average infection time. Hereafter, and without loss of generality, we set the time scale to . Once the system has reached its steady state, we evaluate by selecting random time instants along the process. For each time instant, we select all active links, measure the corresponding values of , and calculate as the average of the hazard rates . The approximate infection rate is calculated from Eq. (12) by integrating numerically the Laplace transforms and with .
Figure 2 shows the prevalence at the steady state as a function of and , for different values of and different network substrates: a two dimensional lattice with periodic boundary conditions, an Erdős Rényi (ER) graph with , a random degree regular (RDR) network with , and a scale-free (SF) network with exponent . 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 . This result is particularly noteworthy since two infection processes with the same average infection time but different forms of are known to behave very differently Van Mieghem and van de Bovenkamp (2013), showing huge differences in the prevalence for the same average infection time. This is particularly true in the case of highly heterogeneous processes, such as the one controlled by , with a very skewed form of the inter-event time distribution , and by , which corresponds to an almost-periodic process.
The curves plotted as functions of the approximate infection rate are also almost indistinguishable from the others, showing that is a very accurate approximation of the exact effective rate, for every underlying network topology. In the SI we also show that is considerably different from the mean-field approximation proposed in Van Mieghem and van de Bovenkamp (2013), and far more accurate in describing extreme cases such as and , see Supplementary Fig. 4. Interestingly, as we show in the SI, Eq. (12) is expected to converge to in the limit of low prevalence and, thus, close to the epidemic threshold, . This implies that the exact critical point 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 , evaluated by means of Eq. (12), and the set of critical exponents , and for a non-Markovian SIS dynamics with and , on top of two-dimensional lattice and degree regular networks, by means of the lifespan method proposed in Ref. Boguñá et al. (2013), see SI and Supplementary Fig. 5. Table 1 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 dynamics 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.
Acknowledgements.We acknowledge support from the James S. McDonnell Foundation; the ICREA Academia prize, funded by the Generalitat de Catalunya; the MINECO projects no. FIS2013-47282-C2-1-P and FIS2016-76830-C2-2-P; Generalitat de Catalunya grant no. 2014SGR608; and Science Foundation Ireland grant no. 11/PI/1026.
- Haight (1967) F. A. Haight, “Handbook of the poisson distribution,” (1967).
- Barrat et al. (2008) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, Cambridge, 2008).
- Barabasi (2005) A.-L. Barabasi, Nature 435, 207 (2005).
- Oliveira and Barabasi (2005) J. G. Oliveira and A.-L. Barabasi, Nature 437, 1251 (2005).
- Gonzalez et al. (2008) M. C. Gonzalez, C. A. Hidalgo, and A.-L. Barabasi, Nature 453, 779 (2008).
- Corral (2004) A. Corral, Phys. Rev. Lett. 92, 108501 (2004).
- Wheatland et al. (1998) M. S. Wheatland, P. A. Sturrock, and J. M. McTiernan, The Astrophysical Journal 509, 448+ (1998).
- Kemuriyama et al. (2010) T. Kemuriyama, H. Ohta, Y. Sato, S. Maruyama, M. Tandai-Hiruma, K. Kato, and Y. Nishida, Biosystems 101, 144 (2010).
- Bratsun et al. (2005) 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 .
- Karsai et al. (2012) M. Karsai, K. Kaski, A.-L. Barabási, and J. Kertész, Sci. Rep. 2 (2012).
- Moinet et al. (2015) A. Moinet, M. Starnini, and R. Pastor-Satorras, Phys. Rev. Lett. 114, 108701 (2015).
- Karsai et al. (2014) M. Karsai, N. Perra, and A. Vespignani, Sci. Rep. 4, srep04001 (2014).
- García-Pérez et al. (2015) G. García-Pérez, M. Boguñá, and M. Á. Serrano, Scientific Reports 5, 9714 EP (2015).
- Jo et al. (2014) H.-H. Jo, J. I. Perotti, K. Kaski, and J. Kertész, Physical Review X 4, 011041 (2014).
- Kiss et al. (2015) I. Z. Kiss, G. Röst, and Z. Vizi, Phys. Rev. Lett. 115, 078701 (2015).
- Van Mieghem and van de Bovenkamp (2013) P. Van Mieghem and R. van de Bovenkamp, Physical Review Letters 110 (2013), 10.1103/PhysRevLett.110.108701.
- Boguñá et al. (2014) M. Boguñá, L. F. Lafuerza, R. Toral, and M. A. Serrano, Phys. Rev. E 90, 042108 (2014).
- Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
- Cator et al. (2013) E. Cator, R. van de Bovenkamp, and P. Van Mieghem, Physical Review E 87, 062816 (2013).
- Iribarren and Moro (2009) J. L. Iribarren and E. Moro, Physical review letters 103, 038702 (2009).
- Iribarren and Moro (2011) J. L. Iribarren and E. Moro, Physical Review E 84, 046116 (2011).
- Lerman (2016) K. Lerman, Future Internet 8, 21 (2016).
- BAILEY (1956) N. T. J. BAILEY, Biometrika 43, 15 (1956), http://biomet.oxfordjournals.org/content/43/1-2/15.full.pdf+html .
- Eichner and Dietz (2003) M. Eichner and K. Dietz, American Journal of Epidemiology 158, 110 (2003), http://aje.oxfordjournals.org/content/158/2/110.full.pdf+html .
- NISHIURA and EICHNER (2006) H. NISHIURA and M. EICHNER, Epidemiology and Infection 135, 1145 (2006).
- Chowell and Nishiura (2014) G. Chowell and H. Nishiura, BMC Medicine 12, 1 (2014).
- Keeling and Rohani (2007) M. Keeling and P. Rohani, Modeling Infectious Diseases in Humans and Animals (Princeton University Press, 2007).
- Cox (1967) D. R. Cox, Renewal Theory (Methuen, London, 1967).
- (29) Note that with this definition, the prevalence of the disease at time is simply given by .
- Boguñá et al. (2013) M. Boguñá, C. Castellano, and R. Pastor-Satorras, Phys. Rev. Lett. 111, 068701 (2013).
Appendix A Supplementary Information
Appendix B Derivation of the differential equation for the point correlation function
The average of Eq (1) of the main text given the state of the system at time , , can be written as
Let us consider a set of nodes . The correlation function between these nodes reads
where the outer average is over the state of the system at time . 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 for different nodes. The first term in Eq (15) can be written, by means of Eq (14) as
where is the set of all subsets of containing nodes. Because of the term , however, the expansion to linear order in is a reduced sum over , and thus
Therefore the point correlation function reads
from which Eq (6) of the main text follows immediately.
Appendix C Time of active link does not depend on the states of other nodes different from and
A critical step in our approach is to prove that
The probability in the left hand side of this equation can be written as
where is the joint probability density, at time , that given that node is susceptible and nodes and are infected, the time elapsed since recovered is and the times elapsed since and became infected are and , respectively. By Bayes’ rule, is the probability density of the time conditioned on the times . However, it is easy to see that since infection events take place in active links independently, once and are fixed, is totally independent of the elapsed times since nodes other than became infected. Therefore,
which directly gives the result in Eq. (20).
Appendix D General formalism for
Using the result in Eq. (22), at the steady state the probability density of the time elapsed since the infection process of node to node started, given that node is susceptible and node is infected, can be written in general as
where is the joint probability that the time elapsed since became infected is equal to and the time elapsed since recovered is equal to . If we assume that the two process are uncorrelated, can be factorized into , and Eq. (23) reduces to
where is the Heaviside step function, is the probability that the time elapsed since became infected is equal to and is the probability that the time elapsed since recovered is equal to . The conditional probability is simply , and
where is the number of infection attempts of node to node , is the probability that the time elapsed since node became infected and the moment of his -th fire is equal to , and is the probability that the time elapsed between the -th fire and the -th fire is greater than . As computed in Eq. (4) of the main text, the probability that the time elapsed since became infected is equal to is simply
The survival probability of recovery events can be written as
where is the inverse Laplace transform of . In Laplace space, the probability distribution has a convenient form, , where is the Laplace transform of . By inserting Eqs. (26) and (27) into Eq. (24) we obtain
By inserting the form of into Eq. (10) of the main text, we obtain an expression for the infection rate
where is the infection hazard rate, and are the survival probabilities of and , respectively and is the convolution between and . At this point, some ansatz regarding the form of , the probability that the time elapsed since recovered is equal to , is needed to continue. We note that if one does not consider the state of node in the probability , which corresponds to inserting into Eq. (28), one obtains
This effective infection rate , already found in Cator et al. Cator et al. (2013) by using a mean field approximation, is now obtained within a more general formalism. A different possibility is to consider equal to an exponential distribution,
with rate . The rate can be written as a simple function of the prevalence , . 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 is correct for large values of the prevalence , but it fails for low prevalence, thus close to the epidemic threshold.
Appendix E Approximations to
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 , which is the join probability that, given that node is infected at the observation time, node is susceptible and the time elapsed since the last infection attempt from to is equal to . Our approximation consists of estimating the probability that node is susceptible at the observation time , which depends on the time instant at which node became infected, this time instant being unknown in principle.
We first consider the case of low prevalence, . Then, let us consider separately the cases in which node attempts at least once to infect node , , and the case of no attempts, . In the first case, at the time of the last infection event from to , , node 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 is in an infected state at time and, thus, the probability that node recovers before the observation time is . If the prevalence is low, the probability that node is subsequently infected by one of its neighbor (other than ) 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 is susceptible at time is simply . If node does not attempt to infect node , we cannot know for certain the state of node . However, given that node became infected at time , one of his neighbors must have infected him. Let us consider that node has degree . If the prevalence is low, it is very unlikely to find more than one neighbor of node infected simultaneously and we assume that only one of his neighbors was infected and infected him. With probability , such infected node is node (See Fig. 3b), so that is infected at time , and the probability that node is then susceptible at the observation time is . With probability , the infected node is a neighbor other than and, thus, we assume that node was susceptible at time and it will remain in this state until time with probability equal to one. Summing up, if node attempts at least once to infect node , then the probability that it is susceptible at the observation time is . Instead, if node does not attempt to infect node , this probability reads . In the limit of low prevalence, we expect this approximation to be exact. In the following, we also approximate the value of by the average degree, .
The probability density can be written as
where again is the probability that the time elapsed since became infected is equal to . The conditional probability is
where the first term accounts for the case in which there are no infection attempts from to , , while the second term accounts for the case . By inserting Eq. (26) into Eq. (32) and integrating over , the probability reads
If we restrict to the case of Markovian recovery, we can use its memoryless property, . By using the convenient Laplacian form of the probability distribution , one can obtain
The normalization of reads
In Fig. 4, we show the steady-state prevalence as a function of the approximate effective infection rate , given by Eq. (12) of the main text, and the mean field effective rate given by Eq. (30) for two extreme values of the exponent controlling the interevent time infection distribution, and . One can see that the curves for do not collapse onto one another, especially for the lattice and SF network substrate.
Appendix F 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 non-Markovian SIS dynamics. For each value of the average infection time and fixed average recovery time , we simulate the non-Markovian SIS dynamics by implementing an algorithm based on a queue of infection and recovery events.
At time , all nodes are in a susceptible state, and a set of randomly chosen nodes, with , change their state to the infected one. In the algorithm, whenever a node changes his state from susceptible to infected at time , he first randomly extracts his recovery time from the distribution , and pushes his recovery event at time to the queue. He also starts independent infection processes to his neighbors. In each infection process to a neighbor , an infection event from node to node is scheduled at time , where is randomly extracted from the distribution , only if , that is if node is still infected at time . A second infection event from node to node is scheduled at time , where is randomly extracted from the distribution , only if , and so on until (with possibly ) 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 , changes his state from infected to susceptible. If the pulled event is an infection event from node to node , and is already in a infected state, nothing happens, otherwise node 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 , set conveniently. In order to measure the prevalence in the steady state and the effective infection rate , we sample time instants uniformly chosen in , with chosen such that the stationary state is reached long before it. For each time instant, we measure the prevalence and the values of for each active link between nodes and , so as to calculate 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 Boguñá et al. (2014), which is much slower, and we obtained identical results for the prevalence and the effective infection rate.
Appendix G 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 and , and two different network substrates, 2D lattice and RDR network. We address the critical properties by means of the lifespan method Boguñá et al. (2013), in which the infection starts with a single infected node. In the lifespan method, each realization is characterized by its lifetime, , and its coverage, , defined as the number of distinct nodes that have become infected at least once. We let each realization run until either the coverage