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

Authors: Hengguo Yu; Min Zhao
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

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


---

## Body

## 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.1)dxdt=r1(t)x(k2-xk1-x)-a1βxyb1+x+c1y-a2xzb2+x+c2z,dydt=(e1a1βxb1+x+c1y-m1(t))y,dzdt=(e2a2xb2+x+c2z-m2(t))z,r1(t)=r1(1+ε1sin⁡(ωt)),m1(t)=m1(1+ε2sin⁡(ωt+φ)),m2(t)=m2(1+ε3sin⁡(ωt+φ)),
where x(t),y(t), and z(t)are the densities of one prey and two predators at time t, respectively, ai(i=1,2) are the cropping rates, ei(i=1,2) denote the efficiency with which preys are converted by new predators, r1(t)k2 is the carrying capacity of prey x, and k1 is the value of the limiting prey. In other words, k1 is the theoretical carrying capacity under ideal conditions if there is no wastage in preys, which is impossible in reality. (k2/k1)(x(t)≤k2<k1)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. bi(i=1,2) are saturation constants, and ci(i=1,2) are scaling factors expressing the impact of predator interference, and β is the relative superiority of predator y. Without loss of generality, it can be assumed that β>1; accordingly, a predator with density dz/dt is an inferior predator, while a predator with density dy/dt is a superior predator because the ratio (dy/dt)/(dz/dt)∝βt increases with time. r1, m1, and m2 are the average values of r1(t), m1(t), and m2(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 r1ε1, m1ε2,and m2ε3 are the magnitudes of the perturbations in r1(t), m1(t), and m2(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 π.


## 2. Mathematical Analysis

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 mx, my, mz, Mx, My, Mz, and T>0 such that mx≤x(t)≤Mx, my≤y(t)≤My, and mz≤z(t)≤Mz for t≥T.Lemma 2.2 (see [26]).
Ifa>0,b>0, (dh(t)/dt)≤(≥)h(t)(b-ah(t)), and h(0)>0, then, for any small constant ε>0, there exists a positive constant T such that h(t)≤(b/a)+ε(h(t)≥(b/a)-ε), for t≥T.Lemma 2.3 (see [27]).
Ifa>0, b>0, (dh(t)/dt)≤(≥)b-ah(t), and h(0)>0, then h(t)≤(≥)(b/a)+(h(0)-(b/a))e-at for t≥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 constantsMx, My, and Mz for any positive solution (x(t), y(t), and z(t)) of system (1.1) with all sufficiently larget.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)dx(t)dt≤r1(t)x(t)k2-x(t)k1-x(t)≤r1(t)x(t)k2-x(t)k1-k2≤r1k2(1+ε1)k1-k2x(t)-r1(1+ε1)k1-k2x2(t).
From Lemma 2.2, we get that, for any small positive constant ε, there exists a constant T1>0 such that x(t)≤k2+ε:=Mx for t≥T1.
Defining thatV(t)=x(t)+(y(t)/e1)+(z(t)/e2), then
(2.2)V̇(t)=ẋ(t)+ẏ(t)e1+ż(t)e2=r1(t)x(t)k2-x(t)k1-x(t)-m1(t)e1y(t)-m2(t)e2z(t)≤r1k2(1+ε1)k1-k2x(t)-r1(1-ε1)k1-k2x2(t)-m1(1-ε2)e1y(t)-m2(1-ε3)e2z(t).
Let0<L1<min{m1(1-ε2),m2(1-ε3)}, then
(2.3)V̇(t)+L1V(t)≤(r1k2(1+ε1)k1-k2+L1)x(t)-r1(1-ε1)k1-k2x2(t)+L1-m1(1-ε2)e1y(t)+L1-m2(1-ε3)e2z(t).
Thus V̇(t)+L1V(t)is bound. Set L2=[r1(1+ε1)+(k1-k2)L1]2/(4r1(1-ε1)(k1-k2)) such that V̇(t)≤-L1V(t)+L2. According to Lemma 2.3, V(t)≤(L2/L1)+(V(0)-(L2/L1))e-L1t. 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 My, Mz(k2+(My/e1)+(Mz/e2)≤[r1(1+ε1)+(k1-k2)L1]2/(4r1L1(1-ε1)(k1-k2))) and T2 such that y(t)≤Myand z(t)≤Mz, for t≥T2. This completes the proof.Theorem 2.5.
IfW1>0, W2>0and W3>0, then system (1.1) is uniformly persistent.Proof.
From Theorem2.4, it is obvious that there exist positive constantsMx, My, Mz, and T3=max{T1,T2}>0 such that x(t)≤Mx, y(t)≤My, z(t)≤Mz for t≥T3.
On the other hand, we get from system (1.1) that
(2.4)dx(t)dt≥r1(t)x(t)k2-x(t)k1-a1βMyb1x(t)-a2Mzb2x(t)≥(r1k2(1-ε1)k1-a1βMyb1-a2Mzb2)x(t)-r1(1+ε1)k1x2(t).
Then it follows from Lemma 2.2 that if W1=(r1k2(1-ε1)/k1)-(a1βMy/b1)-(a2Mz/b2)>0, there is a positive constant T4such that x(t)≥((r1k2b1b2(1-ε1)-a1βb2k1My-a2k1b1Mz)/r1b1b2(1+ε1))-ε:=mxfor t≥T4.
From the second equation of system (1.1), we can obtain that  (dy(t)/dt)≥((e1a1βmx/(b1+k2+c1My))-m1(1+ε2))y(t). It is obvious to know that if W2=((e1a1βmx/(b1+k2+c1My))-m1(1+ε2))>0, there are two positive T5and my such that y(t)≥my for t≥T5.
From the third equation of system (1.1), we can obtain that (dz(t)/dt)≥((e2a2mx/(b2+k2+c2Mz))-m2(1+ε3))z(t). It is obvious to know that if W3=((e2a2mx/(b2+k2+c2My))-m2(1+ε3))>0, there are two positive constants T6 and mz such that z(t)≥mz for t≥T6. Let T=max{T3,T4,T5,T6} then we have mx≤x(t)≤Mx, my≤y(t)≤My, and mz≤z(t)≤Mz 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.


## 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 k2and ω are obtained for three different cases: φ=0, π/2,π. All other parameters are fixed as follows:a1=0.8,a2=0.6, b1=0.25, b2=2.5, e1=0.45, e2=0.35, c1=0.25, c2=0.55, m1=0.1, m2=0.15, r1=0.8, k1=30, β=0.35, ε1=0.5, ε2=0.1, and ε3=0.5.Let us investigate the critical parameterk2. Figure 1 shows the bifurcation diagram and the corresponding largest Lyapunov exponents as a function of k2 in the range 10≤k2≤25 for φ=0, π/2, πand ω=0.3. It is evident from Figure 1 that the system produces complex dynamics for 10≤k2≤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 preyxbiomass level increases with increase of k2, 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 k2. The evidence for the cascade of periodic doubling leading to chaos can be seen in Figure 1 for 11≤k2≤20, 11≤k2≤19.5, and 11≤k2≤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.Bifurcation diagram of system (1.1) with initial conditions x0=(0.1,0.1,0.1). (a) φ=0; (b) φ=π/2; (c) φ=π; (d) the largest Lyapunov exponent (k2 from 10 to 25) for system (1.1).
(a)(b)(c)(d)To investigate more carefully the effect of synchronous variation in the intrinsic growth rate and death rate, three strange attractors are plotted fork2=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 k2=24, but whenφ=π, system (1.1) has a periodic attractor at k2=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) Chaos whenk2=24, φ=0; (b) chaos when k2=24, φ=(π/2); (c) periodic attractor when k2=24, φ=π.
(a)(b)(c)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 k2=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. 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) Strange attractor withφ=0, k2=23; (b) Poincaré section of attractor (a); (c) strange attractor with φ=(π/2), k2=23; (d) Poincaré section of attractor (c); (e) strange attractor φ=π, k2=23; (f) Poincaré section of attractor (e), showing the values of the states x and y  in the plane z=5.
(a)(b)(c)(d)(e)(f)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≤k2≤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 k2=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.Bifurcation diagram of system (1.1) with initial conditions x0=(0.1,0.1,0.1). (a) φ=0; (b) φ=(π/2); (c) φ=π.
(a)(b)(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 (z=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) 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.
(a)(b)(c)(d)(e)(f)Based on the above analysis, it can be seen that the seasonal disturbance can enhance the preyx 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.



---

*Source: 150359-2012-02-15.xml*