This study is concerned with the dynamic behaviors of epidemic spreading in multiplex networks. A model composed of two interacting complex networks is proposed to describe cooperative spreading processes, wherein the virus spreading in one layer can penetrate into the other to promote the spreading process. The global epidemic threshold of the model is smaller than the epidemic thresholds of the corresponding isolated networks. Thus, global epidemic onset arises in the interacting networks even though an epidemic onset does not arise in each isolated network. Simulations verify the analysis results and indicate that cooperative spreading processes in multiplex networks enhance the final infection fraction.

In recent years, epidemic spreading has received considerable attention in multiplex networks. The main focus of studies on this process is on how epidemic spreads in multiplex networks when two spreading processes interact. Epidemic threshold and final infection fraction are the two important parameters underlying any spreading processes, to contrast the two parameters between the interacting networks, and the corresponding isolated networks has practical significance. In this study, a model for cooperative spreading processes on interacting two-layer networks is proposed. This study aims to prove that the epidemic thresholds of interacting two-layer networks can be decreased for cooperative spreading processes, which implies that cooperative spreading processes promote the spread of disease. This study determines that the global epidemic onset arises. Thus, a global epidemic threshold can be used to uniform the two epidemic thresholds. This model can also be used for unidirectional coupled networks. It is revealed that a network with a small epidemic threshold confirms the global epidemic threshold of the multiplex networks.

Epidemic spreading, which is an important dynamic process in complex networks, has attracted great attention for a long time. Many studies initially focused on isolated networks of fixed topology with mean-field approximation.1,2 The spread of sexually transmitted diseases and computer viruses was studied based on scale-free networks.3–6 Several studies also revealed that epidemic processes in scale-free networks do not pose an epidemic threshold because of large connectivity fluctuations in infinite scale-free networks.7–9 Epidemic spreading was investigated in complex random networks with degree correlations.10 Then, reaction-diffusion model in which the nodes contain agents11–14 was used in the process of spreading. The agents can move from a node to another, and spreading processes occurred within agents on the same node. Further, the studies focus on the mobility of individuals.15,16 The movement of individuals among dense crowd in a city has a great effect on the epidemic threshold. A random diffusion model was also used to investigate epidemic spreading among objective traveling people.17,18 The railway leads to a fast increase in the number of infected agents and a high final infection fraction. These models and results are very useful for public health authorities to make effective decisions.

However, many spreading processes have been studied independently in a single network. In the real-world, many spreading processes can occur through multiple routes simultaneously. For example, sexually transmitted diseases can spread both in homosexual and heterosexual networks.19 Avian influenza, such as 2009 H1N1 and 2013 H7N9, has recently made dead-end jumps from poultry to human.20 Thus, these viruses make challenges to human progress and survival. A natural extension is to use the interacting network model in which epidemic can spread from one network to another. Funk and Jansen21 used the bond percolation analysis of two competitive viruses in two-layer networks, and the effects of layer overlapping were investigated. A study on U.S. bisexual men showed that the bisexual men were the medium to connect the pathogens from the network of heterosexual men with the network of homosexual men.22,23 Dickison et al.24 studied the susceptible-infected-refractory (SIR) process and revealed that in strongly coupled networks, epidemics spread from one layer to the other at a critical infection strength βc, below which the disease does not spread. Saumell-Mendiola et al.25 analyzed an epidemic spreading process of two interacting networks. They also developed a heterogeneous mean-field approach and revealed that a global endemic state may arise in the coupled system even though two networks are well below their respective epidemic thresholds. Granell et al.26 used a microscopic Markov chain approach to analyze the interrelation between two processes, namely, the spreading of the epidemic and the spreading of information awareness to prevent infection. They also determined that awareness spread leads to disease suppression. Zhao et al.27 used bond percolation theory and the generating function with SIR model to calculate the epidemic thresholds of multiplex networks. Sahneh and Scoglio28 extended the single-virus spread model to the two exclusive viruses in two-layer networks, and they found analytical expressions that determine extinction, coexistence, and absolute dominance of the viruses. Zhao et al.29 promoted the concept of state-dependent infected rate, and all the different influences between the two epidemics were displayed. Therefore, epidemic spreading processes on top of multiplex networks reveal a rich phase diagram of intertwined effects.

