Multiple peaks patterns of epidemic spreading in multi-layer networks

TL;DR

This paper demonstrates that multi-peak epidemic spreading patterns are characteristic of multi-layer networks, arising from delays and degree differences between layers, with an edge-based theory explaining the phenomena.

Contribution

It reveals the origin of multiple epidemic peaks in multilayer networks due to delays and degree disparities, and develops a theoretical framework to explain these patterns.

Findings

Multiple peaks are linked to delays in spreading between layers.

Degree distribution differences influence epidemic peak patterns.

Edge-based theory accurately predicts numerical results.

Abstract

The study of epidemic spreading on populations of networked individuals has seen recently a great deal of significant progresses. A common point of all past studies is, however, that there is only one peak of infected density in each single epidemic spreading episode. At variance, real data from different cities over the world suggest that, besides a major single peak trait of infected density, a finite probability exists for a pattern made of two (or multiple) peaks. We show that such a latter feature is fully distinctive of a multilayered network of interactions, and reveal that actually a two peaks pattern emerges from different time delays at which the epidemic spreads in between the two layers. Further, we show that essential ingredients are different degree distributions in the two layers and a weak coupling condition between the layers themselves. Moreover, an edge-based theory is developed which fully explains all numerical results. Our findings may therefore be of significance for protecting secondary disasters of epidemics, which are definitely undesired in real life.

Equations

$$S(t)=\sum _{k=0}^{\infty }P(k)\theta (t{)}^k$$

$$\theta (t)={\Phi }^S(t)+{\Phi }^I(t)+{\Phi }^R(t)$$

$${\Phi }^S(t)=\frac{{\sum }_{k^{\prime }}{k}^{\prime }P(k)\theta (t{)}^{k^{\prime }-1}}{⟨k⟩}$$

$$\frac{d{\Phi }^R(t)}{dt}=\mu (1-\beta ){\Phi }^I(t)$$

$$\frac{d(1-\theta (t))}{dt}=\beta {\Phi }^I(t)$$

$$\frac{d\theta (t)}{dt}=-\beta {\Phi }^I(t)$$

$${\Phi }^R(t)=\frac{\mu (1-\theta (t))(1-\beta )}{\beta }$$

$$\frac{d\theta (t)}{dt}=-\beta \theta (t)+\beta \frac{{\sum }_{k^{\prime }}{k}^{\prime }P(k)\theta (t{)}^{k^{\prime }-1}}{⟨k⟩}+\mu (1-\theta (t))(1-\beta )$$

$$\frac{dR(t)}{dt}=\mu I(t),S(t)=\sum _{k=0}^{\infty }P(k)\theta (t{)}^k,I(t)=1-S(t)-R(t)$$

$${\theta }_{1,2}={\Phi }_{1,2}^S+{\Phi }_{1,2}^I+{\Phi }_{1,2}^R$$

$${\Phi }_{1,2}^S$$

$$\frac{d{\Phi }_{1,2}^R}{dt}$$

$$\frac{d{\theta }_{1,2}}{dt}$$

$${\Phi }_{1,2}^R$$

$${\dot{\theta }}_{1,2}$$

$${\dot{\theta }}_{1,1}=-{\beta }_{1,1}{\theta }_{1,1}+{\beta }_{1,1}\frac{{\sum }_{k_1,k_2}{k}_1{P}_1(k_1,k_2){\theta }_{1,1}^{k_1-1}{\theta }_{1,2}^{k_2}}{{\sum }_{k_1,k_2}{k}_1{P}_1(k_1,k_2)}+\mu (1-{\theta }_{1,1})(1-{\beta }_{1,1})$$

$${\dot{\theta }}_{2,1}=-{\beta }_{2,1}{\theta }_{2,1}+{\beta }_{2,1}\frac{{\sum }_{k_1,k_2}{k}_2{P}_1(k_1,k_2){\theta }_{1,1}^{k_1}{\theta }_{1,2}^{k_2-1}}{{\sum }_{k_1,k_2}{k}_2{P}_1(k_1,k_2)}+\mu (1-{\theta }_{2,1})(1-{\beta }_{2,1})$$

