Global stability of a network-based SIS epidemic model with a general nonlinear incidence rate.

In this paper, we develop and analyze an SIS epidemic model with a general nonlinear incidence rate, as well as degree-dependent birth and natural death, on heterogeneous networks. We analytically derive the epidemic threshold R0 which completely governs the disease dynamics: when R0 < 1, the disease-free equilibrium is globally asymptotically stable, i.e., the disease will die out; when R0 > 1, the disease is permanent. It is interesting that the threshold value R0 bears no relation to the functional form of the nonlinear incidence rate and degree-dependent birth. Furthermore, by applying an iteration scheme and the theory of cooperative system respectively, we obtain sufficient conditions under which the endemic equilibrium is globally asymptotically stable. Our results improve and generalize some known results. To illustrate the theoretical results, the corresponding numerical simulations are also given.


(Communicated by Shigui Ruan)
Abstract. In this paper, we develop and analyze an SIS epidemic model with a general nonlinear incidence rate, as well as degree-dependent birth and natural death, on heterogeneous networks. We analytically derive the epidemic threshold R 0 which completely governs the disease dynamics: when R 0 < 1, the disease-free equilibrium is globally asymptotically stable, i.e., the disease will die out; when R 0 > 1, the disease is permanent. It is interesting that the threshold value R 0 bears no relation to the functional form of the nonlinear incidence rate and degree-dependent birth. Furthermore, by applying an iteration scheme and the theory of cooperative system respectively, we obtain sufficient conditions under which the endemic equilibrium is globally asymptotically stable. Our results improve and generalize some known results. To illustrate the theoretical results, the corresponding numerical simulations are also given.

