Suppressing disease spreading by using information diffusion on multiplex networks

Although there is always an interplay between the dynamics of information diffusion and disease spreading, the empirical research on the systemic coevolution mechanisms connecting these two spreading dynamics is still lacking. Here we investigate the coevolution mechanisms and dynamics between information and disease spreading by utilizing real data and a proposed spreading model on multiplex network. Our empirical analysis finds asymmetrical interactions between the information and disease spreading dynamics. Our results obtained from both the theoretical framework and extensive stochastic numerical simulations suggest that an information outbreak can be triggered in a communication network by its own spreading dynamics or by a disease outbreak on a contact network, but that the disease threshold is not affected by information spreading. Our key finding is that there is an optimal information transmission rate that markedly suppresses the disease spreading. We find that the time evolution of the dynamics in the proposed model qualitatively agrees with the real-world spreading processes at the optimal information transmission rate.

The coevolution dynamics on complex networks has attracted much attention in recent years, since dynamic processes, ubiquitous in the real world, are always interacting with each other [1,2].In biological spreading dynamics, two strains of the same disease spread in the same population and interact through cross immunity [3][4][5] or mutual reinforcement [6].In social spreading dynamics, individuals are surrounded by multiple items of information supplied by, e.g., Facebook, Twitter, and YouTube.These sources of information compete with each other for the limited attention-span of users, and the outcome is that only a few items of information survive and become popular [7,8].Recently scholars have become aware of the coevolution or interplay between biological and social spreading dynamics [10][11][12].When a new disease enters a population, if individuals who are aware of its potential spread take preventive measures to protect themselves [9,13] the disease spreading may be suppressed.Our investigation of the intricate interplay between information and disease spreading is a specific example of disease-behavior systems [14].
Studying the micromechanisms of a disease-behavior system can help us understand coevolution dynamics and enable us to develop ways of predicting and controlling the disease spreading [11].In this effort a number of excellent models [15][16][17] have demonstrated the existence of non-trivial phenomena that differ substantially from those when there is independent spreading dynamics [18][19][20][21][22][23][24].Researchers have demonstrated that the outbreak of a disease has a metacritical point [16] that is associated with information spreading dynamics and multiplex network topology and that information propagation is promoted by disease spreading [17].Funk et al. found that the disease threshold is altered once the information and disease evolve simultaneously [15].These mod- * tangminghan007@gmail.com els make assumptions about the coevolution mechanisms of information and disease spreading and do not demonstrate the interacting mechanisms in real-world systems.Because we do not understand the microscopic coevolution mechanisms between information and disease spreading dynamics from realworld disease-behavior systems, we do not have a systematic understanding of coevolution dynamics and do not know how to utilize information diffusion to more effectively suppress the spread of disease.
We present here a systematic investigation of the effects of interacting mechanisms on the coevolution processes of information and disease spreading dynamics.We first demonstrate the existence of asymmetrical interactions between the two dynamics by using real-world data from information and disease systems to analyze the coevolution.We then propose an asymmetric spreading dynamic model on multiplex networks to mimic the coupled spreading dynamics, which will allow us to understand the coevolution mechanics.The results, obtained from both the theoretical analyses and extensive simulations, suggest some interesting phenomena: the information outbreak can be triggered by its own spreading dynamics or the disease outbreak, while the disease threshold is not affected by the information spreading.Our most important finding is that there is an optimal information transmission rate at which the outbreak size of the disease reaches its minimum value, and the time evolution of the dynamics in the proposed model qualitatively agrees with the dynamics of real-world spreading.

