Epidemic threshold has always been a very hot topic for studying epidemic dynamics on complex networks. The previous studies have provided different theoretical predictions of the epidemic threshold for the susceptible-infected-recovered (SIR) model, but the numerical verification of these theoretical predictions is still lacking. Considering that the large fluctuation of the outbreak size occurs near the epidemic threshold, we propose a novel numerical identification method of SIR epidemic threshold by analyzing the peak of the epidemic variability. Extensive experiments on synthetic and real-world networks demonstrate that the variability measure can successfully give the numerical threshold for the SIR model. The heterogeneous mean-field prediction agrees very well with the numerical threshold, except the case that the networks are disassortative, in which the quenched mean-field prediction is relatively close to the numerical threshold. Moreover, the numerical method presented is also suitable for the susceptible-infected-susceptible model. This work helps to verify the theoretical analysis of epidemic threshold and would promote further studies on the phase transition of epidemic dynamics.

Epidemic threshold, which is one of the most important features of the epidemic dynamics, has attracted much attention recently. The existing studies have provided different theoretical predictions for epidemic threshold of the susceptible-infected-recovered (SIR) model on complex networks, while the numerical verification of these theoretical predictions is still lacking. As a result, it is very necessary to develop an effective numerical measure for identifying the SIR epidemic threshold. In this paper, the numerical identification of the SIR epidemic threshold is systematically studied. We present a numerical method by analyzing the peak of the epidemic variability to identify the epidemic threshold. To understand the effectiveness of the variability measure, the distribution of outbreaks sizes is investigated near the epidemic threshold on random regular networks. Based on the analysis of the cutoff hypothesis of the outbreak size distribution, we find that the variability measure can provide an excellent identification of the epidemic threshold. We further use the variability measure to verify the existing theoretical predictions on scale-free and real networks. The results show that the heterogeneous mean-field (HMF) prediction agrees very well with the numerical threshold, except the case that the networks are disassortative, in which the quenched mean-field (QMF) prediction is relatively close to the numerical threshold. The numerical method presented can effectively identify SIR epidemic thresholds on various networks, and could be extended to other dynamical processes such as information diffusion and behavior spreading. This work provides us a deep understanding of epidemic threshold and would promote further studies on phase transition of epidemic dynamics.

Models for disease propagation are at the core of our understanding about epidemic dynamics on complex networks.1,2 Two epidemic models of particular importance are the susceptible-infected-susceptible (SIS) and susceptible-infected-recovered (SIR) models.3 At each time step, an infected node can transmit a disease to each of its susceptible neighbors with probability β. At the same time, the infected nodes become susceptible again in the SIS model or recover in the SIR model with probability μ. In the SIS model, a critical value of the effective transmission rate λ=β/μ separates the absorbing phase with only healthy nodes from the active phase with a stationary density of infected nodes. Differently, no steady state is allowed in the SIR model, but a threshold still exists above which the final fraction of recovered nodes is finite.4 

The traditional theoretical study on the epidemic threshold of the SIS model is based on the heterogeneous mean-field (HMF) theory, which means that all the nodes within a given degree are considered to be statistically equivalent.5,6 According to the HMF theory, the epidemic threshold of SIS model7,8 is given by λcHMF=kk2, where k and k2 are the first and second moments of degree distribution P(k),9 respectively. As the quenched structure of the network and dynamical correlations between the state of adjacent nodes is neglected in the HMF theory,10 researchers have proposed an important improvement over the HMF theory—quenched mean-field (QMF) theory. The QMF theory fully preserves the actual quenched structure of the network described as its adjacency matrix, and the epidemic threshold is predicted to be11–13 

λcQMF=1ΛN,
(1)

where ΛN is the maximum eigenvalue of the adjacency matrix of a given network. Considering that the existing theories more or less have some limitations (e.g., the HMF theory neglects the quenched structure of the network; QMF theory ignores dynamical correlations14), some numerical methods such as the finite-size scaling analysis,15 susceptibility,16 and lifetime measure17 have been proposed to check the accuracies of different theoretical predictions for the SIS model. Among these existing methods, the relatively common one for a network with finite size N is the susceptibility measure

χ=Nρ2ρ2ρ,
(2)

where ρ denotes the outbreak size. In Ref. 18, the susceptibility measure has been shown to be very effective for identifying the SIS epidemic thresholds on various networks.

For another paradigmatic epidemic model, the SIR model, there have been a lot of theoretical studies on its epidemic threshold. The earliest theoretical study on the SIR epidemic threshold is based on the assumption of homogeneous mixing, showing that the SIR epidemic threshold is inversely proportional to the average connectivity k.3 At the HMF level,19 the epidemic threshold of SIR model takes the value

