# A Holling Type II Discrete Switching Host-Parasitoid System with a Nonlinear Threshold Policy for Integrated Pest Management

Authors: Mengqi He; Sanyi Tang; Robert A. Cheke
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

A nonlinear discrete switching host-parasitoid model with Holling type II functional response function, in which the switch is guided by an economic threshold (ET), is proposed. Thus, if the weighted density of two generations of the host population increases and exceeds the ET, then integrated pest management (IPM) measures are enacted, i.e., biological and chemical measures are implemented together, assuming that the chemical immediately precedes the biological inputs to avoid pesticide-induced deaths of the natural enemies. First, the existence and local stability of the equilibria of two subsystems were studied, and the existence and coexistence of several types of equilibria of a nonlinear switching system were analysed. Next, the nonlinear switching system was investigated by numerical simulation, showing that the system exhibits quite complex dynamic behaviour. A two-dimensional bifurcation diagram revealed the existence and coexistence regions of different types of equilibria including regular and virtual equilibria. Moreover, period-adding bifurcations in two-dimensional parameter spaces were found. One-dimensional bifurcation diagrams revealed that the system has periodic, quasiperiodic, and chaotic solutions, Neimark–Sacker bifurcation, multiple coexisting attractors, period-doubling bifurcations, period-halving bifurcations, and so on. Finally, the initial densities of hosts and parasitoids associated with host outbreaks and their biological implications are discussed.


---

## Body

## 1. Introduction

Pest control is important in subjects such as agriculture, fisheries, and ecology [1]. Pest outbreaks often cause huge economic losses and a series of ecological problems. For example, the outbreak of the fall armyworm moth Spodoptera frugiperda in 2019 has caused considerable economic losses in many areas of southern China [2]. Rational control of pests avoids indiscriminate spraying of pesticides, as long-term pesticide usage will not only induce resistance but also cause serious environmental pollution. In order to control pests more effectively and to consider the impact of control measures on the environment, the integrated pest management (IPM) strategy was proposed [3–6], which requires the combination of chemical control and biological control. There are some important concepts in IPM, including the economic injury level (EIL) [7–9], economic threshold (ET) [8–10], and threshold policy control (TPC) [11–13]. The purpose of IPM is not to eliminate pests completely, but to control them at levels below the EIL. Thus, the classical method is to maintain the density of a pest population below the EIL by spraying pesticides and releasing natural enemies and/or other measures such as cultural control once the density of the pest population reaches the ET.Host-parasitoid models including IPM intervention with fixed and unfixed time pulse effects were first proposed and analysed by Tang et al. [8]. Moreover, Yang et al. [14] proposed and analysed a Holling type II host-parasitoid model with IPM intervention as a pulse control strategy, including models with fixed pulses and unfixed pulses. Wang and Zhang [15] studied the Holling type II host-parasitoid model with chemical control at a fixed time. In order to apply such a threshold policy control (TPC) to IPM measures, a switching system (or Filippov system) is proposed and studied. Such systems, based on the ordinary differential equation (ODE) model for the predator-prey interaction, have been investigated in recent publications (see [12, 13, 16–18] for relevant references). In addition, Xiang and his co-authors [19–24] discussed the complex dynamic behaviour of discrete switching systems. However, in the above publications, the authors only considered the relationship between ET and the density of a pest population in generation t, that is, the threshold condition is linear or the switching surface is a straight line.For a pest population with nonoverlapping generations, the best control threshold depends not only on the density of the current generation but also on the density of the next generation. This means that chemical control and biological control applied together (but with the chemical control immediately preceding the biological inputs to avoid pesticide-induced deaths of the natural enemies) should be guided by the densities of both generationt and generation t+1. This means that the threshold condition is nonlinear, which we denote as a nonlinear switching system guided by the weighted density of two generations of the host population. To our knowledge, no previous research has been conducted on such a nonlinear discrete switching system.The structure of this paper is as follows. In Section2, a Holling type II host-parasitoid model with a nonlinear switching threshold strategy is proposed. In Section 3, we discuss the existence and local stability of the equilibria of the two subsystems. Numerical simulations including one-parameter bifurcation and two-parameter bifurcation are applied in Section 4 to reveal the complex and diverse dynamic behaviour of the nonlinear switching system. In Section 5, we study the importance and correlation of initial densities of both host and parasitoid populations in host outbreaks.


## 2. Nonlinear Switching Host-Parasitoid Model

A discrete host-parasitoid model with a Holling type II functional response function was numerically investigated by Tang and Chen [25], who showed that a discrete host-parasitoid model could generate much richer dynamics than those observed in continuous models. With H the host and P the parasitoid populations, this model is [25](1)Ht+1=Htexpr1−HtK−aTPt1+aThHt,Pt+1=Ht1−exp−aTPt1+aThHt,where all parameters are positive, r denotes the intrinsic growth rate of the host population in the absence of the parasitoid population, the carrying capacity is represented by K, a describes the instantaneous search rate (the average encounters per host per unit of search time), assuming that the total time initially available for the search when the hosts are exposed to parasitoids is T, and Th is the handling time (the time between the host being encountered and the search for a host being resumed).In order to take into account IPM strategies for controlling the pest population, Xiang et al. [21] considered that only when the density of generation t of the host population exceeds the ET, is there implementation of biological control measures (natural enemy releases) together with chemical control measures (pesticide spraying, but see comment above on earlier deployment of pesticides), which means that the switching surface is a straight line.But in reality, the optimal control threshold not only depends on the density of the current generation but also depends on the density of the next generation, which implies that the thresholdET will depend on both the host and parasitoid densities. The question is how to establish a mathematical model that uses the relationship between the population density of two adjacent generations and the threshold ET as the switching condition. In order to solve this problem, based on model (1), we consider a more generalized threshold control policy, i.e., once the weighted density of two generations of the host population exceeds the ET, the chemical control is formulated by assuming a proportional killing rate p for the pest population, and the biological control is modelled by releasing a constant number of the natural enemy τ, which obviously generates a nonlinear switching curve. This yields the following model:(2)Ht+1=1−pHtexpr1−HtK−aTPt1+aThHt,Pt+1=Ht1−exp−aTPt1+aThHt+τ,εHt+1−εHt+1≥ET,where 0≤p<1, τ≥0, 0≤ε≤1. In order to control the pest, we assume that if τ>0 and p=0, only biological control is adopted by the release of a sufficient number of the parasitoid to prevent the weighted host population exceeding ET, whereas if 0<p<1 and τ=0, chemical control is adopted to ensure that the weighted host population will be decreased once it reaches the ET. Note that if ε=1, then this special case was investigated by Xiang et al. [21].Therefore, combining model (2) (called a control system) with IPM when the weighted density of two generations of the host population exceeds the ET and model (1) (called a free system) without control measures when the weighted density of two generations of the host population falls below the ET, we have the following nonlinear switching system:(3)Ht+1=Htexpr1−HtK−aTPt1+aThHt,Pt+1=Ht1−exp−aTPt1+aThHt,εHt+1−εHt+1<ET;Ht+1=1−pHtexpr1−HtK−aTPt1+aThHt,Pt+1=Ht1−exp−aTPt1+aThHt+τ,εHt+1−εHt+1≥ET.Obviously, model (3) represents a switching system controlled by a nonlinear threshold policy and it differs from previous studies of a differential switching system in that the model switching surface is nonlinear. Note that the switching threshold between two subsystems depends not only on the density of current generation but also on the next generation, which results in a complex switching curve and could generate new dynamics related to controlling of host population. Moreover, once the weighted parameter ε=0 or 1, then the switching curve reduces to the straight line, which has been studied in literatures [20, 21].


