Fundamental properties of cooperative contagion processes

Contagion processes abound in nature, ranging from the spread of infectious diseases in host populations [1], the spread of information in social networks [2], the adaptation of technology and norms [3, 4], to activation patterns in neural tissue [5, 6], and escape mechanisms from predators in schooling fish [7]. Dynamical computational models are an essential tool for understanding phenomena in all of these contexts. Their application to the spread of infectious diseases has flourished in recent years [8–11], primarily because of the relevance to human health and the spread of human infectious diseases. Dynamical models cover a broad scope in terms of complexity, ranging from qualitative models that focus on universal features of the observed phenomenon [12, 13], network models that account for population structure or host mobility patterns [14–20], to sophisticated, large-scale agent-based models [21, 22] 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 and predicting the trend of ongoing epidemics.

Most epidemic models focus on the transmission dynamics of single, symptomatic pathogenic bacteria or viruses because in most applications it can reasonably be assumed that the phenomena are dominated by host-pathogen interactions. A variety of infectious diseases exist, however, that interact either directly or indirectly e.g. by altering the susceptibility of the host with respect to infection with another pathogen. Furthermore, transmissions of bacterial microorganism between host individuals is not restricted to species that cause disease. The transmission and spread of commensal or mutualistic bacteria as part of the host's microbiome is generic, in fact also often required to sustain a healthy, host specific microbiome. Especially the transmission of bacterial species of the human microbiome has attracted much attention in very recent studies [23, 24]. Microbiotic species are part of a complex microbiological ecosystem of a host, with a densely connected set of metabolic connections [25]. It is reasonable to assume that these interactions, and the presence of particular species in a host's microbiome impacts the propensity of colonization with another. In social science, the adoption of a certain behavioral patterns may impact the propensity to adopt another pattern if exposed to it. Therefore, it is important to understand the basic mechanisms and effects that are generated by interactions of contagion processes in general.

Early network theoretic work focused on competitive coinfection [26–33] with important applications to infection dynamics of virus strains that induce cross immunity. In these systems different pathogens suppress each other, yielding either single pathogen dominance or coexistance solutions. Multiplex network approaches have been applied in this context [34, 35], for situations when contagion processes unfold along different modes of transmission within the same population [26–28, 36]. In [36] a very general framework for describing two concurrent diseases is introduced, relevant factors such as infection type, impact of the underlying networks, positive and negative interactions are analyzed and discussed.

Only recently, cooperative contagion in which infection with one transmissible agent facilitates infection with another was investigated [36–47]. These studies mainly focused on transient dynamics of the generic susceptible-infected-recovery (SIR) model in which individuals acquire immunity after infection. In [39], a simple SIR coinfection model 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. Analytical insights [44] have been obtained that explain the role of network topology in cooperative bond percolation systems, in multiplex systems [45], power-law networks [43], as well as sequential coinfection on Poisson networks [37]. 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 [41]. Two asymmetrically interacting SIR contagion processes were investigated in [47] and backwards bifurcations, i.e. first order phase-transitions, were identified.

Because most of these models focus on transient SIR dynamics they cannot capture situations in which a steady supply of susceptibles permits the existence of a stable endemic state, such as the susceptible-infected-susceptible (SIS) or SIRs or SIR model with vital processes. An exception is [41] where coinfection of two SIS type processes was investigated using pairwise level approximation. Despite of this, co-contagion systems remain poorly understood and some fundamental issues remain: what basic dynamical features can we expect in cooperative contagion processes? To what extent does cooperativity change the classic outbreak scenario, what is the nature of transitions to endemic states? Can we expect multi-stability and multiple thresholds? How does cooperativity impact spatial propagation?

Here we introduce and investigate a model for the dynamics of two transmissible, interacting agents (labeled A and B). The model is based on the well-known SIS model in which host individuals are either susceptible (S) or infected (I). Susceptibles can be infected with either agent. When infected with say A they can transmit A to other susceptibles. Infecteds remain in the infectious state for a typical period after which they recover and susceptible again. The transmission dynamics of agents A and B are governed by agent specific baseline reproduction numbers ${R}_{{A}}$ and ${R}_{{B}}$, respectively that describe the dynamics of an agent in the absence of the other. We incorporate cooperativity by two additional parameters, the cooperativity coefficients ${\xi }_{{A}}$ and ${\xi }_{{B}}$ that capture influence of an infection with A on the subsequent infection with B and vice versa.