Any spreading process has two underlying fundamental parameters: epidemic threshold and final infection fraction. The key issue is the difference in the two parameters between the interacting networks and the corresponding isolated networks. Cozzo et al.30 used perturbation theory to analyze epidemic thresholds of networks, and they revealed that the spectral radius of the whole matrix is always not less than that of each submatrix. Thus, the epidemic threshold of the whole system is always not more than the corresponding isolated networks. Wang et al.31 investigated two types of spreading dynamics: one is the spread of disease, whereas the other is the spread of information about the disease. They used simulation to determine that information spreading can effectively raise the epidemic threshold. Guo et al.32 proposed a local awareness controlled contagion spreading model in multiplex networks, and they determined the emergence of an abrupt transition of epidemic threshold with the local awareness ratio through numerical simulations.

The study on the interaction between two spread processes on two-layer networks has a great significance. In this study, a model composed of two interacting complex networks is proposed to describe the cooperative spreading processes. First, this model extends spreading process in a single network to cooperative spreading in the two-layer networks. It can also be used to analyze unidirectional or bidirectional interaction between two-layer networks. Second, perturbation analysis theory is used to prove that epidemic thresholds of interacting two-layer networks can be decreased in cooperative spreading processes. This theory reveals that the cooperative interaction between two spread processes can promote spread. Epidemic onset arises in one network and that also arises simultaneously in the other network for the interacting networks. Thus, a global epidemic threshold for interacting networks exists. Finally, this theory also displays that cooperative spreading processes in multiplex network enhance the final infection fraction.

The rest of this paper is organized as follows. A model composed of two interacting complex networks and theoretical analysis is proposed in Sec. II. Then, the numerical simulations are used to show the effects of the interactions between two spreading processes in Sec. III. Finally, some discussions and conclusions are given in Sec. IV.

Interacting two-layer networks, including network A=(aij)N×N and B=(bij)N×N, presumably have the same size N with different intra-layer connectivity. Thus, the model is a multiplex networks, as shown in Fig. 1. We focus in the discrete standard susceptible-infected-susceptible (SIS) model1 for both networks with regard to the spreading dynamics. In each subnetwork, the nodes are denoted as susceptible (S) or infected (I), and the links represent the connection along which the infection can propagate. At each time step, susceptible (S) nodes may be infected from intra-layer infected nodes and inter-layer infected nodes simultaneously. On the other hand, infected nodes recover spontaneously. After some transient time, the previous dynamics make the system into a stationary state. For interacting two-layer networks, let β1(β2) be the probability of infection between nodes in network A(B) and γ1(γ2), the probability of infection from a node in B(A) to a node in A(B), μ1(μ2) is the probability of curing for network A(B). The prevalence of infection of A and B, defined as the fraction of infected nodes at a given time, is denoted by p1,i(t) and p2,i(t) for node i, respectively. So, the time evolution processes can be written as

p1,i(t+1)=(1p1,i(t))(1q1,i(t))+(1μ1)p1,i(t)+γ1p2,i(t)(1p1,i(t)),p2,i(t+1)=(1p2,i(t))(1q2,i(t))+(1μ2)p2,i(t)+γ2p1,i(t)(1p2,i(t)),
(1)

for i{1,...,N}, where q1,i(t) and q2,i(t) are the probabilities of node i not being infected by any neighbor in A and B, respectively

q1,i(t)=j=1N(1β1aijp1,j(t)),q2,i(t)=j=1N(1β2bijp2,j(t)).
(2)

In the first line in Eq. (1), the first term on the right-hand side is the probability that node i is susceptible (1p1,i(t)) and is infected (1q1,i(t)) by at least a neighbor, the second term stands for the probability that node i is infected at time t and does not recover, and the last term takes into account the probability that node i is susceptible (1p1,i(t)) and is infected by the counterpart node in other layer, the second line in Eq. (1) has the same definition with the first line. Based on different epidemic parameters, this model can construct various network structures to exhibit rich propagation dynamic behaviors.

FIG. 1.

Example of a two-layer network. The topologies are different per layer, and each node in one layer is connected to its counterparts in the other layer. Dashed and solid lines represent interlayer and intralayer connections, respectively.

FIG. 1.

Example of a two-layer network. The topologies are different per layer, and each node in one layer is connected to its counterparts in the other layer. Dashed and solid lines represent interlayer and intralayer connections, respectively.

Close modal

