1. Introduction

Emerging and reemerging infectious disease, for instance, the 2003 severe acute respiratory syndrome (SARS) [27,19], the 2009 H1N1 influenza epidemic [22,38,44], Ebola virus disease in west Africa [42,23] and Zika virus [2] have become major causes of mortality and morbidity in emergencies, where collapsing health service and disease control programmes, poor access to health care, interrupted supplies and logistics are often the cruel reality. Global public health systems of surveillance and response are used to fighting to curb infectious diseases by controlling it at source [30] or slowing down its development course [20,5]. An efficient way to containing a disease during an outbreak is to provide as much as possible medical resources, including sufficient medicine, health care workers, hospital beds and etc. To explore the impact of limited medical resources on the control of infectious diseases and evaluate the possible control strategies, most of the available compartmental models either assume a constant treatment rate or adopt medical resources dependent continuous functions to represent the treatment rates [41,45]. However, change of the level of available medical resources were not well modeled in these models since usually the allocation and change of medical resources do not occur in a continuous way.

As WHO statistic information system referred, in-patient beds density, as one of the available indicators, can be adopted to access the level of health service delivery [43]. Hence, the number of available hospital beds per 10000 population, termed as hospital bed-population ratio (HBPR), is always used by the public health as an approach of capturing availability of health services delivery. The variations of the level of infected cases between communities and hospitals usually make it difficult for the public health to evaluate and decide the number of hospital beds needed to control the epidemic [1]. Setting sufficient beds to cover peak demands usually leads to quite a lot of beds idle at other times while reducing hospital beds to improving average bed occupancy level results in frequent overflows and congestions within hospital. Hence, quantifying the impact of the number of hospital beds for epidemics control remains a challenge and is therefore an important issue for public health.

Recently, there have been modeling studies dedicated to exploring the impact of the number of hospital beds on the future containing of disease [35,24,33,43,34]. Some have modeled the impact by embedding a continuous function in terms of the number of hospital beds into the recovery rate [35,34]. In fact, the medical resources used to contain the disease are always limited and the policy to intervene the disease varies depending on the number of infected individuals [32]. In particular, treatment rate depends not only on the available resources of health system to the public but also on the total infected individuals seeking treatment. On one hand, the capacity of the hospital settings, especially the number of health workforce and the facilities of the hospital are vital factors for efficient treatment of infection. On the other hand, the proportion of infected individuals getting access of the medical establishment to cure may change. Therefore, a considerable issue is to model and study the impact of these factors on the control of epidemic disease.

In this paper, we choose the number of infected cases as an index to propose a threshold policy [37,44,31,28] (of course, we could adopt the size of population who are exposed to the virus or the total number of susceptible and infected individuals). In this case, the control strategy is formulated by using a piecewise smooth function which depends on the density of infected individuals. The purpose of this paper is then to formulate a novel mathematical model subject to the threshold policy to study a Filippov system [21,39,8,4]. In particular, when the number of infected individuals is below the threshold level, the per capita treatment rate will be a constant, representing the maximum per capita treatment rate; while above the level, a weaker treatment policy will be adopted by incorporating the impact of the capacity and limited resources of the health care system in terms of the HBPR and depends on the number of the infected individuals. Filippov systems have been widely investigated [7,11,3,6], especially the sliding mode bifurcations in generic Filippov systems [17,26,18,29] and the numerical methods for them [15,13]. This type of systems have gained considerable attention in recent years since they not only behave differently from smooth systems but also found very important applications in many fields, such as pest control [37], power control in circuit [9], control of ecological systems [14], evolutionary biology [12,10], forest fire [16] and mechanical systems with friction [25], and a few work found in the disease controlling[32,40]. In [32], two parameters are adopted to model the treatment proportion for selective treatment measure due to the limited medical resources while we will model the hospital bed-population ratio by a parameter in this work.

The main purpose of this paper is to explore the impact of hospital bed-population ratio (i.e. HBPR) on disease control using a Filippov system. Some natural questions to ask include how does the level of infected cases combined with the HBPR affect the existence of the sliding modes, pseudo-equilibrium and sliding bifurcation? Further, how does the threshold policy affect the course of disease transmission and control outcome? To address these questions, we will focus on the sliding mode dynamics as well as the rich sliding bifurcations with variation of the threshold level. The rest of this paper is organized as follows. In Section 2, a piecewise defined treatment programme, which is termed as switching policy in control literature, is proposed to incorporate the level of infected cases and the impact of HBPR into the epidemic model. The sliding mode as well as all possible critical points are examined in Section 3. In Section 4, we study the sliding bifurcation including local sliding bifurcation and global sliding bifurcation. In Section 5, we investigate the impact of some key parameters on the evolution of the epidemic disease. Some concluding remarks and biological interpretations will be presented in the last section.

---

2. SIS Filippov model with switching treatment programme

As mentioned in the introduction, the medical resources can not always meet the increasing treatment demand of the infecteds. As WHO referred, we adopt the hospital bed-population ratio (i.e. HBPR and denoted by $b$) to evaluate the amount of medical resources in the following. Due to the limited medical resources, the treatment policy may vary as the amount of infecteds varies, so it is appropriate to set the number of infected individuals as the index. It is reasonable to assume that when the amount of infecteds is below a certain threshold level (denoted by $I_c$), the maximum per capita treatment rate (denoted by $h_1$) is in effect due to the relatively sufficient medical resources. As the number of infected individuals becomes larger and larger and once it exceeds $I_c$, only a certain ratio of infectious people can get treatment due to the insufficiency of the available medical resource. Note that the available treatment ratio may depend on the HBPR and the number of infected individuals. In particular, the treatment ratio decreases to $h_0$ as the amount of infecteds increases on one hand. On the other hand, the treatment ratio increases as HBPR increases. Thus we adopt a saturated function to represent the treatment policy, i.e. $H_1(b, I) = h_0+(h_1-h_0)b/(b+I)$ for this scenario. It is clear that $H_1(b, I)$ is bounded by the maximum per capita treatment rate $h_1$ and minimum per capita treatment rate $h_0$.

We divide the population into two classes: susceptible ($S$) and infected ($I$) individuals and formulate an SIS epidemic model by considering the effect of disease-induced mortality, which is important and appropriate for some epidemic diseases, such as tuberculosis, AIDS and plague. Then the model takes the following form

