Nonequilibrium phase transitions in metapopulation models of infectious diseases on heterogeneous networks

We study two meta-population models of infectious diseases in heterogeneous networks. We distinguish between asymptomatic and symptomatic infections and these two go through the different courses of infection and recovery. We consider that asymptomatic infections are described by an SIS model and symptomatic infections by an SIR or SIRS model depending on the immunity upon recovery. By introducing the probability of being infected asymptomatically, we combine an SIS model for asymptomatic infections with an SIR or SIRS model for symptomatic infections to obtain the SIS-SIR and SIS-SIRS models. We use a heterogeneous mean-field theory and Monte Carlo simulations to analyze two models and find that both models undergo nonequilibrium continuous phase transitions from the endemic phase to the disease-free phase at certain critical thresholds as we vary the proportion of asymptomatic infections. It suggests that it may be possible to maintain the population in the disease-free phase by controlling the proportion of asymptomatic infections. The SIS-SIRS model shows that asymptomatic infection drives symptomatic infection and vice versa. In addition, the spreading of infections eventually ceases as the population decreases even at a fixed proportion of asymptomatic infections corresponding to the endemic phase. The results provide a theoretical basis for understanding the epidemiological facts that social distancing and reducing asymptomatic infections are important factors in optimizing quarantine measures to prevent the epidemic outbreaks of infectious diseases.

