key: cord-1055649-16nd842j
authors: Van Mieghem, P.; Liu, Qiang
title: Explicit non-Markovian susceptible-infected-susceptible mean-field epidemic threshold for Weibull and Gamma infections but Poisson curings
date: 2019-08-26
journal: Physical review. E
DOI: 10.1103/physreve.100.022317
sha: fa2e2b8922dc8b6343776efca4c1baf2b86cb7dc
doc_id: 1055649
cord_uid: 16nd842j

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.

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 G is unweighted and undirected, containing a set N of N nodes (also called vertices) and a set L of L links (or edges). The topology of the graph G is represented by a symmetric N × N adjacency matrix A with spectral radius λ 1 , which is the largest eigenvalue of the adjacency matrix A of the contact graph. In an SIS epidemic process [2] [3] [4] [5] [6] [7] on the graph G, the viral state of a node i at time t is specified by a Bernoulli random variable X i (t ) ∈ {0, 1}: X i (t ) = 0, when node i is healthy, but susceptible, and X i (t ) = 1, when node i is infected at time t. A node i at time t can only be in one of these two states: infected, with probability Pr[X i (t ) = 1], or healthy but susceptible to the infection, with probability Pr[X i (t ) = 0] = 1 − Pr[X i (t ) = 1]. The recovery or curing process for node i is a Poisson process with rate δ and the infection rate over the link (i, j) is a Poisson process with rate β. Only when node i 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 G. 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 R of a node i is the duration between the time at which node i is just infected and the subsequent time at which node i is again recovered and healthy. Analogously, the infection time T is the time needed for a just infected node to infect * P.F.A.VanMieghem@tudelft.nl † Q.L.Liu@tudelft.nl one of its direct, healthy neighbors. In a Poisson process [8] , interevent times are exponentially distributed so that both recovery time R and infection time T are exponential random variables with mean E [R] = 1 δ and E [T ] = 1 β , respectively. We do not consider (a) heterogeneous epidemics, where each node i can have its own curing rate δ i and each link (i, j) its own infection rate β i j nor (b) time-dependent rates β(t ) and δ(t ) 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 T 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 1/ ln(1 + λ 1 ). Simulations hinted that 1/ ln(1 + λ 1 ) 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 T 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.

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 T and the recovery time R can possess an arbitrary distribution, whereas, in Markov theory as mentioned above, T and R are exponential random variables with mean 1 β and 1 δ , respectively. In the homogeneous setting, each node has the same distribution of the curing or recovery time R and each link transfers the viral item following a same distribution of the infection time T . 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 v i∞ of infection of node i is shown in Ref. [12] to obey

which is surprisingly close to the Markovian equation of the N-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 E [M] of infection events during a healthy period, which is specified in [12] by a (complex) integral

where

are the probability generating function of the infection time T and recovery or curing time R, 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 E [M] > 1 λ 1 , 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 T is exponential or the infection follows a Poisson process, then, for any distribution of the recovery time R, it is shown in Ref. [12] that E [M] = τ . Hence, the epidemic threshold in the non-Markovian SIS epidemics with an arbitrary recovery time R reduces, under the mean-field approximation, precisely to the mean-field epidemic threshold 1 τ (1) c = 1 λ 1 of the Markovian SIS epidemics. The other variant, where the recovery time R is exponential and the infection time T is more attractive and explored further in this paper. As shown in Ref. [12] , the NIMFA epidemic threshold criterion Eq. (4) reduces then to . We consider the two most used distributions for the infection time T 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.

While the curing or recovery process is still Poissonean with rate δ, the infection process at each node infects direct neighbors in a time T Weibull 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 E [T Weibull ] to 1 β , so that

The major theoretical reason to choose the Weibull distribution is that, for α = 1, 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 Var[T Weibull ], while a large α corresponds to almost deterministic infection times with small variance Var[T Weibull ].

The probability generating function (pgf) of a Weibull random variable T Weibull is given in Eq. (A1) in Appendix A,