Empirical analysis of real-world coevolution data.
Information about disease can be obtained in many ways, including face-to-face communication, Facebook, Twitter, and other online tools.Since the growth of the Internet, search engines have enabled anyone to obtain instantaneous information about disease.Patients seek out and analyze prescriptions using search engines in hopes of obtaining a means of rapid recovery.Healthy individuals use search engines to identify protective measures against disease to maintain their good health.
To examine the coevolution of real-world data about information and disease, we use weekly synchronously evolving data on information and disease systems associated with influenza-like illness (ILI) in the US during an approximate 200-week period from 3 January 2010 to 21 September 2013.The ILI dataset records weekly outpatient visits to medical facilities, and Google Flu Trends (GFT) dataset keeps track of week queries in Google search engine about ILI symptoms [25].The GFT is used to analyse the occurrence probability of a disease [26].For simplicity, we assume that the volume of information about the disease is proportional to the GFT volume because any individual can use the Google search engine to gain information about ILI.For a detailed description of the data see Ref. [26].
Figure 1(a) shows the real-data time series of information n G (t) and disease n D (t) indicating that macroscopically the two systems exhibit similar trends and confirming that the GFT effectively predicts disease spreading [26,27]-although some researchers have expressed skepticism [28].To identify the coevolution mechanisms operating between information and disease spreading, we further investigate the time series from a microscopic point of view.Specifically, we study their relative growth rates v G (t) of n G (t) and v D (t) of n D (t) (see definitions in Method Section).Figure 1 To conceptualize the correlations of the growth trends between the two dynamics, we analyze the cross-correlations c(t) between the time series of v G (t) and v D (t) for a given window size w l [29] using the Pearson correlation coefficient c(t) between the two time series {v When c(t) > 0, the growth rates of information and disease share the same trend in the time interval w l .When c(t) < 0, the information and disease have opposite growth trends.Figure 1(c) shows that the positive and negative c(t) are uncovered for w l = 3 and w l = 20, respectively.This may be because individuals tend to search for disease information when they are infected or when someone they know is infected, and thus a disease outbreak promotes the spread of information, i.e., the growth trends of GFT and ILI will be the same.When individuals acquire information about the disease they then take action to protect themselves, and this causes the growth trends of GFT and ILI to go in opposite directions.We thus conclude that there are asymmetric interactions between the dynamics of information and disease spreading, i.e., disease spreading promotes information spreading, but information spreading suppresses disease spreading.Figure 1 The relative growth rate vD(t) (blue dashed line) and vG(t) (red solid line) of nD(t) and nG(t) versus t, respectively.(c) Crosscorrelation c(t) between the two time series of vG(t) and vD(t) for the given window size w l = 3 (blue dashed line) and w l = 20 (red solid line).(d) The fraction of negative correlations fP (blue squares) and positive correlations fN (red circles) as a function of w l .In (a), nG(t) and nD(t) are divided their average values respectively.In (b), the circles and squares denote the relative growth rate at t = 53 and 153, respectively.f P (negative correlations f N ) increases (decreases) with the w l , since individuals taking measures are dependent on the timeliness of the information.Note therefore that asymmetric interactions can only continue over a short period of time.
Coevolution dynamics on multiplex networks.We now propose a novel model based on the coevolution mechanisms in real-world data, i.e., the asymmetric interactions between information and disease spreading.Information spreads through communication networks and disease usually spreads through contact networks.Communication and contact networks usually have different topologies.To describe the distinct transmission topologies of the information and disease we use a multiplex network [30][31][32][33] and construct an artificial communication-contact coupled network without degreedegree correlations in intralayers and interlayers.
We generate uncorrelated two-layer networks A and B with degree distributions P A (k A ) and P B (k B ), where networks A and B represent the communication and contact networks, re- spectively.Nodes are individuals and edges are the interactions among individuals.Each node on layer A is randomly matched one-to-one with a node of layer B. A schematic of the communication-contact coupled networks is shown in Fig. 2(a).
Using the analysis results from real-world data, we construct an asymmetric coevolution information and disease spreading model.In the communication network (layer A) we use the classic susceptible-infected-recovered (SIR) epidemiological model [21,34,35] to describe the spreading of information about the disease.Each node can be in one of three states: susceptible, informed, or recovered.A susceptible individual has not acquired any information about the disease, infected (or informed) individuals are aware of the disease and can transmit their information to their neighbors on the communication layer, and recovered individuals have the information but do not transmit it to their neighbors.At each time step, each informed node transmits their information to each susceptible neighbor on layer A with a probability β A .The informed node recovers with a probability γ A .To include the interacting mechanism between information and disease revealed in the real-world data analysis, i.e., that disease spreading promotes the information spreading, we assume that a susceptible node will become informed when its counterpart in layer B is infected, as shown in Fig. 2(b).
We now introduce a vaccination (V) state into the disease spreading dynamics on the contact network (layer B) and the model becomes SIRV [36,37].The SIR component of the spreading dynamics is the same as the information spreading on layer A and differs only in the infection and recovery rates, β B and γ B , respectively.To introduce the mechanism from our real-world data analysis, i.e., that the spread of information suppresses disease spreading, we assume that an intelligent susceptible individual on layer B is vaccinated with probability p (i) when its counterpart node on layer A is informed and (ii) when the number of its neighbors in the infected state is equal to or greater than a static threshold φ [see Fig. 2(c)].Since immunization is always expensive, condition (i) means that the individual must use the communication network to determine the perniciousness of the disease and condition (ii) means that the individual will adopt immunization measures only when the probability of infection is sufficiently high.
We initiate asymmetrical coupled coevolution dynamics by randomly infecting a tiny fraction of seed nodes on layer B and allowing their counterparts on layer A to become informed.We set the effective information transmission and disease transmission rates to be λ A = β A /γ A and λ B = β B /γ B , respectively.Without lack of generality we set γ A = γ B = 1.
A steady state will be reached when there are no more nodes in the informed or infected state.
Heterogeneous Mean-field theory.To quantify the asymmetrical coevolution dynamics, we develop a heterogeneous mean-field theory.The outbreak threshold and the fraction of infected or informed nodes in the final state are the two quantities that control the outcome.For the information spreading, the densities of susceptible, informed, and recovered nodes with degree k A at time t are denoted by s A kA (t), ρ A kA (t), and r A kA (t), respectively.Analogously, for the disease spreading, the densities of the susceptible, infected, recovered, and vaccinated nodes with degree k B at time t are denoted by s B kB (t), ρ A kB (t), r B kB (t), and v B kB (t), respectively.We first study the time evolution of information spreading on a communication network, i.e., layer A. The evolution equation of the susceptible node with degree k A on layer A can be written where k B is the average degree of layer B, and is the probability that a susceptible node connects to an informed neighbor on uncorrelated layer A (B) (see details in the Supporting Information).The increase in ρ A kA (t) is equal to the decrease in s A kA (t), and thus the evolution equations for ρ A kA (t) and r A kA (t) are and respectively.
We next investigate the evolution of the disease spreading on layer B, the contact network.The time evolution equations for the susceptible, infected, recovered, and vaccinated nodes on layer B are and respectively, where Ψ(k B , t) is the probability that a susceptible node on layer B with degree k B will be vaccinated.More details about the Eqs.( 1)-( 7) can be found in the Supporting Information.
We describe the asymmetrical coevolution dynamics of information and disease spreading using Eqs.( 1)-( 3) and ( 4)- (7), which allow us to obtain the density of each distinct state on layer A and B at time t, i.e., where h ∈ {A, B} and χ ∈ {S, I, R, V }.When t → ∞, in the steady state, the final sizes of information and disease systems are R A and R B , respectively.Initially only a tiny fraction of nodes on layers A and B are informed or infected, and most are susceptible.Thus we have s A kA ≈ 1, s B kB ≈ 1. Linearizing Eqs. ( 2) and ( 5), i.e., neglecting the high order of ρ A kA and ρ B kB , the critical effective information transmission probability is where Λ 1 C is the maximal eigenvalue of matrix from which we obtain where Λ 1 A and Λ 1 B are the maximal eigenvalues of the adjacent matrix of layers A and B, respectively.More details can be found in the Supporting Information.The critical value λ A c separates information spreading dynamics into local and global information regions.When λ A ≤ λ A c , it is in the local information region.When λ A > λ A c , it is in the global information region.In Eq. ( 9) the global information outbreak condition is correlated only with the topologies of layers A and B, i.e., the immunization probability p and threshold φ do not affect the outbreak of information, but increasing the degree heterogeneity of layers A and B increases the information outbreak probability.
When λ > λ A c , immunization can suppress disease spreading on subnetwork B, and thus here immunization process and disease spreading can be treated as competing processes [3].Reference [3] demonstrates that the two competing processes can be treated as one after the other in the thermodynamic limit.When the immunization process spreads more quickly than the disease, it first spreads on layer B and then the disease spreads on the residual network (i.e., the network after immunization).When the disease spreads more quickly than the immunization, the opposite occurs.Using Refs.[3,17] we find that the immunization progresses more quickly than the disease, i.e., λ A λ Bu > λ B λ Au , in which , which are the thresholds for the SIR model on a one-layer network [21], and • • • are the moments of the degree distribution.Because in many real-world scenarios information spreads more quickly than disease, we focus on that case.Thus immunization and disease spreading on layer B can be treated successively and separately.When φ = 0, the approximate disease threshold is which is the same as in Ref. [17].In Eq. ( 11), where V B = pQ A , and Q A is the final density of the informed population without disease spreading obtained using link percolation theory [21].From Eq. ( 11) we can see that, as expected, the threshold is bigger than in the SIR model without vaccination.
When φ ≥ 1 we use competing percolation theory to obtain the approximate disease threshold.The information first spreads on layer A, and then the disease spreads on layer B. Although many nodes on layer A receive the information for large values of λ A , the counterparts of those informed nodes still cannot be immunized when λ B is small.This is the case because according to the proposed model the susceptible nodes that are vaccinated must have authentication from both layers A and B. These informed nodes cannot acquire authentication from layer B when λ B is below the disease threshold.Only for large values of λ B , these informed nodes can obtain authentication simultaneously from layers A and B.Here the immunized nodes are V B ≈ 0 and thus the approximate disease threshold is which is the same as the outbreak threshold of SIR disease [21], i.e., this kind of information-based immunization strategy does not affect the disease outbreak threshold, and this differs from the existing results [16,17].The disease threshold is dependent only on the topology of layer B and is independent of the topology of layer A, the immunization probability p, and the threshold φ.The asymmetrical coevolution mechanisms presented in our model may explain why the disease threshold is not altered in some real-world situations [42][43][44].3) and ( 4)- (7).In (e), the two arrows respectively indicate the numerical disease thresholds for φ ≥ 1 and φ = 0, which are obtained by observing χ.Other dynamical parameters are set to be λB = 0.5 and p = 0.8.

