Filippov Ratio-Dependent Prey-Predator Model with Threshold Policy Control

Abstract

The Filippov ratio-dependent prey-predator model with economic threshold is proposed and studied. In particular, the sliding mode domain, sliding mode dynamics, and the existence of four types of equilibria and tangent points are investigated firstly. Further, the stability of pseudoequilibrium is addressed by using theoretical and numerical methods, and also the local sliding bifurcations including regular/virtual equilibrium bifurcations and boundary node bifurcations are studied. Finally, some global sliding bifurcations are addressed numerically. The globally stable touching cycle indicates that the density of pest population can be successfully maintained below the economic threshold level by designing suitable threshold policy strategies.

1. Introduction

Ordinary differential equation models (ODE models) are widely used to describe the dynamics between predators and their prey, which has long been and will continue to be one of significant fields in mathematical ecology owing to its universal existence and importance [1, 2]. The simplest prey-predator dynamic model is the Lotka-Volterra model [3], which has been modified in many ways since its original and realism formulation in the 1920s.

One important component of the prey-predator relation is predator's functional response which refers to the change in the density of prey attached per unit time per predator as the prey density changes and makes the prey-predator system more realistic. There are several famous functional response types in previous work, which are monotonically increasing and uniformly bounded functions in the first quadrant. Another functional response is the Michaelis-Menten (or Holling-type II) functional response, which is the most common type of functional response among arthropod predators [4, 5]. It takes the form 𝑝(𝑥,𝑦)=𝑎𝑥/(𝑏+𝑥), where 𝑎 and 𝑏 are positive constants that stand for capturing rate, and half capturing saturation constant, respectively.

Considering predators having to search for food, a more suitable general prey-predator theory based on the so-called ratio-dependent theory is involved. It can be roughly comprehended as the per capita predator growth rate should be a function of the ratio of prey to predator abundance. And it is also strongly supported by numerous fields, laboratory experiments, and observations [6–8].

Therefore, we can write the ratio-dependent prey-predator model with Michaelis-Menten functional response as follows: 𝑥̇(𝑡)𝑥=𝑟𝑥1−𝑘−𝑎𝑥𝑦,𝑦̇(𝑡)𝑏𝑦+𝑥=−𝛿𝑦+𝛽𝑎𝑥𝑦,𝑏𝑦+𝑥 (1) where 𝑥 and 𝑦 represent the density of prey (pest) and predator (natural enemy), respectively. The prey is assumed to grow logistically and 𝑘 is the carrying capacity of prey. The positive constants 𝑟 and 𝛿 stand for intrinsic growth rate of prey and mortality rate of predator, respectively. 𝛽 denotes the conversion rate of prey captured by predator.

In population, both ecologist and mathematicians are interested in the ratio-dependent prey-predator model with Michaelis-Menten functional response [2, 7–9]. Hsu et al. [7] resolved a complete classification of the asymptotic behavior of the solutions of ratio-dependent model with the Michaelis-Menten functional response. They also studied the global stability of all equilibria in various cases and reconsidered the uniqueness of limit cycle.

It is well known that pests have been one of the principal threats to crops, important plants, animals, and humans all over the world. Therefore, it is necessary to apply acceptable and effective strategies to control pest outbreak. In practice, it is impossible to eradicate the pests completely, nor is it biologically or economically desirable. Integrated pest management (IPM) is a long term management strategy [10–12], which uses a combination of biological, cultural, and chemical tactics so as to lower cost to the growers, minimize effect on the environment, and maintain pest population below the economic injury level (EIL). On the basis of IPM, biological strategy is useful and effective to suppress pest population, such as releasing beneficial natural enemies; culture strategy makes the environment less favorable to pests, such as catching or harvesting artificially. In most cropping systems, when the above two tactics are unable to keep pest population below the ET, chemical strategy (i.e., insecticide) is still a principal means to control pests and prevent economic loss. Thus, in order to control pest outbreak, we should carry out control strategies when the number of pests reaches or exceeds the ET which is lower than the EIL, and the control strategies should be suspended once the density of pest population falls below the ET, which is the so-called threshold policy control (TPC). Considering IPM strategies, either fixed moment or state-dependent impulsive models with the ratio-dependent or Michaelis-Menten-type response function of prey-predator model have been studied in [13–16].

However, Zhao et al. [17] have stated some disadvantages of the impulsive differential equation models mentioned above. First, in the fixed moment impulsive model, without consideration whether the density of pest reaches the ET or not, control strategies are invariably implemented, which leads to consumption of vast resources. Second, in reality, all kinds of control strategies need some time and cannot be finished instantaneously, but in the state-dependent impulsive differential models, control strategies are carried out instantaneously, which is not reasonable.

Therefore, we use Filippov system which is a vector differential equation with discontinuous right-hand side to describe prey-predator model with both noninstantaneous interventions and the threshold policy. Recently, although Filippov systems have been widely utilized in science and engineering, including harvesting thresholds, oil well drilling, and liquid-gas reaction [18–23]. However, very little is involved that they are used to investigate the ratio-dependent-type predator-prey model with Michaelis-Menten-type functional response. We assume that a proportion of preys are caught or transferred (culture strategy) or killed (chemical strategy), denoted by 𝑞₁; a proportion of predators are released (biological strategy), denoted by 𝑞₂. So we have the following control model for 𝑥>ET: 𝑥̇(𝑡)𝑥=𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝑞₁𝑦̇(𝑡)𝑥,=−𝛿𝑦+𝛽𝑎𝑥𝑦𝑏𝑦+𝑥+𝑞₂𝑦. (2)