## 3. Stability Analysis of Two Subsystems

In order to reveal the dynamics of the whole switching system (3), we first investigate the existence and stability of various equilibria of the control and free systems, providing results that are useful for studying the dynamical behaviour of system (3). For convenience, we denote φZ=εHt+1−εHt+1−ET with vector Z=Ht,PtT, and(4)SG1Z=Htexpr1−HtK−aTPt1+aThHt,Ht1−exp−aTPt1+aThHtT,SG2Z=1−pHtexpr1−HtK−aTPt1+aThHt,Ht1−exp−aTPt1+aThHt+τT.Then, model (3) can be rewritten as the following switching system:(5)Zt+1=SG1Z,Z∈G1,SG2Z,Z∈G2,where(6)G1=Z∈R+2φZ<0,Ht≥0,Pt≥0,G2=Z∈R+2φZ≥0,Ht≥0,Pt≥0.From now on, we call switching system (3) defined in region G1 as subsystem SG1 (i.e., system (1)) and that defined in region G2 as subsystem SG2 (i.e., system (2)).Switching systems have different types of equilibria, and these equilibria play important roles in pest control. In order to show the existence of equilibria of switching system (4), we need the definitions of the regular equilibrium and the virtual equilibrium about switching system (4).Definition 1.
A pointZ∗=H∗,P∗ is called a regular equilibrium of system (4) if Z∗ is an equilibrium of subsystem SG1(SG1Z∗=0) and φH∗,P∗<0 or if Z∗ is an equilibrium of subsystem SG2(SG2Z∗=0) and φH∗,P∗≥0. These equilibria will be denoted by ER1 and ER2. Similarly, a point Z∗ is called a virtual equilibrium of system (4) if Z∗ is an equilibrium of subsystem SG1 and φH∗,P∗≥0 or if Z∗ is an equilibrium of subsystem SG2 and φH∗,P∗<0. These equilibria will be denoted by EV1 and EV2.


### 3.1. Local Stability of Equilibria for SubsystemSG1

IfεHt+1−εHt+1<ET, then the dynamics of system (4) have been determined by subsystem SG1. The equilibrium H∗,P∗ of subsystem SG1 satisfies the following equations:(7)H∗=H∗expr1−H∗K−aTP∗1+aThH∗,P∗=H∗1−exp−aTP∗1+aThH∗.Obviously, subsystemSG1 has the trivial equilibrium E01=0,0.Proposition 1.
The trivial equilibriumE01 is an unstable saddle.Proof.
The stability of0,0 of subsystem SG1 is determined by the eigenvalues of the following Jacobian matrix:(8)J0=er000.
The eigenvalues ofJ0 are given by λ1=er and λ2=0. Since the parameter r is positive, λ1>1,λ2<1. Hence, the trivial equilibrium E01 is an unstable saddle.
In addition, there exists the parasitoid-free equilibriumEb=K,0. It was shown by Kaitala et al. [26] that the boundary equilibrium is locally asymptotically stable when 1−r<1 and aTK/1+aThK<1. Meanwhile, when R0>1 and TR0−1>ThR0lnR0, there exists the unique positive equilibrium E1∗=H1∗,P1∗, where E1∗ satisfies the following equations:(9)H1∗=R0lnR0aTR0−1−ThR0lnR0,P1∗=R0−1lnR0aTR0−1−ThR0lnR0.
Here,R0 denotes the net growth rate of the host per generation with R0=expr1−H1∗/K and reveals that the equilibrium E1∗=H1∗,P1∗ cannot be solved analytically.
In the following, we address the local stability of the unique positive equilibriumE1∗=H1∗,P1∗ of subsystem SG1.Theorem 1.
The unique positive equilibriumE1∗=H1∗,P1∗ of subsystem SG1 is asymptotically stable provided that(10)max−1−trJ1,−1+trJ1<detJ1<1,where(11)J1=∂F1∂H∂F1∂P∂G1∂H∂G1∂HH1∗,P1∗.Proof.
In order to discuss the local stability of the equilibrium, we need to calculate the Jacobian matrix at the equilibrium. Firstly, we can rewrite subsystemSG1 as follows:(12)Ht+1=F1Ht,Pt,Pt+1=G1Ht,Pt.
Therefore, we have(13)J1=∂F1∂H∂F1∂P∂G1∂H∂G1∂PH1∗,P1∗,where(14)∂F1∂HH1∗,P1∗=1−rR0MK+a2TThR0−1R0M21+aThR0M2,∂F1∂PH1∗,P1∗=−aTR0M1+aThR0M,∂G1∂HH1∗,P1∗=1−1+a2TThR0−1R0M21+aThR0M2exp−aTR0−1M1+aThR0M,∂G1∂PH1∗,P1∗=aTR0M1+aThR0Mexp−aTR0−1M1+aThR0M,with(15)M=lnR0aTR0−1−ThR0lnR0,and then we can obtain the characteristic equation about the Jacobian matrix J1 of subsystem SG1:(16)ℱλ=λ2−trJ1λ+detJ1,where(17)trJ1=1−rR0MK+a2TThR0−1R0M21+aThR0M2+aTR0M1+aThR0Mexp−aTR0−1M1+aThR0M,(18)detJ1=1−rR0MKexp−aTR0−1M1+aThR0MaTR0M1+aThR0M.
Conditions (5) ensure that all the inequalities(19)ℱ1>0,ℱ−1>0,detJ1<1,hold, that is, the modulus of all roots of equation (16) is less than 1. By employing the Jury condition, we can derive that the equilibrium point H1∗,P1∗ of subsystem SG1 is local asymptotically stable. The proof is completed.


### 3.2. Local Stability of Equilibria for SubsystemSG2

