# Hopf Bifurcation, Hopf-Hopf Bifurcation, and Period-Doubling Bifurcation in a Four-Species Food Web

Authors: Huayong Zhang; Ju Kang; Tousheng Huang; Xuebing Cong; Shengnan Ma; Hai Huang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Complex dynamics of a four-species food web with two preys, one middle predator, and one top predator are investigated. Via the method of Jacobian matrix, the stability of coexisting equilibrium for all populations is determined. Based on this equilibrium, three bifurcations, i.e., Hopf bifurcation, Hopf-Hopf bifurcation, and period-doubling bifurcation, are analyzed by center manifold theorem, bifurcation theorem, and numerical simulations. We reveal that, influenced by the three bifurcations, the food web can exhibit very complex dynamical behaviors, including limit cycles, quasiperiodic behaviors, chaotic attractors, route to chaos, period-doubling cascade in orbits of period 2, 4, and 8 and period 3, 6, and 12, periodic windows, intermittent period, and chaos crisis. However, the complex dynamics may disappear with the extinction of one of the four populations, which may also lead to collapse of the food web. It suggests that the dynamical complexity and food web stability are determined by the food web structure and existing populations.


---

## Body

## 1. Introduction

The concept and research method of the food web was firstly proposed as early as in 1930s by Elton, who initiated exploring community structure in ecosystems from the viewpoint of network [1]. After that, the food webs have been widely studied in various aspects of characteristics, involving stability, complexity, dynamics, etc. Via analyzing the characteristics of food webs, many research works have played an important role in understanding the properties of natural ecosystems [2–4].Food web systems generally are composed of different combinations of relationships between interacting species, among which predation relationship is a basic one [5, 6]. Understanding the interactions between prey and predator is a significant topic of longstanding interest in biology. The prey-predator models with various functional responses, such as Holling I, II, III, and IV, Leslie-Gower, Michaelis-Menten, ratio-dependent, Beddington-DeAngelis, and so on, have been extensively discussed [7–10]. Many other ecological mechanisms have been also included and investigated in the prey-predator models. For example, Raw et al. proposed and analyzed a tri-trophic prey-predator model with two predators and one prey exhibiting group defense mechanism and found occurrence of complex dynamical behaviors [7].Over past decades, mathematical models have been investigated to study the dynamics of three or more interacting species in food webs. Early studies mainly focused on three-species food web models, demonstrating the complexity of the dynamics for the food webs [6, 11]. The food web systems may converge not only to stable equilibria and periodic orbits, but also to more complex behaviors such as quasiperiodic oscillations and chaos induced by different bifurcations [11]. As studied by Gakkhar and Gupta [6], generalized Hopf bifurcation can exhibit in a food web system with commensalism and Bogdanov-Takens bifurcation point is detected only for a certain value of commensal coefficient. Gupta and Chandra [11] investigated a prey-predator-scavenger model with quadratic harvesting and found that all three species have unpredictable and irregular behavior which may be controlled with the help of an appropriate number of harvesting so that the system can be dynamically balanced.Currently, research on food web dynamics slipped into the cases where four species or more species exist and interact [2, 12, 13]. Rich dynamic behaviors of four-species or more species food webs, such as equilibria, bifurcations, periodic orbits, torus, chaotic attractors, and so on, have been found to be contributing to better comprehending on the complexity of species coexistence and community stability. In Gakkhar et al. [2], a four-species food web model was developed and complex dynamic behaviors including strange attractor and chaos were shown. In the food web studied by Wei [13], Wei proved that the food web system possesses a line segment of degenerate steady states for the parameter values on a bifurcation line in the bifurcation diagram and found diverse routes to chaos.Investigation has demonstrated that the topological structure of food web can play an important role in affecting the food web dynamics. Generally, with the variation of variable diversity and trophic links, the structure of food webs may be different. In Fussmann and Gerd [14], the architecture of a food web with three primary producers, two consumers, one primary predator, one secondary predator, and one tertiary predator was emphasized. Many basic food webs that differ in complexity and in the number of trophic levels and omnivorous links can be generated as substructures of this web. According to the difference of network structure, Kuijper et al. [15] classified the four-species food webs into three types, i.e., competition, omnivory, and food-chain, and found that the architecture of the food web depends on the feeding behaviors of the highest trophic level. Based on the research works of Fussmann and Gerd [14] and Kuijper et al. [15], this research focuses on tri-trophic level food web which includes two-preys, one-middle predator feeding on the two-preys, and one-top predator feeding on all three other populations. The study on the bifurcations and complex dynamics of this food web is still not documented in literature available to the authors.In this research, we modeled the above food web on the basis of the Hasting-Powell [16] model with application of Holling II functional response, and the exploration is arranged as follows. In Section 2, we introduce the four-species food web model and analyze local stability of the coexisting equilibrium. In Section 3, the bifurcation theorem and center manifold theorem are used to determine the conditions for Hopf bifurcation of the food web. Section 4 investigates the Hopf-Hopf bifurcation and performs simulations to demonstrate the quasiperiodic dynamics. Section 5 numerically studies the period-doubling bifurcation and demonstrates the bifurcation diagram in which a route to chaos emerges. In Section 6, discussion and conclusions are provided.


