Synchronized and mixed outbreaks of coupled recurrent epidemics.

Epidemic spreading has been studied for a long time and most of them are focused on the growing aspect of a single epidemic outbreak. Recently, we extended the study to the case of recurrent epidemics (Sci. Rep. 5, 16010 (2015)) but limited only to a single network. We here report from the real data of coupled regions or cities that the recurrent epidemics in two coupled networks are closely related to each other and can show either synchronized outbreak pattern where outbreaks occur simultaneously in both networks or mixed outbreak pattern where outbreaks occur in one network but do not in another one. To reveal the underlying mechanism, we present a two-layered network model of coupled recurrent epidemics to reproduce the synchronized and mixed outbreak patterns. We show that the synchronized outbreak pattern is preferred to be triggered in two coupled networks with the same average degree while the mixed outbreak pattern is likely to show for the case with different average degrees. Further, we show that the coupling between the two layers tends to suppress the mixed outbreak pattern but enhance the synchronized outbreak pattern. A theoretical analysis based on microscopic Markov-chain approach is presented to explain the numerical results. This finding opens a new window for studying the recurrent epidemics in multi-layered networks.

Introduction

Epidemic spreading in complex networks has been well studied and a lot of great progress have been achieved such as the infinitesimal threshold[1-6], reaction-diffusion model[7-10], flow driven epidemic[11-15], and objective spreading etc.[16, 17], see the review refs 18–20. for details. Recently, the attention has been moved to the case of multilayer networks[18, 21–31], which represent the interactions between different real-world networks such as critical infrastructure[32-34], transportation networks[35, 36], living organisms[37-39], and social networks[21, 40] etc. These models enable us to determine how the interplay between network structures influences the dynamic processes taking place on them[41-50]. For instance, a pathogen spreads on a human contact network abetted by global and regional transportation networks[21]. Due to their ubiquitous applications in complex systems[51-55], the understanding of the properties and dynamic processes in multilayered networks carries great practical significance.

Two of the most successful models used to describe epidemic spreading are the susceptible-infected-susceptible (SIS) and susceptible-infected-refractory (SIR) models. Mark et al. used the SIR model to multilayered networks in 2012[41]. Very interesting, they found a mixed phase in weakly coupled networks where an epidemic occurs in one network but does not spread to the coupled network. Saumell-Mendiola et al. used the SIS model to multilayered networks also in 2012[42]. However, they found that such a mixing phase doesn’t exist in both analytic and simulation results. In their work, they mainly focused on the epidemic threshold and studied how epidemics spread from one network to another.

All these studies are focused on the growing aspect of a single epidemic outbreak, no matter it is one layer or multilayered networks. However, in realistic situations, the empirical data shows that epidemic is recurrent, i.e. with outbreaks from time to time[56-60]. Thus, we recently extended the study to the case of recurrent epidemics[57], but limited only to a single network. Considering that the interactions between different networks are ubiquitous, we here recheck several real data of coupled regions or cities such as the General Out-Patient Clinics (GOPC) network and its coupled General Practitioners (GP) network of Hong Kong (see Fig. 1(a) and (b)), the coupled regional networks of California and Nevada (see Fig. 1(c) and (d)), the coupled regional networks of Arizona and California (see Fig. a(a) and (b) in SI), the coupled city networks of Boston and Fall River (see Fig. 1(c) and (d) in SI), and the coupled city networks of Los Angeles and Sacramento (see Fig. 1(e) and (f) in SI). We interestingly find that their recurrent epidemics are closely related to each other. Moreover, we find that the coupled time series of recurrent epidemics can show either synchronized outbreak pattern where outbreaks occur simultaneously in both networks or mixed outbreak pattern where outbreaks occur only in one network but do not in another one. This finding calls our great interest and motivates us to study its underlying mechanism. In this sense, we believe that it is very necessary to further extend the study of recurrent epidemics to the case of multilayered networks.

Figure 1

Time series of recurrent epidemics in two coupled regions or cities. (a) and (b) represent the weekly consultation rates of influenza-like illness (per 1000 consultations) from 1998 to 2014 in Hong Kong for the General Out-Patient Clinics (GOPC) and the General Practitioners (GP), respectively, where the data from 2009/6/13 to 2010/5/23 in (a) are not available and the value of C in (b) is from 0 to 150. (c) and (d) represent the time series of reported weekly measles infective cases I in California and Nevada, respectively.