λcHMF=kk2k.
(3)

On networks with power-law scaling P(k)kγ,9,20 where γ is the degree exponent, the HMF approach predicts a vanishing threshold for scale-free networks with γ3 and the finite threshold for γ>3.4 Mapping the SIR model to a bond percolation process,21 the epidemic threshold coincides with the result of Eq. (3) for a SIR model with unit infection time (e.g., μ = 1), and when the infection times vary among infected nodes (e.g., μ<1), the epidemic threshold is given by4 

λc=kk22k.
(4)

According to the QMF theory,11 the epidemic threshold of the SIR model has the same expression as Eq. (1). However, the QMF result is even qualitatively not correct, as the QMF predicts a vanishing threshold for power-law distributed networks with γ>3 that is in conflict with the visually numerical results.22 

Although the numerical threshold of the SIS model has attracted much attention,15–18 the systematic study of the numerical identification of the SIR epidemic threshold is still insufficient. It is well known that the outbreak size becomes finite above the threshold λc.21 However, as the value of λ increases, the outbreak size continuously changes from an infinitesimal fraction to a finite fraction in the SIR model,23 and thus, it is difficult to determine the value of λ at which the outbreak size turns to be finite. To our knowledge, there has no effective numerical method for identifying the SIR epidemic threshold in previous studies. In this work, we perform extensive numerical simulations of the SIR model on networks with finite size, and present a numerical identification method by analyzing the peak of the epidemic variability24,25 (i.e., the maximal value of the epidemic variability) to identify the epidemic threshold. The effectiveness of the numerical measure is checked on random regular networks (RRN), where the HMF theory is exact. To get a deep understanding of the validity of the numerical method, we investigate the distribution of outbreaks sizes near the epidemic threshold. The robustness of the variability measure is confirmed by the analysis of cutoff hypothesis of the outbreak size distribution. We further employ the variability measure to verify the theoretical predictions on scale-free networks and real-world networks, where the results indicate that the HMF prediction agrees very well with the numerical threshold, except the case that the networks are disassortative, in which the QMF prediction is relatively close to the numerical threshold.

In this section, we give the detailed description of the simulation of SIR model, propose the numerical identification measure of the epidemic threshold, and make deep analysis of the effectiveness of the numerical measure. In the SIR model, each node can be one of three states which are the susceptible state, infected state, and recovered state, respectively. At the beginning, one node is randomly selected as the initial infected (i.e., seed), and all other nodes are susceptible. At each time step t, each susceptible node i becomes infected with probability 1(1β)ni if it has one or more infected neighbors, where ni is the number of its infected neighbors. At the same time, all infected nodes recover (or die) at rate μ and the recovered nodes acquire permanent immunity. Time increases by Δt=1, and the dynamical process terminates when all infected nodes are recovered. In this paper, μ is set to 1, unless otherwise specified.

The susceptibility measure can not only identify an effective SIS epidemic threshold,18 but can also be used to determine the critical point of the percolation process.26 Since the connection between SIR and bond percolation is made through the assimilation of the size of the percolating giant component with the final number of recovered individuals,21 we check the effectiveness of susceptibility measure χ for the SIR model on RRN, where all nodes have exactly the same degree k. On these networks, the HMF prediction λcHMF=1/(k1) (Ref. 5) is accurate for the SIR model. We compare the HMF prediction with the numerical threshold identified by the susceptibility measure λpχ in Fig. 1(a), where the result shows that the numerical threshold of the SIR model identified by χ is significantly larger than 1/(k1).

FIG. 1.

Comparison of theoretical thresholds with numerical thresholds on RRN. (a) The threshold λc vs. k, where N is set to 104. (b) The threshold λc vs. N, where k is set to 10. “Squares,” “circles,” and “triangleups” denote λpχ, λpΔ, and 1/(k1), respectively. Inset: Susceptibility χ and variability Δ as a function of λ. The results are averaged over 102×106 independent realizations on 102 networks.

FIG. 1.

Comparison of theoretical thresholds with numerical thresholds on RRN. (a) The threshold λc vs. k, where N is set to 104. (b) The threshold λc vs. N, where k is set to 10. “Squares,” “circles,” and “triangleups” denote λpχ, λpΔ, and 1/(k1), respectively. Inset: Susceptibility χ and variability Δ as a function of λ. The results are averaged over 102×106 independent realizations on 102 networks.

Close modal

Considering that the fluctuation of the outbreak size is large near the epidemic threshold, we try to identify the epidemic threshold by the variability measure24,25

