A Three-Species Food Chain System with Two Types of Functional Responses

Abstract

In recent decades, many researchers have investigated the ecological models with three and more species to understand complex dynamical behaviors of ecological systems in nature. However, when they studied the models with three species, they have just considered the functional responses between prey and mid-predator and between mid-predator and top predator as the same type. However, in the paper, in order to describe more realistic ecological world, a three-species food chain system with two types of functional response, Holling type and Beddington-DeAngelis type, is considered. It is shown that this system is dissipative. Also, the local and global stability of equilibrium points of the system is established. In addition, conditions for the persistence of the system are found according to the existence of limit cycles. Some numerical examples are given to substantiate our theoretical results. Moreover, we provide numerical evidence of the existence of chaotic phenomena by illustrating bifurcation diagrams of system and by calculating the largest Lyapunov exponent.

1. Introduction

Studying the dynamic interaction between predators and preys has long been one of the main themes in ecological systems. The classical ecological models of interacting populations typically have focused on two-species continuous time systems with one predator and one prey. In general, these models consist of two differential equations and a functional response, which is a function representing the prey consumption per unit time. Many scholars have studied about classical two-species continuous time systems for several functional responses such as Holling-Tanner type (cf. [1–3]), Beddington-DeAngelis type (cf. [4, 5]), and ratio-dependent type (cf. [6, 7]). However, it has been recognized that such classical models with two-species can describe only a small number of the phenomena that are commonly observed in nature. In fact, the models can exhibit only two main patterns: approach to an equilibrium point or to a limit cycle [8]. For the reason, in recent decades, many researchers (cf. [9–17]) have turned their concerns to the ecological models with three and more species to understand complex dynamical behaviors of ecological systems in the real world. They have demonstrated the very complex dynamic phenomena of those models, including cycles, periodic doubling, and chaos. In particular, in [10, 12], the authors showed the occurrence of chaotic dynamics in a simple three-species food chain model with Holling type II functional response. The occurrence of chaos in basic Lotka-Volterra models of four competing species was studied in [16], and, in [14, 17], they investigated the complex dynamical behavior of a three-species food chain model with Beddington-DeAngelis functional response. Thus, in this paper, we consider the following food chain system with three species: 𝑑𝑥(𝑡)𝑑𝑡=𝑥(𝑡)(𝑎−𝑏𝑥(𝑡))−𝑐¹𝑑𝑦(𝑡)𝐹₁(𝑥(𝑡),𝑦(𝑡))𝑦(𝑡),𝑑𝑡=−𝑑¹𝑑𝑦𝑧((𝑡𝑡))+𝑐₂𝐹₁(𝑥(𝑡),𝑦(𝑡))𝑦(𝑡)−𝑐₃𝐹₂(𝑦(𝑡),𝑧(𝑡))𝑧(𝑡),𝑑𝑡=−𝑑²𝑧(𝑡)+𝑐₄𝐹₂(𝑦(𝑡),𝑧(𝑡))𝑧(𝑡),(1.1) where 𝑥(𝑡),𝑦(𝑡), and 𝑧(𝑡) are functions of time representing population densities of the prey and the mid-predator and the top predator, respectively, and all parameters are positive constants. The constant 𝑎 is the intrinsic growth rates of the prey population, 𝑏 is the coefficient of intraspecific competition, 𝑐₁ is the per capita rate of predation of the mid-predator, 𝑑₁ denotes the death rate of the mid-predator, 𝑐₂ is the rate of conversion of a consumed prey to a predator, 𝑐₃ is the per-capita rate of predation of the top predator, 𝑑₂ denotes the death rate of the top predator, 𝑐₄ is the rate of conversion of a consumed prey to a predator, and here the functions 𝐹₁ and 𝐹₂ represent functional responses.

The references mentioned above have though of the functional responses 𝐹₁ and 𝐹₂ as the same type. However, from biological point of view, it is unrealistic. In fact, in real world, predators of different species may feed on preys in different types of consumption ways. For example, consider crops, aphids, and lady beetles as prey, mid-predator, and top predator, respectively. In this case, it is natural to assume that the feeding type of aphids on crops is different from that of lady beetles on aphids. Thus, to describe this phenomenon, two different types of functional responses are needed. So, in this paper, we consider two types of functional responses, the Holling type II functional response (cf. [18]) and the Beddington-DeAngelis functional response (cf. [4]). The former, which can be written as 𝐹₁(𝑥,𝑦)=𝑥/(𝛼₁+𝑥), is adopted to describe the relationship between the prey and the mid-predator, and the latter, which can be written as 𝐹₂(𝑦,𝑧)=𝑦/(𝛼₂+𝑦+𝛽𝑧), is adopted to express the relationship between the mid-predator and the top predator. Thus, the following food chain system with the hybrid type functional responses is considered in this paper: 𝑑𝑥(𝑡)𝑑𝑡=𝑥(𝑡)(𝑎−𝑏𝑥(𝑡))−𝑐¹𝑥(𝑡)𝑦(𝑡)𝑑𝑦(𝑡)𝛼₁+𝑥(𝑡),𝑑𝑡=−𝑑¹𝑦(𝑡)+𝑐²𝑥(𝑡)𝑦(𝑡)𝛼₁+𝑥(𝑡)−𝑐³𝑦(𝑡)𝑧(𝑡)𝑑𝑧(𝑡)𝛼₂+𝑦(𝑡)+𝛽𝑧(𝑡),𝑑𝑡=−𝑑²𝑧(𝑡)+𝑐⁴𝑦(𝑡)𝑧(𝑡)𝛼₂+𝑦(𝑡)+𝛽𝑧(𝑡),(1.2) where 𝛼₁ and 𝛼₂ are the half-saturation constants and the term 𝛽 scales the impact of the predator interference.

