Explicit non-Markovian susceptible-infected-susceptible mean-field epidemic threshold for Weibull and Gamma infections but Poisson curings.

Although non-Markovian processes are considerably more complicated to analyze, real-world epidemics are likely non-Markovian, because the infection time is not always exponentially distributed. Here, we present analytic expressions of the epidemic threshold in a Weibull and a Gamma SIS epidemic on any network, where the infection time is Weibull, respectively, Gamma, but the recovery time is exponential. The theory is compared with precise simulations. The mean-field non-Markovian epidemic thresholds, both for a Weibull and Gamma infection time, are physically similar and interpreted via the occurrence time of an infection during a healthy period of each node in the graph. Our theory couples the type of a viral item, specified by a shape parameter of the Weibull or Gamma distribution, to its corresponding network-wide endemic spreading power, which is specified by the mean-field non-Markovian epidemic threshold in any network.

INTRODUCTION

We confine ourselves here to a particularly simple epidemic model, the susceptible-infected-susceptible (SIS) epidemic process on a graph that, as argued earlier [1], “allows for the highest degree of analytic treatment, which is a major motivation for the continued effort toward its satisfactory understanding.” The graph is unweighted and undirected, containing a set of nodes (also called vertices) and a set of links (or edges). The topology of the graph is represented by a symmetric adjacency matrix with spectral radius , which is the largest eigenvalue of the adjacency matrix of the contact graph. In an SIS epidemic process [2-7] on the graph , the viral state of a node at time is specified by a Bernoulli random variable : , when node is healthy, but susceptible, and , when node is infected at time . A node at time can only be in one of these two states: infected, with probability , or healthy but susceptible to the infection, with probability . The recovery or curing process for node is a Poisson process with rate and the infection rate over the link is a Poisson process with rate . Only when node is infected can it infect its healthy direct neighbors with rate . All Poisson curing and infection processes are independent. This description defines the continuous-time, Markovian homogeneous SIS epidemic process on a graph . The SIS process can be alternatively described in terms of “time” random variables, which is valuable for the remainder of this article. The recovery time of a node is the duration between the time at which node is just infected and the subsequent time at which node is again recovered and healthy. Analogously, the infection time is the time needed for a just infected node to infect one of its direct, healthy neighbors. In a Poisson process [8], interevent times are exponentially distributed so that both recovery time and infection time are exponential random variables with mean and , respectively. We do not consider (a) heterogeneous epidemics, where each node can have its own curing rate and each link its own infection rate nor (b) time-dependent rates and as in Ref. [9]. Perhaps, more importantly, we limit ourselves to a mean-field approximation of the SIS process.

Since real epidemics may not be Markovian, non-Markovian epidemic modeling has been studied for a long time (see, e.g., Refs. [10] and [6], p. 951). For SIS epidemic processes on networks, a non-Markovian infection was reported in Ref. [11] to alter the epidemic threshold significantly. The SIS process with general infection and recovery times has been analyzed in Ref. [12]. In the Weibullian SIS model, as coined in Ref. [13], the exponential distribution of the infection time is extended to a Weibull distribution with probability density function (pdf) [8, p. 56], which is explained in Sec. III and Appendix A. In Refs. [11,12], the Weibullian SIS model is extensively studied, generalizes the Markovian process, and parameterizes the non-Markovian behavior via the shape parameter in Eq. (1). The limiting case of the Weibullian SIS process is analyzed by a mean-field governing equation in Ref. [13], when the Weibull probability density function Eq. (1) reduces to a Dirac delta function and the corresponding epidemic threshold was found to be . Simulations hinted that is the largest epidemic threshold for any Weibull infection time. That result is partly verified in Ref. [14] by an independent simulation, based on Ref. [15] and revisited here in Sec. V.

