key: cord-0796575-7d6qhdv1
authors: Zhang, Yunhu; Xiao, Yanni
title: Global dynamics for a Filippov epidemic system with imperfect vaccination
date: 2020-11-30
journal: Nonlinear Analysis: Hybrid Systems
DOI: 10.1016/j.nahs.2020.100932
sha: 6426b4a29be4f73f35f303b01c270e2fb4142029
doc_id: 796575
cord_uid: 7d6qhdv1

Abstract Given imperfect vaccination we extend the existing non-smooth models by considering susceptible and vaccinated individuals enhance the protection and control strategies once the number of infected individuals exceeds a certain level. On the basis of global dynamics of two subsystems, for the formulated Filippov system, we examine the sliding mode dynamics, the boundary equilibrium bifurcations, and the global dynamics. Our main results show that it is possible that the pseudo-equilibrium exists and is globally stable, or the pseudo-equilibrium, the disease-free equilibrium and the real equilibrium are tri-stable, or the pseudo-equilibrium and the real equilibrium are bi-stable, or the pseudo-equilibrium and disease-free equilibrium are bi-stable, which depend on the threshold value and other parameter values. The global stability of the disease-free equilibrium or pseudo-equilibrium reveals that we may eradicate the disease or maintain the number of infected individuals at a previously given value. Further, the bi-stability and tri-stability imply that whether the number of infected individuals tends to zero or a previously given value or other positive values depends on the parameter values and the initial states of the system. This emphasizes the importance of threshold policy and challenges in the control of infectious diseases if without perfect vaccines.

Mathematical models allow us to simulate the spread of diseases through different kinds of settings, model the possible interventions, and provide useful information to policymakers to quick response to curb the disease spread. There is a number of evidence showing that mathematical models play an essential role in investigating disease spread and identifying the key factors that significantly affect outbreaks [1] [2] [3] [4] [5] [6] [7] . It is well known that mathematical models have successfully helped understand and predict the trend of spread of diseases during the 2003 outbreak of Severe Acute Respiratory Syndrome (SARS) [8] [9] [10] [11] and the 2009 novel A/H1N1 pandemic influenza [12] [13] [14] . In fact, the modeling processes are usually the approximation to the real world, and modelers have been looking for the balance between the formulation of sufficiently fine/realistic models and feasibly manipulated models.

Interventions are usually modeled by continuous processes and disease spread with interventions are often described by the ordinary differential equations. A common assumption for the existing models is that human exploitative activities occur continuously [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] . However, this is not how the thing looks like. Mostly, in the early stages of the outbreak of emergent infectious diseases, the public is less aware of information on the disease itself, the degree of crisis or the necessary measures to protection and control, the disease spreads very quickly. As the disease further spreads or the number of infected individuals exceeds a certain number, the public then begin to implement some strategies to protect or control. Consequently, the threshold policy (TP) provides a natural description of such systems [7, 25] , and consequently Filippov systems become more realistic ones to represent density-dependent control strategies. These systems appear in almost every domain of applied sciences [26] [27] [28] [29] [30] , and in particular, in population dynamics [31] [32] [33] and epidemiology [34] [35] [36] [37] [38] .

