Joint impacts of media, vaccination and treatment on an epidemic Filippov model with application to COVID-19.

A non-smooth SIR Filippov system is proposed to investigate the impacts of three control strategies (media coverage, vaccination and treatment) on the spread of an infectious disease. We synthetically consider both the number of infected population and its changing rate as the switching condition to implement the curing measures. By using the properties of the Lambert W function, we convert the proposed switching condition to a threshold value related to the susceptible population. The classical epidemic model involving media coverage, linear functions describing injecting vaccine and treatment strategies is examined when the susceptible population exceeds the threshold value. In addition, we consider another SIR model accompanied with the vaccination and treatment strategies represented by saturation functions when the susceptible population is smaller than the threshold value. The dynamics of these two subsystems and the sliding domain are discussed in detail. Four types of local sliding bifurcation are investigated, including boundary focus, boundary node, boundary saddle and boundary saddle-node bifurcations. In the meantime, the global bifurcation involving the appearance of limit cycles is examined, including touching bifurcation, homoclinic bifurcation to the pseudo-saddle and crossing bifurcation. Furthermore, the influence of some key parameters related to the three treatment strategies is explored. We also validate our model by the epidemic data sets of A/H1N1 and COVID-19, which can be employed to reveal the effects of media report and existing strategy related to the control of emerging infectious diseases on the variations of confirmed cases.

Introduction

The infectious diseases have been a worldwide threat to public health for many years. The emerging disease outbreaks such as sever acute respiratory syndrome (SARS) in 2003 (Seto et al., 2003, He and Peng, 2004), the A/H1N1 influenza epidemic in 2009 (Khan and Arino, 2009, Jones and Salathe, 2009, Tang et al., 2010), the Ebola virus disease in 2014 (Aylward and Barboza, 2014, Baize and Pannetier, 2014) and the COVID-19 in 2019 (Huang and Wang, 2020, Tang et al., 2020, Chintalapudi et al., 2020) have made strong impacts on human health, global economy, and social behaviors. It is thus extremely important to take timely actions and implement the treatment strategies to relieve and eliminate the risk of these diseases. During the spread of an infectious disease, there are three main therapeutic schedules: reducing the contact rate of the susceptible population with the infected population; injecting vaccine for the susceptible population; and implementing the recovery treatment for the infected population as early as possible (Cui et al., 2008, Tchuenche et al., 2011, Wang and Xiao, 2013, Zhou and Fan, 2012, Huang and Li, 2019).

Most infectious diseases can spread through human contact. Hence, reducing contact rate with the infected individuals will help to prevent the people from getting infected and reduce the risk of epidemic outbreaks. If the disease can spread via air droplets, implementing stay-at-home policy and enforcing the public to wear facial masks in the crowd could be very effective to reduce the rate of infection. In general, if one obtains the latest disease information from the media such as TV or internet, and learns that the outbreak is very serious, she/he may strengthen the self-protection awareness and take actions to reduce the risk of infection. Thus, media coverage plays a significant role in containing the spread of infectious diseases. This motivates us to investigate the impact of the media coverage on the disease transmission. In recent years, many researchers have studied the media coverage impact by establishing and analyzing the corresponding mathematical models (Cui et al., 2008, Cui et al., 2017, Chen et al., 2018, Xiao et al., 2013, Wang and Xiao, 2014, Tchuenche et al., 2011). To describe the reduction of contact rate, the term in Tchuenche et al. (2011) or the saturation incidence rate in Cui et al. (2017) are formulated to investigate the media impact on the epidemic spread, while the exponentially decay functions (Cui et al., 2008), (Wang and Xiao, 2014) and (Liu et al., 2007) (with denoting the exposed, infectious and hospitalized individuals respectively) are proposed to characterize the transmission coefficient. The contact terms presented above are all related to the number of infected population I, where the transmission intensity will be weaken as I increases. In our model, we will adopt the novel exponential decreasing factor as the incidence rate which depends on both I and the changing rate . This term was first proposed in a study of the media effect on the dynamics of a non-smooth system (Xiao et al., 2013). Note that enhancing the media effect will lead to a lower transmission rate during the outbreak of an infectious disease.

In addition to the media coverage, vaccination also plays an important role and has been an effective method to control many infectious diseases such as smallpox, measles and parotitis. There is a rich literature in the study of the impact of vaccine on the spread of infectious diseases via various mathematical models (Buonomo et al., 2008, Zhang et al., 2018, Wang and Xiao, 2013). It has been demonstrated that sufficient vaccine doses will result in disease elimination. However, it takes time to invent vaccine for an emerging infectious disease, and vaccine may not be available during the beginning period of the outbreak; see for example, SARS and A/H1N1. Therefore, it is important to seek for recovery treatment of the infected population, which, in many cases, is the major strategy for controlling the disease. In general, the level of such treatment measure is constrained by the limited medical resources, including medicines, medical workers and hospital beds. The lack of medical resources is a critical issue in underdeveloped regions such as remote rural areas of developing countries, where the medical facilities are insufficient. Hence, optimal allocation of these limited resources is a significant problem for curbing the epidemic spread. In the past decades, there are many studies on the effectiveness of taking recovery treatments (Wang, 2006, Cui et al., 2008, Huang and Li, 2019, Zhou and Fan, 2012, Li et al., 2017, Zhang and Liu, 2008), and it has been shown that the treatment rate indeed relies on the medical resources available for the public health. To describe the limitation of medical resources, Wang (2006) has introduced a piecewise treatment functionwhich assumes that the number of treated patients is proportional to the number of infected population before the capacity of medical resource is reached, and a constant thereafter. The work in Cui et al. (2008) put forward the recovery formwhere denotes the maximum recovery rate, and is a half-saturation constant. Such function is used to detect the impact of limited medical resources on the epidemic controlling (Zhou and Fan, 2012). In Zhang and Liu (2008), a saturation treatment termis proposed to represent the delayed efficiency of the infected cases for being cured. Bifurcation and stability analyses for specific epidemic models involving such a treatment term are explored in Huang and Li, 2019, Li et al., 2017. In this paper, we will employ the same treatment function to describe the limited medical resources.

