key: cord-1015080-744juy8r authors: Qin, Wenjie; Tang, Sanyi; Xiang, Changcheng; Yang, Yali title: Effects of limited medical resource on a Filippov infectious disease model induced by selection pressure date: 2016-06-20 journal: Appl Math Comput DOI: 10.1016/j.amc.2016.02.042 sha: 32ca639226fc62fc3f88258b88d1629e8bd98665 doc_id: 1015080 cord_uid: 744juy8r In reality, the outbreak of emerging infectious diseases including SARS, A/H1N1 and Ebola are accompanied by the common cold and flu. The selective treatment measure for mitigating and controlling the emerging infectious diseases should be implemented due to limited medical resources. However, how to determine the threshold infected cases and when to implement the selective treatment tactics are crucial for disease control. To address this, we derive a non-smooth Filippov system induced by selective treatment measure. The dynamic behaviors of two subsystems have been discussed completely, and the existence conditions for sliding segment, sliding mode dynamics and different types of equilibria such as regular equilibrium, pseudo-equilibrium, boundary equilibrium and tangent point have been provided. Further, numerical sliding bifurcation analyses show that the proposed Filippov system has rich sliding bifurcations. Especially, the most interesting results are those for the fixed parameter set as the bifurcation parameter varies, the sliding bifurcations occur sequentially: crossing → buckling → real/virtual equilibrium → buckling → crossing. The key factors which affect the selective treatment measures and the threshold value of infected cases for emerging infectious disease have been discussed in more detail. In the last few years, frequent outbreaks and quick spread of the emerging infectious disease become a worldwide public healthy problem. It endangers not only people's health but also the stability of the whole society. In 2003, SARS (Severe Acute Respiratory Syndrome) killed 774 people, infected more than 80 0 0 globally and threatened to spread around the world [1, 2] . The A/H1N1 (Hemagglutinin 1, Neuraminidase 1) influenza virus, which caused the 2009 pandemic, continues to circulate in some parts of the world, causing variable levels of disease and outbreaks [3] . By September 14, 2014, a total of 4507 probable and confirmed cases, including 2296 deaths from Ebola virus disease (EVD) had been reported from five countries in West Africa-Guinea, Liberia, Nigeria, Senegal, and Sierra Leone [4] . In the early stages of the emerging infectious disease, it is deficient to recognize the emergence of these infectious systematically and comprehensively, the measures of disinfection and isolation are failed to protect and control infected people, so the patients have little resistance to infection of any kind and emerging infectious spread through the camps like wildfire. The emerging infectious disease including SARS, A/H1N1, Ebola, Dengue fever [5] etc. are often accompanied by other viral diseases such as common flu. Moreover, in the early stages of emerging infectious diseases (taking SARS as an example in the rest of paper) outbreak, there is only a very few number of individuals infected by SARS, and the medical resources are enough at this stage. Meanwhile, the early symptoms are very similar to flu and there is no effective way to identify the infected patients. Thus, patients infected with different viruses can be got treatment in different areas of the hospital at the same time. With the growing numbers of SARS infected cases, those pose a grave threat to public health. Meanwhile, various kinds of control strategies have constraints based on the limited medical resources such as doctors, vaccines, drugs, hospital beds, isolation places, medical devices, and so on, especially in rural areas in many developing countries [6] [7] [8] [9] [10] . The medical resource limitation seriously restricts the prevention and treatment for SARS. At this moment, the department of health or state has to cite the urgency of fighting SARS, adopts green passage policy that speeds for isolation and treatment for SARS, so the doctors have to focus their attentions on the SARS infected cases only. For the patients with common flu, the doctors can only prescribe medicines and advise them to go home for home treatment, which can significantly relieve the pressure of limited medical resources on hospital or doctors. In order to describe the effects of limited medical resource and selection strategy discussed above, the number of the patients infected by SARS in a compartment has been chosen as an index for hospital or doctors to use decisions. That is, if the number of the patients infected by SARS is below the threshold level which can be determined analytically (see main text for more details), there is no limited medical resource and selection pressure; above the threshold, due to the limited resource, and doctors treat SARS only. This type of control strategy is called as threshold control policy [11, 12] , which can be described by Filippov systems [13, 14] . Recently, non-smooth Filippov infectious disease models have been investigated by many researchers [10, [15] [16] [17] [18] . In the present work, a non-smooth Filippov infectious disease model with threshold strategy induced by selective treatment measure is derived. The sliding mode dynamics and the existence of all types equilibria have been discussed. Numerical sliding bifurcation analyses show that the proposed Filippov system has rich sliding bifurcations. The key factors which affect the selective treatment measures and the threshold value of infected cases for emerging infectious disease have been discussed in more detail. Our main results show that reducing the threshold value to an appropriate level could contribute to the efficacy on prevention and treatment of emerging infectious disease, which indicates that the selection pressures can be beneficial to prevent the emerging infectious disease under medical resource limitation. Let S ( t ), I 1 ( t ), I 2 ( t ) and R ( t ) denote the numbers of susceptible, the patients with SARS, common flu and recovered individuals at time t , respectively. For simplification, we assume that the people can only be infected either by SARS virus or by common flu virus. Further, based on the classical infectious disease model with limited capacity for treatment [9, 10, [19] [20] [21] [22] [23] [24] [25] we propose the following SI 1 I 2 R model as the basic model in this study The parameters of model (2.1) are summarized in Table 1 . Obviously, two probabilities p 1 , p 2 ∈ [0, 1]. Note that some very special cases of model (2.1) have been studied in our previous work [10] , and the main purpose in this work is to investigate the generalized cases and reveal the rich dynamics and important biological implications concerning emerging infectious disease control. It follows from model (2.1) that the total recovery rate is the major concern for the doctors once the emerging infectious disease outbreaks. Intuitively, how to choose the treatment proportions p 1 and p 2 for the patients infected by different virus such that the total recovery rate reaches its maximal value? To address this question, we discuss the selective strategies in the following. Conditional upon resource limitation, we first assume that medical treatment service for SARS is more than the patients infected by flu, that is Natural recovery rate for flu c 1 The maximal recovery rate per unit time for the patients with SARS c 2 The maximal recovery rate per unit time for the patients with flu b 1 Delayed effects on the treatment for the patients with SARS b 2 Delayed effects on the treatment for the patients with flu p 1 Probability that doctors treat the patients with SARS p 2 Probability that doctors treat the patients with flu Based on the total recovery rate defined by (2.2) , we define the function F with respect to p 1 and p 2 as follows Taking the derivatives of the function F with respect to p 1 and p 2 respectively, one yields which decides whether the hospital carries out the selective strategy or not. Thus, there are two cases: If I 1 < I c , the function F can obtain the maximum value at p 1 = 1 , p 2 = 1 . If I 1 > I c , the function F can obtain the maximum value at Therefore, taking into account above facts, in the early stages of SARS outbreak ( I 1 < I c ), the patients with flu can be treated simultaneously with SARS, i.e., p 1 = p 2 = 1 . Then model (2.1) becomes With the increasing number of infected cases with SARS ( I 1 > I c ), the doctors have to isolate and treat SARS patients only, i.e., p 1 = 1 , p 2 = 0 . Thus model (2.1) becomes To simplify models (2.6) and (2.7) , we assume that the number of the patients infected by flu each year is a constant, i.e., I 2 = k ∈ Z + . Thus, models (2.6) and (2.7) can be rewritten as the following Filippov system [13, 14] ⎧ ⎨ with (2.9) is a description of the threshold control policy, which is referred to as an on-off control, see [11, 12] for more detailed discussion on Filippov system. Note that the special case (i.e. b 1 = 0) has been investigated recently [10] . Let Then model (2.8) with (2.9) can be rewritten as the following generalized Filippov system [13, 14] ˙ Furthermore, the discontinuity boundary (or manifold) separating two regions G 1 and G 2 is described as The main characteristics of Filippov system (2.10) is that selective strategy is suppressed when the number of patients infected by SARS (i.e., I ( t )) is below a previously chosen threshold policy I c . With this number I ( t ) increases and exceeds the threshold I c , the hospital will treat severe cases only since the shortage of medical resources, that is, the selective strategy is implemented. The following definitions on all types of equilibria of Filippov system (2.10) [26, 27] are necessary throughout the paper. Both the real and virtual equilibria are called regular equilibria. A point Z * is called a pseudo-equilibrium if it is an equilibrium of the sliding mode of system (2.10) , i.e., where ·, · denotes the standard scalar product. If I < I c , then the following system plays a key role in analyzing the Filippov system (2.10) ˙ and the basic reproduction number of subsystem (3.1) reads It is obvious that subsystem (3.1) always has a unique disease-free equilibrium E 0 The endemic equilibria of subsystem (3.1) are solutions of For the simplicity and convenience of exposition, we denote Next, we will study the stability of the endemic equilibria of subsystem (3.1) . The characteristic equation about the endemic equilibrium E i We can obtain the following lemmas from (3.3) . 1 is an unstable node or focus when H(I 1 1 ) < 0 , and system (3.1) has at least one closed orbit in ; E 1 1 is a center of the linear system when H(I 1 is a saddle; and E 1 1 is an unstable node or focus when H(I 1 The methods of proving Lemmas 3.1 -3.3 are similar to those in Refs. [23, 24] , see these references for more details. If I > I c , then the Filippov system (2.10) becomes ˙ (3.5) and it has a unique disease-free equilibrium E 0 2 = (A/μ, 0) , which is globally asymptotically stable if R 2 < 1 . Here, For convenience, we denote The characteristic equation about E i Similar conclusions as Lemmas 3.1 -3.3 for subsystem (3.5) can be obtained, and those are not described here. It follows from Filippov convex method [13, 14] that one can define the sliding vector field as a convex combination of the two vector fields Noting that the control α(Z) = 1 indicates that the flow is governed by F G 1 alone, which must be tangent to the switching surface . Analogously, α(Z) = 0 shows a tangency of flow F G 2 with . Therefore, the sliding subset can be defined as Therefore, the sliding segment of Filippov (2.10) can be defined as Here we employ Utkin's equivalent control method introduced in [14 ] to obtain the differential equation for sliding dynamics defined in the region s . It follows from H = 0 that and solving the above equation with respect to ε yields . According to Utkin's equivalent control method, we can obtain the dynamics defined in s which can be determined by the following scalar differential equation (4.4) For the simplicity and convenience of exposition, we denote real equilibrium as E R , virtual equilibrium as E V , pseudoequilibrium as E P , boundary equilibrium as E B and tangent point as E T , respectively. Regular equilibrium: For convenience, we just consider the subsystem has two endemic equilibria. From Lemma 3.1 , subsystem (3.1) has two endemic equilibria E 1 Pseudo-equilibrium: For the existence of pseudo-equilibrium, we denote E P = (S P , I c ) , according to (4.4) the S P component of the pseudo-equilibrium of sliding flow satisfies the following equation where S P = A/ (μ + βI c ) ∈ s , and S P is a unique positive steady state of (4.4) . That is, if S P = A/ (μ + βI c ) ∈ s holds true, Filippov system (2.10) exists a unique pseudo-equilibrium E P . For the stability of pseudo-equilibrium E P = (S P , I c ) , we rewrite sliding mode equation (4.4) as ˙ it follows from (4.5) that which indicates that the pseudo-equilibrium E P is locally stable in s . Boundary equilibrium: The boundary equilibria of Filippov system (2.10) satisfy equations Hence, there may be two possible tangent points including E 1 The richness of equilibria of Filippov system (2.10) could result in a number of equilibrium bifurcations as the key parameter varies. Thus, in the coming section, we would like to investigate the local and global sliding bifurcations of system (2.10) . By employing numerical methods to investigate the qualitative behaviors of Filippov system (2.10) , the sliding bifurcation analyses including local and global bifurcations of one parameter Filippov system (2.10) are provided in this section. Therefore, we fix all other parameters and choose the maximal recovery rate for the patients with flu (i.e., c 2 ) as a bifurcation parameter, noting that c 2 is directly related to the threshold value I c . Throughout this section, we will investigate the local sliding bifurcation of Filippov system (2.10) . Boundary equilibrium bifurcation, as a type of local sliding bifurcation in Filippov system, is characterized by the collision of pseudo-equilibrium, tangent point, and real equilibrium (or tangent point and real equilibrium) at the discontinuity surface when one parameter passes through a critical value, and includes boundary node, focus and saddle bifurcations which will be addressed in the following. Here the symbol DF G i (E B ) is the characteristic polynomial of subsystems (3.1) or (3.5) about the boundary equilibrium E B ). These bifurcations are classified as boundary saddle, boundary node and boundary focus in [31] . It follows from Section 4.2 that is an saddle (a note, or a focus). Again, similar argument as above will yield for the boundary equilibria E 21 B and E 22 B . Hence, a boundary equilibrium bifurcation occurs at E i j B , i, j = 1 , 2 . Boundary-saddle bifurcation: This type of bifurcation may occur for Filippov system (2.10) if three types of equilibria E R , E T and E P collide together simultaneously as parameter c 2 passes through a critical value [31] . For example, when the parameter c 2 passes through a critical value c t 2 ≈ 3 , a saddle E 22 R , a tangent point E 2 T and a pseudo-equilibrium E P collide together, so the boundary-saddle bifurcation occurs at E 2 B , as shown in Fig. 1 (B) . In this case, the critical value R (an unstable real equilibrium), a stable pseudo-equilibrium E P and an invisible tangent E 2 T can coexist for c 2 < c t 2 , as shown in Fig. 1 (A) with c 2 = 2 . 8 . They collide together simultaneously as c 2 = c t 2 and are substituted by a virtual R and a tangent point E 2 T coexist, as shown in Fig. 3 (A) when c 2 < 0.635. They collide at c 2 = c t 2 (see Fig. 3 (B) with c 2 = 0 . 5 ) and are substituted by a pseudo-equilibrium E P , a tangent point E 2 T and a virtual equilibrium E 21 V when c 2 > 0.635, as shown in Fig. 3 (C) with c 2 = 0 . 7 . Filippov system (2.10) could have standard periodic solutions that lie entirely in regions G 1 or G 2 through a Hopf bifurcation, respectively [31, 32] . Meanwhile, as mentioned in Refs. [26, 31] , Filippov system (2.10) may have additional two types of new periodic solutions: periodic solutions which have a sliding segment in sliding segment s (i.e., sliding periodic solutions) and those which have only isolated points in common with s (i.e., crossing periodic solutions). Noting that a crossing periodic solution can pass through the boundary of the sliding segment s . Accordingly, the orbits corresponding to periodic solutions will be called standard, sliding and crossing cycles. In this section, we focus on the global sliding bifurcations such as grazing bifurcation (i.e., touching bifurcation), buckling bifurcation and crossing bifurcation. Grazing (or touching) bifurcation: It follows from the Refs. [26, 31] , a standard periodic solution can collide with the sliding segments, and this type of bifurcation is called grazing or touching bifurcation. Noting that Filippov system (2.10) has a stable periodic solution in the interior of region G 2 , as shown in Fig. 4 (A) with c 2 = 0 . 5 . At this moment, Filippov system (2.10) has two tangent points E 1 T and E 2 T lying on the boundary of the sliding mode, subsystem (3.5) has an unstable real equilibrium E 2 R while subsystem (3.1) has one unstable virtual equilibrium E 1 V . As the parameter c 2 increases and passes through around 0.85, a grazing or touching bifurcation occurs, as shown in Fig. 4 (B) , which indicates that the standard period solution of Filippov system (2.10) collides with its tangent point E 2 T . As c 2 continues to increase the cycle becomes a sliding cycle, where a piece of sliding segment belongs to the cycle, as shown in Fig. 4 (C) with c 2 = 0 . 95 . Especially, as the bifurcation parameter c 2 is increased to 1.6, the stable periodic cycle is disappeared, and pseudoequilibrium E P appears at c 2 = 1 . 6 . Meanwhile, for subsystem (3.5) , the real equilibrium E 2 R becomes a virtual equilibrium E 2 V , the real/virtual equilibrium bifurcation occurs at c 2 = 1 . 6 , as shown in Fig. 4 (D) . Meanwhile, Fig. 4 (D)also shows that pseudo-equilibrium of Filippov system cannot coexist with the real equilibria. Buckling bifurcation: This type of bifurcation is defined as a standard piece of the cycle starts to pass the invisible quadratic tangent point as the bifurcation parameter varies [31] , as shown in Fig. 5 (C). Fig. 5 (B) clearly shows that Filippov system (2.10) has a stable sliding periodic solution, an invisible quadratic tangent point E 1 T , an unstable regular equilibria E 1 V and E 2 R . As c 2 increases and exceeds 0.7 the piece of the cycle starts to pass the invisible quadratic tangent point E 1 T , and consequently the cycle passes through the whole piece of the sliding segment s , and the sliding cycle does exist in region G 2 and sliding segment s , see Crossing bifurcation: Along with the variation of bifurcation parameter, a stable sliding periodic solution becomes a stable crossing periodic solution, this type bifurcation is called crossing bifurcation. It is interesting to note that a sliding cycle with a single sliding segment ending at E 1 T does exists in both regions G 1 , G 2 and sliding segment s , see Fig. 5 (H). As c 2 increases and reaches 3.45 the sliding cycle only passes the tangent point E 1 T on the sliding segment s . At this moment, the stable sliding periodic solution becomes a stable crossing periodic solution, as shown in Fig. 5 (I) with c 2 = 3 . 45 . Fig. 5 (J) shows that the sliding crossing cycle becomes a crossing cycle as c 2 continues to increase and reach 4. Noting that the sliding segment s is located within the crossing cycle which exists in both regions G 1 and G 2 . It follows from Fig. 5 (H)-(J) that a crossing bifurcation occurs at c 2 = 3 . 45 . Similarly, Filippov system (2.10) must exist one crossing bifurcation from Fig. 5 (A) to (B) . In summary, as bifurcation parameter c 2 increases from 0.5 to 4 and all other parameters are fixed as those in Fig. 5 , Filippov system (2.10) has rich sliding bifurcations, the following local and global sliding bifurcations occur sequentially: crossing → buckling → real/virtual equilibrium → buckling → crossing. Especially, as the bifurcation parameter c 2 changes around 1.6, the stable periodic cycle is disappeared, and pseudo-equilibrium E P appears at c 2 = 1 . 6 , all the orbits tend to the pseudo-equilibrium E P , which is locally asymptotically stable. There exists a real/virtual equilibrium bifurcation for Filippov system (2.10) in such case. Previous analysis indicates that limited medical resources, the basic reproduction numbers R i , i = 1 , 2 and threshold values I c (or the selective strategy) are significant factors affecting the spread of the emerging infectious disease. In this section, we first investigate how the limited medical resources (i.e., b 1 ) affects the dynamic behaviors of Filippov system (2.10) or the spread of SARS. As the medical resource limitation (i.e., b 1 = 0) is taken into account, then the dynamical behaviors of both subsystems (3.1) and (3.5) of Filippov system (2.10) become much more complex. Meanwhile, the dynamic behavior of the Filippov system (2.10) could be dramatically affected by the existence of medical resource limitation. In fact, the sliding mode could change as threshold value I c changes, as shown in Fig. 6 (A) , the length of sliding segment is increased with growing recovery rate c 1 from 0.6 to 1.5. Meanwhile, noting that Filippov system (2.10) does not have pseudo-equilibrium in this case b 1 = 0 . However, with the medical resource limitation (i.e., b 1 = 0), there exists a pseudo-equilibrium E P , which is globally stable with respect to the sliding segment s , as shown in Fig. 6 (B) . Those indicates that the emerging infectious disease will become endemic instead of elimination. Meanwhile, as seen in Fig. 7 , the rate of growth in the number of the patients with SARS is far more quickly before than after controlling. By comparing the red line with the blue one in Fig. 7 , the limited medical resource has been a great influence on the time when the peak of patients appears and how long does the spread of SARS last. In fact, the blue line (with sufficient medical resources, i.e., b 1 = 0 ) shows that there is a peak of patients in rapid and a sharp decline, the spread of SARS lasts only for a short time. The red line (with limited resources, i.e., b 1 = 2 . 5 ) shows there is a delay in appearance of the peak of patients, and the peak value increases obviously, which indicates that the limited medical resources does not facilitate the treatment of infectious disease. Therefore, in order to prevent and control the spread of emerging infectious disease, it is crucial to implement the selective strategy timely to control the basic reproduction numbers R i < 1 , i = 1 , 2 . Then, the key parameters and threshold value I c which affect the basic reproduction numbers R i and selective strategy are investigated, respectively. Obviously, it follows from the expressions of the basic reproduction numbers R i that they are monotonically decreasing functions with respect to c 1 . Meanwhile, we can solve R 01 = 1 with respect to c 1 and derive that c * 1 = (Aβ − μν)(1 + b 2 k ) μ , (6.1) so R 01 (c 1 ) < 1 provided c 1 > c * 1 , as shown in Fig. 8 (A) . An improved the maximum cure rate for SARS will help prevent and control the spread of SARS. Furthermore, the threshold value I c is a monotonically decreasing function with respect to c 1 . Therefore, in order to control the spread of SARS, it is crucial to reduce the threshold value I c to I * c , here That is, the smaller the threshold value I c is, the more beneficial to prevent and control SARS, as shown in Fig. 8 (B) . Those indicate that it is best to timely selective treatment for SARS infected cases so as not to miss the best timing of treatment. In order to understand the effect of the delayed treating ( b 1 = 0) for the patients with emerging infectious disease on control strategy, we have deduced a non-smooth infectious disease model induced by selection pressures. Analysis of this model reveals rich dynamics including local and global sliding bifurcations. Our main results show that reducing the threshold value to an appropriate level could contribute to the efficacy on prevention and treatment of emerging infectious disease, which indicates that the selection pressures can be beneficial to prevent the emerging infectious disease under medical resource limitation. Comparing the results for the model with b 1 = 0 we have studied in Ref. [10] with b 1 = 0, we conclude that the term b 1 makes the dynamical behavior of the system change more interesting and complicated. By using theoretical techniques [26, 27, 31] , the existence conditions for sliding segment, sliding mode dynamics and different types of equilibria such as regular equilibrium, pseudo-equilibrium, boundary equilibrium and tangent point have been provided. Further, numerical sliding bifurcation analyses show that the proposed Filippov system has more rich local and global sliding bifurcations than the case studied in Ref. [10] , as shown in Figs. 1 -5 for more details. Especially, the most interesting results are those for the fixed parameter set as the bifurcation parameter varies, the sliding bifurcations occur sequentially: crossing → buckling → real/virtual equilibrium → buckling → crossing. According to the analyses of key parameters and biological significance, the results indicate that the dynamic behavior of the Filippov system (2.10) could be dramatically affected by the existence of medical resource limitation (i.e., b 1 = 0), see Fig. 6 (B) for more detail, the existence and stability of pseudo-equilibrium shows that the emerging infectious disease will become endemic. That is, the case shown in Figs. 4 (D) and 5 (E), which reveals that there are several hidden factors that can adverse affect the control strategy under medical resource limitation. Therefore, it is very necessary to implement the selective strategy in the control and treatment of SARS, see Fig. 8 . Meanwhile, this points to the urgent need for improvement in medical facilities, access to rapid diagnosis and treatment with more effective drugs for SARS. SARS: A Case Study in Emerging infections Transmission of infectious diseases during commercial air travel Extrapolating from sequence the 2009 h1n1 'swine' influenza virus WHO ebola response team, ebola virus disease in west Africa-the first 9 months of the epidemic and forward projections Dengue fever: Mathematical modelling and computer simulation Optimal vaccine allocation for the early mitigation of pandemic influenza Dissolving polymer microneedle patches for influenza vaccination Vaccine-delivery patch with dissolving microneedles eliminates sharps boosts protection Nonlinear pulse vaccination in an SIR epidemic model with resource limitation The selection pressures induced non-smooth infectious disease model and bifurcation analysis Sliding Modes and Their Applications in Variable Structure Systems Sliding Modes in Control and Optimization Differential Equations with Discontinuous Righthand Sides Sliding Model Control in Electromechanical Systems Sliding mode control of outbreaks of emerging infectious diseases Sliding bifurcation and global dynamics of a filippov epidemic model with vaccination A Filippov system describing media effects on the spread of infectious diseases Global dynamics of a piece-wise epidemic model with switching vaccination strategy Bifurcation in an epidemic model with constant removal rate of the infectives Backward bifurcation of an epidemic model with treatment Backward bifurcation of an epidemic model with standard incidence rate and treatment rate Backward bifurcation of an epidemic model with saturated treatment function Saturation recovery leads to multiple endemic equilibria and backward bifurcation Rich dynamics of an epidemic model with saturation recovery Dynamics of an SIR epidemic model with limited medical resources revisited Bifurcations in nonsmooth dynamical systems Generic bifurcations of low codimension of planar Filippov systems Sur les équations de Hamilton A singular approach to discontinuous vector fields on the plane Canard cycles and Poincar index of non-smooth vector fields on the plane One parameter bifurcations in planar Filippov systems Multiparametric bifurcation analysis of a basic two-stage population model This work is supported by the National Natural Science Foundation of China (NSFCs: 11471201 , 11171199 , 11301320 , 11371030 ), the Fundamental Research Funds for the Central Universities (GK201305010, GK201401004), the Youth Foundation of China Three Gorges University (KJ2015A006), and the Natural Science Foundation of Hubei province ( 2015CFB264 ).