# Bifurcation Analysis of a Discrete Food Chain Model with Harvesting

Authors: Abdul Qadeer Khan; Shahid Mehmood Qureshi
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

We explore existence of fixed points, topological classifications around fixed points, existence of periodic points and prime period, and bifurcation analysis of a three-species discrete food chain model with harvesting. Finally, theoretical results are numerically verified.


---

## Body

## 1. Introduction

Many different types of interactions exist in nature between various species of organisms on this planet Earth and are studied under the discipline of ecology. Ecological interactions are most fundamental part in biology which determines community structure and development. Not all interactions are positive, some are negative also. One of the examples of negative correlation is ammensalism. Ammensalism is a type of ecological interactions between the members of two different species in which one is harmed, destroyed, or inhabited by the member of another species, while the other remains unaffected, neither harmed nor benefitted. It is a type of competitive behavior among different species and is frequently used to refer to asymmetrical competitive association. Research in the field of ecology draws the attention of several mathematicians such as Lotka [1] and Volterra [2]. Nowadays, ecologist and mathematician jointly contributed to the growth of this area of knowledge. Recently, many researchers have investigated the dynamical properties of discrete-time ecological models such as prey-predation, competitions, neutralism, and mutualism by studying fixed points, local and global attractivity, bounded, existence of bifurcation, and many more. For instance, Beddington et al. [3] have explored the behavior of following predator-prey model:(1)xn+1=xner1−xn/K−ayn,yn+1=αxn1−e−ayn.Chen [4] has explored global attractivity and permanence of the following discrete multispecies system:(2)xin+1=xinebin−∑l=1kailnxln−∑l=1kcilnyln,yjn+1=yjne−rjn+∑l=1kdjlnxln−∑l=1kejlnyln.Fang and Chen [5] have explored the permanence of multispecies Lotka–Volterra competition predator-prey system with delays. Furthermore, Fang et al. [6] have explored the dynamics of the following system:(3)xn+1=xnean−bnxn−cnyn/m1n+m2nxn+m3nyn,yn+1=yne−dn−enyn+fnxn/m1n+m2nxn+m3nyn.Agiza et al. [7] have explored chaotic dynamics of the following discrete model with Holling type II:(4)xn+1=axn1−xn−bxnyn1+εxn,yn+1=dxnyn1+εxn.Huo and Li [8] have explored stable periodic solution of the following discrete model:(5)xn+1=xner1n−b1nxn−a1nyn,yn+1=yner2n−a2nyn/xn.Lu and Zhang [9] have studied the permanence and global attractivity of the following discrete system:(6)xn+1=xnean−bnxn−mnyn/An+xn,yn+1=ynedn−enyn/xn.Zhao and Zhang [10] explored the chaos and permanence of the following discrete model:(7)xn+1=1−dxner1−xn/k−ayn+dσxn−1er1−xn−1/k,yn+1=1−dxn1−e−ayn.Zhao et al. [11] have investigated the dynamics of the following discrete model:(8)xn+1=xner1−xn/kfyn,yn−1,yn+1=xn1−fyn,yn−1,where fyn,yn−1=e−a1−dyn+dσyn−1.On the contrary, in recent years, many papers have been published to investigate the bifurcation analysis of certain discrete models by choosing step size as a bifurcation parameter. For example, Salman et al. [12] have explored bifurcation analysis of the following discrete system:(9)xn+1=xn+δ2xn1−xn2−yn,yn+1=yn+δyn−s+cxn,by choosing step size δ as a bifurcation parameter. Liu and Xiao [13] have explored bifurcation analysis of the following discrete system:(10)xn+1=xn+δrxn1−xn−bxnyn,yn+1=yn+δ−d+bxnyn,by choosing step size δ as a bifurcation parameter. Hasan and Hama [14] have explored bifurcation analysis of the following discrete system:(11)xn+1=xn+dxn1−xn−yn1+axn+byn,yn+1=yn+dyncxn1+axn+byn−e,by choosing step size d as a bifurcation parameter. Wu and Zhang [15] have explored bifurcation analysis of the following discrete system:(12)xn+1=xn+δxnK1−α1xn−β12yn−γ1xnyn,yn+1=yn+δynK2−α2yn−β21xn−γ2xnyn,by choosing step size δ as a bifurcation parameter. Rana [16] has explored bifurcation analysis of the following discrete system:(13)xn+1=xn+δxn1−xn−axnynxn+yn,yn+1=xn+δ−dyn+bxnynxn+yn,by choosing step size δ as a bifurcation parameter. Rana and Kulsum [17] have explored bifurcation analysis of the following discrete system:(14)xn+1=xn+δxn1−xn−ynxn2+a,yn+1=yn+δynα−βynxn,by choosing step size δ as a bifurcation parameter. Motivated from the aforementioned studies, in this work, we explore existence of fixed points, topological classifications around fixed points, periodic points, and bifurcation analysis, by choosing step size h as a bifurcation parameter, of the following three species discrete food chain model with harvesting:(15)xn+1=xn+ha11−k1xn−α11xn2,yn+1=yn+ha21−k2yn−α22yn2−α21xnyn,zn+1=zn+ha3zn−α33zn2−α32ynzn,which is a discrete form of the following model:(16)dxdt=a11−k1x−α11x2,dydt=a21−k2y−α22y2−α21xy,dzdt=a3z−α33z2−α32yz,by Euler forward formula, where h is step size and t is customary denoted by n. It is noted that, in model (16), xt, yt, and zt, respectively, denote populations of first, second, and third species. Moreover zt denotes growth rate of first, second, and third species; αiii=1,2,3 denotes the rate of decrease of first, second, and third species due to internal competitions; α21 denotes rate of decrease of second species due to attack of first species; α32 denotes the rate of decrease of third species due to attack of second species; k1 and k2, respectively, denote the harvesting rate of first and second species. It is also important to note that all parameters h,a1,a2,a3,k1,k2,α11,α22,α21,α32, and α33 are positive. In addition, it is important here to mention that we will explore dynamical properties of the discrete-time model (15) instead of the continuous-time model, which is depicted in (16), because discrete-time models governed by difference equations are more realistic and appropriate than the continuous ones in the case where populations have nonoverlapping generations, and moreover, discrete models can also provide efficient computational results for numerical simulations [12, 13].This paper is structured as follows. In Section2, we study the existence of fixed points of model (15) algebraically. The linearized form of model (15) is presented in Section 3. In Section 4, we explored topological classification around fixed points of the model. Existence of periodic points of model (15) is explored in Section 5. In Section 6, we explored detailed analysis of bifurcation around fixed points of model (15). Theoretical results are verified numerically in Section 7. Brief summary of the paper is presented in Section 8.


## 2. Study of Equilibrium Points

Here, we will study the boundary and interior equilibrium points of model (15) as follows.Lemma 1.
Model (15) has atmost eight equilibrium points in ℝ+3. Precisely,(i)
∀h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33>0; model (15) has a trivial equilibrium point: P1=0,0,0.(ii)
∀a3,α33>0; model (15) has boundary equilibrium point: P2=0,0,a3/α33.(iii)
P3=0,1−k2a2/α22,0 is a boundary equilibrium point of (15) if k2<1.(iv)
P4=1−k1a1/α11,0,0 is a boundary equilibrium point of (15) if k1<1.(v)
P5=0,1−k2a2/α22,a3α22−α321−k2a2/α22α33 is a boundary equilibrium point of (15) if a3>α321−k2a2/α22 with k2<1.(vi)
P6=1−k1a1/α11,0,a3/α33 is a boundary equilibrium point of (15) if k1<1.(vii)
P7=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,0 is a boundary equilibrium point of (15) if a2>a11−k1α21/1−k2α11 with k1,k2<1.(viii)
P8=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,a3α11α22−a21−k2α11α32+a11−k1α21α32/α11α22α33 is an interior equilibrium point of (15) if k1<1, a2>a11−k1α21/1−k2α11 and a3>a21−k2α11α32−a11−k1α21α32/α11α22.Proof.
If model (15) has an equilibrium point, P=x,y,z, then(17)x=x+ha11−k1x−α11x2,y=y+ha21−k2y−α22y2−α21xy,z=zn+ha3z−α33z2−α32yz.
The simple computation yields that, for the values ofPii=1,…,7, (17) satisfied identically. So, one can conclude that model (15) has seven boundary points: Pii=1,…,7. In order to find interior point, from (17), one obtains(18)a11−k1−α11x=0,a21−k2−α22y−α21x=0,a3−α33z−α32y=0.
From1st of (18), one obtains(19)x=1−k1a1α11.
From2nd equation of (18) and (19), one obtains(20)y=a21−k2α11−a11−k1α21α11α22.
From3rd equation of (18) and (20), one obtains(21)z=a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33.
From (19)–(21), one can conclude that P8 is an interior equilibrium point of (15) if k1<1, a2>a11−k1α21/1−k2α11, and a3>a21−k2α11α32−a11−k1α21α32/α11α22.


## 3. Linearized Form of Model (15)

The variational matrixJ|P about P under the map:(22)f1,f2,f3↦xn+1,yn+1,zn+1,where(23)f1=x+ha11−k1x−α11x2,f2=y+ha21−k2y−α22y2−α21xy,f3=z+ha3z−α33z2−α32yz,is(24)J|P=1+ha11−k1−2α11x00−hα21y1+ha21−k2−2α22y−α21x00−hα32z1+ha3−2α33z−α32y,with(25)λ1=1+ha11−k1−2α11x,λ2=1+ha21−k2−2α22y−α21x,λ3=1+ha3−2α33z−α32y.


## 4. Dynamical Behavior: Topological Properties of Equilibrium Points

The dynamical behavior about fixed pointsPii=1,…,8 of model (15) is explored in this section.


### 4.1. Dynamical Behavior aboutP1

From (25), eigenvalues of J|P1 about P1 are(26)λ1=1+ha11−k1,λ2=1+ha21−k2,λ3=1+ha3.The dynamical behavior aboutP1 of model (15) is concluded as follows.Lemma 2.
(i)
For all allowed parametric values,h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33>0, P1 is not sink.(ii)
P1 is a source if(27)0<h<min2a1k1−1,2a2k2−1.(iii)
P1 is a saddle if(28)h>max2a1k1−1,2a2k2−1.(iv)
P1 is nonhyperbolic if(29)h=2a1k1−1,or(30)h=2a2k2−1.


### 4.2. Dynamical Behavior aboutP2

