# Neimark–Sacker Bifurcation and Controlling Chaos in a Three-Species Food Chain Model through the OGY Method

Authors: Guo Feng; Ding Yin; Li Jiacheng
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

The dynamics behavior of a discrete-time three-species food chain model is investigated. By using bifurcation theory, it is shown that the equilibrium point of the system loses its stability, and the system undergoes Neimark–Sacker bifurcation, which leads to chaos as the parameter changes. The chaotic motion is controlled on the stable periodic period-1 orbit using the implementation of the hybrid control strategy. The factor affecting the control time of chaos is also studied. Numerical simulations are consistent with the theoretical analysis. The results of this research prove that the chaos control method can be extended to the higher-dimensional biological model and can be realized.


---

## Body

## 1. Introduction

Predation and prey behaviors are common phenomena in the ecosystem. Since the last century when Volterra and Lotka constructed the predator-prey model, the predator-prey model has been concerned by many scholars. Most scholars are deepening the original model to be more practical. They consider many factors such as time delay, functional response, diffusion, and so on. They also discuss the dynamic properties of the system, such as stability, persistence, and the existence of positive periodic solution, and get a lot of good results.Lotka–Volterra model describes interactions between two species in an ecosystem, a predator and a prey. The model was developed independently by Lotka [1] and Voltera [2]. Ali Khan et al. [3] propose a new harvesting strategy for controlling chaotic population in a three-species Hastings and Powell food chain model:(1)dXdT=R0X1−XK0−C1F1XY,dYdT=F1XY−F2YZ−D1Y,dZdT=C2F2YZ−D2Z,where FiU=AiU/Bi+U,i=1,2.Steady state, limit cycle, and period-2 and period-4 behavior of the model are obtained in [3]. This threshold harvesting strategy will be very useful for species conservation and fishery management.Panday et al. [4] propose a three-species food chain model, where the growth rate of middle predator is reduced due to the cost of fear of top predator, and the growth rate of prey is suppressed due to the cost of fear of middle predator:(2)dXdT=R0X1−XK0−C1A1XYB1+X,dYdT=A1XYB1+X−A2YZB2+Y−D1Y,dZdT=C2A2YZB2+X−D2Z.The global stability analysis of the equilibrium points and the chaotic behavior of the system (2) are investigated using Poincaré section and maximum Lyapunov exponent.The Allee effect is an important phenomenon in population biology characterized by positive density dependence, that is, a positive correlation between population density and individual fitness. Parshad et al. [5] consider a ratio-dependent spatially explicit three-species food chain model, where the top predator is subjected to a strong Allee effect. They show the existence of a global attractor for the following system:(3)∂u∂t=d1Δu+u−u2−w1uvu+v,∂v∂t=d2Δv−a2v+w2uvu+v−w3vrv+r,∂r∂t=d3Δr+rr−mc−w4rv+d3.Most of the research is based on continuous model. In reality, for some species with small number, short life cycle, and no cross between mother and son, it is more realistic to establish discrete population model to describe their growth process. Compared with the continuous models, the discrete-time food chain models can better analyze the dynamic behavior of the system. The bifurcation and other dynamic behaviors of the discrete system will be enriched with the change of parameters. The discrete model can make better use of chaos control strategies and computer simulation.In recent years, more and more scholars have studied the dynamic behavior of discrete ecosystems [6–11]. Saeed et al. [12] discuss the qualitative behavior of a four-dimensional discrete-time predator-prey model with parasites.(4)xn+1=xn+hxn1−xnk−xnxn+wn−βxnwn,yn+1=yn+hβθexnwn−γyn−θynzn+wn,zn+1=yn+hxnzn−ynzn−dzn,wn+1=wn+hynzn−d+αwn.In system (4), it is shown that the system shows Neimark–Sacker bifurcation at a positive fixed point. The system shows chaotic dynamics as bifurcation parameter changes. The chaotic motion of the system is controlled through implementation of hybrid control method.Ott et al. [13] proposed a new chaos control strategy. The method is not to use the existing dynamic control strategies or destroy the conditions for the occurrence of chaotic motion. The method is further improved by Romeiras et al. [14]. In 1992, the French physicist Pyragas proposed a time-delayed feedback control method for chaos [15]. In the following ten years, the research on chaos control and synchronization has been booming, which has rapidly become an important hotspot in the field of chaos research. For example, linear state feedback control [16], sliding mode control [17], and adaptive Lyapunov control [18] are widely used. Recently, the chaos control in biological systems has been studied. Some authors have researched that these chaos control methods can be applied to population dynamics and play nontrivial evolutionary roles in [7–9, 19]. Other methods will change the fixed point of the control or change the original properties of the original system. OGY chaos control is a small real-time disturbance to the parameters to achieve the purpose of control and maintain the original characteristics of the original system. Feng [20, 21] studied the dynamics and chaos control of two-dimensional Hassell recruitment population model and Ricker-type recruitment population model using OGY method. The research on chaos and chaos control has been greatly developed and widely used. A lot of valuable research results have been obtained in [22–26].Using the Euler method, the continuous-time system is discretized. Moreover, the Poincaré map was used to discretize a continuous system, to study the Neimark–Sacker bifurcation, and also to control the system by the OGY approach in [27, 28]. Recently, the design of an explicit analytical expression, the Poincaré map, and its use for the control of the chaotic system were achieved in [29–31].In this paper, the chaotic dynamics of a three-species food chain model is investigated and the chaos is controlled. The dynamic behavior of this model is discussed by using the stability theory of nonlinear difference equation and bifurcation theory in [32].(5)F:xn+1=fxn,yn,zn=axn1−xn−bxnyn,yn+1=gxn,yn,zn=cxnyn−bynzn,zn+1=hxn,yn,zn=rynzn.This food chain model describes the insects group of three fully different insects. These insects are a population of preyx, which is predated by a ﬁrst predator with population y. A third species given by a top predator z that predates on the ﬁrst predator y. Each of the species with nonoverlapping generations affects each other’s population dynamics. Parameters a, b, c, d, and r are all positive. Parameter a is a reproduction rate of x; b is the predation rate of predator y to prey x and top predator z to y; c is the growth rate of predators y; and r is the growth rate of predators z due to the consumption of species y.Elsadany [32] has studied the stability or bifurcation of the equilibrium, but has not studied the condition of Neimark–Sacker bifurcation or the control of chaos. In the paper, the conditions of Neimark–Sacker bifurcation are given and proved strictly. By studying chaos control of the model, we can better understand the dynamic behavior of the biological system. We can adjust the value of parameters through chaos control to make the ecological model in a state of equilibrium. Numerical simulation verifies the correctness of our theoretical derivation.For the study of Neimark–Sacker bifurcation, most of the previous studies have used the central manifold theorem to obtain the bifurcation parameter values after tedious derivation. In the paper, a method is proposed to judge the Neimark–Sacker bifurcation (Theorem1), which makes the operation and derivation simple and clear, and gives the bifurcation parameter values.The discrete-time three-species food chain model produces a chaotic attractor whena = 4.26, b =3.7, c =3, d = 3.5, and r = 3.8. The unstable point is controlled on the period-1 orbit. It shows the control times are different as the adjustment values to be selected are different. In biology or ecology, the complex chaotic behavior of this model shows the relationship of the different species, including the number population, reproduction rate, and survival rate, whether they can survive in a balanced state, or makes the population develop in disorder or chaos.