## 2. Mathematical Model and Stability Analysis

In 2002, Fussmann and Gerd [14] have studied the food web complexity and chaotic population dynamics, with proposition of architecture of food web models with two, three, and four to five trophic levels. Their results suggested that natural food webs possess architectural properties that may inherently lower the maybe of chaotic community dynamics [14]. Based on their research work, we choose to study the architecture of food web models with three trophic levels including two preys, a middle predator, and a top predator. With the application of Hasting-Powell [16] model and Holling II functional response, a four-species food web model can be described by a dynamical system of four ordinary differential equations: (1)dXdT=R1X1-XK1-F1XZ-F2XW,dYdT=R2Y1-YK2-F3YZ-F4YW,dZdT=E1F1XZ+E3F3YZ-F5ZZ-I1Z,dWdT=E2F2XW+E4F4YW+E5F5ZW-I2W,with(2)FiU=AiUBiU+1,i=1,…,5.In the above equations, X, Y, Z, and W are population densities of the four interacting species and T is time; R1 and R2 represent the intrinsic growth rates of species X and Y; K1 and K2 are the carrying capacities; Ei(i=1,2,…,5) represents the conversion rate of prey to predator; I1 and I2 are the death rates of species Z and W, respectively. FiU represents a functional response of the Holling-II type; Ai and Bi(i=1,2,…,5) are parameters used in the functional response: Ai describes the effects of capture rate and handling time; Bi determines how fast per capita feeding rate approaches its saturation rate. From the viewpoint of ecology, parameters R1, R2, K1, K2, Ii, Ei, Ai and Bi should all be positive constants.It should also be noticed that population densities are nonnegative. Hence, the state space of the food web system (1) must be strictly restricted in the region R+4=X,Y,Z,W∈R4:X≥0,Y≥0,Z≥0,W≥0. In order to reduce the number of parameters in (1), and to determine which combinations of parameters control the behaviors of the system, system (1) is changed with the application of the following expressions: (3)x=XK1,y=YK2,z=ZA1,w=WA2,t=T.Using (3) to replace the corresponding variables in (1), then we obtain a simpler form of the food web system: (4)dxdt=rx1-x-xza1x+1-xwa2x+1=ψx,y,z,w,dydt=ky1-y-a3yza4y+1-a5ywa6y+1=ξx,y,z,w,dzdt=a7xza1x+1+a8yza4y+1-a9zwa10z+1-a11z=φx,y,z,w,dwdt=a12xwa2x+1+a13ywa6y+1+a14zwa10z+1-a15w=ϕx,y,z,w,in which (5)r=R1,k=R2,a1=B1K1,a2=B2K1,a3=A3A1,a4=B3K2,a5=A4A2,a6=B4K2,a7=A1E1K1,a8=A3E3K2,a9=A5A2,a10=B5A1,a11=I1,a12=A2E2K1,a13=A4E4K2,a14=E5A5A1,a15=I2.The equilibriums of system (4) can be obtained by soling equations dx/dt=0, dy/dt=0, dz/dt=0, and dw/dt=0. Direct calculation yields that system (4) has following six kinds of equilibriums:(i) Trivial equilibriumE0=(0,0,0,0) always exists, describing the absence of populations.(ii) Axial equilibriumsE10=(1,0,0,0) and E11=(0,1,0,0) indicate the prey population x or y increases to carrying capacity in the absence of predation. Likewise, equilibrium E2=(1,1,0,0) shows the case where prey populations x and y increase to the carrying capacity simultaneously in the absence of predation.(iii) When one predator and one prey are absent simultaneously, system (4) degenerates to be a normal Holling II predator-prey system. In such cases, we can find the equilibriums for coexistence of predator and prey, including EXZ, EXW, EYZ, and EYW.(iv) With the consideration of reducing one predator population, system (4) degenerates to be a two-trophic-level food web system with two preys and one predator. The equilibriums of such type of food web have been studied by Misra et al. [17].(v) The equilibriums for the case with absence of one prey population have been investigated by Hsu et al. [18]. One can refer to that literature for detailed calculations of equilibriums.(vi) This research focuses on the positive interior equilibriumE∗=(x∗,y∗,z∗,w∗), which represents the coexistence of the four species in the food web. From system (4), it is easy to know that (x∗,y∗,z∗,w∗) is a positive root of the following nonlinear algebraic equations: (6)r1-x-za1x+1-wa2x+1=0,k1-y-a3za4y+1-a5wa6y+1=0,a7xa1x+1+a8ya4y+1-a9wa10z+1-a11=0,a12xa3x+1+a13ya6y+1+a14za10z+1-a15=0.To determine the local asymptotical stability of the equilibriumE∗, the method of Jacobian matrix is utilized. The Jacobian matrix associated to system (4) at any point E(x,y,z,w) can be described as follows: (7)JE=∂ψ∂x∂ψ∂y∂ψ∂z∂ψ∂w∂ξ∂x∂ξ∂y∂ξ∂z∂ξ∂w∂φ∂x∂φ∂y∂φ∂z∂φ∂w∂ϕ∂x∂ϕ∂y∂ϕ∂z∂ϕ∂w.The characteristic equation of the Jacobian matrix J(E) at E∗ can be written as(8)Pλ=λ4+δ3λ3+δ2λ2+δ1λ+δ0=0,where(9)δ0=c14c23c32c41-c13c24c32c41-c14c22c33c41+c13c22c34c41-c14c23c31c42+c13c24c31c42-c11c24c33c42+c11c23c34c42+c14c22c31c43+c11c24c32c43-c11c22c34c43-c13c22c31c44-c11c23c32c44+c11c22c33c44,δ1=c13c22c31+c11c23c32-c11c22c33+c14c22c41+c14c33c41-c13c34c41+c11c24c42+c24c33c42-c23c34c42-c14c31c43-c24c32c43+c11c34c43+c22c34c43-c11c22c44+c13c31c44+c23c32c44-c11c33c44-c22c33c44,δ2=c11c22-c13c31-c23c32+c11c33+c22c33-c14c41-c24c42-c34c43+c11c44+c22c44+c33c44,δ3=-c11-c22-c33-c44,and cjk is the element of Jacobian matrix J(E∗) at the jth row and kth volume.Define(10)Δ1=δ3,Δ2=δ1δ2-δ1,Δ3=δ1Δ2-δ0δ32,Δ4=δ0Δ3.According to the Routh-Hurwitz criterion, when the following conditions(11)Δ1>0,Δ2>0,Δ3>0,Δ4>0are satisfied simultaneously, the equilibrium E∗ is locally and asymptotically stable. Otherwise, this equilibrium becomes unstable.A numerical example is provided here to illustrate the emergence of stable coexistent equilibrium of the food web system. We set parametric conditions of system (4) as follows:(12)r=1.04,a1=0.5,a2=0.04,a3=0.84,a4=5,a5=3,a6=0.2,a7=0.08,a8=0.1,a9=0.05,a10=0.52,a11=0.02,a12=0.01,a13=0.4,a14=0.004,a15=0.06,k=1.001.Figure1 shows the dynamics of the food web when the parameter values are given in (12). In the graph of time series (Figure 1(a)), variables x, y, z, w oscillate at the beginning with decreasing amplitude and become steady in long-term behaviors. In the x-y-z and y-z-w phase diagrams (Figures 1(b) and 1(c)), it is indicated that the coexistent equilibrium E∗ is a stable focus. It can be verified that (13)Δ1=0.192005,Δ2=0.00153,Δ3=8.9438×10-6,Δ4=1.67489×10-9.From determinative condition (11), we also know that the coexistent equilibrium E∗ is locally and asymptotically stable.Figure 1
(a) Time series of variablesx, y, z, w; (b) and (c) dynamics of the food web shown in x-y-z and y-z-w phase diagrams. The parameter values are provided in (12).
(a) (b) (c)


