# Codimension-2 Bifurcation Analysis and Control of a Discrete Mosquito Model with a Proportional Release Rate of Sterile Mosquitoes

Authors: Qiaoling Chen; Zhidong Teng; Junli Liu; Feng Wang
Journal: Complexity
 (2020)

---

## Abstract

This paper concerns a discrete wild and sterile mosquito model with a proportional release rate of sterile mosquitoes. It is shown that the discrete model undergoes codimension-2 bifurcations with 1 : 2, 1 : 3, and 1 : 4 strong resonances by applying the bifurcation theory. Some numerical simulations, including codimension-2 bifurcation diagrams, maximum Lyapunov exponents diagrams, and phase portraits, are also presented to illustrate the validity of theoretical results and display the complex dynamical behaviors. Moreover, two control strategies are applied to the model.


---

## Body

## 1. Introduction

In recent years, the mosquito-borne diseases, including West Nile virus, dengue, malaria, Japanese encephalitis, and chikungunya, have attracted much attention. Until now, there are no available vaccines for many of these viruses, which constitute a significant threat to global health. To control the transmissions of mosquito-borne diseases, it is necessary to adopt some strategies to reduce the number of wild mosquitoes. One of effective and environment-friendly approaches is sterile insect technique (SIT). Generally, the mass male mosquitoes are exposed to radiation and thereby rendered sexually sterile. After being released into the desired area, the treated male mosquitoes will mate with wide female mosquitoes and their offspring is nonviable. In this process, choosing the number of sterile mosquitoes and suitable releasing strategy are particularly important.To understand the interactive dynamics of wide mosquitoes and sterile mosquitoes, some mathematical models have been introduced, for example, [1–9]. Recently, Cai et al. [7] proposed three mathematical models which represent different releasing strategies of sterile mosquitoes. The existence and stability of equilibria were investigated, and the threshold dynamics were obtained. One of the models considered in [7] is(1)dwdt=aw1+w+g−μ1+ξ1w+gw,dgdt=bw−μ2+ξ2w+gg,where wandg are the numbers of wild mosquitoes and sterile mosquitoes, respectively. All the coefficients are positive constants, a is the number of offspring produced by per mate, bw is the linear release rate of sterile mosquitoes with coefficient b, and μi and ξii=1,2 represent the density independent and dependent death rates of wild and sterile mosquitoes, respectively. It was shown in [10, 11] that continuous model (1) undergoes a sequence of codimensional-1 bifurcations (including saddle-node bifurcation, supercritical and subcritical Hopf bifurcation, and homoclinic bifurcation) and codimensional-2 bifurcation (Bogdanov–Takens bifurcation). When bw in (1) is replaced by bw/1+w, the system becomes a mosquito population model with nonlinear saturated releasing rate of sterile mosquitoes. Some complex bifurcation phenomena, such as Hopf bifurcation of codimension-3, nilpotent cusp of codimension-4, and Bogdanov–Takens bifurcation of codimension-2, were provided in [12] if μ1≠μ2 and ξ1≠ξ2.In mathematical modelling, difference equations (that is, discrete-time models) are also widely used to describe the population dynamics. Compared with continuous models, discrete-time cases have many advantages. The most obvious one is that discrete models are more convenient to formulate and simulate, especially when the populations have nonoverlapping generations [13–15]. Moreover, for several biological models, their discrete-time cases may produce much richer dynamical behaviors than the continuous cases.In this paper, we study the complex dynamics of the following discrete-time model:(2)wn+1=wnexpawn1+wn+gn−μ1+ξ1wn+gn,gn+1=gn+bwn−μ2+ξ2wn+gngn,which is established from (1) by the Euler forward difference method. w,g,a,b,μi,ξii=1,2 are defined as in model (1). By the biological meaning of the model, we only consider system (2) in the region Ω=w,g:w≥0,g≥0 in the w,g-plane. We will study the existence of codimension-2 bifurcations associated with 1 : 2, 1 : 3, and 1 : 4 strong resonances of model (2). The codimension-2 bifurcation diagrams, maximum Lyapunov exponent diagrams, and phase portraits are also presented to illustrate the validity of theoretical results and display the complex dynamical behaviors.The bifurcation and chaotic behavior of discrete models have been widely investigated by many researchers, but most of them focus on codimensional-1 bifurcations, for example, [16–21]. The codimension-2 bifurcations for discrete models are considered very few. The codimension-2 bifurcations of a delayed discrete-time Hopfield neural network, discrete epidemic model with nonlinear incidence rate, and discrete SIS epidemic model with standard incidence are analyzed in [22–24]. A discrete Hindmarsh–Rose model undergoes codimension-2 bifurcations with 1 : 2 and 1 : 4 strong resonances [25], and 1 : 3 strong resonance [26]. 1 : 1 strong resonance bifurcation of a discrete predator-prey model with nonmonotonic functional response has also been studied in [27]. Moreover, Ren and Yu [28] investigated codimension-2 bifurcations of a discrete information diffusion model.The main goal of this paper is to investigate the impact of the mating ratea and the release rate coefficient b on the interactive dynamics of wide mosquitoes and sterile mosquitoes. The rest of this paper is organized as follows. Preliminaries are given in Section 2. In Sections 3–5, some sufficient conditions of codimension-2 bifurcations associated with 1 : 2, 1 : 3, and 1 : 4 strong resonances are established, respectively. The numerical simulations and controller designs are presented in Sections 6 and 7. A brief discussion is given in Section 8.


## 2. Preliminaries