$${\dot{\theta }}_{2,2}=-{\beta }_{2,2}{\theta }_{2,2}+{\beta }_{2,2}\frac{{\sum }_{k_1,k_2}{k}_2{P}_2(k_1,k_2){\theta }_{2,1}^{k_1}{\theta }_{2,2}^{k_2-1}}{{\sum }_{k_1,k_2}{k}_2{P}_2(k_1,k_2)}+\mu (1-{\theta }_{2,2})(1-{\beta }_{2,2})$$

$${\dot{R}}_1=\mu I_1(t),S_1(t)=\sum _{k_1,k_2}^{\infty }{P}_1(k_1,k_2){\theta }_{1,1}^{k_1}{\theta }_{1,2}^{k_2},I_1(t)=1-S_1(t)-R_1(t)$$

$${\dot{R}}_2=\mu I_2(t),S_2(t)=\sum _{k_1,k_2}^{\infty }{P}_2(k_1,k_2){\theta }_{2,2}^{k_2}{\theta }_{2,2}^{k_1},I_2(t)=1-S_2(t)-R_2(t)$$

Full text

{affiliations}

Department of Physics, East China Normal University, Shanghai, 200241, China

(* zhliu@phy.ecnu.edu.cn)

Web Sciences Center, University of Electronic Science and Technology of China, Chengdu 610054, China

School of Information Science and Technology, East China Normal University, Shanghai, 200241, China

CNR-Institute of Complex Systems, Via Madonna del Piano, 10, 50019 Sesto Fiorentino, Florence, Italy

The Embassy of Italy in Tel Aviv, 25 Hamered Street, 68125 Tel Aviv, Israel

Epidemic spreading in networked populations has been studied intensely in the last decade, and a lot of great progresses has been achieved [1, 2, 3, 4] which significantly increased our understanding. This is actually useful for public health authorities to assess situations quickly, to take and enforce informed decisions, and to optimize vaccination and drug delivery policies. Initially, the main attention focused on static networks, where each node represents an immobile agent and the contagion occurs only between neighboring nodes: it was remarkably revealed that scale-free networks display a vanishingly small epidemic threshold in the thermodynamic limit [5, 6]. After that, the focus shifted to reaction-diffusion models [7, 8, 9, 10], flow-driven epidemics [11, 12, 13, 14, 15], objective spreading [16, 17] and adaptive behaviors [18, 19, 20, 21, 22]. Finally, in a third stage, multilayered [23, 24, 25, 26, 27, 28, 29, 32, 30, 31] and temporal [33, 34] networks were assumed to play a critical role on such processes.

A common feature in past studies is the use of phenomenological models which produce a single peak of infected density in each individual epidemic spreading. An interesting question is therefore whether or not all real evolutionary processes are conveniently represented by such a framework. A scrupulous analysis of a large number of real data from different cities over the world surprisingly shows that, besides a major pattern made of a single peak, there is a finite and non negligible probability for a new pattern of epidemic outbreak featuring two (or multiple) peaks. Notice that a two peaks pattern implies two outbreaks in a single spreading period, i.e. a secondary occurrence of the same epidemics, which may in turn produce severe calamities and disasters within unprepared populations. Understanding the underlying mechanism at the basis of this new pattern (with the help of a novel model extracted from real data) is therefore quintessential to properly cope with such life-threatening hazards.

We proceed by making three steps: the first is to build a suitable model for the data, the second is to reproduce the data by the model, and the third is to suggest effective ways for predictions and control of the epidemic spreading. As social activities and interactions occur in network structures, we here consider the epidemics in different geographic regions (or cities) as that occurring in multilayered graphs. Namely, we will take two coupled neighboring regions as an example, and construct a two-layered network model which fully reproduces the observed patterns of epidemics. In particular, we demonstrate that the pattern of two peaks originates from large time delays of epidemic outbreaks between the two layers, which depends in turn on both the difference in the degree distributions of the two layers and a weak coupling condition between them. To better understand the findings, an edge-based theory is developed which perfectly agree with the numerical simulations.