In epidemiology, the infection time is called the generation time [16], which characterizes the infectivity of pathogens and is defined as the time between the infections (or the symptoms onsets) of the primary case and the secondary case infected by the primary case. The generation time is usually obtained by monitoring the first cases and secondary cases in households and follows skewed distributions, which have been fitted by the Gamma, Weibull, or lognormal distribution. For example, Heijne et al. [17] evaluated a norovirus outbreak [18] in Sweden in 1999 by fitting the generation time with a Gamma distribution. Cowling et al. [19] fitted the generation time of influenza in Hong Kong with all the three distributions and indicated that the Weibull distribution performs slightly better than the Gamma and lognormal distributions, based on the Akaike information criterion. The generation time of the Severe Acute Respiratory Syndrome (SARS) in Singapore in 2003 is well fitted by the Weibull distribution [20]. Furthermore, the Weibull distribution was also used to fit [21] the generation time of an Ebola outbreak in Uganda in 2000. In summary, real disease measurements suggest to consider Weibull, Gamma, and lognormal infection times.

Here, we present analytic expressions for the mean-field epidemic threshold of the Weibull and Gamma SIS model on any network. First, we summarize in Sec. II the non-Markovian mean-field method developed earlier in Ref. [12]. In Sec. III, we apply the method to Weibull and Gamma infection times and exponential recovery times. Section IV considers the occurrence time of an infection in a given interval and generalizes a fundamental theorem in Poisson theory to a renewal process setting, which helps to interpret the behavior of the non-Markovian mean-field epidemic threshold for any infection time distribution. We present Lagrange series for the Weibull mean-field epidemic threshold in Sec. VI, which illustrates that a Weibull infection time leads to considerably more complicated computation than a Gamma infection time. Section VII compares the non-Markovian mean-field epidemic threshold for both Weibull and Gamma distribution, whereas the conclusion in Sec. VIII concisely covers the lognormal distribution as well.

BRIEF REVIEW OF NON-MARKOVIAN SIS EPIDEMICS ON NETWORKS

We briefly review the main results in Ref. [12] of the non-Markovian SIS analysis on an arbitrary network, where the infection and recovering process are assumed to be independent (as in the Markovian case). However, both the infection time and the recovery time can possess an arbitrary distribution, whereas, in Markov theory as mentioned above, and are exponential random variables with mean and , respectively. In the homogeneous setting, each node has the same distribution of the curing or recovery time and each link transfers the viral item following a same distribution of the infection time . Only the metastable state of the non-Markovian SIS epidemic process is analyzed in Ref. [12] in a mean-field setting. Assuming that the metastable state exists and invoking renewal theory, the mean-field steady state probability of infection of node is shown in Ref. [12] to obey which is surprisingly close to the Markovian equation of the -intertwined mean-field approximation (NIMFA) [22], Clearly, the Markovian mean-field steady-state regime is transformed into a non-Markovian one (and vice versa), if we replace the effective infection rate by the average number of infection events during a healthy period, which is specified in [12] by a (complex) integral where and are the probability generating function of the infection time and recovery or curing time , respectively.

The analogy with the NIMFA equations in Refs. [22,23] immediately leads to a definition of the mean-field epidemic threshold in non-Markovian SIS epidemics, Thus, if , then the epidemic process is eventually endemic (in the mean-field approximation), in which a nonzero fraction of the nodes remains infected, else the epidemic process dies out after which the network is eventually overall healthy.

If the infection time is exponential or the infection follows a Poisson process, then, for any distribution of the recovery time , it is shown in Ref. [12] that . Hence, the epidemic threshold in the non-Markovian SIS epidemics with an arbitrary recovery time reduces, under the mean-field approximation, precisely to the mean-field epidemic threshold of the Markovian SIS epidemics.

The specification of the actual epidemic threshold , even in the thermodynamic limit when the size of the graph tends to infinity and reduces to one value, is still an open problem in Markovian (and certainly in non-Markovian) epidemics.

The other variant, where the recovery time is exponential and the infection time is more attractive and explored further in this paper. As shown in Ref. [12], the NIMFA epidemic threshold criterion Eq. (4) reduces then to