In this section, we first briefly introduce some known results on the existence of positive equilibria, which have been obtained in [7], and then provide some sufficient conditions for the existence of codimension-2 bifurcations with 1 : 2, 1 : 3, and 1 : 4 strong resonances.Define the function(3)FN=1+Nμ1+ξ1Nb+μ2+ξ2N−aNμ2+ξ2N=ξ1ξ2N3+ξ1b+μ2+ξ2μ1+ξ1−aN2+μ1ξ2+bμ1+ξ1+μ1+ξ1−aμ2N+μ1b+μ2.We can check that model (2) admits a positive equilibrium:(4)w,g=1+N¯μ1+ξ1N¯a,b1+N¯μ1+ξ1N¯aμ2+ξ2N¯,if and only if FN=0 has a positive root N¯, which is also equivalent to the existence of positive solution N¯ to the equation GN=b with(5)GN≔aN−1+Nμ1+ξ1Nμ2+ξ2N1+Nμ1+ξ1N.Assuming thata>μ1+ξ1+2μ1ξ1, we know the equation GN=0 has two positive roots:(6)N1,2=12ξ1a−μ1−ξ1∓a−μ1−ξ12−4μ1ξ1.It follows that there exists a uniqueN∗∈N1,N2 such that G′N∗=0. Let(7)b0≔GN∗,then the equation GN=b has two positive solutions N1,2∗∈N1,N2 with N1∗<N2∗ if and only if b<b0, and only one positive solution N∗∈N1,N2 when b=b0. In summary, we have the following results on the existence of equilibria of model (2).Lemma 1.
Assume thata>μ1+ξ1+2μ1ξ1. Then,(i)
model (2) always has a boundary equilibrium E0,0(ii)
ifb>b0, then model (2) has no any positive equilibrium(iii)
ifb=b0, then model (2) has a unique positive equilibrium E∗w∗,g∗, where(8)w∗,g∗=1+N∗μ1+ξ1N∗a,b01+N∗μ1+ξ1N∗aμ2+ξ2N∗.Moreover,FN∗=0 and F′N∗=0.(iv)
ifb<b0, then model (2) has two positive equilibria E1∗w1∗,g1∗ and E2∗w2∗,g2∗, where(9)w1∗,g1∗=1+N1∗μ1+ξ1N1∗a,b1+N1∗μ1+ξ1N1∗aμ2+ξ2N1∗,w2∗,g2∗=1+N2∗μ1+ξ1N2∗a,b1+N2∗μ1+ξ1N2∗aμ2+ξ2N2∗.Moreover,F′N1∗<0 and F′N2∗>0.Denote byEw,g, the equilibrium of model (2), then the Jacobian matrix at Ew,g is(10)JE=1+wa1+g1+N2−ξ1−waw1+N2+ξ1b−ξ2g1−μ2−ξ2N+g,its trace is TrJE=2+wa1+g/1+N2−ξ1−μ2−ξ2N+g=2+1/aAN, and its determinant is DetJE=1−w/1+NF′N+1/aAN, where(11)AN≔ξ1ξ2−2ξ1N2+aξ1−2ξ2+ξ1ξ2−ξ1+μ1ξ2−3ξ1N+aμ1−μ2+μ1ξ2−ξ1−μ1.If the matrixJE has eigenvalues λ1,2=e±iθ with θ=π, θ=2/3π, or θ=1/2π, then there may present a codimension-2 bifurcation with 1 : 2, 1 : 3, or 1 : 4 strong resonances at Ew,g. For the equilibrium E∗w∗,g∗ in Lemma 1iii, if codimension-2 bifurcation with 1 : 2, 1 : 3, or 1 : 4 strong resonance occurs, then DetJE∗=1 and TrJE∗=−2,−1, or 0, which imply that 1/aAN∗<0 and w∗/1+N∗F′N∗=1/aAN∗<0. It is a contradiction with F′N∗=0. Similarly, since F′N2∗>0, we can also show that there does not exist the above codimension-2 bifurcations at the equilibrium E2∗w2∗,g2∗. Thus, in this paper, we discuss the existence of codimension-2 bifurcation at E1∗w1∗,g1∗.Choosinga and b as bifurcation parameters, if a,b,μ1,μ2,ξ1,ξ2∈F1, where(12)F1≔a,b,μ1,μ2,ξ1,ξ2∈ℝ+6:a>μ1+ξ1+2μ1ξ1,b<b0,w1∗1+N1∗F′N1∗=1aAN1∗=−4,then two eigenvalues of Jacobian matrix JE1∗ are λ1,2=e±πi=−1. This shows that there may present a codimension-2 bifurcation with 1 : 2 strong resonance at E1∗w1∗,g1∗.Ifa,b,μ1,μ2,ξ1,ξ2∈F2, where(13)F2≔a,b,μ1,μ2,ξ1,ξ2∈ℝ+6:a>μ1+ξ1+2μ1ξ1,b<b0,w1∗1+N1∗F′N1∗=1aAN1∗=−3,then two eigenvalues of the Jacobian matrix JE1∗ are λ1,2=e±2π/3i=−1/2±3/2i. This shows that there may present a codimension-2 bifurcation with 1 : 3 strong resonance at E1∗w1∗,g1∗.Furthermore, ifa,b,μ1,μ2,ξ1,ξ2∈F3, where(14)F3≔a,b,μ1,μ2,ξ1,ξ2∈ℝ+6:a>μ1+ξ1+2μ1ξ1,b<b0,w1∗1+N1∗F′N1∗=1aAN1∗=−2,then two eigenvalues of the Jacobian matrix JE1∗ are λ1,2=e±π/2i=±i. This shows that there may present a codimension-2 bifurcation with 1 : 4 strong resonance at E1∗w1∗,g1∗. In the following, we will give the theoretical analysis and numerical simulations.


## 3. 1 : 2 Strong Resonance

