Phase diagrams of interacting spreading dynamics in complex networks

Liming Pan ,1,2 Dan Yang ,2 Wei Wang ,3,2,* Shimin Cai,2,4,5 Tao Zhou,2,4,5 and Ying-Cheng Lai 6 1School of Computer Science and Technology, Nanjing Normal University, Nanjing, Jiangsu 210023, China 2Complex Lab, School of Computer Science and Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China 3Cybersecurity Research Institute, Sichuan University, Chengdu 610065, China 4Institute of Fundamental and Frontier Sciences, University of Electronic Science and Technology of China, Chengdu 610054, China 5Big Data Research Center, University of Electronic Science and Technology of China, Chengdu 610054, China 6School of Electrical, Computer and Energy Engineering, Arizona State University, Tempe, Arizona 85287, USA


I. INTRODUCTION
Spreading dynamics of diseases, behaviors, and information in nature and human society are rarely independent processes but interact with each other in a complex way. Weakened immunity to other viruses due to HIV infection [1,2] and suppression of spreading due to disease-related information exchange on the social media [3] are known examples. To better understand, predict, and control spreading on networks, coevolution of epidemics must be taken into account. In network science, there has been continuous interest in developing interacting epidemic models [1,2,[4][5][6][7][8][9][10][11][12], which can generate surprising behaviors that cannot be predicted by any single-virus epidemic model. For example, spreading of one epidemic can facilitate that of another, leading to a firstorder or explosive transition in the outbreak with significant real-world implications [13]. Many factors can affect the critical behaviors of interacting spreading dynamics, such as self-evolution of each epidemic [3,14], interaction between two epidemics [15,16], and network structure [17,18].
By now, spreading dynamics of a single epidemic on complex networks have been well studied [8,[19][20][21]. For interacting spreading dynamics, the special case of well-mixed populations has been treated [22,23]. A study based on the quenched mean field for two competing pathogens [24] showed that, when simultaneous infection by the two pathogens is not possible (full mutual exclusion), the phase diagram is independent of the spectral radius of the network. There were also theories based on percolation [25], annealed mean field [26], and pair approximations [17] to study the effect of network structure on interacting spreading, leading to a qualitative understanding of critical phenomena. There are difficulties with these theories. For example, the annealed mean-field theory takes into account only the nodal degrees and is not applicable to quenched networks (especially networks with a high clustering coefficient and modularity). For such networks, quenched mean-field theories [27] such as those based on Markov chains [28] and the N-intertwined method [29] is needed. A deficiency of such a theory is that it uses a large number of nonlinear differential equations, with two difficulties: (a) high computational overload for large networks and (b) lack of any analytic insights. Such a theory, due to its heavy reliance on numerics, can lead to inconsistent or even contradicting predictions [30,31]. To our knowledge, a general analytic theory capable of providing a more complete understanding of interacting spreading dynamics is lacking.
In this paper, we develop an analytic theory for interacting spreading dynamics on complex networks through the approach of dimension reduction for complex networks [32][33][34]. From the quenched mean-field equations, we derive a two-dimensional (2D) system that is capable of analytically yielding the full phase diagrams underlying interacting spreading dynamics on any complex network, from which the conditions for various phase transitions can be derived. The analytic model predicts critical phenomena that were previously known numerically, such as the interplay between discontinuous outbreak transitions and hystereses as well as the emergence of tricritical points, providing a solid theoretical foundation for understanding interacting spreading dynamics and articulating optimal control strategies.

A. Model
We consider the susceptible-infected-susceptible (SIS) model of interacting spreading dynamics on complex networks. In the classic SIS model, a single epidemic spreads in the network and a node can be either in the susceptible or in the infected state. Susceptible nodes are infected by their infected neighbors at rate λ and infected nodes recover at rate γ . For interacting SIS dynamics, two epidemics, say 1 and 2, spread simultaneously and interact with each other. Each node infected by a ∈ {1, 2} transmits the infection to neighbors that are susceptible for both epidemics with probability λ a . If a neighbor is susceptible for a but infected by the other epidemic, then the infection will be transmitted with rate λ † a . All the nodes infected by a recover to being susceptible with rate γ a . Without loss of generality, we set γ a = 1 for both a ∈ {1, 2}. In general, the nature of the interacting SIS dynamics depends on the interplay between the rates λ a and λ † a . In particular, for λ † a > λ a , the two epidemics tend to facilitate each other, leading to cooperative SIS dynamics, whereas if λ † a < λ a , infection with one epidemic will suppress infection with the other, giving rise to competitive SIS dynamics. For λ † a > λ a but λ † b < λ b for a, b ∈ {1, 2} with b = a, the interactions are asymmetric.

B. Spectral dimension reduction
The interacting SIS model represents a paradigm to study rich dynamical behaviors such as first-order outbreak transitions and hystereses [35]. The foundation of our study of this model is the quenched mean-field theory (QMF) [36]. While in QMF the dynamical correlations among the neighbors are assumed to be negligible, the theory has been demonstrated to generate reliable prediction of the phase transitions [37]. Since our goal is to analytically map out the complete phase diagram, using the QMF suffices. Let p a,i be the probability that node i ∈ {1, · · · , N} is infected by a ∈ {1, 2} at time t. In the first-order mean-field analysis [27], the evolution of p a,i on a network with adjacency matrix G is governed by for a ∈ {1, 2} and i ∈ {1, · · · , N}. The first term on the right side of Eq. (1) is the rate of recovery from epidemic a for node i, while the second (third) term corresponds to the rate of infection for epidemic a with (without) i already infected by b = a. For a network of size N, the number of equations in Eq. (1) is 2N. To derive an analytic model, we exploit the technique of spectral dimension reduction (SDR) [33] to arrive at an equivalent description of the original system in terms of two macroscopic observables-one for each epidemic. In particular, let α be a vector with nonnegative entries and normalized as i α i = 1. The entries of α represent the nodal weights. We define linear observables as ψ a = α T p a for a ∈ {1, 2}. Since the entries of α are summed to unity, ψ a is a weighted average. The evolution of ψ a is determined by the equation Applying the SDR method, we have that the right-hand side of Eq.
(2) can be written in terms of the macroscopic observables ψ a only.

C. Reduced model
We apply the SDR method to Eq. (2). Microscopic variables p a,i fluctuate about the macroscopic observables ψ a , which can be decomposed as where ρ a , μ a , and ν a are parameters to be determined, and δ p a,i , δ p b,i , and δ p a, j are correction terms. Substituting Eqs. (1) and (3) into Eq. (2) gives (4) whereα = i, j α i G i j and R a is the remainder term that can be decomposed as with R a,1 , R a,2 , and R a,3 containing the first-, second-, and third-order terms in the corrections {δ p a,i }, respectively. Let K be the diagonal matrix with K ii being the degree of node i, the first-order correction R a,1 is given by where δ p a is a vector with δ p a,i in the ith entry and δ p b is defined analogously. The second-order remainder term is and the third-order remainder term is From Eq. (6), the dominant remainder term R a,1 vanishes if the following equations hold where p a is a vector with p a,i in the ith entry and p b is defined similarly.
In general, the equations cannot be satisfied simultaneously. An application of the SDR method in Ref. [33] advocates choosing α as the eigenvector associated with the leading eigenvalue ω of G. For connected undirected networks, the eigenvector associated with the leading eigenvalue ω of G has positive entries. The third equation in Eqs. (9) implies ωψ a =αν a ψ a and, hence, ω = ν aα . Using the definitionα = 1 T Gα = ω, we have ν a = 1. The remaining two parameters, ρ a and μ b , are chosen such that the first two equations in Eqs. (9) are satisfied. The quantities ρ a and μ b can be chosen by minimizing the following squared vector norm which yields A justification of the parameter choices was given in Ref. [33].
With the parameters chosen, R a,1 can be made as small as possible and can be neglected, so can the higher-order terms R a,2 and R a, 3 . Substituting the values of ρ a , μ a , and ν a into Eq. (4), we get The first term on the right side of Eq. (11) accounts for the rate of recovery and the second (third) term represents the rate of infection for epidemic a with (without) being infected by b = a. The quantity R a in Eq. (11) characterizes the fluctuations of the microscopic observables p a,i about the macroscopic observables ψ a , which is small in comparison to other terms on the right side of Eq. (11) due to α's being the leading eigenvector. Since, for finding the phase diagram, it is necessary to analyze the mean-field equations that depend on the macroscopic observables ψ a only, it is justified to drop R a from the analysis. As we will verify numerically, this approximation will not affect the accuracy of the phase diagram as the resulting errors near the phase boundaries are quite insignificant. In Eq. (1), the order parameters are p a,i for a ∈ {1, 2}, where · is the unweighted average over the nodes. In Eq. (11), the order parameters can be chosen to be ψ a = α T p a , a weighted average over the nodes. Since α is the eigenvector associated with the leading eigenvalue of G, its entries are strictly positive. As a result, ψ a = 0 (ψ a > 0) implies p a,i = 0 ( p a,i > 0). When crossing a phase boundary, at least one of ψ a for a ∈ {1, 2} becomes either zero or nonzero, guaranteeing that the corresponding p a,i becomes either zero or nonzero, respectively. We have that p a,i and ψ a give the same phase diagram, which can be obtained analytically through the 2D mean-field system.

III. MAIN RESULT: PHASE DIAGRAM OF REDUCED SYSTEM
The reduced mean-field equations are amenable to analytic treatment. As the derivations are lengthy, we provide a brief sketch of the results from analyzing the reduced system.
The analyses of the 2D mean-field system are performed in the following steps. First, for each point in the parameter space (λ 1 , λ † 1 , λ 2 , λ † 2 ), we determine the equilibrium points of Eq. (11) (Sec. III A) and their stability (Sec. III B). The the equilibrium points have to further satisfy the probability constraint 1 ψ a 1 to be physical meaningful. A detailed analysis of the stability and probability constraints of the equilibrium points leads to the following functions of λ a and λ † a : for a, b ∈ {1, 2} and a = b. Whether an equilibrium point is physical or stable is determined by the signs of these functions.
Calculating the equilibrium points (Sec. III A) and analyzing their stability (Sec. III B) enable us to obtain the full phase diagram of the reduced system and the equations for the phase boundaries (Sec. III C). Based on the numbers of stable and unstable equilibrium points as well as the relationships among them, we can divide the parameter space into distinct regions, where a boundary crossing between two neighboring regions gives rise to a phase transition. A region either can have a unique stable equilibrium point or can have two stable equilibrium points with one unstable point in between, where crossing the latter will result in a hysteresis. We discuss the types of phase transitions crossing the various boundaries and study the interplay between the transitions and the phenomenon of hysteresis (Sec. III D). Finally, we derive the conditions under which a hysteresis can arise (Sec. III E).
The analytical results of the full phase diagram are summarized and discussed in Sec. III F. Readers who are not interested in the technical details of analyzing the 2D mean-field system can skip Secs. III A to III E and check Sec. III F for the results. For convenience, for the rest of the paper, we use the convention that, if variables indexed by a, b ∈ {1, 2} (e.g., ψ a and ψ b ) appear together in an equation or an inequality, then the assumption is a = b.

A. Equilibrium points of the reduced system
The equilibrium points are obtained by setting the right side of Eq. (11) to zero: for a, b ∈ {1, 2} and a = b. Further, the physical solutions have to satisfy the probability constraints 0 ψ a 1. Because of the appearance of terms such as (1 − μψ a ) in Eq. (11), it is necessary to impose the physical condition (1 − μψ a ) 1. It can be proved that μ given by Eq. (10) satisfies μ 1 (see the Appendix for a proof), and it can also be verified that any point with ψ 1 = μ −1 or ψ 2 = μ −1 cannot be an equilibrium point. These, together the probability constraints, imply that all the equilibrium points must satisfy the inequality 0 ψ a < μ −1 for a ∈ {1, 2}.
We are now in a position to discuss the types of equilibrium points of the reduced mean-field equations.
(iii) Partial infection of epidemic 2. Similar to case (ii), we have the equilibrium point: (iv) Coexistence. If ψ a = 0 for both a ∈ {1, 2}, then Eqs. (13) become Substituting this relation into Eq. (17), we get an equation that depends on ψ 1 only. Similarly we can obtain the equation that determines ψ 2 . The two equations for ψ 1 and ψ 2 have the following symmetric form: for a ∈ {1, 2}, where If λ † a = λ a holds for a ∈ {1, 2}, then g a,2 = 0 and Eq. (19) will be quadratic, leading to two solutions, The solutions for a ∈ {1, 2} are paired as If λ † a = λ a for one or both a ∈ {1, 2}, then we have To discuss the probability constraints of the equilibrium points, we consider the following cases.
(iv.1) If g a,2 = 0 for one of a ∈ {1, 2}, i.e., λ † b = λ b , then there is a unique solution given by (24) The probability constraints imply Further, if we have λ † a = λ a , then the two epidemics will become independent of each other with the solution In this case there are two solutions, as shown in Eqs. (22).
Summarizing the above discussions about the equilibrium points, we have that (ψ 1 , ψ 2 ) = (ψ + 1 , ψ + 2 ) is physical for g a,0 < 0, g b,0 < 0 and g b,1 > 0. Note that g b,1 > 0 is implied by the other two. Since we have that g b,1 > 0 always holds given g b,0 < 0. Together, it is sufficient to have g a,0 < 0 and g b,0 < 0. We consider the solution (ψ 1 , . The probability constraints imply which cannot hold since its left side equals μωλ † a . Thus, in this region, no physical solution of (ψ 1 , , we must have g a,1 > 0 and g a,0 < 0 for a ∈ {1, 2}. Similar to the discussions in the case (iv.3), we have that the sufficient condition for an equilibrium point is g a,0 < 0 for a ∈ {1, 2}. The solution (ψ 1 ,

B. Stability analysis
The starting point to study the stability of the equilibrium points is the Jacobian matrix J of the 2D mean-field system, whose entries are We analyze the stability of the different classes of equilibrium points as discussed in Sec. III A.
(ii) Partial infection of epidemic 1. In this case, we have 023233-5 and J 21 = 0, so the Jacobian is upper triangular, whose eigenvalues are simply the diagonal entries: The equilibrium point is stable for J 1,1 < 0 and J 2,2 < 0, i.e., (iii) Partial infection of epidemic 2. For this type of equilibrium point, we have It is stable under the following conditions: (iv) Coexistence. Suppose we have ψ a = 0 for a ∈ {1, 2}. Substituting Eqs. (17) into Eqs. (44), the Jacobian matrix has entries A necessary and sufficient condition for a two-dimensional matrix to have two negative eigenvalues is to have a negative trace (tr(J ) < 0) but a positive determinant (det(J ) > 0). Since μψ a − 1 < 0, the negativity of the trace always holds. The stability of a equilibrium point in this class is fully determined by the determinant. It is stable when det(J ) > 0 and unstable when det(J ) > 0. The stable condition from the determinant is , the inequality can be written as Equations (17) can be rearranged as Multiplying the above two equations, we get where d 2 = 1, Multiplying both sides of Eq. (55) by z and substituting the result into Eq. (53), we obtain For λ † a = λ b for any a ∈ {1, 2}, from Eq. (55), we see that 1/z has two solutions Since we have a pair of solutions for (ψ 1 , ψ 2 ) as in Eq. (22), the following hold: Substituting Eq. (58) into Eq. (57), we have We see that, given d 2 is always unstable. It remains to check the validity of the inequality d 2 1 − 4d 2 d 0 > 0. After some algebra, we have Thus, the inequalities d 2 1 − 4d 2 d 0 > 0 and g 2 a,1 − 4g a,2 g a,0 > 0 are equivalent to each other for a ∈ {1, 2}.

C. Phase diagrams
With full knowledge about the equilibrium points and their stability, we can obtain the phase diagram of the reduced mean-field system. Define the following set of functions: for a, b ∈ {1, 2} and a = b. The distinct phase regions can be defined via various inequalities among these functions.
(ii) Partial infection of epidemic 1. The phase has a stable equilibrium point Combining the probability constraints and the stability analysis, we obtain the phase region as (iii) Partial infection of epidemic 2. The phase is characterized by The phase region is given by (iv) Coexistence. In this region, there is an equilibrium point with both ψ 1 and ψ 2 nonzero, corresponding to the case of double epidemic outbreaks. For cooperative coevolution, i.e., λ † a > λ a for a ∈ {1, 2}, a point in the parameter space belongs to this phase if or s a,0 < 0, s a,2 > 0, for both a ∈ {1, 2}. When coevolution is not cooperative, the coexistence region is given by for both a ∈ {1, 2}. We have verified that the case of λ a = λ † a for one or both a ∈ {1, 2} is well covered by this inequality.
(i ∩ iv). Hysteresis region 1. A hysteresis region appears when there are two stable equilibrium points and one unstable equilibrium point in between. The stability analysis indicates that the solution (ψ + 1 , ψ + 2 ) is always stable while (ψ − 1 , ψ − 2 ) is unstable. In addition to these two equilibrium points, a third stable solution is necessary for a hysteresis to arise. This is only possible when region (iv) overlaps with regions (i), (ii), and (iii). Checking the equilibrium points and their stability, we find that a hysteresis region exists only when the inequality λ † a > λ a holds for a ∈ {1, 2}. The region where (i) and (iv) overlap is s a,0 > 0, s a,1 < 0, s a,2 > 0, s a,3 < 0, s > 0, where the first inequality s a,0 > 0 can in fact be implied by the other inequalities. Since g 1,1 + g 2,1 < 0, we have which further implies ωλ † 1 > 1 and ωλ † 2 > 1. Since s a,3 < 0, we have Altogether, the region is given by s a,1 < 0, s a,2 > 0, s a,3 < 0, s > 0 (73) for a ∈ {1, 2}.

D. Types of phase transition
A phase transition occurs when a point in the parameter space crosses a boundary between two neighboring phase regions. Depending on different combinations of phase-region pairs, the resulting phase transitions can be characteristically distinct. To be concrete, we focus on the phase transitions in the λ 1 -λ 2 plane with fixed values of λ † 1 and λ † 2 . Both continuous and discontinuous phase transitions can arise, as we will show below.
To summarize, if a point in region (ii) has s 2,1 < 0 near the phase boundary s 2,0 = 0, then all the conditions under which the point is in region (ii ∩ iv) hold. Thus, if the point is in the region (iv)\(ii ∩ iv), then we have s 2,1 > 0, which makes the phase transition continuous.
(iii) (iv)\(iii ∩ iv). Following a similar treatment to the preceding case, we have that the phase transition is continuous.
(ii ∩ iv)→(iv). The two phase are separated by the line s 2,0 = 0. As discussed in the case of the (ii) (iv)\(ii ∩ iv) transition, since s 2,1 < 0 holds near the phase boundary, the behavior of the coexistence solution is determined by Eq. (78) when approaching the phase boundary, resulting in a discontinuous phase transition.
The tricritical points that separate the continuous from the discontinuous transition lie in the boundaries of the hysteresis region where q a,1 = 0 holds for either a ∈ {1, 2}. One such point is given by Eq. (80). The second tricritical point can be obtained similarly as

E. Conditions on λ † a for hysteresis
For fixed values of λ † 1 and λ † 2 , a hysteresis can arise in the λ 1 -λ 2 plane. To determine these values, we first note that a hysteresis is possible only when the coevolution dynamics are cooperative, i.e., λ † a > λ a for both a ∈ {1, 2}. A point in the hysteresis region must satisfy the inequality g 1,1 + g 2,1 < 0. Consequently, we have which provides a necessary condition for a hysteresis to arise. We can show that this is also sufficient to guarantee the occurrence of a hysteresis. In particular, suppose inequality Eq. (102) holds. Since s 1,0 = 0 and s 2,0 = 0 intersect at (λ 1 , λ 2 ) = (ω −1 , ω −1 ) in the λ a -λ b plane, there is a neighborhood near (ω −1 , ω −1 ) in which the inequalities s 1,0 > 0 and s 2,0 > 0 hold. At the point (λ a , λ b ) = (ω −1 , ω −1 ), we have We thus have that a hysteresis region exists in the λ 1 -λ 2 plane if and only if Eq. (85) holds.

F. Summary of the phase diagrams
Based on the results of the above analysis, we obtain the structure of the analytically predicted phase diagram, which can be described, as follows.
(ii) Partial infection of epidemic 1. In this phase, an outbreak occurs for the epidemic 1 but not for 2. The equilibrium point is given by which is stable for s 2,0 > 0 and s 1,3 > 0, where the latter gives λ 1 > ω −1 , indicating that epidemic 1 can have an outbreak independently. Similarly, s 2,0 > 0 implies λ 2 < λ 1 and stipulating that λ † 2 cannot be too large, such that the outbreak of epidemic 1 result in an outbreak in epidemic 2.
(iii) Partial infection of epidemic 2. Analogous to (ii), this phase is defined by which is stable for s 1,0 > 0 and s 2,3 > 0.
(iv) Coexistence. In this region the reduced system has a stable equilibrium point with both ψ 1 and ψ 2 nonzero, leading to a simultaneous outbreak of two epidemics. The stable equilibrium points are For cooperative coevolution, i.e., λ † a > λ a for a ∈ {1, 2}, we have s a,2 > 0. A point in the parameter space belongs to this phase if it further satisfies either or s a,0 < 0 (93) for a ∈ {1, 2}. The former case is where region (iv)overlaps with regions (i), (ii), and (iii). As a result, hystereses can arise. For competitive or asymmetric coevolution, the coexistence region is given by s a,0 < 0 for a ∈ {1, 2}. Now we discuss the conditions for observing the coexistence phase in more detail to get an intuitive picture. For the cooperative case, the boundaries of the coexistence region are relatively complex, and we consider the degenerate case of λ 1 = λ 2 and λ † 1 = λ † 2 . For s a,0 < 0, we have ωλ 1 = ωλ 2 > 1 so that both epidemics are able to have an outbreak by themselves. On the contrary, for s a,0 > 0 so that neither epidemic can have an outbreak by itself, we have s a,1 < 0 and λ † 1 > 2λ 1 . The condition s > 0 requires and λ † 1 > 2λ 1 , leading to the necessary condition ωλ † 1 > 2. In this case, the interaction transmission rate must at least double the threshold value of classic SIS outbreak to have coexistence.
For the competitive case, since s a,0 < 0, we have The above inequality implies that at least one of λ a , say λ 1 , must satisfy ωλ 1 > 1. Actually, we must also have ωλ 2 > 1 to observe the coexistence phase. Suppose ωλ 2 < 1, from ωs 1,0 < 0, we have which contradicts with the assumption that the model is competitive. Consequently, to observe the coexistence of two epidemics for the competitive case, a necessary condition is that each of the two epidemics can have an outbreak by itself.
In addition, s a,0 < 0 implies That is to say, to have the coexistence phase, the suppression effect from the other epidemic cannot be too strong. For networks with a larger value of ω, the conditions Eq. (97) and ωλ a > 1 both can be readily satisfied. As a result, the coexistence phase is more likely to be observed in networks with a larger leading eigenvalue ω. It is also interesting to note that, when the two epidemics have fully mutual exclusion (i.e., λ † a = 0), the condition s a,0 < 0 can never be satisfied simultaneously for both a ∈ {1, 2}. In other words, the coexistence phase can never be observed with fully mutual exclusion, and this phenomenon agrees with the prediction in Ref. [38].
(i ∩ iv) Hysteresis region 1. A hysteresis arises when there are two stable equilibrium points and one unstable equilibrium point in between, which occurs when region (iv) overlaps with regions (i), (ii) and (iii). Our analysis reveals that a hysteresis region emerges only for cooperative coevolution. The region where (i) and (iv) overlap is bounded by the inequalities for a ∈ {1, 2}.
(iii ∩ iv) Hysteresis region 3. Similarly, for cooperative coevolution with s a,2 > 0, the region where (iii) and (iv) overlap is bounded by The types of phase transitions that occur when crossing a phase boundary are determined by further checking if the stable solution varies continuously (see Sec. III D). We find all possible phase transitions as a result of crossing a hysteresis region are discontinuous, whereas other transitions are continuous. A result revealed by our analysis of the phase diagrams is that the precursor of a discontinuous transition with an abrupt outbreak of at least one epidemic is a hysteresis. Continuous and discontinuous phase transitions are separated by two tricritical points in the λ 1 -λ 2 plane: The phase diagram also makes it possible to obtain the conditions in the interaction strengths λ † 1 and λ † 2 for a hysteresis to occur. In Sec. III E, we obtain the necessary and sufficient condition where there is a hysteresis region with λ 1 < λ † 1 and λ 2 < λ † 2 .

IV. NUMERICAL VERIFICATION
We provide a numerical illustration of the analytic prediction on the interplay between discontinuous transitions and hystereses with an Erdős-Rényi graph of size N = 100 and average degree k = 4. The inequality Eq. (102) divides the λ † 1 -λ † 2 plane into two regions, as shown in Fig. 1(a). Above the curve defined by a hysteresis region appears while it is absent below. In the limit λ † a → ∞, the curve approaches λ † b = ω −1 , as shown by the orange dashed lines in Fig. 1(a). Note that the curve avoids the dashed lines for finite λ † a . Since ω −1 is also the outbreak threshold of the classic SIS model for a single epidemic, a necessary condition for a hysteresis is that λ † a must be larger than the classic threshold. A special case is λ † 1 = λ † 2 , where Eq. (102) implies that, for a hysteresis to arise, the inequality λ † 1 = λ † 2 > 2ω −1 must hold. That is, the interactive transmission rate must at least twice the classic SIS threshold for a hysteresis to arise, suggesting that networks with a larger leading eigenvalue are more prone to hystereses. Two representative phase diagrams in the λ 1 -λ 2 plane with fixed values of λ † 1 and λ † 2 are shown in Figs. 1(b) and 1(c), corresponding to the points b and c in Fig. 1(a), respectively. For point b, no hysteresis can arise for any values of (λ 1 , λ 2 ) and the phase transitions between different neighboring phase regions are continuous, as indicated by the solid lines in Fig. 1(b). For point c that is slightly above the hysteresis boundary, region (iv) overlaps with regions (i), (ii), and (iii), where a hysteresis can arise. Crossing into region (iv) from any one of the phase regions (i ∩ iv), (ii ∩ iv), and (iii ∩ iv), a discontinuous outbreak transition occurs with some ψ a changing abruptly from zero to a nonzero value. Along the path (i ∩ iv) → (i), (ii ∩ iv) → (ii), and (iii ∩ iv) → (iii), the system displays a discontinuous transition to extinction at which at least one epidemic changes abruptly from a nonzero value to zero. All the phase boundaries with discontinuous transitions are indicated by the dashed lines in Fig. 1(c), where the two tricritical points separating continuous from discontinuous transitions are marked (white circles).
Are the phase diagrams obtained from the reduced mean field equations accurate in comparison with those from the original mean field equations? In the presence of the fluctu-ation terms R a , Eq. (11) are exactly equivalent to Eqs. (1). Consider a system of dimension 2N + 2, which consists of Eqs. (1) and (11). A stable equilibrium point (p 1,i , p 2,i ) 1 i N of the subsystem Eqs. (1) is also one for the 2N + 2 system with ψ a = α T p a . Consider a stable equilibrium point with which neither epidemic has an outbreak. Substituting p a,1 = · · · = p a,N = 0 and ψ a = 0 into the remainder term (with full expression in Sec. II B), we have R a = 0 for a ∈ {1, 2}. In this case the remainder terms can be ignored. Since a zero stable equilibrium point of Eqs. (1) implies the existence of exactly such a point of Eq. (11) (with no remainder terms) and vice versa, any outbreak transition threshold from phase (i) is expected to be exact for Eqs. (1).
There are two cases where the remainder terms do not vanish and can lead to inaccuracies of the analytic prediction. The first case is when Eqs. (1) exhibit a stable equilibrium point at which there is an outbreak for epidemic 1 but extinction for epidemic 2: R 1 = 0 and R 2 = 0. The second case is when Eqs. (1) have a stable equilibrium point with an outbreak for both epidemics: R a = 0 for a ∈ {1, 2}. Since the remainder terms are small by construction, they lead to corrections that can be neglected, which have been verified numerically. Especially, for the Erdős-Rényi network in Fig. 1, we solve Eqs. (1) numerically and compare the solutions with the analytic phase diagram obtained from Eq. (11). The values of ψ a obtained from Eqs. (1) in the λ 1 -λ 2 plane are shown in Figs. 2(a)-2(d), for λ † 1 = 3.5ω −1 and λ † 2 = 2.5ω −1 [so that (102) is satisfied, guaranteeing a hysteresis]. Since in the hysteresis region there are two stable equilibrium points for each ψ a , we plot separately the two solutions for ψ 1 in Figs. 2(a) and 2(b), and those for ψ 2 in Figs. 2(c) and 2(d), respectively. The phase diagram from original Eq. (11) is also shown in Fig. 2 by the orange solid and dashed lines for continuous and discontinuous transitions, respectively. Our analytical phase diagram predicts accurately all outbreak transitions: (i) → (ii), (i) → (iii), (i ∩ iv) → (iv), (ii ∩ iv) → (iv), and (iii ∩ iv) → (iv). However, quantitatively, the predicted extinction transitions (i ∩ iv) → (i), (ii ∩ iv) → (ii), and (iii ∩ iv) → (iii) are less accurate, due to the nonzero remainders R 1 and R 2 Next we consider tests and validation of our analytic prediction from Eq. (11) for a variety of networks, including synthetic networks with strong and weak degree heterogeneity, and real-world networks. For synthetic networks, we have already shown the results for an ER network. Here we also consider networks generated from the uncorrelated configuration model (UCM) with a power-law degree distribution p(k) ∼ k −β . Specifically, we consider three networks with different degree exponents: (1) PL-2.3 with β = 2.3, (2) PL-3 with β = 3, and (3) PL-4 with β = 4. For empirical networks, we have (4) Dolphins [39], a social network of bottle-nose dolphins; (5) HIV [40], a network of sexual contacts between people involved in the early spread of the human immunodeficiency virus (HIV); (6) Highschool [41], a friendship network between boys in a small high school, and (7) Jazz [42], a collaboration network between Jazz musicians. The networks are downloaded from Ref. [43].
Basic features and parameters of the networks considered are listed in Table I. Note that Highschool is a directed and weighted network. Here we simply take it as undirected by assuming that there is an undirected edge between node i and j if there is at least a directed edge between the two nodes in either direction, with the edge weights ignored.
For all the networks, we set λ † 1 = 3.5ω −1 and λ † 2 = 2.5ω −1 to guarantee the emergence of a hysteresis region. To have an idea of the size of the correction terms R a , we show the values of R a at equilibrium. The results of (1) PL-2.

V. DISCUSSION
We have analytically predicted the phase diagram of interacting SIS spreading dynamics using the technique of spectral dimension reduction and provided numerical validation. The analytic phase diagram elucidates the interplay between discontinuous transitions and hystereses as well as the emergence of tricritical points. This method can also be applied to study other interacting epidemic models. For general epidemic models, a one-dimensional description of each epidemics is not sufficient [33]. Determining the number of macroscopic observables required for general epidemic models needs further exploration. Previous theoretical methods for interacting spreading dynamics such as QMF theory [36] employ 2N equations, where N is the network size. For large networks, it is computationally demanding to solve the equations to determine the stability of the fixed points, as this requires manipulating the 2N × 2N Jacobian matrix. It is thus infeasible to use the QMF to map out the phase diagram for interacting spreading dynamics on complex networks, preventing us from gaining a full understanding of the interplay between network topology and the spreading dynamical process as a full phase diagram would reveal. The same difficulty arises with a naive application of the SDR method [33] to obtain the phase diagram for interacting spreading dynamics. In contrast, our approach gives a full picture of the phase diagram on large complex networks with an arbitrary topology through an effective twodimensional system, revealing rich phenomena that have not been systemically investigated. While many previous studies employed the mean-field theory to study different types of spreading dynamics on complex networks [46][47][48], our work is not a simple application of the mean-field theory. In fact, we go way beyond by obtaining, for the first time to our knowledge, a global phase diagram laying out a clear picture of all possible dynamical states and the transitions among them through a comprehensive stability analysis-both at an unprecedented level of details. Taken together, our work gives a full picture of the dependence of phase transition on network topology and spreading parameters for SIS dynamics, and thus lays a foundation for intervening or harnessing this type of interacting spreading processes. For instance, our phase diagram gives possible routes for controlling the type of phase transition through perturbations to the network structure or for controlling one spreading process through manipulating another interacting process. It should be cautioned that, while the SIS model provides phenomenological insights into relatively simply spreading processes and is thus a conceptually useful paradigmatic model, it may be too simplistic to describe spreading processes in the real world which can be significantly more complicated. To apply our analytic approach to irreversible epidemic processes beyond the SIS dynamics is possible but remains to be studied. Since K is positive definite, it can be decomposed as K = K 1/2 K 1/2 , where K 1/2 is a diagonal matrix whose entries are the square root of the degrees. Let y = K 1/2 α. We have α = K −1/2 y. Substituting this back to μ −1 gives which is the Rayleigh quotient of matrix K −1/2 GK −1/2 and, hence, we have μ −1 δ 1 , where δ 1 is the largest eigenvalue of the matrix K −1/2 GK −1/2 . Recall that the symmetric normalized Laplacian matrix of G is defined as which has a smallest eigenvalue ζ n = 0. As a result, we have δ 1 = 1 − ζ n = 1, which gives μ 1.