Susceptible-infected-susceptible model on networks with eigenvector localization

TL;DR

This paper investigates how eigenvector localization in network adjacency matrices influences the SIS epidemic model, revealing a vanishing threshold and highlighting the role of highly connected nodes in sustaining infections.

Contribution

It introduces a critical condition for phase transition based on star subgraph lifespan and elucidates the impact of eigenvector localization on infection dynamics in networks.

Findings

Infection density depends on node degree and star lifespan.

Eigenvector localization enhances infection in hubs and suppresses it in low-degree nodes.

The average star lifespan divergence causes a vanishing epidemic threshold.

Abstract

It is a longstanding debate on the absence of threshold for susceptible-infected-susceptible (SIS) model on networks with finite second order moment of degree distribution. The eigenvector localization of the adjacency matrix for a network gives rise to the inactive Griffiths phase featuring slow decay of the activity localized around highly connected nodes due to the dynamical fluctuation. We show how it dramatically changes our understanding of SIS model, opening up new possibilities for the debate. We derive the critical condition for Griffiths to active phase transition: on average, an infected node can further infect another one in the characteristic lifespan of the star subgraph composed of the node and its nearest neighbors. The system approaches the critical point of avoiding the irreversible dynamical fluctuation and the trap of absorbing state. As a signature of the phase transition, the infection density of a node is not only proportional to its degree, but also proportional to the exponentially growing lifespan of the star. And the divergence of the average lifespan of the stars is responsible for the vanishing threshold in the thermodynamic limit. The eigenvector localization exponentially reinforces the infection of highly connected nodes, while it inversely suppresses the infection of small-degree nodes.

Equations

$$\frac{d{\rho }_i(t)}{dt}=-\bar{\delta }(k_i,\lambda ){\rho }_i(t)+(1-{\rho }_i(t))\sum _{j\ne i}\bar{\lambda }(d_{ij},\lambda ){\rho }_j(t),$$

$$\frac{d{\rho }_k(t)}{dt}=-\bar{\delta }(k,\lambda ){\rho }_k(t)+(1-{\rho }_k(t))\sum _{k^{\prime }}{\lambda }_{k{k}^{\prime }}{\rho }_{k^{\prime }}(t)NP(k^{\prime }),$$

$${\lambda }_{k{k}^{\prime }}\approx \bar{\lambda }(d_{k{k}^{\prime }},\lambda )=\lambda A_{k{k}^{\prime }}^{\ln (1+1/\lambda )/\ln \kappa },$$

$$\begin{array}{rl}{\lambda }_{k{k}^{\prime }}& =\frac{{\sum }_{\ell =1}^{\infty }{\sum }_{deg(i)=k,deg(j)=k^{\prime },d_{ij}=\ell }\bar{\lambda }(d_{ij},\lambda )}{NP(k)NP(k^{\prime })}\\ & =\frac{{\sum }_{\ell =1}^{\infty }k{\kappa }^{\ell -1}P(k^{\prime }|k,\ell )\bar{\lambda }(\ell ,\lambda )NP(k)}{NP(k)NP(k^{\prime })},\end{array}$$

$${\lambda }_{k{k}^{\prime }}=\frac{{\sum }_{\ell =1}^{\infty }\lambda (\mu \kappa {)}^{\ell -1}k{k}^{\prime }}{N⟨k⟩}=\lambda A_{k{k}^{\prime }}\sum _{\ell =1}^{\infty }(\mu \kappa {)}^{\ell -1}=\Pi \lambda A_{k{k}^{\prime }},$$

$$\frac{d{\rho }_k(t)}{dt}=-\bar{\delta }(k,\lambda ){\rho }_k(t)+(1-{\rho }_k(t))\Pi \lambda k\Theta ,$$

$${\rho }_k\simeq \Pi \lambda e^{a(\lambda )k}k\Theta ,$$

$$s(k)=1+a(\lambda )k,$$

$$\sum _{k^{\prime }}k{e}^{a(\lambda )k}\frac{k^{\prime }P(k^{\prime })}{⟨k⟩}{\rho }_{k^{{}^{\prime }}}=\frac{1}{\Pi \lambda }{\rho }_k.$$

$$\sum _{k=k_{min}}^{k_{max}}\frac{kP(k)}{⟨k⟩}\Pi \lambda k{e}^{a(\lambda )k}=1,$$

$${\Pi }_c{\lambda }_c{e}^{a({\lambda }_c)k_x}={\lambda }_c^{HMF},$$

$${\rho }_k\simeq {\lambda }_c^{HMF}k\Theta e^{(k-k_x)a({\lambda }_c)},$$

$$\Pi <\frac{{\lambda }_c^{HMF}}{{\lambda }_c}.$$

$$⟨k⟩=e^{a(\lambda )k_{max}+\ln k_{max}+\ln P(k_{max})+\ln N\ln \lambda }.$$

$$\frac{d{\rho }_i(t)}{dt}=-{\rho }_i(t)+\lambda (1-{\rho }_i(t))\sum _j{a}_{ij}{\rho }_j(t),$$

$$A_{k_i{k}_j}=\frac{P(k_j|k_i)k_i}{NP(k_j)}.$$

$$\frac{d{\rho }_k(t)}{dt}=-{\rho }_k(t)+\lambda (1-{\rho }_k(t))\sum _{k^{\prime }}{A}_{k{k}^{\prime }}{\rho }_{k^{\prime }}(t)NP(k^{\prime }).$$

$$\frac{d{\rho }_k(t)}{dt}=-{\rho }_k(t)+(1-{\rho }_k(t))\lambda k\Theta ,$$

$${\rho }_k\simeq \lambda k\Theta .$$

$$\bar{\lambda }(d_{ij},\lambda )\approx \lambda e^{-b(\lambda )(d_{ij}-1)}=\lambda {\mu }^{(d_{ij}-1)},$$

