Pair quenched mean-field theory for the susceptible-infected-susceptible model on complex networks
We present a quenched mean-field (QMF) theory for the dynamics of the susceptible-infected-susceptible (SIS) epidemic model on complex networks where dynamical correlations between connected vertices are taken into account by means of a pair approximation. We present analytical expressions of the epidemic thresholds in the star and wheel graphs and in random regular networks. For random networks with a power law degree distribution, the thresholds are numerically determined via an eigenvalue problem. The pair and one-vertex QMF theories yield the same scaling for the thresholds as functions of the network size. However, comparisons with quasi-stationary simulations of the SIS dynamics on large networks show that the former is quantitatively much more accurate than the latter. Our results demonstrate the central role played by dynamical correlations on the epidemic spreading and introduce an efficient way to theoretically access the thresholds of very large networks that can be extended to dynamical processes in general.
Networks and genealogical trees Critical point phenomena Dynamics of social systems Nonequilibrium and irreversible thermodynamics
The propagation of an epidemic is an ultimately important issue in a broad collection of systems ranging from the disease dissemination within a population [1, 2] to the virus spreading throughout a computer network [3, 4, 5] among a plenty of examples . The description of the epidemic dynamics has drawn the attention of biologist and mathematicians since a long time  and, more recently, of the statistical physics [3, 4] and computer science communities [7, 8]. The simplest epidemic model is the susceptible-infected-susceptible (SIS) model where individuals in a population can be only in one of two states: infected or susceptible (healthy). Infected individuals become spontaneously healthy at rate (this choice fixes the time scale) while the susceptible ones are infected at rate , where is the number of infected contacts of the individual . The SIS dynamics exhibits a phase transition between a disease-free (absorbing) state and an active stationary phase where a fraction of the population is infected. These regimes are separated by an epidemic threshold , the core issue of many recent works about epidemic spreading [4, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].
The former investigations of the epidemic models relied on the hypothesis of a homogeneous mixing of the population , neglecting the highly heterogeneous structure of the contact network inherent to real systems . Indeed, many bio, socio and technological systems are characterized by heavy tailed distributions of the number of contacts of an individual (the vertex degree, in the network language), for which homogeneity hypothesis is severely violated [20, 6, 21]. Complex networks are, in fact, a framework where the heterogeneity of the contacts can be naturally afforded .
The heterogeneous mean-field (HMF) theory is a benchmark for dynamical process running on the top of complex networks [21, 22]. In a HMF theory, dynamical quantities, as the density of infected individuals in the SIS model, depend only of the vertex degree and do not of their specific location in the network. The HMF epidemic threshold for the SIS dynamics in undirected and uncorrelated networks is given by 
where and is the degree distribution defined as the probability that a randomly chosen vertex has degree . Equation (1) has strong implications since several real networks have a power law degree distribution with exponent in the range . For these distributions, the second moment diverges in the limit of infinite sizes implying a vanishing threshold for the SIS model.
An alternative approach, called of quenched mean-field (QMF) theory , has been developed in parallel to the HMF theory [7, 10, 11, 8, 17, 19, 18]. The QMF theory explicitly takes into account the actual connectivity of the network through its adjacency matrix. The epidemic threshold of the SIS model in a QMF approach is 
where is the largest eigenvalue of the adjacency matrix. The noteworthy differences and similarities between HMF and QMF theories were realized in Ref. . The central point is that diverges for increasing networks with power law degree distributions even when remains finite . Both theories predict vanishing thresholds for despite of different scaling for . However, while HMF predicts a finite threshold for networks with , QMF still predicts a vanishing threshold . Very recently, semi-analytic methods including local and long-range dynamical fluctuations [24, 25] have been proposed as alternatives to HMF and QMF theories.
A numerical investigation of the thresholds of the SIS model on several networks was recently done in Ref. . It was shown that the QMF theory is an improvement of HMF theory but it is still an approximation. In fact, the assumption where the probability that a vertex is infected does not depend of the states of its neighbors is used in both approaches. This approximation neglects all dynamical correlations between vertices, which possibly contribute to the threshold value. Interestingly, the thresholds of some models with absorbing configurations taking place in highly heterogeneous networks are surprisingly well described by a simple homogeneous pair approximation [26, 27, 28], the simplest mean-field theory that considers dynamical correlations between vertices.
In the present work, we develop an extension of the QMF theory for the SIS model, the pair QMF theory, using a heterogeneous pair mean-field approximation. Analogously to the one-vertex QMF theory, our perturbative theory includes the actual connectivity of the network through its adjacency matrix. Analytical expressions are presented for the random regular networks, star and wheel graphs. We also investigated large random networks with a power law degree distribution. The thresholds obtained in pair and one-vertex QMF theories have the same scaling with the system size but the pair QMF theory is quantitatively much more accurate than the one-vertex theory when compared with simulations.
To develop the pair QMF theory, we introduce the following notation: is probability that the vertex is in the state ; is probability that the vertices and are in states and , respectively; is the extension to three vertices; and so forth. The infected state is represented by and the susceptible one by . We introduce the variables and, consequently, , , , , and . Obviously we have that , , and . The following relations hold for any pair of vertices
The dynamical equation for the vertex is
where is the adjacency matrix that, for unweighted and undirected networks, assumes if vertices and are connected and , otherwise. Using a one-vertex approximation where the joint probability is factorized as , one obtains
Performing a linear stability analysis around the trivial fixed point , one has where the Jacobian matrix is , being the Kronecker delta symbol. The transition point is defined when the absorbing phase becomes unstable or, equivalently, when the largest eigenvalue of is null . So, one directly finds the threshold given by Eq. (2).
The dynamical equation for , considering a pair of connected vertices , is given by
where represents the neighborhood of the vertex .
The first three terms represents the reactions involving only the pair ,
that can create or destroy the configuration
Equations (4) and (6) cannot be solved due to the triplets. However, the dynamical equations for triplets will depend on quadruplets, and so forth. So, the hierarchy of clusters must be broken in some point to obtain an approximated solution. In the present work, we apply the standard pair-approximation [30, 31]
and the adjacency matrix in equation (6) to obtain
We now perform a linear stability analysis of Eq. (8) around the fixed point to find
Again, the critical point is obtained when the largest eigenvalue of is null. Equation (11) is the central result of our work. Even though we do not present a closed expression for the threshold in an arbitrary network, we have obtained analytical solutions of transition points for simple networks. These solutions are very important to test the consistency of the results and to unveil the basic mechanisms that sustain an epidemic phase [13, 14]. For the general case, the critical point can be obtained solving Eq. (11) numerically.
An approach similar to ours was developed in Ref. , where a different approximation was used to split the joint probabilities in cluster approximation: instead of the standard pair approximation, Eq. (7), which has been generally accepted in the nonequilibrium statistical community as the most reliable approach .
Let us start the investigation of Eq. (11) for a simple homogeneous network, the random regular network (RRN) illustrated in Fig. 1(a). In this network, all vertices have the same connectivity , , while the connections are done at random, avoiding both self and multiple connections. The largest eigenvalue of is given by
where we used the largest eigenvalue of the given by . The threshold is then obtained as
that corresponds to the threshold of a simple homogeneous pair approximation . Actually, simulations of the SIS model on RRNs reported in Ref.  showed that the thresholds are, in fact, much closer to the homogeneous pair approximation than to the standard QMF theory, a fact captured by the pair QMF theory.
We performed simulations of the SIS model using the quasi-stationary method for dynamical processes with absorbing states suitably adapted for networks . Details of the implementation are given in Ref. . The thresholds in finite networks can be estimated using the peaks of the susceptibility defined, in terms of the density of occupied vertices, as 
Figure 2 shows the susceptibility against infection rate for exhibiting a sharp peak that asymptotically converges to a position that is only 1% above the theoretical value predicted by the pair QMF theory, . So, the pair QMF theory improves a lot the estimate of the thresholds in relation to the simple QMF theory , which errs approximately 22% the peak position for . Moreover, the larger the vertex degree the better the pair approximation. For , the peak converges to that is relatively much farther from the pair QMF threshold, , than in the case . The better accuracy of the pair QMF theory for larger is intuitive since the average distance among vertices decreases as the average degree increases making the mean-field premise a more credible hypothesis.
As an example of simple heterogeneous network, we consider star graph defined as a hub, , connected to leaves, , of degree , as shown in Fig. 1(b). This structure plays a central role in the SIS dynamics since the epidemic activity may be sustained by the sub-graph composed by the most connected vertex and its neighborhood for networks with degree exponent . The adjacency matrix of a star graph is for and , otherwise. Therefore, the elements of the Jacobian are given by
for the diagonal and
for . Thus, the eigenvalue equations for are
which is larger than the prediction of both the standard QMF theory, , and the pair approximation developed in Ref. , . Our result explains very well the pre-factor larger than 1 observed in simulations of the SIS model on star graphs .
Figure 3 shows the susceptibility against the infection rate exhibiting a rounded peak characteristic of SIS model on star graphs . The inset of Fig. 3 compares the thresholds of the simple and the pair QMF theories with the positions of the susceptibility peaks for different network sizes. The last one is assumed as the real threshold of finite networks . An excellent agreement between simulations and pair QMF results is obtained implying a remarkable improvement in relation to the standard QMF approach. However, it is important to emphasize that the pair QMF is still an approximation and that a small discrepancy (less than 5%) is observed for the largest investigated system.
A wheel graph is defined as a regular chain of vertices and periodic boundary conditions and a central vertex connected to all vertices of the chain, Fig. 1(c). This network is interesting in the context of epidemic spreading due to high clustering coefficient that contrasts with the null clustering coefficient of the star graph. Following steps similar to those performed for the star graph, one finds
This threshold is essentially the same obtained for the star graph, which was confirmed in simulations (data not shown). Triangles, characteristic of clustered networks, enhance dynamical correlations and, therefore, are expected to interfere in the thresholds of dynamical processes in general [34, 35]. So, the clustering-independence observed for wheel graphs must be a characteristic of transitions ruled by the star sub-graph centered at the most connected vertex .
For arbitrary random networks, the largest eigenvalue of Eq. (11) can be numerically determined . We considered the uncorrelated configuration model (UCM)  with minimal vertex degree fixed to and a structural cutoff in the degree distribution to investigate SIS dynamics via both Eq. (11) and quasi-stationary simulations. It is known that fluctuations of the degree of the most connected vertex drastically change the position of the threshold for . Therefore, we did our analysis for networks with [13, 16], in order to have representative results from a single sample.
For the range , the QMF and HMF theories are essentially equal  and both theories agree with thresholds estimated from the peaks of the susceptibility curves . We investigated networks with and verified that the pair QMF is very close to QMF and, consequently, to HMF theories, being the pair QMF the one closer to the simulation results. Figure 4 compares the three theories with the peak of the susceptibility (top inset). One can see that pair QMF theory fits better the thresholds for small networks (bottom inset).
Figure 5 shows the thresholds for . In this case, the pair QMF consists in a great improvement in relation to the standard QMF (bottom inset of Fig. 5). The difference is less than 8% while standard QMF errs in approximately 30%. It was proposed that the star sub-graph composed by the most connected vertex and its neighbors is responsible by sustaining activity in random network with . The threshold for the corresponding star sub-graph, , is also included in Fig. 5. One can see that this threshold converges to the pair QMF as network increases confirming that the star sub-graph is actually the structure responsible by sustaining the epidemics in the entire network. Notice that HMF theory does not capture the scaling of threshold against the network size .
The most drastic differences among theories appear for , as illustrated in Fig. 6 where we show the threshold against size for . While HMF theory predicts a finite threshold both simple and pair QMF theories yield asymptotically vanishing thresholds. The vs. curves have a single sharp peak for small sizes but a secondary rounded peak at small emerges for very large sizes. Such doubly peaked susceptibility was interpreted as the competition between two mechanisms triggering the epidemic spreading in the network: The activity in the star sub-graph centered at the most connected vertex against that in the most densely connected component of the network . These competing mechanisms are also associated to the localization/delocalization of the epidemic activity in networks . The peaks at small are well fitted by the pair QMF theory (less than 10% of difference against 40% for QMF theory). The peak at larger is not captured in our theoretical approach. Notice that the threshold of the pair QMF theory quickly converges to that of a star sub-graph . Non-local mean-field approaches are potential candidates to explain the rightmost peaks 
In conclusion, we have investigated the epidemic spreading of the SIS model performing a mean-field pair approximation. Our theory is an extension of the quenched mean-field (QMF) theory [7, 13] that includes the actual network connectivity. The dynamical correlations are introduced by means of a pair approximation . Analytical expressions for the epidemic thresholds are presented for some simple networks while the thresholds for the most interesting case of random networks with a power law degree distribution are obtained from the numerical solution of an eigenvalue problem.
We compared our pair theory with the one-vertex QMF and heterogeneous mean-field (HMF) theories and with the peaks of susceptibility as a function of infection rate obtained in quasi-stationary simulations . We have shown that the thresholds in the pair QMF theory as a function of the network size scale as in the standard QMF theory. However, the pair QMF thresholds are quantitatively much closer to the simulation results than those of standard QMF.
Our theoretical approach represents an important advance in relation to other improvements of the QMF theory, which are that limited to small networks [32, 38] whereas we were able to investigate networks as large as . Despite of the considerable improvement when compared to QMF theory, our approach is still an approximation. Certainly, a -vertex theory with should enhance the accuracy of the threshold determination. Also, the critical exponents associated to the transition are still an open problem in our approach. The pair QMF theory can be extended to other dynamical processes taking place on the complex networks, for which pair approximations have exhibited great improvement in relation to one-site mean-field theories [26, 28, 27, 39].
This work was partially supported by the Brazilian agencies CNPq and FAPEMIG. ASM thanks the financial support from CAPES. Authors thank the enlightening discussions with C. Castellano and R. Pastor-Satorras.
- Spontaneous annihilation events and and the creation in vertex due to , .
- Creation events in or due another vertex , represented by transitions and , respectively, can also destroy/create a configuration .
- \NameAnderson R. M. May R. M. \BookInfectious diseases in humans (Oxford University Press, Oxford) 1992.
- \NameKeeling M. Rohani P. \BookModeling Infectious Diseases in Humans and Animals (Princeton University Press) 2008.
- \NamePastor-Satorras R., Vázquez A. Vespignani A. \REVIEWPhys. Rev. Lett. 872001258701.
- \NamePastor-Satorras R. Vespignani A. \REVIEWPhys. Rev. Lett. 8620013200.
- \NameLloyd A. L. May R. M. \REVIEWScience 29220011316.
- \NameNewman M. \BookNetworks: An Introduction (Oxford University Press, Inc., New York, NY, USA) 2010.
- \NameWang Y., Chakrabarti D., Wang C. Faloutsos C. \BookEpidemic spreading in real networks: An eigenvalue viewpoint in proc. of \Book22nd International Symposium on Reliable Distributed Systems (SRDS’03) (IEEE Computer Society, Los Alamitos, CA, USA) 2003 pp. 25–34.
- \NameGanesh A., Massoulié L. Towsley D. in Proceedings of \Book24th Annual Joint Conference of the IEEE Computer and Communications Societies (IEEE, Piscataway, NJ, 2005), Vol. 2, pp. 1455-1466.
- \NameChatterjee S. Durrett R. \REVIEWAnn. Probab. 3720092332.
- \NameChakrabarti D., Wang Y., Wang C., Leskovec J. Faloutsos C. \REVIEWACM Trans. Inf. Syst. Secur. 1020081.
- \NameDurrett R. \REVIEWProc. Natl. Acad. Sci. USA 10720104491.
- \NameNewman M. E. J. \REVIEWPhys. Rev. E 662002016128.
- \NameCastellano C. Pastor-Satorras R. \REVIEWPhys. Rev. Lett. 1052010218701.
- \NameCastellano C. Pastor-Satorras R. \REVIEWSci. Rep. 22012371.
- \NameVan Mieghem P., Omic J. Kooij R. \REVIEWIEEE ACM T. Network. 1720091 .
- \NameFerreira S. C., Castellano C. Pastor-Satorras R. \REVIEWPhys. Rev. E 862012041125.
- \NameMieghem P. V. \REVIEWEurophys. Lett. 97201248004.
- \NameGómez S., Arenas A., Borge-Holthoefer J., Meloni S. Moreno Y. \REVIEWEurophys. Lett. 89201038009.
- \NameGoltsev A. V., Dorogovtsev S. N., Oliveira J. G. Mendes J. F. F. \REVIEWPhys. Rev. Lett. 1092012128702.
- \NameAlbert R. Barabási A.-L. \REVIEWRev. Mod. Phys. 74200247.
- \NameBarrat A., Barthélemy M. Vespignani A. \BookDynamical Processes on Complex Networks (Cambridge University Press, Cambridge) 2008.
- \NameDorogovtsev S. N., Goltsev A. V. Mendes J. F. F. \REVIEWRev. Mod. Phys. 8020081275.
- \NameChung F., Lu L. Vu V. \REVIEWProc. Natl. Acad. Sci. USA 10020036313.
- \NameLee H. K., Shim P.-S. Noh J. D. \REVIEWPhys. Rev. E 872013062812.
- \NameBoguñá M., Castellano C. Pastor-Satorras R. \REVIEWPhys Rev. Lett.1112013068701.
- \NameFerreira S. C., Ferreira R. S., Castellano C. Pastor-Satorras R. \REVIEWPhys. Rev. E 842011066102.
- \NameJuhász R., Ódor G., Castellano C. Muñoz M. A. \REVIEWPhys. Rev. E 852012066125.
- \NameSander R. S., Ferreira S. C. Pastor-Satorras R. \REVIEWPhys. Rev. E 872013022820.
- \NameHilborn R. C. \BookChaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers (Oxford University Press, New York) 2000.
- \Nameben Avraham D. Köhler J. \REVIEWPhys. Rev. A 4519928358.
- \NameHenkel M., Hinrichsen H. Lübeck S. \BookNon-equilibrium phase transition: Absorbing Phase Transitions (Springer Verlag, Netherlands) 2008.
- \NameCator E. Van Mieghem P. \REVIEWPhys. Rev. E 852012056111.
- \NameFerreira S. C., Ferreira R. S. Pastor-Satorras R. \REVIEWPhys. Rev. E 832011066113.
- \NameRozhnova G. Nunes A. \REVIEWPhys. Rev. E 802009051915.
- \NameFerreira R. S. Ferreira S. C. \REVIEWpreprint arXiv:1307.6186 2013.
- \NamePress W., Teukolsky S., Vetterling W. Flannery B. \BookNumerical Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University Press) 2007.
- \NameCatanzaro M., Boguñá M. Pastor-Satorras R. \REVIEWPhys. Rev. E 712005027103.
- \NameGleeson J. P. \REVIEWPhys. Rev. Lett. 1072011068701.
- \NamePugliese E. Castellano C. \REVIEWEurophys. Lett. 88200958004.