In this paper, we aim to give a detailed analysis of Filippov ratio-dependent prey-predator model with threshold policy control, which describes that control measures are implemented only when the density of pest in a population exceeds the ET. We investigate the sliding mode domain, sliding mode dynamics, the existence of four types of equilibria and tangent point of Filippov system, regular/virtual equilibrium bifurcation, and boundary node bifurcations. In addition, the local stability of pseudoequilibrium implies global stability in our numerical simulations. Globally touching bifurcation especially indicates that the density of pest can be successfully maintained below the ET by designing suitable threshold policy strategy. Therefore our control objective can be achieved in the above two cases, which are desired situations in crop, livestock sectors and forestry.

The organization of this paper is as follows: in Section 2, we give some basic results and preliminaries for ODE system and Filippov system. In Section 3, the existence of sliding segments and sliding mode dynamics for the Filippov system (3) are addressed. In Section 4, we give the null-isoclines and equilibria. Based on those results, in Section 5, we consider the bifurcation sets of equilibria and sliding bifurcation analyses. Then, in the last section, we give some discussions.

2. The ODE System and Filippov System

2.1. The Basic Preliminaries and Results for ODE System

The ODE model (1) has been well studied in [2, 7], and a complete classification of the asymptotic behavior of the solutions of the ratio-dependent model with Michaelis-Menten functional response has been proposed. In the following, we present some primary results in the following Lemma which are useful for this study.

Lemma 1. System (1) includes three equilibria (0,0), (𝑘,0), and a unique positive equilibrium 𝐸^∗(𝑥^∗,𝑦^∗) if and only if the following two conditions are true: 𝛽𝑎−𝛿>0 and 𝑟>(𝛽𝑎−𝛿)/𝛽𝑏, where 𝐸^∗=(𝑥^∗,𝑦^∗)=(𝑘(𝛽𝑏𝑟+𝛿−𝛽𝑎)/𝛽𝑏𝑟,𝑘(𝛽𝑎−𝛿)(𝛽𝑏𝑟+𝛿−𝛽𝑎)/𝛽𝑏²𝑟𝛿). If 𝛽𝑎−𝛿>0 and 𝛽𝑎−𝛿<𝑟≤(𝑎−𝑏𝛿)/𝑏, (0,0) is globally asymptotically stable; if 𝛽𝑎−𝛿<0 and 𝑟≥(𝑎−𝑏𝛿)/𝑏, (𝑘,0) is globally asymptotically stable; if 𝛽𝑎−𝛿>0 and (𝛽𝑎−𝛿)/𝛽𝑏<𝑟≤𝛽𝑎−𝛿 or 𝛽𝑎−𝛿>0, 𝑟>𝛽𝑎−𝛿, and 𝑟≥(𝑎−𝑏𝛿)/𝑏, the positive equilibrium of 𝐸^∗=(𝑥^∗,𝑦^∗) is globally asymptotically stable. However if the positive equilibrium is locally asymptotically stable, then the system (1) has no nontrivial positive periodic solutions. If 𝛽𝑎−𝛿>0 and max{𝛽𝑎−𝛿,(𝛽𝑎−𝛿)/𝛽𝑏}<𝑟<(𝑎−𝑏𝛿)/𝑏 hold true, then the system (1) has at most one stable limit cycle.

2.2. Filippov Ratio-Dependent Prey-Predator Model and Preliminaries

By now, based on IPM strategies and TPC, the models (1) and (2) can be incorporated and rewritten in the following form: 𝑥̇(𝑡)𝑥=𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝜀𝑞₁𝑦̇(𝑡)𝑥,=−𝛿𝑦+𝛽𝑎𝑥𝑦𝑏𝑦+𝑥+𝜀𝑞₂𝑦, (3) with 𝜀=0,𝑥<ET,1,𝑥>ET. (4)

We first introduce some useful properties and definitions on Filippov system according to [24, 25], so that we can investigate the model (3) in more detail. Let 𝐻(𝑍)=𝑥−ET with vector 𝑍=(𝑥,𝑦)^𝑇, and 𝐹_𝑆1(𝑍)=𝑥𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥,−𝛿𝑦+𝛽𝑎𝑥𝑦𝑏𝑦+𝑥^𝑇,𝐹_𝑆2(𝑍)=𝑥𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝑞₁𝑥,−𝛿𝑦+𝛽𝑎𝑥𝑦𝑏𝑦+𝑥+𝑞₂𝑦^𝑇. (5) Then the system (3) can be rewritten as the following Filippov system: 𝑍̇(𝑡)=𝐹_𝑆1(𝑍),𝑍∈𝑆₁,𝐹_𝑆2(𝑍),𝑍∈𝑆₂. (6)

In addition, we define the discontinuity boundary set (or the switching line) Σ={𝑍∈𝑅²₊∣𝐻(𝑍)=0}, which divides 𝑅²₊ into two regions 𝑆₁ and 𝑆₂, where 𝑆₁=𝑍∈𝑅²₊(𝑍),𝑆∣𝐻<0₂=𝑍∈𝑅²₊∣𝐻(𝑍)>0. (7)

From now on, Filippov system (3) in different regions 𝑆₁ or 𝑆₂ is named as system 𝑆₁ (i.e., system (1)) or system 𝑆₂ (i.e., system (2)) correspondingly.

Denote 𝜎(𝑍)=𝐻_𝑍(𝑍),𝐹_𝑆1(𝑍)𝐻_𝑍(𝑍),𝐹_𝑆2(𝑍), (8) where 𝐻_𝑍 is a nonvanishing gradient of the smooth scale function 𝐻 on Σ, and ⟨⋅⟩ denotes the standard scalar product; then the sliding mode domain can be defined as Σ_𝑆={𝑍∈Σ∣𝜎(𝑍)≤0}. (9)

