Physical Review. E
American Physical Society
image
Limitations of discrete-time approaches to continuous-time contagion dynamics
Volume: 94, Issue: 5
DOI 10.1103/PhysRevE.94.052125
Abstract

Continuous-time Markov process models of contagions are widely studied, not least because of their utility in predicting the evolution of real-world contagions and in formulating control measures. It is often the case, however, that discrete-time approaches are employed to analyze such models or to simulate them numerically. In such cases, time is discretized into uniform steps and transition rates between states are replaced by transition probabilities. In this paper, we illustrate potential limitations to this approach. We show how discretizing time leads to a restriction on the values of the model parameters that can accurately be studied. We examine numerical simulation schemes employed in the literature, showing how synchronous-type updating schemes can bias discrete-time formalisms when compared against continuous-time formalisms. Event-based simulations, such as the Gillespie algorithm, are proposed as optimal simulation schemes both in terms of replicating the continuous-time process and computational speed. Finally, we show how discretizing time can affect the value of the epidemic threshold for large values of the infection rate and the recovery rate, even if the ratio between the former and the latter is small.

Fennell, Melnik, and Gleeson: Limitations of discrete-time approaches to continuous-time contagion dynamics

I.. INTRODUCTION

A feature of our environment is the existence of networks, from real-life human contact networks, to virtual networks such as online social networks, to functional and technological networks such as transport networks and the Internet [1]. Networks form a medium for contagions, which spread from node to node through the links of the networks. Contagions can be physical [2,3], cultural [4,5], societal [6–8], financial [9–11], and the modeling of such contagions [12–16]—along with the understanding of the suitability of various modeling approaches [17–19]—is vital for matters of the utmost public importance [20–22].

A common modeling paradigm for studying contagions is the framework of continuous-time Markov processes [23–25], where events (such as the infection of a susceptible individual by an infected individual) occur at certain rates. The most well known of these models are epidemiological compartment models [2], which, although introduced as models of disease spread [26], are also widely used as models of social contagions such as the diffusion of information and innovations [27–29]. Continuous-time Markov process models can provide valuable insights into contagion processes, and have real value in both predicting and controlling contagious outbreaks [30–32].

One avenue to study continuous-time Markov process models is by using discrete-time approximations [32–48]. Such approaches can be either numerical (i.e., synchronous updating Monte Carlo simulations) or theoretical. In a discrete-time approach, time is discretized into time steps of length Δt (which usually takes the value Δt=1), and events occur with certain probabilities. These probabilities are known as the state transition probabilities, and are simply the product of the corresponding rate and the time step Δt.

Although discrete-time approaches correspond to their continuous-time counterpart in the limit Δt0, they can significantly differ in the case that Δt is finite. Allen, in her work [33], shows that discrete-time susceptible-infected-susceptible (SIS) and susceptible-infected-recovered (SIR) models can produce complex behavior such as period doubling and chaotic effects for sufficiently large values of the time step and/or contact rate. This behavior is not possible in the continuous-time SIS and SIR models, and is thus no more than an artifact of discretizing time. Similarly, Gomez et al.  [39] observe that differences between continuous and discrete-time SIS dynamics are substantial when an arbitrary time step of Δt=1 is employed. An understanding of the discrepancies introduced as a result of discretizing time is thus important, allowing us to gauge the validity of discrete-time approaches and when they may accurately be employed.

In this paper, we show the limitations of discrete-time approaches when used to study continuous-time contagion dynamics. Our message is clear–that the accuracy of such methods will be poor if state transition probabilities are too large, leading to deviations from the underlying continuous-time process. The repercussions of this are manifold. Discrete-time theoretical approaches can be significantly inaccurate for large values of the contagion parameter values (such as infection and recovery rates), and thus the analysis of such approaches will not be valid. Furthermore, discrete-time Monte Carlo simulations—often used as a gold standard [34,49,50]—can be inaccurate for large parameter values, and such inaccurate simulations can lead to misleading conclusions. We illustrate this latter point with an example from the literature in Sec. IV. Our work highlights the consequences of erroneous approaches to studying continuous-time contagion dynamics, which has important implications not only for the academic study of these dynamics [34–48] but also for the implementation of such dynamics within large-scale simulators for real contagions [31,32].

II.. CONTINUOUS-TIME CONTAGION DYNAMICS AND THE DISCRETE-TIME APPROXIMATION

To begin, we describe in some detail both continuous and discrete-time Markov processes to illustrate mathematically the difference between the two. In continuous-time Markov processes, events are described by rates λ, while events in the discrete-time analog are described by transition probabilities λ~, where λ~=λΔt. In the course of our analysis we focus on the specific example of SIS dynamics; however, our analysis holds for any continuous-time Markovian dynamics, with the core message being the limitations on the size of the transition probabilities λ~ for which discrete-time approaches are accurate.

A.. Continuous-time SIS dynamics

