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.
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 . 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 , which, although introduced as models of disease spread , 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 (which usually takes the value ), 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 .
Although discrete-time approaches correspond to their continuous-time counterpart in the limit , they can significantly differ in the case that is finite. Allen, in her work , 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.  observe that differences between continuous and discrete-time SIS dynamics are substantial when an arbitrary time step of 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].
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 . 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.
Consider SIS dynamics taking place on a network of nodes. This is a continuous-time Markov process where at any time each node in the network has a corresponding state , which is either susceptible () or infected () [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
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 the state of the network at time , and by the probability that the network is in state , then the master equation is given by
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., .
Lemma 2. The probability that the next node in the network to change state will be node is .
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 changes its state, where and are random numbers drawn according to Lemmas 1 and 2 (Fig. 1). Stochastic simulations give the opportunity to construct 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.
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 . Instantaneous transition rates are replaced by transition probabilities. In a single time interval, susceptible nodes become infected through their infected neighbors with probability per infected neighbor, while infected nodes recover with probability . Note that is often assumed to take the value , 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.,
Figure 2 shows the actual probability along with the discrete-time probability , 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 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 and when . 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 .
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 ) 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 should be small while if , 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., . 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
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 . We show that discrete-time dynamics can accurately reproduce continuous-time dynamics for small values of , but that they incur a breakdown in accuracy as increases. Further to this, we show in Sec. IV that when the time step is fixed to the value , 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.
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 . We do this by carrying out synchronous updating simulations for various values of and comparing them against exact results obtained from the master equation. Numerical simulations are carried out in c ++ and the code is available online .
As our example, we consider SIS dynamics on a complete graph of 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
We present the results for SIS dynamics with running on a complete graph with 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 . For the numerical simulations, we performed realizations and obtained the corresponding by taking the fraction of realizations in which there are infected nodes at time . For the synchronous updating simulations, we consider time steps of and 1. From Fig. 2, it is clear that these values of with will give a comprehensive range on which to judge the accuracy of the discrete-time approach, while noting that for and for these and parameter values the system is deterministic.
We consider the SIS process at time at which stage the expected fraction of infected nodes has reached a metastable state (Fig. 3) . At we empirically construct from the synchronous updating simulations and compare it to 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 this accuracy can fully break down when 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 the results are highly inaccurate with all of the probability concentrated on , i.e., . Even the case , while fairly accurate, shows discrepancies in both the probability distribution and the expected fraction of infected nodes [Fig. 3(a)]. Considering that the error between () and () is less than 0.005 for and (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 constructed empirically from the Gillespie algorithm to calculated from the master equation. The Gillespie algorithm is extremely accurate and matches the exact 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 realizations for the Gillespie algorithm and for synchronous updating with various values of . For —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.
To summarize, the accuracy of discrete-time approximations to continuous-time dynamics depends highly on the size of the discrete-time step 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 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.
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
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. , 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
However, the comparison of discrete-time and continuous-time formulations in this manner is biased. Synchronous updating with a time step 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 . The example is SIS dynamics on an Erdős-Rényi network of 1000 nodes and mean degree . Figure 4 shows various numerical simulations of these dynamics. Again, the computer code used to perform the simulations is available from Ref. . Included in Fig. 4 are synchronous updating simulations with a time step , as in Ref. , along with synchronous updating simulations with a small time step and Gillespie algorithm simulations. In Fig. 4(c) of Ref. , where and , 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 . The KW theory, which in Ref.  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.  should be reversed: the KW model is more accurate than the NDLS model.
For fixed , the accuracy of the discrete-time approach decreases as and increase. In the example above, when is decreased from to the discrete-time simulations match relatively closer to the continuous-time simulations [Fig. 4(a)], while when is decreased further to 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 [Fig. 4(a)] to [Fig. 4(b)], the fraction of infected nodes in the metastable state (at, for example, ) 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 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.
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, 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. ). Note however that the epidemic threshold predicted by steady-state analysis [i.e., setting in Eq. (11)], such as in Ref. , 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 and where and are chosen so that the epidemic threshold lies between them, i.e., . Thus when is small (large), will be small (large) so that their ratio lies in the range . We perform standard synchronous updating simulations (with ) and obtain the critical value 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 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 nodes and mean degree , 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 . From Fig. 5 we see that when is small (), the epidemic threshold predicted by synchronous updating simulations corresponds to this value . 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 when to when , 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 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 by forward iterating the system from an initial condition . 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.
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.
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.).