$$\begin{array}{rl}& \frac{d{\rho }_k(t)}{dt}=-\bar{\delta }(k,\lambda ){\rho }_k(t)+\\ & (1-{\rho }_k(t))\frac{{\sum }_{deg(i)=k}{\sum }_{j\ne i}\bar{\lambda }(d_{ij},\lambda ){\rho }_j(t)}{NP(k)}.\end{array}$$

$$\begin{array}{rl}& \frac{{\sum }_{deg(i)=k}{\sum }_{j\ne i}\bar{\lambda }(d_{ij},\lambda ){\rho }_j(t)}{NP(k)}\\ & =\frac{{\sum }_{k^{\prime }}{\sum }_{deg(i)=k,deg(j)=k^{\prime }}\bar{\lambda }(d_{ij},\lambda ){\rho }_j(t)}{NP(k)}\\ & =\frac{{\sum }_{k^{\prime }}{\rho }_{k^{\prime }}(t){\sum }_{deg(i)=k,deg(j)=k^{\prime }}\bar{\lambda }(d_{ij},\lambda )}{NP(k)}\\ & =\sum _{k^{\prime }}NP(k^{\prime }){\rho }_{k^{\prime }}(t)\frac{{\sum }_{deg(i)=k,deg(j)=k^{\prime }}\bar{\lambda }(d_{ij},\lambda )}{NP(k)NP(k^{\prime })}\\ & =\sum _{k^{\prime }}NP(k^{\prime }){\rho }_{k^{\prime }}(t){\lambda }_{k{k}^{\prime }},\end{array}$$

$${\lambda }_{k{k}^{\prime }}=\frac{{\sum }_{deg(i)=k,deg(j)=k^{\prime }}\bar{\lambda }(d_{ij},\lambda )}{NP(k)NP(k^{\prime })}$$

$$\begin{array}{rl}& \frac{d{\rho }_k(t)}{dt}=-\bar{\delta }(k,\lambda ){\rho }_k(t)+\\ & \lambda N{\left[\frac{k}{N⟨k⟩}\right]}^{\frac{b(\lambda )}{\ln \kappa }}\sum _{k^{\prime }}{k}^{\prime \frac{b(\lambda )}{\ln \kappa }}P(k^{\prime }){\rho }_{k^{\prime }}(t)[1-{\rho }_k(t)].\end{array}$$

$${\rho }_k\simeq e^{a(\lambda )k}\lambda N{\left[\frac{k}{N⟨k⟩}\right]}^{\frac{b(\lambda )}{\ln \kappa }}\sum _{k^{\prime }}{k}^{\prime \frac{b(\lambda )}{\ln \kappa }}P(k^{\prime }){\rho }_{k^{\prime }}$$

$${\rho }_k\propto k^{\beta }{e}^{a(\lambda )k},$$

$$\frac{d{\rho }_i(t)}{dt}=-\bar{\delta }(k_i,\lambda ){\rho }_i(t)+(1-{\rho }_i(t))\sum _{\ell }\sum _j\bar{\lambda }(\ell ,\lambda )b_{ij}(\ell ){\rho }_j,$$

$$\sum _{\ell }\sum _{k^{\prime }}\bar{\lambda }(\ell ,\lambda )B_{k{k}^{\prime }}(\ell ){\rho }_{k^{\prime }}(t)NP(k^{\prime }),$$

$$\begin{array}{rl}& \sum _{\ell }\sum _{k^{\prime }}\bar{\lambda }(\ell ,\lambda )B_{k{k}^{\prime }}(\ell ){\rho }_{k^{\prime }}(t)NP(k^{\prime })\\ & =\sum _{\ell }\bar{\lambda }(\ell ,\lambda )k{\kappa }^{\ell -1}\sum _{k^{\prime }}P(k^{\prime }|k,\ell ){\rho }_{k^{\prime }}(t)\\ & =\lambda k\Theta \sum _{\ell }(\mu \kappa {)}^{\ell -1}\\ & =\Pi \lambda k\Theta ,\end{array}$$

$$\Pi =\sum _{\ell }(\mu \kappa {)}^{\ell -1}=\frac{1}{1-\mu \kappa }.$$

Full text

I Introduction

Networks network ; cn ; netS essentially capture the structure of individual interactions, contact and mobility patterns, through which information, innovation, fads, epidemics, and human behaviors spread across us mobp ; behavior ; metap ; hidden ; universal ; voter . Burst of investigations have been devoted to SIS model on networks to understand, predict and control diverse recurrent contagion phenomena, such as influenza-like diseases, computer virus, memes etc. model ; process ; unifi . There is a threshold $\lambda_{c}$ for the spreading rate corresponding to the nonequilibrium absorbing state phase transition npt ; npt1 . The contagion will outbreak above $\lambda_{c}$, and die out below it. The heterogeneity of real-world networks characterized by a power law degree distribution $P(k)\sim k^{-\gamma}$ has a nontrivial impact on the threshold. As a seminal theory to understand SIS model, the heterogeneous mean field (HMF) theory predicts that for uncorrelated random networks $\lambda_{c}^{HMF}=\langle k\rangle/\langle k^{2}\rangle$, where $\langle k\rangle$ and $\langle k^{2}\rangle$ are the first and second moments of $P(k)$ hmfP . For $2<\gamma<3$, $\langle k^{2}\rangle$ is divergent in the thermodynamic limit, thus, $\lambda_{c}$ will converge to zero. While for $\gamma>3$, $\lambda_{c}>0$, since $\langle k^{2}\rangle$ is finite.