GENERAL INFECTION TIME AND EXPONENTIAL RECOVERY TIME

To compare different distributions for the infection time and recovery time , we require that they all have the same average and , so that the effective infection rate in Markovian SIS epidemics equals .

We consider the two most used distributions for the infection time in real diseases, a Weibull and Gamma distribution, whose main difference lies in the tail behavior. Physically, both distributions provide the same insight, although they are computationally very different. We will show that the Gamma distribution possesses a simple analytic mean-field epidemic threshold Eq. (11) below, while the corresponding computation for the Weibull (Sec. VI) is involved. However, the computation of the time occurrence of an infection event in a given time interval (see Fig. 1 below) is more elegant for the Weibull random variable, because its distribution function has an closed analytic form (see Appendix A) in contrast to the Gamma random variable.

FIG. 1.

The pdf of the occurrence of a Weibull event in the interval for various values of .

The infection time has a Weibull distribution

While the curing or recovery process is still Poissonean with rate , the infection process at each node infects direct neighbors in a time with a Weibull pdf Eq. (1), mean , and variance Eq. (A6), computed in Appendix A, To compare the Weibull with the exponential distribution, we fix the average infection time to , so that

The major theoretical reason to choose the Weibull distribution is that, for , the Weibull distribution reduces to the exponential distribution and, hence, to Markovian SIS epidemics. A small shape parameter in Eq. (1) corresponds to heavy tails and a large variance , while a large corresponds to almost deterministic infection times with small variance .

The probability generating function (pgf) of a Weibull random variable is given in Eq. (A1) in Appendix A, where and (in seconds) is chosen as the unit of time. The NIMFA epidemic threshold criterion Eq. (5) of the non-Markovian SIS process with Weibullian infection time is the solution for in After inversion of Eq. (6), the NIMFA epidemic threshold is equal to The bounds Eq. (A7) on the inverse function of show that the NIMFA mean-field epidemic threshold is bounded by In Sec. VI, we will derive exact Lagrange series Eqs. (20) and (19) for the Weibull NIMFA epidemic threshold , in which the respective first terms are precisely equal to the above bounds.

The infection time has a Gamma distribution

Instead of a Weibull distribution, we also consider a Gamma distribution [8, pp. 45–46] for the infection time , with mean , variance , and with corresponding pgf Similar to the Weibull distribution, the Gamma distribution reduces for to an exponential distribution. After fixing the average infection to , the value of . The NIMFA epidemic threshold is the solution for in the criterion Eq. (5), which translates with to from which the NIMFA epidemic threshold follows as

If is an integer, then the Gamma random variable equals the sum of independent and identically distributed exponential random variables [8, pp. 45–46]. Thus, the SIS model with a Gamma infection time can be interpreted as a dose-infection process: Each infected node can infect each healthy neighbor via a Poisson process with rate , but only a small dose of infection is transmitted. A healthy node needs to receive continuous doses of infection from an infected neighbor to become infected. The infection time in this interpretation follows a Gamma distribution with and . The overall effective infection rate . Thus, there exists a dose threshold such that, if , then the Gamma SIS process is below the epidemic threshold, while if , then it is above the threshold. Here, the dose threshold can be a real number. Equating in Eq. (11) with , we obtain the following dose threshold: Equation (12) shows that the dose threshold increases logarithmically with the largest eigenvalue of the contact graph, that can be interpreted as a “dynamic” average nodal degree [24]. When , then leading to for sufficiently large : the dose threshold increases approximately linearly with the transmission rate of each dose of infection.

In summary, the Weibull analysis is characterized by , whereas the key parameter for the Gamma infection time is . Appendix B investigates their relation as well as whether can be transformed into for some transform so that leads to .

OCCURRENCE OF INFECTION EVENTS IN AN INTERVAL