Δ=ρ2ρ2ρ,
(5)

which is a standard measure to determine the critical point in equilibrium phase on magnetic system.27 The inset of Fig. 1(a) shows that the variability Δ exhibits a peak over a wide range of λ, so we estimate the epidemic threshold from the position of the peak of the variability λpΔ. On RRN with different values of k, we find that λpΔ is always consistent with the HMF prediction. When the degree k is given, we further consider the relationship between the epidemic threshold and the network size N in Fig. 1(b), where the numerical thresholds λpχ and λpΔ do not change with N. λpΔ is very close to the HMF prediction, while there is an obvious gap between λpχ and HMF prediction. From the above, we know that the variability Δ can identify an effective SIR epidemic threshold, while the epidemic threshold identified by the susceptibility χ is overestimated on RRN. Thus, a new problem has arisen: Why the variability Δ performs well but the susceptibility χ goes awry for the SIR model?

Next, we make an analysis of the effectiveness of the numerical identification measures above, by investigating the distribution of outbreak sizes which has the strong heterogeneity near the epidemic threshold. On RRN with k = 10, where the epidemic threshold λc=1/(k1)=1/9, Fig. 2(a) shows the distribution of outbreak sizes near λc. The outbreak sizes follow approximately an exponential distribution at λ=0.1 that is smaller than λc. At λ=λc, the outbreak sizes follow a power-law distribution P(ρ)ρα with a cutoff at some values, where α1.5.28–30 Since the disease may die out quickly or infect a subset of nodes when λ>λc, the distribution of outbreak sizes is bimodal,31,32 with two peaks occurring at ρ=1/N and ρ0.2 for λ=0.12, respectively.

FIG. 2.

Analysis of effectiveness and robustness of numerical methods on RRN. (a) Numerical distribution of outbreak sizes in SIR model for λ=0.10 (circles), λ=1/9 (triangles), and λ=0.12 (squares), where blue solid, red short dash, and black dot lines represent the theoretical distributions given by Eq. (A4). (b) χ vs. λ, where only the small outbreak sizes with ρrc are considered to numerically count χ. (c) Δ vs. λ, where the lump is assumed to be located at rc in the theoretical distribution diagram when λ>λc. In (b) and (c), “triangles,” “circles,” and “diamonds” denote cutoff values rc = 0.05, 0.2 and 0.4, respectively. The parameters are chosen as N=104 and k = 10. The results are averaged over 102×106 independent realizations on 102 networks.

FIG. 2.

Analysis of effectiveness and robustness of numerical methods on RRN. (a) Numerical distribution of outbreak sizes in SIR model for λ=0.10 (circles), λ=1/9 (triangles), and λ=0.12 (squares), where blue solid, red short dash, and black dot lines represent the theoretical distributions given by Eq. (A4). (b) χ vs. λ, where only the small outbreak sizes with ρrc are considered to numerically count χ. (c) Δ vs. λ, where the lump is assumed to be located at rc in the theoretical distribution diagram when λ>λc. In (b) and (c), “triangles,” “circles,” and “diamonds” denote cutoff values rc = 0.05, 0.2 and 0.4, respectively. The parameters are chosen as N=104 and k = 10. The results are averaged over 102×106 independent realizations on 102 networks.

Close modal

Moreover, the theoretical distribution of outbreak sizes (see Appendix for details) is compared with the results obtained by numerical simulations in Fig. 2(a), where the theoretical probability from Eq. (A4) is consistent with the numerical results for relatively small outbreak sizes (ρ<0.05). At the epidemic threshold, the theoretical results also show that the outbreak sizes obey a power-law distribution with the exponent of about −1.5. When λ>λc, some large outbreak sizes constitute a lump in the numerical scattergram, but the probability of large outbreak sizes cannot be solved from Eq. (A4). We thus speculate that the non-ignorable lump may affect the numerical identification of SIR epidemic threshold for the susceptibility measure.

To verify the speculation, Fig. 2(b) investigates the effectiveness of the susceptibility measures under some cutoff hypotheses. We set the cutoff value of the outbreak size as rc, which means that only the outbreak sizes with ρrc are used to numerically count the susceptibility χ under the cutoff hypothesis. Three kinds of rc are considered, where rc=0.05 corresponds to the maximum value of small outbreak size before the lump appears in the numerical scattergram, rc=0.2 means that the numerical scattergram consists of a part of the lump, and rc=0.4 means that there is a complete lump in the numerical scattergram. As shown in Fig. 2(b), the susceptibility measure can indeed give a quite effective estimate of the SIR epidemic threshold when the whole lump is ignored (i.e., rc=0.05). With the increase of rc, the position of peak value of susceptibility χ gradually shifts to the right for large outbreak sizes are considered. This shows that the susceptibility χ loses its effectiveness in identifying the SIR epidemic threshold due to the existence of the lump.