Simulation results.
We perform extensive stochastic simulations to study the proposed asymmetrically interacting spreading dynamics on multiplex networks.In the simulations the network sizes and average degrees are set at N A = N B = 10 4 and k A = k B = 8, respectively.We use the uncorrelated configuration model to generate layers A and B according to the given degree distributions [45].For each multiplex network, we perform the dynamics 10 4 times and measure the average final fraction of information size R A , disease size R B , and immunization size V B with five randomly selected seeds in layer B. We then average these results over 100 network realizations.
To understand the coevolution dynamics of information and disease, we use Erdős-Rényi (ER) networks to represent the communication and contact networks.The degree distributions of layer A and layer B are Figure 3 shows how the immunization threshold φ affects the final information, disease, and vaccination sizes.For the information spreading on layer A, we find that R A increases with λ A and λ B [see Figs.3(a) and (d)].In addition, R A increases with φ because the individuals in layer B need a large φ value to guide their immunization decisions [see Figs.3(c) and (f)], which causes R B to increase with φ [see Figs.3(b)  and (e)].As a result, the information spreading increases as disease spreading increases.
Figures 3(b) and (e) show that R B increases with φ, since individuals are increasingly reluctant to be immunized as φ increases, and this causes V B to decrease with φ [see Figs.3(c) and (f)].Note that R B and V B as a function of λ A have a nonmonotonic shape for φ = 2 and 4, that R B (V B ) first decreases (increases) with λ A and then increases (decreases) with λ A .Thus there is an optimal information transmission rate λ O A at which R B (V B ) reaches its minimum (maximum) value.Qualitatively this is because a node on layer B will be immunized only (i) when its counterpart on layer A is informed, and (ii) when the number of its infected neighbors n B I is larger than φ.For a given λ B , condition (i) is difficult to fulfill when λ A is small and the spread of the information is slow.Increasing λ A allows more nodes to fulfill condition (i) and allows V B (R B ) to increase (decrease) with λ A .When the value of λ A is very large the information spreads so rapidly that condition (ii) can no longer be satisfied.Thus V B decreases with λ A , which enhances the spread of disease.The optimal phenomenon is not qualitatively affected by the recovery rates of information and disease.As shown in Fig. 3(e), R B versus λ B displays a non-monotonic shape for φ = 2 and 4, i.e., R B first increases with λ B and then decreases.When λ A = 0.5 the information spreading is rapid.Increasing λ B allows more nodes to fulfill the second immunization condition and to be immunized [see Fig. 3(f)], and further leads to the decrease (φ = 2) or saturation (φ = 4) of R B with λ B .The theoretical predictions of our heterogeneous mean-field theory agree with the simulation predictions.The differences between the theoretical predictions and the simulations are caused by the dynamic correlations among the states of the neighbors and by finite-size network effects [17].The dynamic correlations are produced when the information (disease) transmission events to one node in layer A (B) coming from two distinct neighbors are correlated [41].In the case of coevolution dynamics, the dynamic correlations are also induced by the counterparts of susceptible nodes [4].
For the disease spreading on layer B, the disease threshold λ B c for φ = 0 is clearly larger than the threshold λ B c0 = 1/ k B , which is the disease threshold without immunization (i.e., p = 0) [see the right arrow in Fig. 3(e)].We can determine the numerical disease threshold by measuring the susceptibility [39] or variability [40] (see details in Method).Note that the disease threshold λ B c for φ ≥ 1 is the same as λ B c0 , which is consistent with the theoretical prediction [see Eq. ( 12) and the left arrow in Fig. 3(e)].This occurs because individuals choose immunization only when the number of their infected neighbors is equal to or greater than φ.The asymmetrical coevolution mechanisms proposed in our model may explain why choosing to be immunized during disease spreading does not affect the disease threshold [42][43][44].
We use φ = 2 to measure the final information and disease sizes (see Fig. 4).According to Eq. ( 12), the disease threshold is λ B c = 1/ k B = 0.125.When λ B = 0.2, 0.5, and 0.8, any value of λ A can cause an information outbreak due to an outbreak of disease on layer B [see Fig. 4(a)].Thus the information outbreak threshold λ A c is zero.