In spite of the rather simplicity of SIS model, considerable controversies and confusions arise from a number of competing theories. The focus is that $\lambda_{c}>0$ for $\gamma>3$ seems to be a well established fact, however, different voices never stop. The following is a rapid review. The quenched mean field (QMF) theory replaces the annealed adjacency matrix in HMF by the quenched one, predicting $\lambda_{c}^{QMF}=1/\Lambda_{1}$, where $\Lambda_{1}\approx\max\{\sqrt{k_{max}},\langle k^{2}\rangle/\langle k\rangle\}$ is the largest eigenvalue of the adjacency matrix eigenv ; eigenv1 . QMF readily gives the same threshold as HMF for $2<\gamma<2.5$. While for $\gamma>2.5$, $\Lambda_{1}$ shifts from $\langle k^{2}\rangle/\langle k\rangle$ to $\sqrt{k_{max}}$. Then $\lambda_{c}^{QMF}\simeq 1/\sqrt{k_{max}}$, which is zero for any networks with divergent maximum degree qmf0 , just at odds with HMF theory. The node with maximum degree which independently determines $\Lambda_{1}$ was deemed the self-sustainable source. Goltsev et al. applied the spectral approach to QMF, only to find that it causes the eigenvector localization (EL) local , for which most of the weight of eigenvector concentrates around a small number of nodes with highest degree el1 ; el2 ; el3 . It immediately means that there are only a finite number of active nodes localized around hubs. Therefore, $\lambda_{c}^{QMF}\simeq 1/\sqrt{k_{max}}$ is not a genuine threshold. The true threshold appears at the inverse of upper delocalized eigenvalue, which is slightly smaller than $\lambda_{c}^{HMF}=\langle k\rangle/\langle k^{2}\rangle$. Once again, it is a positive value for $\gamma>3$.

Lee et al. sparked a big stir lee . It was showed that when $\lambda_{c}^{QMF}<\lambda<\lambda_{c}$, the activity localized around hubs in Ref. local in fact belongs to the Griffiths phase grif , featuring slow relaxation dynamics. Owing to the irreversible dynamical fluctuation, which is seldom considered in most mean-field methods, these active domains eventually fall into absorbing state with exceedingly long relaxation time. Moreover, they conjectured that the Griffiths to active phase transition is triggered by the percolation of the active domains through direct connection of hubs, which enables the mutual reinfection of hubs. Depending on the corresponding percolation threshold, the epidemic threshold $\lambda_{c}$ is also positive for $\gamma>3$. However, this compelling mechanism is criticized by Boguñá et al.nature . They invoked Durrett’s et al. work durret1 ; durret2 , pointing out that the direct connection of hubs is not a necessary condition for reinfection. The null threshold is retrieved when considering the possibility of reinfection between any pair of nodes.

Here we argue that these studies involve the key ingredients toward resolving the controversies, including EL, dynamical fluctuation and Griffiths effect, the reinfection between distant nodes. In particular, EL interpolates the Griffiths phase between absorbing and active phase gp , which may dramatically change our understanding of SIS model. The dynamical evolution in this situation cannot be described by those conventional theories neglecting these points such as HMF or QMF. There are other smart studies without considering above ingredients but pairwise dynamical correlations, such as the dynamical message passing approach dmp , effective degree edegree and effective branching factor methods ebf . These studies give higher threshold than $\lambda_{c}^{HMF}$, which are not in agreement with simulations on networks with EL sample ; qs ; nature .

Ref. nature seems to put an end to the ongoing debate on the vanishing threshold. Due to the rather rough annealed mean-field approximation, however, it provides few results other than the upper bound of threshold, and their theory lacks thorough and convincing examination. Moreover, the role of EL and Griffiths effect have not been explicitly incorporated in their theory. And it remains to be uncovered how the system avoids the irreversible dynamical fluctuation, triggering Griffiths to active phase transition. A consistent and comprehensive investigation of above theoretical ingredients is absent. We develop a new mean-field equation based on Ref. nature for this attempt. It gives us the leverage to examine the theory presented in Ref. nature in several aspects. It also allows us to study the impact of EL in depth beyond the framework of QMF. Our results suggest that SIS model on networks with EL is markedly different from the delocalized case.

II Degree-based mean-field approximation

Let us begin by recalling the theory of Ref. nature . The core is to take into account that an infected node $i$ may infect any other node $j$ in the network separated by distance $d_{ij}$ at a long enough time scale at the rate $\bar{\lambda}(d_{ij},\lambda)$. “On long time scales node $i$ is considered as susceptible only when the node and all of its nearest neighbors in the original graph are susceptible” nature . Therefore, the corresponding recovery rate is replaced by the inverse of the characteristic lifespan of the star $\bar{\delta}(k_{i},\lambda)=\tau^{-1}(k_{i},\lambda)\approx e^{-a(\lambda)k_{i}}$. The evolution of SIS model is depicted by

[TABLE]

where $\bar{\lambda}(d_{ij}\lambda)\approx\lambda\mu^{(d_{ij}-1)}$, and $\mu=\lambda/(1+\lambda)$ is the probability that an infected node infects its nearest neighbor before returning to susceptible state. Compared with QMF working on the quenched adjacency matrix, Eq. (1) is defined on a fully connected graph, where $k_{i}$ and $d_{ij}$ reflect the structure of the original network. See Appendix A for detailed introduction about Eq. (1).

In the following, we develop the degree-based mean-field equation based on Eq. (1). HMF equation is just a result of the degree-based mean-field approximation of QMF equation. Nonetheless, similar mean-field calculation regarding Eq. (1) seems sophisticated. Incorporating the definition $\rho_{k}(t)\equiv\sum_{deg(i)=k}\rho_{i}(t)/(NP(k))$ with Eq. (1), we obtain

[TABLE]

where $\lambda_{kk^{\prime}}\equiv\sum_{deg(i)=k,deg(j)=k^{\prime}}\bar{\lambda}(d_{ij},\lambda)/NP(k)NP(k^{\prime})$, see Eq. (21)-Eq. (23) in Appendix A for the derivations.

A very important mean-field approximation was made in Ref. nature , which merely considers the dynamical influence of nodes at the average distance while omits others:

[TABLE]

