Diversity of hysteresis in a fully cooperative coinfection model

TL;DR

This paper introduces a fully cooperative coinfection model demonstrating three novel types of hysteresis in disease spread, with detailed analysis and numerical validation on network structures.

Contribution

It identifies and characterizes three new types of hysteresis in a fully cooperative coinfection model, expanding understanding of epidemic transitions.

Findings

Three types of hysteresis ($C$, $S_l$, $S_r$) identified in the model.

Discontinuous outbreak and eradication transitions characterized.

Hysteresis types analyzed in well-mixed and network-based approaches.

Abstract

We propose a fully cooperative coinfection model in which singly infected individuals are more likely to acquire a second disease than those who are susceptible, and doubly infected individuals are also assumed to be more contagious than those infected with one disease. The dynamics of such fully cooperative coinfection model between two interacting infectious diseases is investigated through well-mixed and network-based approaches. We show that the former approach exhibits three types of hysteresis, namely, $C$, $S_l$ and $S_r$ types, where the last two types have not been identified before. The first (resp., the second and the third) type exhibits (resp., exhibit) discontinuous outbreak transition from the disease free (resp., low prevalence) state to the high prevalence state when a transmission rate crosses a threshold from the below. Moreover, the third (resp., the first and the second) type possesses (resp., possess) discontinuous eradication transition from the high prevalence state to the low prevalence (resp., disease free) state when the transmission rate reaches a threshold from the above. Complete characterization of these three types of hysteresis in term of parameters measuring the uniformity of the model is also provided. Finally, we assess numerically this epidemic dynamics in random networks.

Tables

Table 1. Table 1: Basic reaction scheme for the disease spreading in the well-mixed model and stochastic network models.

Transmission	Well-mixed           (rate)	Stochastic network      (probability)
$S\stackrel{𝐴}{\to }A$	$\alpha$	$p$
$S\stackrel{𝐵}{\to }B$	$\alpha$	$p$
$S\stackrel{AB}{\to }A$	$a\alpha$	$1-{\left(1-p\right)}^a$
$S\stackrel{AB}{\to }B$	$a\alpha$	$1-{\left(1-p\right)}^a$
$A\stackrel{𝐵}{\to }AB$	$\beta \left(=C\alpha \right)$	$1-{\left(1-p\right)}^C$
$B\stackrel{𝐴}{\to }AB$	$\beta \left(=C\alpha \right)$	$1-{\left(1-p\right)}^C$
$A\stackrel{AB}{\to }AB$	$b\beta \left(=bC\alpha \right)$	$1-{\left(1-p\right)}^{bC}$
$B\stackrel{AB}{\to }AB$	$b\beta \left(=bC\alpha \right)$	$1-{\left(1-p\right)}^{bC}$

Equations

$$\dot{S}(t)$$

$${\dot{I}}_A(t)$$

$${\dot{I}}_B(t)$$

$${\dot{I}}_{AB}(t)$$

$$I_A^{\ast }=I_B^{\ast }$$

$$I_{AB}^{\ast }$$

$$f(S,\alpha ):=$$

$$+[\alpha a(a-b)C-abC+2a-1]{\alpha }^2{S}^2$$

$$+(\alpha abC-2a+2)\alpha S-1=0.$$

$$T=\alpha S.$$

$$\bar{\alpha }(S)$$

$$=:\tilde{\alpha }(T).$$

$$C_2=\frac{2}{a}.$$

$$\text{mbox}{C}_1$$

$$C_{12}$$