In Fig. 2(c), the robustness of the variability Δ is further checked in theory. As the numerical distribution of the large outbreak sizes is concentrated, we assume that there is a lump located at rc with P(rc)=1Σρ<rcP(ρ) in the theoretical probability distribution diagram of outbreak sizes while λ>λc. Based on such theoretical distribution, we plot the variability measure as a function of λ for different values of rc in Fig. 2(c). Since the variability Δ measures the heterogeneity of the outbreak size distribution, which is strongest at the epidemic threshold,28–30 the peak position of the variability measure does not change with the position of the lump.

From the above analysis, we can conclude that the variability Δ is effective in identifying the epidemic threshold of SIR model, while the bimodal distribution of outbreak sizes for λ>λc leads to the overestimation of the SIR epidemic threshold when using the susceptibility χ.

We further verify the theoretical predictions on scale-free and real networks, by comparing them with the numerical thresholds from the variability Δ.

We build scale-free networks (SFN) with degree distribution P(k)kγ based on the configuration model.9 The so-called structural cutoff kmaxN1/2 and natural cutoff kmaxN1/γ1 (Ref. 33) are considered to constrain the maximum possible degree kmax on SFN. Differently, the degree-degree correlations vanish on scale-free networks with structural off, while the disassortative degree-degree correlations exist when γ<3 for scale-free networks with natural cutoff,33 because high degree vertices connect preferably to low degree ones in this case. We consider the SIR model on SFN with structural cutoff in Figs. 3(a) and 3(c), where the SIR epidemic threshold increases monotonically with the degree exponent γ, and the variation of epidemic threshold with network size N is approximately linear in logarithmic relationship.16 The HMF prediction λcHMF is very close to the numerical threshold λpΔ, while there is an obvious difference between the QMF prediction λcQMF and λpΔ.

FIG. 3.

Comparison of theoretical thresholds with numerical thresholds on SFN. The threshold λc vs. γ on SFN with structural cutoff (a) and natural cutoff (b), where N is set to 104. λc vs. N on SFN with structural cutoff (c) and natural cutoff (d), where solid and empty symbols denote γ=2.25 and 3.50, respectively. “Squares,” “circles,” and “triangles” denote λcQMF,λcHMF, and λpΔ, respectively. The results are averaged over 102×106 independent realizations on 102 networks.

FIG. 3.

Comparison of theoretical thresholds with numerical thresholds on SFN. The threshold λc vs. γ on SFN with structural cutoff (a) and natural cutoff (b), where N is set to 104. λc vs. N on SFN with structural cutoff (c) and natural cutoff (d), where solid and empty symbols denote γ=2.25 and 3.50, respectively. “Squares,” “circles,” and “triangles” denote λcQMF,λcHMF, and λpΔ, respectively. The results are averaged over 102×106 independent realizations on 102 networks.

Close modal

The SFN with natural cutoff are considered in Figs. 3(b) and 3(d), where the variations of epidemic threshold with γ and N are similar to the results on SFN with structural cutoff. When γ>3, the HMF prediction is close to the numerical threshold, while there is a gap between the QMF prediction and the numerical threshold. Since the disassortative degree-degree correlations exist for γ<3, there is a slight difference between λcHMF and λpΔ. In Fig. 3(d), the distinction between λcHMF and λpΔ becomes large with the increase of N for γ=2.25, while in such case the QMF prediction is always close to the numerical threshold since the principle eigenvector is delocalized when 2<γ5/2.34 It can been seen from the above analysis that the numerical method presented provides quantitive indexes for the observations in Ref. 22, where the HMF theory is relatively accurate for the SIR model.

To check the performance of the variability Δ on real-world networks,35–42 Fig. 4 depicts Δ as a function of λ for Hamsterster full network,35 which contains friendships and family links between users of the website hamsterster.com, and Facebook (NIPS) network,41 which contains Facebook user-user friendships. The numerical results intuitively show that the variability Δ always reaches a maximum value near the value of λ above which the outbreak size ρ becomes finite. The theoretical predictions of the HMF theory and of the QMF theory are quite close to the numerical threshold identified by Δ on Hamsterster full network, but they become poor on Facebook (NIPS) network. Considering the difference in the behaviors of the theoretical predictions on the two networks above, the detailed comparisons between the numerical and theoretical thresholds on other real networks are presented in Table I. The results indicate that although the HMF prediction and the numerical threshold λpΔ(SIR) are nearly the same for assortative networks, there is an obvious difference between them for the networks showing significant disassortative mixing. The QMF prediction is relatively worse than the HMF prediction for assortative networks, but the former is close to λpΔ(SIR) for some disassortative networks (e.g., Router views,37 CAIDI,37 and email contacts42). These findings numerically show the accuracies of existing theoretical predictions on real-world networks.

