The role of time scale in the spreading of asymmetrically interacting diseases

Diseases and other contagion phenomena in nature and society can interact asymmetrically, such that one can benefit from the other, which in turn impairs the first, in analogy with predator-prey systems. Here, we consider two models for interacting disease-like dynamics with asymmetric interactions and different associated time scales. Using rate equations for homogeneously mixed populations, we show that the stationary prevalences and phase diagrams of each model behave differently with respect to variations of the relative time scales. We also characterize in detail the regime where transient oscillations are observed, a pattern that is inherent to asymmetrical interactions but often ignored in the literature. Our results contribute to a better understanding of disease dynamics in particular, and interacting processes in general, and could provide interesting insights for real-world applications, most notably, the interplay between the dynamics of fact-checked and fake news.


I. INTRODUCTION
Spreading processes are ubiquitous in nature and society. A primary example of these phenomena is the propagation of diseases in human and animal populations [1,2]. Despite many recent advances in the theoretical, computational and data-driven modeling of diseases, there are still many scientific problems that remain open. Two of such problems are currently of high relevance. On the one hand, we have many limitations when it comes to understand the effects of mobility restrictions and human behavioral changes on the evolution of a disease [3]. Secondly, there are many diseases that rarely evolve in isolation, on the contrary, diverse viral strains or different pathogens can either compete for the susceptible population or establish cooperative interactions, both at a population level or inside a host's organism. This paper focuses on studying spreading processes with the aim of shedding some light into the second kind of challenge.
Early works on interacting diseases [4][5][6][7][8] have modeled the dynamics of pathogens or strains that interact competitively, showing that their coexistence is not always possible. In the last years, these models have been adapted to networks [9][10][11][12][13][14][15], including single layer and multiplex networks [2,16], as well as metapopulations [13]. These works have showed that the network organization plays a fundamental role on the evolution of the dynamical processes, affecting the epidemic threshold and prevalence (e.g. [17]). Recently, the focus has also been placed on collaborative contagion, in which there is a positive feedback between diseases [18][19][20][21][22][23][24]. In these models, discontinuous phase transitions, in which the prevalence goes from zero to a finite value abruptly, have been observed, as well as simultaneous stability of two epidemic states. It is worth remarking that although the network topology can generate important phenomena, some essential properties can still be observed in homogeneously mixed populations, as demonstrated by Chen and collaborators [24].
General models for interacting diseases have also been developed recently. This is the case of [25], where the authors introduced a generic model for two interacting dis-eases in multiplex networks, comprising the cases of competitive, collaborative and asymmetrical interactions, for both SIS (Susceptible-Infected-Susceptible) and SIR (Susceptible-Infected-Removed) compartmental models. They calculated the respective epidemic thresholds and studied the prevalence for both the competitive and the collaborative scenarios. However, the asymmetrical case has remained unexplored.
Asymmetrical interactions occur in many complex systems. For example, there are reported cases of HIV viral load suppression by the presence of some other pathogens [26][27][28], which suggests an asymmetrical interaction between HIV and many diseases. Asymmetrical interactions have been studied mostly in the case of the interplay between epidemic and awareness [17,[29][30][31][32][33][34], malicious computer worms and spreading countermeasures [35][36][37][38], and antigens and immune system agents [39][40][41], among other examples. In the case of disease and awareness, the epidemic stimulates information awareness, which in turn tends to reduce the exposure to disease and therefore also in disease prevalence, configuring an asymmetrical interaction between the two processes. It has been shown that, in general, the awareness can effectively reduce the disease prevalence and increase its threshold [29,30]. The epidemic, in turn, can sustain an information outbreak even when the latter is bellow its "independent" threshold [17,31]. Correlations and structural properties of the underlying multiplex network can also influence this propagation process [17,29,31].
As mentioned above, some previous works have examined asymmetrically interacting spreading phenomena. However, this has been done for specific cases, and a general description of these processes is, to the best of our knowledge, still missing. Of special interest is the fact that the interacting dynamics can do so at different time-scales, whose effect has not been deeply investigated yet. Changes in the relative timescale (i.e., the possibility that one process has a different intrinsic clock concerning the other) typically do not change the qualitative behavior of interacting systems, but often lead to relevant quantitative effects. For example, Karrer and Newman [10] showed that, for competing SIR diseases, the slower pathogen may take advantage, and even win the competition due to its inherent time scale. Oliveira and Dickman [43] also showed that the slower species can win a competition in a contact process (CP), including scenarios with periodic and stochastic environmental variations. In the context of asymmetrical interactions, Wu and collaborators [41] and Poletto et al. [42] demonstrated that the epidemic threshold is affected by changes in the recovery rate when the reproductive number (the ratio between the spreading and the recovery rates) is kept constant, which is equivalent to a variation of the time scale. Equivalently, it has also been shown that a slower information awareness may have more impact on the epidemic prevalence than a faster one [33,34]. Therefore, the development of a more general understanding about the influence of the relative time scales of asymmetrically interacting disease-like models is of great interest.
In this paper, we study two SIS-like models for interacting diseases when they are in an asymmetrically interacting regime over homogeneously mixed populations and continuous-time evolution. We rely on the most simple setup that doesn't involve structured populations to extract the intrinsic properties of the models. In each model, we assign a parameter π that allows us to control the relative time scale between the two dynamics while keeping constant the epidemic forces of each disease. We adopt a descriptive analysis of each model, determining whether different aspects are or aren't influenced by the time scale, including the respective phase diagrams, coexistence and oscillatory behavior. We also compare both models and discuss our results in the light of some previous works.

II. EPIDEMIC MODELS WITH ASYMMETRICAL INTERACTION
In what follows, regardless of the model considered, we assume that there are two diseases that interact. Moreover, we consider that disease I is the "prey" (it is impaired by the other disease) and disease II is the "predator" (it is benefited from the presence of the first disease), in an analogy with the asymmetrical interaction of predator-prey systems. A parameter, π, controls the relative time scale between the two diseases. It is implemented as follows: the rates of all processes promoted by disease I are multiplied by 1 − π, whereas the rates of disease II processes are multiplied by π. Making π to range between 0 and 1, we sweep through scenarios in which disease I is faster (π < 0.5) or slower (π > 0.5) than disease II, as well as the balanced case (π = 0.5). Although π does not add a new degree of freedom of the model, the advantages of using this approach are: (i) it allows us to control the relative time scale with a single parameter; (ii) the "overall rate" of the system, which can be regarded as ∼ (1 − π) + π, is kept constant when varying π; and (iii) plotting variables as a function of π is simple because it is limited between 0 and 1. Let us now describe each of the models scrutinized in the rest of the paper.

A. Model A: interacting diseases through susceptibility change
In this variant, we consider that the presence of one disease impairs the spreading of the second one by changing the individual's susceptibility to catch the other disease. It has been used to describe the dynamics of competing pathogens with partial cross immunity [8] and for collaborative contagion [24,25], but its asymmetrical version is still largely unexplored. The latter regime is a prototypical model to describe processes in which an epidemic coexists with an information dynamics in which awareness regarding the disease plays a role in the spreading of it. Indeed, our Model A is mathematically similar to models used for that purpose [34].
In this model, each individual can either be susceptible to both diseases (S 1 S 2 ), infected by one disease and susceptible to the second one (I 1 S 2 ) or (S 1 I 2 ), or infected by both diseases (I 1 I 2 ). Note that we denote as I 1 and I 2 individuals that are, respectively, infected by disease I and II, regardless of their state with respect to the other disease. For completely susceptible individuals S 1 S 2 , the baseline contagion rate of disease I (II) when in contact with an individual infected by disease I (II) is β 1 (β 2 ). For individuals already infected by disease I but susceptible to disease II (I 1 S 2 ), the contagion rate for disease II is Γ 2 · β 2 , i.e., it is multiplied by a factor Γ 2 . The same holds for S 1 I 2 individuals, for which the contagion rate by disease I is changed to Γ 1 · β 1 . The healing rates from disease I and II are respectively µ 1 and µ 2 , and are not affected by the other disease. Figure 1 represents all the possible transitions for this model, with their respective time scale factors as explained before.
The asymmetrical interaction between the two diseases can be obtained by setting 0 ≤ Γ 1 < 1 and Γ 2 > 1. This means that individuals that hold disease II are less susceptible to catch disease I in comparison to fully susceptible individuals. On the other hand, individuals infected by disease I are more likely to catch disease II. Therefore, disease I enhances the propagation of disease II, whereas disease II impairs the propagation of disease I. We represent the density of individuals in a given state X by ρ x , with X being either a composite state (like I 1 S 2 ) or a simple state (like I 2 ). Based on the diagram from figure 1, the time evolution of the composite state densities is given by the following equations: The densities are yet subject to the normalization constraint State transitions allowed in model A. The baseline infection and healing rates of disease I (II) are respectively β1 (β2) and µ1 (µ2). Γ1 (Γ2) represents the modification to the baseline transmission rate of of disease I (II) due to the presence of the other disease in the host. Besides, each rate is multiplied by its corresponding time scale factor: 1 − π for processes of disease I and π for disease II.
This constraint makes the system effectively three-dimensional. One can reduce the number of equations and simplify the notation by using the following variable change (as done in [24]): For which the dynamical equations are:

B. Model B: competing diseases with superinfection
Model B constitutes a modification of models of competing strains, in which a host cannot have the two diseases at the same time. This could be achieved from model A by setting Γ 1 = Γ 2 = 0. However, we also allow the in-host disease replacement via superinfection: if an individual infected by disease I contacts another one infected by disease II, the first S I 1 FIG. 2. State transitions of model B. Infection and healing rates of disease I (II) are respectively β1 (β2) and µ1 (µ2). I1 individuals can be (super)infected by disease II with rate αβ2. Each rate is multiplied by its corresponding time scale factor: 1 − π for disease I processes and π for disease II processes.
can also become infected by disease II, which immediately replaces disease I in the host. The other way around, from disease II to I, is not possible. Superinfection is a phenomenon that is claimed to occur for some diseases such as HIV [44][45][46] and bacterial pathogens [47], although it does not necessarily leads to in-host replacement of the first infection. There is a considerable amount of works about epidemic models with superinfection, both with homogeneous populations [48][49][50] and complex networks [39][40][41] and including other realistic aspects such as demography. Beyond epidemiological examples, the present model could abstract the interaction between computer viruses and the spreading of anti-malware, or the dynamics of fake news and fact-checking messages. In such scenarios, either the anti-malware or the fact-checked message replaces the previous "infectious agent". We also note that ours is a particular case of the in [41], and can also be interpreted as a generalized predator-prey model [51]. Specifically, in model B, we represent susceptible individuals simply by S, and infected individuals of diseases I and II, respectively, by I 1 and I 2 . The transmission rates are β 1 and β 2 , the healing rates are µ 1 and µ 2 , and the rate at which I 1 individuals are "superinfected" by disease II when exposed to I 2 individuals is given by a modified term α · β 2 . As in model A, each term is also multiplied by the corresponding time scale factor ((1 − π) for disease I and π for disease II). The transitions are schematically represented in figure 2. Despite the competitive aspect of the model, one can generate an asymmetrical interaction by setting α > 1. This is because I 1 individuals, in this case, are more easily infected by disease II than susceptible ones, meaning that the spreading of disease II is enhanced by the presence of disease I, which in turn is in disadvantage due to the competition with disease II.
Using the same notation as before, ρ x , for the density of individuals in state X, we write the dynamical equations: where the normalization constraint is ρ s +ρ i1 +ρ i2 = 1. As in model A, we make the change of variables u = ρ i1 , v = ρ i2 to obtain the reduced set of dynamical equations: III. RESULTS

A. Phase diagrams
In the asymmetrically interacting regime, models A and B share a common feature: for any given set of parameters, there is exactly one stable fixed point within the region of the phase portrait that corresponds to physically possible solutions. Therefore, unlike what has been reported for mutually competitive or cooperative scenarios [24,48], our models do not present bistability [52]. This in intrinsic to the positivenegative feedback of the asymmetric interaction between the diseases. Both models A and B present four phases, separated by transcritical bifurcations: (I) no disease, (II) disease I only, (III) disease II only and (IV) coexistence of both diseases. Here, we investigate the λ 1 ×λ 2 phase diagrams, setting the other parameters to fixed values. Figure 3 shows the phase diagrams for both models.
Through a stability analysis, we can obtain the phase transition curves of both models. The boundary between phases (II) and (IV) in model A is expressed as: whereas, by the symmetry of the model, the boundary between phases (III) and (IV) is given by: The other two bifurcations are trivial and are given by λ 2 = 1 for 0 ≤ λ 1 ≤ 1 (boundary between (I) and (III)) and λ 1 = 1 for 0 ≤ λ 2 ≤ 1 ((I) and (II)). It is worth noticing that, for model A, none of the phase transition curves depends on the time scale parameter π.
For model B, the boundary between (II) and (IV) is given by: which is similar to that in equation 6, only replacing Γ 2 by α.
The boundary between (III) and (IV) is expressed as: where we define χ (interpreted as the time scale ratio between diseases II and I) as: The trivial boundaries of region (I) with regions (II) and (III) are the same as in model A. Notice, however, that the parameter χ depends increasingly on the time scale parameter π, and so does the critical λ 1 value from equation 9. This means that, as π increases (i.e., when disease II propagates on shorter time scales), the phase transition between (III) and (IV) moves to lower values of λ 2 , contracting the region of coexistence. Therefore, in this model, a faster clock for disease II makes it more effective to suppress disease I. This is an interesting result that might be used to control the prevalence of disease (or any other "infectious agent") in the host population. The π-dependence of the phase diagrams of model B constitutes an important difference with respect to model A, and was already reported by Wu and collaborators [41] as a dependency on the recovery rate when keeping the ratios λ/µ constant.

B. Behavior of the stationary prevalence
In the phase of coexistence (region (IV) of the phase diagrams), both models present a single stable fixed point, for which u, v = 0. One can use the dynamical equations (eqs. 1 to 3 for model A and 4 and 5 for model B) to derive analytical expressions for the prevalences at the coexistence fixed point. The derivation and the final expressions are shown in appendix A. As expected, the stationary prevalence of each disease is a non-decreasing function of its reproduction ratio λ = β/µ, when considering other parameters as fixed. However, the dependence of the prevalences with the relative time scale parameter π is not trivial, and is different in each model. Figure 4 shows the basic behavior of the fixed point prevalences with π for models A and B. While the prevalence v of disease II decreases with π for both models, the prevalence u of disease I has opposite behaviors in each of them. In model A, the prevalence w of coinfection increases with π. This variant shows the same behavior (if equating disease II with the information) as recently reported in some works on epidemic with awareness [33,34] in complex networks, namely, a faster relative clock of the information induces an increase on the disease prevalence. For model B, however, the behavior is the opposite: the prevalence of disease I decreases with π. Considering also how the time scale parameter distorts the phase diagram of model B (see figure 3), we see that a faster clock of the "predator" process (disease II) effectively decreases (and possibly leads to the extinction of) the spreading of the "prey" process. Therefore, the relationship between the prevalences and the relative time scale in asymmetrically interacting spreading phenomena is a feature that depends on the specific shape of the considered model. In figures 5 and 6, we further analyze the behavior of the prevalences with π and the parameters that control the interactions: Γ 1 , Γ 2 for model A and α for model B. We first notice that the increasing or decreasing trends with respect to π, as observed in figure 4, are not changed for different values of the interaction parameters: disease I increases while disease II decreases with π for model A (figure 5), whereas both prevalences decrease with π for model B ( figure 6). To rule out the possibility that there is a region of the parameter space in which the reported behaviors with respect to π might be different, we show in appendix B that this is not possible, i.e., that the behavior with π is always the same for each model in the coexistence region.
We can also analyze how the prevalences vary with changes in the interaction parameters. For model A (figure 5), we see that an increase in Γ 1 increases both the prevalences of disease I (a) and II (b). This is expected, as a greater Γ 1 value means a weaker impairing to the propagation of disease I, which is beneficial for both diseases (as disease II benefits from disease I). On the other hand, an increase in Γ 2 causes a decrease in disease I (c) and an increase in disease II (d). This is also expected, as a larger value of Γ 2 means a greater benefit to disease II, which in turn is detrimental to disease I. Therefore, for model A, the effect of the two interaction parameters in each disease's prevalence is intuitive and predictable. The effects are also numerically influenced by the time scale π, yet not qualitatively changed. However, for model B, which has a single interaction parameter α, the behavior of the prevalence is not trivial. From figure 6, we see that, while disease I prevalence (a) is always reduced with an increase in α, the behavior of the prevalence of disease II (b) with α is not uniform, and may have an optimal value that depends on π. This happens because the superinfection transition, controlled by α, is simultaneously beneficial to disease II and detrimental to disease I. Thus, as seen from model A, an increase in α certainly reduces the prevalence of disease I, but has a "conflicting" effect to the prevalence of disease II. While a value of α close to 1 means almost no benefit to disease II from disease I, a large value α 1 means an excessive "predation" from disease II, therefore existing an optimal intensity of the interaction α. This is further illustrated in Figure 7,where we show the prevalence of disease II in model B as a function of α, for different values of π. As it can be seen, there is an optimal value of α that maximizes the prevalence, for fixed values of the other parameters.
C. Damped oscillations: node vs spiral point Some predator-prey systems, which naturally have an asymmetrical relationship between two processes, are known to present stable closed orbits (i.e., sustained oscillations) [51]. For both epidemic models addressed in our work, there are no closed orbits in the physical region of the phase portrait [53]. However, in the coexistence phase (region (IV)), the stable fixed point can be either a node or a spiral point. In the second case, the transient dynamics of the system towards the fixed point may present some damped oscillations. The presence of such local oscillations is determined by the imaginary part of the eigenvalues of the model's Jacobian matrix, calculated at the stable fixed point. For model B, which is two-dimensional, the Jacobian's eigenvalues σ 1 , σ 2 can either be both real or complex conjugate to each other. For model A, which is three-dimensional, the 3x3 Jacobian matrix can either have none or two non-real conjugate eigenvalues. For both models A and B, one can use the quantity: to measure the "quality factor" of the oscillations when Re(σ i ) = 0 ({σ i } are the eigenvalues of the Jacobian at the fixed point, d is the dimensionality of the system). This is because the imaginary part is responsible of the oscillations, and the real part of the damping, thus the imaginary-to-real part ratio measures the propensity of the system to oscillate around the fixed point. In Figure 3, together with the regular phases of the model, we show in (color-coded) green the numerically calculated values of Q. The coexistence phase can thus be subdivided according to the existence of a non-real Jacobian eigenvalue: within the white regions, the eigenvalues are real and the fixed point is a node whereas in the green areas, the fixed point is a spiral point and there may occur oscillations around it. Notice that, for both models A and B, the shape of the spiral point region depends considerably on the time scale parameter π. Greater values of π seem to shrink the oscillatory region and reduce the values of Q for model A, but the opposite appears to happen with model B. Furthermore, we show the difference between a node and a spiral point in figure 8, in which the time evolution of the prevalences of model B is shown for two situations: one with Q = 0 (hence no local oscillations (left)), and another with Q = 3.65 (thus there are damped oscillations before reaching the steady state). Interestingly, we only had to change the relative time scale parameter π to switch between the two situations.
An important fact is that damped oscillations can only occur when the interaction between the diseases is asymmetrical, for both models A and B. For model A, we numerically check this by observing that the oscillatory region of the phase diagram shrinks and disappears as Γ 1 or Γ 2 leaves the region for which interactions are asymmetric. For model B, it can be shown that the Jacobian's eigenvalue equation (which is quadratic) can only assume non-real solutions if α > 1. Finally, we note that despite that our deterministic formulation predicts that the oscillations are always damped, stochasticity -which is intrinsic to real-world systems, and are inherent when performing Monte Carlo simulations -could cause such oscillations to last for the long term, as small perturbations to the prevalences could recover the oscillatory pattern.

IV. CONCLUSIONS
In this work, we have studied two minimalist models for interacting diseases in the asymmetrical regime, for homogeneously mixed populations and with continuous-time evolution. We focus on the influence of the relative time scale between the two diseases. The simplicity of our framework not only provides us with analytical tractability, but also reveals the fundamental properties of asymmetrically interacting contagion. The models are simple enough to be applied to different situations, and the choice for two models is justified by our goal of achieving a more general understanding about interacting processes.
Model A is a mathematical prototype for epidemic with awareness and the possible asymmetrical interaction between HIV and some specific diseases [26][27][28]. Our results show that, despite the epidemic thresholds are not affected by the time scale parameter π, the prevalence is influenced in a nonintuitive way: the "prey" disease (I) has greater prevalence when the "predator" disease (II) has a faster time scale. If we take the HIV as an example, which naturally has a long time scale for its development, detection and treatment, this means that the information awareness (which is more quickly transmitted and forgotten) may not be as efficient as it could due to its shorter time scale.
Model B, which is inspired in situations such as computer viruses and spreading countermeasures [35][36][37][38] or fake vs fact-checked news, has a different behavior with the relative time scale: if the "predator" disease (II) has a faster clock, both prevalences of diseases I and II decay, with the possibility of disease I eradication (see figure 4.b). The interpretation then varies according to the situation being modeled. For example, if the main "goal" of disease II is to fight the other disease (such as a fact-checked against fake news), then a shorter time scale (higher π) is desirable. If however the goal is to maximize disease II prevalence by taking advantage of disease I (one can imagine a computer virus that makes use of vulnerabilities generated by another one), then a longer time scale for disease II is preferable.
Another interesting phenomenology revealed is the behavior of model B with its superinfection parameter α, which controls the interaction between diseases I and II in both directions. The interpretation again depends on the situation: if the goal is to minimize disease I, then a greater value of α is always better, whereas for the optimization of the prevalence of disease II an intermediate value of α (which depends on π) should be sought. Finally, we have shown that oscillatory behavior, which is a well-known feature of predator-prey systems, is also present in asymmetrically interacting epidemic models, and that its expression crucially depends on the relative time scale. Although the theory predicts that they are strongly damped, real world systems could persistently display such oscillations due to random fluctuations.
To round off, we note that apart from the above connections to real dynamical systems, we expect that our work helps future studies of asymmetrically interacting spreading phenomena by providing general guidelines. By studying two different models, we provide a more general understanding of such systems, but we recognize that the results may have some bias due to the specific choice of models. An interesting extension of our work would be to consider a broader family of models, using more generic formulations, and determine for example the conditions that make the prevalences to depend, in one way or another, on the time scale. Finally, we also expect some variations when these models are studied on top of structured populations. For such scenarios, our results provide a baseline of which qualitative and quantitative patterns are not to be associated to a network effect.  The coexistence phase is characterized by the long-term presence of both diseases, thus having u, v > 0 in the steady state. In this appendix, we show how to analytically calculate the prevalences for the coexistence fixed points of both models. For simplicity of notation, throughout the whole appendix section we use u, v and w to denote the fixed point values of the prevalences (and not their time-dependent functions).

Model A
We find the fixed points by settingu =v =ẇ = 0 in the dynamical equations 1 to 3. Using the fact that u, v > 0 at the coexistence, we divide equation 1 by (1 − π)µ 1 u, equation 2 by πµ 2 v and equation 3 by (1 − π)µ 1 , obtaining the following reduced system of equations: in which only equation A3 is non-linear. The disease II time scale factor χ is the same as defined in Eq. 10. Using also the definition: where s i is the solution for a non-interacting SIS model with reproduction ratio λ i , we can further simplify equations A1 and A2 respectively to: from which it is intuitive to see that the interaction between the diseases (represented by Γ 1 and Γ 2 ) causes a deviation from the non-interacting solution.
At this point, we split the solution in two cases: (a) the simpler case Γ 1 = 0 (i.e., when disease II completely inhibits infection by disease I) and (b) the more general case Γ 1 > 0.
a. Case Γ1 = 0 In this case, equation A5 simplifies to: Combining this equation with A6, we get a 2x2 system for v and (u − w), for which the solution is: Now using that Γ 1 = 0 in equation A3, we obtain an expression for w in terms of known variables: Plugging equations A8 and A9 in the above equation, and then using it back to equation A9, we get the full expressions for the fixed points when Γ 1 = 0: where we have also replaced λ 2 = 1/(1 − s 2 ). With some analysis of the above expressions, noticing that χ/(1 − χ) is an increasing function of χ > 0, we can infer that u and w increase with χ (and thus with π), whereas v does not depend on the time scale factor. This is in agreement with the plots in figure 5 (a) and (b).
b. Case Γ1 > 0 In this case, the fact that equation A3 is quadratic on its variables cannot be avoided. Our strategy is to use equations A1 and A2 to write u, v, (u − w) and (v − w) as functions of the variable w and the model parameters. This can be used in equation A3 to find the solution for w and, consequently, for the other variables.
We can isolate v in equation A6 and apply it to equation A5, obtaining: where we simplified the notation using the definitions We can also manipulate equations A14 and A15 to obtain: Now we plug the above expressions into equation A3 which, after redistributing the terms and dividing them by Γ 1 Γ 2 m 2 (m > 0 for asymmetrical interactions and, by hypothesis, Γ 1 > 0), becomes the quadratic equation aw 2 + bw + c = 0, where: For asymmetrical interactions, it is possible that a = 0 in the coexistence region. Thus we write the coexistence fixed point of the system as: From the above expressions, it is difficult to extract the behavior of the prevalences with respect to the time scale parameter χ (or π). However, in appendix B we present an alternative argument that reinforces the behavior observed in Fig.  5.

Model B
Model B has a much simpler procedure to find analytical expressions for the coexistence fixed point u and v, for general values of the parameters. Knowing that u, v > 0, one can divide eq. 4 by (1 − π)µ 1 u and eq. 5 by πµ 2 v and set their left-hand sides to 0, obtaining: which is a 2x2 linear system in u and v. Passing convenient terms to the left side of each equation, one can write the system in terms of s 1 and s 2 as defined in A4: where we define φ as: The solution of this 2x2 system is: Appendix B: Behavior of the prevalences with π By plotting the values of the fixed point prevalences, analytically derived in appendix A, we could analyze the behavior of such prevalences with the time scale parameter π for models A and B. In this section, we provide analytical support for the observed behaviors in both models.

Model A
A possible approach to determine the slope of the fixed point prevalences u, v, w with π is to differentiate the expressions A24 to A26 with respect to χ (notice, from the definition in equation 10, that χ is an increasing function of π). This procedure, however, can be very tedious and provides little or no analytical insight. An alternative approach is to implicitly differentiate the reduced equations A1 to A3 w.r.t. χ, obtaining more insightful expressions.
Let us define here u = ∂u/∂χ, v = ∂v/∂χ and w = ∂w/∂χ. Implicit differentiation of equations A1 to A3 yield: With some rearrangement, the above expressions can be written as a linear system given by: where − → x = (u , v , w ) T , and Moreover, the vector of independent coefficients is Using Cramer's rule, we can obtain the χ-derivatives as functions of the model parameters and the prevalences: Thus, in the asymmetrically interacting regime (0 ≤ Γ 1 < 1 and Γ 2 > 1), u has the same sign as w (thus u and w have the same slope with χ and π), whereas v (v) has opposite sign (slope). By showing that the ratio b w / det A is positive, we could demonstrate that u and w actually increase with π, while v decreases. Although we could not mathematically determine the signals of b w and det A for arbitrary model parameters, we collected robust numerical evidence that b w and det A are both negative in the coexistence phase (region IV in figure 3) for a wide set of model parameters, and so the ratio b w / det A is positive. This suggests that the behaviors shown in figures 4.a) and 5 are robust and should remain for the whole coexistence region.

Model B
For model B, one can extract the dependence of the prevalences with π (or χ) directly from their analytical expressions (eqs. A32 and A33), noticing that φ = 1 + αχλ 2 /λ 1 is an increasing function of χ. The prevalence v of disease II, as in equation A33, is clearly a decreasing function of φ for α > 1 (which is the asymmetrically interacting case). From equation A32, we can also directly infer that the prevalence u of disease I also decreases with φ for s 2 > 0 (or equivalently, λ 2 ≥ 1). However, the coexistence phase (region (IV)) also comprehends a region at which λ 2 < 1, for which we should check the behavior with φ more carefully. Taking the partial derivative of u with respect to φ, we get: From the above expression, the condition ∂u ∂φ < 0 can be simplified as: which is equivalent to the condition that λ 2 is above its critical value for coexistence, expressed by equation 8. This means that ∂u ∂φ < 0 in the whole coexistence region and, therefore, both prevalences v and u are decreasing functions of φ, χ and π. This is consistent with the observed behaviors in figures 4.b) and 6.