Recently, a number of mathematical models have been proposed to investigate the effect of threshold policy on disease dynamics [31, [34] [35] [36] [37] [39] [40] [41] [42] . Existing approaches to model the impact of threshold policy have assumed that the threshold policy only influences on a certain type of population such that they implement control measures or change their behaviors. For example, only susceptible individuals enhance the protection strategies to let the incidence rate be reduced by a rate once the number of infected individuals exceeds a threshold level [35, 36] . However, not only susceptible but also (imperfect) vaccinated individuals enhance protection and control measures to avoid to be infected. How effective the threshold policy implemented by not only one compartment remains, is therefore an issue of great importance for epidemics control, and quantifying this policy through a mathematical modeling framework falls within the scope of this study.

The purpose of this paper is to investigate the effect of threshold policy on disease dynamics, and further to study what particular threshold level can be used to guide to eradicate the infectious disease. We formulate a general piecewise susceptible-infectious-vaccinated (SIV) type of model with nonlinear incidence to investigate the effects of TP implemented by both susceptible and vaccinated individuals. In Section 2, we examine two subsystems and their own dynamical behaviors. In the following section, the sliding mode dynamics is proposed to find out the pseudo-equilibrium. We prove that the pseudo-equilibrium is stable when it is feasible. In Section 4, we investigate the types and stability of all possible equilibria, and global dynamic behaviors by considering several scenarios. Further, the boundary equilibrium bifurcation is examined in Section 5. Finally concluding marks are given in the last section.

First we divide the total population N into three compartments: susceptible (S), infected (I) and vaccinated (V ) to establish the SIV model with threshold policy. In our model we assume that a fraction φ of the susceptible class is vaccinated per unit time. Although vaccination can reduce the infection by protecting susceptibles, it may not be completely effective. We let a factor of parameter σ (0 < σ < 1) to measure the efficiency of the vaccine as a multiplier to the infection rate, here σ = 0 means that the vaccine is perfectly effective, while σ = 1 means the vaccine has no effect at all. Given that the vaccination is not permanent we let θ be the rate of loss of immunity, then 1 θ is the immunity period after an effective vaccination. We divide each variable by N to normalize the variables, and hence in the following S, I, and V represent the proportions of three compartments, then we get

here H(I) is a given threshold function which may depend on the number of infected individuals, ϵ(t) is a piecewise function which is dependent on the sign of H(I). Constant µ is the natural birth (death) rate of population, β is the basic transmission coefficient, and γ represents the rate of recovery/remove from infected class. All the parameters are assumed to be positive constants. Model (1) with (2) is a description of the normalized threshold policy (TP), which is referred to as an on-off control. We denote the structure without intervention (ϵ(t) = 0) by the free system (S 1 ) and the structure with intervention (ϵ(t) = 1) by the control system (S 2 ). Here ϵ(t) = 1 means that when the number of infected individuals exceeds a certain number (I(t) > I c ), people begin to enhance some precaution and control measures such as wearing face masks, hand-washing and avoiding crowded places, and consequently the baseline transmission rate β is then reduced to be β(1 − f ). Positive constant f represents the intensity level of the implemented precaution and control strategies.

There is a variable structure with two distinct structures with their own equilibrium points due to the threshold policy.

The total population is supposed to be constant here, so we have (S + I + V ) ′ = 0. Substituting the variable 1 − I − V for S to yield the following two dimensional system

It is easy to get the invariant region by using the similar method in [43] 

Let H(U) = I − I c with the vector U = (V , I) for convenience. Then the discontinuous switching surface Σ can be defined

which divides the plane R 2 + into two separate parts

Here we characterize the right side of the equation in G i (i = 1, 2) to the form of column vector F i (U) with the elements F ij for convenient, therefore

so we can formulate the Filippov system by the following forṁ

Firstly, we will do the preliminary work to show the respective dynamics behavior for two subsystems.

Free system S 1 gives the following model equations

By employing the well-known next generation method in [44] , we can derive the basic reproduction number

It is easy to get the disease-free equilibrium E 1 0 = ( φ µ+θ +φ , 0). The endemic equilibria are solutions of the algebraic

One can get the quadratic equation for I

with

, we can easily obtain I 1 1 > 0 and I 1 2 < 0. When R 1 0 < 1 and B 1 < 0, we have I 1 1 > I 1 2 > 0. Then the second equation of (7) gives 

In fact, system S 1 has been carefully studied in [43] , and here we only introduce main results without giving detailed calculations. We easily know that there is a backward bifurcation at R 1 0 = 1 when the conditions B 1 < 0 and B 2 1 > 4A 1 C 1 are satisfied. There are two positive endemic equilibria when the parameter β changes from β c to β 1 , which are corresponding to B 1 = −2 √ A 1 C 1 and R 1 0 = 1 respectively. By computing the equation 

The disease-free equilibrium is E 2

The basic reproduction number for control system S 2 is

Using the same method as free system S 1 , we can derive the following quadratic equation for I

with

It follows that

Hence, for subsystem S 2 , we get the two equilibria E 2 1 (V 2

and only one feasible equilibrium E 2 1 when ∆ 2 > 0 and R 2 0 > 1 hold. Using the same method as free system S 1 , we get Lemma 2. The disease-free equilibrium E 2 0 is a LAS node when R 2 0 < 1 and unstable when R 2 0 > 1, while the endemic

there are precisely two endemic equilibria for subsystem (6) provided B 2 < 0 and ∆ 2 > 0, then we conclude that E 2 1 is a locally stable node, while E 2 2 is an unstable saddle point. When R 2 0 < R 1 c , the disease-free equilibrium E 2 0 is GAS. It is worth noting that both the two subsystems may have two positive equilibria, we can order them by the elementary of variable I as I 1 2 < I 2 2 < I 2 1 < I 1 1 , provided they are feasible (see detailed calculation in Appendix). Further, we note that R 2 c and R 1 c , which are independent of transmission coefficient β, are the same, so we denote them by R c , that is, we use

Now we briefly recall the definitions of sliding segment and crossing segment, then calculate the pseudo-equilibrium and give the sufficient conditions for existence of the pseudo-equilibrium. By using the Filippov convex method [45] , we let

where ⟨·, ·⟩ denotes the standard scalar product, H U (U) = ( ∂H ∂V , ∂H ∂I ) is the nonvanishing gradient on the discontinuity boundary Σ. It is well-known that the crossing segment Σ c ⊂ Σ is defined as

which means that the vectors F 1 (U) and F 2 (U) have the same sign with nontrivial normal components [45] . A sliding mode ensures that both trajectories in the vicinity of vector fields along I c are directing toward each other.

When there are some subregions of the line Σ such that both vectors of two subsystems S 1 and S 2 are directed towards each other, then the sliding segment appears on the boundary. We employ the equivalent control method as in [26] to examine the sliding domain and sliding mode dynamics of Filippov system (3) . So the sliding domain can be defined as

Here H = I − I c , which means H U (U) = (0, 1) T . Since F 1 (U) = (F 11 (U), F 12 (U)) T , F 2 (U) = (F 21 (U), F 22 (U)) T . When h(U) = ⟨H U (U), F 1 (U)⟩⟨H U (U), F 2 (U)⟩ = F 12 (U)F 22 (U) ≤ 0, note that F 12 (U) > F 22 (U) for any f > 0, it is impossible to be F 12 ≤ 0 and F 22 ≥ 0. Therefore, (11) is equivalent to the following

) denote the endpoints of the sliding segment. So we have the sliding domain of the Filippov system as follows 

with

Hence, there are two real roots for Eq. (12):

We now claim V 2

Then, the right hand side of (13)

Secondly, we want to find the sufficient conditions under which V 1

As A 2 > 0, the parabolic curve of G 2 is opening upward, then when

we have G 2 (I c ) > 0, and hence ψ(

As A 1 > 0, then when

it follows that G 1 (I c ) < 0, and consequently we have ψ(V R ) < 0.

In order to make sure V 1

, we need the sufficient condition to get ψ(V L ) > 0 and ψ(V R ) < 0. Therefore, when I c satisfies (14) and (15), i.e., I 1 2 < I c < I 2 2 or I 2

, then we get

which means that E 1 p is LAS, that is to say, the pseudo-equilibrium is locally stable when it is feasible. As E 2 p does not locate in the sliding segment, may not be named as a pseudo-equilibrium, we now use E p to denote the only possible pseudo-equilibrium (when it exist) in the following (i.e., E p = E 1 p ).

In this section, we examine the global dynamics of the piecewise system. We initially analyze the nonexistence of all possible limit cycles, that is, we exclude the existence of three possible limit cycles, the canard limit cycles, which contain part of the sliding segment, the limit cycles that totally in the region G i (i = 1, 2) and the limit cycles that surround the sliding segment. By combining the type (i.e., regular/sliding equilibria) and local stability of equilibria, we can obtain the reachability of orbits theoretically and numerically. Then we can conclude the global stability of the equilibria. To the end, we classify all the types of equilibria, and discuss the global dynamics of the proposed system in terms of relations of R 1 0 , R 2 0 , R c and 1. According to the previous conclusion, we get R 1

In such case we easily know that the endemic equilibria E 1 1 and E 2 1 are GAS for subsystem S 1 and S 2 according to Lemmas 1 and 2, respectively. In this case, I i 1 (i = 1, 2) is monotonously increase with β (the detailed proof is given in Appendix A.1), so we have I 1 Fig. A.1(A) ). According to Section 3, when the threshold value I c is between I 2 1 and I 1 1 , the sliding mode with a pseudo-equilibrium appears. In the following, we will examine the global analysis of the switching system (3) with (2). Here we take case (a) as an example to show non-existence of limit cycle.

Note that the existence of limit cycle of other cases are similar to do, hence we omit the detailed proof for other cases. 

By using the method in Lemma 1 [25] , we conclude that there is no limit cycle totally located in G 1 for the subsystem S 1 when R 1 0 > 1. Similarly, we get that no limit cycle for subsystem S 2 exists in region G 2 for R 2 0 > 1. Therefore, we have the following Lemma 3.

There is no closed orbit which is totally located in the region G i (i = 1, 2) for Case (a).

We now examine the limit cycle which contains part of the sliding segment or surround the sliding domain, and have the following two Lemmas.

There is no limit cycle which contains part of the sliding segment AB for Case (a).

Proof. By Section 3, we know that the sliding mode does exist for I 2 1 < I c < I 1 1 and the pseudo-equilibrium E p is feasible and LAS. This implies that trajectories which hit the sliding line will tend to the pseudo-equilibrium E p , hence there is no limit cycle containing part of the sliding line for I 2

We now need to prove this Lemma for the case I c < I 2 1 or I c > I 1 1 . Without loss of generality, we only prove the former.

on the sliding segment. So the trajectories move from the right to the left along the sliding line. Note that four isoclinic lines and the switching line I = I c divide the whole region into seven parts, in which the general trend of trajectories are identified (shown in Fig. 1(A) ). Trajectories hitting the sliding line will move from right to left and pass through the left endpoint A, as E 2 1 is the only real stable node in G 2 and A is a visible tangent point (see the definition in Section 5), so they move towards up and left in the domain delimited by the horizontal isocline g 2 I (red line), the vertical isocline g 2 V (green line) and the switching line I = I c (black dashed line), without hitting the switching line again. Then we know that no closed orbit containing part of the sliding segment exists for I c < I 2 1 . We can also use the similar method to prove the non-existence of limit cycle for I c > I 1 1 .

There is no closed orbit surrounding the sliding segment AB in Case (a).

Proof. Suppose there is a limit cycle L = L 1 + L 2 surrounding AB, which intersects with switching line I = I c at point P and Q , as shown in Fig. 1(B) . We add two auxiliary lines I = I c − ϵ and I = I c + ϵ in the form of a dotted line (∀ϵ > 0). The two lines intersect with L at points P 1 , Q 1 and P 2 , Q 2 respectively. Denote the region delimited by L 1 and P 1 Q 1 by K 1 in the lower part, while the region delimited by L 2 and P 2 Q 2 by K 2 in the higher part. The directed boundarẏ

Moreover, suppose that the abscissas of the points P, Q , 2, 3, 4) . Let the Dulac function be B = 1 VI , we can get the following results by the Green's theorem.