Hybrid dynamic systems including impulsive system, non-smooth continuous system and Filippov system have been constructed to investigate the impacts of various treatment strategies on the dynamical behaviour of the epidemic models; see Zhang et al., 2018, Xiao et al., 2013, Chen et al., 2018, Mu et al., 2019, Wang and Xiao, 2014. In this paper, we will propose a Filippov SIR model to incorporate the aforementioned three control strategies, analyze the model dynamics and explore the biological significance. Our main purpose is to maintain the infected population at a desired level by choosing appropriate implementation time and strength for the above three protection measures. We will find effective switching containment strategies for the infectious diseases.

The rest of our paper is organized as follows. In Section 2, we propose an epidemic Filippov model with two subsystems and present some preliminaries related to the proposed system. In Section 3, we explore the dynamics of two subsystems in detail. In Section 4, we examine the dynamics on the sliding domain of the Filippov system. In Section 5, we conduct local and global sliding bifurcation analysis, and investigate the global stability of the proposed model. In Section 6, some key parameters involving the control measures are chosen to examine the effect of threshold policy on the disease transmission. Our model is also validated by the epidemic data of A/HAN1 and COVID-19 outbreaks. The conclusion and discussion on the biological significance of the proposed model will be given in the last section.

A Filippov SIR model with media coverage, vaccination and treatment

Filippov epidemic model is widely used to describe when and how to take control measures to inhibit disease transmission. In many cases, the prevention measures are implemented once the size of infected population exceeds a certain critical value (Chen et al., 2018, Mu et al., 2019, Wang and Xiao, 2014). However, during a disease outbreak, not only does the number of infected population play an important role, but also the changing rate of infected population is also a significant factor. A small size of infected population with a rapid growth rate presents a strong threat to public health and should trigger a serious reaction. It is thus important to consider both factors in the switching condition of the Filippov model (Xiao et al., 2013). We propose the threshold value for our Filippov system; namely, the control measures for the infectious diseases involving media coverage, the maximum vaccination rate for susceptible individuals and the maximum recovery treatment rate for infected population are implemented once , otherwise, the natural state with only vaccination and treatment presented by the saturation functions is introduced; see more details in the following content.

First, we divide the population into three classes: susceptible (S), infected (I) and recovered (R). To consider the joint impacts of media coverage, vaccination and treatment, we propose the following Filippov epidemic systemwithandwhere represents the recruitment rate, is the basic transmission rate, denotes the natural death rate, stands for the disease-induced death rate. The saturation function represents the vaccination treatment, denotes the maximal vaccination rate and describes the effect of the delay when injecting the vaccine. Moreover, represents the recovery treatment term, denotes the maximal recovery rate and describes the effect of medical resource limitation on the treatment.

In addition, is the changing rate of infected population for subsystem (see more details about the definitions of subsystems in later discussion). The term indicates the media coverage may reduce the contact rate, and the non-negative constants and characterize the dependence of media impact on the size of infected population and its changing rate, respectively. The implicit equation defines the switching manifold (denoted by ) for the Filippov model (1) with (2). In general, this equation is difficult to solve. Using a similar method as in Xiao et al. (2013), we rewrite asin terms of Lambert W function (Corless et al., 1996), here . Now, the switching condition is equivalent to the switching lineMoreover, it can be shown that () is equivalent to (). To ensure the positivity of , we assume throughout this paper. Note that the recovered population does not affect the variation of susceptible and infected, then system (1) can be transformed into the following equationswithwhere and are defined in Eqs. (4), (5). System (6) with (7) describes a SIR Filippov model where the vaccination and recovery treatment strategies are saturation functions when and linear functions when .

For convenience, we denote and defineHence system (6) with (7) can be rewritten as the following Filippov system (Filippov, 1988)Define . Then the switching manifold becomes and it divides into two separate regionsFor each , we denote by the subsystem of (8) in region . It is easily seen that is an attractive region of system (8). We adopt the following definitions of three types of equilibria and the tangent point of Filippov systems (Kuznetsov et al., 2003, di Bernardo et al., 2008).

An equilibrium of Filippov system (8) is called real if it is an equilibrium of subsystem () and lies in the region ; i.e., a real equilibrium of () satisfies (). On the other hand, an equilibrium of subsystem which belongs to another region is defined as virtual; i.e., a virtual equilibrium of () satisfies (). We use and to denote the real equilibrium and virtual equilibrium, respectively.

An equilibrium of Filippov system (8) is defined as a boundary equilibrium if it is an equilibrium of subsystem and lies in the switching boundary ; i.e., it satisfies (or ).

Denotewhere denotes the standard scalar product, and (with ) is the Lie derivative of H with respect to the vector field at the point Z. Hence the sliding domain of system (8) can be defined as , while the crossing domain is . Moreover, we define the sliding mode differential equation on aswith

An equilibrium of Filippov system (8) is called a pseudo-equilibrium if it is an equilibrium on the sliding domain; i.e., it satisfies .

A point of Filippov system (8) is called a tangent point if the orbit initiating from this point is tangent to ; i.e., it satisfies . Further, the tangent point is visible if the orbit of () initiating from the tangent point belongs to region for all sufficiently small neighbourhood of , otherwise it is called invisible.

Dynamics of two subsystems

In this section, we will examine the dynamical behaviour of two subsystems , including all possible equilibria and their stabilities. For subsystem , only one disease-free equilibrium exists, denoted by , whereAccording to Definition 1, is a real (virtual) equilibrium if (). The basic reproduction number of can be defined as and the Jacobian matrix of linearized about an equilibrium isDenote