In the following, a simplifying assumption should be made. The probability of infection of inter-layer is much smaller than the probability of infection of intra-layer,30 since the probability of infection of inter-layer describes spreading process from one species to another species.20 So, one obtains:

γ1β1,γ2β2.
(3)

The fixed point iteration method is used to calculate nontrivial stationary solution of Eq. (1)

p1,i=(1p1,i)(1q1,i)+(1μ1)p1,i+γ1p2,i(1p1,i),p2,i=(1p2,i)(1q2,i)+(1μ2)p2,i+γ2p1,i(1p2,i),
(4)

where

q1,i=j=1N(1β1aijp1,j)q2,i=j=1N(1β2bijp2,j).
(5)

p1,i and p2,i are the probability of infection of node i at the stationary state for A and B, respectively. Thus, the final infection fractions ρ1 and ρ2 for the two interacting networks A and B are computed as

ρi=1Nj=1Npi,j,i=1,2.
(6)

When the values of μi and γi are fixed, there are epidemic thresholds β1c and β2c for A and B, respectively, thus ρ1=0 if β1<β1c, and ρ2=0 if β2<β2c. ρ1>0 if β1>β1c, and ρ2>0 if β2>β2c. When β1β1c and β2β2c, the probabilities p1,i1 and p2,i1, so from Eq. (2), one obtains

q1,i1β1j=1Naijp1,j,q2,i1β2j=1Nbijp2,j.
(7)

Inserting Eq. (7) into Eq. (4) and neglecting second-order terms, one obtains

β1j=1Naijp1,jμ1p1,i+γ1p2,i=0,β2j=1Nbijp2,jμ2p2,i+γ2p1,i=0,
(8)

thus, we obtain the following equation:

β1AP1μ1P1+γ1P2=0,β2BP2μ2P2+γ2P1=0,
(9)

where P1=[p1,1,p1,2,...,p1,N],P2=[p2,1,p2,2,...,p2,N]. From the second line of Eq. (9), we obtain

(Bμ2β2I)P2=γ2β2P1,
(10)

when μ2β2 is not the eigenvalue of matrix B, then matrix (Bμ2β2I) is reversible, note that irreversible measure of (Bμ2β2I) is zero, since the number of eigenvalues for matrix B is limited. One obtains

P2=γ2β2(Bμ2β2I)1P1.
(11)

By the same methods, one obtains

P1=γ1β1(Aμ1β1I)1P2.
(12)

Inserting Eqs. (11) and (12) into Eq. (9), we obtain

(Aγ1γ2β1β2(Bμ2β2I)1μ1β1I)P1=0,(Bγ1γ2β1β2(Aμ1β1I)1μ2β2I)P2=0,
(13)

after symbol substitution in Eq. (13), one gets

(A¯μ1β1I)P1=0,(B¯μ2β2I)P2=0,
(14)

where

A¯=Aγ1γ2β1β2(Bμ2β2I)1,B¯=Bγ1γ2β1β2(Aμ1β1I)1.
(15)

When Eq. (14) has nonzero solution (P1>0,P2>0), if and only if μ1/β1 is the eigenvalue of matrix A¯ and μ2/β2 is the eigenvalue of matrix B¯.10,33 Looking for the onset of spread, the lowest values of β1 and β2 satisfying Eq. (14) are written as

β1c=μ1Λmax(A¯),β2c=μ2Λmax(B¯),
(16)

where Λmax is the largest eigenvalue of the matrix.

Considering the two isolated networks A and B, probabilities of infection are denoted by p1,i(t) and p2,i(t) for node i, respectively

p1,i(t+1)=(1p1,i(t))(1q1,i(t))+(1μ1)p1,i(t),p2,i(t+1)=(1p2,i(t))(1q2,i(t))+(1μ1)p2,i(t).
(17)

Inserting Eq. (2) into Eq. (17) and neglecting second-order terms, we can calculate nontrivial stationary solution for isolated networks by fixed point iteration method, one obtains

(Aμ1β1I)P1*=0,(Bμ2β2I)P2*=0.
(18)

P1* and P2* are the probability of infection vectors for all nodes at the stationary state for A and B, respectively. Looking for the onset of spread (P1*>0,P2*>0), the lowest values of infected rate β1 and β2 for isolated networks satisfying Eq. (18) are

β1*=μ1Λmax(A)β2*=μ2Λmax(B).
(19)