## 3. Hopf Bifurcation

In this section, Hopf bifurcation theorem and center manifold theorem [19] are used to determine the occurrence conditions of Hopf bifurcation for the four-species food web model. To study the Hopf bifurcation, bifurcation parameter should be chosen at first. Among all parameters of system (4), parameter r is basic and represents the growth rate of prey X, controlling one of energy inputs of the food web and therefore we choose it as the bifurcation parameter.According to the Hopf bifurcation theorem, when system (4) undergoes Hopf bifurcation, the characteristic equation Pλ=0 should have a pair of conjugate pure imaginary roots, i.e.,(14)λ1,2=±iωi=-1.Substituting the eigenvalues (14) into the characteristic equation (8), and separating the real and imaginary parts, then we have(15)δ0+ω4-δ2ω2=0,(16)δ3ω3-δ1ω=0.Reducing the unknown quantity ω, the above equations (15) and (16) turn to be (17)δ0δ32+δ12-δ1δ2δ3=0.Equation (17) determines the relationship of the parameters which exactly leads to the emergence of Hopf bifurcation. With the other parameter values given, we can calculate the threshold value of the bifurcation parameter r. In the following, we denote this threshold value of Hopf bifurcation point as r=r0.When the value of parameterr changes around the Hopf bifurcation point, the two conjugate imaginary eigenvalues ±iω will change to be complex eigenvalues, written as λ=υ+iω~. Substituting λ=υ+iω~ into (8), then we obtain the following [12]:(18)υ4+δ3υ3+δ2-6ω~2υ2+δ1-3δ3ω~2υ+δ0+ω~4-δ2ω~2=0,(19)δ3+4υω~3-3δ3υ2+3δ2υ+4υ3+δ1ω~=0.Through (19), we find ω~ can be expressed by υ. Replacing such expression of ω~ into (18) results in(20)υ4+δ3υ3+δ2-6fυυ2+δ1-3δ3fυυ+δ0+f2υ-δ2fυ=0,where f(υ)=(3δ3υ2+3δ2υ+4υ3+δ1)/(δ3+4υ). Differentiating υ with respect to r and putting r=r0, we get(21)dυdrr=r0=δ32dδ0/dr+2δ1-δ2δ3dδ1/dr-δ1δ3dδ2/dr+2δ0δ3-δ1δ2dδ3/dr3δ22δ3+2δ1δ32-2δ1δ2-8δ0δ3r=r0.If system (4) undergoes Hopf bifurcation around the interior equilibrium E∗, we should have dυ/drr=r0≠0.System (4) undergoing Hopf bifurcation at r=r0 should be also under the satisfaction of the last determinative value. We consider two conjugate imaginary eigenvalues λ1,2=±iω and two conjugate complex eigenvalues λ3,4=σ±iη for the characteristic equation (8). To calculate the determinative value, we further transform the equations of system (4) by applying the following translation:(22)x1=x-x∗,x2=y-y∗,x3=z-z∗,x4=w-w∗.Substituting (22) into the equations of system (4) and separating the linear part and nonlinear part. Then system (4) can be rewritten as(23)X˙=JE∗r=r0X+D,where(24)X=x1,x2,x3,x4T,D=D1x1,x2,x3,x4,D2x1,x2,x3,x4,D3x1,x2,x3,x4,D4x1,x2,x3,x4T,D1x1,x2,x3,x4=∑2≤j1+j2+j3+j4≤3cj1j2j3j4x1j1x2j2x3j3x4j4+Ox1+x2+x3+x44,D2x1,x2,x3,x4=∑2≤j1+j2+j3+j4≤3dj1j2j3j4x1j1x2j2x3j3x4j4+Ox1+x2+x3+x44,D3x1,x2,x3,x4=∑2≤j1+j2+j3+j4≤3mj1j2j3j4x1j1x2j2x3j3x4j4+Ox1+x2+x3+x44,D4x1,x2,x3,x4=∑2≤j1+j2+j3+j4≤3nj1j2j3j4x1j1x2j2x3j3x4j4+Ox1+x2+x3+x44,in which Ox1+x2+x3+x44 is a fourth-order polynomial functions in regard to independent variable x1,x2,x3,x4, and (25)cj1j2j3j4=1j1+j2+j3+j4!∂ψj1+j2+j3+j4∂x1j1∂x2j2∂x3j3∂x4j4,dj1j2j3j4=1j1+j2+j3+j4!∂ξj1+j2+j3+j4∂x1j1∂x2j2∂x3j3∂x4j4,mj1j2j3j4=1j1+j2+j3+j4!∂φj1+j2+j3+j4∂x1j1∂x2j2∂x3j3∂x4j4,nj1j2j3j4=1j1+j2+j3+j4!∂ϕj1+j2+j3+j4∂x1j1∂x2j2∂x3j3∂x4j4.Via introducing an invertible transformation, i.e.,(26)X=PY,where(27)P=1011p21p22p23p24p31p32p33p34p41p42p43p44,Y=y1,y2,y3,y4T,where the elements pl1l2l1,l2=1,⋯,4 of the matrix P are provided in Appendix, then the normal form of system (23) can be given by(28)Y˙=P-1JE∗PY+D~,where(29)P-1JE∗P=0-ω00ω00000σ-η00ησ,D~y1,y2,y3,y4T=P-1Dx1,x2,x3,x4T,in which c~j1j2j3j4, d~j1j2j3j4, m~j1j2j3j4, and n~j1j2j3j4 are coefficients of the polynomial function D~y1,y2,y3,y4.Next, the center manifoldWC(0,0,0,0) of system (28) at the origin is calculated. On the basis of center manifold theorem, a center manifold WC(0,0,0,0) exists and can be described as follows: (30)WC0,0,0,0=y1,y2,y3,y4∈R4∣y3=hy1,y2,y4=gy1,y2,h0,0=0,g0,0=0,Dh0,0=0,Dg0,0=0,where h(y1,y2) and g(y1,y2) are assumed as (31)hg=∑1≤j+k≤2hjky1y2+Oy1+y23∑1≤j+k≤2gjky1y2+Oy1+y23.By computation, we get(32)h10=h01=h11=0,h20=m~0200ω2,h02=-m~2000ω2,g10=g01=g11=0,g20=m~0200ησ-m~2000ηω2+n~2000σω2+n~0200ω4ω2ω4-σ2,g02=m~0200ηω2-m~2000ησ+n~0200σω2+n~2000ω4ω2ω4-σ2.Accordingly, the dynamics restricted to the center manifoldWC(0,0,0,0) can now be considered. Combining with (31), system (28) is simplified as(33)y˙1y˙2=G1y1,y2G2y1,y2,where(34)G1y1,y2=-ωy2+μ1y12+μ2y12y2+μ3y1y22+μ4y13,G2y1,y2=ωy1+γ1y22+γ2y12y2+γ3y1y22+γ4y23,in which(35)μ1=c~2000,μ2=0,μ3=c~1001g02+c~1010h02,μ4=c~3000+c~1001g20+c~1010h20,γ1=d~0200,γ2=d~0101g20+d~0110h20,γ3=0,γ4=d~0300+d~0101g02+d~0110h02.We introduce the following partial derivative notation:(36)∂G1∂y1yr0=Gy11,∂3G1∂y12∂y2yr0=Gy1y1y21,∂2G2∂y22yr0=Gy2y22,⋯,in which subscript y1 demonstrates the partial derivative of the first variable; subscript y2 is the partial derivative of the second variable. On account of (34), we have Gy11=0, Gy21=-ω≠0, Gy12=ω≠0, Gy22=0 and Gy21Gy12<0. Moreover, system (33) has two eigenvalues as ±iGy21Gy12. Consequently, the final determination for the system (4) experiencing Hopf bifurcation is calculated by the following [20]: (37)Kr0=116Gy21Gy121/2Gy1y1y11+Gy12Gy211/2Gy1y2y21+Gy21Gy121/2Gy1y1y22+Gy12Gy211/2Gy2y2y22+116Gy12Gy12Gy211/2Gy1y11Gy1y21+Gy12Gy21Gy2y21Gy1y21+Gy12Gy21Gy2y21Gy2y22+116Gy21Gy21Gy121/2Gy2y22Gy1y22+Gy21Gy12Gy1y12Gy1y22+Gy21Gy12Gy1y11Gy1y12.Direct calculation yields(38)Kr0=116c~1001g02+c~1010h02+c~3000+c~1001g20+c~1010h20+d~0101g20+d~0110h20+d~0300+d~0101g02+d~0110h02.Summarizing the above analysis, if condition (38) is satisfied, then system (4) undergoes Hopf bifurcation at equilibrium (xr0∗,yr0∗,zr0∗,wr0∗). Furthermore, when dυ/drr=r0<0 establishes, if Kr0<0 and 0<r-r0<<1, the system experiences supercritical bifurcation, and a stable orbit generates with radius as about υr0/Kr0; if Kr0>0 and 0<r0-r<<1, subcritical Hopf bifurcation occurs and results in the emergence of an unstable orbit with radius as about υr0/Kr0.For illustrating the occurrence of Hopf bifurcation in the food web system, an example is provided here. First, parameterr is still set as the bifurcation parameter and the other parameter values are given in (12). Under such parametric conditions, the bifurcation critical value of r for occurrence of Hopf bifurcation in system (4) is determined as r0=1.1621. Simultaneously, the coexistent equilibrium at the bifurcation critical point is E∗(0.2412,0.1422,0.8133,0.1575) and the corresponding eigenvalues are λ1,2=±0.1642i, λ3=-0.1260, and λ4=-0.0557, which show two pairs of imaginary numbers. Moreover, we can derive dυ/drr=r0=-0.190021<0 and the determinative value Kr0 is calculated as Kr0=-5.62119<0. Consequently, we know that system (4) experiences supercritical Hopf bifurcation at the coexistence equilibrium E∗, i.e., when r<r0, the dynamics of the food web converg to the coexistent equilibrium; when r=r0, the coexistent equilibrium is weakly attracting; when r>r0, the dynamics of the food web are attracted to a limit cycle (Figures 2 and 3).Figure 2
(a) Time series of variablesx, y, z, w; (b) and (c) dynamics of the food web in x-y-z and y-z-w diagrams. r=1.2, and the other parameter values are provided in (12).
(a) (b) (c)Figure 3
Hopf bifurcation diagrams corresponding to the variation of parameterr. The parametric conditions are given in (12).
(a) (b)Figure2 is plotted to show the dynamics of the food web when r=1.2. By Routh-Hurwitz criterion, we obtain (39)Δ1=0.1536,Δ2=-0.0001,Δ3=-4.54×10-7,Δ4=-8.56×10-11,suggesting that the coexistent equilibrium is unstable. In this case, x, y, z, and w will keep periodically oscillating with long-term behavior (Figure 2(a)). In the phase diagram, the periodically oscillating behavior is clearly manifested as a limit cycle (Figures 2(b) and 2(c)).Figure3 further demonstrates the Hopf bifurcation diagrams with the change of parameter r. As the value of parameter r increases, the system dynamics vary from stable equilibrium to limit cycle, forming bowl-shaped surfaces in Figure 3, which vividly and clearly demonstrate the occurrence of the Hopf bifurcation in the food web system.