In this work, we present a two-layered network model of coupled recurrent epidemics to reproduce the synchronized and mixed outbreak patterns. To guarantee the appearance of recurrent outbreaks, we choose the susceptible-infected-refractory-susceptible (SIRS) model for each node of network and let the infectious rate be time dependent, symbolizing the larger annual variation of environment. By this model, we find that the average degrees of both the intra- and inter-networks play key roles on the emergence of synchronized and mixed outbreak patterns. The synchronized outbreak pattern tends to be triggered in two coupled networks with the same average degree while the mixed outbreak pattern is likely to show for the case with different average degrees. Further, we show that the increasing of coupling strength, i.e. either the inter-layer infection rate or inter-layer average degree, will tend to suppress the mixed outbreak pattern but enhance the synchronized outbreak pattern. A theoretical analysis based on microscopic Markov-chain approach is presented to explain the numerical results. This finding may be of significance to the long-term prediction and control of recurrent epidemics in multi-areas or cities.

Results

The synchronized and mixed outbreak patterns of recurrent epidemics in real data

Monitoring epidemic spreading is vital for us to prevent and control infectious diseases. For this purpose, Hong Kong Department of Health launched a sentinel surveillance system to collect data of infectious diseases, aiming to analyze and predict the trend of infectious spreading in different regions of Hong Kong. In this system, there are about 64 General Out-Patient Clinics (GOPC) and 50 General Practitioners (GP), which form two sentinel surveillance networks of Hong Kong[56], respectively. By these two networks, we can obtain the weekly consultation rates of influenza-like illness (per 1000 consultations). Figure 1(a) and (b) show the collected data from 1998 to 2014 for the GOPC and GP, respectively, where the data from 2009/6/13 to 2010/5/23 in (a) was not collected by Hong Kong Department of Health and the value of C in (b) is from 0 to 150. This weekly consultation rates of influenza-like illness can well reflect the overall influenza-like illness activity in Hong Kong. From the data in Fig. 1(a) and (b) we easily find that there are intermittent peaks, marking the recurrent outbreaks of epidemics. By a second check on the data in Fig. 1(a) and (b) we interestingly find that some peaks occur simultaneously in the two networks at the same time, indicating the appearance of synchronized outbreak pattern. While other peaks appear only in Fig. 1(a) but not in Fig. 1(b), indicating the existence of mixed outbreak pattern (see the light yellow shaded areas).

Is this finding of synchronized and mixed outbreak patterns in recurrent epidemics a specific phenomenon only in Hong Kong? To figure out its generality, we have checked a large number of other recurrent infectious data in different pathogens and in different states and cities in the United States and found the similar phenomenon. Figure 1(c) and (d) show the data of weekly measles infective cases I in the states of California and Nevada, respectively, which were obtained from the USA National Notifiable Diseases Surveillance System as digitized by Project Tycho[58-60]. As California is adjacent to Nevada in west coast of the United States, their climatic conditions are similar. Thus, they can be also considered as two coupled networks. Three more these kinds of examples have been shown in Fig. 1 of SI where the coupled networks are based on states-level influenza data and cities-level measles data, respectively. Therefore, the synchronized and mixed outbreak patterns are general in recurrent epidemics. We will explain their underlying mechanisms in the next subsection.

A two-layered network model of coupled recurrent epidemics

Two classical models of epidemic spreading are the Susceptible-Infected-Susceptible (SIS) model and Susceptible-Infected-Refractory (SIR) model[18]. In the SIS model, a susceptible node will be infected by an infected neighbor with rate β. In the meantime, each infected node will recover with a probability γ at each time step. After the transient process, the system reaches a stationary state with a constant infected density ρ . Similarly, in the SIR model, each node can be in one of the three states: Susceptible, Infected, and Refractory. At each time step, a susceptible node will be infected by an infected neighbor with rate β and an infected node will become refractory with probability γ. The infection process will be over when there is no infected I. These two models have been widely used in a variety of situations. However, it was pointed out that both the SIS and SIR models are failed to explain the recurrence of epidemics in real data[57, 61].

To reproduce the synchronized and mixed outbreak patterns, we here present a two-layered network model of coupled recurrent epidemics, shown in Fig. 2 where (a) represents its schematic figure of network topology and (b) denotes the epidemic model at each node. In Fig. 2(a), two networks are coupled through some interconnections between them, which form the inter-network . For the sake of simplicity, we let the two networks have the same size N  = N . We let 〈k 〉, 〈k 〉, and 〈k 〉 represent the average degrees of the networks , respectively, see Methods for details. In Fig. 2(b), the epidemic model is adopted from our previous work[57] by two steps. In step one, we extend the SIR model to a Susceptible-Infected-Refractory-Susceptible (SIRS) model where a refractory node will become susceptible again with probability δ. In step two, we let each susceptible node have a small probability p 0 to be infected, which represents the natural infection from environment. Moreover, we let the infectious rate β(t) be time dependent, representing its annual and seasonal variations etc. To distinguish the function of the interconnections from that of those links in , we let β be the inter-layer infectious rate. In this way, the interaction between can be described by the tunable parameter β .