For subsystem , the disease-free equilibrium is a locally asymptotically stable node for and a saddle for .

The Jacobian matrix at isIt is easy to calculate that . If , we have which indicates that is a saddle. On the other hand, since and for , and , it follows that is a stable node.

For subsystem , an endemic equilibrium satisfies the following equationswhereSince there exist at most three positive roots for the first equation in (12), there are at most three endemic equilibria for . Denote , and let be the three roots (counting multiplicity) of . It follows from Vieta theorem thatIf , then the derivative function has two distinct real rootsDenoteOn account of Cardano’s formula (Cardano, 1968), there are three cases for the roots of : (A) three distinct real roots when ; (B) three real roots and at least two of them are equal when ; (C) one real root and two imaginary roots when . We will investigate these three cases respectively.

Case (A): .

In order to investigate the sign of three distinct real roots of , the parameters sets are established to explore the positivity of (), which yields four subcases: (A1) one positive root and two negative roots or three positive roots for ; (A2) one positive root and two negative roots for ; (A3) two positive roots and one negative root for ; (A4) two positive roots and one negative root or three negative roots for . Denote , then two roots of may appear withwhere and are presented in (13). What we concern is the appearance of positive roots for . Hence, three cases will be focused: (i) one positive real root; (ii) two positive real roots; (iii) three positive real roots. Note that two subcases (A1) and (A2) involve one positive real root, and the difference between the occurrence of one positive root and three positive roots in (A1) is the sign of , as shown in Fig. 1 (a)–(b). By elaborate analyses for the signs of and , we can conclude that the inequalities (18)–(20) are set to ensure the existence of one positive real root. Two positive real roots exist under subcases (A3) and (A4), and the difference between the existence of two positive roots and none positive root in (A4) is the sign of depicted in Fig. 1(c)–(d), furthermore inequalities (21)–(23) can assure case (ii) hold. Only in subcase (A1) may three positive real roots occur, which demands and . Thus, the following inequalities can be obtained.Hence, there exist at most three endemic equilibria for (denoted by , if two endemic equilibria appear, then we denote them as and , further is set for only one endemic equilibrium existing), moreover is real (virtual) if (). Consider the stability of endemic equilibrium , which is locally asymptotically stable if , and a stable node for while a stable focus for .

Fig. 1

Eight categories illustrate the appearance of the positive real roots for in Cases (A)–(C).

Case (B): .

It can be found that two real roots (one of them is a double root) or one real root (which is a triple root) of will occur for , in addition, the double root ( or ) and triple root () are presented in Fig. 1(e)–(g). Only and in (16) can assure the triple root exist for , otherwise a double root occurs. Thus the positive triple root happens only when the formula (24) and () hold true. For the double root, denote the corresponding endemic equilibrium as (the collision of and ) or (the collision of and ), while is set for the triple root.

Case (C): .

Only one real root of arises in such case, and two types of one positive root existing are depicted in Fig. 1(h). For convenience to discuss, denote both positive roots of these two curves as . Then it is crucial to analyse the sign of , and we can yield that for while for . Hence, one positive equilibrium occurs if . Based on the elaborate discussion in the above three Cases (A)-(C), the existence of the endemic equilibria for subsystem are addressed in Table 1 .

Table 1

Conditions for the number of endemic equilibria for .

	One	Two	Three
N<0	Eq. (18) or (19) or (20)	Eq. (21) or (22) or (23)	Eq. (24)
N=0	Eq. (24), n_i=0(i=1,2)	Eq. (24), n_1≠0 or n_2≠0	None
	Eq. (21) or (22) or (23), f(I^′′)=0
N>0	β>βb_0′′	None	None

The subsystem readswhere . The disease-free equilibrium of is with . The basic reproduction number of is defined as . It is readily seen that is stable if and unstable if . In addition, is a real (virtual) equilibrium if (). The endemic equilibrium of (if exists) is denoted by , whereNote that is a necessary condition for the existence of , and (if exists) is always a real equilibrium because . The following two lemmas are coming from Xiao et al. (2013).

No limit cycle occurs for subsystem .

For subsystem , the disease-free equilibrium is globally asymptotically stable if , while the endemic equilibrium is globally asymptotically stable if .

The aforementioned discussions reveal the complex dynamical behaviour of two subsystems in Filippov system (8). Since , the two disease-free equilibria and can not be both real. Note that the endemic equilibrium is real once it exists. Hence, we only consider the real and virtual equilibrium bifurcations for . Here, we choose and as the bifurcation parameters and assume to investigate how different parameter values of () affect existence of the equilibria for . See more details in Fig. 2 and Table 2 .

Fig. 2

(a) Real and virtual equilibrium bifurcations for subsystem . Here, we select and as the bifurcation parameters and fix all other parameters as follows: . (b) The magnification for the interval and in Fig. 2(a). (c)–(d) Two categories for the occurrence of and .

Table 2

Illustration of real and virtual equilibrium characteristics in Fig. 2(b).

Region (1)	E10v and E11r coexist	Region (2)	E10v and E11v coexist
Region (3)	E10v and E1ir(i=1,2,3) coexist	Region (4)	E1iv (i=0,3) and E1jr(j=1,2) coexist
Region (5)	E1iv (i=0,2,3) and E11r coexist	Region (6)	E1iv (i=0,1,2,3) coexist
Region (7)	E1ir (i=0,1,2) coexist	Region (8)	E10v and E1ir (i=1,2) coexist
Region (9)	E1iv (i=0,2) and E11r coexist	Region (10)	E1iv (i=0,1,2) coexist
Region (11)	E10r exists	Region (12)	E10v exists