We distinguish the following regions on Σ:(i)Σ₁∈Σ is the escaping region if ⟨𝐻_𝑍(𝑍),𝐹_𝑆1(𝑍)⟩<0 and ⟨𝐻_𝑍(𝑍),𝐹_𝑆2(𝑍)⟩>0 on Σ₁;(ii)Σ₂∈Σ is the sliding region if ⟨𝐻_𝑍(𝑍),𝐹_𝑆1(𝑍)⟩>0 and ⟨𝐻_𝑍(𝑍),𝐹_𝑆2(𝑍)⟩<0 on Σ₂;(iii)Σ₃∈Σ is the sewing region if ⟨𝐻_𝑍(𝑍),𝐹_𝑆1(𝑍)⟩⟨𝐻_𝑍(𝑍),𝐹_𝑆2(𝑍)⟩>0 on Σ₃.

The following definitions about all types of equilibria for Filippov system are necessary throughout the paper, so we list them as follows.

Definition 2. A point 𝑍^∗ is called a regular equilibrium of system (3) if 𝐹_𝑆1(𝑍^∗)=0, with 𝐻(𝑍^∗)<0, or 𝐹_𝑆2(𝑍^∗)=0, with 𝐻(𝑍^∗)>0. A point 𝑍^∗ is called a virtual equilibrium of system (3), if 𝐹_𝑆1(𝑍^∗)=0, with 𝐻(𝑍^∗)>0, or 𝐹_𝑆2(𝑍^∗)=0, with 𝐻(𝑍^∗)<0.

Definition 3. A point 𝑍^∗ is called a pseudoequilibrium if it is an equilibrium of the sliding mode of system (3); that is, (1−𝜆)𝐹_𝑆1(𝑍^∗)+𝜆𝐹_𝑆2(𝑍^∗)=0, 𝐻(𝑍^∗)=0, and 0<𝜆<1, where 𝐻𝜆=^𝑍(𝑍),𝐹_𝑆1(𝑍)𝐻_𝑍(𝑍),𝐹_𝑆1(𝑍)−𝐹_𝑆2(𝑍). (10)

Definition 4. A point 𝑍^∗ is called a boundary equilibrium of system (3) if 𝐹_𝑆1(𝑍^∗)=0, with 𝐻(𝑍^∗)=0, or 𝐹_𝑆2(𝑍^∗)=0, with 𝐻(𝑍^∗)=0.

Definition 5. A point 𝑍^∗ is called a tangent point of system (3) if 𝑍^∗∈Σ_𝑆 and ⟨𝐻_𝑍(𝑍^∗),𝐹_𝑆1(𝑍^∗)⟩=0, or ⟨𝐻_𝑍(𝑍^∗),𝐹_𝑆2(𝑍^∗)⟩=0.

The details and knowledge about the Filippov system, such as the concepts of Filippov solution, sliding mode solution, and bifurcation can be found in reference [24].

3. Sliding Region and Sliding Mode Dynamics

A sliding mode exists if there are regions in the discontinuity boundary Σ, where the vectors for the two subsystems of the system (3) are directed towards each other. It is well known that two basic methods the so-called Filippov convex method [24] and Utkin equivalent control method [18], are developed for the sliding mode and its domains, which are shown in the appendix.

3.1. Sliding Segment and Region

Based on the appendix, we have 𝜆(𝑍)=𝑟(1−ET/𝑘)−𝑎𝑦/𝑏𝑦+ET𝑞₁, (11) because the sliding mode regions can be determined by solving the inequalities; that is, 𝜆(𝑍)≥0 and 𝜆(𝑍)≤1. In order to solve the above two inequalities with respect to 𝑦, we need to consider the following two algebraic equations: 𝑟1−ET𝑘−𝑎𝑦𝑟𝑏𝑦+ET=0,1−ET𝑘−𝑎𝑦𝑏𝑦+ET−𝑞₁=0. (12)

Solving the above two algebraic equations with respect to 𝑦 yields two roots, denoted by 𝑦₁=𝑟ET(𝑘−ET)𝑎𝑘−𝑏𝑟(𝑘−ET),𝑦₂=𝑟(𝑘−ET)−𝑞₁𝑘ET𝑎𝑘−𝑏𝑟(𝑘−ET)−𝑞₁𝑘. (13)

Based on the relations between 𝑦₁ and 𝑦₂, there exist two cases for the existence of sliding segments of Filippov system (3). By simply calculating and arranging, we have the following results. (i)When 𝑟<𝑎𝑘/𝑏(𝑘−ET), the sliding segment can be described as Σ¹_𝑆=𝑥,𝑦∣max0,𝑦₂≤𝑦≤𝑦₁,𝑥=ET. (14)(ii)When max{𝑞₁𝑘/(𝑘−ET),𝑎𝑘/𝑏(𝑘−ET)}<𝑟<min{𝑘(𝑎+𝑏𝑞₁)/𝑏(𝑘−ET),𝑎(𝑘+𝑏𝑞₁)/𝑏(𝑘−ET)}, the sliding segment can be described as Σ²_𝑆=𝑥,𝑦∣0≤𝑦≤𝑦₂,𝑥=ET. (15)

3.2. Sliding Mode Dynamics

Filippov system (3) only has one piece of sliding segment, and the solutions defined in it can be obtained from the sliding mode dynamics, which can be determined by employing the Utkin equivalent control method (see the appendix).

From 𝐻=0, we get that 𝐻_𝑍=𝑥^′𝑥=𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝜀𝑞₁𝑥=0,𝑥=ET. (16) And solving the above equations with respect to 𝜀 yields 𝑟()𝜀=1−ET/𝑘−𝑎𝑦/𝑏𝑦+ET𝑞₁. (17)