Based on this model, we show that cooperativity between contagion processes generates a variety of interesting properties that are absent in single agent dynamics. For sufficiently strong cooperativity, increasing the baseline reproduction number of one or both agents yields abrupt, discontinuous outbreak transitions and multi-stability (i.e. the coexistence of different stable asymptotic states). Furthermore, cooperativity exhibits dynamical hysteresis, a consequence of the split of the ordinary epidemic threshold into two separate thresholds (an inception and extinction threshold). We derive these features analytically in a deterministic well-mixed model. Their robustness is corroborated by numerical simulations of analogous stochastic dynamical processes on both lattices and generic network systems. Finally we investigate cooperative contagion in spatially extended systems. We show that the interplay of different thresholds and hysteresis yields a rich set of wavefront dynamics and invasion dynamics, e.g. accelerated propagation in certain parameter regimes, stable heterogeneous patterns as well as negative wavefront speeds (receding wavefronts).

Our model is an extension of the generic SIS compartmental model: host individuals are either susceptible (S) or infected (I) and change state by two reactions, the transmission of the infection $S+I\to 2I$ and recovery $I\to S$ at rates α and β, respectively. In a well mixed, large, and conserved population the fraction of infected individuals u(t) can be described by $\dot{u}=\alpha u(1-u)-\beta u$. The basic reproduction ratio is defined by $R=\alpha /\beta$. For $R\lt 1$ the trivial state u = 0 is globally stable. If R is increased beyond the critical threshold ${R}_{c}=1$ the system exhibits a transcritical bifurcation, u = 0 becomes unstable and $u=1-{R}^{-1}$ is the stable endemic state in which transmission and recovery events balance. The SIS system thus exhibits a continuous transition as R crosses the critical threshold ${R}_{c}=1$. Analogous stochastic lattice models in which lattice sites can transmit to neighboring sites and recover exhibit the same type of threshold behavior and a continuous phase transition. Here, we consider a generalization of the SIS model that captures the dynamics of two interacting transmissible agents: A and B. A host can be in one of four states S, A, B, and AB, corresponding to susceptible, infected with A but not B, infected with B but not A, and infected with both A and B, respectively, see figure 1. Transmissions in this system occur by interactions of host individuals in these four different states and can be summarized as follows:

$$\begin{eqnarray}&&{A}\vee {AB}+S\mathop{\longrightarrow }\limits^{{\alpha }_{{A}}}{A}\vee {AB}+{A},\\ &&\quad {B}\vee {AB}+S\mathop{\longrightarrow }\limits^{{\alpha }_{{B}}}{B}\vee {AB}+{B},\\ &&\quad {A}\vee {AB}+{B}\mathop{\longrightarrow }\limits^{{\alpha }_{{BA}}}{A}\vee {\text{AB}}+{AB},\\ &&\quad {B}\vee {AB}+{A}\mathop{\longrightarrow }\limits^{{\alpha }_{{AB}}}{B}\vee {\text{AB}}+{AB},\end{eqnarray} \tag{ 1 }$$

where e.g. $A\vee {AB}$ represents an individual that is either in state A or in state AB such that e.g. the first reaction represents the transmission of agent A to a susceptible individual. The system is defined by four different transmission rates ${\alpha }_{{A}}$, ${\alpha }_{{B}}$, ${\alpha }_{{BA}}$, and ${\alpha }_{{AB}}$ that correspond to transmission of A to a susceptible, of B to a susceptible, of A to an individual carrying B, of B to an individual carrying A, respectively. For simplicity we assume uniform recovery rates:

$$\begin{eqnarray}&&{AB}\mathop{\to }\limits^{\beta }A\vee B,\qquad A\vee B\mathop{\to }\limits^{\beta }S.\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Because we focus on cooperative contagion we restrict the transmission rates:

$$\begin{eqnarray}&&{\alpha }_{{BA}}={\xi }_{{B}}{\alpha }_{{A}}\geqslant {\alpha }_{{A}},\qquad {\alpha }_{{AB}}\equiv {\xi }_{{A}}{\alpha }_{{B}}\geqslant {\alpha }_{{B}},\end{eqnarray} \tag{ 3 }$$