Figures 4(b)-(c) show the optimal information transmission rate λ O
A at which R B (V B ) reaches its minimum (maximum) value.When λ A = 0.2, 0.5, and 0.8, R A increases with λ B because of the increase in the disease [see Fig. 4(d)].Note that λ B c is not affected by λ A [see the arrow in Fig. 4(e)].As shown in Fig. 4(e), R B versus λ B first increases and then decreases for large λ A = 0.5 and 0.8.This phenomenon can be understood in the same way with Fig. 3(e).There is again good agreement between the theoretical and numerical results.
Figure 5 shows the effects of λ A and λ B on the final steady state for R A , R B , and V B for φ = 2 and shows the phase diagrams for the final sizes as a function of λ A and λ B .A is slightly smaller than λ B because whether information induces an individual to be vaccinated depends on the infection level of their neighbors.Our heterogeneous mean-field theory describes this phenomenon very well.
Thus we know that for a given disease transmission rate there is an optimal information transmission rate at which the disease spreading is markedly reduced.In order to determine the coevolution characteristics of information and disease spreading when the information reaches its optimal transmission, we first look at the macroscopic coevolution of the two dynamics under different information transmission rates as shown in Fig. 6.We denote the fraction of nodes on layer A informed by their neighbors or by their counterpart nodes using ρ A A (t) and ρ B A (t), respectively.Here ρ A (t) [ρ B (t)] is the fraction of nodes obtaining the information (disease) on layer A (B) at time t.For small λ A = 0.13 below λ O A [see Fig. 6 neously.Note that ρ B (t) is larger than ρ A A (t) and very close to ρ B A (t), which means that the spread of information is primarily induced by the disease outbreak.At λ O A = 0.22, we find that ρ A A (t), ρ B A (t), and ρ B (t) reach their peaks simultaneously, and that ρ B (t) is closer to ρ A A (t) than to ρ B A (t). Thus the information and disease have a similar spreading velocity.For a large value of λ A = 0.4, the information spreads more quickly than the disease.Our results suggest that information and disease spreading have a similar macroscopic coevolution characteristic when the information transmission rate is at its optimal value.
Figure 7 shows the microscopic coevolution characteristics of the two dynamics at the optimal information transmission rate.Figure 7(a) shows the time evolution of information and disease in three independent dynamical realizations that have similar trends in their macroscopic coevolution of information spreading and disease spreading.Figure 7(b) shows the relative growth rates of information v I (t) and disease v D (t).As in the real-world case in Fig. 1(b), the same and opposite growth trends are observed.Figure 7(c) shows the calculated cross-correlations between the two time series of v D (t) and v I (t).Both positive and negative cross-correlations exist when the window size is small [see Fig. 7(d)].Note that Fig. 7 agrees well with the real-world situation shown in Fig. 1.Through extensive simulations, we find that heterogeneous networks display a similar phenomenon.Thus the coevolution between information and disease can become optimal in which the macroscopic and microscopic coevolution characteristics of information and disease exhibit similar trends and the information diffusion greatly suppresses the spread of disease.To examine how topology affects multiplex systems, we next simulate different possible heterogeneities in the communication and contact networks (see Fig. 8).We generate scale-free (SF) networks with a power-law degree distribution  P (k) ∼ k −γD by using an uncorrelated configuration model [45,46] in which γ D is the degree exponent.Through extensive simulations we find that the values of γ D do not qualitatively affect the results.Without loss of generality we set γ D = 3.0.Note that there is an optimal information transmission rate at which the disease is significantly suppressed [see Figs.8(b)-(c)], and thus heterogeneity in network topology does not qualitatively affect this optimal phenomenon.We also find that the multiplex networks with a homogeneous communication layer and a heterogeneous contact layer have a greater optimal information transmission rate.As the information (disease) spreads more (less) widely on homogeneous (heterogeneous) networks for a large transmission rate, R B is further reduced.Figure 8(e) shows that the disease threshold λ B c is determined only by the topology of layer B, and that the topology of layer A does not affect λ B c .For information spreading on layer A as shown in Fig. 8(a), R A decreases with the degree heterogeneity of layer B, since a homogeneous contact network facilitates the spread of disease for large λ B = 0.5 [20].In Figs.c the heterogeneity of layer A does not increase information diffusion, but promotes disease spreading because nodes are less likely to be immunized.We examine the effects of the heterogeneity of layer B and find that R A and R B increase (decrease) with the degree heterogeneity of layer B for small (large) λ B .When the degree heterogeneity of layer B is increased, the network has a large number of individuals with very small degrees and more individuals with large degrees.When λ B is small there are more hubs in heterogeneous networks that facilitate disease spreading because they are more likely to be infected, and this increases information diffusion.When λ B is large, however, there are many smalldegree nodes with a low probability of being infected, and this produces smaller values of R B , which causes smaller values of R A .

DISCUSSION
We have systematically investigated the coevolution dynamics of information and disease spreading on multiplex networks.We first discover indications of asymmetrical interactions between the two spreading dynamics by analyzing real data, i.e., the weekly time series of information spreading and disease spreading in the form of influenza-like illness (ILI) evolving simultaneously in the US during an approximate 200-week period from 3 January 2010 to 10 December 2013.Using these interacting mechanisms observed in real data, we propose a mathematical model for describing the coevolution spreading dynamics of information and disease on multiplex networks.We investigate the coupled dynamics using heterogeneous mean-field theory and stochastic simulations.We find that information outbreaks can be triggered by the spreading dynamics within a communications network and also by disease outbreaks in the disease contact network, but we also find that the disease threshold is not affected by information spreading, i.e., that the outbreak of disease is solely dependent on the topology of the contact network.More important, for a given rate of disease transmission we find that there is an optimal information transmission rate that decreases the disease size to a minimum value, and the modeled evolution of information and disease spreading is consistent with realworld behavior.We also verify that heterogeneity in network topology does not invalidate the results.In addition, we find that when information diffuses slowly, the degree heterogeneity of the communication network has a trivial impact on disease spreading.The homogeneity of the communication network can enhance the vaccination size and thus prevent disease spreading more effectively when the spread of information is rapid.
The asymmetrical interacting mechanism we discover by analyzing real-world data provides solid evidence supporting the basic assumptions of previous researches [16,17].Our data-driven model also reveals some fundamental coevolution mechanisms in the coevolution dynamics.Using these coevolution dynamics of information and disease we are able to identify phenomena that differ qualitatively from those found in previous research on disease-behavior systems.Our results enable us to quantify the optimal level of information transmission that suppresses disease spreading.The coevolution mechanisms also enable us to better understand why the disease threshold is unchanged even when information spreading in some real-world situations undergoes coevolution.
Further research on disease-behavior systems promises to discover additional real-world mechanisms that can be used to refine models of coevolution spreading dynamics.Developing a more accurate theoretical method is full of challenges because it is difficult to describe the strong dynamic correlations among the states of neighboring nodes in a network.If we take dynamical correlations into account, we may be able to use such advanced theoretical methods as dynamic messagepassing [47,48] or pair approximation [49,50].

Relative growth rates. We define the relative growth rates
and where R h is the final information size R A or disease size R B , and • • • is the ensemble averaging.The value of χ exhibits a peak at the critical point at which the thresholds can be computed.3) and ( 4)- (7).In (e), the two arrows respectively indicate the numerical disease thresholds for φ ≥ 1 and φ = 0, which are obtained by observing χ.Other dynamical parameters are set to be λ B = 0.5 and p = 0.8.
(b) shows the evolution of v G (t) and v D (t).Note that the same and opposite growth trends of v G (t) and v D (t) coexist.For example, at week 53 (week 153), v G (53) > 0 [v G (153) > 0] and v D (53) < 0 [v D (153) > 0].Thus the GFT and ILI show the opposite (the same) growth trends.

