Table of Contents Author Guidelines Submit a Manuscript
Abstract and Applied Analysis
Volume 2011, Article ID 934569, 16 pages
http://dx.doi.org/10.1155/2011/934569
Research Article

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

1Department of Mathematics, Kyungpook National University, Daegu 702-701, Republic of Korea
2Department of Mathematics Education, Catholic University of Daegu, Gyeongsan-si, Gyeongbuk 712-702, Republic of Korea

Received 20 July 2010; Revised 11 December 2010; Accepted 19 January 2011

Academic Editor: H. B. Thompson

Copyright © 2011 Younghae Do et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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. [13]), 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. [917]) 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 ( 𝑥 ( 𝑡 ) , 𝑦 ( 𝑡 ) ) 𝑦 ( 𝑡 ) , 𝑑 𝑦 ( 𝑡 ) 𝑑 𝑡 = 𝑑 1 𝑦 ( 𝑡 ) + 𝑐 2 𝐹 1 ( 𝑥 ( 𝑡 ) , 𝑦 ( 𝑡 ) ) 𝑦 ( 𝑡 ) 𝑐 3 𝐹 2 ( 𝑦 ( 𝑡 ) , 𝑧 ( 𝑡 ) ) 𝑧 ( 𝑡 ) , 𝑑 𝑧 ( 𝑡 ) 𝑑 𝑡 = 𝑑 2 𝑧 ( 𝑡 ) + 𝑐 4 𝐹 2 ( 𝑦 ( 𝑡 ) , 𝑧 ( 𝑡 ) ) 𝑧 ( 𝑡 ) , ( 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, 𝑐 1 is the per capita rate of predation of the mid-predator, 𝑑 1 denotes the death rate of the mid-predator, 𝑐 2 is the rate of conversion of a consumed prey to a predator, 𝑐 3 is the per-capita rate of predation of the top predator, 𝑑 2 denotes the death rate of the top predator, 𝑐 4 is the rate of conversion of a consumed prey to a predator, and here the functions 𝐹 1 and 𝐹 2 represent functional responses.

The references mentioned above have though of the functional responses 𝐹 1 and 𝐹 2 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 𝐹 1 ( 𝑥 , 𝑦 ) = 𝑥 / ( 𝛼 1 + 𝑥 ) , is adopted to describe the relationship between the prey and the mid-predator, and the latter, which can be written as 𝐹 2 ( 𝑦 , 𝑧 ) = 𝑦 / ( 𝛼 2 + 𝑦 + 𝛽 𝑧 ) , 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 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝛼 1 , + 𝑥 ( 𝑡 ) 𝑑 𝑦 ( 𝑡 ) 𝑑 𝑡 = 𝑑 1 𝑐 𝑦 ( 𝑡 ) + 2 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝛼 1 𝑐 + 𝑥 ( 𝑡 ) 3 𝑦 ( 𝑡 ) 𝑧 ( 𝑡 ) 𝛼 2 , + 𝑦 ( 𝑡 ) + 𝛽 𝑧 ( 𝑡 ) 𝑑 𝑧 ( 𝑡 ) 𝑑 𝑡 = 𝑑 2 𝑧 𝑐 ( 𝑡 ) + 4 𝑦 ( 𝑡 ) 𝑧 ( 𝑡 ) 𝛼 2 , + 𝑦 ( 𝑡 ) + 𝛽 𝑧 ( 𝑡 ) ( 1 . 2 ) where 𝛼 1 and 𝛼 2 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 3 + = { ( 𝑥 , 𝑦 , 𝑧 ) 𝑇 𝑥 0 , 𝑦 0 , 𝑧 0 } . In fact, they are Lipschitzian on 3 + 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 3 + is an invariant domain of system (1.2).

A system is said to be dissipative if all population initiating in 3 + 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 𝐶 = ( 𝑎 𝑏 𝑥 0 ) / 𝑥 0 , which implies that 𝑥 ( 𝑡 ) 𝑎 / 𝑏 for sufficiently large 𝑡 . Define 𝑉 ( 𝑡 ) = ( 𝑐 2 / 𝑐 1 ) 𝑥 ( 𝑡 ) + 𝑦 ( 𝑡 ) + ( 𝑐 3 / 𝑐 4 ) 𝑧 ( 𝑡 ) . Then 𝑑 𝑉 ( 𝑡 ) / 𝑑 𝑡 𝑎 ( 𝑎 + 1 ) 𝑐 2 / 𝑏 𝑐 1 𝑚 𝑉 ( 𝑡 ) , where 𝑚 = m i n { 1 , 𝑑 1 , 𝑑 2 } . So, by comparison theorem, we obtain that 𝑉 ( 𝑡 ) 𝑎 ( 𝑎 + 1 ) 𝑐 2 / 𝑏 𝑐 1 𝑚 ( 𝑎 ( 𝑎 + 1 ) 𝑐 2 / 𝑏 𝑐 1 𝑚 ) 𝑒 𝑚 𝑡 for 𝑡 0 . Thus, ( 𝑐 2 / 𝑐 1 ) 𝑥 ( 𝑡 ) + 𝑦 ( 𝑡 ) + ( 𝑐 3 / 𝑐 4 ) 𝑧 ( 𝑡 ) 𝑎 ( 𝑎 + 1 ) 𝑐 2 / 𝑏 𝑐 1 𝑚 for sufficiently large 𝑡 , which means that all species are uniformly bounded for any initial value in 3 + . 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 𝑑 1 < 𝑎 𝑐 2 𝑏 𝛼 1 . + 𝑎 ( 2 . 1 )

Proof. From the second equation of system (1.2), we get 𝑑 𝑦 ( 𝑡 ) 𝑑 𝑡 = 𝑑 1 𝑐 𝑦 ( 𝑡 ) + 2 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝛼 1 𝑐 + 𝑥 ( 𝑡 ) 3 𝑦 ( 𝑡 ) 𝑧 ( 𝑡 ) 𝛼 2 + 𝑦 ( 𝑡 ) + 𝛽 𝑧 ( 𝑡 ) 𝑑 1 𝑐 𝑦 ( 𝑡 ) + 2 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝛼 1 + 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝑑 1 + 𝑎 𝑐 2 𝑏 𝛼 1 + 𝑎 ( B y T h e o r e m 2 . 1 ) . ( 2 . 2 ) Then we have 𝑦 ( 𝑡 ) 𝑦 0 𝑒 𝐴 𝑡 , where 𝐴 = 𝑑 1 + 𝑎 𝑐 2 / ( 𝑏 𝛼 1 + 𝑎 ) . Thus, for 𝐴 < 0 , l i m 𝑡 𝑦 ( 𝑡 ) = 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: 𝑑 𝑥 ( 𝑡 ) 𝑐 𝑑 𝑡 = 𝑥 ( 𝑡 ) ( 𝑎 𝑏 𝑥 ( 𝑡 ) ) 1 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝛼 1 , + 𝑥 ( 𝑡 ) 𝑑 𝑦 ( 𝑡 ) 𝑑 𝑡 = 𝑑 1 𝑐 𝑦 ( 𝑡 ) + 2 𝑥 ( 𝑡 ) 𝑦 ( 𝑡 ) 𝛼 1 . + 𝑥 ( 𝑡 ) ( 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 < 1 𝑑 1 𝑐 2 𝑑 1 < 𝑎 𝑏 . ( 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 = ( 0 , 0 ) always exists and is a saddle point.(II)The equilibrium point 𝐸 0 1 = ( 𝑎 / 𝑏 , 0 ) always exists and is a saddle point.(III)The positive equilibrium point 𝐸 0 2 = ( ̃ 𝑥 , ̃ 𝑦 ) exists, where 𝛼 ̃ 𝑥 = 1 𝑑 1 𝑐 2 𝑑 1 𝛼 , ̃ 𝑦 = ( 𝑎 𝑏 ̃ 𝑥 ) 1 + ̃ 𝑥 𝑐 1 , ( 3 . 3 ) and it is a locally asymptotically stable point if the following condition holds: 𝑑 1 > 𝑐 2 𝑎 𝑏 𝛼 1 𝑎 + 𝑏 𝛼 1 . ( 3 . 4 ) Moreover, the solution to system (3.1) approaches to a stable limit cycle if 𝑑 1 < 𝑐 2 ( 𝑎 𝑏 𝛼 1 ) / ( 𝑎 + 𝑏 𝛼 1 ) .

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 , 0 ) and one species equilibrium point 𝐸 1 = ( 𝑎 / 𝑏 , 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 𝐸 2 = ( ̃ 𝑥 , ̃ 𝑦 , 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 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) exists in the interior of the first octant if and only if there exists a positive solution to the following algebraic nonlinear simultaneous equations: 𝑓 1 𝑐 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑎 𝑏 𝑥 1 𝑦 𝛼 1 𝑓 + 𝑥 = 0 , 2 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑑 1 + 𝑐 2 𝑥 𝛼 1 𝑐 + 𝑥 3 𝑧 𝛼 2 𝑓 + 𝑦 + 𝛽 𝑧 = 0 , 3 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑑 2 + 𝑐 4 𝑦 𝛼 2 + 𝑦 + 𝛽 𝑧 = 0 . ( 3 . 5 ) By applying elementary calculation to (3.5), we obtain that 𝑦 = 𝑎 𝑏 𝑥 𝛼 1 + 𝑥 𝑐 1 , 𝑧 = 𝑐 4 𝑑 2 𝑦 𝑑 2 𝛼 2 𝑑 2 𝛽 , ( 3 . 6 ) where 𝑥 is a positive solution of the quadratic equation 𝑝 𝑥 2 + 𝑞 𝑥 + 𝑟 = 0 , 𝑝 = 𝑏 𝑐 2 𝑐 4 𝛽 + 𝑏 𝑐 4 𝑑 1 𝛽 + 𝑏 𝑐 3 𝑐 4 𝑏 𝑐 3 𝑑 2 , 𝑞 = 𝑎 𝑐 2 𝑐 4 𝛽 + 𝑏 𝑐 4 𝑑 1 𝛼 1 𝛽 𝑎 𝑐 4 𝑑 1 𝛽 + 𝑎 𝑐 3 𝑑 2 + 𝑏 𝑐 3 𝑐 4 𝛼 1 𝑎 𝑐 3 𝑐 4 𝑏 𝑐 3 𝑑 2 𝛼 1 , 𝑟 = 𝑎 𝑐 4 𝑑 1 𝛼 1 𝛽 + 𝑐 1 𝑐 3 𝑑 2 𝛼 2 𝑎 𝑐 3 𝑐 4 𝛼 1 + 𝑎 𝑐 3 𝑑 2 𝛼 1 . ( 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: 𝑞 2 4 𝑝 𝑟 0 , 0 < 𝑥 < 𝑎 𝑏 𝛼 , 0 < 2 𝑑 2 𝑐 4 𝑑 2 < 𝑦 . ( 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 𝑣 𝑉 ( 𝑥 , 𝑦 , 𝑧 ) = 1 1 𝑣 1 2 𝑣 1 3 𝑣 2 1 𝑣 2 2 𝑣 2 3 𝑣 3 1 𝑣 3 2 𝑣 3 3 , ( 3 . 9 ) where 𝑣 1 1 = 𝑥 𝜕 𝑓 1 𝜕 𝑥 + 𝑓 1 𝑐 = 𝑎 2 𝑏 𝑥 1 𝛼 1 𝑦 𝛼 1 + 𝑥 2 , 𝑣 1 2 = 𝑥 𝜕 𝑓 1 𝑐 𝜕 𝑦 = 1 𝑥 𝛼 1 + 𝑥 , 𝑣 1 3 = 𝑥 𝜕 𝑓 1 𝑣 𝜕 𝑧 = 0 , 2 1 = 𝑦 𝜕 𝑓 2 = 𝑐 𝜕 𝑥 2 𝛼 1 𝑦 𝛼 1 + 𝑥 2 , 𝑣 2 2 = 𝑦 𝜕 𝑓 2 𝜕 𝑦 + 𝑓 2 = 𝑑 1 + 𝑐 2 𝑥 𝛼 1 𝑐 + 𝑥 3 𝑧 𝛼 2 + 𝛽 𝑧 𝛼 2 + 𝑦 + 𝛽 𝑧 2 , 𝑣 2 3 = 𝑦 𝜕 𝑓 2 𝑐 𝜕 𝑧 = 3 𝑦 𝛼 2 + 𝑦 𝛼 2 + 𝑦 + 𝛽 𝑧 2 , 𝑣 3 1 = 𝑧 𝜕 𝑓 3 𝜕 𝑥 = 0 , 𝑣 3 2 = 𝑧 𝜕 𝑓 3 = 𝑐 𝜕 𝑦 4 𝑧 𝛼 2 + 𝛽 𝑧 𝛼 2 + 𝑦 + 𝛽 𝑧 2 , 𝑣 3 3 = 𝑧 𝜕 𝑓 3 𝜕 𝑧 + 𝑓 3 = 𝑑 2 + 𝑐 4 𝑦 𝛼 2 + 𝑦 𝛼 2 + 𝑦 + 𝛽 𝑧 2 , ( 3 . 1 0 ) and 𝑓 1 , 𝑓 2 , and 𝑓 3 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 𝐸 0 is a hyperbolic saddle point. In fact, near 𝐸 0 = ( 0 , 0 , 0 ) the prey population is increasing, while both of the predators populations are decreasing. ( I I ) The equilibrium point 𝐸 1 = ( 𝑎 / 𝑏 , 0 , 0 ) is locally stable if 𝑑 1 > 𝑎 𝑐 2 / ( 𝑏 𝛼 1 + 𝑎 ) . However, under the Kolmogorov condition (3.2), that is, 𝑑 1 < 𝑎 𝑐 2 / ( 𝑏 𝛼 1 + 𝑎 ) , the point 𝐸 1 is a saddle point with locally stable manifold in 𝑥 𝑧 plane and with locally unstable manifold in 𝑦 -direction. ( I I I ) Clearly, the equilibrium point 𝐸 2 = ( ̃ 𝑥 , ̃ 𝑦 , 0 ) has the same stability behavior as 𝐸 0 2 = ( ̃ 𝑥 , ̃ 𝑦 ) in the interior of positive coordinate 𝑥 𝑦 plane. However, the stability of the point 𝐸 2 is determined by the positive direction orthogonal to the 𝑥 𝑦 plane, that is, 𝑧 -direction, depending on whether the eigenvalue ̃ 𝜆 3 = 𝑑 2 + 𝑐 4 ̃ 𝑦 / ( 𝛼 2 + ̃ 𝑦 ) is negative or positive, respectively. ( I V ) Let 𝑉 = ( 𝑣 𝑖 , 𝑗 ) be the variational matrix at the equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) . Then we have 𝑣 1 1 = 𝑥 ( 𝑏 + 𝑐 1 𝑦 / ( 𝛼 1 + 𝑥 ) 2 ) , 𝑣 2 2 = 𝑐 3 𝑦 𝑧 / ( 𝛼 2 + 𝑦 + 𝛽 𝑧 ) 2 ( > 0 ) and 𝑣 3 3 = 𝑐 4 𝛽 𝑦 𝑧 / ( 𝛼 2 + 𝑦 + 𝛽 𝑧 ) 2 ( < 0 ) because 𝑓 1 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑓 2 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑓 3 ( 𝑥 , 𝑦 , 𝑧 ) = 0 and 𝑣 1 2 < 0 , 𝑣 2 1 > 0 , 𝑣 2 3 < 0 , and 𝑣 3 2 > 0 . Moreover, the characteristic equation of 𝑉 is 𝜆 3 + 𝐴 𝜆 2 + 𝐵 𝜆 + 𝐶 = 0 , where 𝑣 𝐴 = 1 1 + 𝑣 2 2 + 𝑣 3 3 , 𝐵 = 𝑣 1 1 𝑣 3 3 + 𝑣 2 2 𝑣 3 3 + 𝑣 1 1 𝑣 2 2 𝑣 1 2 𝑣 2 1 𝑣 2 3 𝑣 3 2 , 𝑣 𝐶 = 1 2 𝑣 2 1 𝑣 1 1 𝑣 2 2 𝑣 3 3 + 𝑣 1 1 𝑣 2 3 𝑣 3 2 . ( 3 . 1 1 ) By the Routh-Hurwitz criterion [20], 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) is locally asymptotically stable if and only if 𝐴 , 𝐶 , and 𝐴 𝐵 𝐶 are positive.

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

Theorem 3.1. Suppose that the positive equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) exists in the interior of the positive octant. Then 𝐸 3 is locally asymptotically stable if the following conditions hold: 𝑏 𝑥 > 𝑦 𝑐 1 𝑥 𝜑 2 + 𝑐 3 𝑧 𝜙 2 , ( 3 . 1 2 ) 𝑏 𝑐 3 𝑧 𝜑 3 < 𝑐 1 𝑐 2 𝛼 1 𝜙 2 + 𝑐 3 𝑦 𝑧 𝜑 , 𝑐 4 𝛽 𝑐 3 > 0 , ( 3 . 1 3 ) where 𝜑 = 𝛼 1 + 𝑥 and 𝜙 = 𝛼 2 + 𝑦 + 𝛽 𝑧 .

Proof. Under the condition (3.12), it is easy to show that 𝐴 , 𝐶 > 0 . By expanding 𝐴 𝐵 and calculating 𝐴 𝐵 𝐶 , we get 𝑣 𝐴 𝐵 𝐶 = 1 1 + 𝑣 2 2 𝑣 1 2 𝑣 2 1 𝑣 1 1 𝑣 2 2 + 𝐴 𝑣 3 3 + 𝑣 2 3 𝑣 3 2 𝑣 2 2 + 𝑣 3 3 . ( 3 . 1 4 ) From the conditions (3.13), it is proven that 𝑣 1 2 𝑣 2 1 𝑣 1 1 𝑣 2 2 < 0 and 𝑣 2 2 + 𝑣 3 3 < 0 . Therefore, from the signs of variational matrix elements 𝑣 𝑖 , 𝑗 ( 𝑖 , 𝑗 = 1 , 2 , 3 ) , we get 𝐴 𝐵 𝐶 > 0 . Hence, 𝐸 3 is locally asymptotically stable.

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

Theorem 3.2. Suppose that the positive equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) is locally asymptotically stable. Then it is a globally asymptotically stable if the following condition is satisfied: 𝑐 2 𝑎 + 𝑏 𝑥 2 4 𝑏 𝑐 1 + 𝑑 1 𝑦 + 𝑐 3 𝑑 2 𝑐 4 𝑧 𝑐 + 𝑀 2 𝑥 𝛼 1 + 𝑐 3 𝑦 𝛼 2 < 𝑎 𝑐 2 𝑐 1 𝑥 𝑐 + 𝑀 2 𝑦 𝛼 1 + 𝑐 + 𝑀 3 𝑧 𝛼 2 , + ( 1 + 𝛽 ) 𝑀 ( 3 . 1 5 ) where 𝑀 = 𝑎 ( 𝑎 + 1 ) 𝑐 2 / 𝑏 𝑐 1 𝑚 and 𝑚 = m i n { 1 , 𝑑 1 , 𝑑 2 } .

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 𝑊 ( 𝑥 , 𝑦 , 𝑧 ) = ( 𝑐 2 / 𝑐 1 ) 𝑊 1 ( 𝑥 , 𝑦 , 𝑧 ) + 𝑊 2 ( 𝑥 , 𝑦 , 𝑧 ) + ( 𝑐 3 / 𝑐 4 ) 𝑊 3 ( 𝑥 , 𝑦 , 𝑧 ) in the interior of the positive octant, where 𝑊 1 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑥 𝑥 𝑥 l n 𝑥 , 𝑊 2 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑦 𝑦 𝑦 l n 𝑦 , 𝑊 3 ( 𝑥 , 𝑦 , 𝑧 ) = 𝑧 𝑧 𝑧 l n 𝑧 . ( 3 . 1 6 ) Note that 𝑑 𝑊 1 = 𝑐 𝑑 𝑡 𝑎 𝑏 𝑥 1 𝑦 𝛼 1 + 𝑥 𝑥 𝑥 , 𝑑 𝑊 2 = 𝑑 𝑡 𝑑 1 + 𝑐 2 𝑥 𝛼 1 𝑐 + 𝑥 3 𝑧 𝛼 2 + 𝑦 + 𝛽 𝑧 𝑦 𝑦 , 𝑑 𝑊 3 = 𝑑 𝑡 𝑑 2 + 𝑐 4 𝑦 𝛼 2 + 𝑦 + 𝛽 𝑧 𝑧 𝑧 . ( 3 . 1 7 ) Thanks to Theorem 2.1, without loss of generality, we may assume that there exists a constant 𝑀 = 𝑎 ( 𝑎 + 1 ) 𝑐 2 / 𝑏 𝑐 1 𝑚 satisfying 𝑥 ( 𝑡 ) , 𝑦 ( 𝑡 ) , 𝑧 ( 𝑡 ) < 𝑀 , where 𝑚 = m i n { 1 , 𝑑 1 , 𝑑 2 } . Thus, we have 𝑑 𝑊 = 𝑐 𝑑 𝑡 2 𝑐 1 𝑑 𝑊 1 + 𝑑 𝑡 𝑑 𝑊 2 + 𝑐 𝑑 𝑡 3 𝑐 4 𝑑 𝑊 3 𝑐 𝑑 𝑡 2 𝑐 1 𝑎 𝑥 𝑏 𝑥 2 𝑎 𝑥 + 𝑏 𝑥 𝑥 + 𝑐 2 𝑥 𝑀 𝛼 1 𝑑 1 𝑦 + 𝑑 1 𝑦 𝑐 2 𝑦 𝑀 𝛼 1 + 𝑐 + 𝑀 3 𝑦 𝑀 𝛼 2 + 𝑐 3 𝑐 4 𝑑 2 𝑧 + 𝑑 2 𝑧 𝑐 3 𝑧 𝑀 𝛼 2 + ( 1 + 𝛽 ) 𝑀 𝑏 𝑐 2 𝑐 1 𝑥 𝑎 + 𝑏 𝑥 2 𝑏 2 𝑑 1 𝑐 𝑦 3 𝑑 2 𝑐 4 𝑐 𝑧 + 𝑀 2 𝑥 𝛼 1 𝑐 2 𝑦 𝛼 1 + 𝑐 + 𝑀 3 𝑦 𝛼 2 𝑐 3 𝑧 𝛼 2 + 𝑐 + ( 1 + 𝛽 ) 𝑀 2 𝑎 + 𝑏 𝑥 2 4 𝑏 𝑐 1 𝑎 𝑐 2 𝑐 1 𝑥 + 𝑑 1 𝑦 + 𝑐 3 𝑑 2 𝑐 4 𝑧 . ( 3 . 1 8 ) 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 𝐸 3 . Hence, the equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) is globally asymptotically stable.