$\left\{ \begin{aligned} &\frac{dS}{dt} = \Lambda-\mu S-\frac{\beta SI}{S+I}+H(I, \epsilon)I, \\ &\frac{dI}{dt} = \frac{\beta SI}{S+I}-\mu I-\nu I-H(I, \epsilon)I \end{aligned} \right.$

(1)

with

$H(I, \epsilon) = \epsilon H_1(b, I)+(1-\epsilon)h_1$

(2)

and

$\epsilon = \left\{ \begin{aligned} &0, \sigma(S, I, I_c)<0, \\ &1, \sigma(S, I, I_c)>0, \end{aligned} \right.$

(3)

where $\Lambda$ represents the recruitment rate, $\mu$ is the natural death rate, $\beta$ denotes the basic transmission rate, $\nu$ stands for disease-induced mortality rate, the function $H(I, \epsilon)$ represents the per capita treatment rate. $\sigma(S, I, I_c)$ is a given function which may depend on the number of susceptible individuals, infected individuals and the threshold value $I_c$. But here we choose a particular form for function $\sigma(S, I, I_c)$, i.e. $\sigma(S, I, I_c) = \sigma(I, I_c) = I-I_c$, which indicates the amount of infected individuals is the index to guide which treatment policy works for different cases. Further, $H_1(b, I)$ is increasing as a function of the parameter $b$, where $b$ represents the HBPR and is used to characterize the impact of the capacity as well as the limited hospital resources. A detailed examination of the property on the treatment rate $H_1(b, I)$ can be found in [35].

In particular, the function $H(I, \epsilon)$ can be rewritten as the following form

$H(I, \epsilon) = \left\{ \begin{align} &{{h}_{1}},\,\,\,\,\, I<{{I}_{c}}, \\ &{{h}_{0}}+({{h}_{1}}-{{h}_{0}})\frac{b}{b+I},\,\,\,\,\, I>{{I}_{c}}. \\ \end{align} \right.$

Model system (1) with (2) and (3) is a piecewise smooth dynamical system (PWS)[3]. In particular, it subjects to a threshold value and is indeed a so-called Filippov system [21], i.e. systems of ordinary differential equations (ODES) with non-smooth right-hand sides.

Denote $Z = (S, I)$ and

$\begin{array}{rcl} X_1(S, I)& = &\displaystyle\bigg(\Lambda-\mu S-\frac{\beta SI}{S+I}+h_1I, \frac{\beta SI}{S+I}-\mu I-\nu I-h_1I\bigg)^T\\ &\dot{ = }& (f_{11}(S, I), f_{12}(S, I))^T, \\ X_2(S, I)& = &\displaystyle\bigg(\Lambda-\mu S-\frac{\beta SI}{S+I}+ \bigg[h_0+(h_1-h_0)\frac{b}{b+I}\bigg]I, \\ &&\displaystyle\frac{\beta SI}{S+I}-\mu I-\nu I-\bigg[h_0+(h_1-h_0)\frac{b}{b+I}\bigg]I\bigg)^T\\ &\dot{ = }&(f_{21}(S, I), f_{22}(S, I))^T, \end{array}$

then system (1) with (2) and (3) can be written as the following Filippov system

$\frac{\mathrm{d}Z}{\mathrm{d}t} = \left\{ \begin{array}{l} X_1(S, I), \sigma(I, I_c)<0\\ X_2(S, I), \sigma(I, I_c)>0. \end{array} \right.$

(4)

The attraction region of system (4) is $\Omega = \{(S, I)\in R^2_+:S+I\leq \Lambda/\mu\}$. For system (4), $\sigma = 0$ defines the switching manifold

$\Sigma = \{(S, I):I = I_c, (S, I)\in R^2_+\},$

which splits $R^2_+$ into two parts:

$\begin{aligned} &G_1(S, I) = \{(S, I):I< I_c, (S, I)\in R^2_+\}, \\ &G_2(S, I) = \{(S, I):I> I_c, (S, I)\in R^2_+\}. \end{aligned}$

For convenience, system (4) restricted in region $G_i$ is denoted as system $S_{G_i}, i = 1, 2$. Different from the smooth counterparts, system (4) possesses more than one equilibria [29,26], which play an important role in determining complex dynamics of the system.

Definition 2.1. Regular equilibrium of system (4) refers to those equilibria located in the region $G_i(i = 1, 2)$, which consists of real equilibrium and virtual equilibrium.

A real equilibrium $Z_*$ of Filippov system (4) (denoted by $Z_*^r$) refers to one which is an equilibrium of subsystem $S_{G_i}(i = 1, 2)$ and locates in the corresponding region $G_i$, i.e., it satisfies

$X_1(Z_*) = 0, \ \sigma(Z_*)<0$

or

$X_2(Z_*) = 0, \ \sigma(Z_*)>0.$

A virtual equilibrium $Z_*$ of Filippov system (4) (denoted by $Z_*^v$) refers to one which is an equilibrium of subsystem $S_{G_i}(i = 1, 2)$ but locates in the opposite region $G_j(j = 1, 2\ \mbox{and}\ j\neq i)$, i.e., it satisfies

$X_1(Z_*) = 0, \ \sigma(Z_*)>0$

or

$X_2(Z_*) = 0, \ \sigma(Z_*)<0.$

Definition 2.2. A boundary equilibrium $Z_*$ of Filippov system (4) refers to one which is an equilibrium of subsystem $S_{G_i}(i = 1, 2)$ but locates on the discontinuity boundary $\Sigma$, i.e. it satisfies

$X_1(Z_*) = 0, \ \sigma(Z_*) = 0$

or

$X_2(Z_*) = 0, \ \sigma(Z_*) = 0.$

Definition 2.3. A tangency point $Z_*$ of system (4) refers to a point $Z_*\in \Sigma$ (i.e. $\sigma(Z_*) = 0$) and $X_1\sigma(Z_*) = 0$ or $X_2\sigma(Z_*) = 0$, where $X_i\sigma(Z)(i = 1, 2)$ stands for the Lie derivative of $\sigma$ with respect of the vector field $X_i$ at the point $Z$, i.e.

$X_i\sigma(Z) = \langle\sigma(Z), X_i(Z)\rangle.$

We examine the dynamics of the subsystems $S_{G_1}$ and $S_{G_2}$ separately in the following.

Dynamics of subsystem $S_{G_1}$. Let $r_i = \mu+\nu+h_i(i = 1, 2)$. The Jacobian of system $S_{G_1}$ reads

$J_1(S, I) = \left( \begin{array}{cc} \displaystyle -\mu-\frac{\beta I^2}{(S+I)^2}& \displaystyle-\frac{\beta S^2}{(S+I)^2}+h_1\\ \displaystyle\frac{\beta I^2}{(S+I)^2}&\displaystyle\frac{\beta S^2}{(S+I)^2}-r_1 \end{array} \right).$

If we denote

$\delta_1 = (\mbox{tr}(J_1(S_1, I_1)))^2-4\mbox{det}(J_1(S_1, I_1))$

and the dynamics of system $S_{G_1}$ can be summarized in the following theorem.

Theorem 2.4. There is a disease-free equilibrium $E_0 = (\Lambda/\mu, 0)$ for system $S_{G_1}$, which is a stable node for $R_0<1$ and a saddle for $R_0>1$ with $R_0 = \beta/r_1$. The endemic equilibrium $E_1 = (S_1, I_1)$ exists if $R_0>1$ with

$S_1 = \frac{r_1I_1}{\beta-r_1}, \ I_1 = \frac{(\beta-r_1)\Lambda}{\mu r_1+(\mu+\nu)(\beta-r_1)},$

which is a stable node for $\delta_1\geq 0$ and a stable focus for $\delta_1<0$. No endemic equilibrium exists for $R_0 = 1$.

Dynamics of subsystem $S_{G_2}$. Denote

$a_0 = \beta\mu+\nu(\beta-r_0), a_1 = b\beta(\mu+\nu)+\Lambda(r_0-\beta)-b\nu r_1, \ a_2 = b\Lambda(r_1-\beta)$

and

$C_0 = [\beta b(\mu+\nu)-\Lambda(r_0-\beta)-b\nu r_1]^2+4\beta\mu b\Lambda(h_0-h_1).$

One can verify that there are two possible endemic equilibria of system $S_{G_2}$ if $C_0>0$ and $a_0>0$, denoted by $E_i = (S_i, I_i), i = 2, 3,$ with

$S_i = \frac{\Lambda-(\mu+\nu)I_i}{\mu}, \ I_2 = \frac{-a_1+\sqrt{C_0}}{2a_0}, \ I_3 = \frac{-a_1-\sqrt{C_0}}{2a_0}.$

Only one endemic equilibrium $E_* = (S_*, I_*)$ or $E_4 = (S_4, I_4)$ exists if $C_0 = 0$ or $a_0 = 0$, where

$\begin{aligned} &S_* = \frac{\Lambda-(\mu+\nu)I_*}{\mu}, I_* = \frac{-a_1}{2a_0}, \\ &S_4 = \frac{\Lambda-(\mu+\nu)I_4}{\mu}, I_4 = \frac{b\Lambda(\beta-r_1)}{b\nu(r_0-r_1)+\Lambda(r_0-\beta)}. \end{aligned}$

We address the existence of equilibria for system $S_{G_2}$ in detail according to whether the basic reproduction number $R_0$ is greater than unity and summarize the result as follows.

Theorem 2.5. The disease-free equilibrium $E_0$ always exists for system $S_{G_2}$ while the existence of possible endemic states is given in Table 1 .

The disease-free equilibrium $E_0$ is a node for $R_0<1$ and it is a saddle for $R_0>1$. To determine the type of the endemic equilibria, we initially compute the Jacobian of system $S_{G_2}$ as follows

$J_2(S, I) = \left( \begin{array}{cc} \displaystyle -\mu-\frac{\beta I^2}{(S+I)^2}& \displaystyle-\frac{\beta S^2}{(S+I)^2}+\frac{(h_1-h_0)b^2}{(b+I)^2}+h_0\\ \displaystyle\frac{\beta I^2}{(S+I)^2} & \displaystyle\frac{\beta S^2}{(S+I)^2}-\frac{(h_1-h_0)b^2}{(b+I)^2}-r_0 \end{array} \right).$

For any endemic equilibrium $\overline{E} = (\overline{S}, \overline{I})$, we denote

$\begin{aligned} &\mbox{tr}(J_2(\overline{S}, \overline{I})) = \frac{(-\beta+r_0-\mu)\overline{I}^2+ 2b(r_1-\beta-\mu)\overline{I}+(-\mu-\beta+r_1)b^2}{(b+\overline{I})^2}, \\ &\mbox{det}(J_2(\overline{S}, \overline{I})) = \frac{\mu\Big\{[(\beta-r_0)\Lambda+\nu b(r_1-r_0)]\overline{I}^2 +2\Lambda b(\beta-r_1)\overline{I}+\Lambda b^2(\beta-r_1) \Big\}} {(\Lambda-\nu \overline{I})(b+\overline{I})^2}. \end{aligned}$

One can get that the unique endemic equilibrium $E_2$ is a stable focus or node if $R_0>1$.

For the coexistence of $E_2$ and $E_3$, note that $\mbox{det}(J_2(S_*, I_*)) = 0$ and

$\mbox{sgn}\Big(\mbox{det}\big(J_2(\overline{S}, \overline{I})\big)\Big) = \mbox{sgn} \Big(\mathcal{A}\overline{I}^2+ 2\Lambda b(\beta-r_1)\overline{I}+\Lambda(\beta-r_1)b^2\Big)$

with

$\mathcal{A} = \big((\beta-r_0)\Lambda+\nu b(\beta-r_0)-\nu b(\beta-r_1)\big)>0,$

then

$\mbox{det}\big(J_2(S_2, I_2)\big)>0, \ \mbox{det}\big(J_2(S_3, I_3)\big)<0.$

This indicates that $E_3$ is a saddle while $E_2$ and $E_*$ are anti-saddles. Further discussion yields

$\mbox{tr}\big(J_2(S_2, I_2)\big) =\frac{(\nu+h_0-\beta)I_2^2+2b(\nu+h_1-\beta)I_2+(\nu+h_1-\beta)b^2} {(b+I_2)^2}.$

It follows that

$\begin{aligned} &\mbox{sgn}\Big(\mbox{tr}\big(J_2(S_2, I_2))\Big)\\ = &\mbox{sgn}\big((\nu+h_0-\beta)I_2^2+2b(\nu+h_1-\beta)I_2+(\nu+h_1-\beta)b^2\big)\\ = &\mbox{sgn}\Big\{\Big[a_1\big(2a_0bp_1-a_1p_0\big)+2a_0\big(a_2p_0-a_0b^2p_1\big)\Big] -\big(2a_0bp_1-a_1p_0\big)\sqrt{C_0}\Big\}, \end{aligned}$

where $p_0 = \beta-\nu-h_0, p_1 = \beta-\nu-h_1$. Denote

$\mathcal{B} = a_1\big(2a_0bp_1-a_1p_0\big)+2a_0\big(a_2p_0-a_0b^2p_1\big)$

and

$\delta_2 = \Big[\mbox{tr}(J_2(S_2, I_2))\Big]^2-4\mbox{det}(J_2(S_2, I_2)).$

Then one can verify that $E_2$ is a node for $\delta_2\geq 0$ while it is a focus for $\delta_2<0$. The further conclusion is as following:

(ⅰ) $E_2$ is stable if one of the following conditions holds

$\bullet$ $\beta\geq \nu+h_1$;

$\bullet$$\nu+h_0<\beta<\nu+h_1, \displaystyle\frac{p_1}{p_0}>\frac{a_1}{2a_0b}, (2a_0bp_1-a_1p_0)^2C_0>\mathcal{B}^2$.

(ⅱ)$E_2$ is unstable if $p_1/p_0\neq a_1/(2a_0b), (2a_0b p_1-a_1p_0)^2C_0\neq \mathcal{B}^2$ and both of the above conditions do not hold.

For the existence of endemic equilibrium $E_4$, it follows from

$\begin{aligned} \mbox{det}\big(J_2(S_4, I_4)\big)& = \frac{\mu(\beta-r_1)\mathcal{C}} {(b\nu+\Lambda)^2(\beta-r_0)^2(h_1-h_0)}>0, \\ \mbox{tr}\big(J_2(S_4, I_4)\big)& = (\beta-r_0-\mu)- \frac{\big[b\nu(h_0-h_1)+\Lambda(r_0-\beta)\big]^2}{(h_1-h_0)(\Lambda+b\nu)^2}<0, \end{aligned}$

where

$\mathcal{C} = -\mu\nu\beta b^2(h_1-h_0)^2+2b\mu\beta\Lambda(r_0-\beta)(h_1-h_0)-\Lambda^2(r_0-\beta)^3,$

we conclude that $E_4$ is a stable node or stable focus.

Theorem 2.6. (i) If $R_0>1$, the endemic equilibrium $E_2$ is a stable node or focus for subsystem $S_{G_2}$.

(ⅱ)If there are two endemic equilibria $E_2$ and $E_3$, then $E_3$ is a saddle while $E_2$ is an anti-saddle.

(ⅲ) If only one endemic equilibrium $E_*$ exists for the system, it appears as a stable node or stable focus.

By implementing similar discussion as in [35], we obtain the rich dynamics of system $S_{G_2}$ as the parameters vary. Here we omit the details and summarize the main result as follows.

Theorem 2.7. (i) If there are two endemic equilibria $E_2$ and $E_3$ for system $S_{G_2}$, Hopf bifurcation occurs at $E_2$ if

$\nu+h_0<\beta<\nu+h_1, \displaystyle\frac{p_1}{p_0}>\frac{a_1}{2a_0b}, (2a_0bp_1-a_1p_0)^2C_0 = \mathcal{B}^2.$

(ⅱ) Bogdanov-Takens bifurcation occurs if $R_0<1, a_0>0, a_1<0, C_0 = 0, -a_1/ (2a_0)<\Lambda/(\mu+\nu)$ and

$\nu+h_0<\beta<\nu+h_1, \displaystyle\frac{p_1}{p_0}>\frac{a_1}{2a_0b}, (2a_0bp_1-a_1p_0)^2C_0 = \mathcal{B}^2.$

---

3. Sliding mode and equilibria

To establish the dynamics given by the Filippov system (4), the first step is to rigorously define the local trajectory through a point $Z\in R^2_+$, especially those trajectories initiating from $Z\in \Sigma$ and it becomes the focus of this section. Those local trajectories initiating from the subregions $G_1$ and $G_2$ are defined by the subsystems $S_{G_1}$ and $S_{G_2}$ as usual. To extend the definition of a trajectory to the discontinuous boundary $\Sigma$, it is necessary to split $\Sigma$ into three distinct parts depending on whether the vector field points towards it [21,25].

$\bullet$ Attracting sliding mode $\Sigma_s$, where the vector fields points toward $\Sigma_s$ in both $G_1$ and $G_2$.

$\bullet$ Repulsing sliding mode $\Sigma_e$, where the vector fields points away from $\Sigma_e$ in either of $G_1$ and $G_2$.

$\bullet$ Transversal sliding mode $\Sigma_c$, where the vector fields points toward $\Sigma_c$ in $G_1(G_2)$ and away from $\Sigma_c$ in $G_2(G_1)$.

Existence of sliding mode. Indeed, no attracting sliding mode region exists for Filippov system (4) due to the mechanisms of our model. Solving $X_i\sigma(I, I_c) = 0\ (i = 1, 2)$ yields

$S_{c1} = \frac{(\mu+\nu+h_1)I_c}{\beta-(\mu+\nu+h_1)}, \ S_{c2} = \frac{[(b+I_c)(\mu+\nu+h_0)+b(h_1-h_0)]I_c}{(b+I_c) (\beta-(\mu+\nu+h_0))-b(h_1-h_0)}.$

The repulsing sliding mode region takes the form

$\Sigma_e = \{(S, I_c):S_{c2}\leq S\leq S_{c1}\}\dot{ = }\Sigma_{e1}$

for $R_0>1$ and

$\Sigma_e = \{(S, I_c):S\geq S_{c2}\}\dot{ = }\Sigma_{e2}$

for $R_0<1, \beta>r_0$ and $I_c>I_{c0}$, where $I_{c0} = b(r_1-\beta)/(\beta-r_0)$. For $R_0 = 1$, the sliding mode region is also $\Sigma_{e2}$. On the sliding mode region $\Sigma_{e_1}$ or $\Sigma_{e_2}$, it follows from Filippov's convex method [21,39] that the sliding mode dynamics reads

$\frac{\mathrm{d}S}{\mathrm{d}t} = -\mu\bigg[ S-\frac{\Lambda-(\mu+\nu)I_c}{\mu}\bigg]\dot{ = }F_s(S, I_c).$

(5)

Equation (5) defines a one dimensional dynamical system on $\Sigma_{e1}$ or $\Sigma_{e2}$.

We will define the type of equilibrium on the repulsing sliding mode $\Sigma_e$ [29], which is special for Filippov system and plays an important role in the analysis of global behavior.

Definition 3.1. A point $Z_*\in \Sigma_e$ is called a pseudo-equilibrium of Filippov system (4) if it satisfies $F_s(Z_*) = 0$, which suggests a pseudo-equilibrium is indeed an equilibrium of the sliding mode dynamics. We call pseudo-saddle to any pseudo-equilibrium $Z_*\in \Sigma_e$ such that $F_s^{\prime}(Z_*)<0$.

The following definition gives a more detailed characterization of a tangency point.

Definition 3.2. A point $Z_*\in\Sigma$ is termed as a fold tangency point if one of the following conditions hold.

(a) $X_1\sigma(Z_*) = 0,\,\,\,\, X_1^2\sigma(Z_*)\neq 0,\,\,\,\, X_2\sigma(Z_*)\neq 0.$

(b) $X_2\sigma(Z_*) = 0,\,\,\,\, X_2^2\sigma(Z_*)\neq 0,\,\,\,\, X_1\sigma(Z_*)\neq 0.$

Furthermore, if $X_2\sigma(Z_*)>0$ in possibility (a) or $X_1\sigma(Z_*)<0$ in possibility (b), then $Z_*\in \partial\Sigma_e$ and it appears as a repulsing sliding fold tangency.

For case (a), we call a fold tangency point visible (invisible) if the trajectories of subsystem $S_{G_1}$ initiating from $Z_*$ remains in region $G_1(G_2)$ for any sufficiently small period. It follows that a fold tangency point is visible (invisible) if we further have $X_1^2\sigma(Z_*)<0 (X_1^2\sigma(Z_*)>0)$ in possibility (a). Similar definitions hold for case (b).

An equilibrium $S_s = \big(\Lambda-(\mu+\nu)I_c\big)/\mu$ exists for the sliding mode dynamics (5), which suggests a possible pseudo-equilibrium $E_s = (S_s, I_c)$ exists for system (4) by Definition 3.1. The feasibility of pseudo-equilibrium for system (4) in different cases depends on whether $E_s$ lies in the domain of repulsing sliding mode region $\Sigma_{e1}$ or $\Sigma_{e2}$.

Note that

$S_s-S_{c1} = \frac{(I_1-I_c)\big[\beta\mu+\nu(\beta-r_1)\big]} {\mu(\beta-r_1)},$

we have

$S_s\leq S_{c1}\Longleftrightarrow I_c\geq I_1$

for $R_0>1$. It follows from

$S_s-S_{c2} = \frac{\big[\beta\mu+\nu(\beta-r_0)\big]I_c^2+ \big[b\beta(\mu+\nu)-\Lambda(\beta-r_0)-b\nu r_1\big]I_c-b\Lambda(r_1-\beta)} {-\mu(\beta-r_0)(I_c-I_{c0})}$

that $S_s\geq S_{c2}$ provided one of the following set of inequalities satisfies.

$\bullet$ $R_0>1, I_c\leq I_2$;

$\bullet$ $R_0<1, \beta>r_0, I_c>I_{c0}, I_3\leq I_c \leq I_2$;

$\bullet$ $R_0 = 1, I_c <I_{c0}$.

Then we conclude that the pseudo-equilibrium $E_s$ is feasible if one of the following conditions holds.

(a) $R_0>1, I_1\leq I_c\leq I_2$;

(b) $R_0<1, \beta>r_0, I_c>I_{c0}, I_3\leq I_c\leq I_2$;

(c) $R_0 = 1, I_c < I_{c0}$.

Equilibria. We address all possible critical points for system (4) as follows.

Regular Equilibrium. In terms of Definition 2.1, the disease-free equilibrium $E_0$ is real and independent of the size of threshold level $I_c$; the endemic equilibrium $E_1$ is real provided $I_c>I_1$ while it is virtual provided $I_c<I_1$; the endemic equilibrium $E_i(E_*)$ is real for $I_c<I_i(I_c<I_*)$ while $E_i(E_*)$ is virtual for $I_c>I_i(I_c>I_*)$, $i = 2, 3, 4$.

Boundary Equilibrium. There are up to five equilibria of Filippov system (4) colliding with the switching boundary $\Sigma$ at distinct moments, so five possible boundary equilibria exist according to Definition 2.2, i.e.

$E_1^b = \bigg(\frac{r_1I_c}{\beta-r_1}, I_{c1}\bigg), \ E_i^b = \bigg(\frac{\Lambda-(\mu+\nu)I_c}{\mu}, I_{ci}\bigg), i = 2, 3, 4, 5$

if $I_c = I_{ci}, I_{ci} = I_i$ for $E_i^b$ with $i = 1, 2, 3, 4$ and $I_{c} = I_{c5}, I_{c5} = I_*$ for $E_5^b$.

Pseudo-equilibrium. Only a pseudo-equilibrium $E_s$ is feasible if one of the following conditions satisfies.

(a) $R_0>1, I_1\leq I_c\leq I_2$;

(b) $R_0<1, \beta>r_0, I_c>I_{c0}, I_3\leq I_c\leq I_2$;

(c) $R_0 = 1, I_c < I_{c0}$.

Furthermore, the pseudo-equilibrium $E_s$ is a pseudo-saddle if it is defined well on the sliding mode region $\Sigma_{e1}$ or $\Sigma_{e2}$.

Tangency point. By Definition 2.3, two tangency points coexist for system (4) if $R_0>1$, represented by

$E_1^t = (S_{c1}, I_c), E_2^t = (S_{c2}, I_c),$

and there is one tangency point $E_2^t$ if $R_0<1, \beta>r_0, I_c>I_{c0}$ or $R_0 = 1$.

According to Definition 3.2, fold tangency points can be classified as visible and invisible ones, so we examine the type of tangency points $E_i^t (i = 1, 2)$ as follows.

For $E_1^t$, it follows from

$\left\{ \begin{aligned} &\mbox{sgn}\big(X_2\sigma(E_1^t)\big) = \mbox{sgn}\bigg( \frac{\beta S_{c1}}{S_{c1}+I_c}-\frac{b(h_1-h_0)}{b+I_c}-r_0\bigg) = 1, \\ &X_1^2\sigma(E_1^t) =\frac{\beta I_{c}^2}{(S_{c1}+I_{c})^2}f_{11}(S_{c1}, I_{c})+ \Bigg[\frac{\beta S_{c1}^2}{\big(S_{c1}+I_{c}\big)^2}-r_1\Bigg]f_{12}(S_{c1}, I_{c}), \end{aligned} \right.$

that

$\mbox{sgn}\big(X_1^2\sigma(E_1^t)\big) = \mbox{sgn}\Bigg( \frac{\Lambda(\beta-r_1)-\mu r_1I_{c}-(\mu+\nu)(\beta-r_1)I_{c}}{\beta-r_1}\Bigg),$

so if $I_c < I_1$, one gets $X_1^2\sigma(E_1^t)>0$ and then $E_1^t$ is invisible in such case; if $I_c>I_1$, one has $X_1^2\sigma(E_1^t)<0$ and then $E_1^t$ is visible. Furthermore, since $X_2\sigma(E_1^t)>0$, the fold tangency point $E_1^t$ is a repulsing fold tangency.

For the tangency point $E_2^t$, similar argument gives

$\left\{ \begin{aligned} &\mbox{sgn}\big(X_1\sigma(E_2^t)\big) = \mbox{sgn} \big((\beta-r_1)S_{c2}-r_1I_c\big) = -1, \\ &X_2^2\sigma(E_2^t) = \!\frac{\beta I_c^2}{(\!S_{c2}\!+\!I_c\!)^2}f_{21}(S_{c2}, I_c)\!\!+ \!\!\Bigg[\!\!\frac{\beta S_{c2}^2}{\big(\!S_{c2}\!+\!I_c\!\big)^2}\!-\! \frac{(\!h_1\!-\!h_0\!)b^2}{(b+I_c)^2}\!-\!(\!\mu\!+\!\nu\!+\!h_0\!)\!\Bigg] \!f_{22}(S_{c2}, \!I_c), \end{aligned} \right.$

so

$\begin{aligned} &\mbox{sgn}\Big(X_2^2\sigma(E_2^t)\Big)\\ = &\mbox{sgn}\Bigg\{ \frac{-\big[\mu\beta+\nu(\beta-r_0)\big]I_c^2+ \big[\Lambda(\beta-r_0)-b(\beta\mu+\beta\nu-r_1\nu)\big]I_c+\Lambda b(\beta-r_1)} {(\beta-r_0)I_c+b(\beta-r_1)}\Bigg\}. \end{aligned}$

(6)

There are two possibilities according to whether $R_0$ is greater or less than unity. If $R_0>1$, one gets that $E_2^t$ is visible for $I_c<I_2$ and invisible for $I_c>I_2$. If $R_0<1$, since the sliding mode region exists for $\beta>r_0$ and $I_c>I_{c0}$, we focus our examination on such case and there are the following three subcases to consider in terms of the existence of endemic equilibrium for system $S_{G_2}$. For the coexistence of two endemic equilibria $E_2$ and $E_3$, one gets

$I_3 < I_c<I_2\Rightarrow X_2^2\sigma(E_2^t)>0$

according to (6), so $E_2^t$ is visible; while we get $X_2^2\sigma(E_2^t)<0$ if $I_c>I_2$ or $I_c<I_3$, so $E_2^t$ is invisible. (ⅱ) For the existence of endemic equilibrium $E_*$, if $I_c\neq I_*$, we have $X_2^2\sigma(E_2^t)<0$, so $E_2^t$ is invisible. (ⅲ) For the nonexistence of endemic equilibrium, we derive $X_2^2\sigma(E_2^t)<0$, so $E_2^t$ is invisible. Moreover, it follows from $X_1\sigma(E_2^t)<0$ that the fold tangency point $E_2^t$ is a repulsing fold independent of the level of threshold $I_c$.

Once we have examined the sliding dynamics as well as all possible critical points for system (4), we turn to discuss the discontinuity-induced bifurcation in the following section, which is special for Filippov system and refer to a type of bifurcations involving structural changes in the sliding mode domain.

---

4. Sliding bifurcation

In this section, we focus our attention on sliding bifurcation, in which some sliding segment on the switching boundary is involved. According to [29,17], all sliding bifurcations can be classified as local and global bifurcations. The local sliding bifurcations include boundary equilibrium bifurcation and double tangency bifurcation. Those sliding bifurcations are called global if they involve nonvanishing cycles. Global bifurcations include grazing (touching) bifurcation, sliding homoclinic bifurcation to pseudo-saddle, or pseudo-saddle-node or saddle, and the sliding heteroclinic bifurcation between pseudo-saddles, etc. In the following, we choose the threshold value $I_c$ as a bifurcation parameter and fix all other parameters unchanged.

---

4.1. Local sliding bifurcation

We initially examine the local sliding bifurcation for system (4), which include boundary equilibrium bifurcation and double tangency bifurcation. There are four distinguished boundary equilibrium bifurcations, i.e. boundary node (BN), boundary focus (BF), boundary saddle (BS) and boundary saddle-node bifurcation (BSN) [29,26]. Further, if the boundary equilibrium bifurcation involves changes of the real/virtual character, it is also termed as real/virtual equilibrium bifurcation. Double tangency bifurcation is a type of bifurcations triggered by the collision of two tangency points. Since the tangency points $E^t_1$ and $E^t_2$ are distinct and impossible to collide for $I_c>0$, no double tangency bifurcation will occur. Thus, in the following we are devoted to the boundary equilibrium bifurcation, which is induced by the collision of a hyperbolic equilibrium with the discontinuity boundary $\Sigma$. For such purpose, we introduce the following definition [17].

Definition 4.1. A boundary equilibrium bifurcation occurs at a critical point $Z_*\in\Sigma$ if $X_1(Z_*) = 0$ (or $X_2(Z_*) = 0$) and $DX_i(Z_*)$ is invertible (i.e. $\mbox{det}(DX_i(Z_*))\neq 0$) for $i = 1, 2$.

When a boundary equilibrium bifurcation occurs, we can observe two cases. One is persistence, where a real equilibrium becomes a pseudo-equilibrium; the other is non-smooth fold, where a real equilibrium disappears after colliding together with a pseudo-equilibrium. Implementing a similar discussion as in [17], one gets the conditions to distinguish between the two possible types of unfolding of a boundary equilibrium as the parameter $I_c$ is perturbed from the critical value. Without loss of generality, we take the boundary equilibrium $E_1^b$ as an example and state the conditions as following.

Definition 4.2. For system (4), if

$\begin{array}{rcl} \mbox{det}(J_1(S_{ic1}, I_{c1}))&\neq&0\\ \displaystyle\frac{\partial \sigma(I_{c1}, I_c)}{\partial I_c}-J_3(S_{ic1}, I_{c1})J_1^{-1}(S_{ic1}, I_{c1}) J_{4}(S_{ic1}, I_{c1})&\neq&0\\ J_3(S_{ic1}, I_{c1})J_1^{-1}(S_{ic1}, I_{c1})\big(X_2(S_{ic1}, I_{c1})- X_1(S_{ic1}, I_{c1})\big)&\neq&0, \end{array}$

a non-smooth fold is observed for

$J_3(S_{ic1}, I_{c1})J_1^{-1}(S_{ic1}, I_{c1})\big(X_2(S_{ic1}, I_{c1})-X_1(S_{ic1}, I_{c1})\big)<0$

while persistence is derived at the boundary equilibrium bifurcation point for

$J_3(S_{ic1}, I_{c1})J_1^{-1}(S_{ic1}, I_{c1})\big(X_2(S_{ic1}, I_{c1})-X_1(S_{ic1}, I_{c1})\big)>0,$

where $J_3(S, I)$ stands for the Jacobian of $\sigma(I, I_c)$ with respect to the the state variable $(S, I)$, $J_4(S, I)$ represents the Jacobian of $X_1$ with respect of the parameter $I_c$, $S_{ic1}$ denotes the abscissa of the boundary equilibrium $E_1^b$.

The similar condition can be presented for other boundary equilibria and we omit them here.

For the case that only one endemic equilibrium $E_2$ exists for system $S_{G_2}$, there are two boundary equilibria $E_1^b$ and $E_2^b$ for Filippov system (4). According to section 2, $\mbox{det}(J_2(E_2^b))>0$ for $R_0>1$. Further calculation yields

$\begin{aligned} \mbox{det}(J_1(E_1^b))& = \frac{\mu r_1I_1}{S_1+I_1}+\frac{(\mu+\nu)\beta I_1^2}{(S_1+I_1)^2}, \\ \mbox{det}(J_1(E_2^b))& = \frac{\beta(\mu+\nu)I_2^2}{(S_2+I_2)^2}+ \mu\bigg[\frac{\beta S_2I_2}{S_2+I_2}+\frac{(h_1-h_0)I_2}{b+I_2}\bigg], \\ \mbox{det}(J_2(E_1^b))& = \frac{\big[\nu(r_1+\beta)^2+\beta\mu(\beta-r_1)\big] \big[(\beta-r_1)(\Lambda+b\nu)+b\beta\mu\big]^2-\mathcal{C}} {\beta\big[(\beta-r_1)(\Lambda+b\nu)+b\beta\mu\big]^2} \end{aligned}$

with

$\mathcal{C} = \beta\mu(r_1-r_0)\big[(\beta-r_1)^2(\Lambda^2+2b\nu\Lambda)+ 2b\beta\mu\Lambda(\beta-r_1)\big],$

so a boundary focus (or node) bifurcation occurs for $R_0>1$ when $I_c$ passes through the critical value $I_{c2}$. If $\mbox{det}J_2(E_1^b)\neq 0$, a boundary focus (node) bifurcation happens for $R_0>1$ when $I_c$ passes through another critical value $I_{c1}$. Indeed, a boundary equilibrium bifurcation is triggered by the collision of the regular equilibrium, pseudo-equilibrium and tangency point (or two of them) as the threshold value $I_c$ passes through the critical level.

For boundary equilibrium $E_1^b$, the first three conditions in Definition 4.2 hold, so according to Definition 4.2, the type of unfolding depends on the sign of the following formula

$\begin{aligned} &\mbox{sgn}\Big(J_3(S_{ic1}, I_{c1})J_1^{-1}(S_{ic1}, I_{c1}) \big(X_2(S_{ic1}, I_{c1})-X_1(S_{ic1}, I_{c1})\big)\Big)\\ = &\mbox{sgn}\bigg\{ \frac{(h_1-h_0)I_{c1}}{b+I_{c1}}\bigg[h_1-\mu-\frac{\beta I_{c1}^2}{(S_{ic1}+I_{c1})^2}- \frac{\beta S_{ic1}^2}{(S_{ic1}+I_{c1})^2}\bigg] \bigg\}\\ = &\mbox{sgn}\bigg\{ \frac{(h_1-\mu-\beta)r_1^2+2(h_1-\mu)r_1(\beta-r_1)+(h_1-\mu-\beta)(\beta-r_1)^2}{(\beta-r_1)^2} \bigg\}\\ = &\mbox{sgn}\big(-2r_1^2+\beta(r_1+\nu+2h_1-\beta)\big). \end{aligned}$

Thus, a non-smooth fold occurs if

$-2r_1^2+\beta(r_1+\nu+2h_1-\beta)>0,$

(7)

while persistence is derived if

$-2r_1^2+\beta(r_1+\nu+2h_1-\beta)<0.$

(8)

Implementing the similar process for boundary equilibrium $E_2^b$, we derive that the type of unfolding of $E_2^b$ is determined by

$\begin{aligned} &\mbox{sgn}\Big(J_3(S_{ic2}, I_{c2})J_2^{-1}(S_{ic2}, I_{c2}) \big(X_1(S_{ic2}, I_{c2})-X_2(S_{ic2}, I_{c2})\big)\Big)\\ = &\mbox{sgn}\bigg\{ \frac{(h_1-h_0)I_{c2}}{b+I_{c2}}\bigg[-h_0+\mu-\frac{(h_1-h_0)b^2}{(b+I_{c2})^2}+ \frac{\beta I_{c2}^2}{(S_{ic2}+I_{c2})^2}+ \frac{\beta S_{ic2}^2}{(S_{ic2}+I_{c2})^2}\bigg]\bigg\}\\ = &\mbox{sgn}\bigg\{ (\mu-h_0)-\frac{b^2(h_1-h_0)}{(b+I_{c2})^2}+\frac{\beta\big[ \Lambda-(\mu+\nu)I_{c2}\big]^2+\beta\mu^2I_{c2}^2}{(\Lambda-\nu I_{c2})^2} \bigg\}\\ = &\mbox{sgn}\bigg\{\!\! (\Lambda\!-\!\nu I_{c2})^2\big[\!(-h_0\!+\!\beta\!+\!\mu)(b\!+\!I_{c2})^2\!-\!b^2(h_1\!-\!h_0)\!\big]\!-\! 2\mu\beta I_{c2}(b\!+\!I_{c2})^2\big[\!\Lambda\!-\!(\mu\!+\!\nu)I_{c2}\!\big] \!\!\bigg\}, \end{aligned}$

where $S_{ic2}$ is the abscissa of boundary equilibrium $E_2^b$. It follows that the persistence happens for $\beta+\mu\geq h_0$ and

$\frac{2\beta(b+I_{c2})^2}{b^2(h_1-h_0)-(b+I_{c2})^2(\beta+\mu-h_0)} <\frac{(\Lambda-\nu I_{c2})^2}{\mu I_{c2}\big[\Lambda-(\mu+\nu)I_{c2}\big]}$

(9)

while a non-smooth fold appears for $\beta+\mu\geq h_0$ and

$\frac{2\beta(b+I_{c2})^2}{b^2(h_1-h_0)-(b+I_{c2})^2(\beta+\mu-h_0)} >\frac{(\Lambda-\nu I_{c2})^2}{\mu I_{c2}\big[\Lambda-(\mu+\nu)I_{c2}\big]}.$

(10)

If we choose the parameters as $\Lambda = 8, \mu = 0.1, \beta = 1.8, h_0 = 0.2, h_1 = 0.8, b = 5, \nu = 0.6, I_c = 8.2934 \mbox{(a)}, I_c = 6.6667 \mbox{(b)}, I_c = 9.2934\mbox{(c)}, I_c = 9.7934\mbox{(d)}$ and keep $I_c$ as a bifurcation parameter, Fig. 1 shows the boundary node bifurcation occurs when the threshold parameter $I_c$ passes through the critical value $I_{c1}$ or $I_{c2}$, where $I_{c1} = 6.6667, I_{c2} = 9.2934, \mbox{det}(J_2(E_1^b)) = 0.011$. In Fig. 1, the red solid circle points (i.e. $E_1^r, E_2^r$) stand for the real endemic equilibria while the red hollow ones (i.e. $E_2^v$) denote the virtual equilibria. The hollow square dots (i.e. $E_0$) represent the unstable disease-free equilibria and the solid diamond dot (i.e. $E_s$) denote the pseudo-equilibria. The stable (unstable) manifolds of the pseudo-saddle $E_s$ are depicted by blue (green) lines. The grey thick solid lines are used to show the repulsing sliding mode region while the grey thin dashed lines are to represent the transversal sliding mode region. All the black solid lines stand for the trajectories of system (4). In the rest of this work, unless otherwise stated, we always use these notations.

The pseudo-saddle $E_s$, stable node $E_1^r$ and visible tangency point $E_1^t$ coexist for $I_c>I_{c1}$ and $I_c = 8.2934$ as shown in Fig. 1(a). They collide together simultaneously when $I_c$ decreases and reaches the critical value $I_{c1}$, denoted by $E_1^b$, as shown in Fig. 1(b). In this case, $E_1^b$ is a pseudo-saddle with an incoming sliding orbit. Further, since (7) holds ture, a non-smooth fold occurs, which indicates the real equilibrium $E_1^r$ disappears after it collides with the pseudo-equilibrium $E_s$. However, as $I_c$ increases and reaches the critical value $I_{c2}$, another boundary node bifurcation occurs, which is due to the collision of real endemic equilibrium $E_2^r$, the visible tangency point $E_2^t$ and the pseudo-equilibrium $E_s$, denoted by $E_2^b$, as shown in Fig. 1(c). In such scenario, $E_2^b$ is also a pseudo-saddle with an incoming sliding orbit. It is worth noting that (10) is satisfied, which implies the occurrence of non-smooth fold in this process. In fact, as $I_c$ increases further, the pseudo-saddle $E_2^b$ is superseded by an invisible tangency point $E_2^t$ and a virtual equilibrium $E_2^v$ for $I_c>I_{c2}$ and $I_c = 9.7934$ as shown in Fig. 1(d). This shows how a catastrophic disappearance of a stable equilibrium occurs.

Note that the endemic equilibrium $E_3$ is a saddle provided it is feasible, so a boundary saddle bifurcation may occur once $E_3$ collides with the tangency point at the critical parameter value $I_c = I_{c3}$ as shown in Fig. 2. The type of unfolding of the boundary equilibrium $E_3^b$ is determined by

$\begin{aligned} &\mbox{sgn}\Big(J_3(S_{c3}, I_{c3})J_2^{-1}(S_{c3}, I_{c3}) \big(X_1(S_{c3}, I_{c3})-X_2(S_{c3}, I_{c3})\big)\Big)\\ = &\mbox{sgn}\bigg\{\!\! (\Lambda\!-\!\nu I_{c3})^2\big[\!(h_0\!-\!\beta\!-\!\mu)(b\!+\!I_{c3})^2\!+\!b^2(h_1\!-\!h_0)\!\big]\!+\! 2\mu\beta I_{c3}(b\!+\!I_{c3})^2\big[\!\Lambda\!-\!(\mu\!+\!\nu)I_{c3}\!\big] \!\!\bigg\}, \end{aligned}$

where $S_{c3}$ represents the abscissa of $E_3^b$. Then a non-smooth fold is observed if (9) is true and persistence is derived if (10) holds. In this scenario, (10) is true and so persistence occurs. In Fig. 2(a), the virtual equilibrium $E_3^v$, stable focus $E_2^r$, visible tangency point $E_2^t$, pseudo-saddle $E_s$ and the disease-free equilibrium $E_0$ coexist for $I_c>I_{c3}$ with $I_{c3} = 1.1463$ and $I_c = 2$. As $I_c$ decreases and passes through the critical value $I_{c3}$, the virtual equilibrium $E_3^v$, visible tangency point $E_2^t$ and pseudo-saddle $E_s$ collide together simultaneously, denoted by $E_3^b$, which results in a boundary saddle bifurcation as shown in Fig. 2(b). As $I_c$ further decreases and $I_c<I_{c3}$, the boundary saddle $E_3^b$ is replaced by the visible tangency point $E_2^t$ and the standard saddle $E_3^r$, as shown in Fig. 2(c) with $I_c = 1$. It is worth mentioning that different monotonicity levels of $I_c$ yield different phenomena if we take $I_c$ as a variable. In fact, if $I_c$ is decreasing, we get the following variation order, i.e. Fig. 2$(a)\rightarrow (b)\rightarrow (c)$, which shows how a pseudo-saddle becomes a standard saddle. If $I_c$ is increasing, we get the change order Fig. 2$(c)\rightarrow (b)\rightarrow (a)$, which indicates the real equilibrium $E_3^r$ turns into the pseudo-equilibrium $E_s$ in this case.

It follows from Theorem 2.2 that a unique endemic equilibrium $E_*$ exists for system $S_{G_2}$ if $R_0<1, a_0>0, a_1<0, C_0 = 0$. In fact, $E_*$ is a saddle node point of system $S_{G_2}$ and is actually the result of the collision of two endemic equilibria $E_2$ and $E_3$. A boundary saddle node bifurcation occurs when $E_*$ collides with the repulsing sliding mode region $\Sigma_{e2}$, according to [26], as shown in Fig. 3(b) with $I_c = I_{c5}$ and $I_{c5} = 2.5661$. Saddle node point $E_*^r$ and invisible tangency point $E_2^t$ coexist for $I_c<I_{c5}$, as shown in Fig. 3(a). The saddle node point $E_*^r$ collides with the invisible tangency point $E_2^t$, which results in occurrence of the boundary saddle node point $E_2^b$ with a sliding stable manifold as shown in Fig. 3(b). If $I_c>I_{c5}$, $E_2^b$ is replaced by invisible tangency point $E_2^t$ and virtual equilibrium $E_*^v$, as shown in Fig. 3(c).

For other cases with the existence of endemic equilibrium $E_2$, $E_4$ or $E_*$, boundary node (or focus) bifurcation occurs once the threshold parameter $I_c$ passes through the critical value $I_{c2}, I_{c4}$ or $I_{c5}$ as shown in Fig. 4 (g)-(h).

---

4.2. Global sliding bifurcation

Now we turn to examine the bifurcations of limit cycles of (4). Basically, there are three types of limit cycles for Filippov system (4) [29], i.e. standard periodic cycles, sliding periodic cycles and crossing periodic cycles. Standard periodic solutions refer to the cycles lying entirely in region $G_1$ or $G_2$; sliding periodic cycles are those cycles having a sliding segment in $\Sigma_e$; crossing periodic cycles represent the periodic solutions having only isolated points in common with $\Sigma_e$. In what follows we address grazing bifurcation and pseudo-homoclinic bifurcations for system (4) by choosing the threshold level $I_c$ as bifurcation parameter and specifying other parameters as shown in Fig. 4.

Grazing bifurcation. The bifurcation that occurs once the standard piece of a periodic cycle touches the switching boundary $\Sigma$ is said to be grazing bifurcation (or touching bifurcation) [29]. It follows from Fig. 4(c) that an unstable standard cycle exists in region $G_2$ with $I_c = 4.95$.In such scenario, there are five critical points, i.e. the real endemic equilibrium $E_2^r$, virtual endemic equilibrium $E_3^v$, stable disease-free equilibrium $E_0^r$, pseudo-equilibrium $E_s$ and visible tangency point $E_2^t$, two stable manifolds of the pseudo-saddle $E_s$, denoted by the blue solid lines, and two unstable manifolds of it, denoted by the green solid lines. As the parameter $I_c$ increases and reaches the value $I_c = 5.03$, grazing bifurcation occurs as shown in Fig. 4(d), where the closed orbit is tangent to the repulsing sliding mode region $\Sigma_{e2}$ at the tangency point $E_2^t$, which is a so-called crossing cycle. As $I_c$ further increases, the crossing cycle becomes a sliding cycle, which contains a piece of sliding segment, as shown in Fig. 4(e) with $I_c = 5.2$.

It is worth noting that local sliding bifurcation happens at the critical threshold value $I_{c3}$, i.e. real/virtual equilibrium bifurcation or boundary saddle bifurcation, before the occurrence of grazing bifurcation, as shown in Fig. 4(a)-(c) with $I_{c3} = 4.7844$. A real endemic equilibrium $E_3^r$ and invisible fold tangency point $E_2^t$ coexist as shown in Fig. 4(a) with $I_c<I_{c3}$ and $I_c = 4$, where $E_3^r$ is a saddle. With the increasing of $I_c$ and $I_c = I_{c3}$, the saddle point $E_3^r$ collides with tangency point $E_2^t$ and becomes boundary saddle $E_3^b$, which results in boundary saddle bifurcation, where $E_3^b$ possesses a stable sliding manifold as shown in Fig. 4(b). If we further increase $I_c$, the boundary saddle disappears and is replaced by the tangency point $E_2^t$, pseudo-saddle $E_s$ and virtual equilibrium $E_3^v$. This process shows how a real equilibrium becomes a virtual equilibrium.

Bifurcation of a sliding homoclinic orbit to pseudo-saddle. If a pseudo-equilibrium of system (4) is a pseudo-saddle, it can have a sliding trajectory which initiates and returns back to it at certain threshold parameter values. This is indeed the so-called sliding homoclinic orbit. If a sliding cycle collides with such a pseudo-saddle, a sliding homoclinic bifurcation occurs as shown in Fig. 4(e)-(g). It follows from Fig. 4(e) that a sliding cycle, two pieces of stable manifolds and unstable manifolds of the pseudo-saddle $E_s$ coexist with $I_c = 5.2$. As $I_c$ increases and reaches about $6.009782$, the trajectory departing from the tangency point $E_2^t$ becomes the stable manifold of the pseudo-saddle and the sliding cycle collides with the pseudo-saddle $E_s$ as shown in Fig. 4(f), where the sliding cycle is replaced by a sliding homoclinic orbit. As $I_c$ further increases, the sliding homoclinic orbit breaks and no periodic orbit exists as shown in Fig. 4(g). If $I_c$ continues to increase and passes through $I_{c2}$ with $I_{c2} = 6.8556$, a local sliding bifurcation (i.e. boundary node bifurcation) occurs as shown in Fig. 4(h).

It follows from Fig. 4 that as the threshold level increases from 4 to 6.8556, Filippov system (4) exhibits the interesting local and global sliding bifurcations sequentially, i.e. boundary saddle bifurcation $\rightarrow$ grazing bifurcation $\rightarrow$ sliding homoclinic bifurcation $\rightarrow$ boundary node bifurcation. This illustrates a series of rich dynamics and suggests that the control outcomes may be sensitive to the threshold value.

---

5. Impact of interventions of Filippov system (4)

One of our purposes in this study is to seek better strategies for curbing the spread of diseases, and to examine when and how to implement the control measure to contain the number of infected individuals to be less than some acceptable level if it is almost impossible to eradicate an infectious disease. An efficient way is to reduce the level of infected cases in steady state as small as possible. To realize this goal, we will examine the impact of parameters associated with interventions on the dynamics, especially the equilibrium level of infected individuals of Filippov system (4) in this section. Due to the biological significance of each parameter and their effect on disease control, we focus only on the impact of threshold parameter $I_c$, HBPR parameter $b$, maximum and minimum per capita treatment rate $h_1$ and $h_0$ in this section.

Impact of threshold level $I_c$. We start from the effect of the variation of the threshold level $I_c$ on the the sliding mode region. Fig. 5 shows the possible sliding modes and pseudo-equilibrium for different values of parameter $I_c$. We select $I_c = 0, 1, 2, \cdots, 20$ in each subplot. The black solid lines represent the vertical isoclines for system $S_{G_1}$ and system $S_{G_2}$, denoted by $g_1^s$ and $g_2^s$, respectively. The grey thin solid lines represent the horizontal isoclines for system $S_{G_1}$ and system $S_{G_2}$, denoted by $g_1^i$ and $g_2^i$, respectively. The grey thick solid lines stand for the sliding mode regions with various threshold levels $I_c$. The red cycle points are the regular endemic equilibria while the red diamond points denote the pseudo-equilibria of the system (4). Subplot (a) and (b) show the case $R_0>1$ and the case $R_0<1$, respectively.

It is clear that the pseudo-equilibria and sliding modes are sensitive to the threshold level $I_c$. If the basic reproduction number is greater than unity, Fig. 5(a) indicates that the sliding mode region enlarges as $I_c$ increases. In particular, if we let $I_c<6.6667$, for instance, $I_c = 6$, the endemic equilibrium of system $S_{G_1}$ (i.e.$E_1$) is virtual, that of system $S_{G_2}$ (i.e. $E_2$) is real and no pseudo-equilibrium exists. If $I_c$ increases and reaches $I_c = 6.6667$, a boundary equilibrium $E_1^b$ appears. If $I_c$ further increases such that $6.6667<I_c<9.2934$, for example $I_c = 7$, a pseudo-equilibrium $E_s$ and two real endemic equilibria appears. If $I_c$ increases to $I_c = 9.2934$, another boundary equilibrium $E_2^b$ is feasible. Note that the pseudo-equilibrium disappears while virtual endemic equilibrium $E_2$ and real endemic equilibrium $E_1$ coexist for $I_c>9.2934$, for example $I_c = 10$. If $R_0<1$, it follows from Fig. 5(b) that the repulsing sliding mode region is feasible only when the threshold level satisfies $I_c>I_{c0}$. Similar argument to that on Fig. 5(a) yields the evolution of sliding mode regions as well as equilibria.

It follows from the above discussion that we can increase the threshold level $I_c$ such that only endemic equilibrium $E_1 = (S_1, I_1)$ of system $S_{G_1}$ is real for $R_0>1$, where $I_1$ represents a relatively low equilibrium level of infection. If $R_0<1$, the infectious disease cannot be eradicated due to the existence of multiple endemic equilibria of system $S_{G_2}$. However, Fig. 5(b) shows that both endemic equilibria of system $S_{G_2}$ can be virtual if $I_c$ is increased to a high level in such scenario, so neither of them can be the attractor of the Filippov system. This indicates that increasing the threshold level $I_c$ can lead to the nonexistence of real endemic equilibria, i.e. only the disease free equilibrium is real in such scenario.

Impact of HBPR parameter $b$. Due to the significance of hospital bed-population ratio (HBPR) on the control of disease, it is essential to examine the impact of HBPR parameter $b$ on disease infection. To this purpose, we numerically conduct the analysis about the effect of $b$ on the variation of sliding mode and equilibria, as shown in Fig. 6. The black circle points represent the endemic equilibria of system $S_{G_2}$ with higher endemicity while the square ones denote the endemic equilibria of system $S_{G_1}$ (shown in Fig. 6(a)) or the endemic equilibrium of system $S_{G_2}$ with lower endemicity (shown in Fig. 6(b)). The solid ones of those are real equilibria while the hollow ones are virtual equilibria. We select the parameter $b$ as a bifurcation parameter and all other parameters are fixed as indicated in Fig. 6, with $R_0>1$ for subplot (a) and $R_0<1$ for subplot (b). Fig. 6 shows that the sliding mode region shrinks as the parameter $b$ increases for both $R_0>1$ and $R_0<1$. For $R_0>1$, Fig. 6(a) indicates that two endemic equilibria ($E_1, E_2$) coexist for system (4), one (i.e. $E_1$) for system $S_{G_1}$ and the other (i.e. $E_2$) for system $S_{G_2}$. Further, $E_1$ is real and independent of the size of parameter $b$ while a real/virtual equilibrium bifurcation occurs at $E_2$, i.e. the real equilibrium $E_2$ becomes a virtual one when the parameter $b$ passes through a critical value around $b = 12$, where the pseudo-equilibrium $E_s$ disappears. In fact, $E_s$ is feasible for $b$ less than the critical value $b = 12$. It follows from Fig. 6(b) that the sliding mode region is feasible only when the parameter $b$ satisfies $b < b_c$ with $b_c = (\beta-r_0)I_c/(r_1-\beta)$, and here we have $b_c = 8$. In such scenario, both endemic equilibrium $E_2$ and $E_3$ exist for the parameter $b$ less than a critical value about $b = 3.4315$. The one with lower endemicity (i.e. $E_3$) is virtual provided it is feasible while a real/virtual equilibrium bifurcation is triggered at $E_2$ when the parameter $b$ passes through a certain critical value around $2.6$, below which the pseudo-equilibrium is feasible.

The above discussion demonstrates that only endemic equilibrium of free system exists for Filippov system (4) by increasing the HBPR (i.e. $b$) for $R_0>1$, which indicates $I_1$ stands for the equilibrium level of infected individuals in this case. If $R_0<1$, a backward bifurcation is triggered for system $S_{G_2}$ when the parameter $b$ passes through the critical value $b = 3.4315$. In particular, if $R_0<1$ and $b>3.4315$, no stable endemic equilibrium exists for system $S_{G_2}$; for $b<3.4315$, there are two endemic equilibria, which collide together and are substituted by a unique endemic equilibrium for $b = 3.4315$. Note that only those real stable equilibria can be the attractor of system (3) and endemic equilibrium $E_2$ is real for $R_0<1, b<2.6$, so the coexistence of stable disease free equilibrium and endemic equilibrium occurs at $b = 2.6$ for system (3). This indicates that no endemic equilibrium exists besides the disease-free equilibrium $E_0$ when the number of hospital beds is enough large, so $E_0$ becomes the unique attractor in such scenario.

Note that here increasing $I_c$ or $b$ makes the Filippov system work as system $S_{G_1}$, and consequently this bifurcation result is rather obvious. Moreover, we refer backward to the co-existence of a stable disease free equilibrium and a stable endemic equilibrium when the basic reproductive number is less than unity, i.e., the backward bifurcation leads to bistability for $R_0<1$ as a parameter varies. In particular, the real endemic equilibrium $E_2$, which can be stable, begins to exist for system (3) besides the disease free equilibrium when the parameter $b$ passes through the value $b = 2.6$, as addressed above, so a backward bifurcation occurs.

Impact of the maximum and minimum treatment rate. In this section, we focus on the impact of maximum and minimum per capita treatment rate (i.e. $h_1$ and $h_0$) on the equilibrium level of infected individuals for system $S_{G_2}$ in the following.

Fig. 7 is to show the variation of the number of infected individuals in steady states as the maximum per capita treatment rate $h_1$ varies. The black solid (dashed) lines represent the number of infected individuals with high (low) endemicity in steady states, i.e. $I_2$ ($I_3$). The grey solid line describes the variation of basic reproduction number $R_0$ with the maximum treatment rate $h_1$. Solving $R_0 = 1$ and $C_0 = 0$ with respect to $h_1$ yields $h_{11} = \beta-(\mu+\nu)$ and

$\begin{aligned} h_{12}& = \frac{4\beta\mu b\Lambda-2\mathcal{D}_1-\sqrt{\big(4\beta\mu b\Lambda-2\mathcal{D}_1\big)^2-4b^2\nu^2\big(\mathcal{D}_1^2+4\beta\mu br_0\Lambda\big)}}{2b^2\nu^2}-(\mu+\nu), \\ \mathcal{D}_1& = b\nu\big[\Lambda(r_0-\beta)-\beta b(\mu+\nu)\big], \end{aligned}$

respectively. Denote

$\begin{aligned} L_1& = \big\{(S, I)\in R^2_+:\ S = h_{11}\big\}, L_2 = \big\{(S, I)\in R^2_+:\ S = h_{12}\big\}. \end{aligned}$

Then $L_1$ divides the whole domain into two parts, the one to the left of which is for $R_0>1$ while the one to the right is for $R_0<1$. Line $L_2$ also divides the whole region into two parts, the part to the left of which stands for $C_0>0, (-a_1+\sqrt{C_0})/(2a_0)<\Lambda/(\mu+\nu)$ while the part to the right of which is for $C_0\leq0$ or $(-a_1+\sqrt{C_0})/(2a_0)>\Lambda/(\mu+\nu))$. So, there is a unique endemic equilibrium $E_2$ for subsystem $S_{G_2}$ if $h_1<h_{11}$ (i.e. $R_0>1, C_0>0$), which is stable in some cases; two endemic equilibria $E_2$ and $E_3$ coexist for subsystem $S_{G_2}$ if $h_{11} < h_1<h_{12}$ (i.e. $R_0<1, C_0>0, (-a_1+\sqrt{C_0})/(2a_0)<\Lambda/(\mu+\nu)$), and $E_2$ is stable in some cases while $E_3$ is unstable; no endemic equilibrium is feasible for $h_1>h_{12}$ (i.e. $R_0<1$ and $C_0<0$ or $(-a_1+\sqrt{C_0})/(2a_0)<\Lambda/(\mu+\nu)$); a unique endemic equilibrium $E_*$ exists for $h_1 = h_{12}$ (i.e. $R_0<1, C_0 = 0$). Further discussion yields that a backward bifurcation occurs when the parameter $h_1$ passes through the critical value $h_{12}$ on one hand. On the other hand, a saddle-node bifurcation occurs at $h_1 = h_{12}$, which implies $E_*$ is a saddle-node as the result of the collision of $E_2$ and $E_3$.

It follows from Fig. 7 that $I_2$ decreases as $h_1$ increases for $h_1 < h_{12}$, which implies that the number of infected cases decreases as the maximum treatment rate increases. If we select the maximum treatment rate properly such that $h_1>h_{12}$, only the disease free equilibrium $E_0$ is the attractor.

By choosing the minimum per capita treatment rate $h_0$ as a bifurcation parameter and fixing all other parameters, the evolution of the number of infected cases is clearly depicted in Fig. 8. Solving $C_0 = 0$ in terms of $h_0$ gives

$\begin{aligned} h_{01}& = \frac{2\mathcal{D}_2\Lambda-4\beta\mu b\Lambda-\sqrt{\big(4\beta\mu b\Lambda-2\mathcal{D}_2\Lambda\big)^2-4\Lambda^2\big(\mathcal{D}_2^2-4\beta\mu b\Lambda r_1\big)}}{2\Lambda^2}-(\mu+\nu), \\ \mathcal{D}_2& = \beta\Lambda-b\nu r_1+b\beta(\mu+\nu). \end{aligned}$

Denote

$\begin{aligned} L_4& = \big\{(S, I)\in R^2_+:\ S = h_{01}\big\}. \end{aligned}$

According to Fig. 8(a), the equilibrium number of infected cases for system $S_{G_2}$ can be contained at a low level by enhancing the minimum per capita treatment rate $h_0$ for $R_0>1$. Similar argument yields that increasing the minimum treatment rate can lead to the existence of only disease free equilibrium $E_0$ for system $S_{G_2}$.

---

6. Discussion and conclusions

It has been observed that treatment programme plays a significant role in controlling the emerging and reemerging infectious diseases, such as HIV [36] and Ebola [23], and usually the control programme is implemented only when the number of infected individuals reaches or exceeds the threshold level $I_c$. In this study, we propose a Filippov epidemic model by using a piecewise smooth function to reflect the treatment rate, which incorporates the impact of number of hospital beds. By taking the number of infected cases as a critical index, we define a threshold policy as follows: the treatment policy, depending on HBPR and the number of infected cases, is implemented when the number of infected individuals is above the threshold level; below this level, a larger treatment ratio, which is indeed the maximum per capita treatment ratio the health care system can provide, is adopted. Our modeling and results extend the recent work [35] on the impact of threshold policy on the control of infectious diseases.

Based on the dynamics of two subsystems of the Filippov system, we investigate the long-term dynamical behavior, which reveals much more complex dynamics compared to those for the continuous counterpart. In particular, we have examined the sliding mode, pseudo-equilibrium, multiple attractors in Section 3 and Section 4, respectively. As the threshold value varies, bifurcation analysis of piecewise smooth systems [17,29] on the proposed Filippov epidemic model yields the following local sliding bifurcations theoretically, i.e. regular/virtual equilibrium bifurcation, boundary equilibrium bifurcation including BN (see Fig. 1), BF, BS (see Fig. 2) and BSN bifurcations (see Fig. 3), and global sliding bifurcations including grazing bifurcation (see Fig. 4) and sliding homoclinic bifurcation (see Fig. 4). Our results demonstrate that variety of threshold level can give rise to diversity of long-term dynamical behavior.

It is worth emphasizing that understanding of impact of threshold policy will lead to the development of effective control programmes for public health, so we have examined the impact of some key parameters related to the control measure in Section 5. According to Fig. 5(a), if the basic reproduction number is greater than unity, we can choose proper threshold level such that the real equilibrium of system $S_{G_1}$, or that of system $S_{G_2}$, or both of them together with the unique pseudo-equilibrium act as the attractors for Filippov system (4), which implies that different threshold levels and initial states lead to different levels of total number of infected cases. It follows from Fig. 5 and Fig. 6 that we can select a proper threshold value or the number of hospital beds such that the endemic equilibrium with the lowest endemicity or the disease free equilibrium appears as the attractor (i.e. the endemic one of system $S_{G_1}$ for $R_0>1$ while disease free equilibrium for $R_0<1$). This observation indicates that the course of an outbreak as well as the control outcome are sensitive to the threshold level and the available HBPR. If the infectious disease is on the course of becoming endemic ($R_0>1$), it is essential to increase the medical resources, i.e, to increase the number of hospital beds therefore the possibility of implementing maximum treatment ratio to control the infection to a relatively low level. For the case when $R_0<1$, the infection can still become endemic due to the possible existence of backward bifurcation of system $S_{G_2}$, increasing the medical resources therefore be essential to control the spreading of infection.

In most cases, it is impossible to eradicate an epidemic disease from the population. We have investigated how the maximum and minimum treatment rate affect the equilibrium level of infected individuals in Section 5. Fig. 7 and Fig. 8 show that the number of infected individuals at the steady state is closely related to the maximum and minimum per capita treatment rate. The maximum and minimum per capita treatment rate remain effective for reducing transmission during the outbreak. In particular, if $R_0>1$, strengthening the basic medical conditions and increasing medical resources, i.e. increasing the minimum and maximum treatment ratio, will result in a relative low level of infected cases, as shown in Fig. 8 (a) and Fig. 7. Increasing maximum treatment rate $h_1$ (or minimum treatment rate $h_0$ if $R_0<1$) such that $h_1>h_{12}$ (or $h_0>h_{01}$) will aid in eradicating the disease, as shown in Fig. 7 and Fig. 8(b).

It is worth noting that we choose the size of infected individuals as an index. In fact, the number of susceptible individuals can also affect on the implementation of control measures, especially on those immunization policies, so it is more appropriate for modeling the impact of limited resources on vaccination policies. Moreover, it is more natural to choose the total number of susceptible individuals and infected individuals as an index, especially for the combined control measure based on treatment and vaccination policy, which may be difficult and interesting, and we leave this for future work.

In conclusion, we propose a Filippov epidemic model to study the impact of HBPR on the disease control by incorporating a piecewise defined treatment programme in this paper. The main results illustrate the significant role of switching treatment programme in response to the call of controlling the epidemic disease. Indeed, to illustrate the main idea, we formulate a simple SIS model without considering the effect of hospitalized infection. Generally, those hospitalized has no contact with susceptible individuals, so no transmission occurs. An SIH (susceptible-infective-hospitalized) model of three dimension could then be more natural, which we will consider in the future work.

---

Acknowledgments

The first author was supported by Scientific research plan projects of Shaanxi Education Department (16JK1047). Xiao was supported by the National Natural Science Foundation of China (NSFC, 11571273 and 11631012) and Fundamental Research Funds for the Central Universities(GK 08143042). Zhu was supported by NSERC of Canada. This research was finished when Wang visited Lamps and Department of Mathematics and Statistics, York University.

---