where w = bz and b (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 T Weibull is the solution for τ in

After inversion of Eq. (6), the NIMFA epidemic threshold

The bounds Eq. (A7) on the inverse function w = ϕ −1 T Weibull (y; α) of y = ϕ T Weibull (w; α) 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 τ (1) c;W (α), in which the respective first terms are precisely equal to the above bounds.

Instead of a Weibull distribution, we also consider a Gamma distribution [8, pp. 45-46] for the infection time T ,

with mean E [T Gamma ] = b ξ , variance Var[T Gamma ] = ξ b 2 , and with corresponding pgf

Similar to the Weibull distribution, the Gamma distribution reduces for ξ = 1 to an exponential distribution. After fixing the average infection E [T ] to 1 β , the value of b = 1 βξ . The NIMFA epidemic threshold τ (1) c; (ξ ) is the solution for τ in the criterion Eq. (5), which translates with b = 1 βξ to

from which the NIMFA epidemic threshold follows as

If ξ = k 1 is an integer, then the Gamma random variable equals the sum of k 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 r, but only a small dose of infection is transmitted. A healthy node needs to receive k continuous doses of infection from an infected neighbor to become infected. The infection time T in this interpretation follows a Gamma distribution with ξ = k and b = 1/r. The overall effective infection rate τ = E [R] E [T ] = r kδ . Thus, there exists a dose threshold k c such that, if k > k c , then the Gamma SIS process is below the epidemic threshold, while if k < k c , then it is above the threshold. Here, the dose threshold k c can be a real number. Equating τ (1) c; (ξ ) = r k c δ in Eq. (11) with ξ = k c , we obtain the following dose threshold:

Equation (12) shows that the dose threshold k c increases logarithmically with the largest eigenvalue λ 1 of the contact graph, that can be interpreted as a "dynamic" average nodal degree [24] . When δ r > 1, then

δ r for sufficiently large δ r > 1: the dose threshold k c increases approximately linearly with the transmission rate r 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 τ (1) c; (ξ ) can be transformed into τ (1) c;W (α) for some transform g

The theorem [8, p. 146 ] "Given that exactly one event of a Poisson process has occurred in the interval [0, t], the time of occurrence of this event is uniformly distributed over [0, t]" 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 y(s, t ) that a renewal event happens between s and s + ds, given that precisely one renewal occurs in the interval [0, t], where s t. We denote the interarrival times in the renewal process by τ 1 , τ 2 , . . . and the time that the kth renewal occurs by W k . The interarrival times are i.

dx . The number of renewal events N (t ) at time t is related to the waiting time of W n by the fundamental equivalence

Instead of a fixed length t, we present the derivation for a random time interval and we replace t by a non-negative

dt . We assume that the length V 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 [0, V ] is

Invoking the renewal fundamental equivalence Pr[N (t ) n] = Pr[W n t],

and the law of total probability [8] , p. 23],

we obtain, taking into account that τ and V are independent,

where the convolution (see Ref. [8, p. 160 ]) is

Thus, we find that

has not yet occurred so that its interarrival time must be larger than V − s. Since renewals are independent (in nonoverlapping intervals), it holds that

With the law of total probability,

Combining all yields the probability y V (s) that a renewal event occurs between s and s + ds, given that there is precisely one renewal in an interval [0, V ] with random length V ,

We verify from Eq. (13) that 

which is, unfortunately, demanding to evaluate numerically. Therefore, we proceed with the simplest case where the time interval is fixed V = t and f V (u) = δ D (t − u) is the Dirac function. In that case, y V (s) in Eq. (13) simplifies to the probability y(s, t ) that a renewal event occurs between s and s + ds, given that there is precisely one renewal in an interval

For example, if the interarrival time τ is exponentially dis-

which, indeed, reflects that the occurrence of a Poisson event is uniformly distributed over [0, t] and that y(s, t ) is independent of s, meaning that each time s ∈ [0, t] is equally 

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 σ + dσ , given that there is precisely one renewal in an interval [0, ρ] versus normalized time σ ∈ [0, ρ] for ρ = 1 time unit, expressed in b seconds, and various values of α. Figure 1 indicates that the regime α < 1 models a different behavior than the regime α > 1 (as also follows from the functional equation (A5) of the probability generating function). For α < 1, a Weibull event occurs increasingly likely at smaller times, the smaller α is. For α > 1, 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 t, 2t, 3t, . . ., 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 V , a similar interpretation holds, but the recovery times for each node occur at different times V 1 , V 2 , . . . (almost surely) and these differences per node complicate network-wide prognoses of the evolutions of the SIS epidemic.

c;W (α) EQUALS lim ξ→∞ τ (1) c; (ξ)

Apart from the analysis in Ref. [13] that established

and claimed that any mean-field SIS epidemic threshold τ (1) c 1 ln (1+λ 1 ) , 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. (1) c;W (α) for large α,