where $d_{kk^{\prime}}=1+\ln(N\langle k\rangle/kk^{\prime})/\ln\kappa$ is the average distance between nodes of degree $k$ and $k^{\prime}$ for random networks, and $\kappa=\langle k^{2}\rangle/\langle k\rangle-1$ is the average branching factor, $A_{kk^{\prime}}=kk^{\prime}/(N\langle k\rangle)$ is the annealed adjacency matrix for random networks.

The upper bound of threshold was estimated by plugging Eq. (3) into Eq. (2) in Ref. nature . Though, the rather rough mean-field approximation made in Eq. (3) impedes further convincing examination of Eq. (1). In particular, the reinfection of a node based on a fully connected graph characterized by the sum $\sum_{j\neq i}\bar{\lambda}(d_{ij},\lambda)\rho_{j}(t)$ lacks faithful validation. Here we circumvent above problem by making a complete consideration of the reinfections of nodes at each distance

[TABLE]

where $P(k^{\prime}|k,\ell)$ is the probability of finding a neighboring node with degree $k^{\prime}$ at distance $\ell$ from a given node of degree $k$. For sparse uncorrelated random networks, $P(k^{\prime}|k,\ell)=k^{\prime}P(k^{\prime})/\langle k^{\prime}\rangle$. Insert it into Eq. (4), we have

[TABLE]

where $\Pi=\sum_{\ell=1}^{\infty}(\mu\kappa)^{\ell-1}$ is a parameter implying the reinfection patterns. Numerical simulations show that $\lambda_{c}<\lambda_{c}^{HMF}$ for $\gamma>2.5$ (except for networks with very small size) nature ; qs , hence, $\mu\kappa<1$. Then we have $\Pi=1/(1-\mu\kappa)$. And we see the ratio of Eq. (5) to Eq. (3), $\Pi A_{kk^{\prime}}^{1-\ln(1+1/\lambda)/\ln\kappa}\gg 1$.

Instead of Eq. (3), substitute Eq. (5) into Eq. (2), we obtain

[TABLE]

where $\Theta=\sum_{k^{\prime}}\rho_{k^{\prime}}(t)k^{\prime}P(k^{\prime})/\langle k^{\prime}\rangle$ is the probability of a randomly selected link connected to an infected node. It demonstrates a markedly different physical picture compared with HMF or QMF, which can be validated by quasi-stationary simulation limit ; sample (see Appendix B for technical details).

It is worth noting that each term in $\Pi\lambda k\Theta=\sum_{\ell=1}^{\infty}\lambda k\Theta(\mu\kappa)^{\ell-1}$ denoted as $\psi_{k}(\ell)=\lambda k\Theta(\mu\kappa)^{\ell-1}$ describes the reinfections between a node of degree $k$ and those nodes separated by distance $\ell$. Let us make it clear. Obviously, the sum in Eq. (1) $\sum_{j\neq i}\bar{\lambda}(d_{ij},\lambda)\rho_{j}(t)=\sum_{\ell}\sum_{j}\bar{\lambda}(\ell,\lambda)b_{ij}(\ell)\rho_{j}$, where $b_{ij}(\ell)=\delta_{d_{ij},\ell}$, and $\delta_{d_{ij},\ell}$ is Kronecker function. The degree-based mean-field approximation assumes that all the nodes with degree $k$ are statistically equivalent process . The state of each node of degree $k$ is specified by $\rho_{k}(t)$, and the quenched network structure is replaced by the average of network ensembles. Thus, the degree-based mean-field approximation for $\sum_{j\neq i}\bar{\lambda}(d_{ij},\lambda)\rho_{j}(t)$ is $\sum_{\ell}\sum_{k^{\prime}}\bar{\lambda}(\ell,\lambda)B_{kk^{\prime}}(\ell)\rho_{k^{\prime}}(t)NP(k^{\prime})$. Here $B_{kk^{\prime}}(\ell)=P(k^{\prime}|k,\ell)k\kappa^{\ell-1}/(NP(k^{\prime}))$ is the annealed counterpart of $b_{ij}(\ell)$. As a result, $\psi_{k}(\ell)=\bar{\lambda}(\ell,\lambda)k\kappa^{\ell-1}\sum_{k^{\prime}}P(k^{\prime}|k,\ell)\rho_{k^{\prime}}(t)=\lambda k\Theta(\mu\kappa)^{\ell-1}$, which decays exponentially with distance. Detailed derivations for above results are showed by Eq. (27)-Eq. (30) in Appendix A. As the impact of reinfection at each distance is explicitly taken into account, Eq. (6) enables us to faithfully verify Eq. (1) from a broad perspective.

III The exotic critical behavior of $\rho_{k}$ under the influence of eigenvector localization

First, we study the critical behavior of $\rho_{k}$ around the critical point. It is not difficult to derive from the HMF that $\rho_{k}\propto k$. It seems to be a widely accepted claim. For example, it was used to estimate the average infection density $\rho$ lee . Indeed, it is true for $2<\gamma<2.5$. A justified example is shown in Fig. 1 (a). For $\gamma>2.5$, where EL arises, however, it is striking that it no longer holds. As shown in Fig. 1 (b)-(d), with increasing of $k$, the deviations from the benchmark get more and more remarkable. Behind the significant discrepancy lies the manifest failure of the HMF equation. It calls for new theory. If Eq. (3) is plugged into Eq. (2), we obtain $\rho_{k}\propto k^{\ln(1+1/\lambda)/\ln\kappa}e^{a(\lambda)k}$ (see Eq. (24)-Eq. (26) in Appendix A). The derivative of this prediction under double log plot $s(k)\simeq\ln(1+1/\lambda)/\ln\kappa+a(\lambda)k\gg 1$ can explain the huge deviation for the highly connected nodes, yet fails to predict that $s(k)\simeq 1$ for the small-degree nodes shown in Fig. 1.