From (25), eigenvalues of J|P2 about P2 are(31)λ1=1+ha11−k1,λ2=1+ha21−k2,λ3=1−ha3.The dynamical behavior aboutP2 is concluded as follows:Lemma 3.
(i)
P2 is a sink if(32)h>max2a1k1−1,2a2k2−1,0<h<2a3.(ii)
P2 is a source if(33)0<h<min2a1k1−1,2a2k2−1,h>2a3.(iii)
P2 is a saddle if(34)h<min2a1k1−1,2a2k2−1,0<h<2a3.(iv)
P2 is nonhyperbolic if(35)h=2a3,or(36)h=2a1k1−1,or(37)h=2a2k2−1.


### 4.3. Dynamical Behavior aboutP3

From (25), eigenvalues of J|P3 about P3 are(38)λ1=1+ha11−k1,λ2=1−ha21−k2,λ3=1+ha3α22−a2α321−k2α22.The dynamical behavior aboutP3 is concluded as follows.Lemma 4.
(i)
P3 is a sink if(39)h>max2a1k1−1,2α22a2α321−k2−a3α22,0<h<2a21−k2.(ii)
P3 is a source if(40)0<h<min2a1k1−1,2α22a2α321−k2−a3α22,h>2a21−k2.(iii)
P3 is a saddle if(41)h>max2a1k1−1,2a21−k2,0<h<2α22a2α321−k2−a3α22.(iv)
P3 is nonhyperbolic if(42)h=2a21−k2,or(43)h=2a1k1−1,or(44)h=2α22a2α321−k2−a3α22.


### 4.4. Dynamical Behavior aboutP4

From (25), eigenvalues of J|P4 about P4 are(45)λ1=1−ha11−k1,λ2=1+ha21−k2,λ3=1+ha3.The dynamical behavior aboutP4 is concluded as follows.Lemma 5.
(i)
For all allowed parametric values,h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33>0, P4 is not sink.(ii)
P4 is a source if(46)h>max2a11−k1,2a2k2−1.(iii)
P4 is a saddle if(47)0<h<min2a11−k1,2a2k2−1.(iv)
P4 is nonhyperbolic if(48)h=2a11−k1,or(49)h=2a2k2−1.


### 4.5. Dynamical Behavior aboutP5

From (25), eigenvalues of J|P5 about P5 are(50)λ1=1+ha11−k1,λ2=1−ha21−k2,λ3=1+h−a3+α321−k2a2α22.The dynamical behavior aboutP5 is concluded as follows.Lemma 6.
(i)
P5 is a sink if(51)h>max2a1k1−1,2α22a3α22−a2α321−k2,0<h<2a21−k2.(ii)
P5 is a source if(52)0<h<min2a1k1−1,2α22a3α22−a2α321−k2,h>2a21−k2.(iii)
P5 is a saddle if(53)h>max2a1k1−1,2α22a3α22−a2α321−k2,h>2a21−k2.(iv)
P5 is nonhyperbolic if(54)h=2a1k1−1,or(55)h=a3α22−2α22a2α321−k2,or(56)h=2a21−k2.


### 4.6. Dynamical Behavior aboutP6

From (25), eigenvalues of J|P6 about P6 are(57)λ1=1−ha11−k1,λ2=1+ha21−k2−a1α211−k1α11,λ3=1−ha3.The dynamical behavior aboutP6 is concluded as follows.Lemma 7.
(i)
P6 is a sink if(58)0<h<min2a11−k1,2a3h>2α11a1α211−k1−a2α111−k2.(ii)
P6 is a source if(59)h>max2a11−k1,2a3,0<h<2α11a1α211−k1−a2α111−k2.(iii)
P6 is a saddle if(60)0<h<min2a11−k1,2a3,0<h<2α11a1α211−k1−a2α111−k2.(iv)
P6 is nonhyperbolic if(61)h=2a11−k1,or(62)h=2a3,or(63)h=2α11a1α211−k1−a2α111−k2.


### 4.7. Dynamical Behavior aboutP7

From (25), eigenvalues of J|P7 about P7 are(64)λ1=1−ha11−k1,λ2=1−ha21−k2−2a1α211−k1α11,λ3=1+ha3−a21−k2α32α22+a1α21α321−k1α11α22.The dynamical behavior aboutP7 is concluded as follows:Lemma 8.
(i)
P7 is a sink if(65)0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,h>2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(ii)
P7 is a source if(66)h>max2a11−k1,2α11a2α111−k2−2a1α211−k1,0<h<2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iii)
P7 is a saddle if(67)0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,0<h<2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iv)
P7 is nonhyperbolic if(68)h=2a11−k1,or(69)h=2α11a2α111−k2−2a1α211−k1,or(70)h=2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.


### 4.8. Dynamical Behavior aboutP8

From (25), eigenvalues of J|P8 about P8 are(71)λ1=1−ha11−k1,λ2=1−ha21−k2−2a1α211−k1α11,λ3=1−ha3−a21−k2α32α22+a1α21α321−k1α11α22.The dynamical behavior aboutP8 is concluded as follows.Lemma 9.
(i)
P8 is a sink if(72)!0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(ii)
P8 is a source if(73)!h>max2a11−k1,2α11a2α111−k2−2a1α211−k1,2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iii)
P8 is a saddle if(74)0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,h>2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iv)
P8 is nonhyperbolic if(75)h=2a11−k1,or(76)h=2α11a2α111−k2−2a1α211−k1,or(77)h=2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.


## 4.1. Dynamical Behavior aboutP1

From (25), eigenvalues of J|P1 about P1 are(26)λ1=1+ha11−k1,λ2=1+ha21−k2,λ3=1+ha3.The dynamical behavior aboutP1 of model (15) is concluded as follows.Lemma 2.
(i)
For all allowed parametric values,h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33>0, P1 is not sink.(ii)
P1 is a source if(27)0<h<min2a1k1−1,2a2k2−1.(iii)
P1 is a saddle if(28)h>max2a1k1−1,2a2k2−1.(iv)
P1 is nonhyperbolic if(29)h=2a1k1−1,or(30)h=2a2k2−1.


## 4.2. Dynamical Behavior aboutP2

From (25), eigenvalues of J|P2 about P2 are(31)λ1=1+ha11−k1,λ2=1+ha21−k2,λ3=1−ha3.The dynamical behavior aboutP2 is concluded as follows:Lemma 3.
(i)
P2 is a sink if(32)h>max2a1k1−1,2a2k2−1,0<h<2a3.(ii)
P2 is a source if(33)0<h<min2a1k1−1,2a2k2−1,h>2a3.(iii)
P2 is a saddle if(34)h<min2a1k1−1,2a2k2−1,0<h<2a3.(iv)
P2 is nonhyperbolic if(35)h=2a3,or(36)h=2a1k1−1,or(37)h=2a2k2−1.


## 4.3. Dynamical Behavior aboutP3

From (25), eigenvalues of J|P3 about P3 are(38)λ1=1+ha11−k1,λ2=1−ha21−k2,λ3=1+ha3α22−a2α321−k2α22.The dynamical behavior aboutP3 is concluded as follows.Lemma 4.
(i)
P3 is a sink if(39)h>max2a1k1−1,2α22a2α321−k2−a3α22,0<h<2a21−k2.(ii)
P3 is a source if(40)0<h<min2a1k1−1,2α22a2α321−k2−a3α22,h>2a21−k2.(iii)
P3 is a saddle if(41)h>max2a1k1−1,2a21−k2,0<h<2α22a2α321−k2−a3α22.(iv)
P3 is nonhyperbolic if(42)h=2a21−k2,or(43)h=2a1k1−1,or(44)h=2α22a2α321−k2−a3α22.


## 4.4. Dynamical Behavior aboutP4

From (25), eigenvalues of J|P4 about P4 are(45)λ1=1−ha11−k1,λ2=1+ha21−k2,λ3=1+ha3.The dynamical behavior aboutP4 is concluded as follows.Lemma 5.
(i)
For all allowed parametric values,h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33>0, P4 is not sink.(ii)
P4 is a source if(46)h>max2a11−k1,2a2k2−1.(iii)
P4 is a saddle if(47)0<h<min2a11−k1,2a2k2−1.(iv)
P4 is nonhyperbolic if(48)h=2a11−k1,or(49)h=2a2k2−1.


## 4.5. Dynamical Behavior aboutP5

From (25), eigenvalues of J|P5 about P5 are(50)λ1=1+ha11−k1,λ2=1−ha21−k2,λ3=1+h−a3+α321−k2a2α22.The dynamical behavior aboutP5 is concluded as follows.Lemma 6.
(i)
P5 is a sink if(51)h>max2a1k1−1,2α22a3α22−a2α321−k2,0<h<2a21−k2.(ii)
P5 is a source if(52)0<h<min2a1k1−1,2α22a3α22−a2α321−k2,h>2a21−k2.(iii)
P5 is a saddle if(53)h>max2a1k1−1,2α22a3α22−a2α321−k2,h>2a21−k2.(iv)
P5 is nonhyperbolic if(54)h=2a1k1−1,or(55)h=a3α22−2α22a2α321−k2,or(56)h=2a21−k2.


## 4.6. Dynamical Behavior aboutP6

From (25), eigenvalues of J|P6 about P6 are(57)λ1=1−ha11−k1,λ2=1+ha21−k2−a1α211−k1α11,λ3=1−ha3.The dynamical behavior aboutP6 is concluded as follows.Lemma 7.
(i)
P6 is a sink if(58)0<h<min2a11−k1,2a3h>2α11a1α211−k1−a2α111−k2.(ii)
P6 is a source if(59)h>max2a11−k1,2a3,0<h<2α11a1α211−k1−a2α111−k2.(iii)
P6 is a saddle if(60)0<h<min2a11−k1,2a3,0<h<2α11a1α211−k1−a2α111−k2.(iv)
P6 is nonhyperbolic if(61)h=2a11−k1,or(62)h=2a3,or(63)h=2α11a1α211−k1−a2α111−k2.


## 4.7. Dynamical Behavior aboutP7

From (25), eigenvalues of J|P7 about P7 are(64)λ1=1−ha11−k1,λ2=1−ha21−k2−2a1α211−k1α11,λ3=1+ha3−a21−k2α32α22+a1α21α321−k1α11α22.The dynamical behavior aboutP7 is concluded as follows:Lemma 8.
(i)
P7 is a sink if(65)0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,h>2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(ii)
P7 is a source if(66)h>max2a11−k1,2α11a2α111−k2−2a1α211−k1,0<h<2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iii)
P7 is a saddle if(67)0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,0<h<2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iv)
P7 is nonhyperbolic if(68)h=2a11−k1,or(69)h=2α11a2α111−k2−2a1α211−k1,or(70)h=2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.


## 4.8. Dynamical Behavior aboutP8

