Seasonally Perturbed Prey-Predator Ecological System with the Beddington-DeAngelis Functional Response

Abstract

On the basis of the theories and methods of ecology and ordinary differential equation, a seasonally perturbed prey-predator system with the Beddington-DeAngelis functional response is studied analytically and numerically. Mathematical theoretical works have been pursuing the investigation of uniformly persistent, which depicts the threshold expression of some critical parameters. Numerical analysis indicates that the seasonality has a strong effect on the dynamical complexity and species biomass using bifurcation diagrams and Poincaré sections. The results show that the seasonality in three different parameters can give rise to rich and complex dynamical behaviors. In addition, the largest Lyapunov exponents are computed. This computation further confirms the existence of chaotic behavior and the accuracy of numerical simulation. All these results are expected to be of use in the study of the dynamic complexity of ecosystems.

1. Introduction

Population communities are embedded in periodically varying environments. Therefore, it is appropriate to identify the functional role that seasons play in the behavior of population communities [1]. This seasonal variation can cause changes in the dynamics of an ecological system. The study of ecological systems which are subjected to seasonal perturbations is important for both theoretical and experimental ecologists. This study should also take into account the intrinsic nature of environmental and seasonal perturbations [2]. The basic problem is to understand the relationship between the magnitude of the seasonal variations and the complexity of the system. Numbers of studies have been performed on the interactions between the seasons and the internal biological rhythms of simple prey-predator ecosystems [3–6]. These studies have shown that these interactions can have spectacular consequences, such as multiplicity of attractors, catastrophes, and chaos [7].

In recent decades, it has been demonstrated that complex dynamics can appear in continuous-time models with three or more species [8–12], and specifically that nonlinear dynamics, including cycles, quasicycles, and chaos, can occur in such biological systems. In this context, continuous-time and discrete-time systems [13–21] have been discussed extensively by a number of researchers. The Beddington-DeAngelis functional response is introduced by Beddington [22] and DeAngelis et al. [23], independently. The main difference of this functional response from other functional responses is that it contains an extra term presenting mutual interference by predators [24]. Although a direct link between the predators and preys cannot be established unless quantitative methods are used, the precious works clearly show that the amount of three species are often related, and a change in one species can cause a change in the others, especially predator. Thus, we apply Beddington-DeAngelis functional response to describe their relationship with sufficient accuracy in this paper.

This paper considers a seasonally perturbed prey-predator system with the Beddington-DeAngelis functional response, which can be described by the following differential equations:

𝑑𝑥𝑑𝑡=𝑟₁(𝑡)𝑥𝑘₂−𝑥𝑘₁−𝑎−𝑥¹𝛽𝑥𝑦𝑏₁+𝑥+𝑐₁𝑦−𝑎₂𝑥𝑧𝑏₂+𝑥+𝑐₂𝑧,𝑑𝑦=𝑒𝑑𝑡¹𝑎₁𝛽𝑥𝑏₁+𝑥+𝑐₁𝑦−𝑚₁(𝑡)𝑦,𝑑𝑧=𝑒𝑑𝑡²𝑎₂𝑥𝑏₂+𝑥+𝑐₂𝑧−𝑚₂(𝑡)𝑟𝑧,₁(𝑡)=𝑟₁1+𝜀₁(,𝑚sin𝜔𝑡)₁(𝑡)=𝑚₁1+𝜀₂,𝑚sin𝜔𝑡+𝜑₂(𝑡)=𝑚₂1+𝜀₃sin𝜔𝑡+𝜑, (1.1) where 𝑥(𝑡),𝑦(𝑡), and 𝑧(𝑡)are the densities of one prey and two predators at time 𝑡, respectively, 𝑎_𝑖(𝑖=1,2) are the cropping rates, 𝑒_𝑖(𝑖=1,2) denote the efficiency with which preys are converted by new predators, 𝑟₁(𝑡)𝑘₂ is the carrying capacity of prey 𝑥, and 𝑘₁ is the value of the limiting prey. In other words, 𝑘₁ is the theoretical carrying capacity under ideal conditions if there is no wastage in preys, which is impossible in reality. (𝑘₂/𝑘₁)(𝑥(𝑡)≤𝑘₂<𝑘₁)expresses the efficiency of nutrient utilization by a species and has a value between zero and one. If the ratio approaches unity, the efficiency is high; a lower ratio indicates that population increase is quickly restricted by the limiting prey. 𝑏_𝑖(𝑖=1,2) are saturation constants, and 𝑐_𝑖(𝑖=1,2) are scaling factors expressing the impact of predator interference, and 𝛽 is the relative superiority of predator 𝑦. Without loss of generality, it can be assumed that 𝛽>1; accordingly, a predator with density 𝑑𝑧/𝑑𝑡 is an inferior predator, while a predator with density 𝑑𝑦/𝑑𝑡 is a superior predator because the ratio (𝑑𝑦/𝑑𝑡)/(𝑑𝑧/𝑑𝑡)∝𝛽^𝑡 increases with time. 𝑟₁, 𝑚₁, and 𝑚₂ are the average values of 𝑟₁(𝑡), 𝑚₁(𝑡), and 𝑚₂(𝑡), respectively, and 𝜔 is the angular frequency of the fluctuations caused by seasonality. The parameters 𝜀_𝑖(𝑖=1,2,3)(0≤𝜀_𝑖≤1) represent the degree of seasonality, and 𝑟₁𝜀₁, 𝑚₁𝜀₂,and 𝑚₂𝜀₃ are the magnitudes of the perturbations in 𝑟₁(𝑡), 𝑚₁(𝑡), and 𝑚₂(𝑡), respectively [25]. Finally, the parameter 𝜑, where 0≤𝜑≤2𝜋, can be interpreted as the difference in phase angle between the seasonality of the above three parameters. Clearly, when𝜑=0, there is synchronous variation in the intrinsic growth rate and death rate, while when 𝜑=𝜋, the variation is antisynchronous. In this paper, only three values of phase angles will be considered: 𝜑=0,𝜋/2, and 𝜋.

2. Mathematical Analysis

Definition 2.1. System (1.1) is said to be uniformly persistent if, for any positive solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (1.1), there exist positive constants 𝑚_𝑥, 𝑚_𝑦, 𝑚_𝑧, 𝑀_𝑥, 𝑀_𝑦, 𝑀_𝑧, and 𝑇>0 such that 𝑚_𝑥≤𝑥(𝑡)≤𝑀_𝑥, 𝑚_𝑦≤𝑦(𝑡)≤𝑀_𝑦, and 𝑚_𝑧≤𝑧(𝑡)≤𝑀_𝑧 for 𝑡≥𝑇.

Lemma 2.2 (see [26]). If 𝑎>0,𝑏>0, (𝑑ℎ(𝑡)/𝑑𝑡)≤(≥)ℎ(𝑡)(𝑏−𝑎ℎ(𝑡)), and ℎ(0)>0, then, for any small constant 𝜀>0, there exists a positive constant 𝑇 such that ℎ(𝑡)≤(𝑏/𝑎)+𝜀(ℎ(𝑡)≥(𝑏/𝑎)−𝜀), for 𝑡≥𝑇.

Lemma 2.3 (see [27]). If 𝑎>0, 𝑏>0, (𝑑ℎ(𝑡)/𝑑𝑡)≤(≥)𝑏−𝑎ℎ(𝑡), and ℎ(0)>0, then ℎ(𝑡)≤(≥)(𝑏/𝑎)+(ℎ(0)−(𝑏/𝑎))𝑒^−𝑎𝑡 for 𝑡≥0.

Now we can obtain the threshold expression of some critical parameters under the condition of all species persistence.