cooperativity thus implies that ${\xi }_{{A}},{\xi }_{{B}}\geqslant 1$. For example, a value ${\xi }_{{A}}=5$ means that transmission of B to an individual already carrying A is 5-fold the transmission compared to the baseline transmission to an S individual. Based on the above reactions one can obtain a set of ordinary differential equations for the fraction of individuals in each state. The reactions above, however, suggest a more suitable set of compartments ${ \mathcal S }$, ${{ \mathcal A }}^{+}$, ${{ \mathcal B }}^{+}$ and ${ \mathcal A }{ \mathcal B }={{ \mathcal A }}^{+}\cap {{ \mathcal B }}^{+}$ with the corresponding dynamical variables s, u, v, and w: the fractions of susceptibles, individuals infected with A (including those that are also infected with B), individuals infected with B (including those that are infected with A), and individuals infected with both A and B, respectively, see figure 1. In the limit of a large, well-mixed host population the dynamics is described by

$$\begin{eqnarray}\dot{u} & = & {R}_{{A}}{su}+{\xi }_{{B}}{R}_{{A}}(v-w)u-u,\\ \dot{v} & = & {R}_{{B}}{sv}+{\xi }_{{A}}{R}_{{B}}(u-w)v-v,\\ \dot{w} & = & {\xi }_{{A}}{R}_{{B}}(u-w)v+{\xi }_{{B}}{R}_{{A}}(v-w)u-2w,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&s=1-u-v+w,\end{eqnarray} \tag{ 5 }$$

where ${R}_{{A}}={\alpha }_{{A}}/\beta$, ${R}_{{B}}={\alpha }_{{B}}/\beta$ and time is measured in units of ${\beta }^{-1}$. For cooperativity coefficients ${\xi }_{{A}}={\xi }_{{B}}=1$ the above system describes two independent contagion processes: if ${R}_{{A}},{R}_{{B}}\gt 1$ the stable endemic state is given by ${u}^{\star }=1-{R}_{{A}}^{-1}$, ${v}^{\star }=1-{R}_{{B}}^{-1}$, ${w}^{\star }=(1-{R}_{{A}}^{-1})(1-{R}_{{B}}^{-1})$, and ${s}^{\star }={({R}_{{A}}{R}_{{B}})}^{-1}.$

We now consider the effect of cooperativity. In the following and in analogy to the labels used to identify the state of an individual in the population, it is useful to assign the same labels S, A, B, and AB to the potential asymptotic states of the entire host population. We say, e.g., that the system is in state A if only agent A is present in the population, the contagion free state is S, etc. We begin with a symmetric system in which ${\xi }_{{A}}={\xi }_{{B}}=\xi$ and identical baseline reproduction ratios ${R}_{{A}}={R}_{{B}}=R$. In this case the above system reduces to:

$$\begin{eqnarray}\dot{u} & = & {Rsu}+\xi R(v-w)u-u,\\ \dot{v} & = & {Rsv}+\xi R(u-w)v-v,\\ \dot{w} & = & \xi R[2{uv}-(u+v)w]-2w.\end{eqnarray} \tag{ 6 }$$

Figure 2 illustrates the bifurcation analysis of the system. When $1\leqslant \xi \lt 2$ the system exhibits a behavior similar to independent contagion processes: at R = 1 we observe a transcritical bifurcation yielding a stable endemic population state AB for $R\gt 1$. This means that even when cooperativity amplifies transmission rates by up to a factor of two, we see no qualitative dynamical difference.

Zoom In Zoom Out Reset image size

However, when cooperativity exceeds a critical magnitude, i.e. for $\xi \gt {\xi }_{c}=2$, a different bifurcation behavior emerges. As R is increased and before the conventional critical point ${R}_{c}=1$ is reached, a saddle-node bifurcation emerges at ${R}_{e}=2\sqrt{\xi -1}/\xi \lt 1$. When ${R}_{e}\lt R\lt {R}_{c}$, in addition to the trivial stable state S, two AB stationary states exist:

$$\begin{eqnarray}{u}_{\pm }^{\star }={v}_{\pm }^{\star } & = & \displaystyle \frac{\xi R-2\pm \sqrt{{\xi }^{2}{R}^{2}-4\xi +4}}{2\xi R},\\ {w}_{\pm }^{\star } & = & \displaystyle \frac{{{Ru}}_{\pm }^{\star }(\xi -2)+R-1}{(\xi -1)R},\end{eqnarray} \tag{ 7 }$$