From (25), eigenvalues of J|P8 about P8 are(71)λ1=1−ha11−k1,λ2=1−ha21−k2−2a1α211−k1α11,λ3=1−ha3−a21−k2α32α22+a1α21α321−k1α11α22.The dynamical behavior aboutP8 is concluded as follows.Lemma 9.
(i)
P8 is a sink if(72)!0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(ii)
P8 is a source if(73)!h>max2a11−k1,2α11a2α111−k2−2a1α211−k1,2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iii)
P8 is a saddle if(74)0<h<min2a11−k1,2α11a2α111−k2−2a1α211−k1,h>2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.(iv)
P8 is nonhyperbolic if(75)h=2a11−k1,or(76)h=2α11a2α111−k2−2a1α211−k1,or(77)h=2α11α22a3α11α22−a2α11α321−k2+a1α21α321−k1.


## 5. Periodic Points

We will prove thatPii=1,…,8 of model (15) are periodic points of period n.Theorem 1.
Equilibrium pointsPii=1,…,8 of model (15) are periodic points of prime period 1.Proof.
From (15), define(78)Fx,y,z≔f1,f2,f3,where f1,f2, and f3 are represented in (23). From (78), the computation yields(79)F|P1=0,0,0=P1,F|P2=0,0,a3/α33=P2,F|P3=0,1−k2a2/α22,0=P3,F|P4=1−k1a1/α11,0,0=P4,F|P5=0,1−k2a2/α22,a3α22−α321−k2a2/α22α33=P5,F|P6=1−k1a1/α11,0,a3/α33=P6,F|P7=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,0=P7,F|P8=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,a3α11α22−a21−k2α11α32+a11−k1α21α32/α11α22α33=P8.
Hence, from (79), we can say that equilibrium points Pii=1,…,8 of three species model (15) are periodic points of prime period 1.Now, it is proved that equilibrium pointsPii=1,…,8 are period points of period n.Theorem 2.
P1 of model (15) is a periodic point of period n.Proof.
From (78), the following computation yields the required statement:(80)F2=f1+ha11−k1f1−α11f12,f2+ha21−k2f2−α22f22−α21f1f2,f3+ha3f3−α33f22−α32f2f3⇒F2|P1=P1,F3=f12+ha11−k1f12−α11f122,f22+ha21−k2f22−α22f222−α21f12f22x,y,z,f32+ha3f32−α33f222−α32f22f32⇒F3|P1=P1,⋮Fn=f1n+ha11−k1f1n−α11f1n2,f2n+ha21−k2f2n−α22f2n2−α21f1nf2n,f3n+ha3f3n−α33f2n2−α32f2nf3n⇒Fn|P1=P1.Theorem 3.
P2 of model (15) is a periodic point of period n.Proof.
Utilizing the computation as we have done in (80), one gets the following required statement:(81)F2|P2=0,0,a3/α33=P2,F3|P2=0,0,a3/α33=P2,⋮Fn|P2=0,0,a3/α33=P2.Theorem 4.
P3 of model (15) is a periodic point of period n.Proof.
In view of (80), one gets the following required statement:(82)F2|P3=0,1−k2a2/α22,0=P3,F3|P3=0,1−k2a2/α22,0=P3,⋮Fn|P3=0,1−k2a2/α22,0=P3.Theorem 5.
P4 of model (15) is a periodic point of period n.Proof.
In view of (80), one gets the following required statement:(83)F2|P4=1−k1a1/α11,0,0=P4,F3|P4=1−k1a1/α11,0,0=P4,⋮Fn|P4=1−k1a1/α11,0,0=P4.Theorem 6.
P5 of model (15) is a periodic point of period n.Proof.
From (80), one obtains(84)F2|P5=0,1−k2a2/α22,a3α22−α321−k2a2/α22α33=P5,F3|P5=0,1−k2a2/α22,a3α22−α321−k2a2/α22α33=P5,⋮Fn|P5=0,1−k2a2/α22,a3α22−α321−k2a2/α22α33=P5.Theorem 7.
P6 of model (15) is a periodic point of period n.Proof.
From (80), one obtains(85)F2|P6=1−k1a1/α11,0,a3/α33=P6,F3|P6=1−k1a1/α11,0,a3/α33=P6,⋮Fn|P6=1−k1a1/α11,0,a3/α33=P6.Theorem 8.
P7 of model (15) is a periodic point of period n.Proof.
From (80), one obtains(86)F2|P7=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,0=P7,F3|P7=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,0=P7,⋮Fn|P7=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,0=P7.Theorem 9.
P8 of model (15) is a periodic point of period n.Proof.
From (80), one obtains(87)F2|P8=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,a3α11α22−a21−k2α11α32+a11−k1α21α32/α11α22α33=P8,F3|P8=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,a3α11α22−a21−k2α11α32+a11−k1α21α32/α11α22α33=P8,⋮Fn|P8=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,a3α11α22−a21−k2α11α32+a11−k1α21α32/α11α22α33=P8.


## 6. Analysis of Bifurcation

In this section, we give analysis of bifurcation about fixed pointsPii=1,…,8 of model (15) by bifurcation theory [18, 19].


### 6.1. Analysis of Bifurcation atP1

Here, we will study analysis of bifurcation atP1 of model (15). From (26), the simple computation yields λ1|29=−1, but λ2,3|29=1−2a21−k2/a11−k1,1−2a3/a11−k1≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P1 if Ω=h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33 passes the curve:(88)ℱ|P1=Ω:h=2a1k1−1.However, flip bifurcation cannot occur by computation, soP1 is degenerated with high co-dimension as Ω∈ℱ|P1.


### 6.2. Analysis of Bifurcation atP2

We will study analysis of bifurcation atP2 of model (15). From (26), the simple computation yields λ3|35=−1, but λ1,2|35=1+2a11−k1/a3,1+2a21−k2/a3≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P2 if Ω passes the curve:(89)ℱ|P2=Ω:h=2a3.The proof of following theorem shows that model (15) undergoes flip bifurcation around P2 if Ω∈ℱ|P2.Theorem 10.
Model (15) undergo flip bifurcation around P2 if Ω∈ℱ|P2.Proof.
It is noticed that three-species model (15) is invariant with respect to x=y=0. Thus, we restrict (15) on x=y=0, to determine the bifurcation, where it takes the form(90)zn+1=zn+ha3zn−α33zn2.
From (90), define(91)fz≔z+ha3z−α33z2.
Now, one denotesh=h∗=2/a3andz=z∗=a3/α33. The computation yields(92)fz|h=h∗=2/a3,z=z∗=a3/α33=−1,(93)fzz|h=h∗=2/a3,z=z∗=a3/α33=−4α33a3≠0,(94)fh|h=h∗=2/a3,z=z∗=a3/α33=−a32α33≠0.
From (92)–(94), it can be concluded that the model undergoes flip bifurcation around P2 if Ω∈ℱ|P2.


### 6.3. Analysis of Bifurcation atP3

From (38), the computation yields λ2|42=−1, but λ1,3|42=1+2a11−k1/a21−k2,1+2/a21−k2a3α22−a2α321−k2≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P3 if Ω passes the curve:(95)ℱ|P3=Ω:h=2a21−k2.The proof of following theorem shows that model (15) undergoes flip bifurcation around P3 if Ω∈ℱ|P3.Theorem 11.
Model (15) undergoes flip bifurcation around P3 if Ω∈ℱ|P3.Proof.
It is noticed that,w.r.tx=z=0, model (15) is invariant. So, one restricts model (15) on x=z=0, where it becomes(96)yn+1=yn+ha21−k2yn−α22yn2.
From (96), define(97)fy≔y+ha21−k2y−hα22y2.
Denoteh=h∗=2/a21−k2andy=y∗=a21−k2/α22. By computation, one obtains(98)fy|h=h∗=2/a21−k2,y=y∗=a21−k2/α22=−1,(99)fyy|h=h∗=2/a21−k2,y=y∗=a21−k2/α22=−4α22a21−k2≠0,(100)fh|h=h∗=2/a21−k2,y=y∗=a21−k2/α22=a21−k2α22≠0.
So, model (15) undergoes flip bifurcation by (98)–(100) if Ω∈ℱ|P3.


### 6.4. Analysis of Bifurcation atP4

From (45), the computation yields λ1|48=−1, but λ2,3|48=1+2a21−k2/a11−k1,1+2a3/a11−k1≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P4 if Ω passes the curve:(101)ℱ|P3=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P4 if Ω∈ℱ|P4.Theorem 12.
Model (15) undergoes flip bifurcation around P4 if Ω∈ℱ|P4.Proof.
It is noticed that,w.r.ty=z=0, model (15) is invariant. So, one restricts model (15) on y=z=0, where it becomes(102)xn+1=xn+ha11−k1xn−α11xn2.
From (102), define(103)fx≔x+ha11−k1x−hα11x2.
Denoteh=h∗=2/a11−k1,x=x∗=a11−k1/α11. By computation, one obtains(104)fx|h=h∗=2/a11−k1,x=x∗=a11−k1/α11=−1,(105)fxx|h=h∗=2/a11−k1,x=x∗=a11−k1/α11=−4α11a11−k1≠0,(106)fh|h=h∗=2/a11−k1,x=x∗=a11−k1/α11=a11−k1α11≠0.
So, model (15) undergoes flip bifurcation by (104)–(106) if Ω∈ℱ|P4.


### 6.5. Analysis of Bifurcation atP5