IfεHt+1−εHt+1≥ET, then the dynamical behaviour of system (4) has been determined by subsystem SG2. Similarly, we can discuss the existence and stability of equilibria of subsystem SG2. Obviously, SG2 has a host-free equilibrium Eb′=0,τ if τ>0. Moreover, if(20)R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,then subsystem SG2 has a unique positive equilibrium E2∗=H2∗,P2∗, where(21)H2∗=1−pR1ln1−pR1−τaTaT1−pR1−T−Th1−pR1ln1−pR1,P2∗=1−pR1−1ln1−pR1−τaTaT1−pR1−T−Th1−pR1ln1−pR1+τ.Similarly,R1 represents the net growth rate of the host per generation with R1=expr1−H2∗/K. Similarly, the equilibrium E2∗=H2∗,P2∗ cannot be analytically solved. In the following theorems, the local stability of the host-free equilibrium 0,τ and the positive equilibrium H2∗,P2∗ of subsystem SG2 are proved.Theorem 2.
The host-free equilibrium0,τ of subsystem SG2 is asymptotically stable provided that(22)1−per−τaT<1.Proof.
The stability of0,τ of subsystem SG2 is determined by the eigenvalues of the following Jacobian matrix:(23)J=1−per−τaT01−e−τaT0.
Accordingly, we find that the eigenvalues areλ1=1−per−τaT,λ2=0. From this, 0,τ is asymptotically stable only when 1−per−τaT<1. The proof is completed.
For the stability of interior equilibriumE2∗=H2∗,P2∗, we denote the Jacobian matrix as follows:(24)J2=∂F2∂H∂F2∂P∂G2∂H∂G2∂PH2∗,P2∗,with(25)∂F2∂HH2∗,P2∗=1−r1−pR1QK+a2TTh1−pR1Q1−pR1−1Q+τ1+aTh1−pR1Q2,∂F2∂PH2∗,P2∗=−aT1−pR1Q1+aTh1−pR1Q,∂G2∂HH2∗,P2∗=1−1+a2TTh1−pR1Q1−pR1−1Q+τ1+aTh1−pR1Q2exp−aT1−pR1−1Q+τ1+aTh1−pR1Q,∂G2∂PH2∗,P2∗=aT1−pR1Q1+aTh1−pR1Qexp−aT1−pR1−1Q+τ1+aTh1−pR1Q,Q=ln1−pR0−τaTaT1−pR0−T−Th1−pR0ln1−pR0.Theorem 3.
The unique positive equilibriumE2∗=H2∗,P2∗ of subsystem SG2 is asymptotically stable if and only if(26)max−1−trJ2,−1+trJ2<detJ2<1.The detailed mathematical expressions oftrJ2 and detJ2 can be given as follows:(27)trJ2=1−r1−pR1QK+a2TTh1−pR1Q1−pR1−1Q+τ1+aTh1−pR1Q2+aT1−pR1Q1+aTh1−pR1Qexp−aT1−pR1−1Q+τ1+aTh1−pR1Q,(28)detJ2=1−r1−pR1QKexp−aT1−pR1−1Q+τ1+aTh1−pR1QaT1−pR1Q1+aTh1−pR1Q.The proof of Theorem3 is similar to Theorem 1, and so we omit it.


## 3.1. Local Stability of Equilibria for SubsystemSG1

IfεHt+1−εHt+1<ET, then the dynamics of system (4) have been determined by subsystem SG1. The equilibrium H∗,P∗ of subsystem SG1 satisfies the following equations:(7)H∗=H∗expr1−H∗K−aTP∗1+aThH∗,P∗=H∗1−exp−aTP∗1+aThH∗.Obviously, subsystemSG1 has the trivial equilibrium E01=0,0.Proposition 1.
The trivial equilibriumE01 is an unstable saddle.Proof.
The stability of0,0 of subsystem SG1 is determined by the eigenvalues of the following Jacobian matrix:(8)J0=er000.
The eigenvalues ofJ0 are given by λ1=er and λ2=0. Since the parameter r is positive, λ1>1,λ2<1. Hence, the trivial equilibrium E01 is an unstable saddle.
In addition, there exists the parasitoid-free equilibriumEb=K,0. It was shown by Kaitala et al. [26] that the boundary equilibrium is locally asymptotically stable when 1−r<1 and aTK/1+aThK<1. Meanwhile, when R0>1 and TR0−1>ThR0lnR0, there exists the unique positive equilibrium E1∗=H1∗,P1∗, where E1∗ satisfies the following equations:(9)H1∗=R0lnR0aTR0−1−ThR0lnR0,P1∗=R0−1lnR0aTR0−1−ThR0lnR0.
Here,R0 denotes the net growth rate of the host per generation with R0=expr1−H1∗/K and reveals that the equilibrium E1∗=H1∗,P1∗ cannot be solved analytically.
In the following, we address the local stability of the unique positive equilibriumE1∗=H1∗,P1∗ of subsystem SG1.Theorem 1.
The unique positive equilibriumE1∗=H1∗,P1∗ of subsystem SG1 is asymptotically stable provided that(10)max−1−trJ1,−1+trJ1<detJ1<1,where(11)J1=∂F1∂H∂F1∂P∂G1∂H∂G1∂HH1∗,P1∗.Proof.
In order to discuss the local stability of the equilibrium, we need to calculate the Jacobian matrix at the equilibrium. Firstly, we can rewrite subsystemSG1 as follows:(12)Ht+1=F1Ht,Pt,Pt+1=G1Ht,Pt.
Therefore, we have(13)J1=∂F1∂H∂F1∂P∂G1∂H∂G1∂PH1∗,P1∗,where(14)∂F1∂HH1∗,P1∗=1−rR0MK+a2TThR0−1R0M21+aThR0M2,∂F1∂PH1∗,P1∗=−aTR0M1+aThR0M,∂G1∂HH1∗,P1∗=1−1+a2TThR0−1R0M21+aThR0M2exp−aTR0−1M1+aThR0M,∂G1∂PH1∗,P1∗=aTR0M1+aThR0Mexp−aTR0−1M1+aThR0M,with(15)M=lnR0aTR0−1−ThR0lnR0,and then we can obtain the characteristic equation about the Jacobian matrix J1 of subsystem SG1:(16)ℱλ=λ2−trJ1λ+detJ1,where(17)trJ1=1−rR0MK+a2TThR0−1R0M21+aThR0M2+aTR0M1+aThR0Mexp−aTR0−1M1+aThR0M,(18)detJ1=1−rR0MKexp−aTR0−1M1+aThR0MaTR0M1+aThR0M.
Conditions (5) ensure that all the inequalities(19)ℱ1>0,ℱ−1>0,detJ1<1,hold, that is, the modulus of all roots of equation (16) is less than 1. By employing the Jury condition, we can derive that the equilibrium point H1∗,P1∗ of subsystem SG1 is local asymptotically stable. The proof is completed.


## 3.2. Local Stability of Equilibria for SubsystemSG2