The main object of this paper is to investigate the dynamic properties and behaviors of system (1.2). In this context, the paper is organized as follows. In Section 2, we show the dissipativeness of system (1.2) and find a necessary condition for the mid-predator to survive. The local stabilities of the equilibrium points of system (1.2) are examined in Section 3, and the conditions for persistence of system (1.2) are found out according to the existence of limit cycles in Section 4 and some numerical examples are given to substantiate our theoretical results. Moreover, in Section 5, we provide numerical evidence of the existence of chaotic phenomena by illustrating bifurcation diagrams of system (1.2) and by calculating the largest Lyapunov exponent. Finally, conclusions are given in Section 6.

2. Dissipativeness

Obviously, the right-hand sides of system (1.2) are continuous and have continuous partial derivatives on the state space ℝ³⁺={(𝑥,𝑦,𝑧)^𝑇∣𝑥≥0,𝑦≥0,𝑧≥0}. In fact, they are Lipschitzian on ℝ³⁺ and then the solution of system (1.2) with nonnegative initial condition exists and is unique, as the solution of system (1.2) initiating in the nonnegative octant is bounded. Moreover, from [19], it is easy to show that ℝ³⁺ is an invariant domain of system (1.2).

A system is said to be dissipative if all population initiating in ℝ³⁺ are uniformly limited by their environment [20]. Thus, the dissipativeness of system (1.2) is carried out in the following theorem.

Theorem 2.1. System (1.2) is dissipative.

Proof. From the first equation of system (1.2), we get 𝑑𝑥(𝑡)/𝑑𝑡≤𝑥(𝑡)(𝑎−𝑏𝑥(𝑡)). So, by comparison theorem, 𝑥(𝑡)≤𝑎/(𝑏+𝐶𝑒^−𝑎𝑡) for all 𝑡≥0, where 𝐶=(𝑎−𝑏𝑥₀)/𝑥₀, which implies that 𝑥(𝑡)≤𝑎/𝑏 for sufficiently large 𝑡. Define 𝑉(𝑡)=(𝑐₂/𝑐₁)𝑥(𝑡)+𝑦(𝑡)+(𝑐₃/𝑐₄)𝑧(𝑡). Then 𝑑𝑉(𝑡)/𝑑𝑡≤𝑎(𝑎+1)𝑐₂/𝑏𝑐₁−𝑚𝑉(𝑡), where 𝑚=min{1,𝑑₁,𝑑₂}. So, by comparison theorem, we obtain that 𝑉(𝑡)≤𝑎(𝑎+1)𝑐₂/𝑏𝑐₁𝑚−(𝑎(𝑎+1)𝑐₂/𝑏𝑐₁𝑚)𝑒^−𝑚𝑡 for 𝑡≥0. Thus, (𝑐₂/𝑐₁)𝑥(𝑡)+𝑦(𝑡)+(𝑐₃/𝑐₄)𝑧(𝑡) ≤𝑎(𝑎+1)𝑐₂/𝑏𝑐₁𝑚 for sufficiently large 𝑡, which means that all species are uniformly bounded for any initial value in ℝ³⁺. Therefore, system (1.2) is dissipative.

The following proposition provides a necessary condition for survival of the mid-predator in system (1.2).

Proposition 2.2. A necessary condition for the mid-predator species 𝑦 to survive is 𝑑₁<𝑎𝑐²𝑏𝛼₁+𝑎.(2.1)

Proof. From the second equation of system (1.2), we get 𝑑𝑦(𝑡)𝑑𝑡=−𝑑¹𝑦(𝑡)+𝑐²𝑥(𝑡)𝑦(𝑡)𝛼₁+𝑥(𝑡)−𝑐³𝑦(𝑡)𝑧(𝑡)≤−𝑑₁𝛼₂+𝑦(𝑡)+𝛽𝑧(𝑡)𝑦(𝑡)+𝑐²𝑥(𝑡)𝑦(𝑡)≤𝑦(𝑡)−𝑑₁𝛼₁+𝑥(𝑡)+𝑎𝑐²𝑏𝛼₁+𝑎(ByTheorem2.1).(2.2) Then we have 𝑦(𝑡)≤𝑦₀𝑒^𝐴𝑡, where 𝐴=−𝑑₁+𝑎𝑐₂/(𝑏𝛼₁+𝑎). Thus, for 𝐴<0, lim_𝑡→∞𝑦(𝑡)=0. Hence, (2.1) is a necessary condition for the survival of the mid-predator 𝑦.

3. Stability Analysis