Based on the dynamical analysis for subsystem , there are supposed to exist fourteen cases to explain the real and virtual equilibrium bifurcations for . The two cases not presenting in Table 2 are: (i) () coexist; and (ii) () coexist. Consider the appearance of and without other endemic equilibrium of , it can be classified into two kinds as shown in Fig. 2(c)–(d) by exploring the null-isoclines of . The star point denotes , the solid circle point stands for , and the black (dashed) curve is the vertical (horizonal) isocline of .

Denote the intersection point of horizonal isocline and S-axis by with . Note that () are irrelevant to . Thus, we can examine the real or virtual properties of by varying . Letting , then the maximum switching value follows for , i.e., . Hence holds under the case described in Fig. 2(c), which implies is virtual invariably. For another category in Fig. 2(d), where the endemic equilibrium is discussed above in Case (B) and . Then the two equilibria () can be both real when the value belongs to the interval (). Consequently, the case () coexisting occurs only when the parameters satisfy the condition in Fig. 2(d). Consider the second case () coexisting, a similar analysis shows that . Therefore, the four equilibria (three endemic equilibria and one disease-free equilibrium) of subsystem can not be all real.

Sliding domain and pseudo-equilibrium

By Definition 4, the tangent point of two subsystems for system (8) satisfiesThen we have the tangent points of () as , whereNote that . The condition guarantees the existence and positivity of two tangent points . The sliding domain satisfies andIt is easily seen that is equivalent to , which corresponds to the sliding domainThe sliding mode differential equation followsSubstituting the formulas of and into Eq. (30) and letting , we obtainwhere . Since , the equation has a positive rootHence the pseudo-equilibrium occurs if . It can be demonstrated that for while for , which implies that is unstable. For the two tangent points in (28), we haveIt follows from Eq. (5) that , which implies due to . Hence, exists only if . Note that is equal to the sign of , where is the horizonal isocline of . Furthermore, such curve is monotonic decreasing with S in the first quadrant. Thus will hold only if the tangent point lies above the curve . According to (27), lies in the curve (the vertical null-isoclinic curve of ), which is also monotonic decreasing with S for . In addition, noting that the switching value possesses the maximum size for , we obtain the following existence criterion of by examining two null-isoclinic curves of .

Consider the number of endemic equilibria for subsystem , the existence criterion of the pseudo-equilibrium described with the threshold value for Filippov system (8) can be divided into the following seven categories.

[] for one endemic equilibrium existing;

[] for two endemic equilibria and existing;

[] or for three endemic equilibria ( ) existing;

[] for the endemic equilibria and existing;

[ ] or for the endemic equilibria and existing;

[] for the endemic equilibrium existing;

[ ] None pseudo-equilibrium happens for no one endemic equilibrium or only the endemic equilibrium existing.

Sliding bifurcation

The bifurcation involving sliding domain is a significant property for Filippov systems, and it can produce much complex dynamical behaviour. In this section, we focus on the sliding bifurcation of system (8), including local and global sliding bifurcations. Note that is related to the threshold and the media coverage in the proposed model, here we choose as the bifurcation parameter and let to examine the joint impact of such two parameters on the dynamics of (8).

Local sliding bifurcation

Note that the local sliding bifurcation includes two types of bifurcations: boundary equilibrium and tangency bifurcations (Kuznetsov et al., 2003, di Bernardo et al., 2008). The tangency bifurcation can occur only when the two tangent points () collide. Hence, no tangency bifurcation occurs for system (8) due to in (28). Consider the boundary equilibrium bifurcation, which contains four different types, i.e., boundary focus, boundary node, boundary saddle and boundary saddle-node bifurcations. On account of the variety of endemic equilibria for subsystem , we mainly focus on the boundary equilibrium bifurcation with respect to (). It is similar to examine the bifurcation properties of and we omit the detail here.

Boundary focus (node) bifurcation: First, the boundary focus bifurcation related to the endemic equilibrium is illustrated in Fig. 3 , where the black dotted line denotes the switching manifold and the thick line stands for the sliding domain . Further, the red (blue) solid circle point is the tangent point of (), the red (blue) star point is the endemic equilibrium of () and the black solid circle point is the pseudo-equilibrium . Without loss of generality, we will employ these notations in the rest analysis of this paper. Note that in Fig. 3(a), two real stable endemic equilibria and coexist which induce the bistable properties, accompanied with two tangent points and an unstable pseudo-equilibrium . Choosing as and collide at one point, denoted by , which implies a boundary focus bifurcation occurring as shown in Fig. 3(b). Moreover, is a pseudo-focus with incoming orbits, and becomes invisible. In this case, the orbits reaching except all tend to (no displaying in Fig. 3(b)–(c)) ultimately.

Fig. 3

Boundary focus bifurcation for Filippov system (8), here we select as the bifurcation parameter and fix all other parameters as follows: .

It is observed in Fig. 3(c) that disappears and turns into virtual (). Similarly, those orbits reaching the sliding domain eventually run to . It can be shown that there exists only one endemic equilibrium () which is globally asymptotically stable for and the disease-free equilibrium is virtual on the parameters space in Fig. 3. Hence is globally asymptotically stable for Filippov system (8) in Fig. 3(c). The variation in Fig. 3 shows how a catastrophic disappearance of a stable equilibrium happens. Fig. 4 expresses a boundary node bifurcation, whose properties is similar to the boundary focus bifurcation. Furthermore, denotes the boundary node bifurcation point induced by the collision of and the real node point as shown in Fig. 4(b).

Fig. 4

Boundary node bifurcation for Filippov system (8), here we select as the bifurcation parameter and fix all other parameters as follows: .