Theorem 2.4. There exist positive constants 𝑀_𝑥, 𝑀_𝑦, and 𝑀_𝑧 for any positive solution (𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡)) of system (1.1) with all sufficiently large𝑡.

Proof. Assume that (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) is an arbitrary positive solution of system (1.1), then the first equation of system (1.1) can yield 𝑑𝑥(𝑡)𝑑𝑡≤𝑟₁(𝑡)𝑥(𝑡)𝑘₂−𝑥(𝑡)𝑘₁−𝑥(𝑡)≤𝑟₁(𝑡)𝑥(𝑡)𝑘₂−𝑥(𝑡)𝑘₁−𝑘₂≤𝑟₁𝑘₂1+𝜀₁𝑘₁−𝑘₂𝑥(𝑡)−𝑟₁1+𝜀₁𝑘₁−𝑘₂𝑥²(𝑡). (2.1) From Lemma 2.2, we get that, for any small positive constant 𝜀, there exists a constant 𝑇₁>0 such that 𝑥(𝑡)≤𝑘₂+𝜀∶=𝑀_𝑥 for 𝑡≥𝑇₁.
Defining that 𝑉(𝑡)=𝑥(𝑡)+(𝑦(𝑡)/𝑒₁)+(𝑧(𝑡)/𝑒₂), then 𝑉̇(𝑡)=𝑥̇(𝑡)+𝑦̇(𝑡)𝑒₁+𝑧̇(𝑡)𝑒₂=𝑟₁(𝑡)𝑥(𝑡)𝑘₂−𝑥(𝑡)𝑘₁(𝑡)−𝑚−𝑥¹(𝑡)𝑒₁𝑦(𝑡)−𝑚₂(𝑡)𝑒₂𝑧(𝑡)≤𝑟₁𝑘₂1+𝜀₁𝑘₁−𝑘₂𝑥(𝑡)−𝑟₁1−𝜀₁𝑘₁−𝑘₂𝑥²(𝑡)−𝑚₁1−𝜀₂𝑒₁𝑦(𝑡)−𝑚₂1−𝜀₃𝑒₂𝑧(𝑡). (2.2)
Let 0<𝐿₁<min{𝑚₁(1−𝜀₂),𝑚₂(1−𝜀₃)}, then 𝑉̇(𝑡)+𝐿₁𝑉(𝑡)≤𝑟₁𝑘₂1+𝜀₁𝑘₁−𝑘₂+𝐿₁𝑥(𝑡)−𝑟₁1−𝜀₁𝑘₁−𝑘₂𝑥²(𝑡)+𝐿₁−𝑚₁1−𝜀₂𝑒₁𝑦(𝑡)+𝐿₁−𝑚₂1−𝜀₃𝑒₂𝑧(𝑡). (2.3) Thus 𝑉̇(𝑡)+𝐿₁𝑉(𝑡)is bound. Set 𝐿₂=[𝑟₁(1+𝜀₁)+(𝑘₁−𝑘₂)𝐿₁]²/(4𝑟₁(1−𝜀₁)(𝑘₁−𝑘₂)) such that 𝑉̇(𝑡)≤−𝐿₁𝑉(𝑡)+𝐿₂. According to Lemma 2.3, 𝑉(𝑡)≤(𝐿₂/𝐿₁)+(𝑉(0)−(𝐿₂/𝐿₁))𝑒^−𝐿1𝑡. Therefore 𝑉(𝑡) is ultimately bounded, and it follows that each positive solution of system (1.1) is uniformly ultimately bounded. Hence, there exist two positive 𝑀_𝑦, 𝑀_𝑧(𝑘₂+(𝑀_𝑦/𝑒₁)+(𝑀_𝑧/𝑒₂)≤[𝑟₁(1+𝜀₁)+(𝑘₁−𝑘₂)𝐿₁]²/(4𝑟₁𝐿₁(1−𝜀₁)(𝑘₁−𝑘₂))) and 𝑇₂ such that 𝑦(𝑡)≤𝑀_𝑦and 𝑧(𝑡)≤𝑀_𝑧, for 𝑡≥𝑇₂. This completes the proof.