one of which $({u}_{+}^{\star },{v}_{+}^{\star },{w}_{+}^{\star })$ is stable, the other unstable. Thus, when ${R}_{e}\lt R\lt {R}_{c}$ sufficiently small perturbations to the S state will have no effect as S is stable. However, when perturbations are sufficiently large, the system will approach the endemic AB state with ${u}_{+}^{\star }={v}_{+}^{\star },{w}_{+}^{\star }$ . Furthermore, when R is increased beyond the critical value ${R}_{c}=1$, state S loses stability and any arbitrarily small perturbation will yield a discontinuous jump to the endemic state, reminiscent of a first order phase-transition. For example when $R={R}_{c}+\varepsilon$ with $0\lt \varepsilon \ll 1$ the stable endemic state is ${u}_{+}^{\star }={v}_{+}^{\star }\approx (\xi -2)/\xi$, ${w}_{+}^{\star }\approx {(\xi -2)}^{2}/\xi (\xi -1)$. So if say $\xi =10$ this yields an endemic state in which 71% of the population is in state AB, immediately after R_c is crossed. Cooperative contagion also exhibits hysteresis: starting with $R\gt {R}_{c}$ and state AB, decreasing R across the critical value R_c from above will not result in immediate extinction. A high endemic state is maintained until the eradication threshold R_e is reached, which can be substantially smaller than the ordinary epidemic threshold R_c. Decreasing R below R_e will then yield a sudden collapse of AB into S. These newly observed dynamics is also termed as backward bifurcation, and its ocurrence mathematically could come from quite different aspects [47–50]. It is also worth to note that complex contagions through reinforced infection via different neighboring nodes but of single agent can also generate a similar phase transitions [51, 52].

Equation (6) capture the symmetric special case of the more general system defined by equation (4), the latter of which has four parameters, ${R}_{{A}}$, ${R}_{{B}}$, ${\xi }_{{A}}$, and ${\xi }_{{B}}$. Figure 3 illustrates the phase diagram for a more general choice of these parameters. Fixing the cooperativity coefficients to ${\xi }_{{A}}={\xi }_{{B}}=5$ we investigated the phases of the system as a function of baseline reproduction ratios ${R}_{{A}}$ and ${R}_{{B}}$. Apart from the expected stable states we observe a rich variety of bistable states in the region in which baseline reproduction is near unity, as is illustrated in figure 4. For example, when ${R}_{{B}}=1.1$ and starting with ${R}_{{A}}\approx 0$ the system is initially in state B. Increasing ${R}_{{A}}$ further a saddle-node bifurcation occurs and the system enters a regime in which B and AB are both stable. When ${R}_{{B}}\lt 1$ , e.g. ${R}_{{B}}=0.57$, increasing ${R}_{{A}}$ first yields an ordinary transcritical bifurcation to the A state, followed by a second bifurcation into a regime in which A and AB are stable, and finally, a third bifurcation to into the regime in which only AB is stable, see also figure 3. A key property of the system is that the complexity of transitions is only observed for baseline reproduction ratios near unity. If one baseline reproduction is too low or two high, only a single ordinary transitions and no state coexistence is observed. This is an interesting property from an evolutionary point of view. When new strains of transmissible agents emerge, typically they are not adapted to the host and possess baseline reproduction not significantly larger than unity, or even smaller. Cooperativity with other transmissible agents and coexistence of stable endemic states may present an opportunity for developing a species rich system with higher evolvability. This type of complexity is expected to increase dramatically when more than two transmissible agents are involved, yielding a potentially rich space of stable states and an increasing complexity in phase separation manifolds in parameter space.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The deterministic model discussed above cannot account for fluctuations or population heterogeneities. An important question is therefore whether the observed phenomena prevail in a more complex scenario in which transmissions and recovery events are stochastic, the host population is finite and not every host interacts with every other host at equal rates. Typically, stochastic effects in a well-mixed system are modeled by birth-death type stochastic processes equivalent to the reactions depicted in figure 1 and generating solutions to the corresponding master-equation for a fixed but finite population size N. Population heterogeneities are typically addressed by modeling these processes on fixed network topologies or lattices in which host individuals only interact with the neighbors defined by the network links. In order to address the robustness of effects and properties derived for the deterministic system of equation (4) we investigated cooperative contagion dynamics in a stochastic 2d-lattice system and and Erdős-Rényi (ER) network with equal mean degree and number of nodes. The results are compiled in figure 5. In both cases, we observe hysteresis, and a separation of extinction and outbreak thresholds for large cooperativity coefficients ξ. Interestingly, the extinction transition is continuous in the lattice, a consequence of the local coupling of the system. The ER network exhibits discontinuous transitions, as predicted by the above mean field treatment. This is not surprising as the ER network is topologically more similar to the well-mixed scenario, and lattice of low dimension like 2D is physically far from that. Based on these observations, we believe that the key features of cooperative contagion can be expected also in more realistic, structured populations⁶ .