Boundary saddle bifurcation: In Fig. 5 (a), a real stable saddle exists and the green thick line denotes the stable manifold of while the magenta one describes the unstable manifold. Note that there exists no pseudo-equilibrium. Hence, the orbits reaching the sliding domain finally run to . As increases to the critical value in Fig. 5(b), and collide which results in a boundary saddle bifurcation, denoted by . Furthermore, the boundary saddle disappears and occurs; see in the upper right graph of Fig. 5(c). Meanwhile, the invisible tangent point turns into visible and the real endemic equilibrium becomes virtual (). In such parameters space, possesses one virtual disease-free equilibrium () and two endemic equilibria ( and ). Moreover, it can be demonstrated that is unstable for and there is no limit cycle for . Consequently, attracts all orbits except for those initiating from the stable manifold of ; see Fig. 5.

Fig. 5

Boundary saddle bifurcation for Filippov system (8), where we select as the bifurcation parameter and fix all other parameters as follows: .

Boundary saddle-node bifurcation: It follows from Case (B) in Section 3 that a saddle-node for can occur for the collision of saddle and node points happening. Denote the real saddle-node as presented in Fig. 6 (a), and the definitions of green and magenta thick lines are similar to in Fig. 5(a). Note that and collide in Fig. 6(b), which induces a boundary saddle-node () bifurcation occurring, and converts into a virtual equilibrium when . There exists a real stable equilibrium which is not depicted in Fig. 6, in such case, is globally asymptotically stable for the proposed system except the orbits initiating from the stable manifold of in Fig. 6(a); is globally asymptotically stable except in Fig. 6(b); is globally asymptotically stable in Fig. 6(c).

Fig. 6

Boundary saddle-node bifurcation for Filippov system (8), where we select as the bifurcation parameter and fix all other parameters as follows: .

Global sliding bifurcation

The global sliding bifurcation occurs accompanied with the variation of non-vanishing cycles. Hence, it is significant to examine the bifurcation involving limit cycles of model (8). In general, there are three types of periodic solutions (i.e., limit cycles) for Filipppov systems (Kuznetsov et al., 2003), including standard, sliding and crossing periodic solutions. Further, the standard periodic solution denotes the periodic solution which entirely lies in region or ; the sliding periodic solution refers to such periodic solution owing a piece of sliding segment in and the periodic solution passing through the switching boundary with only isolated intersection points in is defined as crossing periodic solution. Accordingly, such three kinds of periodic solutions are called standard, sliding and crossing limit cycles. In what follows we will address grazing (i.e., touching) bifurcation and sliding homoclinic bifurcation to a pseudo-saddle for (8), where we still employ as the bifurcation parameter.

Grazing bifurcation: It follows that once the standard limit cycle touches the sliding domain, then it is called grazing bifurcation. Note that a standard stable limit cycle (red orbit) of , the unstable real endemic equilibrium and the stable real endemic equilibrium coexist in Fig. 7 (a). The magnification near the sliding domain can be seen clearly in the upper right graph, where a pseudo-saddle emerges. For , the stable limit cycle becomes tangent to on the tangent point seen from the upper right in Fig. 7(b), then a grazing bifurcation occurs. Furthermore, the orbits reaching the sliding domain above (below) will run to the limit cycle ().

Fig. 7

Global bifurcation for Filippov system (8), where we select as the bifurcation parameter and fix all other parameters as follows: .

Sliding homoclinic bifurcation to a pseudo-saddle: If there exists a pseudo-saddle , then the orbit initiating from and returning to is called a pseudo–homoclinic. In other words, the sliding pseudo–homoclinic bifurcation can exist. When reaches to 0.504 in Fig. 7(c), there occurs a sliding pseudo–homoclinic orbit with which implies a sliding homoclinic bifurcation to a pseudo-saddle arising.

The global sliding bifurcation in Fig. 7 illustrates the variation of a stable limit cycle of to a stable sliding cycle, and finally to a homoclinic orbit involving a pseudo-saddle, from which we can see the bistable properties of model (8). A stable limit cycle (or stable sliding cycle) and a real stable endemic equilibrium can coexist. The proper parameters will be chosen to investigate the global sliding bifurcation involving crossing limit cycle, and we still select as the bifurcation parameter.

Crossing bifurcation: A stable crossing limit cycle (blue orbit) appears in Fig. 8 (a), accompanied with an unstable real equilibrium , and the details near the sliding domain are magnified in Fig. 8(h). It is not hard to demonstrate that there exist an endemic equilibrium and an unstable disease-free equilibrium for subsystem in such parameters space. Hence the orbits initiating from the external region of the crossing cycle will all tend to such cycle. On the other hand, three kinds of orbits initiating form region in the interior space of the crossing cycle are exhibited here. It can be seen from Fig. 8(h) that the mid orbit traverses the switching line and runs into , then reaches the sliding domain below and tends to at last. However, the two other orbits run to the crossing cycle eventually. In conclusion, for the orbits initiating from the interior space of the crossing cycle, once they transverse the switching line, then such cycle is the limit set if the orbits reach above the pseudo-equilibrium, while other orbits all tend to .

Fig. 8

Global bifurcation for Filippov system (8), here we select as the bifurcation parameter and fix all other parameters as follows: .

As increases to , a stable crossing limit cycle with an isolated point on occurs as shown in Fig. 8(b), and the stable cycle involves a piece of sliding segment when . The orbits depicted in Fig. 8(a)–(c) illustrate the crossing bifurcation for system (8). Enlarging to 0.141 in Fig. 8(d), another sliding homoclinic orbit to comes into existence. In addition, the orbits initiating from in the interior space of the pseudo–homoclinic orbit will tend to once they traverse the switching line and enter into . No cycle happens for , and the real equilibrium is globally asymptotically stable except . It can be seen in Fig. 8(f) that the orbit starting from will not touch the switching line, which is different from the trend in Fig. 8(e). Furthermore, the orbit possessing the same initial value (orange one) will tend to directly once it traverses into for . Fig. 8(e)–(g) illustrate that increasing properly will expedite the convergence of infected population size to the stable state.