The theorem [8, p. 146] “Given that exactly one event of a Poisson process has occurred in the interval , the time of occurrence of this event is uniformly distributed over ” is of a remarkable simplicity. That theorem associates the Poisson process to the uniform distribution and indicates that a Poisson process can be viewed as the “most” random process, void of any correlation. Here, we investigate the generalization of this theorem to a renewal process [8, Chapter 8], where events still occur independent of the others, but where the interarrival time of renewal events has a general distribution, instead of an exponential distribution as in the Poisson process. In particular, we compute the probability that a renewal event happens between and , given that precisely one renewal occurs in the interval , where . We denote the interarrival times in the renewal process by and the time that the renewal occurs by . The interarrival times are i.i.d with distribution and probability density function (pdf) . The number of renewal events at time is related to the waiting time of by the fundamental equivalence .

Instead of a fixed length , we present the derivation for a random time interval and we replace by a non-negative random variable , with distribution and pdf . We assume that the length and the interarrival time of the renewal process are independent (as in the non-Markovian SIS epidemic). Applying the formula for the conditional probability, the occurrence of precisely one renewal event in the interval is Invoking the renewal fundamental equivalence , and the law of total probability [8], p. 23], we obtain, taking into account that and are independent, where the convolution (see Ref. [8, p. 160]) is

Thus, we find that The event is equivalent to , because the second renewal event has not yet occurred so that its interarrival time must be larger than . Since renewals are independent (in nonoverlapping intervals), it holds that With the law of total probability, Combining all yields the probability that a renewal event occurs between and , given that there is precisely one renewal in an interval with random length , We verify from Eq. (13) that , because becomes, after reversal of the - and -integration, equal to the denominator. For exponential distributions, and , we find from Eq. (13) that . For a Weibull distribution Eq. (1) and in Eq. (15) becomes which is, unfortunately, demanding to evaluate numerically.

Therefore, we proceed with the simplest case where the time interval is fixed and is the Dirac function. In that case, in Eq. (13) simplifies to the probability that a renewal event occurs between and , given that there is precisely one renewal in an interval , For example, if the interarrival time is exponentially distributed as , then so that which, indeed, reflects that the occurrence of a Poisson event is uniformly distributed over and that is independent of , meaning that each time is equally likely. The uniform spread of a Poisson event over can be regarded as a reflection of “complete,” unbiased and uncorrelated, randomness of the Poisson process. For a Weibull distribution Eq. (1), we find from Eq. (15) or from Eq. (14) Let and measure time in units of , then with normalized time and , we have For Gamma distribution Eq. (9), Eq. (15) reduces to which is less attractive than the Weibull case. Thus, we confine ourselves further to the Weibullian infection time.

Figure 1 illustrates the scaled probability that a Weibullian renewal event occurs between and , given that there is precisely one renewal in an interval versus normalized time for time unit, expressed in seconds, and various values of . Figure 1 indicates that the regime models a different behavior than the regime (as also follows from the functional equation (A5) of the probability generating function). For , a Weibull event occurs increasingly likely at smaller times, the smaller is. For , the occurrence of a Weibull event is more evenly distributed over the entire interval with preference at later times. Figure 1 may also be consulted in practice to characterize a disease, if the occurrence time of an infection can be measured for an SIS process (with repeated infections and recoveries for each node). Indeed, if the recovery of each node is deterministic at , then Fig. 1 shows that Weibullian infections with small shape parameter occur most likely soon after the recovery moment, implying that the node is infected most of the time. A large , however, reflects that a node is with high probability healthy for a long time, but becomes infected, just before the recovery of the next cycle takes place. For a random time interval , a similar interpretation holds, but the recovery times for each node occur at different times (almost surely) and these differences per node complicate network-wide prognoses of the evolutions of the SIS epidemic.

THE EPIDEMIC THRESHOLD EQUALS

Apart from the analysis in Ref. [13] that established and claimed that any mean-field SIS epidemic threshold , we present here different derivations and an asymptotic result for the Weibull and Gamma infection time, that support our earlier claim.

First, the limit case in Eq. (16) is immediate from Eqs. (5) and (A3),

