A Bistable Phenomena Induced by a Mean-Field SIS Epidemic Model on Complex Networks: A Geometric Approach

In this paper, we propose a degree-based mean-field SIS epidemic model with a saturated function on complex networks. First, we adopt an edge-compartmental approach to lower the dimensions of such a proposed system. Then we give the existence of the feasible equilibria and completely study their stability by a geometric approach. We show that the proposed system exhibits a backward bifurcation, whose stabilities are determined by signs of the tangent slopes of the epidemic curve at the associated equilibria. Our results suggest that increasing the management and the allocation of medical resources effectively mitigate the lag effect of the treatment and then reduce the risk of an outbreak. Moreover, we show that decreasing the average of a network sufficiently eradicates the disease in a region or a country.


INTRODUCTION
Mathematical modeling plays a crucial role in fighting against large scale infectious disease such as Tuberculosis, HIV, COVID-19, etc., Compartment models have been used to anticipate the progression of diseases and evaluate the effect of interventions on disease spread [1]. One of such models separates the total population into two distinct categories with respect to disease status. People who have not gotten the disease are labeled "susceptibles"; while those who have been infected by a certain disease are called "infectives". This kind of compartment model is denoted by "an SIS epidemic model" [2][3][4][5][6][7], which has been extensively used to address the dynamics of those diseases, describing an individual infected by a disease as having no immunity, thus becoming a susceptible again.
Most of the existing models assume that all the individuals are well-mixed and they have homogeneous mixing of surfaces, which implies that each individual has the same probability to contact other individuals and ignores the degree of social heterogeneity induced by age, household, spatial structures, and social spheres, etc. Generally, the social interactions of individuals generate a certain pattern based on social preferences, which contributes to transmission heterogeneity. Indeed, such factors may play a decisive role in the disease transmission and they also may help health policymakers to take more effective control measures for curbing the disease spread [7,8]. Epidemic models on complex networks incorporate such contact heterogeneity and take account for how the structures of the networks affect the disease prevalence. A popular degree-based SIS epidemic model has been built [9] and it exhibits threshold dynamics [4,5]. Since then, many factors including vector-borne [10,11], infective media [12], awareness reaction [13], and the gene diversity of pathogens [14], etc., have been incorporated into the study of the co-evolution of networks and epidemics.
In epidemiology, the basic reproduction number, R 0 , denotes the average number of secondary cases produced by a typical disease carrier during their infectious period [15][16][17]. Generally, if R 0 < 1, the infection will be eradicated as each infected individual suffers a disease, on average, to less than one other individual; if it is greater than one, then the outbreak will grow exponentially. For a degree-based SIS epidemic model, the expression of the basic reproduction number is given by with transmission rate β and treatment rate c, which serves a vital role in mitigating a disease transmission. However, when people face an emerging infectious disease, such as COVID-19, SARS, H1N1 etc., it may lead to a sharp increase in the number of patients and may cause serious runs of medical resources. A delay phenomena occurs since some infected individuals can not be instantly treated in hospitals. Such properties can be addressed by a Hill function by where α denotes a hysteretic effect due to medical limitations. Apparently, T(I) → 0 when I is small enough and it tends towards to c/α when I is large enough. The constant c/α represents the capacity for treatment. Some scholars proposed a type of step function [18] and a Verhulst-type function [19], which have similar properties to T(I).
In view of such epidemiological models incorporating a saturated treatment function, it is not hard to find that most of them enable such models, essentially changing their dynamics. Once a saturated function has been introduced, there always exists a backward bifurcation, which implies that even if some certain control measures make R 0 < 1, it does not guarantee the eradication of a disease in a region or a county. The presence of a backward bifurcation indicates that there exists a bistable phenomena, i.e, a stable endemic equilibria and a stable disease-free equilibrium coexist when R 0 < 1. The epidemic curve converging to zero or a positive constant is definitely determined by initial sizes of an infection. Some other epidemiological mechanisms, including partial immunity by vaccination [20], death caused by the disease [21], and susceptible heterogeneity [22] etc., have been identified to produce a backward bifurcation. However, the existing results are usually shown in homogeneous models. Up to date, there are few results to study backward bifurcation phenomena in epidemic models on complex networks from theoretical view of points. Recently, Li and Yousef proposed two degree-based SIR and SIS models with a saturated function to study a backward bifurcation phenomena. They showed a sufficient condition for the occurrence of a backward bifurcation [23,24]. Indeed, since network models have higher dimensions than those of well-mixed models, studying the complex dynamics of such models is a challenging issue from a theoretical view of points.
There are three main contributions in this paper. First, we propose a degree-based SIS epidemic model with a saturated function to study its long-term behaviors. Second, to overcome the difficulty of high dimensions for a network, we adopt an edgebased compartmental approach to lower the dimensions of an SIS epidemic model. Such an approach changes a complete degreebased model to a degree-edge-mixed model, and hence, it lowers the dimension of such a model from 2n to n + 1. Third, we proposed a geometric method to completely characterize the local stability of each equilibrium.
The organization of this paper is as follows: In Section 2, a degree-based SIS epidemic model on complex networks with a saturated function is proposed. Furthermore, we adopt an edgecompartmental approach to rewrite it as a degree-edge-mixed model. Section 3 gives a geometric approach to study the local stability of each equilibrium. In Section 4, we conducted some numerical simulations to illustrate our theoretical results. We give a brief discussion in the last section.