Figure 2

Schematic figure of the epidemic model to reproduce the synchronized and mixed outbreak patterns. (a) Schematic figure of the two-layered network, where the “black”, “blue” and “red” lines represent the links of the networks and the inter-network , respectively. β denotes the infectious rate through one interconnection between . (b) Schematic figure of the extended SIRS model for each node in , where the symbols S, I and R represent the susceptible, infectious, and refractory states, respectively, and the parameters β, γ and δ represent the infectious, refractory and recovery rates, respectively. p 0 denotes the probability for a susceptible node to be naturally infected by environment or other factors.

In numerical simulations, we let both the networks be the Erdös-Rényi (ER) random networks[62]. To guarantee a recurrent epidemics in each of , we follow ref. 57 to let β(t) be the truncated Gaussian distribution and choose p 0 = 0.01. Figure 3(a) and (b) show the evolutions of the infected density ρ in , respectively, where the parameters are taken as 〈k 〉 = 6.5, 〈k 〉 = 1.5, 〈k 〉 = 1.0, and β  = 0.09. It is easy to observe that some peaks of ρ appear simultaneously in , indicating the synchronized outbreak pattern. We also notice that some peaks of ρ in do not have corresponding peaks in (see the light yellow shadowed areas in Fig. 3(a) and (b)), indicating the mixed outbreak pattern. Moreover, we do not find the contrary case where there are peaks of ρ in but no corresponding peaks in , which is also consistent with the empirical observations in Fig. 1(a–d).

Figure 3

Reproduced time series of recurrent epidemics in the two-layered network model. (a) and (b) represent the evolutions of the infected density ρ in the two-layered networks , respectively, where the parameters are taken as 〈k 〉 = 6.5, 〈k 〉 = 1.5, 〈k 〉 = 1.0, , β  = 0.09, γ = 0.2, δ = 0.02, p 0 = 0.01, and N  = N  = 1000.

Mechanism of the synchronized and mixed outbreak patterns

To understand the phenomenon of synchronized and mixed outbreak patterns better, we here study their underlying mechanisms. A key quantity for the phenomenon is the outbreak of epidemic, i.e. the peaks in the time series of Fig. 1. Notice that a peak is usually much higher than its background oscillations. To pick out a peak, we need to define its background/baseline first. As the distributions of both the real data in Fig. 1 and numerical simulations in Fig. 3 are approximately satisfied the normal distribution (see Fig. 2 SI), we define the baseline as μ + 3σ with μ and σ being the mean and standard deviation, respectively, which contains about 99.7% data in the normal distribution[63]. Then, we can count the number of outbreaks in a measured time t. Let n be the average number of outbreaks in realizations of the same evolution t. Larger n implies more frequent outbreaks. Let be the difference of outbreak numbers between the networks . Larger Δn implies more frequent emergence of the mixed outbreak pattern. In particular, the mixed outbreak pattern will disappear when Δn = 0.

We are interested in how the average degrees and coupling influence the numbers n and Δn. Figure 4(a) and (b) show the dependence of n on the average degree 〈k 〉 of network for β  = 0.09 and 0.3, respectively, where the average degree of network is fixed as 〈k 〉 = 6.5. It is easy to see from Fig. 4(a) and (b) that the number of network will gradually increase with the increase of 〈k 〉, while the number of network keeps approximately constant, indicating that a larger 〈k 〉 is in favor of the recurrent outbreaks. Specifically, will reach when 〈k 〉 is increased to the value of 〈k 〉 = 〈k 〉 = 6.5, see the insets in Fig. 4(a) and (b) for the minimum of Δn. For details, Fig. 3 in SI shows the evolution of infected densities ρ for the cases of 〈k 〉 = 〈k 〉, confirming the result of Δn = 0 in Fig. 4(a) and (b). On the other hand, comparing the two insets in Fig. 4 (a) and (b), we find that Δn in Fig. 4(a) is greater than the corresponding one in Fig. 4(b), indicating that a larger β is in favor of suppressing the mixed outbreak pattern. These results can be also theoretically explained by the microscopic Markov-chain approach, see Methods for details. The solid lines in Fig. 4 (a) and (b) represent the theoretical results from Eqs (6) and (7). It is easy to see that the theoretical results are consistent with the numerical simulations very well.

