We investigate the effects of cooperation between two interacting infectious diseases that spread and stabilize in a host population. We propose a model in which individuals that are infected with one disease are more likely to acquire the second disease, both diseases following the susceptibleinfected-susceptible reaction scheme. We analyze cooperative coinfection in stochastic network models as well as the idealized, well-mixed mean field system and show that cooperative mechanisms dramatically change the nature of phase transitions compared to single disease dynamics. We show that, generically, cooperative coinfection exhibits discontinuous transitions from the disease free to high prevalence state when a critical transmission rate is crossed. Furthermore, cooperative coinfection exhibits two distinct critical points, one for outbreaks the second one for eradication that can be substantially lower. This implies that cooperative coinfection exhibits hysteresis in its response to changing effective transmission rates or equivalently the basic reproduction number. We compute these critical parameters as a function of a cooperativity coefficient in the well-mixed mean field system. We finally investigate a spatially extended version of the model and show that cooperative interactions between diseases change the general wave propagation properties of conventional spreading phenomena of single diseases. The presented work may serve as a starting and reference point for a more comprehensive understanding of interacting diseases that spread in populations. 
The proliferation and prevalence of infectious diseases continue to pose a significant threat to human society [1] . The combination of increasing population densities and global mobility facilitates the emergence of new pathogens that can quickly spread on a global scale, exemplified by outbreaks in the last decade, e.g. the 2003 SARS epidemic [2, 3] , the 2009 H1N1 pandemic [4, 5] , the emergence of the novel MERS-CoV in 2013 [6] , and the 2014 Ebola crisis [7] .
Mathematical, dynamical, statistical, and computational models have become an essential tool for understanding disease dynamics and for designing effective mitigation strategies [8] . Dynamical models cover a broad scope in terms of complexity, ranging from parsimonious and conceptual models [9, 10] , network models that account for population structure and/or host mobility patterns [11] [12] [13] [14] [15] [16] [17] , to sophisticated, large-scale computational models [18, 19] that incorporate high resolution data on multi-scale transportation, demographics, epidemiological factors, and behavioral response rules. State-of-the-art computational models have become remarkably successful in reproducing observed patterns * Electronic address: ChenL@rki.de † Electronic address: fakhteh@pks.mpg.de ‡ Electronic address: dirk.brockmann@hu-berlin.de and predicting the trend of ongoing epidemics. The majority of epidemic models are designed to capture epidemiological dynamics of single pathogen species in a host population, characterized by epidemiological parameters associated with a single virus or bacterium. However, bacterial and viral pathogens generically interact, either directly or indirectly by means of their effects on the host. Prominent examples are HIV infections that through immune suppression increase the host's susceptibility towards other infections [20] , including hepatitis [21] , malaria [22] , syphilis, herpes virus, tuberculosis [23, 24] and more. This phenomenon, i.e. being infected with more than one pathogen is known as syndemics or coinfection [25] . Although, coinfections have been well documented and received considerable attention by epidemiologist, the dynamical implication of disease interactions in epidemic progression is unclear and the generic dynamical behavior is poorly understood.
Although some models have been developed recently to describe coinfection dynamics, the majority of them consider competitive interactions between diseases [26] [27] [28] [29] [30] [31] [32] , i.e. infection with one disease induces a host's cross immunity towards another infection. In these systems different phases such as extinction, coexistence, and onestrain dominance have been identified. In one class of models coinfection is captured by a multiplex network approach [33] [34] [35] , where each infectious agent is transmitted along a different sets of links in the same population [26] [27] [28] . Recently, researchers also considered cooperative pathogen interactions [36] [37] [38] 40] tion with one disease facilitating infection with another pathogen and vice versa. In Ref. [38] , a simple coinfection model of susceptible-infected-recovered (SIR) type was investigated within the framework of cooperative bond percolation. This model exhibits avalanche-like outbreak scenarios, depending on the level of cooperation and the structure of the underlying transmission network. Furthermore, it has been found that highly clustered structures in population aid the proliferation of coinfections, contrary to the effect observed in single disease dynamics [40] . Nevertheless, some fundamental issues remain elusive: What basic dynamical features can we expect in multiple-epidemic spreading with disease-disease interactions? To what extent does coinfection change the classic outbreak scenario of single diseases? How does coinfection dynamics shape the generic features of spatial propagation?
Here we introduce and investigate a straightforward network model for the contagion process of two interacting epidemics (labeled A and B) in which the interaction is explicitly incorporated. The kinetics of both diseases follows the generic susceptible-infected-susceptible (SIS) reaction scheme, in which infected individuals transmit the disease to susceptibles at a specified rate. Infected individuals remain infectious for a typical infectious period before they become susceptible again. Based on this model, we show that cooperative coinfections possess a variety of important features that are absent in single disease dynamics. For sufficiently strong cooperation, increasing the infection rate of both diseases yields abrupt, discontinuous outbreak transitions. Furthermore, cooperative coinfection generically exhibits dynamical hysteresis, a consequence of different threshold for disease outbreak and disease eradication. Finally we show that in spatially extended systems, wavefront velocities are considerably accelerated in certain parameter regimes, while the propagation wave can be considerably slowed down when the system approaches the eradication threshold, even progress in a backwards direction.
We consider a stochastic network in which each node represents a host individual, e.g. a person, in a population and links are potential transmission channels. In the conventional SIS system each node is either susceptible or infected. If the transmission probability across a link is p, a susceptible node can acquire an infection from one of its m infected neighbors with probability 1−(1−p) m . An infected node recovers from infection and becomes susceptible again with probability r. The coinfection scenario is modeled analogously. Having two types of infections, we denote the infected states by one of the diseases by A and B, respectively and coinfection by AB (see Fig. 1 ). We assume that susceptibles are infected at probability p, whereas individuals in class A and B are infected with the complementary disease at probability q. Cooperative coinfection corresponds to a scenario with q > p, i.e. being infected by one pathogen increases the infection probability with the other. In contrast, the condition q < p corresponds to antagonistic effects. When q = p the system corresponds to two non-interacting infections. For simplicity we consider identical recovery probabilities (r = 0.5) for both infections [41] .
To understand the general behavior of the system and identify the nature of phase transitions, we investigated the dynamics on 2d square lattices and random networks while monitoring the fraction of non-susceptible nodes, 1 − ρ S in equilibrium, i.e. the fraction of nodes that are either infected by A or B or both, as a function of transmission rate parameters p and q. 1−ρ S can be considered the order parameter of the system. At fixed q, we distinguish forward (i.e. increasing p) and backward (i.e. decreasing p) direction, yielding an outbreak phase transition (OPT) and an eradication phase transition (EPT), respectively. These transitions are identical in the ordinary SIS model. To initiate an epidemic, two nodes are randomly chosen and seeded with A and B, respectively.
We first study the system on 2d square lattice. In the limit of strong cooperation [42] p q = 1, we find a discontinuous OPT. When p increases beyond the critical point p o c ≈ 0.181, the system exhibits a discontinuous jump to a large endemic state, see Fig. 2a , contrary to the single disease case (conventional SIS dynamics with a threshold p c ). However, the critical point is identical in both systems (i.e., p o c ≈ p c ). This is understandable, as for A and B seeds to generate cooperative effects, each single disease must be sustained for some time in order for single disease fronts to meet and initiate cooperative effects. Once the two wave fronts touch, the amplified coinfection establishes and acquires access to a consider-able fraction of the population yielding a discontinuous phase transition [43] . This type of behavior was also observed in cooperative bond percolation [38] .
In addition to the discontinuous OPT we find that eradication of an endemic state occurs at a substantially reduced threshold p e c (p e c < p o c ) in the cooperative regime [43] . This means that even a very small basic transmission probability p is sufficient to sustain the coinfection in the population. Note that although the transition is very steep, it is continuous in the 2d-lattice. Cooperative coinfection therefore exhibits a clear hysteresis in its order parameter. Once a coinfection outbreak has established an endemic state in the population, a large effort is required to eradicate the cooperating diseases. The magnitude of this hysteresis effect increases with increasing degree of cooperation, i.e. increasing q, see Fig. 2c . The hysteresis loop disappears if q < q c ≈ 0.167, where these two thresholds coincide.
Qualitatively, we observe the same behavior in Erdős-Réyni (ER) random networks (Figure 2b and 2d) . A prominent difference, however, is that here the eradication phase transition is also discontinuous. This implies that the infection suddenly disappears as the eradication threshold is crossed from above. This different behavior at the EPT in the ER-network is due to the net- work's lack of short loops that are present in the lattice and that are required to sustain coinfection. In the ERnetwork, single infection pertrusions have to move along much longer paths to meet and incite a coinfection spark. Fig. 3a and 3b show typical endemic patterns in the 2dlattice system for independent, non-cooperative (p = q) and cooperative (q = 1) parameter values, both being slightly beyond the critical point. In the first case (Fig. 3a) , the pattern is loosely distributed without any significant spatial structure, with a small fraction of AB infected sites. In contrast, in the cooperative case (Fig. 3b) , the infected clusters exhibit spatial structure, in the sense that the backbones of the infected islands always consist of essential AB individuals that sustain the coinfection. This morphological difference is further evidence that in the cooperative case, it is the disease interaction and a significant fraction of coinfected AB individuals that support the survival in the regime p < p c (see Fig. 2 ). This is evident in Fig. 3c that depicts the fraction individuals that are infected by both diseases normalized by the fraction of infected individuals of all types, i.e. ρ AB /(ρ A +ρ B +ρ AB ). Without cooperation this fraction continuously decreases until the critical points is reached. When cooperation is acting, the fraction of AB individuals remains almost constant until the eradication threshold is reached below which the endemic state col-lapses.
To better understand the nature of phase transitions, we provide the full phase diagrams for the two topologies in Fig. 4a and 4b . In either case, the outbreak transition point is at p o c = p c , irrespective to q. Tricritical points [44] are identified, which separate the discontinuous transition (for q > q c ) and the continuous transition (for q < q c ). Coincidentally, this tricritical point of outbreak transition is also the tricritical point of the eradication phase transition for ER networks, but not for the 2d lattice, since its eradication phase transitions are all continuous.
Apart from the outbreak size a characteristic feature of the stochastic dynamics is the outbreak probability P , i.e. the chance that an outbreak occurs in response to a small infected seed. In the coinfection dynamical system we distinguish two outbreak probabilities. First we define P (A ∪ B) as the probability that an initial seed yields an endemic state for either single infection A, B or coinfection AB. We also define P (A ∩ B) as for generating coinfection endemic where both diseases must survive. The results are shown in Fig. 4c and 4d . Although the outbreak size transition is discontinuous, the two probability transitions are continuous and occur at the same threshold p c . Note that, ER network and 2dlattice exhibit clearly different P (A ∩ B) profiles. In the random network P (A ∩ B) increases more slowly with p although both topologies have identical average degree. Again, a plausible explanation here is that the lattice system exhibits a higher likelihood of short loops that facilitate the outbreak of coinfection.
The similarities as well as differences in the nature of phase transitions in the lattice system on one hand and the ER-network on the other hand pose the question which properties of the phase transitions are observed in a well mixed particle kinetic system that is unaffected by network structure and fluctuations in the system. To address this question we analyzed the mean field dynamical system that corresponds to cooperative SIS dynamics. The basic reaction scheme underlying the dynamical process is given by as a function of infection rate α for different levels of cooperativity C. Solid and dashed lines correspond to stable and unstable stationary states, respectively. (b) Typical bifurcation types for continuous (C = 1), tri-critical (C = 2), and discontinuous (C = 5) transitions. All show transcritical bifurcation at αc except for tri-critical case, which is of pitchfork type instead. Furthermore, a saddle-node bifurcation precedes at α e c in the discontinuous case. As a consequence, a bistable region is found for the case with C > 2 , which is bounded by αc and α e c , corresponding to critical outbreak and eradication transition rate, respectively.
where S, I A , I B , I AB correspond to the fraction of S, A, B and AB individuals in the population. As in the previous section we assume for simplicity that the diseases are characterized by identical parameters. This assumption can be relaxed without changing the qualitative behavior but eases the analysis. We can simplify the system by assuming that I A (t) = I B (t) which follows for symmetric initial conditions. However, even for nonsymmetric conditions the system approaches this symmetric state. Defining new variables P = I A = I B and X = I A + I AB = I B + I AB , the above system simplifies tȯ
where we have absorbed the recovery rate γ into time, other rates are expressed in units of γ. We can conveniently define C = β/α as the cooperativity coefficient.
The epidemic interaction is cooperative if C > 1 and antagonistic for C < 1. The special cases with C = 1 corresponds to non-interacting diseases (equivalent to p = q in Sec. II). Figure 5 depicts the asymptotic behavior of the system for different levels of cooperativity C. A bifurcation analysis shows that for weak cooperation (C < 2), the system exhibits a continuous transition, at α c = 1, the same threshold as the single SIS case. For C > 2, a bistable region emerges, with a different eradication threshold α e c = 2 √ C − 1/C. Hence, a hysteresis loop is formed on the interval (α e c , α c ). Both bifurcations at α c are transcritical. The key bifurcation is the additional saddle-node bifurcation at α e c that generates a bistable region. The case of C = 2 has the critical level of cooperation which connects continuous and discontinuous bifurcations yielding a pitchfork bifurcation. This analysis also implies that in the bistable regime, an unstable stationary state exists between both stable states and that only for sufficiently substantial perturbations above this unstable stationary state will invoke a non-trivial endemic state. Overall, these results qualitatively are in good agreement with what we observed in numerical experiments in Sec. II, especially in the ER-network system, which is naturally more similar to the well-mixed population structure as required by the mean field treatment.
The defining feature of the bistable regime in cooperative coinfection is the fact that two stable stationary states co-exist: i) the disease free state and ii) the macroscopic endemic state, separated by an unstable intermediate stationary state. Because of this, a non-infinitesimal perturbation of the disease free state is required to push the system into the attractor of the endemic state. Therefore, it is reasonable to assume that the standard properties of wave propagations of pulled-fronts (diffusion into an unstable state) in a spatially extended system is dif- ferent for cooperative coinfection. In order to address these question we investigate the generic generalization of the above model to include spatial coordinates and allow individuals to diffuse in space. For simplicity, we first consider the corresponding reaction-diffusion system in one dimension:
where D in last terms is the diffusion coefficient. We assume that diffusion is constant and independent of the state of a host individual. Here we use also the symmetric parameters for the two diseases and set γ = 1. For sufficiently localized initial conditions (and only one disease) the above system exhibits a propagating wave asymptotically with a wave front velocity of V = 2 αD(1 − 1/α) for α > 1, these types of wave fronts are usually referred to as "pulled fronts" [46] . Interestingly, very different propagation modes arise in the cooperative coinfection system (red solid line in Fig. 6a ). The first observation is that the velocity is larger than the single disease case for a given α, which can be explained by cooperative effect and is less surprising. However, a negative wavefront velocity V is found for values of α e c < α < α st , equivalent to a wave that is propagating backwards and shrinking the coinfection bulk. Coinfec-tion cannot be sustained in this range of α. Note however, that at these values of α the diffusion free system possesses a stable coinfection state. Therefore, the diffusion in this regime dilutes the coinfected individuals to a level where the endemic state vanishes. This can only occur, however, in a system that requires a finite threshold (in coinfection the unstable stationary state) above which a concentration needs to grow in order to reach the stable endemic state. The effect becomes more prominent when α → α e c , because the attractive basin of the non-outbreak branch expands to substantial fraction of phase space. A marginal scenario is observed when diffusion and reaction forces balance and yield a stable "standing wave" (V = 0) that we also observe in numerical experiments. We denote this point by α st , which marks the separation of the two moving directions. Figure 7 illustrating the three typical propagation modes in 2d space are depicted (4)). Bottom panels (i-l): backward propagation (α = 0.61). The infected fraction IA+IB +IAB is color coded. The sequence (from left to right) of panels depicts the time course of the infected regions at time t = 0, 100, 200, 300, respectively. Here αst is slightly different from the value in 1d space, up to the dimension correction. Initial conditions start from some randomly infected round regions with random radius as shown in the first column. Other parameters: C = 10, D = 1.
(see the movies in the Supplemental Material [47] ).
One situation of particularly practical interest is when a disease invades a region that has already been infected by the other one. This happens when the two diseases enter a region at different time so that the secondary infection acts in an established region by the first infection wave. When this happens, the spread of the secondary disease (blue dashed-dotted) is found not only much faster than its usual speed of single disease (black dashed), but is also faster than the wave of coinfection propagation (red solid).
We presented a class of structurally parsimonious models of two cooperating epidemics that spread within the same structured population. We have shown that new dynamical properties emerge as a consequence of cooperative coinfection that are absent in single disease dynamics or antagonistic coinfection dynamics. The most prominent property is the existence of discontin-uous transitions to an endemic state when the disease specific outbreak threshold is crossed and that outbreak and eradication transitions occur at different critical parameters. This is an important insight for designing eradication protocols for pathogens that are known to interact with other diseases in a cooperative way. Because of the universality of this property of intrinsic hysteresis, we believe that with increasing data on interactions between pathogens in traditional epidemiology the simple single threshold concept may have to be revised and extended.
The proposed model can serve as a valuable reference point for dynamical systems in which not only two but more contagion processes cooperate, potentially with a variety of disease specific parameters. For example, an entire network of diseases with different degrees of cooperation or conflict may yield a complex system of possible stationary state and much more complex transitions as parameter thresholds are crossed.
Furthermore, because the simple SIS dynamical system finds applications outside of epidemiology, for example in adoption dynamics, our results may also find applications in situations when technologies are adopted in social systems.