MODEL FORMULATION
In this paper, we focus on the complete stability of each equilibrium by a novel approach. Let us assume that the maximum contact number of an individual is n and then the degree set is N {1, 2, /, n}. The total subpopulation N k (t) is categorized into two statuses: susceptible and infective. Let S k (t) and I k (t) represent the densities of susceptible nodes and infected nodes at time t and degree k(k ∈ N), respectively. Inspired by the modeling approach [9], we account for the limitation of medical resources and propose a mean-field SIS epidemic model as follows: , where β denotes the transmission rate and c represents the treatment rate; α stands for the lag effect of the treatment due to the limitation of medical resources. From epidemiological view of points, the term denotes the probability of a given node connecting to an infected node at time t. Hence, it can be considered as a density of an [SI] edge at time t. Following the steps [10], model (1) can be rewritten as a degree-edge-mixed mode.
PROOF Observing the (n + 1) equation of system (2), we solve it and obtain where Obviously, Θ(t) > 0 if Θ(0) > 0. Now, we will prove the positivity of S k (t) for k ∈ N and t ∈ R + . By the continuity of the solution, if we claim that there exists a t 0 > 0 such that S k (t 0 ) 0 for some k ∈ N. Then S k (t) > 0 for t ∈[0, t 0 ) and S k ′(t 0 ) < 0. However, from the first n equation of (2), it follows that This leads to a contradiction with the claim. Therefore, S k (t) > 0 for k ∈ N + and t ∈ R + . Remark 2.1. Lemma 2.1 ensures that the solution of system (2) is strictly positive if Θ 0 > 0, which supports that I k (t) > 0 for all t ∈ R + and k ∈ N. The strong connectivity of a network guarantees the positivity of the solution as long as I k (0) ≠ 0 or there exists at least k 0 ∈ N such that I k 0 > 0.Noting that N k (t) S k (t) + I k (t) 1, we have from Lemma 2.1 that S k (t) ≤ 1, and I k (t) < 1, for all k ∈ N and t ∈ R + . Employing the expression of Θ, we assert that Therefore, for all t ∈ R + and k ∈ N, the set is positively invariant associate with system (2) and ϕ (S k0 , Θ 0 ).
In what follows, we focus on the long term behaviors of system (2) taking the initial data from Γ.