In this section, by applying the bifurcation theory in [29], we analyze the codimension-2 bifurcation associated with 1 : 2 strong resonance of model (2) at the equilibrium E1∗w1∗,g1∗.We choosea and b as bifurcation parameters and assume that a1,b1,μ1,μ2,ξ1,ξ2∈F1. Consider model (2) with a=a1+a1∗ and b=b1+b1∗, where a1∗ and b1∗ are small perturbations of the parameters a1 and b1. Let w^n=wn−w1∗ and g^n=gn−g1∗, and model (2) becomes into the following form which has the equilibrium 0,0:(15)w^n+1=w^n+w1∗expaw^n+w1∗1+w^n+g^n+w1∗+g1∗−μ1+ξ1w^n+g^n+w1∗+g1∗−w1∗,g^n+1=g^n+bw^n+w1∗−μ2+ξ2w^n+g^n+w1∗+g1∗g^n+g1∗.By expanding the right-hand side of (15) into the Taylor series at 0,0, we have(16)w^n+1=1+w1∗αw^n−w1∗βg^n+∑2≤i+j≤3kija,bw^ing^jn+Ow^n+g^n4,g^n+1=b−ξ2g1∗w^n+1−μ2+ξ2g1∗+ξ2N1∗g^n−ξ2w^ng^n−ξ2g^2n+Ow^n+g^n4,where(17)α=a1+g1∗1+N1∗2−ξ1,β=aw1∗1+N1∗2+ξ1,k20a,b=α+12w1∗α2−w1∗1+N1∗α+ξ1,k11a,b=−β−w1∗αβ−aw1∗1+N1∗2+2aw1∗21+N1∗3,k02a,b=12w1∗β2+aw1∗21+N1∗3,k30a,b=12α2+16w1∗α3−w1∗1+N1∗αα+ξ1−1aα+ξ12,k03a,b=−16w1∗β3−aw1∗21+N1∗3β−aw1∗21+N1∗4,k21a,b=−αβ−12w1∗α2β+aw1∗1+g1∗1+N1∗3β−aw1∗2+g1∗−3w1∗21+N1∗3α−a1+N1∗2+4aw1∗1+N1∗3−3aw1∗21+N1∗4,k12a,b=aw1∗1+g1∗−w1∗1+N1∗3β+aw1∗21+N1∗3α+12β2+12w1∗αβ2+aw1∗2+2N1∗−3w1∗1+N1∗4.Fora1,b1,μ1,μ2,ξ1,ξ2∈F1, the Jacobian matrix at E1∗ is(18)JE1∗=1+w1∗α1−w1∗β1b1−ξ2g1∗1−μ2−ξ2N1∗+g1∗,where(19)α1≔a11+g1∗1+N1∗2−ξ1,β1≔a1w1∗1+N1∗2+ξ1.Take two linearly independent (generalized) eigenvectors ofJE1∗,(20)p1=β1w1∗2+α1w1∗1,p2=1b1−ξ2g1∗0,and that of the transposed matrix JTE1∗,(21)q1=b1−ξ2g1∗−2+α1w1∗,q2=01.It is easy to check thatp1,q2=p2,q1=1,p1,q1=p2,q2=0, where x,y≔xTy=x1y1+x2y2 for any x=x1,x2T,y=y1,y2T∈ℝ2. Therefore, any vector w^n,g^nT∈ℝ2 can be represented in the form(22)w^n,g^nT=xnp1+ynp2,which implies that(23)xn=w^n,g^nT,q2,yn=w^n,g^nT,q1.Under new coordinatesxn,yn, map (16) can be written as(24)xn+1=−1+Θ1a,bxn+1−12Θ1a,byn−ξ21+β1w1∗2+α1w1∗x2n−ξ2b1−ξ2g1∗xnyn+Oxn+yn4,yn+1=Θ2a,bxn+−1−12Θ2a,byn+∑2≤i+j≤3rija,bxinyjn+Oxn+yn4,where(25)Θ1a,b=21−b−ξ2g1∗b1−ξ2g1∗,Θ2a,b=−22+w1∗α−2+w1∗α1b−ξ2g1∗b1−ξ2g1∗,r20=b1−ξ2g1∗k11β1w1∗2+α1w1∗+k20β1w1∗2+α1w1∗2+k02+2+α1w1∗ξ2+β1w1∗ξ2,r02=k20b1−ξ2g1∗,r11=k11+2k20β1w1∗2+α1w1∗+2+α1w1∗ξ2b1−ξ2g1∗,r30=b1−ξ2g1∗k30β1w1∗2+α1w1∗3+k21β1w1∗2+α1w1∗2+k12β1w1∗2+α1w1∗+k03,r21=3k30β1w1∗2+α1w1∗2+2k21β1w1∗2+α1w1∗+k12,r12=3k30β1w1∗2+α1w1∗b1−ξ2g1∗+k211b1−ξ2g1∗,r03=k301b1−ξ2g1∗2.By introducing a coordinate transformation(26)xnyn=1−12Θ1a,b0−Θ1a,b1x^ny^n,map (24) can be reduced into the following form:(27)x^n+1=−x^n+y^n+k^20x^2n+k^11x^ny^n+Ox^n+y^n4,y^n+1=Θ2x^n+−1+Θ3y^n+∑2≤i+j≤3r^ija,bx^iny^jn+Ox^n+y^n4,where(28)Θ3a,b=Θ1a,b−12Θ2a,b,k^11k,r=−ξ2b1−ξ2g1∗,r^20a,b=Θ1k^20+1−Θ122r20+Θ12r02−1−Θ12Θ1r11,k^20k,r=−1−Θ12ξ21+β1w1∗2+α1w1∗+Θ1ξ2b1−ξ2g1∗,r^02a,b=r02,r^11a,b=Θ1k^11+1−Θ12r11−2Θ1r02,r^30a,b=1−Θ123r30−Θ13r03−Θ11−Θ122r21+Θ121−Θ12r12,r^21a,b=3Θ12r03+1−Θ122r21−2Θ11−Θ12r12,r^12a,b=1−Θ12r12−3Θ1r03,r^03a,b=r03.Next, we aim to eliminate all quadratic terms and some cubic terms in (27). To achieve it, we introduce the following transformation:(29)x^n=un+∑2≤i+j≤3φija,buinvjn,y^n=vn+∑2≤i+j≤3ψija,buinvjn,where φij and ψij2≤i+j≤3 will be determined later.By applying transformation (29) and its inverse, map (27) becomes(30)un+1=−un+vn+∑2≤i+j≤3ηija,buinvjn+Oun+vn4,vn+1=Θ2a,bun+−1+Θ3a,bvn+∑2≤i+j≤3ζija,buinvjn+Oun+vn4,where(31)η20a,b=k^20+ψ20−2φ20−Θ22φ02+Θ2φ11,η11a,b=k^11+ψ11−2Θ21+Θ3φ02+Θ3−Θ2φ11+2φ20,η02a,b=ψ02−1+1+Θ32φ02+1+Θ3φ11−φ20,ζ20a,b=r^20−Θ22ψ02+Θ2ψ11+Θ2+Θ3ψ20,ζ11a,b=r^11−2Θ21+Θ3ψ02+2−Θ2+2Θ3ψ11+2ψ20+Θ2φ11,ζ02a,b=r^20−Θ31+Θ3ψ02−1+Θ3ψ11−ψ20+Θ2φ02.Assume that the following vanishing conditions hold:(32)η20=η11=η02=ζ20=ζ11=ζ02=0,and we can obtain all the coefficients φij and ψij satisfying i+j=2. Moreover, by taking η30=η12=η21=η03=ζ12=ζ03=0, we can determine the coefficients φij and ψij for i+j=3, see more details in [24, 25]. Thus, model (15) can be transformed into the normal form for 1 : 2 strong resonance bifurcation as follows:(33)un+1=−un+vn+Oun+vn4,vn+1=Θ2a,bun+−1+Θ3a,bvn+Ca,bu3n+Da,bu2nvn+Oun+vn4,where Ca,b and Da,b satisfy(34)Ca1,b1=r^30a1,b1+k^20a1,b1r^20a1,b1+12r^202a1,b1+12r^20a1,b1r^11a1,b1,Da1,b1=r^21a1,b1+12k^20a1,b1r^11a1,b1+54r^20a1,b1r^11a1,b1+r^20a1,b1r^02a1,b1+3k^202a1,b1+52k^20a1,b1r^20a1,b1+52k^11a1,b1r^20a1,b1+r^202a1,b1+12r^112a1,b1.Applying the results established in [29], we can get the existence of 1 : 2 strong resonance bifurcation.Theorem 1.
Assume thatCa1,b1≠0 and Da1,b1+3Ca1,b1≠0. Then, model (2) has a 1 : 2 strong resonance bifurcation at E1∗w1∗,g1∗ when the parameters aandb vary in a sufficiently small neighborhood of F1 with fixed μ1,μ2,ξ1,andξ2. If we further assume Ca1,b1<0, then E1∗w1∗,g1∗ is a saddle; if Ca1,b1>0, then E1∗w1∗,g1∗ is elliptic. Da1,b1+3Ca1,b1 determines the bifurcation scenario under small perturbations. Moreover, model (2) has the following bifurcation behaviors:(i)
There is a pitchfork bifurcation curveP=θ1,θ2:θ1=0, and there exist nontrivial equilibria for θ1<0(ii)
There is a nondegenerate Neimark–Sacker bifurcation curveH=θ1,θ2:θ1=−θ2+Oθ1+θ22,θ1<0(iii)
There is a heteroclinic bifurcation curveHL=θ1,θ2:θ1=−5/3θ2+Oθ1+θ22,θ1<0