In order to investigate the stability of the equilibrium points of system (1.2), first, we consider the following two-dimensional dynamical system: 𝑑𝑥(𝑡)𝑑𝑡=𝑥(𝑡)(𝑎−𝑏𝑥(𝑡))−𝑐¹𝑥(𝑡)𝑦(𝑡)𝑑𝑦(𝑡)𝛼₁+𝑥(𝑡),𝑑𝑡=−𝑑¹𝑦(𝑡)+𝑐²𝑥(𝑡)𝑦(𝑡)𝛼₁+𝑥(𝑡).(3.1)

It is well known that the Kolmogorov theorem is applicable in two-dimensional dynamical system and guarantees the existence of either a stable equilibrium point or stable limit cycle behavior in the positive quadrant of phase space of the system, provided certain conditions are satisfied (cf. [20, 21]). Such conditions ensure that the parametric values are biologically relevant.

Now, it is observed that subsystem (3.1) is a Kolmogorov system under the condition 0<𝛼¹𝑑₁𝑐₂−𝑑₁<𝑎𝑏.(3.2) From now on, we assume that subsystem (3.1) satisfies the condition (3.2). By applying the local stability analysis to a Kolmogorov system (3.1) we have the following results [11].(I)The equilibrium point 𝐸₀₀=(0,0) always exists and is a saddle point.(II)The equilibrium point 𝐸₀₁=(𝑎/𝑏,0) always exists and is a saddle point.(III)The positive equilibrium point 𝐸₀₂=(𝑥̃,𝑦)̃ exists, where 𝑥̃=𝛼¹𝑑₁𝑐₂−𝑑₁,𝑦̃=(𝑎−𝑏𝑥̃)𝛼¹+𝑥̃𝑐₁,(3.3) and it is a locally asymptotically stable point if the following condition holds: 𝑑₁>𝑐²𝑎−𝑏𝛼₁𝑎+𝑏𝛼₁.(3.4) Moreover, the solution to system (3.1) approaches to a stable limit cycle if 𝑑₁<𝑐₂(𝑎−𝑏𝛼₁)/(𝑎+𝑏𝛼₁).

Now, we will study the dynamic behavior of the solution of system (1.2). First, we think over the stability of equilibrium points of system (1.2). In fact, there are at most four nonnegative equilibrium points of system (1.2). The existence conditions of them are mentioned as follows.(I)The trivial equilibrium point 𝐸₀=(0,0,0) and one species equilibrium point 𝐸₁=(𝑎/𝑏,0,0) always exist. However, the predators die out in the absence of the prey. Thus the equilibrium points (0,𝑦_𝑐,0) and (0,0,𝑧_𝑐) with 𝑦_𝑐,𝑧_𝑐>0 do not exist.(II)Two-species equilibrium point 𝐸₂=(𝑥̃,𝑦,̃0)exists in the interior of positive quadrant of 𝑥𝑦 plane under the Kolmogorov condition (3.2), where 𝑥̃ and 𝑦̃ are given in (3.3). On the other hand, the absence of the mid-predator causes no equilibrium point in the 𝑥𝑧 plane. Moreover, if there exists no prey, then neither 𝑦 nor 𝑧 can survive, which means that there is no equilibrium point in the 𝑦𝑧 plane.(III)The positive equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗) exists in the interior of the first octant if and only if there exists a positive solution to the following algebraic nonlinear simultaneous equations: 𝑓₁(𝑥,𝑦,𝑧)=𝑎−𝑏𝑥−𝑐¹𝑦𝑓𝛼₁+𝑥=0,₂(𝑥,𝑦,𝑧)=−𝑑₁+𝑐²𝑥𝛼₁+𝑥−𝑐³𝑧𝑓𝛼₂+𝑦+𝛽𝑧=0,₃(𝑥,𝑦,𝑧)=−𝑑₂+𝑐⁴𝑦𝛼₂+𝑦+𝛽𝑧=0.(3.5) By applying elementary calculation to (3.5), we obtain that 𝑦^∗=𝑎−𝑏𝑥^∗𝛼₁+𝑥^∗𝑐₁,𝑧^∗=𝑐⁴−𝑑₂𝑦^∗−𝑑₂𝛼₂𝑑₂𝛽,(3.6) where 𝑥^∗ is a positive solution of the quadratic equation 𝑝𝑥²+𝑞𝑥+𝑟=0, 𝑝𝑞==𝑎−𝑐𝑏₂𝑐₂𝑐₄𝛽+𝑏𝑐₄𝑑₁𝛽+𝑏𝑐₃𝑐₄−𝑏𝑐₃𝑑₂,𝑟=−𝑎𝑐𝑐₄₄𝛽+𝑏𝑐₄𝑑₁𝛼₁𝛽−𝑎𝑐₄𝑑₁𝛽+𝑎𝑐₃𝑑₂+𝑏𝑐₃𝑐₄𝛼₁−𝑎𝑐₃𝑐₄−𝑏𝑐₃𝑑₂𝛼₁,𝑑₁𝛼₁𝛽+𝑐₁𝑐₃𝑑₂𝛼₂−𝑎𝑐₃𝑐₄𝛼₁+𝑎𝑐₃𝑑₂𝛼₁.(3.7) Therefore, sufficient conditions for the existence of the positive equilibrium point in the interior of the first octant are easily obtained as follows: 𝑞²−4𝑝𝑟≥0,0<𝑥^∗<𝑎𝑏,0<𝛼²𝑑₂𝑐₄−𝑑₂<𝑦^∗.(3.8)