In contrast, the linear approximation of Eq. (6) gives

[TABLE]

which is perfectly in agreement with simulation, see Fig. 1 (b)-(d). The derivative of $\rho_{k}$,

[TABLE]

is a coherent explanation for the deviation from benchmark in Fig. (1). Assume that $a(\lambda)\approx a\lambda^{2}$ and $a\sim\mathcal{O}(1)$ nature , the deviation $\Delta_{s}=a(\lambda)k$ is pretty tiny for small-degree nodes, as it only slightly deviates from $s(k)=1$. Though, with increasing of $k$, $\Delta_{s}$ becomes larger and larger. For highly connected nodes, e.g., $\lambda^{2}k\sim 1$, the deviation is so huge that by no means can be neglected. We argue that the exponential term $e^{a(\lambda)k}$ in Eq. (7) responsible for the increasing deviation originates from EL.

For random networks, $\Lambda_{1}$ is largely determined by two special subgraphs in isolation: the ultramost dense $K$-core core and the star subgraph composed of the hub with degree $k_{max}$ and its immediate neighbors. In specific, $\Lambda_{1}\approx\max\{\Lambda_{1}^{core},\Lambda_{1}^{star}\}$, where $\Lambda_{1}^{core}$ and $\Lambda_{1}^{star}$ are the largest eigenvalues for $K$-core and the star respectively eigenv1 . This insight suggests that the dynamics are determined by either $K$-core or the star. When $2<\gamma<2.5$, $\Lambda_{1}\approx\Lambda_{1}^{core}$. In this case, the component of the leading eigenvector is proportional to the degree of the corresponding node. The dense core, which has a large proportion of the weight of eigenvector el2 ; el3 , plays the role of sustainable spreading source. And the spreading can be seen as a branching process.

When $\gamma>2.5$, it turns out that $\Lambda_{1}\approx\Lambda_{1}^{star}$. More generally, the $n$th largest eigenvalue of the quenched adjacent matrix $\Lambda_{n}$ is almost solely determined by the star centered in the node with $n$th largest degree $k_{n}$. For sufficiently large $k_{n}$, $\Lambda_{n}\approx\sqrt{k_{n}}$ eigenv . And the corresponding eigenvector is localized around the center of the star nlocal . Therefore, each star centered in the hub with degree, e.g., $k\gtrsim 1/\lambda^{2}$, is approximately dynamically quasi-independent with its own activation threshold $\lambda_{n}\approx 1/\sqrt{k_{n}}$ and lifespan $\tau(k,\lambda)\approx e^{a(\lambda)k}$ just like an isolated graph. Hence, the assumption of $\bar{\delta}(k_{i},\lambda)$ is conditional on EL. As shown in Fig. 1, Eq. (7) does not hold for networks with delocalized eigenvector when $\gamma<2.5$. As a consequence, Eq. (1) is invalid when $\gamma<2.5$.

Eq. (7) can be detailed as

[TABLE]

Based on this result, we define the spreading matrix S with element $S_{kk^{\prime}}=ke^{a(\lambda)k}k^{\prime}P(k^{\prime})/\langle k\rangle$. Thereby Eq. (9) is equivalent to $\textbf{S}\hat{\rho}=\Lambda\hat{\rho}$. According to the Perron-Frobenius theorem, if $\Lambda$ is equal to the spectral radius, then $\hat{\rho}$ is a positive eigenvector, and vice versa. Hence, $\Lambda=1/\Pi\lambda$ is a function of the genuine threshold, marking the onset of endemic state, where a finite fraction of nodes are infected. Substitute Eq. (7) into Eq. (9), we obtain a definite critical condition for Griffiths to active phase transition.

[TABLE]

where $k_{min}$ denotes the minimum degree.

It is interesting that Eq. (10) means that $\textmd{Tr}(\textbf{S})=\Lambda_{1}$, where $\textmd{Tr}(\cdot)$ is the trace operator. Similar result also holds for the annealed adjacency matrix $A_{kk^{\prime}}$, that is, $\textmd{Tr}(\textbf{A})=\langle k^{2}\rangle/\langle k\rangle$. It is not surprising. Owing to $[\textmd{Tr}(\textbf{S})]^{2}=\textmd{Tr}(\textbf{S}^{2})$ and $[\textmd{Tr}(\textbf{A})]^{2}=\textmd{Tr}(\textbf{A}^{2})$, there is only one non-zero eigenvalue for the matrix.

Since $e^{a(\lambda)k_{min}}\sum_{k}k^{2}P(k)<\sum_{k}k^{2}P(k)e^{a(\lambda)k}<e^{a(\lambda)k_{max}}\sum_{k}k^{2}P(k)$, there must exist $k_{x}\in(k_{min},k_{max})$ making $\sum_{k}k^{2}P(k)e^{a(\lambda)k}=\langle k^{2}\rangle e^{a(\lambda)k_{x}}$. With this result, Eq. (10) is equivalent to,

[TABLE]

where $\Pi_{c}$ is the corresponding value of $\Pi$ at the critical point. Combine Eq. (11) with Eq. (7), in the vicinity of critical point, the later can be written as

[TABLE]

which has appealing implications. EL exponentially reinforces the infection of large-degree nodes ($k>k_{x}$), while inversely suppresses the infection of small-degree nodes ($k<k_{x}$), see Fig. 1.

IV Examining the value of $\Pi$ and estimation of the threshold for Griffiths to active phase transition

There is a big drawback for the rough mean-field approximation made in Eq. (3). It is unable to confirm the validity of the sum in Eq. (1) characterizing the reinfection patterns defined on a fully connected graph. Eq. (6) suggests that $\Pi$ can serve as an effective indicator. It can be best verified by checking whether $\Pi=\sum_{\ell=1}^{\infty}(\mu\kappa)^{\ell-1}=1/(1-\mu\kappa)$. We provide three approaches to examine this theoretical prediction.