The asymptotic expression Eq. (17) illustrates that τ (1) c (α) < 1 ln (1+λ 1 ) , which supports the claim in Ref. [13] . The corresponding limiting epidemic threshold lim ξ →∞ τ (1) c; (ξ ) for Gamma infection times follows from Eq. (11), after substitution x = 1 ξ , as

where de l'Hospital's rule has been used. Since τ (1) c; (ξ ) is increasing in ξ , the claim τ (1) c; (ξ ) < 1 ln (1+λ 1 ) 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 B k , from which an expansion for the Gamma epidemic threshold τ (1) c; (ξ ), analogous to the asymptotic series Eq. (17) for the Weibull epidemic threshold τ (1) c;W (α), follows as

which converges provided that ln (1+λ 1 ) ξ < 2π .

For each α, the epidemic threshold τ (1) c;W (α) in Eq. (7) is expressed in terms of the inverse function ϕ −1 T Weibull (y; α). In Appendix A 5, we derive the Lagrange series for ϕ −1 T Weibull (y; α). In particular, for α 1, the epidemic threshold τ (1) c;W (α) follows after substitution of the Lagrange series Eq. (A20) into Eq. (7) as 1 τ (1) c;W (α)

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 c n (α), defined by Eq. (A21), in Appendix A 5 and mention the interesting result that c n (0) = 1 n c n ( 1 α ) 1 = c n (1) for α 1. As shown in Appendix A 5, both extremes have well-known series. The geometric series 1 1−z = ∞ n=0 z n indicates that c;W (α) 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 d-dimensional lattices (and the worst for a path, a one-dimensional lattice), which explains the larger deviation for a rectangular grid. are positive. Consequently, τ (1) c,W (α; K ) > τ (1) c;W (α), for α > 1. Moreover, c n (α) decreases with α ∈ [1, ∞), which implies that τ (1) c;W (α) increases with α towards τ (1) c;W (∞) = 1 ln (1+λ 1 ) . For large realistic graphs, the spectral radius λ 1 can be large so that ( 1 α +1) → 1 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 z = 1. In fact, for λ 1 → ∞, we find that τ (1) c;W (α) → 0 and the Lagrange series Eq. (19) , valid for α 1, 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 K = 30 terms in the truncated series of Eq. (19) lead to a two-digit accuracy for τ (1) c,W (α; 30) for α ∈ [1, 100], 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) .

When α 1, we find from Eq. (A22) that

Since 1 1+λ 1 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 τ (1) c;W (α) reduces to the estimate in Ref. [12] , Eq. (11)] for a large spectral radius λ 1 ,

More accurate expressions for α < 1 than the above estimate are immediate from Eq. 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 α > 10.

Furthermore, the interpretation of Fig. 1 above intuitively explains the behavior of τ (1) c;W (α) 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 τ (1) c;W (α), is sufficient to cause endemic behavior. As a consequence, τ (1) c;W (α) 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 τ (1) c;W (α). The one-to-one relation between α and ξ of the infection time distributions with the corresponding thresholds τ (1) c;W (α) and τ (1) c; (ξ ), immediately couples the type of viral item k (via a specific value of α k or ξ k or of another distribution) to its endemic impact τ (1) c;W (α k ) and τ (1) c; (ξ k ) 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. (1) c; (ξ ) 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.

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 τ (1) c (α) 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 τ (1) c;W (α) and τ (1) c; (ξ ) increases with α and ξ , respectively, from 0 for α, ξ = 0 to [ln (1 + λ 1 )] −1 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 E [T lognormal ] = exp (μ) exp ( σ 2 2 ) and pgf

The non-Markovian SIS epidemic threshold criterion Eq. (5) with lognormal infection time T and exponential recovery time R, normalized by E [T lognormal ] = 1 β , which allows us to eliminate μ = − ln β − σ 2 2 , is the solution τ (1) c;l (σ ) for τ in

Since (1) c; (ξ ), which further supports the claim that 1 ln (λ 1 +1) 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 τ (1) c;l (σ −1 ) as for τ (1) c;W (α) and τ (1) c; (ξ ). 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 τ (1) c;W (α), τ (1) c; (ξ ) and τ (1) c;l (σ −1 ) may be more suitable than the classical reproduction number R 0 , which is critically reconsidered in Ref. [29] .

We are very grateful to Eric Cator and Caterina Scoglio for their constructive input.