Generally, probability of infection between layers is smaller than the probability of infection in the layer. Following the assumptions in Eq. (3), we obtain γ1γ2β1β2. By comparing Eq. (13) with Eq. (18), we can set matrix γ1γ2β1β2(Bμ2β2I)1 as a disturbance of the matrix A, and γ1γ2β1β2(Aμ1β1I)1 as a disturbance of the matrix B. Therefore, in order to contrast the epidemic thresholds β1c and β2c of interacting networks with β1* and β2* of the corresponding isolated networks, we can use perturbation analysis method to analyze the thresholds of isolated networks. The perturbed solutions to both infected thresholds β1* and β2* and infection rates P1* and P2* of isolated networks are proposed as

β1c=β1*+ϵ1β1*+O(ϵ12)β2c=β2*+ϵ2β2*+O(ϵ22)P1=P1*+ϵ3P1*+O(ϵ32)P2=P2*+ϵ4P2*+O(ϵ42).
(20)

Inserting Eq. (20) into Eq. (9), using Eq. (19) and neglecting second-order terms, Eq. (9) after some algebra yields

ϵ1β1*AP1*+γ1(1+ϵ4)P2*=0,ϵ2β2*BP2*+γ2(1+ϵ3)P1*=0,
(21)

since |ϵ3|1 and |ϵ4|1, A and B are adjacency matrixes, the results show that ϵ1<0 and ϵ2<0, one obtains β1c<β1* and β2c<β2*. Therefore, we can conclude that the epidemic threshold of the interacting two-layer networks can be decreased for two cooperative spreading processes, which means that cooperative spreading promotes the spreading processes. To obtain more accurate results, we consider the second order approximation of Eq. (5) as follows:

q1,i1β1j=1Naijp1,j+β12j<laijailp1,jp1,l,q2,i1β2j=1Nbijp2,j+β22j<lbijbilp2,jp2,l.
(22)

The second-order corresponds to reinfections and multiple infections.

The conclusion can be extended from duplex networks to multiplex networks. Take a multiplex network composed of three layers, for example. First, we consider a network which contains two layers of sub-networks. Based on the conclusion, we obtain that the epidemic threshold of this duplex network is decreased compared with the two corresponding isolated networks. Second, we take the duplex network as the first layer, and a third sub-network as the second layer. Similarly, based on the conclusion, we can conclude that the epidemic threshold of the multiplex network composed of three layers is decreased compared with the three corresponding isolated networks. Analogously, the conclusion can be generalized to multiplex networks composed of more layers of sub-networks.

It is observed that a global endemic activity arises in the two interacting complex networks, the reason being that the state ρ10 and ρ2=0 is not a fixed point of the dynamics shown in Eq. (9). This conclusion is proved as below. The state ρ10 and ρ2=0 is presumably a fixed point. Thus, the probability is P1>0, and P2=0. When P2=0, one obtains P1=0 after substitution into Eq. (9), wherein contradiction occurs. The state ρ1=0 and ρ20 is also not a fixed point, as can be revealed with the same method. Thus, the fixed point is at the state ρ10 and ρ20. Thus, when β2β2c, then β1β1c simultaneously. This finding indicates that if epidemic activity arises in one network, it will also arise in the other coupled network. That is to say, if the epidemic spreads in one of the coupled networks, it will spread to the whole system.

Generalized outer synchronization in unidirectional interaction of two networks has attracted great attention for a long time.34 However, the unidirectional interaction of two networks for an epidemic has not been considered thus far. The proposed model can also be used to analyze the unidirectional interaction of two networks by setting γ1=0 or γ2=0. In the unidirectional interaction of two-layer networks with γ1=0, network A affects network B, but the inverse is not true. When the epidemic threshold of isolated network A is smaller than that of isolated network B, we call network A as the dominant network and network B as the nondominant network. The epidemic threshold of nondominant network B is proved smaller than that of corresponding isolated network with perturbation theory. Based on Eq. (9) with γ1=0, the following can be obtained:

β1AP1μ1P1=0,β2BP2μ2P2+γ2P1=0.
(23)

The perturbed solution is proposed to the infected threshold β2* and the infection rate P2* of the isolated network B

β2c=β2*+ϵ2β2*+O(ϵ22),P2=P2*+ϵ4P2*+O(ϵ42).
(24)