## 2. The Equilibrium Point and Local Stability Analysis of the Food Chain Model

In order to find out the equilibria of this system and study its dynamic properties. Equation (5) is written as(6)fx,y,z=x,gx,y,z=y,hx,y,z=z,where fx,y,z=ax1−x−bxy,gx,y,z=cxy−byz,hx,y,z=ryz. Equation (5) has four equilibria:(i)
E1=0,0,0 is the origin; all species are extinct(ii)
E2=a−1/a,0,0 is the axial fixed point in the absence of midlevel species y = 0 and top-level predator z = 0 that exists for a>1(iii)
E3=1/c,a/b1−1/c−1/b,0 is the axial fixed point in the absence of top-level predator z = 0 that exists for a>c−1/c(iv)
The interior (positive) fixed pointE4=ra−1−b/ar,1/r,rac−a−c−cb/adrAssumeEix∗,y∗,z∗,i=1,2,3,4 is positive and denote equilibrium of equation (5). By the theorem in [20], the equilibria are classified according to the eigenvalues of linearized matrix of equation (5). To carry out linear stability analysis, the Taylor series expansion of equation (5) may be written as(7)xn+1yn+1zn+1=∂f∂x∂f∂y∂f∂z∂g∂x∂g∂y∂g∂z∂h∂x∂h∂y∂h∂zx∗,y∗,z∗xnynzn,where (8)∂f∂x|x∗,y∗,z∗=a1−2x∗−bx∗+y∗,∂f∂y|x∗,y∗,z∗=−bx∗,∂f∂z|x∗,y∗,z∗=0,∂g∂x|x∗,y∗,z∗=cy∗,∂g∂y|x∗,y∗,z∗=cbx∗−bz∗,∂g∂z|x∗,y∗,z∗=−by∗,∂h∂x|x∗,y∗,z∗=0,∂h∂y|x∗,y∗,z∗=rz∗,∂h∂z|x∗,y∗,z∗=ry∗.Consider the matrix(9)J=P11P12P13P21P22P23P31P32P33,where(10)P11=∂f∂x|x∗,y∗,z∗,P12=∂f∂y|x∗,y∗,z∗,P13=∂f∂z|x∗,y∗,z∗=0,P21=∂g∂x|x∗,y∗,z∗,P22=∂g∂y|x∗,y∗,z∗,P23=∂g∂z|x∗,y∗,z∗,P31=∂h∂x|x∗,y∗,z∗=0,P32=∂h∂y|x∗,y∗,z∗,P33=∂h∂z|x∗,y∗,z∗.The characteristic equation is(11)detP11−λP12P13P21P22−λP23P31P32P33−λ=0,which may be rewritten in the form(12)λ3+c1λ2+c2λ+c3=0,where(13)c1=−P11+P22+P33,c2=P11P22−P12P21+P11P33+P22P33−P23P32,c3=−P11P22P33+P11P23P32+P12P21P33.There exist different topological types ofEix∗,y∗,z∗,i=1,2,3,4 for all permissible values of parameters.(i)
The eigenvalues of matrix (9) are λ11=a, λ12=λ13=0 at E1, with eigenvectors, respectively, given by d1 = (1,0,0), d2 = (0,1,0), and d3 = (0,0,3). When a < 1, the equilibrium point E1 is a stable point; a > 1, the equilibrium E1 is an unstable point; otherwise, when a = 1, E1 is called nonhyperbolic point.(ii)
The eigenvalues of matrix (9) are λ21=2−a, λ22=ac−c/a, λ23=0 at E2. When 1 < a < 3 and c/c−1<a<c/c+1, λ21<1,λ22<1, so the equilibrium E2 is a sink. If a > 3 and a>c/c+1, λ21>1,λ22>1, E2 is a saddle point; if a = 3 or a=c/c+1, λ21=1,λ22=1, E2 is a nonhyperbolic point.(iii)
The eigenvalues of matrix (9) are λ31=abrc−1/bc, λ32=1−a/2c+a/c+22−4a, and λ33=1−a/2c−a/c+22−4a at E3. If a/c+22−4a>0, λ31,λ32,λ33 are real, and they are all inside the unit circle when cr−b/brc−1<a<cr+b/brc−1 and a/2a−1<c<2a/a−1, the equilibrium E3 is asymptotically stable point. The interior equilibrium will be locally asymptotically stable if the coefficients of the characteristic equation (12) c1,c2, and c3 satisfy the Routh–Hurwitz [33] stability criterion, i.e., c1>0,c3>0 and c1c2−c3>0.(iv)
The eigenvaluesλ41,λ42,λ43 of matrix (9) at E4 are all inside the unit circle, when the coefficients c1, c2, and c3 of equation (12) are satisfied:(14)c1=2ar3+a2br−a2r3+2a2r2a2r2,c2=ab2c+acr2+2abcr+2a2r3−a2cr2−4ar3−a2bcr−2a2bra2r2,c3=abcr2+2a2bcr+3acr3−2abcr−2cr3−2bcr2−a2cr3−2ab2ca2r2.Meanwhile, the equilibriumE4 is a stable point, if and only if the coefficients c1, c2, and c3 of equation (12) are satisfied:(15)3+c1−c2−2c3>0,1−c2+c3c1−c3>0,1−c1+c2−c3>0.


