Dynamics of interacting diseases

Current modeling of infectious diseases allows for the study of complex and realistic scenarios that go from the population to the individual level of description. However, most epidemic models assume that the spreading process takes place on a single level (be it a single population, a meta-population system or a network of contacts). In particular, interdependent contagion phenomena can only be addressed if we go beyond the scheme one pathogen-one network. In this paper, we propose a framework that allows describing the spreading dynamics of two concurrent diseases. Specifically, we characterize analytically the epidemic thresholds of the two diseases for different scenarios and also compute the temporal evolution characterizing the unfolding dynamics. Results show that there are regions of the parameter space in which the onset of a disease's outbreak is conditioned to the prevalence levels of the other disease. Moreover, we show, for the SIS scheme, that under certain circumstances, finite and not vanishing epidemic thresholds are found even at the thermodynamic limit for scale-free networks. For the SIR scenario, the phenomenology is richer and additional interdependencies show up. We also find that the secondary thresholds for the SIS and SIR models are different, which results directly from the interaction between both diseases. Our work thus solve an important problem and pave the way towards a more comprehensive description of the dynamics of interacting diseases.


I. INTRODUCTION
The problem of how diseases spread has been studied for a long time [1][2][3][4][5][6][7][8]. In the past few years, due to the increasing availability of data about transmission patterns, contact networks, and population mobility, it has become apparent that we are in a position to develop theoretical and computational frameworks that will ultimately allow the forecast of epidemic outbreaks [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. The latter is a consequence of the increasing availability of new models that are capable of providing better predictions of real epidemic scenarios [27]. In particular, we have been able to identify what ingredients are key to characterizing the unfolding of a disease outbreak. This is the case of the structure of the population in the form of networks of contacts at a local scale [9,17,19,28], or of individual mobility patterns that facilitate the spreading of diseases in wider geographical areas [29]. The theoretical implications of considering these aspects have been thoroughly addressed [30][31][32][33], and at the same time, the amount and quality of real data that are relevant to epidemic spreading are constantly increasing [34,35].
In this context, an emergent field of research is the modeling of coupled spreading phenomena, whether they are two pathogens or multiple strains of the same disease that propagate concurrently on the same population [33,[36][37][38][39]. Focusing on the two-pathogen scenario, the complexity of the problem increases because now the natural history of one of the diseases is affected by the presence of the second one, typically as a consequence of the modification of the host's immune response after infectionwith a plethora of possibilities as given by different interaction schemes. In addition, the networks of contacts through which the pathogens spread can vary from one disease to another. Typical examples of these coupled spreading phenomena are given by the interaction between HIV infection and the spreading of certain opportunist Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. pathogens [40,41], or the strain-strain interactions of a viral pathogen like the flu virus [42].
From a theoretical point of view, one of the first works that considered the above problem from a networked point of view is from Funk and Jensen [37]. In their work, they exploited the analogy of the spreading process with a bondpercolation problem so as to address the issue of epidemic thresholds in a susceptible-infected-removed (SIR) model. However, this first work did not consider a framework in which the temporal evolution of the epidemics could be studied. Following a similar bond-percolation approach, Mills et al. [43], motivated by the interaction between AIDS and tuberculosis, studied two coupled SIR models where the diseases can spread in both layers at different rates. Also in this case, due to the static nature of the bondpercolation approach, the temporal evolution of the diseases cannot be studied. Recently, Marceau et al. presented a new model [39] aimed at studying the latter aspect, also within the SIR scheme, using "on-the-fly graphs," a network generation model they had previously introduced [44]. In this work, although the temporal evolution of the system can be monitored, the modeling approach does not provide information about epidemic thresholds. In addition, the model cannot be trivially generalized so as to cover other classical models like susceptible-infected-susceptible (SIS). Other works have considered the problem by studying the spreading of the diseases as a Markovian process. In Ref. [45], Sahneh and Scoglio, using two coupled Markov processes, studied the competition of two epidemics and found different regions of the parameter space where the two diseases can coexist. However, the model in Ref. [45] represents only the specific scenario of two mutually excluding diseases. It is also worth mentioning that, using a similar approach, the so-called microscopic Markov-Chain [25], the specific case in which a disease and an "awareness" diffusion process coexist has been studied in Ref. [46].
In this work, we propose a modeling framework, based on a heterogeneous mean-field (HMF) approximation, for the spreading of two concurrent diseases that interact with each other. With this purpose, we analyze in detail the dynamics of a two-layer system that represents the networks of contacts on top of which either two coupled SIS or SIR processes spread. Although we use as a convenient representation of the system under study a topology that has recently attracted a lot of attention, i.e., the so-called multilayer networks [47,48], we stress that the main focus is to determine what the effects are of the complex interplay between the different interaction mechanisms considered and the temporal and topological scales onto the dynamical behavior of the diseases. Specifically, we explicitly derive the epidemic threshold of each disease in terms of the parameters that characterize their evolution and the topology of the networks of contacts. Remarkably, we find that the epidemic thresholds in the SIR case are different from those of the SIS case, which is a consequence of the emergence of new mechanisms of interaction and of the transitory nature of the epidemic outbreaks in the SIR case. Additionally, different scenarios and limiting cases of both theoretical and epidemiological interest are scrutinized. Finally, results from numerical simulations are presented to validate our analytical results. Our work, thus, shows how the different interaction mechanisms considered can give rise to new phenomenological insights regarding the dynamics of interacting diseases.

II. MODELING FRAMEWORK
Our purpose in this work is to explore epidemiological scenarios in which two different infectious diseases simultaneously spread through the same host population, and whose dynamical parameters (i.e., infectiousnesses and recovery rates) depend, at the level of single individuals, on whether the agents involved have caught only one of the diseases or both. In addition, we consider that neither the mechanisms behind each disease evolution nor the contact networks through which they spread have to be the same, and, thus, are considered independently for each disease. Let us then consider that we have two diseases (disease 1 and disease 2) that spread over two different networks of contacts: disease 1 propagates over network 1, which has a mean degree hki ¼ P k;l Pðk; lÞk, while disease 2 propagates over network 2, whose mean connectivity is equal to hli ¼ P k;l Pðk; lÞl. The composed degree distribution Pðk; lÞ gives the proportion of nodes (individuals) having k and l links in networks 1 and 2, respectively [49]. In what follows, we consider both the SIS and SIR schemes independently.

A. SIS scenario
As a first step, we consider the baseline scenario in which the isolated dynamics of each disease, when the second is absent, is described by a simple SIS scheme. Each individual belonging to a composed connectivity class ðk; lÞ can be in four different dynamical states: susceptible with respect to both diseases SSðk; lÞ, infected of both IIðk; lÞ, and infected with the first [second] one and still susceptible to catch the second [first] disease, ISðk; lÞ [SIðk; lÞ], SSðk; lÞ; ISðk; lÞ; SIðk; lÞ; IIðk; lÞ represent the proportion of individuals in each disease state within the composed degree class ðk; lÞ. Thus, we have that SSðk; lÞ þ ISðk; lÞ þ SIðk; lÞ þ IIðk; lÞ ¼ 1 ∀ðk; lÞ. In addition, regardless of the connectivities of the nodes involved, we have eight possible contagion transitions after a contact (four for each disease). On the other hand, four other possible transitions correspond to the cases in which infected individuals go back to the susceptible class. This amounts to a total of 12 elementary transitions, schematized as follows (see also Fig. 1): The model thus contains two basic infection probabilities, λ 1 and λ 2 , as well as two basic recovery rates, μ 1 and μ 2 , one for each disease. In addition, infection probabilities are affected by scaling factors-the four β's and combinations of them-in the same way that the recovery rates are affected by the η's coefficients, as we explain in detail in Table I. These parameters describe the interaction of the diseases through three different effects that are concurrently taken into account. The first effect is the variation of the susceptibility of healthy individuals to get infected with one disease as a consequence of being infected with another. This mechanism is described by β a 1 for the variation of the infection risk of disease 1 caused by disease 2, and β a 2 for the symmetric case. The second effect is the change of the spreading capabilities of double-infected individuals with respect to single-infected ones, which is described by the parameters β b 1 and β b 2 for diseases 1 and 2, respectively. The last effect is the variation of the infectious periods of double-infected individuals also with respect to singleinfected ones, described by η 1 and η 2 .
Therefore, these parameters exhaustively describe all the ways in which two diseases can interact according to a SIS scheme and allow us to isolate the different effects of one disease on the spreading of the other by making infectious individuals more efficient spreaders or by making susceptible individuals more prone to get sick. Once we have defined the whole set of parameters in Table I, we have all possible transitions between dynamical states well defined (see Fig. 1).

Epidemic thresholds
In order to analyze the most relevant dynamical properties of the system, we look for the values [IS Ã ðk; lÞ; SI Ã ðk; lÞ; II Ã ðk; lÞ] that define the stationary state: ½ _ IS Ã ðk; lÞ; _ SI Ã ðk; lÞ; _ II Ã ðk; lÞ ¼ ð0; 0; 0Þ ∀ðk; lÞ for the system Eq. (8). In order to do that, we are forced to consider σ 1 and σ 2 as additional parameters, although these two quantities are linear combinations of the other variables and are ultimately responsible for the coupling among all connectivity classes. As commonly happens in this kind of model [6,17], there exists a trivial fixed point of the dynamics in which there are no infected individuals in the system: ½II Ã 1 ðk;lÞ;II Ã 2 ðk;lÞ;II Ã ðk;lÞ ¼ ð0;0;0Þ ∀ðk; lÞ. This fixed point represents the absorbing state of our model. In addition, there are other possible fixed points for which the densities of infected individuals could be written as a function of the σ parameters: ½II Ã 1 ðk; l; σ 1 ; σ 2 Þ; II Ã 2 ðk; l; σ 1 ; σ 2 Þ; II Ã ðk; l; σ 1 ; σ 2 Þ ≠ ð0; 0; 0Þ. It is thus possible to get self-consistent equations for the variables σ as Pðk; lÞk½ISðk; l; σ 1 ; σ 2 Þ þ β b 1 IIðk; l; σ 1 ; σ 2 Þ; ð9Þ Pðk; lÞl½SIðk; l; σ 1 ; σ 2 Þ þ β b 2 IIðk; l; σ 1 ; σ 2 Þ: ð10Þ The condition σ 1 ¼ f 1 ðσ 1 ; σ 2 Þ > 0 implies the existence of a stable active state for the dynamics of disease 1, Baseline recovery rate of disease 1 μ 2 Baseline recovery rate of disease 2 β a 1 Variation of disease 1 infectiousness due to the fact that the susceptible individual exposed to disease 1 is infected with disease 2 β a 2 Variation of disease 2 infectiousness due to the fact that the susceptible individual exposed to disease 2 is infected with disease 1 β b 1 Variation of disease 1 infectiousness due to the fact that the spreader is also infected with disease 2 β b 2 Variation of disease 2 infectiousness due to the fact that the spreader is also infected with disease 1 η 1 Variation of disease 1 recovery rate for individuals also infected with disease 2 η 2 Variation of disease 2 recovery rate for individuals also infected with disease 1 Probability that a given link of network 2 points to an II node i.e., a state in which disease 1 becomes endemic in the population. For this situation to take place for disease 2, the condition σ 2 ¼ f 2 ðσ 1 ; σ 2 Þ > 0 must be fulfilled. Given the symmetry between the two expressions in Eq. (10), we focus only on the analysis of the first equation. In fact, it can be shown that ∂ 2 fðσ 1 ;σ 2 Þ ∂σ 2 1 < 0 always, which means that, if we think of the graphical solution of the equation σ 1 ¼ f 1 ðσ 1 ; σ 2 Þ, for this nontrivial solution to exist-given that, obviously, fðσ 1 ¼ 0; σ 2 ¼ 0Þ ¼ 0-it must be verified that ½ ∂fðσ 1 ;σ 2 Þ ∂σ 1 σ 1 ¼0 > 1, as in Ref. [17]. After some algebraic operations, this condition yields the following expression: which allows us to derive the epidemic threshold as Looking at the latter expression-which contains the underlying topologies in a more intricate way than for the uncoupled, classical case-the threshold dependence on disease 2's prevalence via σ 2 becomes explicit. If we evaluate λ c 1 ðσ 2 ¼ 0Þ, we recover the classical result λ c 1 ¼ μ 1 hki=hk 2 i [9,17]. Therefore, in the following we refer to this baseline case as primary threshold, λ c 1 ð0Þ, whereas the more general case is referred to as the secondary threshold, λ c 1 ðσ 2 Þ (with σ 2 > 0). Obviously, the same stands for the primary [λ c 2 ð0Þ] and secondary [λ c 2 ðσ 1 Þ] thresholds of the second disease. A particular case for the topologies on top of which both diseases are spreading corresponds to the homogeneous mean-field version of the system, i.e., Pðk; lÞ ¼ δðk − k o Þδðl − l o Þ, for which the last expression can be rewritten as An independent derivation of this expression can be obtained by analyzing the Jacobian matrix of the homogeneous mean-field system analogous to Eq. (8), as shown in Appendix A.

Phase diagrams
In order to explore the quality of the threshold prediction of our model, we designed a Monte Carlo simulation scheme in which a single state transition is allowed per individual per time step. First, infected individuals will eventually spread the disease(s) that they carry. As double events are not allowed at a single time step, forbidden double transitions from SS to II are resolved by choosing the disease that an individual will catch proportionally to the status of their infected neighbors: the more neighbors one individual has, say, carrying disease 1, the more likely the individual will become infected with disease 1 rather than with disease 2. After the spreading loop is completed for both diseases, infected nodes that have not suffered any contagion at the present time step eventually get back to the susceptible state of the disease(s) they carry. To avoid forbidden double transitions from the II class to the SS class, in those cases the only disease the individual is going to recover from is also chosen stochastically, according to the probabilities p 1 ¼ η 1 μ 1 =ðη 1 μ 1 þ η 2 μ 2 Þ for the first disease and p 2 ¼ 1 − p 1 for the second disease.
We first explore a simple scenario in which we assume that the dynamical effects of one disease on the other are totally symmetric. In terms of the parameters of the model, this implies that β a In this case, we focus on two opposite scenarios: (i) mutual enhancement, η < 1 and β > 1, and (ii) mutual impairment, η > 1 and β < 1. In the case of mutual enhancement, individuals who are infected with the second disease spread and become infected with disease 1 more easily than those who are not (this is because β > 1). In addition, since also η < 1, individuals infected with disease 1 remain so for longer times if they are also infected with the second disease. These two effects imply that the appearance of disease 2 in the system enhances the spreading capabilities of disease 1. The reciprocal situation is also true, as the interaction between the two diseases is symmetric. Finally, in the case of mutual impairment, the effects on the infectiousness and recovery rates are the opposite, and so the appearance of one of the diseases at a certain prevalence impairs the spreading of the other disease. In Figs. 2 and 3, we represent the prevalence of each disease as a function of the baseline infectiousnesses ðλ 1 ; λ 2 Þ for a given set of parameters after the introduction of infection seeds in the order shown. The networks through which diseases spread are, for this first case, two uncorrelated Erdös-Renyi (ER) graphs.
The color maps represent the prevalence levels of disease 1 (upper panels) and disease 2 (lower panels) at different stages of the Monte Carlo simulations. As in Fig. 2, we successively introduce three infectious seeds ðIS; SI; ISÞ ¼ ð0.005; 0.005; 0.005Þ and plot the stationary prevalences after each fluctuation in the three columns of the figure. As can be seen, the reintroduction of the third seed of infection 1 in the system does not affect the prevalence levels, as global stability is reached after the second stage. For the case of mutual enhancement (Fig. 2), given our set of parameters, the analytically obtained curves for the secondary threshold remain below the primary thresholds, leading to the appearance of two regions in the plane ðλ 1 ; λ 2 Þ, for which it is verified that λ c 1 ðσ 2 Þ < λ 1 < λ c 1 ð0Þ and λ c 2 ðσ 1 Þ < λ 2 < λ c 2 ð0Þ, respectively. The dynamical relevance of these regions is that within them the appearance of an outbreak of one of the diseases is conditional to the previous installation of the other infection in the system. In this way, we observe that, after an initial seed of IS individuals, disease 1 does not become endemic in the region λ c Fig. 2(a)], but then, after the outbreak of the second disease in the network, the same seed leads disease 1 to become endemic in that same region [ Fig. 2(c)], as predicted by our model. Regarding the conjugate region in which λ c 2 ðσ 1 Þ < λ 2 < λ c 2 ð0Þ, we see in Fig. 2(e) that disease 2 directly becomes endemic after the introduction of an infection seed SI due to the fact that, previously, disease 1 was already introduced in the system.

FIG. 2. Reciprocally enhanced spreading. Dynamical parameters
The situation for the mutual impairment case is the opposite, and the secondary thresholds remain, in this case, above the primary ones. So, in this scenario, we have two more relevant regions in which λ c 1 ð0Þ < λ 1 < λ c 1 ðσ 2 Þ and λ c 2 ð0Þ < λ 2 < λ c 2 ðσ 1 Þ, respectively. In Fig. 3, we show the behavior of the system under these conditions. In Fig. 3(a), we see how disease 1 becomes endemic after the introduction of an initial seed above its primary threshold. Then, after introducing a seed of disease 2, as shown in Fig. 3(b) for the area comprised between λ c 1 ð0Þ and λ c 1 ðσ 2 Þ, the prevalence of disease 1 vanishes. In other words, in that region, the introduction of disease 2 makes it possible for the system to recover from disease 1. If we look at the behavior of the second disease in the region in which λ c 2 ð0Þ < λ 2 < λ c 2 ðσ 1 Þ, we see how the disease is unable to become endemic as a consequence of the fact that the first disease has already been introduced in the system. This situation suggests that, as has already been addressed in the context of computational sciences [50], the introduction of an infectious agent designed to immunize its host with respect to another, more harmful infection might be a conceptually feasible option to reduce the prevalence of the latter, or even to eradicate it. This has also been recently reported in the context of multistrain diseases, in which more than one strain of the same disease competes for the host population [33,45].
Once the dynamics of the model has been exhaustively characterized when the diseases spread over homogeneous networks, we move on and explore the influence of degree heterogeneity on the dynamics. To this end, we also perform intensive numerical simulations in an analogous way, but using scale-free (SF) graphs of the same size as before with different exponents. In , also at the final state. Scale-free networks are generated using the uncorrelated configuration model with N 1 ¼ N 2 ¼ 5000 agents, hki ¼ 4.00, and hli ¼ 5.11. The figure represents the final prevalence of (a),(b) disease 1 and (c),(d) disease 2 in each case.
configurations, secondary thresholds are closer to primary ones than in the case of homogeneous networks.

System sizes and epidemic thresholds: General case
In the previous sections, we described the baseline cases in which none of the dynamical parameters vanish, and both networks are, in each case, of the same kind-ER or uncorrelated SF graphs. A relevant theoretical question remains unanswered, i.e., how the epidemic thresholds behave when the system size grows and eventually reaches the thermodynamic limit. Regarding this question, we present an exhaustive analysis in Appendix B, in which we show the conditions that lead to having vanishing small secondary thresholds as a function of the underlying topologies and some of the dynamical parameters; note, however, that this question is of interest from a theoretical viewpoint, as, strictly speaking, all real systems are finite and, thus, an effective epidemic threshold exits. In addition, our approach is based on a HMF approximation, and that means that influences of dynamical correlations on the analysis are neglected.
The results show that the secondary epidemic threshold associated with any of the diseases that is spreading over a scale-free network with 2 < γ ≤ 3 vanishes at the thermodynamic limit, regardless of the topology of the network on top of which its conjugate disease propagates and regardless of the values of the dynamical parameters. In an analogously robust way, power laws with γ > 3-or homogeneous degree distributions-yield finite, nonvanishing secondary thresholds at the thermodynamic limit, regardless of the conjugate topologies or parameter values, with some exceptions.
The relevance of this result relies on the fact that the model predicts the same behavior for the epidemic thresholds in heterogeneous and homogeneous networks as compared with HMF models of uncoupled (single) diseases that spread over simple networks. In addition, our analysis shows that the eventual vanishing of the epidemic thresholds for infinite systems is determined only by the topology of the network under consideration rather than by any possible coupling with another disease that spreads over any other possible conjugate network within our model framework. However, there is an exception to this general behavior that is meaningful from an epidemiological viewpoint. This is the case when both diseases spread over two highly and positively correlated scale-free networks with composed degree distribution Pðk; lÞ ¼ δðk − lÞαk −γ , where δ stands for the Kronecker δ function and α is a normalization constant. In that situation, if we focus, for example, on disease 1, there exist two different interaction schemes of interest for which we recover finite epidemic thresholds λ c 1 > 0 at the thermodynamic limit even for scale-free graphs with 2 < γ ≤ 3: (i) Case 1: Individuals infected with disease 2 become immune to infection by disease 1. β a If β a 1 ≠ 0, individuals infected with both diseases cannot cause contagion of disease 1: β b 1 ¼ 0. In addition, disease 1 cannot cause total immunity to disease 2: β a 2 ≠ 0. β b 2 is a free parameter.
To illustrate this situation, we take as an example a particular case of the first scheme, a mutual cross-immunity scenario given by β a 1 ¼ β a 2 ¼ 0. In order to point out the role of interlayer degree correlations on this effect, we can directly compare the expression for the threshold when both networks are totally correlated with the analogous expression derived from an uncorrelated combined degree distribution: In Fig. 5, we represent the values predicted by these expressions for different network sizes, when γ ¼ Γ ¼ 2.5.
As we see in Fig. 5(a), for uncorrelated networks, regardless of the value of σ 2 , the threshold continuously decreases as we increase network sizes. The result thus shows that the existence of a coupling with another disease present in the system with a certain prevalence proportional to σ 2 does not play any role, since the degree heterogeneities are still the main reason leading to the vanishing of the threshold at the thermodynamic limit. This picture turns out to be remarkably different when we introduce positive, strong correlations between the two networks. In that case, as we see in Fig. 5(b), the appearance of the second disease, characterized by a certain prevalence level σ 2 > 0, implies a sudden change in the behavior of the threshold, which no longer vanishes, even when N → ∞.
The influence of degree correlations between networks for this case of full cross-immunity becomes evident also at finite sizes, since the differences between primary and secondary thresholds, as seen in Fig. 6, are also greatly amplified. Another eventually relevant effect that can be observed in Fig. 6 is that the transition that takes place at the epidemic threshold is much sharper in the case of correlated networks. All the previous results point out that the worst scenario for the spreading of a disease when it interacts with a second one that confers immunity to the former corresponds to the case in which there is a correlation between heterogeneous networks of contacts. This finding is essentially equivalent to what was found previously in Ref. [37]; what show here is that the effect comes to revert the vanishing threshold at the thermodynamic limit for a disease spreading on top of scale-free networks. In addition, it is not needed that the second disease causes total immunity against the first for this effect to be present.
Remarkably, all these "pathological" cases can be identified without abandoning HMF descriptions. Beyond HMF, the behavior of epidemic thresholds varies with respect to mean-field predictions as a consequence of dynamical correlations [20,51]. Here, however, we identify that the interactions between diseases can modify the size scaling of thresholds from the classical mean-field picture without the need of dynamical correlations, which, for certain cases, remarkably modify the whole picture at the thermodynamic limit, even in annealed networks in which adjacency matrices are only fixed on average, and so, dynamical correlations do not exist.
FIG. 6. Effect of degree correlations between the two scale-free networks on the steady prevalence levels. Dynamical parameters:

B. SIR scenario
In the previous sections, we presented the behavior of systems of interacting diseases that spread according to a SIS scheme, thus leading, above the epidemic threshold, to stationary, endemic states with a prevalence greater than zero. In the following, we explore the case of transient, interacting epidemic phenomena.
In order to do so, we can extend the framework and describe the dynamics of two SIR epidemics interacting among them. In classical, noninteracting systems-either homogeneous or heterogeneous-the resemblance of both types of models translates into a strong mathematical symmetry between them that yields identical expressions for the epidemic thresholds under mean-field descriptions (i.e., when neglecting the effects of dynamical correlations). In this section, we see the way in which part of this symmetry is broken as a consequence of the interacting nature of the epidemic processes.

Mathematical description
In this case, the first new aspect to note is that not just the condition of being infected with one disease can modify subjects' susceptibility to the conjugate infection, but also whether they have been infected and have recovered. For instance, this might represent situations in which some kind of immunological memory is acquired after the first infection-for example, partial immunity in front of other strains is gained after suffering an influenza infection, especially if both are phylogenetically close enough [52]. This new phenomenology translates into the need of introducing new parameters (see Table III where all the parameters and variables that were present in the SIS formulation retain their original meaning, with the nuance that now σ 1 and σ 2 have an extra term related to the appearance of classes IR and RI:  Variation of disease 1 infectiousness due to the fact that the susceptible individual exposed to disease 1 has been infected and recovered from disease 2 ϕ a 2 Variation of disease 2 infectiousness due to the fact that the susceptible individual exposed to disease 2 has been infected and recovered from disease 1 ϕ b 1 Variation of disease 1 infectiousness due to the fact that the spreader has been infected and recovered from disease 2 ϕ b 2 Variation of disease 2 infectiousness due to the fact that the spreader has been infected and recovered from disease 1 ζ 1 Variation of disease 1 recovery rate for individuals that have been infected and recovered from disease 2 ζ 2 Variation of disease 2 recovery rate for individuals that have been infected and recovered from disease 1 where θ IR 1 is the probability of a link of network 1 pointing to an IR node and θ RI 2 is the probability of a link of network 2 pointing to a node in the RI state.

Epidemic thresholds
In order to characterize the epidemic threshold of disease 2, we need to study the dynamics of the system around any point in which no infected individual of disease 2 has yet been introduced and so ½SIðk; lÞ; IIðk; lÞ; RIðk; lÞ ¼ ð0; 0; 0Þ ∀ðk; lÞ, as well as θ SI Around such a disease-free point, either fixed or not, it is trivial to see that all the partial derivatives of classes SIðk; lÞ; IIðk; lÞ and RIðk; lÞ with respect to the rest of dynamic classes vanish. This makes the subset of variables fSIðk; lÞ; IIðk; lÞ; RIðk; lÞg locally autonomous around that point, and so, its linearization might serve us to address the stability inversion yielding the emergence of the epidemic threshold.
In addition, the dimensionality of the system can be greatly reduced if we write the equations driving the time evolution of the probabilities θ SI 2 , θ II 2 , and θ RI 2 as follows: where we substitute hklSIi by hkliθ SI 2 , an approximation that is valid around the point ½SIðk; lÞ; IIðk; lÞ; RIðk; lÞ ¼ ð0; 0; 0Þ∀ðk; lÞ. Obviously, as happened for fSIðk; lÞ; IIðk; lÞ; RIðk; lÞg, all of the partial derivatives of θ SI 2 , θ II 2 , and θ RI 2 with respect to variables unrelated to θ SI 2 , θ II 2 , and θ RI 2 vanish, which allows us to study the stability of the system with respect to disease 2 by linearizing the The corresponding Jacobian matrix J can be calculated this way, and, by evaluating the condition J ¼ 0 for stability shift, the epidemic threshold takes its final value: which depends on the initial prevalence of disease 1 via σ 1 , hl 2 ISi, hl 2 RSi [and hl 2 SSi]. In the case of having nonconcurrent outbreaks, we have that the second disease arrives at the system after the outbreak of the first disease has come to an end. In that case, we have that ISðk; lÞ ¼ 0 and RSðk; lÞ ¼ R 1 ∞;k;l ∀ðk; lÞ, where R 1 ∞;k;l is the fraction of recovered individuals at the end of an outbreak of disease 1 alone, in the composed degree class ðk; lÞ. In such a case, the problem is much simpler, as σ 1 ¼ hl 2 ISi ¼ 0 and hl 2 SSi ¼ hl 2 i − hl 2 R 1 ∞;k;l i, and the threshold reads as λ c 2 ðhl 2 i; hl 2 R 1 ∞;k;l iÞ ¼ Obviously, if disease 1 has not yet appeared in the system, the threshold for disease 2 becomes λ c 2 ¼ μ 2 hli hl 2 i , as in the classical, noninteracting HMF case [19]. As is done for the SIS model, we refer to the latter expression as the primary threshold of disease 2, in counterposition to the secondary threshold of Eq. (26). A remarkable property of the threshold for the SIR scenario is that its dependence on the dynamical state of the conjugate disease is more complex than in the SIS case. Specifically, once an outbreak of one influencing infection is unfolding, the epidemic thresholds of the other disease may vary with time in nontrivial ways, depending on the nature and intensity of the different interaction mechanisms present. Figure 8 shows results from numerical simulations illustrating the previous phenomenology and the agreement with the analytical thresholds. Specifically, each panel represents the case in which the infection seed of the second disease is introduced in different moments for each topology, coinciding with different stages of the outbreak of disease 1: early phases [ Fig. 8(a), SF; Fig. 8(c), ER), the outbreak's peak [ Fig. 8(b), SF], and at the end of the outbreak [ Fig. 8(d), ER]. As we see, regardless of the topology, the parameters' values or the moment of appearance of the infection seed, our model adequately foresees the epidemic threshold and its variations with time. As in the SIS case, the influence of the interactions on the epidemic threshold is smaller in the case of scale-free networks, due to the smaller values of the primary thresholds in that case [in fact, primary thresholds have not been represented in Figs. 8(a) and 8(b) for the sake of clarity, because its values λ 1 2c ¼ 0.00316 virtually overlap the secondary thresholds].

III. CONDITIONS FOR DISEASE ENHANCEMENT AND IMPAIRMENT
Finally, it is interesting to analyze which combinations of parameters give raise to situations in which there is either enhancement or impairment of the diseases, as well as some other cases in which both effects are possible. This can be done for both the SIS and the SIR scenarios by studying the sign of the difference between the primary and the secondary thresholds.
For the SIS case, if we focus, for example, on the second disease, we have  ) and in an Erdos-Renyi (ER, right) network: I 1 ¼ P k;l Pðk; lÞðIS k;l þ II k;l þ IR k;l ÞðtÞ, R 1 ¼ P k;l Pðk; lÞðRS k;l þ RI k;l þ RR k;l ÞðtÞ. Disease 1 enhances the spreading of disease 2 in both cases, but still is not influenced by the second infection. Lower panels: Total number of individuals affected by disease 2 as a function of λ 2 in each network: R ∞ 2 ¼ lim t→∞ P k;l Pðk; lÞðSR k;l þ IR k;l þ RR k;l ÞðtÞ. Dynamical parameters: (a),(b) scale-free network: N ¼ 5000, hki ¼ 4.00,hli ¼ 5.11, As before, in both cases, disease 1 is not influenced by disease 2, and so β a where Cðk; lÞ is always positive. In this sense, λ c 2 ðσ 1 Þ − λ c 2 ð0Þ > 0 implies that disease 1 is reducing the epidemic threshold of disease 2, thus, enhancing its spreading. In the opposite case, λ c 2 ðσ 1 Þ − λ c 2 ð0Þ < 0, disease 1 causes the secondary threshold for disease 2 to be bigger than the primary one, hence, impairing its spreading. As we see, the condition yielding one or the other case involves a complex combination of the parameters. However, it is trivial to show that, provided that β a 2 > 1, β b 2 > 1, and η 2 < 1, disease 1 enhances disease 2 spreading for any value of σ 1 greater than zero. We call this scenario coherent enhancement, because all the interaction mechanisms contribute to enhance the spreading of disease 2. In a similar way, if β a 2 < 1, β b 2 < 1, and η 2 > 1, the interaction has the opposite sign regardless of σ 1 , and we have a situation of coherent impairment. Noticeably, if none of these conditions is fulfilled, there exists the possibility that the sign of the influence that disease 1 exerts on the spreading of disease 2 depends on its prevalence via σ 1 (see below). In Fig. 9, we represent λ c 2 ðσ 1 Þ − λ c 2 ð0Þ for different parameter combinations that cover all possible phenomenologies.
For the SIR case, the proliferation of dynamical classes and the possibility that other mechanisms (i.e., those including R individuals) carry the interaction between both diseases makes the equivalent expression more complex, but still derivable as λ c 2 ðhl 2 SSi; hl 2 ISi; hl 2 RSi; σ 1 Þ − λ c 2 ðhl 2 i; 0; 0; 0Þ where C 0 is always positive too. As noted previously, in the interacting SIR model, once an outbreak of one disease has started, the temporal dependence of the epidemic threshold of the second infection depends on the dynamic state of the system as well as on the parameter values that account for the mechanisms of interaction that are present in the system and their intensities. In such a case, if the interaction between the diseases is mediated by the class R (i.e., being recovered from one disease is what causes subjects' dynamics with respect to the other infection to be altered, and so β and η parameters are equal to 1), the density of individuals belonging to the R class is a monotonically increasing function of time. This trivially implies that, in general, the closer the epidemic process for one disease is to its end, the deeper its impact on the conjugate one. However, if it is the I class, the one that carries the interaction, it is easy to see that an interaction will take place if and only if the outbreak of disease 2 happens within a time window that, at least partially, overlaps with that characterizing the outbreak of disease 1. In such a situation, the interaction between the two diseases is a transitory property of the system, and the shift in the epidemic threshold crucially depends on the temporal co-occurrence of the outbreaks of both diseases.
A similar time-dependent shift of the epidemic thresholds can also take place in the SIS model, for example, for the thresholds of disease 2, if the difference λ c 2 ðσ 1 Þ − λ c 2 ð0Þ changes its sign during an outbreak of disease 1, as a consequence of the evolution of σ 1 . This is the scenario given by the green curves and the red curves in Fig. 9. However, it is worth noting that such combinations of parameters are somehow pathological, as they imply the appearance of diverse noncoherent mechanisms of FIG. 9. Relative variation of epidemic threshold of disease 2 as a function of σ 1 for different parameter sets. Black line: Coherent enhancement of disease 1 over disease 2: The networks are two Erdös-Renyi graphs: interaction (i.e., one disease simultaneously enhances subjects' susceptibility to disease but, at the same time, it reduces their infectiousness) that are epidemiologically less plausible. On the contrary, in the SIR framework, this transitoriness is an intrinsic property of the system when the infected class is responsible for influencing the dynamics of the conjugate disease. This transitory nature of the interactions constitutes an essential difference between the double SIR model with respect to the SIS one, which is at the root of the remarkable difference between their corresponding thresholds and the richer phenomenology of two interacting SIR-like processes.

IV. CONCLUSIONS
In this work, we propose a composed HMF model to describe the simultaneous spreading of two diseases over the same population but driven by independent mechanisms of transmission taking place on different networks of contacts. Within this framework, we thoroughly study extensions of the SIS and SIR models, where the parameters defining the infectious and recovery transitions of one disease depend on the state of each node with respect to the conjugate disease, in this way establishing a coupling between both diffusion processes. Our modeling approach presents different advantages with respect to previous models [36][37][38][39], as it simultaneously allows analytical derivations of the epidemic thresholds and an approximate description of the temporal evolution of the system, in addition to providing a way to isolate the effects on spreading dynamics of each possible interaction mechanism, such as variations of infectivity, susceptibility, or infectious periods. In addition, it enables us to solve the two paradigmatic modeling scenarios (SIS and SIR), identifying relevant differences between the two cases that arise as a consequence of disease interactions.
The model we present and analyze, even if it is based on a HMF and so constitutes a first approximation to the problem that neglects dynamical correlations among nodes, is able to identify novel phenomena qualitatively different from what is found on HMF descriptions of noninteracting systems. First, our model foresees different thresholds for SIR and SIS models, which arises from the asymmetry between the interaction mechanisms that take place under both models. Additionally, in what regards size scaling of epidemic thresholds, at least in the SIS case, we identify some relevant situations that yield threshold dependences on the network's size that are different from what is found in HMF models for single diseases. At this point, it is worth mentioning that asymmetries between SIS and SIR critical properties [20,51], or a different behavior with respect to that predicted by the HMF for epidemic thresholds [24], have also been identified in the context of single diseases spreading over one network, if the HMF approach is abandoned.
The situation here is qualitatively different, as all these divergent results with respect to classical HMF models of single diseases do not arise as a consequence of considering any dynamical correlation but as a consequence of disease-disease interactions. This has two relevant implications. On the one hand, deviations from HMF results on epidemic thresholds on single diseases identified as a consequence of considering dynamical correlations are of quantitatively residual relevance, precisely because, at the epidemic threshold, those correlations tend to vanish and, there, HMF models perform consistently well [24]. Instead, in our case, differences between SIS and SIR thresholds, for example, are completely different, as they may even depend on conceptually different interaction mechanisms. On the other hand, HMF is exact when dealing with annealed networks [20], and so, on these type of networks, which lack, by definition, dynamical correlations, both SIS-SIR asymmetries and eventually anomalous threshold dependences with size that we have identified in this work constitute phenomena not found before.
In addition, for the first time, the modeling framework we propose here allows one to isolate the independent effects of the different mechanisms of interaction that can determine the critical properties of the model. On the one hand, this allows us to foresee that, in a SIS model, certain interaction schemes may yield to the effects of one disease on another disease of different sign as a function of the prevalence levels of the former. On the other hand, for the SIR model, we also discuss how the interaction between the two diseases and the different dynamical classes give rise to a richer phenomenology. In particular, we show that, due to the transitory nature of the SIR spreading processes, the moment at which the interaction of the two diseases is made effective might greatly determine the values of the epidemic threshold of the disease whose course is modified by the other disease.
Further advances in the field should address the influences of dynamical correlations between nodes on the spreading of interacting epidemics. This constitutes a truly conceptual challenge as, for example, pairwise descriptions of the disease dynamics [51] should transform into quadruplet-based descriptions, as, in this case, dynamical correlations go beyond the dynamical state of neighbors in a network, but comprise both dynamical states on the two networks involved, which multiplies the complexity of the description.

ACKNOWLEDGMENTS
The authors thank Maitane Gutierrez for assistance with design software to produce some of the figures. J. S. was supported by MINECO through the FPU program. This work has been partially supported by the National Natural In the particular case in which both networks are regular graphs, an independent derivation of the epidemic threshold can be obtained by a linear stability analysis. In this case, for which Pðk; lÞ ¼ δðk − k o Þδðl − l o Þ, the dynamics can be described by the following system of four equations: one of which is linearly dependent of the rest. Thus, we analyze only the system constituted by the last three equations and use s ¼ 1 − IS − SI − II. In order to perform a linear stability analysis, we first linearize the system around the equilibrium point and we calculate the Jacobian to get which, taking advantage of the fact that the Jacobian itself is just the product of the elements in the diagonal, yields the stability conditions detailed in Table IV. Looking carefully at the results shown in Table IV, we first recognize the appearance of a couple of critical values for the infectiousnesses λ c 1i and λ c 2i which are referred to in the following as the primary thresholds of their respective diseases, just as in the general case described in the main text. These threshold values stand for the minimum values of the infectiousnesses that yield epidemic outbreaks after the introduction of infinitesimal seeds of infected individuals of each disease in an initially healthy population. Therefore, the condition λ 1 > λ c 1i must be verified in order to have an epidemic outbreak for the first disease once an infinitesimal seed ðIS; SI; IIÞ ¼ ðϵ; 0; 0Þ has been introduced on a system being at the disease-free fixed point ðIS; SI; IIÞ ¼ ð0; 0; 0Þ. On the other hand, the condition λ c 2i plays an equivalent role for the second disease with respect to a seed ðIS; SI; IIÞ ¼ ð0; ϵ; 0Þ.
The fixed point ðIS; SI; IIÞ ¼ ð0; 0; 0Þ is, strictly speaking, not the only possible disease-free fixed point in our model. Another two partially disease-free fixed points can exist: a first fixed point in which disease 1 is installed in the system at a certain prevalence π 1 while disease 2 is absent ðIS;SI;IIÞ ¼ ðπ 1 ;0;0Þ and its cognate ðIS; SI; IIÞ ¼ ð0; π 2 ; 0Þ, for which there is no individual infected with disease 1. As we show, the stability of these fixed points depends on the prevalence fractions π 1 and π 2 , which are the stationary proportions of sick individuals of each disease: Considering Eqs. (A2) and (A3), we study the stability of the first fixed point ðIS; SI; IIÞ ¼ ðπ 1 0Þ as a function of λ 1 : the Jacobian, around this fixed point, takes the following form: By solving the equation derived from imposing that J ¼ 0, we obtain the values of the parameters leading to stability inversion. As can be seen by the naked eye, ðμ 1 − λ 1 k o Þ is the eigenvalue ξ 1 associated with the eigenvectorψ 1 ¼ ð1; 0; 0Þ, and the condition ξ 1 ¼ 0 again yields the same condition λ 1 ¼ μ 1 =k o , thus, defining the threshold of the classical SIS model. Regarding the other two eigenvalues, ξ 2 and ξ 3 , the result is more cumbersome. Despite that, the vanishing of the 2 × 2 determinant of the right-hand inferior corner of the Jacobian matrix [Eq. (A4)] yields to which the general expression presented in the main text for the epidemic threshold reduces when Pðk; lÞ ¼ δðk − k o Þδðl − l o Þ. The agreement between numerical simulations and the analytic expression of the threshold presented here is accurate only when λ 1 k o ≃ μ 1 , and the reason is easily understood. We compare the reaction of the system to the introduction of a small seed of SI individuals when both diseases are absent [ðIS; SI; IIÞ ¼ ð0; 0; 0Þ, case 1] or when the first disease was already installed in the system [ðIS; SI; IIÞ ¼ ðπ 1 ; 0; 0Þ, case 2]. As we argued in the preceding sections, the epidemic threshold is different in each case, and the reason is simply the presence, in the second case, of a fraction π 1 of IS of individuals for which the infectiousness and recovery rates for disease 2 are different with respect to the rest of the individuals. As a consequence, the mean values of the dynamical parameters averaged over the whole population hλ 2 i and hμ 2 i are different in both cases and will yield different values for the epidemic thresholds. Thus, in order to accurately evaluate the secondary threshold for disease 2, it is essential to know with enough precision the prevalence π 1 corresponding to a single SIS model, as a function of λ 1 . The problem arises from the fact that, precisely, the derivation of Eq. (A5) explicitly assumes that the bijection between π 1 and λ 1 is governed by the mean-field stationary expression [Eq. (A2)], which is precise only when λ 1 k o ≃ μ 1 [25]. In fact, a way to rebuild a more accurate secondary threshold curve can be achieved if, from Eq. (A2), we substitute λ 1 as a function of π 1 , and then introduce the so-obtained expression into Eq. (A5) to get which allows one to directly evaluate the threshold as a function of π 1 rather than of λ 1 . Thus, by introducing in Eq. (A6) the π 1 values obtained from the simulations instead of the theoretical prediction of the mean field [Eq. (A2)],, we recover the curves for the secondary threshold represented with red lines in Figs. 2 and 3, in quantitative agreement with results from simulations. Obviously, the same arguments stand for the secondary threshold of the first disease, which obeys the following expression: When the networks are heterogeneous, this reformulation of the threshold curves as a function of π 1 or π 2 cannot be done straightforwardly, as no analytical bijection λ 1 ðπ 1 Þ or λ 2 ðπ 2 Þ can be reached. The best we can do is to use a numerically builtup relationship λ 1 ðθ 1 Þ [or λ 2 ðθ 2 Þ] and introduce it into the general expression for the threshold. Although the accuracy of the curves for the secondary threshold is very satisfactory in Fig. 4, we identify this effect to be the source of the slight divergence between the analytical and numerical secondary thresholds shown in Fig. 6.

APPENDIX B: VANISHING CONDITIONS FOR EPIDEMIC THRESHOLDS
Here, we provide a systematic analysis of the behavior of the secondary threshold in our model for large scale-free networks. In particular, we inspect whether the coupling between the spreading of the two diseases in the terms described in our model can modify the classical, singledisease scheme, and, if so, under what conditions. The epidemic threshold for the first disease in our model reads as and we have to address its behavior in the limit N → ∞. Substituting sums by integrals in Eq. (B1), one gets To study the behavior of the threshold in the thermodynamic limit, we present an analysis that is essentially based on the following result, whose proof is an exercise of elementary algebra that we present for the sake of completeness in the Appendix C. Given two polynomials PðkÞ, QðkÞ, we have that lim Focusing on scale-free connectivity distributions, we distinguish two different scenarios in this section: uncorrelated layers and totally correlated layers. In the first case, both networks present a scale-free distribution in which the connectivity of a node in a layer is essentially independent of its degree on the other layer, in such a way that the composed connectivity distribution verifies Pðk; lÞ ¼ C o k −γ l −Γ . Instead, in the second case, although both layers are also scale-free networks, the degree of any given node in both layers is forced to be the same. Therefore, nodes that are hubs in a layer are also hubs in the other layer, and the composed degree distribution is Pðk; lÞ ¼ C o δðk − lÞk −γ . In addition, we assume that both γ and Γ exponents are rational; hence, we can write γ ¼ w=x and Γ ¼ y=z, with ðw; x; y; zÞ ∈ N. By addressing these two opposite scenarios, our aim is to characterize the difference, in terms of the spreading dynamics, between the coupling of two diseases that spread by independent means-that is, through different networks of contactsand the coupling of two related diseases-or variations of the same disease-that spread following the very same mechanisms and, thus, they do so over highly correlated networks of contacts.

Uncorrelated scale-free layers
If we substitute Pðk; lÞ ¼ C o k −γ l −Γ into Eq. (B2), we can factorize the double integrals into independent terms: The condition for the numerator to diverge is γ ≤2∨Γ≤1. These conditions are verified by networks whose mean connectivity diverges in the thermodynamic limit. For this reason, networks with those small exponents are neither reasonable systems to be used to describe any information diffusion process over them nor are they found in real systems. Taking this into account, we restrict our analysis to the more epidemiologically relevant scenario in which γ > 2 and Γ > 2, although our reasoning is generalizable to any value of the exponents. Thus, in the scenario γ > 2 and Γ > 2, the numerator always remains finite and the only phenomenon of interest that could be found is an eventual vanishing of the threshold due to a divergence in the denominator-R Dðk; lÞdkdl in what follows. Dðk; lÞ can be factorized to give Dðk; lÞ ¼ D 1 ðkÞD 2 ðlÞ, which allows us to factorize the integral, as can be seen from Eq. (B4). The integral in k, R D 1 ðkÞdk will diverge if and only if γ ≤ 3, but the integral in l, R D 2 ðlÞdl, might independently diverge under some conditions. That situation would suppose the vanishing of the threshold of the first disease as a consequence of its coupling to the second disease rather than to internal, dynamical, or topological features. To find out the conditions for R D 2 ðlÞdl to diverge, we make the following change of variable: so as to change the dependence of the denominator Dðk; lÞ ¼ D 1 ðkÞD 2 ðlÞ in Eq. (B4) to Dðk; mÞ ¼ D 1 ðkÞD 2 ðmÞ: min fm 2z σ 2 2 β a 2 β a 1 β b 1 þm z σ 2 ½η 2 μ 2 β a 1 þβ b 1 ðβ a 1 μ 1 þβ a 2 μ 2 Þþμ 2 ðη 1 μ 1 þη 2 μ 2 Þgzm z−1 m y ½m 2z σ 2 2 β a 2 η 1 þm z σ 2 ðη 1 μ 1 þη 2 μ 2 þβ a 2 η 1 μ 2 Þþμ 2 ðη 1 μ 1 þη 2 μ 2 Þ dm : ðB6Þ The last expression has the advantage that the argument of the integral in m, D 2 ðmÞ is just the quotient of two polynomials, and thus, its behavior in the limit l max → ∞ is governed by Eq. (B3) and depends only on the difference of degrees of the denominator and the numerator. If min fm 2z σ 2 2 β a 2 β a 1 β b 1 þ m z σ 2 ½η 2 μ 2 β a 1 þ β b 1 ðβ a 1 μ 1 þ β a 2 μ 2 Þ þ μ 2 ðη 1 μ 1 þ η 2 μ 2 Þgzm z−1 m y ½m 2z σ 2 2 β a 2 η 1 þ m z σ 2 ðη 1 μ 1 þ η 2 μ 2 þ β a 2 η 1 μ 2 Þ þ μ 2 ðη 1 μ 1 þ η 2 μ 2 Þ dm ¼ Therefore, in the region of interest γ > 2 ∧ Γ > 2, the integral in m: R D 2 ðmÞdm does not diverge, and so, the threshold does not vanish. It is worth noticing that the condition Γ < 1 does not guarantee the vanishing of the threshold, as it also causes the numerator to diverge. Expression (B8) is valid for the case in which none of the β and η parameters of the model vanish. If this is not the case, i.e., if some of the infectiousness variations β do vanishthe rest of the parameters cannot do that within a realistic epidemiological framework, the degrees of the numerator and the denominator in Eq. (B7) could vary. In Table V, we systematically address all possible combinations of null parameters and the composed conditions yielding a vanishing threshold for each case.
Once again, we see that, whatever the interacting scheme between both diseases is, the denominator always remains finite provided that Γ > 2. In conclusion, if no degree correlation is introduced between layers, and for realistic systems characterized by double power laws verifying γ > 2 ∧ Γ > 2, the behavior of the thresholds in the thermodynamic limit is essentially the same as the behavior of the uncoupled systems: the threshold associated with the first disease vanishes if and only if the exponent of its own network verifies γ ≤ 3, whatever the exponent of the second network. The coupling will introduce, in general, only a finite prefactor. The symmetric situation obviously stands for the threshold of the second disease.