Results

0.1 Patterns of epidemic outbreak with two or multiple peaks in real data.

Monitoring the potential outbreaks of an epidemic spreading is of extreme importance for protection of our society. Based on the detected trend of spreading of infections such as SARS (Severe Acute Respiratory Syndrome), H1N1 (Swine Influenza), H5H1 (Avian Influenza), and Ebola, one can indeed attempt to enforce suitable measures able to reduce the epidemic at its maximum extent. For this purpose, many countries have established their sentinel surveillance systems to collect epidemic data. For instance, Hong Kong Department of Health has organized a surveillance system, with the aim of collecting empirical data of infectious diseases, and of analyzing and predicting the trend of the infection. In such a system, there are about $64$ General Out-Patient Clinics (GOPC) and $50$ General Practitioners (GP), which form two distinct sentinel surveillance networks of the city [35, 36]. In these two networks, one obtains for instance the weekly consultation rates of influenza-like illness (per $1,000$ consultations), which reflect the overall influenza-like illness activity in Hong Kong.

Fig. 1(a) shows the collected data from $1998/1/3$ to $2014/8/2$ in GP, while the corresponding data of GOPC is shown in Fig. 1 of the Supplementary Information (SI). From Fig. 1(a) one see that there are many events of epidemic spreading, and the intervals between two consecutive events are not regular, indicating non-periodic outbreaks of recurrent epidemics [36]. On the other hand, one notices from Fig. 1(a) that most of the outbreaks correspond to a single peak of infected density, which is the pattern well described by the classical susceptible-infected-refractory (SIR) models. However, one notices also that there is a finite probability for a novel pattern of epidemic outbreak which features, instead, two or multiple peaks. Fig. 1(b) shows one of such patterns (with two peaks) occurring at around $2005/6$, which indicates that an infectious disease raised two times during that epidemic period in Hong Kong. Such unexpected phenomenon also exists in the data from GOPC (as one can see in Fig. 1 of the SI).

These multiple peaks patterns are actually occurring generically, and are not limited to a specific geographical region. Remrkably, indeed, epidemic data from other sources and cities display ubiquitously patterns similar to that reported in Fig. 1(b). For instance, Fig. 1(c) shows the data of the weekly measles infective cases (WMICs) from 1908 to 1937 in Boston [37, 38], and a two peaks pattern is shown in Fig. 1(d) at around 1915/4. In addition, two or multiple peaks of infected density may characterize the outbreaks in the total number of WMICs of two neighboring cities, also when uni-modal patterns are actually observed in each individual city. For example, Fig. 1(e) and (f) report the data of WMICs in Bristol and Newcastle [39]: the yellow line denotes the total number of WMICs, whereas the blue and green lines represent the data in Bristol and Newcastle, respectively. In Fig. 1(f) one can well appreciate that the pattern of two peaks occurs only in the total number of WMICs. Similarly, Fig. 1(g) and (h) show the case of Bristol and Sheffield [39], and once again a typical two peaks pattern [Fig. 1(h)] occurs.

0.2 The two-layered network model.

To capture the underlying mechanism, we introduce a model of a two-layered network, where the two layers represent actually two interconnected regions or cities. Fig. 2 is a sketch of the model: $\mathcal{A}$ and $\mathcal{B}$ are the two layers, which are coupled through the inter-network $\mathcal{AB}$. For the sake of simplicity, we let the two networks $\mathcal{A}$ and $\mathcal{B}$ have the same size $N_{a}=N_{b}$. Furthermore, $\langle k_{a}\rangle$, $\langle k_{b}\rangle$, and $\langle k_{ab}\rangle$ represent the average degrees of networks $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{AB}$, respectively (see Methods for details).