IfεHt+1−εHt+1≥ET, then the dynamical behaviour of system (4) has been determined by subsystem SG2. Similarly, we can discuss the existence and stability of equilibria of subsystem SG2. Obviously, SG2 has a host-free equilibrium Eb′=0,τ if τ>0. Moreover, if(20)R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,then subsystem SG2 has a unique positive equilibrium E2∗=H2∗,P2∗, where(21)H2∗=1−pR1ln1−pR1−τaTaT1−pR1−T−Th1−pR1ln1−pR1,P2∗=1−pR1−1ln1−pR1−τaTaT1−pR1−T−Th1−pR1ln1−pR1+τ.Similarly,R1 represents the net growth rate of the host per generation with R1=expr1−H2∗/K. Similarly, the equilibrium E2∗=H2∗,P2∗ cannot be analytically solved. In the following theorems, the local stability of the host-free equilibrium 0,τ and the positive equilibrium H2∗,P2∗ of subsystem SG2 are proved.Theorem 2.
The host-free equilibrium0,τ of subsystem SG2 is asymptotically stable provided that(22)1−per−τaT<1.Proof.
The stability of0,τ of subsystem SG2 is determined by the eigenvalues of the following Jacobian matrix:(23)J=1−per−τaT01−e−τaT0.
Accordingly, we find that the eigenvalues areλ1=1−per−τaT,λ2=0. From this, 0,τ is asymptotically stable only when 1−per−τaT<1. The proof is completed.
For the stability of interior equilibriumE2∗=H2∗,P2∗, we denote the Jacobian matrix as follows:(24)J2=∂F2∂H∂F2∂P∂G2∂H∂G2∂PH2∗,P2∗,with(25)∂F2∂HH2∗,P2∗=1−r1−pR1QK+a2TTh1−pR1Q1−pR1−1Q+τ1+aTh1−pR1Q2,∂F2∂PH2∗,P2∗=−aT1−pR1Q1+aTh1−pR1Q,∂G2∂HH2∗,P2∗=1−1+a2TTh1−pR1Q1−pR1−1Q+τ1+aTh1−pR1Q2exp−aT1−pR1−1Q+τ1+aTh1−pR1Q,∂G2∂PH2∗,P2∗=aT1−pR1Q1+aTh1−pR1Qexp−aT1−pR1−1Q+τ1+aTh1−pR1Q,Q=ln1−pR0−τaTaT1−pR0−T−Th1−pR0ln1−pR0.Theorem 3.
The unique positive equilibriumE2∗=H2∗,P2∗ of subsystem SG2 is asymptotically stable if and only if(26)max−1−trJ2,−1+trJ2<detJ2<1.The detailed mathematical expressions oftrJ2 and detJ2 can be given as follows:(27)trJ2=1−r1−pR1QK+a2TTh1−pR1Q1−pR1−1Q+τ1+aTh1−pR1Q2+aT1−pR1Q1+aTh1−pR1Qexp−aT1−pR1−1Q+τ1+aTh1−pR1Q,(28)detJ2=1−r1−pR1QKexp−aT1−pR1−1Q+τ1+aTh1−pR1QaT1−pR1Q1+aTh1−pR1Q.The proof of Theorem3 is similar to Theorem 1, and so we omit it.


## 4. Complex Dynamical Behaviour Analysis

In contrast to two subsystems (1) and (2), model (4) characterizes the nonlinear switching threshold policy, and its dynamical behaviour may be more complex compared with results obtained in the literature [25]. The nonlinear interaction between the host and the parasitoid populations and the nonlinear switching threshold is very difficult or even impossible to analyse theoretically, so we will study the dynamical behaviour of system (4) through numerical investigations.


### 4.1. Equilibria Bifurcation for Nonlinear Switching System (4)

According to Definition1, nonlinear switching system (4) has different types of equilibria, and these equilibria play key roles in disease and pest control. For example, Xiao et al. [27] investigated the control of emerging infectious disease and Tang et al. [16] addressed pest control with an ET. In the following, we focus on how the system with a nonlinear switching threshold affects the distribution of regular and virtual equilibria.By a simple calculation, we conclude that ifR0>1,(29)TR0−1>ThR0lnR0,a>R0lnR0ETTR0−1−ThR0lnR0,then there exists a unique regular equilibrium for subsystem SG1, that is, ER1=H1∗,P1∗. If R0>1,(30)TR0−1>ThR0lnR0,a<R0lnR0ETTR0−1−ThR0lnR0,then ER1 becomes a virtual equilibrium denoted by EV1.Similarly, if(31)R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,a<1−εpR1ln1−pR1−τaTETT1−pR1−T−Th1−pR1ln1−pR1,then there exists a unique regular equilibrium for subsystem SG2, that is, ER2=H2∗,P2∗. Moreover, if(32)R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,a<1−εpR1ln1−pR1−τaTETT1−pR1−T−Th1−pR1ln1−pR1,then ER2 becomes a virtual equilibrium denoted by EV2.Further, ifR0>1,(33)TR0−1>ThR0lnR0,R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,then the two regular equilibria ER1 and ER2 can coexist only when (q≐1−p):(34)R0lnR0ETTR0−1−ThR0lnR0<a<1−εpR1lnqR1−τaTETTqR1−T−ThqR1lnqR1.Finally, ifR0>1,(35)TR0−1>ThR0lnR0,R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,then the two regular equilibria EV1 and EV2 can coexist only when(36)1−εpR1lnqR1−τaTETTqR1−T−ThqR1lnqR1<a<R0lnR0ETTR0−1−ThR0lnR0.The interesting question is how the types of the interior equilibria of model (4) change as the parameters vary. Theoretically, it is impossible to answer this question due to the positions related to the switching curve of all interior equilibria of the two subsystems which cannot be analytically determined. Therefore, in the following, we reveal the existence of regular and virtual equilibria and their coexistence by two-parameter bifurcation analyses.We first choosea and ET as bifurcation parameters and fix the other parameters as given in Figure 1. Letting a increase from 0.001 to 0.1 and with ET increasing from 0.01 to 1, the results show that the parameter space is divided into eight regions. The results shown in Figure 1 reveal that the variations of two parameters a and ET can significantly affect the types of equilibria of system (4), which are crucial for pest control and host outbreaks. For example, when the instantaneous search rate a∈0.001,0.0038, there is no interior equilibrium in the two subsystems SG1 and SG2, as shown in region I. If the instantaneous search rate a in 0.0038,0.011, the parameter space can be divided into two regions: region II-1 and region II-2. In region II-1, only EV2 exists, while in region II-2, only ER2 exists. In particular, if we let a=0.008, then according to region II-2 of Figure 1 with ET=0.6, the unique interior equilibrium of subsystem SG2 is regular (ER2 here), but no interior equilibrium of subsystem SG1 exists. If ET further increases and goes beyond a certain threshold value ET=0.76, the regular equilibrium ER2 becomes a virtual equilibrium EV2 as shown in region II-1. If the instantaneous search rate a is in 0.011,0.069, the parameter space is divided into three regions which are regions III-1, III-2, and III-3. In region III-1, ER1 and EV2 coexist; in region III-2, EV1 and EV2 coexist; in region III-3, EV1 and ER2 coexist. When the instantaneous search rate a increases and exceeds a certain threshold value a=0.069, that is, a∈0.069,0.1, then the parameter space is divided into two regions indicated as IV-1 and IV-2. In this case, the interior equilibrium of subsystem SG2 disappears, and only the interior equilibrium of subsystem SG1 exists. In region IV-1, only ER1 exists, while in region IV-2, only EV1 exists.Figure 1
Bifurcation diagrams for the existence and coexistence of regular and virtual equilibria of system (4) with respect to parameters a and ET. The different regions refer to different types of equilibria which are clarified in the text. The other parameters are fixed as follows: r=3.5,p=0.1,T=100,Th=1,τ=0.4,K=1,andε=0.3.The stability and types of equilibria of free system are important for host control. For example, it is finally free of control if the equilibrium of the free system is globally stable. Therefore, in order to design the best control policies to prevent host outbreaks, one of the possible ways is to choose a desirable switching curve such that all equilibria of subsystemSG1 become regular and all equilibria of subsystem SG2 become virtual. Thus, the parameters a and ET should be carefully chosen such that the interior equilibria of the two subsystems SG1 and SG2 are in region III-1. Similarly, we can address the effects of all other parameters on the regular and virtual equilibria of system (4).