From (50), the computation yields λ1|54=−1, but λ2,3|54=1+2a21−k2/a11−k1,1−2/a11−k1−a3+α321−k2a2/α22≠1or−1. This suggests that model (15) could undergo flip bifurcation around P5 if Ω passes the curve:(107)ℱ|P5=Ω:h=2a1k1−1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P5 if Ω∈ℱ|P5.Theorem 13.
Model (15) undergoes flip bifurcation around P5 if Ω∈ℱ|P5.Proof.
Recall that ifΩ∈ℱ|P5, then λ1|54=−1, but λ2,3|54=1+2a21−k2/a11−k1,1−2/a11−k1−a3+α321−k2a2/α22≠1or−1. So, hereafter, detailed flip bifurcation is explored if Ω varies in the nbhd of h, i.e.,h=h+ε, by assuming h≠2α22/a3α22−a2α321−k2,2/a21−k2. Let(108)un=xn,vn=yn−1−k2a2α22,wn=zn−a3α22−α321−k2a2α22α33.
Then, (15) gives(109)un+1vn+1wn+1=1+ha11−k100−hα211−k2α22a21−ha21−k200−hα32a3α32−α321−k2a2α22α331+h−a3+α321−k2a2α22unvnwn+hα11un2+a11−k1εun−α11εun2ha21−k22α22−α22vn+1−k2a2α222+α21unvn+εa21−k2vn+a21−k2α22−α22vn+1−k2a2α222−α21unvn+a21−k2α22ha3a3α32−α321−k2a2α22α33−α33wn+a3α22−α321−k2a2α22α332−α32vnwn−1−k2a2α22a3α32−α321−k2a2α22α33+εa3wn+a3α32−α321−k2a2α22α33−α33wn+a3α22−α321−k2a2α22α332−α32vn+1−k2a2α22wn+a3α22−α321−k2a2α22α33.
By using transformation,(110)unvnwn=a1100a21a220111xnynzn.
(109) takes the form(111)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+Fxn,yn,zn,εGxn,yn,zn,εHxn,yn,zn,ε,where(112)a11=−a1−a2+a1k1+a2k2α22a1α22+a3α22−a1k1α22−a2α32+a2k2α32α33a2−1+k2−a2+a3+a2k2α21α322,a21=−a1α22+a3α22−a1k1α22−a2α32+a2k2α32α33−a2+a3+a2k2α322,a22=−−a2α22+a3α22+a2k2α22−a2α32+a2k2α32α33−a2+a3+a2k2α322,F=xna11−k1ε−xn2a11α11+hxn2a11α11−xn2a11α11ε,G=−a1a111−k1ε−xn2a112α11+hxn2a112α11−xn2εa112α11a11+ha221−k22/α22+xna11α21xna21+ynα22−α22xna21+ynα221−k2/α222a21+εa21−k2xna21+ynα22+a21−k2/α22−xnεa11α21xna21+ynα22+a21−k2/α22a21−εα22xna21+ynα22+a21−k2/α112a21,H=h−xnxn+yn+zna11α32−a21−k2a3α32−a2α321−k2α222α33+a3a3α32−a2α321−k2α22α33−xn+yn+zn+a3α32−a2α321−k2α22α332α33−ha221−k22/α22+xna11α21xna21+ynα22−α22xna21+ynα221−k2/α222a21+εa21−k2xna21+ynα22+a21−k2/α22−xnεa11α21xna21+ynα22+a21−k2/α22a21−εα22xna21+ynα22+a21−k2/α112a21εa3xn+yn+zn+a3α32−a2α321−k2α22α33−εxna21+ynα22+a21−k2α22α32xn+yn+zn+a3α32−a21−k2α32α22α33−xn+yn+zn+a3α32−a2α321−k2α22α332α33.
Now, consider (111) on the center manifold, i.e.,(113)Wc0=xn,yn,zn|yn,zn=χ1xn,χ2xn,χi0=0,Dχi0=0,i=1,2,where(114)χixn=aix2+bix3+Ox4,fori=1,2.
From (111) and (113), one obtains(115)χ1−xn+Fxn,χ1,χ2=λ2χ1xn+Gxn,χ1,χ2,χ2−xn+Fxn,χ1,χ2=λ3χ2xn+Hxn,χ1,χ2.
From (115), computation yields a1=b1=a2=b2=0. Finally, map (111); restrict to Wc0 as(116)fxn=−xn+xna11−k1ε−xn2a11α11+hxn2a11α11−xn2a11α11ε+Oxn+ε3.
For the model to undergo flip bifurcation, the following should be nonzero:(117)Ω1=∂2f∂xn∂ε+12∂f∂ε∂2f∂xn2|O=a11−k1≠0,Ω2=16∂3f∂xn3+12∂2f∂xn22|O=a11α11−ha11α112>0.
From (117), one can say that about P5 model (15) undergoes flip bifurcation if Ω∈ℱ|P5. Moreover, period-2 points bifurcating from P5 are stable since Ω2=a11α11−ha11α112>0.


### 6.6. Analysis of Bifurcation atP6

From (57), the computation yields λ1|61=−1, but λ2,3|61=1+2/a11−k1a21−k2−a1α211−k1/α11,1−2a3/a11−k1≠1or−1. This suggests that model (15) could undergo flip bifurcation around P6 if Ω passes the curve:(118)ℱ|P6=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P6 if Ω∈ℱ|P6.Theorem 14.
Model (15) undergoes flip bifurcation around P6 if Ω∈ℱ|P6.Proof.
Recall that ifΩ∈ℱ|P6, then λ1|61=−1, but λ2,3|61=1+2/a11−k1a21−k2−a1α211−k1/α11,1−2a3/a11−k1≠1or−1. So, hereafter, detailed flip bifurcation is explored if Ω varies in the nbhd of h, i.e.,h=h+ε, by assuming h≠2/a3,2α11/a1α211−k1−a2α111−k2. Let(119)un=xn−1−k1a1α11,vn=yn,wn=zn−a3α33.
Then, (15) gives(120)un+1vn+1vn+1=1−ha11−k10001+ha21−k2−a1α211−k1α1100−hα32a3α321−ha3unvnwn+ha11−k12α11−hα11un+a11−k1α112+εa11−k1un+a11−k1α11−εα11un+a11−k1α112−hα22vn2−hα22unvn+εa21−k2vn−εα22vn2−εα21vnun+a11−k1α11ha32α33−α33wn+a3α332−α32vnwn+εa3wn+a3α33−εα33wn+a3α332−εα32vnwn+a3α33.
Using transformation,(121)unvnwn=01000b23101xnynzn.
(120) becomes(122)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+F1xn,yn,zn,εG1xn,yn,zn,εH1xn,yn,zn,ε,where(123)b23=α33−a2α11−a3α11+a2k2α11+a1α21−a1k1α21a3α11α32,F1=−zna21−k2ε+znyn+a11−k1α11α21ε+hznα22+hynznα22+zn2b23α22ε+a3xn+yn+a3α33ε−znb23α32xn+yn+a3α33ε−xn+yn+a3α332α33ε+h−xn+ynznb23α32+a32α33−xn+yn+a3α332α33,G1=a11−k1yn+a11−k1α11ε−hyn+a11−k1α112+ha121−k12α11−yn+a11−k1α112α11ε,H1=−zna21−k2ε+znyn+a11−k1α11α21ε+hznα22+hynznα22+zn2b23α22ε.
Now, from model (122) on the center manifold,(124)Wc0=xn,yn,zn|yn,zn=χ3xn,χ4xn,χi0=0,Dχi0=0,i=3,4,where(125)χixn=aix2+bix3+Ox4,fori=3,4.
From (122) and (124), one has(126)χ3−xn+F1xn,χ3,χ4=λ2χ3xn+G1xn,χ3,χ4,χ4−xn+F1xn,χ3,χ4=λ3χ4xn+H1xn,χ3,χ4.
From (126), the calculation yields: a3=b3=a4=b4=0. Thus, map (122); restrict to Wc0 as(127)fxn=−xn+a3xn+a3α33ε−xn+a3α332α33ε+ha32α33−xn+a3α332α33.
From (117) and (127), the computation yields: Ω1=3a3−ha32/α33≠0 and Ω2=h2α332>0. This implies that about P6 model (15) undergoes flip bifurcation if Ω∈ℱ|P6. Moreover, period-2 points bifurcating from P6 are stable since Ω2=h2α332>0.


### 6.7. Analysis of Bifurcation atP7

From (64), the computation yields λ1|68=−1, but λ2,3|68=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1+2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. This suggests that model (15) could undergo flip bifurcation around P7 if Ω passes the curve:(128)ℱ|P7=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P7 if Ω∈ℱ|P7.Theorem 15.
Model (15) undergoes flip bifurcation around P7 if Ω∈ℱ|P7.Proof.
Recall that ifΩ∈ℱ|P7, then λ1|68=−1, but λ2,3|68=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1+2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. So, in the following, flip bifurcation is explored by assuming h≠2α11/a2α111−k2−2a1α211−k1,2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1. Let(129)un=xn−1−k1a1α11,vn=yn−a21−k2α11−a11−k1α21α11α22,wn=zn.
Then, (15) gives(130)un+1vn+1vn+1=1−ha11−k100−hα21a21−k2α11−a11−k1α21α11α221−ha21−k2−2a1α211−k1α110001+ha3−a21−k2α32α22+a1α21α321−k1α11α22unvnwn+ha11−k12α11−α11un+a11−k1α112+εa11−k1un+a11−k1α11−α11un+a11−k1α112ha21−k2a21−k2α11−a11−k1α21α11α22+α22vn+a21−k2α11−a11−k1α21α11α222−α21a11−k1α11a21−k2α11−a11−k1α21α11α22+εa21−k2vn+a21−k2α11−a11−k1α21α11α22−α22vn+a21−k2α11−a11−k1α21α11α222−α21un+a11−k1α11vn+a21−k2α11−a11−k1α21α11α22h−α33wn2−α32vnwn+εa3wn−α33wn2−α32wnvn+a21−k2α11−a11−k1α21α11α22.
Using transformation,(131)unvnwn=c1100110001xnynzn.
(130) gives(132)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+F2xn,yn,zn,εG2xn,yn,zn,εH2xn,yn,zn,ε,where(133)c11=α11α22ha11−k1−ha21−k2−2ha11−k1α21/α11hα21a21−k2α11−a11−k1α21,F2=ha11−k12/α11−xnc11+a11−k1/α112α11+εa11−k1xnc11+a11−k1/α11−xnc11+a11−k1/α112α11c11,G2=ha21−k2a2a111−k2−a11−k1α21α11α22−a11−k1α21a2a111−k2−a11−k1α21α112α22+xn+yn+a2a111−k2−a11−k1α21α11α222α22ha11−k12/α11−xnc11+a11−k1/α112α11+εa11−k1xnc11+a11−k1/α11−xnc11+a11−k1/α112α11c11+εa21−k2xn+yn+a2a111−k2−a11−k1α21α11α22+a11−k1α21xn+yn+a2a111−k2−a11−k1α21/α11α223xnc11a2a111−k2−a11−k1α21/α11α22α22a11,H2=h−xn+ynznα32−z2α33+εzna3−znxn+yn+a2a111−k2−a11−k1α21α11α22α32−zn2α33.
Now, using system (132) on the center manifold,(134)Wc0=xn,yn,zn|yn,zn=χ5xn,χ6xn,χi0=0,Dχi0=0,i=5,6,where(135)χixn=aix2+bix3+Ox4,fori=5,6.
In view of (132) and (134), we obtain(136)χ5−xn+F2xn,χ5,χ6=λ2χ5xn+G2xn,χ5,χ6,χ6−xn+F2xn,χ5,χ6=λ3χ6xn+H2xn,χ5,χ6.
From (136), one gets: a5=b5=a6=b6=0. Finally, map (132); restrict to Wc0 as(137)fxn=−xn+ha11−k12/α11−xnc11+a11−k1/α112α11+εa11−k1xnc11+a11−k1/α11−xnc11+a11−k1/α112α11c11.
From (117) and (137), the computation yields: Ω1=a11−k11−2c11≠0. Moreover, Ω2=hα11c11+2α11c11ε2>0. This implies that about P7 model (15) undergoes flip bifurcation if Ω∈ℱ|P7. Moreover, period-2 points bifurcating from P7 are stable since Ω2=hα11c11+2α11c11ε2>0.