Each node is a unit of the SIR model, where $S,I$ and $R$ represent the susceptible, infected and refractory phases of individuals, respectively. At each time step, a susceptible node will be infected by an infected neighbor with rate $\beta$, and an infected node will become refractory with probability $\mu$. The infectious process will be considered terminated when no more infected nodes exist. While $\mu$ is taken to be the same for all networks, we let $\beta_{a}$, $\beta_{b}$ and $\beta_{ab}$ be the infectious rates of networks $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{AB}$, respectively.

A pattern of two peaks appears in the numerical simulations of our model. We choose $\mathcal{A}$ to be a scale-free (SF) network with degree distribution $P_{A}(k)\sim k^{-\gamma}$ [40], and $\mathcal{B}$ a random regular (RR) network with a constant degree $k_{b}$ [41]. The inter-network $\mathcal{AB}$ is constructed by randomly adding links between $\mathcal{A}$ and $\mathcal{B}$ until an averaged degree $\langle k_{ab}\rangle$ is attained. In simulations, we fix $N_{a}=N_{b}=10\,000$, and set initially $0.1\%$ of the individuals in $\mathcal{A}$ to be infected. The yellow circles in Fig. 3(a) report the evolution of the infected density $\rho_{I}$ in the whole network with the parameters $\gamma=2.1$, $\langle k_{a}\rangle=\langle k_{b}\rangle=6$, $\langle k_{ab}\rangle=1.0$, $\beta_{a}=\beta_{b}=0.05$, $\beta_{ab}=0.005$ and $\mu=0.1$. One can easily differentiate two peaks of $\rho_{I}$, indicating that the empirical observations in Fig. 1 can be fully reproduced.

To gather a deeper understanding, we also measure the evolutions of $\rho_{I}^{A}$ and $\rho_{I}^{B}$ in both layers $\mathcal{A}$ and $\mathcal{B}$, and report them as green triangles and blue squares in Fig. 3(a), respectively. It is easy to see that the times at which the maximum infected density is obtained in layers $\mathcal{A}$ and $\mathcal{B}$ are different, implying that the pattern of two peaks is likely triggered by the time difference of epidemic outbreaks in the two layers.

The next step is to focus on the key factors that determines the occurrence of such two peaks pattern. For this purpose, we first concentrate on the role of the average degree $\langle k_{ab}\rangle$ of the inter-network $\mathcal{AB}$. Fig. 3(b) shows the evolution of the infected density $\rho_{I}$ with different $\langle k_{ab}\rangle$ (with triangles, squares and circles denoting the cases of $\langle k_{ab}\rangle=3.0,1.0$ and $0.5$, respectively). One notices that the pattern of $\rho_{I}$ is uni-modal when $\langle k_{ab}\rangle$ is large, but bimodal when $\langle k_{ab}\rangle$ is sufficiently small, indicating that $\langle k_{ab}\rangle$ is a key factor for the appearance of a bimodal pattern: a smaller $\langle k_{ab}\rangle$ favours the appearance of the two peaks pattern, indicating that the two main networks $\mathcal{A}$ and $\mathcal{B}$ should be only weakly coupled among them.

As a second step, we study the influence of the infectious rate $\beta_{ab}$ on the pattern. Fig. 3(c) reports the results obtained for different $\beta_{ab}$ (with triangles, squares and circles representing the cases of $\beta_{ab}=0.05,0.005$ and $0.001$, respectively). Once again one may notice that the condition of a weak coupling is essential for a bimodal pattern: $\rho_{I}$ is indeed uni-modal when $\beta_{ab}$ is large, while the two peaks appear when $\beta_{ab}$ is small. Finally, we study the influence of the exponent $\gamma$ of the SF network. Fig. 3(d) shows that the bimodal feature is reduced with the increase of $\gamma$. As a larger $\gamma$ means a smaller difference between the structures of the SF and RR networks, one can infer that also the heterogeneity between the two layers is a key factor for the appearance of two peaks patterns.

All the numerical results are fully confirmed by an edge-based compartmental theory, see our theoretical Eqs. (18) and (19) in Methods. The solid curves in Fig. 3(a)-(d) show the corresponding theoretical results.