On one hand, there are two parameters $\Pi$ and $a$, the values of which remain to be determined. On the other hand, we have two equations (Eq. (7) and Eq. (10)) involving the two parameters. This mathematically allows us to solve the two parameters with assistance of the simulation data of $\lambda_{c}$ and $\rho_{k}$, see Appendix C for how this is done. As shown in Fig. 2, the corresponding numerical solutions suggest that $\Pi<1$, which is significantly smaller than $\Pi=\sum_{\ell=1}^{\infty}(\mu\kappa)^{\ell-1}>1$.

Besides, Eq. (11) provides another two methods to verify the theoretical value of $\Pi$. Plug $\Pi=1/(1-\mu\kappa)$ and the simulated value of $\lambda_{c}$ into Eq. (11), $e^{-a(\lambda_{c})k_{x}}=\lambda_{c}\Pi/\lambda_{c}^{HMF}$. It unexpectedly produces an invalid result $e^{-a(\lambda_{c})k_{x}}>1$, see Fig. 2 for demonstration. In contrast, $e^{-a(\lambda_{c})k_{x}}<1$ always hold if $\Pi<1$ and $\lambda_{c}<\lambda_{c}^{HMF}$ .

Utilizing Eq. (11) and the strong constraint $\lambda_{c}<\lambda_{c}^{HMF}$, we obtain the upper bound of $\Pi$. Since $e^{-a(\lambda)k_{x}}<1$, we have

[TABLE]

Since $\lambda_{c}$ can be faithfully obtained from simulation, this upper bound is useful for finite size systems. As displayed in Fig. 2, the theoretical prediction $\Pi=1/(1-\mu\kappa)$ is substantially bigger than this upper bound, which thus violates the inequality of Eq. (13). However, the simulation value $\Pi_{s}<1$ well obeys this inequality .

Now we are prepared to analyze the threshold $\lambda_{c}$. Similar to Ref. nature , the upper bound of $\lambda_{c}$ can be estimated by merely considering the reinfections from those nodes at the average distance, since $\Pi\lambda\gg\bar{\lambda}(\langle\ell\rangle,\lambda)$, where $\langle\ell\rangle\sim\ln N$ is the average distance between any two nodes in a small-world network. The upper bound of $\lambda_{c}$ is obtained by substituting $\bar{\lambda}(\langle\ell\rangle,\lambda)$ for $\Pi\lambda$ in Eq. (10). Due to $\sum_{k}k^{2}P(k)e^{a(\lambda)k}/\langle k\rangle>\sum_{k}kP(k)e^{a(\lambda)k}/\langle k\rangle=\langle\tau\rangle$, the upper bound can be estimated by $\bar{\lambda}(\langle\ell\rangle,\lambda)\langle\tau\rangle=1$, or $1/\bar{\lambda}(\langle\ell\rangle,\lambda)=\langle\tau\rangle$. Here, $1/\bar{\lambda}(\langle\ell\rangle,\lambda)$ is just the average time it takes for an infected node to infect for the first time a node at the average distance Hence, it also provides a condition that on average a star with an active center is capable of propagating the infection to nodes as far as at the average distance before spontaneously recovering. It requires $\langle\tau\rangle\rightarrow\infty$ in the thermodynamic limit. Without compromising the upper bound, $\bar{\lambda}(\langle\ell\rangle,\lambda)\langle\tau\rangle=1$ can be reduced to

[TABLE]

If $k_{max}$ grows faster than $\ln N$, and $P(k)$ decays slower than exponential, with increasing of $N$ while keeping $\lambda$ fixed, then the exponential term will become bigger than $\langle k\rangle$. As a result, $\lambda$ must converge to zero in the thermodynamic limit. Assume that $a(\lambda)\propto\lambda^{2}$, we conclude that $\lambda_{c}\rightarrow 0$ at a speed much slower than $\lambda_{c}^{QMF}\simeq 1/\sqrt{k_{max}}$ for $\gamma>2.5$.

Quasi-stationary simulations have already shown that $\lambda_{c}^{QMF}<\lambda_{c}<\lambda_{c}^{HMF}$ for $\gamma>2.5$ sample ; qs ; nature . Here we provide closer bounds for $\lambda_{c}$ by substituting the corresponding bound of $\Pi$ into Eq. (10). The vanishing threshold means $\mu\kappa\rightarrow 0$ for $\gamma>3$, while the simulated value $\Pi_{s}\rightarrow 1$. As a result, for large-size networks, $\mu\kappa<\Pi_{s}<1$ for $\gamma>3$. For $2.5<\gamma<3$, Fig. 2 shows that it is true as well. Therefore, a closer upper bound for $\lambda_{c}$ can be obtained if $\Pi=\mu\kappa$, and $\Pi=1$ accordingly leads to a closer lower bound. In contrast, inserting $\Pi=1/(1-\mu\kappa)$ into Eq. (10) yields a threshold much smaller than simulation, see Fig. 3.

V Discussion and conclusion

Our results allow us to evaluate to what extent Ref. nature is valid. First, the introduction of $\delta(k_{i},\lambda)$ is closely associated with EL. It is invalid when $\gamma<2.5$. However, Ref. nature claimed that Eq. (1) can be applied to any network. The second concern is the remarkable discrepancy between the value of $\Pi$ obtained from simulation and the mean-field expectation $\Pi=1/(1-\mu\kappa)$, which stems from the sum in Eq. (1) defined on a fully connected graph portraying the pattern of mutual reinfections between distant nodes accounting for the null threshold. This remarkable contradiction underscores that Eq. (1) is inadequate to characterize the dynamical process. Since it is difficult to reconcile this confusing contradiction, We leave it as an open question.