Example 3.3. Let 𝑎 = 1 . 5 , 𝑏 = 1 , 𝑐 1 = 0 . 6 , 𝑐 2 = 0 . 6 , 𝑐 3 = 0 . 7 , 𝑐 4 = 0 . 7 , 𝑑 1 = 0 . 1 , 𝑑 2 = 0 . 4 5 , 𝛼 1 = 1 , 𝛼 2 = 1 . 1 , and 𝛽 = 1 . 1 . Then it follows from (3.6) and Theorem 3.1 that the unique positive equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) ( 0 . 3 4 3 5 , 2 . 5 8 6 3 , 0 . 3 0 7 9 ) 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 𝐸 3 seems to be globally stable.

fig1
Figure 1: 𝑎 = 1 . 5 , 𝑏 = 1 , 𝑐 1 = 0 . 6 , 𝑐 2 = 0 . 6 , 𝑐 3 = 0 . 7 , 𝑐 4 = 0 . 7 , 𝑑 1 = 0 . 1 , 𝑑 2 = 0 . 4 5 , 𝛼 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 𝑧 ( 𝑡 ) .
fig2
Figure 2: 𝑎 = 1 . 5 , 𝑏 = 1 , 𝑐 1 = 0 . 6 , 𝑐 2 = 0 . 6 , 𝑐 3 = 0 . 7 , 𝑐 4 = 0 . 7 , 𝑑 1 = 0 . 1 , 𝑑 2 = 0 . 4 5 , 𝛼 1 = 1 , 𝛼 2 = 1 . 1 , and 𝛽 = 1 . 1 and an initial condition ( 1 0 0 , 1 0 0 , 1 0 0 ) . (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 , 𝑐 1 = 0 . 7 , 𝑐 2 = 0 . 6 , 𝑐 3 = 0 . 0 1 , 𝑐 4 = 1 . 4 , 𝑑 1 = 0 . 3 , 𝑑 2 = 0 . 4 5 , 𝛼 1 = 1 . 3 , 𝛼 2 = 0 . 6 , and 𝛽 = 1 . 0 . From Theorem 3.2, the positive equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) ( 1 . 3 5 0 2 , 1 . 7 4 1 3 , 3 . 0 7 6 2 ) is globally stable as shown in Figure 3.