## 4. 1 : 3 Strong Resonance

In this section, we prove that there exists a codimension-2 bifurcation with 1 : 3 strong resonance of model (2) at E1∗w1∗,g1∗. We choose a and b as bifurcation parameters of model (2) and assume that a2,b2,μ1,μ2,ξ1,ξ2∈F2.Similar to Section3, we can also transform E1∗w1∗,g1∗ into 0,0 and then expand the right-hand side of the map into Taylor series. For a2,b2,μ1,μ2,ξ1,ξ2∈F2, the Jacobian matrix at E1∗w1∗,g1∗ is(35)JE1∗=1+w1∗α2−w1∗β2b2−ξ2g1∗1−μ2−ξ2N1∗+g1∗,with(36)α2≔a21+g1∗1+N1∗2−ξ1,β2≔a2w1∗1+N1∗2+ξ1,and it has two complex eigenvalues λ1,2=−1/2±3/2i. The corresponding eigenvector ζ1∈ℂ2 and adjoint eigenvector ζ2∈ℂ2 of JE1∗ are defined as(37)ζ1=x1y1≔w1∗β232+w1∗α2−32i,ζ2=x2y2≔3+33+2w1∗α2i6w1∗β2−33i.Any vectorZn=w^n,g^nT∈ℝ2 can be represented as Zn=znζ1+z¯nζ¯1. Thus, for a2,b2,μ1,μ2,ξ1,ξ2∈F2, model (16) can be written in the complex form as follows:(38)zn+1=λ1zn+∑2≤i+j≤3σiji!j!zinz¯jn+Ozn4,where λ1=−1/2+3/2i and(39)σ20=2x¯2k20a2,b2x12+k11a2,b2x1y1+k02a2,b2y12−2ξ2y¯2x1y1+y12,σ02=2x¯2k20a2,b2x¯12+k11a2,b2x¯1y¯1+k02a2,b2y¯12−2ξ2y¯2x¯1y¯1+y¯12,σ11=x¯22k20a2,b2x1x¯1+k11a2,b2x1y¯1+x¯1y1+2k02a2,b2y1y¯1−ξ2y¯2x1y¯1+x¯1y1+2y1y¯1,σ30=6x¯2k30a2,b2x13+k21a2,b2x12y1+k12a2,b2x1y12+k03a2,b2y13,σ03=6x¯2k30a2,b2x¯13+k21a2,b2x¯12y¯1+k12a2,b2x¯1y¯12+k03a2,b2y¯13,σ21=2x¯23k30a2,b2x12x¯1+k21a2,b2x12y¯1+2k21a2,b2x1x¯1y1+k12a2,b2y12x¯1+2x1y1y¯1+3k03a2,b2y12y¯1,σ12=2x¯23k30a2,b2x1x¯12+k21a2,b2x¯12y1+2k21a2,b2x1x¯1y¯1+k12a2,b2y¯12x1+2x¯1y1y¯1+3k03a2,b2y1y¯12.Next, we shall establish the normal form for 1 : 3 strong resonance bifurcation. We first eliminate some quadratic terms in (38) through the transformation(40)zn=ϕn+h202ϕ2n+h11ϕnϕ¯n+h022ϕ¯2n,where the coefficients hij satisfying i+j=2 will be determined later.Using transformation (40) and its inverse transformation,(41)ϕn=zn−h202z2n−h11znz¯n−h022z¯2n+12h202+h11h¯02z3n+h112+32h11h20+12h022z2nz¯n+h02h¯11+h112+12h11h¯20+12h20h02znz¯2n+12h11h02+h02h¯20z¯3n+Ozn4,map (41) is changed into the following form:(42)ϕn+1=λ1ϕn+∑2≤i+j≤3σ^iji!j!ϕinϕ¯jn+Oϕn4,where(43)σ^20=3ih20+σ20,σ^11=−32+32ih11+σ11,σ^02=σ02,σ^30=9−33i2σ20h20+3σ11h¯02+σ30+9+33i2h202−3λ1σ¯02h11,σ^21=2σ11h¯11+σ11h20+σ02h¯02+σ21+1−3iσ11h20+1+3i2σ¯02h02+−3+3ih112+−3+53i2h11h20+1−3ih11σ¯11+5+3i2σ20h11,σ^12=3+3iσ11h11+σ11h¯20+2σ02h¯11+σ20h02+σ12+1−3iσ02h20−3+3ih112−3+3i2h11h¯20+1−3i2h11σ¯20−23ih02h¯11+1+3iσ¯11h02,σ^03=3σ11h02+3σ02h¯20+σ03+31+3i2σ02h11+9−33i2h02h¯20+31+3i2σ¯20h02.Take(44)h20=−σ203i,h11=2σ113−3i,h02=0,then transformation (40) is determined and map (42) becomes(45)ϕn+1=λ1ϕn+σ022ϕ¯2n+∑i+j=3σ^iji!j!ϕinϕ¯jn+Oϕn4.In order to eliminate some cubic terms, we take the transformation(46)ϕn=χn+16h30χ3n+12h21χ2nχ¯n+12h12χnχ¯2n+16h03χ¯3n.It is easy to check that the inverse transformation of (46) is(47)χn=ϕn−16h30ϕ3n−12h21ϕ2nϕ¯n−12h12ϕnϕ¯2n−16h03ϕ¯3n+Oϕn4.By applying (46) and (47), map (45) becomes(48)χn+1=λ1χn+σ022χ¯2n+∑i+j=3σ˜iji!j!χinχ¯jn+Oχn4,where(49)σ˜30=−3−3i2h30+σ^30,σ˜21=σ^21,σ˜12=3ih12+σ^12,σ˜03=−3−3i2h03+σ^03.Take(50)h30=2σ^303−3i,h12=3i3σ^12,h03=2σ^033−3i,h21=0,then transformation (46) is determined and map (48) reduces into the following form:(51)χn+1=λ1χn+σ022χ¯2n+σ^212χn2χ¯n+Oχn4.We defineCa2,b2≔−3/43i+1σ02 and Da2,b2≔−3/4σ022−3/43i+1σ^21. From the 1 : 3 strong resonance bifurcation analysis established in [29], we can get the following theorem.Theorem 2.
Assume thatCa2,b2≠0 and ReDa2,b2≠0. Then, model (2) undergoes a 1 : 3 strong resonance bifurcation at E1∗w1∗,g1∗. ReDa2,b2≠0 determines the stability of the bifurcation invariant closed curve. Moreover, model (2) admits a number of the following complex codimension-1 bifurcation curves:(i)
There is a nondegenerate Neimark–Sacker bifurcation at the trivial fixed pointE0 of (51)(ii)
There is a saddle cycle of periodic three corresponding to saddle fixed pointsEkk=1,2,3 of map (51)(iii)
There is a homoclinic structure formed by the stable and unstable invariant manifolds of the period three cycle intersecting transversally in an exponentially narrow parameter region