Now, in order to investigate the stabilities of the equilibrium points, we consider the variational matrix 𝑉(𝑥,𝑦,𝑧) of system (1.2). Thus, we get the matrix𝑉(𝑥,𝑦,𝑧)=⎛⎜⎜⎝𝑣¹¹𝑣₁₂𝑣₁₃𝑣₂₁𝑣₂₂𝑣₂₃𝑣₃₁𝑣₃₂𝑣₃₃⎞⎟⎟⎠,(3.9) where𝑣₁₁=𝑥𝜕𝑓¹𝜕𝑥+𝑓¹=𝑎−2𝑏𝑥−𝑐¹𝛼₁𝑦𝛼₁+𝑥²,𝑣₁₂=𝑥𝜕𝑓¹𝜕𝑦=−𝑐¹𝑥𝛼₁+𝑥,𝑣¹³=𝑥𝜕𝑓¹𝑣𝜕𝑧=0,₂₁=𝑦𝜕𝑓²𝜕𝑥=𝑐²𝛼₁𝑦𝛼₁+𝑥²,𝑣₂₂=𝑦𝜕𝑓²𝜕𝑦+𝑓²=−𝑑₁+𝑐²𝑥𝛼₁+𝑥−𝑐³𝑧𝛼₂+𝛽𝑧𝛼₂+𝑦+𝛽𝑧²𝑣₂₃,=𝑦𝜕𝑓²𝜕𝑧=−𝑐³𝑦𝛼₂+𝑦𝛼₂+𝑦+𝛽𝑧²,𝑣₃₁=𝑧𝜕𝑓³𝜕𝑥=0,𝑣³²=𝑧𝜕𝑓³𝜕𝑦=𝑐⁴𝑧𝛼₂+𝛽𝑧𝛼₂+𝑦+𝛽𝑧²𝑣₃₃,=𝑧𝜕𝑓³𝜕𝑧+𝑓³=−𝑑₂+𝑐⁴𝑦𝛼₂+𝑦𝛼₂+𝑦+𝛽𝑧²,(3.10) and 𝑓₁,𝑓₂, and 𝑓₃ are in (3.5).

Using the variational matrix 𝑉(𝑥,𝑦,𝑧), the local stability of system (1.2) near the equilibrium points are obtained as follows.(I)The trivial equilibrium point 𝐸₀ is a hyperbolic saddle point. In fact, near 𝐸₀=(0,0,0) the prey population is increasing, while both of the predators populations are decreasing.(II)The equilibrium point 𝐸₁=(𝑎/𝑏,0,0) is locally stable if 𝑑₁>𝑎𝑐₂/(𝑏𝛼₁+𝑎). However, under the Kolmogorov condition (3.2), that is, 𝑑₁<𝑎𝑐₂/(𝑏𝛼₁+𝑎), the point 𝐸₁ is a saddle point with locally stable manifold in 𝑥𝑧 plane and with locally unstable manifold in 𝑦-direction.(III)Clearly, the equilibrium point 𝐸₂=(𝑥̃,𝑦,̃0) has the same stability behavior as 𝐸₀₂=(𝑥̃,𝑦)̃ in the interior of positive coordinate 𝑥𝑦 plane. However, the stability of the point 𝐸₂ is determined by the positive direction orthogonal to the 𝑥𝑦 plane, that is, 𝑧-direction, depending on whether the eigenvalue 𝜆̃₃=−𝑑₂+𝑐₄𝑦/̃(𝛼₂+𝑦)̃ is negative or positive, respectively.(IV)Let 𝑉^∗=(𝑣^∗𝑖,𝑗) be the variational matrix at the equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗). Then we have 𝑣^∗11=𝑥^∗(−𝑏+𝑐₁𝑦^∗/(𝛼₁+𝑥^∗)²), 𝑣^∗22=𝑐₃𝑦^∗𝑧^∗/(𝛼₂+𝑦^∗+𝛽𝑧^∗)² (>0) and 𝑣^∗33= −𝑐₄𝛽𝑦^∗𝑧^∗/(𝛼₂+𝑦^∗+𝛽𝑧^∗)² (<0) because 𝑓₁(𝑥^∗,𝑦^∗,𝑧^∗)=𝑓₂(𝑥^∗,𝑦^∗,𝑧^∗)=𝑓₃(𝑥^∗,𝑦^∗,𝑧^∗)=0 and 𝑣^∗12<0, 𝑣^∗21>0, 𝑣^∗23<0, and 𝑣^∗32>0. Moreover, the characteristic equation of 𝑉^∗ is 𝜆³+𝐴𝜆²+𝐵𝜆+𝐶=0, where 𝐴𝐵==𝑣−_∗11𝑣^∗11+𝑣^∗22+𝑣^∗33,𝐶=𝑣𝑣^∗33_∗12+𝑣^∗22𝑣^∗33+𝑣^∗11𝑣^∗22−𝑣^∗12𝑣^∗21−𝑣^∗23𝑣^∗32,𝑣^∗21−𝑣^∗11𝑣^∗22𝑣^∗33+𝑣^∗11𝑣^∗23𝑣^∗32.(3.11) By the Routh-Hurwitz criterion [20], 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗) is locally asymptotically stable if and only if 𝐴, 𝐶, and 𝐴𝐵−𝐶 are positive.