These first numerics point that both a weak coupling and a difference in heterogeneity between the two layers are necessary conditions for the emergence of the new pattern. An interesting question is whether the observed behavior corresponds to a critical phenomenon. To figure out the answer, we let $\tau$ be the time interval between the two peaks of $\rho_{I}$. In particular, the two peaks will merge into a single one when $\tau=0$. Similarly, we let $\delta_{t}=\mid t^{B}_{max}-t^{A}_{max}\mid$ be the time delay between the two peaks in layers $\mathcal{A}$ and $\mathcal{B}$, where $t^{A}_{max}$ and $t^{B}_{max}$ are the times of occurrence of the peak in $\rho_{I}^{A}$ and $\rho_{I}^{B}$, respectively. The trivial situation would be that for which $\tau=\delta_{t}$, but our numerical simulations show that this condition is attained only when $\tau$ is large, whereas one has $\tau<\delta_{t}$ when $\tau$ is sufficiently small. And, in particular, one has $\delta_{t}>0$ when $\tau=0$.

Fig. 4(a) and (b) show the dependence of $\tau$ and $\delta_{t}$ on $\beta_{ab}$ for fixed $\gamma=2.1$ and different $\langle k_{ab}\rangle$, respectively. From Fig. 4(a), one sees that when $\langle k_{ab}\rangle$ is small, $\tau$ will decrease monotonically and be non-vanishing with the increasing of $\beta_{ab}$, indicating that the event of two peaks always exists in the pattern. However, when $\langle k_{ab}\rangle$ is increased, $\tau$ will decrease rapidly to zero, implying that there is a critical $\beta_{ab}^{c}$ for different $\langle k_{ab}\rangle$. When $\beta_{ab}<\beta_{ab}^{c}$, the epidemic spreading event occurs through a pattern of two peaks in the infected density, while for $\beta_{ab}>\beta_{ab}^{c}$ it occurs via the traditional pattern with a single peak. On its turn, Fig. 4(b) shows that $\delta_{t}$ decreases monotonically with the increase of $\beta_{ab}$ for all the three cases of $\langle k_{ab}\rangle$, and it never vanishes. This is because the spreading speed is different in homogeneous and heterogeneous networks. Generally speaking, epidemic spreading is faster in heterogeneous network than in homogeneous network [42].

We then move to investigate the influence of the heterogeneity in degree distribution on the occurrence of the two peaks pattern. Fig. 4(c) and (d) show the dependence of $\tau$ and $\delta_{t}$ on $\beta_{ab}$ for different $\gamma$ and fixed $\langle k_{ab}\rangle=1.0$. While the network heterogeneity decreases with the increasing of $\gamma$, it is easy to see that when $\gamma$ is large, $\tau$ and $\delta_{t}$ decrease more prominently with $\beta_{ab}$. Specifically, the difference of spreading speeds between the two layers is not distinctive for large $\gamma$, resulting in the disappearance of the two peaks. Therefore, increasing the coupling strength (i.e. $\langle k_{ab}\rangle$ and $\beta_{ab}$) or decreasing the heterogeneity of network topology between the networks $\mathcal{A}$ and $\mathcal{B}$ will decrease the time delay of epidemic outbreak and then suppress the pattern of two peaks.

We have also confirmed all these numerical results in Fig. 4(a)-(d) by the theoretical Eqs. (18) and (19) in Methods. For each set of parameters in Fig. 4(a)-(d), we first produce the infected densities $\rho_{I},\rho_{I}^{A}$, and $\rho_{I}^{B}$ from Eqs. (18) and (19), as done in Fig. 3(a)-(d), and then measure the corresponding $\tau$ and $\delta_{t}$. The solid curves in Fig. 4(a)-(d) show the theoretical $\tau$ and $\delta_{t}$. One can easily see that the theoretical results are fully consistent with the numerical results.

Discussion