Figure 4

(a) and (b) Dependence of n on the average degree 〈k 〉 of network for β  = 0.09 and 0.3, respectively, where the average degree of network is fixed as 〈k 〉 = 6.5, 〈k 〉 = 1.0, and the insets show the dependence of on 〈k 〉. The solid lines represent the theoretical results from Eqs (6) and (7). (c) Dependence of n on β with 〈k 〉 = 6.5 and 〈k 〉 = 1.5 where the “squares” and “circles” represent the case of 〈k 〉 = 0.6 for the networks , respectively, and the “up triangles” and “down triangles” represent the case of 〈k 〉 = 1.2 for the networks , respectively. (d) Dependence of Δn on β with 〈k 〉 = 6.5 and 〈k 〉 = 1.5 where the “squares” and “circles” represent the cases of 〈k 〉 = 0.6 and 1.2, respectively. The other parameters are set as , γ = 0.2, δ = 0.02, p 0 = 0.01, and N  = N  = 1000. All the results are averaged over 1000 independent realizations with the simulation steps t = 20000.

Figure 4(c) and (d) show the influences of the coupling parameters such as β and k for the outbreak number n and the difference of outbreak number Δn, respectively, where 〈k 〉 = 6.5 and 〈k 〉 = 1.5. From Fig. 4(c) we see that for the case of k  = 0.6, is an approximate constant and gradually increase with β . While for the case of , both and change with β , indicating that both and β take important roles in the synchronized and mixed outbreak patterns. From Fig. 4(d) we see that the case of k  = 1.2 decreases faster than the case of k  = 0.6, indicating that both the larger k and larger β are in favor of the synchronized outbreak pattern. That is, the stronger coupling will suppress the mixed outbreak pattern but enhance the synchronized outbreak pattern. On the contrary, the weaker coupling is in favor of the mixed epidemic outbreak pattern but suppress the synchronized outbreak pattern. For details, Fig. 4 in SI shows the evolution of infected densities ρ for different β , confirming the above results.

Coupling induced correlation between the epidemics of the two networks

The coupling between the two layers is represented by the pair of variables (β , 〈k 〉). Qualitatively, larger β and 〈k 〉 represent stronger coupling. To quantitatively represent the effects of β and 〈k 〉, we here measure the cross-correlation, r, defined in Eq. (19), which can show some information beyond the synchronized and mixed outbreak patterns. By Eq. (19) we first calculate the coefficient r between the two time series of GOPC and GP and find r = 0.66, indicating that these two data are highly correlated. Figure 5(a) shows the correspondence between the two time series of GOPC and GP. Then, we check the influence of coupling on the coefficient r. Figure 5(b) shows the dependence of r on β for 〈k 〉 = 0.2, 0.5, 1.0 and 2.0, respectively. We see that r increases with β for a fixed 〈k 〉 and also increases with 〈k 〉 for a fixed β , indicating the enhanced correlation by the coupling strength. Very interesting, we find that for the case of 〈k 〉 = 1.0 in Fig. 5(b), the point of r = 0.66 corresponds to 〈β 〉 = 0.09 (see the purple “star”), implying that the coupling in Fig. 5(a) is equivalent to the case of 〈k 〉 = 1.0 and 〈β 〉 = 0.09 in Fig. 5(b). In this sense, we may draw a horizontal line passing the purple “star” in Fig. 5(b) (see the dashed line) and its crossing points with all the curves there will also represent the equivalent coupling in Fig. 5(a). However, we notice that the horizontal line has no crossing point with the curve of 〈k 〉 = 0.2, indicating that there is a threshold of 〈k 〉 for the appearance of r = 0.66 in Fig. 5(a).

Figure 5

Correlation between the two coupled layers. (a) Correlation between the two time series of GOPC and GP. (b) Dependence of the correlation coefficient r on β for 〈k 〉 = 0.2, 0.5, 1.0 and 2.0, respectively, where the purple “star” and its related dashed line represent the point of r = 0.66. Other parameters are the same as in Fig. 3(a) and (b).