Sum up the both sides of the two equations above, we obtain the following

On the switching line, we have I = I c and dI dt = 0, which mean

Therefore, we have

On the other hand, − β V − φ(1−I)

This is an obvious contradiction with Eq. (16). So we have the conclusion that there is no limit cycle surrounding the sliding segment AB, the proof is then completed.

In the following, we will investigate the global stability of endemic equilibria in terms of relationship of I 1 1 , I 2 1 and I c . Here, three subcases are considered according to this relationship. Subcase (a1): assume I c > I 1 1 . In such subcase, the endemic equilibrium E 1 1 is real, while E 2 1 is virtual. We know that the sliding mode does exist but there is not a pseudo-equilibrium according to the previous derivation. We conclude that E 1 1 is GAS in the following. Proof. According to Lemma 1, the endemic equilibrium E 1 1 is real and is a LAS node for subsystem S 1 . And we know that any trajectory, once touching the sliding mode, moves from the left to the right along the sliding segment AB by the procedure of proof of Lemma 4. Note that there is no limit cycle totally located in the region G 1 or G 2 according to Lemma 3. Moreover, we know that there is no limit cycle which contains part of the sliding line or surrounding the sliding line AB by Lemmas 4 and 5. Hence, trajectories initiating from region G 1 will either go toward E 1 1 directly or hit the sliding line and slide along this line from the left side to the right endpoint B, then go toward E 1 1 ultimately. Trajectories starting from region G 2 will either cross over the switching line to get into the region G 1 or touch the sliding domain and then move along this line from left side to right in order to approach E 1 1 (shown in Fig. 2(B) ). Thus, all trajectories will tend to equilibrium E 1 1 ultimately. So the endemic equilibrium E 1 1 is GAS. This completes the proof.