Theorem 2.5. If 𝑊₁>0, 𝑊₂>0and 𝑊₃>0, then system (1.1) is uniformly persistent.

Proof. From Theorem 2.4, it is obvious that there exist positive constants𝑀_𝑥, 𝑀_𝑦, 𝑀_𝑧, and 𝑇₃=max{𝑇₁,𝑇₂}>0 such that 𝑥(𝑡)≤𝑀_𝑥, 𝑦(𝑡)≤𝑀_𝑦, 𝑧(𝑡)≤𝑀_𝑧 for 𝑡≥𝑇₃.
On the other hand, we get from system (1.1) that 𝑑𝑥(𝑡)𝑑𝑡≥𝑟₁(𝑡)𝑥(𝑡)𝑘₂−𝑥(𝑡)𝑘₁−𝑎₁𝛽𝑀_𝑦𝑏₁𝑥(𝑡)−𝑎₂𝑀_𝑧𝑏₂𝑥(𝑡)≥𝑟₁𝑘₂1−𝜀₁𝑘₁−𝑎₁𝛽𝑀_𝑦𝑏₁−𝑎₂𝑀_𝑧𝑏₂𝑥(𝑡)−𝑟₁1+𝜀₁𝑘₁𝑥²(𝑡). (2.4) Then it follows from Lemma 2.2 that if 𝑊₁=(𝑟₁𝑘₂(1−𝜀₁)/𝑘₁)−(𝑎₁𝛽𝑀_𝑦/𝑏₁)−(𝑎₂𝑀_𝑧/𝑏₂)>0, there is a positive constant 𝑇₄such that 𝑥(𝑡)≥((𝑟₁𝑘₂𝑏₁𝑏₂(1−𝜀₁)−𝑎₁𝛽𝑏₂𝑘₁𝑀_𝑦−𝑎₂𝑘₁𝑏₁𝑀_𝑧)/𝑟₁𝑏₁𝑏₂(1+𝜀₁))−𝜀∶=𝑚_𝑥for 𝑡≥𝑇₄.
From the second equation of system (1.1), we can obtain that (𝑑𝑦(𝑡)/𝑑𝑡)≥((𝑒₁𝑎₁𝛽𝑚_𝑥/(𝑏₁+𝑘₂+𝑐₁𝑀_𝑦))−𝑚₁(1+𝜀₂))𝑦(𝑡). It is obvious to know that if 𝑊₂=((𝑒₁𝑎₁𝛽𝑚_𝑥/(𝑏₁+𝑘₂+𝑐₁𝑀_𝑦))−𝑚₁(1+𝜀₂))>0, there are two positive 𝑇₅and 𝑚_𝑦 such that 𝑦(𝑡)≥𝑚_𝑦 for 𝑡≥𝑇₅.
From the third equation of system (1.1), we can obtain that (𝑑𝑧(𝑡)/𝑑𝑡)≥((𝑒₂𝑎₂𝑚_𝑥/(𝑏₂+𝑘₂+𝑐₂𝑀_𝑧))−𝑚₂(1+𝜀₃))𝑧(𝑡). It is obvious to know that if 𝑊₃=((𝑒₂𝑎₂𝑚_𝑥/(𝑏₂+𝑘₂+𝑐₂𝑀_𝑦))−𝑚₂(1+𝜀₃))>0, there are two positive constants 𝑇₆ and 𝑚_𝑧 such that 𝑧(𝑡)≥𝑚_𝑧 for 𝑡≥𝑇₆. Let 𝑇=max{𝑇₃,𝑇₄,𝑇₅,𝑇₆} then we have 𝑚_𝑥≤𝑥(𝑡)≤𝑀_𝑥, 𝑚_𝑦≤𝑦(𝑡)≤𝑀_𝑦, and 𝑚_𝑧≤𝑧(𝑡)≤𝑀_𝑧 for 𝑡≥𝑇. This completes the proof.