A sufficient condition for the local stability of 𝐸₃ is given in the following theorem.

Theorem 3.1. Suppose that the positive equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗) exists in the interior of the positive octant. Then 𝐸₃ is locally asymptotically stable if the following conditions hold: 𝑏𝑥^∗>𝑦^∗𝑐¹𝑥^∗𝜑²+𝑐³𝑧^∗𝑏𝑐₃𝜙²,(3.12)𝑧^∗𝜑³<𝑐𝑐₁₄𝑐₂𝛼₁𝜙²+𝑐₃𝑦^∗𝑧^∗𝜑,𝛽−𝑐₃>0,(3.13) where 𝜑=𝛼₁+𝑥^∗ and 𝜙=𝛼₂+𝑦^∗+𝛽𝑧^∗.

Proof. Under the condition (3.12), it is easy to show that 𝐴,𝐶>0. By expanding 𝐴𝐵 and calculating 𝐴𝐵−𝐶, we get 𝐴𝐵−𝐶=𝑣^∗11+𝑣^∗22𝑣^∗12𝑣^∗21−𝑣^∗11𝑣^∗22+𝐴𝑣^∗33+𝑣^∗23𝑣^∗32𝑣^∗22+𝑣^∗33.(3.14) From the conditions (3.13), it is proven that 𝑣^∗12𝑣^∗21−𝑣^∗11𝑣^∗22<0 and 𝑣^∗22+𝑣^∗33<0. Therefore, from the signs of variational matrix elements 𝑣^∗𝑖,𝑗(𝑖,𝑗=1,2,3), we get 𝐴𝐵−𝐶>0. Hence, 𝐸₃ is locally asymptotically stable.

In the following theorem, the global stability of the equilibrium point 𝐸₃ is investigated.

Theorem 3.2. Suppose that the positive equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗) is locally asymptotically stable. Then it is a globally asymptotically stable if the following condition is satisfied: 𝑐₂𝑎+𝑏𝑥^∗²4𝑏𝑐₁+𝑑₁𝑦^∗+𝑐³𝑑₂𝑐₄𝑧^∗+𝑀𝑐²𝑥^∗𝛼₁+𝑐³𝑦^∗<𝑎𝑐₂𝛼₂𝑐₁𝑥^∗+𝑀𝑐²𝑦^∗𝛼₁+𝑀+𝑐³𝑧^∗𝛼₂+(1+𝛽)𝑀,(3.15) where 𝑀=𝑎(𝑎+1)𝑐₂/𝑏𝑐₁𝑚 and 𝑚=min{1,𝑑₁,𝑑₂}.

Proof. The proof can be reached by using a Lyapunov stability theorem which gives a sufficient condition. Now, let us consider a positive definite function 𝑊(𝑥,𝑦,𝑧)=(𝑐₂/𝑐₁)𝑊₁(𝑥,𝑦,𝑧)+𝑊₂(𝑥,𝑦,𝑧)+(𝑐₃/𝑐₄)𝑊₃(𝑥,𝑦,𝑧) in the interior of the positive octant, where 𝑊₁(𝑥,𝑦,𝑧)=𝑥−𝑥^∗ln𝑥𝑊₂𝑥^∗,(𝑥,𝑦,𝑧)=𝑦−𝑦^∗ln𝑦𝑊₃𝑦^∗,(𝑥,𝑦,𝑧)=𝑧−𝑧^∗ln𝑧𝑧^∗.(3.16) Note that 𝑑𝑊₁𝑑𝑡=𝑎−𝑏𝑥−𝑐¹𝑦𝑑𝑊₂𝛼₁+𝑥𝑥−𝑥^∗,𝑑𝑡=−𝑑¹+𝑐²𝑥𝛼₁+𝑥−𝑐³𝑧𝑑𝑊₃𝛼₂+𝑦+𝛽𝑧𝑦−𝑦^∗,𝑑𝑡=−𝑑²+𝑐⁴𝑦𝛼₂+𝑦+𝛽𝑧𝑧−𝑧^∗.(3.17) Thanks to Theorem 2.1, without loss of generality, we may assume that there exists a constant 𝑀=𝑎(𝑎+1)𝑐₂/𝑏𝑐₁𝑚 satisfying 𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)<𝑀, where 𝑚=min{1,𝑑₁,𝑑₂}. Thus, we have 𝑑𝑊𝑑𝑡=𝑐𝑐²₁𝑑𝑊₁𝑑𝑡+𝑑𝑑𝑊𝑡²+𝑐𝑐³₄𝑑𝑊₃≤𝑐𝑑𝑡₂𝑐₁𝑎𝑥−𝑏𝑥²−𝑎𝑥^∗+𝑏𝑥𝑥^∗+𝑐²𝑥^∗𝑀𝛼₁−𝑑₁𝑦+𝑑₁𝑦^∗−𝑐²𝑦^∗𝑀𝛼₁+𝑀+𝑐³𝑦^∗𝑀𝛼₂+𝑐³𝑐₄−𝑑₂𝑧+𝑑₂𝑧^∗−𝑐³𝑧^∗𝑀≤−𝑏𝑐₂𝛼₂+(1+𝛽)𝑀𝑐₁𝑥−𝑎+2𝑏𝑏𝑥^∗²−𝑑₁𝑦−𝑐³𝑑₂𝑐₄𝑧+𝑀𝑐²𝑥^∗𝛼₁−𝑐²𝑦^∗𝛼₁+𝑀+𝑐³𝑦^∗𝛼₂−𝑐³𝑧^∗+𝑐𝛼₂+(1+𝛽)𝑀₂𝑎+𝑏𝑥^∗²4𝑏𝑐₁−𝑎𝑐²𝑐₁𝑥^∗+𝑑₁𝑦^∗+𝑐³𝑑₂𝑐₄𝑧^∗.(3.18) It is easy to see that 𝑑𝑊/𝑑𝑡<0 under the condition (3.15). Therefore, the function 𝑊 in the interior of the positive octant is a Lyapunov function with respect to 𝐸₃. Hence, the equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗) is globally asymptotically stable.