Subcase (a2): assume I 2 1 < I c < I 1 1 . In such subcase both E 1 1 and E 2 1 are virtual, the sliding mode does exist, and the pseudo-equilibrium E p does exist and is LAS. Now we want to show that E p is GAS in the following.

Proof. In such case we easily know that the virtual equilibrium E 1 1 is in region G 2 while E 2 1 is in region G 1 . So trajectories go toward the opposite regions in order to approach their own equilibrium. Consequently, all the trajectories collide with the switching line. In view of Section 3, E p is locally stable. Hence, some trajectories which intersect with the switching line at the sliding segment move to the pseudo-equilibrium E p along the switching line, while others which intersect with the switching line at the crossing segment enter into the opposite region. Further, there is no limit cycle according to Lemmas 3-5, thus all the trajectories touch the sliding segment finally (shown in Fig. 2(C) ). Hence, we obtain that the pseudo-equilibrium E p is GAS. Subcase (a3): assume 0 < I c < I 2 1 . Note that this subcase is opposite of subcase (a1), and in such subcase the endemic equilibrium E 1 1 is virtual while E 2 1 is real, and there is no pseudo-equilibrium here. As there is no limit cycle, while the equilibrium E 2 1 of subsystem S 2 is globally stable in region G 2 (shown in Fig. 2(D) ). Then we get the following conclusion:

Theorem 3. The endemic equilibrium E 2 1 = (V 2 1 , I 2 1 ) is a GAS node if 1 < R 2 0 < R 1 0 and I c < I 2 1 .

In such case, there is only one globally stable endemic equilibrium E 1 1 for free subsystem S 1 , while there are two endemic equilibria E 1 2 and E 2 2 for control subsystem S 2 . In view of Lemmas 1 and 2, E 1 1 and E 2 1 are LAS nodes, while E 2 2 is a unstable saddle point, and the disease-free equilibrium E 0 is unstable for subsystem S 1 as R 1 0 > 1. There is a backward bifurcation at R 2 0 = 1 (i.e., C 2 = 0) for control system S 2 , then we have the order I 2 2 < I 2 1 < I 1 1 according to Fig. A.1(B) , and see detailed calculation in the Appendix. Again, we discuss the asymptotic behavior in terms of order of I 1 1 , I 2 1 , I 2 2 and I c .

Subcase (b1): assume I c > I 1 1 . The equilibrium E 1 1 is real and LAS, while E 2 1 is virtual but stable in this subcase. It is mentionable that the sliding mode exist but there is no pseudo-equilibrium. As we can preclude the existence of limit cycle using the same method as case (a), by analyzing the trend of any trajectory (similar to case (a)) we can obtain that E 1 1 is GAS (shown in Fig. 3(A) ).

Subcase (b2): assume I 2 1 < I c < I 1 1 . All the equilibria E 1 1 , E 2 1 and E 2 2 are virtual. So all trajectories initiating from G 1 go upward in order to approach to E 1 1 , while trajectories initiating from G 2 go downward in order to approach the equilibrium E 2 1 . As there is no limit cycle, the sliding mode happens and the pseudo-equilibrium appears (shown in Fig. 3(B) ), then the pseudo-equilibrium E p is GAS. Subcase (b3): assume I 2 2 < I c < I 2 1 . In this subcase, E 2 1 is a real and a stable equilibrium, while E 1 1 and E 2 2 are virtual. The sliding mode without pseudo-equilibrium exists. On the basis of non-existence of limit cycle for the whole system, we can confirm that E 2 1 is the only GAS equilibrium (shown in Fig. 3(C) ).

. Both E 2 1 and E 2 2 are real equilibria, while E 1 1 is virtual. Trajectories starting above Γ 2 (the stable manifold of E 2 2 , the pink line in Fig. 3(D) ) in region G 2 tend to equilibrium E 2 1 , while trajectories initiating below Γ 2 in G 2 will tend to equilibrium E 0 . Further, when trajectories starting from G 1 intent do go upward to approach the equilibrium E 1 1 , they collide with the switching line, and consequently the locally stable pseudo-equilibrium E p appears. Hence, both E 2 1 and E p are bi-stable (shown in Fig. 3(D) ). We summarize the above conclusion as the following:

Suppose R c < R 2 0 < 1 < R 1 0 , we conclude that the endemic equilibrium E 1 1 is GAS for I c > I 1 1 ; the pseudo-equilibrium E p is GAS for I 2 1 < I c < I 1 1 ; E 2 1 is GAS for I 2 2 < I c < I 2 1 ; and both E 2 1 and E p are bi-stable for I c < I 2 2 .

In such case endemic equilibrium E 1 1 of system S 1 is GAS in G 1 , while system S 2 has no positive equilibrium. We consider the following two subcases in terms of relations of I c and I 1 1 .

Subcase (c1): assume I c > I 1 1 . In such subcase, E 1 1 is a real equilibrium and is stable in G 1 , while the disease-free equilibrium E 1 0 is unstable for subsystem S 1 . Since subsystem S 2 has no endemic states, its disease-free equilibrium E 2 0 is stable in G 2 . Any trajectory initiating from G 2 intends to approach to the disease-free equilibrium, and then go down to the region G 1 by crossing switching line I = I c . As the limit cycles are ruled out, we can derive that the equilibrium E 1 1 is GAS (shown in Fig. 4(A) ). Subcase (c2): assume 0 < I c < I 1 1 . In this subcase, E 1 1 is in G 2 and is virtual. Thus, all trajectories initiating from G 1 intend to go upward in order to approach E 1 1 , while trajectories initiating from G 2 intend go downward in order to approach the equilibrium E 0 (E 2 0 ). Then two types of trajectories collide at the switching line, the sliding mode with a pseudo-equilibrium appears. Using the same method as before we rule out the existence of limit cycles, we can easily derive that the pseudo-equilibrium E p is GAS (shown in Fig. 4(B) ). Then we have the following conclusion:

0 , we conclude that the endemic equilibrium E 1 1 is GAS for I c > I 1 1 ; while the pseudo-equilibrium E p is GAS for I c < I 1 1 .

In such case, both the two subsystems have backward bifurcations for R i 0 = 1(i = 1, 2), respectively. As the order of the four biological roots are determined (see Fig. A.1(C) ), that is, I 1 2 < I 2 2 < I 2 1 < I 1 1 , we can classify all the types of phase portraits in term of the relation of I 1 1 , I 1 2 , I 2 1 , I 2 2 and I c .

Subcase (d1): assume I c > I 1 1 . In such subcase, the switching line I = I c is higher than all the endemic equilibria, E 1 1 and E 1 2 are real, while E 2 1 and E 2 2 are virtual. When trajectories starting from G 2 intend to go downward to approach E 2 1 , they touch the switching line, and either slide along the switching line or enter the region G 1 . While trajectories initiating from G 1 approach either E 1 1 or E 0 , depending on the starting points. In particular, when the starting point is above the stable manifold of E 1 2 (the pink line Γ 1 in Fig. 5(A) ), the trajectories move upward to tend to E 1 1 , when the starting point is below the stable manifold of E 1 2 , the trajectories go downward to E 0 . Hence both E 1 1 and E 0 are LAS in their own attraction domain (shown in Fig. 5(A) ), and then they are bi-stable in this subcase. Subcase (d2): assume I 2 1 < I c < I 1 1 . Here we know E 1 2 is real, while E 1 1 , E 2 1 and E 2 2 are virtual. Trajectories starting from G 2 intend to go downward to approach E 2 1 , after a slide or refraction, then enter into the region G 1 . While the trajectories initiating above the stable manifold of E 1 2 in G 1 intend to go upward to approach to E 1 1 . Then two types of trajectories touch the switching line, and slide to the pseudo-equilibrium. The trajectories initiating below the stable manifold Γ 1 of E 1 2 in region G 1 intend to go downward to tend to E 0 locally. Hence, both E p and E 0 are bi-stable (shown in Fig. 5(B) ).

Subcase (d3): assume I 2 2 < I c < I 2 1 . The endemic equilibriums E 1 1 and E 2 2 are virtual, while E 2 1 and E 1 2 are real. Then trajectories starting from G 2 intend to tend to E 2 1 , while the trajectories starting from the region above Γ 1 (the stable manifold of E 1 2 ) in G 1 intend to go upward to tend to E 1 1 , after a slide or refraction at the switching line, move upward to E 2 1 following the vector field for S 2 system. Trajectories starting from the region below the stable manifold Γ 1 of E 1 2 in G 1 go downward to approach E 0 locally (shown in Fig. 5(C) ). Therefore, the real equilibrium E 2 1 and the disease-free equilibrium E 0 are bi-stable in this subcase.