fig3
Figure 3: 𝑎 = 1 , 𝑏 = 0 . 4 , 𝑐 1 = 0 . 7 , 𝑐 2 = 0 . 6 , 𝑐 3 = 0 . 0 1 , 𝑐 4 = 1 . 4 , 𝑑 1 = 0 . 3 , 𝑑 2 = 0 . 4 5 , 𝛼 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 ̃ 𝜆 3 0 ( 4 . 1 ) and the sufficient condition for the persistence of system (1.2) is ̃ 𝜆 3 > 0 , ( 4 . 2 ) where ̃ 𝜆 3 = 𝑑 2 + 𝑐 4 ̃ 𝑦 / ( 𝛼 2 + ̃ 𝑦 ) = ( 𝛼 2 𝑑 1 𝑐 1 ( 𝑐 2 𝑑 1 ) 2 + 𝛼 1 𝑐 2 ( 𝑐 4 𝑑 2 ) ( 𝑎 𝑐 2 𝑎 𝑑 1 𝑏 𝛼 1 𝑑 1 ) ) / ( 𝛼 2 𝑐 1 ( 𝑐 2 𝑑 1 ) 2 + 𝛼 1 𝑐 2 ( 𝑎 𝑐 2 𝑎 𝑑 1 𝑏 𝛼 1 𝑑 1 ) ) .