By using the theoretical analysis, we obtain the threshold expression of some critical parameters under the condition of all species persistence, which in turn provides a theoretical basis for the numerical simulation.

3. Simulation Analysis and Results

To study complex dynamics of the prey-predator system with seasonality, the solution of the system (1.1) with initial conditions is obtained numerically for the biologically feasible range of parametric values. In order to investigate the interplay among seasonal perturbation and some critical parameters, the diagrams as a function of 𝑘₂and 𝜔 are obtained for three different cases: 𝜑=0, 𝜋/2,𝜋. All other parameters are fixed as follows:𝑎₁=0.8,𝑎₂=0.6, 𝑏₁=0.25, 𝑏₂=2.5, 𝑒₁=0.45, 𝑒₂=0.35, 𝑐₁=0.25, 𝑐₂=0.55, 𝑚₁=0.1, 𝑚₂=0.15, 𝑟₁=0.8, 𝑘₁=30, 𝛽=0.35, 𝜀₁=0.5, 𝜀₂=0.1, and 𝜀₃=0.5.

Let us investigate the critical parameter𝑘₂. Figure 1 shows the bifurcation diagram and the corresponding largest Lyapunov exponents as a function of 𝑘₂ in the range 10≤𝑘₂≤25 for 𝜑=0, 𝜋/2, 𝜋and 𝜔=0.3. It is evident from Figure 1 that the system produces complex dynamics for 10≤𝑘₂≤25 when 𝜑=0, 𝜋/2 and 𝜋, such as chaotic band, period-doubling bifurcation, and chaotic band with wide or narrow periodic windows. Nonetheless, it is interesting to notice that the prey𝑥biomass level increases with increase of 𝑘₂, but it can be found that when 𝜑=𝜋/2, the prey 𝑥 biomass level is relatively high. Comparison of Figures 1(a), 1(b), and 1(c) suggests that seasonality may still give a region of periodic solutions for 𝑘₂. The evidence for the cascade of periodic doubling leading to chaos can be seen in Figure 1 for 11≤𝑘₂≤20, 11≤𝑘₂≤19.5, and 11≤𝑘₂≤21.5. However, Figure 1(a) shows chaotic bands with a narrow periodic window, and Figures 1(b) and 1(c) show chaotic bands with wide or narrow periodic windows. It is obvious that the seasonal disturbance not only enhances the prey 𝑥 biomass level, but also changes the complex dynamics when 𝜑=𝜋/2.

(a)

(b)

(c)

(d)

Figure 1 

Bifurcation diagram of system (1.1) with initial conditions 𝑥0=(0.1,0.1,0.1). (a) 𝜑=0; (b) 𝜑=𝜋/2; (c) 𝜑=𝜋; (d) the largest Lyapunov exponent (𝑘2 from 10 to 25) for system (1.1).

To investigate more carefully the effect of synchronous variation in the intrinsic growth rate and death rate, three strange attractors are plotted for 𝑘₂=24 and three different values of the synchronous variation parameter: 𝜑=0, 𝜋/2, 𝜋. Figure 2 shows that when 𝜑 is 0 or (𝜋/2), system (1.1) has a chaotic attractor at 𝑘₂=24, but when𝜑=𝜋, system (1.1) has a periodic attractor at 𝑘₂=24. It can be observed from Figure 2 that seasonality in the intrinsic growth rate and death rate can change the chaotic behavior of the system at various values of the synchronous variation𝜑.

(a)

(b)

(c)

Figure 2 

(a) Chaos when 𝑘2=24, 𝜑=0; (b) chaos when 𝑘2=24, 𝜑=(𝜋/2); (c) periodic attractor when 𝑘2=24, 𝜑=𝜋.