The following is obtained by inserting Eq. (24) into Eq. (23), the second line in Eq. (19), and by neglecting the second-order terms:

ϵ2β2*BP2*+γ2P1*=0.
(25)

Given that B is adjacency matrix, the results show that ϵ2<0. Then, the following is obtained:

β2c<β2*.
(26)

This finding reveals that the epidemic threshold of B is well below that of the corresponding isolated network. The dominant network A with a small epidemic threshold confirms the global epidemic threshold of the multiplex networks, given that it induces a shift of the epidemic threshold of the network B to small values. In other words, the multiplex nature of the system leads to an earlier endemic activity also in the nondominant network, given that its epidemic threshold is smaller than that for the isolated networks.

In numerical simulations, two different networks formed by 1000 nodes each are generated, which follow the power law degree distribution with exponent 2.5, given that the epidemic thresholds are found at the limit for scale-free networks. Figs. 2–9 show results averaged over 20 runs of randomly generated two-layer models.

FIG. 2.

Comparison of the final infection fraction ρ1 as a function of β1 for cooperative network A with the corresponding isolated network.

FIG. 2.

Comparison of the final infection fraction ρ1 as a function of β1 for cooperative network A with the corresponding isolated network.

Close modal
FIG. 3.

Comparison of the final infection fraction ρ2 as a function of β2 for cooperative network B with the corresponding isolated network.

FIG. 3.

Comparison of the final infection fraction ρ2 as a function of β2 for cooperative network B with the corresponding isolated network.

Close modal
FIG. 4.

Comparison of the epidemic thresholds for two-layer interacting networks A and B. The x-axis represents the number of runs, and the y-axis represents the epidemic thresholds obtained by simulations for the two-layer interacting networks.

FIG. 4.

Comparison of the epidemic thresholds for two-layer interacting networks A and B. The x-axis represents the number of runs, and the y-axis represents the epidemic thresholds obtained by simulations for the two-layer interacting networks.

Close modal
FIG. 5.

Expected infection ratios ρ=(ρ1+ρ2)/2 for two-layer interacting networks as a function of probability of infection β1 and β2 using first order approximation.

FIG. 5.

Expected infection ratios ρ=(ρ1+ρ2)/2 for two-layer interacting networks as a function of probability of infection β1 and β2 using first order approximation.

Close modal
FIG. 6.

Expected infection ratios ρ=(ρ1+ρ2)/2 for two-layer interacting networks as a function of probability of infection β1 and β2 by using second order approximation.

FIG. 6.

Expected infection ratios ρ=(ρ1+ρ2)/2 for two-layer interacting networks as a function of probability of infection β1 and β2 by using second order approximation.

Close modal
FIG. 7.

Expected infection ratios ρ=(ρ1+ρ2)/2 for corresponding two isolated networks as a function of probability of infection β1 and β2.

FIG. 7.

Expected infection ratios ρ=(ρ1+ρ2)/2 for corresponding two isolated networks as a function of probability of infection β1 and β2.

Close modal
FIG. 8.

Comparison of the final infection fraction ρ1 as a function of β1 for cooperative network A with the corresponding isolated network for different values of μ1.

FIG. 8.

Comparison of the final infection fraction ρ1 as a function of β1 for cooperative network A with the corresponding isolated network for different values of μ1.

Close modal
FIG. 9.

Comparison of the final infection fraction ρ1 as a function of β1 for cooperative network A with the corresponding isolated network for different values of γ.

FIG. 9.

Comparison of the final infection fraction ρ1 as a function of β1 for cooperative network A with the corresponding isolated network for different values of γ.

Close modal

Monte Carlo simulations are employed to obtain the epidemic thresholds. The initial fraction of the infected nodes is set as 0.05, and the probability of curing are μ1=0.4 and μ2=0.4. The values for interlayer probability of infection are γ1=0.05 and γ2=0.05. For simplicity, we always set μ1=μ2 and γ1=γ2. Note that two additional types of topologies, two interacting random networks, and two interacting small-world networks are also used for simulations. The numerical results which are not shown here for brevity agree with those of the two interacting scale-free networks.

