Multistage onset of epidemics in heterogeneous networks

These scale-free static networks in this paper can be generated using an algorithm proposed by Catanzaro et al. Catanzaro et al. (2005), which is called the uncorrelated configuration model. For each node we first assign it a degree $k_{i}$ according to the prescribed degree distribution, whose value is restricted to the range of [$k_{\min}$, $\sqrt{N}$], where $k_{\min}$ is the minimum degree of the node in the network and $N$ is the network size. And then we create a set of $k_{i}$ stubs that represent each of these edges with only a single tail connected to a node. If there is an uneven number of stubs, a random individual is given an extra stub. Finally, we construct the network by connecting pairs of these stubs chosen uniformly at random to make complete edges respecting the preassigned degrees. Self connections or duplicate edges between nodes are not allowed in the generation process. We note that the are several subtleties when implementing the configuration model by random stub-matching Fosdick et al. (2018), and the method we use might introduce a small bias. This might explain some of the discrepancies between analytical and numerical results but should not invalidate the approach.

The SIS dynamic process in static networks can be simulated as follows Ferreira et al. (2012); Cota and Ferreira (2017): First, we build a list $\nu$ for infected nodes and calculate the total transition rate $\omega=\sum_{i\in I}(g+rk_{i})$ at the initial state, where $k_{i}$ is the degree of node $i$. The list $\nu$ and the total transition rate $\omega$ are constantly updated in one simulation, and the time is incremented by $dt=1/\omega(t)$. And then, there are two possible events happening at time $t+dt$: (1) With probability $gI/\omega$, an infected node $i$ is chosen with equal probability from $\nu$ and recovered. (2) With probability $\sum_{i\in I}(rk_{i})/\omega$, an infected node $i$ is chosen from $\nu$ at random and accepted with probability $k_{i}/k_{\max}$, which is repeated until one choice is accepted. Finally, a neighbor of $i$ is chosen randomly. If the neighbor is infected, nothing happens; otherwise it becomes infected. We iterate the whole process longer than it typically takes for the system to reach the steady-state.

The simulations in this paper were performed using the quasi-stationary (QS) method Ferreira et al. (2012), in which every time the system tries to visit the absorbing state it jumps to one of the pre-stored active configurations. These pre-stored active configurations are updated constantly, i.e. , a pre-stored configuration is chosen randomly and replaced by the present active configuration with a probability $p_{r}dt$. After a relaxation time $t_{r}$, the QS probability $\overline{P}(n)$ that the system has $n$ infected individuals is computed during an averaging time $t_{a}$. From the distribution $\overline{P}(n)$, the moments of the activity distribution can be computed as $\left\langle\rho^{k}\right\rangle=\sum_{n}(n/N)^{k}\overline{P}(n)$. And then, the epidemic threshold can be obtained by the maximum value of the susceptibility $\chi$, whose value is defined as $\chi=N\left\langle\rho^{2}\right\rangle-\left\langle\rho\right\rangle^{2})/\left\langle\rho\right\rangle$. In this paper, the number of active configurations is set as 100, $p_{r}=0.02$, $t_{a}=3t_{r}$, and $t_{r}$ depends on $N$ and $\lambda$.

We denote $p_{k}$ and $q_{k}$, respectively, as the probabilities of reaching an arbitrary infected individual by following a randomly chosen edge from a susceptible and infected individual of degree $k$ Cai et al. (2016). To calculate the $p_{k}$ and $q_{k}$ on scale-free networks, our analysis takes five steps.

First, we choose a node $i$ whose degree is $k$, which is the only known condition. And then, we divide the whole scale-free network into two subnetworks. One is a star subnetwork (denoted star) where the center node is node $i$. The other is the scale-free network minus the star subnetwork (denoted minus-star). Note that we do not ignore triangles. For any triangle $[ijj^{\prime}]$ where the nodes $j$ and $j^{\prime}$ are the neighbors of $i$, the edges $[ij]$ and $[ij^{\prime}]$ belong to the former subnetwork, and the edge $[jj^{\prime}]$ belongs to the latter subnetwork.