It follows from Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8 that the variation of could induce rich dynamics for system (8), different values of can lead the infected population to a stable state, periodicity or even vanish, which suggests that suitable control measures for a certain epidemic is very significant. See term (5) that decreases with the increase of , then bigger value of can prolong the time to implement containment strategies for the epidemic. Note that describes the media impact. Hence, we can implement proper media coverage to curb the disease transmission effectively.

Impacts of key parameters on the dynamics of Filippov system (8)

Note that the model (8) involves media coverage, vaccination for susceptible and recovery treatment for infected. Our purpose is to find an effective strategy to regulate the infected population less than an expectant value or to be eliminated. It is crucial to implement the control strategy by choosing appropriate threshold value . We choose and to be the key parameters and investigate how these parameters affect the dynamics of (8). In this section, we will present detailed analysis and numerical calculation into two subsections: time series analysis for the infected cases and data fitting for the epidemics A/H1N1, COVID-19.

Time series analysis for the infected cases

Based on the purpose of model (8), without implementing control strategies, the disease transmission follows the natural status, i.e., the dynamics of uncontrolled subsystem , and the red dashed curve in Fig. 9 (a) denotes the variation of I with time without control measures. Furthermore, the controlled subsystem is employed for , which implies the infected will vary following the dynamics of in the initial stage of epidemic spread. That is because the quantity of susceptible is relatively large at the beginning during the epidemic outbreak.

Fig. 9

Time series of the infected cases for Filippov system (8), the parameters are: (a). In (b)-(d), all the parameters are the same as presented in (a) except and , which are depicted in each subfigure orderly.

The solid curve consisting of two parts in Fig. 9(a) describes the controlled results, where the blue part follows the dynamics of and the red part obeys the dynamics of . Comparing the two curves in Fig. 9(a), note that once the control measures have been implemented, the peak value of infected cases will reduce to a lower level and tend to the disease-free state faster, which indicates that the containment strategies in our model is very effective for protecting individuals from the infectious diseases. Note that and describe the strengths of three control measures: media coverage, vaccination for susceptible and recovery treatment for infected respectively, which demands us to alter the value of such three parameters to investigate the variation of infected population. Enlarging from 0.15 to 0.3 can reduce the peak value of I as shown in Fig. 9(b), which implies adding more media impact at the beginning of disease transmission can weaken the spread outbreak. Similar tendencies appear in Fig. 9(c)–(d) for increasing and , which implies that strengthening the vaccination and recovery treatment can be helpful as well.

It is necessary and crucial to investigate the joint impact of two parameters on the variation of infected population. Here we select and to examine how the infected size changes by implementing the media coverage and vaccination strategies. The other parameters can be analysed in the same way. In order to control the disease transmission better, it is aimed to reduce the maximum and final size of infected. Especially, the convergence of infected population size to zero (i.e., the epidemic diseases be eliminated) is our main focus.

Note that the maximum and final infected size will both reduce if one enlarges the parameters and ; see Fig. 10 (a)–(b) (the blue region denotes lower level size while the yellow region represents higher level size). Moreover, it is observed from Fig. 10(b) that the final infected population is more sensitive to , where the final size tends to zero rapidly by increasing when keeps fixed. Another objective is to study the elapsed time for the final infected population reducing to zero. Hence, the interval in Fig. 10(b) is selected and the details are expressed in Fig. 10(c). It reads that enlarging the two parameters and can expedite the infected numbers to zero. Note that the basic reproduction numbers of two subsystems are independent of . Thus, two parameters and are chosen to investigate the variation of (); see Fig. 10(d). Note that increasing and can reduce to be less than 1, and hence eliminate the disease.

Fig. 10

The monotonicity of (a) the maximum infected population, (b) the final infected population and (c) the elapsed time to the final infected population with respect to and for the Filippov system (8). Here the parameters are: . (d) The monotonicity of basic reproduction number () with respect to and , where the blue (yellow) surface denotes the variation of ().

Examination of data fitting for epidemics A/H1N1 and COVID-19

The proposed system (8) can provide us an effective switching strategy induced by the media coverage to contain the epidemic spread. Here we validate our model by the data of the laboratory-confirmed cases during the A/H1N1 influenza pandemic in Shaanxi province of China. The dataset is obtained from the Province’s Public Health Information System (Tang et al., 2010, Xiao et al., 2015). The variation of the daily hospital notifications (the magenta circle points) of A/H1N1 flu are depicted in Fig. 11 (a), and the least squares (LS) criterion is applied to explore the well data fitting.

Fig. 11

Fitting Filippov system (8) to the data of A/H1N1 (a) and COVID-19 (b), here the magenta circle points denote the reported daily number of hospital notifications of A/H1N1 flu in Shaanxi Province from September 3 to October 12, 2009 in (a); the global reported daily new laboratory-confirmed COVID-19 cases from April 8 to May 17, 2020 in (b). The curve consisting of blue and red parts presents the data fitting results, and the dashed line in (a) stands the switching time.

For simplicity, we assume that the total population at the outbreak period is a constant and the numbers of recruitment and natural death are ignored; i.e., . Moreover, there is no reported cases for the disease-related death before mid October 2009, and no vaccine was injected until the end of November 2009 (Xiao et al., 2015). Thus, .

A well fitting curve is shown in Fig. 11(a), and all estimated parameters are listed in Table 3 . At the beginning of the epidemic spread, the variation of infected population follows the dynamics of (blue curve), and the subsystem is applied when it reaches the switching time (dashed line). It is observed that the number of infected cases continues to reduce (red curve). Our estimation results () reveal that the recovery treatment plays an important role in containing A/H1N1 epidemic spread. In addition, for the media impact proposed in our model, note that the individual behaviour variations are more sensitive to the changing rate of infected cases than the infected population size ().

