Backward bifurcation and stability analysis of a network-based SIS epidemic model with saturated treatment function

Physica A 527 (2019) 121407

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Backward bifurcation and stability analysis of a network-based SIS epidemic model with saturated treatment function Yi-Jie Huang, Chun-Hsien Li

∗

Department of Mathematics, National Kaohsiung Normal University, Yanchao District, Kaohsiung City 82444, Taiwan

highlights • • • •

We present a study on a network-based SIS epidemic model with a saturated treatment function. A threshold which decides the stability of the disease-free equilibrium is obtained. Backward bifurcation occurs under certain conditions based on the threshold. We provide the results from numerical simulations to validate the proposed model.

article

info

Article history: Received 1 July 2018 Received in revised form 6 April 2019 Available online 11 May 2019 MSC: 05C82 34D23 37B25 92D30 Keywords: Complex network Epidemic model Saturated treatment function Backward bifurcation Stability

a b s t r a c t In this paper, we present a study on a network-based SIS epidemic model with a saturated treatment function to characterize the saturation phenomenon of limited medical resources. In this model, we first obtain a threshold value R0 , which is the threshold condition for the stability of the disease-free equilibrium. We show that a backward bifurcation occurs under certain conditions. More precisely, the saturated treatment function leads to a such bifurcation. In this case, R0 < 1 is not sufficient to eradicate the disease from the population. Furthermore, we also study the stability of the endemic equilibrium and the corresponding stability condition is given. Numerical experiments are conducted and their results validate the theoretical results. © 2019 Elsevier B.V. All rights reserved.

1. Introduction In recent years, infectious diseases have become a serious problem and have increased the morbidity and mortality of afflicted individuals around the world. Thus, studying the transmission mechanism of various infectious diseases and developing appropriate control strategies are relevant. To better predict and capture the factors affecting the spread of infectious diseases, mathematical modeling has been widely applied to perform various types of analyses. The commonly used models in epidemiology are generally called compartmental models. In such types of models, the population is divided into compartments and individuals in the same compartment have the same epidemiological states. To study the dynamics of each compartment with respect to time, a coupled system of differential equations is formulated and studied. For various compartmental epidemic models, we refer the reader to [1,2]. ∗ Corresponding author. E-mail address: [email protected] (C.-H. Li). https://doi.org/10.1016/j.physa.2019.121407 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

It is well-known that interactions among individuals in a population can have specific effects on the transmission of infectious diseases. The concept of a complex network is therefore incorporated into epidemic models to describe the effect of contact heterogeneity. Generally, a complex network is composed of a set of nodes and edges. Each node represents an individual in the population associated with an epidemic status, and each edge between two different nodes indicates an interaction that may enable disease transmission. One of the most popular epidemic models using complex networks was proposed by Pastor-Satorras and Vespignani [3,4]. This model is a network-based SIS model and the authors showed that there exists an epidemic threshold for the transmission rate. If the transmission rate is below the aforementioned threshold, the disease will be eradicated. Otherwise, the infection spreads and becomes persistent. Subsequently, many epidemic models have been extensively studied in the literature. For example, the SIR model and SIRS model were studied in [5] and [6], respectively. To study the effects of quarantine on the spread of diseases, SIQRS epidemic models were proposed and analyzed in [7,8]. For the diseases transmitted by a vector, network-based models with infectious media were investigated in [9–11]. In addition, multi-strain models on complex networks were proposed in [12–14], and recent studies [15,16] have considered optimal control theory into to study epidemic network models. As we have seen in the literature, many efforts have been made to obtain a threshold value (e.g., the basic reproduction number) that would determine the spreading dynamics [3–8,11,15]. Specifically, most epidemic models have shown that there is always a disease-free equilibrium and an endemic equilibrium exists when the threshold value is greater than unity. In this case, the aforementioned forward bifurcation occurs. However, there are epidemic models in which the disease transmission undergoes another type of bifurcation, known as backward bifurcation, in which a stable endemic equilibrium co-exists with a stable disease-free equilibrium when the threshold is less than unity. In other words, the threshold is less than unity is insufficient to eradicate the disease. This has serious consequences for disease control. Thus, studying the mechanisms that cause the phenomenon of backward bifurcation in epidemic models is important in epidemiology. In [17], epidemiological mechanisms that can induce the phenomenon of backward bifurcation in some epidemic models are identified. A model incorporating a nonlinear incidence rate is considered in [18]; an imperfect vaccine is provided to individuals and the conditions that ensure backward bifurcation are derived. In [19], an SIS epidemic model with the standard incidence rate and saturated treatment function was studied. The saturated treatment function was used to describe the effects of limited resources for treatment on the spread of diseases. It is shown that insufficient capacity for treatment is one of the factors that can cause backward bifurcation. Recently, the SIRS model with a saturated treatment function was investigated in [20], and the authors showed that this model also exhibits backward bifurcation. An SEIR model with a saturated incidence rate and treatment is presented in [21]. The existence of backward bifurcation depending upon some model parameters is reported therein. Thus far, to the best of our knowledge, limited efforts have been made towards studying backward bifurcation in epidemic network models. Therefore, in this paper, we present our study on a network-based SIS epidemic model with a saturated treatment function. As stated earlier, this type of treatment function can be used to describe the effects of limited medical resources on the spread of diseases. We will show that incorporating this treatment function does indeed lead to the occurrence of backward bifurcation. More precisely, we will derive the necessary and sufficient condition for the occurrence of backward bifurcation. The remainder of this paper is organized as follows. In Section 2, we formulate our model and present some prerequisites. In Section 3, we present our study on the stability of the disease-free equilibrium and determine the conditions that ensure backward bifurcation. Stability and existence of the endemic equilibrium are discussed in Section 4, and in Section 5, numerical experiments are used to illustrate the theoretical results. Finally, conclusions and future work are discussed in Section 6. 2. Model formulation The network-based SIS model proposed by Pastor-Satorras and Vespignani [3,4] takes the following form:

{

Sk′ (t) = −β kSk (t)Θ (t) + γ Ik (t), Ik′ (t) = β kSk (t)Θ (t) − γ Ik (t),

(1)

k = 1, 2, . . . , n,

where Sk (t) represents the densities of susceptible nodes with node degree k at time t, Ik (t) represents the densities of infected nodes with node degree k at time t, β > 0 is the transmission rate, γ > 0 is the recovery rate of infected nodes, and n is the maximum degree of the network. The term Θ (t) is the probability that a given edge is connected to an infected node. Throughout this paper, we assume that the connectivity of nodes in the network is uncorrelated, then

Θ (t) =

n 1 ∑

⟨k ⟩

hP(h)Ih (t),

h=1

∑n

∑n

where P(h) > 0 is the probability that a node has degree h and consequently, h=1 P(h) = 1; ⟨k⟩ = h=1 hP(h) denotes the mean degree of the network. Let Nk (t) = Sk (t) + Ik (t) for all k = 1, 2, . . . , n. It follows from (1) that Nk′ (t) = 0. Thus, Nk (t) is a constant for each k. If we incorporate the normalization condition, namely, Sk (0) + Ik (0) = 1 for all k, then system (1) can be reduced to Ik′ (t) = β k(1 − Ik (t))Θ (t) − γ Ik (t),

k = 1, 2, . . . , n.

(2)

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

3

It is worth noting that, in system (1) or (2), the recovery from the infected class per unit of time is assumed to be proportional to the densities of infected nodes, that is, γ Ik . However, in some situations, this approximation may not be reasonable. For example, when a severe infectious disease takes the world by surprise, patients will emerge in large numbers for a short period of time, and subsequently, it may temporarily paralyze public services and economic productivity. Although treatment is an important and effective method to prevent and control the spread of such severe infectious disease, limited resources for treatment of patients will restrict the maximum number of recovered patients. In other words, the recovery from the infective class will reach a saturated level. Thus, we must investigate how the saturation recovery affects the spreading dynamics. In classical epidemic models, Wang [22] considered an SIR model with a treatment rate described as

{ T (I) =

rI , rI0 ,

if 0 ≤ I ≤ I0 , if I > I0 .

This means that the treatment rate is proportional to the number of infected individuals when the capacity of treatment is not exceeded; otherwise, it takes the maximal capacity. It is shown that a backward bifurcation occurs if the capacity is small and that there exists a bistable endemic equilibria if the capacity is low. In [23], the treatment function takes the form T (I) =

cI b+I

,

where c > 0 represents the maximum recovered number from the infective class, and b > 0 is half-saturation constant, which measures the efficiency of the medical resource supply in the sense that, if b is small, then the efficiency is high. Related studies regarding this type of treatment function are reported in [19,24]. Moreover, Zhang and Liu [25] introduced the following saturated treatment function T (I) =

rI 1 + αI

to describe the effect of the infected individuals being delayed for treatment. The parameter r > 0 is the cure rate and

α ≥ 0 is used to measure the extent of the effect of the infected being delayed for treatment.

To summarize, there has been relatively little research on epidemic network models with a saturated treatment function thus far. Hence, in this paper, we present our study on the following network-based SIS epidemic model with a saturated treatment function: ⎧ r ⎪ ⎨Sk′ (t) = −β kSk (t)Θ (t) + Ik (t) + γ Ik (t), 1 + α Θ (t) (3) r ⎪ Ik (t) − γ Ik (t), k = 1, 2, . . . , n, ⎩Ik′ (t) = β kSk (t)Θ (t) − 1 + α Θ (t) where γ > 0 and γ Ik is the recovery from the kth infected group without treatment; in the present paper, the parameters r and α have the same meaning as those in [25], and hence, the term rIk /(1 + α Θ ) represents the recovery from the kth infected group with treatment. Thus, the probability that a given edge connected to a node that has recovered from infected individuals with treatment is g(Θ ) =

n 1 ∑

⟨k⟩

h=1

hP(h)

r 1 + α Θ (t)

Ih (t) =

r Θ (t) 1 + α Θ (t)

.

This equation is an analogous form of the treatment function given in [25]. Besides, the parameters r and α can be used to determine the saturation levels (cf. Fig. 1). Now, we assume that the initial conditions for system (3) satisfy the normalization condition; hence, the following derivation can be made: Sk (t) + Ik (t) = 1 for all t ≥ 0 and for all k = 1, 2, . . . , n. Consequently, system (3) can be reduced to the following form: Ik′ (t) = β k(1 − Ik (t))Θ (t) −

r 1 + α Θ (t)

Ik (t) − γ Ik (t),

k = 1, 2, . . . , n.

(4)

Throughout this paper, we will focus on the dynamics of the solutions of system (4) in the following bounded region:

Ω = {(I1 , I2 , . . . , In ) : 0 ≤ Ik ≤ 1, 1 ≤ k ≤ n}. We can show that system (4) is positively invariant in Ω . The proof of this result is analogous to the proof of Lemma 3.1 in [26]; hence, we shall omit this proof in this paper. 3. Disease-free equilibrium and backward bifurcation In this section, we will discuss the stability of the disease-free equilibrium and find a condition that ensures the occurrence of backward bifurcation for system (4). To find the equilibria of system (4), we set Ik′ (t) = 0, and it follows

4

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

Fig. 1. The graphs of the treatment function g(Θ ) with r = 0.8 and various values of α .

from (4) that Ik =

β kΘ γ+

r 1 + αΘ

+ β kΘ

.

(5)

Substituting (5) into Θ , we obtain

⎛

⎞ n 2 ∑ β h P(h) ⎜ 1 ⎟ Θ⎝ − 1⎠ = 0. (6) r ⟨k⟩ + β hΘ h=1 γ + 1 + αΘ Clearly, Θ = 0 is always a solution of Eq. (6). It follows from Eq. (5) that there is always a disease-free equilibrium E0 = (0, 0, . . . , 0) for system (4). To study the stability of E0 , we introduce the following lemma to facilitate the stability analysis. Lemma 3.1 ([13]). For a real n × n matrix A = (aij ) where aij = δij σi + pi qj (pi , qj ≥ 0, i, j = 1, 2, . . . , n) and δij is the Kronecker symbol. The determinant of A is given by det(A) = σ1 σ2 · · · σn + p1 q1 σ2 σ3 · · · σn + σ1 p2 q2 σ3 · · · σn + · · · + σ1 σ2 · · · σn−1 pn qn .

(7)

If σi ̸ = 0, i = 1, 2, . . . , n, then

( det(A) =

1+

n ∑ pi qi i=1

)

σi

n ∏

σi ;

(8)

i=1

if σk = σ for all k, then

( det(A) = σ

n−1

σ+

n ∑

) pi qi

.

(9)

i=1

Now, we will discuss the stability of the disease-free equilibrium E0 . Theorem 3.1.

Define the threshold value

β⟨k2 ⟩ , (γ + r)⟨k⟩ ∑n 2 where ⟨k2 ⟩ = h=1 h P(h). If R0 < 1, then E0 is asymptotically stable; if R0 > 1, it is unstable. R0 :=

Proof. The Jacobian matrix of system (4) evaluated at an equilibrium (I1 , I2 , . . . , In ) is given by A = (aij ), where

⎧ r −r α Ii jP(j) β ijP(j) ⎪ ⎪ (1 − Ii ) − β iΘ − + , ⎨−γ + ⟨k⟩ 1 + αΘ ⟨k⟩(1 + α Θ )2 aij = β ijP(j) −r α Ii jP(j) ⎪ ⎪ ⎩ (1 − Ii ) + , ⟨k⟩ ⟨k⟩(1 + α Θ )2

if i = j, if i ̸ = j.

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

5

To obtain the stability of the disease-free equilibrium, we evaluate A at E0 :

⎧ β ijP(j) ⎪ ⎪ , ⎨−(γ + r) + ⟨k⟩ aij = β ijP(j) ⎪ ⎪ ⎩ , ⟨k ⟩

if i = j, if i ̸ = j.

In order to assess the eigenvalues of the Jacobian matrix A by Lemma 3.1, we denote the (i, j)-entries of matrix A − λI be −(γ + r + λ)δij + β ijP(j)/⟨k⟩. Set σi = −(γ + r + λ), pi = β i, and qj = jP(j)/⟨k⟩. From (9), the following derivation can be made:

( det(A − λI) = (−1)

n−1

(γ + r + λ )

n−1

−(γ + r + λ) +

n 1 ∑

⟨k⟩

) β i P(i) . 2

i=1

We can see that λ = −(γ + r) is a negative characteristic root with multiplicity n − 1. Thus, the stability of the disease-free equilibrium completely depends on the sign of the root of

−(γ + r + λ) +

n 1 ∑

⟨k⟩

β i2 P(i) = 0,

i=1

which implies that

λ = (γ + r)(R0 − 1). Accordingly, if R0 < 1, then all eigenvalues of the Jacobian matrix A are negative. Hence, E0 is asymptotically stable. Conversely, if R0 > 1, then E0 is unstable because there exists a positive characteristic root. This completes the proof. □ Theorem 3.2.

ˆ R0 = (

Denote

β⟨k2 ⟩

.

)

r

γ+

1+α

⟨k⟩

If ˆ R0 < 1, then E0 is globally asymptotically stable. Proof. Given that R0 ≤ ˆ R0 < 1, it suffices to show that E0 is globally attractive. According to system (4), we have the following:

Θ ′ (t) =

n 1 ∑

⟨k⟩

(

hP(h) β h(1 − Ih (t))Θ (t) −

h=1

r 1 + α Θ (t)

)

Ih (t) − γ Ih (t)

n

1 ∑ 2 r Θ (t) β⟨k2 ⟩ Θ (t) − β h P(h)Ih (t)Θ (t) − − γ Θ (t) ⟨k⟩ ⟨k⟩ 1 + α Θ (t) h=1 ( ) β⟨k2 ⟩ r ≤ Θ (t) − −γ ⟨k⟩ 1 + α Θ (t) ( ) β⟨k2 ⟩ r ≤ Θ (t) − −γ ⟨k⟩ 1+α ( ) r = γ+ Θ (t)(ˆ R0 − 1). 1+α

=

Given that ˆ R0 < 1, we have

( ς := γ +

r 1+α

)

(ˆ R0 − 1) < 0.

This implies that Θ ′ (t) ≤ ς Θ (t), and that one can obtain Θ (t) → 0 as t → ∞. Consequently, we can conclude that Ik (t) → 0 as t → ∞. This proves that the disease-free equilibrium of system (4) is globally attractive if ˆ R0 < 1. By applying Theorem 3.1 to this result, we can conclude that the disease-free equilibrium is globally asymptotically stable if ˆ R0 < 1. This completes the proof. □ Remark 1. From the result of Theorem 3.2, we cannot show that the disease-free equilibrium is globally asymptotically stable when R0 < 1. This motivates us to consider whether the infection can persist even if R0 < 1.

6

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

Now, we will determine whether there exists a backward bifurcation at R0 = 1. More precisely, we will derive a necessary and sufficient condition on the system parameters, such that backward bifurcation occurs. To this end, let us consider the endemic equilibria, that is, Ik > 0 for all k = 1, 2, . . . , n. Given that Ik > 0 for all k, we have Θ > 0. It follows from Eq. (6) that the endemic equilibria should satisfy the following equation: n 1 ∑

⟨k⟩

h=1

β h2 P(h) γ+

r 1 + αΘ

+ β hΘ

= 1.

Now, we want to express the above equation in terms of R0 and Θ . To develop this expression, if both denominator and numerator are multiplied by ⟨k2 ⟩/((γ + r)⟨k⟩), then we can obtain n 1 ∑

⟨k⟩

h=1

R0 h2 P(h)

⟨k2 ⟩ (γ + r)⟨k⟩

( γ+

r

)

1 + αΘ

= 1.

(10)

+ hR0 Θ

We may consider Eq. (10) as the definition of an endemic equilibrium curve in the (R0 , Θ ) positive quadrant. If we consider locally Θ as a function of R0 , then the backward bifurcation is characterized by a negative slope at (R0 , Θ ) = (1, 0), that is, the derivative at the critical value (R0 , Θ ) = (1, 0) is negative:

⏐ ∂ Θ ⏐⏐ < 0. ∂ R0 ⏐(R0 ,Θ )=(1,0) We can use this inequality to derive a necessary and sufficient condition for the system parameters such that backward bifurcation occurs. To obtain this condition, we compute the derivative ∂ Θ /∂ R0 at (R0 , Θ ) = (1, 0) by differentiating Eq. (10) implicitly. Thus, differentiating Eq. (10) implicitly with respect to R0 , we can obtain the following: n 1 ∑

⟨k⟩

h=1

I1 − I2

{

2

(

⟨k ⟩ r γ+ (γ + r)⟨k⟩ 1 + αΘ

)

+ hR0 Θ

}2 = 0,

(11)

where I1 = h2 P(h)

{

( ) } ⟨k2 ⟩ r γ+ + hR0 Θ (γ + r)⟨k⟩ 1 + αΘ

and

{

I2 = h2 P(h)R0 −

} ∂Θ ⟨k2 ⟩ rα ∂Θ . + h Θ + hR 0 (γ + r)⟨k⟩ (1 + α Θ )2 ∂ R0 ∂ R0

Now, if R0 = 1 and Θ = 0, then I1 = h2 P(h)

⟨k2 ⟩ ⟨k⟩

and

I2 = −h2 P(h)

r α⟨k2 ⟩

∂Θ

(γ + r)⟨k⟩ ∂ R0

+ h3 P(h)

Therefore, Eq. (11) becomes

( ) ⟨k2 ⟩ r α⟨k2 ⟩ ∂ Θ ∂Θ − − + h n 1 ∑ 2 ⟨k⟩ (γ + r)⟨k⟩ ∂ R0 ∂ R0 h P(h) = 0, ( 2 )2 ⟨k⟩ ⟨k ⟩ h=1 ⟨k⟩ which implies that

−

n

1 ∑

⟨k⟩

h2 P(h)

h=1

r α⟨k2 ⟩

∂Θ ∂Θ +h (γ + r)⟨k⟩ ∂ R0 ∂ R0 = 1. ( 2 )2 ⟨k ⟩ ⟨k⟩

Hence, we can obtain

⏐ ∂ Θ ⏐⏐ = ∂ R0 ⏐(R0 ,Θ )=(1,0)

1

−

rα

γ +r

+

⟨k⟩⟨k3 ⟩ ⟨k 2 ⟩2

.

∂Θ . ∂ R0

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

7

This derivative can be negative if and only if the denominator is negative. Hence, we obtain the following necessary and sufficient condition for backward bifurcation:

−

rα

+

γ +r

⟨k⟩⟨k3 ⟩ < 0, ⟨k2 ⟩2

which leads to

α>

(γ + r)⟨k⟩⟨k3 ⟩ r ⟨k2 ⟩2

.

To conclude, we have the following theorem: Theorem 3.3.

α>

System (4) has backward bifurcation at R0 = 1 if and only if

(γ + r)⟨k⟩⟨k3 ⟩ r ⟨k2 ⟩2

where ⟨k3 ⟩ =

∑n

h=1

,

(12)

h3 P(h).

Remark 2. From Theorem 3.3, we can see that the nonlinear treatment function does indeed play a key role in causing backward bifurcation. More precisely, if the delayed effect for treatment is stronger than some level, i.e., α is large enough so that condition (12) is met, then a backward bifurcation will take place. Otherwise, there is no backward bifurcation when such effect is weak. 4. Endemic equilibrium and its stability In this section, we discuss the existence of endemic equilibria and their stability. To find the endemic equilibrium, namely, Ik > 0 (k = 1, 2, . . . , n), we know from Eq. (6) that Ik should satisfy n 1 ∑

⟨k⟩

β h2 P(h) γ+

h=1

r 1 + αΘ

+ β hΘ

= 1.

(13)

Define the function f (Θ ) :=

n 1 ∑

⟨k⟩

h=1

β h2 P(h) γ+

r 1 + αΘ

+ β hΘ

.

Then we have f (0) = R0 and f (1) < 1. The following theorem states the result of the existence of an endemic equilibrium: Theorem 4.1.

The following results hold.

(1) If ˆ R0 < 1, then no endemic equilibrium exists. (2) If R0 > 1 and 0 ≤ r α ≤ β , then a unique endemic equilibrium exists. Proof. Given that ˆ R0 < 1, we can establish that f (Θ ) < ˆ R0 < 1 for all Θ ∈ [0, 1]. Thus, there is no Θ ∈ [0, 1] such that Eq. (13) holds. This implies that there exists no endemic equilibrium, and part (1) is proved. To prove part (2), we first compute d dΘ

f (Θ ) =

n 1 ∑

⟨k⟩

h=1

rα

β h2 P(h) (

(1 + α Θ )2

γ+

r 1 + αΘ

− βh + β hΘ

)2 =

n 1 ∑

⟨k⟩

β h2 P(h)

r α − β h(1 + α Θ )2 Dh2 (Θ )

h=1

,

where Dh (Θ ) = γ (1 + α Θ ) + r + β h(1 + α Θ ) > 0 for all h = 1, 2, . . . , n. If r α = 0, then df (Θ )/dΘ < 0 for Θ ∈ [0, 1]. Given that f (0) = R0 > 1 and f (1) < 1, we see that there is a unique Θ ∗ ∈ (0, 1) such that f (Θ ∗ ) = 1. From Eq. (5), we obtain a unique endemic equilibrium. We now turn to the case 0 < r α ≤ β . Define Nh (Θ ) = r α − β h(1 + α Θ )2 .

It is clear that Nh (Θ ) > Nh+1 (Θ ) for h = 1, 2, . . . , n − 1 and for each fixed Θ . Thus, one can derive the following:

( 0 = N1

−1 +

√

r α/β

α

)

> N2

(

−1 +

√

r α/β

α

)

> · · · > Nn

(

−1 +

√

α

r α/β

)

.

(14)

8

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

√

If 0 < r α ≤ β , then (√ −1 + r α/β )/α ≤ 0. Because N1 (Θ ) is a decreasing function with respect to Θ , it follows that N1 (Θ ) ≤ 0 for (−1 + r α/β )/α ≤ Θ ≤ 1. Combining this result with relation (14), we have d dΘ

n 1 ∑

f (Θ ) =

⟨k⟩

β h2 P(h)

h=1

N h (Θ ) Dh2 (Θ )

< 0.

As in the case of r α = 0, we can conclude that there is a unique Θ ∗ ∈ (0, 1) such that f (Θ ∗ ) = 1. Consequently, only one endemic equilibrium exists. This completes the proof. □ The following theorem is a stability result for the unique endemic equilibrium in Theorem 4.1 (2). If R0 > 1 and 0 ≤ r α ≤ β , the unique endemic equilibrium is asymptotically stable.

Theorem 4.2.

Proof. Let E ∗ = (I1∗ , I2∗ , . . . , In∗ ) be the unique equilibrium. The Jacobian matrix of system (4) evaluated at E ∗ is now given by A = (aij ), where

⎧ r β ijP(j) −r α Ii jP(j) ⎪ ⎨−γ + (1 − Ii∗ ) − β iΘ ∗ − + , if i = j, ∗ ⟨k⟩ 1 + αΘ ⟨k⟩(1 + α Θ ∗ )2 aij = ∗ −r α Ii jP(j) β ijP(j) ⎪ ⎩ (1 − Ii∗ ) + , if i ̸ = j, ⟨k⟩ ⟨k⟩(1 + α Θ ∗ )2 ∑n and Θ ∗ = ⟨k⟩−1 h=1 hP(h)Ih∗ . To assess the eigenvalues of the Jacobian matrix A by Lemma 3.1, we denote the entries of A − λI by (A − λI)ij = δij σi + pi qj , where ( ) r α Ii∗ r jP(j) ∗ σi = − γ + β i Θ ∗ + + λ , p = β i(1 − I ) + , and qj = . i i ∗ ∗ 2 1 + αΘ (1 + α Θ ) ⟨k⟩ According to (7), σi ̸ = 0 for each i, otherwise there exists i ̸ = j such that σi = σj , and this is impossible when the expression for σi is considered. From (8), det(A − λI) = 0 is equivalent to n 1 ∑

⟨k⟩

r α Ih∗

β h(1 − Ih∗ ) +

(1 + α Θ ∗ )2 = 1, ∆(h)

hP(h)

h=1

(15)

where

∆(h) = γ +

r

+ β hΘ ∗ + λ.

1 + αΘ ∗

Due to the identity ∗

Ih

(

γ+

β h(1 − Ih∗ ) =

)

r 1 + αΘ ∗

,

Θ∗

Eq. (15) becomes the following: n ∑

hP(h)Ih∗

h=1

T (Θ ∗ )∆(h)

=

n ∑ hP(h)I ∗ h

1

h=1

,

(16)

where T (Θ ∗ ) =

(

γ+

r 1 + αΘ ∗

+

)−1

r αΘ ∗ (1 + α Θ ∗ )2

.

Observe that Eq. (16) establishes n complex roots λν for ν = 1, 2, . . . , n. Let aν = Re(λν ) and bν = Im(λν ). Substituting each λν into Eq. (16), we obtain

( hP(h)Ih γ +

n ∑ h=1

r

∗

{( ∗ T (Θ ) γ+

1 + αΘ ∗ r

1 + αΘ ∗

+ β h Θ + aν ∗

+ β h Θ ∗ + aν

) }=

)2 + b2ν

n ∑ hP(h)I ∗ h

h=1

1

and

−bν

n ∑ h=1

hP(h)Ih∗

{( T (Θ ∗ ) γ+

r 1 + αΘ ∗

+ β h Θ + aν ∗

} = 0.

)2 +

b2 ν

(17)

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

9

Eq. (17) implies that bν = 0 for all ν . Thus, the eigenvalues λν are all real numbers. We next claim that λν < 0 for all ν . If we prove that T (Θ ) ∗

(

γ+

r 1 + αΘ ∗

+ β hΘ

∗

)

> 1 for h = 1, 2, . . . , n,

then the assertion follows from Eq. (16). For this purpose, observe that r ) ( γ+ + β hΘ ∗ r 1 + αΘ ∗ ∗ ∗ + β h Θ = . T (Θ ) γ + r r αΘ ∗ 1 + αΘ ∗

γ+

1 + αΘ ∗

+

(18)

(1 + α Θ ∗ )2

Thus, it is sufficient to show that

β hΘ ∗ >

r αΘ ∗

for h = 1, 2, . . . , n

(1 + α Θ ∗ )2

(19)

whenever 0 ≤ r α ≤ β . Clearly, inequalities (19) are fulfilled if r α = 0. For 0 < r α ≤ β , we have

β hΘ ∗ −

r αΘ ∗ (1 + α Θ ∗ )2

> β hΘ ∗ − r α Θ ∗ ≥ (β − r α )Θ ∗ ≥ 0.

Therefore, we can conclude that all eigenvalues of the Jacobian matrix A are negative. Hence, the endemic equilibrium E ∗ is asymptotically stable. This completes the proof. □ Remark 3. Now, we will discuss how the network structure affects the epidemic dynamics. We will analyze the sensitivity of ˆ R0 with respect to the parameter of the network structure by computing the elasticity of ˆ R0 . The elasticity of ˆ R0 with respect to the structure parameter ⟨k2 ⟩/⟨k⟩ is given by ˆ R0

ε ⟨k2 ⟩ = ⟨k⟩

⟨k2 ⟩ ⟨k⟩

∂ˆ R0

2 ˆ R 0 ∂ ⟨k ⟩ ⟨k⟩

= 1.

Note that the positive value of the elasticity means ˆ R0 increases with the increase in ⟨k2 ⟩/⟨k⟩. Thus, it follows from 2 Theorem 3.2 that decreasing ⟨k ⟩/⟨k⟩ could help us to control the disease, whereas increasing ⟨k2 ⟩/⟨k⟩ could increase the difficulty in controlling the disease. 5. Numerical experiments In this section, we will present a number of numerical simulations to validate our theoretical results. In the following examples, we consider scale-free ∑ networks with degree distribution P(h) = ηh−2.4 for h = 1, 2, . . . , 100, where the 100 constant η is chosen to maintain h=1 P(h) = 1. Under the aforementioned conditions, we can obtain ⟨k⟩ = 1.9607, 2 3 ⟨k ⟩ = 18.3468, and ⟨k ⟩ = 722.2562. Example 5.1. We first consider the stability of the disease-free equilibrium E0 . In this example, we choose β = 0.008, γ = 0.2, r = 0.8, and α = 6. Thus, R0 = 0.0749 < 1 and ˆ R0 = 0.2382 < 1 can be established. It follows from Theorem 3.2 that the disease-free equilibrium E0 is globally asymptotically stable. We use 10 different initial conditions to plot the trajectories of Ik (t) for k = 25, 50, 75, and 100. From Fig. 2, it can be seen that all the trajectories converge to 0. This supports the global stability of the disease-free equilibrium E0 . Example 5.2. We now move to the case for backward bifurcation. The parameters are chosen as β = 0.1, γ = 0.2, r = 0.8, and α = 10. Thus, we have R0 = 0.9357 < 1 and ˆ R0 = 3.4310 > 1. Furthermore, we can also check that

α = 10 >

γ + r ⟨k⟩⟨k3 ⟩ = 5.2589. r ⟨k2 ⟩2

It therefore follows from Theorem 3.3 that backward bifurcation occurs at R0 = 1. The bifurcation diagrams are given in Fig. 3. In the left panel, the parameter α = 5. Thus, condition (12) is not met. This implies that only forward bifurcation occurs at R0 = 1. In the right panel, we choose α = 10; consequently, the backward bifurcation can be observed. The trajectories of Ik (t) for k = 25, 50, 75, and 100 with 10 different initial conditions are depicted in Fig. 4. It can be seen that all trajectories converge to a positive constant despite R0 < 1. Remark 4. We should mention here that the condition ˆ R0 < 1 in Theorem 3.2 is quite conservative. In fact, it can be seen from the right panel of Fig. 3 that, there is a critical value of R0 , denoted as R∗0 (R∗0 ≈ 0.8935 in Example 5.2), such that if R0 < R∗0 , then system (4) has only the disease-free equilibrium. Furthermore, we may conjecture that if R0 < R∗0 , then the disease-free equilibrium is globally asymptotically stable. However, we cannot currently validate this claim because it is difficult to establish the theoretical form of R∗0 .

10

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

Fig. 2. Time evolutions of Ik (t) for k = 25, 50, 75, and 100 in Example 5.1. The 10 different initial conditions are given by Ik (0) = 0.05i for i = 1, 2, . . . , 10.

Fig. 3. Bifurcation diagram in the R0 Θ -plane. The model parameters are β = 0.1, γ = 0.2, and r = 0.8. In the left panel, we choose α = 5 and forward bifurcation occurs at R0 = 1. In the right panel, we set α = 10; hence, backward bifurcation occurs at R0 = 1.

Example 5.3. In the final example, we consider the stability of the endemic equilibrium E ∗ in Theorem 4.2. Let β = 0.1, γ = 0.3, r = 0.01, and α = 5. Then, it can be derived that R0 = 3.0185 > 1 and r α = 0.05 < β = 0.1. By Theorem 4.2, the endemic equilibrium E ∗ is asymptotically stable. The trajectories of Ik (t) for k = 25, 50, 75, and 100 are plotted in

Fig. 5. We can observe that E ∗ is indeed asymptotically stable. 6. Conclusions

In this paper, we proposed a network-based SIS epidemic model with saturation recovery from infected individuals to understand the effect of limited medical resources on disease transmission. We first obtained a threshold value R0 , which determines the local stability of the disease-free equilibrium. Furthermore, we showed that there exists a backward bifurcation at R0 = 1 provided condition (12) is met. In this case, one cannot eradicate the disease unless the value of R0 decreases such that R0 < R∗0 < 1 for some critical value R∗0 . The result discussed in Theorem 3.2 indicates that the disease does indeed die out if R0 < ˆ R0 < 1. For the existence of endemic equilibria, it can be seen from our numerical

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

11

Fig. 4. Time evolutions of Ik (t) for k = 25, 50, 75, and 100 in Example 5.2. The 10 different initial conditions are given by Ik (0) = 0.05i for i = 1, 2, . . . , 10.

Fig. 5. Time evolutions of Ik (t) for k = 25, 50, 75, and 100 in Example 5.3. The 10 different initial conditions are given by Ik (0) = 0.05i for i = 1, 2, . . . , 10.

12

Y.-J. Huang and C.-H. Li / Physica A 527 (2019) 121407

simulation that there exists a unique endemic equilibrium if R0 > 1 and there are two equilibria if R∗0 < R0 < 1 (see Fig. 3 (right)). However, we cannot currently prove this because the model is complicated and only partial results are provided in Theorems 4.1 and 4.2. We have also performed numerical experiments to validate the theoretical results. The unsolved problems for the endemic equilibria and epidemic dynamics on other types of complex networks [27,28] will be addressed in the future. Acknowledgments The authors would like to thank two anonymous referees for their valuable comments and suggestions that improved the presentation of this paper. This work was partially supported by the Ministry of Science and Technology of Taiwan under the grant MOST 106-2115-M-017-001. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, second ed., Springer, New York, 2012. M. Martcheva, An Introduction to Mathematical Epidemiology, Springer, New York, 2015. R. Pastor-Satorras, A. Vespignani, Epidemic spreading in scale-free networks, Phys. Rev. Lett. 86 (2001) 3200–3203. R. Pastor-Satorras, A. Vespignani, Epidemic dynamics and endemic states in complex networks, Phys. Rev. E 63 (2001) 066117. Y. Moreno, R. Pastor-Satorras, A. Vespignani, Epidemic outbreaks in complex heterogeneous networks, Eur. Phys. J. B 26 (2002) 521–529. C. Li, C. Tsai, S. Yang, Analysis of epidemic spreading of an SIRS model in complex heterogeneous networks, Commun. Nonlinear Sci. Numer. Simul. 19 (2014) 1042–1054. S. Huang, F. Chen, L. Chen, Global dynamics of a network-based SIQRS epidemic model with demographics and vaccination, Commun. Nonlinear Sci. Numer. Simul. 43 (2017) 296–310. T. Li, Y. Wang, Z. Guan, Spreading dynamics of a SIQRS epidemic model on scale-free networks, Commun. Nonlinear Sci. Numer. Simul. 19 (2014) 686–692. J. Juang, Y.-H. Liang, Analysis of a general SIS model with infective vectors on the complex networks, Physica A 437 (2015) 382–395. Y. Wang, Z. Jin, Z. Yang, Z. Zhang, T. Zhou, G. Sun, Global analysis of an SIS model with an infective vector on complex networks, Nonlinear Anal. RWA 13 (2012) 543–557. M. Yang, G. Chen, X. Fu, A modified SIS model with an infective medium on complex networks and its global stability, Physica A 390 (2011) 2408–2413. Q. Wu, X. Fu, M. Yang, Epidemic thresholds in a heterogenous population with competing strains, Chin. Phys. B 20 (2011) 046401. Q. Wu, M. Small, H. Liu, Superinfection behaviors on scale-free networks with competing strains, J. Nonlinear Sci. 23 (2013) 113–127. J. Yang, C.-H. Li, Dynamics of a competing two-strain SIS epidemic model on complex networks with a saturating incidence rate, J. Phys. A 49 (2016) 215601. L. Chen, J. Sun, Global stability and optimal control of an SIRS epidemic model on heterogeneous networks, Physica A 410 (2014) 196–204. L. Chen, J. Sun, Optimal vaccination and treatment of an epidemic network model, Phys. Lett. A 378 (2014) 3028–3036. A.B. Gumel, Causes of backward bifurcations in some epidemiological models, J. Math. Anal. Appl. 395 (2012) 355–365. B. Buonomo, D. Lacitignola, On the backward bifurcation of a vaccination model with nonlinear incidence, Nonlinear Anal.-Model Control 16 (2011) 30–46. J. Wei, J. Cui, Dynamics of SIS epidemic model with the standard incidence rate and saturated treatment function, Int. J. Biomath. 5 (2012) 1260003. Y. Gao, W. Zhang, D. Liu, Y. Xiao, Bifurcation analysis of an SIRS epidemic model with standard incidence reate and saturated treatment function, J. Appl. Anal. Comput. 7 (2017) 1070–1094. M.A. Khan, Y. Khan, S. Islam, Complex dynamics of an SEIR epidemic model with saturated incidence rate and treatment, Physica A 493 (2018) 210–227. W. Wang, Backward bifurcation of an epidemic model with treatment, Math. Biosci. 201 (2006) 58–71. J. Cui, X. Mu, H. Wan, Saturation recovery leads to multiple endemic equilibria and backward bifurcation, J. Theoret. Biol. 254 (2008) 275–283. J. Wang, S. Liu, B. Zhang, Y. Takeuchi, Qualitative and bifurcation analysis using an SIR model with a saturated treatment function, Math. Comput. Modelling 55 (2012) 710–722. X. Zhang, X. Liu, Backward bifurcation of an epidemic model with saturated treatment function, J. Math. Anal. Appl. 348 (2008) 433–443. A. Lajmanovich, J. Yorke, A deterministic model for gonorrhea in a nonhomogeneous population, Math. Biosci. 28 (1976) 221–236. M. Sun, H. Zhang, H. Kang, G. Zhu, Z. Fu, Epidemic spreading on adaptively weighted scale-free networks, J. Math. Biol. 74 (2017) 1263–1298. J. Zhang, J. Sun, Exponential synchronization of complex networks with continuous dynamics and Boolean mechanism, Neurocomputing 378 (2018) 146–152.