FIG. 4.

Variability Δ and outbreak size ρ as a function of λ on Hamsterster full network (a) and Facebook (NIPS) network (b). “Triangleup” and “triangledown” denote λcHMF=k/[k2k] and λcQMF=1/ΛN, respectively. The variability Δ for each λ is normalized by the maximal variability Δmax. The results are averaged over 106 independent realizations on each network.

FIG. 4.

Variability Δ and outbreak size ρ as a function of λ on Hamsterster full network (a) and Facebook (NIPS) network (b). “Triangleup” and “triangledown” denote λcHMF=k/[k2k] and λcQMF=1/ΛN, respectively. The variability Δ for each λ is normalized by the maximal variability Δmax. The results are averaged over 106 independent realizations on each network.

Close modal
TABLE I.

Characteristics and epidemic thresholds of real-world networks. N is the network size, kmax is the maximum degree, r is the correlation coefficient of the degrees, λcHMF and λcQMF are the HMF and QMF predictions for SIR, respectively, λpΔ(SIR) denotes the numerical threshold of SIR identified by Δ, and λpΔ(SIS) and λpχ(SIS) represent the numerical thresholds of SIS identified by Δ and χ, respectively.

NetworkNkmaxrλcHMFλcQMFλpΔ(SIR)λpΔ(SIS)λpχ(SIS)
Hamsterster full35  2000 273 0.023 0.023 0.020 0.023 0.025 0.025 
Brightkite36  56739 1134 0.010 0.016 0.010 0.014 0.012 0.012 
arXiv astro-ph37  17903 504 0.201 0.015 0.011 0.012 0.012 0.012 
Pretty good privacy38  10680 206 0.239 0.056 0.024 0.053 0.033 0.033 
US power grid39  4941 19 0.003 0.348 0.134 0.446 0.261 0.264 
Euroroad40  1039 10 0.090 0.479 0.249 0.498 0.331 0.331 
Facebook (NIPS)41  2888 769 −0.668 0.004 0.036 0.075 0.079 0.077 
Route views37  6474 1458 −0.182 0.006 0.022 0.037 0.034 0.036 
CAIDA37  26475 2628 −0.195 0.004 0.014 0.019 0.019 0.019 
email contacts42  12625 576 −0.387 0.009 0.02 0.027 0.024 0.025 
NetworkNkmaxrλcHMFλcQMFλpΔ(SIR)λpΔ(SIS)λpχ(SIS)
Hamsterster full35  2000 273 0.023 0.023 0.020 0.023 0.025 0.025 
Brightkite36  56739 1134 0.010 0.016 0.010 0.014 0.012 0.012 
arXiv astro-ph37  17903 504 0.201 0.015 0.011 0.012 0.012 0.012 
Pretty good privacy38  10680 206 0.239 0.056 0.024 0.053 0.033 0.033 
US power grid39  4941 19 0.003 0.348 0.134 0.446 0.261 0.264 
Euroroad40  1039 10 0.090 0.479 0.249 0.498 0.331 0.331 
Facebook (NIPS)41  2888 769 −0.668 0.004 0.036 0.075 0.079 0.077 
Route views37  6474 1458 −0.182 0.006 0.022 0.037 0.034 0.036 
CAIDA37  26475 2628 −0.195 0.004 0.014 0.019 0.019 0.019 
email contacts42  12625 576 −0.387 0.009 0.02 0.027 0.024 0.025 

In summary, we have studied the numerical identification of SIR epidemic threshold on complex networks with finite size. We have checked the effectiveness of the susceptibility measure for SIR on RRN, where the HMF is exact. The results showed an obvious gap between the numerical threshold identified by the susceptibility χ and the HMF prediction. Then, we proposed the numerical identification method by analyzing the peak of the epidemic variability, and found that the numerical threshold identified by the variability measure Δ agrees very well with the HMF prediction on RRN.