## 5. 1 : 4 Strong Resonance

In this section, we prove that there exists a codimension-2 bifurcation with 1 : 4 strong resonance of model (2) at E1∗w1∗,g1∗. We choose a and b as bifurcation parameters of model (2) and assume that a3,b3,μ1,μ2,ξ1,ξ2∈F3.Similar to Section3, we can also transform E1∗w1∗,g1∗ into the origin and then expand the right-hand side of the map into Taylor series. For a3,b3,μ1,μ2,ξ1,ξ2∈F3, the Jacobian matrix at E1∗w1∗,g1∗ is(52)JE1∗=1+w1∗α3−w1∗β3b3−ξ2g1∗1−μ2−ξ2N1∗+g1∗,with(53)α3≔a31+g1∗1+N1∗2−ξ1,β3≔a3w1∗1+N1∗2+ξ1,and it has two complex eigenvalues λ1,2=±i. The corresponding eigenvector ζ˜1∈ℂ2 and adjoint eigenvector ζ˜2∈ℂ2 of JE1∗ are defined as(54)ζ˜1=x˜1y˜1≔w1∗β31+w1∗α3−i,ζ˜2=x˜2y˜2≔1+1+w1∗α3i2w1∗β3−12i.Any vectorZn=x^n,y^nT∈ℝ2 can be represented as Zn=znζ˜1+z¯nζ˜¯1. Thus, for a3,b3,μ1,μ2,ξ1,ξ2∈F3, model (16) can be written in the following complex form:(55)zn+1=λ1zn+∑2≤i+j≤3σiji!j!zinz¯jn+Ozn4,where λ1=i and σij2≤i+j≤3 are defined as in (39) with x1,y1,x2,y2,kija2,b2 replaced by x˜1,y˜1,x˜2,y˜2,kija3,b3, respectively.Similar to Section4, by (40) and (41) we can transform (55) into the following form:(56)ϕn+1=λ1ϕn+∑2≤i+j≤3σ^iji!j!ϕinϕ¯jn+Oϕn4,where(57)σ^20=i+1h20+σ20,σ^11=i−1h11+σ11,σ^02=i+1h02+σ02,σ^30=31−iσ20h20+3σ11h¯02+σ30−3i−1h202−3i+1h11h¯02−3iσ¯02h11,σ^21=2σ11h¯11+1−2iσ11h20+2+iσ20h11+σ02h¯02+σ21+3i+1h20h11+2i−1h112−2ih11σ¯11+i+1h022+iσ¯02h02,σ^12=2i+1σ11h11+σ11h¯20+2σ02h¯11+σ20h02+σ12−i−1h20h02−2iσ02h20−2i+1h112−i+1h11h¯20−ih11σ¯20−2i−1h02h¯11+2iσ¯11h02,σ^03=3σ11h02+3σ02h¯20+σ03+3i−1h11h02+3iσ02h11+3i+1h02h¯20+3iσ¯20h02.Take(58)h20=−σ201+i,h11=σ111−i,h02=−σ021+i,then (56) becomes into the following form:(59)ϕn+1=λ1ϕn+∑i+j=3σ^iji!j!ϕinϕ¯jn+Oϕn4.By (46) and (47), we can transform (59) into the following form as (51):(60)χn+1=λ1χn+σ^212χn2χ¯n+σ^036χ¯3n+Oχn4.LetCa3,b3=−2iσ^21 and Da3,b3=−2/3iσ^03. If Da3,b3≠0, we denote B=Ca3,b3/Da3,b3. From the 1 : 4 strong resonance bifurcation analysis established in [29], we can obtain the following results.Theorem 3.
Assume thatDa3,b3≠0, ReB≠0, and ImB≠0. Then, model (2) has a 1 : 4 strong resonance bifurcation at E1∗w1∗,g1∗. B determines the bifurcation scenario near E1∗. Near the equilibrium E1∗, there exist two-parameter families of equilibria of order four bifurcation from E1∗, and one of them is unstable, or attractive or repelling invariant closed curves, which depends on the choices of parameters a and b. Moreover, there are many complex codimension-1 bifurcation curves of map (60) in a sufficiently small neighborhood of a3,b3.