Let us remark that the weak coupling condition predicted by us for the occurrence of the novel epidemic pattern is actually consistent with the cases of the data of Fig. 1. As it is well known, indeed, Hong Kong in Fig. 1(a) consists of islands, and the movement of individuals between different islands is not as convenient as that within each single island, and thus the coupling between neighboring islands can be considered to be weak. At the same time, the population distribution in Hong Kong central island is significantly different from that characterizing the surrounding islands, confirming the presence of the second ingredient predicted by our theory: i.e. a difference in the heterogeneity of the layers’ structures.

In Boston, a river separates the city into two parts, which (to all extent) can be considered as equivalent to two islands. The same reasoning applies to the neighboring cities of Bristol and Newcastle and the neighboring cities of Bristol and Sheffield. As Bristol and Newcastle, Newcastle and Sheffield are separated regions in the United Kingdom, and they can therefore be considered as a pair of weakly coupled networks.

Our predictions were obtained on coupled SF-RR networks, and it is legitimate to seek for generality of the two peaks pattern phenomenon, by means of investigating coupled networks with other topological structures. For this purpose, we have also studied the case of SF-SF and RR-RR networks, respectively. Very interestingly, one finds that the pattern of two peaks can be still observed by adjusting the coupling strength between the coupled layers (see Fig. 7 in SI for details). On the other hand, extension to three-layered model was also considered, and it was found that there is a small probability to produce a pattern of three peaks (see Fig. 8 in SI for details). Therefore, while in principle one can expect a multi-peaks pattern to occur in a multilayered network, the majority of unusual cases (i.e. cases in which the epidemic event is not happening with a single maximum of infected density) will be characterized by just two peaks, in full consistency with the data of Fig. 1.

In summary, epidemic spreading has been well studied in the past decades but mainly focused on outbreaks corresponding to patterns with a single peak of infected density. We here reported (from real data) the evidence that also a pattern of two peaks in a single epidemic period is possible. We pointed out that such a pattern is a genuine product of a multi-layered interaction structure, and we have introduced a proper model able to fully capture the mechanisms for its occurrence. Our model, together with reproducing the classical pattern of a single peak, can generate the pattern with two peaks when proper conditions on weak coupling between the layers and difference in heterogeneity of the layers’ structures are satisfied.

Methods

0.3 A two-layered network epidemic model.

We consider a two-layered network model with coupling between its two layers, i.e. the networks $\mathcal{A}$ and $\mathcal{B}$ in Fig. 2. We let the two networks have the same size $N_{a}=N_{b}=N$ and their degree distributions $P_{A}(k)$ and $P_{B}(k)$ be different. Each node has two kinds of links, i.e. intra-connection (within $\mathcal{A}$ or $\mathcal{B}$) and interconnection between $\mathcal{A}$ and $\mathcal{B}$. The former consists of the degree distributions of $P_{A}(k)$ and $P_{B}(k)$ while the latter gives rise to the interconnection network $\mathcal{AB}$. We let $\langle k_{a}\rangle$, $\langle k_{b}\rangle$, and $\langle k_{ab}\rangle$ represent the average degrees of networks $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{AB}$, respectively. In details, we first generate two separated networks $\mathcal{A}$ and $\mathcal{B}$ with the same size $N$ and different degree distributions $P_{A}(k)$ and $P_{B}(k)$, respectively. Then, we add links between $\mathcal{A}$ and $\mathcal{B}$. That is, we randomly choose two nodes from $\mathcal{A}$ and $\mathcal{B}$ and then connect them if they are not connected yet. The process is repeated until all the needed specifications are attained.

In the above way, one obtains a uncorrelated two-layered network.

To discuss epidemic spreading in such a framework, we let each node represent a SIR model. In this model, a susceptible node has two ways to be infected. One is from contacting with infected individuals in the network $\mathcal{A}$ (or $\mathcal{B}$), represented by $\beta_{a}$ (or $\beta_{b}$). The other is from the coupled network $\mathcal{AB}$, represented by $\beta_{ab}$ (see Fig. 2). Thus, a susceptible node will be infected with a probability $1-(1-\beta_{a})^{k^{inf}}(1-\beta_{ab})^{k_{ab}^{inf}}$, where $k^{inf}$ is the infected neighbors in the same network and $k_{ab}^{inf}$ is the infected neighbors in the coupled network. At the same time, an infected node will become refractory by a probability $\mu$.