STABILITY OF EQUILIBRIA
In this section, we will consider the local stability of system (2) by a novel-geometric approach, which resolves such a matter once and for all. First, we try to give the basic reproduction number R 0 in an explicit form by the approach in [25]. Obviously, system (2) always has a disease-free equilibrium E 0 (1, 0). Linearising the n + 1 equation in system (2) yields to Solving Eq. 5 by a constant variation method, one drives a renew equation So that the basic reproduction number is calculated in form of The epidemiological meaning of R 0 is the average number of secondary infected edges produced by a typical [SI] edge during its infectious period. Now, let us account for the local stability of disease-free equilibrium E 0 . THEOREM 3.1. If R 0 < 1, then the disease-free equilibrium E 0 is locally asymptotically stable for any ϕ ∈ Γ.
PROOF Linearizing system (2) around the disease-free equilibrium E 0 results in Let us assume that system (8) has the solution with exponential forms, i.e, s k (t) s k0 e λt and θ(t) θ 0 e λt . Plugging them into (8), we have that If θ 0 ≠ 0, we cancel θ 0 on both sides of system (8) to derive that λ c(R 0 − 1). Otherwise, θ 0 0. From the first equation of (8), it follows that λ −c < 0. From what has been discussed, we conclude that if R 0 < 1, then the disease-free equilibrium E 0 is locally stable for any ϕ ∈ Γ.
PROOF Let us candidate a Lyapunov function by Differentiating V along the solution of system (2) leads to dV(t) dt (2) dΘ(t) dt here we have used the fact that Θ(t) < 1 and S k (t) ≤ 1 for all The equality holds if and only if Θ(t) 0. Hence, the largest invariant set M {ϕ ∈ Γ _ V(t) 0} only contains a singleton point E 0 . Consequently, the Lyapunov-Kasovskii-LaSalle theorem ensures that E 0 is globally asymptotically stable if R 0 < 1.Now, we are in a position to study the existence of the endemic steady state, whose components of that feasible equilibrium satisfy From the first equation of (10), we get Substituting Eq. 11 into the second equation of (10) and canceling cΘ * 1+αΘ * , we have that Apparently, if F(Θ) 0 has a solution Θ > 0, then system (10) has an endemic equilibrium. Note that F(Θ) has the following properties: Consequently, we conclude the following theorem on the existence of system (2). THEOREM 3.3 If R 0 > 1, then F(Θ) has at least one positive solution Θ * . In turn, system (2) has at least one positive endemic equilibrium E * .
From Theorem 3.3, we assert that system (2) has the existence of the endemic equilibrium, but it does not guarantee the uniqueness of the positive solution when R 0 > 1. However, if α is too small or large enough, it is easy to find that the function F is a decreasing function associated with Θ. Then we can claim the uniqueness of the endemic equilibrium whereas R 0 > 1. However, for R 0 < 1, we do not catch the existence of such equilibria due to the complex structure of F. To overcome that difficulty, we select a sensitive analysis for model parameters to investigate the existence of a backward bifurcation of that model for R 0 < 1. Now, we are concerned with the endemic curve which bifurcates backwards at (R 0 , Θ * ) (1, 0). If a backward bifurcation takes place, there always exist some additional conditions except R 0 < 1. To assess such a phenomena, we express β R 0 c 〈k〉 〈k 2 〉 from R 0 . Substituting this expression into Eq. 12, one arrives at Calculating the derivative of Eq. 14 with respect to R 0 by the Implicit Function Theorem, we obtain where After a simple computation, we have that LEMMA 3.4. Suppose R 0 < 1. A backward bifurcation occurs if the following inequality holds: Proof This is a direct result from Section 2.3 [26].
Next, let us move our attention on to the local stability of equilibria if they exist, which is a challenge issue for a degreebased epidemic model due to the complex structure. We will propose a geometric approach to deal with such an issue.Linearising system (2) around E * and taking the exponential perturbation solution similar to Theorem 3.1, we obtain If λ ≠ − βkΘ * + c 1+αΘ * , solving Eq. 18 yields to Replacing s k0 in Eq. 19 by Eq. 20 and canceling θ 0 , one admits Recall that β 〈k〉 Alternatively, taking the derivative of F with respect to Θ * leads to Plugging Eq. 23 into Eq. 22, we have that If λ − βkΘ * + c 1+αΘ * , then θ 0 0. From Eq. 19, it follows the positivity of β, k, p(k), and 〈k〉 that s k0 0. This case is impossible since (s k0 , θ 0 ) is an eigenvector of the eigenvalue λ − βkΘ * + c 1+αΘ * . From what has been discussed, we give some local stability of endemic equilibria. THEOREM 3.5. Let E * (S * k , Θ * ) be any feasible endemic equilibrium. The following statements are valid.
PROOF To address the stability of case (1), we rewrite H(λ) defined by Eq. 21 in the form of where L(k, Θ) βkΘ + c 1+αΘ , k ∈ N. Suppose (6) has a solution λ * with Reλ ≥ 0. Noting that Hence, if F′(Θ * ) > 0. Then H(λ) > 0 contradicts with Eq. 6. This implies that Eq. 6 has no solution with nonnegative real parts. Consequently, the endemic equilibrium E * is locally stable if F′(Θ * ) > 0.On the contrary, if F′(Θ * ) < 0, it follows from Eq. 21 that This, together with the Intermediate Value Theorem, ensures that Eq. 21 has at least one positive real solution. Hence, the endemic equilibrium E * is unstable. This completes the proof. If a backward bifurcation takes place, Lemma 3.4 and Theorem 3.5 give the stability of two positive endemic equilibria. THEOREM 3.6. If R 0 < 1, and α > 〈k〉〈k 3 〉 〈k 2 〉 2 , system (2) has two endemic equilibria. The one with a smaller quantitative of infected nodes is unstable; while the other, with a higher value of infected nodes, is locally asymptotically stable. THEOREM 3.7. Suppose R 0 > 1. If α ≪ 1, or α is large enough, then the endemic equilibrium E * is globally asymptotically stable if ϕ ∈ Γ\{E 0 }. PROOF Let us pick up a candidate Lyapunov function by where g(x) x − 1 − lnx with x > 0. Taking the derivative of V S along the trajectory of Eq. 2, we obtain dV S (t) dt On the contrary, differentiating V Θ along the solution of Eq. 2 leads to Adding Eqs 30, 31 together, one derives that Frontiers in Physics | www.frontiersin.org June 2021 | Volume 9 | Article 681268 If α ≪ 1 or it is large enough, then the term 1 + αΘ(t) 1 + αΘ * → 1, here we have used the fact that Θ * and Θ are both smaller than one. Therefore, the first and the last terms move towards to zero. Apparently, dV[S,Θ](t) dt ≤ 0. The equality holds if and only if Consequently, the largest invariant set of M {ϕ ∈ Γ _ V[S, Θ](t) 0} contains only a singleton point E * . From LaSalle's invariance principle, it follows that the endemic equilibrium E * is globally asymptotically stable.