We study two meta-population models of infectious diseases in heterogeneous networks. We distinguish between asymptomatic and symptomatic infections and these two go through the different courses of infection and recovery. We consider that asymptomatic infections are described by an SIS model and symptomatic infections by an SIR or SIRS model depending on the immunity upon recovery. By introducing the probability of being infected asymptomatically, we combine an SIS model for asymptomatic infections with an SIR or SIRS model for symptomatic infections to obtain the SIS-SIR and SIS-SIRS models. We use a heterogeneous mean-field theory and Monte Carlo simulations to analyze two models and find that both models undergo nonequilibrium continuous phase transitions from the endemic phase to the disease-free phase at certain critical thresholds as we vary the proportion of asymptomatic infections. It suggests that it may be possible to maintain the population in the disease-free phase by controlling the proportion of asymptomatic infections. The SIS-SIRS model shows that asymptomatic infection drives symptomatic infection and vice versa. In addition, the spreading of infections eventually ceases as the population decreases even at a fixed proportion of asymptomatic infections corresponding to the endemic phase. The results provide a theoretical basis * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The evolution of infectious diseases and their control have been the significant focus of concern through history and, especially, in the last couple of decades. Various diseases outbreaks including the SARS epidemic of 2002-2003, the H5N1 influenza epidemic of 2005, the H1N1 influenza pandemic of 2009, the ebola outbreak of 2014, and COVID-19 pandemic of 2020 were the focus of concern for physicists and mathematicians looking for some sort of mathematical modeling in order to predict the evolution of these infectious diseases [1][2][3][4][5][6]. The basic compartmental model, an SIR model, to describe the transmission of infectious diseases are first introduced by Kermack and McKendrick in the series of seminal works [7][8][9]]. An SIR model and its variations (an SIS model and an SIRS model) are standard models in the study of infectious diseases in epidemiology [4,6]. In these compartmental models, the population is divided into compartments labeled S, I, R and so on. S, I, and R denote susceptible individuals, infected individuals, and individuals who have been infected and recovered with immunity, respectively.
Infected individuals may show severe symptoms (symptomatic infection) or no symptoms in the whole course of infection (asymptomatic infection). Recent studies of COVID-19 reported that asymptomatically infected individuals are also infectious as well as symptomatically infected individuals. Because asymptomatically infected individuals are unaware of their infections, they may do a lot of daily activities, through which they can spread diseases unconsciously in the population. Hence asymptomatic infection can act as a silent driver of the pandemic of diseases [10]. In this respect, it is of great concern to understand the role of asymptomatic infection in the spreading of diseases in a unified theoretical framework for the control and the prevention of diseases spreading.
Motivated by the epidemiological characteristics of COVID-19, we distinguish the course of infection for the asymptomatically infected from the course of infection for the symptomatically infected. We consider that the asymptomatically infected follow the dynamics of the SIS model assuming that asymptomatic infection does not evoke immunity and the symptomatically infected follow the dynamics of the SIR or SIRS model with lifelong immunity (the SIR model) or temporary immunity (the SIRS model) upon recovery. By combining two courses of infections, we propose two metapopulation models (the SIS-SIR model and the SIS-SIRS model) and investigate these models in complex networks which represent the irregular structures of connected regions such as towns, cities, and countries.
In the SIS model for asymptomatic infections, the susceptible class is further divided into S and S A classes. The S class includes susceptible individuals being never infected once, whereas the S A class includes susceptible individuals who recovered from asymptomatic infections without immunity. I class includes asymptomatically infected individuals denoted as A.
In the SIR model for symptomatic infections, the S class only includes susceptible individuals being never infected once due to the absence of reinfection processes. I class includes infected individuals with symptoms denoted as B, while R class includes recovered individuals from symptomatic infections with immunity. The SIRS model includes a new susceptible class S B in addition to the classes of the SIR model. S B class includes susceptible individuals who recovered from symptomatic infections but lost immunity. Hence the susceptible class is divided into S and S B classes.
The combined SIS-SIR and SIS-SIRS models are defined by incorporating the SIS model with the SIR and the SIRS model as follows. For the incorporation, a parameter p is introduced as the probability of a susceptible individual in the S class being infected asymptomatically given that an individual is infected. With the probability 1 − p, symptomatic infection occurs. Hence the state of an individual in the S class changes as S → A with p and S → B with 1 − p, respectively. The probability p is then the proportion of asymptomatic infections in total infections. The schematic flow diagrams of the dynamics of the SIS-SIR and the SIS-SIRS models are depicted in figures 1 and 5. We study two models separately, because the SIR and SIRS models show different stationary flows. The SIR model shows an irreversible stationary flow S → B → R, whereas the SIRS model shows a cyclic stationary flow We consider the metapopulation version of the SIS-SIR and the SIS-SIRS model in quenched complex networks, irregular structures consisting of nodes connected by links. The degree k of a node is the number of links of a node, and a degree distribution P(k) characterizes the irregularity of node connections in networks [11][12][13]. In the limit of infinite network size, the second moment of degree ⟨k 2 ⟩ is finite for random networks of Poissonian P(k), but it diverges for SFNs of P(k) ∼ k −γ with γ ⩽ 3. The diverging ⟨k 2 ⟩ means that the networks are highly irregular in node connections because there are hubs having a significant portion of total degree unlike random networks without such hubs. From a metapopulation point of view, nodes can be thought of as stores, towns, cities, and even countries as one scales up the regions represented by nodes. Then corresponding links can be roads, highways, and airline routes [11].
In the metapopulation models, we consider that individuals randomly move to connected nodes after all processes are completed in every node. Recently, Zheng et al presented the effects of the diffusion, delay, and driving factors on the stationary behavior of epidemic model on the small-size network with 100 nodes by calculating the eigenvalues of the Laplacian matrix of the network [14][15][16]. Infection processes are assumed to be frequency-dependent, which means that infection rate does not depend on the population [4]. It is assumed that infections take place among individuals in the same node only, whereas infections between nodes are forbidden.
We analytically study both models using heterogeneous mean-field (HMF) approximation [17,18] and numerically confirm the analytic results via Monte Carlo simulations in the random networks and the scale-free networks (SFNs) with γ = 2.5.
In the limit of infinite population, we find that both models undergo nonequilibrium continuous phase transitions from the endemic phase to the disease-free phase at certain critical thresholds by varying p and the critical thresholds are independent of P(k). The SIS-SIRS model shows that symptomatic infections are proportional to asymptomatic infections. In addition, for both models with finite populations, the critical threshold depends on population density ρ and shifts to the inside of the endemic phase as ρ decreases, which means that the lower the population density is, the weaker the infection becomes.
The outline of the paper is as follows. In section 2, we introduce the SIS-SIR model and present the analytic results obtained from the HMF theory. The simulation results for the SIS-SIR model are presented in section 3. In section 4, we introduce the SIS-SIRS model and present the analytic results obtained from the HMF theory. The simulation results for the SIS-SIRS model are presented in section 5. The dependence of the critical threshold on population density is discussed in section 6 and we conclude with a discussion in section 7.