In our numerical simulations, we use scale-free (SF) and regular random (RR) graphs as networks $\mathcal{A}$ and $\mathcal{B}$, respectively. The network size is $N_{a}=N_{b}=10\,000$, the average degree $\langle k_{a}\rangle=\langle k_{b}\rangle=6$, and initially $0.1\%$ of individuals of network $\mathcal{A}$ are chosen to be infected.

0.4 Edge-based compartmental theory for a single network.

Let us first illustrate the edge-based compartmental theory for a single network, by following the methods and tools introduced in Refs.[43, 44, 45, 46, 47, 48, 49, 50, 51, 52].

For an uncorrelated, large and sparse network, the SIR model can be described in terms of the quantities $S(t)$, $I(t)$, and $R(t)$, which represent the densities of the susceptible, infected, and recovered nodes at time $t$, respectively. Let $\theta(t)$ be the probability that a neighbor $v$ of $u$ has not transmitted the disease to $u$ along the edge connecting them up to time $t$. Then, the node $u$ with degree $k$ is susceptible at time $t$ as $s(k,t)=\theta(t)^{k}$. Averaging over all $k$, the density of susceptible nodes at time $t$ is given by

[TABLE]

where $P(k)$ is the degree distribution of the network. In order to solve for $S(t)$, one needs to know $\theta(t)$. Since a neighbor $v$ of node $u$ may be susceptible, infected, or recovered, $\theta(t)$ can be expressed as

[TABLE]

where $\Phi^{S}(t),\Phi^{I}(t),\Phi^{R}(t)$ is the probability that the neighbor $v$ is in the susceptible, infected, recovery state, respectively, and has not transmitted the disease to node $u$ through their connection. Once these three parameters can be derived, we will get the density of susceptible nodes at time $t$ by substituting them into Eq. (2) and then into Eq. (1). To this purpose, in the following, we will focus on how to solve them.

To find $\Phi^{S}(t)$, we now consider a randomly chosen node $u$, and assume this node is in the cavity state, which means that it cannot transmit any disease to its neighbors $v$ but can be infected by its neighbors. In this case, the neighbor $v$ can only get the disease from its other neighbors except the node $u$. Thus, node $v$ with degree $k^{\prime}$ is susceptible with probability $\theta(t)^{k^{\prime}-1}$ at time $t$ . For uncorrelated networks, the probability that one edge from node $u$ connects with a node $v$ with degree $k^{\prime}$ is $k^{\prime}P(k)/\langle k\rangle$. Summing over all possible $k^{\prime}$, one obtains

[TABLE]

According to the SIR spreading process, the growth of $\Phi^{R}(t)$ includes two consecutive events: first, an infected neighbor has not transmitted the infection to node $u$ via with probability $1-\beta$; second, the infected neighbor has been recovered with probability $\mu$. Combining these two events, the $\Phi^{I}(t)$ to $\Phi^{R}(t)$ flux is $\mu(1-\beta)\Phi^{I}(t)$. Thus, one gets

[TABLE]

Once the infected neighbor $v$ transmits the disease to $u$ successfully (with probability $\beta$), the $\Phi^{I}(t)$ to $1-\theta(t)$ flux will be $\beta\Phi^{I}(t)$, which means

[TABLE]

That is

[TABLE]

Combining Eqs. (4) and (5), and considering (as initial conditions) $\theta(0)=1$ and $\Phi^{R}(0)=0$, one obtains

[TABLE]

Substituting Eqs. (3) and (6) into Eq. (2), one gets an expression for $\Phi^{I}(t)$ in terms of $\theta(t)$, and then one can rewrite Eq. (5) as

[TABLE]