## 6. Numerical Simulations

In this section, we give a simulation illustration for codimension-2 bifurcations with 1 : 2, 1 : 3, and 1 : 4 strong resonances. The codimension-2 bifurcation diagrams, maximum Lyapunov exponents, and phase portraits of (2) are presented to illustrate the validity of theoretical results and display the complex dynamical behaviors.Case 1.
Choosing the values of the parameters asμ1=8,μ2=2,ξ1=50,andξ2=11, we can numerically obtain a1=230.060722andb1=12.87486 and equilibrium E1∗w1∗,g1∗=0.603987,0.535162. The corresponding eigenvalues of JE1∗ are λ1,2=−1. By further calculations, we have Ca1,b1=2.1132 and Da1,b1=4.2543. From Theorem 1, map (2) is equivalent to the normal form for 1 : 2 strong resonance bifurcation and E1∗ is a 1 : 2 strong resonance bifurcation point.
Figure1(a) shows 2-dimensional bifurcation diagram when b=b1=12.87486 and a varies in 230.060715,230.06075, which is a neighborhood of a1=230.060722. Figure 1(b) shows 3-dimensional bifurcation diagram when a,b vary in a neighborhood of a1,b1=230.060722,12.87486. Figures 2(a)–2(c) show the phase portraits of (2) near equilibrium E1∗ for different parameters a and b. From the phase portraits, we know that there are periodic orbits and chaos at the corresponding parameter values.
The maximum Lyapunov exponents corresponding to the bifurcation diagram Figure1(a) are presented in Figure 2(d). From Figure 2(d), we can see that some maximum Lyapunov exponents are positive and some ones are negative, which means that map (2) has chaotic behaviors near the 1 : 2 strong resonance point E1∗.Figure 1
1 : 2 strong resonance bifurcation diagram at (0.603987, 0.535162): (a) ina−w plane with b=12.87486 and (b) in b−a−w space.
(a)(b)Figure 2
(a) Phase portraits fora=230.060692 and b=12.87486. (b) Phase portraits for a=230.060731 and b=12.87486. (c) Local amplification corresponding to (b). (d) The maximum Lyapunov exponents corresponding to (a).
(a)(b)(c)(d)Case 2.
Choosing the values of the parameters asμ1=8,μ2=2,ξ1=50,andξ2=11, we can numerically obtain a2=248.97816andb2=15.1575017 and equilibrium E1∗w1∗,g1∗=0.61636576,0.605198. The corresponding eigenvalues of JE1∗ are λ1,2=−1/2±3/2i. By further calculations, we have Ca2,b2=1.8958−1.1471i and Da2,b2=−8.362+5.4278i. From Theorem 2, map (2) is equivalent to the normal form for 1 : 3 strong resonance bifurcation and E1∗ is a 1 : 3 strong resonance bifurcation point.
Figure3(a) shows 2-dimensional bifurcation diagram when b=b2=15.1575017 and a varies in 248.97811,248.97822, which is a neighborhood of a2=248.97816. Figure 3(b) shows 3-dimensional bifurcation diagram when a,b vary in a neighborhood of a2,b2=248.97816,15.1575017. Figures 4(a)–4(d) show phase portraits of map (2) near E1∗ for different a and b. From Figure 4(a), we can see that there is a homoclinic structure formed by the stable and unstable invariant manifolds of the period three cycle intersecting transversally in a narrow parameter region.
The maximum Lyapunov exponents corresponding to the bifurcation diagrams Figures3(a) and 3(b) are presented in Figures 5(a) and 5(b), respectively. We can see that some maximum Lyapunov exponents are positive and some ones are negative, which means that map (2) has chaotic behaviors near the 1 : 3 strong resonance point E1∗.Figure 3
1 : 3 strong resonance bifurcation diagram at (0.61636576, 0.605198): (a) ina−w plane with b=15.1575017 and (b) in b−a−w space.
(a)(b)Figure 4
(a) Phase portraits fora=248.97819 and b=15.1575. (b) Phase portraits for a=248.97816 and b=15.157509. (c)-(d) Local amplifications corresponding to (b).
(a)(b)(c)(d)Figure 5
Maximum Lyapunov exponents of map (2) near (0.61636576, 0.605198) as a and b varying: (a) in a−MLE plane with b=15.1575017 and (b) in b−a−MLE space.
(a)(b)Case 3.
Choosing the values of the parameters asμ1=8,μ2=2,ξ1=50,andξ2=11, we can numerically obtain a3=267.20079andb3=17.42246 and equilibrium E1∗w1∗,g1∗=0.627048,0.671077. The corresponding eigenvalues of JE1∗ are λ1,2=±i. By further calculations, we have Ca3,b3=−1.3731−9.4851i and Da3,b3=−4.2121+2.2817i. From Theorem 3, map (2) is equivalent to the normal form for 1 : 4 strong resonance bifurcation and E1∗ is a 1 : 4 strong resonance bifurcation point.
Figure6(a) shows 2-dimensional bifurcation diagram when b=b3=17.42246 and a varies in 267.2004,267.2015, which is a neighborhood of a3=267.20079. From Figure 6(a), we can see that there are flip bifurcations occur after Neimark–Sacker bifurcation near the 1 : 4 strong resonance point. Figure 6(b) shows 3-dimensional bifurcation diagram when a,b vary in a neighborhood of a3,b3=267.20079,17.42246. Figures 7(a)–7(d) show phase portraits of map (2) near E1∗ for different a and b. When a=267.2007andb=17.422455, the equilibrium is locally stable. When aandb increase, map (2) has chaotic behaviors.
The maximum Lyapunov exponents corresponding to the bifurcation diagrams (Figures6(a) and 6(b)) are presented in Figures 8(a) and 8(b), respectively. We can see that some maximum Lyapunov exponents are positive and some ones are negative, which means that map (2) has chaotic behaviors near the 1 : 4 strong resonance point E1∗.Figure 6
1 : 4 strong resonance bifurcation diagram at (0.627048, 0.671077): (a) ina−w plane with b=17.42246 and (b) in b−a−w space.
(a)(b)Figure 7
Phase portraits: (a)a=267.2007 and b=17.422455; (b) a=267.2008 and b=17.422455; (c) a=267.20079 and b=17.42246; (d) a=267.20076 and b=17.42246.
(a)(b)(c)(d)Figure 8
Maximum Lyapunov exponents of map (2) near (0.627048, 0.671077) as a and b varying: (a) in a−MLE plane with b=17.42246 and (b) in b−a−MLE space.
(a)(b)