SIS-SIR model
The metapopulation SIS-SIR model is defined in quenched undirected networks with a degree distribution P(k) as follows.
An undirected network of size V is represented by the V × V adjacency matrix A whose element a ij = a ji is 1 if two nodes i and j are connected by a link, and otherwise a ij = a ji = 0. Conventionally self-connections are not allowed, so a ii = 0. In quenched networks, the connections between nodes do not vary in time, so the matrix A is invariant in time. The degree k i of a node i is given as k i = ∑ j a ij . The probability T ji of random hopping from node i to j solely depends on the degree of the departure node i: T ji = 1/k i if a ij = 1 and T ji = 0 if a ij = 0 As defined in the previous section, the population is divided into 5 classes in the SIS-SIR model. In what follows, S, S A , A, B, and R denote either the symbol of individuals of each class or the number of individuals of each class depending on the context. We consider the closed population of fixed size N, i.e. S + S A + I + R = N, where I = A + B is the total density of infected individuals.
For the metapopulation dynamics, we adopt the discrete-time two-step updating scheme presented in the study of the metapopulation SIS model by Colizza et al [19]. So all rates of reactions are probabilities which are shown in figure 1. We assume the homogeneous mixing of individuals occupying the same node, and all infections occur among individuals on the same node only.
We consider the frequency-dependent disease transmission (standard incidence), which reflects the situation where on average an individual makes a constant number (λ) of contacts to transmit infection in unit time [4,19]. Then in a node i, because the contact rate is λ, the probability that a random contact by an infective with a susceptible is S i /N i , and the number of infectives is I i , the rate of new infections is In each node, susceptible S's are infected asymptomatically with rate pλ and symptomatically with (1 − p)λ, respectively. Asymptomatically infected A recovers and becomes susceptible S A with rate µ, and S A is infected again with rate λ. Hence after S is asymptomatically infected, the SIS dynamics characterizes a circulation between A and S A , S A ⇆ A. Symptomatically infected B recovers with rate r, and R is not infected again due to the immunity upon recovery. In this work, we will assume that the duration of stay in each compartment is exponentially distributed so that our models will be systems of differential equations.
All reactions are depicted in figure 1. After completing reactions in all nodes, every individual randomly moves to one of the connected nodes with rate D. Then time t is increased by τ . Diffusion rate may depend on classes in general, but we numerically confirmed that classdependent diffusion rates give the same results. Hence, for simplicity, we consider the same diffusion rate D for all classes.
We consider the set Ω k of size V k consisting of nodes with the same degree k and define the densities of all classes on the set Ω k such as S k = ∑ i∈Ω k S i /V k and similarly to the others. The population of a node with degree k is given by Then the network-averaged density S is given as S = ∑ k S k P(k) and similarly to the others. In HMF approximations, we neglect the statistical fluctuations of S k and the others, for instance, S j = S k , ∀j ∈ Ω k . By setting up the discrete-time equations governing the evolutions of the densities of a node with degree k ∈ [k min , k max ], where k min and k max are the minimum and the maximum degree of a given network of size V, and taking the continuous-time limit of τ → 0, we obtain the differential equation for the density S k (t) as where the reaction kernel Θ k is defined as Θ k = I k S k /ρ k and represents the infection probability of S k by I k according to the frequency-dependent transmission. S and Θ are the average of S k and Θ k over P(k). The second term in the LHS of equation (1) denotes the outgoing individuals from the node of degree k to k nodes connected after the infection process so that it is proportional to (S k − λΘ k ) because the number of the susceptibles at the node of degree k is reduced from S k by the size of the infection λΘ k . Similarly, the third term denotes the incoming individuals from the nodes of all kinds of degree k ′ after the infection process. The hopping probability from a node of degree k ′ to a node of degree k is denoted as T kk ′ and given by T kk ′ = 1/k ′ for random hopping. The P(k ′ |k) denotes the conditional probability of a node with degree k being connected with a node with degree k ′ . For uncorrelated networks, [11][12][13], where ⟨k⟩ is the average degree defined as ⟨k⟩ = ∑ k kP(k). Similarly, with Θ Ak = I k S Ak /ρ k , we obtain the equation for S Ak as The equations for A k , B k , and R k are written as By averaging the above equations over P(k), we obtain the equations for network-averaged densities and find that the incoming flow of individuals is the same as the outgoing flow of individuals so that all diffusion terms are canceled out. The diffusion process does not affect the stationary behavior of the metapopulation model in the thermodynamic limit and does not play any role in the phase transition between two homogeneous phases, the disease-free phase and the endemic phase. However, the diffusion process may affect the stationary behavior of the epidemic model in the small-size network under some specific control parameters. Zheng et al found periodic outbreaks via Turing instabilities caused by the diffusion, delay, and driving factors [14][15][16]. The equations for network-averaged densities are written as The reaction kernels Θ and Θ A can be expressed in terms of S, S A , A, and B. by using the fact that when individuals perform random walks without any reactions, the stationary density of total individuals in a node with degree k is given by ρ k = ρk/⟨k⟩, where ρ is the total density of walkers [27]. In the SIS-SIR model, the stationary density of each class in a node with degree k is given by the product of k/⟨k⟩ and the network-averaged density of each class so that S k = kS/⟨k⟩, A k = kA/⟨k⟩ and similarly to the others. Then the reaction kernels in the steady state are written as