Thus, the equation of the system comes out to be

[TABLE]

In fact, Eq. (7) does not depend on Eq. (8), so the system is governed by the single ordinary differential equation (7). Although the resulting equation are simpler than those found by other methods, it can be proven to exactly predict the disease dynamics in the large-population limit for different network topologies[49, 53].

0.5 The theory for two-layered networks.

When one assumes that the population is made up of two interacting networks, then $P_{j}(k_{1},k_{2})$ denote the probability that a node of network $j$ has $k_{1}$ degree in network $1$ and $k_{2}$ in network $2$. For the sake of simplicity, one can name the two networks $\mathcal{A}$ and $\mathcal{B}$ as $1$ and $2$. Let $\beta_{j,l}$ be the rate of transmission across an edge from network $l$ to network $j$, and let us define $\mu$ to be the recovery rate of a node in any network.

$\theta_{j,l}$ can be defined to be the probability that an edge to a test node $u$ in network $j$ ($j=1,2$) is coming from network $l$ ($l=1,2$), and has not transmitted the infection.

Now, $\theta_{1,2}$ can be solved as in the case of a single network. Since a neighbor $v$ in network 2 of node $u$ in network 1 may be susceptible, infected, or recovered, $\theta_{1,2}$ can be expressed as

[TABLE]

where $\Phi^{S}_{1,2}$, $\Phi^{I}_{1,2}$, $\Phi^{R}_{1,2}$ is the probability that the neighbor $v$ is in the susceptible, infected, recovery state, and has not transmitted the disease to node $u$ through their connection.

Similarly, to find $\Phi^{S}_{1,2}$, the neighbor $v$ in network 1 can only get the disease from its other neighbors except the node $u$ in network 2. Thus, the node $v$ with degree $k_{1}$ in network 1 and degree $k_{2}$ in network 2 is susceptible with probability $\theta_{2,1}^{k_{1}-1}\theta_{2,2}^{k_{2}}$ at time $t$. For uncorrelated networks, the probability that one edge from node $u$ connects with a node $v$ with degree $(k_{1},k_{2})$ is $\frac{k_{1}P_{2}(k_{1},k_{2})}{\sum_{k_{1},k_{2}}k_{1}P_{2}(k_{1},k_{2})}$. Thus, one has

[TABLE]

It is easily to know that the growth of $\Phi^{R}_{1,2}$ includes two consecutive events: first, an infected neighbor has not transmitted the infection to node $u$ via with probability $1-\theta_{1,2}$; second, the infected neighbor has been recovered with probability $\mu$. Combining these two events, the $\Phi^{I}_{1,2}$ to $\Phi^{R}_{1,2}$ flux is $\mu(1-\theta_{1,2})\Phi^{I}_{1,2}$. Thus, one gets

[TABLE]

Once the infected neighbor $v$ in network 1 transmits the disease to node $u$ in network 2 successfully (with probability $\beta_{1,2}$), the $\Phi^{I}_{1,2}$ to $1-\theta_{1,2}$ flux will be $\beta_{1,2}\Phi^{I}_{1,2}$, which means

[TABLE]

Combining Eqs. (11) and (12), and considering (as initial conditions) $\theta_{1,2}(0)=1$ and $\Phi^{R}_{1,2}(0)=0$, one obtains

[TABLE]

So, one gets

[TABLE]

Similarly, one can write down $\theta_{1,1}$, $\theta_{2,1}$ and $\theta_{2,2}$ as follows

[TABLE]

With Eqs. (14-17) on hand, the densities associated with each distinct state can be obtained by

[TABLE]

Eqs. (18) and (19) are the main theoretical results from which the theoretical curves in Figs. 3 and 4 are calculated. Furthermore, we find that the threshold for the whole network to show epidemic outbreak can be theoretically figured out by the Jacobian matrix J of Eqs. (14-17) (see Figs. 5 and 6 in SI for details). Especially, when the coupling between the two layers is very weak, the obtained threshold will be consistent with the previous findings [54, 25].