Table 3

Estimation of initial values and parameters for Filippov system (8).

Parameters	Description	Estimated value for A/H1N1	Estimated value for COVID-19
S_0	Initial value of susceptible population	3.606×10^4	3.458×10^5
I_0	Initial value of infected population	1	8.735×10^3
Λ	Recruitment rate of S	0	9.27×10^4
β	Probability of transmission from I to S per contact	1.218×10^-4	2.438×10^-5
μ	Natural death rate	0	0.033
δ	Disease-related rate	0	0.252
γ_1	Maximum vaccination rate	0	0
ω_1	Time delayed for vaccinating	-	-
γ_2	Maximum recovery treatment rate	3.903	8.092
ω_2	Limitation of Medical resources	5.102×10^-6	3.867×10^-6
a_1	Weight of media effect to the number of I	6.633×10^-6	6.877×10^-8
a_2	Weight of media effect to the changing rate of I	0.016	1.813×10^-6

The novel coronavirus disease (COVID-19) induced by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has become a worldwide healthy concern since its outbreak. It has spread so fast that more than ten million people around the world are infected by July 10, 2020 and almost all countries are suffering (World Health Organization, 2020a). WHO declared the COVID-19 as a pandemic on March 11, 2020 (World Health Organization, 2020b). In the meantime, some countries released stringent strategies to contain the increasing confirmed infected cases, such as lockdown in Italy on March 9, 2020 (COVID-19 Health System Response Monitor, 2020c). We exhibit the reported daily new laboratory-confirmed COVID-19 cases by WHO from March 1 to July 10 (magenta circle points) in Fig. 12 (a). During this time period the epidemic was getting worse worldwide. It seems that there was a plateau period from April to May, where the new infected population fluctuates around 80,000 cases. Our main purpose is to explore the potential factors that could generate such fluctuation period by fitting the Filippov model (8). More accurately, the data from April 8 to May 17 (depicted by green curve in Fig. 12(a)) are chosen for further discussion.

Fig. 12

(a) The global reported daily new laboratory-confirmed COVID-19 cases (magenta circle points) from March 1 to July 10, 2020, the blue curve is the fitting results and the green part is the plateau period during April 8 to May 17. (b)–(m) The data from April 8 to May 17 in twelve typical countries.

The daily new confirmed COVID-19 cases from April 8 to May 17 in twelve countries where the epidemic was more serious are depicted in Fig. 12(b)–(m), revealing the variations of evolution trends. Note that the epidemic in Spain, Italy, France and Germany had been contained well, where the confirmed cases declined gradually during this period. The infected population in the US and the UK had been fluctuated all the time, while the situation in Brazil, Peru, Russia, Saudi Arabia and India had been worsening. Different from other countries, the epidemic in Iran recrudesced in early May because of the relaxation of restricting strategies (Devi, 2020).

It can be found that the first confirmed COVID-19 case in all twelve territories arose during the end of January to the early March, although the time interval of the first case happening in these countries is not too far, the trends of epidemic spread are various. It follows from Fig. 12(f)–(i) that the infected cases in Spain, Italy, France and Germany are prevented well, which is induced by the effective measures implementing at the initial stage during the epidemic spread, such as lockdown in some districts or even the whole country, cancelling the public activities and closing schools as well as non-essential stores (Chintalapudi et al., 2020, COVID-19 Health System Response Monitor, 2020c, COVID-19 Health System Response Monitor, 2020d, COVID-19 Health System Response Monitor, 2020a, COVID-19 Health System Response Monitor, 2020b). Thus the peak value of new cases had already reached before April 8th. The epidemic situation in severe countries such as the US, Brazil, Russia is mainly caused by the lack of attention to epidemic by the government, the unimplemented stringent control measures timely, or the inadequacy of medical facilities (van Dorn et al., 2020, Lancet, 2020, Lotta et al., 2020). Hence the stringent protective polices implemented in several countries like Italy and unenforced in other territories like the US gave rise to the plateau period in global.

To examine the applicability of system (8) to the COVID-19 epidemic, we still employ the least squares method to detect the fitting results for the data belonging to the plateau. Note that the vaccine of COVID-19 has been researched all the time and not put into use formally yet (Leidy and Garca, 2020), which indicates that . It follows from the fitting curve in Fig. 11(b) that the controlled system ( with the blue curve) and the uncontrolled system ( with the red curve) play an alternate role in the disease transmission, such curve is well fitted to the real data (magenta circle points). The estimated values in Table 3 show that the number of initial susceptible cases () holds a small proportion in the world’s total population, however, the large recruitment rate for susceptible () will induce the exacerbation of this epidemic. During the therapeutic process for COVID-19, the recovery treatment for confirmed cases is the main strategy to implement, which is in accordance with our fitted parameter (). Compared to the number of confirmed cases, its changing rate possesses a bigger effect on the variation of individuals (). The fitting results reveal that the control measures implemented once in a while can generate the fluctuation of infected cases, which is consistent with the aforementioned discussion. It is because the prevention strategies can not be enforced sternly so that the outbreak of COVID-19 epidemic happened repeatedly. The successful experience like that in Italy alerts us that the stringent strategies implemented timely and sustainedly is very vital to contain the disease transmission.

The well fitness depicted in Fig. 11 indicates that the proposed model (8) can be applicable to describe the existing strategies for each country, then we can analyse the impact of such strategies on the epidemic fluctuation as well as the second outbreak. The successful prevention like that in Italy implies that the control measures should be implemented as soon as possible, which accords exactly with the intention of our model. It can be found that the controlled system is employed at the initial stage of epidemic outbreak, where the susceptible cases are very large. Furthermore, the main aim of the stringent strategies in the early phase is to reduce the contact rate from infected to susceptible cases, which demands the government for enhancing the media impact to warn people to make much account of the disease severity and take self-protection.