1.
Introduction. In the traditional epidemiology, mathematical models become important tools in understanding epidemic dynamics and making strategies to control disease (see [1,4,9,15,18] and the references therein). All the above researches are mainly based on the homogeneous mixing approximation, that is, each susceptible individual has the same rate of disease-causing contacts. However, in reality, each individual has limited contact with those who can spread disease. For better handling the effects of contact heterogeneity, the disease transmission should be modeled over complex heterogenous networks [11,27]. On networks, nodes stand for individuals and an edge connecting two nodes describes the interaction between individuals, in which the infection may spread. In this way, the node with more edges has a higher possibility of being infected. Recently, considerable concern has arisen over the study of epidemic models on complex heterogeneous networks (see [2,3,5,7,8], [11]− [14], [16,17], [19]− [27]) To deal with the heterogeneity of contact patterns, one needs to consider the difference of node degree. For epidemic spreading of SIS process, each node in the network can be either susceptible (S) or infected (I) at any time. Since a real network is composed of a finite number of nodes, we Let N be the total nodes. Then N = S + I. We classify all the nodes into groups based on the numbers of edges emanating from a node. That is, each node in the k−th group has the same edges (i.e., the same degree), say k, for k = 1, 2, · · · , n. Here n is the maximum node degree of the finite-size network (1 ≤ n ≤ N ) [7]. We let S k (t) and I k (t) be the densities of susceptible and infected nodes with a given degree k at time t, respectively, and let N k (t) be the number of nodes with degree k at time t, that is, N k (t) = S k (t) + I k (t). Since different groups may have different fertility levels, we let the degree-dependent parameter b k > 0 be the number of newly born nodes with degree k per unit time. And we assume that each newly born node is susceptible. Based on the above assumptions and the models in [17,26], we have the following dynamics model: where λ > 0 is the transmission rate; the natural deaths are proportional to the densities of nodes with death rate µ > 0; γ > 0 is the recovery rate of the infected nodes. According to [5,7,22], the probability Θ(t) that any given edge emanating from a node with degree k connects to an infected node can be written as where the factor 1/i accounts for the probability that one of the infected neighbors of a node, with degree i, will contact this node at the present time step. P (i|k) is the conditional probability that a node with degreek is connected to a node with degreei. ϕ(k) represents the infectivity of a node with degree k, i.e., ϕ(k) denotes the average number of occupied edges from which a node with degree k can transmit the disease [8]. This means that ϕ(k) ≤ k. It should be noted that various types of the infectivity ϕ(k) have been studied, such as ϕ(k) = k [11,12,14,17,20,27]; ϕ(k) = A [21]; ϕ(k) = min{A, σk}, 0 < σ ≤ 1 [3]; ϕ(k) = k σ [2] and ϕ(k) = ak σ /(1 + νk σ ), a > 0, ν ≥ 0 [22]. That is, the function ϕ(i) in (2) can take any of the above forms according to the degree of real networks. For simplicity, in this paper, we assume that the connectivity of nodes is uncorrelated, that is P (i|k) = iP (i)/ < k > [13], where < k >= n k=1 kP (k) is the average degree of the network; P (k) is the probability that a randomly chosen node has degree k, which is also named as the degree distribution. For convenience, we define < u(k) >:= n k=1 u(k)P (k), where u(k) is a function of the variable k. Model (1) is used to describe a kind of diseases, such as Tuberculosis, which is persistent and can last for a individual's lifetime. Some special cases of model (1) were studied, such as b k = µ = 0, γ = 1, ϕ(k) = k and N k (0) = 1 in [17,20]; b k = µ = 0, γ = 1, ϕ(k) = min{A, σk} and N k (0) = 1 in [3]; and b k = µ = 0, γ = 1 and N k (0) = 1 in [22]. And if b k = µN k (t) (i.e., deaths are balanced by births) and N k (0) = 1, then model (1) becomes model (2) with λ(k) = λk in [26].
However, there are some reasons for using nonlinear incidence rates such as saturating and nearly bilinear in the process of disease modelling. For example, to consider the instinctive reaction of people, Zhang and Sun [23] introduced a nonlinear incidence rate λkS k (t)(1 − αΘ(t))Θ(t) into their network-based SIS epidemic model, where the positive parameter α is called 'fear factor'. To describe the psychological effect of certain diseases spread in a contact network, Li [11] proposed a nonmonotone incidence rate λkS k (t)Θ(t)/(1 + αΘ 2 (t)). Enlightened by [1,9], we will study a general nonlinear incidence rate given by λkS k (t)Θ(t)/f (Θ(t)). Obviously it is a bilinear incidence rate as f (Θ(t)) = 1. We assume that the function f (Θ(t)) satisfies (H 1 ): Some of the specific forms of f (Θ(t)) appearing in the literature and satisfying (H 1 ) and (H 2 ) are f (Θ(t)) = 1 + αΘ 2 (t) with α > 0 [11] andf (Θ(t)) = 1/(1 − αΘ(t)) with 0 < α < 1 [23]. One can note from (H 1 ) that, for Θ(t) small enough, the bilinear term dominates. What's more, it follows from (H 2 ) that 1/f (Θ) is increasing when Θ(t) is small and decreasing when Θ(t) is large (for example f (Θ) = 1 + αΘ(t), α ≥ 0). In this case, the function 1/f (Θ) can be used to interpret the psychological or inhibition effect from the behavior change of the susceptible individuals. This is because the number of contacts with the infected individuals (or the force of infection) may decrease, when the probability that one may contact with infected individuals increases (i.e., Θ(t) is becoming large). As pointed out by Zhang and Sun [23], people will consciously reduce the number of contacts with others during the period of the diseases prevalence. The larger the probability that one may contact with infected individuals is, the more cautious the people will become, and the more number of contacts will be reduced in everyday life.
Motivated by the above consideration, we propose the following network-based SIS epidemic model with a general nonlinear incidence rate, as well as degreedependent birth and natural death: where the parameters and variables are the same as aforementioned. Note that when f (Θ(t)) = 1, the nonlinear incidence rate becomes the bilinear one and model (3) can be simplified to model (1). When f (Θ(t)) = 1/(1 − αΘ(t)), 0 < α < 1, b k = µN k (t), N k (0) = 1 and ϕ(k) = k, model (3) becomes model (4) in [23]. When f (Θ(t)) = 1 + αΘ 2 (t), α > 0, b k = µ = 0, N k (0) = 1 and ϕ(k) = k, model (3) becomes model (4) in [11]. A fundamental problem in epidemiology is to study the global dynamics of epidemic spreading. The global behaviors of network-based SIS epidemic models are well studied in [20,26], but they are only based on the bilinear incidence rate. To date, there has still been relatively little research studied on network-based epidemic models with nonlinear incidence rate. In [14], Liu and Ruan introduced an SIS model with a generalized nonlinear incidence rate on scale-free networks. They derived the basic reproduction number and studied the stability of the disease-free equilibrium, but they neglected to give the stability of the endemic equilibrium. As already mentioned, Li [11], Zhang and Sun [23] analyzed the dynamics of their network-based models with different nonlinear incidence rates, respectively. In [23], Zhang and Sun obtained the globally asymptotical stability of the disease-free equilibrium and the local stability of the endemic equilibrium. Later, in [24], they further studied an SIS model with a generalized feedback mechanism on weighted networks and obtained the similar results. In [11], Li proved that if the transmission rate is greater than the threshold value, the disease is permanent; otherwise, the disease-free equilibrium is globally attractive. By the numerical simulations, all of [11,23,24] observed that the endemic equilibrium is globally asymptotically stable. However, the authors have not found a strict mathematical proof of this conclusion in the literatures, which is a very challenging issue.
The aim of this paper is to investigate the global dynamics of system (3). The rest of this paper is organized as follows. In Section 2, we reveal some properties of the solutions and obtain the epidemic threshold. In Section 3, the globally asymptotical stability of the disease-free equilibrium and the permanence of epidemic are showed. In Section 4, the globally asymptotical stability of the endemic equilibrium is discussed. In Section 5, numerical simulations are given to support the theoretical analysis, and then the paper concludes with a brief discussion in Section 6.
2. Positivity, boundedness and equilibria. Before going into details, let us simply system (3). For each k, summing the two equations in (3), it follows that . Thus, we assume that the initial value is N k (0) = S k (0) + I k (0) = b k µ =: δ k , for k = 1, 2, · · · , n, in order to have a population of constant size (i.e., S k (t)+I k (t) ≡ δ k , for all t ≥ 0, k = 1, 2, · · · , n). Then system (3) becomes the following form: Considering the uncorrelated networks, it follows from (2) that In order to investigate the global stability of system (3), we only need to study the global stability of system (4). From a practical perspective, only the case of P (k) > 0, for k = 1, 2, · · · , n, is considered, and the initial conditions for system (3) (or system (4)) satisfy: In the following lemma, some properties of the solutions are obtained.
Proof. First, it follows from (4) and (5) that Θ(t) satisfies which implies that Multiplying the above inequality by exp (µ+γ)t+λk f (Θ(s)) ds and integrating from 0 to t, we get Next, it can be verified that the function δ k − I k (t) satisfied the equation Similarly proof shows that Now, we discuss all biologically feasible equilibria of System (4). One can easily find that there exists a zero equilibrium I k = 0 (k = 1, 2, · · · , n), which is corresponding to the disease-free equilibrium of system (3). Let dI k (t) dt = 0, then it follows from system (4) that Then From the above discussion, we have the following result.
Lemma 2.2. If and only if R 0 > 1, system (4) has a unique positive equilibrium I * k , k = 1, 2, · · · , n, which is corresponding to the endemic equilibrium of system (3) and satisfies (1) Lemma 2.2 shows that the existence of the endemic equilibrium depends on the epidemic threshold R 0 , which is determined by the model parameters and network structure.
(2) It is seen that the decrease of the transmission rate λ and the increase of the recovery rate γ can deduce the decrease of R 0 . Thus it will be easier for us to control the disease.
(3) More interestingly, the epidemic threshold R 0 bears no relation to the functional form of the nonlinear incidence rate and degree-dependent birth b k . In other words, the nonlinear incidence rate and degree-dependent birth do not affect the epidemic threshold R 0 .
3. Stability of the disease-free equilibrium and permanence of the disease. In this section, the global behavior of the disease-free equilibrium and the permanence of the disease are discussed.
Theorem 3.1. If R 0 < 1, then the disease-free equilibrium of system (4) is globally asymptotically stable, i.e., the disease fades out.
Obviously, (12) can be written as where y(t) = (y 1 (t), y 2 (t), · · · , y n (t)) T , The following characteristic polynomial can be calculated by mathematical induction method: Specially, if λ = −A i , i = 1, 2, · · · , n, then Since the function L( λ) is continuous, there exits at least one root in (−A i+1 , −A i ), for i = 1, 2, · · · , n − 1. Namely, there are n − 1 negative roots in (−A n , −A 1 ). It is clear that < 0, and according to (10), Consequently, there exits a negative root in (−A 1 , 0). Up to now, it is proven that all the eigenvalues of matrix J are negative, that is, the endemic equilibrium of system (4) is locally asymptotically stable. The proof is completed.
Next, applying a novel monotone iterative technique in [5,20,26,27], we obtain sufficient conditions for the global asymptotical stability of the endemic equilibrium of system (4).
Turning back to system (4), one has 1 . Consequently, for any enough small constant 1 , there exits a T such that

