Investigation on dynamics of an impulsive predator–prey system with generalized Holling type IV functional response and anti-predator behavior

In the present article, we propose and analyze a new mathematical model for a predator–prey system including the following terms: a Monod–Haldane functional response (a generalized Holling type IV), a term describing the anti-predator behavior of prey populations and one for an impulsive control strategy. In particular, we establish the existence condition under which the system has a locally asymptotically stable prey-eradication periodic solution. Violating such a condition, the system turns out to be permanent. Employing bifurcation theory, some conditions, under which the existence and stability of a positive periodic solution of the system occur but its prey-eradication periodic solution becomes unstable, are provided. Furthermore, numerical simulations for the proposed model are given to confirm the obtained theoretical results.


Introduction
Natural phenomena such as porous medium problems and periodically oscillatory waves that we usually observe in science, engineering or some other fields can be explained via mathematical models which are often formulated by means of a differential equation or a system of differential equations [1][2][3]. However, a relationship between predators and preys in an ecological system is considered as one of the most interesting natural phenomena because it is a type of the most fundamental biological systems in nature. Due to its universal existence and importance, the dynamical behaviors of predator-prey systems have been studied by a great number of scholars [4][5][6][7]. Interactions between predator and prey populations are extensively exploited for arranging economically damaging prey (pest) species in an eco-friendly community. In 1925-1926, Lotka and Volterra described the relationship between predators and their preys in an ecological system by establishing a system of two autonomous ordinary differential equations, namely the Lotka-Volterra model. The Lotka-Volterra system is the simplest model of predator-prey interactions based on linear per capita growth rates [8,9]. Some recent developments of complex predator-prey systems, which have been investigated in depth by many researchers to obtain their dynamics, are as follows. For example, Tiwari et al. [10] proposed and analyzed a mathematical model for a predator-prey system with multiple Allee effect acting on the growth rate of the prey population. Haldar et al. [11] formulated a delay-induced ecoepidemic model using a reconstructed Leslie-Gower-type growth rate and then analyzed the dynamics of the system with infection in the prey population. A predator-prey model with the additive Allee effect and the fear effect in the prey was presented and investigated by Lai et al. [12]. The bifurcation analysis of the model was explored by considering on the impact of these two effects. To further grasp other crucial predator-prey mathematical models in detail, one can refer to [13][14][15][16][17][18][19].
It has been found that the prey adopts several mechanisms for overcoming the predation pressure in their interactions. Two such techniques, which are used in the predator-prey models for an ecological system but not well scrutinized in the literature, are group defense and anti-predator behavior in prey [20]. Firstly, group defense [21][22][23][24] is a term used to describe the phenomenon whereby predation is decreased or even prevented altogether because of the enlarged ability of the prey to maintain its own survival via better defending or disguising itself when a sufficiently large number of prey appear. The explorations for predator-prey models with group defense were considered by many researchers such as in [24][25][26][27]. Particularly, Freedman and Wolkowicz introduced a general model of predatorprey interactions in which the prey behaves a group defense and there is no mutual interference among predators. Their model consisting of autonomous ordinary differential equations of generalized Gauss-type can be expressed as follows: x (t) = xg(x, k)yp(x), y (t) = y -d + q(x) , where x(t) and y(t) represent the density of prey and predator populations, respectively. The functions g, p and q are assumed to be continuously differentiable and the parameters d and k are positive constants. More details as regards system (1) can be found in [21]. Secondly, preys can sometimes jeopardize their predators. This leads to role reversals in predators and prey, namely anti-predator behaviors. Several authors studied the complicated dynamics of predator-prey systems with anti-predator behaviors such as in [20,[28][29][30][31][32][33]. In particular, an anti-predator behavior of prey [20] is a counterattacking approach in which adult preys attack the juvenile predators to decrease their future predation circumstances. Biological experiments demonstrate that anti-predator behaviors of prey populations are characterized as (i) morphological changes or behavior changes [34,35], or (ii) attacking their predators [29,36]. Normally, juvenile preys escaping from predation and becoming adult can counterattack and kill juvenile predators, but do not consume them so that a number of preys probably killed in the future are reduced [37]. Therefore, this cyclic dominance is very significant for predator-prey interactions because there are adverse effects on the biological control strategy. Here are some recent examples of studies of predator-prey models with anti-predator behaviors. In 2017, Sun et al. [38] introduced a predator-prey model equipped with anti-predator behaviors occurring only when the population size of the prey is greater than a threshold. They concluded that a large anti-predator rate induces the prey population to persist but whether we have coexistence between prey and predator populations depends upon the threshold activating anti-predator behavior. In 2019, Prasad et al. [20] analyzed the predator-prey model with an anti-predator behavior in prey to overcome the loss in predator populations via catering the predators with some additional food. However, we are quite interested in studying the effect of anti-predator behavior on the predator-prey system proposed by Tang and Xiao [30]. Their mathematical model is expressed as in which the functional response Ψ (x) is specifically chosen to be the simplified Monod-Haldane functional response Ψ (x) = βx a+x 2 , where β is the capture rate of the predator and a is the reciprocal of group defense in prey. The remaining parameters in system (2) are described as follows: r is the intrinsic growth rate of the prey, k is the prey carrying capacity, μ is the conversion fraction of prey into predator and d is the natural death rate of the predator population. The parameter η is the rate of anti-predator behavior of prey to the predator population in such a way that this behavior reduces the growth of the predator population because the prey population is not mainly fed on the predator population. All of the parameter values in system (2) are assumed to be positive. Actually, there are many types of functional response Ψ (x) indicating the change in the density of prey attacked per unit of time per predator when the prey density changes. Some interesting functional responses, widely used in predator-prey models, are Holling type, sigmoid type and Beddington-DeAngelis functional responses. More details as regards the functional responses can be found in [23,[39][40][41]. However, the generalized Holling type IV functional response or the Monod-Haldane function [5,[42][43][44][45] where β, a, b are positive constants, will be selected to use in our work due to its nonmonotonicity. Next, we adopt the predator-prey model with anti-predator behavior (2) via replacing Ψ (x) with the generalized Holling type IV functional response in Eq. (3). As mentioned, we obtain the following autonomous system: System (4) can be thought of as adding the anti-predator term, -ηxy, to the predator-prey model with prey group defense in Eq. (1) from which the functions g, p and q are identified as The functions g, p and q in (5) satisfy the conditions (2)-(4) given in [21], respectively, Using Descartes' rule of signs, if aη < βμ, then the equation Q (x) = 0 has exactly one positive real root denoted by x = M * . It is not difficult to see that Q(x), x ∈ [0, ∞) has a global maximum at M * . Now, we must assume for system (4) that Q(M * ) > d because of otherwise the predator cannot survive on the prey at any density. The predator-prey system (4) will be a part of our impulsive problem proposed later on.
The remaining parts of this paper are organized as follows. In Sect. 2, an impulsive differential system based on the predator-prey model (4) is formulated. The meaning of the state variables and the parameters of the system is given there. Section 3 provides the theoretical outcomes of the impulsive model including the basic results, the extinction and permanence and the bifurcation and existence of positive periodic solution obtained using the Floquet theory of an impulsive system and bifurcation theory. Numerical results showing that the system can exhibit several interesting phenomena are given in Sect. 4. We discuss and conclude our results in the last section.

Impulsive model formulation
Over the last few decades, impulsive differential equations have been found in almost every region of applied science and have been studied in many fields such as impulsive birth [46], impulsive vaccination [47], chemotherapeutic treatment of disease [48] and population ecology [49,50]. In general, they are exploited to describe the phenomena depending upon steep or instantaneous changes which make the differential system more intractable. A large number of researchers have suggested to apply an impulsive control strategy [50][51][52][53][54] to a modified predator-prey model for obtaining some desired results such as a reduction in pest populations in a farm. This approach is known as integrated pest management (IPM), which is a suitable method for manipulating pests (or preys) with a combination of biological, cultural, physical and chemical tools so that one can minimize economic, health and environmental risks [51,55]. Periodically releasing natural enemies or spraying pesticides at a different fixed time is regarded as a certain example of the IPM.
The present section is devoted to formulating an impulsive mathematical model for predator-prey model (4) in which the generalized Holling type IV functional response and the anti-predator behavior are inserted. Our impulsive differential system with a fixed moment can be expressed as and T is the impulsive period for releasing predators in order to remove target pests, prevent harmless pests from extinction and drive target pests to extinction, or control target pests at acceptably low level to avoid an increase of pest population due to an economic loss. The notations x(t + ) and y(t + ) are the right limits representing the values of x and y immediately after a pulse at time t. The number n ∈ N = {1, 2, 3, . . .} and p 1 , p 2 , p 3 are nonnegative constants with 0 ≤ p 1 , p 2 < 1, p 3 ≥ 0. Here, p 1 and p 2 represent the proportional constants of the prey and predator, respectively which are eradicated by using chemical poisons in agriculture or by harvesting. The constant p 3 is the density of predators released each time. Consequently, the last two equations of Eq. (7) can be considered as the IPM of the system. The values of the parameters r, k, β, μ, a, b, d, η are assumed to be positive. The descriptive list of the variables and parameters presented in system (7) is shown in Table 1. To the best of the authors' knowledge, the impulsive system (7) has never been proposed and investigated before so the dynamics of the system remains unclear, and falls within the scope of our exploration.

Analysis of the model
With a view on the use of the generalized Holling type IV functional response, the antipredator behavior and the impulsive IPM appearing in the impulsive predator-prey system (7), we are interested in analyzing some important properties of the system in this section.

Fundamental results
In the following section, we will show that a solution of system (7) is nonnegative and bounded above for sufficient large t. In addition, the system has a prey-eradication periodic solution.
denote the map defined by the right hand of system (7). Then a function V : R + × R 2 + → R + is said to belong to class V 0 if [56]: (i) V is continuous in (t, z) ∈ (nT, (n + 1)T] × R 2 + and for each z ∈ R 2 + , n ∈ N, exists.
(ii) V is locally Lipschitzian in z.
From the first equation of system (7), we have x (t) = 0, which is a contradiction. Furthermore, we also obtain where t ∈ (0, nT]. Hence, The proof for y(t) can be conducted in a like manner.
Next, we show that all solutions of system (7) are uniformly ultimately bounded. (7) for sufficient large t.

Lemma 3.2 There exists a constant M > 0 such that x(t) ≤ M and y(t) ≤ M for each solution z(t) = (x(t), y(t)) of system
and using the upper right derivative of V (t, z(t)) along a solution of system (7), we obtain for t = nT Hence, for t ∈ (nT, (n + 1)T] and t k = kT , k = 1, 2, . . . , n, Lemma 2.2 of [58] implies that Therefore, V (t, z) is ultimately bounded by M > 0. In consequence, x(t) ≤ M and y(t) ≤ M when t is sufficiently large. In other words, the positive solution z(t) of system (7) is uniformly ultimately bounded. This completes the proof.
The following comparison theorem on impulsive differential systems [57] will be used for system (7).
where g : Let R(t) be the maximal solution of the scalar impulsive differential equation is any solution of system (7). Let r(t) be the minimal solution of (16) existing on [0, ∞), and assume that the inequalities in (15) are reversed and h n is nonincreasing. Then V (t, z(t)) ≥ r(t) for all t ≥ 0.
Note that if we have some smoothness conditions on g to ensure the existence and uniqueness of solutions for (16), then R(t) and r(t) are exactly the unique solution of (16). Consider the following system: which is a subsystem of system (7). Some important properties of subsystem (17) are given as follows. It is not difficult to see that is a positive solution of system (17) for t ∈ (nT, (n + 1)T], n ∈ N with Since is the solution of system (17) for t ∈ (nT, (n + 1)T] where n ∈ N, we have the following lemma. (17) has a positive periodic solutionỹ(t) with a period T. For every solution y(t) of (17) with y 0 ≥ 0, we have |y(t) -ỹ(t)| → 0 as t → ∞. Therefore, system (7) has a prey-eradication periodic solution

Extinction and permanence
Before studying the permanence of system (7), we will give the condition ensuring the locally asymptotical stability of the prey-eradication periodic solution (0,ỹ(t)). Next, the condition for the permanence of the system will be given.
, y(t)) be any solution of system (7). Then (0,ỹ(t)) is locally asymptotically stable provided that Proof The local stability of the periodic solution (0,ỹ(t)) in Eq. (21) may be determined by scrutinizing the behavior of small amplitude perturbations of the solution. Defining then we may express where (t) is the fundamental solution matrix satisfying and (0) = I, the 2 × 2 identity matrix. Hence, The term ( * ) will not be used in further calculation; then there is no need to calculate the exact expression for ( * ) in the above equation.
The linearization of the last two equations of system (7) becomes According to the Floquet theory [58] of an impulsive differential equation, if both eigenvalues of the matrix have absolute values less than one, then the periodic solution (0,ỹ(t)) is locally asymptotically stable. Let λ 1 and λ 2 denote the eigenvalues of M . We have and It is obvious that |λ 1 | < 1. However, |λ 2 | < 1 if and only if rT - Hence, the solution (0,ỹ(t)) of system (7) is locally asymptotically stable when the condition (31) holds. This completes the proof. (7) is said to be permanent if there exist positive constants m, M and T 0 such that each positive solution z(t) = (x(t), y(t)) of system (7) with the initial values x(0 + ) > 0 and y(0 + ) > 0 satisfies m ≤ x(t) ≤ M and m ≤ y(t) ≤ M for all t ≥ T 0 . Note that m, M are independent of the initial values but T 0 may depend on the initial values. (7) is permanent if the condition rT -
Proof Suppose that z(t) = (x(t), y(t)) is a solution of system (7) with x(0 + ) > 0, y(0 + ) > 0. Then Lemma 3.2 guarantees that z(t) is bounded above when t is large enough. Consequently, there is a constant M > ra β > 0 such that x(t) ≤ M and y(t) ≤ M when t is sufficiently large.
Since μβxy a+bx 1 +x 2 1 ≥ 0, system (7) implies that where D 0 ≡ d + ηM > 0. Consider We see that, for t ∈ (nT, (n + 1)T], n ∈ N, is a positive solution of the comparison system (34) with Therefore, the solution of the comparison system (34) is for t ∈ (nT, (n + 1)T], n ∈ N and s(t) →s(t) as t → ∞. Then we have for some ε 2 > 0 and for all sufficiently large t. According to [57], we have y(t) ≥ s(t) by Lemma 3.3 and hence, when t is large enough, Then we only need to find a constant m 1 > 0 such that x(t) ≥ m 1 for all t large enough. Following the two steps below, we will obtain the desired result.
Step 1. From condition (32), we can equivalently get Based on condition (40), one can choose m 3 > 0 and ε 1 > 0 sufficiently small such that We will prove that x(t) < m 3 cannot hold for all t ≥ 0. Otherwise, from the second and fourth equations of system (7) and the fact that x 2 + bx ≥ 0 and ηx ≥ 0 due to x(t) ≥ 0, we consequently have From Lemmas 3.3 and 3.4, we obtain y(t) ≤ w(t) and w(t) →w(t) as t → ∞, where w(t) is the solution of the following system: , t ∈ nT, (n + 1)T .
Therefore, there exists a bound T 1 > 0 such that and for t ≥ T 1 . Let k ∈ N and kT ≥ T 1 . Integrating (46) on (nT, (n + 1)T], n ≥ k, we obtain where is expressed in Eq. (41). Then x((n + k)T) > x(nT) k → ∞ as k → ∞, which contradicts the boundedness of x(t). Hence, there exists a time t 1 > 0 such that x(t 1 ) ≥ m 3 .
Step 2. If x(t) ≥ m 3 for all t ≥ t 1 , then our aim is attained. Otherwise, let t * = inf t>t 1 {t : x(t) < m 3 }. Next, we consider two possible cases for t * .
From the second condition of (49) along with the inequalities in (54) and (57), we have which contradicts the assumption that x(t) ≤ m 3 for all t ∈ (t * , t * + T ]. Thus, there must exist t 2 ∈ (t * , t * + T ] so that x(t 2 ) > m 3 . Next, we lett = inf t>t * {t : x(t) > m 3 }. Since x(t) is left continuous and x(t + ) = (1p)x(t) ≤ x(t) when t = nT , we can conclude that x(t) ≤ m 3 for t ∈ (t * ,t) and x(t) = m 3 . For t ∈ (t * ,t), we assume that t ∈ (t * + (l -1)T, t * + lT] where l ∈ N and l ≤ n 2 + n 3 . From (55) we have Hence, we can conclude that x(t) ≥ m 1 for t ∈ (t * ,t). For t >t, the same arguments can be continued since x(t) ≥ m 3 . In consequence, when t is sufficiently large, we obtain x(t) ≥ m 1 > 0.

Remark 3.1 Let
Since g(0) = -ln ( 1 1-p 1 ) < 0, g(T) → ∞ as T → ∞ and g(T) = 0 has a unique positive root, denoted by T max . From Theorems 3.1 and 3.2 we know that T max is a threshold. If T < T max , then the prey-eradication periodic solution (0,ỹ(t)) is asymptotically stable. However, if T > T max , then system (7) is permanent.

Bifurcation and stability of positive periodic solution
In this section, we will investigate the existence of a nontrivial periodic solution to system (7) near the prey-eradication periodic solution (0,ỹ(t)) via bifurcation. Such a positive periodic solution is stable when the prey-eradication loses its stability. To do this, we apply Theorem 2 of [59] to system (7) by interchanging the state variables x(t), y(t) and giving a new variable for the period T. Letting x 1 (t) = y(t), x 2 (t) = x(t) and τ = T, system (7) consequently becomes All notations used in this section are the same as those in [59]. Then whereỹ is defined in Eq. (18) with T = τ and where τ 0 is the root ∂x 2 ) (τ 0 ,X 0 ) = 0 in which = ( 1 , 2 ) T is the flow associated to system (68). It will be notified in a later step that τ 0 = T max which is the positive root of g(T) = 0, where g(T) is defined in Eq. (66).
As defined above, we have Next, we obtain the following terms involving the flow = ( 1 , 2 ) T : , Consider It is not difficult to see from (73) and (74) that ∂ 2 2 (τ 0 ,X 0 ) is true.
In addition, we have Then we compute If d 0 = 0, then there is only one positive value of τ 0 satisfying the equation Furthermore, we obtain Since b 0 < 0 provided that the condition (75) holds.

Numerical simulations
The algorithm for solving impulsive differential equations is based on well-known numerical schemes [60][61][62] such as the spline approximation method, the θ -method, the multistep method and the Runge-Kutta method. However, the adaptive step-size control should be incorporated as a part of such methods because it can handle an impulsive differential equation, considered as a stiff problem involving rapidly changing components together with slowly changing ones, quite well. Here, numerical simulations for system (7) are obtained using the built-in function ode15s [63] in MATLAB which is a variableorder, variable step size and multistep solver based on the backward differentiation formulas [64]. The entire time of investigation is partitioned into equal subintervals. Each of the subintervals has length T (the impulsive period). The behavior of a solution for system (7) is characterized right before a pulse and immediately after the pulse by exploiting ode15s to numerically solve each subinterval of time. The information, evaluated at the jump point or the last time of the current subinterval, is used to compute the right limit of such a point. The right limit is then used as the initial condition of the next consecutive subinterval. The process is repeated in a similar manner for the rest subintervals and is terminated when we apply ode15s to perform numerical results for the last subinterval. The obtained numerical simulations in this section will be carried out to confirm our theoretical outcomes reported in Sect. 3. Studying the consequences of the impulsive effect to the dynamics of the impulsive prey-predator population model (7) will be done by the numerical experiments through adjusting the values of p 1 , p 2 , p 3 and T. Some parametric values for system (7) are exactly taken from the literature [23,30] but the others are appropriately estimated in order to satisfy the relevant conditions. The numerical simulations are as follows. Figure 1 displays the simulation results of the impulsive system of Eq. (7) with the parameter values r = 3.1, k = 2.067, β = 1.083, a = 1.031, b = 0.001, μ = 0.85, d = 0.3, η = 0.01. Figure 1(a) shows the coexisting circumstance in which the prey and predator densities are the periodic oscillations starting at (x(0), y(0)) = (0.5, 0.5). This is because the value of the parameters p 1 , p 2 , p 3 is all set to be zero, that is, any pest-management strategy is not yet applied to the system. The following numerical experiments confirm the result of Theorem 3.1 indicating that the solution of system (7) locally asymptotically converges to the prey-eradication periodic solution (0,ỹ(t)) as T < T max . When we simulate the system with (x(0), y(0)) = (1, 1) and p 1 = p 2 = 0, p 3 = 0.8 (in other words, there is only a periodic releasing of the predators but without giving pesticide), then the solution trajectory tends toward the oscillatory solution (0,ỹ(t)) when time passes long enough and T = 0.8 < T max = 0.9. This phenomenon can be seen in Fig. 1(b). It can be observed from Fig. 1(c) that the predator density y(t) oscillates in a stable periodic solution but the prey population x(t) quickly reduces to zero. This prey-eradication periodic solution occurs when we set p 1 = 0.85, p 2 = 0.2, p 3 = 0.9, x(0) = y(0) = 1 and T = 1 < T max = 1.335 in the simulation. According to the used parameter values for the last two simulations, one can verify that condition (22) in Theorem 3.1 is satisfied.
In Fig. 2, the numerical simulation beginning with (x(0), y(0)) = (0.1, 0.5) of the impulsive system of Eq. (7) locally asymptotically converges to the prey-eradication periodic solution (0,ỹ(t)) when time sufficiently increases and T = 2 < T max = 2.47. The parameter values used for this computation are r = 0.8, k = 0.5, β = 0.7, a = 0.3, b = 0.05, μ = 0.85, d = 0.85, η = 0.1, p 1 = p 2 = p 3 = 0.5. With this IPM, the density of prey x(t) converges asymptotically to zero but the density of predators y(t) tends to a periodic oscillation of period T = 2. This means that when t is large enough, the prey population tends to extinction, however, the number of the predators changes in terms of the periodic oscillation. These obtained results are in good agreement with Theorem 3.1 since the condition (22) in Theorem 3.1 holds.
Next, we utilize the same parameter values and initial condition as used in the simulation for Fig. 2 except using T = 8 > T max = 2.47 to simulate the numerical solution of the impulsive system (7). In consequence, the impulsive system with the use of this IPM is permanent as shown in Fig. 3. In other words, the prey and predator densities of the system are positively bounded as predicted in Theorem 3.2 because the condition (32) holds. In order to verify Theorem 3.3 applied for system (7), we must select the values of T, k, β, a, b, μ, η to satisfy all of the conditions in Theorem 3.3. Specifically, one can choose the impulsive period T such that T > T max and T is close to T max . In the final simulation, we here use the same parameter values and initial condition as used in the simulation for Fig. 2 except using T = 2.7 > T max = 2.47. The numerical solution plotted in Fig. 4 shows that the prey-eradication periodic solution (0,ỹ(t)) becomes unstable and then the preys and predators can coexist, which this characterization is exhibited in terms of a sustained positive periodic oscillation when time passes and T = 2.7 > T max = 2.47. This phenomenon allows the system to have a supercritical bifurcation, which can be used to control the magnitude of the prey population.

Discussion and conclusions
We have studied the dynamical behaviors of the mathematical impulsive model for the predator-prey system as proposed in Eq. (7). The impulsive system consists of the predator-prey model in which the generalized Holling type IV functional response and the anti-predator term are included and the impulsive control strategy concerning with giving pesticide and releasing of predators. Firstly, we have shown that system (7) has a positive solution which is uniformly ultimately bounded. Secondly, we have proved that the prey-eradication periodic solution (0,ỹ(t)) is locally asymptotically stable provided Figure 4 The numerical simulation demonstrates the occurrence of a stable positive periodic solution of the impulsive system in Eq. (7). Here, the parameter values and the initial condition are exactly the same as used in Fig. 2 except T = 2.7 > T max = 2.47: (a) the time series simulation of the prey density x(t); (b) the time series simulation of the predator density y(t); (c) the trajectory simulation plotted on the xy-plane that the impulsive period T is less than the critical value T max , or equivalently the condition (22) holds. Practically, if the impulsive control strategy in system (7) is selected to push the preys to extinction, then we can determine the impulsive period T based on the influences of spraying chemical pesticide and releasing natural predators such that T < T max . Of course, an absolute eradication of the prey population is not suitable for a real ecological system, so the effective prey-control practice should decrease the number of preys in the system to acceptable levels. In contrast, system (7) is permanent (i.e., the prey and predator populations are positively bounded) if T > T max . Profoundly, the system can lose the stability of the prey-eradication periodic solution and have a positive periodic solution when T > T max and T is close to T max . From the numerical experiments, we can observe that the narrower duration of T results in the lower number of preys. Thirdly, system (7) can undergo either supercritical or subcritical bifurcations via choosing the values of the parameters in such a way that the sign of the discriminant BC in Lemma 3.5 changes. Finally, changing the functional response, the anti-predator term or the IPM in system (7) could result in obtaining interesting and complicated dynamical behaviors of the preys and predators.