The probability generating function (pgf) ϕ T (z) = E [e −zT ] of a continuous random variable T is defined as the double-sided Laplace transform [8, p. 20] of the probability density function f T (x). Clearly, ϕ T (0) = 1. In this section, we explore properties of the pgf of the Weibull distribution with pdf, defined in Eq. (1), but we replace T Weibull by T to shorten the notation,

and pgf

Let x = u b , then we obtain with w = bz and explicitly expressing the dependence on the "shape" parameter α,

which illustrates that the pgf of a Weibull distribution consists of two parameters, w = zb (which is a complex number) and the real nonnegative number α. In fact, and to be precise, the substitution w = bz means that

which indicates that ϕ T (w; α) should be written as ϕ T b (w; α). Assuming that all moments of T exists, by expanding the entire function e z into a Taylor series that converges everywhere in the complex z-plane, we obtain

which relates all moments E [T k ] to the pgf ϕ T b (w). In the sequel, we simplify the notation, use T as a dimensionless random variable instead of T b and remember that the infection time T is measured in units of b seconds.

Further, differentiating Eq. (A1),

demonstrates, since the integrand is always nonnegative for real w, that ϕ T (w; α) monotonously decreases with w along the real w axis from ϕ T (0; α) = 1 towards lim w→∞ ϕ T (w; α) = 0. After substituting u = x α , we find for all α > 0 that

Differentiating Eq. (A3),

indicates that ϕ T (w; α) is decreasing for real, positive w and all real α > 0 and that | dϕ T (w;α)

Partial integration of the integral in Eq. (A1) yields

Let u = wx, then ϕ T (w; α) = 1 − ∞ 0 e −w −α u α e −u du and comparison with Eq. (A3) leads to a functional equation for the pgf of a Weibull random variable T ,

which illustrates the central role of α = 1.