Hence, the dynamics on the sliding mode Σ_𝑆 can be determined by the following scalar differential equation: 𝑦̇(𝑡)=−𝛿𝑦+𝛽𝑎𝑦ET+𝑟()𝑏𝑦+ET1−ET/𝑘−𝑎𝑦/𝑏𝑦+ET𝑞₁𝑞₂𝑦,=𝑃₁𝑦𝑃₂𝑦+𝑃₃≜𝜙𝑦, (18) where 𝑦∈Σ¹_𝑆 or Σ²_𝑆,  𝑃₁(𝑦)=𝑦/𝑞₁(𝑏𝑦+ET), 𝑃₂=−𝛿𝑏𝑞₁+𝑟𝑏𝑞₂(1−ET/𝑘)−𝑎𝑞₂, and 𝑃₃=−𝛿𝑞₁ET+𝛽𝑎𝑞₁ET+𝑟𝑞₂(1−ET/𝑘)ET.

4. Null-Isoclines and Equilibria

4.1. The Null-Isoclines of Filippov System (3)

Null-isoclines of both systems 𝑆₁ and 𝑆₂ are related to the existence of equilibria and are useful for analysis of sliding dynamic.

Null-isoclines 𝑥̇(𝑡)=0 and 𝑦(̇𝑡)=0 for both systems 𝑆₁ and 𝑆₂ can be determined as follows. For the system 𝑆₁, solving the equation of the null-isocline 𝑥̇(𝑡)=0, yeilds 𝑓_𝑆1(𝑥)=𝑟𝑥(1−𝑥/𝑘)(),𝑎−𝑏𝑟1−𝑥/𝑘 (19) and null-isocline 𝑦(̇𝑡)=0 gives 𝑔_𝑆1(𝑥)=𝛽𝑎−𝛿𝑥.𝑏𝛿 (20)

For the system 𝑆₂, solving the equation of the null-isocline 𝑥̇(𝑡)=0, yields 𝑓_𝑆2(𝑥)=𝑟𝑥(1−𝑥/𝑘)−𝑞₁𝑥𝑎+𝑏𝑞₁(),−𝑏𝑟1−𝑥/𝑘 (21) and null-isocline 𝑦(̇𝑡)=0 gives 𝑔_𝑆2(𝑥)=𝛽𝑎−𝛿−𝑞₂𝑥𝑏𝛿−𝑞₂. (22)

4.2. The Equilibria of Filippov System (3)

Because the solutions of Filippov system (3) are composed of connecting standard solutions in subsystems 𝑆₁, 𝑆₂ and sliding mode solutions on Σ¹_𝑆 or Σ²_𝑆. The definitions of all types of equilibria of Filippov system have been provided in Section 2, which are important to bifurcation analysis. There may be several types of equilibria for Filippov system (3) which include regular equilibrium (denoted by 𝐸_𝑅), virtual equilibrium (denoted by 𝐸_𝑉), pseudoequilibrium (denoted by 𝐸_𝑃), boundary equilibrium (denoted by 𝐸_𝐵), and one type of special point named as tangent point (denoted by 𝐸_𝑇). Detailed definitions of these equilibria and the tangent point can be found in the literature [24].

From Lemma 1, we know that the subsystem 𝑆₁ with 𝑥<ET has three possible equilibria: (0,0), (𝑘,0), and interior equilibrium, 𝐸^∗_𝑆1=𝑥^∗_𝑆1,𝑦^∗_𝑆1=𝑘𝛽𝑏𝑟+𝛿−𝛽𝑎,𝑘𝛽𝑏𝑟𝛽𝑎−𝛿𝛽𝑏𝑟+𝛿−𝛽𝑎𝛽𝑏².𝑟𝛿 (23)

If 𝛽𝑎−𝛿>0 and 𝑟>(𝛽𝑎−𝛿)/𝛽𝑏, then there is a positive interior equilibrium for the system 𝑆₁.

Equilibria for the subsystem 𝑆₂ with 𝑥>ET satisfy 𝑥𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝑞₁𝑥=0,−𝛿𝑦+𝛽𝑎𝑥𝑦𝑏𝑦+𝑥+𝑞₂𝑦=0. (24)

Solving the above equations with respect to 𝑥 and 𝑦 yields three possible equilibria, that is, (0,0), (𝑘(𝑟−𝑞₁)/𝑟,0), and interior equilibrium, 𝐸^∗_𝑆2=𝑥^∗_𝑆2,𝑦^∗_𝑆2=𝑘𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎,𝑘𝛽𝑏𝑟𝛽𝑎−𝛿+𝑞₂𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎𝛽𝑏²𝑟𝛿−𝑞₂, (25) which indicates that, if 𝑟>(𝛽𝑎−𝛿+𝑞₂)/𝛽𝑏+𝑞₁ and 𝑞₂<𝛿<𝛽𝑎+𝑞₂, then there is a positive interior equilibrium for the system 𝑆₂.

4.2.1. Regular Equilibria

For the subsystem 𝑆₁ with 𝑥<ET, (0,0) is a regular equilibrium, while (𝑘,0) is a virtual equilibrium. In addition, according to the coordinate of equilibrium 𝐸^∗_𝑆1, we have the following results. If 𝑟>(𝛽𝑎−𝛿)/𝛽𝑏>0 and 𝑘𝛽𝑏𝑟+𝛿−𝛽𝑎𝛽𝑏𝑟<ET, (26) then it is a regular equilibrium for the system 𝑆₁, denoted by 𝐸¹_𝑅. If 𝑟>(𝛽𝑎−𝛿)/𝛽𝑏>0 and 𝑘𝛽𝑏𝑟+𝛿−𝛽𝑎𝛽𝑏𝑟>ET, (27) then the equilibrium 𝐸^∗_𝑆1 becomes a virtual equilibrium, denoted by 𝐸¹_𝑉.