## 3. Neimark–Sacker Bifurcation Analysis

Neimark–Sacker bifurcation is a critical value of a parameter where the system stability changes and periodic solution arises. In the next theorem, the existence of Neimark–Sacker bifurcation is given, and we take the reproduction rate parametera of the prey x as a bifurcation parameter. In the second part, we have analyzed the local stability near the positive equilibrium in detail. Here, without losing generality, we give the Neimark–Sacker bifurcation conclusion around the equilibrium E3 and prove it.Theorem 1.
When the reproduction rate parameter a of the prey x crosses a critical valuea0, then the system shows Neimark–Sacker bifurcation around the positive equilibrium E3 if the following conditions hold:(16)c1a0>0,c3a0>0,c1a0c2a0−c3a0=0,c1a0c2a0′≠c3′a0.Proof.
The position equilibriumE3 is locally asymptotically stable and the system (6) will lose its stability when some parameter values are changed. Hence, we take the reproduction rate parameter a of the prey x as a bifurcation parameter. If there exists a critical value, a critical value a0 such that(17)c1a0c2a0−c3a0=0.
Whena=a0, the characteristic equation (14) can be of the form(18)λ2a0+c2a0λa0+c1a0=0.
The above equation has three roots−c1a0,ic2a0 and −ic2a0. To show that Neimark–Sacker bifurcation occurs at a0, we need to satisfy the condition(19)dReλada|a=a0≠0.
For alla0, the roots are in general of the form(20)λ1a=μa+iνa,λ2a=μa−iνa,λ3a=−c1a.
Now, we will verify the condition(21)dReλkada|a=a0≠0,k=1,2.
Substitutingλka=μa±iνa into equation (18) and calculating the derivative, we have(22)Kaμ′a−Laν′a+Ma=0,Kaμ′a+Laν′a+Na=0,where(23)Ka=3μ2a+2c1aμa+c2a−3ν2a,La=6μaνa+2c1aνa,Ma=μ2ac1′a+c2′aμa+c3′a−c1′aν2a,Na=2μaνac1′a+c2′aνa.
Notice thatμa0=0,νa0=λ2a0; then, we have(24)Ka0=−2c1aμa0,La0=2c1a0λ2a0,Ma0=λ3′a0−c1′a0c2a0,Na0=c2′a0λ2a0.
Solvingμ′a0 from equation (22), we get(25)dReλkada|a=a0=μ′aa=a0=−La0Na0+Ka0Ma0L2a0+K2a0=12λ3′a0−c1a0c2a0′c12a0+c2a0≠0,k=1,2.
Ifc1a0c2a0′≠c3′a0, and c3a0=−λ1a0, thus the conditions hold and Neimark–Sacker bifurcation occurs at a=a0. From numerical simulation, we observed that the Neimark–Sacker bifurcation occurs with respect to the reproduction rate parameter a at the critical value a0=3.01 around the equilibrium 0.5062,0.2632,0.1743.