. Then we get E 1 1 is virtual, while E 1 2 , E 2 2 and E 2 1 are real. We divide the whole region into four parts according to the attraction domain separatrix and the switching line. The trajectories starting above Γ 2 (the stable manifold of E 2 2 ) in G 2 of subsystem 2 intend to tend to E 2 1 as the real equilibrium E 2 1 is locally stable, while trajectories starting below the stable manifold of E 2 2 in G 2 intend to go downward to E 0 . Further, trajectories starting above Γ 1 (the stable manifold of E 1 2 ) in G 1 intend to go upward to approach E 1 1 . When these two types of trajectories touch the switching line, the sliding mode with locally stable pseudo-equilibrium appears. Trajectories initiating from the region below Γ 1 in G 1 intend to go downward to E 0 , then E 0 is locally stable in this situation. As is shown in Fig. 5(D) , the real equilibrium E 2 1 , the pseudo-equilibrium E p and the disease-free equilibrium E 0 are locally stable. Hence three equilibria are tri-stable in such subcase.

Here we get that E 1 1 and E 1 2 are virtual, while E 2 1 and E 2 2 are real. Trajectories starting above Γ 2 (the stable manifold of E 2 2 ) in G 2 intend to tend to equilibrium E 2 1 . Trajectories initiating below Γ 2 in G 2 will tend to E 0 , and it is notable that when they go across the switching line, there is a slide or refraction on this line, and then these trajectories move to E 0 according to the vector field for subsystem S 1 . Meanwhile, it is obviously that trajectories starting from G 1 tend to E 0 directly. Hence, both E 2 1 and E 0 are bi-stable (shown in Fig. 5(E) ). Then we conclude the following: 

In such case, subsystem S 1 has two endemic equilibria, while subsystem S 2 have no biological equilibrium. So we consider the following three subcases in term of the relation of I 1 1 , I 1 2 and I c .

Subcase (e1): assume I c > I 1 1 . Under this condition, both E 1 1 and E 1 2 are real, E 0 is real too. All trajectories initiating above the stable manifold Γ 1 of E 1 2 in G 1 intend to go upward to E 1 1 , while trajectories initiating below Γ 1 in G 1 go downward to E 0 . And trajectories starting from G 2 intend to go downward to approach E 0 , then they collide at the switching line.

As there is no closed orbit, these trajectories will either slide along the switching line or refract at this line, and then continue to move according to vector field of system S 1 . Therefore, both E 1 1 and E 0 are bi-stable (shown in Fig. 6(A) ).

Subcase (e2): assume I 1 2 < I c < I 1 1 . Here we have E 1 1 is virtual, while E 1 2 is real. All trajectories starting from G 2 intend to go downward to approach E 0 , trajectories initiating above the stable manifold Γ 1 of E 1 2 in G 1 go upward to approach E 1 1 . These two types of trajectories then collide at the switching line, the sliding mode with pseudo-equilibrium appears. We know the pseudo-equilibrium E p is LAS. Meanwhile, trajectories initiating below Γ 1 in G 1 go downward to E 0 , indicating E 0 is LAS. Hence equilibria E 0 and E p are bi-stable (shown in Fig. 6(B) ). Subcase (e3): assume 0 < I c < I 1 2 . Both E 1 1 and E 1 2 are virtual, while E 0 is real. It is obviously that all the trajectories intend to go downward to tend to E 0 , and moreover we note that trajectories initiating from G 2 slide from left to right side along the switching line, or refract at the line. Therefore, E 0 is GAS (shown in Fig. 6(C) ). Then we conclude the following: 

In such case, there is no biological endemic state for subsystem S 1 or S 2 , the disease-free equilibrium E 0 = E 1 0 = E 2 0 is the only equilibrium for two subsystems. It is obvious that E 0 is always located in the G 1 region. Thus it is always real for subsystem S 1 but virtual for subsystem S 2 . As there is no other locally stable equilibrium or limit cycle, simply analyzing the trend of orbits can yield that the disease-free equilibrium is GAS for this scenario (shown in Fig. 6(D) ). In summary, we discuss the types and stability of all possible equilibria for Filippov system (3) , and obtain the global asymptotic behavior for all cases. We then summarize the stability of all possible equilibria in Table 1 . It is worth noting that the idea of proving global stability of Filippov system (3) is based on the local stability of equilibria, reachability of orbits (i.e., regular/sliding equilibria) and nonexistence of limit cycles. For this purpose we carefully analyze the types of various equilibria and trend of orbits in phase plane under different cases.

Now we investigate the boundary equilibrium bifurcation of the switching system. The readers can find the detailed definitions for the boundary equilibrium and the tangent point in Appendix.

Let E B be a boundary equilibrium, namely, E B satisfies the following equations:

So we can get two possible boundary equilibria by solving the above equations in (17) 

The equilibrium does not exist; GAS: globally asymptotically stable.

which are corresponding to ϵ = 0 or ϵ = 1, respectively. Here I c satisfies G 1 (I c ) = 0 or G 2 (I c ) = 0, so we have I c = I 1 1 or I 1 2 for G 1 (I c ) = 0, and I c = I 2 1 or I 2 2 for G 2 (I c ) = 0. Then we get four possible boundary equilibria.

Let T be a tangent point, according to the definition we can get the following equations:

It follows that the possible tangent points are