It is maybe helpful to understand the influence of (β , 〈k 〉) on epidemics from the aspect of purely coupling in network structure. We know that a critical problem in coupled networks is when they behave as separate networks vs when they behave as a solid single network[64, 65]. In our case here, the synchrony of epidemic in the two layers corresponds the case of strong coupling while the mixed pattern represents the case of weak coupling. Therefore, the synchrony and mixed patterns also reflect the influence of network structure on its dynamics. On the other hand, we notice from Fig. 5(b) that there is a finite value of r ≈ 0.43 when β  = 0, implying that the same infection probability β(t) in the two networks do have a contribution to the correlation r. However, the further increase of r in Fig. 5(b) reflects the influence from the coupling parameters of (β , 〈k 〉). From Fig. 5(b) we see that for a fixed β , r increases with 〈k 〉, indicating the influence of 〈k 〉. While for a fixed 〈k 〉, r increases with β , indicating the influence of β . Therefore, the further increase of correlation from r ≈ 0.43 do come from the coupling between the two time series.

Discussion

The dependence of the synchronized and mixed outbreak patterns on the average degrees of networks may be also understood from the aspect of their epidemic thresholds. By the theoretical analysis in Methods we obtain the epidemic thresholds as in Eq. (18), when the two networks are weakly coupled. For the case of 〈k 〉 > 〈k 〉, we have . When β satisfies , the epidemics cannot survive in each of networks . Thus, the infected fraction will be approximately zero, i.e. no epidemic outbreak in both and . It should be noticed that this result just holds for the case of weakly coupling. When coupling is strong, it is possible for epidemic to occur in the coupled network even when β is below the epidemic threshold of each layer[24]. When β satisfies , the epidemics will survive in both networks , indicating an outbreak will definitely occur in both of them. These two cases are trivial. However, when β is in between , some interesting results may be induced by coupling. When coupling is weak, it is possible for epidemic to outbreak only in network but not in network . When coupling is increased slightly, it will be also possible for epidemic to outbreak sometimes in network , i.e. resulting a mixed outbreak pattern. Once coupling is further increased to large enough, a synchronized outbreak pattern will be resulted.

So far, the reported results are from the ER random networks . We are wondering whether it is possible to still observe the phenomenon of the synchronized and mixed outbreak patterns in other network topologies. For this purpose, we here take the scale-free network[62] as an example. Very interesting, by repeating the above simulation process in scale-free networks we have found the similar phenomenon as in ER random networks, see Figs 5–7 in SI for details. Therefore, the synchronized and mixed outbreak patterns are a general phenomena in multi-layered epidemic networks.

In sum, the epidemic spreading has been well studied in the past decades, mainly focused on both the single and multi-layered networks. However, only a few works focused on the aspect of recurrent epidemics, including both the models in homogeneous population[61] and our recent model in a single network[57]. We here report from the real data that the epidemics from different networks are in fact not isolated but correlated, implying that they should be considered as a multi-layered network. Motivated by this discovery, we present a two-layered network model to reproduce the correlated recurrent epidemics in coupled networks. More importantly, we find that this model can reproduce both the synchronized and mixed outbreak patterns in real data. The two-layered network favors to show the synchronized pattern when the average degrees of the two coupled networks have a large difference and shows the mixed pattern when their average degrees are very close. Besides the degree difference between the two networks, the coupling strength between the two layers has also significant influence to the synchronized and mixed outbreak patterns. We show that both the larger β and larger 〈k 〉 are in favor of the synchronized pattern but suppress the mixed pattern. This finding thus shows a new way to understand the epidemics in realistic multi-layered networks. Its further studies and potential applications in controlling the recurrent epidemics may be an interesting topic in near future.

Methods

A two-layered network model of recurrent epidemic spreading

We consider a two-layered network model with coupling between its two layers, i.e. the networks in Fig. 2(a). We let the two networks have the same size N  = N  = N and their degree distributions P (k) and P (k) be different. We may imagine the network as a human contact network for one geographic region and the network for a separated region. Each node in the two-layered network has two kinds of links, i.e. intra-links within or and the interconnection between . The former consists of the degree distributions P (k) and P (k) while the latter the interconnection network. We let 〈k 〉, 〈k 〉, and 〈k 〉 represent the average degrees of networks and interconnection network , respectively. In details, we firstly generate two separated networks with the same size N and different degree distributions P (k) and P (k), respectively. Then, we add links between . That is, we randomly choose two nodes from and then connect them if they are not connected yet. Repeat this process until the steps we planned. In this way, we obtain an uncorrelated two-layered network.