FIG. 1 .
FIG. 1. (Color online) Weekly outpatient visits and Google Flu Trends (GFT) of influenza-like illness (ILI) from 3 January 2010 to and 21 September 2013 in the United States.(a) The relative number of outpatient visits nD(t)/ nD(t) (blue dashed line) and relative search queries aggregated in GFT nG(t)/ nG(t) (red solid line) versus t, where nD(t) = tmax t=1 nD(t)/tmax and nG(t) = tmax t=1 nG(t)/tmax, and tmax is the number of weeks.(b) The relative growth rate vD(t) (blue dashed line) and vG(t) (red solid line) of nD(t) and nG(t) versus t, respectively.(c) Crosscorrelation c(t) between the two time series of vG(t) and vD(t) for the given window size w l = 3 (blue dashed line) and w l = 20 (red solid line).(d) The fraction of negative correlations fP (blue squares) and positive correlations fN (red circles) as a function of w l .In (a), nG(t) and nD(t) are divided their average values respectively.In (b), the circles and squares denote the relative growth rate at t = 53 and 153, respectively.

FIG. 2 .
FIG. 2. (Color online) Illustration of asymmetrical mechanisms of information and disease on multiplex networks.(a) A multiplex network is used to represent communication and contact networks, which are denoted as layer A and layer B, respectively.Each layer has 5 nodes.(b) The promotion of information spreading by disease.If node 5 on layer B is infected, its counterpart on layer A becomes informed.(c) The suppression of disease spreading by information diffusion.Node 3 in layer B becomes vaccination only when: (1) its counterpart on layer A is in the informed state and (2) the number of its infected neighbors on layer B is equal to the threshold φ = 2.