Proof. Note that the boundedness of system (1.2) is shown in Theorem 2.1 and ̃ 𝜆 3 is the eigenvalue, which gives the stability of the equilibrium point 𝐸 2 = ( ̃ 𝑥 , ̃ 𝑦 , 0 ) in the positive direction orthogonal to the 𝑥 𝑦 plane, that is, 𝑧 -direction. Since there are no nontrivial periodic solutions in the 𝑥 𝑦 plane, so 𝐸 0 2 becomes a stable equilibrium point under the Kolmogorov condition (3.2). Therefore, if (4.1) does not hold (i.e., ̃ 𝜆 3 < 0 ), then there is an orbit in the positive cone, which approaches to 𝐸 2 . 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 𝑓 1 , 𝑓 2 , and 𝑓 3 of system (1.2) in (3.5). Then the following four conditions are satisfied: ( C 1 ) 𝜕 𝑓 1 / 𝜕 𝑦 < 0 , 𝜕 𝑓 1 / 𝜕 𝑧 = 0 , 𝜕 𝑓 2 / 𝜕 𝑥 > 0 , 𝜕 𝑓 3 / 𝜕 𝑥 = 0 , 𝜕 𝑓 3 / 𝜕 𝑦 > 0 , 𝑓 2 ( 0 , 𝑦 , 𝑧 ) < 0 , and 𝑓 3 ( 0 , 0 , 𝑧 ) < 0 , ( C 2 ) the prey population grows up to its carrying capacity 𝑏 / 𝑎 in the absence of predators, that is, 𝑓 1 ( 0 , 0 , 0 ) = 𝑎 > 0 , 𝑓 1 ( 𝑏 / 𝑎 , 0 , 0 ) = 0 , and ( 𝜕 𝑓 1 / 𝜕 𝑥 ) ( 𝑥 , 0 , 0 ) < 0 , ( C 3 ) there are no equilibrium points in the 𝑥 𝑦 and 𝑦 𝑧 planes, ( C 4 ) the mid-predator can survive on its prey in the absence of the top predator. In other words, there exists an equilibrium point 𝐸 2 = ( ̃ 𝑥 , ̃ 𝑦 , 0 ) in the 𝑥 𝑦 plane such that 𝑓 3 ( ̃ 𝑥 , ̃ 𝑦 , 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 𝑓 3 ( 𝑢 ( 𝑡 ) , 𝑣 ( 𝑡 ) , 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 ) = 𝑢 ( 𝑡 ) 𝜕 𝑓 1 𝜕 𝑥 + 𝑓 1 𝑢 ( 𝑡 ) 𝜕 𝑓 1 0 𝜕 𝑦 𝑣 ( 𝑡 ) 𝜕 𝑓 2 𝜕 𝑥 𝑣 ( 𝑡 ) 𝜕 𝑓 2 𝜕 𝑦 + 𝑓 2 𝑣 ( 𝑡 ) 𝜕 𝑓 2 𝜕 𝑧 0 0 𝑓 3 , ( 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 ( 𝑢 0 , 𝑣 0 , 𝑤 0 ) sufficiently close to the limit cycle. It is from the variational matrix 𝑉 ( 𝑢 ( 𝑡 ) , 𝑣 ( 𝑡 ) , 0 ) that 𝜕 𝑧 / 𝜕 𝑤 0 is a solution of system 𝑑 𝑧 / 𝑑 𝑡 = 𝑓 3 ( 𝑢 ( 𝑡 ) , 𝑣 ( 𝑡 ) , 0 ) 𝑧 with 𝑧 ( 0 ) = 1 . So we have 𝜕 𝑧 𝜕 𝑤 0 𝑡 , 𝑢 0 , 𝑣 0 , 𝑤 0 = e x p 𝑡 𝑜 𝑓 3 ( . 𝑢 ( 𝑠 ) , 𝑣 ( 𝑠 ) , 0 ) 𝑑 𝑠 ( 4 . 5 ) Therefore, by using Taylor expansion, we have 𝑧 𝑡 , 𝑢 0 , 𝑣 0 , 𝑤 0 𝑧 𝑡 , 𝑢 0 , 𝑣 0 , 0 e x p 𝑡 𝑜 𝑓 3 ( 𝑤 𝑢 ( 𝑠 ) , 𝑣 ( 𝑠 ) , 0 ) 𝑑 𝑠 0 . ( 4 . 6 ) Thus, the value of 𝑧 increases or decreases according to the sign of that of e x p ( 𝑡 𝑜 𝑓 3 ( 𝑢 ( 𝑠 ) , 𝑣 ( 𝑠 ) , 0 ) 𝑑 𝑠 ) . Since the equilibrium point 𝐸 3 = ( 𝑥 , 𝑦 , 𝑧 ) 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 𝑑 1 and 𝑑 2 of the mid-predator and the top predator, respectively, are considered. In order to accomplish our purpose, fix the parameters except 𝑑 1 and 𝑑 2 as follows: 𝑐 𝑎 = 1 . 5 , 𝑏 = 1 , 1 = 1 , 𝑐 2 = 0 . 9 , 𝑐 3 = 1 , 𝑐 4 𝛼 = 0 . 8 , 1 = 0 . 5 , 𝛼 2 = 0 . 5 , 𝛽 = 0 . 1 . ( 5 . 1 )

First, fix 𝑑 2 = 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 < 𝑑 1 < 0 . 5 .

fig4
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.

fig5
Figure 5: Magnification of Figure 4 when 0 . 3 5 5 < 𝑑 1 < 0 . 3 7 5 .

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 𝑑 1 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.

fig6
Figure 6: Period-doubling bifurcation: (a) attractor with period-1 when 𝑑 1 = 0 . 1 5 ; (b) attractor with period-2 when 𝑑 1 = 0 . 2 2 ; (c) attractor with period-2 when 𝑑 1 = 0 . 2 5 ; (d) chaotic attractor when 𝑑 1 = 0 . 3 3 .
fig7
Figure 7: Period-halving bifurcation: (a) attractor with period-2 when 𝑑 1 = 0 . 4 2 ; (b) attractor with period-1 when 𝑑 1 = 0 . 4 7 .

Now, in order to investigate the dynamical behaviors with respect to the death rate 𝑑 2 of the top predator via bifurcation diagrams, fix the death rate 𝑑 1 = 0 . 3 5 , 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 𝑑 2 are plotted in Figure 8.

fig8
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 𝑑 2 varies between 0 and 0.43. However, a stable limit cycle can be shown in the range 𝑑 2 > 0 . 4 3 . As is the case of 𝑑 1 , 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 𝑑 1 and 𝑑 2 . 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 ( 𝑑 1 from 0.24 to 0.45 and 𝑑 2 from 0.02 to 0.42) to system (1.2) are calculated in Figure 9 by using the method in [23].

fig9
Figure 9: The largest Lyapunov exponent of system (1.2): (a) 0 . 2 4 < 𝑑 1 < 0 . 4 5 ; (b) 0 . 0 2 < 𝑑 2 < 0 . 4 2 .

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).

References

  1. P. A. Braza, “The bifurcation structure of the Holling-Tanner model for predator-prey interactions using two-timing,” SIAM Journal on Applied Mathematics, vol. 63, no. 3, pp. 889–904, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. C. Cosner, D. L. Deangelis, J. S. Ault, and D. B. Olson, “Effects of spatial grouping on the functional response of predators,” Theoretical Population Biology, vol. 56, no. 1, pp. 65–75, 1999. View at Publisher · View at Google Scholar · View at Scopus
  3. S. Ruan and D. Xiao, “Global analysis in a predator-prey system with nonmonotonic functional response,” SIAM Journal on Applied Mathematics, vol. 61, no. 4, pp. 1445–1472, 2001. View at Publisher · View at Google Scholar
  4. M. Fan and Y. Kuang, “Dynamics of a nonautonomous predator-prey system with the Beddington-DeAngelis functional response,” Journal of Mathematical Analysis and Applications, vol. 295, no. 1, pp. 15–39, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  5. S. Gakkhar and R. K. Naji, “Seasonally perturbed prey-predator system with predator-dependent functional response,” Chaos, Solitons & Fractals, vol. 18, no. 5, pp. 1075–1083, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  6. R. Arditi and L. R. Ginzburg, “Coupling in predator-prey dynamics: ratio-dependence,” Journal of Theoretical Biology, vol. 139, no. 3, pp. 311–326, 1989. View at Google Scholar · View at Scopus
  7. H. I. Freedman and R. M. Mathsen, “Persistence in predator-prey systems with ratio-dependent predator influence,” Bulletin of Mathematical Biology, vol. 55, no. 4, pp. 817–827, 1993. View at Publisher · View at Google Scholar · View at Scopus
  8. L. A. Segel, Modeling Dynamic Phenomena in Molecular and Cellular Biology, Cambridge University Press, Cambridge, UK, 1984.
  9. S. Gakkhar and R. K. Naji, “Order and chaos in predator to prey ratio-dependent food chain,” Chaos, Solitons & Fractals, vol. 18, no. 2, pp. 229–239, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. A. Hastings and T. Powell, “Chaos in a three-species food chain,” Ecology, vol. 72, no. 3, pp. 896–903, 1991. View at Google Scholar · View at Scopus
  11. S.-B. Hsu, T.-W. Hwang, and Y. Kuang, “A ratio-dependent food chain model and its applications to biological control,” Mathematical Biosciences, vol. 181, no. 1, pp. 55–83, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. A. Klebanoff and A. Hastings, “Chaos in three-species food chains,” Journal of Mathematical Biology, vol. 32, no. 5, pp. 427–451, 1994. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  13. S. Lv and M. Zhao, “The dynamic complexity of a three species food chain model,” Chaos, Solitons & Fractals, vol. 37, no. 5, pp. 1469–1480, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  14. R. K. Naji and A. T. Balasim, “Dynamical behavior of a three species food chain model with Beddington-DeAngelis functional response,” Chaos, Solitons & Fractals, vol. 32, no. 5, pp. 1853–1866, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  15. C. Shen, “Permanence and global attractivity of the food-chain system with Holling IV type functional response,” Applied Mathematics and Computation, vol. 194, no. 1, pp. 179–185, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. J. A. Vano, J. C. Wildenberg, M. B. Anderson, J. K. Noel, and J. C. Sprott, “Chaos in low-dimensional Lotka-Volterra models of competition,” Nonlinearity, vol. 19, no. 10, pp. 2391–2404, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  17. M. Zhao and S. Lv, “Chaos in a three-species food chain model with a Beddington-DeAngelis functional response,” Chaos, Solitons & Fractals, vol. 40, no. 5, pp. 2305–2316, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  18. B. Liu, Z. Teng, and L. Chen, “Analysis of a predator-prey model with Holling II functional response concerning impulsive control strategy,” Journal of Computational and Applied Mathematics, vol. 193, no. 1, pp. 347–362, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  19. F. Cao and L. Chen, “Asymptotic behavior of nonautonomous diffusive Lotka-Volterra model,” Systems Science and Mathematical Sciences, vol. 11, no. 2, pp. 107–111, 1998. View at Google Scholar · View at Zentralblatt MATH
  20. H. I. Freedman and P. Waltman, “Persistence in models of three interacting predator-prey populations,” Mathematical Biosciences, vol. 68, no. 2, pp. 213–231, 1984. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. R. M. May, Stablitiy and Complexity in Model Ecosystem, Princeton University Press, Princeton, NJ, USA, 1973.
  22. M. T. Rosenstein, J. J. Collins, and C. J. De Luca, “A practical method for calculating largest Lyapunov exponents from small data sets,” Physica D, vol. 65, no. 1-2, pp. 117–134, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  23. J. C. Sprott, Chaos and Time-Series Analysis, Oxford University Press, New York, NY, USA, 2003.