Steady-state solutions
To find the stationary density of each class, we set LHS of equations (6)- (10) to zero and substitute equations (11) and (12) for the reaction kernels. First, we consider B. From equation (10), we find B = 0, which yields from equations (11) and (12) Then equations (6)-(9) become Equations (14) and (17) give the same result, SA = 0. Because we know that the stationary solution for S is 0 from figure 1, we can have the endemic phase with A > 0 or the diseasefree phase with A = 0. Then equations (15) and (16) become equivalent and yield A = 0 or Next, for R, we combine equations (9) and (10) to obtain By integrating equation (18), we obtain the final size relation We consider the initial situation that one individual is infected symptomatically in the population consisting of all individuals belonging to the S class and spreads a disease in the population. Hence the initial conditions are S(0) = ρ, B(0) = 1/(ρV), and R(0) = 0. With these initial conditions and the conservation of ρ = S + S A + I + R, we derive the stationary-state solutions from equation (19).

The endemic phase.
Because Thus we find stationary-state solutions for S A , A, R, S, and B in the endemic phase as  (19), we find S A = ρ − R = pρ. Thus we find stationary-state solutions for S A , R, S, A, and B in the absorbing phase as The behaviors of the SIS-SIR model in the absorbing phase is quite similar to the ordinary SIR model [4].
The endemic phase exists only if µ/λ < p, otherwise only absorbing disease-free phase of A = 0 is possible. Therefore nonequilibrium continuous phase transition occurs at the critical threshold p c = µ/λ. This type of transitions is known as the absorbing phase transition (APT) in the literature of nonequilibrium critical phenomena [28,29]. This result agrees with the result of the ordinary SIS model of p = 1 in which the endemic phase exists for µ/λ < 1 [4,19].
We observe the absorbing phase transition in the SIS-SIR model with the critical value p c = µ/λ and For p < p c , only the absorbing phase of A = B = 0 is possible and thus the SIS-SIR model is reduced to the ordinary SIR model. On the other hand, if p is fixed and ψ = µ/λ is chosen as an external parameter, the APT is observed at ψ c = p and the phase for ψ < p is endemic.
We compare the phases of the SIS-SIR model with those of the SIS model (p = 1) and the SIR model (p = 0). The basic reproduction number for the SIS model (p = 1) is R 0 = λ/µ. If R 0 < 1 the population is in the disease-free phase, while if R 0 > 1 in the endemic phase. For the SIR model (p = 0), the disease will die out eventually (I(∞) = 0) and the population will be in the disease-free phase for all values of the parameters λ and r. By combining these two models via the probability p of being asymptomatically infected, we find the phase diagram of the SIS-SIR model as follows. For λ/µ < 1, the population is in the disease-free phase for all p and r. For λ/µ > 1, if λ/µ < 1/p it is in the disease-free phase, while if λ/µ > 1/p in the endemic phase for all r.