$$\alpha_{e}(C)\left\{\begin{array}[]{cl}>1&\text{if }C<C_{12},\\[1.0pt] =1&\text{if }C=C_{12},\\[1.0pt] <1&\text{if }C>C_{12}.\\[1.0pt] \end{array}\right.$$

$$\text{sgn}\left({\left[\frac{d}{dT}\tilde{\alpha }(T)\right]}_{T=t}\right)=\text{sgn}\left({\left[\frac{d}{dS}\bar{\alpha }(S)\right]}_{S=\frac{t}{\tilde{\alpha }(t)}}\right).$$

$${\left[\frac{d}{dS}\bar{\alpha }(S)\right]}_{S=\frac{t}{\tilde{\alpha }(t)}}={\left[\frac{\frac{d}{dT}\tilde{\alpha }(T)}{\frac{d}{dT}\left[\frac{T}{\tilde{\alpha }(T)}\right]}\right]}_{T=t}.$$

$$\frac{d}{dT}\left[\frac{T}{\tilde{\alpha }(T)}\right]=\frac{N_1(T,a,b)}{aCT[\tilde{\alpha }(T){]}^2},$$

$$N_1(T,a,b)$$

$$+\{T(1-T)[2b(1-T)+(3-T)]\}a$$

$$+b(2-T)(T-1{)}^2$$

$$>0.$$

$${\tilde{\alpha }}^{\prime }(T)$$

$${\tilde{\alpha }}^{\prime \prime }(T)$$

$${\tilde{\alpha }}^{\prime \prime \prime }(T)$$

$$M_1(T,a,b,C)$$

$$+[a{b}^{2}C-2a^2+2a-b]T^2$$

$$+2(b-a)T-b,$$

$$M_2(T,a,b)$$

$$+3(a-b{)}^2{T}^2+3b(a-b)T+b^2,$$

$$M_3(T,a,b)$$

$$-4(a-b{)}^3{T}^3-6b(a-b{)}^2{T}^2$$

$$-4b^2(a-b)T-b^3.$$

$$M_3(T,a,b)$$

$$=b^2\left[c_0{b}^{-2}+c_1{b}^{-1}+c_2\right]+c_3{b}^3,$$

$$\tilde{\alpha}^{\prime\prime}(T)\left\{\begin{array}[]{cl}>0&\text{if }T\in(0,\gamma),\\ =0&\text{if }T=\gamma,\\ <0&\text{if }T\in(\gamma,1].\\ \end{array}\right.$$

$${\tilde{\alpha }}^{\prime }(\gamma )={\tilde{\alpha }}^{\prime \prime }(\gamma )=0.$$

$$\tilde{\alpha }(T)=1(=\tilde{\alpha }(1))$$

$$\text{sgn}({\tilde{\alpha }}^{\prime }(T))=\text{sgn}(M_1(T,a,b,C))$$

$$M_1(T,a,b,C)$$

Full text

Coinfection, bistability, hysteresis, multiple epidemic equilibria

I INTRODUCTION

Mathematical modelling has become a fundamental tool to understand the basic features of disease spreading in a group Hethcote2000 . We can deal with this problem with two approaches: considering a well-mixed system and approximating it with an interaction all-to-all, i.e., the mean-field approximation, or taking into account a topology that specifies who interacts with whom. Hence, the theory of both dynamical systems and complex networks Newman2003 is fundamental for this kind of studies. Two main models have been proposed for the dynamics of a spreading disease: one, considering three kinds of individuals, where susceptibles (S), when exposed to the disease, become infective (I), and afterwards become recovered (R), getting immunized against the disease, is known as SIR model KermacK1927 ; the other model represents diseases that are spread by pathogens that frequently mutate, such that individuals do not get immunized, but come back to susceptible state after the infective phase (SIS model) Anderson1992Book . In fact, those models have been proposed not only in the context of disease spreading, but also for other spreading processes, such as computer viruses epidemics Pastor-Satorras2001A .

Most studies on disease spreading focus on one disease infecting a population Pastor-Satorras2001B ; Moreno2002 ; Newman2002 . However, different pathogens may interact with each other on the host. For instance, there are cases in which infection with one disease strengthens a host’s immunity against other infection Newman2005 ; Ahn2006 ; Funk2010C ; Marceau2011 ; Karrer2011 ; Poletto2013 ; Sahneh2014 . In those competitive systems, regimes such as extinction, coexistence and one strain dominance have been identified. On the other hand, it has also been reported that people with HIV infections are more likely to get other infections such as hepatitis or tuberculosis Alter2006 ; Abu-Raddad2006 ; Pawlowski2012 ; WHO2014 . Such cooperative disease spreading models have recently been studied in Refs. Hethcote2000 ; Pastor-Satorras2001A ; Newman2005 ; Ahn2006 ; Alter2006 ; Abu-Raddad2006 ; Colizza2006 ; Barrat2008Book ; Funk2010C ; Marceau2011 ; Karrer2011 ; Pawlowski2012 ; Poletto2013 ; Chen2013 ; WHO2014 ; Sahneh2014 ; Sanz2014 ; Cai2015 ; Grassberger2016 ; Chen2016 . Interesting dynamics, like the appearance of abrupt transitions, in contrast to the typical second order phase transition found in single disease spreading model, have been reported Cai2015 ; Grassberger2016 ; Chen2016 . Specifically, for SIS dynamics, there is hysteresis behaviour: the models exhibit discontinuous outbreak transition from disease free state to the high prevalence state as a transmission rate crosses a threshold $\alpha_{o}$ from the below, while in the reverse scenario, they exhibit discontinuous eradication transition from the high prevalence state to disease free state as the transmission rate reaches a threshold $\alpha_{e}\,(<\alpha_{o})$ from the above is also observed. The above described phase transition is to be called a hysteresis of $C$ type. At $\alpha=\alpha_{o}$, we see that the size of the infected population has a discontinuous jump, indicating the outbreak of the disease. The difference in $\alpha_{e}$ and $\alpha_{o}$ points out that, once the outbreak of the disease occurs, driving down the transmission rate to $\alpha_{o}$ is not enough to eradicate the diseases. In fact, we need to drive the rate further down to $\alpha_{e}$ for the epidemic to die out, which requires more effort and results in greater economic costs. Note that, in the cooperative coinfection models considered in Refs. Chen2013 ; Cai2015 ; Grassberger2016 ; Chen2016 , it was only assumed that a singly infected individual is more likely to get infected with a second disease than a susceptible individual.

In this work, we consider a fully cooperative coinfection model for two diseases, in which we further assume that a doubly infected individual is more contagious than a singly infected individual; see Fig. 1a. Hence, diseases cooperate not only when one individual is getting infected, but also when they meet in one host and start infecting other. We observe two additional types of hysteresis. In particular, the model has discontinuous outbreak transition from the low prevalence state to the high prevalence state as the transmission rate crosses a threshold $\alpha_{o}$ from the below. Moreover, a discontinuous eradication transition from the high prevalence state to the disease free state or the lower prevalence state can both be observed as the transmission rate reaches a threshold $\alpha_{e}\,(<\alpha_{o})$ from the above. They are to be called hysteresis of $S_{l}$ and $S_{r}$ types, respectively; see Fig. 1b. It is worthwhile to point out that, for model exhibiting hysteresis of $S_{l}$ or $S_{r}$ types, it must exist two stable endemic fixed points in some range of transmission rates as contrast to the existence of at most one stable endemic fixed point for all range of transmission rates in the case of hysteresis of $C$ type.

To depict the above mentioned dynamics, we first consider the model through the well-mixed approach. Three positive quantities $C$, $a$, $b$, called the relative infection factor, relative coinfection factors for susceptible class and singly infected class, respectively, are to be introduced to measure the uniformity of the model; see Fig. 1a. When $C=a=b=1$, the model is said to be uniform and corresponds to noninteracting disease spreading model. The larger deviation of these quantities from $1$, the more nonuniform the model is. The case that $\min\{C,a,b\}>1$ (resp., $\max\{C,a,b\}<1$) corresponds to a fully cooperative (resp., competitive) coinfection model. If $\min\{C,a,b\}<1<\max\{C,a,b\}$, it is called a mixed type. In Ref. Chen2016 , the authors consider a partially cooperative model, i.e., $C>1$ and $a,b=1$, where the hysteresis of type $C$ was found provided that $C$ is sufficiently large. In this work, for sufficiently large $C$ and any $a,b\geq 1$, the model behaves like the one considered in Ref. Chen2016 . In case that $b-a>1/2$ and $C$ in an intermediate range, hysteresis of types $S_{l}$ and $S_{r}$ emerge. Subsequently, we investigate this problem through the network-based approach.

The organization of the paper is as follows. In Section II, the well-mixed approach for our fully cooperative coinfection model is introduced and analytically studied. In Section III, we describe our stochastic network model and the numerical simulation results are also provided. The summary of our obtained results and the future work are given in Section IV. All proofs of our results are recorded in Appendix A.

II Well-Mixed approach

We study a well-mixed particle system for the fully cooperative coinfection dynamics with two diseases, A and B. We define $S$, $I_{A}$, $I_{B}$ and $I_{AB}$ as, respectively, the fractions of susceptible individuals and the ones infected with disease A, B and AB in the population. Susceptible individuals get a primary infection from a singly infected with a rate $\alpha$. The cooperation between diseases on individuals that are singly infected is translated into a higher rate for getting secondary infection, $\beta>\alpha$. Moreover, the cooperation makes the doubly infected individuals more virulent, such that they infect susceptibles with a rate $a\alpha$, $a\geq 1$ and primary infected with a rate $b\beta$, $b\geq 1$. The case $a=b=1$ has been previously reported Chen2016 , implying cooperation for getting a secondary infection. Our case is more general, translating cooperation between the diseases not only for receiving, but also for transmitting diseases. Then the corresponding mean field equations for this dynamical system are:

[TABLE]

Here the recovery rates for diseases A and B, without loss of generality, are set equal to the time scale, while the basic reaction scheme underlying the dynamical process is given in the first and second columns of Table 1; see also Fig. 1a. More precisely, $\alpha$ and $\alpha a$ (resp., $\beta$ and $\beta b$) are infection rates for susceptible class (resp., singly infected class) who get infected by singly infected and doubly infected class (resp., the complementarily infected class and doubly infected class), respectively. The quantity $a$ (resp., $b$) is the ratio between the infection rates of the doubly infected class and singly infected class for infecting the susceptible class (resp., singly infected class). The larger $a$ and $b$ are, the more contagious the doubly infected class AB is. We call $a$ and $b$ the relative coinfection factors for susceptible class and singly infected class, respectively. The quantity $C:=\beta/\alpha$, which is termed the relative infection factor, is the ratio of the infection rates for singly infected class and susceptible class. The larger $C$ is, the easier singly infected class gets infected with a secondary infection. Note that $C$, $a$ and $b$ measure the degree of uniformity in fully cooperative coinfection models. The larger those quantities deviate from $1$, the more nonumiform the model is.

The dynamics of this model is depicted by its fixed points. Assuming a symmetric behaviour between the diseases A and B, the endemic fixed point $E^{*}:=(S^{*},I_{A}^{*},I_{B}^{*},I_{AB}^{*})$ of (1) satisfies

[TABLE]

where $S^{*}$, satisfying $0\leq S^{*}<1$ and $0\leq\alpha S^{*}\leq 1$, is a solution of

[TABLE]

For any $a,b\geq 1$ and $C>0$, Eq. (3) defines a solution curve $(S,\bar{\alpha}(S))$, $0<S\leq 1$, satisfying $f(S,\bar{\alpha}(S))=0$. In particular, for any fixed $a,b\geq 1$, $C>0$ and $0<S\leq 1$, there exists a unique $\bar{\alpha}\,(=\bar{\alpha}(S))$ on $(0,1/S]$ such that $f(S,\bar{\alpha}(S))=0$ (see Proposition A.1 in Appendix A). Our goal next is to study the graph the function $\bar{\alpha}(S)$, $0<S\leq 1$. We address it by making the change of variable

[TABLE]

Clearly, by (3),

[TABLE]

Here $0<T\leq 1$. Using Proposition A.2 in Appendix A, we have that the monotonicity structure of functions $\bar{\alpha}(S)$ and $\tilde{\alpha}(T)$ on the interval $(0,1]$ is identical. Hence, to study the dependence of the endemic fixed points on the parameters, it suffices to graph the function of $\tilde{\alpha}(T)$. The detailed analysis of the graph of $\tilde{\alpha}(T)$ is provided in Propositions A.3, A.4, A.5 and A.6 in Appendix A. We summarize the results in the following theorem.

Theorem 2.1**.**

For parameters $a$, $b$, $C\geq 1$, $\lim\limits_{T\rightarrow 0^{+}}\tilde{\alpha}(T)=\infty$ and $\tilde{\alpha}(1)=1$. Define $C_{2}\,(=C_{2}(a,b))$ by

[TABLE]

Then the following assertions (I) and (II) hold.

1. (I)

(See Figs. 2a and 3a) For $b\leq b_{c}\,(:=a+1/2)$, $\tilde{\alpha}(T)$ is concave upward on $(0,1]$. Moreover, the following two assertions hold.

1. (I-i)

For $C\leq C_{2}$, $\tilde{\alpha}(T)$ is strictly decreasing on $(0,1]$. Consequently, the endemic fixed point curve, denoted by $\tilde{S}$, on the $(\alpha,1-S)$ plane is strictly increasing. 2. (I-ii)

For $C>C_{2}$, $\tilde{\alpha}(T)$ has exactly one critical point at some point $T_{*}$ on $(0,1)$, which corresponds to a local minimum $\alpha_{e}\,(=\alpha_{e}(C))$. Moreover, $\alpha_{e}<\tilde{\alpha}(1)\,(=1)$. Consequently, the corresponding $\tilde{S}$, on the $(\alpha,1-S)$ plane is $C$-shaped. 2. (II)

(See Figs. 2b and 3b) For $b>b_{c}$, $\tilde{\alpha}(T)$ has exactly one inflection point $\gamma$ on $(0,1)$. Define $C_{1}\,(=C_{1}(a,b))$ and $C_{12}\,(=C_{12}(a,b))$, respectively, by

[TABLE]

Then the following three assertions hold.

1. (II-i)

For $C\leq C_{1}$, $\tilde{\alpha}(T)$ is strictly decreasing on $(0,1]$. Consequently, the corresponding $\tilde{S}$, on the $(\alpha,1-S)$ plane is strictly increasing. 2. (II-ii)

For $C_{1}<C<C_{2}$, $\tilde{\alpha}(T)$ has a unique local minimum $\alpha_{e}\,(=\alpha_{e}(C))$ and maximum $\alpha_{o}\,(=\alpha_{o}(C))$ at $T_{*}$ and $T^{*}$, respectively, on $(0,1)$ with $T_{*}<T^{*}$ such that $\tilde{\alpha}(T)$ is strictly decreasing on $(0,T_{*}]\cup[T^{*},1]$ while strictly increasing on $(T_{*},T^{*})$. Moreover,

[TABLE]

Consequently, the corresponding $\tilde{S}$, on the $(\alpha,1-S)$ plane is $S$-shaped. 3. (II-iii)

For $C>C_{2}$, $\tilde{\alpha}(T)$ has exactly one critical point at some point $T_{*}$ on $(0,1)$, which corresponds to a local minimum $\alpha_{e}\,(=\alpha_{e}(C))$. Moreover, $\alpha_{e}<\tilde{\alpha}(1)\,(=1)$. Consequently, the corresponding $\tilde{S}$, on the $(\alpha,1-S)$ plane is $C$-shaped.

The graphs of $\tilde{\alpha}(T)$ for various values of $C$ with $b\leq b_{c}$ and $b>b_{c}$ are given in Figs. 2a and 2b, respectively. From Fig. 3, we see clearly that $\tilde{S}$ is $C$-shaped (resp., $S$-shaped) as long as $C\geq C_{2}$ (resp., $b>b_{c}$ and $C_{1}<C<C_{2}$). Otherwise, it is strictly increasing. For $C_{1}<C<C_{12}$ (resp., $C_{12}<C<C_{2}$), the corresponding $\tilde{S}$ has the property that the discontinuous eradication transition is from the high prevalence state to the low prevalence state (resp., the disease free state). Hence, to further distinguish these two types of $S$-shaped curves, we shall term the corresponding $\tilde{S}$ for $C\in(C_{1},C_{12})$ and $C\in(C_{12},C_{2})$ $S_{r}$ and $S_{l}$-shaped, respectively. From our summary above, we see that nonuniformity in the model is the primary reason for the occurrence of the hysteresis. To see this, in Fig. 4, we classify the shape of $\tilde{S}$ by the uniformity measurements $a,b,C$ on the $(a,C)$ and $(b,C)$ planes with the other parameter being fixed. In Figs. 5 and 6, we consider a coinfection model of the mixed type. Specifically, we consider the parameters $a$, $b$ and $C$ satisfying $0<C<1$ and $a,b>1$. From Figs. 3, 4, 5 and 6, we arrive at the following conclusions:

1. (A)

As mentioned earlier, $C$, $a$, $b$ are the measurement for the degree of uniformity in the coinfection model. The assumptions that $b\leq b_{c}$ and $C\leq C_{2}$, or $b>b_{c}$ and $C\leq C_{1}$ indicate that the model under consideration is relatively uniform. Model exhibits a continuous transition at $\alpha=1$, the same threshold as the single disease spreading according to SIS model. In short, for such system, it acts as non interacting. Mathematically, the system undergoes a transcritical bifurcation at $\alpha=1$. 2. (B)

If the relative infection factor $C$ is sufficiently large, i.e., $C>C_{2}$, then the model exhibits the hysteresis of type $C$ regardless of what relative coinfection factors $a$ and $b$ are. In addition to a transcritical bifurcation at $\alpha=1$, a saddle-node bifurcation occurs at $\alpha=\alpha_{e}$ that generates a bistable region. A hysteresis loop is found on the interval $[\alpha_{e},1]$. The case $C=C_{2}$ connecting continuous and discontinuous transition from the disease free state has a pitchfork bifurcation at $\alpha=1$, which is supercritical as $b\leq b_{c}$ while subcritical as $b>b_{c}$. 3. (C)

If the difference in $b$ and $a$ becomes greater than $1/2$, i.e, $b>b_{c}$, then hysteresis of types $S_{r}$ and $S_{l}$ occur for $C\in(C_{1},C_{2})$. Two saddle-node bifurcations occur at $\alpha=\alpha_{e}$ and $\alpha=\alpha_{o}$. The hysteresis loop for both types is then found on the interval $[\alpha_{e},\alpha_{o}]$. As mentioned in the introductory section, the occurrence of any type of hysteresis presents greater challenge for health departments. The implication of our new results is that if the relative coinfection factors differ by less than $1/2$, hysteresis of any type will never occur in an intermediate range of $C$, i.e., $C_{1}<C<C_{2}$. Otherwise, hysteresis of types $S_{l}$ and $S_{r}$ will take place in such range. It is worthwhile to point out that the existence of hysteresis of $S_{l}$ or $S_{r}$ types guarantees the existence of two stable endemic fixed points in some range of transmission rates while only one stable endemic fixed point is needed for hysteresis of $C$. 4. (D)

If the devastation caused by the diseases is measured by how large the size of infected population jumps in the outbreak and by how much effort needed to eradicate the diseases, then hysteresis of type $C$ (resp., $S_{r}$) is the most (resp., least) devastating. 5. (E)

For models of a mixed type, we see, via Figs. 5 and 6, that three types of hysteresis still exist provided that the model is extremely nonuniform. More precisely, with $(b,C)=(25.5,0.082)$ by increasing the nonuniformity measurement in $a$, we see that the model goes through the following phase transition: continuous phase transition to hysteresis of $S_{r}$ type to continuous phase transition to hysteresis of $S_{r}$ type to hysteresis of $S_{l}$ type to hysteresis of $C$ type. 6. (F)

The most interesting case, where $b-a>1/2$, exhibits four different regimes in the parameter space, see Theorem 2.1, Figs. 4a and 4c. Those behaviours can be understood taking into account the microscopic features of the dynamics. In this situation, the cooperation mainly comes through state AB, which are considerably more virulent (Table 1). First, in the case that the value of $C$ is small, there are not enough AB states at the transition point in which the disease breaks out, so the transition is continuous. Moreover, enough AB states appear only when the system is already in a state with most of nodes infected. This makes the effect on cooperation small and hence we have a continuous branch. If the value of $C$ is further increased, the appearance of more AB states is stimulated. Eventually, in the case that enough AB states appear after a threshold $\alpha_{o}$, which is greater than the point for the outbreak of diseases, is reached, cooperation for doubly infected AB infecting other states starts to be effective, and then an abrupt transition emerges. Finally, if we keep on increasing $C$, the cooperation for secondary infection dominates the dynamics and the model matches the results reported in Ref. Chen2016 . 7. (G)

In the case that the relative infection factor $a$ for susceptible class is relatively higher than that of $b$ for singly infected class, or equivalently, $b\leq a+1/2$, there exist only two different regimes in the parameter space. The reason for the disappearance of $S_{r}$ and $S_{l}$ shape is due to the fact that a higher $a$ coupled with a larger $C$ would simulate enough states AB at the transition point where the disease breaks out. Consequently, an abrupt transmission emerges.

III Stochastic network approach

We consider a random network where diseases spread through the links that specify who interacts with whom. If the transmission probability across a link $i$ is $\bar{p}_{i}$, then the probability that a node is infected through a link with one of its $m$ infected neighbors is $1-\prod\limits_{i=1}^{m}(1-\bar{p}_{i})$. An infected node recovers from every infection with a uniform probability $r=0.5$. Having two different diseases, we denote by, respectively, the first type, the second type and the doubly infected type A, B and AB. We assume that susceptible nodes are infected by diseases A and B due to its interaction with an individual from class A or B with probability $p$. We also assume that a susceptible node is infected by disease A or B due to its interaction with an individual from class AB with probability $1-(1-p)^{a}$. We further assume that individuals in class A (resp., B) are infected with the complementary disease by individuals from classes B and AB (resp., A and AB) with probabilities $1-(1-p)^{C}$ and $1-(1-p)^{bC}$. The basic reaction scheme underlying the stochastic model, comparing with the well-mixed model, is given in Table 1.

For $a=b=C=1$, the system corresponds to two noninteracting infections. The case with $C>1$, $a=b=1$, corresponding to a partially cooperative coinfection model, has been addressed in Ref. Chen2016 , where cooperative coinfection exhibits hysteresis of $C$ type.

We simulate the stochastic dynamics in an Erdos-Renyi network of size $N=2^{14}$ with average degree $\langle k\rangle=4$. In our simulations, the primary infection probability $p$ is varied between $p_{0}$ and $p_{f}$, both with increasing and decreasing $p$. The initial condition for the simulations at $p=p_{0}$ when increasing $p$ and at $p=p_{f}$ when decreasing $p$ is set to be the same. Specifically, we set it with a fraction $\epsilon=10^{-4}$ doubly infected that is in state AB and the rest of the nodes in susceptible state. When the stationary value is reached for every value of $p$, we measure the fraction of infected nodes $1-S$ and then we use that final state for $p$ as an initial condition for next value $p+\Delta p$. However, this approach has absorbing states for individual diseases, which are states from which you cannot leave when one of the diseases has no infected nodes. In that case, for the next value of $p$, we keep final condition for the active disease as the corresponding initial condition and set the initial condition for the disease that reached the absorbing state with a fraction $\epsilon$.

We set the values $b=2$ and $a=1$ for our parameters, assessing how the dynamics is influenced by the changes in $C$. To avoid including in our statistics realizations in which one disease is in the absorbing state, which could bias our results towards the results of one single disease spreading, we do not include points in which this happens. First of all, as $C$ is increased, a discontinuous transition appears, where the transition point is smaller, with a higher transition gap, for higher values of $C$; see Fig. 7. Secondly, the jump for increasing $p$ happens at a higher value of $p$ than in the case of decreasing it; see Fig. 7 inset. This leads to hysteresis of $C$ type, while we do not find the other two types reported for the well-mixed dynamics.

IV CONCLUSION

In this paper, we propose a fully cooperative coinfection model, which is investigated through well-mixed and network-based approaches. For the former approach, we prove analytically that the model exhibits three types of hysteresis provided that the model is relatively nonuniform in its parameters. However, we were able to identify just one type of hysteresis in the simulations on random networks.

It is of interest for future work to study the effect of vaccination and awareness on such fully cooperative coinfection model. Moreover, finding the topologies in which $S_{r}$ and $S_{l}$ hysteresis appear will underline the microscopic mechanisms that underscore this dynamics.

Appendix A Proofs of Some Claims in Section II

In Appendix A, we provide the rigorous proof for some claims in the main content. We begin with proving that the solutions $(S,\alpha)$ of $f(S,\alpha)=0$ with $0<S<1$ implicitly defines a function $\bar{\alpha}(s)$ from to $(0,1/S]$ such that $f(S,\bar{\alpha}(S))=0$.

Proposition A.1**.**

Fix $a,b\geq 1$ and $C>0$. Then, for any $0<S\leq 1$, there exists a unique $\bar{\alpha}\,(=\bar{\alpha}(S))$ on $(0,1/S]$ such that $f(S,\bar{\alpha}(S))=0$.

Proof.

Fix $a,b\geq 1$ and $C>0$. We first consider the case $0<S<1$. Writing $f(S,\alpha)$ as $g_{S}(\alpha)$, we have that $g_{S}(0)=-1$ and $g_{S}(1/S)=a^{2}C(1-S)/S>0$ for $0<S<1$. Hence, there must exist a $\bar{\alpha}$ on $(0,1/S)$ such that $g_{S}(\bar{\alpha})=0$. Suppose that there exist two roots $\alpha_{1}$ and $\alpha_{2}$ of $g_{S}(\alpha)$ on $(0,1/S)$. Then, since $g_{S}(0)=-1<0$ and $g_{S}^{\prime}(0)=-2(a-1)S\leq 0$, there exists a critical point on $[0,\alpha_{1})$. Moreover, another critical point should exist on $(\alpha_{1},\alpha_{2})$. Noting $g_{S}(1/S)>0$ for $0<S<1$, we also conclude that a third critical point exists on $(\alpha_{2},1/S)$, a contradiction to the fact that $g_{S}(\alpha)$ is a cubic polynomial. Suppose there exist three roots on $(0,1/S)$. Upon using the fact that $g_{S}(0)<0$, we conclude that $g_{S}(\alpha)$ must be increasing on $(-\infty,0)$, a contradiction to the fact that $g_{S}^{\prime}(0)\leq 0$. Therefore, $g_{S}(\alpha)$ has a unique root on $(0,1/S)$ for $0<S<1$. For $S=1$, the assertion of the proposition can be verified directly. ∎

We next show that the graphs of $\bar{\alpha}(S)$ and $\tilde{\alpha}(S)$ defined in (5) have the same shape on the interval of $(0,1]$ in the following sense. Note that the transformation given in (4) defines a bijective mapping between $(S,\bar{\alpha}(S))$ and $(T,\tilde{\alpha}(T))$.

Proposition A.2**.**

For any $0<t<1$, we have that

[TABLE]

Here function $\text{sgn}(x)$ takes the sign of $x$.

Proof.

By the chain rule, we have that, for any $0<t<1$,

[TABLE]

To complete the proof, it suffices to show that $\frac{d}{dT}\left[\frac{T}{\tilde{\alpha}(T)}\right]>0$ for all $0<T<1$. Indeed, some direct computation yield that

[TABLE]

where

[TABLE]

We have just completed the proof of the proposition. ∎

To see the shape of the graph of $\tilde{\alpha}(T)$, we first compute the derivatives of $\tilde{\alpha}(T)$ at its two end points, [math] and $1$.

Proposition A.3**.**

The following assertions hold.

1. (i)

$\lim_{T\rightarrow 0^{+}}\tilde{\alpha}(T)=\infty$ and $\tilde{\alpha}(1)=1$. 2. (ii)

$\lim_{T\rightarrow 0^{+}}\tilde{\alpha}^{\prime}(T)=-\infty$ and $\tilde{\alpha}^{\prime}(1)=1-\frac{2}{aC}$. 3. (iii)

$\lim_{T\rightarrow 0^{+}}\tilde{\alpha}^{\prime\prime}(T)=\infty$ and $\tilde{\alpha}^{\prime\prime}(1)=-\frac{4}{a^{2}C}[b-a-1/2]$. 4. (iv)

$\tilde{\alpha}^{\prime\prime\prime}(T)<0$ for all $0<T\leq 1$ and $a\geq 1$.

Proof.

Some direct computations yield that

[TABLE]

where

[TABLE]

Then the assertions of parts (i), (ii) and (iii) in the proposition can be verified directly. To see (iv), we rewrite $M_{3}(T,a,b)$ as the polynomial form of $b$. Specifically, we have that

[TABLE]

where $c_{3}=-(1-T)^{4}\,(<0)$, $c_{2}=-2aT[(a-2)T^{3}+6T^{2}-6T+2]$, $c_{1}=a^{2}T^{2}[(4a-5)T^{2}+12T-6]$, $c_{0}=-2a^{3}T^{3}[(a-1)T+2]\,(<0)$. Hence, to conclude the assertion (iv) in the proposition, we only need to show that $c_{1}^{2}-4c_{0}c_{2}<0$ for all $0<T\leq 1$ and $a\geq 1$. Indeed, we compute that $c_{1}^{2}-4c_{0}c_{2}=a^{4}T^{4}[(8a-7)(1-T)^{4}-12(1-T)^{3}-6(1-T)^{2}-4(1-T)+(1-8a)]\leq a^{4}T^{4}[(8a-7)+(1-8a)]=-6a^{4}T^{4}<0$. The proof is now complete. ∎

Proposition A.4**.**

For any $a\geq 1$ and $b\leq b_{c}(:=a+1/2)$, the following two assertions hold.

1. (i)

If $C\leq C_{2}$, then $\tilde{\alpha}(T)$ is a strictly decreasing function of $T$ on $(0,1]$. 2. (ii)

If $C>C_{2}$, then $\tilde{\alpha}(T)$ has a unique minimum on $(0,1)$. In particular, the graph of $\tilde{\alpha}(T)$ is $C$-shaped.

Here $C_{2}=\frac{2}{a}$ is define as in (6).

Proof.

Since $b\leq b_{c}$, we have that $\tilde{\alpha}^{\prime\prime}(1)\geq 0$ by Proposition A.3(iii). Consequently, $\tilde{\alpha}^{\prime\prime}(T)>0$ for all $0<T<1$ by Proposition A.3(iv), which implies that $\tilde{\alpha}(T)$ is concave up on $(0,1]$. In addition, $\tilde{\alpha}(T)$ is strictly decreasing on $(0,1]$ if $\tilde{\alpha}^{\prime}(1)\leq 0$, while first decreasing and then increasing on $(0,1]$ if $\tilde{\alpha}^{\prime}(1)>0$. Above cases correspond to $C\leq C_{2}$ or $C>C_{2}$, respectively, by Proposition A.3(ii). Hence, the proof is complete. ∎

For $a\geq 1$ and $b>b_{c}$, the graph of $\tilde{\alpha}(T)$ is more complicated since the concavity of $\tilde{\alpha}(T)$ changes from the upward to the downward on $(0,1)$ by Proposition A.3(iii) and (iv). More precisely, $\tilde{\alpha}(T)$ has an inflection point $\gamma$ on $(0,1)$ such that

[TABLE]

To classify its shape, we first find $C$ so that the corresponding $\tilde{\alpha}(T)$ has the some particular properties.

Proposition A.5**.**

For any $a\geq 1$ and $b>b_{c}\,(:=a+1/2)$, define $\gamma$ as in (A12). Then the following three assertions hold.

1. (i)

There exists a unique positive $C_{1}$, defined as in (7), such that

[TABLE] 2. (ii)

There exists a unique positive $C_{12}$, defined as in (8), such that the equation

[TABLE]

has a unique solution on $(0,1)$ when $C=C_{12}$. Moreover, (A14) has no solution on $(0,1)$ when $C<C_{12}$. 3. (iii)

It holds that $C_{1}<C_{12}<C_{2}$. Here $C_{2}=\frac{2}{a}$ is define as in (6).

Proof.

Fix $a\geq 1$ and $b>b_{c}(:=a+1/2)$. We first show Proposition A.5(i). Since $\gamma$ is the inflection point of $\tilde{\alpha}$, it holds that $\tilde{\alpha}^{\prime\prime}(\gamma)=0$. We have, via (A11), that $M_{2}(\gamma,a,b)=0$. Moreover, since function $M_{2}$ is independent of $C$, we have that $\gamma$ is also independent of $C$. On the other hand, note that

[TABLE]

Clearly, $\alpha_{1}(T)$ and $\alpha_{0}(T)$ are positive for all $0<T<1$. By choosing $C_{1}=\frac{\alpha_{0}(\gamma)}{\alpha_{1}(\gamma)}$, we have that

[TABLE]

In fact, we can obtain the explicit form of $C_{1}$, as given in (7) by using the Maple.

We next show Proposition A.5(ii). By solving $\tilde{\alpha}(T)=1$ directly, we have that $T=1$ or

[TABLE]

where

[TABLE]

Hence, for equation $\tilde{\alpha}(T)=1$ to have exactly one solution on $\mathbb{R}$, we must have $L(a,b,C)=0$ or, equivalently,

[TABLE]

We are then able to derive the following results.

1. (a)

For $C=C_{-}$, we have that $\hat{T}=\frac{(abC_{-}-2a+1)}{2a(b-a)C_{-}}<0$ by some tedious computation. It implies that $\tilde{\alpha}(T)$ is either greater or less than $1$ on the entire interval $(0,1)$. In fact, $\tilde{\alpha}(T)>1$ on $(0,1)$ since $\lim_{T\rightarrow 0^{+}}\tilde{\alpha}(T)=\infty$ by Proposition A.3(i). 2. (b)

For $C<C_{-}$, we have, via direct computation, that

[TABLE]

whenever $0<T<1$. This implies that $\left[\tilde{\alpha}(T)\right]_{C}>\left[\tilde{\alpha}(T)\right]_{C=C_{-}}>1$ for $C<C_{-}$. 3. (c)

For $C\in(C_{-},C_{+})$, since $L(a,b,C)<0$, equation $\tilde{\alpha}(T)=1$ has no solution on $\mathbb{R}$ and hence $[\tilde{\alpha}(T)]_{C}>1$ on $(0,1)$. 4. (d)

For $C=C_{+}$, we have that $\hat{T}=\frac{(abC_{+}-2a+1)}{2a(b-a)C_{+}}\in(0,1)$ by some tedious computation. Moreover, since $[\tilde{\alpha}(T)]_{C}>1$ on $(0,1)$ for $C<C_{+}$ as claimed above, we have that $[\tilde{\alpha}(T)]_{C_{+}}\geq 1$ on $(0,1)$ and the equality holds if and only if $T=\hat{T}$.

Hence, the proof of Proposition A.5(ii) is complete by choosing $C_{12}=C_{+}$.

Finally, we show Proposition A.5(iii). As in the proof of part (ii), we have that $\tilde{\alpha}(T)>1$ on $(0,1)$ whenever $C<C_{12}$; while $\tilde{\alpha}(T)\geq 1$ on $(0,1)$ and there exists $T_{*}$ on $(0,1)$ such that $\tilde{\alpha}(T_{*})=1$ whenever $C=C_{12}$; $\tilde{\alpha}(T)<1$ for some point on $(0,1)$ whenever $C>C_{12}$ by (A14) and (A17). Then, by (A12) and (A13), we have that, when $C=C_{1}$, $\tilde{\alpha}(T)$ is strictly decreasing on $(0,1)$ and hence $\tilde{\alpha}(T)>\tilde{\alpha}(1)\,(=1)$ on $(0,1)$. So $C_{1}<C_{12}$. On the other hand, since $\tilde{\alpha}^{\prime}(1)=0$ by Proposition A.3(ii) and $\tilde{\alpha}^{\prime\prime}(1)<0$ by (A12), we have that $\tilde{\alpha}^{\prime}(T)<0$ and hence $\tilde{\alpha}(T)<\tilde{\alpha}(1)\,(=1)$ whenever $T$ is sufficiently close to $1^{-}$. So $C_{2}>C_{12}$.

We have just completed the proof of Proposition A.5. ∎

Proposition A.6**.**

For any $a\geq 1$ and $b>b_{c}(:=a+1/2)$, the following two assertions hold.

1. (i)

If $C<C_{1}$, then $\tilde{\alpha}(T)$ is strictly decreasing. 2. (ii)

If $C\in(C_{1},C_{12})$, then $\tilde{\alpha}(T)$ has a unique local minimum and a unique local maximum on $(0,1)$. The graph of such $\tilde{\alpha}(T)$ is to be called an $S$-shaped. Moreover, we have that $\tilde{\alpha}(T)>1$ for all $0<T<1$. 3. (iii)

If $C\in(C_{12},C_{2})$, then the graph of $\tilde{\alpha}(T)$ is also $S$-shaped with $\tilde{\alpha}(T)<1$. for some $T$ on $(0,1)$. 4. (iv)

If $C>C_{2}$, then $\tilde{\alpha}(T)$ has a unique critical point. The graph of such $\tilde{\alpha}(T)$ is to be termed a $C$-shaped.

Proof.

We first prove Proposition A.6(i). By (A12) and (A13), we have that, when $C=C_{1}$, $\tilde{\alpha}^{\prime}(T)\leq 0$ on $(0,1)$. Then, by (A15), we have that $\tilde{\alpha}^{\prime}(T)<0$ on $(0,1)$ if $C<C_{1}$. So the assertion of part (i) holds.

We next prove Proposition A.6LABEL:{item:case_b_larger_b_c_and_C_in_C_1_C_12} and (iii). For $C\in(C_{1},C_{2})$, we have that $\tilde{\alpha}^{\prime}(\gamma)>0$ by (A16) and $\tilde{\alpha}^{\prime}(1)<0$ by Proposition A.3(ii). Then by (A12), we conclude that $\tilde{\alpha}(T)$ has a unique local minimum and a unique local maximum on $(0,1)$. The remaining parts of the proof have been claimed in the proof of Proposition A.5(iii). So the assertions of part LABEL:{item:case_b_larger_b_c_and_C_in_C_1_C_12} and (iii) hold.

Finally, we prove Proposition A.6(iv). In fact, the assertion holds by (A12) and since $\tilde{\alpha}^{\prime}(\gamma)>0$ by (A16) and $\tilde{\alpha}^{\prime}(1)>0$ by Proposition A.3(ii).

We have just completed the proof of Proposition A.6. ∎

Acknowledgements.

J.P.R. acknowledges support from Ministry of Science and Technology (Taiwan) and the FPU program of MECD (Spain), and thanks Víctor M. Eguíluz and Fakhteh Ghanbarnejad for comments and discussion. Y.H.L. and J.J. acknowledge support from the Ministry of Science and Technology of Taiwan under grant No. Most 105-2115-M-009-002-MY2.