A second derivation interprets the general Eq. (4) directly, without resorting to the integral representation in Eq. (3). For , the average number of infection events during a healthy period is computed as follows. Without loss of generality, we assume that a node is infected at time . The infected node infects its neighbors at times , until node is recovered at time . The recovery time follows an exponential distribution with expectation . Thus, if the recovery time falls into , then the infected node infects each of its neighbor times, and we obtain Equating in Eq. (4) again leads to Eq. (16).

Direct application of Eq. (A17) for and leads to the asymptotic expression of the epidemic threshold for large , The asymptotic expression Eq. (17) illustrates that , which supports the claim in Ref. [13].

The corresponding limiting epidemic threshold for Gamma infection times follows from Eq. (11), after substitution , as where de l'Hospital's rule has been used. Since is increasing in , the claim is again demonstrated, now based on the Gamma distribution.

Finally, we rewrite Eq. (11) as and recognize the right-hand side as the famous generating function [25, Sec. 23] of the Bernoulli numbers , from which an expansion for the Gamma epidemic threshold , analogous to the asymptotic series Eq. (17) for the Weibull epidemic threshold , follows as which converges provided that .

LAGRANGE SERIES APPROACH FOR

For each , the epidemic threshold in Eq. (7) is expressed in terms of the inverse function . In Appendix A 5, we derive the Lagrange series for .

In particular, for , the epidemic threshold follows after substitution of the Lagrange series Eq. (A20) into Eq. (7) as The companion series of Eq. (19) for the Gamma infection time is given in the Appendix in Eq. (B1). We explicitly listed the first seven coefficients of , defined by Eq. (A21), in Appendix A 5 and mention the interesting result that for . As shown in Appendix A 5, both extremes have well-known series. The geometric series indicates that indeed agreeing with the Markovian NIMFA threshold , while the Taylor series , obtained by integrating the geometric series, shows that which agrees with our earlier result Eq. (16). When only the first terms are computed, the truncated series at the first terms is a lower bound for the infinite series in Eq. (19), because all coefficients are positive. Consequently, , for . Moreover, decreases with , which implies that increases with towards .

For large realistic graphs, the spectral radius can be large so that for sufficiently large and the right-hand side Lagrange series converges slowly with a comparable convergence rate as the geometric series toward its pole at . In fact, for , we find that and the Lagrange series Eq. (19), valid for , diverges! However, in these limit regimes, either the asymptotic expression Eq. (17) is applicable or the series needs to be transformed, e.g., by the Euler transform [26]. Since terms in the truncated series of Eq. (19) lead to a two-digit accuracy for for , which is indistinguishable from the exact computation on a plot as Fig. 2, we content ourselves to ignore further numerical considerations of the Lagrange series Eq. (19).

FIG. 2.

The mean-field epidemic threshold of a Weibullian SIS process versus the shape parameter in different type of graphs. Simulations of the precise Weibullian SIS process are compared with mean-field theory and with a mean-field asymptotic approximation Eq. (17). It is known [27] that the SIS mean-field approximation is reasonably accurate for ER graphs, less accurate for scale-free graphs, but inaccurate for -dimensional lattices (and the worst for a path, a one-dimensional lattice), which explains the larger deviation for a rectangular grid.

When , we find from Eq. (A22) that Since is small for realistic graphs, the series in Eq. (20) converges amazingly fast and only a few terms are sufficient. Indeed, the first term in the series Eq. (20) for the epidemic threshold reduces to the estimate in Ref. [12], Eq. (11)] for a large spectral radius , More accurate expressions for than the above estimate are immediate from Eq. (20) by incorporating a few higher order terms in . Just in the regime, where a few terms () in the series Eq. (20) already provide very precise values, accurate simulations to deduce are very difficult, which underlines the utility of the series Eq. (20) in practice.

NUMERICAL EVALUATION OF AND