which are the solutions of Eq. (18) corresponding to ϵ = 0 and ϵ = 1, respectively. With the variation of the threshold value I c , the boundary equilibrium bifurcation occurs when the regular equilibrium collides with the tangent point and the boundary equilibrium. Fig. 7 illustrates a series of the boundary equilibrium bifurcations for case (d), in which each subsystem has two positive equilibria: a stable node and a saddle point. The real and stable node E 1 1 coexists with the visible tangent point T 2 for I c > I 1 1 (shown in Fig. 7(A) ). As I c decreases from I c > I 1 1 to I 1 1 , E 1 1 collides with T 2 (shown in Fig. 7(B) ). As threshold I c continues to decrease to I 2 1 < I c < I 1 1 , the stable pseudo-equilibrium E p appears and T 2 becomes an invisible tangent point (shown in Fig. 7(C) ). This bifurcation shows how a stable pseudo-equilibrium appears. Moreover, another boundary bifurcation occurs when I c passes through the critical value I 2 1 . The tangent point T 1 , the real node E 2 1 and the pseudo-equilibrium collide when I c = I 2 1 (shown in Fig. 7(D) ), then the pseudo-equilibrium E p disappears, and stable node E 2 1 becomes the locally attractor (shown in Fig. 7(E) ). When I c continuously decreases to I 2 2 , the third boundary bifurcation occurs, the visible tangent point T 1 collides with the saddle point E 2 2 (shown in Fig. 7(F) ). When I c further passes through I 2 2 to I 1 2 < I c < I 2 2 , the locally stable pseudo-equilibrium E p appears (Fig. 7(G) ), and the tri-stable phenomenon (E p , the disease-free equilibrium E 0 and real node E 2 1 ) emerges. The fourth boundary equilibrium occurs when I c passes through I 1 2 , the tangent point T 2 collides with saddle point E 2 2 and the pseudo-equilibrium E p for I c = I 2 2 (Fig. 7(H) ). When the threshold continues to decrease to be lower than I 1 2 , the pseudo-equilibrium E p disappears and tangent point T 2 becomes invisible ( Fig. 7(I)) , and consequently the disease-free equilibrium E 0 and the node point E 2 1 are bi-stable in this situation. Moreover, it is notable that the first two boundary equilibrium bifurcations are boundary node bifurcations, while the last two are boundary saddle bifurcations.

It has been observed that interventions, such as quarantine, isolation, treatment and vaccination, play a significant role in controlling the emerging and reemerging infectious diseases. However, not all interventions are implemented from the beginning of the outbreak or during the entire outbreak. Then the threshold policy can be naturally described by the Filippov system [35, 37, 39] . In this study, we extend the existing SIV-type model by including the extra control strategies once the number of infected individuals exceeds a certain level. In particular, based on the SIV-type model we consider susceptible and imperfect vaccinated individuals to enhance protection and control measures conditional upon relatively high prevalence.

The global dynamics of the Filippov system is fully investigated for different cases. It is interesting to note that two subsystems can undergo backward bifurcations, and consequently one endemic state and disease-free equilibrium for free or control system may be bistable when its specific reproduction number is less than one. The switching system may stabilize at regular equilibrium E 1 1 or E 2 1 , or the pseudo-equilibrium E p , which depends on the threshold value I c . The global stability of the pseudo-equilibrium indicates that when we choose the threshold value properly, the number of infected individuals can stabilize at a previously given value. These results are similar to those obtained by Wang and Xiao [25, 37] . For cases (a), (b) and (c), we observed some new phenomena. The disease-free equilibrium E 0 and E 1 1 , or E 0 and E 2 1 , or E 0 and E p may be bistable (as shown in Fig. 5(B,C,D) ). Especially, the disease-free equilibrium E 0 , E p and E 2 1 may be tri-stable, and the disease-free equilibrium E 0 may be GAS. These results imply that the dynamic behaviors of the system depend not only on the parameter values but also on the initial values of the system. Moreover, the threshold policy makes the whole system show more complex dynamical behaviors. Further, we analyze the boundary equilibrium bifurcation of the case (d). The four boundary bifurcations are illustrated with the variation of threshold value I c . All the evolutionary processes are displayed when I c passes across the four particular values I 1 1 , I 2 1 , I 2 2 and I 1 2 . It is worth noticing that the idea of proving global stability of the Filippov system (3) is based on the local stability of equilibria, reachability of orbits (i.e., regular/sliding equilibria) and nonexistence of limit cycles. For this purpose, we carefully analyze the types of various equilibria and trend of orbits in phase plane under different cases. Consequently, the proof of the global stability is a bit descriptive instead of using mathematical notations, which has also been intensely applied to investigating the global dynamics for other Filippov systems [27, 37, 38, 42] . Note that the proof of global stability can also be given by using mathematical notations like attractor, limit set and etc. We choose this way of description seems more readable and simple. The proposed switching system exhibits rich dynamics, implying that three equilibria including pseudo-equilibrium may be tri-stable or bi-stable. In particular, the pseudo-equilibrium or disease-free equilibrium may be GAS for particular conditions and suitable threshold level I c , which means that the disease can be eradicated or stabilize to a previously given level. Hence the modeling approach and main results give us a new understanding to realize the disease eradication and epidemic control. Definition 4 (Tangent Point). The tangent point T , which refers to the point at which F i (T ) ̸ = 0(i = 1, 2), and the trajectory of S i (i = 1, 2) tangents to the switching line Σ. That is to say, the tangent points satisfy ⟨F 1 (T ), H U (T )⟩ = 0 or ⟨F 2 (T ), H U (T )⟩ = 0.