### 6.8. Analysis of Bifurcation atP8

From (71), the computation yields λ1|75=−1, but λ2,3|75=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1−2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. This suggests that model (15) could undergo flip bifurcation around P8 if Ω passes the curve:(138)ℱ|P8=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P8 if Ω∈ℱ|P8.Theorem 16.
Model (15) undergoes flip bifurcation around P8 if Ω∈ℱ|P8.Proof.
Recall that ifΩ∈ℱ|P8, then λ1|75=−1, but λ2,3|75=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1−2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. So, in the following, flip bifurcation is explored by assuming h≠2α11/a2α111−k2−2a1α211−k1,2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1. Let(139)un=xn−1−k1a1α11,vn=yn−a21−k2α11−a11−k1α21α11α22,wn=zn−a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33.
Then, (15) becomes(140)un+1vn+1vn+1=1−ha11−k100−hα21a21−k2α11−a11−k1α21α11α221−ha21−k2−2a1α211−k1α1100−hα32a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α331−ha3−a21−k2α32α22+a1α21α321−k1α11α22unvnwn+F3xn,yn,zn,εG3xn,yn,zn,εH3xn,yn,zn,ε,where(141)F3=ha11−k12α11−α11un+a11−k1α112+εa11−k1un+a11−k1α11−α11un+a11−k1α112,G3=ha21−k2a21−k2α11−a11−k1α21α11α22−α21vn+a21−k2α11−a11−k1α21α11α222−α21unvn−α21a11−k1α11a21−k2α11−a11−k1α21α11α22+εa21−k2vn+a21−k2α11−a11−k1α21α11α22−α22vn+a21−k2α11−a11−k1α21α11α222−α21un+a11−k1α11vn+a21−k2α11−a11−k1α21α11α22,H3=ha3a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−α32vnwn−α33wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332−α32a21−k2α11−a11−k1α21α11α22a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33+εa3wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−α33wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332−α32vn+a21−k2α11−a11−k1α21α11α22ε×wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33.
Now, by utilizing transformation,(142)unvnwn=d1100d21d220111xnynzn,gives(143)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+F2∗xn,yn,zn,εG2∗xn,yn,zn,εH2∗xn,yn,zn,ε,where(144)d11=α112ha1−ha1k1−ha21−k2−2ha11−k1α21/α11α222α33ha1−ha1k1−ha3−a21−k2α32/α22+a11−k1α21α32/α11α22h2α21a21−k2α11−a11−k1α21α32a3α11α22−a21−k2α11α32+a11−k1α21α32,d21=α11α22ha1−ha1k1−ha3−a21−k2α32/α22+a11−k1α21α32/α11α22α33hα32a3α11α22−a21−k2α11α32+a11−k1α21α32,d22=−α33−a2α11α22+a3α11α22+a2k2α11α22−2a1α21α22+2a1k1α21α22−a2α11α32+a2k2α11α32+a1α21α32−a1k1α21α32α32a3α11α22−a2α11α32+a2k2α11α32+a1α21α32−a1k1α21α32,F2∗=ha11−k12/α11−xnd11+a11−k1/α112α11+εa11−k1xnd11+a11−k1/α11−xnd11+a11−k1/α112α11d11,G2∗=−d21ha11−k12/α11−xnd11+a11−k1/α112α11+εa11−k1xnd11+a11−k1/α11−xnd11+a11−k1/α112α11d11d22+hd22−xnd11d21xn+d22ynα21−α21d21xn+d22yn+a2a111−k2−a11−k1α21α11α222+a21−k2a2a111−k2−a11−k1α21α11α22−a11−k1α21a2a111−k2−a11−k1α21α112α22+εd22a21−k2d21xn+d22yn+a2a111−k2−a11−k1α21α11α22−a11−k1α21d11x+a11−k1/α11d21xn+d22yn+a2a111−k2−a11−k1α21/α11α22α11−d21xn+d22yn+a2a111−k2−a11−k1α21α11α222α22,H2∗=h−yn+2znd21xn+d22ynα32−α32xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332+a3a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−a2a111−k2−a11−k1α21α32a3α11α22−a21−k2α11α32+a11−k1α21α32α112α222α33+d21−d22ha11−k12/α11−d11xn+a11−k1/α112α11+ϵa11−k1d11xn+a11−k1/α11−d11xn+a11−k1/α112α11d11d22−hd22−d11xnd21xn+d22ynα21−α21d21xn+d22yn+a2a111−k2−a11−k1α21α11α222+a21−k2a2a111−k2−a11−k1α21α11α22−a11−k1α21a2a111−k2−a11−k1α21α112α22+εd22a21−k2d21xn+d22yn+a2a111−k2−a11−k1α21α11α22−a11−k1α21d11xn+a11−k1/α11d21xn+d22yn+a2a111−k2−a11−k1α21/α11α22α11−d21xn+d22yn+a2a111−k2−a11−k1α21α11α222α22+εa3xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−α32d21xn+d22yn+a2a111−k2−a11−k1α21α11α22×xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332α33.
Now, from system (143) on the center manifold(145)Wc0=xn,yn,zn|yn,zn=χ7xn,χ8xn,χi0=0,Dχi0=0,i=7,8,where(146)χixn=aix2+bix3+Ox4,fori=7,8.
From (143) and (145), we obtain(147)χ7−xn+F2∗xn,χ7,χ8=λ2χ7xn+G2∗xn,χ7,χ8,χ8−xn+F2∗xn,χ7,χ8=λ3χ8xn+H2∗xn,χ7,χ8.
From (147), computation yields: a7=b7=a8=b8=0. Finally, map (143); restrict to Wc0 as(148)f1∗xn=−xn+ha11−k12/α11−xnd11+a11−k1/α112α11+εa11−k1xnd11+a11−k1/α11−xnd11+a11−k1/α112α11d11.
From (117) and (148), the computation yields: Ω1=a11−k11−2d11≠0. Moreover, Ω2=hα11d11+2α11d11ε2>0. This implies that about P8 model (15) undergoes flip bifurcation if Ω∈ℱ|P8. Moreover, period-2 points bifurcating from P8 are stable since Ω2=hα11d11+2α11d112>0.


## 6.1. Analysis of Bifurcation atP1

Here, we will study analysis of bifurcation atP1 of model (15). From (26), the simple computation yields λ1|29=−1, but λ2,3|29=1−2a21−k2/a11−k1,1−2a3/a11−k1≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P1 if Ω=h,a1,a2,a3,k1,k2,α11,α22,α21,α32,α33 passes the curve:(88)ℱ|P1=Ω:h=2a1k1−1.However, flip bifurcation cannot occur by computation, soP1 is degenerated with high co-dimension as Ω∈ℱ|P1.


## 6.2. Analysis of Bifurcation atP2

We will study analysis of bifurcation atP2 of model (15). From (26), the simple computation yields λ3|35=−1, but λ1,2|35=1+2a11−k1/a3,1+2a21−k2/a3≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P2 if Ω passes the curve:(89)ℱ|P2=Ω:h=2a3.The proof of following theorem shows that model (15) undergoes flip bifurcation around P2 if Ω∈ℱ|P2.Theorem 10.
Model (15) undergo flip bifurcation around P2 if Ω∈ℱ|P2.Proof.
It is noticed that three-species model (15) is invariant with respect to x=y=0. Thus, we restrict (15) on x=y=0, to determine the bifurcation, where it takes the form(90)zn+1=zn+ha3zn−α33zn2.
From (90), define(91)fz≔z+ha3z−α33z2.
Now, one denotesh=h∗=2/a3andz=z∗=a3/α33. The computation yields(92)fz|h=h∗=2/a3,z=z∗=a3/α33=−1,(93)fzz|h=h∗=2/a3,z=z∗=a3/α33=−4α33a3≠0,(94)fh|h=h∗=2/a3,z=z∗=a3/α33=−a32α33≠0.
From (92)–(94), it can be concluded that the model undergoes flip bifurcation around P2 if Ω∈ℱ|P2.


## 6.3. Analysis of Bifurcation atP3

From (38), the computation yields λ2|42=−1, but λ1,3|42=1+2a11−k1/a21−k2,1+2/a21−k2a3α22−a2α321−k2≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P3 if Ω passes the curve:(95)ℱ|P3=Ω:h=2a21−k2.The proof of following theorem shows that model (15) undergoes flip bifurcation around P3 if Ω∈ℱ|P3.Theorem 11.
Model (15) undergoes flip bifurcation around P3 if Ω∈ℱ|P3.Proof.
It is noticed that,w.r.tx=z=0, model (15) is invariant. So, one restricts model (15) on x=z=0, where it becomes(96)yn+1=yn+ha21−k2yn−α22yn2.
From (96), define(97)fy≔y+ha21−k2y−hα22y2.
Denoteh=h∗=2/a21−k2andy=y∗=a21−k2/α22. By computation, one obtains(98)fy|h=h∗=2/a21−k2,y=y∗=a21−k2/α22=−1,(99)fyy|h=h∗=2/a21−k2,y=y∗=a21−k2/α22=−4α22a21−k2≠0,(100)fh|h=h∗=2/a21−k2,y=y∗=a21−k2/α22=a21−k2α22≠0.
So, model (15) undergoes flip bifurcation by (98)–(100) if Ω∈ℱ|P3.


## 6.4. Analysis of Bifurcation atP4

From (45), the computation yields λ1|48=−1, but λ2,3|48=1+2a21−k2/a11−k1,1+2a3/a11−k1≠1or−1. This suggests that model (15) could undergo a flip bifurcation around P4 if Ω passes the curve:(101)ℱ|P3=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P4 if Ω∈ℱ|P4.Theorem 12.
Model (15) undergoes flip bifurcation around P4 if Ω∈ℱ|P4.Proof.
It is noticed that,w.r.ty=z=0, model (15) is invariant. So, one restricts model (15) on y=z=0, where it becomes(102)xn+1=xn+ha11−k1xn−α11xn2.
From (102), define(103)fx≔x+ha11−k1x−hα11x2.
Denoteh=h∗=2/a11−k1,x=x∗=a11−k1/α11. By computation, one obtains(104)fx|h=h∗=2/a11−k1,x=x∗=a11−k1/α11=−1,(105)fxx|h=h∗=2/a11−k1,x=x∗=a11−k1/α11=−4α11a11−k1≠0,(106)fh|h=h∗=2/a11−k1,x=x∗=a11−k1/α11=a11−k1α11≠0.
So, model (15) undergoes flip bifurcation by (104)–(106) if Ω∈ℱ|P4.