## 4. Numerical Simulations

The dynamic behavior of this system is studied by numerical analysis. The bifurcation diagrams of equation (5) are shown in Figures 1(a)–1(c) for each species x, y, and z with b = 3.7, c = 3, d = 3.5, and r = 3.8, as the parameter a changes from 2 to 4.3. We now discuss the salient features of the bifurcation diagram. Because of the similarity of the bifurcation diagrams, only Figure 1(b) is shown magnified in Figure 2 to analyze the dynamic properties.Figure 1
Bifurcation diagrams of equation (5) (b = 3.7, c = 3, d = 3.5, r = 3.8, and 2 < a < 4.3).
(a)(b)(c)Figure 2
Magnification of part of Figure1(b).
(a)(b)When the reproduction ratea changes between 2 and 4.3, equation (5) generates complicated features. The equilibrium point E3 is a positive fixed point that is asymptotically stable for 2.94 < a < 2.99. The phase portraits of various a corresponding to Figure 2(a) are plotted in Figures 3(a) and 3(b). When the parameter a changes from 2.98 to 3.01, the equilibrium point E3 becomes unstable, and at a = 3.01, a Neimark–Sacker bifurcation appears. The size and smoothness of a closed curve around the equilibrium point change continuously with the change of parameters for 3.01 < a < 3.6, and finally the closed curve converges to the equilibrium point E3 again at a = 3.6. The phase portraits of various a corresponding to Figure 2(a) are plotted in Figures 3(c)–3(e).Figure 3
The phase diagrams of equation (5) ((a, b) 2.94 < a < 2.99, (c) a = 3.01, (d, e) 3.01 < a < 3.6, (f) a = 3.87, (g) a = 3.904, and (h–j) 3.904 < a < 3.998).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)When the parametera passes through the range (3.6, 3.87), the system has an unstable period-1 orbit, and a Neimark–Sacker bifurcation occurs again at a = 3.87. There is also a closed curve whose shape and size are constantly changing for 3.87 < a < 3.904. At a = 3.904, this closed curve breaks, and the system shows period-17 orbits. With the parameter a changing from 3.904 to 3.998, the system shows a cascade of period-doubling bifurcations leading to a chaotic attractor finally. The phase portraits of various a corresponding to Figure 2(b) are plotted in Figures 3(f)–3(j).