Simulation results of SIS-SIR model
To confirm the results of the HMF analysis, we perform Monte Carlo simulations in a random network of Poissonian P(k) and an SFN of P(k) ∼ k −γ with γ = 2.5. The average degree for the total links K and the size of network V is given by ⟨k⟩ = 2K/V, and we set ⟨k⟩ = 8 and V = 2 × 10 4 . The minimum and maximum degree, k min and k max , are set to k min = 4 for all networks, and k max = 50 for the random network and k max = 100 for the SFN.
The construction of the networks follows the algorithm of [24,30]. The number of nodes with degree k, V(k), is deterministically calculated by denotes the integer part of a real number x. To assign the degree k to V(k) nodes, V(k) number of nodes are randomly selected among nodes that are not yet assigned degrees. Then the degree k is assigned to the selected V(k) nodes. By repeating these steps for all integer k values (k ∈ [k min , k max ]), every node is assigned with a degree and then the degree sequence Finally, to construct a quenched undirected network, two nodes are randomly selected, which satisfy the condition that the two nodes are not connected and their current degrees do not exceed their assigned degrees. If the conditions are satisfied, then the two nodes are connected. Otherwise randomly select another pair of nodes satisfying the conditions. This connecting processes are repeated until all nodes are connected to each other according to the degree sequence K V .
Because the algorithm is deterministic, all quenched networks generated in this way have the same K V . Only difference is the connections between nodes. Random rewiring V nodes merely changes the indices of nodes linked to one node and doses not affect P(k). Hence networks with the same K V give the same network-average of an observable X as ⟨X⟩ = ∑ k XP(k). So we perform simulations in a single network.
We consider the invasion of a disease in the population consisting of all individuals belonging to the S class and set the total density ρ = 100. As initial conditions, we select one node randomly, and select one individual in the selected node. Then the state of the selected individual is changed into B.
The metapopulation SIS-SIR model is simulated according to the two-step parallel updating processes. For reactions in a node i, each S individual is infected with probability [19], where I i = A i + B i , and then becomes A with probability p and B with (1 − p) respectively. Each S A individual is also infected with probability (1 − [1 − λ/N i ] Ii ) and becomes A with unit probability. A becomes S A with probability µ, while B becomes R with probability r (figure 1). After completing all reactions in all nodes, every individual randomly moves to one of the connected nodes with probability D. Then time t increases by one unit (τ = 1).
In the present simulations, we set r = 0.1 and D = 1 for simplicity. We confirm that the simulation results in the steady state are not significantly changed for other values of D. For several values of ψ = µ/λ, we estimate the critical p c using the steady-state density for A. Then we compare the results with the HMF prediction, p HMF c = µ/λ. First, we present the simulation results for λ = 0.6 and µ = 0.3 in the random network and the SFN of γ = 2.5, for which the HMF theory predicts p HMF c = 1/2. Densities of all classes are measured up to 2 × 10 3 Monte Carlo time steps (t), and averaged over 20 samples starting from the independent initial-conditions. Figure 2 shows the plot of densities against t in the SFN for p = 0.4 and 0.56. As shown, the density A(t) saturates to a steady-state value for p = 0.56 corresponding to the endemic phase, but decreases to zero for p = 0.4 corresponding to the absorbing phase. So the p c locates in between. In both phases, B(t) increases at the beginning, but decreases to zero after reaching a maximum, which is the typical behavior of the density of infected individuals in the ordinary SIR model [4]. As expected by the HMF theory, the density S(t) also decays to zero in both phases.
To locate p c for λ = 0.6 and µ = 0.3 in both networks, we measure the steady-state density of A by increasing p from 0.51 to 0.58. Since A is expected to linearly scale with p as A ∼ (p − p c ) in equation (20), we estimate p c as the value p * yielding unit slope of the least-square fitting line in double-logarithmic plot of A against δ(= p − p * ). When one varies p * either larger or smaller than a true p c , A(δ) decreases nonlinearly in the double-logarithmic plot as δ decreases. The results are shown in figure 3 and p c is estimated as p c = 0.50(1), where the number in the   Figure 4 shows the phase diagram, the plot of the estimates of p c against ψ. As shown, the estimates of p c are independent of P(k), and agree well with p HMF c = µ/λ. The endemic and the absorbing phases correspond to the regions of p > ψ and p < ψ, respectively.
The SIS-SIR model shows that the endemic phase of asymptomatic infection can exist, whereas symptomatic infection always ceases in the steady state. In the next section, we show from the analysis of the SIS-SIRS model that asymptomatic infection causes symptomatic infection and vice versa.