Second, from the star, we can get the probability $A$ that any node $j$ (a neighbor of $i$) is infectious. However, the calculation of this probability currently does not consider the node $j$’s neighbors except node $i$.

Third, for a large $k$, node $i$ has many neighbors. For a small $k$, there are many nodes with degree $k$ on a scale-free network. That is, there are many nodes $j$ regardless of $k$. For the minus-star, we can, by mean-field theory, get the average probability $B(k^{\prime})$ that any node $j^{\prime\prime}$ with degree $k^{\prime}$ is infectious. We can get the joint probability $C$ that the node $j^{\prime\prime}$ links the node $i$ and the node $j^{\prime\prime}$ is infected.

Fourth, we do not need to calculate these exact results $A$ and $C$, because they will affect each other on the whole scale-free network. Here, we assume approximately $C$ converges to $x$. And then, we preset that a fraction $x$ of the nodes are infected nodes all the time (they are fixed in the infection state) on a star network where the degree of its center node is $k$.

Fifth, $p_{k}$ can be calculated on the star network with permanently infected nodes (that we mentioned in Step 4), but it still contains an unknown $x$. We solve the $x$ by using the other global relation, i.e. Eq. (13). As a result, the probability $p_{k}$ (or $q_{k}$) can be divided into two parts: one from the nearest-neighbor network structure of the chosen node itself; the other representing the dynamic correlation coming from the whole network except the chosen node. The latter can be approximated as a fixed value on uncorrelated networks. Suppose all neighbors of one center node are identical, we consider a star subnetwork with a fraction $x$ of permanently infected nodes to calculate the probabilities $p_{k}$ and $q_{k}$, where $k$ is the degree of the star. Note that $x$ represents a joint probability, which cannot be replaced by the prevalence $\rho$ of scale-free networks.

In Fig. 1, we show a schematic illustration of one state cycle of the center node. We denote the number of infected leaf nodes at time $t$ by $I(t)$. To make it easier for the reader to understand our model, we summarize our ideas with simple words before the formula derivation begins: First, we derive the relationship between $I(t)$ and the joint probability $x$ on a star subnetwork with nonzero permanently infected leaf nodes at any given time (see Eqs. (1)–(3)); Second, according to the definition of $p_{k}$ (i.e. the relationship between $I(t)$ and $p_{k}$), we can estimate the relationship between $p_{k}$ and $x$ (see Eqs. (III.1)–(11)); Third, we solve $p_{k}$ and $x$ by introducing another global law, Eq. (14).

Since both $I(0)$ and $I(\tau_{1})$ have peaked distributions (see Appendix B), we use their expected values as approximations. From Fig. 1, we assume that the center node is infected at time $0$, and there are $n$ infected leaf nodes, i.e., $I(0)=n$. After a time interval $\tau_{1}$, the center node is cured but the number of infected leaf nodes rises to $m$ (so $I(\tau_{1})=m$). Then, $I(t)$ in the recovery process is determined by

where the second and third terms on the right hand side are the number of newly infected, and recovered, leaf nodes in the time interval $[0,t]$, respectively. After a time $\tau_{2}$, the center node becomes reinfected, and the number of infected leaf nodes reduces to $I(\tau_{1}+\tau_{2})=n$. $I(t)$ in the infection process is

where the second term on the right hand side is the number of newly recovery leaf nodes for times in $[\tau_{1},\tau_{1}+t]$. Combining Eqs. (1) and (2) gives