. (16)
Accordingly, it follows from system (4) and (15) that 2 . So, for any enough small constant there exits a T 2 .
Continuously, i = 3, 4, · · · , for any enough small constant And for any enough small constant i ,

Remark 2.
From the conditions of Theorem 4.2, especially the condition λ > µ+γ, we can conclude that the disease has the potential to become endemic when the transmission rate λ is greater than the natural death rate µ and the recovered rate γ. Therefore, in order to control the disease, we should reduce the transmission rate λ and increase the recovered rate γ, which is in accord with Remark 1(2).
Similar to the analysis of Theorem 4.2, we have the following corollaries.
A natural and interesting question now arises about whether the same result holds for a larger parameter α. In other words, shall we ignore the conditions f (Θ) ≤ 1 and λ > µ + γ ? To solve this question, we turn to the following results from the theory of cooperative system.
is called cooperative in an open set X ⊂ R n , if ∂Fi ∂xj (x) ≥ 0 for i = j and for all x ∈ X.
Lemma 4.4. [6] Suppose that X = R n , or IntR n + , or [[p, q]]. Then the cooperative system (28) has a globally asymptotically stable equilibrium if and only if the following conditions hold in X: (a) every forward semi-obit has compact closure; and (b) there is not more than one equilibrium.
Remark 4. Corollary 2 and Theorem 4.5 show that the outstanding issues in [11,23] are partly solved. So the results we obtained improve and complement those of [11,23].
5. Numerical simulations. In this section, some numerical simulations are given to illustrate the theoretical analysis. All the simulations are based on a finite scalefree networks, where the degree distribution is P (k) = Ck −τ , 2 < τ ≤ 3 and the constant C is chosen to satisfy n k=1 P (k) = 1. Let I(t) = n k=1 P (k)I k (t) and S(t) = n k=1 P (k)S k (t) be the global average densities of the two epidemic classes. Since S k (t) + I k (t) ≡ δ k , the variables I k (t) (k = 1, 2, · · · , n) are only considered. Now we study the dynamical behaviors of system (3) with τ = 2.6. We assume ϕ(k) = ak σ 1+bk σ . In Fig.1−3 and Fig.5, we chose n = 500, a = 0.3, σ = 0.75, ν = 0.02, then β = ϕ(k) / k = 0.2299. Figure 1 displays the time series I(t) with different incidence rates. The initial value of Fig.1(a) is I(0) = 80. The parameters are chosen as: λ = 0.2, µ = 0.06, γ = 0.25 and b k = 6. Then the epidemic threshold R 0 = 0.6739 < 1. In Fig.1(b), the initial value is I(0) = 5, and the parameters are listed as follows: λ = 0.5,   µ = 0.05, γ = 0.2 and b k = 4. Then R 0 = 2.0891 > 1. It can be seen from Fig.1, regardless of the functional form of the nonlinear incidence rate, when R 0 < 1, the disease will disappear; when R 0 > 1, the epidemic disease is permanent on the network. In the following Fig.2 and Fig.3, we only show f (Θ(t)) = 1 + αΘ 2 (t) (α = 1.75) on behalf of other forms of the function f (Θ(t)). To further study the detailed outcome of system (3), we should examine the time series of those nodes with different degree. In Fig.2(a) and Fig.2(b), the initial value and the parameters are the same as those of Fig.1(a) and Fig.1(b), respectively. Figure 2 also verifies that when R 0 < 1, the disease-free equilibrium is globally asymptotically stable; when R 0 > 1 and λ > µ+γ, the number of the infected with different degree will converge to a positive constant, respectively. Figure 3 depicts the relevance I k (t) versus t with different initial values. Here we choose k = 50 on behalf of other degrees. It should be noted that the time evolution of the infected nodes with other degrees are analogous. The parameters in Fig.3(a) and Fig.3(b) are the same as those in Fig.1(a) and Fig.1(b), respectively. One can observe from Fig.3 that, no matter how many the initial values of the infected nodes are, the density function I k (t) (k = 1, 2, · · · , n) tends to 0 and approaches to a positive stationary level according to above two cases, respectively. The numerical results mentioned above coincide with our theoretical analysis.    From Theorem 4.2 and Theorem 4.5, we know that the endemic equilibrium of system (3) (for example f (Θ(t)) = 1 + αΘ 2 (t)) is globally asymptotically stable under some conditions (that is R 0 > 1, λ > µ + γ and 0 < α ≤ 1/(2β); or R 0 > 1 and 0 < α ≤ 1/β 2 ), which is shown in Fig.1(b), Fig.2(b) and Fig.3(b). However, in Fig.4, we choose n = 100, a = 0.85, σ = 0.75, ν = 0.01, α = 3, λ = 0.02, µ = 0.01, γ = 0.02. Then R 0 = 1.7184 > 1 and β = ϕ(k) / k = 0.6814. It is clear that λ < µ + γ. Through simple computation, we have α − 1/(2β) = 1.0162 > 0 and α − 1/β 2 = 0.8463 > 0. It follows from Fig.4 that the endemic equilibrium of system (3) is also globally asymptotically stable only when R 0 > 1, though the rigorous analysis does not present in this paper. In addition, it is also found in Fig.4 and Fig.2(b) that if R 0 > 1 and b k is a monotone increasing function of degree k or a degree-independent constant, the larger the degree number is, the higher the endemic level will be, although degree-dependent birth b k cannot change the epidemic threshold R 0 .
Finally, we investigate the effect of the nonlinear incidence rate on the spread of a disease. Without loss of generality, we choose the saturation incidence rate, such as f (Θ(t)) = 1 + αΘ(t), where the parameter α > 0 describes the inhibition effect of the general public toward the infectivities (the same meaning as [18]). In this case, we only study the effect of the parameter α on the disease transmission. In Fig.5, the parameters are chosen as: λ = 0.3, µ = 0.01, γ = 0.2 and b k = 4, then R 0 = 1.4922 > 1. An interesting discovery shown in Fig.5 is that, when the disease is endemic, the larger the value of parameter α is, the lower the endemic level will be, although the parameter α cannot affect the epidemic threshold R 0 . This result is consistent with that in [11,23].
6. Discussion. The purpose of this paper is to study the global dynamics of a newly proposed SIS epidemic model which incorporates a general nonlinear incidence rate, as well as degree-dependent birth and natural death, on heterogeneous networks. Some special cases of this model were studied in [3,11,17,20,22,23]. We analytically derive the expression for the epidemic threshold R 0 , which determines not only the existence of endemic equilibrium but also the global dynamics of the model. Interestingly, the epidemic threshold R 0 is not dependent on the functional form of the nonlinear incidence rate, but our simulations show that the nonlinear incidence rate does affect the epidemic dynamics.
By constructing Lyapunov function, we show that when R 0 < 1, the disease-free equilibrium of system (3) is globally asymptotically stable, i.e., the disease will die out. And when R 0 > 1, the disease will persist on the network. Furthermore, by applying an iteration scheme and the theory of cooperative system respectively, we obtain sufficient conditions which ensure the globally asymptotical stability of the endemic equilibrium of system (3) and offer partial answers to the open problems in [11,23]. We believe that the idea here can also be applied to study the global dynamics of the model in [24]. In particular, the endemic equilibrium of system (1) with bilinear incidence rate (i.e., f (Θ(t)) = 1) is globally asymptotically stable when R 0 > 1, which enriches the related result in [20,26]. Hence, our results are more general and richer.
At the same time, we have considered the effect of degree-dependent birth b k on the epidemic spreading. The simulations illustrate that when the disease is endemic and b k is a monotone increasing function of degree k or a degree-independent constant, the larger the degree number is, the higher the endemic level will be, although degree-dependent birth b k cannot change the epidemic threshold R 0 . In addition, if b k is a degree-independent constant, from the expression of Θ(t), it follows that the larger the value of b k , the lower the probability Θ(t) will be. And we can note from (H 1 ) that for Θ(t) small enough, the bilinear incidence rate dominates. Therefore, according to Corollary 1, if b k is a degree-independent constant and becomes large enough, the endemic equilibrium of system (3) is globally asymptotically stable when R 0 > 1. From Theorem 4.2 and Theorem 4.5, one can see that in order to obtain the globally asymptotical stability of the endemic equilibrium of system (3), besides the epidemic threshold R 0 > 1, some additional conditions are required. That is, if R 0 > 1, we showed that the endemic equilibrium of system (3) is globally asymptotically stable provided that (C 1 ): λ > µ + γ and f (Θ) ≤ 1 or (C 2 ): f (Θ) ≥ Θf (Θ). Obviously, if f (Θ) ≤ 1, it follows from (H 2 ) that f (Θ) ≥ f (Θ), thus (C 2 ) holds. In this case, the condition λ > µ+γ in (C 1 ) can be ignored. What's more, from the numerical simulation results, we can obtain the endemic equilibrium of system (3) is globally asymptotically stable only when R 0 > 1. However, we could not give a rigorous proof for the problem in this paper. We leave this for future work.