SIS-SIRS model
The metapopulation SIS-SIRS model is defined by adding another process R → S B to the dynamics of the SIS-SIR model. The schematic flow diagram of the SIS-SIRS dynamics is shown in figure 5. The reaction R → S B occurs with probability h and a new susceptible class S B is introduced. The S B class includes susceptible individuals who recovered from symptomatic infections but lost immunity due to the temporary immunity upon recovery. We suppose that S B can be repeatedly infected with symptoms because they have already experienced symptomatic infections unlike S A . Of course, S A (S B ) can be symptomatically (asymptomatically) infected, but we neglect this possibility in the present study.
The dynamics of the SIS-SIRS model also follow the discrete-time two-step parallel updating scheme, and the disease transmission is frequency-dependent as for the SIS-SIR dynamics. Then the HMF analysis of the SIS-SIRS model is straightforward in the quenched networks defined in section 3.

HMF analysis
Similarly to the HMF analysis of the SIS-SIR model, we confirm that diffusion terms are canceled out after taking network average and thus irrelevant to the network-averaged densities. Thus we directly write down the differential equations for the network-averaged densities as follows.
The new reaction kernel Θ B is defined as The last equality holds in the steady state.

Steady-state solutions
In the steady state, we obtain the following equations for S, S A , and S B from equations (22)-(24) Because the stationary solution for S is 0 from figure 5, equation (29) is satisfied automatically. Thus equations (26) and (27) become identical in the steady state because Θ = 0 from equation (22) and λΘ B = hR from equation (24). We have an equation rB = hR from equations (26) and (27). For more informations, we add equations (26) and (27) using equations (22) and (24) to obtain the equation We again consider the invasion of a disease in the population of S's so that the initial conditions are B(0) = 1/(ρV), R(0) = S B (0) = 0, and S(0) = ρ. By integrating equation (32) together with S(∞) = 0, we obtain the final size relation neglecting 1/(ρV) term in the limit V → ∞.
From equation (33) and the conservation of the total density ρ = A + B + R + S + S A + S B , we obtain For A, we substitute S A obtained from equation (30) into equation (34) and find Similarly, by using equations (31), (33), and rB = hR, we obtain Next, we analyze equations (35) and (36). The density of infected individuals I (= A + B) should simultaneously satisfy these two equations, which can be rearranged for I as follows By equating equations (37) and (38), we obtain B as a function of A Hence B is proportional to A, and both A and B are positive in the endemic phase unlike the SIS-SIR model in which B is always zero in the steady state.
To find a critical threshold, we consider the scaling regime nearby p c , where A decreases to zero in power of δ(= p c − p). In this scaling regime, we neglect A in the denominator in equation (39) and approximate B as B ≈ µ(1 − p)A/rp. Then, we obtain A + B in the limit A → 0 as By equating equation (37) to equation (40), we obtain A nearby p c as We rearrange equation (41) for both λ > r and µ > r. The critical threshold p c can be either p c > 1 or 0 < p c < 1 depending on a ratio µ/λ. The condition for p c > 1 gives the condition of λ > µ. In this case, the absorbing phase is not possible because two conditions, p ⩽ 1 and p c > 1, means that A and B are always positive. The system is always in the endemic state for any p ∈ (0, 1) due to p c > 1, and no absorbing phase ransition takes place.
If µ > λ, p c given in equation (43) becomes less than 1 and the SIS-SIRS model undergoes the absorbing phase transition at the p c as we vary the external parameter p. When p > p c , the system is in the disease-free phase with A = B = 0. As we decrease p smaller than p c , the system starts to show the endemic behavior. In the endemic phase, A and B are proportional to each other and linearly go to zero with δ(= p c − p) as p → p − c . Equations (39) and (42) show that the SIS-SIRS model becomes the SIS model of A class for p = 1 and the SIRS model of B class for p = 0 as expected.
From an epidemiological point of view, the occurrence of the phase transition with changes in p means that it can be possible to keep the population in the disease-free phase by controlling the proportion of asymptomatic infections. In addition, asymptomatic infections drive symptomatic infections and vice versa as shown by equation (39), which supports the empirical evidence that asymptomatic infection can be acting as a silent driver of the global pandemic of COVID-19 [10].