Although Eq. (6) does not have complete theoretical power as the correct expression for $\Pi$ remains unclear, it has already produced a series of novel results, answering some important questions of broad interest. Eq. (10) portrays the critical condition for Griffiths to active phase transition: on average, a node getting infected can further infect another node within the lifespan of the star centered in itself. Thereby, the star is able to avoid getting trapped in the absorbing state, marking the onset of the phase transition. As a hallmark of the transition, we show how the critical behavior of $\rho_{k}$ is dramatically changed by EL. Instead of the conventional belief $\rho_{k}\propto k$, it is striking that $\rho_{k}\propto ke^{a(\lambda)k}$ for $\gamma>2.5$, where the exponential term originates from EL. As a consequence, HMF and QMF are not valid when $\gamma>2.5$.

On one hand, Eq. (14) reproduces the conclusion of Ref. nature regarding the vanishing threshold. On the other hand, it suggests that whether $\lambda_{c}\rightarrow 0$ depends on the competition between the damping long-range reinfection rate $\bar{\lambda}(\langle\ell\rangle,\lambda)$ and the growing average lifespan $\langle\tau\rangle$. The divergence of $\langle\tau\rangle$ is a crucial condition for the effective mutual reinfection between remote nodes and the null threshold in the thermodynamic limit. Eq. (11) makes a direct connection between $\lambda_{c}$ and $\lambda_{c}^{HMF}$. It could be seen as a quantitative criterion explicitly separating highly connected nodes ($k>k_{x}$) from small-degree nodes ($k<k_{x}$). Together with Eq. (12), we conclude that owing to EL, the infection of highly connected nodes is exponentially reinforced, while the infection of small-degree nodes is inversely suppressed.

Acknowledgement. We thank S. C. Ferreira for the discussion of quasistationary simulation. We are also grateful to Jia-Rong Xie for discussing the manuscript. This work is funded by the National Natural Science Foundation of China (Grant No. 71874172).

Appendix A Some derivation details for the degree-based mean-field equation and threshold

If the dynamical correlations are neglected, the quenched mean field (QMF) equation takes the form:

[TABLE]

where $\rho_{i}(t)$ is the infection density of node $i$, and $a_{ij}$ is the element of adjacency matrix. Next, we develop the degree-based mean-field equation based on QMF. The degree-based mean-field approximation basically assumes that all the nodes with the same degree are statistically equivalent. This assumption has two implications. All the nodes of degree $k$ share the same sate $\rho_{k}(t)$. And the quenched adjacency matrix $a_{ij}$ is replaced by the average of network ensembles.

[TABLE]

Applying above two principles to Eq. (15), the resulting degree-based mean-field equation is

[TABLE]

Plug $A_{kk^{\prime}}=kk^{\prime}/(N\langle k\rangle)$ for uncorrelated random networks into Eq. (17), we obtain the familiar form of HMF equation:

[TABLE]

where $\Theta=\sum_{k^{\prime}}\rho_{k^{\prime}}(t)k^{\prime}P(k^{\prime})/\langle k^{\prime}\rangle$ refers to the probability of a randomly selected link pointing to an infected node. The linear approximation of Eq. (18) around the critical point leads to

[TABLE]

The QMF equation only involves the influence of the nearest neighbors. In contrast, Boguñá et al. took a leap by considering that at a long enough time scale, a node suffers the infections of all the nodes in the network nature . The spreading rate $\lambda$ is replaced by the effective rate

[TABLE]

where $b(\lambda)=\ln(1+1/\lambda)$, $\mu=\lambda/(1+\lambda)$. It denotes the rate at which an infected node $i$ can infect a node $j$ separated by distance $d_{ij}$. “On long time scales node $i$ is considered as susceptible only when the node and all of its nearest neighbors in the original graph are susceptible” nature , the recovery rate of node $i$ is then replaced by the inverse of the characteristic lifetime of star subgraph composed of the node and its nearest neighbors: $\bar{\delta}(k_{i},\lambda)=\tau^{-1}(k_{i},\lambda)\approx e^{-a(\lambda)k_{i}}$, where $a(\lambda)\propto\lambda^{2}$. Based on these points, they suggested a time coarse-grained equation, i.e., Eq. (1). Combining $\rho_{k}(t)\equiv\sum_{deg(i)=k}\rho_{i}(t)/(NP(k))$ with Eq. (1), we obtain

[TABLE]

The sum in this equation can be further calculated as follows:

[TABLE]

where

[TABLE]

The mean-field calculation of $\lambda_{kk^{\prime}}$ becomes crucial. Ref.nature coped with this problem by making an important mean field approximation which merely keeps the contribution of nodes at the average distance while omits others: $\lambda_{kk^{\prime}}\approx\bar{\lambda}(d_{kk^{\prime}},\lambda)$ where $d_{kk^{\prime}}=1+\ln(N\langle k\rangle/kk^{\prime})/\ln\kappa$ is the average distance between nodes of degree $k$ and $k^{\prime}$ for small-world random networks. It leads to a rather sophisticated equation reported in Ref.nature :

[TABLE]

In the vicinity of critical point, Eq. (24) yields a prediction significantly different from HMF theory (see Eq. (19)):

[TABLE]

In short,

[TABLE]

where $\beta=b(\lambda)/\ln\kappa$, and $\kappa=\langle k^{2}\rangle/\langle k\rangle-1$ is the average branching factor.

We propose a different degree-based mean-field approximation of Eq. (1), which can be rewritten as

[TABLE]

where $b_{ij}(\ell)=\delta_{d_{ij},\ell}$. According to the principles of degree-based mean-field approximation, the sum in above equation is replaced by

[TABLE]

where $B_{kk^{\prime}}(\ell)=P(k^{\prime}|k,\ell)k\kappa^{\ell-1}/(NP(k^{\prime}))$ is the annealed counterpart of $b_{ij}(\ell)$, and $P(k^{\prime}|k,\ell)$ is the probability of finding a neighboring node with degree $k^{\prime}$ at distance $\ell$ from a given node of degree $k$. For sparse uncorrelated random networks, $P(k^{\prime}|k,\ell)=k^{\prime}P(k^{\prime})/\langle k^{\prime}\rangle$. As a result,