Only for a few values of α, the pgf ϕ T (w; α) can be analytically evaluated. For α = 1, the Weibull distribution reduces to an exponential and Eq. (A1) is ϕ T (w; 1) = 1 1+w . For α = 2, we find that ϕ T (w; 2) = 1 − we w 2 4 ∞ w 2 e −u 2 du and the functional equation (A5) gives us ϕ T (w; 1 2 

Invoking the duplication formula [25, 6.1.18] of the Gamma function, (2z) = 2 2z 2 √ π (z) (z + 1 2 ) leads to

where the factor

is monotonously decreasing from ∞ at α = 0 towards 0 when α → ∞. Further, r(1) = 1 and r(α) > 1 for α < 1, while r(α) < 1 for α > 1. Thus, the variance Var[T ] is always larger than the mean E [T ] for α < 1, while smaller than the mean in the regime α > 1, which demonstrates, besides the functional equation (A5), the fundamental difference between the two α-regimes. The difficulty to express the shape parameter α as a simple function of the moments of T or probabilistic measures such as the mean and variance, underlines the importance of α as the basic characteristic of the Weibull distribution.

we find from Eq. (A4) the lower bound

The inversion y = ϕ T (w; α) indicates that

An upper bound is immediate from

Combining the lower and upper bound for the inverse function ϕ −1 T (y; α) yields, for any α > 0, an inequality

Invoking the mean-value theorem [30, p. 322] , the bounds on ϕ T (w; α) can be turned into equalities,

for some nonnegative numbers q 1 and q 2 . Inversion yields

and if we are satisfied to use the inverse 2 function W λ (z) of ze −λz (see a series expansion in Ref. [31, pp. 36-37] ), then

The constants q 1 and q 2 , which depend on both y and α, may be fitted from simulations or theory. For example, when α = 1, then ϕ −1 T (y; 1) = 1 y − 1 = e −q 1 y , from which q 1 = − ln (1 − y).

Taking the Mellin transform of ϕ T (w; α),

valid for 0 < Re (s) < α. The inverse Mellin transform [32] gives

(A9)

We present two series for the integral in Eq. (A1) by expanding first e −wx and next e −x α in a Taylor series around x = 0. Thus, first we have

which converges 4 for all w, only when α 1. The limit lim α→∞ ϕ T (w; α) = e −w . The second series is derived analogously, 

showing that the dominant factor k ( 1 α −1)k only tends to zero with k if α 1, for any w.

which converges for α 1. After letting k → k − 1 and using the functional equation of the Gamma function (z + 1) = z (z), we arrive at

from which lim α→0 ϕ T (w; α) = 1 − e −1 . Comparing both series Eqs. (A10) and (A11),

again verifies the functional equation (A5). By closing the contour 5 in Eq. (A9) over the half-plane where lim s→∞ (1 − s α ) (s) = 0, we again find the series Eqs. (A10) and (A11).

Next, we expand ( k α + 1) in a Taylor series around z = 1,

which converges for k a < 1. Introduced into the series Eq. (A10) yields

We assume that the summations can be reversed and obtain

Invoking characteristic coefficients [31] , we can show that 5 By moving the line of integration to the left at −k < c < −k + 1, we encounter k poles of (s) at negative integers and Cauchy's residue theorem [33] shows that

where S (k) n are the Stirling numbers of the second kind [25, 21.1.4] . Hence, we arrive at

which converges if the polynomial m j=0 S ( j) m w j α m for large m and which is useful as an asymptotic expansion 6 for large α.

We deduce an approximate inversion of ϕ T (w; α) = y for large α. Limiting Eq. (A14) to the first few values of m,

where the Euler constant equals γ = 0.577216 . . . and d 2 (z) dz 2 | z=1 = 1.97811, we consider ln ϕ T (w) = ln y and take the logarithm of the last order estimate,

Since α is large, we may expand ln (1

Invoking the Taylor expansion Eq. (A12) for (1 + 1 α ), we obtain 6 Indeed, reversing the summation yields

With the bounds on the Stirling numbers of the second kind, 

After expanding the square root and choosing the minus sign, as the only possible physically realistic case, the above simplifies, for large α, to

− γ 2 ln 2 y

The logarithm of Weierstrass's infinite product for the Gamma function 7

where ζ (z) is the Riemann-Zeta function [35] and γ = 0.557216 . . . is the Euler constant. By differentiation of the Taylor expansion Eq. (A16), we find d 2 (z) dz 2 | z=1 − γ 2 = ζ (2) = π 2 6 = 1.64493 and arrive at the solution of ϕ T (w) = y for large α, w = ϕ −1 T (y; α) = − ln y 1+ 1 α + 1 12 π 2 ln 2 y α 2 3 

(A17) The general inversion of ϕ T (w; α) = y for all α is presented in Eq. (A 5).

The equation ϕ T (w; α) = y is formally solved as w = ϕ −1 T (y; α), where f −1 (.) is the inverse function of f (.) obeying

Since we possess the Taylor series Eq. (A10) in w and Eq. (A11) in 1 w α of ϕ T (w; α), we can directly apply our characteristic coefficients to compute the Lagrange series for the inverse function ϕ −1 T (y; α). Since |ϕ T (w; α)| 1 for Re (w) 0, a solution for positive Re (w) is only possible provided |y| 1.

The Lagrange series around z 0 for the inverse function f −1 (z) of a general function f (z) is given [31] in terms of the Taylor coefficients f k (z 0 ) of the Taylor series of f (z) = 1 τ (1) c; (ξ )

Invoking the binomial series (1 + x) q = ∞ n=0 ( q n )x n , convergent for all q provided |x| < 1, shows that should holds for all n > 2, which is, unfortunately, not the case. However, the transform Eq. (B2) with w = 1.83 seems quite accurate in the sense that τ (1) c; (ξ ) ≈ τ (1) c;W (α) (in the range α ∈ [1, 10] ) and the difference between right-and left-hand side of the higher-order coefficient is positive and smaller than 0.035 for 2 < n < 10 and α = 1, 2, . . . 10. Figure 6 illustrates the correspondence for α max = 100 with w = 1.64 in Eq. (B2) on the same graphs considered in Figs. 2-4. The goodness of an approximate ξ -axis scaling is not so surprising, because (a) τ (1) c; (ξ ) = τ (1) c;W (α) for precisely three points: both end-points 0 and ∞ and α = ξ = 1, and (b) we impose that the two distributions have the same mean and variance.

Universality of the SIS prevalence in networks

The Mathematical Theory of Infectious Diseases and its Applications

Infectious Diseases of Humans: Dynamics and Control

Epidemic Modeling: An Introduction

Mathematical Tools for Understanding Infectious Disease Dynamics

Epidemic processes in complex networks

Mathematics of Network Epidemics: From Exact to Approximate Models

Performance Analysis of Complex Networks and Systems

SIS epidemics with time-dependent rates describing ageing of information spread and mutation of pathogens

Mathematics of Epidemics on Networks: From Exact to Approximate Models

Non-Markovian Infection Spread Dramatically Alters the Susceptible-Infected-Susceptible Epidemic Threshold in Networks

Susceptible-infected-susceptible epidemics on networks with general infection and curing times

Burst of virus infection and a possibly largest epidemic threshold of non-Markovian susceptibleinfected-susceptible processes on networks

Simulating non-Markovian stochastic processes

A note on generation times in epidemic models

Enhanced hygiene measures and norovirus transmission during an outbreak

Clinical spectrum and transmission characteristics of infection with norwalk-like virus: Findings from a large community outbreak in sweden

Estimation of the serial interval of influenza

Transmission dynamics and control of severe acute respiratory syndrome

Transmission dynamics and control of Ebola virus disease (EVD): A review

The N-Intertwined SIS epidemic network model

Virus spread in networks

Graph Spectra for Complex Networks

Handbook of Mathematical Functions

Divergent Series

Accuracy criterion for the mean-field approximation in SIS epidemics on networks

Lognormal infection times of online information spread

Measurability of the epidemic reproduction number in data-driven contact networks

A Course of Pure Mathematics

The asymptotic behavior of queueing systems: Large deviations theory and dominant pole approximation, Queueing Syst

Introduction to the Theory of Fourier Integrals

The Theory of Functions

Die Gammafunktion: Band I. Handbuch der Theorie der Gammafunktion und Band II

The Theory of the Zeta-function

∞ k=0 f k (z 0 )(z − z 0 ) k around z 0 aswhere s * [k, n](z 0 ) = s[k, n](z 0 )| ∀m : f m (z 0 )→ f m+1 (z 0 ) and the characteristic coefficient s[k, m] is defined asand obeys the recursionApplying Eq. (A18) to the Taylor series Eq. (A10) forwhere the coefficients are defined as c 1 (α) = 1 (α+1) and, for n > 1,The Langrange series Eq. (A20) converges fast if y is close to 1, but slowly for small y. In fact, the pgf ϕ T (z) = E [e −zT ] is ϕ T (0; α) = 1 and ϕ T (∞; α) = 0, thus ϕ −1 T (0; α) = ∞. All coefficients c n ( 1 α ) of the series in Eq. (A20) are positive and smaller than 1. If α = 1, then ϕ −1 T (y; 1) = − ln (y). Since 1When we apply Eq. (A18) to the Taylor series Eq. (A11) for α 1 around z 0 = 0, written with v = 1 w α asThe solution of g(v) = y equals v = g −1 (y) and, since w = v − 1 a , we arrive at w = [g −1 (y)] − 1 α and the Langrange seriesWe list the first coefficients c n (α) in the Lagrange series Eq. (A22), up to n = 7,,, − 7{3 (5a + 1) (2a + 1) + 5 (3a + 1) (4a + 1)} 720 2 (a + 1) + (6a + 1) 720 (a + 1) , c 7 (α) = 33 6 (2a + 1) 16 6 (a + 1) − 55 (3a + 1) 4 (2a + 1) 16 5 (a + 1) + 5{ (4a + 1) 3 (2a + 1) + 2 2 (3a + 1) 2 (2a + 1)} 8 4 (a + 1) − {27 (5a + 1) 2 (2a + 1) + 90 (3a + 1) (4a + 1) (2a + 1) + 20 3 (3a + 1)} 360 3 (a + 1) + {4 (6a + 1) (2a + 1) + 5 2 (4a + 1) + 8 (3a + 1) (5a + 1)} 720 2 (a + 1) − (7a + 1) 5040 (a + 1) .Both Langrange series Eqs. (A20) and (A22) can be transformed into each other after α → 1 α and (1 − y) → y and indicate that the coefficients c n (α) are only needed for α ∈ [0, 1]. Alternatively, the functional equation (A5) indicates that y = ϕ T (w; α) or w = ϕ −1 T (y; α) as well asHence, for all α 0, we find the functional equation for the inverse pgf of a Weibull random variable T ,The upper and lower bound in Eq. (A7) are the first terms in the Langrange series Eqs. (A20) and (A22), respectively, and can be deduced from those Lagrange series [because c n (α) is positive]. We know already that ϕ −1 T (y; 1) = 1 y − 1 and that ϕ −1 T (y; ∞) = − ln (y). Figure First, let x = y q in Eq. (9), thenand comparison with the Weibull pdf Eq. (1) shows thatwhich illustrates that both pdfs can be mapped into one another by a "power law" transform x = y α . Next, we rewrite the Gamma epidemic threshold in Eq. (11) as