Simulation results of SIS-SIRS model
By using the same networks of V = 2 × 10 4 as in section 3, we present the simulation results of the SIS-SIRS model for ρ = 200, r = 0.2, and h = 0.1. In order to observe the absorbing phase transition at p c ∈ (0, 1), equations (42) and (44) restrict the values of λ and µ for the given h and r values. Equation (44) requires λ > r, µ > r, and λ < µ for 0 < p c < 1, and the First, we present the simulation results in the SFN of γ = 2.5 for λ = 0.25 and µ = 0.45, which gives p HMF c = 0.36. We carry out simulations starting from the same initial conditions used in section 3, and measure all densities up to 2 × 10 4 time steps depending on p and use more than 10 2 samples. Figure 6 shows A and B densities for p = 0.42 in the panel (a) and all densities for p = 0.3 in the panel (b). In figure 6(b), we take the logarithm of y-axis to show all densities clearly. As expected from figure 6, p = 0.42 and 0.3 correspond to the absorbing phase and the endemic phase, respectively. Thus p c locates between 0.3 and 0.42.
To locate p c , we measure the steady-state densities by increasing p by 0.02 from 0.26 to 0.32, and estimate p c using the same method described in section 3. Unlike the SIS-SIR model, the equation (42) for the steady-state of A shows the prefactor depending on p. The p-dependent prefactor may play a strong correction to the linear scaling of A ∼ δ. For this reason, we use the scaled density A/ϕ with ϕ = p/(µ(1 − p) + rp). However, the p-dependence of the prefactor of B is not strong, which can be confirmed by substituting equation (42) into the approximated expression of B. Figure 7 shows the plot of A/ϕ and B against δ with the estimate of p c = 0.354(4).  Similarly, for the same p values ranging from 0.26 to 0.32, we estimate p c in the random network for λ = 0.25 and µ = 0.45. Figure 8 shows the plot of A/ϕ and B against δ with the estimate of p c = 0.357 (4). As a result, the estimates of p c for the SFN and the random networks agree well with p HMF c = 0.36, which confirms the independence of p c on P(k) distributions. In this way, we estimate p c for four values of ψ = µ(λ−r) λ(µ−r) with µ = 0.375, 0.4, 0.425, 0.45, and λ = 0.7 − µ for both networks. Figure 9 shows the resultant phase diagram. As shown, the estimates of p c are independent of degree distribution P(k), and agree well with the prediction p HMF c = ψ from the HMF theory. The endemic and the absorbing phase correspond to the regions of p < ψ and p > ψ, respectively.

Critical thresholds for low populations
In the HMF analysis, we assume the infinite limit of network size V for a fixed ρ, which means the infinite limit of population N with ρ = N/V fixed. The limit of N → ∞ can be also achieved from the limit of ρ → ∞ with a finite V, in which the HMF theory is valid as well.
In simulations, we always deal with finite V and ρ. Since the limit of V → ∞ cannot be achieved practically due to the limitations of computing resources, it is practically better to take the limit of ρ → ∞ for a finite V. In the simulations for both models, we use V = 2 × 10 4 with ρ = 100 for the SIS-SIR model and ρ = 200 for the SIS-SIRS model, and confirm that the simulation results agree well with the predictions of the HMF analysis. Thus the chosen ρ values are high enough to realize the limit of ρ → ∞.
At this point, it is natural to question whether the simulation results for low ρ agree with the theoretical predictions. In this section, we investigate the effects of the finite population not large enough to achieve the limit of infinite population.
We perform simulations by increasing ρ from 20 to 200 in the random networks of V = 2 × 10 4 for the SIS-SIR model with λ = 0.6 and µ = 0.3, and for the SIS-SIRS model with λ = 0.25 and µ = 0.45. Simulation results show that the HMF analysis is valid for ρ ⩾ 100 for the SIS-SIR model and ρ ⩾ 150 approximately for the SIS-SIRS model ( figure 10).
For the fixed V, on the other hand, when ρ is not high enough, the critical threshold depends on ρ and approaches to p HMF c as ρ increases as shown in figure 10. The ρ-dependence of p c results from the fact that the initial population of S class, S(0) = N, is not large enough to sustain infections for a given p corresponding to the endemic phase, i.e. p > p HMF c for the SIS-SIR and p < p HMF c for the SIS-SIRS model. This argument can be confirmed by the following simulation results. Figure 11 shows   These results mean that infection in the endemic phase becomes weaker as ρ decreases. Eventually, the spreading ceases at sufficiently low ρ, which leads to the shift of p c to the inside of the endemic phase. For the SIS-SIRS model, for instance at p = 0.3, A for ρ = 20 is zero, whereas A for ρ = 200 is non-zero ( figure 11). On the other hand, the equations of A, equations (20) and (42), also show that A decreases with ρ in both models, but it does not mean the change of p c . The ρ-dependence of p c results from the effects of small populations, which are not captured by the present HMF analysis.
One of the well-known quarantine measures to slow down the spreading of infectious diseases is the social distancing or lowering the number of people gathering in places, which corresponds to lowering ρ for a fixed V in the present study.