We define $P(\tau_{1})$ as the probability density of the duration $\tau_{1}$ in the recovery process. Since the total recovery rate is $g\tau_{1}$ during the time interval $[0,\tau_{1}]$, the center node’s recovery probability can be calculated as $1-e^{-g\tau_{1}}$ in this interval Ferreira et al. (2016). Then, we obtain $P(\tau_{1})=ge^{-g\tau_{1}}$, by the definition of cumulative distributions. We define $P(\tau_{2})$ as the probability density of the duration $\tau_{2}$ in the infection process. We can calculate the infection probability of the center node as $1-\exp^{-rI(t)dt}$ in a time interval $[t,t+dt]$. The cumulative infection probability of time interval $[\tau_{1},\tau_{1}+\tau_{2}]$ can be calculated as $1-\prod_{t=\tau_{1}}^{\tau_{1}+\tau_{2}}\exp^{-rI(t)dt}$. Combining this with Eq. (3), we can obtain the probability density $P(\tau_{2})$,

Moreover, combining Eq. (3) and Eq. (III.1), we can further find the relation below,

Eq. (5) applies to any star subnetwork with non-zero permanently infected leaf nodes, which implies that the disease in Fig. 1 will never become extinct.

Next, we average both sides of Eq. (2) by $\tau_{2}$ and set $t=\tau_{1}+\tau_{2}$. Considering Eq. (5), the expected decrease of $I(t)$ is $g/r-xkg\langle\tau_{2}\rangle$ in the infection process. Thus, we have

where $\langle\tau_{2}\rangle=\int\tau_{2}P(\tau_{2})d\tau_{2}$ is the average duration of infection. Combining Eq. (3) with Eq. (6) gives the number of infected leaf nodes at time $t=\tau_{1}$

Finally, according to the definition of $p_{k}$ and $q_{k}$,

we can get

where $\langle 1-e^{-g\tau_{2}}\rangle=\int(1-e^{-g\tau_{2}})P(\tau_{2})d\tau_{2}$. Here, we are able to rewrite $p_{k}$ to a simpler formula $1/(rk\langle\tau_{2}\rangle)$ by using the Eq. (5), but this formula is sensitive to the upper limit of the integral, where the upper limit is set as $\tau_{2}=50$ in this paper.

Moreover, we can obtain the following inequality between $p_{k}$ and $q_{k}$ by combining Eq. (7) and Eq. (9),

Ref. Cai et al. (2016) also presents this relation, but without a derivation from an underlying theory.

As we are modeling a stochastic process, there might be no infected leaf nodes other than the permanently infected ones. This situation becomes typical for small $k$ or small $\lambda$, which means that $I(t)<xk$. When this situation occurs, Eq. (6) becomes inaccurate because of the average loss of $I(t)$ might not be $g/r-xkg\langle\tau_{2}\rangle$. We can circumvent this problem by using the lower boundary $I(0)=xk$ for small $k$ or small $\lambda$. In particular, by replacing the Eqs. (6)–(7) with $n=xk$, we obtain new estimates for $p_{k}$ and $q_{k}$:

Here, Eq. (11) is not appropriate when $x=0$, because of the high chance of extinction.

Now, let us turn to the SIS dynamics on scale-free networks. When a scale-free network is in its epidemic state, we assume, for theoretical purposes, that there is a fixed fraction $x>0$ of infected nodes. As long as the epidemic process is in its quasi-stationary Nåsell (2011) state, the total recovery rate of $I_{k}$ must be equal to the total infection rate of $S_{k}$ Pastor-Satorras and Vespignani (2001). That means

where $N_{k}=S_{k}+I_{k}$ and $S_{k}$ ($I_{k}$) is the number of susceptible (infected) individuals with degree $k$. In the framework of DC Cai et al. (2016), we still require that the expected variations in the total recovery rate and total infection rate are equal in the steady state:

where $\omega=\sum_{k}\sum_{j}(S_{k,j}jr+I_{k,j}g)$ is the total rate and $S_{k,j}$ ($I_{k,j}$) is the number of susceptible (infected) nodes which have $j$ infected nodes among the total $k$ neighbors. Following Ref. Cai et al. (2016), this reduces to

