Neimark–Sacker Bifurcation and Controlling Chaos in a Three-Species Food Chain Model through the OGY Method
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.
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  and Voltera . Ali Khan et al.  propose a new harvesting strategy for controlling chaotic population in a three-species Hastings and Powell food chain model:where .
Steady state, limit cycle, and period-2 and period-4 behavior of the model are obtained in . This threshold harvesting strategy will be very useful for species conservation and fishery management.
Panday et al.  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:
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.  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:
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.  discuss the qualitative behavior of a four-dimensional discrete-time predator-prey model with parasites.
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.  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. . In 1992, the French physicist Pyragas proposed a time-delayed feedback control method for chaos . 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 , sliding mode control , and adaptive Lyapunov control  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 .
This food chain model describes the insects group of three fully different insects. These insects are a population of prey x, 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  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 (Theorem 1), 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 when a = 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 aswhere . Equation (5) has four equilibria:(i) is the origin; all species are extinct(ii) is the axial fixed point in the absence of midlevel species y = 0 and top-level predator z = 0 that exists for (iii) is the axial fixed point in the absence of top-level predator z = 0 that exists for (iv)The interior (positive) fixed point
Assume is positive and denote equilibrium of equation (5). By the theorem in , 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 aswhere
Consider the matrixwhere
The characteristic equation iswhich may be rewritten in the formwhere
There exist different topological types of for all permissible values of parameters.(i)The eigenvalues of matrix (9) are , 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 , , at E2. When 1 < a < 3 and , , so the equilibrium E2 is a sink. If a > 3 and , , E2 is a saddle point; if a = 3 or , , E2 is a nonhyperbolic point.(iii)The eigenvalues of matrix (9) are , , and at E3. If , are real, and they are all inside the unit circle when and , the equilibrium E3 is asymptotically stable point. The interior equilibrium will be locally asymptotically stable if the coefficients of the characteristic equation (12) , and satisfy the Routh–Hurwitz  stability criterion, i.e., and .(iv)The eigenvalues of matrix (9) at E4 are all inside the unit circle, when the coefficients c1, c2, and c3 of equation (12) are satisfied:
Meanwhile, the equilibrium E4 is a stable point, if and only if the coefficients c1, c2, and c3 of equation (12) are satisfied:
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 parameter a 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 value , then the system shows Neimark–Sacker bifurcation around the positive equilibrium E3 if the following conditions hold:
Proof. The position equilibrium E3 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 such thatWhen , the characteristic equation (14) can be of the formThe above equation has three roots and . To show that Neimark–Sacker bifurcation occurs at , we need to satisfy the conditionFor all , the roots are in general of the formNow, we will verify the conditionSubstituting into equation (18) and calculating the derivative, we havewhereNotice that ; then, we haveSolving from equation (22), we getIf , and , thus the conditions hold and Neimark–Sacker bifurcation occurs at . From numerical simulation, we observed that the Neimark–Sacker bifurcation occurs with respect to the reproduction rate parameter a at the critical value around the equilibrium .
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.
When the reproduction rate a 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).
When the parameter a 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 Figure 4, 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.
The chaotic motion of the three-species food chain model is controlled on the stable periodic period-1 orbit using the method in . 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: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 , at b = 3.7, c = 3, d = 3.5, and r = 3.8. By equation (6), the equilibrium E3 is .
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.
Let denotes the unstable point E3; by first-order Taylor expansion, equation (5) can be described aswhere A is the derivative matrix of F(, a) to the variable , and is the variable set of (x, y, z), ,and B is the derivative matrix of F(, a) to the variable a, ,
Put the equilibrium point E3 into the matrixes A and B:
The time-dependent control parameter a is in the form of a linear function with respect to the variable:
If the modulus of the matrix eigenvalues is less than 1, then the equilibrium point is stable. According to Ogata , we can find the matrix is controllable matrix:
The solution of pole placement is found out through the matrix , which can stabilize the chaotic motion to a stable periodic orbit, where and T = CQ, Q is a matrix of order 3.(i = 1, 2, 3) are the coefficients of the characteristic polynomial of the matrix A; that is,
Put matrix A into the following equation:
So, . Suppose are the coefficients of the characteristic polynomial of the characteristic polynomial :
The eigenvalues of A, which can also be obtained at the equilibrium point , are , and . Assume are the coefficients of the characteristic polynomial of the matrix . The eigenvalues are called adjustment values:
According to the relation between roots and coefficients, we obtain
From equation (33), we know that the matrix is not unique. Let ; then, can be obtained by and equation (35). When goes into this region, whose width is , the parameter is slightly perturbed; otherwise, the parameter is not perturbed. The controlling equation for the parameter a iswhere is a piecewise function and is a small positive number whose value affects the control time.
When , the chaotic motion can be controlled on the period-1 orbit at n = 2000 at (as shown in Figures 5(a)–5(c)). The chaotic motion can be controlled on the period-1 orbit at n = 750 at (as shown in Figures 5(d)–5(f)). According to the numerical simulation, when different values of , , and are taken, there is a great difference in the time of chaos control. When the values of , , and are not selected properly, chaos control may not be realized.
When , the chaotic motion can be controlled on the period-1 orbit at n = 990 at (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.
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 parameters a, 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.
No data, models, or codes were generated or used during the study.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
This work was supported by the National Natural Science Foundation of China (no. 12075043) and the High-Level Scientific Research Project Cultivation Fund of Shandong Women’s University.
A. J. Lotka, Elements of Mathematical Biology, Williams & Wilkins, Baltimore, MD, USA, 1925.
V. Voltera, Opere Matematiche: Mmemorie e Note, Accademia Nazionale dei Lincei, Rome, Italy, 1962.
P. V. Ivanchikov and L. V. Nedorezov, “About a modification of May model of parasite-host system dynamics,” Computational Ecology and Software, vol. 2, no. 1, pp. 42–52, 2012.View at: Google Scholar
H. Li, “Hopf Bifurcation of delayed density-dependent dredator-drey model,” Acta Mathematica Scientia, Series A, vol. 39, no. 2, pp. 358–371, 2019.View at: Google Scholar
E. Ott, C. Grebogi, and J. A. Yorke, “Controlling chaos,” Physical Review Letters, vol. 64, no. 11, pp. 1196–1199, 1990.View at: Google Scholar
A. A. Gomes, E. Manica, and M. C. Varriale, “Applications of chaos control techniques to a three-species food chain,” Chaos, Solitons and Fractals, vol. 36, no. 4, pp. 1097–1107, 2006.View at: Google Scholar
G. Feng, “Controlling chaos of the Ricker population model,” American Journal of Bioscience and Bioengineering, vol. 16, no. 3, pp. 424–431, 2020.View at: Google Scholar
S. Hashemi, M. Ali Pourmina, S. Mobayen, and M. R. Alagheband, “Design of a secure communication system between base transmitter station and mobile equipment based on finite-time chaos synchronisation,” International Journal of Systems Science, Taylor & Francis Journals, vol. 51, no. 11, pp. 1969–1986, 2020.View at: Publisher Site | Google Scholar
S. Mobayen, J. Ma, G. Pujol-Vazquez, L. Acho, and Q. Zhu, “Adaptive finite-time stabilization of chaotic flow with a single unstable node using a nonlinear function-based global sliding mode,” Iranian Journal of Science and Technology, Transactions of Electrical Engineering, vol. 43, no. 1, pp. 339–347, 2019.View at: Publisher Site | Google Scholar
S. Mobayen, K. V. Christos, S. Kaçar, Ü. Çavuşoğlu, and B. Vaseghi, “A chaotic system with infinite number of equilibria located on an exponential curve and its chaos-based engineering application,” International Journal of Bifurcation and Chaos, vol. 28, no. 9, Article ID 1850112, 2018.View at: Publisher Site | Google Scholar
K. Ogata, Controlling Engineering, Prentice-Hall, Englewood Cliffs, NJ, USA, Second edition, 1990.