For further analysis of the effects of seasonality, three chaotic attractors as well as a Poincaré section (the Poincaré section is a classical technique for analyzing dynamic systems [27, 28]) have been obtained for 𝑘₂=23 and three different values of synchronous variation: 𝜑=0, 𝜋/2, 𝜋. In this paper, the plane (𝑧=6) was chosen, and values of 𝑦 were plotted against 𝑥. The resulting attractors and Poincaré section are shown in Figure 3. Figure 3(a) shows a chaotic attractor at 𝜑=0, and Figure 3(b) shows a Poincaré section for the chaotic attractor shown in (a). Three wide strips can be seen in Figure 3(b). Figure 3(c) shows a chaotic attractor at 𝜑=𝜋/2, and Figure 3(d) shows a Poincaré section for the chaotic attractor shown in (c). One wide strip can be seen in Figure 3(d). Figure 3(e) shows a chaotic attractor at 𝜑=𝜋, and Figure 3(f) shows a Poincaré section for the chaotic attractor shown in (e). Four wide strips can be seen in Figure 3(f). By comparison of the three Poincaré sections, it is apparent that the three chaotic attractors have different structures in the three different cases 𝜑=0, 𝜋/2, 𝜋. However, they might be better described as stranger chaotic attractors and fractals. These results show that seasonal perturbation can change chaotic attractor structure. In a word, seasonal perturbation is sufficient to change the chaotic attractor trajectory in system (1.1).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

(a) Strange attractor with 𝜑=0, 𝑘2=23; (b) Poincaré section of attractor (a); (c) strange attractor with 𝜑=(𝜋/2), 𝑘2=23; (d) Poincaré section of attractor (c); (e) strange attractor 𝜑=𝜋, 𝑘2=23; (f) Poincaré section of attractor (e), showing the values of the states 𝑥 and 𝑦  in the plane 𝑧=5.

Moreover, convincing evidence for deterministic chaos has come from several recent experiments [29–34]. The results of chaos studies have confirmed the importance of detecting and exploring chaos. Here, the largest Lyapunov exponents are considered because they have proved to offer the most useful diagnostics for a chaotic system [35–40]. The Lyapunov exponents take into account the average exponential rates of divergence or convergence of nearby orbits in phase space [29]. For a chaotic attractor, the values of the largest Lyapunov exponent must be positive; if they are always negative, chaos cannot be observed, and a stable state or a period attractor will instead appear. In the bifurcation diagram for resource 𝑥 shown in Figure 1(b), the corresponding largest Lyapunov exponent 10≤𝑘₂≤25 can be calculated for the system described by (1.1), which is shown in Figure 1(d). It should be stressed that the result is consistent with Figure 1(b), which shows the accuracy and effectiveness of numerical simulation. Moreover, using the simulation of the largest Lyapunov exponents, the existence of chaotic behavior in system (1.1) can be further confirmed.

The other critical parameter to be investigated is 𝜔. Figure 4 depicts the local long-term dynamics maxima for 𝜑=0, 𝜋/2, 𝜋 and the range 0≤𝜔≤0.4 and 𝑘₂=20. As 𝜔 increases from 0 to 0.14, Figure 4(a) shows that the system has rich dynamics, including periodic windows, period doubling, chaos, and a chaotic crisis. However, only chaos, chaotic crisis, and periodic windows are observed in Figures 4(b) and 4(c). Figure 4(b) shows that when 𝜔 is slightly increased beyond 0.14, the size of the chaotic attractor abruptly changes, constituting a type of attractor crisis. As 𝜔 continues to increase, a cascade of period-doubling bifurcation occurs, the period-2 attractor changes to a period-4 attractor, then a period-8 attractor and period-doubling bifurcation finally leads system (1.1) into chaos. However, only chaos and periodic windows are observed in Figures 4(a) and 4(c). In the range (0.18,0.4), the solution is chaotic with intermittent periodic windows. Similar behavior is observed in Figures 4(a) and 4(c). Hence, it is obvious that the seasonality has a strong effect on the dynamical behaviors of system (1.1) with the synchronous variation 𝜑 varying.