For a given pair of $r$ and $g$ on a scale-free network, numerical calculation is divided into two major steps. First, we get $p_{k}$ for a preset $x$ by using Eqs. (III.1)–(11). Then, combining $p_{k}$ with the Eq. (14), we confirm whether the preset $x$ is a solution. In particular, we are able to calculate the probability $p_{k}$ in the steady state as follows: (i) Give any $x^{\prime}\in(0,1]$, the values of $m$ and $n$ are calculated by Eqs. (III.1)–(7); (ii) The values of various $p_{k}$ are calculated by Eq. (9) or Eq. (11); (iii) Plugging the set of $p_{k}$ into Eq. (14), and the value $x$ is the solution when Eq. (14) is valid. (iv) Inserting the $x$ into Eq. (9) and Eq. (11) to solve $p_{k}$.

In Fig. 2, we plot $p_{k}$ as a function of degree for several infection rates $r$ for the SIS model on a scale-free network with $\gamma=4.5$, $N=10^{7}$, and $k_{\max}=314$. As $r$ increases, the $k_{\max}$-star subnetwork (the star subnetwork of the largest-degree node) will become active before the other large-degree star networks. As $r$ increases, more star subnetworks will become active until the system reaches the global epidemic threshold ($x\to 0$). We show that our theoretical method captures $p_{k}$’s properties beyond the global threshold in Fig. 2(c)–Fig. 2(e) and near it in Fig. 2(b). It only fails for the smallest infection rates, see Fig. 2(a).

For sufficiently low infection rates, the epidemic is localized around one or a few high-degree star subnetworks in Fig. 2(a). Our theoretical method only captures two points, which implies that these two stars may be approximately continuously activated. By the actual SIS dynamics, the epidemics at any star network would go extinct after a average lifetime $T(k,r,g)$ Boguñá et al. (2013); Castellano and Pastor-Satorras (2020). So, there are many long intervals on star subnetworks between deaths and reactivations, which leads to $p_{k}\to 0$. Meanwhile, the mean-field treatment in Eqs. (1)–(3) is no longer appropriate in localized spreading. These contrast with when $r$ is larger, and many star subnetworks are active. Then the infected nodes are, as assumed by mean-field theory, more evenly distributed across the network. This explanation also implies HMF, PHMF, and DC are not suitable for localized spreading. They mistakenly give a finite threshold on scale-free networks with $\gamma>3$, but their threshold expressions are still valuable on networks with no localized spreading, such as random regular networks, Erdős-Rényi networks, and scale-free networks with $\gamma<3$.

Finally, combining these $p_{k}$ with Eq. (12), we can obtain the more accurate prevalence $\rho$ in population-wide outbreaks ($x>0$), that is,

where $\lambda=r/g$ is the effective infection rate and $P(k)$ is the degree distribution of the network. In Fig. 3, we show numerical and approximate results of the SIS dynamics on a scale-free network with $N=10^{7}$ and $\gamma=4.5$ and $3.5$. We find that our theoretical approach’s estimations of prevalence match those from stochastic simulations well on scale-free networks and are much better than the results predicted by the heterogeneous mean-field theory and the dynamic correlation method Pastor-Satorras and Vespignani (2001); Cai et al. (2016). Here, we only compare HMF and DC in Fig. 3 because they do not require the iterations of a large body of equations with respect to time (where PQMF, PHMF, ED does) to obtain the epidemic prevalence, which is the same as our theory. In Fig. 3(e), we show the ratio between the epidemic prevalence $\rho$ and the joint probability $x$ as a function of infection rate $r$. Clearly, $x$ fails to approximate $\rho$ for low enough $r$.

The successive activation of star subnetworks with increasing $r$ is a natural explanation for the multi-peak phenomenon of the susceptibility $\chi$ observed in the literature Mata and Ferreira (2015); Cota et al. (2016). If this scenario is correct, the separation of the hubs should affect the location of the peaks. Comparing with the double random regular network of Ref. Mata and Ferreira (2015) which is formed by two random regular networks connected by a single edge, we consider a network consisting of two stars where the two center nodes are separated by $l$ edges in Fig. 4(a).