About regular equilibria for the subsystem 𝑆₂ with 𝑥>ET,  (0,0) is a virtual equilibrium, while (𝑘(𝑟−𝑝₁)/𝑟,0) may be a regular or virtual equilibrium which depends on the parameter space. Moreover, according to the coordinate of equilibrium 𝐸^∗_𝑆2, we have the following conclusions. If 𝑟>(𝛽𝑎−𝛿+𝑞₂)/𝛽𝑏+𝑞₁, 𝑞₂<𝛿<𝛽𝑎+𝑞₂ and 𝑘𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎𝛽𝑏𝑟>ET, (28) then it is a regular equilibrium for the system 𝑆₂, denoted by 𝐸²_𝑅. Note that if 𝑟>(𝛽𝑎−𝛿+𝑞₂)/𝛽𝑏+𝑞₁, 𝑞₂<𝛿<𝛽𝑎+𝑞₂ and 𝑘𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎𝛽𝑏𝑟<ET, (29) then the equilibrium 𝐸^∗_𝑆2 becomes a virtual equilibrium, denoted by 𝐸²_𝑉. Based on the above analysis, if 𝑘𝛽𝑎−𝛿()𝑘𝛽𝑏𝑘−ET<𝑟<𝛽𝑎−𝛿+𝑞₂+𝛽𝑏𝑞₁(),𝛽𝑏𝑘−ET (30) then the two virtual equilibria 𝐸¹_𝑉 and 𝐸²_𝑉 can coexist.

4.2.2. Pseudoequilibrium

For the existence of pseudoequilibrium 𝐸_𝑃=(ET,𝑦_𝑃), 𝑦_𝑃 component of the pseudoequilibrium of sliding flow satisfies the following equation: 𝑃₁𝑦𝑃₂𝑦+𝑃₃=0, (31) where 𝑦∈Σ¹_𝑆 or Σ²_𝑆. Solving the above equation with respect to 𝑦 yields two possible roots denoted by 𝑦_𝑃1=0 and 𝑦_𝑃2=𝛿𝑞₁−𝛽𝑎𝑞₁−𝑟𝑞₂(1−ET/𝑘)ET−𝛿𝑏𝑞₁+𝑟𝑏𝑞₂(1−ET/𝑘)−𝑎𝑞₂. (32)

Further, if 𝑘𝑞₁𝛿−𝛽𝑎𝑞₂()𝑘𝑘−ET<𝑟<𝑎𝑞₂+𝛿𝑏𝑞₁𝑏𝑞₂(𝑘−ET) (33) holds true, then 𝑦_𝑃2 is a positive root. Note that if 𝐸_𝑃1=(ET,𝑦_𝑃1) or 𝐸_𝑃2=(ET,𝑦_𝑃2) lies in the sliding region Σ¹_𝑆 or Σ²_𝑆, then the model has pseudoequilibrium. To do this, we consider the following two cases.

Case 1 (sliding segment defined by Σ¹_𝑆). If the inequality 𝑟<𝑞₁𝑘/(𝑘−ET) holds true, then 𝐸_𝑃1=(ET,𝑦_𝑃1)∈Σ¹_𝑆 is a pseudoequilibrium of the Filippov system (3). Note that −𝛿𝑏𝑞₁+𝑟𝑏𝑞₂(1−ET/𝑘)−𝑎𝑞₂<0 (i.e., 𝑟<𝑘(𝑎𝑞₂+𝛿𝑏𝑞₁)/𝑏𝑞₂(𝑘−ET)) is well defined in this case. Therefore if the inequalities max𝑘𝑞₁𝛿−𝛽𝑎𝑞₂(),𝑘𝑘−ET𝛽𝑎−𝛿()𝑘𝛽𝑏𝑘−ET<𝑟<min𝛽𝑎−𝛿+𝑞₂+𝛽𝑏𝑞₁(),𝑘𝛽𝑏𝑘−ET𝑎𝑞₂+𝛿𝑏𝑞₁𝑏𝑞₂(𝑘−ET) (34) hold true, then 𝐸_𝑃2=(ET,𝑦_𝑃2)∈Σ¹_𝑆 is a positive pseudoequilibrium of the Filippov system (3).

Case 2 (sliding segment defined by Σ²_𝑆). In this case, 𝐸_𝑃1=(ET,𝑦_𝑃1)∈Σ²_𝑆 is a pseudoequilibrium of the Filippov system (3). If the inequalities 𝑘max𝛽𝑎−𝛿+𝑞₂+𝛽𝑏𝑞₁(),𝛽𝑏𝑘−ET𝑘𝑞₁𝛿−𝛽𝑎𝑞₂()<𝑘𝑘−ET≤𝑟𝑎𝑞₂+𝛿𝑏𝑞₁𝑏𝑞₂(𝑘−ET) (35) hold true, 𝐸_𝑃2=(ET,𝑦_𝑃2)∈Σ²_𝑆 is a positive pseudoequilibrium of the Filippov system (3).

4.2.3. Boundary Equilibrium

The boundary equilibria of Filippov system (3) satisfy 𝑥𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝜀𝑞₁𝑥=0,−𝛿𝑦+𝛽𝑎𝑥𝑦𝑏𝑦+𝑥+𝜀𝑞₂𝑦=0,𝑥=ET, (36) with 𝜀=0 or 1, which indicate that, if 𝑘𝛽𝑏𝑟+𝛿−𝛽𝑎𝑘𝛽𝑏𝑟=ETor𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎𝛽𝑏𝑟=ET, (37) then we have the boundary equilibria 𝐸¹_𝐵=ET,𝛽𝑎−𝛿ET,𝐸𝑏𝛿_2𝐵=ET,𝛽𝑎−𝛿+𝑞₂ET𝑏𝛿−𝑞₂. (38)

4.2.4. Tangent Point

According to the definition of tangent point, we can see that the tangent point 𝐸_𝑇=(ET,𝑦_𝑇) on sliding segment Σ_𝑆 satisfies 𝑥𝑟𝑥1−𝑘−𝑎𝑥𝑦𝑏𝑦+𝑥−𝜀𝑞₁𝑥=0,𝑥=ET. (39)