### 4.2. Bifurcation Analysis about Key Parameters

In order to investigate the complex dynamic behaviour of system (4), we illustrate the one-dimensional parameter diagrams of the following parameters: the instantaneous search rate a, the intrinsic growth rate r, and the killing rate p. Moreover, we also give the two-dimensional parameter bifurcation plane of r×a. The purpose of these numerical analyses is to reveal the complex dynamics and the types of attractors and their changes as the parameters change.We first choose the instantaneous search ratea as the bifurcation parameter and fix the other parameters as shown in Figure 2. It follows from Figure 2 that nonlinear switching system (4) has more complex and interesting dynamic behaviour as the instantaneous search rate increases. When a increases from 0.03 to 0.06, we can find periodic, quasiperiodic, and chaotic oscillations and a Neimark–Sacker bifurcation [28–30]. As the parameter a increases from 0.03 to 0.033, system (4) has a stable solution (see Figure 3(a) for details when a=0.032). When a is further increased to 0.03305, an attracting invariant curve generates a Neimark–Sacker bifurcation, which is the discrete-time analogue of a Hopf bifurcation (see Figures 3(b) and 3(c) for details). Moreover, the quasiperiodic attractor appears suddenly at a=0.0345 and a=0.0455, at which the quasiperiodic attractors in the phase plane give one and eight closed curves (see Figures 3(d) and 3(e)). However, a chaotic solution emerges abruptly at a=0.055 (see Figure 3(f) for details).Figure 2
Bifurcation diagrams of system (4) with respect to parameter a. The other parameters are fixed as follows: r=3.2,p=0.5,T=100,Th=1,τ=0.1,K=1,ET=0.8,andε=0.5.
(a)(b)Figure 3
Phase trajectories of system (4) with different a and initial values. (a) a=0.032. (b) a=0.03305. (c) a=0.0335. (d) a=0.0345. (e) a=0.0455. (f) a=0.055. The initial values H0,P0 in (a)–(c) are 0.6,0.4 and in (d)–(f) are 0.7,0.31. The other parameters are fixed as follows: r=3.2,p=0.5,T=100,Th=1,τ=0.1,K=1,ET=0.8,andε=0.5.
(a)(b)(c)(d)(e)(f)In addition, when we choose the intrinsic growth rater and the killing rate p as the bifurcation parameters, system (4) also shows complex dynamic behaviour as shown in Figures 4–6. It follows from the bifurcation diagram with respect to parameter r (see Figure 4) that there is period-doubling bifurcation for a range of parameter values, that is, a solution with period-16 to period-32 by the period-doubling bifurcation and then the solution with period-32 to period-16 through the period-halving phenomena when r∈2.99,3.14. Moreover, we can also observe the coexistence of multiple attractors, for example, r∈2.63,2.66. The details will be discussed in the following section.Figure 4
Bifurcation diagrams of system (4) with respect to parameter r. The other parameters are fixed as follows: a=0.041,T=100,Th=1,p=0.5,τ=0.1,K=1,ET=0.8,andε=0.6.
(a)(b)Figure 5
Periods of period point cycles of system (4) with respect to parameters r and a, where different color intensities represent different periodic values. The other parameters are fixed as follows: p=0.4,T=100,Th=1,τ=0.6,K=1,ET=0.8,andε=0.6.Figure 6
Bifurcation diagrams of system (4) with respect to parameter p. The other parameters are fixed as follows: r=1.9,a=0.033,T=100,Th=1,τ=0.2,K=1,ET=0.8,andε=0.6.
(a)(b)To investigate the effect of different parameter values on the periodicity of solutions for both host and parasitoid populations, according to published methods [31–34], a two-dimensional bifurcation diagram is developed in Figure 5 using the intrinsic growth rate r and the instantaneous search rate a as bifurcation parameters, in which the periods of period point cycles have been marked with different colors. With parameter increases, system (4) experiences period-adding bifurcation and there are many periodic windows embedded in the chaotic areas. The periodic solutions with different periods are plotted in different colors and marked by the corresponding numbers, in particular, the number 2 represents a period 2-point cycle, the number 14 represents period 14-point cycle, and the chaotic region has been marked when the number is equal or greater than 30. The variations of complex parameter spaces related to the periods of solutions of system (4) reveal that small changes of key parameters could significantly influence the oscillation patterns of both host and parasitoid populations. Consequently, this could result in difficulties for successful pest control.The bifurcation diagram with respect to parameterp shown in Figure 6 reveals a phenomenon such that if the killing rate of the pesticides is too low or too large (e.g., p∈0.01,0.22 and 0.81,0.98), the host population will have outbreaks. From the ecological point of view, if we choose a pesticide with low efficacy, i.e., p∈0.01,0.22, it may not be able to kill pests and eventually lead to outbreaks; if we choose a pesticide with high efficacy, i.e., p∈0.81,0.98, it may kill natural enemies while killing pests, but in this case, the pest population can still have an outbreak. All these results reveal that when using high-dose insecticides to kill hosts in large quantities, the parasitoid population will die in large numbers later due to the lack of hosts, which will eventually lead to the resurgence of the host population. Moreover, the main results shown in Figure 6 clearly reveal how to choose the appropriate killing rate to make the whole system stabilize at the attractor of subsystem SG1. In particular, if we choose the killing rate p∈0.23,0.80, then the pest population can stabilize in subsystem SG1 so as to prevent the outbreak of pests.


## 4.1. Equilibria Bifurcation for Nonlinear Switching System (4)

