† Corresponding author. E-mail:

Project supported by the National Natural Science Foundation of China (Grant Nos. 11135001, 11375066, and 11405059) and the National Basic Key Program of China (Grant No. 2013CB834100).

Epidemic spreading has been studied for a long time and is currently focused on the spreading of multiple pathogens, especially in multiplex networks. However, little attention has been paid to the case where the mutual influence between different pathogens comes from a fraction of epidemic propagators, such as bisexual people in two separated groups of heterosexual and homosexual people. We here study this topic by presenting a network model of two layers connected by impulsive links, in contrast to the persistent links in each layer. We let each layer have a distinct pathogen and their interactive infection is implemented by a fraction of propagators jumping between the corresponding pairs of nodes in the two layers. By this model we show that (i) the propagators take the key role to transmit pathogens from one layer to the other, which significantly influences the stabilized epidemics; (ii) the epidemic thresholds will be changed by the propagators; and (iii) a reverse-feeding effect can be expected when the infective rate is smaller than its threshold of isolated spreading. A theoretical analysis is presented to explain the numerical results.

Epidemic spreading has been a challenging problem for a long time as it usually causes major public health threats owing to their potential mortalities and substantial economic impacts. Its recent examples are the outbreaks of global infectious diseases including SARS (Severe Acute Respiratory Syndrome), H1N1 (Swine Influenza), H5H1 (Avian Influenza), and Ebola, etc. To understand how epidemics spread in a system, which is a crucial step to preventing and controlling outbreaks, several models have been proposed such as the susceptible-infected-susceptible (SIS) and susceptible-infected-refractory (SIR) models, etc.,^{[1]} where *S*, *I*, and *R* represent the susceptible, infected, and refractory individuals, respectively. Based on these models, it is revealed that the epidemic can be spread out only when its infection rate is greater than a threshold.

The study of epidemics has been a very hot research topic in recent decades, due to the fast progress in network science.^{[2,3]} It is revealed that the stabilized epidemic depends on many factors such as the network structure, adopted immunization, received information, and chosen traveling, etc. For example, it was found that the epidemic spreading on scale-free networks has no epidemic threshold in the thermodynamic limit.^{[4–11]} Different approaches of immunization will result in different effects of controlling an epidemic, such as the targeted immunization^{[12–14]} and acquaintance immunization,^{[15,16]} etc. Received information will make individuals take self-protective actions to protect themselves.^{[17–21]} Epidemic spreading on one layer will asymmetrically interact with the communication on another layer.^{[22,23]} Objective traveling of choosing the shortest path will accelerate the epidemic spreading,^{[24–28]} in contrast to the random diffusion model where individuals randomly choose one of its neighboring nodes for diffusion.^{[29–33]}

Currently, the study of epidemics is focusing on the multi-layered networks where a variety of dynamics can be investigated.^{[34–36]} A great deal ofattention has been paid to the case of multiple pathogens, where two or more diseases spread in the same population and each person can be infected by two or more of them.^{[37–42]} In fact, this case is more realistic. For example, it was revealed that the persons infected by HIV have a higher risk to be infected by other pathogens, including hepatitis, tuberculosis, and malaria.^{[43,44]} It was also found that the incidence of tuberculosis was increased during the 1918–1919 Spanish flu.^{[45,46]} These evidences show that infection by one pathogen can effectively influence the infection process by the other pathogen, which may be either competitive or cooperative for the same population of hosts. Focusing on the two-pathogen scenario, it has been found that the onset of a disease’s outbreak is conditioned to the prevalence levels of the other disease. For the case of competitive (cooperative) interaction between two pathogens, the presence of one pathogen makes the other less (more) likely to spread.^{[37,42,44,47–50]}