Solving the above equations with respect to 𝑥 and 𝑦 yields two tangent points, including 𝐸¹_𝑇=()ET,𝑟ET1−ET/𝑘(),𝐸𝑎−𝑏𝑟1−ET/𝑘_2𝑇=ET𝑟(1−ET/𝑘)−𝑞₁𝑎+𝑞₁().𝑏−𝑏𝑟1−ET/𝑘 (40)

If we fix all parameters, then the relations among null-isoclines, regular/virtual equilibria, pseudoequilibrium, and sliding segment are provided in Figure 1. Note that virtual equilibria of both systems 𝑆₁ and 𝑆₂ imply the existence of pseudoequilibrium, and we will prove this general result later.

Figure 1 

Notations for null-isoclines, regular/virtual equilibria, pseudoequilibrium, and sliding segment. The parameter values are fixed as follows: 𝑎=0.7, 𝑏=0.6, 𝑟=1.2, 𝑘=12, 𝛽=0.5, 𝛿=0.25, 𝑞1=0.25, 𝑞2=0.1, and ET=5.

4.3. The Stability of Pseudoequilibrium

In the process of pest control, we should apply all kinds of control strategies so as to prevent multiple pest outbreaks or make sure that the total density of the pest stabilizes at a desired level of ET. If the unique positive equilibrium of system 𝑆₁ and system 𝑆₂ is virtual simultaneously, then the sliding flow has a unique pseudoequilibrium. In order to realize this goal, we can choose a set of parameters such that all the equilibria of subsystems 𝑆₁ and 𝑆₂ are virtual equilibria, and the pseudoequilibria are globally stable, which have been widely used in pest control [26–28]. For example, if we fixed all parameter values as those in Figure 2, then both virtual equilibria coexist. Therefore, we address the stability of pseudoequilibrium 𝐸_𝑃2 in the following which is important to control pest.

Figure 2 

Regular/virtual equilibrium. The parameter values are fixed as follows: 𝑎=0.7, 𝑏=0.6, 𝑘=12, 𝛽=0.5, 𝛿=0.25, 𝑞1=0.25, and 𝑞2=0.1.

Theorem 6. Either the inequalities (34) or (35) hold true or the two virtual equilibria 𝐸¹_𝑉 and 𝐸²_𝑉 can coexist; then Filippov system (3) contains a positive pseudoequilibrium 𝐸_𝑃2. Regardless of which cases would occur, the pseudoequilibrium 𝐸_𝑃2 is locally stable with respect to sliding mode domain.

Proof. According to the conditions of Theorem 6 we see that the inequalities (30) hold true, which implies that the inequalities (34) are true. Based on the discussions about the existence of pseudoequilibrium in Section 4.2, we have concluded that if inequalities (34) or (35) hold true, then the system (3) contains a positive pseudoequilibrium 𝐸_𝑃2.
It follows from (18) that we have𝑑𝜙𝑦𝑑𝑦_{𝑦=𝑦𝑃2}=−𝛿𝑏²𝑘𝑞₁+𝑟𝑏²𝑞₂(𝑘−ET)−𝑎𝑏𝑘𝑞₂𝑦²_𝑃2𝑘𝑞₁𝑏𝑦_𝑃2+ET²+−2𝛿𝑏𝑘𝑞₁ET−2𝑎𝑘𝑞₂ET+2𝑟𝑏𝑞₂ET(𝑘−ET)𝑦_𝑃2𝑘𝑞₁𝑏𝑦_𝑃2+ET²+−𝛿𝑘𝑞₁ET²+𝑎𝛽𝑘𝑞₁ET²+𝑟𝑘𝑞₂ET²−𝑟𝑞₂ET³𝑘𝑞₁𝑏𝑦_𝑃2+ET²=𝑘𝑞₁𝛿−𝛽𝑎−𝑟𝑞₂(𝑘−ET)𝛿𝑏𝑘𝑞₁−𝑟𝑏𝑞₂(𝑘−ET)+𝑎𝑘𝑞₂𝑎𝑞₁𝑘²𝑏𝛽𝑞₁+𝑞₂. (41)

Therefore, if the inequalities (33) hold true, then 𝑑𝜙(𝑦)/𝑑𝑦<0 and 𝑦_𝑃2 is a positive root of (31). Note that the inequalities (34) and (35) indicate the inequalities (33). Thus, the positive pseudoequilibrium 𝐸_𝑃2 is locally stale if it exists.

It is difficult to directly prove the global stability of pseudoequilibrium 𝐸_𝑃2 in this case. Because we cannot employ the classical Bendixson-Dulac theorem due to the discontinuity of vector fields. However, if there is not crossing cycle surrounding sliding segment, then pseudoequilibrium is globally stable [29]. From the analysis of global bifurcation for the system (3) in the following section, the system just has touching bifurcation and there is nonexistence of a sliding cycle which surrounds the 𝐸_𝑃2. By using the similar methods, we have that pseudoequilibrium 𝐸_𝑃2 is globally stable. That is to say, the local stability of pseudoequilibrium 𝐸_𝑃2 with respect to sliding mode domain indicates its global stability in the first quadrant (shown in Figures 3(a), 3(b), 4(b), and 4(c)). In practice, in order to control pest outbreak, we should choose the desirable ET at first, so that all equilibria of each system such as system 𝑆₁ and system 𝑆₂ become virtual; then pseudoequilibrium not only exists but also is globally stable. In other words, the density of pest can be stable at the ET. When the density of pest reaches or exceeds the ET, we should carry out control strategies (e.g., releasing natural enemy, etc.), until it falls below the ET. In this way, our control goal can be realized fully.

(a) ET=8.3

(b) ET=9.3333

(c) ET=9.8

Figure 3 