According to Definition1, nonlinear switching system (4) has different types of equilibria, and these equilibria play key roles in disease and pest control. For example, Xiao et al. [27] investigated the control of emerging infectious disease and Tang et al. [16] addressed pest control with an ET. In the following, we focus on how the system with a nonlinear switching threshold affects the distribution of regular and virtual equilibria.By a simple calculation, we conclude that ifR0>1,(29)TR0−1>ThR0lnR0,a>R0lnR0ETTR0−1−ThR0lnR0,then there exists a unique regular equilibrium for subsystem SG1, that is, ER1=H1∗,P1∗. If R0>1,(30)TR0−1>ThR0lnR0,a<R0lnR0ETTR0−1−ThR0lnR0,then ER1 becomes a virtual equilibrium denoted by EV1.Similarly, if(31)R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,a<1−εpR1ln1−pR1−τaTETT1−pR1−T−Th1−pR1ln1−pR1,then there exists a unique regular equilibrium for subsystem SG2, that is, ER2=H2∗,P2∗. Moreover, if(32)R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,a<1−εpR1ln1−pR1−τaTETT1−pR1−T−Th1−pR1ln1−pR1,then ER2 becomes a virtual equilibrium denoted by EV2.Further, ifR0>1,(33)TR0−1>ThR0lnR0,R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,then the two regular equilibria ER1 and ER2 can coexist only when (q≐1−p):(34)R0lnR0ETTR0−1−ThR0lnR0<a<1−εpR1lnqR1−τaTETTqR1−T−ThqR1lnqR1.Finally, ifR0>1,(35)TR0−1>ThR0lnR0,R1>expτaT1−p,T1−pR1−1>Th1−pR1ln1−pR1,then the two regular equilibria EV1 and EV2 can coexist only when(36)1−εpR1lnqR1−τaTETTqR1−T−ThqR1lnqR1<a<R0lnR0ETTR0−1−ThR0lnR0.The interesting question is how the types of the interior equilibria of model (4) change as the parameters vary. Theoretically, it is impossible to answer this question due to the positions related to the switching curve of all interior equilibria of the two subsystems which cannot be analytically determined. Therefore, in the following, we reveal the existence of regular and virtual equilibria and their coexistence by two-parameter bifurcation analyses.We first choosea and ET as bifurcation parameters and fix the other parameters as given in Figure 1. Letting a increase from 0.001 to 0.1 and with ET increasing from 0.01 to 1, the results show that the parameter space is divided into eight regions. The results shown in Figure 1 reveal that the variations of two parameters a and ET can significantly affect the types of equilibria of system (4), which are crucial for pest control and host outbreaks. For example, when the instantaneous search rate a∈0.001,0.0038, there is no interior equilibrium in the two subsystems SG1 and SG2, as shown in region I. If the instantaneous search rate a in 0.0038,0.011, the parameter space can be divided into two regions: region II-1 and region II-2. In region II-1, only EV2 exists, while in region II-2, only ER2 exists. In particular, if we let a=0.008, then according to region II-2 of Figure 1 with ET=0.6, the unique interior equilibrium of subsystem SG2 is regular (ER2 here), but no interior equilibrium of subsystem SG1 exists. If ET further increases and goes beyond a certain threshold value ET=0.76, the regular equilibrium ER2 becomes a virtual equilibrium EV2 as shown in region II-1. If the instantaneous search rate a is in 0.011,0.069, the parameter space is divided into three regions which are regions III-1, III-2, and III-3. In region III-1, ER1 and EV2 coexist; in region III-2, EV1 and EV2 coexist; in region III-3, EV1 and ER2 coexist. When the instantaneous search rate a increases and exceeds a certain threshold value a=0.069, that is, a∈0.069,0.1, then the parameter space is divided into two regions indicated as IV-1 and IV-2. In this case, the interior equilibrium of subsystem SG2 disappears, and only the interior equilibrium of subsystem SG1 exists. In region IV-1, only ER1 exists, while in region IV-2, only EV1 exists.Figure 1
Bifurcation diagrams for the existence and coexistence of regular and virtual equilibria of system (4) with respect to parameters a and ET. The different regions refer to different types of equilibria which are clarified in the text. The other parameters are fixed as follows: r=3.5,p=0.1,T=100,Th=1,τ=0.4,K=1,andε=0.3.The stability and types of equilibria of free system are important for host control. For example, it is finally free of control if the equilibrium of the free system is globally stable. Therefore, in order to design the best control policies to prevent host outbreaks, one of the possible ways is to choose a desirable switching curve such that all equilibria of subsystemSG1 become regular and all equilibria of subsystem SG2 become virtual. Thus, the parameters a and ET should be carefully chosen such that the interior equilibria of the two subsystems SG1 and SG2 are in region III-1. Similarly, we can address the effects of all other parameters on the regular and virtual equilibria of system (4).


## 4.2. Bifurcation Analysis about Key Parameters

In order to investigate the complex dynamic behaviour of system (4), we illustrate the one-dimensional parameter diagrams of the following parameters: the instantaneous search rate a, the intrinsic growth rate r, and the killing rate p. Moreover, we also give the two-dimensional parameter bifurcation plane of r×a. The purpose of these numerical analyses is to reveal the complex dynamics and the types of attractors and their changes as the parameters change.We first choose the instantaneous search ratea as the bifurcation parameter and fix the other parameters as shown in Figure 2. It follows from Figure 2 that nonlinear switching system (4) has more complex and interesting dynamic behaviour as the instantaneous search rate increases. When a increases from 0.03 to 0.06, we can find periodic, quasiperiodic, and chaotic oscillations and a Neimark–Sacker bifurcation [28–30]. As the parameter a increases from 0.03 to 0.033, system (4) has a stable solution (see Figure 3(a) for details when a=0.032). When a is further increased to 0.03305, an attracting invariant curve generates a Neimark–Sacker bifurcation, which is the discrete-time analogue of a Hopf bifurcation (see Figures 3(b) and 3(c) for details). Moreover, the quasiperiodic attractor appears suddenly at a=0.0345 and a=0.0455, at which the quasiperiodic attractors in the phase plane give one and eight closed curves (see Figures 3(d) and 3(e)). However, a chaotic solution emerges abruptly at a=0.055 (see Figure 3(f) for details).Figure 2
Bifurcation diagrams of system (4) with respect to parameter a. The other parameters are fixed as follows: r=3.2,p=0.5,T=100,Th=1,τ=0.1,K=1,ET=0.8,andε=0.5.
(a)(b)Figure 3
Phase trajectories of system (4) with different a and initial values. (a) a=0.032. (b) a=0.03305. (c) a=0.0335. (d) a=0.0345. (e) a=0.0455. (f) a=0.055. The initial values H0,P0 in (a)–(c) are 0.6,0.4 and in (d)–(f) are 0.7,0.31. The other parameters are fixed as follows: r=3.2,p=0.5,T=100,Th=1,τ=0.1,K=1,ET=0.8,andε=0.5.
(a)(b)(c)(d)(e)(f)In addition, when we choose the intrinsic growth rater and the killing rate p as the bifurcation parameters, system (4) also shows complex dynamic behaviour as shown in Figures 4–6. It follows from the bifurcation diagram with respect to parameter r (see Figure 4) that there is period-doubling bifurcation for a range of parameter values, that is, a solution with period-16 to period-32 by the period-doubling bifurcation and then the solution with period-32 to period-16 through the period-halving phenomena when r∈2.99,3.14. Moreover, we can also observe the coexistence of multiple attractors, for example, r∈2.63,2.66. The details will be discussed in the following section.Figure 4
Bifurcation diagrams of system (4) with respect to parameter r. The other parameters are fixed as follows: a=0.041,T=100,Th=1,p=0.5,τ=0.1,K=1,ET=0.8,andε=0.6.
(a)(b)Figure 5
Periods of period point cycles of system (4) with respect to parameters r and a, where different color intensities represent different periodic values. The other parameters are fixed as follows: p=0.4,T=100,Th=1,τ=0.6,K=1,ET=0.8,andε=0.6.Figure 6
Bifurcation diagrams of system (4) with respect to parameter p. The other parameters are fixed as follows: r=1.9,a=0.033,T=100,Th=1,τ=0.2,K=1,ET=0.8,andε=0.6.
(a)(b)To investigate the effect of different parameter values on the periodicity of solutions for both host and parasitoid populations, according to published methods [31–34], a two-dimensional bifurcation diagram is developed in Figure 5 using the intrinsic growth rate r and the instantaneous search rate a as bifurcation parameters, in which the periods of period point cycles have been marked with different colors. With parameter increases, system (4) experiences period-adding bifurcation and there are many periodic windows embedded in the chaotic areas. The periodic solutions with different periods are plotted in different colors and marked by the corresponding numbers, in particular, the number 2 represents a period 2-point cycle, the number 14 represents period 14-point cycle, and the chaotic region has been marked when the number is equal or greater than 30. The variations of complex parameter spaces related to the periods of solutions of system (4) reveal that small changes of key parameters could significantly influence the oscillation patterns of both host and parasitoid populations. Consequently, this could result in difficulties for successful pest control.The bifurcation diagram with respect to parameterp shown in Figure 6 reveals a phenomenon such that if the killing rate of the pesticides is too low or too large (e.g., p∈0.01,0.22 and 0.81,0.98), the host population will have outbreaks. From the ecological point of view, if we choose a pesticide with low efficacy, i.e., p∈0.01,0.22, it may not be able to kill pests and eventually lead to outbreaks; if we choose a pesticide with high efficacy, i.e., p∈0.81,0.98, it may kill natural enemies while killing pests, but in this case, the pest population can still have an outbreak. All these results reveal that when using high-dose insecticides to kill hosts in large quantities, the parasitoid population will die in large numbers later due to the lack of hosts, which will eventually lead to the resurgence of the host population. Moreover, the main results shown in Figure 6 clearly reveal how to choose the appropriate killing rate to make the whole system stabilize at the attractor of subsystem SG1. In particular, if we choose the killing rate p∈0.23,0.80, then the pest population can stabilize in subsystem SG1 so as to prevent the outbreak of pests.