Zoom In Zoom Out Reset image size

An important aspect of contagion processes is their spatial propagation. When simple contagion processes with $R\gt 1$ are seeded in a spatially homogeneous susceptible host population and contagion dynamics is combined with diffusive dispersal of the host these systems typically exhibit propagating wavefronts that travel at constant speeds. The endemic state invades the unstable S domain. The generic SIS contagion process, e.g., can be described by:

$$\begin{eqnarray}&&{\partial }_{t}u=R(1-u)u-u+D{\partial }_{x}^{2}u,\end{eqnarray} \tag{ 8 }$$

where $u=u({\bf{x}},t)$ is the density of infected individuals at location ${\bf{x}}$ at time t. The combination of local initial exponential growth (for $R\gt 1$) and diffusion yields a front-speed depending on the basic reproduction ratio and diffusion coefficient D: $c=2\sqrt{(R-1)D}$. This is a generic feature of processes that exhibit pulled fronts [53]. Given the more complex nature of cooperative coinfection, especially the dynamical bistability for intermediate baseline reproduction ratio ${R}_{e}\lt R\lt {R}_{c}$ and large cooperativity coefficient ξ, we can expect a richer set of phenomena when cooperative contagion processes expand in space. To account for a diffusing host we extend equation (4) and consider the corresponding reaction–diffusion system:

$$\begin{eqnarray}{\partial }_{t}u & = & {f}_{u}(u,v,w)+D{\partial }_{x}^{2}u,\\ {\partial }_{t}v & = & {f}_{v}(u,v,w)+D{\partial }_{x}^{2}v,\\ {\partial }_{t}w & = & {f}_{w}(u,v,w)+D{\partial }_{x}^{2}w,\end{eqnarray} \tag{ 9 }$$

where D in last terms is the diffusion coefficient and the functions f_u, f_v, and f_w are the same as on the rhs of equation (6). The dynamical variables are function of time t and space ${\bf{x}}$, e.g. $u=u({\bf{x}},t)$. We assume that the diffusion coefficient is constant and independent of the state of a host individual primarily focusing on contagion processes that do not affect the host's dispersal behavior. We also consider a constant overall density which implies that $s=1-u-v+w$ at every position ${\bf{x}}$. As before we use labels S, A, B and AB to refer to region that are contagion free, only affected by A, only affected by B, and both A and B, respectively.

The system defined by equation (9) exhibits a range of front velocities, each one corresponding to different states invading regions in a different state. For example a localized A-patch invades an S-region at a different speed than a uniform B-region (turning the latter into an ${AB} \mbox{-}$region). A localized AB-patch invades an S-region differently than an A-region. To understand the asymptotics and transients of the system we first consider a uniform population in state S, with the exception of two localized patches, each being in state A and B respectively and separated by some distance, see figure 6. When $R\gt {R}_{c}$ cooperative contagion plays no role at the beginning, each patch will expand at a constant front speed of ${c}_{0}\propto \sqrt{(R-1)D}$. Once these growing patches touch, cooperativity kicks in at the A–B interface. The emerging AB-nucleus has interfaces to the A and B regions as well as to the S region. For $\xi \gt 1$ the invasion of AB into the S-region is faster at a speed ${c}_{{AB}\to S}\gt {c}_{0}$ as expected. Interestingly, the invasion of the ${AB} \mbox{-}$region into A-region (and B-region) occurs at an even higher speeds ${c}_{{AB}\to A,B}\gt {c}_{{AB}\to S}$. Using a propagating wave ansatz $u=u(x-{ct})$ for a 1-d spatial support (analogously for variables v and w) one can compute a lower bound ${c}_{{AB}\to A,B}\approx \sqrt{\xi (R-1)D}=\sqrt{\xi }{c}_{0}$. Because ${c}_{{AB}\to A,B}\gt {c}_{{AB}\to S}$ the system will eventually converge to a uniform AB-region that spreads at speed ${c}_{{AB}\to S}$. Regions affected only by one agent will not persist. This effect of enhanced wave-front speed might be particularly relevant in situations in which a covert, unknown and commensal agent is endemic in some region and a known process with known baseline reproduction ratio expands somewhere else in the system at a speed that is computed based on its baseline reproduction ratio. If this front enters the region in which the unknown but highly cooperative covert agent prevails, a sudden but potentially unexpected boost in the proliferation of the initial spreading process could occur.