The comparison of the final infection fractions ρ1 and ρ2 of infected nodes for cooperative interaction networks with the corresponding isolated networks is shown in Figs. 2 and 3. Using first order approximation and second order approximation for cooperative spreading processes in numerical simulations, Figs. 2 and 3 show that the difference between first order approximation and second order approximation is trivial. It is also obvious that the epidemic thresholds are decreased for cooperative spreading processes compared with those of the corresponding isolated networks. This observation indicates that cooperative interacting networks promote propagation process. Figs. 2 and 3 also show that the final infection fraction for cooperative interaction networks is larger than that for isolated networks whether we use the first or second order approximation.

Further, the endemic activity of the two interacting complex networks arises simultaneously. Thus, the epidemic thresholds of the two interacting complex networks are always the same. Fig. 4 shows the comparison of epidemic thresholds for interacting networks A and B, where the x-axis represents the number of runs. In particular, red stars represent epidemic threshold β1c for network A, and black circles represent epidemic threshold β2c for network B. Fig. 4 shows that the epidemic thresholds β1c are consistent with epidemic thresholds β2c for 20 runs.

We call the sub-network with the smallest epidemic threshold as the dominant network. Numerical analysis is used to verify the conclusion indicating that spreading process in the dominant network can promote the spreading process in a nondominant network with unidirectional interaction. We assume that network A is dominant and network B is nondominant, and A affects B although the inverse is not the case. The simulation result which is not shown here for brevity is similar to Fig. 3. Thus, it reveals that the epidemic threshold of B in cooperative unidirectional interaction networks is well below the isolated network B. Therefore, we obtain the result with Eq. (26).

We provide three colormaps to illustrate the overall behavior. Figs. 5–7 show that the expected infection ratios ρ=(ρ1+ρ2)/2 for multiplex networks as a function of parameters β1 and β2. The white lines in the lower left corner are used to separate the epidemic survival region from epidemic die out region. ρ = 0 in this lower left smaller region means that the epidemic is eventually dying out. ρ>0 in the upper right region means that the epidemic is going to be eventually persistent in the population. Thus, horizontal and vertical white lines represent the epidemic thresholds of complex network A and complex network B, respectively. Comparing the horizontal and vertical white lines in Figs. 5–7, we show that epidemic thresholds are decreased for cooperative multiplex networks using both first order approximation and second order approximation as compared to the epidemic thresholds of corresponding isolated networks. At the same time, Figs. 5 and 6 show that the die out regions differ trivially between first and second order approximation. Therefore, the simulation results are consistent with that of the previous simulations.

Finally, we use simulations for sensitivity analysis of parameters. The second order approximation for cooperative spreading processes is used in numerical simulations. Considering the interacting network A, the comparison of ρ1 for different values of μ1 with γ1=0.04 is shown in Fig. 8, which reveals that larger μ1 leads to a smaller final infection fraction. Further, the comparison of ρ1 for different values of γ1 with μ1=0.4 is shown in Fig. 9. Simultaneously, we can also obtain similar results which were not displayed for interacting network B. The same observation for different parameters is obtained, which reveals that the conclusion is not sensitive to parameters.

In conclusion, a model for cooperative spreading processes on interacting two-layer networks is proposed. In particular, the epidemic thresholds of interacting two-layer networks can be decreased for cooperative spreading processes, thereby implying that cooperative spreading processes promote the spread of the disease. In theory, the global epidemic onset arises in the interacting networks simultaneously. Thus, a global epidemic threshold can be used to uniform the two epidemic thresholds. This model can also be used for unidirectional coupled networks. The results can provide hints for public health authorities to make effective measures for disease control and prevention.

This work was supported by the National Natural Science Foundation of China under Grant Nos. 61273215, 61573262, and 41301442.