In order to get a deep understanding of the effectiveness of the two numerical measures above, we have analyzed the distribution of outbreak sizes near the epidemic threshold λc. The outbreak sizes follow approximately an exponential distribution when λ<λc. At the epidemic threshold, the outbreak sizes follow a power-law distribution with the exponent of about −1.5. When λ>λc, the numerical distribution of outbreak sizes is bimodal with two peaks occurring at ρ=1/N and O(1), respectively. The probability of small outbreak sizes in theory is consistent with that obtained by numerical simulations, but the probability of large outbreak sizes that constitute a lump in the numerical scattergram cannot be obtained theoretically. Based on the analysis of the cutoff hypothesis of the outbreak size distribution, we found that the susceptibility measure can give a fairly effective SIR epidemic threshold when the lump is ignored. Since the variability measure reflects the heterogeneity of the outbreak size distribution, it is always effective in identifying the epidemic threshold, where the distribution of outbreak sizes has the very strong heterogeneity.

We further employed the variability measure to verify the theoretical predictions on scale-free and real networks. The HMF prediction is close to the numerical threshold on most of the networks, but on SFN with natural cutoff and degree exponent γ<5/2, it becomes poor due to the existence of disassortative mixing. Similarly, the HMF prediction agrees well with the numerical method on real networks with assortative mixing, while it becomes very poor for disassortative networks, where the QMF prediction is relatively close to the numerical threshold. These findings provide quantitive indexes for the accuracies of existing theoretical predictions from the perspective of simulation.

As part of the discussion, we have considered the epidemic threshold for μ<1 in Fig. 5. The results on RRN and SFN all show that the numerical thresholds for μ=0.1 are a little larger than those for μ=0.5. As shown in the inset, μ0 leads to an epidemic threshold close to Eq. (4), while μ1 leads to an epidemic threshold close to Eq. (3). It should be pointed out that the numerical threshold of SIR model is inclined to the theoretical prediction λc=kk2+(μ2)k from the edge-based compartmental theory43,44 when 0<μ<1. These findings could be complementary to some existing results.4 

FIG. 5.

Comparison of theoretical thresholds with numerical thresholds for μ<1. The threshold λc vs. k for μ=0.1 (a) and μ=0.5 (b) on RRN. The threshold λc vs. γ for μ=0.1 (c) and μ=0.5 (d) on SFN, where solid and empty symbols denote SFN with structural cutoff kmaxN1/2 and natural cutoff kmaxN1/γ1, respectively. “Circles” and “squares” denote λcHMF and λpΔ, respectively. Inset: λc as a function of μ on RRN with N=104 and k = 6. The results are averaged over 102×106 independent realizations on 102 networks.

FIG. 5.

Comparison of theoretical thresholds with numerical thresholds for μ<1. The threshold λc vs. k for μ=0.1 (a) and μ=0.5 (b) on RRN. The threshold λc vs. γ for μ=0.1 (c) and μ=0.5 (d) on SFN, where solid and empty symbols denote SFN with structural cutoff kmaxN1/2 and natural cutoff kmaxN1/γ1, respectively. “Circles” and “squares” denote λcHMF and λpΔ, respectively. Inset: λc as a function of μ on RRN with N=104 and k = 6. The results are averaged over 102×106 independent realizations on 102 networks.

Close modal

Moreover, we have tried applying the variability measure to the identification of the SIS epidemic threshold. As shown in Fig. 6 and Table I, the numerical threshold λpΔ from the variability measure agrees very well with the threshold λpχ identified by the susceptibility measure, whose validity for the SIS model has been confirmed in Ref. 18. This shows that the variability measure can also provide an effective estimate of the SIS epidemic threshold.

FIG. 6.

Comparison of theoretical thresholds with numerical thresholds for SIS model on RRN. (a) The threshold λc vs. k, where N is set to 104. (b) The threshold λc vs. N, where k is set to 10. “Squares” and “circles” denote λpχ and λpΔ, respectively. The results are averaged over 102×106 independent realizations on 102 networks.

FIG. 6.

Comparison of theoretical thresholds with numerical thresholds for SIS model on RRN. (a) The threshold λc vs. k, where N is set to 104. (b) The threshold λc vs. N, where k is set to 10. “Squares” and “circles” denote λpχ and λpΔ, respectively. The results are averaged over 102×106 independent realizations on 102 networks.

Close modal

We have put forward a numerical method for identifying the epidemic threshold for SIR model, which is also suitable for the SIS model. This method can effectively identify epidemic thresholds on various networks, and could be extended to other dynamical processes such as information diffusion and behavior spreading. Further work should be done to check the effectiveness of this method on more complicated networks (e.g., temporal networks45 and multilayer networks46). Besides, the accurate analytic approximation of the epidemic threshold for general networks remains an important problem. This work helps to verify theoretical analysis of epidemic threshold and would promote further studies on phase transition of epidemic dynamics.