Zoom In Zoom Out Reset image size

In the bistable region ${R}_{e}\lt R\lt {R}_{c}$, isolated islands of A nor B cannot persist. If we initialize the system with A and B patches that share a small overlapping region in the AB state cooperativity can yield the survival of the AB state while the homogeneous A and B states fade. The remaining AB patch then proliferates at speed ${c}_{{AB}\to S}$. Interestingly, we observe negative propagation speed in this regime, ${c}_{{AB}\to S}\lt 0$, which implies a receding AB-region. This behavior is caused by the dispersal of A or B affected individuals into the S-region in which the S-state is also stable. Once individuals enter this state, they have a higher likelihood of becoming susceptible than being colonized by both agents. The wavefront acts as a drain for infected agents and a competition exists between the supply of new agents of type A or B and the diffusive dilution of their concentration. For a critical choice of parameters, e.g. the baseline reproduction ratio we observe a stationary heterogeneous solution, an immobile front, that separates S from AB regions. Figure 7 illustrating the three typical propagation modes in 2d space are depicted (see also the movies in the supplemental materials⁷ , available online at: stacks.iop.org/NJP/19/103041/mmedia).

Zoom In Zoom Out Reset image size

The emergence of new front wave propagation modes is dynamically rooted in the multi-stability of coinfection dynamics. The next step would be an analytical derivation e.g. of the dependence of R_st on system parameters and physical quantities. This remains an open question, and some previous studies offer a good starting point [54–57]. For example, [55] studied a cooperative reaction–diffusion system analytically, in which one species diffuses while the other is immobile.

We present a reaction kinetic model for the dynamics of cooperative contagion of the SIS class. The most prominent property of the model is the existence of discontinuous transitions to an endemic state when the traditional outbreak threshold is crossed and a separation of outbreak and extinction thresholds, the magnitude of which depends on the degree of cooperativity. Although we derive the key properties analytically and numerically in a deterministic model suitable for large, well-mixed populations, we observe the key features of discontinuous transitions also in a stochastic network variant of the model. The system of two cooperating agents that proliferate in a host population exhibits diverse properties when spatial diffusion is incorporated yielding different types of transients and spreading speeds.

The proposed model and presented results can be employed and adapted to understand realistic systems, e.g. pneumonia in which with bacterium like Streptococcus pneumoniae interacts with viral respiratory infections (e.g. influenza viruses), where one pathogen increases the susceptibility towards the other up to 100-fold [58, 59]. Another prominent example are human immunodeficiency virus (HIV) syndemics [60], where the suppressed immune system of the hosts greatly increases the susceptibility towards secondary infections: hepatitis, malaria, syphilis, herpes virus and tuberculosis. In the latter case, the cooperative interactions are mutual, as hosts with tuberculosis are found also more likely to acquire HIV [61, 62]. These new complexities uncovered in our study suggest that we may need to devise new containment strategies, e.g. by vaccination programs [63], to combat the epidemic spreading in some more realistic circumstances.

Although we discussed the model predominantly in the context of the spread of transmissible diseases, the model is sufficiently generic to be applied to other transmissible contagion processes that influence each other cooperatively, e.g. the adoption of one technology may increase the infection of a user with another type of technology which may then occur explosively or at different speeds than expected. While in the context of epidemic spreading, the implication of our results is bad news for the infection and is detrimental, but in the majority of other contexts such as social contagions, where a high prevalence is usually desired, our results are good news indeed.

Considering a model for only two interacting agents is a foundation that can easily be generalized to a larger number, potentially a network, of interacting agents. If the baseline reproduction of a family of transmissible agents is in the critical regime, we expect in such a system a diverse set of stable configurations and we believe that the model presented here is a helpful starting point for investigating these more general systems.

LC and FG would like to thank W Cai and P Grassberger for valuable comments and discussions. DB acknowledges fruitful discussion with Osamah Hamouda and Imbish Mortimer and Annegret Tierkrper. LC would also like to thank the hospitality of Beijing Computational Science Research Center (CSRC) for a pleasant stay during the final revision stage.

LC, FGh, DB conceived research, LC and DB developed dynamical model, LC, DB theoretical, analytical treatment of mean field and spatial model, LC performed numerical analysis. DB and LC wrote manuscript.