[TABLE]

where

[TABLE]

Since the influence of nodes at each distance is completely taken into account, we derive the least approximated degree-based mean-field equation for Eq. (1), i.e., Eq. (6). It enables us to understand Griffiths to active phase transition.

Appendix B Quasistationary simulation

In this paper, we use the configuration model to generate uncorrelated random networks for quasistationary simulation ucm . The degree sequence $\{k_{i}\}$ is generated as follows. First, the degree distribution $P(k)=ck^{-\gamma}$ is normalized as $\sum_{k_{min}}^{k_{c}}P(k)dk=1$, where $k_{c}\leq\sqrt{N}$ is the structural degree cutoff that can avoid inherent degree correlation ucm , and $k_{min}$ refers to the minimum degree. Second, the number of nodes with degree $k$ is $n(k)=NP(k)$.

The principle of quasistationary simulation is to calculate the statistical average over possible evolutionary configurations by certain sampling strategy without being trapped in the absorbing state. The discrete-time approach, which synchronously updates the states of all the nodes with a fixed time interval such as $\Delta t=1$, can lead to incorrect result, especially when the recovery rate is unit limit . Instead, the continuous-time approach is employed in this paper.

It is implemented according to the following algorithm. At each time step, the number of infected nodes $N_{i}$ and the sum of the degree of these nodes $E_{i}$ are computed respectively. With probability $N_{i}/(N_{i}+\lambda E_{i})$, a randomly picked infected node returns to susceptible state. While with probability $\lambda E_{i}/(N_{i}+\lambda E_{i})$, an edge connected to an infected node is chosen at random. If another endpoint of the edge is a susceptible node, the node is set as infected. Otherwise, the simulation just continues. After this operation, recompute the values of $N_{i}$ and $E_{i}$, and the total time increases by the amount $\Delta t=1/(N_{i}+\lambda E_{i})$. Above procedures are iterated. The threshold is determined by the position of the peak in the susceptibility curve. And the susceptibility is defined as

[TABLE]

Once the system visits absorbing state, how to select the seeds to restart the evolution is crucial. It determines whether the threshold can be accurately captured. The standard updating scheme is to randomly select a historical evolution configuration from a costly maintained memory. Such scheme is not a good one, since it does not produce a susceptibility curve with sharp peak, but multiple ambiguous peaks. In this case, it is difficult to confirm the accurate position of the transition point. That is why Ref. error provided misleading results. Another scheme which fixes the hub node as restarting seed also suffers such a drawback sample . Instead, we employ the so called reflection boundary condition, which simply selects the last active node as the restarting seed sample . This method is able to overcome above problem. See Fig. 4 for examples.

Appendix C How to determine the values of $\Pi$ and $a$

The sum over infections of all the other nodes in Eq. (1) is defined on a fully connected graph, which highlights its difference compared with QMF theory. The validity of this theoretical proposition remains to be confirmed. It can be effectively addressed by examining whether $\Pi=\sum_{\ell=1}^{\infty}(\mu\kappa)^{\ell-1}=1/(1-\mu\kappa)$. The quasi-empirical parameter $a(\lambda)$ in Eq. (7) makes it a nontrivial job. Though, the combination of Eq. (7) and Eq. (10) involving the two parameters mathematically allows us to determine their values using the simulation data of $\lambda_{c}$ and $\rho_{k}$. By fitting the data of $\rho_{k}$ with Eq. (7), we can roughly estimate their values. Meanwhile, only the values of $\Pi$ and $a$ satisfying Eq. (10) can be accepted. It is not difficult to find the valid values by iterative computations. In the following, let us clarify how this works.

Eq. (7) suggests that the derivative of $\rho_{k}$ under double log plot is $s(k)=1+a(\lambda)k$. Assume that $a(\lambda)\approx a\lambda^{2}$ nature , and $a\sim\mathcal{O}(1)$,

[TABLE]

which is well in agreement with simulated results. Note that if $a\gg 1$, $s(k)$ seriously deviates from the benchmark $s(k)=1$ for the small-degree nodes. And if $a\rightarrow 0$, the deviation from the benchmark for highly connected nodes will be so tiny that can be neglected . Both results do not coincide with simulation. Therefore, it is reasonable to assume that $a\sim\mathcal{O}(1)$. See the main text for detailed illustration. For small degree nodes, since $a\lambda^{2}k\ll 1$, Eq. (7) can be approximated as

[TABLE]

In this case, the value of $a$ has little influence to estimate $\Pi$. Thus, it is plausible to just try a specific value from $0<a\leq 1$, and $\Pi$ is estimated by fitting the data of $\rho_{k}$ constrained by $k<k_{u}$ with Eq. (7). Figure. 5 shows that both the values of $a$ and $k_{u}$ do have little influence to estimate $\Pi$. With increasing of $a$, as shown in Fig. 5 (a), $\Pi$ decreases slowly in a linear way.

The values of $\Pi$ and $a$ can be further determined in a consistent way. Fit the data of $\rho_{k}$ with a trial value of $a$ denoted as $a_{t}$ to yield the value of $\Pi$. Eq. (10) puts a strong constraint on $\Pi$ and $a$. Hence, substitute this value of $\Pi$ into Eq. (10), and obtain $a_{p}$ referring to the predicted value of $a$. Once $a_{p}$ is worked out, give this value to $a_{t}$. Iterate above computation procedure for several times. It quickly converges to $a_{p}=a_{t}$, see Fig. 6 (a). And $\Pi$ is accordingly determined as well. Eq. (7) is then specified by substituting the values of $a$ and $\Pi$. As shown in Fig. 6 (b), the specified expression of Eq. (7) fits fairly well with the simulated data. Finally, the fidelity of above computation is enough to come to the conclusion $\Pi<1$.