Discussion
We present two metapopulation models, the SIS-SIR and SIS-SIRS models, to investigate the effects of asymptomatic infections on the spreading of infectious diseases such as the SARS, the H5N1 influenza, the H1N1 influenza, the ebola virus, COVID-19 and so on. For our two models, we consider the frequency-dependent disease transmission and study the invasion of an infectious disease causing asymptomatic infections into the population consisting of individuals being never infected once in random networks and SFNs by using the HMF theory and Monte Carlo simulations.
The results of the HMF analysis show that both models undergo nonequilibrium absorbing phase transitions from the endemic phase to the disease-free phase at certain critical thresholds by varying the rate p of asymptomatic infection, and also the critical thresholds are independent of the degree distribution P(k) of the networks. The predictions of the HMF theory are numerically confirmed via Monte Carlo simulations.
The SIS-SIRS model shows that asymptomatic infections drive symptomatic infections and vice versa, which supports the empirical evidence that asymptomatic infection can be acting as a silent driver of the ongoing pandemic of COVID-19 [10]. More importantly, the occurrence of the phase transition with changes in p suggests that it can be possible to keep the population in the disease-free phase by controlling the proportion of asymptomatic infections.
In addition, the simulation results for low population densities confirm that the lower the population density is, the weaker the infection becomes. Eventually, the spreading of infections ceases for a sufficiently low density for a fixed p value corresponding to the endemic phase, which results in the shift of the critical threshold to the inside of the endemic phase. Hence these results support the epidemiological facts that the social distancing and restricting social gatherings of people are effective quarantine measures, and also infectious diseases spread more easily in big cities with large populations.
The frequency-dependent disease transmission means that the number of contacts between people is independent of population, so the effective infection rate is the same for all nodes representing spatial regions with different populations. As a result, the degree distribution P(k) does not affect the critical threshold, which should be the same for random and SF networks. This P(k) independence of critical threshold was also reported in the study of a metapopulation SIS model in SF networks [19].
We adopt the parallel updates in simulations instead of sequential updates for the states of individuals in simulations. Hence it can be questioned whether the discrete-time simulation results differ from those of the continuous-time rate equations used in the HMF analysis. We confirm that the discrete-time equations of the SIS-SIR and the SIS-SIRS model also show the independence of diffusion after taking network average and give the same expressions of A and B in the steady state. It was also shown in the study of the metapopulation SIS model that the discrete-time and the continuous-time formulations give the same results for the frequencydependent disease transmission [19,20].
On the other hand, for the density-dependent transmission reflecting the situation that the contact rate is proportional to population [4], the two formulations were shown to give different results [20]. In epidemiology, the density-dependent transmission is usually considered in the spreading of animal or plant diseases, whereas the frequency-dependent transmission is appropriate for the spreading of human diseases [4]. Hence we are convinced that the frequencydependent transmission in the present study properly describe the spreading processes of COVID-19 and similar types of infectious diseases.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).