FIG. 3 .
FIG. 3. (Color online) With immunization thresholds φ being the parameter of interest, the final sizes of information, disease and vaccination on two layer ER-ER multiplex networks.(a) The final information size RA, (b) the final disease size RB, and (c) the final vaccination size VB versus information transmission rate λA for different values of immunization threshold φ with λB = 0.5.For different values of φ, (d) RA, (e) RB and (f) VB as a function of λB at λA = 0.5.The symbols represent the simulation results and the lines are the theoretical predictions obtained by numerically solving Eqs.(1)-(3) and (4)-(7).In (e), the two arrows respectively indicate the numerical disease thresholds for φ ≥ 1 and φ = 0, which are obtained by observing χ.Other dynamical parameters are set to be λB = 0.5 and p = 0.8.

FIG. 4 .
FIG. 4. (Color online) With disease transmission rate λB being the parameter of interest, the asymmetrically interacting dynamics spreads on ER-ER networks.(a) The final information size RA, (b) the final disease size RB, and (c) the vaccination size VB versus the information transmission rate λA for the disease transmission rate λB = 0.2, 0.5 and 0.8.For λA = 0.2, 0.5 and 0.8, (d) RA, (e) RB and (f) VB as a function of λB.In the figures, symbols are the simulation results and the lines are the theoretical predictions.In (e), the arrow indicates the numerical disease threshold.We set other parameters to be φ = 2 and p = 0.8.
FIG. 5. (Color online) Asymmetrically interacting dynamics on ER-ER networks.The final density in each state relating the parameters λA and λB: (a) the final information size RA, (b) the final disease size RB and (c) the vaccination size VB.In (a), the horizontal and vertical dashed lines separate the λA − λB plane into local and global information outbreak regions, which are denoted as regions I and II.In (b), the vertical dashed line divides the plane into a local (region I) and a global (region II) disease outbreak regions.In (b), the blue circles (λA = 0.13, λB = 0.3), green up triangle (λA = 0.22, λB = 0.3) and gray diamond (λA = 0.4, λB = 0.3) represent λA being below, at and above λ OA , respectively (see more discussions in Fig.6).The black squares (black lines) in (b) and (c) represent the optimal information transmission rate λ O A versus λB.Other parameters are set to be φ = 2 and p = 0.8.