This work was partially supported by National Natural Science Foundation of China (Grant Nos. 11105025 and 91324002), China Postdoctoral Science Special Foundation (Grant No. 2012T50711), the Program of Outstanding Ph. D. Candidate in Academic Research by UESTC (Grant No. YXBSZC20131033), and Open Foundation of State key Laboratory of Networking and Switching Technology (Beijing University of Posts and Telecommunications) (SKLNST-2013-1-18). Y. Do was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2013R1A1A2010067).

For the case of the SIR model and similar models with no steady-state, the static properties (e.g., the final outbreak size and the epidemic threshold) of the epidemic outbreak can be mapped into a suitable bond percolation problem. In this framework, the distribution of occupied cluster sizes is related to the distribution of outbreak sizes. To get the distribution of small outbreak sizes in the SIR model with a fixed value of λ when recovery rate μ = 1, we will present the derivation of the distribution of small occupied cluster sizes in bond percolation with bond occupation probability λ.21 

After the percolation process on a general network with arbitrary degree distribution pk, the average degree of the occupied network A1, which composes of vertices and occupied edges, is kT=λk, where k is the average degree of the original network A0. And, the size distribution of the small subgraphs of network A1 is

πs=kT(s1)![ds2dzs2[g1(z)]s]z=0,
(A1)

where s is the small subgraphs size and g1(z) is the generating function of the excess degree of network A1. In addition, the generating function of degree distribution of A1 is

g0(z)=k=0pk(1λ+zλ)k,

and we thus have

g1(z)=g0(z)g0(1).

In a random regular network, which has an unique degree k with pk = 1, we can easily obtain that

g0(z)=[1+(z1)λ]k,
(A2)

and

g1(z)=[1+(z1)λ]k1.
(A3)

Substituting Eq. (A3) into Eq. (A1), we can obtain the distribution of small outbreak sizes

πs=kΓ(a2)Γ(a0)Γ(a1)λs1(1λ)s(k1)(s2),
(A4)

where Γ(x+1)=x!,a0=(s2),a1=s(k1)(s1), and a2=s(k1)1.