## 4. Hopf-Hopf Bifurcation

Beyond the Hopf bifurcation, the food web system can still further undergo Hopf-Hopf bifurcation with parameter variation. Hopf-Hopf bifurcation is a complicated codimension-two bifurcation, which only appears in four-species food web system or higher-species food web systems. Hopf-Hopf bifurcation can lead to the emergence of a torus solution, appearing like tyre. Hopf-Hopf bifurcation in physical engineering field was fully studied [19, 21, 22]. However, the research on Hopf-Hopf bifurcations of food web system is still not documented in the literature available to the authors.For studying the Hopf-Hopf bifurcation, we choose two parametersk and r as the bifurcation parameters. Depending on the variation of bifurcation parameters, the food web system may undergo Hopf-Hopf bifurcation. At the bifurcation point, equilibrium E∗ has two pairs of different purely imaginary eigenvalues, denoted as(40)λ1,2=±iω1,λ3,4=±iω2,where ω1>ω2>0.Then (8) can be rewritten as(41)Pλ=λ2+ω12λ2+ω22.Balancing the right side of both (8) and (41), we easily obtain(42)δ3=0,δ2=ω12+ω22,δ1=0,δ0=ω12ω22.From the conditionsδ1r,k=0 and δ3r,k=0, we can get the threshold values of the two bifurcation parameters as r=r0 and k=k0. Moreover, if the conditions δ2r0,k0>0, δ0r0,k0>0, and δ22r0,k0-4δ0r0,k0>0 are satisfied, the values of ω1 and ω2 can be also obtained, i.e., (43)ω1=δ2r0,k0+δ22r0,k0-4δ0r0,k02,ω2=δ2r0,k0-δ22r0,k0-4δ0r0,k02.In this critical case, the food web system undergoes Hopf-Hopf bifurcation at the equilibrium E∗. To demonstrate the Hopf-Hopf bifurcation numerically, we set parameter values of the food web system as follows: (44)a1=0.86,a2=0.55,a3=0.59,a4=2.277,a5=2.5,a6=2.556,a7=2.5,a8=0.05,a9=1.9,a10=0.83,a11=0.01,a12=0.0001,a13=0.25,a14=0.03,a15=0.05.Under the parametric conditions of (44), the threshold values of bifurcation parameters at the Hopf-Hopf bifurcation point are determined as r0=1.1985, k0=1.171. Meanwhile, the coexisting equilibrium is obtained as E∗=(0.2746,0.2474,0.6053,0.4371) and its eigenvalues are λ1,2=±0.4465i and λ3,4=±0.0790i, which verify the occurrence of the Hopf-Hopf bifurcation in the four-species food web system. As shown in Figure 4, the dynamics of the system converge to a stable invariant torus in such case, which represents quasiperiodic dynamics of the food web system.Figure 4
The torus solution generated by Hopf-Hopf bifurcation of system (4). r=1.28, k=1.2001, and the others parameter values are provided in (44).
(a) (b) (c) (d)Figure5 demonstrates the dynamic transition around the Hopf-Hopf bifurcation point with the change of parameter k. As the value of parameter k increases, the system dynamics vary from limit cycle to torus (Figure 5(a)). In the Poincaré map, this phenomenon reflects as a dynamic variation from a stable fixed point to a limit cycle, as shown in Figure 5(b).Figure 5
Dynamic transition around the Hopf-Hopf bifurcation point corresponding to the variation of parameterk. The parametric conditions are provided in (44).
(a) (b)Figure6 further exhibits a bifurcation diagram when the food web system experiences Hopf-Hopf bifurcation. The bifurcation diagram is plotted with the successive maximum of z as the parameter k varies in the range [1.08,1.32]. As shown in Figure 6, with the k value increasing around the critical value k0, the system dynamics experience transition from a stable limit cycle to a torus. After that the food web system comes into the zone of quasiperiodic dynamics, in which a few periodic windows can still take place, such as periodic windows with period 7, 8 orbits (Figures 7(b) and 7(d)). These periodic dynamics are intermittent and appear and disappear abruptly, resulting in alteration of periodic and quasiperiodic dynamics with the increase of bifurcation parameter. At the end of the quasiperiodic zone, the system dynamics turn to be periodic, as demonstrated by period 4 orbit in Figure 7(f).Figure 6
Bifurcation diagram of the food web system corresponding to the variation of parameterk when Hopf-Hopf bifurcation occurs. r=1.28, and the other parameter values are given in equation (44).Figure 7
Periodic and quasiperiodic dynamics in Poincaré maps, corresponding to the dynamic variation in the bifurcation diagram of Figure6.
(a)
k = 1.12 (b)
k = 1.148 (c)
k = 1.15 (d)
k = 1.175 (e)
k = 1.23 (f)
k = 1.27