Figure 5 (
b) shows how region I and region II are separated by λ B c (see vertical white dashed line).For the minimum value of R B in region II, λ O A increases linearly with λ B , as shown in Fig. 5(b) [see black lines and symbols in (b) and (c)].At the optimal λ O A , R B (V B ) reaches its minimum (maximum) value, as shown in Fig. 5(b) [Fig.5(c)].Note that λ O

FIG. 7 .
FIG. 7. (Color online) Asymmetrically interacting spreading dynamics on coupled ER-ER networks at the optimal information transmission rate.(a)The fractions of nodes in the informed state ρA(t) (red solid line) and infected state ρB(t) (blue dashed line) versus t.(b) The relative growth rates vD(t) (blue dashed line) and vI (t) (red solid line) of ρB(t) and ρA(t) versus t, respectively.(c) Cross-correlations c(t) between vI (t) and vD(t) for the given window size w l = 3 (blue dashed line) and w l = 5 (red solid line).(d) The fractions of negative correlations fP (blue squares) and positive correlations fN (red circles) as a function of w l .We set other parameters to be λA = 0.22, λB = 0.3 and p = 0.8, respectively.

FIG. 8 .
FIG. 8. (Color online) Effect of degree heterogeneity on coevolution dynamics.(a) The final information size RA, (b) the final disease size RB and (c) the vaccination size VB versus the information transmission rate λA on ER-ER, ER-SF, SF-ER and SF-SF coupled networks with λB = 0.5.For ER-ER, ER-SF, SF-ER and SF-SF networks with λA = 0.5, (d) RA, (e) RB and (f) VB as a function of λB.Other parameters are set to be φ = 2, p = 0.8 and kA = kB = 8.
FIGURE LEGENDS tmax t=1 n D (t)/t max and n G (t) =