Example 3.3. Let 𝑎=1.5, 𝑏=1, 𝑐₁=0.6, 𝑐₂=0.6, 𝑐₃=0.7, 𝑐₄=0.7, 𝑑₁=0.1, 𝑑₂=0.45, 𝛼₁=1, 𝛼₂=1.1, and 𝛽=1.1. Then it follows from (3.6) and Theorem 3.1 that the unique positive equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗)≈(0.3435,2.5863,0.3079) is locally stable as shown in Figure 1. However, even if these parameters do not satisfy the conditions of Theorem 3.2, Figure 2 exhibits that the positive equilibrium point 𝐸₃ seems to be globally stable.

Figure 1 

𝑎=1.5, 𝑏=1, 𝑐1=0.6, 𝑐2=0.6, 𝑐3=0.7, 𝑐4=0.7, 𝑑1=0.1, 𝑑2=0.45, 𝛼1=1, 𝛼2=1.1, and 𝛽=1.1 and an initial condition (0.3,2.5,0.3). (a) The trajectory of system (1.2). (b)–(d) Time series for 𝑥(𝑡),𝑦(𝑡), and 𝑧(𝑡).

Figure 2 

𝑎=1.5, 𝑏=1, 𝑐1=0.6, 𝑐2=0.6, 𝑐3=0.7, 𝑐4=0.7, 𝑑1=0.1, 𝑑2=0.45, 𝛼1=1, 𝛼2=1.1, and 𝛽=1.1 and an initial condition (100,100,100). (a) The trajectory of system (1.2). (b)–(d) Time series for 𝑥(𝑡),𝑦(𝑡), and 𝑧(𝑡).

Example 3.4. Consider system (1.2) with the parameters 𝑎=1, 𝑏=0.4, 𝑐₁=0.7, 𝑐₂=0.6, 𝑐₃=0.01, 𝑐₄=1.4, 𝑑₁=0.3, 𝑑₂=0.45, 𝛼₁=1.3, 𝛼₂=0.6, and 𝛽=1.0. From Theorem 3.2, the positive equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗)≈(1.3502,1.7413,3.0762) is globally stable as shown in Figure 3.

Figure 3 

𝑎=1, 𝑏=0.4, 𝑐1=0.7, 𝑐2=0.6, 𝑐3=0.01, 𝑐4=1.4, 𝑑1=0.3, 𝑑2=0.45, 𝛼1=1.3, 𝛼2=0.6, and 𝛽=1.0 and an initial condition (1,1,1). (a) The trajectory of the system (1.2). (b)–(d) Time series for 𝑥(𝑡),𝑦(𝑡), and 𝑧(𝑡).

4. Persistence of System (1.2)

The term persistence is given to systems in which strict solutions do not approach the boundary of the nonnegative cones as time goes to infinity. Therefore, for the continuous biological system, survival of all interacting species and the persistence are equivalent. In the following, we will find out some persistence conditions of the food chain system (1.2) with two-type functional responses.

Theorem 4.1. Suppose that there are no nontrivial periodic solutions in the 𝑥𝑦 plane. Then the necessary condition for the persistence of system (1.2) is 𝜆̃₃≥0(4.1) and the sufficient condition for the persistence of system (1.2) is 𝜆̃₃>0,(4.2) where 𝜆̃₃=−𝑑₂+𝑐₄𝑦/̃(𝛼₂+𝑦)̃= (−𝛼₂𝑑₁𝑐₁(𝑐₂−𝑑₁)²+𝛼₁𝑐₂(𝑐₄−𝑑₂)(𝑎𝑐₂−𝑎𝑑₁−𝑏𝛼₁𝑑₁))/(𝛼₂𝑐₁(𝑐₂−𝑑₁)²+𝛼₁𝑐₂(𝑎𝑐₂−𝑎𝑑₁−𝑏𝛼₁𝑑₁)).