In order to deal with localized spreading, we use the local condition (separated distance $l$) instead of those global conditions in Eq. (12) and Eq. (14). For $l=1$, the $x$ of $k_{2}$-star approximately equal to $1/k_{2}$ which is a finite value when the infection rate is slightly larger than the threshold of $k_{1}$-star. According to the Eq. (5), $k_{2}$-star is also active when $k_{1}$-star is active, which means that there is just one peak. When we set $x=0$ in Eqs. (6)–(7), the threshold of $k_{1}$-star can be easily predicted as an isolated star network, that is

In Fig. 4(b), we can see that the Eq. (16) captures the first peak very well.

For $l\geq 2$, the $x$ of $k_{2}$-star is $\left(r/g\right)^{l-2}p_{k_{1}}/k_{2}$ which is very small and vanish for large $l$. After setting $x=0$, the Eq. (5) can be rewritten as

where $P(\tau_{2})=rm\exp^{-g\tau_{2}}\exp^{-\frac{rm}{g}(1-e^{-g\tau_{2}})}$. The Eq. (17) means that the isolated star network will die with a probability $Y=\left(1+\frac{rm}{g}\right)\exp^{-\frac{rm}{g}}$. And then, its average lifetime $T(k,r,g)$ can be calculated as

where $\langle\tau_{1}\rangle=1/g$ and $m(k)=(kr^{2}-g^{2})/[r(g+r)]$. In Fig. 5, we can see that the Eq. (III.2) is a more accurate prediction of average lifetime on $k$-star network than that given from Ref. Boguñá et al. (2013); Castellano and Pastor-Satorras (2020). Note that $T(k,r,g)$ depends on the individual values of $r$ and $g$, not the ratio $r/g$.

We denote $\eta(k)$ as the probability that any leaf node is in the infected state on an active $k$-star network, that is

The upper boundary of Eq. (19) is $\eta(k)<r/(g+r)$ by using the Eq. (10). And then, we can calculate the probability of $k_{2}$-star being active per unit time on the two star network in Fig. 4(a),

Note that the recovery process of the leaf nodes of $k_{1}$-star has been considered by $\eta(k_{1})$. The $f(k_{1},k_{2},r,g)\geq 1$ means that $k_{2}$-star subnetwork can continuously be active. Combining $f(k_{1},k_{2},r,g)=1$ and $\langle\tau_{1}\rangle\gg\langle\tau_{2}\rangle$ with Eqs. (III.2)–(20), we can predict the activation point of $k_{2}$-star subnetwork (the lower peak of $\chi$) as follows,

In Fig. 4(b), we can see that the activation point of the $k_{2}$-star subnetwork strongly depends on the distance $l$ between the center nodes, and the Eq. (21) captures them very well. As $l$ decreases, the lower peak moves toward small $r$ until extinction ($l=1$), which all can be predicted quantitatively by the localized analysis of our theory.

Finally, based on the localized analysis above, we can predict the localized prevalence as follows,

where $\Omega$ contains all star subnetworks, $\overline{I(k)}=\frac{\langle\tau_{1}\rangle}{\langle\tau_{1}\rangle+\langle\tau_{2}(k)\rangle}+k\eta(k)$ is the average number of infected individuals when the $k$-star subnetwork is active, the $F(k)H[\lambda-\lambda_{c_{1}}(k)]$ represents the activation probability of the $k$-star subnetwork. The $F(k)$ is a ramp function whose value is $\min\{f(k_{\max},k,r,g),1\}$ and the $H[\lambda-\lambda_{c_{1}}(k)]$ is a Heaviside step function whose value is zero for negative argument and one for positive argument. In Fig. 4(c), we can see that the estimations in Eq. (22) match those from stochastic simulations quite well.

Since epidemic thresholds are defined only in the $N\rightarrow\infty$ limit, the peaks of $\chi$ from our quasi-stationary simulations are, technically speaking, not true thresholds [$\lambda_{c_{1}}(\text{$k_{\max}$})<\lambda<\lambda_{c}^{\text{global}}$]. For uncorrelated static scale-free networks $\lambda_{c_{1}}(\text{$k_{\max}$})$ is related to $\lambda_{c}^{\text{global}}$ Castellano and Pastor-Satorras (2020) via