A typical but interesting case of multiple pathogens is sexually transmitted diseases, which can propagate in both heterosexual and homosexual networks of sexual contacts.^{[51,52]} This problem has been recently considered by Sanz *et al.*^{[53]} and our team.^{[54]} Sanz *et al.* thoroughly studied extensions of the SIS and SIR models, based on the condition of constant infection rates. While we presented a unified framework to study both the cooperative and competitive interactions between two pathogens in multiplex networks, based on a new concept of a state-dependent infectious rate.^{[54]} However, all these considerations ignore the effect of bisexual men, which is evidenced to increase the AIDS-related risk.^{[55,56]} A study on the number of partners among U.S. bisexual men shows that bisexual men are the medium to connect the pathogens in the network of heterosexual men with the pathogens in the network of homosexual men.^{[55–59]} Based on these evidences, we here study this problem by focusing on the role of bisexual individuals, i.e., connecting the heterosexual network and homosexual network. Thus, we present a two-layered network model without fixed links between the two layers, but with dynamic propagators to connect the two layers. This model is not only limited for sexually transmitted diseases, but also for other pathogens caused by traveling between different cities. In detail, we let each layer have a distinct pathogen and let the propagators have a chance to jump between the two layers at each time step, see Fig.

The rest of this paper is organized as follows. In Section 2, we introduce the model to study the effect of the propagators. Then in Section 3, we make numerical simulations to show how the propagators influence the spreading of the two pathogens, especially the epidemic thresholds. After that, we give a brief theoretical analysis to explain the numerical results in Section 4. Finally, we give discussions and conclusions in Section 5.