## 5. Initial Sensitivities and Coexistence of Multiple Attractors

It is well known that the initial densities of host and parasitoid populations will have a significant impact on the dynamic behaviour of the system, so they are key to biological control, and in this section, we will investigate how they affect outbreak patterns or final states of the host population.


### 5.1. Initial Sensitivities

In order to analyse the relationship between initial densities of populations and host outbreak patterns, Figure7 illustrates the effect of initial densities on the outbreak frequencies related to the nonlinear switching curve. In Figure 7(a), the initial densities of the host and parasitoid populations are 0.4,0.5, respectively, and the numerical results indicate that the density of the host population never reaches the switching surface, which shows that the solutions initiating from (0.4, 0.5) are free from IPM control measures. If we choose the initial densities as 0.4,0.65 and 0.5,0.56, the simulation results indicate that the host population will not have an outbreak after applying one or two IPM strategies (see Figures 7(b) and 7(c) for details). However, if we set the initial densities to be 0.5,0.6, in order to control the weighted density of the host population below ET, the IPM measures should be implemented multiple times, as shown in Figure 7(d).Figure 7
Illustrating the switching effects of initial densities of the host and parasitoid populations of system (4) related to the switching curve. The parameters are fixed as follows: r=1.9,a=0.033,T=100,Th=1,τ=0.2,K=1,ET=0.8,andε=0.6. The initial values H0,P0 in (a)–(d) are 0.4,0.5, 0.4,0.65, 0.5,0.56, and 0.5,0.6. The red lines in each subplot represent the switching curves.
(a)(b)(c)(d)In order to reveal the effects of initial values of both host and parasitoid populations on the outbreak patterns of the host population in more detail, the basin of attraction of the host outbreak frequencies with respect to initial densities is shown in Figure8, and the regions I, II, III, IV, and V with different colors clearly clarify the initial sensitivities. An outbreak of the weighted host population could not occur in region I, where the system is always stabilized at the attractor of subsystem SG1. In regions II, III, and IV, the weighted density of two generations of the host density can be lower than the given threshold ET after implementing the IPM strategy one, two, and three times, while the weighted host population will have outbreaks several times if the initial densities lie in region V before the system stabilizes in subsystem SG1. All these results confirm that the initial densities of both populations are crucial for the outbreak patterns and frequencies of host population.Figure 8
Dependence of the host outbreak frequencies on the initial valuesH0,P0 of system (4). The parameters are fixed as follows: r=1.9,p=0.7,a=0.033,T=100,Th=1,τ=0.2,K=1,ET=0.8,andε=0.6.


### 5.2. Coexistence of Multiple Attractors

The initial values not only influence the host outbreak frequencies but can also affect the host final stable states, i.e., multiple attractors could coexist. It follows from the bifurcation analyses that multiple attractors can coexist for a wide range of parameters (see Figure4 as an example). To further confirm this and address their biological implications, we set all parameters in Figure 9 and choose different initial densities to depict some typical trajectories. Especially, three host outbreak attractors coexist when r=2.288, as shown in Figure 9, which show different periods. If we set the initial densities to be 0.4,0.3, then the solution of system (4) approaches the periodic attractor with period 9 as shown in Figure 9(a). If we choose the initial values as 0.4,0.1, then the solution of system (4) is a periodic attractor of period 14 (see Figure 9(b)). Moreover, when the initial value is 0.5,0.3, the third attractor becomes quite complex (Figure 9(c)).Figure 9
Three coexisting attractors of system (4) with different initial values. The parameters are fixed as follows: r=2.288,a=0.05,T=100,Th=1,p=0.5,τ=0.3,K=1,ET=0.75,andε=0.6. (a) H0,P0=0.4,0.3. (b) H0,P0=0.4,0.1. (c) H0,P0=0.5,0.3.
(a)(b)(c)In order to illustrate the initial sensitivities more specifically, the basin of attraction with respect to three different host outbreak solutions of coexistence is shown in Figure10. A basin of attraction exists for three attractors: the blue, magenta, and green regions. These are related to the periodic solutions which are shown in Figures 9(a)–9(c). We can also see that the curve εHt+1−εHt+1=ET divides the attraction regions into two distinct areas, each of which has a significant difference in dynamic behaviour. Moreover, the sensitivity and dependence of initial values are very obvious, and the final attractors of host and parasitoid populations depend on their initial densities. Note that all these results are confirmed by basins of attraction of initial densities.Figure 10
Basins of attraction of three coexisting attractors of system (4). The range of the left one is 0.01≤H0≤1,0.01≤P0≤1. The blue, magenta, and green regions are attracted to the attractors shown in Figure 9 from left to right. The right one is an enlargement of the basins of attraction with range 0.3≤H0≤1,0.01≤P0≤0.4. The parameters are fixed as follows: r=2.288,a=0.05,T=100,Th=1,p=0.5,τ=0.3,K=1,ET=0.75,andε=0.6. The white lines in each subplot represent the switching curves.
(a)(b)


## 5.1. Initial Sensitivities