## 7. Controller Design

In this section, we apply two control strategies (tracking control and hybrid control) to (2). By using tracking control, we can annihilate wide mosquitoes and then prevent the transmission of mosquito-borne diseases. By applying hybrid control, we can advance the codimension-2 bifurcation and extend the range of parameter a from the equilibrium to chaos.


### 7.1. Tracking Control

In this subsection, we design a tracking controller such that the number of wide mosquitoes vanishes gradually as in [23]. Consider a new controlled system(61)wn+1=wnexpawn1+wn+gn−μ1+ξ1wn+gn,gn+1=gn+bwn−μ2+ξ2wn+gngn+Rn,where R is an exterior control signal. When R=0, map (61) has chaos or period orbit under suitable parameter conditions. Let w∗=0 be our control aim, and we shall design a tracking controller R such that wn⟶w∗=0 as the number of iterations is sufficiently large.Theorem 4.
If the controller of controlled map (61) is(62)Rn=−lnϑn+μ1+ξ1+lnϑn+μ1+ξ12+4ξ1aΦ2n−lnϑn−μ12ξ1−Φ2n+gn+bwn−μ2+ξ2wn+gngn,where(63)ϑn=k1exp−awn1+wn+gn+μ1+ξ1wn+gn+k2,Φ2n=wnexpawn1+wn+gn−μ1+ξ1wn+gn,and the constants k1andk2 are positive constants satisfying all roots of z2−k2z−k1=0 lie in the unit circle, then wn⟶w∗=0 as number of iterations tends to infinity.Proof.
Take the following transformation(64)Φ1n=wn,Φ2n=wnexpawn1+wn+gn−μ1+ξ1wn+gn,then map (61) becomes(65)Φ1n+1=wn+1=Φ2n,Φ2n+1=Ψn,where(66)Ψn=wn+1expawn+11+wn+1+gn+1−μ1+ξ1wn+1+gn+1=Φ2nexpaΦ2n1+Φ2n+gn+bΦ1n−μ2+ξ2Φ1n+gngn+Rn−μ1+ξ1Φ2n+gn+bΦ1n−μ2+ξ2Φ1n+gngn+Rn.
Taking (62) into (66), we can obtain the following map:(67)Φ1n+1=Φ2n,Φ2n+1=k1Φ1n+k2Φ2n.
By choosing suitable positive constantsk1andk2 such that all roots of z2−k2z−k1=0 lie in the unit circle, map (61) after feedback is asymptotically stable. That is, Φ2n⟶0 as the number of iterations tends to infinity. Therefore, the trajectory of w tracks an ideal state w∗=0. This completes the proof.
We choose the values of the parameters asμ1=8,μ2=2,ξ1=50,ξ2=11,a=267.2014, and b=17.42246, and take k1=0.75 and k2=0.125 in Rn. When the number of iterations is less than 2000, the chaotic behavior of wn for uncontrolled system (2) is presented in Figure 9(a). When n>2000, the controller is put on (2). From Figure 9(b), we can see that the trajectory of wn tracks ultimately an ideal state w∗=0.Figure 9
(a) Chaotic behavior ofw for uncontrolled system (2). (b) Time responses for the state w of controlled system (61). (c) Bifurcation diagram of controlled system (68).
(a)(b)(c)


### 7.2. Hybrid Control

In this section, we apply hybrid control, which is used in [28, 30], to (2). For any 0<τ<1, consider a new controlled system(68)wn+1=τwnexpawn1+wn+gn−μ1+ξ1wn+gn+1−τwn,gn+1=τgn+bwn−μ2+ξ2wn+gngn+1−τgn,which has the same equilibria with (2). The Jacobian matrix at equilibrium Ew,g is(69)JE=1+τwa1+g1+N2−ξ1−τwaw1+N2+ξ1τb−ξ2g1−τμ2−τξ2N+g,its trace is TrJE=2+τwa1+g/1+N2−ξ1−μ2−ξ2N+g=2+τ/aAN, and determinant is DetJE=1−w/1+Nτ2F′N+τ/aAN. By the analysis in Section 2, we know that only E1∗w1∗,g1∗ may be a 1 : 2, 1 : 3, or 1 : 4 strong resonance point.Here, we only consider the codimension-2 bifurcation with 1 : 4 strong resonance. Takingμ1=8,μ2=2,ξ1=50,ξ2=11, and τ=0.8 in (68), we can numerically obtain a4=258.9192andb4=16.3854 and equilibrium E1∗w1∗,g1∗=0.6234,0.64176, which is a 1 : 4 strong resonance point. The bifurcation diagram of controlled system (68) for fixed b=b4=16.3854 is presented in Figure 9(c) and is very similar to Figure 6(a) of the uncontrolled system. Recall that for the uncontrolled system, we have a3=267.20079andb3=17.42246. Thus, by using hybrid control, we advance the corresponding codimension-2 bifurcation with 1 : 4 strong resonance. Compared with these two figures, we can also find the range of parameter a from 1 : 4 strong resonance point to chaos in Figure 9(c) is about three times larger than that of Figure 6(a).