## 6.5. Analysis of Bifurcation atP5

From (50), the computation yields λ1|54=−1, but λ2,3|54=1+2a21−k2/a11−k1,1−2/a11−k1−a3+α321−k2a2/α22≠1or−1. This suggests that model (15) could undergo flip bifurcation around P5 if Ω passes the curve:(107)ℱ|P5=Ω:h=2a1k1−1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P5 if Ω∈ℱ|P5.Theorem 13.
Model (15) undergoes flip bifurcation around P5 if Ω∈ℱ|P5.Proof.
Recall that ifΩ∈ℱ|P5, then λ1|54=−1, but λ2,3|54=1+2a21−k2/a11−k1,1−2/a11−k1−a3+α321−k2a2/α22≠1or−1. So, hereafter, detailed flip bifurcation is explored if Ω varies in the nbhd of h, i.e.,h=h+ε, by assuming h≠2α22/a3α22−a2α321−k2,2/a21−k2. Let(108)un=xn,vn=yn−1−k2a2α22,wn=zn−a3α22−α321−k2a2α22α33.
Then, (15) gives(109)un+1vn+1wn+1=1+ha11−k100−hα211−k2α22a21−ha21−k200−hα32a3α32−α321−k2a2α22α331+h−a3+α321−k2a2α22unvnwn+hα11un2+a11−k1εun−α11εun2ha21−k22α22−α22vn+1−k2a2α222+α21unvn+εa21−k2vn+a21−k2α22−α22vn+1−k2a2α222−α21unvn+a21−k2α22ha3a3α32−α321−k2a2α22α33−α33wn+a3α22−α321−k2a2α22α332−α32vnwn−1−k2a2α22a3α32−α321−k2a2α22α33+εa3wn+a3α32−α321−k2a2α22α33−α33wn+a3α22−α321−k2a2α22α332−α32vn+1−k2a2α22wn+a3α22−α321−k2a2α22α33.
By using transformation,(110)unvnwn=a1100a21a220111xnynzn.
(109) takes the form(111)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+Fxn,yn,zn,εGxn,yn,zn,εHxn,yn,zn,ε,where(112)a11=−a1−a2+a1k1+a2k2α22a1α22+a3α22−a1k1α22−a2α32+a2k2α32α33a2−1+k2−a2+a3+a2k2α21α322,a21=−a1α22+a3α22−a1k1α22−a2α32+a2k2α32α33−a2+a3+a2k2α322,a22=−−a2α22+a3α22+a2k2α22−a2α32+a2k2α32α33−a2+a3+a2k2α322,F=xna11−k1ε−xn2a11α11+hxn2a11α11−xn2a11α11ε,G=−a1a111−k1ε−xn2a112α11+hxn2a112α11−xn2εa112α11a11+ha221−k22/α22+xna11α21xna21+ynα22−α22xna21+ynα221−k2/α222a21+εa21−k2xna21+ynα22+a21−k2/α22−xnεa11α21xna21+ynα22+a21−k2/α22a21−εα22xna21+ynα22+a21−k2/α112a21,H=h−xnxn+yn+zna11α32−a21−k2a3α32−a2α321−k2α222α33+a3a3α32−a2α321−k2α22α33−xn+yn+zn+a3α32−a2α321−k2α22α332α33−ha221−k22/α22+xna11α21xna21+ynα22−α22xna21+ynα221−k2/α222a21+εa21−k2xna21+ynα22+a21−k2/α22−xnεa11α21xna21+ynα22+a21−k2/α22a21−εα22xna21+ynα22+a21−k2/α112a21εa3xn+yn+zn+a3α32−a2α321−k2α22α33−εxna21+ynα22+a21−k2α22α32xn+yn+zn+a3α32−a21−k2α32α22α33−xn+yn+zn+a3α32−a2α321−k2α22α332α33.
Now, consider (111) on the center manifold, i.e.,(113)Wc0=xn,yn,zn|yn,zn=χ1xn,χ2xn,χi0=0,Dχi0=0,i=1,2,where(114)χixn=aix2+bix3+Ox4,fori=1,2.
From (111) and (113), one obtains(115)χ1−xn+Fxn,χ1,χ2=λ2χ1xn+Gxn,χ1,χ2,χ2−xn+Fxn,χ1,χ2=λ3χ2xn+Hxn,χ1,χ2.
From (115), computation yields a1=b1=a2=b2=0. Finally, map (111); restrict to Wc0 as(116)fxn=−xn+xna11−k1ε−xn2a11α11+hxn2a11α11−xn2a11α11ε+Oxn+ε3.
For the model to undergo flip bifurcation, the following should be nonzero:(117)Ω1=∂2f∂xn∂ε+12∂f∂ε∂2f∂xn2|O=a11−k1≠0,Ω2=16∂3f∂xn3+12∂2f∂xn22|O=a11α11−ha11α112>0.
From (117), one can say that about P5 model (15) undergoes flip bifurcation if Ω∈ℱ|P5. Moreover, period-2 points bifurcating from P5 are stable since Ω2=a11α11−ha11α112>0.


## 6.6. Analysis of Bifurcation atP6

From (57), the computation yields λ1|61=−1, but λ2,3|61=1+2/a11−k1a21−k2−a1α211−k1/α11,1−2a3/a11−k1≠1or−1. This suggests that model (15) could undergo flip bifurcation around P6 if Ω passes the curve:(118)ℱ|P6=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P6 if Ω∈ℱ|P6.Theorem 14.
Model (15) undergoes flip bifurcation around P6 if Ω∈ℱ|P6.Proof.
Recall that ifΩ∈ℱ|P6, then λ1|61=−1, but λ2,3|61=1+2/a11−k1a21−k2−a1α211−k1/α11,1−2a3/a11−k1≠1or−1. So, hereafter, detailed flip bifurcation is explored if Ω varies in the nbhd of h, i.e.,h=h+ε, by assuming h≠2/a3,2α11/a1α211−k1−a2α111−k2. Let(119)un=xn−1−k1a1α11,vn=yn,wn=zn−a3α33.
Then, (15) gives(120)un+1vn+1vn+1=1−ha11−k10001+ha21−k2−a1α211−k1α1100−hα32a3α321−ha3unvnwn+ha11−k12α11−hα11un+a11−k1α112+εa11−k1un+a11−k1α11−εα11un+a11−k1α112−hα22vn2−hα22unvn+εa21−k2vn−εα22vn2−εα21vnun+a11−k1α11ha32α33−α33wn+a3α332−α32vnwn+εa3wn+a3α33−εα33wn+a3α332−εα32vnwn+a3α33.
Using transformation,(121)unvnwn=01000b23101xnynzn.
(120) becomes(122)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+F1xn,yn,zn,εG1xn,yn,zn,εH1xn,yn,zn,ε,where(123)b23=α33−a2α11−a3α11+a2k2α11+a1α21−a1k1α21a3α11α32,F1=−zna21−k2ε+znyn+a11−k1α11α21ε+hznα22+hynznα22+zn2b23α22ε+a3xn+yn+a3α33ε−znb23α32xn+yn+a3α33ε−xn+yn+a3α332α33ε+h−xn+ynznb23α32+a32α33−xn+yn+a3α332α33,G1=a11−k1yn+a11−k1α11ε−hyn+a11−k1α112+ha121−k12α11−yn+a11−k1α112α11ε,H1=−zna21−k2ε+znyn+a11−k1α11α21ε+hznα22+hynznα22+zn2b23α22ε.
Now, from model (122) on the center manifold,(124)Wc0=xn,yn,zn|yn,zn=χ3xn,χ4xn,χi0=0,Dχi0=0,i=3,4,where(125)χixn=aix2+bix3+Ox4,fori=3,4.
From (122) and (124), one has(126)χ3−xn+F1xn,χ3,χ4=λ2χ3xn+G1xn,χ3,χ4,χ4−xn+F1xn,χ3,χ4=λ3χ4xn+H1xn,χ3,χ4.
From (126), the calculation yields: a3=b3=a4=b4=0. Thus, map (122); restrict to Wc0 as(127)fxn=−xn+a3xn+a3α33ε−xn+a3α332α33ε+ha32α33−xn+a3α332α33.
From (117) and (127), the computation yields: Ω1=3a3−ha32/α33≠0 and Ω2=h2α332>0. This implies that about P6 model (15) undergoes flip bifurcation if Ω∈ℱ|P6. Moreover, period-2 points bifurcating from P6 are stable since Ω2=h2α332>0.


## 6.7. Analysis of Bifurcation atP7