1.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. E
63
,
066117
(
2001
).
2.
M. E. J.
Newman
,
Phys. Rev. E
66
,
016128
(
2002
).
3.
A. E.
Motter
,
T.
Nishikawa
, and
Y. C.
Lai
,
Phys. Rev. E
66
,
065103
(
2002
).
4.
K. T. D.
Eames
and
M. J.
Keeling
,
Math. Biosci.
189
,
115
(
2004
).
5.
F.
Liljeros
,
C. R.
Edling
, and
L. A. N.
Amaral
,
Microbes Infect.
5
,
189
(
2003
).
7.
R.
Pastor-Satorras
and
A.
Vespignani
,
Phys. Rev. Lett.
86
,
3200
(
2001
).
8.
V. M.
Eguíluz
and
K.
Konstantin
,
Phys. Rev. Lett.
89
,
108701
(
2002
).
9.
M.
Boguna
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
Phys. Rev. Lett.
90
,
028701
(
2003
).
10.
M.
Boguñá
and
R.
Pastor-Satorras
,
Phys. Rev. E
66
,
047104
(
2002
).
11.
V.
Colizza
and
A.
Vespignani
,
Phys. Rev. Lett.
99
,
148701
(
2007
).
12.
V.
Colizza
,
R.
Pastor-Satorras
, and
A.
Vespignani
,
Nat. Phys.
3
,
276
282
(
2007
).
13.
V.
Colizza
and
A.
Vespignani
,
J. Theor. Biol.
251
,
450
(
2008
).
14.
M.
Catanzaro
,
M.
Boguñá
, and
R.
Pastor-Satorras
,
Reaction-Diffusion Processes in Scale-Free Networks
(
Springer
,
Heidelberg, Berlin
,
2008
).
15.
A.
Buscarino
,
L.
Fortuna
,
M.
Frasca
, and
V.
Latora
,
Europhys. Lett.
82
,
38002
(
2008
).
16.
J.
Zhou
and
Z.
Liu
,
Physica A
388
,
1228
(
2009
).
17.
M.
Tang
,
Z.
Liu
, and
B.
Li
,
Europhys. Lett.
87
,
18005
(
2009
).
19.
F.
Liljeros
,
C. R.
Edling
,
L. A.
Amaral
,
H. E.
Stanley
, and
Y.
Aberg
,
Nature
411
,
907
(
2001
).
20.
D. M.
Morens
,
G. K.
Folkers
, and
A. S.
Fauci
,
Nature
430
,
242
(
2004
).
21.
S.
Funk
and
V. A. A.
Jansen
,
Phys. Rev. E
81
,
036118
(
2010
).
22.
J.
Tao
,
Y.
Ruan
,
L.
Yin
,
S. H.
Vermund
,
B. E.
Shepherd
,
Y.
Shao
, and
H. Z.
Qian
,
Aids Patient Care STD
27
,
524
(
2013
).
23.
G. C.
Van
,
K.
Vongsaiya
,
C.
Hughes
,
R.
Jenkinson
,
A. L.
Bowring
,
A.
Sihavong
,
C.
Phimphachanh
,
N.
Chanlivong
,
M.
Toole
, and
M.
Hellard
,
AIDS Educ. Prev.
25
,
232
(
2013
).
24.
M.
Dickison
,
S.
Havlin
, and
H. E.
Stanley
,
Phys. Rev. E
85
,
066109
(
2012
).
25.
A.
Saumell-Mendiola
,
M.
Serrano
, and
M.
Boguna
,
Phys. Rev. E
86
,
026106
(
2012
).
26.
C.
Granell
,
S.
Gomez
, and
A.
Arenas
,
Phys. Rev. Lett.
111
,
128701
(
2013
).
27.
D.
Zhao
,
L.
Li
,
H.
Peng
,
Q.
Luo
, and
Y.
Yang
,
Phys. Lett. A
378
,
770
(
2014
).
28.
F. D.
Sahneh
and
C.
Scoglio
,
Phys. Rev. E
89
,
062817
(
2014
).
29.
Y.
Zhao
,
M.
Zheng
, and
Z.
Liu
,
Chaos
24
,
043129
(
2014
).
30.
E.
Cozzo
,
R. A.
Baños
,
S.
Meloni
, and
Y.
Moreno
,
Phys. Rev. E
88
,
050801
(
2013
).
31.
W.
Wang
,
M.
Tang
,
H.
Yang
,
Y.
Do
,
Y. C.
Lai
, and
G.
Lee
,
Sci. Rep.
4
,
5097
(
2014
).
32.
Q.
Guo
,
X.
Jiang
,
Y.
Lei
,
M.
Li
,
Y.
Ma
, and
Z.
Zheng
,
Phys. Rev. E
91
,
012822
(
2015
).
33.
S.
Gomez
,
A.
Arenas
,
J.
Borge-Holthoefer
,
S.
Meloni
, and
Y.
Moreno
,
Europhys. Lett.
89
,
38009
(
2010
).
34.
X.
Wu
,
W. X.
Zheng
, and
J.
Zhou
,
Chaos
19
,
013109
(
2009
).