## 7.1. Tracking Control

In this subsection, we design a tracking controller such that the number of wide mosquitoes vanishes gradually as in [23]. Consider a new controlled system(61)wn+1=wnexpawn1+wn+gn−μ1+ξ1wn+gn,gn+1=gn+bwn−μ2+ξ2wn+gngn+Rn,where R is an exterior control signal. When R=0, map (61) has chaos or period orbit under suitable parameter conditions. Let w∗=0 be our control aim, and we shall design a tracking controller R such that wn⟶w∗=0 as the number of iterations is sufficiently large.Theorem 4.
If the controller of controlled map (61) is(62)Rn=−lnϑn+μ1+ξ1+lnϑn+μ1+ξ12+4ξ1aΦ2n−lnϑn−μ12ξ1−Φ2n+gn+bwn−μ2+ξ2wn+gngn,where(63)ϑn=k1exp−awn1+wn+gn+μ1+ξ1wn+gn+k2,Φ2n=wnexpawn1+wn+gn−μ1+ξ1wn+gn,and the constants k1andk2 are positive constants satisfying all roots of z2−k2z−k1=0 lie in the unit circle, then wn⟶w∗=0 as number of iterations tends to infinity.Proof.
Take the following transformation(64)Φ1n=wn,Φ2n=wnexpawn1+wn+gn−μ1+ξ1wn+gn,then map (61) becomes(65)Φ1n+1=wn+1=Φ2n,Φ2n+1=Ψn,where(66)Ψn=wn+1expawn+11+wn+1+gn+1−μ1+ξ1wn+1+gn+1=Φ2nexpaΦ2n1+Φ2n+gn+bΦ1n−μ2+ξ2Φ1n+gngn+Rn−μ1+ξ1Φ2n+gn+bΦ1n−μ2+ξ2Φ1n+gngn+Rn.
Taking (62) into (66), we can obtain the following map:(67)Φ1n+1=Φ2n,Φ2n+1=k1Φ1n+k2Φ2n.
By choosing suitable positive constantsk1andk2 such that all roots of z2−k2z−k1=0 lie in the unit circle, map (61) after feedback is asymptotically stable. That is, Φ2n⟶0 as the number of iterations tends to infinity. Therefore, the trajectory of w tracks an ideal state w∗=0. This completes the proof.
We choose the values of the parameters asμ1=8,μ2=2,ξ1=50,ξ2=11,a=267.2014, and b=17.42246, and take k1=0.75 and k2=0.125 in Rn. When the number of iterations is less than 2000, the chaotic behavior of wn for uncontrolled system (2) is presented in Figure 9(a). When n>2000, the controller is put on (2). From Figure 9(b), we can see that the trajectory of wn tracks ultimately an ideal state w∗=0.Figure 9
(a) Chaotic behavior ofw for uncontrolled system (2). (b) Time responses for the state w of controlled system (61). (c) Bifurcation diagram of controlled system (68).
(a)(b)(c)


## 7.2. Hybrid Control

In this section, we apply hybrid control, which is used in [28, 30], to (2). For any 0<τ<1, consider a new controlled system(68)wn+1=τwnexpawn1+wn+gn−μ1+ξ1wn+gn+1−τwn,gn+1=τgn+bwn−μ2+ξ2wn+gngn+1−τgn,which has the same equilibria with (2). The Jacobian matrix at equilibrium Ew,g is(69)JE=1+τwa1+g1+N2−ξ1−τwaw1+N2+ξ1τb−ξ2g1−τμ2−τξ2N+g,its trace is TrJE=2+τwa1+g/1+N2−ξ1−μ2−ξ2N+g=2+τ/aAN, and determinant is DetJE=1−w/1+Nτ2F′N+τ/aAN. By the analysis in Section 2, we know that only E1∗w1∗,g1∗ may be a 1 : 2, 1 : 3, or 1 : 4 strong resonance point.Here, we only consider the codimension-2 bifurcation with 1 : 4 strong resonance. Takingμ1=8,μ2=2,ξ1=50,ξ2=11, and τ=0.8 in (68), we can numerically obtain a4=258.9192andb4=16.3854 and equilibrium E1∗w1∗,g1∗=0.6234,0.64176, which is a 1 : 4 strong resonance point. The bifurcation diagram of controlled system (68) for fixed b=b4=16.3854 is presented in Figure 9(c) and is very similar to Figure 6(a) of the uncontrolled system. Recall that for the uncontrolled system, we have a3=267.20079andb3=17.42246. Thus, by using hybrid control, we advance the corresponding codimension-2 bifurcation with 1 : 4 strong resonance. Compared with these two figures, we can also find the range of parameter a from 1 : 4 strong resonance point to chaos in Figure 9(c) is about three times larger than that of Figure 6(a).


## 8. Discussion

In this paper, we analyzed codimension-two bifurcations with 1 : 2, 1 : 3, and 1 : 4 strong resonances of a discrete wild and sterile mosquito model, which is established from continuous case (1) proposed in [7] by the Euler forward difference method. In [10, 11], we have known that continuous model (1) undergoes a sequence of bifurcations such as saddle-node bifurcation, supercritical and subcritical Hopf bifurcation, homoclinic bifurcation, and Bogdanov–Takens bifurcation. Through the theoretical analysis and numerical simulations given in this paper, we find that corresponding discrete mosquito model (2) also has many complex dynamical behaviors, including codimension-1, codimension-2 bifurcations, and chaos. From the impact of the mating rate a and the release rate coefficient b on the dynamics of wide mosquitoes, we know that the release strategy plays an important role in controlling the number of wide mosquitoes. Moreover, two control strategies are applied to model (2). However, whether there is codimension-2 bifurcation with 1 : 1 strong resonance is an interesting problem worth considering.



---

*Source: 3075487-2020-08-17.xml*