Figure ^{[2]} with sizes *N*_{1} and *N*_{2} and average degrees ⟨*k*_{1}⟩ and ⟨*k*_{2}⟩, respectively. We assume that there are totally *m* propagators in the network and they are initially half in the up layer and the other in the down layer. The propagators have a double identity. They will behave only as homosexual (heterosexual) people when they are at the homosexual (heterosexual) layer. At each time step, every propagator will have a probability *p* to jump from its layer to the other layer. To smoothly implement this jumping process, we randomly choose *m* nodes from each layer and let them be fixed as the host nodes for the propagators. When the propagators are staying at these host nodes, their connections to neighbors will be active (see the “red” nodes in Fig. *m*/2 occupied nodes for the propagators and *m*/2 “hole” nodes at each layer. We will always have *N*_{1} − *m* non-jumping nodes at the up layer and *N*_{2} − *m* non-jumping nodes at the down layer. We let a propagator from one layer and a hole node from the other layer form a corresponding pair of nodes, see those pairs of nodes connected by the dotted “red” lines between the two layers in Fig.

Then, we assume two pathogens in the network, i.e., pathogen-1 on the up layer and pathogen-2 on the down layer. Without propagators, the two pathogens will independently spread on their individual layers, respectively. With propagators, the two pathogens will spread on both the two layers. Therefore, there will be a dynamical interaction between the two pathogens once they spread to the same nodes. To understand how this dynamical interaction influences epidemic spreading, we here take the susceptible-infected-susceptible (SIS) model for both the two pathogens as an example. In the SIS model, a susceptible node may be infected by an infected neighbor at rate *β* when it contacts one infected node. At the same time, each infected node will become susceptible again at rate *μ* at each time step.^{[1,4,5,7,60–63]} In the case of complex networks, a susceptible node may have *k*_{inf} infected neighbors, thus it will become infected with probability 1 − (1 − *β*)^{kinf}. While the present case is more complicated as it has two interacting pathogens. To describe it, we let the set (*S*_{1}, *I*_{1}) represent pathogen-1 and the set (*S*_{2}, *I*_{2}) represent pathogen-2. In this case, a node will have 4 states: *S*_{1}*S*_{2} — the node is susceptible to both pathogen-1 and pathogen-2; *I*_{1}*S*_{2} — the node is infected by the first pathogen-1 but susceptible to pathogen-2; *S*_{1}*I*_{2} — the node is susceptible to the first pathogen-1 but infected by the second pathogen-2; and *I*_{1}*I*_{2} — the node is infected by both pathogen-1 and pathogen-2.^{[53,54]}

Figure *S*_{1}*S*_{2} and *I*_{1}*I*_{2} is ignored as a person will usually not be infected by two pathogens at exactly the same time. Let *β*_{1} (*β*_{2}) and *μ*_{1} (*μ*_{2}) be the infective and recovery rates of pathogen-1 (pathogen-2), respectively. Then, an *S*_{1}*S*_{2} node will be infected into an *I*_{1}*S*_{2} status by a probability 1 − (1 − *β*_{1})^{kinf1}, where *k*_{inf1} is the number of its total infected neighbors with either status *I*_{1}*S*_{2} or status *I*_{1}*I*_{2}. At the same time, each *I*_{1}*S*_{2} node will recover to the *S*_{1}*S*_{2} status by a probability *μ*_{1}. Similarly, an *S*_{1}*S*_{2} node will become the *S*_{1}*I*_{2} status by a probability 1 − (1 − *β*_{2})^{kinf2}, where *k*_{inf2} is the number of its total infected neighbors with either status *S*_{1}*I*_{2} or status *I*_{1}*I*_{2}. At the same time, each *S*_{1}*I*_{2} node will recover to the *S*_{1}*S*_{2} status by a probability *μ*_{2}.

For an infected node, its secondary infection by another pathogen will be different from its first infection. It has been revealed that an infection by malaria is shown to increase an individual’s attractiveness to mosquitos,^{[64]} which in turn may increase the chance to be infected by other strains. Generally, 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. On the other hand, when an infected person takes medicine, he will be less likely to get the other, which corresponds to the case of competition between two pathogens. We here introduce a factor to reflect this kind of cooperation or competition. In detail, we let *f*_{1} be the infection factor and *f*_{2} be the recovery factor. Thus, an *I*_{1}*S*_{2} node will become the *I*_{1}*I*_{2} status by a probability 1 − (1 − *f*_{1}*β*_{2})^{kinf2}. At the same time, each *I*_{1}*I*_{2} node will recover to the *S*_{1}*I*_{2} status by a probability *f*_{2}*μ*_{2}. Similarly, an *S*_{1}*I*_{2} node will become the *I*_{1}*I*_{2} status by a probability 1 − (1 − *f*_{1}*β*_{1})^{kinf1}. At the same time, each *I*_{1}*I*_{2} node will recover to the *S*_{1}*I*_{2} status by a probability *f*_{2}*μ*_{1}.

We are here mainly focused on how the impulsive jumping between the two networks influences the epidemic spreading, as it is what happens in reality and has not been paid attention so far. It is well known that the network topology will seriously influence the epidemic spreading.^{[2,3]} To reduce the influence from the aspect of network topology, we limit the jumping probability to a very small range, i.e., *p* ≪ 1. In this condition, the main influence between the two layers will come from the effect of jumping. We will use the random Erdös–Rényi networks as an example, as it is a good candidate for the mean-field approach. We find that the obtained results also work for other networks such as the scale-free networks.

In numerical simulations, we first construct two independent random Erdös–Rényi networks^{[2]} as the up and down layers of Fig. *N*_{1} = 1000 nodes and the down network have *N*_{2} = 4000 nodes. To reflect the different ways of communication in the homosexual and heterosexual communities, we let the average degrees of the up network be ⟨*k*_{10}⟩ = 8 and that of the down network be ⟨*k*_{20}⟩ = 3. We let the number of total propagators be *m* = 400, i.e., there are initially *m*/2 = 200 propagators in each of the up and down layers. The non-jumping nodes will be 600 in the up layer and 3600 in the down layer. When the jumping probability *p* is zero, the up and down networks become independent. Notice that the links of the hole nodes in Fig. *k*_{10}⟩ *m*/2 inactive links in the up network and ⟨*k*_{20}⟩ *m*/2 inactive links in the down network. In general, it is possible for the hole nodes to also be neighbors of each other, indicating that the average degree of the up layer ⟨*k*_{1}⟩ will be less than ⟨*k*_{10}⟩ but greater than (*N*_{1} ⟨*k*_{10}⟩ − 2 ⟨*k*_{10}⟩ *m*/2)/(*N*_{1}−*m*/2) = ⟨*k*_{10}⟩ (*N*_{1} − *m*)/(*N*_{1} − *m*/2). Similarly, the average degree of the down layer ⟨*k*_{2}⟩ will be less than ⟨*k*_{20}⟩ but greater than ⟨*k*_{20}⟩ (*N*_{2} − *m*)/(*N*_{2} − *m*/2). When the hole nodes are non-neighbors, we have ⟨*k*_{1}⟩ = ⟨*k*_{10}⟩ (*N*_{1} − *m*)/(*N*_{1} − *m*/2) = 6 and ⟨*k*_{2}⟩ = ⟨*k*_{20}⟩ (*N*_{2} − *m*)/(*N*_{2} − *m*/2) ≈ 2.84. For a homogeneous random network, we have the critical infection rate *β*_{c} = *μ*/⟨*k*⟩.^{[65,66]} Taking *μ*_{1} = *μ*_{2} = 0.01 as an example, we will have the thresholds *β*_{1c0} ≈ *μ*/⟨*k*_{1}⟩ ≈ 0.00167 for the up layer and *β*_{2c0} ≈ *μ*/⟨*k*_{2}⟩ ≈ 0.00352 for the down layer when the two layers are isolated with *p* = 0.

We now consider epidemic spreading by fixing *μ*_{1} = *μ*_{2} = 0.01 in this paper. We choose the factors *f*_{1} = 1.4 and *f*_{2} = 1, unless otherwise specified. Initially, we randomly choose a few nodes of the up (down) layer of Fig. *ρ*_{I}, which includes both the first infection *I*_{1}*S*_{2} or *S*_{1}*I*_{2} and the secondary infection *I*_{1}*I*_{2}. We let *ρ*_{I} on the jumping probability *p* with the “squares”, “circles”, “uptriangles”, and “downtriangles” being *β*_{1} = 0.003 and *β*_{2} = 0.005, (b) the case of *β*_{1} = 0.003 and *β*_{2} = 0.003, (c) the case of *β*_{1} = 0.0014 and *β*_{2} = 0.005, and (d) the case of *β*_{1} = 0.0014 and *β*_{2} = 0.003. The parameters in Fig. *β*_{1} > *β*_{1c0} and *β*_{2} > *β*_{2c0} in panel (a), *β*_{1} > *β*_{1c0} and *β*_{2} < *β*_{2c0} in panel (b), *β*_{1} < *β*_{1c0} and *β*_{2} > *β*_{2c0} in panel (c), and *β*_{1} < *β*_{1c0} and *β*_{1c0} < *β*_{2} < *β*_{2c0} in panel (d). A common point in the four panels of Fig. *p* = 0, there is at most one pathogen of *ρ*_{I} > 0 for each layer, confirming the independence of the two layers at *p* = 0.

The chosen values of *β*_{1} > *β*_{1c0} and *β*_{2} > *β*_{2c0} in Fig. *p* = 0. When *p* > 0, the infected seeds of pathogen-1 (pathogen-2) will be transmitted to the down (up) layer and then cause its spreading there. In the up layer of the case of Fig. *β*_{1} and *β*_{2} are greater than *β*_{1c0} and thus they can both spread out, provided that the propagators can transmit the infected seeds of pathogen-2 from the down layer to the up layer. The final *β*_{1} − *β*_{1c0} and *β*_{2} − *β*_{1c0}, which results in their approximate constant values in Fig. *β*_{1} < *β*_{2c0}, thus the seeds of pathogen-1 cannot spread out but can only stay there for a transient time. As the recovery rate *μ*_{2} = 0.01 < 1, a transmitted seed of pathogen-1 will keep its status for some time steps before it jumps back to the up layer. In this sense, *p* > 0 will provide persistent jumping seeds of pathogen-1 into the down layer and thus the value of *p*. On the other hand, for an individual seed, its effect of epidemic spreading depends not only on the jumping speed but also on the contacting time,^{[25]} resulting in a nonlinear dependence of *p*, as in Fig.

We now turn to Fig. *β*_{1} < *β*_{1c0} and *β*_{2} < *β*_{2c0}. When *p* = 0, there are neither nonzero pathogen-1 nor nonzero pathogen-2 in the system. However, *p* > 0 induces an interesting phenomenon that *β*_{2} > *β*_{1c0}, the transmitted seeds of pathogen-2 will spread out on the up layer. Then, the propagators will send pathogen-2 back to the down layer and thus sustain the nonzero

To understand the underlying mechanism of the nonlinear dependence of *ρ*_{I} on *p* in the four panels of Fig. *β*_{2} = 0, where panel (a) represents the case of *β*_{1} = 0.003 > *β*_{1c0} and panel (b) the case of *β*_{1} = 0.0014 < *β*_{1c0}. From Fig. *p*, *β*_{1} = 0.003 < *β*_{2c0}, the epidemic cannot spread out at the down layer, indicating that *p* will result in an increase of *β*_{1} = 0.0014 being smaller than both *β*_{1c0} and *β*_{2c0}.

Figure *β*_{1} = 0, where panel (c) represents the case of *β*_{1} = 0.005 > *β*_{2c0} and (d) the case of *β*_{2} = 0.003 < *β*_{2c0}. From Fig. *p*, both *β*_{2} = 0.005 is greater than both *β*_{2c0} and *β*_{1c0}, the epidemic can spread out at both the up and down layers. In this situation, *β*_{2}–*β*_{1c0} and *β*_{2} − *β*_{2c0}. Thus, the influence of *p* can be ignored. Figure *p* but with *β*_{2} = 0.003 is smaller than *β*_{2c0} but greater than *β*_{1c0}. The condition of *β*_{2} < *β*_{2c0} will make the epidemic not be spread out, indicating no epidemic at both the up and down layers when *p* = 0. However, when *p* > 0, it is possible for pathogen-2 to be transmitted to the up layer by the propagators. Once it happens, pathogen-2 will spread out at the up layer because *β*_{2} > *β*_{1c0}. Then the propagator will bring the infected ones back to the down layer and thus result in a sustained *β*_{2c0} and *β*_{1c0} caused by the different average degrees between the up and down layers.

Secondly, we consider the aspect of mutual influence of two pathogens, where *f*_{1} and *f*_{2} are the two key parameters to determine the feature of interaction, i.e. cooperation or competition. To clearly see its effect, we focus on the influence of *f*_{1} and *f*_{2} on the epidemic thresholds. We let *β*_{1c} and *β*_{2c} be the epidemic thresholds of the up and down layers with *p* > 0, respectively. We first check the effect of *f*_{1} by fixing *f*_{2} = 1 and *p* = 0.005. Figure *β*_{1c} depends on *f*_{1} where the “squares” and “uptriangles” represent the cases of *β*_{2} = 0.005 and 0.003, respectively. From Fig. *β*_{1c} decreases with the increase of *f*_{1} for the two cases but is faster in the case of *β*_{2} = 0.005 than that of *β*_{2} = 0.003. This result can be understood as follows. As *β*_{2} = 0.005 is greater than both *β*_{1c0} and *β*_{2c0}, pathogen-2 can spread out in both the up and down layers, implying an interaction between the two pathogens in the up layer. Notice that an increase of *f*_{1} means the enhanced cooperation between the two pathogens. Therefore, the increase of *f*_{1} will result in an enhanced spreading of *β*_{1c} on *f*_{1}. While in the case of *β*_{2} = 0.003, its *β*_{2} = 0.005, indicating its interaction with *β*_{2} = 0.005 and thus will result in a slower decreased *β*_{1c} on *f*_{1}.

Figure *β*_{2c} depends on *f*_{1} where the “squares” and “uptriangles” represent the cases of *β*_{1} = 0.003 and 0.0014, respectively. From Fig. *β*_{2c} decreases with the increase of *f*_{1} for the case of *β*_{1} = 0.003 but keeps a constant for the case of *β*_{1} = 0.0014. This result is easy to understand. As *β*_{1} = 0.003 is greater than *β*_{1c0} but smaller than *β*_{2c0}, there will be a small quantity of *β*_{2c} on *f*_{1}. While *β*_{1} = 0.0014 is smaller than both *β*_{1c0} and *β*_{2c0}, there will be no *β*_{2c}.

We now turn to check the effect of *f*_{2} by fixing *f*_{1} = 1 and *p* = 0.005. Figure *β*_{1c} depends on *f*_{2} where the “squares” and “uptriangles” represent the cases of *β*_{2} = 0.005 and 0.003, respectively. Figure *β*_{2c} depends on *f*_{2} where the “squares” and “uptriangles” represent the cases of *β*_{1} = 0.003 and 0.0014, respectively. Comparing Figs. *f*_{1} and *f*_{2} is the opposite, confirming the role of the infection factor *f*_{1} and the role of the recovery factor *f*_{2}.

Finally, we discuss how the parameter *m* influences the thresholds *β*_{1c} and *β*_{2c}. We fix the total numbers of nodes *N*_{1} and *N*_{2} and let *p* = 0.005, *f*_{1} = 1.4, and *f*_{2} = 1. When *m* increases, the non-jumping nodes in both the up and down layers will decrease by *N*_{1} − *m* and *N*_{2} − *m*, respectively, which results in the decrease of ⟨*k*_{1}⟩ and ⟨*k*_{2}⟩ and then the increase of *β*_{1c} and *β*_{2c}. Figure *β*_{2} = 0.005 and 0.003, respectively, and the “squares” and “uptriangles” in panel (b) represent the cases of *β*_{1} = 0.003 and 0.0014, respectively. From Fig. *β*_{2} has a stronger influence on the threshold *β*_{1c}. We see the same result in Fig. *β*_{1} also has a stronger influence on the threshold *β*_{2c}.

We now conduct a theoretical analysis on the above numerical simulations. For convenience, we introduce

We divide the dynamics process into two steps, i.e. reaction and jumping. For the former, we use the mean-field approach to figure out the densities *ρ*_{SS}, *ρ*_{IS}, *ρ*_{SI}, and *ρ*_{II}. According to Fig.

*t*. Obviously, they satisfy the conservation laws

In the jumping process, we have

*N*

^{u}=

*N*

_{1}−

*m*/2 and

*α*

_{1}=

*mp*/2

*N*

^{u}represents the density variation caused by the jumping propagators. Similarly, we have

*N*

^{d}=

*N*

_{2}−

*m*/2 and

*α*

_{2}=

*mp*/2

*N*

^{d}. When the system arrives at a stationary state, we have

*X*= SS, IS, SI, and II, respectively, and

The nonlinear algebraic equations (*β*_{2} = 0. From Eq. (

*β*

_{2}= 0 we have

*β*

_{1c}=

*μ*

_{1}/⟨

*k*

_{1}⟩, which is exactly the threshold in an isolated SIS random network.

^{[66]}

For a general solution of the nonlinear algebraic equations (*k*_{1}⟩ and ⟨*k*_{2}⟩ by counting the total links of the active nodes and making average on time. Then, we substitute them into Eqs. (*β*_{1c} and *β*_{2c}. The “circles” and “down triangles” in Fig. *β*_{1c} and *β*_{2c} for different *m*, see the “circles” and “downtriangles” in Figs.

The present study is for the interaction between two SIS pathogens, but it can be easily extended to other cases such as two SIR pathogens, one SIS and one SIR pathogens, or even more than two pathogens. Further, our model can also be extended to the case of traveling between different cities where the travelers take the role of propagators. Although this work is the primary try to consider the role of propagators, its further studies may provide guiding to human mobility patterns and thus may help public health services understand and control human infectious diseases.

In conclusion, we have presented a model of two interactive pathogens in two-layered networks where the connection between the two layers is implemented by a small fraction of propagators. We interestingly find a reverse-feeding phenomenon that the stabilized epidemic spreading on one layer can be sustained by the persistent jumping of propagators from the other layer. When focusing on the thresholds with *p* > 0, we find that the jumping and the mutual interaction between the two pathogens will reduce both *β*_{c1} and *β*_{c2}. We further show that the reduced epidemic thresholds come mainly from two aspects. One is the different average degrees between the two layers. The other is the modified immune system of the secondary infection, which can be characterized by the factors of *f*_{1} and *f*_{2}. A mean-field theory is provided to explain the underlying mechanism of the epidemic spreading in the model. This work provides significant implications for controlling epidemic spreading in realistic cases.

**Reference**

1 | |

2 | Rev. Mod. Phys. 74 47 |

3 | Rev. Mod. Phys. 80 1275 |

4 | Phys. Rev. Lett. 86 3200 |

5 | Phys. Rev. Lett. 89 108701 |

6 | Phys. Rev. E 66 016128 |

7 | Phys. Rev. Lett. 90 028701 |

8 | Phys. Rev. Lett. 96 208701 |

9 | Phys. Rev. Lett. 104 258701 |

10 | Phys. Rev. Lett. 105 218701 |

11 | Phys. Rev. E 82 036112 |

12 | Phys. Rev. E 65 036104 |

13 | Phys. Rev. E 65 055103 |

14 | Phys. Rev. E 67 031911 |

15 | Phys. Rev. Lett. 91 247901 |

16 | |

17 | |

18 | |

19 | Physics 3 17 |

20 | Phys. Rev. E 86 036117 |

21 | Chin. Phys. B 23 118904 |

22 | |

23 | |

24 | Europhys. Lett. 87 18005 |

25 | Phys. Rev. E 81 016110 |

26 | Eur. Phys. J. B 86 13 |

27 | Eur. Phys. J. B 86 149 |

28 | J. Stat. Mech. 2013 04004 |

29 | Nat. Phys. 3 276 |

30 | Phys. Rev. Lett. 99 148701 |

31 | J. Theor. Biol. 251 450 |

32 | Phys. Rev. E 78 016111 |

33 | Phys. Rev. E 79 016108 |

34 | Phys. Rev. E 89 052813 |

35 | Chaos, Solitons and Fractals 80 7 |

36 | Chaos, Solitons and Fractals 68 72 |

37 | Europhys. Lett. 104 50001 |

38 | J. Theor. Biol. 320 47 |

39 | Phys. Rev. E 81 036118 |

40 | Phys. Rev. E 85 066109 |

41 | Phys. Rev. E 86 026106 |

42 | Phys. Rev. Lett. 95 108701 |

43 | |

44 | Science 314 1603 |

45 | Emerg. Infect. Dis. 14 1193 |

46 | Comput. Math. Methods Med. 2012 124861 |

47 | SIAM J. Appl. Math. 66 843 |

48 | Phys. Rev. E 84 026105 |

49 | Phys. Rev. E 87 060801 |

50 | PLoS ONE 8 e71321 |

51 | Theor. Popul. Biol. 72 389 |

52 | Nature 411 907 |

53 | |

54 | Chaos 24 043129 |

55 | Amer. J. Public Health 92 203 |

56 | Sexu. Trans. Diseases 33 585 |

57 | Perspect. Sex. Reprod Health 43 151 |

58 | J. Sex Research 44 278 |

59 | Science 342 47 |

60 | Phys. Rev. E 70 030902 |

61 | Europhys. Lett. 72 315 |

62 | Chaos 22 013101 |

63 | |

64 | PLoS Biol. 3 e298 |

65 | Phys. Rev. E 86 026115 |

66 | Phys. Rev. E 63 066117 |