## 5. Period-Doubling Bifurcation

Based on the Hopf bifurcation, we find that the food web can further undergo period-doubling bifurcation. Under the influence of the period-doubling bifurcation, the food web system may experience period-doubling cascade, eventually leading to the occurrence of complex periodic behaviors and chaotic dynamics.In order to numerically display the period-doubling bifurcation and corresponding system dynamics, we give the parametric conditions as follows:(45)a1=0.86,a2=0.05,a3=0.5,a4=4,a5=4,a6=0.556,a7=2.2,a8=0.005,a9=1.5,a10=0.0005,a11=0.002,a12=0.002,a13=1.8,a14=0.001,a15=0.4,k=1.0001.Here, we also choose parameter r as the bifurcation parameter, and the initial condition is set as (0.5, 0.5, 0.5, 0.5). The numerical simulation results are shown in Figures 7–11. With the change of parameter r, the dynamics of all species exhibit period-doubling cascade, leading to dynamical transition from period to chaos. When r<1.16, the system dynamics exhibit simple periodic oscillation, which appears as a limit cycle in phase portrait (Figure 8(a)) and a single point in Poincaré map (Figure 9(a)). When the value of parameter r increases to r∈1.16,1.316, the period of the periodic oscillation doubles, i.e., period 2 orbit. This reflects as two cycles in Figure 8(b) and two points in the Poincaré map of Figure 9(b). Consequently with the increase of r value, we find periodic orbits with period 4, period 8, and period 16 for r∈1.316,1.342, r∈1.342,1.448, and r∈1.448,1.498, respectively. In the phase portrait and Poincaré map, these periodic dynamical behaviors represent as 4, 8, and 16 cycles and points (see Figures 8(c)-8(e) and 9(c)-9(e)). As a result of the period-doubling cascade, chaotic attractors are discovered in long-term behaviors when r>1.3498, as shown in Figures 8(f) and 9(f).Figure 8
Period-doubling cascade and chaotic dynamics of the food web system as shown in thex-y-z phase portraits. The other parameter values are provided in (45).
(a)
r = 1.15 (b)
r = 1.25 (c)
r = 1.3425 (d)
r = 1.3452 (e)
r = 1.3495 (f)
r = 1.45Figure 9
Periodic and chaotic dynamics corresponding to Figure8 as shown in the Poincaré maps.
(a)
r = 1.15 (b)
r = 1.25 (c)
r = 1.3425 (d)
r = 1.3452 (e)
r = 1.3495 (f)
r = 1.45Figure 10
Sensitivity analysis on the dynamics of food web corresponding to the parametric conditions of Figure8(f).
(a) (b)Figure 11
Bifurcation diagrams of the food web system. (a) Maxima for the population density of middle predatorz; (b) and (c) magnifications of (a) when the r value ranges in [1.41,1.435] and [1.51,1.57]. The parameter values are provided in (45).
(a) (b) (c)To further determine the chaotic behaviors of the food web, sensitivity analysis is performed via applying two close initial conditions and then observing the system dynamics in long-term behaviors. In Figure10, red and blue waves show the dynamical behaviors with reference to two close initial conditions, namely,x0 = 0.5,y0 = 0.5,z0 = 0.5, w0 = 0.5 andx0 = 0. 5001,y0 = 0.5,z0 = 0.5, w0 = 0.5. We find that a tiny change (0.0001 in magnitude ofx0) in initial conditions can result in a large difference ofz and w dynamics as the time progresses. It indicates that the solutions under such parametric conditions are highly sensitive to initial conditions, suggesting the emergence of chaotic dynamics for the food web.To explicitly demonstrate the dynamic transition between nonchaos and chaos with the increase of bifurcation parameter, bifurcation diagrams are plotted, as shown in Figure11. Firstly, we find that the period-doubling bifurcation starts a route to chaos (Figure 11(a)). Secondly, when the food web goes into the chaotic zone, periodic windows can repeatedly take place. In the periodic windows, period-doubling cascade also emerges and leads to multiperiodic behaviors (Figure 11(b)). Figure 12 displays the system dynamics in a periodic window with period 3 orbit and period-doubling cascade in orbits of period 3, 6, and 12. Thirdly, when the periodic windows end, the food web may experience chaos crisis, which leads to an abrupt transition from chaotic subsets to a chaotic band (Figure 11(c)). Lastly, the chaos disappears when parameter r increases to higher values (such as r>1.75), and the food web dynamics return to be periodic.Figure 12
Period-doubling cascade in orbits of period 3, 6, and 12, as shown in the Poincaré maps, corresponding to the periodic window in Figure11(b).
(a)
r = 1.418 (b)
r = 1.424 (c)
r = 1.4242