From (64), the computation yields λ1|68=−1, but λ2,3|68=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1+2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. This suggests that model (15) could undergo flip bifurcation around P7 if Ω passes the curve:(128)ℱ|P7=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P7 if Ω∈ℱ|P7.Theorem 15.
Model (15) undergoes flip bifurcation around P7 if Ω∈ℱ|P7.Proof.
Recall that ifΩ∈ℱ|P7, then λ1|68=−1, but λ2,3|68=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1+2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. So, in the following, flip bifurcation is explored by assuming h≠2α11/a2α111−k2−2a1α211−k1,2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1. Let(129)un=xn−1−k1a1α11,vn=yn−a21−k2α11−a11−k1α21α11α22,wn=zn.
Then, (15) gives(130)un+1vn+1vn+1=1−ha11−k100−hα21a21−k2α11−a11−k1α21α11α221−ha21−k2−2a1α211−k1α110001+ha3−a21−k2α32α22+a1α21α321−k1α11α22unvnwn+ha11−k12α11−α11un+a11−k1α112+εa11−k1un+a11−k1α11−α11un+a11−k1α112ha21−k2a21−k2α11−a11−k1α21α11α22+α22vn+a21−k2α11−a11−k1α21α11α222−α21a11−k1α11a21−k2α11−a11−k1α21α11α22+εa21−k2vn+a21−k2α11−a11−k1α21α11α22−α22vn+a21−k2α11−a11−k1α21α11α222−α21un+a11−k1α11vn+a21−k2α11−a11−k1α21α11α22h−α33wn2−α32vnwn+εa3wn−α33wn2−α32wnvn+a21−k2α11−a11−k1α21α11α22.
Using transformation,(131)unvnwn=c1100110001xnynzn.
(130) gives(132)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+F2xn,yn,zn,εG2xn,yn,zn,εH2xn,yn,zn,ε,where(133)c11=α11α22ha11−k1−ha21−k2−2ha11−k1α21/α11hα21a21−k2α11−a11−k1α21,F2=ha11−k12/α11−xnc11+a11−k1/α112α11+εa11−k1xnc11+a11−k1/α11−xnc11+a11−k1/α112α11c11,G2=ha21−k2a2a111−k2−a11−k1α21α11α22−a11−k1α21a2a111−k2−a11−k1α21α112α22+xn+yn+a2a111−k2−a11−k1α21α11α222α22ha11−k12/α11−xnc11+a11−k1/α112α11+εa11−k1xnc11+a11−k1/α11−xnc11+a11−k1/α112α11c11+εa21−k2xn+yn+a2a111−k2−a11−k1α21α11α22+a11−k1α21xn+yn+a2a111−k2−a11−k1α21/α11α223xnc11a2a111−k2−a11−k1α21/α11α22α22a11,H2=h−xn+ynznα32−z2α33+εzna3−znxn+yn+a2a111−k2−a11−k1α21α11α22α32−zn2α33.
Now, using system (132) on the center manifold,(134)Wc0=xn,yn,zn|yn,zn=χ5xn,χ6xn,χi0=0,Dχi0=0,i=5,6,where(135)χixn=aix2+bix3+Ox4,fori=5,6.
In view of (132) and (134), we obtain(136)χ5−xn+F2xn,χ5,χ6=λ2χ5xn+G2xn,χ5,χ6,χ6−xn+F2xn,χ5,χ6=λ3χ6xn+H2xn,χ5,χ6.
From (136), one gets: a5=b5=a6=b6=0. Finally, map (132); restrict to Wc0 as(137)fxn=−xn+ha11−k12/α11−xnc11+a11−k1/α112α11+εa11−k1xnc11+a11−k1/α11−xnc11+a11−k1/α112α11c11.
From (117) and (137), the computation yields: Ω1=a11−k11−2c11≠0. Moreover, Ω2=hα11c11+2α11c11ε2>0. This implies that about P7 model (15) undergoes flip bifurcation if Ω∈ℱ|P7. Moreover, period-2 points bifurcating from P7 are stable since Ω2=hα11c11+2α11c11ε2>0.


## 6.8. Analysis of Bifurcation atP8

From (71), the computation yields λ1|75=−1, but λ2,3|75=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1−2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. This suggests that model (15) could undergo flip bifurcation around P8 if Ω passes the curve:(138)ℱ|P8=Ω:h=2a11−k1.The proof of the following theorem shows that model (15) undergoes flip bifurcation around P8 if Ω∈ℱ|P8.Theorem 16.
Model (15) undergoes flip bifurcation around P8 if Ω∈ℱ|P8.Proof.
Recall that ifΩ∈ℱ|P8, then λ1|75=−1, but λ2,3|75=1−2/a11−k1a21−k2−2a1α211−k1/α11, 1−2/a11−k1a3−a21−k2α32/α22+a1α21α321−k1/α11α22≠1or−1. So, in the following, flip bifurcation is explored by assuming h≠2α11/a2α111−k2−2a1α211−k1,2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1. Let(139)un=xn−1−k1a1α11,vn=yn−a21−k2α11−a11−k1α21α11α22,wn=zn−a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33.
Then, (15) becomes(140)un+1vn+1vn+1=1−ha11−k100−hα21a21−k2α11−a11−k1α21α11α221−ha21−k2−2a1α211−k1α1100−hα32a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α331−ha3−a21−k2α32α22+a1α21α321−k1α11α22unvnwn+F3xn,yn,zn,εG3xn,yn,zn,εH3xn,yn,zn,ε,where(141)F3=ha11−k12α11−α11un+a11−k1α112+εa11−k1un+a11−k1α11−α11un+a11−k1α112,G3=ha21−k2a21−k2α11−a11−k1α21α11α22−α21vn+a21−k2α11−a11−k1α21α11α222−α21unvn−α21a11−k1α11a21−k2α11−a11−k1α21α11α22+εa21−k2vn+a21−k2α11−a11−k1α21α11α22−α22vn+a21−k2α11−a11−k1α21α11α222−α21un+a11−k1α11vn+a21−k2α11−a11−k1α21α11α22,H3=ha3a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−α32vnwn−α33wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332−α32a21−k2α11−a11−k1α21α11α22a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33+εa3wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−α33wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332−α32vn+a21−k2α11−a11−k1α21α11α22ε×wn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33.
Now, by utilizing transformation,(142)unvnwn=d1100d21d220111xnynzn,gives(143)xn+1yn+1zn+1=−1000λ2000λ3xnynzn+F2∗xn,yn,zn,εG2∗xn,yn,zn,εH2∗xn,yn,zn,ε,where(144)d11=α112ha1−ha1k1−ha21−k2−2ha11−k1α21/α11α222α33ha1−ha1k1−ha3−a21−k2α32/α22+a11−k1α21α32/α11α22h2α21a21−k2α11−a11−k1α21α32a3α11α22−a21−k2α11α32+a11−k1α21α32,d21=α11α22ha1−ha1k1−ha3−a21−k2α32/α22+a11−k1α21α32/α11α22α33hα32a3α11α22−a21−k2α11α32+a11−k1α21α32,d22=−α33−a2α11α22+a3α11α22+a2k2α11α22−2a1α21α22+2a1k1α21α22−a2α11α32+a2k2α11α32+a1α21α32−a1k1α21α32α32a3α11α22−a2α11α32+a2k2α11α32+a1α21α32−a1k1α21α32,F2∗=ha11−k12/α11−xnd11+a11−k1/α112α11+εa11−k1xnd11+a11−k1/α11−xnd11+a11−k1/α112α11d11,G2∗=−d21ha11−k12/α11−xnd11+a11−k1/α112α11+εa11−k1xnd11+a11−k1/α11−xnd11+a11−k1/α112α11d11d22+hd22−xnd11d21xn+d22ynα21−α21d21xn+d22yn+a2a111−k2−a11−k1α21α11α222+a21−k2a2a111−k2−a11−k1α21α11α22−a11−k1α21a2a111−k2−a11−k1α21α112α22+εd22a21−k2d21xn+d22yn+a2a111−k2−a11−k1α21α11α22−a11−k1α21d11x+a11−k1/α11d21xn+d22yn+a2a111−k2−a11−k1α21/α11α22α11−d21xn+d22yn+a2a111−k2−a11−k1α21α11α222α22,H2∗=h−yn+2znd21xn+d22ynα32−α32xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332+a3a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−a2a111−k2−a11−k1α21α32a3α11α22−a21−k2α11α32+a11−k1α21α32α112α222α33+d21−d22ha11−k12/α11−d11xn+a11−k1/α112α11+ϵa11−k1d11xn+a11−k1/α11−d11xn+a11−k1/α112α11d11d22−hd22−d11xnd21xn+d22ynα21−α21d21xn+d22yn+a2a111−k2−a11−k1α21α11α222+a21−k2a2a111−k2−a11−k1α21α11α22−a11−k1α21a2a111−k2−a11−k1α21α112α22+εd22a21−k2d21xn+d22yn+a2a111−k2−a11−k1α21α11α22−a11−k1α21d11xn+a11−k1/α11d21xn+d22yn+a2a111−k2−a11−k1α21/α11α22α11−d21xn+d22yn+a2a111−k2−a11−k1α21α11α222α22+εa3xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−α32d21xn+d22yn+a2a111−k2−a11−k1α21α11α22×xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α33−xn+yn+zn+a3α11α22−a21−k2α11α32+a11−k1α21α32α11α22α332α33.
Now, from system (143) on the center manifold(145)Wc0=xn,yn,zn|yn,zn=χ7xn,χ8xn,χi0=0,Dχi0=0,i=7,8,where(146)χixn=aix2+bix3+Ox4,fori=7,8.
From (143) and (145), we obtain(147)χ7−xn+F2∗xn,χ7,χ8=λ2χ7xn+G2∗xn,χ7,χ8,χ8−xn+F2∗xn,χ7,χ8=λ3χ8xn+H2∗xn,χ7,χ8.
From (147), computation yields: a7=b7=a8=b8=0. Finally, map (143); restrict to Wc0 as(148)f1∗xn=−xn+ha11−k12/α11−xnd11+a11−k1/α112α11+εa11−k1xnd11+a11−k1/α11−xnd11+a11−k1/α112α11d11.
From (117) and (148), the computation yields: Ω1=a11−k11−2d11≠0. Moreover, Ω2=hα11d11+2α11d11ε2>0. This implies that about P8 model (15) undergoes flip bifurcation if Ω∈ℱ|P8. Moreover, period-2 points bifurcating from P8 are stable since Ω2=hα11d11+2α11d112>0.


## 7. Numerical Simulations