1.
A.
Barrat
,
M.
Barthélemy
, and
A.
Vespignani
,
Dynamical Processes on Complex Networks
(
Cambridge University Press
,
New York
,
2008
).
2.
A.
Vespignani
,
Nat. Phys.
8
,
32
(
2012
).
3.
R. M.
Anderson
and
R. M.
May
,
Infections Diseases in Humans
(
Oxford University Press
,
Oxford
,
1992
).
4.
R.
Pastor-Satorras
,
C.
Castellano
,
P. V.
Mieghem
, and
A.
Vespignani
, e-print arXiv:1408.2701v1 (
2014
).
5.
S. N.
Dorogovtsev
,
A. V.
Goltsev
, and
J. F. F.
Mendes
,
Rev. Mod. Phys.
80
,
1275
(
2008
).
6.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. E
63
,
066117
(
2001
).
7.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. Lett.
86
,
3200
(
2001
).
8.
M.
Boguñá
and
R.
Pastor-Satorras
,
Phys. Rev. E
66
,
047104
(
2002
).
9.
M. E. J.
Newman
,
Networks: An Introduction
(
Oxford University Press
,
Oxford
,
2010
).
10.
O.
Givan
,
N.
Schwartz
,
A.
Cygelberg
, and
L.
Stone
,
J. Theor. Biol.
288
,
21
(
2011
).
11.
D.
Chakrabarti
,
Y.
Wang
,
C.
Wang
,
J.
Leskovec
, and
C.
Faloutsos
,
ACM Trans. Inf. Syst. Secur.
10
,
1
(
2008
).
12.
P.
Van Mieghem
,
J.
Omic
, and
R.
Kooij
,
IEEE/ACM Trans. Network
17
,
1
(
2009
).
13.
S.
Gómez
,
A.
Arenas
,
J.
Borge-Holthoefer
,
S.
Meloni
, and
Y.
Moreno
,
Europhys. Lett.
89
,
38009
(
2010
).
14.
J. P.
Gleeson
,
Phys. Rev. Lett.
107
,
068701
(
2011
).
15.
J.
Marro
and
R.
Dickman
,
Nonequilibrium Phase Transitions in Lattice Models
(
Cambridge University Press
,
Cambridge
,
1999
).
16.
K.
Binder
and
D. W.
Heermann
,
Monte Carlo Simulation in Statistical Physics
, 5th ed. (
Springer-Verlag
,
Berlin
,
2010
).
17.
M.
Boguñá
,
C.
Castellano
, and
R.
Pastor-Satorras
,
Phys. Rev. Lett.
111
,
068701
(
2013
).
18.
S. C.
Ferreira
,
C.
Castellano
, and
R.
Pastor-Satorras
,
Phys. Rev. E
86
,
041125
(
2012
).
19.
M.
Barthélemy
,
A.
Barrat
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
Phys. Rev. Lett.
92
,
178701
(
2004
).
20.
R.
Albert
and
A.-L.
Barabási
,
Rev. Mod. Phys.
74
,
47
(
2002
).
21.
M. E. J.
Newman
,
Phys. Rev. E
66
,
016128
(
2002
).
22.
C.
Castellano
and
R.
Pastor-Satorras
,
Phys. Rev. Lett.
105
,
218701
(
2010
).
23.
C.
Castellano
and
R.
Pastor-Satorras
,
Sci. Rep.
2
,
371
(
2012
).
24.
P.
Crépey
,
F. P.
Alvarez
, and
M.
Barthélemy
,
Phys. Rev. E
73
,
046131
(
2006
).
25.
P.
Shu
,
M.
Tang
,
K.
Gong
, and
Y.
Liu
,
Chaos
22
,
043124
(
2012
).
26.
N.
Bastas
,
K.
Kosmidis
, and
P.
Argyrakis
,
Phys. Rev. E
84
,
066112
(
2011
).
27.
S. C.
Ferreira
,
R. S.
Ferreira
,
C.
Castellano
, and
R.
Pastor-Satorras
,
Phys. Rev. E
84
,
066102
(
2011
).
28.
E.
Ben-Naim
and
P. L.
Krapivsky
,
Phys. Rev. E
69
,
050901(R)
(
2004
).
29.
E.
Ben-Naim
and
P. L.
Krapivsky
,
Eur. Phys. J. B
85
,
145
(
2012
).
30.
D. A.
Kessler
and
N. M.
Shnerb
,
Phys. Rev. E
76
,
010901(R)
(
2007
).
31.
D. H.
Zanette
,
Phys. Rev. E
64
,
050901(R)
(
2001
).
32.
A.
Khalleque
and
P.
Sen
,
J. Phys. A: Math. Theor.
46
,
095007
(
2013
).
33.
M.
Boguñá
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
Eur. Phys. J. B
38
,
205
(
2004
).
34.
A. V.
Goltsev
,
S. N.
Dorogovtsev
,
J. G.
Oliveira
, and
J. F. F.
Mendes
,
Phys. Rev. Lett.
109
,
128702
(
2012
).
35.
J.
Kunegis
, see http://konect.uni-koblenz.de/networks/petster-hamster for Hamsterster full (accessed 5th March
2014
).
36.
E.
Cho
,
S. A.
Myers
, and
J.
Leskovec
, in
Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining
(
San Diego, California
,
2011
), pp.
1082
1090
.
37.
J.
Leskovec
,
J.
Kleinberg
, and
C.
Faloutsos
,
ACM Trans. Knowl. Discovery Data
1
,
1
(
2007
).
38.
M.
Boguñá
,
R.
Pastor-Satorras
,
A.
Diaz-Guilera
, and
A.
Arenas
,
Phys. Rev. E
70
,
056122
(
2004
).
39.
D. J.
Watts
and
S. H.
Strogatz
,
Nature
393
,
440
(
1998
).
40.
L.
Šubelj
and
M.
Bajec
,
Eur. Phys. J. B
81
,
353
(
2011
).
41.
J.
McAuley
and
J.
Leskovec
, in
Advances in Neural Information Processing Systems
(
Lake Tahoe
,
Nevada
,
2012
), pp.
539
547
.
42.
M.
Kitsak
,
L. K.
Gallos
,
S.
Havlin
,
F.
Liljeros
,
L.
Muchnik
,
H. E.
Stanley
, and
H. A.
Makse
,
Nat. Phys.
6
,
888
(
2010
).
43.
44.
W.
Wang
,
M.
Tang
,
H. F.
Zhang
,
H.
Gao
,
Y.
Do
, and
Z. H.
Liu
,
Phys. Rev. E
90
,
042803
(
2014
).
45.
P.
Holme
and
J.
Saramäki
,
Phys. Rep.
519
,
97
(
2012
).
46.
S.
Boccaletti
,
G.
Bianconi
,
R.
Criado
,
C. I.
del Genio
,
J.
Gómez-Gardeñes
,
M.
Romance
,
I.
Sendiña-Nadai
,
Z.
Wang
, and
M.
Zanin
,
Phys. Rep.
544
,
1
(
2014
).