To discuss epidemic spreading in the two-layered network, we use the extended SIRS model, see Fig. 2(b) for its schematic illustration. In this model, a susceptible node has three ways to be infected. The first one is the natural infection from environment or unknown reasons, represented by a small probability p 0. The second one is the infection from contacting with infected individuals in the network (or ), represented by β(t). And the third one is the infection from another network, represented by β (see Fig. 2(a)). Thus, a susceptible node will become infected with a probability where k is the infected neighbors in the same network and is the infected neighbors in another network. At the same time, an infected node will become refractory by a probability γ and a refractory node will become susceptible again by a probability δ.

In numerical simulations, the dependence of β(t) on time is implemented as follows[57]: we divide the total time t into multiple segments with length T and let T = 52, corresponding to the 52 weeks in one year. We let β(t) be a constant in each segment, which is randomly chosen from the truncated Gaussian distribution . Once a β(t) < 0 or β(t) > 1 is chosen, we discard it and then choose a new one. At the same time, we fix γ = 0.2 and δ = 0.2 and set β as the tunable parameter.

A theoretical analysis based on microscopic Markov-chain approach

Let be the probabilities for node i in network to be in one of the three states of S, I and R at time t, respectively. Similarly, we have in network . They satisfy the conservation law

According to the Markov-chain approach[43, 57, 66–68], we introduce where represent the densities of susceptible, infected, and refractory individuals at time t in network , respectively. Similarly, we have in network .

Let, be the transition probabilities from the state S to I, I to R, and R to S in network , respectively. By the Markov chain approach[66, 68] we havewhere Λ represents the neighboring set of node i in network . The term (1 − p 0) in Eq. (2) represents the probability that node i is not infected by the environment. The term is the probability that node i is not infected by the infected neighbors in network . While the term is the probability that node i is not infected by the infected neighbors in another network. Thus, is the probability for node i to be in the susceptible state. Similarly, for the node in network , we have

Based on these analysis, we formulate the following difference equations

The first term on the right-hand side of the first equation of Eq. (4) is the probability that node i is remained as susceptible state. The second term stands for the probability that node i is changed from refractory to susceptible state. Similarly, we have the same explanation for the other equations of Eqs (4) and (5). Substituting Eqs (2) and (3) into Eqs (4) and (5), we obtain the microscopic Markov dynamics as follows

Instead of getting the analytic solutions of Eqs (6) and (7), we solve them by numerical integration. We set the initial conditions as . To conveniently compare the solutions with the numerical simulations in the section Results, we use the same set of β(t) for both the integration and numerical simulations. Figure 6 shows the results where the left and right panels are for the networks and , respectively. In Fig. 6(b–d) and (f–h), the solid curves represent the theoretical solutions while the “circles” represent the numerical simulations. It is easy to see that the theoretical solutions are consistent with the numerical simulations very well.

Figure 6

Comparison between the theoretical solutions and numerical simulations. The left and right panels are for the networks , respectively. All the parameters are the same as in Fig. 3(a) and (b). (a) and (e) β(t) versus t; (b) and (f) ρ versus t; (c) and (g) ρ versus t; (d) and (h) ρ versus t. In (b–d) and (f–h), the solid curves represent the theoretical solutions while the “circles” represent the numerical simulations.

A theoretical analysis based on mean field theory

To go deeper into the mechanism of the synchronized and mixed outbreak patterns, we try another theoretical analysis on mean field equations. Let s (t), i (t) and r (t) represent the densities of susceptible, infected, and refractory individuals at time t in network , respectively. Similarly, we have s (t), i (t) and r (t) in network . Then, they satisfy

According to the mean-field theory, we have the following ordinary differential equations

Specifically, we consider a case of single epidemic outbreak with extremely weak coupling, i.e. p 0 = 0 and β  ≈ 0. In the steady state, we havewhich gives

Substituting Eq. (16) into Eq. (8) we obtain

The critical point can be given by letting i (t) and i (t) in Eq. (17) change from zero to nonzero, which givesFor the considered networks with 〈k 〉 > 〈k 〉, we have .

Cross-correlation measure

In statistics, the Pearson correlation coefficient is a measure of the linear correlation between two variables. If two time series {X } and {Y } have the mean values , we can define the correlation coefficient r as follows

To analyze the correlations of the growth trends between the two time series, we investigates their cross-correlation r(t) in a given window w [69, 70], i.e. . {X } and {Y } will share the same trend in the time interval w when r(t) > 0 and the opposite growth trend when r(t) < 0. In this work, we let the whole time series be one window, i.e. with w being the total evolutionary time t.