Figure 2 shows the epidemic threshold of the Weibullian SIS process for different values of the shape parameter on (a) an Erdős-Réyni graph on nodes with link density , (b) a rectangular grid with nodes and links, and (c) a scale-free, Barabási-Albert graph with and links. The simulation results are obtained by averaging over realizations with curing rate and each realization runs for 50 time units. Since the Weibullian infection time may cause oscillations in the prevalence (as illustrated in Fig. 1 in Ref. [13]), the threshold is chosen as the value of the effective infection rate , which leads to the maximum prevalence around 0.001 at the last oscillating period.

The simulations agree with the Lagrange series Eqs. (19) and (20) and with numerical inversion of Eq. (6). Figure 2, shows, interestingly, that the asymptotic expansion Eq. (17) is already accurate for reasonably small .

Furthermore, the interpretation of Fig. 1 above intuitively explains the behavior of versus the shape parameter in Figure 2. In the small regime, infections occur predominantly early during the healthy period for each node, so that the larger part of the nodes is longer infected. To cause a network-wide persistence of the epidemic, the viral item only needs a small “push” to infect a substantial part of the nodes long enough to enter the endemic state. In other words, a weak strength of the viral item to infect nodes, which is related to the effective infection rate , is sufficient to cause endemic behavior. As a consequence, is small for small , which agrees with Figure 2. The other regime for large is understood analogously. Most nodes are likely most of the time healthy and infected only shortly when is large. Hence, the strength of the viral item to cause the network-wide persistence of the epidemic must be high, resulting in a high effective infection rate .

The one-to-one relation between and of the infection time distributions with the corresponding thresholds and , immediately couples the type of viral item (via a specific value of or or of another distribution) to its endemic impact and on any contact network (assuming a mean-field approximation). Both Fig. 2 for the Weibull and Fig. 3 for the Gamma SIS threshold look very similar.

FIG. 3.

The NIMFA epidemic threshold of a Gamma SIS process versus the parameter for the same graphs as in Fig. 2. The theory in Eq. (11) is also added in full line.

CONCLUSION

A mathematical analysis of the Weibull probability generating function and its inverse function in Appendix A have led to analytical expressions for the mean-field epidemic threshold of the Weibullian SIS process on a graph as a function of the shape parameter . Similar results for a Gamma infection time are deduced. The mean-field epidemic threshold and increases with and , respectively, from 0 for to for . The Lagrange series Eqs. (19) and (20), and perhaps the asymptotic expansion Eq. (17), may be used to rapidly compute the epidemic threshold of a real epidemic, once its shape parameter is specified.

The Weibull and Gamma distribution were focal here, while digital spread of information, for instance on Twitter [28], has an infection time close to a lognormal [8, p. 60-64], with mean and pgf The non-Markovian SIS epidemic threshold criterion Eq. (5) with lognormal infection time and exponential recovery time , normalized by , which allows us to eliminate , is the solution for in Since , where is the Dirac function, the threshold criterion Eq. (5) with lognormal infection time becomes Thus, it holds that , which further supports the claim that is a maximum possible SIS epidemic threshold. The solution for in Eq. (23) poses a similar difficulty as for the Weibull infection time and is placed on the agenda for future research, because simulations in Fig. 4 exhibit a similar behavior for the lognormal epidemic threshold as for and . The criterion Eq. (4) for the epidemic threshold shows that the analysis can be repeated for any other distribution of the infection time and other recovery times than exponential. Thus, the presented Weibull and Gamma analysis may serve as a guideline for computations with other distributions. Since the SIR epidemic threshold is slightly higher than the SIS epidemics (due to re-infections that potentially lead to more infected nodes), the SIS non-Markovian threshold can be regarded as a lower bound for SIR, so enlarging the scope of the presented theory. Finally, as the influence of the underlying contact graph is incorporated, the non-Markovian thresholds and may be more suitable than the classical reproduction number , which is critically reconsidered in Ref. [29].

FIG. 4.

The NIMFA epidemic threshold of a lognormal SIS process versus the parameter for the same graphs as in Fig. 2. The theory (full line) is obtained by numerical evaluation of Eq. (23).