Proof. Note that the boundedness of system (1.2) is shown in Theorem 2.1 and 𝜆̃₃ is the eigenvalue, which gives the stability of the equilibrium point 𝐸₂=(𝑥̃,𝑦,̃0) in the positive direction orthogonal to the 𝑥𝑦 plane, that is, 𝑧-direction. Since there are no nontrivial periodic solutions in the 𝑥𝑦 plane, so 𝐸₀₂ becomes a stable equilibrium point under the Kolmogorov condition (3.2). Therefore, if (4.1) does not hold (i.e., 𝜆̃₃<0), then there is an orbit in the positive cone, which approaches to 𝐸₂. Hence, the condition (4.1) is one of the necessary conditions for the persistence of system (1.2). Now, we will use the abstract theorem of Freedman and Waltman [20] for a sufficient condition of the persistence of system (1.2). For this, consider the growth functions 𝑓₁,𝑓₂, and 𝑓₃ of system (1.2) in (3.5). Then the following four conditions are satisfied:(C1)𝜕𝑓₁/𝜕𝑦<0, 𝜕𝑓₁/𝜕𝑧=0, 𝜕𝑓₂/𝜕𝑥>0, 𝜕𝑓₃/𝜕𝑥=0, 𝜕𝑓₃/𝜕𝑦>0, 𝑓₂(0,𝑦,𝑧)<0, and 𝑓₃(0,0,𝑧)<0,(C2)the prey population grows up to its carrying capacity 𝑏/𝑎 in the absence of predators, that is, 𝑓₁(0,0,0)=𝑎>0, 𝑓₁(𝑏/𝑎,0,0)=0, and (𝜕𝑓₁/𝜕𝑥)(𝑥,0,0)<0,(C3)there are no equilibrium points in the 𝑥𝑦 and 𝑦𝑧 planes,(C4)the mid-predator can survive on its prey in the absence of the top predator. In other words, there exists an equilibrium point 𝐸₂=(𝑥̃,𝑦,̃0) in the 𝑥𝑦 plane such that 𝑓₃(𝑥̃,𝑦,̃0)>0 by (4.2). However, the top predator cannot survive on the prey 𝑥 in the 𝑥𝑧 plane.Therefore, by Freedman and Waltman theorem [20], system (1.2) persists if condition (4.2) holds.

Theorem 4.2. Suppose that condition (4.2) holds and there are a finite number of limit cycles in 𝑥𝑦 plane. Then, for each limit cycle (𝑢(𝑡),𝑣(𝑡)) in 𝑥𝑦 plane, system (1.2) is persistent if the following condition holds: ^𝑇0𝑓₃(𝑢(𝑡),𝑣(𝑡),0)𝑑𝑡>0,(4.3) where 𝑇 is the time period of the limit cycle.

Proof. Note that the variational matrix 𝑉 about the limit cycle (𝑢(𝑡),𝑣(𝑡),0) is as follows: 𝑉(𝑢(𝑡),𝑣(𝑡),0)=⎛⎜⎜⎜⎜⎝𝑢(𝑡)𝜕𝑓¹𝜕𝑥+𝑓¹𝑢(𝑡)𝜕𝑓¹𝑣(𝑡)𝜕𝑓₂𝜕𝑦0𝜕𝑥𝑣(𝑡)𝜕𝜕𝑓𝑦²+𝑓²𝑣(𝑡)𝜕𝑓²00𝑓₃𝜕𝑧⎞⎟⎟⎟⎟⎠,(4.4) where all the partial derivatives and 𝑓_𝑖(𝑖=1,2,3) are computed at the limit cycle (𝑢(𝑡),𝑣(𝑡),0). Now, consider a solution of system (1.2) with positive initial condition (𝑢₀,𝑣₀,𝑤₀) sufficiently close to the limit cycle. It is from the variational matrix 𝑉(𝑢(𝑡),𝑣(𝑡),0) that 𝜕𝑧/𝜕𝑤₀ is a solution of system 𝑑𝑧/𝑑𝑡=𝑓₃(𝑢(𝑡),𝑣(𝑡),0)𝑧 with 𝑧(0)=1. So we have 𝜕𝑧𝜕𝑤₀𝑡,𝑢₀,𝑣₀,𝑤₀=exp^𝑡𝑜𝑓₃(𝑢(𝑠),𝑣(𝑠),0)𝑑𝑠.(4.5) Therefore, by using Taylor expansion, we have 𝑧𝑡,𝑢₀,𝑣₀,𝑤₀−𝑧𝑡,𝑢₀,𝑣₀,0≈exp^𝑡𝑜𝑓₃(𝑢(𝑠),𝑣(𝑠),0)𝑑𝑠𝑤₀.(4.6) Thus, the value of 𝑧 increases or decreases according to the sign of that of exp(∫^𝑡𝑜𝑓₃(𝑢(𝑠),𝑣(𝑠),0)𝑑𝑠). Since the equilibrium point 𝐸₃=(𝑥^∗,𝑦^∗,𝑧^∗) and these limit cycles are the only possible limit in the 𝑥𝑦 plane of trajectories with positive initial condition, these trajectories go away from the 𝑥𝑦 plane if the conditions (4.2) and (4.3) hold.

5. Bifurcation Analysis via Numerical Simulations

In this section, we will investigate various dynamic behaviors of system (1.2) by using numerical simulations. In particular, we will focus on exploring the possibility of chaotic behavior for the food chain with two different types of functional responses.