Boundary node bifurcation for Filippov system (3). The parameter values are fixed as follows: 𝑎=0.7, 𝑏=0.6, 𝑟=1.5, 𝑘=12, 𝛽=0.5, 𝛿=0.25, 𝑞1=0.25, and 𝑞2=0.1.

(a) ET=4

(b) ET=4.6667

(c) ET=5.8

Figure 4 

Boundary node bifurcation for Filippov system (3). The parameter values are fixed as follows: 𝑎=0.7, 𝑏=0.6, 𝑟=1.5, 𝑘=12, 𝛽=0.5, 𝛿=0.25, 𝑞1=0.25, and 𝑞2=0.1.

5. Equilibria and Sliding Bifurcation Set

5.1. Regular/Virtual Equilibrium Bifurcation

According to the above discussions, it is obvious that 𝑟 and ET are primary factors in determining the existence of the above different types of equilibria of the system (3). So we define four curves about parameters 𝑟 and ET as follows: 𝐿₁=(𝑟,ET)∣𝑟=𝛽𝑎−𝛿,𝐿𝛽𝑏₂=(𝑟,ET)∣𝑟=𝛽𝑎−𝛿+𝑞₂𝛽𝑏+𝑞₁,𝐿₃=)𝑘(𝑟,ET∣ET=𝛽𝑏𝑟+𝛿−𝛽𝑎,𝐿𝛽𝑏𝑟₄=)𝑘(𝑟,ET∣ET=𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎.𝛽𝑏𝑟 (42)

The four curves (i.e., 𝐿₁, 𝐿₂, 𝐿₃, and 𝐿₄) divide the 𝑟 and ET parameter space into six regions, and the existence or coexistence of regular or virtual equilibria is indicated in each region. The boundary equilibria 𝐸¹_𝐵 and 𝐸²_𝐵 can appear on the lines 𝐿₃ and 𝐿₄ accordingly. In particular, it follows from Figure 2 that the two virtual equilibria 𝐸¹_𝑉 and 𝐸²_𝑉 can coexist, which is very important to pest control. However, the two regular equilibria 𝐸¹_𝑅 and 𝐸²_𝑅 cannot coexist.

5.2. Boundary Node Bifurcations

This type of bifurcation may occur for Filippov system (3) once equilibria 𝐸_𝑃, 𝐸_𝑇, and 𝐸_𝑅 or 𝐸_𝑇 and 𝐸_𝑅 collide together simultaneously when ET passes through a critical value. In this part, we choose ET as bifurcation parameter, and all other parameters are fixed as those in Figure 3. Note that once the parameter ET passes through the first critical value ET_𝑐1=9.3333, the regular equilibrium 𝐸¹_𝑅, tangent point 𝐸_𝑇, and pseudoequilibrium 𝐸_𝑃 collide together (see Figure 3(b)), where ET_𝑐1 is determined by ET_𝑐1=𝑘𝛽𝑏𝑟+𝛿−𝛽𝑎.𝛽𝑏𝑟 (43)

A virtual equilibrium 𝐸¹_𝑉, an invisible tangent point 𝐸¹_𝑇 and a pseudoequilibrium 𝐸_𝑃 coexist, as shown in Figure 3(a), when ET<ET_𝑐1. They collide at ET=ET_𝑐1 and are substituted by a visible tangent point 𝐸¹_𝑇, as shown in Figure 3(c), when ET>ET_𝑐1.

Similarly, another boundary node bifurcation of Filippov system (3) occurs at ET_𝑐2=4.6667 (see Figure 4(b)), where ET_𝑐2 is determined by ET_𝑐2=𝑘𝛽𝑏𝑟−𝑞₁+𝛿−𝑞₂−𝛽𝑎.𝛽𝑏𝑟 (44)

A stable regular equilibrium 𝐸²_𝑅 and a visible tangent point 𝐸²_𝑇 coexist as shown in Figure 4(a) when ET<ET_𝑐2. They collide with a pseudoequilibrium 𝐸_𝑃 at ET=ET_𝑐2 and are substituted by an invisible tangent point 𝐸²_𝑇 as shown in Figure 4(c) when ET>ET_𝑐2.

5.3. Global Sliding Bifurcation

Global sliding bifurcations involve nonvanishing cycles, which include sliding disconnection, touching (or grazing) bifurcation, buckling bifurcation, crossing bifurcation, bifurcation of a sliding homoclinic orbit, and heteroclinic orbit [22]. Touching (or grazing) bifurcation implies that a positive period solution can collide with the sliding segments. From the work of Kuang and Beretta and Hsu et al. [2, 7], they concluded that when the value of the parameter 𝑎 passed slightly through the bifurcation value, a stable limit cycle bifurcates from the unstable positive interior equilibrium for the system (1). According to numerical simulations, the system (3) just has touching bifurcation. For example, if we choose ET as bifurcation parameter and fix all other parameters as shown in Figure 5, when the value of parameter ET varies, touching bifurcation occurs at the critical value ET_𝑐3=2.59 for the system (3). Note that extensive numerical simulations indicate that, no matter what the value of ET is in touching bifurcation, the whole periodic solution lies in the region 𝑆₁ (shown in Figure 5). This implies that the density of pest can be successfully maintained below ET by designing suitable threshold policy strategies. So our control objective can be fully realized, which is a desired situation in crop, livestock sectors and forestry.

(a) ET=1.5

(b) ET=2.59

(c) ET=3.5

Figure 5 

Globally touching bifurcation for Filippov system (3). The parameter values are fixed as follows: 𝑎=0.51, 𝑏=0.4, 𝑟=1.2, 𝑘=6, 𝛽=0.6863, 𝛿=0.05, 𝑞1=0.4, and 𝑞2=0.1.

6. Discussion

Recently, the threshold policy and IPM strategies have attracted great attention in agriculture, forestry, animal husbandry, and so on [10–12]. In the process of pest control, IPM strategies would be used only when the density of pest reaches or exceeds the ET. In addition, Filippov system provides a natural and rational framework for those real world problems, so it has been widely used in different fields such as in science and engineering [19–21, 23]. In this paper, we employ the sliding analysis of Filippov system to describe and investigate the long term dynamical behavior of the ratio-dependent-type prey-predator model with Michaelis-Menten-type functional response, so that we can use the Filippov system to model intervention of pest control policy. Firstly, we investigate the sliding mode domain, sliding mode dynamics. Secondly, the null-isoclines and the existence of four types of equilibria for Filippov system, including regular, virtual, boundary, and pseudoequilibrium and the tangent points are discussed in detail. Moreover, the stability of pseudoequilibrium is also studied. Thirdly, we have investigated the local sliding bifurcations including regular/virtual equilibrium bifurcation, and boundary node bifurcations. Further, global touching bifurcation is also studied by numerical techniques.

In the process of pest control, we should apply all kinds of control strategies so as to prevent multiple pest outbreaks or make sure that the density of pest stabilizes at a desired level of ET. In order to realize these goals, on the one hand, we can determine a set of parameters such that not only all the equilibria of subsystems 𝑆₁ and 𝑆₂ are virtual equilibria, which ensures that the pseudoequilibria exist, but also the pseudoequilibrium is globally stable. From Theorem 6, we provided the conditions of the existence of pseudoequilibrium 𝐸_𝑃2 for the system (3) and showed that it is locally stable if it exists. Further, we showed that the local stability of pseudoequilibrium 𝐸_𝑃2 implies global stability by numerical simulation. On the other hand, globally touching bifurcation indicates that the density of pest can be successfully maintained below the ET by designing suitable threshold policy strategies (shown in Figure 5). Therefore our control objective can be achieved fully in the above two cases, which can be used to pest control in crop, livestock sectors and forestry.

Although impulsive prey-predator models with the ratio-dependent- or Michaelis-Menten-type response function have been studied in [13–16], the Filippov systems have many advantages in describing interventions including spraying pesticides and releasing natural enemies compared with impulsive models. In this work, the number of natural enemy to be released is proportional to its number. It is interesting that releasing number of natural enemy can be described by constant, independently of the existing numbers of pest and natural enemy. Moreover, in practice, considering environmental energy resource finiteness, we should choose the total number of both populations as a guide to switch the system, which is called the weighted escapement policy (WEP) [19, 30–32]. Therefore, in the future work, we will focus on the above two cases, which could result in richer dynamics.

Appendix

Methods for Analysis of the Sliding Solution

Filippov Convex Method. The Filippov method associates the following convex combination 𝐹_𝑆(𝑍) of the two vectors 𝐹_𝑆1(𝑍) and 𝐹_𝑆2(𝑍) with each nonsingular sliding point 𝑍∈Σ₂; that is,𝐹_𝑆(𝑍)=(1−𝜆(𝑍))𝐹_𝑆1(𝑍)+𝜆(𝑍)𝐹_𝑆2(𝑍), (A.1) where 𝜆(𝑍)=𝐻_𝑍,𝐹_𝑆1𝐻_𝑍,𝐹_𝑆1−F_𝑆2, (A.2)𝐹_𝑆(𝑍) is tangent to Σ₂ and 0≤𝜆(𝑍)≤1.

Thus, the sliding mode dynamics can be determined by 𝑍̇=𝐹_𝑆(𝑍),𝑍∈Σ₂, (A.3) which is smooth on a one-dimensional sliding interval of Σ₂. The solution of the above equation is the sliding solution.

The equation 𝜆(𝑍)=0 (or 𝜆(𝑍)=1) indicates that the flow is determined by 𝐹_𝑆1 (or 𝐹_𝑆2) alone. Therefore the sliding mode domain can be defined as Σ_𝑆=𝑍=𝑥,𝑦∈Σ∣0≤𝜆(𝑍)≤1, (A.4) which is equivalent to Σ_𝑆={𝑍∈Σ∣𝜎(𝑍)≤0}. Denote Σ⁺_𝑆={𝑍∈Σ∣𝜆(𝑍)=1},Σ⁻_𝑆={𝑍∈Σ∣𝜆(𝑍)=0}, (A.5) by the boundary of the sliding mode domain, so the vector fields are tangent to boundary Σ⁺_𝑆 or Σ⁻_𝑆.

Utkin Equivalent Control Method. It follows from Utkin's works [18] on sliding mode dynamics along the manifold Σ; we note that Filippov system (3) can be also rewritten as follows: 𝑍̇=𝐹𝑍,𝑈_𝐻, (A.6) where the control 𝑈_𝐻 is defined as𝑈_𝐻=(𝑍)𝑈0,𝐻<0,₁(𝑍,𝑡),𝐻(𝑍)>0, (A.7) where 𝑈₁ is a continuous function. In the controlled system (i.e., system 𝑆₂), the control 𝑈_𝐻=𝑈₁ is applied, and in the free system (i.e., system 𝑆₁) 𝑈_𝐻=0 is applied. Assume that a sliding mode exists on manifold Σ; that is, Σ₂ is nonempty. Solving the following algebraic equation: 𝐻̇=𝜕𝐻𝐹𝜕𝑍𝑍,U_𝐻=0 (A.8) with respect to 𝑈_𝐻 on this manifold, gives the solution denoted by 𝑈^∗(𝑍,𝑡), which is referred to as equivalent control. Substituting for 𝑈_𝐻 in system (A.6) yields 𝑍̇=𝑓𝑍,𝑈^∗(𝑍,𝑡),𝑍∈Σ₂ (A.9) which determines the sliding mode dynamics of the Filippov system (3).

Acknowledgments

This work is supported by the National Natural Science Foundation of China (NSFC, 11171199) and by the Fundamental Research Funds for the Central Universities (GK201003001).