Definition 5 (Visible (Invisible) Tangent Point). If the orbit ofU = F 1 (U) starting at T belongs to S 1 for all sufficiently small |t| ̸ = 0, we say that this tangent point is visible. While if this orbit belongs to S 2 , this tangent point is invisible. Similar definitions hold for the vector field F 2 (U).

We now give detailed calculation for the order of arrangement of the four roots I 1 1 , I 2 1 , I 2 2 and I 1 2 . Our purpose is to show the order of the four biological roots are determined when all the four or two or three of them are feasible.

Firstly, we want to show G 2 (I 1 1 ) > 0 when R 1 0 > 1. 

To simplify the expression, let M = (µ + θ + σ φ) + σ (µ + γ ), N = (µ + γ )(µ + θ + φ), K = µ + θ + σ φ, then we need to show G 2 (I 1

The both two sides of the inequation are multiplied by 

divided by σ 2 N 2 at each side of the inequation, one obtains

The left side is negative, while the right side is positive, the inequality is evidently valid so that (A.1) holds. Summing up the above two cases, we have G 2 (I 1 A.2. R 1

Secondly, when R 2 0 < 1, C 2 > 0, so I 2 2 > 0. We want to prove the inequality G 1 (I 2 1 ) < 0 holds when R 1 0 > 1 > R 2 0 . G 1 (I 2 1 ) = A 1 (I 2 1 ) 2

Consequently, it is obviously G 1 (I 2 1 ) < 0. Using the same method to get G 1 (I 2 2 ) < 0. Combine these two inequalities with G 1 (I 1 1 ) = 0, we have the relation I 1 1 > I 2 1 > I 2 2 (see Fig. A.1(B) ). Besides, when R 1 0 < 1, C 1 > 0, we can get the relation I 1 1 > I 2 1 > I 2 2 > I 1 2 by using the similar method (see sketch map in Fig. A.1(C) ).

Modeling and Dynamics of Infectious Diseases

Infectious Diseases of Humans

Bifurcation analysis of a discrete SIS model with bilinear incidence depending on new infection

The cost-effectiveness of oral HIV pre-exposure prophylaxis and early antiretroviral therapy in the presence of drug resistance among men who have sex with men in San Francisco

Global hopf bifurcation of a delayed equation describing the lag effect of media impact on the spread of infectious disease

Vaccination threshold size and backward bifurcation of SIR model with state-dependent pulse control

Effects of limited medical resource on a Filippov infectious disease model induced by selection pressure

Abu-Raddad, Transmission dynamics of the etiological agent of SARS in Hong Kong: impact of public health interventions

The outbreak pattern of SARS cases in China as revealed by a mathematical model

Mathematical model and prediction of epidemic trend of SARS

Analysis of the efficiency of the preventing and Isolating Treatments of SARS based on mathematical model

Modelling strategies for controlling H1N1 Outbreaks in China

Campus quarantine (Fengxiao) for curbing emergent infectious diseases: Lessons from mitigating A/H1N1 in Xi'an, China

Community-based measures for mitigating the 2009 H1N1 Pandemic in China

Seasonality and mixed vaccination strategy in an epidemic model with vertical transmission

A feedback control model of comprehensive therapy for treating immunogenic tumours

A delay SIR epidemic model with pulse vaccination and incubation times

Transmission dynamics of an influenza model with vaccination and antiviral treatment

Global results for an SIRS model with vaccination and isolation

Dynamics of an infectious disease where media coverage influences transmission

Media/psychological impact on multiple outbreaks of emerging infectious diseases

An SIS infection model incorporating media coverage

The impact of media on the spreading and control of infectious disease

The impact of media coverage on the transmission dynamics of human influenza

A Filippov system describing media effects on the spread of infectious diseases

Sliding Modes in Control and Optimization

Modelling the effects of cutting off infected branches and replanting on fire-blight transmission using Filippov systems

Sliding Mode Control in Electro-Mechanical Systems

Control of DC electric motor

Using Sliding modes in control theory

Media impact switching surface during an infectious disease outbreak

Threshold policies in the control of Predator-Prey Models

Optimal Foraging and Predator-Prey dynamics, Theor

Non-smoooth plant disease models with economic thresholds

Sliding mode control of outbreaks of emerging infectious diseases

Sliding bifurcation and global dynamics of a Filippov epidemic model with vaccination

A threshold policy to interrupt transmission of West Nile Virus to birds

A filippov model describing the effects of media coverage and quarantine on the spread of human influenza

Dynamics of an infectious diseases with media/psychology induced non-smooth incidence

Dynamics of a Filippov Epidemic model with limited Hospital Beds

Global dynamics of a piece-wise epidemic model with switching vaccination strategy

Multiscale system for environmentally-driven infectious disease with threshold control strategy

Backward bifurcations in simple vaccination models

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

One-parameter Bifurcations in Planar Filippov Systems

The work was supported by the National Natural Science Foundation of China (NSFC 11631012(YX)).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

In this appendix, we initially give some definitions of the various type of equilibria for the switching system [25, 45] Definition 1 (The Real (Virtual) Equilibrium). Let U * be such that F i (U * ) = 0 with threshold policy (2) . Then U * is called a real equilibrium of system (5) with (4) if it belongs to G i , and a virtual equilibrium if it belongs to G j , where i, j = 1, 2, (i ̸ = j). Pseudo-equilibrium) . A point U * is said to be a pseudo-equilibrium of system (5)