## 6. Discussion and Conclusion

In this research, we have studied the dynamics of a four-species food web with Holling II functional response. Three bifurcations, Hopf bifurcation, Hopf-Hopf bifurcation, and period-doubling bifurcation, are found. Various dynamical behaviors are observed, including stable equilibrium, periodic oscillations, quasiperiodic oscillations, period-doubling cascade, periodic window, chaotic interior crisis, and chaotic attractors. With the variation of bifurcation parameter, the food web system may even experience the transitions between equilibrium and period, period and quasiperiod, or nonchaos and chaos.The complex dynamics occurring in the food web result from the interactions among all populations. If one of the populations is extinct, possibly due to catastrophic shifts of environmental conditions, the system dynamics may experience dramatic change. In order to exhibit the contribution of populations on the dynamics of the food web system, we here choose three cases of typical dynamical behaviors, period (Figure2), quasiperiod (Figure 4), and chaos (Figure 8(f)) and then remove one population and observe the dynamics change. Figure 13 shows the situation where the prey population x is extinct. We find that with the absence of one prey, the middle predator or top predator is going to be extinct in the three cases of Figure 13 and that the food web system degenerates a predator-prey system. Similar phenomenon also occurs in the situation when prey y or predator w is extinct. As shown in Figure 14, it is discovered that, with the absence of predator w, the prey population x always tends to extinction and that populations y and z keep periodically oscillating in long-term behavior for all three cases. This result demonstrates that the extinction of one population can lead to instability and collapse of the food web.Figure 13
System dynamics in the situation with the absence of prey populationx, corresponding to the three cases in Figures 2, 4, and 8(f).
(a) (b) (c)Figure 14
System dynamics in the situation with the absence of predator populationw, corresponding to the three cases in Figures 2, 4, and 8(f).
(a) (b) (c)When we consider that the middle predatorz is extinct, the periodic, quasiperiodic, and chaotic behaviors of the food web system all change to be a stable equilibrium (Figure 15). Although the left three populations are still coexisting, the dynamical complexity of the food web is reduced in these cases. The results in Figures 14 and 15 suggest that the omnivorous population can balance the population abundance and enhance the stability of the food web system. Moreover, Figure 15 also implies that the food web has unpredictable and irregular dynamic behaviors which may be controlled with the help of appropriate number of predation so that the food web can be dynamical equilibrium.Figure 15
System dynamics in the situation with the absence of predator populationz, corresponding to the three cases in Figures 2, 4, and 8(f).
(a) (b) (c)



---

*Source: 8394651-2018-09-27.xml*