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

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.


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: the magnitudes of the perturbations in r 1 t , m 1 t , and m 2 t , 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 π.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 M x , M y , and M z for any positive solution x t , y t , and z t of system 1.1 with all sufficiently large t.

Mathematical Analysis
Proof.Assume that x t , y t , z t is an arbitrary positive solution of system 1.1 , then the first equation of system 1.1 can yield

2.1
From Lemma 2.2, we get that, for any small positive constant ε, there exists a constant . Therefore V t is ultimately bounded, and it follows that each positive solution of system 1.1 is uniformly ultimately bounded.Hence, there exist two positive This completes the proof.
Proof.From Theorem 2.4, it is obvious that there exist positive constants M x , M y , M z , and On the other hand, we get from system 1.1 that

2.4
Then it follows from Lemma 2.
From the second equation of system 1.1 , we can obtain that there are two positive T 5 and m y such that y t ≥ m y for t ≥ T 5 .
From the third equation of system 1.1 , we can obtain that dz t /dt there are two positive constants T 6 and m z such that z t ≥ m z for t ≥ T 6 .Let T max{T 3 , T 4 , T 5 , T 6 } then we have m x ≤ x t ≤ M x , m y ≤ y t ≤ M y , and m z ≤ z t ≤ M z for t ≥ T .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.

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 Let us investigate the critical parameter k 2 .Figure 1 shows the bifurcation diagram and the corresponding largest Lyapunov exponents as a function of k 2 in the range 10 ≤ k 2 ≤ 25 for ϕ 0, π/2, π and ω 0.3.It is evident from Figure 1 that the system produces complex dynamics for 10 ≤ k 2 ≤ 25 when ϕ 0, π/2 and π, such as chaotic band, perioddoubling bifurcation, and chaotic band with wide or narrow periodic windows.Nonetheless, it is interesting to notice that the prey x biomass level increases with increase of k 2 , but it can be found that when ϕ π/2, the prey x 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 k 2 .The evidence for the cascade of periodic doubling leading to chaos can be seen in Figure 1 for 11 ≤ k 2 ≤ 20, 11 ≤ k 2 ≤ 19.5, and 11 ≤ k 2 ≤ 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 x biomass level, but also changes the complex dynamics when ϕ π/2.
To investigate more carefully the effect of synchronous variation in the intrinsic growth rate and death rate, three strange attractors are plotted for k 2 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 k 2 24, but when ϕ π, system 1.1 has a periodic attractor at k 2 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 ϕ.
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 k 2 23 and three different values of synchronous variation: ϕ 0, π/2, π.In this paper, the plane z 6 was chosen, and values of y were plotted against x.The resulting attractors and Poincaré section are shown in Figure 3.    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 .
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 x shown in Figure 1 b , the corresponding largest Lyapunov exponent 10 ≤ k 2 ≤ 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 longterm dynamics maxima for ϕ 0, π/2, π and the range 0 ≤ ω ≤ 0.4 and k 2 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.
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 z 5 .The attractors and Poincaré sections are shown in Figure 5.    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.
Based on the above analysis, it can be seen that the seasonal disturbance can enhance the prey x 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.

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

Figure 3 :
Figure 3: a Strange attractor with ϕ 0, k 2 23; b Poincaré section of attractor a ; c strange attractor with ϕ π/2 , k 2 23; d Poincaré section of attractor c ; e strange attractor ϕ π, k 2 23; f Poincaré section of attractor e , showing the values of the states x and y in the plane z 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

Figure 5 :
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 x and y in the plane z 5.
r 1 , m 1 , and m 2 are the average values of r 1 t , m 1 t , and m 2 t , respectively, and ω is the angular frequency of the fluctuations caused by seasonality.The parameters ε i i 1, 2, 3 0 ≤ ε i ≤ 1 represent the degree of seasonality, and r 1 ε 1 , m 1 ε 2 , and m 2 ε 3 are 1.1 where x t , y t , and z t are the densities of one prey and two predators at time t, respectively, a i i 1, 2 are the cropping rates, e i i 1, 2 denote the efficiency with which preys are converted by new predators, r 1 t k 2 is the carrying capacity of prey x, and k 1 is the value of the limiting prey.In other words, k 1 is the theoretical carrying capacity under ideal conditions if there is no wastage in preys, which is impossible in reality.k 2 /k 1 x t ≤ k 2 < k 1 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 Definition 2.1.System 1.1 is said to be uniformly persistent if, for any positive solution x t , y t , z t of system 1.1 , there exist positive constants m x , m y , m z , M x , M y , M z , and T > 0 such that m x ≤ x t ≤ M x , m y ≤ y t ≤ M y , and m z ≤ z t ≤ M z for t ≥ T .