Numerical simulations of three-species model (15) are performed in this section to check previous theoretical findings and to show rich dynamical behaviors. In this regard, following eight cases are presented to address the accuracy of theoretical results obtained about fixed points for model (15):Case I: if a1=4.1,a2=4.2,a3=0.08,k1=2.2,k2=0.6,α11=0.04,α22=0.4,α21=0.09,α32=0.09,andα33=0.4, then, from (29), one gets h=0.4065040650406504. From (27), if h=0.01<min0.4065040650406504,0.20703933747412007 and starting from 0.9,0.1,0.2, then Figure 1(a) indicates that P1 of (15) is a source. However, if h=0.5>max0.406504065040,0.2070393374741, then Figure 1(b) indicates that P1 of (15) is a saddle. Hence, theoretical results obtained in Lemma 2 coincide with numerical simulations.Case II: if a1=4.1,a2=4.2,a3=2.7,k1=2.2,k2=0.6,α11=1.2,α22=1.4,α21=1.9,α32=0.9,andα33=0.4, then, from (35), one gets h=0.7407407407407407. Figure 2(a) indicates if h=0.01<0.7407407407407407, then P2 of (15) is a sink. However, if h=0.95>0.7407407407407407, then Figure 2(b) indicates that P2 is unstable. Moreover, if h=0.7407407407407407, then P2 exchanges the stability, and in fact, flip bifurcation takes place by Theorem 10. Therefore, the flip bifurcation diagrams are presented in Figure 3. Finally, maximum Lyapunov exponents corresponding to Figure 3 are drawn in Figure 4.Case III: if a1=4.1,a2=4.2,a3=2.7,k1=2.2,k2=1.6,α11=1.2,α22=1.4,α21=1.9,α32=0.9,andα33=0.4, then, from (42), one gets h=0.7936507936507935. Hence, P3 is stable if h<0.7936507936507935, and exchange stability is h=0.7407407407407407, and in fact, flip bifurcation takes place by Theorem 11. Therefore, the flip bifurcation diagrams with initial value 0,0.3,0.1 are presented in Figure 5. Finally, maximum Lyapunov exponents corresponding to Figure 5 are drawn in Figure 6.Case IV: if a1=5.1,a2=4.2,a3=2.7,k1=0.09,k2=1.6,α11=1.2,α22=1.4,α21=1.9,α32=0.9,andα33=0.4, then, from (48), one gets h=0.43094160741219567. Hence, P4 is stable if h<0.43094160741219567, and exchange stability is h=0.43094160741219567, and in fact, flip bifurcation takes place by Theorem 12. Therefore, the flip bifurcation diagrams with initial value 0.1,0.2,0.1 are presented in Figure 7. Finally, maximum Lyapunov exponents corresponding to Figure 7 are drawn in Figure 8.Case V: if a1=7.1,a2=4.2,a3=2.7,k1=1.9,k2=1.6,α11=1.2,α22=1.4,α21=1.9,α32=0.9,andα33=0.4, then, from (54), one gets h=0.31298904538341166. Hence, P5 is stable if h<0.31298904538341166, and exchange stability is h=0.31298904538341166, and in fact, flip bifurcation takes place by Theorem 13. Therefore, the flip bifurcation diagrams with initial value 0,0.1,0.1 are presented in Figure 9. Finally, maximum Lyapunov exponents corresponding to Figure 9 are drawn in Figure 10.Case VI: if a1=7.1,a2=4.2,a3=2.7,k1=0.0009,k2=1.6,α11=1.2,α22=1.4,α21=1.9,α32=0.9,α33=0.4 then from (61) one gets: h=0.2819438903463822. Hence P6 is stable if h<0.2819438903463822, and exchange stability if h=0.2819438903463822 and infact flip bifurcation takes place by Theorem 14. Therefore the flip bifurcation diagrams with initial value 0.1,0,0.1 are presented in Figure 11. Finally maximum lypunov exponents corresponding to Figure 11 are drawn in Figure 12.Case VII: If a1=9.1,a2=4.2,a3=2.7,k1=0.23,k2=4.6,α11=1.2,α22=1.4,α21=0.099,α32=5.9,andα33=2.4, then, from (68), one gets h=0.2854288568574283. Hence, P7 is stable if h<0.2854288568574283, and exchange stability is h=0.2854288568574283, and in fact, flip bifurcation takes place by Theorem 15. Therefore, the flip bifurcation diagrams with initial value 0.2,0.1,0 are presented in Figure 13. Finally, maximum Lyapunov exponents corresponding to Figure 13 are drawn in Figure 14.Case VIII: if a1=9.1,a2=1.2,a3=2.7,k1=0.9,k2=0.006,α11=4.2,α22=1.4,α21=4.9,α32=1.9,andα33=0.4, then, from (75), one gets h=2.1978021978021984. Hence, P7=0.21666666666666,0.09366666666666,6.305083333333 is stable if h<2.1978021978021984, and exchange stability is h=2.1978021978021984, and in fact, flip bifurcation takes place by Theorem 16. Therefore, flip bifurcation diagrams with 0.2,0.2,0.4 are presented in Figure 15 which indicates that period-2 points bifurcate from P8 are stable, since Ω2=hα11d11+2α11d112=310.844023668>0. Finally, maximum Lyapunov exponents corresponding to Figure 15 are drawn in Figure 16.Figure 1
Phase portrait aboutP1.
(a)(b)Figure 2
Phase portrait aboutP2.
(a)(b)Figure 3
Flip bifurcation diagrams atP2, where h∈ [0.1, 6.4].
(a)(b)(c)Figure 4
Maximum Lyapunov exponents corresponding to Figure3.Figure 5
Flip bifurcation diagrams atP3, where h∈ [0.1, 1.4].
(a)(b)(c)Figure 6
Maximum Lyapunov exponents corresponding to Figure5.Figure 7
Flip bifurcation diagrams atP4, where h∈ [0 : 1; 1 : 4].
(a)(b)(c)Figure 8
Maximum Lyapunov exponents corresponding to Figure7.Figure 9
Flip bifurcation diagrams atP5, where h∈ [0.1, 1.4].
(a)(b)(c)Figure 10
Maximum Lyapunov exponents corresponding to Figure9.Figure 11
Flip bifurcation diagrams atP6, where h∈ [0.1, 1.4].
(a)(b)(c)Figure 12
Maximum Lyapunov exponents corresponding to Figure11.Figure 13
Flip bifurcation diagrams atP7, where h∈ [0.1, 1.0].
(a)(b)(c)Figure 14
Maximum Lyapunov exponents corresponding to Figure13.Figure 15
Flip bifurcation diagrams atP8, where h∈ [0.1, 4.2].
(a)(b)(c)Figure 16
Maximum Lyapunov exponents corresponding to Figure15.


## 8. Conclusion

The work is about the existence of fixed points, topological classifications around fixed points, periodic points, and bifurcations of a three-species discrete food chain model with harvesting in the region:ℝ+3=x,y,z:x,y,z≥0. We proved that, for all parametric values h, a1, a2, a3, k1, k2, α11, α22, α21, α32, and α33, model (15) has trivial fixed point: P1=0,0,0; boundary fixed points: P2=0,0,a3/α33∀a3,α33>0; P3=0,1−k2a2/α22,0ifk2<1; P4=1−k1a1/α11,0,0ifk1<1; P5=0,1−k2a2/α22,a3α22−α321−k2a2/α22α33ifa3>α321−k2a2/α22withk2<1; P6=1−k1a1/α11,0,a3/α33ifk1<1; P7=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,0ifa2>a11−k1α21/1−k2α11withk1,k2<1. We also proved that if k1<1, a2>a11−k1α21/1−k2α11 and a3>a21−k2α11α32−a11−k1α21α32/α11α22; then, P8=1−k1a1/α11,a21−k2α11−a11−k1α21/α11α22,a3α11α22−a21−k2α11α32+a11−k1α21α32/α11α22α33 is an interior equilibrium point of (15). Furthermore, we studied the local stability with different topological classifications around each fixed points whose main findings are presented in Table 1. Next, for under consideration model (15), we also studied existence of periodic points by existing theory. Furthermore, we explored the existence of possible bifurcations about each fixed points in order to understand dynamics of model (15) deeply. It is proved that (i) around P1 model undergoes no flip bifurcation if Ω∈ℱ|P1=Ω:h=2/a1k1−1, (ii) around P2 model undergoes flip bifurcation if Ω∈ℱ|P2=Ω:h=2/a3, (iii) around P3 model undergoes flip bifurcation if Ω∈ℱ|P3=Ω:h=2/a21−k2, (iv) around P4 model undergoes flip bifurcation if Ω∈ℱ|P4=Ω:h=2/a11−k1, (v) around P5 model undergoes flip bifurcation if Ω∈ℱ|P5=Ω:h=2/a1k1−1, (vi) around P6 model undergoes flip bifurcation if Ω∈ℱ|P6=Ω:h=2/a11−k1, (vii) around P7 model undergoes flip bifurcation if Ω∈ℱ|P7=Ω:h=2/a11−k1, and (viii) around P8 model undergoes flip bifurcation if Ω∈ℱ|P8=Ω:h=2/a11−k1. Finally, obtained results are verified numerically. This research can provide a framework for theoretical basis and help for the research in different aspects of biology specifically in the field of ecology.Table 1
Dynamical classifications around fixed points of model (15).
Fixed pointsCorresponding behaviorP1Not sink; source if0<h<min2/a1k1−1,2/a2k2−1; saddle ifh>max2/a1k1−1,2/a2k2−1; nonhyperbolic if h=2/a1k1−1 or h=2/a2k2−1P2Sink ifh>max2/a1k1−1,2/a2k2−1and0<h<2/a3Source if 0<h<min2/a1k1−1,2/a2k2−1andh>2/a3Saddle if 0<h<min2/a1k1−1,2/a2k2−1and0<h<2/a3Nonhyperbolic if h=2/a3 or h=2/a1k1−1 or h=2/a2k2−1P3Sink ifh>max2/a1k1−1,2α22/a2α321−k2−a3α22and0<h<2/a21−k2Source if 0<h<min2/a1k1−1,2α22/a2α321−k2−a3α22andh>2/a21−k2Saddle if h>max2/a1k1−1,2/a21−k2and0<h<2α22/a2α321−k2−a3α22Nonhyperbolic if h=2/a21−k2 or h=2/a1k1−1 or h=2α22/a2α321−k2−a3α22P4Not sink; source ifh>max2/a11−k1,2/a2k2−1Saddle if 0<h<min2/a11−k1,2/a2k2−1Nonhyperbolic if h=2/a11−k1 or h=2/a2k2−1P5Sink ifh>max2/a1k1−1,2α22/a3α22−a2α321−k2and0<h<2/a21−k2Source if 0<h<min2/a1k1−1,2α22/a3α22−a2α321−k2andh>2/a21−k2Saddle if h>max2/a1k1−1,2α22/a3α22−a2α321−k2andh>2/a21−k2Nonhyperbolic if h=2/a1k1−1 or h=a3α22−2α22/a2α321−k2 or h=2/a21−k2P6Sink if0<h<min2/a11−k1,2/a3andh>2α11/a1α211−k1−a2α111−k2;source if h>max2/a11−k1,2/a3and0<h<2α11/a1α211−k1−a2α111−k2;saddle if 0<h<min2/a11−k1,2/a3and0<h<2α11/a1α211−k1−a2α111−k2;nonhyperbolic if h=2/a11−k1 or h=2/a3 or h=2α11/a1α211−k1−a2α111−k2.P7Sink if0<h<min2/a11−k1,2α11/a2α111−k2−2a1α211−k1 andh>2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1Source if h>max2/a11−k1,2α11/a2α111−k2−2a1α211−k1 and0<h<2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1Saddle if 0<h<min2/a11−k1,2α11/a2α111−k2−2a1α211−k1 andh<2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1Nonhyperbolic if h=2/a11−k1 or h=2α11/a2α111−k2−2a1α211−k1 orh=2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1P8Sink if0<h<min2/a11−k1,2α11/a2α111−k2−2a1α211−k1,2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1Source if h>max2/a11−k1,2α11/a2α111−k2−2a1α211−k1,2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1Saddle if 0<h<min2/a11−k1,2α11/a2α111−k2−2a1α211−k1 and h>2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1 Nonhyperbolic if h=2/a11−k1 or h=2α11a2α111−k2−2a1α211−k1 or h=2α11α22/a3α11α22−a2α11α321−k2+a1α21α321−k1



---

*Source: 5987435-2021-10-14.xml*