2 .Fig. 3 .
tmax t=1 n G (t)/t max , and t max is the number of weeks.(b) The relative growth rate v D (t) (blue dashed line) and v G (t) (red solid line) of n D (t) and n G (t) versus t, respectively.(c) Cross-correlation c(t) between the two time series of v G (t) and v D (t) for the given window size w l = 3 (blue dashed line) and w l = 20 (red solid line).(d) The fraction of negative correlations f P (blue squares) and positive correlations f N (red circles) as a function of w l .In (a), n G (t) and n D (t) are divided their average values respectively.In (b), the circles and squares denote the relative growth rate at t = 53 and 153, respectively.FIG.Illustration of asymmetrical mechanisms of information and disease on multiplex networks.(a) A multiplex network is used to represent communication and contact networks, which are denoted as layer A and layer B, respectively.Each layer has 5 nodes.(b) The promotion of information spreading by disease.If node 5 on layer B is infected, its counterpart on layer A becomes informed.(c) The suppression of disease spreading by information diffusion.Node 3 in layer B becomes vaccination only when: (1) its counterpart on layer A is in the informed state and (2) the number of its infected neighbors on layer B is equal to the threshold φ = 2.With immunization thresholds φ being the parameter of interest, the final sizes of information, disease and vaccination on two layer ER-ER multiplex networks.(a) The final information size R A , (b) the final disease size R B , and (c) the final vaccination size V B versus information transmission rate λ A for different values of immunization threshold φ with λ B = 0.5.For different values of φ, (d) R A , (e) R B and (f) V B as a function of λ B at λ A = 0.5.The symbols represent the simulation results and the lines are the theoretical predictions obtained by numerically solving Eqs.(1)-(

FIG. 4 .
With disease transmission rate λ B being the parameter of interest, the asymmetrically interacting dynamics spreads on ER-ER networks.(a) The final information size R A , (b) the final disease size R B , and (c) the vaccination size V B versus the information transmission rate λ A for the disease transmission rate λ B = 0.2, 0.5 and 0.8.For λ A = 0.2, 0.5 and 0.8, (d) R A , (e) R B and (f) V B as a function of λ B .In the figures, symbols are the simulation results and the lines are the theoretical predictions.In (e), the arrow indicates the numerical disease threshold.We set other parameters to be φ = 2 and p = 0.8.FIG.5.Asymmetrically interacting dynamics on ER-ER networks.The final density in each state relating the parameters λ A and λ B : (a) the final information size R A , (b) the final disease size R B and (c) the vaccination size V B .In (a), the horizontal and vertical dashed lines separate the λ A − λ B plane into local and global information outbreak regions, which are denoted as regions I and II.In (b), the vertical dashed line divides the plane into a local (region I) and a global (region II) disease outbreak regions.In (b), the blue circles (λ A = 0.13, λ B = 0.3), green up triangle (λ A = 0.22, λ B = 0.3) and gray diamond (λ A = 0.4, λ B = 0.3) represent λ A being below, at and above λ O A , respectively (see more discussions in Fig. 6).The black squares (black lines) in (b) and (c) represent the optimal information transmission rate λ O A versus λ B .Other parameters are set to be φ = 2 and p = 0.8.FIG.6.On ER-ER coupled networks, the time evolution of each type of nodes.The time evolution of ρ A A (t), ρ B A (t), ρ A (t) and ρ B (t) for (a) λ A = 0.13, (b) λ A = 0.22 and (c) λ A = 0.40.Other parameters are set to be λ B = 0.3, φ = 2 and p = 0.8.

FIG. 7 .FIG. 8 .
FIG.7.Asymmetrically interacting spreading dynamics on coupled ER-ER networks at the optimal information transmission rate.(a)The fractions of nodes in the informed state ρ A (t) (red solid line) and infected state ρ B (t) (blue dashed line) versus t.(b) The relative growth rates v D (t) (blue dashed line) and v I (t) (red solid line) of ρ B (t) and ρ A (t) versus t, respectively.(c) Cross-correlations c(t) between v I (t) and v D (t) for the given window size w l = 3 (blue dashed line) and w l = 5 (red solid line).(d) The fractions of negative correlations f P (blue squares) and positive correlations f N (red circles) as a function of w l .We set other parameters to be λ A = 0.22, λ B = 0.3 and p = 0.8, respectively.
) [n D (t)] shows an increasing trend at time t.If not, n G (t) [n D (t)] shows a decreasing trend at time t.Variability measure.The variability χ [40] is χ