Conclusion and discussion

Media coverage, vaccination for susceptible population, and treatment for infected population are three main control strategies to protect human from the infectious diseases. It has been investigated via various mathematical models that implementing these three measures could be effective in controlling epidemic outbreaks (Chen et al., 2018, Wang and Xiao, 2013, Zhang et al., 2018). However, these works focused only on one type of the three control strategies and explored its effect on the transmission of an infectious disease. In this paper, we propose a novel Filippov SIR model (8) which takes into consideration all of these three control strategies. We use an implicit formula related to the number of infected cases and its changing rate to address the switching condition (Xiao et al., 2013). Moreover, this formula can be converted into an explicit equation which is a straight line () and represents the switching manifold for the proposed model. For our proposed epidemic Filippov model, if , the natural growing status for the individuals describing with two saturation functions related to vaccination for susceptible and recovery treatment for infected is presented in the uncontrolled subsystem ; on the other hand, if the susceptible cases exceed the threshold value (), the media coverage, linear functions involving vaccination and recovery treatment control strategies are implemented to decide the controlled subsystem .

Compared to the smooth systems, Filippov systems can exhibit very rich dynamics (Kuznetsov et al., 2003, di Bernardo et al., 2008) due to its non-smooth and discontinuous properties. Our proposed Filippov model involves exponential and non-linear functions, which induce many difficulties for theoretical analysis. We also present some interesting and illustrative numerical simulations as depicted in each figure. In comparison to the models which only focus on one type of control measures such as media coverage, vaccination for susceptible population, or recovery treatment for infected population (Buonomo et al., 2008, Wang and Xiao, 2013, Zhang and Liu, 2008), our model has much richer dynamical behaviour and indicates more biological significance for the joint effects of these three strategies. To investigate the dynamics of the proposed model, we study the properties of two subsystems and the bifurcations involving the sliding domain. We also examine how the three strategies presented in our model affect the epidemic spread, and whether our threshold policy is effective in disease control. Our results suggest that it is necessary and effective to implement control strategies as soon as possible at the early stage of a disease outbreak. It is also demonstrated that, by choosing the switching policy properly, one could reduce the infected size to a very low level or even close to zero.

In this paper, we first incorporate the Filippov system (1) with (2) and the implicit threshold policy (3). We use the Lambert W function to convert (3) into an explicit formula (5), which results in the main Filippov system (8). Note that the two non-linear functions presented in induce complex dynamics. It is observed that the subsystem possesses only one disease-free equilibrium which is stable when and unstable when . Moreover, there are at most three endemic equilibria of , and the existence theory of theses three equilibria is given in Cases (A)-(C) and displayed in Table 1. The parameters and are chosen to examine the real and virtual equilibrium bifurcations for ; see Fig. 2. The dynamics of subsystem has been studied in Xiao et al. (2013), and we state their main conclusions in Lemma 2, Lemma 3. A further analysis for the sliding domain of the proposed model is given in Section 4. It is demonstrated that only one sliding segment appears and holds all the time. Hence, the tangent point is above . Only one pseudo-equilibrium may appear on the sliding segment and it is unstable; see Theorem 4.

Note that the sliding bifurcation including local and global sliding bifurcations exhibits very rich dynamics as shown in Fig. 3, Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8. In addition, the local sliding bifurcation related to the boundary equilibrium bifurcation is explored for four types: boundary focus, boundary node, boundary saddle and boundary saddle-node bifurcations; see Fig. 3, Fig. 4, Fig. 5, Fig. 6 respectively. The global sliding bifurcation involving limit cycle shows some interesting dynamics. The grazing (touching) bifurcation, sliding homoclinic bifurcation to the pseudo-saddle and crossing bifurcation are presented in Fig. 7, Fig. 8.

We also examine the impact of the key parameters and related to the aforementioned three control strategies respectively on the dynamics of the model (8). Time series analysis for the infected cases reveals that enlarging theses three parameters are very effective to control the infectious disease transmission. Furthermore, implementing the containment measures can reduce the peak value of infected population and elapsed time to the final state. The joint impact for and on the dynamics of the proposed system has also been examined, where the maximum and final infected sizes as well as the elapsed time to the final state are chosen as indexes. The monotone decreasing of these three indexes with increasing and can be seen in Fig. 10(a)–(c). Furthermore, the basic reproduction numbers of two subsystems related to the parameters and are examined in Fig. 10(d), which indicates that increasing such two parameters may eliminate the infectious diseases. To validate our model by the real data, we fit the Filippov model (8) to the daily number of hospital notifications of A/H1N1 flu in Shaanxi Province from September 3 to October 12, 2009 and the global reported daily new laboratory-confirmed COVID-19 cases from April 8 to May 17, 2020 respectively. The well fitting results in Fig. 11 indicate that our proposed model can be applicable to describe the epidemic spread trend, and the switching strategy induced by the media coverage reflects that reinforcing the media report as early as possible can suppress the epidemic outbreak effectively. In addition, the treatment of infected individuals plays an important role in controlling an epidemic outbreak.

To address the multiple prevention and control measures for the infectious diseases, we have proposed a novel epidemic Filippov model. A detailed theoretical analysis together with numerical simulations reveal how the switching strategies in the proposed system affect the epidemic spread. This provides an effective treatment policy to protect individuals from infection. The results concluded in our work imply that the final outbreak size can be confined in a desired level if the control strategies are implemented appropriately.

CRediT authorship contribution statement

Jiawei Deng: Conceptualization, Software, Writing - original draft. Sanyi Tang: Methodology, Supervision, Funding acquisition. Hongying Shu: Formal analysis, Writing - review & editing.

Declaration of Competing Interest

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.