In fact, the dynamic behavior of system (1.2) depends on eleven independent parameters. So it is very difficult to study the system for complete range of parametric space. For the reason, two bifurcation parameters, the death rates 𝑑₁ and 𝑑₂ of the mid-predator and the top predator, respectively, are considered. In order to accomplish our purpose, fix the parameters except 𝑑₁ and 𝑑₂ as follows: 𝑐𝑎=1.5,𝑏=1,₁=𝛼1,𝑐₂=0.9,𝑐₃=1,𝑐₄=0.8,₁=0.5,𝛼₂=0.5,𝛽=0.1.(5.1)

First, fix 𝑑₂=0.2 and initial value (1,1,1). The bifurcation diagrams of system (1.2) for the successive maxima of the prey population 𝑥 and the mid-predator population 𝑦 and the top predator population 𝑧 are plotted in Figure 4 in the range of 0<𝑑₁<0.5.

(a)

(b)

(c)

Figure 4 

Bifurcation diagrams of system (1.2) with respect to 𝑑1: (a) maxima for the prey; (b) maxima for the mid-predator; (c) maxima for the top predator.

Figure 5 is the magnified parts of Figure 4, and the windows of periodic behaviors are more visible.

(a)

(b)

(c)

Figure 5 

Magnification of Figure 4 when 0.355<𝑑1<0.375.

The resulting bifurcation diagrams clearly show that system (1.2) has rich dynamics including stable limit cycles, period-doubling bifurcation, chaotic bands, periodic windows, and period-halving bifurcation. Especially, the values of 𝑑₁ between 0 and 0.2 result in a stable limit cycle, then a cascade of period-doubling bifurcation leads system (1.2) to a chaotic region (see Figure 6). And Figure 7 illustrates periodic halving phenomena.

(a)

(b)

(c)

(d)

Figure 6 

Period-doubling bifurcation: (a) attractor with period-1 when 𝑑1=0.15; (b) attractor with period-2 when 𝑑1=0.22; (c) attractor with period-2 when 𝑑1=0.25; (d) chaotic attractor when 𝑑1=0.33.

(a)

(b)

Figure 7 

Period-halving bifurcation: (a) attractor with period-2 when 𝑑1=0.42; (b) attractor with period-1 when 𝑑1=0.47.

Now, in order to investigate the dynamical behaviors with respect to the death rate 𝑑₂ of the top predator via bifurcation diagrams, fix the death rate 𝑑₁=0.35, and keep other parameters as given in (5.1). The bifurcation diagrams of system (1.2) for the successive maxima of the prey population 𝑥 and the mid-predator population 𝑦 and the top predator population 𝑧 as a function of 𝑑₂ are plotted in Figure 8.

(a)

(b)

(c)

Figure 8 

Bifurcation diagrams of system (1.2) with respect to 𝑑2: (a) maxima for the prey; (b) maxima for the mid-predator; (c) maxima for the top predator.

According to these bifurcation diagrams the dynamic behaviors of the solutions of system (1.2) also have different types of attracting sets including periodic and chaotic as 𝑑₂ varies between 0 and 0.43. However, a stable limit cycle can be shown in the range 𝑑₂>0.43. As is the case of 𝑑₁, these figures show the evidence of the route to chaos through the cascade of periodic doubling.

Thus, from the above facts, we can induce that the dynamic behavior of system (1.2) is extremely sensitive to the death rates 𝑑₁ and 𝑑₂. Moreover, in both cases, chaotic attractors can be easily observed as shown in Figure 6. In order to provide a quantitative measure of the degree of chaotic motion, the largest Lyapunov exponent is considered, which is the average exponential rates of divergence or convergence of nearby orbits in phase space (cf. [22]). In fact, a trajectory with the positive largest Lyapunov exponent is chaotic provided that it is not asymptotic to an unstable periodic solution. Or if the largest Lyapunov exponent of a trajectory is negative, then it is stable. Reviewing the bifurcation diagrams in Figures 4 and 8, the corresponding largest Lyapunov exponents (𝑑₁ from 0.24 to 0.45 and 𝑑₂ from 0.02 to 0.42) to system (1.2) are calculated in Figure 9 by using the method in [23].

(a)

(b)

Figure 9 

The largest Lyapunov exponent of system (1.2): (a) 0.24<𝑑1<0.45; (b) 0.02<𝑑2<0.42.

6. Conclusion

In this paper, we have proposed an ecological system with the hybrid type of functional responses, Holling type, and Beddington-DeAngelis type and have studied their dynamic behaviors. In the context, we have shown the dissipative nature of the system and have given sufficient conditions for the local and global stability of the equilibrium points of the system, as well as for permanence of the system. Furthermore, by illustrating bifurcation diagrams, we have shown that a food chain system with two kinds of functional responses may have a chaotic phenomenon. Based on these facts, we have concluded that three-species ecological models with hybrid functional responses can also have complex dynamical phenomena including chaotic behaviors.

Acknowledgments

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper. This work was supported by WCU (World Class University) program through the Korea Science and Engineering Foundation funded by the Ministry of Education, Science and Technology (Grant no. R32-2009-000-20021-0).