where we have used the natural degree cut-off Dorogovtsev and Mendes (2002), $k_{\mathrm{cut}}/k_{\min}\sim N^{1/(\gamma-1)}$ as an estimate of $k_{\max}$. Equation (23) means that $\lambda_{c_{1}}(\text{$k_{\max}$})$ decreases faster than $\lambda_{c}^{\text{global}}$ with increasing $N$, which leaves a finite range of $\lambda$ where epidemic localization can emerge so that the outbreak survives in the $k_{\max}$-star subnetwork, but not in the rest of the network. We also note that is $\gamma>3$, hubs will be separated in large networks since the probability that nodes of degrees $k_{1}$ and $k_{2}$ are connected is proportional to $k_{1}k_{2}/N\langle k\rangle$ which goes zero with $N$.

For networks with more homogeneous degree distributions, such as random regular networks and Erdős-Rényi networks, the largest degree is not large enough to ensure $\lambda_{c_{1}}(\text{$k_{\max}$})<\lambda_{c}^{\mathrm{DC}}$, where $\lambda_{c}^{\mathrm{DC}}=\langle k\rangle/(\langle k^{2}\rangle-\langle k\rangle)$ is the epidemic threshold predicted by DC or PHMF Cai et al. (2016); Mata et al. (2014). Thus, a $k_{\max}$-star subnetwork would not become active until $\lambda$ exceeds the true threshold. Similarly, for the scale-free network with $2<\gamma<2.5$, we also have the relations of $\lambda_{c}^{\mathrm{QMF}}<\lambda_{c_{1}}(\text{$k_{\max}$})$ and $\lambda_{c}^{\mathrm{DC}}<\lambda_{c_{1}}(\text{$k_{\max}$})$ (see Appendix C). Thus, also for these networks, we expect one susceptibility peak, not because the largest degree is too small, but because the largest-degree nodes are connected.

Note that the discussion of the scale-free network with $2<\gamma<2.5$ mentioned above is from the uncorrelated configuration model. For a degree-degree disassortative correlations network, low degree nodes act as bridges linking the hubs, i.e. , the situation similar to Fig. 4(a) can easily arise. Interestingly, we may see the multi-peak phenomenon on the degree-degree disassortative correlations of the scale-free network with $2<\gamma<2.5$.

We return the epidemic threshold prediction on a scale-free network with $\gamma>3$, and summarize the results in Fig. 6. First, we can predict the global peak by using the Eq. (14) with the global coupling $x\to 0$, and show that our results are less than the prediction of Castellano and Pastor-Satorras (2020) where the latter tends to zero in the thermodynamic limit. Second, we can predict the first localized peak by using the Eq. (16) and $k=k_{\max}$. Note that they almost coincide between our results (black dashes) and the prediction of PQMF (red lines). Third, we can predict the second localized peak by using the Eq. (21) whose parameters are set as $k_{1}=k_{\max}$, $k_{2}\sim k_{\max}$, and $l=\langle l_{i}\rangle$. Compared with $l$, the value of $k_{2}$ has less influence on $f(k_{1},k_{2},r,g)$ in Eq. (20), which is why we use an approximation $k_{2}\sim k_{\max}$ in here. The $l_{i}$ is the distance between a star subnetwork with center node $i$ and the $k_{\max}$-star subnetwork. We average the $l_{i}$ of nodes whose degree is greater than $k^{\prime}$, where $k^{\prime}$ satisfy $\sum_{k=k^{\prime}}^{k_{\max}}kP(k)N\sim 0.001N$. This value of $0.001N$ is selected to ensure that the number of continuously activated subnetworks is sufficient to generate fluctuations of prevalence, but not reach the global spreading. Although the prediction we gave has many rough approximations, it is able to roughly hold the position of the second localized peak in Fig. 6.

Finally, we would like to point out that the nodes with a small degree still cannot be continuously active in Fig. 2(b), though its infection rate is higher than the second peak of Fig. 6(a).