NUMERICAL SIMULATION
In this section, we will proceed with some numerical experiments to validate our theoretical results. We account for an epidemic spreading on a scale-free network. Hence, we assume that the degree distribution of that network is p(k) ξk −2.5 and maximum degree n 100. Then ξ 1/ 100 k 1 p(k) 0.7458. Hence, the average degree, the second movement and the third movement of this network are 〈k〉 1.7995, 〈k 2 〉 13.8643, and 〈k 3 〉 500.7832.
First, we fix α 0.01. If we take β 0.032, we calculate R 0 0.9862 < 1 and R 0 0.9960 < 1. Theorem 3.2 ensures that the disease-free equilibrium E 0 is globally asymptotically stable. Let us pick up twenty different initial conditions for system (2). We find that the densities of infected edges associated with those initial data converge to zero (see Figure 1A). Enlarging β 0.33, we estimate R 0 1.071 > 1. Since α 0.01 ≪ 1, Theorem 3.7 ensures that the endemic equilibrium E * is globally asymptotically stable. Figure 1B displays that all the trajectories of system (2) move towards to a positive constant.
Second, we fix β 0.033 and hence R 0 1.071. We verify α and consider the existence and the uniqueness of the endemic equilibria. If we choose α 100 and it is large enough, Theorem 3.7 shows that the endemic equilibrium E * is globally asymptotically stable. Figure 2A depicts that the densities of infected edges from differential initial data converge to an endemic equilibrium Θ * 0.9546. If we pick up α 5, 10, 15, 20, 25, which are mediated values, all the trajectories of system (2) converge to their corresponding positive equilibria. From Figure 2B we find that increasing values of α gradually advances the arrival times and increases the sizes of the associated endemic equilibria. This indicates that enriching the adequate medical resources effectively reduce the risk of an infection.
Third, if we fix the structure of the network, then a key value α 0 〈k〉〈k 3 〉 〈k 2 〉 2 4.6883 determines whether or not a backward bifurcation happens when R 0 < 1. If we select α 6 > α 0 and β 0.032, then R 0 < 1. From Lemma 3.4, system exhibits a backward bifurcation at (R 0 , Θ * ) (1, 0). From Figure 3, it follows that the trajectories of system (2) partially converge to a positive level or partially move towards to zero. However, when we take α 4 < α 0 , however, R 0 4.9308 > 1 does not satisfy the conditions of Theorem 3.2. Figure 3B shows that the disease-free equilibrium E 0 is still asymptotically stable.
Finally, we want further insight into the existence of endemic equilibria when R 0 < 1. We take the parameters similar to case one except α. In this case, R 0 0.9862 < 1. As we know, the endemic Θ * is a function with respect to α and R 0 . If we change values of α from 4.8 to 6.8 with step 0.5. Figure 4 suggests that there always exists a backward bifurcation; In addition, the lag effect α controls the depth of the backward bifurcation; the larger α, the bigger the depth of the bifurcation. α is bigger in case of strength of lag effects of medical resources. Figure 4B displays the distributions of the positive equilibria in a given region (R 0 , α). In the blue region, system (2) has only unique disease-free equilibrium E 0 , which is globally asymptotically stable; In the grey region, a backward bifurcation occurs, i.e, system (2) exhibits a bistable phenomena, one endemic equilibria with a large value is stable and the other one with a small value is unstable; Moreover, the disease-free equilibrium E 0 is locally stable. In the green region, system (2) has a unique endemic equilibrium and it is stable.