Consider SIS dynamics taking place on a network of N nodes. This is a continuous-time Markov process where at any time t each node i in the network has a corresponding state Xti, which is either susceptible (Xti=S) or infected (Xti=I[23,51,52]. The states of each node in the network change dynamically over time. Susceptible nodes become infected through each of their infected neighbours at a rate β per infected neighbor, while infected nodes recover at a rate μ. “Rate” here refers to instantaneous transition rates, which in continuous-time dynamics define the transitions between states; these are defined in terms of probabilities as

μ=limΔt0PXt+Δti=S|Xti=IΔt,
β=limΔt0PXt+Δti=Iviaj|Xti=S,Xtj=IΔt,
where {Xt+Δti=Iviaj} is the event that a susceptible node i became infected through an infected neighbor j. The fraction terms in the right hand sides of Eqs. (1) and (2) are the probabilities of state changes per unit time and taking the limit of these fractions as Δt0 leads to the concept of transition rates. In general, we can define ri as the rate at which a node i changes from its current state to the opposite state; this is given by
ri=βmi,tifXti=SμifXti=I,
where mi,t is the number of infected neighbors of node i at time t.

The evolution of the dynamics in the network can be fully described by the master equation for the Markov process [24,53]. If we denote by Yt={Xti}i=1N the state of the network at time t, and by p(y,t) the probability that the network is in state Yt=y, then the master equation is given by

ddtp(y,t)=yy(p(y,t)ryyp(y,t)ryy),
with initial conditions p(y,0)=p0(y). Here ryy is the instantaneous rate at which the network changes from state y to y and is fully determined by the network structure and the transition rates μ and β. While the master equation is the gold standard—exactly describing the evolution of SIS dynamics—the dimension of its sample space Ω is 2N, which in general is prohibitively large for analytical or numerical studies. One way to tackle this problem is to study the dynamics as a series of individual transitions between states. In continuous-time dynamics nodes change state one at a time, or asynchronously [54]. Given the state of the network, the probability distributions governing both the length of time until the next state change and the node which will change state can be constructed. These are given by the following lemmas (of which rigorous derivations can be found in the literature, e.g., Ref. [23], chapter 10):

Lemma 1. Let τ be the holding time of the network, the length of time that the network remains in its current state before changing to the next state. Then τ is an exponentially distributed random variable and the parameter of the distribution is the sum of the individual node transition rates, i.e., i=1Nri.

Lemma 2. The probability that the next node in the network to change state will be node i is ri/j=1Nrj.

Lemmas 1 and 2 describe how the network probabilistically evolves from one state to another. They are the basis of continuous-time stochastic simulation methods such as the well-known Gillespie algorithm, also known as the stochastic simulation algorithm or kinetic Monte Carlo [55–58]. Such simulations are often referred to as event-based simulations because the time intervals are not fixed but rather correspond to the time between consecutive state-changes in the system. At each step in such algorithms, time advances by an amount τ and node i changes its state, where τ and i are random numbers drawn according to Lemmas 1 and 2 (Fig. 1). Stochastic simulations give the opportunity to construct p(y,t) empirically by running multiple realizations of the stochastic process and aggregating over an ensemble of realizations. Such simulations are statistically exact as they are fully based on Lemmas 1 and 2, which are derived without approximation from the axioms of the Markov process.

Schematics of both (a) the Gillespie algorithm and (b) the synchronous updating scheme. Vertical ticks on the t axis indicate the moments through which the simulation advances—in synchronous updating the interval between these moments is a fixed time step with value Δt while in the Gillespie algorithm the interval is a random variable τ given by Lemma 1. The light green and dark red circles are nodes in the network, which are in the susceptible and infected states, respectively. A square around a node means that the node has been chosen for updating at a certain moment and may change its state. In the Gillespie algorithm a node is chosen according to Lemma 2 and will always change its state while in the synchronous updating scheme every node has the chance to change state and will do so with a probability that depends on their state and the states of their nearest neighbors.
FIG. 1.
Schematics of both (a) the Gillespie algorithm and (b) the synchronous updating scheme. Vertical ticks on the t axis indicate the moments through which the simulation advances—in synchronous updating the interval between these moments is a fixed time step with value Δt while in the Gillespie algorithm the interval is a random variable τ given by Lemma 1. The light green and dark red circles are nodes in the network, which are in the susceptible and infected states, respectively. A square around a node means that the node has been chosen for updating at a certain moment and may change its state. In the Gillespie algorithm a node is chosen according to Lemma 2 and will always change its state while in the synchronous updating scheme every node has the chance to change state and will do so with a probability that depends on their state and the states of their nearest neighbors.

B.. Discrete-time approach

In a discrete-time framework, time is no longer treated as a continuous variable but rather takes the form of a discrete variable, which advances in time intervals of length Δt. Instantaneous transition rates are replaced by transition probabilities. In a single time interval, susceptible nodes become infected through their infected neighbors with probability β~=βΔt per infected neighbor, while infected nodes recover with probability μ~=μΔt. Note that Δt is often assumed to take the value Δt=1, but even in this case it should be included in the expression for β~ and μ~ to clarify that a rate needs to be multiplied by a time step before it can be expressed as a probability.

The discretization of time in this manner leads to two deviations from the continuous-time process. These deviations arise through both the transition probabilities, which are used in place of transition rates, as well as the parallel (synchronous) state changes in discrete-time systems that are uncharacteristic of continuous-time dynamics. To understand the roots of the deviations introduced through the transition probabilities we can examine the definitions of μ and β as rates given in Eqs. (1) and (2). These equations can be rearranged to give transition probabilities in terms of these rates, i.e.,

PXt+Δti=S|Xti=I=μΔt
PXt+Δti=Iviaj|Xti=S,Xtj=I=βΔt,
where in this case Δt is an infinitesimally small length of time. In the case that Δt is not infinitesimally small, Eqs. (5) and (6) become approximations. In a time interval of length Δt in the continuous-time Markov process, the exact probability that an infected node will recover is 1eμΔt while the probability that a susceptible node will become infected by a given infected neighbor is 1eβΔt. The transition probabilities μ~ and β~, the right-hand side of Eqs. (5) and (6), are approximations to 1eμΔt and 1eβΔt, respectively, and an important question then arises of the effect these approximations have on the dynamics.

Figure 2 shows the actual probability 1eλΔt along with the discrete-time probability λ~=λΔt, where we use the parameter λ to represent either μ or β. We also plot the error ε, which is defined as the difference between the discrete-time probability and the actual probability. When λ~<0.1,ε<0.01 and so the approximation is fairly accurate in this range. For larger values of the state transition probability λ~, however, the approximation differs significantly from the true values. At λ~=0.5,ε0.1 and when λ~=1,ε0.37. These individual errors can accumulate and have significant implications on the dynamics as a whole; indeed we show empirically in Secs. III and IV that although discrete-time approaches can be very accurate when μ~ and β~ are very small they begin to lose accuracy when μ~ and β~ are of the order of magnitude of 101.

The actual probability (blue solid line) that a rate λ event will occur in a time step of length Δt plotted along with the approximate probability λ~ (black dash-dotted line) as used in discrete-time formalisms. The error ε is defined as the absolute distance between the two.
FIG. 2.
The actual probability (blue solid line) that a rate λ event will occur in a time step of length Δt plotted along with the approximate probability λ~ (black dash-dotted line) as used in discrete-time formalisms. The error ε is defined as the absolute distance between the two.

Second, we comment on the synchronous updating nature of discrete-time approaches. This is in contrast to the continuous-time process where nodes change state asynchronously and the change of state of one node immediately affects the transition rates of the other nodes (Fig. 1). The strength of effect will depend on the transition probabilities, as the values μ~ and β~ dictate the number of state changes that take place in each time step and thus the propensity of multiple nodes to change state at the same time.

Thus, we arrive at a simple conclusion: the values of μ~ and β~ (and thus μ,β, and Δt) used in discrete-time approaches should be controlled so that these approaches are accurate representations of the continuous-time process. For large values of μ or β, the time step Δt should be small while if Δt=1, as in the case of the majority of discrete-time approaches, the values of μ and β should be relatively small. Throughout the rest of this paper we will give empirical evidence of this conclusion.

Finally, we comment on discrete-time numerical simulation schemes that are used to stochastically simulate SIS dynamics. A commonly used simulations scheme is synchronous updating, also referred to as rejection sampling (Fig. 1[34,49,57,59]. In this case, time advances in steps of one time unit, i.e., Δt=1. In a single time unit, a susceptible node will become infected by its infected neighbors with probability β~ per infected neighbor while infected nodes become susceptible with probability μ~. Synchronous updating simulations are statistically exact realizations of the discrete-time dynamics; these dynamics are fully described by the discrete-time master equation

p(y,t+1)=yp(y,t)qyy,
where qyy is the probability that the network changes from state y to state y in a time step of length Δt=1 and is fully determined by the network structure and the transition probabilities μ~ and β~ [24]. Because synchronous updating simulations exactly mimic the discrete-time dynamics and master equation they will be used throughout this paper to gauge the accuracy of the discrete-time approach.

In the remainder of the paper, we show how the approximations introduced in discrete-time approaches can lead to misrepresentation of the actual continuous-time dynamics. We begin in the next section by examining the discrete-time approximations of Eqs. (5) and (6) for fixed μ and β and various values of Δt. We show that discrete-time dynamics can accurately reproduce continuous-time dynamics for small values of Δt, but that they incur a breakdown in accuracy as Δt increases. Further to this, we show in Sec. IV that when the time step is fixed to the value Δt=1, as in much of the literature, discrete-time approaches break down in accuracy when the transition rates (μ and β) are too large. This limits the range of parameters that can be studied with discrete-time approaches. We illustrate this with an example from the literature, also showing how synchronous updating simulation schemes can favour discrete-time formalisms leading to biased conclusions when comparing against continuous-time theories. Finally, in Sec. V we show that overly large values of β~ and μ~ can affect the value of the epidemic threshold, even if the effective transition rate defined as γ=β/μ=β~/μ~ is small.

III.. EFFECT OF THE TIME STEP Δt WHEN DISCRETIZING TIME

In this section, we analyze the discrete-time approximations introduced in Sec. II B as a function of the size of the discrete time step Δt. We do this by carrying out synchronous updating simulations for various values of Δt and comparing them against exact results obtained from the master equation. Numerical simulations are carried out in c ++ and the code is available online [60].

As our example, we consider SIS dynamics on a complete graph of N nodes, i.e., a graph where every pair of nodes is connected. On such a graph, the SIS dynamics are defined by the rate functions

ri=βZtifXti=SμifXti=I,
where Zt is the number of infected nodes at time t and β and μ are the infection rate and recovery rate respectively, consistent with Eq. (3) for the complete graph. We choose the complete graph because on such a graph the master equation given in Eq. (4) can be reduced from a system of 2N equations for p(y,t) to a system of N+1 equations for p(n,t), the probability that there are n infected nodes in the graph at time t [24]. This reduced system is given by
ddtp(n,t)=(μn+βn(Nn))p(n,t)+μ(n+1)p(n+1,t)+β(n1)×(Nn+1)p(n1,t),
for 0nN with initial conditions p(n,0)=p0(n). For small values of N, this system can easily be solved using standard differential-equation solvers, giving us a gold standard against which to compare the discrete-time simulations. We also perform Gillespie algorithm simulations to illustrate the accuracy and the speed of such simulations and thus their efficacy in simulating continuous-time dynamics.

We present the results for SIS dynamics with β=1,μ=1 running on a complete graph with N=10 nodes in Fig. 3. We plot the solution of Eq. (9) as well as the numerical results given by the Gillespie algorithm and synchronous updating schemes with different time steps Δt. For the numerical simulations, we performed 106 realizations and obtained the corresponding p(n,t) by taking the fraction of realizations in which there are n infected nodes at time t. For the synchronous updating simulations, we consider time steps of Δt=0.01,0.1 and 1. From Fig. 2, it is clear that these values of Δt with μ=β=1 will give a comprehensive range on which to judge the accuracy of the discrete-time approach, while noting that for Δt=1 and for these μ and β parameter values the system is deterministic.

SIS dynamics with β=μ=1 on a complete graph with N=10 nodes. (a) Time evolution of the expected fraction of infected nodes ηt for both the exact master equation Eq. (9) (solid line) and synchronous updating simulations with time steps Δt=0.1 and Δt=1 (dashed lines). The Gillespie algorithm and Exact curves are indistinguishable. (b) and (c) Histograms showing the exact probability mass function p(z,t) at t=3—calculated from numerical integration of Eq. (9)—and the probability mass functions obtained empirically from 106 simulation realizations for both (b) synchronous updating simulations with time steps Δt=0.01,Δt=0.1, and Δt=1 and (c) the Gillespie algorithm, respectively.
FIG. 3.
SIS dynamics with β=μ=1 on a complete graph with N=10 nodes. (a) Time evolution of the expected fraction of infected nodes ηt for both the exact master equation Eq. (9) (solid line) and synchronous updating simulations with time steps Δt=0.1 and Δt=1 (dashed lines). The Gillespie algorithm and Exact curves are indistinguishable. (b) and (c) Histograms showing the exact probability mass function p(z,t) at t=3—calculated from numerical integration of Eq. (9)—and the probability mass functions obtained empirically from 106 simulation realizations for both (b) synchronous updating simulations with time steps Δt=0.01,Δt=0.1, and Δt=1 and (c) the Gillespie algorithm, respectively.

We consider the SIS process at time t=3 at which stage the expected fraction of infected nodes ηt=n=010np(n,t) has reached a metastable state (Fig. 3[61]. At t=3 we empirically construct p(n,3) from the synchronous updating simulations and compare it to p(n,3) calculated from the master equation (9). The histogram of Fig. 3(b) shows this comparison. From this histogram, it is clear that while the discrete-time simulations are quite accurate for small Δt this accuracy can fully break down when Δt is too large. The accuracy of the probability distribution in the metastable state depends highly on the value of the time step used to reach the metastable state. In the synchronous updating simulations with Δt=1 the results are highly inaccurate with all of the probability concentrated on n=9, i.e., p(n,3)=δn,9. Even the case Δt=0.1, while fairly accurate, shows discrepancies in both the probability distribution p(n,3) and the expected fraction of infected nodes ηt [Fig. 3(a)]. Considering that the error ε between μ~ (β~) and 1eμΔt (1eβΔt) is less than 0.005 for μ=β=1 and Δt=0.1 (Fig. 2), we conclude that these discrepancies are due to the simultaneous state changes in synchronous updating, which are uncharacteristic of the continuous-time process.

In the histogram of Fig. 3(c) we compare p(n,3) constructed empirically from the Gillespie algorithm to p(n,3) calculated from the master equation. The Gillespie algorithm is extremely accurate and matches the exact p(n,t) to a high degree of precision. Furthermore, this algorithm is computationally rapid. We performed a short comparison of the simulation algorithms in terms of speed, showing in Table I the simulation run times for the 106 realizations for the Gillespie algorithm and for synchronous updating with various values of Δt. For Δt=0.01—corresponding to the simulations which most closely match the accuracy of the Gillespie simulations—the Gillespie algorithm is an order of magnitude faster. This computational speed, along with its natural precision of the algorithm, make the Gillespie algorithm an optimal algorithm for simulating continuous-time dynamics.

TABLE I.
Time T (in seconds) taken to carry out the Gillespie algorithm and synchronous updating numerical simulations for the examples described in Sec. III and Sec. IV. The simulation code was written in c++ and the simulations were run on a single-CPU contemporary desktop computer.
Complete network (Sec. III)
GillespieSynchronous
ΔtTΔtT
-5.670.0165.62
  0.112.46
  14.48
Erdős-Rényi network (Sec. IV)
GillespieSynchronous
ΔtTΔtT
-42.10.012555.8
  181.4

To summarize, the accuracy of discrete-time approximations to continuous-time dynamics depends highly on the size of the discrete-time step Δt at which the system evolves. This has extremely important implications for real-world predictive models of epidemic spreads that are discrete-time based [31,32], as overly large time steps can affect the prediction of both the expected evolution of a contagion [Fig. 3(a)] as well as variance or confidence intervals around the expected evolution [Figs. 3(b) and 3(c)].

In the next section, we fix the time step at Δt=1 and show how the accuracy breaks down when the infection and recovery rates are too large, showing that discrete-time formalisms using this approach are limited in the ranges of the rate parameters that they can study and thus their ability to match continuous-time dynamics.

IV.. LIMITATIONS ON RANGE OF PARAMETER VALUES WHEN Δt=1

As mentioned in Sec. II B, synchronous updating has the same characteristics of discrete-time systems, which are characterized by transition probabilities and difference equations of the form

p(y,t+Δt)=f({p(y,t)}yΩ),
where p(y,t+Δt)—the probability that the system is in state y at time t+Δt—is a function of the probabilities p(y,t) for all possible states y in the sample space Ω. On the other hand, continuous-time systems are characterized by transition rates and differential equations of the form given by the master equation (4). Although the discrete-time formulation coincides with the continuous-time one in the limit Δt0, the dynamics will differ for noninfinitesimal Δt. Issues then arise when comparing discrete and continuous-time systems and the choice of numerical scheme becomes important. We illustrate this now with an example from the literature, while also showing how the accuracy of discrete-time approaches with Δt=1 can be insufficient for large values of the transition rates.

A prominent current strand of research is the behavior of the SIS model on infinite networks with power-law degree distributions [36,49,62–66]. In Ref. [36], Chakrabarti et al. introduced the nonlinear dynamical systems (NLDS) theory, a discrete-time approach to SIS modeling with a set of mean-field difference equations of the form

pi,t+1=f{pj,t}j=1N,
for 0iN, where pi,t+1 is the probability that node i is infected at time t. They compare their results to two continuous-time formulations, the heterogeneous mean-field (HMF) approach of Pastor-Sattoras and Vespignani [49] and the Kephart-White (KW) [67] approach. The bases of the comparison are synchronous updating numerical simulations with a time step Δt=1 and it is found (see for example Fig. 4 of Ref. [36]) that the NLDS theory is much closer to the numerical simulations than both the HMF and KW theories.
Time evolution of the expected fraction of infected nodes ηt for SIS dynamics with infection rate β=0.2 and recovery rates (a) μ=0.48 and (b) μ=0.72. The network here is an Erdős-Rényi network with 1000 nodes and average degree 〈k〉=4. Each trajectory is averaged over 104 realisations. Synchronous updating with a time step of Δt=1—as in [36]—is given by the circular symbols (SΔt=1). This significantly deviates from Gillespie algorithm simulations (GA, triangular symbols) and synchronous updating simulations with a small time step of Δt=0.01 (SΔt=0.01, square symbols) if the transition rates μ and/or β are too large.
FIG. 4.
Time evolution of the expected fraction of infected nodes ηt for SIS dynamics with infection rate β=0.2 and recovery rates (a) μ=0.48 and (b) μ=0.72. The network here is an Erdős-Rényi network with 1000 nodes and average degree k=4. Each trajectory is averaged over 104 realisations. Synchronous updating with a time step of Δt=1—as in [36]—is given by the circular symbols (SΔt=1). This significantly deviates from Gillespie algorithm simulations (GA, triangular symbols) and synchronous updating simulations with a small time step of Δt=0.01 (SΔt=0.01, square symbols) if the transition rates μ and/or β are too large.

However, the comparison of discrete-time and continuous-time formulations in this manner is biased. Synchronous updating with a time step Δt=1 is the correct procedure for numerically simulating discrete-time dynamics. On the other hand, to simulate continuous-time dynamics either synchronous updating with a vanishingly small time step or a continuous-time simulation scheme such as the Gillespie algorithm should be used.

To illustrate the difference resulting from the use of the different updating methods we reproduce an example from [36]. The example is SIS dynamics on an Erdős-Rényi network of 1000 nodes and mean degree k=4. Figure 4 shows various numerical simulations of these dynamics. Again, the computer code used to perform the simulations is available from Ref. [60]. Included in Fig. 4 are synchronous updating simulations with a time step Δt=1, as in Ref. [36], along with synchronous updating simulations with a small time step Δt=0.01 and Gillespie algorithm simulations. In Fig. 4(c) of Ref. [36], where μ=0.72 and β=0.2, it can be seen that the fraction η¯ of infected nodes in the metastable state given by the NLDS theory matches very closely the synchronous updating numerical simulations. However, as can be seen in Fig. 4(b) here, these synchronous updating simulations differ quite significantly from continuous-time simulations, which plateau at the metastable state with η¯0.04. The KW theory, which in Ref. [36] is rejected as being inaccurate, actually converges to a value much closer to the continuous-time simulations than the NLDS theory. Thus, using the correct simulation technique, the conclusions in Ref. [36] should be reversed: the KW model is more accurate than the NDLS model.

For fixed Δt=1, the accuracy of the discrete-time approach decreases as μ and β increase. In the example above, when μ is decreased from μ=0.72 to μ=0.48 the discrete-time simulations match relatively closer to the continuous-time simulations [Fig. 4(a)], while when μ is decreased further to μ=0.24 the discrepancy between the two simulations in negligible. Chakrabarti et al. state that their model “outperforms (the KW model) when μ is high.” However the opposite is the case: their discrete-time approach breaks down in accuracy (as an approximation to the continuous-time process) as μ increases.

We conclude with an observation to motivate the next section. As μ is increased from μ=0.48 [Fig. 4(a)] to μ=0.72 [Fig. 4(b)], the fraction of infected nodes in the metastable state η¯ (at, for example, t=100) decreases for both the continuous-time simulations and discrete-time simulations. However, η¯ decreases quicker for the continuous-time simulations and so it would seem that the critical value μc at which η¯ first becomes zero will be different depending on whether a discrete-time or continuous-time approach is used. This has implications for the epidemic threshold, which is the focus of the next section.

V.. EFFECT ON EPIDEMIC THRESHOLD

A characteristic of SIS dynamics is the occurrence of phase transitions as the effective transition rate γ is varied. Recall that the effective transition rate is defined as the ratio of the infection rate to the recovery rate, i.e., γ=β/μ. Depending on the structure of the network and whether the network is finite or infinite, the critical point, or epidemic threshold, γc between different phases can vary. As mentioned in Sec. IV, there are still remaining questions about the steady-state behavior of the SIS model—particularly the value of the epidemic threshold on such networks—and so a good understanding of how different approximations affect the value of the epidemic threshold is important. In this section, we show that although the epidemic threshold is defined in terms of the ratio γ=β/μ=β~/μ~, the individual values of the transition probabilities β~ and μ~ used in discrete-time approaches affect the value of the epidemic threshold when it is calculated by (i) performing discrete-time numerical simulations or (ii) iterating a discrete-time system [such as Eq. (10)] from a set of initial conditions (as, for example, in Ref. [37]). Note however that the epidemic threshold predicted by steady-state analysis [i.e., setting pt+1=pt in Eq. (11)], such as in Ref. [36], is completely valid.

We show how μ~ and β~ affect the value of the epidemic threshold in the following manner. For a given network, we fix the value of μ~ and vary β~ so that the effective transition rate γ varies between γmin and γmax where γmin and γmax are chosen so that the epidemic threshold lies between them, i.e., γminγcγmax. Thus when μ~ is small (large), β~ will be small (large) so that their ratio lies in the range γminβ~/μ~γmax. We perform standard synchronous updating simulations (with Δt=1) and obtain the critical value γc as the smallest value of γ such that the fraction of nodes in the metastable state is nonzero. If the epidemic threshold depends only on the ratio γ=β~/μ~ and is independent of the individual values of μ~ and β~, then γc should be the same regardless of the value of μ~, which is fixed. However, we find that this is not the case.

We perform this experiment on an Erdős-Rényi network with N=104 nodes and mean degree k=4, similar to the network used in the example of Sec. IV (Fig. 5) but with a greater number of nodes. On such a network, the epidemic threshold is predicted by steady-state analysis of both the NDLS and HMF theories as γc=0.25. From Fig. 5 we see that when μ~ is small (μ~=0.1), the epidemic threshold predicted by synchronous updating simulations corresponds to this value γc=0.25. However, as μ~ (and thus β~) increases, the accuracy of the discrete-time approach breaks down and both the fraction of infected nodes in the metastable state and the epidemic threshold deviate from the true values. The epidemic threshold decreases from γc=0.25 when μ~=0.1 to γc=0.2 when μ~=1, even though the ratio γ=β~/μ~ remains in the same range. Thus, in discrete-time formalisms the steady-state behavior is not fully determined by the effective transition rate γ but also depends on μ~ and β~. From our analysis in Sec. III (Fig. 3) it is clear that the metastable state reached iteratively from an initial condition depends on the single-step transition probabilities μ~ and β~. If these are too large, the errors introduced in the discrete-time approximation become significant, affecting the metastable state and the value of the epidemic threshold.

The fraction of infected nodes η∞ in the metastable state in an Erdős-Rényi network with 104 nodes and mean degree 〈k〉=4 for various values of the transition probabilities μ~ and β~. Along each of the five curves shown from left to right, μ~ is fixed at the values 1, 0.775, 0.55, 0.325, and 0.1 respectively and β~ varies so that γ=β~/μ~ varies between 0.15 and 0.35. The vertical dashed line indicates the epidemic threshold γc=0.25 predicted by both the NLDS and HMF theories.
FIG. 5.
The fraction of infected nodes η in the metastable state in an Erdős-Rényi network with 104 nodes and mean degree k=4 for various values of the transition probabilities μ~ and β~. Along each of the five curves shown from left to right, μ~ is fixed at the values 1, 0.775, 0.55, 0.325, and 0.1 respectively and β~ varies so that γ=β~/μ~ varies between 0.15 and 0.35. The vertical dashed line indicates the epidemic threshold γc=0.25 predicted by both the NLDS and HMF theories.

The results of this section have important implications for discrete-time approaches. First, they show that the epidemic threshold calculated empirically using synchronous updating simulations can be incorrect if μ~ and β~ are too large, even if the ratio between them is small. Second, they have implications for calculating the epidemic threshold from discrete-time systems of the form pt+1=f(pt) by forward iterating the system from an initial condition [37]. If the transition probabilities μ~ and β~ used in such systems are too large then the metastable state will be affected, possibly leading to a miscalculation of the epidemic threshold.

VI.. CONCLUSIONS

In this paper, we have provided conclusive evidence of the limitations of discrete-time approaches as approximations to continuous-time contagion processes. When the state transition probabilities are too large, such approaches become inaccurate and misrepresentative of the underlying continuous-time processes, thus compromising their utility and their applicability to prediction and analysis.

Our message is clear: Due care needs to be taken when implementing discrete-time methods as approximations to continuous-time dynamical processes. Being constructive, we have briefly discussed alternatives. For simulations of continuous-time processes on networks, event-based simulations such as the Gillespie algorithm are more favorable than synchronous updating schemes both in terms of accuracy and speed. For theoretical analysis, continuous-time analogs [63,68] of discrete-time approaches should be employed because they are unconstrained in the range of dynamics parameter values that can be studied.

ACKNOWLEDGMENTS

This work has been supported by Science Foundation Ireland, Grants No. 11/PI/1026 and No. 12/PI/1683 (P.G.F, S.M, J.P.G.), the James S. McDonnell Foundation postdoctoral fellowship (P.G.F.), and the European Commission Seventh Framework Programme for the FET-Proactive project PLEXMATH (FP7-ICT-2011-8; Grant No. 317614) (P.G.F., S.M., J.P.G.).

REFERENCES

[1] 

M. E. Newman, Networks: An Introduction (Oxford University Press, Oxford, 2010).

[2] 

R. M. Anderson and R. M. May, Nature (London)280, 361 (1979)., doi: 10.1038/280361a0

[3] 

R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys.87, 925 (2015).RMPHAT0034-6861

[4] 

R. Axelrod, J. Conflict Resolution41, 203 (1997)., doi: 10.1177/0022002797041002001

[5] 

D. J. Watts, Proc. Natl. Acad. Sci. USA99, 5766 (2002).PNASA60027-8424

[6] 

C. Castellano, S. Fortunato, and V. Loreto, Rev. Mod. Phys.81, 591 (2009).RMPHAT0034-6861

[7] 

D. Centola, Science329, 1194 (2010).SCIEAS0036-8075

[8] 

D. J. O'Sullivan, G. J. O'Keeffe, P. G. Fennell, and J. P. Gleeson, Front. Phys.3, 71 (2015).2296-424X, doi: 10.3389/fphy.2015.00071

[9] 

J. P. Gleeson, T. Hurd, S. Melnik, and A. Hackett, in Advances in Network Analysis and its Applications (Springer, Berlin, 2012), pp. 27–56.

[10] 

P. Gai and S. Kapadia, Proc. Roy. Soc. A466, 2401 (2010)., doi: 10.1098/rspa.2009.0410

[11] 

R. M. May and N. Arinaminpathy, J. Roy. Soc. Int.7, 823 (2010).1742-5689, doi: 10.1098/rsif.2009.0359

[12] 

J. P. Gleeson, D. Cellai, J.-P. Onnela, M. A. Porter, and F. Reed-Tsochas, Proc. Natl. Acad. Sci. USA111, 10411 (2014).PNASA60027-8424

[13] 

S. Melnik, J. A. Ward, J. P. Gleeson, and M. A. Porter, Chaos23, 013124 (2013).CHAOEH1054-1500

[14] 

J. P. Gleeson, J. A. Ward, K. P. O'Sullivan, and W. T. Lee, Phys. Rev. Lett.112, 048701 (2014).PRLTAO0031-9007

[15] 

J. P. Gleeson, Phys. Rev. Lett.107, 068701 (2011).PRLTAO0031-9007

[16] 

J. P. Gleeson, Phys. Rev. X3, 021004 (2013).2160-3308, doi: 10.1103/PhysRevX.3.021004

[17] 

J. P. Gleeson, S. Melnik, J. A. Ward, M. A. Porter, and P. J. Mucha, Phys. Rev. E85, 026106 (2012).PLEEE81539-3755

[18] 

S. Melnik, A. Hackett, M. A. Porter, P. J. Mucha, and J. P. Gleeson, Phys. Rev. E83, 036112 (2011).PLEEE81539-3755

[19] 

A. Faqeeh, S. Melnik, and J. P. Gleeson, Phys. Rev. E91, 052807 (2015).PLEEE81539-3755

[20] 

M. F. Gomes, A. Piontti, L. Rossi, D. Chao, I. Longini, M. E. Halloran, and A. Vespignani, PLoS: Currents Outbreaks1 (2014).

[21] 

S. González-Bailón, J. Borge-Holthoefer, A. Rivero, and Y. Moreno, Sci. Rep.1, 197 (2011)., doi: 10.1038/srep00197

[22] 

D. Lazer, A. S. Pentland, L. Adamic, S. Aral, A. L. Barabási, D. Brewer, N. Christakis, N. Contractor, J. Fowler, M. Gutmann, Science323, 721 (2009).SCIEAS0036-8075

[23] 

P. Van Mieghem, Performance Analysis of Complex Networks and Systems (Cambridge University Press, Cambridge, 2014).

[24] 

A. Barrat, M. Barthelemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, Cambridge, 2008).

[25] 

D. R. Cox and H. D. Miller, The Theory of Stochastic Processes (CRC Press, Boca Raton, 1977).

[26] 

W. O. Kermack and A. G. McKendrick, Proc. Roy. Soc. A115, 700 (1927)., doi: 10.1098/rspa.1927.0118

[27] 

J. L. Toole, M. Cha, and M. C. González, PloS One7, e29528 (2012).1932-6203, doi: 10.1371/journal.pone.0029528

[28] 

M. Kitsak, L. K. Gallos, S. Havlin, F. Liljeros, L. Muchnik, H. E. Stanley, and H. A. Makse, Nature Phys.6, 888 (2010).1745-2473, doi: 10.1038/nphys1746

[29] 

D. Daley and D. G. Kendall, IMA J. Appl. Math.1, 42 (1965).IJAMDM0272-4960

[30] 

N. M. Ferguson, D. A. Cummings, C. Fraser, J. C. Cajka, P. C. Cooley, and D. S. Burke, Nature (London)442, 448 (2006).NATUAS0028-0836

[31] 

M. L. C. Degli Atti, S. Merler, C. Rizzo, M. Ajelli, M. Massari, P. Manfredi, C. Furlanello, G. S. Tomba, and M. Iannelli, PLoS One3, e1790 (2008).1932-6203, doi: 10.1371/journal.pone.0001790

[32] 

D. Balcan, B. Gonçalves, H. Hu, J. J. Ramasco, V. Colizza, and A. Vespignani, J. Comput. Sci.1, 132 (2010)., doi: 10.1016/j.jocs.2010.07.002

[33] 

L. J. Allen, Math. Biosci.124, 83 (1994).MABIAR0025-5564

[34] 

Y. Wang, D. Chakrabarti, C. Wang, and C. Faloutsos, in Proceedings of the 22nd International Symposium on Reliable Distributed Systems, 2003 (IEEE, Piscataway, 2003), pp. 25–34.

[35] 

C. C. Zou, D. Towsley, and W. Gong, IEEE Transactions on Dependable and Secure Computing,5, 105 (2007).1545-5971, doi: 10.1109/TDSC.2007.1001

[36] 

D. Chakrabarti, Y. Wang, C. Wang, J. Leskovec, and C. Faloutsos, ACM Transactions on Information and System Security (TISSEC)10, 1 (2008).1094-9224, doi: 10.1145/1284680.1284681

[37] 

S. Gómez, A. Arenas, J. Borge-Holthoefer, S. Meloni, and Y. Moreno, Europhys. Lett.89, 38009 (2010)., doi: 10.1209/0295-5075/89/38009

[38] 

D. Trpevski, W. K. S. Tang, and L. Kocarev, Phys. Rev. E81, 056102 (2010)., doi: 10.1103/PhysRevE.81.056102

[39] 

S. Gómez, J. Gómez-Gardenes, Y. Moreno, and A. Arenas, Phys. Rev. E84, 036105 (2011)., doi: 10.1103/PhysRevE.84.036105

[40] 

E. Estrada, F. Kalala-Mutombo, and A. Valverde-Colmeiro, Phys. Rev. E84, 036110 (2011)., doi: 10.1103/PhysRevE.84.036110

[41] 

D. Smilkov and L. Kocarev, Phys. Rev. E85, 016114 (2012)., doi: 10.1103/PhysRevE.85.016114

[42] 

M. De Domenico, A. Lima, P. Mougel, and M. Musolesi, Sci. Rep.3, 2980 (2013)., doi: 10.1038/srep02980

[43] 

X. Wei, N. C. Valler, B. A. Prakash, I. Neamtiu, M. Faloutsos, and C. Faloutsos, IEEE J. Sel. Areas Commun.31, 1049 (2013)., doi: 10.1109/JSAC.2013.130607

[44] 

C. Granell, S. Gómez, and A. Arenas, Phys. Rev. Lett.111, 128701 (2013)., doi: 10.1103/PhysRevLett.111.128701

[45] 

E. Cozzo, R. A. Banos, S. Meloni, and Y. Moreno, Phys. Rev. E88, 050801 (2013)., doi: 10.1103/PhysRevE.88.050801

[46] 

C.-R. Cai, Z.-X. Wu, and J.-Y. Guan, Phys. Rev. E90, 052803 (2014)., doi: 10.1103/PhysRevE.90.052803

[47] 

E. Valdano, L. Ferreri, C. Poletto, and V. Colizza, Phys. Rev. X5, 021005 (2015)., doi: 10.1103/PhysRevX.5.021005

[48] 

J. Gómez-Gardeñes, L. Lotero, S. Taraskin, and F. Pérez-Reche, Sci. Rep.6, 19767 (2016)., doi: 10.1038/srep19767

[49] 

R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett.86, 3200 (2001)., doi: 10.1103/PhysRevLett.86.3200

[50] 

M. Barthélemy, A. Barrat, R. Pastor-Satorras, and A. Vespignani, Phys. Rev. Lett.92, 178701 (2004)., doi: 10.1103/PhysRevLett.92.178701

[51] 

N. T. Bailey, The Mathematical Theory of Infectious Diseases and its Applications (Charles Griffin & Company, London, 1975).

[52] 

T. M. Liggett, Stochastic Interacting Systems: Contact, Voter and Exclusion processes, Vol. 324 (Springer, Berlin, 1999).

[53] 

C. W. Gardiner, Handbook of Stochastic Methods, Springer Series in Synergetics (Springer, Berlin, 1985).

[54] 

T. M. Liggett, Interacting Particle Systems (Springer, Berlin, 1985).

[55] 

D. T. Gillespie, J. Phys. Chem.81, 2340 (1977)., doi: 10.1021/j100540a008

[56] 

S. C. Ferreira, C. Castellano, and R. Pastor-Satorras, Phys. Rev. E86, 041125 (2012)., doi: 10.1103/PhysRevE.86.041125

[57] 

C. L. Vestergaard and M. Génois, PLoS Comput. Biol.11, e1004579 (2015)., doi: 10.1371/journal.pcbi.1004579

[58] 

A. F. Voter, in Radiation Effects in Solids (Springer, Berlin, 2007), pp. 1–23.

[59] 

R. M. Neal, Ann. Stat.31, 705 (2003)., doi: 10.1214/aos/1056562461

[61] 

We refer to the metastable state as the state when the expected fraction of infected nodes has reached a plateau, which very slowly decreases towards 0.

[62] 

C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett.105, 218701 (2010)., doi: 10.1103/PhysRevLett.105.218701

[63] 

A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F. Mendes, Phys. Rev. Lett.109, 128702 (2012)., doi: 10.1103/PhysRevLett.109.128702

[64] 

P. Van Mieghem, Europhys. Lett.97, 48004 (2012)., doi: 10.1209/0295-5075/97/48004

[65] 

M. Boguñá, C. Castellano, and R. Pastor-Satorras, Phys. Rev. Lett.111, 068701 (2013)., doi: 10.1103/PhysRevLett.111.068701

[66] 

A. S. Mata and S. C. Ferreira, Phys. Rev. E91, 012816 (2015)., doi: 10.1103/PhysRevE.91.012816

[67] 

J. O. Kephart and S. R. White, in Proceedings of the 1991 IEEE Computer Society Symposium on Research in Security and Privacy (IEEE, Piscataway, 1991), pp. 343–359.

[68] 

P. Van Mieghem, J. Omic, and R. Kooij, IEEE/ACM Transactions on Networking17, 1 (2009)., doi: 10.1109/TNET.2008.925623

https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1103/PhysRevE.94.052125&title=Limitations of discrete-time approaches to continuous-time contagion dynamics&author=Peter G. Fennell,Sergey Melnik,James P. Gleeson,&keyword=&subject=Articles,Statistical Physics,