## 5. The Controlling Chaos of the Three-Species Food Chain Model

As shown in Figure4, for b = 3.7, c = 3, d = 3.5, and r = 3.8, when a = 4.26, the dynamic behavior of the system is chaotic. At the same time, there is a chaotic attractor, which is the closure of the unstable manifolds of the saddle points. And there is an infinite number of unstable periodic orbits in the chaotic attractor. There is an unstable period-1 orbit embedding in the chaotic attractor.Figure 4
The chaotic attractor of the food chain model with an unstable fixed point (b = 3.7, c = 3, d = 3.5, r = 3.8, and a = 4.26).The chaotic motion of the three-species food chain model is controlled on the stable periodic period-1 orbit using the method in [18]. Control parameter a is perturbed slightly with time. When the mapping point wanders to the neighborhood of the periodic-1 orbit, the control parameter is perturbed.Assume the system is written in the following form:(26)Wi+1=FWi,a,Wi∈R3,a∈R,F is sufficiently smooth, and a is an externally adjustable real parameter. Let the control parameter a be a variable parameter near the rated value a¯=4.26, at b = 3.7, c = 3, d = 3.5, and r = 3.8. By equation (6), the equilibrium E3 is x∗,y∗,z∗=0.5062,0.2632,0.1743.Now, the aim is to change the parameter such that the chaotic attractor involves almost all of the initial conditions, so that the dynamic behavior of the system converges to the desired periodic orbit. Due to the ergodicity of the chaos dynasmics, when the state trajectory enter the vicinity of the unstable periodic orbit to be stabilized, a feedback control law is applied to control the trajectory to move to the desired unstable periodic orbit.LetW∗a¯ denotes the unstable point E3; by first-order Taylor expansion, equation (5) can be described as(27)Wi+1−W∗a¯=AWi−W∗a¯+Ba−a¯,where A is the derivative matrix of F(w, a) to the variable w, and w is the variable set of (x, y, z), A=DwFw,a,(28)A=a1−2x−by−bx0cycx−dz−dy0rzry,and B is the derivative matrix of F(w, a) to the variable a, B=DaFw,a,(29)B=x1−x−by00.Put the equilibrium pointE3x∗,y∗,z∗=0.5062,0.2632,0.1743 into the matrixes A and B:(30)A=−1.027−1.87300.78960.9086−0.921200.66231.000,B=−0.723900.The time-dependent control parametera is in the form of a linear function with respect to the variable:(31)a−a¯=−KTWi−W∗a¯.Substituting equation (32) into equation (28),(32)Wi+1−W∗a¯=A−BKTWi−W∗a¯.If the modulus of the matrixA−BKT eigenvalues is less than 1, then the equilibrium point W∗a¯ is stable. According to Ogata [34], we can find the matrix C3×3 is controllable matrix:(33)C=BABA2B=−0.72390.74230.30710−0.57160.067700−0.3786.The solution of pole placement is found out through the matrixKT, which can stabilize the chaotic motion to a stable periodic orbit, where KT=α3−a3α2−a2α1−a1T−1 and T = CQ, Q is a matrix of order 3.(34)Q=a2a11a110100,ai(i = 1, 2, 3) are the coefficients of the characteristic polynomial of the matrix A; that is,(35)detλI−A=λ3+a1λ2+a2λ+a3.Put matrixA into the following equation:(36)detλI−A=detλ+1.0271.8730−0.7896λ−0.90860.92120−0.6623λ−1.000=−λ3+0.8816λ2−1.038λ−0.0808.So,a1=0.8816,a2=−1.038,a3=−0.0808. Suppose α1,α2,α3 are the coefficients of the characteristic polynomial of the characteristic polynomial detA−BKT:(37)detsI−A−BKT=s3+α1s2+α2s+α3,T=CQ=−0.72390.74230.30710−0.57160.067700−0.3786−1.0380.881610.881610100=1.71290.1041−0.7239−0.4362−0.57160−0.378600,T−1=00−2.64130−1.74952.0157−1.3814−0.2516−5.9601.The eigenvalues ofA, which can also be obtained at the equilibrium point x∗,y∗,z∗=0.506,0.263,0.174, are λs=−0.073,λ1u=0.477+0.937i, and λ2u=0.477−0.937i. Assume α1,α2,α3 are the coefficients of the characteristic polynomial of the matrix A−BKT. The eigenvalues μ1,μ2,μ3 are called adjustment values:(38)detsI−A−BKT=s3−μ1+μ2+μ3s2+μ1μ2+μ2μ3+μ1μ3s−μ1μ2μ3.According to the relation between roots and coefficients, we obtain(39)α1=−μ1+μ2+μ3,α2=μ1μ2+μ2μ3+μ1μ3,α3=−μ1μ2μ3.From equation (33), we know that the matrix KT is not unique. Let α1,α2,α3=−λs,−λs,0; then, (40)KT=0−a3−λs−a2−λs−a1T−1=0.0808,0.965,0.954600−2.64130−1.74952.0157−1.3814−0.2516−5.9601=−1.3814−1.9284−3.9578,KTWi−W∗a¯<δ can be obtained by a−a¯<δ and equation (35). When Wi goes into this region, whose width is 2δ/KT, the parameter is slightly perturbed; otherwise, the parameter is not perturbed. The controlling equation for the parameter a is(41)a−a¯=−KTWi−W∗a¯×uδ−KTZi−Z∗a¯,where u is a piecewise function and δ is a small positive number whose value affects the control time.(42)uα=0,α<0,1,α>0.Whenδ=0.036, the chaotic motion can be controlled on the period-1 orbit at n = 2000 at α1=−0.073,α2=−0.073,α3=0 (as shown in Figures 5(a)–5(c)). The chaotic motion can be controlled on the period-1 orbit at n = 750 at α1=−0.056,α2=0.433,α3=−0.025 (as shown in Figures 5(d)–5(f)). According to the numerical simulation, when different values of α1, α2, and α3 are taken, there is a great difference in the time of chaos control. When the values of α1, α2, and α3 are not selected properly, chaos control may not be realized.Figure 5
The phase diagrams of equation (5) ((a, b) 2.94 < a < 2.99, (c) a = 3.01, (d, e) 3.01 < a < 3.6, (f) a = 3.87, (g) a = 3.904, and (h–j) 3.904 < a < 3.998).
(a)(b)(c)(d)(e)(f)Whenδ=0.063, the chaotic motion can be controlled on the period-1 orbit at n = 990 at α1=−0.073,α2=−0.073,α3=0 (as shown in Figures 6(a)–6(c)). According to the numerical simulation, when different values of δare taken, there is a great difference in the time of chaos control.Figure 6
Controlling of period-1 of the food chain modelδ=0.063.
(a)(b)(c)


## 6. Conclusion

It is shown that the equilibrium point of the system loses its stability and Neimark–Sacker bifurcation appears leading to chaos for some parameter values. The chaotic food chain is controlled on the stable periodic period-1 to obtain steady-state orbit using implementation of hybrid control strategy, which is based on feedback control methodology and parameter perturbation. Numerical simulations are presented to illustrate our results with the theoretical analysis and show the effect of the control method. The results of this research prove that the chaos control method can be extended to the higher-dimensional biological model and can be realized. In the study of bifurcation, it is easier to obtain bifurcation parameters than the central manifold theorem. Compared with other methods, chaos control will save more control time and be easier to realize chaos control.This paper only studies the dynamic behavior of the system in mathematics. In ecology and biology, how to control the parametersa, b, c, d, and r properly to achieve an ecological balance, so as not to make a species extinct or to cause damage to the ecological environment will be a very meaningful and challenging topic, which will be paid more attention in the future research. When the parameter is in a certain range, there should be subcritical or supercritical bifuraction. Whether there will be flip bifurcation is also a question worthy of discussion. I will discuss this problem in the future research work.



---

*Source: 6316235-2021-06-28.xml*