Physical Review. E

American Physical Society

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.

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 $\mathrm{\Delta}t$ (which usually takes the value $\mathrm{\Delta}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 $\mathrm{\Delta}t$.

Although discrete-time approaches correspond to their continuous-time counterpart in the limit $\mathrm{\Delta}t\to 0$, they can significantly differ in the case that $\mathrm{\Delta}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 $\mathrm{\Delta}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].

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 $\lambda $, while events in the discrete-time analog are described by transition probabilities $\stackrel{~}{\lambda}$, where $\stackrel{~}{\lambda}=\lambda \mathrm{\Delta}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 $\stackrel{~}{\lambda}$ for which discrete-time approaches are accurate.

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 ${X}_{t}^{i}$, which is either susceptible (${X}_{t}^{i}=S$) or infected (${X}_{t}^{i}=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 $\beta $ per infected neighbor, while infected nodes recover at a rate $\mu $. “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

$\begin{array}{cc}\hfill \mu & =\underset{\mathrm{\Delta}t\to 0}{lim}\frac{P\left({X}_{t+\mathrm{\Delta}t}^{i}=S|{X}_{t}^{i}=I\right)}{\mathrm{\Delta}t},\hfill \end{array}$

$\begin{array}{cc}\hfill \beta & =\underset{\mathrm{\Delta}t\to 0}{lim}\frac{P\left({X}_{t+\mathrm{\Delta}t}^{i}=I\phantom{\rule{4.pt}{0ex}}\text{via}\phantom{\rule{4.pt}{0ex}}j|{X}_{t}^{i}=S,{X}_{t}^{j}=I\right)}{\mathrm{\Delta}t},\hfill \end{array}$

where {${X}_{t+\mathrm{\Delta}t}^{i}=I\phantom{\rule{4.pt}{0ex}}\text{via}\phantom{\rule{4.pt}{0ex}}j$} 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 $\mathrm{\Delta}t\to 0$ leads to the concept of transition rates. In general, we can define ${r}_{i}$ as the rate at which a node $i$ changes from its current state to the opposite state; this is given by
${r}_{i}=\left\{\begin{array}{cc}\hfill \beta {m}_{i,t}& \hfill \text{if}\phantom{\rule{4.pt}{0ex}}{X}_{t}^{i}=S\\ \hfill \mu & \hfill \text{if}\phantom{\rule{4.pt}{0ex}}{X}_{t}^{i}=I\end{array}\right.,$

where ${m}_{i,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 ${\mathbf{Y}}_{t}={\left\{{X}_{t}^{i}\right\}}_{i=1}^{N}$ the state of the network at time $t$, and by $p(\mathbf{y},t)$ the probability that the network is in state ${\mathbf{Y}}_{t}=\mathbf{y}$, then the master equation is given by

$\frac{d}{dt}p(\mathbf{y},t)=\sum _{{\mathbf{y}}^{\prime}\ne \mathbf{y}}(p({\mathbf{y}}^{\prime},t){r}_{{\mathbf{y}}^{\prime}\to \mathbf{y}}-p(\mathbf{y},t){r}_{\mathbf{y}\to {\mathbf{y}}^{\prime}}),$

with initial conditions $p(\mathbf{y},0)={p}_{0}\left(\mathbf{y}\right)$. Here ${r}_{\mathbf{y}\to {\mathbf{y}}^{\prime}}$ is the instantaneous rate at which the network changes from state $\mathbf{y}$ to ${\mathbf{y}}^{\prime}$ and is fully determined by the network structure and the transition rates $\mu $ and $\beta $. While the master equation is the gold standard—exactly describing the evolution of SIS dynamics—the dimension of its sample space $\mathrm{\Omega}$ is ${2}^{N}$, 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 $\tau $ 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 $\tau $ is an exponentially distributed random variable and the parameter of the distribution is the sum of the individual node transition rates, i.e., ${\sum}_{i=1}^{N}{r}_{i}$.

*Lemma 2.* The probability that the next node in the network to change state will be node $i$ is ${r}_{i}/{\sum}_{j=1}^{N}{r}_{j}$.

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 $\tau $ and node $i$ changes its state, where $\tau $ and $i$ are random numbers drawn according to Lemmas 1 and 2 (Fig. 1). Stochastic simulations give the opportunity to construct $p(\mathbf{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.

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 $\mathrm{\Delta}t$. Instantaneous transition rates are replaced by transition probabilities. In a single time interval, susceptible nodes become infected through their infected neighbors with probability $\stackrel{~}{\beta}=\beta \mathrm{\Delta}t$ per infected neighbor, while infected nodes recover with probability $\stackrel{~}{\mu}=\mu \mathrm{\Delta}t$. Note that $\mathrm{\Delta}t$ is often assumed to take the value $\mathrm{\Delta}t=1$, but even in this case it should be included in the expression for $\stackrel{~}{\beta}$ and $\stackrel{~}{\mu}$ 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 $\mu $ and $\beta $ as rates given in Eqs. (1) and (2). These equations can be rearranged to give transition probabilities in terms of these rates, i.e.,

$\begin{array}{cc}& P\left({X}_{t+\mathrm{\Delta}t}^{i}=S|{X}_{t}^{i}=I\right)=\mu \mathrm{\Delta}t\hfill \end{array}$

$\begin{array}{cc}& P\left({X}_{t+\mathrm{\Delta}t}^{i}=I\phantom{\rule{4.pt}{0ex}}\text{via}\phantom{\rule{4.pt}{0ex}}j|{X}_{t}^{i}=S,{X}_{t}^{j}=I\right)=\beta \mathrm{\Delta}t,\hfill \end{array}$

where in this case $\mathrm{\Delta}t$ is an infinitesimally small length of time. In the case that $\mathrm{\Delta}t$ is not infinitesimally small, Eqs. (5) and (6) become approximations. In a time interval of length $\mathrm{\Delta}t$ in the continuous-time Markov process, the exact probability that an infected node will recover is $1-{e}^{-\mu \mathrm{\Delta}t}$ while the probability that a susceptible node will become infected by a given infected neighbor is $1-{e}^{-\beta \mathrm{\Delta}t}$. The transition probabilities $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$, the right-hand side of Eqs. (5) and (6), are approximations to $1-{e}^{-\mu \mathrm{\Delta}t}$ and $1-{e}^{-\beta \mathrm{\Delta}t}$, respectively, and an important question then arises of the effect these approximations have on the dynamics.Figure 2 shows the actual probability $1-{e}^{-\lambda \mathrm{\Delta}t}$ along with the discrete-time probability $\stackrel{~}{\lambda}=\lambda \mathrm{\Delta}t$, where we use the parameter $\lambda $ to represent either $\mu $ or $\beta $. We also plot the error $\epsilon $, which is defined as the difference between the discrete-time probability and the actual probability. When $\stackrel{~}{\lambda}<0.1,\phantom{\rule{0.16em}{0ex}}\epsilon <0.01$ and so the approximation is fairly accurate in this range. For larger values of the state transition probability $\stackrel{~}{\lambda}$, however, the approximation differs significantly from the true values. At $\stackrel{~}{\lambda}=0.5,\phantom{\rule{0.16em}{0ex}}\epsilon \approx 0.1$ and when $\stackrel{~}{\lambda}=1,\phantom{\rule{0.16em}{0ex}}\epsilon \approx 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 $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ are very small they begin to lose accuracy when $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ are of the order of magnitude of ${10}^{-1}$.

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 $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ 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 $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ (and thus $\mu ,\phantom{\rule{0.16em}{0ex}}\beta $, and $\mathrm{\Delta}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 $\mu $ or $\beta $, the time step $\mathrm{\Delta}t$ should be small while if $\mathrm{\Delta}t=1$, as in the case of the majority of discrete-time approaches, the values of $\mu $ and $\beta $ 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., $\mathrm{\Delta}t=1$. In a single time unit, a susceptible node will become infected by its infected neighbors with probability $\stackrel{~}{\beta}$ per infected neighbor while infected nodes become susceptible with probability $\stackrel{~}{\mu}$. Synchronous updating simulations are statistically exact realizations of the discrete-time dynamics; these dynamics are fully described by the discrete-time master equation

$p(\mathbf{y},t+1)=\sum _{{\mathbf{y}}^{\prime}}p({\mathbf{y}}^{\prime},t){q}_{{\mathbf{y}}^{\prime}\to \mathbf{y}},$

where ${q}_{{\mathbf{y}}^{\prime}\to \mathbf{y}}$ is the probability that the network changes from state ${\mathbf{y}}^{\prime}$ to state $\mathbf{y}$ in a time step of length $\mathrm{\Delta}t=1$ and is fully determined by the network structure and the transition probabilities $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ [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 $\mu $ and $\beta $ and various values of $\mathrm{\Delta}t$. We show that discrete-time dynamics can accurately reproduce continuous-time dynamics for small values of $\mathrm{\Delta}t$, but that they incur a breakdown in accuracy as $\mathrm{\Delta}t$ increases. Further to this, we show in Sec. IV that when the time step is fixed to the value $\mathrm{\Delta}t=1$, as in much of the literature, discrete-time approaches break down in accuracy when the transition rates ($\mu $ and $\beta $) 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 $\stackrel{~}{\beta}$ and $\stackrel{~}{\mu}$ can affect the value of the epidemic threshold, even if the effective transition rate defined as $\gamma =\beta /\mu =\stackrel{~}{\beta}/\stackrel{~}{\mu}$ 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 $\mathrm{\Delta}t$. We do this by carrying out synchronous updating simulations for various values of $\mathrm{\Delta}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

${r}_{i}=\left\{\begin{array}{cc}\hfill \phantom{\rule{0.28em}{0ex}}\beta {Z}_{t}& \hfill \text{if}\phantom{\rule{4.pt}{0ex}}{X}_{t}^{i}=S\\ \hfill \phantom{\rule{0.28em}{0ex}}\mu & \hfill \text{if}\phantom{\rule{4.pt}{0ex}}{X}_{t}^{i}=I\end{array}\right.,$

where ${Z}_{t}$ is the number of infected nodes at time $t$ and $\beta $ and $\mu $ 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 ${2}^{N}$ equations for $p(\mathbf{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
$\begin{array}{ccc}\hfill \frac{d}{dt}p(n,t)& =& -(\mu n+\beta n(N-n\left)\right)p(n,t)\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}\mu (n+1)p(n+1,t)+\beta (n-1)\hfill \\ & & \times \phantom{\rule{0.16em}{0ex}}(N-n+1)p(n-1,t),\hfill \end{array}$

for $0\le n\le N$ with initial conditions $p(n,0)={p}_{0}\left(n\right)$. 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 $\beta =1,\phantom{\rule{0.16em}{0ex}}\mu =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 $\mathrm{\Delta}t$. For the numerical simulations, we performed ${10}^{6}$ 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 $\mathrm{\Delta}t=0.01,0.1$ and 1. From Fig. 2, it is clear that these values of $\mathrm{\Delta}t$ with $\mu =\beta =1$ will give a comprehensive range on which to judge the accuracy of the discrete-time approach, while noting that for $\mathrm{\Delta}t=1$ and for these $\mu $ and $\beta $ parameter values the system is deterministic.

We consider the SIS process at time $t=3$ at which stage the expected fraction of infected nodes ${\eta}_{t}={\sum}_{n=0}^{10}np(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 $\mathrm{\Delta}t$ this accuracy can fully break down when $\mathrm{\Delta}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 $\mathrm{\Delta}t=1$ the results are highly inaccurate with all of the probability concentrated on $n=9$, i.e., $p(n,3)={\delta}_{n,9}$. Even the case $\mathrm{\Delta}t=0.1$, while fairly accurate, shows discrepancies in both the probability distribution $p(n,3)$ and the expected fraction of infected nodes ${\eta}_{t}$ [Fig. 3(a)]. Considering that the error $\epsilon $ between $\stackrel{~}{\mu}$ ($\stackrel{~}{\beta}$) and $1-{e}^{-\mu \mathrm{\Delta}t}$ ($1-{e}^{-\beta \mathrm{\Delta}t}$) is less than 0.005 for $\mu =\beta =1$ and $\mathrm{\Delta}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 ${10}^{6}$ realizations for the Gillespie algorithm and for synchronous updating with various values of $\mathrm{\Delta}t$. For $\mathrm{\Delta}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.

To summarize, the accuracy of discrete-time approximations to continuous-time dynamics depends highly on the size of the discrete-time step $\mathrm{\Delta}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 $\mathrm{\Delta}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.

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(\mathbf{y},t+\mathrm{\Delta}t)=f\left({\left\{p({\mathbf{y}}^{\prime},t)\right\}}_{{\mathbf{y}}^{\prime}\in \mathrm{\Omega}}\right),$

where $p(\mathbf{y},t+\mathrm{\Delta}t)$—the probability that the system is in state $\mathbf{y}$ at time $t+\mathrm{\Delta}t$—is a function of the probabilities $p({\mathbf{y}}^{\prime},t)$ for all possible states ${\mathbf{y}}^{\prime}$ in the sample space $\mathrm{\Omega}$. 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 $\mathrm{\Delta}t\to 0$, the dynamics will differ for noninfinitesimal $\mathrm{\Delta}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 $\mathrm{\Delta}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

${p}_{i,t+1}=f\left({\left\{{p}_{j,t}\right\}}_{j=1}^{N}\right),$

for $0\le i\le N$, where ${p}_{i,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 $\mathrm{\Delta}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.
However, the comparison of discrete-time and continuous-time formulations in this manner is biased. Synchronous updating with a time step $\mathrm{\Delta}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 $\langle k\rangle =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 $\mathrm{\Delta}t=1$, as in Ref. [36], along with synchronous updating simulations with a small time step $\mathrm{\Delta}t=0.01$ and Gillespie algorithm simulations. In Fig. 4(c) of Ref. [36], where $\mu =0.72$ and $\beta =0.2$, it can be seen that the fraction $\overline{\eta}$ 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 $\overline{\eta}\approx 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 $\mathrm{\Delta}t=1$, the accuracy of the discrete-time approach decreases as $\mu $ and $\beta $ increase. In the example above, when $\mu $ is decreased from $\mu =0.72$ to $\mu =0.48$ the discrete-time simulations match relatively closer to the continuous-time simulations [Fig. 4(a)], while when $\mu $ is decreased further to $\mu =0.24$ the discrepancy between the two simulations in negligible. Chakrabarti *et al.* state that their model “outperforms (the KW model) when $\mu $ 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 $\mu $ increases.

We conclude with an observation to motivate the next section. As $\mu $ is increased from $\mu =0.48$ [Fig. 4(a)] to $\mu =0.72$ [Fig. 4(b)], the fraction of infected nodes in the metastable state $\overline{\eta}$ (at, for example, $t=100$) decreases for both the continuous-time simulations and discrete-time simulations. However, $\overline{\eta}$ decreases quicker for the continuous-time simulations and so it would seem that the critical value ${\mu}_{c}$ at which $\overline{\eta}$ 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 $\gamma $ is varied. Recall that the effective transition rate is defined as the ratio of the infection rate to the recovery rate, i.e., $\gamma =\beta /\mu $. Depending on the structure of the network and whether the network is finite or infinite, the critical point, or epidemic threshold, ${\gamma}_{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 $\gamma =\beta /\mu =\stackrel{~}{\beta}/\stackrel{~}{\mu}$, the individual values of the transition probabilities $\stackrel{~}{\beta}$ and $\stackrel{~}{\mu}$ 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 ${p}_{t+1}={p}_{t}$ in Eq. (11)], such as in Ref. [36], is completely valid.

We show how $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ affect the value of the epidemic threshold in the following manner. For a given network, we fix the value of $\stackrel{~}{\mu}$ and vary $\stackrel{~}{\beta}$ so that the effective transition rate $\gamma $ varies between ${\gamma}_{\mathrm{min}}$ and ${\gamma}_{\mathrm{max}}$ where ${\gamma}_{\mathrm{min}}$ and ${\gamma}_{\mathrm{max}}$ are chosen so that the epidemic threshold lies between them, i.e., ${\gamma}_{\mathrm{min}}\ll {\gamma}_{c}\ll {\gamma}_{\mathrm{max}}$. Thus when $\stackrel{~}{\mu}$ is small (large), $\stackrel{~}{\beta}$ will be small (large) so that their ratio lies in the range ${\gamma}_{\mathrm{min}}\ll \stackrel{~}{\beta}/\stackrel{~}{\mu}\ll {\gamma}_{\mathrm{max}}$. We perform standard synchronous updating simulations (with $\mathrm{\Delta}t=1$) and obtain the critical value ${\gamma}_{c}$ as the smallest value of $\gamma $ such that the fraction of nodes in the metastable state is nonzero. If the epidemic threshold depends only on the ratio $\gamma =\stackrel{~}{\beta}/\stackrel{~}{\mu}$ and is independent of the individual values of $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$, then ${\gamma}_{c}$ should be the same regardless of the value of $\stackrel{~}{\mu}$, 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={10}^{4}$ nodes and mean degree $\langle k\rangle =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 ${\gamma}_{c}=0.25$. From Fig. 5 we see that when $\stackrel{~}{\mu}$ is small ($\stackrel{~}{\mu}=0.1$), the epidemic threshold predicted by synchronous updating simulations corresponds to this value ${\gamma}_{c}=0.25$. However, as $\stackrel{~}{\mu}$ (and thus $\stackrel{~}{\beta}$) 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 ${\gamma}_{c}=0.25$ when $\stackrel{~}{\mu}=0.1$ to ${\gamma}_{c}=0.2$ when $\stackrel{~}{\mu}=1$, even though the ratio $\gamma =\stackrel{~}{\beta}/\stackrel{~}{\mu}$ remains in the same range. Thus, in discrete-time formalisms the steady-state behavior is not fully determined by the effective transition rate $\gamma $ but also depends on $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$. 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 $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$. 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 $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ 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 ${p}_{t+1}=f\left({p}_{t}\right)$ by forward iterating the system from an initial condition [37]. If the transition probabilities $\stackrel{~}{\mu}$ and $\stackrel{~}{\beta}$ 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.).

Citing articles via

Tweets

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,

© 2020 Newgen KnowledgeWorks | Privacy & Cookie Policy | Powered by: Nova