In order to analyse the relationship between initial densities of populations and host outbreak patterns, Figure7 illustrates the effect of initial densities on the outbreak frequencies related to the nonlinear switching curve. In Figure 7(a), the initial densities of the host and parasitoid populations are 0.4,0.5, respectively, and the numerical results indicate that the density of the host population never reaches the switching surface, which shows that the solutions initiating from (0.4, 0.5) are free from IPM control measures. If we choose the initial densities as 0.4,0.65 and 0.5,0.56, the simulation results indicate that the host population will not have an outbreak after applying one or two IPM strategies (see Figures 7(b) and 7(c) for details). However, if we set the initial densities to be 0.5,0.6, in order to control the weighted density of the host population below ET, the IPM measures should be implemented multiple times, as shown in Figure 7(d).Figure 7
Illustrating the switching effects of initial densities of the host and parasitoid populations of system (4) related to the switching curve. The parameters are fixed as follows: r=1.9,a=0.033,T=100,Th=1,τ=0.2,K=1,ET=0.8,andε=0.6. The initial values H0,P0 in (a)–(d) are 0.4,0.5, 0.4,0.65, 0.5,0.56, and 0.5,0.6. The red lines in each subplot represent the switching curves.
(a)(b)(c)(d)In order to reveal the effects of initial values of both host and parasitoid populations on the outbreak patterns of the host population in more detail, the basin of attraction of the host outbreak frequencies with respect to initial densities is shown in Figure8, and the regions I, II, III, IV, and V with different colors clearly clarify the initial sensitivities. An outbreak of the weighted host population could not occur in region I, where the system is always stabilized at the attractor of subsystem SG1. In regions II, III, and IV, the weighted density of two generations of the host density can be lower than the given threshold ET after implementing the IPM strategy one, two, and three times, while the weighted host population will have outbreaks several times if the initial densities lie in region V before the system stabilizes in subsystem SG1. All these results confirm that the initial densities of both populations are crucial for the outbreak patterns and frequencies of host population.Figure 8
Dependence of the host outbreak frequencies on the initial valuesH0,P0 of system (4). The parameters are fixed as follows: r=1.9,p=0.7,a=0.033,T=100,Th=1,τ=0.2,K=1,ET=0.8,andε=0.6.


## 5.2. Coexistence of Multiple Attractors

The initial values not only influence the host outbreak frequencies but can also affect the host final stable states, i.e., multiple attractors could coexist. It follows from the bifurcation analyses that multiple attractors can coexist for a wide range of parameters (see Figure4 as an example). To further confirm this and address their biological implications, we set all parameters in Figure 9 and choose different initial densities to depict some typical trajectories. Especially, three host outbreak attractors coexist when r=2.288, as shown in Figure 9, which show different periods. If we set the initial densities to be 0.4,0.3, then the solution of system (4) approaches the periodic attractor with period 9 as shown in Figure 9(a). If we choose the initial values as 0.4,0.1, then the solution of system (4) is a periodic attractor of period 14 (see Figure 9(b)). Moreover, when the initial value is 0.5,0.3, the third attractor becomes quite complex (Figure 9(c)).Figure 9
Three coexisting attractors of system (4) with different initial values. The parameters are fixed as follows: r=2.288,a=0.05,T=100,Th=1,p=0.5,τ=0.3,K=1,ET=0.75,andε=0.6. (a) H0,P0=0.4,0.3. (b) H0,P0=0.4,0.1. (c) H0,P0=0.5,0.3.
(a)(b)(c)In order to illustrate the initial sensitivities more specifically, the basin of attraction with respect to three different host outbreak solutions of coexistence is shown in Figure10. A basin of attraction exists for three attractors: the blue, magenta, and green regions. These are related to the periodic solutions which are shown in Figures 9(a)–9(c). We can also see that the curve εHt+1−εHt+1=ET divides the attraction regions into two distinct areas, each of which has a significant difference in dynamic behaviour. Moreover, the sensitivity and dependence of initial values are very obvious, and the final attractors of host and parasitoid populations depend on their initial densities. Note that all these results are confirmed by basins of attraction of initial densities.Figure 10
Basins of attraction of three coexisting attractors of system (4). The range of the left one is 0.01≤H0≤1,0.01≤P0≤1. The blue, magenta, and green regions are attracted to the attractors shown in Figure 9 from left to right. The right one is an enlargement of the basins of attraction with range 0.3≤H0≤1,0.01≤P0≤0.4. The parameters are fixed as follows: r=2.288,a=0.05,T=100,Th=1,p=0.5,τ=0.3,K=1,ET=0.75,andε=0.6. The white lines in each subplot represent the switching curves.
(a)(b)


## 6. Discussion

In this paper, we discussed a nonlinear switching discrete host-parasitoid model with Holling type II functional response function, where the switching strategies and control measures are guided by the weighted density of two generations of the pest population. The important innovation of the model lies in the threshold control strategy. Here, we consider the weighted density of two generations of the pest population as the threshold to judge whether the control strategy is implemented or not, which not only brings new challenges in theory but also makes pest control more operable. The most interesting theoretical analyses are to determine the types of all possible equilibria of the whole switching system and their stabilities because if the real equilibrium of the free subsystem is stable for the whole switching system, then the purpose of pest control can be realized easily. Further, the sensitivity analyses related to the key parameters and initial densities of both host and parasitoid populations can help us to design suitable control measures.We first investigated the existence and local stability of equilibria of two subsystems and briefly analysed the existence and coexistence of the equilibria of the whole nonlinear switching system. We used thea×ET two-parameter bifurcation diagram to reveal the distribution regions of different types of equilibria including regular and virtual equilibria, as shown in Figure 1. From the perspective of pest control, the parameters a and ET should be carefully chosen such that the interior equilibria of the two subsystems SG1 and SG2 are in region III-1. Moreover, we provided the r×a two-parameter bifurcation diagram that reveals the existence of period-adding bifurcation in this nonlinear switching system, as shown in Figure 5. This suggests that we should choose appropriate parameters to avoid the solution of system (4) falling into the chaotic region for successful implementation of an IPM strategy.On the other hand, the one-parameter bifurcation diagrams which were derived from system (4) reveal the complex dynamical behaviour. For example, the bifurcation diagram of the instantaneous search rate a shows that system (4) may have very complex dynamical behaviour such as periodic, quasiperiodic, and chaotic solutions and a Neimark–Sacker bifurcation (Figures 2 and 3). The bifurcation diagram of the intrinsic growth rate r illustrates that system (4) may have multiple coexisting attractors, period-doubling bifurcation, and period-halving bifurcation (Figure 4). Moreover, it follows from the bifurcation diagram of the killing rate p that the host population can stabilize in subsystem SG1 (Figure 6), which emphasises the importance of choosing the appropriate concentration of pesticide.In addition, the relationship between initial densities and pest control was also studied. The results show that the initial densities of the host and parasitoid populations will affect the outcome of an IPM strategy, and the final stable states of the populations depend on their initial densities (Figures7–10). In those figures, the basins of initial attractors and outbreak frequencies of the host population have been revealed, which clarify that the outbreak patterns of the host population could become markedly different as parameters vary and change.Compared with the basic model and the main published results [21], we conclude that it is more flexible for long-term monitoring and successful control of the host population if we employ the weighted density of two generations of its population as the standard for deciding whether the control strategy is implemented or not. In particular, in a certain parameter space, the nonlinear switching curve can ensure that the stable region of the real equilibrium of the free subsystem increases, thus increasing the controllability of the host population. Moreover, the one-parameter bifurcation and two-parameter bifurcation analyses reveal that the dynamic behaviour of system (4) is more variable, including the appearance of Neimark–Sacker bifurcations and period-adding bifurcations. Obviously, the nonlinear switching curve divides the basin of attraction of outbreak frequencies and coexistence attractors of the host and parasite populations into two distinct patterns, as shown in Figures 8 and 10, respectively.



---

*Source: 9425285-2020-06-09.xml*