(a)

(b)

(c)

Figure 4 

Bifurcation diagram of system (1.1) with initial conditions 𝑥0=(0.1,0.1,0.1). (a) 𝜑=0; (b) 𝜑=(𝜋/2); (c) 𝜑=𝜋.

In order to study the effect of seasonality, three chaotic attractors as well as a Poincaré section have been obtained for 𝜔=0.26 and three different values of the synchronous variation parameter: 𝜑=0, 𝜋/2, 𝜋 and the plane (𝑧=5). The attractors and Poincaré sections are shown in Figure 5. Figure 5(a) shows a chaotic attractor at 𝜑=0, and Figure 5(b) shows a Poincaré section for the chaotic attractor shown in (a). Two wide point regions can be seen in Figure 5(b). Figure 5(c) shows a chaotic attractor at 𝜑=𝜋/2, and Figure 5(d) shows a Poincaré section for the chaotic attractor shown in (c). Two point regions which are wider than in the Poincaré section shown in (b) can be seen in Figure 5(d). Figure 5(e) shows a chaotic attractor at𝜑=𝜋, and Figure 5(f) shows a Poincaré section for the chaotic attractor shown in (e). Two point regions which are narrower than in the other two Poincaré sections can be seen in Figure 5(f). By comparison of the three Poincaré sections, it is apparent that the three chaotic attractors have different structures in the three different cases 𝜑=0, 𝜋/2, 𝜋. However, they are all fractals. These results show that the angular frequency of the fluctuations caused by seasonality can change chaotic attractor structure with the synchronous variation 𝜑 varying.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

(a) Strange attractor with 𝜑=0, 𝜔=0.26; (b) Poincaré section of attractor (a); (c) strange attractor with 𝜑=(𝜋/2), 𝜔=0.26; (d) Poincaré section of attractor (c); (e) strange attractor with 𝜑=𝜋, 𝜔=0.26; (f) Poincaré section of attractor (e), showing the values of the states 𝑥 and 𝑦 in the plane 𝑧=5.

Based on the above analysis, it can be seen that the seasonal disturbance can enhance the prey 𝑥 biomass level for𝜑=𝜋/2, in which result is agreed with some results in reality. Further, it is also interesting to point out that the seasonality in three different parameters can come into rich and complex dynamical behaviors, but these dynamical behaviors are different. Moreover, the use of mathematical model with seasonality is considered to investigate some biological problems, and the numerical simulation provides an approximation of the real biological system behaviors; hence, these results can promote the study of ecological dynamics.

4. Conclusions and Remarks

In this paper, a prey-predator model with Beddington-DeAngelis functional response and seasonal perturbation has been studied analytically and numerically. Mathematical theoretical works have been pursuing the investigation of uniformly persistent, which depicts the threshold expression of some critical parameters and in turn provides a theoretical basis for the numerical simulation. Numerical analysis indicates that the seasonality has a strong effect on the dynamical complexity and species biomass using bifurcation diagrams and Poincaré sections. Bifurcation diagrams show that the seasonality in three different parameters can give rise to rich and complex dynamical behaviors, including periodic, period-doubling, period-halving, chaotic crises, and chaos. By comparing the Poincaré sections of various chaotic attractors, it can be determined that these chaotic attractors have different structures with the seasonality in three different parameters. Using numerical simulation of the largest Lyapunov exponents, the existence of chaotic behavior in system (1.1) and the accuracy and effectiveness of numerical simulation can be further confirmed. All these results are expected to be of use in the study of the dynamic complexity of ecosystems.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 30970305 and Grant no. 31170338) and also by the Zhejiang Provincial Natural Science Foundation of China (Grant no. Y505365).