A unified framework of mutual influence between two pathogens in multiplex networks Yanping Zhao, Muhua Zheng, and Zonghua Liu Citation: Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 043129 (2014); doi: 10.1063/1.4902254 View online: http://dx.doi.org/10.1063/1.4902254 View Table of Contents: http://scitation.aip.org/content/aip/journal/chaos/24/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in The impact of vaccine failure rate on epidemic dynamics in responsive networks Chaos 25, 043116 (2015); 10.1063/1.4919245 Epidemic spreading in time-varying community networks Chaos 24, 023116 (2014); 10.1063/1.4876436 Evolution of N-species Kimura/voter models towards criticality, a surrogate for general models of accidental pathogens AIP Conf. Proc. 1479, 1331 (2012); 10.1063/1.4756401 The impact of awareness on epidemic spreading in networks Chaos 22, 013101 (2012); 10.1063/1.3673573 The barrier method: A technique for calculating very long transition times J. Chem. Phys. 133, 124103 (2010); 10.1063/1.3485285
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
CHAOS 24, 043129 (2014)
A unified framework of mutual influence between two pathogens in multiplex networks Yanping Zhao, Muhua Zheng, and Zonghua Liua) Department of Physics, East China Normal University, Shanghai 200062, China
(Received 13 September 2014; accepted 11 November 2014; published online 20 November 2014) There are many evidences to show that different pathogens may interplay each other and cause a variety of mutual influences of epidemics in multiplex networks, but it is still lack of a framework to unify all the different dynamic outcomes of the interactions between the pathogens. We here study this problem and first time present the concept of state-dependent infectious rate, in contrast to the constant infectious rate in previous studies. We consider a model consisting of a two-layered network with one pathogen on the first layer and the other on the second layer, and show that all the different influences between the two pathogens can be given by the different range of parameters in the infectious rates, which includes the cases of mutual enhancement, mutual suppression, and even initial cooperation (suppression) induced final suppression (acceleration). A theoretical C 2014 AIP Publishing LLC. analysis is present to explain the numerical results. V [http://dx.doi.org/10.1063/1.4902254]
Epidemic spreading has been studied for a long time and has recently gotten great attention on complex networks. The main focus is on how the network topology influences epidemic spreading and how epidemic spreads on multilayered networks. Recently, people begin to consider the case of two interacting pathogens. In these studies, the infectious rate is always supposed to be a constant. However, in realistic situation, two pathogens in a network dynamically interact each other and the influence from one to another depends on their individual dynamical status. For example, an infected pathogen-I will make the host weaker and then enhance the spreading of the pathogen-II. While a refractory pathogen-I will make the host be immunized and then suppress the spreading of the pathogen-II. Therefore, the constant infectious rate is not correct in the case of two interacting pathogens. In this paper, we try to consider the dynamical interaction by presenting a concept of state-dependent infectious rate. We find that based on the varying infectious rate, all the different kinds of interaction between two pathogens can be unified into a common framework.
I. INTRODUCTION
Epidemic spreading on complex networks has been well studied in the past decade and many significant results have been achieved.1 The study was first focused on the static networks where each node represents an immobile agent and the contagion occurs only between the neighboring nodes through links. It was revealed that for scale-free networks, there is no epidemic threshold in the thermodynamic limit.2–9 Then, the study was moved to the reaction-diffusion model where agents can move to their neighboring nodes by a possibility.10–14 In this model, the contagious process takes place only within the agents on the same node and there is a)
Electronic mail: [email protected]
1054-1500/2014/24(4)/043129/7/$30.00
no contagion between the neighboring nodes, i.e., the links are used only for diffusion. Later on, the random diffusion of epidemics was replaced by the objective traveling of human being.15–19 After that, some adaptive networks were proposed to discuss the interplay between epidemic spreading and network structure,20–23 where adaptive rewiring of connections among agents is included to represent the selfprotection of human beings. Furthermore, the study on epidemic spreading was extended to the case of traffic-driven epidemic spreading17,24–27 and the case of epidemic spreading at links.18 These quantitative, efficient, and predictive models are very useful for public health authorities to assess situations quickly, make informed decisions, and optimize vaccination and drug delivery plans. Recently, the study is moving to the case of two-layered networks.28–31 In this case, epidemic spreading may be transmitted from one layer (e.g., animals) to the other (e.g., human beings)32 such as the SARS epidemic of 2003, the 2009 H1N1 influenza pandemic, and 2013 H7N9 avian influenza, i.e., a process that continued throughout the twentieth century and up to the present. A two-layered network may also consist of one layer for physical network and the other for online information network.33–40 In this case, besides the infectious contact, the obtained information will jointly influence the epidemic spreading. In conclusion, all these studies on two-layered networks are focused only on the epidemic spreading of one pathogen. However, in reality, it is possible for two or more pathogens to coexist. For example, in south-east Asia, the type-I of adenoviruses could exist for 15 months at one end of a village, whereas the type-II predominated at the other end. In the center of the village, both types coexisted.41,42 Newman discussed the threshold effect for two pathogens spreading on a single network by the bond percolation approach,43 where two diseases compete for the same population of hosts. There are also evidences to show that the two pathogens may stay in different networks but influence each other.
24, 043129-1
C 2014 AIP Publishing LLC V
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-2
Zhao, Zheng, and Liu
For example, for the sexually transmitted diseases, it can propagate both in heterosexual and homosexual networks of sexual contacts. These two networks are not completely isolated due to the existence of bisexual individuals. A study on the number of partners among U.S. bisexual men shows that the bisexual men are the medium to connect the pathogens in the network of heterosexual men with the pathogens in the network of homosexual men.44–48 More examples may come from the interplay between two different pathogens such as flu virus and sexual diseases or Tuberculosis and AIDS, where the former spreads on a respiratory contact network while the latter spreads on a sexual contact network. Thus, it is interesting to study the interaction between two pathogens spreading on two-layered networks. Current studies on this topic are focused on the bidirectional interactions and it has been found that for the case of competitive (cooperative) interaction between two pathogens, the presence of one pathogen makes the other less (more) likely to spread.43,49–55 Of course, it is also possible for the interaction between two pathogens to be unidirectional, which has not been considered so far. All of these cases have corresponding examples in reality. For examples, an infected person will become weakened and thus can be more likely to be attacked by another pathogen, which corresponds to the case of cooperation between two pathogens. When an infected person takes medicine or stay at home, he will be less likely to get the other, which corresponds to the case of competition between two pathogens. And when the recovery from one pathogen leaves an individual weakened, he will be more easy to acquire another pathogen, which corresponds to a new type of interaction between two pathogens, and so on. Therefore, it is necessary to systematically study the interactions. We here study this problem and present a model to unify all the kinds of interaction into a unified framework. The model incorporates the mutual influences between two pathogens into a state-dependent infectious rate, in contrast to the constant infectious rate in previous studies. We here consider the pathogen-I with the susceptible-infected-susceptible (SIS) model and the pathogen-II with the susceptible-infected-refractory (SIR) model. We find that the stabilized epidemics depends on the range of parameters in the infectious rate and may result in an interesting phenomenon, i.e., initial cooperation (suppression) induced final suppression (acceleration). We also discuss the effect of time delay. The numerical results have been explained by a theoretical analysis. The rest of this paper is organized as follows. In Sec. II, we present a model to study the mutual influences between two pathogens on a two-layered network, where a set of parameters is introduced to represent different kinds of interaction. Then in Sec. III, we make numerical simulations to show the effects of the interactions between two pathogens. After that, we give a brief theoretical analysis to the obtained numerical results in Sec. IV. Finally, we give discussions and conclusions in Sec. V. II. MODEL
We first construct two independent random Erd}os-Renyi networks A and B and let them have the same size N and the
Chaos 24, 043129 (2014)
same average degree hki,56 thus the nodes on the two networks will have one-to-one correspondence. We randomly match two nodes from the networks A and B as a pair. For simplicity, we let each pair of the corresponding nodes represent an identical host, i.e., each node has two set of links with one on the network A and the other on the network B. Then, we assume that two pathogens spread on the two networks, respectively. The interaction between the two pathogens occurs dynamically once they spread to the same nodes. Figure 1 shows the schematic figure of multiple epidemic spreading on a two-layered network where the dashed lines imply that each pair of connected nodes represent an identical host and the “blue,” “red,” and “black” circles denote the status of susceptible, infected and refractory, respectively. To understand how the interaction between the two pathogens influences their spreading, instead of considering the same kind of pathogens on the two networks as before, we here assume that the pathogen-I on the network A satisfies the SIS model and the pathogen-II on the network B satisfies the SIR model, where S, I, and R represent the susceptible, infected and refractory individuals, respectively. Other choices such as SIS to SIS or SIR to SIR are also fine. The SIS model is for diseases that cannot be immunized, and an infected person will recover and become susceptible again, as in tuberculosis and gonorrhea.2,3,5,57–59 In an isolated SIS model, a susceptible node may be infected by an infected neighbor at rate b1. When a susceptible node has kinf1 infected neighbors, it will become infected with probability 1 ð1 b1 Þkinf 1 . At the same time, each infected node will become susceptible at rate l at each time step. The SIR model is for diseases that can be immunized such as parotitis, measles, chickenpox, pertussis, and influenza.4,57,60,61 In an isolated SIR model, a susceptible node may be infected by an infected neighbor at rate b2 and will also become
FIG. 1. Schematic figure of multiple epidemic spreading on different networks where A and B represent two different networks with the same size N and the same average degree hki, respectively, and the dashed lines imply that each pair of connected nodes represent an identical host. Two different pathogens are spreading on the two networks, i.e., A and B, respectively, where the “blue,” “red,” and “black” circles denote the status of susceptible, infected and refractory, respectively. a1(t) (a2(t)) represents the influence of the network B (A) on the network A (B).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-3
Zhao, Zheng, and Liu
Chaos 24, 043129 (2014)
infected with the probability 1 ð1 b2 Þkinf 2 with kinf2 being the number of its infected neighbors. At the same time, the infected node will decay into a refractory one with probability l. Without loss of generality, we set l ¼ 1 since it only affects the definition of the time scale of epidemic spreading. In Fig. 1, the two networks are not isolated but connected by the dashed lines. Once a node is infected by the pathogen-I, its immunity against the pathogen-II will be significantly changed, and vice-versa. For the same reason, the refractory status of the pathogen-II will also influence the spreading of the pathogen-I. Considering that the status of a node will change with time, the interplay between the two pathogens will depend on time. Thus, for the SIS model on the layer A, the infection probability 1 ð1 b1 Þkinf 1 of an isolated network has to be replaced by 1 ð1 a1 ðtÞb1 Þkinf 1 with a1(t) representing the influence from the pathogen-II on the layer B. In detail, we let a1(t) ¼ a11, a12 and a13 when the pathogen-II is in the status of S, I, and R, respectively. Obviously, we always have a11 ¼ 1. The values of a12 > 1 and 1 or 1, a13 < 1, and a22 > 1, which represent the effects of acceleration by the infected status and suppression by the refractory status. We let the network size be N ¼ 2000 and the average degree be hki ¼ 6, with ki being the degree of node i. An infectious process is initialized by having 0.5% of the nodes randomly infected for both the two pathogens on the two networks. To measure the epidemic spreading, we let qI1 represent the density of infected nodes in the network A, and let qI2 and qR2 represent the densities of infected and refractory nodes in the network B, respectively. Figure 2 shows the evolution of the two epidemics on the two networks with b1 ¼ b2 ¼ 0.25, where the solid lines represent the case of having state-dependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2. To compare with the spreading of a single
FIG. 2. Evolution of the two epidemics on the two networks with b1 ¼ b2 ¼ 0.25, where the solid lines represent the case of having statedependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2, the dashed lines represent the case without interaction between the two pathogens with the set of parameters a1(t) ¼ a2(t) ¼ 1, and the dash-dotted lines represent the case of having constant interaction between the two pathogens with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. (a) qI1 versus t, where the right inset represents the log plot and the left inset represents the amplification of qI1 for t 2 [7, 9]. (b) qI2 and qR2 versus t.
disease, we have calculated the case without interaction between the two pathogens, i.e., a1(t) ¼ a2(t) ¼ 1, shown in Fig. 2 by the dashed lines. To compare with the previous studies of having constant interaction between the two pathogens, we have also calculated the case of constant a1(t) and a2(t), shown in Fig. 2 by the dash-dotted lines with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. In this case, we let a node in the network A have a1(t) ¼ 1.2 once it is (or was) infected in the network B. That is, its value of a1(t) ¼ 1.2 will be kept after its first infection by the second disease. In the same way, we let a node in the network B have a2(t) ¼ 1.2 once it is (or was) infected in the network A. Fig. 2(a) represents the evolution of qI1 , where the right inset represents the log plot and the left inset represents the amplification of qI1 for t 2 [7, 9]. Fig. 2(b) represents the evolution of qI2 and qR2 . From Fig. 2, we see that (i) both qI1 and qR2 show a rapid increase at a short time followed by reaching a steady-state value, indicating that a stabilized state has been reached. And (ii) there is an approximately exponential growth with exponent c 0.48 in the initial evolutionary process of the pathogen-I, which will be explained in the theoretical part. Moreover, we surprisingly find from the left inset of Fig. 2(a) that the solid line is higher than the dashed line in the initial stage but become the contrary in stabilized stage, indicating that an initial cooperative interaction between the two pathogens results in a final suppressive effect! This contradictory behavior can be explained by combining Fig. 2(a) with 2(b). The initial cooperative interaction results in more qI2 and thus more qR2 . In the stabilized state, we have qI2 ¼ 0 and thus the influence of the SIR pathogen on the SIS pathogen is only through qR2 . More qR2 will make less qI1 , which is just what we have observed in the stabilized state of Fig. 2(a). To measure the influence of the parameter b1 on epidemic spreading, we focus on the steady status. Fig. 3(a)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-4
Zhao, Zheng, and Liu
FIG. 3. Influence of the parameters b1 and b2 on epidemic spreading where the solid lines represent the case of having state-dependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2, the dashed lines represent the case lacking interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1, and the dash-dotted lines represent the case of having constant interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. All the results are obtained by taking the average on 100 different realizations. (a) qI1 and qR2 versus b1 for fixed b2 ¼ 0.25; (b) qR2 and qI1 versus b2 for fixed b1 ¼ 0.25.
shows the dependence of the stabilized qI1 and qR2 on b1 for b2 ¼ 0.25, where the solid lines represent the case of having state-dependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2, the dashed lines represent the case lacking interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1, and the dash-dotted lines represent the case of having constant interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. All the results are obtained by taking the average on 100 different realizations. From the three lines of qI1 in Fig. 3(a), we interestingly find that the threshold b1c (the point where qI1 changes from zero to nonzero) for the case with constant interaction is less than that without interaction, while the value of b1c for the case with statedependent interaction is larger than that without interaction, indicating that the case with constant interaction will speed up the epidemic spreading while the case with statedependent interaction will slow down the epidemic spreading. Similar to Fig. 2(a), this abnormal phenomenon may be understood as follows. Before the steady status, the existence of qI2 does enhance qI1 . But in the steady status, qI2 will become zero and thus the influence of the pathogen-II on the pathogen-I comes only from qR2 , which will increase nodes’ immunity against the pathogen-I and thus reduce qI1 and increase b1c. We will discuss this problem further in Fig. 6. From the solid lines of Fig. 3(a), we find that the existence of qI1 has resulted in the increase of qR2 through accelerating qI2 , indicating the coinfections between the two pathogens. Similarly, Fig. 3(b) shows the dependence of qR2 and qI1 on b2 for fixed b1 ¼ 0.25. From Fig. 3(b), we find that the threshold b2c for the case of having state-dependent interaction is smaller than that of having constant interaction and is further smaller than that without interaction, indicating the accelerating effect of the SIR disease in the network B by
Chaos 24, 043129 (2014)
FIG. 4. Effect of time delay with the set of parameters a11 ¼ 1, a12 ¼ 2, a13 ¼ 0.6 and a21 ¼ 1, a22 ¼ 2, where the three lines represent the evolution of qI1 ; qI2 , and qR2 , respectively. (a) The case for a time delay s ¼ 10 on the network B. (b) The case for a time delay s ¼ 10 on the network A.
the SIS disease in the network A. From the solid line of qI1 in Fig. 3(b), we find that when b2 > b2c, qI1 decreases with the increase of b2, confirming again the suppressive influence from the pathogen-II. Considering that the two epidemics may not be initiated at the same time, it is necessary to discuss the effect of time delay. Figure 4(a) shows the situation that the epidemic spreading on the network B is delayed for a time s ¼ 10, where the three lines represent the evolution of qI1 ; qI2 , and qR2 , respectively. Comparing Fig. 4(a) with Fig. 2(a), we see that the time delay s has slightly influenced the stabilized value of qI1 . Similarly, Fig. 4(b) shows the situation that the epidemic spreading on the network A is delayed for a time s ¼ 10, where the three lines represent the evolution of qI1 ; qI2 , and qR2 , respectively. Comparing Fig. 4(a) with 4(b) we see that their qR2 is different, indicating that different time delays may have different effects. To show the detailed influence, Fig. 5 shows the dependence of stabilized qI1 and qR2 on the delay s where (a) and (b) represent the case that
FIG. 5. Dependence of stabilized qI1 and qR2 on the delay s with the set of parameters a11 ¼ 1, a12 ¼ 2, a13 ¼ 0.6 and a21 ¼ 1, a22 ¼ 2. (a) and (b): The case for a time delay s on the network B. (c) and (d): The case for a time delay s on the network A.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-5
Zhao, Zheng, and Liu
the epidemic spreading on the network B is delayed s while (c) and (d) represent the case that the epidemic spreading on the network A is delayed s. From Figs. 5(a) and 5(b), we see that with the increase of s, qI1 will slightly decrease while qR2 will slightly increase. From Figs. 5(c) and 5(d), we see that with the increase of s, qI1 will slightly increase while qR2 will significantly decrease. Therefore, the effects of time delay from the two pathogens are asymmetric. Based on the conditions of a12 > 1, a13 < 1, and a22 > 1, we now examine how the concrete values of a12, a13, and a22 influence the two epidemics. For this purpose, Figs. 6(a)–6(f) show the dependence of qI1 and qR2 on a12, a13, and a22, respectively, where the parameters are taken as b1 ¼ b2 ¼ 0.25, and a13 ¼ 0.6 and a22 ¼ 2 in (a) and (b), a12 ¼ 2 and a22 ¼ 2 in (c) and (d), and a12 ¼ 2 and a13 ¼ 0.6 in (e) and (f). From Figs. 6(a)–6(f) we find that qI1 may either increase or decrease, depending on the different influences from the pathogen-II, while qR2 is always enhanced. We now turn to other cases of the parameters a12, a13, and a22. Without loss of generality, we here fix a12 ¼ a22. Considering that a13 > 1 and < 1 represent the effects of acceleration and suppression, respectively, we choose two typical a13 and study the relationship between qI1 and a12. Fig. 7(a) shows the results where the “squares” and “circles” represent the cases of a13 ¼ 0.5 and 2, respectively. We see that qI1 decreases with the increase of a12 for the case of a13 ¼ 0.5 but increases with a12 for the case of a13 ¼ 2. This result can be explained as follows. For the former, the increase of a12 will result in the increase of qR2 and then suppress qI1 . For the latter, the increase of qR2 and the accelerating effect of a13 will compete and then results in the slow increase of qI1 . Similarly, Fig. 7(b) shows the relationship between qI1 and a13 where the “squares” and “circles” represent the cases of a12 ¼ 0.5 and 2, respectively. It is interesting to notice that the two curves have a cross point at a13 ¼ 1. Before the cross point, qI1 for a12 ¼ 0.5 will be greater than that for a12 ¼ 2. While after the cross point, qI1 for a12 ¼ 0.5 will be less than that for a12 ¼ 2. The reason is
FIG. 6. Dependence of qI1 and qR2 on a12, a13, and a22, respectively, where the parameters are taken as a11 ¼ 1, a13 ¼ 0.6 and a21 ¼ 1, a22 ¼ 2 in (a) and (b), a11 ¼ 1, a12 ¼ 2 and a21 ¼ 1, a22 ¼ 2 in (c) and (d), and a11 ¼ 1, a12 ¼ 2, a13 ¼ 0.6 and a21 ¼ 1 in (e) and (f).
Chaos 24, 043129 (2014)
FIG. 7. Comparison of the effects between acceleration and suppression with fixed a12 ¼ a22 and b1 ¼ b2 ¼ 0.25. (a) qI1 versus a12 where the “squares” and “circles” represent the cases of a13 ¼ 0.5 and 2, respectively. (b) qI1 versus a13 where the “squares” and “circles” represent the cases of a12 ¼ 0.5 and 2, respectively.
that before the cross point, a13 behaves as suppression and thus larger a12 will result in less qI1 . While after the cross point, a13 behaves as acceleration and thus larger a12 will result in more qI1 . A more comprehensive study is shown in Fig. 8, where qI1 is plotted on both the parameters a12 and a13. We see that qI1 will be zero when a12 is larger and a13 is small, decrease with the increase of a12 when a13 < 1, increase with a12 when a13 > 1, and reach saturation when both a12 and a13 are larger. Therefore, the different parameter regions in Fig. 8 represent different interactions between the two pathogens such as the cooperation, suppression and so on, implying that all of them can be unified into the same framework of Eq. (1). IV. A BRIEF THEORETICAL ANALYSIS
We now conduct a theoretical analysis to the above numerical simulations. In the framework of mean-field, the evolutionary equations of the two epidemics can be given as
FIG. 8. Dependence of qI1 on a12 and a13 with fixed a12 ¼ a22 and b1 ¼ b2 ¼ 0.25.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-6
Zhao, Zheng, and Liu
Chaos 24, 043129 (2014)
dqI1 ¼ qI1 þ a1 ðtÞb1 qI1 qS1 hki; dt dqI2 ¼ qI2 þ a2 ðtÞb2 qI2 qS2 hki; dt dqR2 ¼ qI2 ; dt
From the second and third equations of Eq. (2), we have dqS2 =dt ¼ a2 ðtÞb2 qI2 qS2 hki. In the initial condition b2 hki of Ð t qR2 ð0Þ ¼ 0 and qS2 ð0Þ 1, we have qS2 ðtÞ ¼ e 0 a2 ðtÞqI2 ðtÞdt. Thus, in the stabilized state, we have 0
qR2 ð1Þ ¼ 1 ea2 b2 hkiqR2 ð1Þ ;
(2)
with the conservation laws qS1 þ qI1 ¼ 1 and qS2 þ qI2 þqR2 ¼ 1. In the mean-field level, the coefficient a1(t) and a2(t) can be given as
where a02 is defined by 0 a2 ðtÞqI2 ðtÞdt a02 qR2 ð1Þ and can be considered as an average of a2(t). From Eq. (3) we know a2(t) > 1, which gives a02 > 1. By the condition of a0 b hkiq
dð1e 2 2 R2 dqR2 ð1Þ
a1 ðtÞ ¼ a11 qS2 ðtÞ þ a12 qI2 ðtÞ þ a13 qR2 ðtÞ; a2 ðtÞ ¼ a21 qS1 ðtÞ þ a22 qI1 ðtÞ;
qI1 ðtÞ ¼
qI1 ð0ÞðC 1ÞeðC1Þt : C 1 þ CqI1 ð0ÞeðC1Þt
(4)
As qI1 ð0Þ 0, we ignore the part of CqI1 ð0ÞeðC1Þt in the denominator and thus obtain qI1 ðtÞ qI1 ð0Þe
ðC1Þt
;
(5)
with approximately exponential growth rate (C 1). For the parameters in Fig. 2, we have C 1 ¼ a11 b1 hki 1 ¼ 0:5, which is approximately consistent with the slope c ¼ 0.48 in the right inset of Fig. 2(a). Notice from Fig. 2 that there is a time delay between the setting up of steady-state of the two diseases, i.e., the SIR model in Fig. 2(b) will reach its steady-state at about t ¼ 20, while the SIS model in Fig. 2(a) will reach its steady-state at about t > 25. Thus, the evolution of the SIS model for t > 20 is not a co-evolution of two pathogens but a single pathogen evolution. In this sense, we cannot naturally obtain the solution of steady-state of Eq. (2). Therefore, we here focus on the critical points. In the stabilized state, we have qI2 ¼ 0 and qI1 and qR2 constant. Thus, a1 ¼ a11 ð1 qR2 Þ þ a13 qR2 and a2 ¼ a21 ð1 qI1 Þ þ a22 qI1 are both constants. When b1 b1c, we have dqI1 =dt ¼ 0 in the first equation of Eq. (2), which gives qI1 ¼ 1
1 : a1 b1 hki
(6)
The critical point b1c can be obtained by letting qI1 change from zero to nonzero, which gives b1c ¼
1 : a1 hki
ð1Þ
Þ
jqR
2
ð1Þ¼0
> 1, we obtain the critical point
(3)
with a11 ¼ a21 ¼ 1. To get the solution of Eq. (2), we take the initial conditions of qS1 ð0Þ 1; qI1 ð0Þ 0; qS2 ð0Þ 1; qI2 ð0Þ 0 and qR2 ð0Þ ¼ 0. In the initial steps, we have a1(t) a11. From the first equation of Eq. (2), we have dqI1 =dt ¼ qI1 þ CqI1 ð1 qI1 Þ ¼ ðC 1ÞqI1 Cq2I1 with C ¼ a11 b1 hki. Its solution is
(7)
It will return to the threshold bc ¼ 1=hki in an isolated SIS random network.62 For the parameters in Fig. 3, we have qR2 0:57 and thus obtain a1 ¼ 0.77. Substituting a1 into Eq. (7) we have b1c ¼ 0.21, which is approximately consistent with the threshold b1c ¼ 0.2 in Fig. 3(a) and thus confirms the numerical result.
(8)
Ð1
b2c ¼
1 ; a02 hki
(9)
which depends on the first pathogen through a02 or a2(t). b2c will also return to the threshold bc ¼ 1=hki in an isolated SIR random network. For the parameters in Fig. 3, we have qI1 0:28 and thus a2 ¼ 1.28, which gives a02 1:14. Substituting it into Eq. (9) we have b2c ¼ 0.146, which is consistent with the threshold b2c ¼ 0.14 in Fig. 3(b) and thus confirms the numerical result. V. DISCUSSIONS AND CONCLUSIONS
The present study is for the interactions between the two pathogens with the SIS and SIR models, respectively, but it can be easily extended to the case between two SIS pathogens or the case between two SIR pathogens and so on. The framework of Eq. (1) can be further extended to the case of more than two pathogens. On the other hand, the considered mutual influences in this work are only for the steady status. From Fig. 2, we see that in the initial evolution process, both qI1 and qI2 are enhanced, which is one of the main features to cause crowd panicked. Thus, the control of epidemic spreading in the cases of multiple pathogens should be paid more attention to the initial evolution process and thus deserved to be studied further. In conclusions, we have presented a unified framework to study the cases of two pathogens on two networks and figured out that the key point is the time dependent infectious rate, in contrast to the previous constant infectious rate. Through this framework we find that the final infectious depends on the range of the set of parameters a12, a13, and a21 and also on the types of pathogens such as SIS or SIR. This framework thus opens a new window to study the epidemic spreading on multiple networks with multiple pathogens and may stimulate further studies. ACKNOWLEDGMENTS
This work was partially supported by the NNSF of China under Grant Nos. 11135001 and 11375066, Joriss project under Grant No. 78230050, and 973 Program under Grant No. 2013CB834100. 1
S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008). 2 R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001). 3 V. M. Eguiluz and K. Klemm, Phys. Rev. Lett. 89, 108701 (2002).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-7
Zhao, Zheng, and Liu
Chaos 24, 043129 (2014)
4
32
5
33
M. E. J. Newman, Phys. Rev. E 66, 016128 (2002). M. Boguna, R. Pastor-Satorras, and A. Vespignani, Phys. Rev. Lett. 90, 028701 (2003). 6 T. Gross, C. J. D. D’Lima, and B. Blasius, Phys. Rev. Lett. 96, 208701 (2006). 7 R. Parshani, S. Carmi, and S. Havlin, Phys. Rev. Lett. 104, 258701 (2010). 8 C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. 105, 218701 (2010). 9 Y. Schwarzkopf, A. Rakos, and D. Mukamel, Phys. Rev. E 82, 036112 (2010). 10 V. Colizza, R. Pastor-Satorras, and A. Vespignani, Nat. Phys. 3, 276 (2007). 11 V. Colizza and A. Vespignani, Phys. Rev. Lett. 99, 148701 (2007). 12 V. Colizza and A. Vespignani, J. Theor. Biol. 251, 450 (2008). 13 A. Baronchelli, M. Catanzaro, and R. Pastor-Satorras, Phys. Rev. E 78, 016111 (2008). 14 M. Tang, L. Liu, and Z. Liu, Phys. Rev. E 79, 016108 (2009). 15 M. Tang, Z. Liu, and B. Li, Europhys. Lett. 87, 18005 (2009). 16 Z. Liu, Phys. Rev. E 81, 016110 (2010). 17 Z. Ruan, P. M. Hui, H. Q. Lin, and Z. Liu, Eur. Phys. J. B 86, 13 (2013). 18 Z. Ruan, M. Tang, and Z. Liu, Eur. Phys. J. B 86, 149 (2013). 19 P. Grassberger, J. Stat. Mech.: Theory Exp. 2013, P04004. 20 D. H. Zanette and S. Risau-Gusman, J. Biol. Phys. 34, 135 (2008). 21 V. Marceau, P.-A. Noel, L. Hebert-Dufresne, A. Allard, and L. J. Dube, Phys. Rev. E 82, 036116 (2010). 22 I. B. Schwartz and L. B. Shaw, Physics 3, 17 (2010). 23 C. Lagorio, M. Dickison, F. Vazquez, L. A. Braunstein, P. A. Macri, M. V. Migueles, S. Havlin, and H. E. Stanley, Phys. Rev. E 83, 026102 (2011). 24 A. Vazquez, B. Racz, A. Lukacs, and A. Barabasi, Phys. Rev. Lett. 98, 158702 (2007). 25 F. C. Coelho, O. G. Cruz, and C. T. Codeco, Source Code Biol. Med. 3, 3 (2008). 26 S. Meloni, A. Arenas, and Y. Moreno, Proc. Natl. Acad. Sci. 106, 16897 (2009). 27 D. Balcan, V. Colizza, B. Goncalves, H. Hu, J. J. Ramasco, and A. Vespignani, Proc. Natl. Acad. Sci. 106, 21484 (2009). 28 S. Funk and V. A. A. Jansen, Phys. Rev. E 81, 036118 (2010). 29 M. Dickison, S. Havlin, and H. E. Stanley, Phys. Rev. E 85, 066109 (2012). 30 A. Saumell-Mendiola, M. A. Serrano, and M. Boguna, Phys. Rev. E 86, 026106 (2012). 31 F. D. Sahneh, C. Scoglio, and P. Van Mieghem, IEEE/ACM Trans. Networking 21, 1609 (2013).
D. M. Morens, G. K. Folkers, and A. S. Fauci, Nature 430, 242 (2004). C. Granell, S. Gomez, and A. Arenas, Phys. Rev. Lett. 111, 128701 (2013). 34 G. Su, Z. Ruan, S. Guan, and Z. Liu, Europhys. Lett. 103, 48004 (2013). 35 Z. Ruan, M. Tang, and Z. Liu, Phys. Rev. E 86, 036117 (2012). 36 S. Funk, E. Gilad, C. Watkins, and V. A. A. Jansen, Proc. Natl. Acad. Sci. 106, 6872 (2009). 37 F. Bagnoli, P. Lio, and L. Sguanci, Phys. Rev. E 76, 61904 (2007). 38 Q. Wu, X. Fu, M. Small, and X. Xu, Chaos 22, 013101 (2012). 39 H. Zhang, J. Zhang, C. Zhou, M. Small, and B. Wang, New J. Phys. 12, 023015 (2010). 40 A. Perisic and C. T. Bauch, PLoS Comput. Biol. 5, e1000280 (2009). 41 K. Dietz, J. Math. Biol. 8, 291 (1979). 42 F. B. Bang, Int. J. Epidemiol. 4, 337 (1975). 43 M. E. J. Newman, Phys. Rev. Lett. 95, 108701 (2005). 44 W. L. Jeffries, Perspect. Sexual Reprod. Health 43, 151 (2011). 45 W. L. Jeffries and B. Dodge, J. Sex Res. 44, 278 (2007). 46 C. Goodenow, J. Netherland, and L. Szalacha, Am. J. Public Health 92, 203 (2002). 47 L. B. Hightow et al., Sexually Transmitted Dis. 33, 585 (2006). 48 C. T. Bauch and A. P. Galvani, Science 342, 47 (2013). 49 L. J. Abu-Raddad, P. Patnaik, and J. G. Kublin, Science 314, 1603 (2006). 50 M. Martcheva and S. S. Pilyugin, SIAM J. Appl. Math. 66, 843 (2006). 51 V. Marceau, P.-A. Noel, L. Hebert-Dufresne, A. Allard, and L. J. Dube, Phys. Rev. E 84, 026105 (2011). 52 J. C. Miller, Phys. Rev. E 87, 060801(R) (2013). 53 M. E. J. Newman and C. R. Ferrario, PLoS One 8, e71321 (2013). 54 L. Chen, F. Ghanbarnejad, W. Cai, and P. Grassberger, Europhys. Lett. 104, 50001 (2013). 55 H. L. Mills, A. Ganesh, and C. Colijn, J. Theor. Biol. 320, 47–57 (2013). 56 R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47 (2002). 57 R. M. Anderson and R. M. May, Infectious Diseases of Humans: Dynamics and Control (Oxford University Press, Oxford, New York, 1991). 58 R. Olinky and L. Stone, Phys. Rev. E 70, 030902 (2004). 59 Z. Liu and B. Hu, Europhys. Lett. 72, 315 (2005). 60 Y. Moreno, R. Pastor-Satorras, and A. Vespignani, Eur. Phys. J. B 26, 521 (2002). 61 Y. Moreno, M. Nekovee, and A. Vespignani, Phys. Rev. E 69, 055101 (2004). 62 R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 63, 066117 (2001).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
CHAOS 24, 043129 (2014)
A unified framework of mutual influence between two pathogens in multiplex networks Yanping Zhao, Muhua Zheng, and Zonghua Liua) Department of Physics, East China Normal University, Shanghai 200062, China
(Received 13 September 2014; accepted 11 November 2014; published online 20 November 2014) There are many evidences to show that different pathogens may interplay each other and cause a variety of mutual influences of epidemics in multiplex networks, but it is still lack of a framework to unify all the different dynamic outcomes of the interactions between the pathogens. We here study this problem and first time present the concept of state-dependent infectious rate, in contrast to the constant infectious rate in previous studies. We consider a model consisting of a two-layered network with one pathogen on the first layer and the other on the second layer, and show that all the different influences between the two pathogens can be given by the different range of parameters in the infectious rates, which includes the cases of mutual enhancement, mutual suppression, and even initial cooperation (suppression) induced final suppression (acceleration). A theoretical C 2014 AIP Publishing LLC. analysis is present to explain the numerical results. V [http://dx.doi.org/10.1063/1.4902254]
Epidemic spreading has been studied for a long time and has recently gotten great attention on complex networks. The main focus is on how the network topology influences epidemic spreading and how epidemic spreads on multilayered networks. Recently, people begin to consider the case of two interacting pathogens. In these studies, the infectious rate is always supposed to be a constant. However, in realistic situation, two pathogens in a network dynamically interact each other and the influence from one to another depends on their individual dynamical status. For example, an infected pathogen-I will make the host weaker and then enhance the spreading of the pathogen-II. While a refractory pathogen-I will make the host be immunized and then suppress the spreading of the pathogen-II. Therefore, the constant infectious rate is not correct in the case of two interacting pathogens. In this paper, we try to consider the dynamical interaction by presenting a concept of state-dependent infectious rate. We find that based on the varying infectious rate, all the different kinds of interaction between two pathogens can be unified into a common framework.
I. INTRODUCTION
Epidemic spreading on complex networks has been well studied in the past decade and many significant results have been achieved.1 The study was first focused on the static networks where each node represents an immobile agent and the contagion occurs only between the neighboring nodes through links. It was revealed that for scale-free networks, there is no epidemic threshold in the thermodynamic limit.2–9 Then, the study was moved to the reaction-diffusion model where agents can move to their neighboring nodes by a possibility.10–14 In this model, the contagious process takes place only within the agents on the same node and there is a)
Electronic mail: [email protected]
1054-1500/2014/24(4)/043129/7/$30.00
no contagion between the neighboring nodes, i.e., the links are used only for diffusion. Later on, the random diffusion of epidemics was replaced by the objective traveling of human being.15–19 After that, some adaptive networks were proposed to discuss the interplay between epidemic spreading and network structure,20–23 where adaptive rewiring of connections among agents is included to represent the selfprotection of human beings. Furthermore, the study on epidemic spreading was extended to the case of traffic-driven epidemic spreading17,24–27 and the case of epidemic spreading at links.18 These quantitative, efficient, and predictive models are very useful for public health authorities to assess situations quickly, make informed decisions, and optimize vaccination and drug delivery plans. Recently, the study is moving to the case of two-layered networks.28–31 In this case, epidemic spreading may be transmitted from one layer (e.g., animals) to the other (e.g., human beings)32 such as the SARS epidemic of 2003, the 2009 H1N1 influenza pandemic, and 2013 H7N9 avian influenza, i.e., a process that continued throughout the twentieth century and up to the present. A two-layered network may also consist of one layer for physical network and the other for online information network.33–40 In this case, besides the infectious contact, the obtained information will jointly influence the epidemic spreading. In conclusion, all these studies on two-layered networks are focused only on the epidemic spreading of one pathogen. However, in reality, it is possible for two or more pathogens to coexist. For example, in south-east Asia, the type-I of adenoviruses could exist for 15 months at one end of a village, whereas the type-II predominated at the other end. In the center of the village, both types coexisted.41,42 Newman discussed the threshold effect for two pathogens spreading on a single network by the bond percolation approach,43 where two diseases compete for the same population of hosts. There are also evidences to show that the two pathogens may stay in different networks but influence each other.
24, 043129-1
C 2014 AIP Publishing LLC V
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-2
Zhao, Zheng, and Liu
For example, for the sexually transmitted diseases, it can propagate both in heterosexual and homosexual networks of sexual contacts. These two networks are not completely isolated due to the existence of bisexual individuals. A study on the number of partners among U.S. bisexual men shows that the bisexual men are the medium to connect the pathogens in the network of heterosexual men with the pathogens in the network of homosexual men.44–48 More examples may come from the interplay between two different pathogens such as flu virus and sexual diseases or Tuberculosis and AIDS, where the former spreads on a respiratory contact network while the latter spreads on a sexual contact network. Thus, it is interesting to study the interaction between two pathogens spreading on two-layered networks. Current studies on this topic are focused on the bidirectional interactions and it has been found that for the case of competitive (cooperative) interaction between two pathogens, the presence of one pathogen makes the other less (more) likely to spread.43,49–55 Of course, it is also possible for the interaction between two pathogens to be unidirectional, which has not been considered so far. All of these cases have corresponding examples in reality. For examples, an infected person will become weakened and thus can be more likely to be attacked by another pathogen, which corresponds to the case of cooperation between two pathogens. When an infected person takes medicine or stay at home, he will be less likely to get the other, which corresponds to the case of competition between two pathogens. And when the recovery from one pathogen leaves an individual weakened, he will be more easy to acquire another pathogen, which corresponds to a new type of interaction between two pathogens, and so on. Therefore, it is necessary to systematically study the interactions. We here study this problem and present a model to unify all the kinds of interaction into a unified framework. The model incorporates the mutual influences between two pathogens into a state-dependent infectious rate, in contrast to the constant infectious rate in previous studies. We here consider the pathogen-I with the susceptible-infected-susceptible (SIS) model and the pathogen-II with the susceptible-infected-refractory (SIR) model. We find that the stabilized epidemics depends on the range of parameters in the infectious rate and may result in an interesting phenomenon, i.e., initial cooperation (suppression) induced final suppression (acceleration). We also discuss the effect of time delay. The numerical results have been explained by a theoretical analysis. The rest of this paper is organized as follows. In Sec. II, we present a model to study the mutual influences between two pathogens on a two-layered network, where a set of parameters is introduced to represent different kinds of interaction. Then in Sec. III, we make numerical simulations to show the effects of the interactions between two pathogens. After that, we give a brief theoretical analysis to the obtained numerical results in Sec. IV. Finally, we give discussions and conclusions in Sec. V. II. MODEL
We first construct two independent random Erd}os-Renyi networks A and B and let them have the same size N and the
Chaos 24, 043129 (2014)
same average degree hki,56 thus the nodes on the two networks will have one-to-one correspondence. We randomly match two nodes from the networks A and B as a pair. For simplicity, we let each pair of the corresponding nodes represent an identical host, i.e., each node has two set of links with one on the network A and the other on the network B. Then, we assume that two pathogens spread on the two networks, respectively. The interaction between the two pathogens occurs dynamically once they spread to the same nodes. Figure 1 shows the schematic figure of multiple epidemic spreading on a two-layered network where the dashed lines imply that each pair of connected nodes represent an identical host and the “blue,” “red,” and “black” circles denote the status of susceptible, infected and refractory, respectively. To understand how the interaction between the two pathogens influences their spreading, instead of considering the same kind of pathogens on the two networks as before, we here assume that the pathogen-I on the network A satisfies the SIS model and the pathogen-II on the network B satisfies the SIR model, where S, I, and R represent the susceptible, infected and refractory individuals, respectively. Other choices such as SIS to SIS or SIR to SIR are also fine. The SIS model is for diseases that cannot be immunized, and an infected person will recover and become susceptible again, as in tuberculosis and gonorrhea.2,3,5,57–59 In an isolated SIS model, a susceptible node may be infected by an infected neighbor at rate b1. When a susceptible node has kinf1 infected neighbors, it will become infected with probability 1 ð1 b1 Þkinf 1 . At the same time, each infected node will become susceptible at rate l at each time step. The SIR model is for diseases that can be immunized such as parotitis, measles, chickenpox, pertussis, and influenza.4,57,60,61 In an isolated SIR model, a susceptible node may be infected by an infected neighbor at rate b2 and will also become
FIG. 1. Schematic figure of multiple epidemic spreading on different networks where A and B represent two different networks with the same size N and the same average degree hki, respectively, and the dashed lines imply that each pair of connected nodes represent an identical host. Two different pathogens are spreading on the two networks, i.e., A and B, respectively, where the “blue,” “red,” and “black” circles denote the status of susceptible, infected and refractory, respectively. a1(t) (a2(t)) represents the influence of the network B (A) on the network A (B).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-3
Zhao, Zheng, and Liu
Chaos 24, 043129 (2014)
infected with the probability 1 ð1 b2 Þkinf 2 with kinf2 being the number of its infected neighbors. At the same time, the infected node will decay into a refractory one with probability l. Without loss of generality, we set l ¼ 1 since it only affects the definition of the time scale of epidemic spreading. In Fig. 1, the two networks are not isolated but connected by the dashed lines. Once a node is infected by the pathogen-I, its immunity against the pathogen-II will be significantly changed, and vice-versa. For the same reason, the refractory status of the pathogen-II will also influence the spreading of the pathogen-I. Considering that the status of a node will change with time, the interplay between the two pathogens will depend on time. Thus, for the SIS model on the layer A, the infection probability 1 ð1 b1 Þkinf 1 of an isolated network has to be replaced by 1 ð1 a1 ðtÞb1 Þkinf 1 with a1(t) representing the influence from the pathogen-II on the layer B. In detail, we let a1(t) ¼ a11, a12 and a13 when the pathogen-II is in the status of S, I, and R, respectively. Obviously, we always have a11 ¼ 1. The values of a12 > 1 and 1 or 1, a13 < 1, and a22 > 1, which represent the effects of acceleration by the infected status and suppression by the refractory status. We let the network size be N ¼ 2000 and the average degree be hki ¼ 6, with ki being the degree of node i. An infectious process is initialized by having 0.5% of the nodes randomly infected for both the two pathogens on the two networks. To measure the epidemic spreading, we let qI1 represent the density of infected nodes in the network A, and let qI2 and qR2 represent the densities of infected and refractory nodes in the network B, respectively. Figure 2 shows the evolution of the two epidemics on the two networks with b1 ¼ b2 ¼ 0.25, where the solid lines represent the case of having state-dependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2. To compare with the spreading of a single
FIG. 2. Evolution of the two epidemics on the two networks with b1 ¼ b2 ¼ 0.25, where the solid lines represent the case of having statedependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2, the dashed lines represent the case without interaction between the two pathogens with the set of parameters a1(t) ¼ a2(t) ¼ 1, and the dash-dotted lines represent the case of having constant interaction between the two pathogens with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. (a) qI1 versus t, where the right inset represents the log plot and the left inset represents the amplification of qI1 for t 2 [7, 9]. (b) qI2 and qR2 versus t.
disease, we have calculated the case without interaction between the two pathogens, i.e., a1(t) ¼ a2(t) ¼ 1, shown in Fig. 2 by the dashed lines. To compare with the previous studies of having constant interaction between the two pathogens, we have also calculated the case of constant a1(t) and a2(t), shown in Fig. 2 by the dash-dotted lines with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. In this case, we let a node in the network A have a1(t) ¼ 1.2 once it is (or was) infected in the network B. That is, its value of a1(t) ¼ 1.2 will be kept after its first infection by the second disease. In the same way, we let a node in the network B have a2(t) ¼ 1.2 once it is (or was) infected in the network A. Fig. 2(a) represents the evolution of qI1 , where the right inset represents the log plot and the left inset represents the amplification of qI1 for t 2 [7, 9]. Fig. 2(b) represents the evolution of qI2 and qR2 . From Fig. 2, we see that (i) both qI1 and qR2 show a rapid increase at a short time followed by reaching a steady-state value, indicating that a stabilized state has been reached. And (ii) there is an approximately exponential growth with exponent c 0.48 in the initial evolutionary process of the pathogen-I, which will be explained in the theoretical part. Moreover, we surprisingly find from the left inset of Fig. 2(a) that the solid line is higher than the dashed line in the initial stage but become the contrary in stabilized stage, indicating that an initial cooperative interaction between the two pathogens results in a final suppressive effect! This contradictory behavior can be explained by combining Fig. 2(a) with 2(b). The initial cooperative interaction results in more qI2 and thus more qR2 . In the stabilized state, we have qI2 ¼ 0 and thus the influence of the SIR pathogen on the SIS pathogen is only through qR2 . More qR2 will make less qI1 , which is just what we have observed in the stabilized state of Fig. 2(a). To measure the influence of the parameter b1 on epidemic spreading, we focus on the steady status. Fig. 3(a)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-4
Zhao, Zheng, and Liu
FIG. 3. Influence of the parameters b1 and b2 on epidemic spreading where the solid lines represent the case of having state-dependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2, the dashed lines represent the case lacking interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1, and the dash-dotted lines represent the case of having constant interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. All the results are obtained by taking the average on 100 different realizations. (a) qI1 and qR2 versus b1 for fixed b2 ¼ 0.25; (b) qR2 and qI1 versus b2 for fixed b1 ¼ 0.25.
shows the dependence of the stabilized qI1 and qR2 on b1 for b2 ¼ 0.25, where the solid lines represent the case of having state-dependent interaction between the two pathogens with the set of parameters a12 ¼ 2, a13 ¼ 0.6, and a22 ¼ 2, the dashed lines represent the case lacking interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1, and the dash-dotted lines represent the case of having constant interaction between them with the set of parameters a1(t) ¼ a2(t) ¼ 1.2. All the results are obtained by taking the average on 100 different realizations. From the three lines of qI1 in Fig. 3(a), we interestingly find that the threshold b1c (the point where qI1 changes from zero to nonzero) for the case with constant interaction is less than that without interaction, while the value of b1c for the case with statedependent interaction is larger than that without interaction, indicating that the case with constant interaction will speed up the epidemic spreading while the case with statedependent interaction will slow down the epidemic spreading. Similar to Fig. 2(a), this abnormal phenomenon may be understood as follows. Before the steady status, the existence of qI2 does enhance qI1 . But in the steady status, qI2 will become zero and thus the influence of the pathogen-II on the pathogen-I comes only from qR2 , which will increase nodes’ immunity against the pathogen-I and thus reduce qI1 and increase b1c. We will discuss this problem further in Fig. 6. From the solid lines of Fig. 3(a), we find that the existence of qI1 has resulted in the increase of qR2 through accelerating qI2 , indicating the coinfections between the two pathogens. Similarly, Fig. 3(b) shows the dependence of qR2 and qI1 on b2 for fixed b1 ¼ 0.25. From Fig. 3(b), we find that the threshold b2c for the case of having state-dependent interaction is smaller than that of having constant interaction and is further smaller than that without interaction, indicating the accelerating effect of the SIR disease in the network B by
Chaos 24, 043129 (2014)
FIG. 4. Effect of time delay with the set of parameters a11 ¼ 1, a12 ¼ 2, a13 ¼ 0.6 and a21 ¼ 1, a22 ¼ 2, where the three lines represent the evolution of qI1 ; qI2 , and qR2 , respectively. (a) The case for a time delay s ¼ 10 on the network B. (b) The case for a time delay s ¼ 10 on the network A.
the SIS disease in the network A. From the solid line of qI1 in Fig. 3(b), we find that when b2 > b2c, qI1 decreases with the increase of b2, confirming again the suppressive influence from the pathogen-II. Considering that the two epidemics may not be initiated at the same time, it is necessary to discuss the effect of time delay. Figure 4(a) shows the situation that the epidemic spreading on the network B is delayed for a time s ¼ 10, where the three lines represent the evolution of qI1 ; qI2 , and qR2 , respectively. Comparing Fig. 4(a) with Fig. 2(a), we see that the time delay s has slightly influenced the stabilized value of qI1 . Similarly, Fig. 4(b) shows the situation that the epidemic spreading on the network A is delayed for a time s ¼ 10, where the three lines represent the evolution of qI1 ; qI2 , and qR2 , respectively. Comparing Fig. 4(a) with 4(b) we see that their qR2 is different, indicating that different time delays may have different effects. To show the detailed influence, Fig. 5 shows the dependence of stabilized qI1 and qR2 on the delay s where (a) and (b) represent the case that
FIG. 5. Dependence of stabilized qI1 and qR2 on the delay s with the set of parameters a11 ¼ 1, a12 ¼ 2, a13 ¼ 0.6 and a21 ¼ 1, a22 ¼ 2. (a) and (b): The case for a time delay s on the network B. (c) and (d): The case for a time delay s on the network A.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-5
Zhao, Zheng, and Liu
the epidemic spreading on the network B is delayed s while (c) and (d) represent the case that the epidemic spreading on the network A is delayed s. From Figs. 5(a) and 5(b), we see that with the increase of s, qI1 will slightly decrease while qR2 will slightly increase. From Figs. 5(c) and 5(d), we see that with the increase of s, qI1 will slightly increase while qR2 will significantly decrease. Therefore, the effects of time delay from the two pathogens are asymmetric. Based on the conditions of a12 > 1, a13 < 1, and a22 > 1, we now examine how the concrete values of a12, a13, and a22 influence the two epidemics. For this purpose, Figs. 6(a)–6(f) show the dependence of qI1 and qR2 on a12, a13, and a22, respectively, where the parameters are taken as b1 ¼ b2 ¼ 0.25, and a13 ¼ 0.6 and a22 ¼ 2 in (a) and (b), a12 ¼ 2 and a22 ¼ 2 in (c) and (d), and a12 ¼ 2 and a13 ¼ 0.6 in (e) and (f). From Figs. 6(a)–6(f) we find that qI1 may either increase or decrease, depending on the different influences from the pathogen-II, while qR2 is always enhanced. We now turn to other cases of the parameters a12, a13, and a22. Without loss of generality, we here fix a12 ¼ a22. Considering that a13 > 1 and < 1 represent the effects of acceleration and suppression, respectively, we choose two typical a13 and study the relationship between qI1 and a12. Fig. 7(a) shows the results where the “squares” and “circles” represent the cases of a13 ¼ 0.5 and 2, respectively. We see that qI1 decreases with the increase of a12 for the case of a13 ¼ 0.5 but increases with a12 for the case of a13 ¼ 2. This result can be explained as follows. For the former, the increase of a12 will result in the increase of qR2 and then suppress qI1 . For the latter, the increase of qR2 and the accelerating effect of a13 will compete and then results in the slow increase of qI1 . Similarly, Fig. 7(b) shows the relationship between qI1 and a13 where the “squares” and “circles” represent the cases of a12 ¼ 0.5 and 2, respectively. It is interesting to notice that the two curves have a cross point at a13 ¼ 1. Before the cross point, qI1 for a12 ¼ 0.5 will be greater than that for a12 ¼ 2. While after the cross point, qI1 for a12 ¼ 0.5 will be less than that for a12 ¼ 2. The reason is
FIG. 6. Dependence of qI1 and qR2 on a12, a13, and a22, respectively, where the parameters are taken as a11 ¼ 1, a13 ¼ 0.6 and a21 ¼ 1, a22 ¼ 2 in (a) and (b), a11 ¼ 1, a12 ¼ 2 and a21 ¼ 1, a22 ¼ 2 in (c) and (d), and a11 ¼ 1, a12 ¼ 2, a13 ¼ 0.6 and a21 ¼ 1 in (e) and (f).
Chaos 24, 043129 (2014)
FIG. 7. Comparison of the effects between acceleration and suppression with fixed a12 ¼ a22 and b1 ¼ b2 ¼ 0.25. (a) qI1 versus a12 where the “squares” and “circles” represent the cases of a13 ¼ 0.5 and 2, respectively. (b) qI1 versus a13 where the “squares” and “circles” represent the cases of a12 ¼ 0.5 and 2, respectively.
that before the cross point, a13 behaves as suppression and thus larger a12 will result in less qI1 . While after the cross point, a13 behaves as acceleration and thus larger a12 will result in more qI1 . A more comprehensive study is shown in Fig. 8, where qI1 is plotted on both the parameters a12 and a13. We see that qI1 will be zero when a12 is larger and a13 is small, decrease with the increase of a12 when a13 < 1, increase with a12 when a13 > 1, and reach saturation when both a12 and a13 are larger. Therefore, the different parameter regions in Fig. 8 represent different interactions between the two pathogens such as the cooperation, suppression and so on, implying that all of them can be unified into the same framework of Eq. (1). IV. A BRIEF THEORETICAL ANALYSIS
We now conduct a theoretical analysis to the above numerical simulations. In the framework of mean-field, the evolutionary equations of the two epidemics can be given as
FIG. 8. Dependence of qI1 on a12 and a13 with fixed a12 ¼ a22 and b1 ¼ b2 ¼ 0.25.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-6
Zhao, Zheng, and Liu
Chaos 24, 043129 (2014)
dqI1 ¼ qI1 þ a1 ðtÞb1 qI1 qS1 hki; dt dqI2 ¼ qI2 þ a2 ðtÞb2 qI2 qS2 hki; dt dqR2 ¼ qI2 ; dt
From the second and third equations of Eq. (2), we have dqS2 =dt ¼ a2 ðtÞb2 qI2 qS2 hki. In the initial condition b2 hki of Ð t qR2 ð0Þ ¼ 0 and qS2 ð0Þ 1, we have qS2 ðtÞ ¼ e 0 a2 ðtÞqI2 ðtÞdt. Thus, in the stabilized state, we have 0
qR2 ð1Þ ¼ 1 ea2 b2 hkiqR2 ð1Þ ;
(2)
with the conservation laws qS1 þ qI1 ¼ 1 and qS2 þ qI2 þqR2 ¼ 1. In the mean-field level, the coefficient a1(t) and a2(t) can be given as
where a02 is defined by 0 a2 ðtÞqI2 ðtÞdt a02 qR2 ð1Þ and can be considered as an average of a2(t). From Eq. (3) we know a2(t) > 1, which gives a02 > 1. By the condition of a0 b hkiq
dð1e 2 2 R2 dqR2 ð1Þ
a1 ðtÞ ¼ a11 qS2 ðtÞ þ a12 qI2 ðtÞ þ a13 qR2 ðtÞ; a2 ðtÞ ¼ a21 qS1 ðtÞ þ a22 qI1 ðtÞ;
qI1 ðtÞ ¼
qI1 ð0ÞðC 1ÞeðC1Þt : C 1 þ CqI1 ð0ÞeðC1Þt
(4)
As qI1 ð0Þ 0, we ignore the part of CqI1 ð0ÞeðC1Þt in the denominator and thus obtain qI1 ðtÞ qI1 ð0Þe
ðC1Þt
;
(5)
with approximately exponential growth rate (C 1). For the parameters in Fig. 2, we have C 1 ¼ a11 b1 hki 1 ¼ 0:5, which is approximately consistent with the slope c ¼ 0.48 in the right inset of Fig. 2(a). Notice from Fig. 2 that there is a time delay between the setting up of steady-state of the two diseases, i.e., the SIR model in Fig. 2(b) will reach its steady-state at about t ¼ 20, while the SIS model in Fig. 2(a) will reach its steady-state at about t > 25. Thus, the evolution of the SIS model for t > 20 is not a co-evolution of two pathogens but a single pathogen evolution. In this sense, we cannot naturally obtain the solution of steady-state of Eq. (2). Therefore, we here focus on the critical points. In the stabilized state, we have qI2 ¼ 0 and qI1 and qR2 constant. Thus, a1 ¼ a11 ð1 qR2 Þ þ a13 qR2 and a2 ¼ a21 ð1 qI1 Þ þ a22 qI1 are both constants. When b1 b1c, we have dqI1 =dt ¼ 0 in the first equation of Eq. (2), which gives qI1 ¼ 1
1 : a1 b1 hki
(6)
The critical point b1c can be obtained by letting qI1 change from zero to nonzero, which gives b1c ¼
1 : a1 hki
ð1Þ
Þ
jqR
2
ð1Þ¼0
> 1, we obtain the critical point
(3)
with a11 ¼ a21 ¼ 1. To get the solution of Eq. (2), we take the initial conditions of qS1 ð0Þ 1; qI1 ð0Þ 0; qS2 ð0Þ 1; qI2 ð0Þ 0 and qR2 ð0Þ ¼ 0. In the initial steps, we have a1(t) a11. From the first equation of Eq. (2), we have dqI1 =dt ¼ qI1 þ CqI1 ð1 qI1 Þ ¼ ðC 1ÞqI1 Cq2I1 with C ¼ a11 b1 hki. Its solution is
(7)
It will return to the threshold bc ¼ 1=hki in an isolated SIS random network.62 For the parameters in Fig. 3, we have qR2 0:57 and thus obtain a1 ¼ 0.77. Substituting a1 into Eq. (7) we have b1c ¼ 0.21, which is approximately consistent with the threshold b1c ¼ 0.2 in Fig. 3(a) and thus confirms the numerical result.
(8)
Ð1
b2c ¼
1 ; a02 hki
(9)
which depends on the first pathogen through a02 or a2(t). b2c will also return to the threshold bc ¼ 1=hki in an isolated SIR random network. For the parameters in Fig. 3, we have qI1 0:28 and thus a2 ¼ 1.28, which gives a02 1:14. Substituting it into Eq. (9) we have b2c ¼ 0.146, which is consistent with the threshold b2c ¼ 0.14 in Fig. 3(b) and thus confirms the numerical result. V. DISCUSSIONS AND CONCLUSIONS
The present study is for the interactions between the two pathogens with the SIS and SIR models, respectively, but it can be easily extended to the case between two SIS pathogens or the case between two SIR pathogens and so on. The framework of Eq. (1) can be further extended to the case of more than two pathogens. On the other hand, the considered mutual influences in this work are only for the steady status. From Fig. 2, we see that in the initial evolution process, both qI1 and qI2 are enhanced, which is one of the main features to cause crowd panicked. Thus, the control of epidemic spreading in the cases of multiple pathogens should be paid more attention to the initial evolution process and thus deserved to be studied further. In conclusions, we have presented a unified framework to study the cases of two pathogens on two networks and figured out that the key point is the time dependent infectious rate, in contrast to the previous constant infectious rate. Through this framework we find that the final infectious depends on the range of the set of parameters a12, a13, and a21 and also on the types of pathogens such as SIS or SIR. This framework thus opens a new window to study the epidemic spreading on multiple networks with multiple pathogens and may stimulate further studies. ACKNOWLEDGMENTS
This work was partially supported by the NNSF of China under Grant Nos. 11135001 and 11375066, Joriss project under Grant No. 78230050, and 973 Program under Grant No. 2013CB834100. 1
S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008). 2 R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001). 3 V. M. Eguiluz and K. Klemm, Phys. Rev. Lett. 89, 108701 (2002).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
043129-7
Zhao, Zheng, and Liu
Chaos 24, 043129 (2014)
4
32
5
33
M. E. J. Newman, Phys. Rev. E 66, 016128 (2002). M. Boguna, R. Pastor-Satorras, and A. Vespignani, Phys. Rev. Lett. 90, 028701 (2003). 6 T. Gross, C. J. D. D’Lima, and B. Blasius, Phys. Rev. Lett. 96, 208701 (2006). 7 R. Parshani, S. Carmi, and S. Havlin, Phys. Rev. Lett. 104, 258701 (2010). 8 C. Castellano and R. Pastor-Satorras, Phys. Rev. Lett. 105, 218701 (2010). 9 Y. Schwarzkopf, A. Rakos, and D. Mukamel, Phys. Rev. E 82, 036112 (2010). 10 V. Colizza, R. Pastor-Satorras, and A. Vespignani, Nat. Phys. 3, 276 (2007). 11 V. Colizza and A. Vespignani, Phys. Rev. Lett. 99, 148701 (2007). 12 V. Colizza and A. Vespignani, J. Theor. Biol. 251, 450 (2008). 13 A. Baronchelli, M. Catanzaro, and R. Pastor-Satorras, Phys. Rev. E 78, 016111 (2008). 14 M. Tang, L. Liu, and Z. Liu, Phys. Rev. E 79, 016108 (2009). 15 M. Tang, Z. Liu, and B. Li, Europhys. Lett. 87, 18005 (2009). 16 Z. Liu, Phys. Rev. E 81, 016110 (2010). 17 Z. Ruan, P. M. Hui, H. Q. Lin, and Z. Liu, Eur. Phys. J. B 86, 13 (2013). 18 Z. Ruan, M. Tang, and Z. Liu, Eur. Phys. J. B 86, 149 (2013). 19 P. Grassberger, J. Stat. Mech.: Theory Exp. 2013, P04004. 20 D. H. Zanette and S. Risau-Gusman, J. Biol. Phys. 34, 135 (2008). 21 V. Marceau, P.-A. Noel, L. Hebert-Dufresne, A. Allard, and L. J. Dube, Phys. Rev. E 82, 036116 (2010). 22 I. B. Schwartz and L. B. Shaw, Physics 3, 17 (2010). 23 C. Lagorio, M. Dickison, F. Vazquez, L. A. Braunstein, P. A. Macri, M. V. Migueles, S. Havlin, and H. E. Stanley, Phys. Rev. E 83, 026102 (2011). 24 A. Vazquez, B. Racz, A. Lukacs, and A. Barabasi, Phys. Rev. Lett. 98, 158702 (2007). 25 F. C. Coelho, O. G. Cruz, and C. T. Codeco, Source Code Biol. Med. 3, 3 (2008). 26 S. Meloni, A. Arenas, and Y. Moreno, Proc. Natl. Acad. Sci. 106, 16897 (2009). 27 D. Balcan, V. Colizza, B. Goncalves, H. Hu, J. J. Ramasco, and A. Vespignani, Proc. Natl. Acad. Sci. 106, 21484 (2009). 28 S. Funk and V. A. A. Jansen, Phys. Rev. E 81, 036118 (2010). 29 M. Dickison, S. Havlin, and H. E. Stanley, Phys. Rev. E 85, 066109 (2012). 30 A. Saumell-Mendiola, M. A. Serrano, and M. Boguna, Phys. Rev. E 86, 026106 (2012). 31 F. D. Sahneh, C. Scoglio, and P. Van Mieghem, IEEE/ACM Trans. Networking 21, 1609 (2013).
D. M. Morens, G. K. Folkers, and A. S. Fauci, Nature 430, 242 (2004). C. Granell, S. Gomez, and A. Arenas, Phys. Rev. Lett. 111, 128701 (2013). 34 G. Su, Z. Ruan, S. Guan, and Z. Liu, Europhys. Lett. 103, 48004 (2013). 35 Z. Ruan, M. Tang, and Z. Liu, Phys. Rev. E 86, 036117 (2012). 36 S. Funk, E. Gilad, C. Watkins, and V. A. A. Jansen, Proc. Natl. Acad. Sci. 106, 6872 (2009). 37 F. Bagnoli, P. Lio, and L. Sguanci, Phys. Rev. E 76, 61904 (2007). 38 Q. Wu, X. Fu, M. Small, and X. Xu, Chaos 22, 013101 (2012). 39 H. Zhang, J. Zhang, C. Zhou, M. Small, and B. Wang, New J. Phys. 12, 023015 (2010). 40 A. Perisic and C. T. Bauch, PLoS Comput. Biol. 5, e1000280 (2009). 41 K. Dietz, J. Math. Biol. 8, 291 (1979). 42 F. B. Bang, Int. J. Epidemiol. 4, 337 (1975). 43 M. E. J. Newman, Phys. Rev. Lett. 95, 108701 (2005). 44 W. L. Jeffries, Perspect. Sexual Reprod. Health 43, 151 (2011). 45 W. L. Jeffries and B. Dodge, J. Sex Res. 44, 278 (2007). 46 C. Goodenow, J. Netherland, and L. Szalacha, Am. J. Public Health 92, 203 (2002). 47 L. B. Hightow et al., Sexually Transmitted Dis. 33, 585 (2006). 48 C. T. Bauch and A. P. Galvani, Science 342, 47 (2013). 49 L. J. Abu-Raddad, P. Patnaik, and J. G. Kublin, Science 314, 1603 (2006). 50 M. Martcheva and S. S. Pilyugin, SIAM J. Appl. Math. 66, 843 (2006). 51 V. Marceau, P.-A. Noel, L. Hebert-Dufresne, A. Allard, and L. J. Dube, Phys. Rev. E 84, 026105 (2011). 52 J. C. Miller, Phys. Rev. E 87, 060801(R) (2013). 53 M. E. J. Newman and C. R. Ferrario, PLoS One 8, e71321 (2013). 54 L. Chen, F. Ghanbarnejad, W. Cai, and P. Grassberger, Europhys. Lett. 104, 50001 (2013). 55 H. L. Mills, A. Ganesh, and C. Colijn, J. Theor. Biol. 320, 47–57 (2013). 56 R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47 (2002). 57 R. M. Anderson and R. M. May, Infectious Diseases of Humans: Dynamics and Control (Oxford University Press, Oxford, New York, 1991). 58 R. Olinky and L. Stone, Phys. Rev. E 70, 030902 (2004). 59 Z. Liu and B. Hu, Europhys. Lett. 72, 315 (2005). 60 Y. Moreno, R. Pastor-Satorras, and A. Vespignani, Eur. Phys. J. B 26, 521 (2002). 61 Y. Moreno, M. Nekovee, and A. Vespignani, Phys. Rev. E 69, 055101 (2004). 62 R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 63, 066117 (2001).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.81.226.78 On: Sun, 17 May 2015 08:56:24
