also studied the household model, in which each vertex has a graded state corresponding to the number of patients in a household [26] . They showed that the epidemic thresholds are proportional to k / k 2 also in these models. However, how general the universality of the k / k 2 -scaling is has not been addressed.
In this paper, we analyze various multi-state epidemic processes originally proposed for lattices or populations with homogeneous contact rates. We concentrate on endemic models in which the phase boundary between the coexistence state, where both susceptible and infected individuals can survive, and the disease-free state with the susceptible individuals only is concerned. Of course, we then disregard the other important class of epidemic models including the SIR model, in which possibility of one-shot waves of outbreaks rather than coexistence is concerned (for this distinction, see [5, 16] ). Instead of covering both classes, we investigate phase boundaries of endemic models with more than two states and possibly with more than two phases. We show that the critical infection rates are largely extinct as k 2 → ∞. However, there are a couple of important exceptions to this rule: models with mutation or cyclic interaction. In Sec. II, we start with studying two models that follow the k / k 2 -scaling. Sections III, IV, V, and VI are devoted to the analysis of representative models with different characteristics and outcomes. Some of them do not follow the k / k 2 -scaling. The results are summarized and discussed in Sec. VII.
In the following, we investigate various contagion models with different transition rules.
Because we can never be exhaustive, our strategy is to study known multiple-state models.
It will turn out that the extinction of epidemic thresholds as k 2 / k goes to infinity is not a universal story, as summarized in Sec. VII. What determines the nature of phase boundaries is considered to be the gross organization of the transition pathways, some of which are independent of the state of the neighbors (solid lines in Figs. 1, 2, 3, 5, 7) and others are neighbor-dependent (dashed lines). For the introductory purpose, we analyze in this section two three-state models whose epidemic thresholds are proportional to k / k 2 .
A. Two-stage Contact Process
The state-transition rules of the two-stage contact process [22] are schematically shown in Fig. 2(A) . It differs from the SIRS model ( Fig. 1(C) ) in that the transition rate from 0 to 1 is proportional to the density of the neighbors with state 2 (not 1) of a state-0 vertex.
This model corresponds to dynamics with two life stages. The states can be interpreted as 0: vacant, 1: occupied by young individuals, and 2: occupied by adults. Only adults are reproductive and generate offsprings in neighboring vacant sites at a birth rate equal to λ. In other words, a birth event occurs at a vacant site at a rate proportional to λ and the number of neighboring adults. Youngs (state 1) spend random time of mean r before becoming adult (state 2). They are also subject to random death events at a rate of δ. Adults die at a rate of 1, which gives normalization of the entire model. Alternatively, we can interpret the three states as 0: vacant, 1: partially occupied, and 2: fully occupied colonies. Then, only fully occupied colonies are potent enough to colonize vacant lands. In the disease analogue, 0 is susceptible, 1 is infected but in the incubation period (not infectious), and 2 is infected and infectious. We note that the model corresponds to the SEIRS model without the recovery state R.
The ordinary contact process ( Fig. 1(A) ) is recovered as r = ∞. In this case, a state-1 vertex immediately turns into state 2. The meanfield approximation and the rigorous results on regular lattices conclude that there are two phases: {0} (the steady state that consists only of state 0) and {0, 1, 2} (positive probability of steady coexistence of 0, 1, and 2). Naturally, large λ or r promotes survival [22] . Particularly, in the meanfield approximation, these two phases are divided by rλ = δ + r.
To derive meanfield dynamics for populations with heterogeneous contact rates, let us denote by p k the probability that a vertex has degree k. Obviously, ∞ k=1 p k = 1. We denote by ρ i,k (i = 1, 2) the probability that a vertex with degree k takes state i. The probability that a vertex with degree k takes state 0 is equal to 1 − ρ 1,k − ρ 2,k . The probability that a neighbor located at the end of a randomly chosen edge takes state i is denoted by Θ i . This edge-conditioned probability, particularly Θ 2 in this model, is used to determine the effective infection rate [36, 37] .
The dynamics are given bẏ
We note that the effective birth rate (the first term in the first equation of Eq. (1)) is
proportional to kΘ 2 , which is the average number of state-2 vertices in the neighborhood of a degree-k vertex.
When we choose an arbitrary edge, the probability that a specific vertex is connected to this edge is proportional to its degree k [9, 10, 33, 36, 37] . Therefore, we obtain
The steady state is given by
where * indicates the steady state. Plugging Eq. (3) into Eq. (2) leads to
Equation (4) is satisfied when Θ * 2 = 0, corresponding to the disease-free state {0}. When 0 < Θ * 2 < 1, state 2 survives. Equation (3) implies ρ * 1,k , ρ * 2,k > 0 in this situation. Accordingly, Θ * 2 > 0 is equivalent to the {0, 1, 2} phase. Because the RHS of Eq. (4) is smaller than the LHS when Θ * 2 = 1, the condition for the {0, 1, 2} phase is
This yields
When p k ∝ k −γ (γ ≤ 3), the critical infection rate vanishes. This agrees with the results for the percolation, the contact process, and the SIR model.
Model A next example is the tuberculosis model shown in Fig. 2 (B) [45] . In the tuberculosis, a majority of the infected individuals does not become infectious before recovery. Accordingly, state 0 corresponds to healthy, 1 to infected but not infectious, and 2 to infected and infectious. This model, which is an extension of the two-stage contact process, has two important features. One is that we have two state-transition routes that depend on neighbors' states.
Both of them depend on the density of infectious individuals (state 2). Antoher feature is that disease progression from state 1 to 2 is possible through two parallel routes. An early-stage patient (state 1) can develop the disease on its own at a rate of r. At the same time, infectious patients (state 2) in the neighborhood increases the possibility of disease progression for a state-1 individual.
When r = 0, the {0, 1, 2} phase is possible on regular lattices but impossible in perfectly mixed populations [45] . When r > 0, the {0} phase and the {0, 1, 2} phase appear for both lattices and meanfield populations, depending on values of other parameters. According to the ordinary meanfield analysis, these phases are divided by λr = (δ + r), which is the condition identical to that for the two-stage contact process.
In heterogeneous populations, the dynamics becomė
The steady state is given by
and
The {0, 1, 2} phase appears when Θ * 2 > 0, or
This is identical to Eq. (6).
In many epidemics, multiple strains or pathogens with different transmissivity, virulence, and mobility coexist. They compete with each other by preying on common susceptible individuals. On top of that, transitions between infected states with different pathogens or stages of the disease can occur owing to mutation and disease progression. Coexistence of two distinct diseases in a population also introduces competition between them.
In terms of schematic diagrams, this means that multiple neighbor-dependent infection pathways (dashed lines) emanate from a susceptible state. We analyze three models with competiting pathogens and show that epidemic thresholds are not entirely governed by the degree-dependent factors in two models.
Let us consider a minimal model for competing pathogens whose transition rules are shown in Fig. 3 (A). State 1 and state 2 correspond to patients infected with strain 1 and strain 2, respectively. The two strains compete with each other by feeding on the common resource, namely, susceptible individuals (state 0), at rates β 1 and β 2 . Any patient spontaneously recovers at rate 1. Unidirectional and spontaneous mutations occur (1 → 2) at rate r [43] . As the mutation rate increases, strain 2 becomes more viable than strain 1. In the poupulation ecology context, the three states can be also interpreted as empty (0), species A (1), and species B (2).
Exploiting the fact that the absence of state-1 vertices is equivalent to the contact process, Schinazi showed for regular lattices that there are three phases {0}, {0, 2}, and {0, 1, 2} separated by nontrivial critical lines [43] . Small r and large β 1 support the {0, 1, 2} state.
The meanfield analysis predicts that the boundary between {0} and {0, 2} is β 2 = 1, and {0, 1, 2} appears when both β 1 > β 2 (r + 1) and β 1 > r + 1 are additionally satisfied.
For heterogeneous populations, the dynamics are given bẏ
The boundary between {0} and {0, 2} is obtained easily. Application of the contact-process result implies that {0, 2} emerges when
which extends the condition β 2 > 1 for homogeneous populations.
We next analyze the boundary between {0, 2} and {0, 1, 2}. With the steady state
we obtain
and
In the {0, 1, 2} phase, Θ 1 > 0 and Θ 2 > 0 are simultaneously satisfied. The condition Θ * 2 > 0 is satisfied when
Equation (16) extends one of the conditions for the {0, 1, 2} phase in the homogeneous population: β 1 > r + 1. Imposing Θ * 1 > 0 leads to β 1 > β 2 (r + 1), a condition independent of heterogeneity in contact rates. Because of this condition, k 2 → ∞ means survival of state 2 but not necessarily survival of state 1. Increasing k 2 / k makes the survival of state 1 and state 2 more or less equally likely because they feed on the susceptibles in the same manner. Then, the strengths of 1 and 2 must be balanced so that one does not completely devour the other. This constrains the range of the mutation rate in a degree-independent manner. Figure 4 shows numerically obtained steady densities of state 1 (A) and 2 (B) in networks with 10000 vertices and k = 12. Initially, each vertex takes one of the three states independently with probability 1/3. We set β 1 = 0.5 and β 2 = 0.25, which predicts the phase boundary: r = β 1 /β 2 − 1 = 1.
Scale-free networks with the degree distribution p k ∝ k γ with γ = 2.5 (thickest solid lines), γ = 3.0 (moderate solid lines), and γ = 4.0 (thinnest solid lines) are produced using a static model [14] . Even though produced networks can be disconnected in general, more than 95 % of the vertices constitute one component in every run of our simulations. The results for the random graph, whose p k follows the Poisson distribution, are also shown (dashed lines). Regardless of networks, the numerical results support the existence of finite thresholds in terms of r. The values of critical r are slightly smaller than r = 1 but do not vary so much for different p k , as our theory predicts.
In many diseases such as tuberculosis, drug-resistant strains emerge by the mutation of a wild strain [8] . Misuse of antibiotics is a major cause of appearance of such stronger strains. A four-state model shown in Fig. 3(B) represents spreads of a drug-resistant strain [42] . The states 0, 1, 2, and 3 represent empty, susceptible, infected with the wild strain, and infected with the drug-resistant strain, respectively. The wild strain (state 2) and the drug-resistant strain (state 3) compete with each other to devour the susceptible (state 1).
Only the drug-resistant strain is supposed to be fatal with death rate 1. When infected by the wild strain, which is not lethal, an individual recovers at rate r owing to the drug. The mutation rate is denoted by φ.
The model without state 2 (or β 2 = 0) is equivalent to the SIRS model in Fig. 1 (C) with δ = 0. On regular lattices, {0, 1, 2, 3} appears when r + φ is small enough. Otherwise, state 2 is extinguished [42] . Based on the meanfield analysis, the {0, 1, 2, 3} phase ensues when
If any of the conditions in Eq. (17) 
In heterogeneous populations, we obtaiṅ
The steady densities are given by
where
Equation (19) assures that ρ * 1,k > 0. In addition,
provided that β 2 /(r + φ) > β 3 . Under this condition, Θ * 2 > 0 is equivalent to the {0, 1, 2, 3} phase. Combining Eqs. (2), (19) and (21) results in
Since ∆ is of the form: ∆ = C 3 + C 4 kΘ * 2 , Eq. (22) is essentially the same as Eq. (4).
Even for large r + φ, sufficiently heterogeneous p k with large k 2 / k makes Eq. (23) valid.
However, for existence of state 2, β 2 (r + φ) > β 3 must be also satisfied. This condition is independent of the degree distribution and identical to one of the conditions for the {0, 1, 2, 3} phase in the ordinary meanfield case (see Eq. (17)). As is the case for the previous model (Sec. III A), this condition regulates the rates of mutation and spontaneous recovery so that the strengths of the wild strain and the drug-resistant strain are roughly balanced.
If β 2 /(r + φ) ≤ β 3 , then Θ * 2 = 0. In this case, Eq.
which yields
As a result, {0, 1, 3} and {1} are separated by
In summary, the {0, 1, 2, 3} phase dominates even for infinitesimally small infection rates when k 2 = ∞ and β 2 /(r + φ) > β 3 . When the latter condition is violated, {0, 1, 3} appears for a tiny infection rate.
Superspreaders, namely, patients that infect many others in comparison with normal patients, are identified in the outbreaks in, for example, gonorrhea [16] , HIV [4, 11, 29] , and SARS, [1, 23] . Such heterogeneous infection rates are sometimes ascribed to heterogeneous social or sexual contact rates as specified by p k [4, 16, 24, 29, 31, 46] . However, superspreaders exist even in diseases whose associated networks are considered to have relatively homogeneous k, such as SARS [27] . Then, an alternative way to introduce superspreaders is to consider two types of patients with different infection rates. Then, competition occurs between superspreaders and normal patients.
We assume that superspreaders and normally infectious patients are not distinguished by the strain types. Accordingly, an infection event originating from a superspreader can cause a normal patient and vice versa. There is no mutation, and competition occurs indirectly.
That is how the superspreader model considered here is essentially different from the models in Secs. III A and III B.
We analyze a three-state model with different infection rates, which is schematically shown in Fig. 3 (C) [18, 44] . The states 0, 1, and 2 mean susceptible, patients of type 1 with infection rate λ 1 and recovery rate 1, and patients of type 2 with infection rate λ 2 ≥ λ 1 and recovery rate δ, respectively. State 2 corresponds to the superspreader, and in practical appliations, we put λ 1 ≪ λ 2 . Upon infection, a state-0 vertex changes its state into 1 (2) with probability 1 − p (p). There are two neighbor-dependent infection routes in Fig. 3(C) . 
Otherwise, {0} results. Slightly weaker but qualitatively similar results have been proved for lattices [44] .
In heterogeneous populations, dynamics are given bẏ
The steady state is:
which leads to
Therefore, Θ * 1 > 0, or equivalently, Θ * ≡ λ 1 Θ * 1 + λ 2 Θ * 2 > 0, is the condition for the {0, 1, 2} phase. With Eq. (29), we obtain
and {0, 1, 2} emerges when
Equation (32) suggests that an arbitrary small infection rate λ 1 or λ 2 allows the endemic state as k 2 / k → ∞. In contrast to the models in Secs. III A and III B, degree-independent factor does not play a role in this model. This is due to the absence of mutation from 1 to 2 or from 2 to 1.
For a realistic network model of populations with natural birth, we have to take into account that birth can occur only at empty sites. Then, a minimal model such as the contact process must be extended to a three-state process in which states 0, 1, and 2 mean empty, susceptible, and infected, respectively. A susceptible gives birth to an offspring in a neighboring empty vertex, which is operationally similar to a contagion event. An infected preys on a susceptible in its neighborhood. Then there are two neighbor-dependent transition rates, namely, 0 → 1 transition at a rate proportional to n 1 and 1 → 2 transition at a rate proportional to n 2 . A similar situation arises in the context of double infection in which states 0, 1, and 2 correspond to susceptible, infected, and infected by another pathogen, respectively. A second pathogen (state 2) targets individuals infected by a first pathogen (state 1).
Keeping these interpretations in mind, we analyze the model shown in Fig. 5 [15, 39] .
State 2 may more virulent than state 1, and evolution of strains has been investigated with a more complicated version of this model [15] . The model extends the contact process and the SIR model, which correspond to µ = 0 and λ = δ = r = 0, respectively. 
Another paradoxical behavior occurs regarding to population density. For example, an increase in λ lessens the number of state-1 vertices on lattices because a larger λ creates more preys to devour for state-2 vertices [52] . Similarly, too large µ drives state 2 to perish due to excess mortality [15, 39] , which defines another {0, 1} phase that we do not examine here. This {0, 1} phase is distinct from the {0, 1} phase revealed by the meanfield equations and mathematical analysis on lattices.
This model is complex also in the sense that the phase diagram in the parameter space is even qualitatively unknown for lattices. The meanfield solution indicates that the two critical lines approach each other as µ → ∞. On the other hand, they may cross at finite λ and µ, as supported by the improved pair approximation ansatz [15, 39] .
In heterogeneous populations, the dynamics reaḋ
The steady state is given by
which leads to
Clearly, (Θ * 1 , Θ * 2 ) = (0, 0), which corresponds to the {0} phase, solves Eq. (36) . To explore other phases, let us set
and
Then Eq. (36) is equivalent to Θ * 1 f 1 (Θ * 1 , Θ * 2 ) = 0 and Θ * 2 f 2 (Θ * 1 , Θ * 2 ) = 0. First, we identify the boundary between {0} and {0, 1}. In the {0, 1} phase, we have Θ * 2 = 0. Under this condition, we look for 0 < Θ * 1 < 1 that satisfies f 1 (Θ * 1 , 0) = 0. By substituting Θ * 2 = 0 into Eq. (37), we obtain
Because
Eq. (39) is satisfied when f 1 (0, 0) > 0, that is,
This conclusion complies with the results for the contact process (Sec. I).
Second, we examine the boundary between {0, 1} and {0, 1, 2}. The {0, 1, 2} phase implies Θ * 1 > 0 and Θ * 2 > 0 such that f 1 (Θ * 1 , Θ * 2 ) = 0 and f 2 (Θ * 1 , Θ * 2 ) = 0. Let us suppose that Eq. (41) is satisfied because {0} results otherwise. Then, there is a 0 < Θ 1 ≤ 1 that satisfies Eq. (39) . Since ∂f 1 /∂Θ * 1 < 0 and ∂f 1 /∂Θ * 2 < 0, there is a curve f 1 (Θ * 1 , Θ * 2 ) = 0 in the Θ * 1 -Θ * 2 space that looks like the solid line (I) (when f 1 (0, 1) ≤ 0), or (II) (when f 1 (0, 1) > 0) in Fig. 6(A) .
Because ∂f 2 /Θ * 1 > 0 and ∂f 2 /∂Θ * 2 < 0, f 2 monotonically increases in Θ * 1 on lines Θ * 2 = 0 and Θ * 1 + Θ * 2 = 1. If the isocline f 2 (Θ * 1 , Θ * 2 ) = 0 is located as the dashed line in Fig. 6(B) , it nontrivially crosses f 1 (Θ * 1 , Θ * 2 ) = 0 (one of the solid lines) to yield the {0, 1, 2} phase. Because f 2 (0, Θ * 2 ) = −1 (0 ≤ Θ * 2 ≤ 1), f 2 (1, 0) > 0 is necessary for the dashed line to exist. Then, there is a unique 0 < Θ < 1 so that f 2 (Θ, 1 − Θ) = 0. Finally, f 1 (Θ * 1 , Θ * 2 ) = 0 and f 2 (Θ * 1 , Θ * 2 ) = 0 cross when (i) f 2 (Θ 1 , 0) > 0 and (ii) f 1 (Θ, 1 − Θ) < 0 are both satisfied. Using
the condition (i) is equivalent to
that is,
which is an extension of the ordinary meanfield solution (Eq. (33)). We note that (i) implies f 2 (1, 0) > 0, which we required beforehand.
To show (ii), define
We obtain
To summarize, both {0} and {0, 1} disappear to be replaced by {0, 1, 2} as k 2 → ∞.
V. ROCK-SCISSORS-PAPER GAME
Winnerless cyclic competition among different phenotypes abounds in nature. For example, real microbial communities of Escherichia coli [19, 20] and color polymorphisms of natural lizards [48] have cyclically dominating three states. The evolutionary public-good game with volunteering (choice of not joining the game) also defines a three-state population dynamics with cyclic competition [47, 49] . A typical consequence of these dynamics is oscillatory population density with each phenotype alternatively dominant.
Let us consider a simple rock-scissors-paper game with cyclic competition (Fig. 7(A) ). The ordinary meanfield analysis yields a single phase: a neutrally stable periodic orbit in which the densities of states 0, 1, and 2 are cyclically dominant [17] . The coexistence equilibrium inside the limit cycle is stabilized on regular lattices [13, 51] and trees [40] . Consequently, convergence to {0, 1, 2} with a damped oscillation occurs on these graphs. We focus on the steady states only in the following analysis.
The dynamics in heterogeneous populations are given bẏ
which yields
Noting that Eq. (48) is independent of k, we obtain
for any k. The degree distribution does not affect the steady state.
Similarly, it is straightforward to show that degree distributions are irrelevant in the voter model ( Fig. 7(B) ) and cyclic interaction models with more than three states [41, 50] . A lesson from these examples is that heterogeneous contact rates do not influence the equilibrium population density if there is no neighbor-independent state transition. However, models with at least one neighbor-independent transitions show extinction of the epidemic threshold as k 2 / k → ∞. For example, variants of the rock-scissors-paper game supplied with spontaneous death or mutation rates [52] are expected to belong to the class analyzed in Sec. IV.
In many infectious diseases of humans, such as malaria, yellow feber, and dengue feber, transmission are mediated by other hosts such as mosquitoes. In these diseases, direct human-human or mosquito-mosquito infection is absent, and two distinct populations transmit diseases between each other. As an example, we analyze a model for malaria spreads (Anderson and May, 1991 , Ch. 14). Following the original notation, the meanfield dynamics are represented byẏ
where y and ψ are the proportions of infected humans and mosquitoes, respectively. The size of the human population is denoted by N, and Ω is the size of the female mosquito population; only female mosquitoes infect humans. In addition, a is the bite rate, b and c are infection rates, γ is the recovery rate of the humans, and µ is the death rate of the mosquitoes. Infected humans and infected mosquitoes survive simultaneously if they do.
Based on Eq. (50), this occurs when
Otherwise, the malaria is eradicated.
How does this condition change by introducing heterogeneous contact rates? In network terminology, relevant networks are bipartite, with each part representing the human population and the mosquito population. Let us denote the degree distributions of humans and mosquitoes by {p y,k } and {p ψ,k }, respectively. Real data do not suggest that {p y,k } or {p ψ,k } is scale-free. However, we treat general degree distributions because the obtained results may apply to other multi-population contagion processes.
The dynamics are given byẏ
where y k (ψ k ) is the probability that humans (mosquitoes) with degree k are infected, and
The steady state is calculated as
which is solved trivially by (Θ * y , Θ * ψ ) = (0, 0). To explore nontrivial solutions, let us eliminate Θ * ψ from Eqs. (53) and (54) to obtain
The RHS of Eq. (55) is less than 1 when Θ * y = 1. The endemic state results if
which extends Eq. (51) . Equation (56) indicates that divergence of just either k 2 y or k 2 ψ is sufficient for the epidemic threshold to disappear. Even if both moments are finite, their effects are multiplicative. This is consistent with the results for two-sex models of heterosexual HIV transmissions [29] and the bond percolation on bipartite graphs [31, 33] .
We have analyzed various models of endemic infectious diseases in populations with heterogeneous contact rates, or on complex networks. The effects of heterogeneity on epidemic thresholds are summarized in Tab. I. In many models, diverging second moments of the degree distribution extinguishes the epidemic thresholds, as reported previously for the percolation [2, 9, 10, 33] , the contact process [16, 36, 37] , the SIR model [4, 29, 30, 32] , the SIRS model [25] , and the household moodel [26] . On scale-free networks, which underly sexually transmitted diseases and computer viruses (see Sec. I for references), introduction of a tiny amount of virus to a population can cause an endemic state even in models with more complex transition rules.
However, the models with competing pathogens and mutation (Secs. III A and III B) show different behavior. The heterogeneity in contact rates equally boosts the infection strengths of competing pathogens in these models. Whether multiple strains survive or one overwhelms the others depends on network-independent mutation rates that modulate relative strengths of pathogens. The rock-scissors-paper game and the voter model stipulate another class of models in which heterogeneity does not affect the equilibrium population density at all. For the heterogeneity to take effects, there must be at least one neighbor-independent transition rate.
For a fixed graph or a contact rate distribution, one can create new models that are not convered by this paper. Our analysis has not been exhaustive. However, we consider that what essentially matters is gross arrangements of contagion pathways. If there is at least one transmission route independent of the neighbors' states, which is not the case for the rocksrissors-paper game, the equilibrium population density will depend on degree distributions as does the contact process. A care must be paid to cases of competing pathogens with mutation. We hope that our results give prescription for understanding complex-network consequences of other models.
A last note is on the stability of solutions. For the contact process on complex networks, the coexistence phase is stable if it exists [5, 35] . Similarly, the stability of the simplest nontrivial phase is assured for models analyzed in this paper. However, stability analysis of other phases seems mathematically difficult. In addition, the stability of coexistence solutions of the rock-scissors-paper game shows network dependence [28] . These topics are warranted for future work. 