DISCUSSION
In this paper, we considered a mean-field degree-based SIS epidemic model with a saturated treatment function. First, we adopted an edge-compartmental approach to simplify a pure degree-based model to a degree-edge-mixed model. Second, we proposed a novel method-a geometric approach to completely study the stability of each feasible equilibrium. The proposed model exhibits a backward bifurcation, i.e, R 0 < 1 does not sufficiently guarantee the eradication of an outbreak.
Compared with the results in [5,25], a degree-based SIS epidemic model on complex networks shows a threshold dynamic, in the sense that, if R 0 < 1, then the disease-free equilibrium E 0 is globally asymptotically stable; Otherwise, the unique endemic equilibrium E * is globally asymptotically stable. However, in our system, a saturated treatment function radically altered such a threshold property. Lemma 3.4 suggests that the hysteretic effect leads to the occurrence of backward bifurcation.
The basic reproduction number R 0 has no relationship with the lag effect of α. However, it is a vital value which determines whether or not a backward bifurcation occurs. In view of most existing results about that phenomena, they usually showed the existence and did not point out the stability of each equilibrium except some simple models with lower dimensions [2,3,18,19]. In this paper, we proposed a geometric approach to completely solve such an issue. Our results suggest that whereas R 0 < 1 and α > 〈k〉〈k 3 〉 〈k 2 〉 2 , one positive equilibrium with a larger value is stable; the other with small data is unstable. Their stabilities definitely depend on those signs of tangent slopes of the epidemic curve F(Θ). To establish it, we build a bridge between the local stability and the derivative of the epidemic curve, which is a universal result for any system even if it has a more complex structure [23].
Generally, the contact magnitude of an outbreak is characterized by the average degree 〈k〉 and its heterogeneity is measured by the two movement degrees 〈k 2 〉 of a network. To assess how the heterogeneity of a network affects the transmission of an infection, we fixed all the parameters as in part one of Section 4. Then we take the degree distribution function p(k) k −r with r −2.49, −2.5, −2.505, respectively. Then 〈k〉 1.840, 1.7995 and 1.7924. As 〈k〉 increases, the system undergoes three different phenomena: the extinction of an outbreak, the occurrence of a backward bifurcation, and the persistence of the disease (see Figure 5). This expounds that the system (2) undergoes a phase transition with the change of the network structure.
From an epidemiological viewpoint, the occurrence of a backward bifurcation implies that those control measures enabling R 0 < 1 do not efficiently ensures the eradication of an infection. In this case, the evolution of the disease heavily depends on the initial numbers of infected individuals. If it is small, the disease will die out; Otherwise, it will maintain a positive level and it cannot be radically eradicated from a region or a county. Hence, enhancing the management and allocations of medical resources plays an important role in slowing down a disease prevalence. Those control measures are in favor of effectively reducing the lag effect due to the limitations of medical resources as long as people face an emerging disease.
However, there are some limitations of this paper. First, we do not incorporate the population demography into the modeling process because introducing birth and death nodes  essentially changes the topology of a network [27,28]. This makes the model become too complex and then it has been become an unresolved issue to analyze its dynamical behaviors from mathematical view of points. Second, we do not couple individual contact data with some reported data for a realistic disease to study its evolutionary behaviors [29]. Third, we do not consider the convolution of information spread and epidemic transmission on multi-layered networks [30]. To carry out such a project, it may provide some valuable control suggestions for policymakers and public health government. We leave these for our future works.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.