Stability and Bifurcation Analysis of Time Delayed Prey-Predator System with Holling Type-III Response Function

Interaction between prey and predator is a recurring event that occurs continuously and the presence of both can mutually affect each other’s population. This paper discusses the stability and bifurcation analysis of time delayed prey-predator system with Holling type-III response function incorporating a prey refuge. Holling type-III response function is used because the availability of the prey in nature is decreasing. Time delay represents the time for predators to find another prey population when the available population is decreasing. Dynamic analysis is used to study the influence of a given response function. The equilibrium point and stability analysis of the model with and without time delay has been calculated and analyzed. Model analysis under the influence of time delay and coefficient of competition among predators shows an indication of Hopf bifurcation in the neighborhood of the co-existing equilibrium point.


I. INTRODUCTION
P REDATOR is an organism that is looking for, hunting, and eating other organisms, while prey is an organism that is hunted and eaten by predators. The dynamic relationship between prey and their predator is a recurring event and will continue to be one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance. Alfred Lotka (1925) and Vito Volterra (1927) in Beals et al. [1] develop differential equations that describe the phenomenon of prey-predator for the first time. A pair of differential equations known as a model of Lotka-Volterra. There are several assumptions on the Lotka-Volterra models: the prey population grows exponentially in the absence of the predator, the predator population would starve without prey population, predators can consume an infinite number of preys, and no complexity environment. Furthermore, Beals et al. [1] state that one of the shortcomings of the Lotka-Volterra models is the reliance on unrealistic assumptions. Similar opinion was also expressed by Gasull et al. [2] that the Lotka-Volterra models is unrealistic because the prey population can grow without limit, mostly in the absence of the predators. In addition to that view, we begin to develop some models that are the modification of the Lotka-Volterra models, one of them is Holling-Tanner. Gasull et al. [2] revealed that the model of Holling-Tanner gives an overview of the competition going on Manuscript received July 29, 2019; accepted October 19, 2020. The authors are with the Department of Mathematics, Institut Teknologi Sepuluh Nopember, Kampus ITS Sukolilo-Surabaya 60111, Indonesia. E-mail: nurainamaziun@yahoo.com among the preys in a high density. In a high density, the prey will compete for their resources.
We have studied prey-predator phenomena. Among the references is a study conducted by Tapan Kumar Kar [3]. He analyzed the dynamic behavior of prey-predator models with Holling type-II function response incorporation of prey protection and assess the harvesting business as a control to prove that it is possible to came into effect on the cyclic behavior of the system and directs it to the required state. Jana et al. [4] studied the global and bifurcation analysis system of prey-predator combining time, of protection, using Holling response function of type-II. From the results, they obtained that the presence of protection has a significant effect on the coexistence of predator and prey population. Huang et al. [5] discuss the prey-predator models with Holling type-III function response combining prey protection. They have analyzed the model and discussed some of the significant qualitative results from a biological point of view.
Liu and Zhang [6] investigated the spatio-temporal dynamics of a delayed diffusive prey-predator model with nonlinear predator harvesting. Through mathematical analysis, they obtained the conditions for Turing and Hopf bifurcation. Ma et al. [7] presented a predator-prey system with Holling type function response incorporating prey refuge. By applying the analytical approaches, the dynamic behavior of the considered system is investigated, including stability, limit cycle and bifurcation. Balilo and Collera [8] considered delayed threespecies predator-prey model with non-monotonic functional response where two predator populations feed on a single prey population. Lajmiri et al. [9] studied the bifurcation and stability of a ratio-dependent predator-prey model with nonconstant predator harvesting rate. The analysis is carried out both analytically and numerically. Meng and Wu [10] proposed a differential algebraic predator-prey model including two delays, Beddington-DeAngelis functional response and nonlinear predator harvesting. Without considering time delay, the existence of singularity induced bifurcation is analyzed by regarding economic interest as a bifurcation parameter.
Based on these studies, we discuss the stability analysis and bifurcation of a predator-prey system incorporating a prey refuge in the time delay using Holling response function of type-III considering the prey migration from predators area to the refuge and vice versa.

A. Formulation of the Model
Suppose b 1 and b 2 are prey birth rate in refuge and predatory regions, respectively. While d 1 and d 2 are respectively the prey death rate in refuge and predatory regions and then x(t) and y(t) are respectively prey populations in the refuge and predatory region at time t. Parameter σ 1 is the migration rate of the prey population in the refuge region to the predatory region and σ 2 is the migration rate from predatory region to the refuge region. By taking r = b 1 − d 1 and s = b 2 − d 2 , we then obtained r and s as a respective intrinsic growth rate for the prey populations x and y. Prey growth in each region without the presence of predators is assumed logistics are also assumed that there are predators population in the predator area. Suppose P(t) is the predator population, where α is the maximum predator per capita consumption rate, i.e. the maximum number of prey populations from predatory area that can be eaten by a predator in each time unit α > 0, a is Michaelis-Menten constant and γ is the intra-specific competition coefficient or the density dependent mortality rate of predator population. Notation K and L respectively are the environmental carrying capacity for the prey populations x and y. The predator population consumes prey at the rate β y (we assumed 0 < β ≤ α). We assumed that predatory predators on prey only occurs in regions of prey, and the refuge region is predator-free, so that predation does not influence the growth of the predator population directly. Despite this, it is assumed that predators compete against each other to survive and response function used in the following equation is Holling type-II. It is assumed that all metabolism energy from only the food is used specifically for growth, and ultimately increase the predator population. Predator population consumes prey population at a constant rate, but the reproduction of predators after consuming prey is not directly done, thus will be put a delay required for the reproduction of prey. Suppose that the time interval between the individual and the number of dead prey individuals added to the predator population is considered as a time delay τ [4].
Form αyP a + y in the equation of Holling type II response function, according to [5] it usually describes the uptake of substrates by the microorganisms in microbial dynamics kinetics. If the predator is an invertebrate, it is always the case. Form αy 2 P a 2 + y 2 which is Holling type III suits the vertebral predator. In this paper, we will use Holling type III response function so the equation will be dx dt The initial conditions for the system are x(0) > 0, y(0) > 0, and P(0) > 0.

B. Domain, Terms, and Limitations of Early Settlement
To analyze system (2), we first discuss the boundedness criteria of the system (2). From the first two equations of system (2), we have K + L as the total carrying capacity of the total prey population. So we have x + y ≤ K + L + ε as t → ∞. Thus, we may take x ≤ K + ε 1 as t → ∞ and y ≤ L + ε 2 as t → ∞. where ε, ε 1 , ε 2 are three positive numbers.
and integrated both sides, we obtained the settlement area boundary, Therefore, the equations

C. Equilibria
We study all possible equilibria of system (2): 3) The interior (positive) equilibrium isB (x,ȳ,P) wherex is the positive root of the equation

D. Dynamic Behavior
In this subsection, we discuss the stability properties of the equilibria B 0 , B 1 and B 2 . The Jacobian for non delay and delay parameter of the system (2) is given by (3) and The characteristic equation is given by Clearly, for r − σ 1 > max rσ 2 s , σ 2 − s all the eigenvalues of system (2) become negative at B 0 (0, 0, 0) thus the system is locally asymptotically stable around B 0 (0, 0, 0).
2) Local Stability of Equilibrium Point B 1 (x 1 , y 1 , 0): Jacobian matrix for B 1 (x 1 , y 1 , 0) is given by The characteristic equation in the case of non delay (τ = 0) is given by Thus, using the Routh Hurwitz criterion if r 2 , r 4 and r 5 all positive, we can conclude that all the eigenvalues of the characteristic equation (5) contains negative real part. Hence system (2) is asymptotically stable around the equilibrium point. The characteristic equation in the presence of delay (τ = 0) is given by Thus, using the Routh Hurwitz criterion if r 2 , r 5 positive and β 2 1 − d 2 3 ≤ 0, we can conclude that all the eigenvalues of the characteristic equation (6) contains negative real part. Hence the system (2) is asymptotically stable around the equilibrium point.
It is interpreted that the presence of time delay does not change the stability from B 1 (x 1 , y 1 , 0) because, time delay required to seek other prey when the prey were previously reduced or discharged. In the equilibrium point B 1 (x 1 , y 1 , 0), there is no predator population (P = 0), so the equilibrium point B 1 (x 1 , y 1 , 0) will remain stable for τ ≥ 0 if and only if r 2 > 0, r 5 > 0 and β 2 1 − d 2 3 ≤ 0. 3) Local Stability of Equilibrium Point B 2 (x,ȳ,P): Jacobian matrix for B 2 (x,ȳ,P) is given by (7) and The characteristic equation is given by where (a 2 +ȳ 2 ) 2 ; m 6 = αȳ 2 a 2 +ȳ 2 ; 2βPȳ a 2 +ȳ 2 ; m 7 = σ 1 σ 2 ; for non delay (τ = 0), Equation (8) is given by where We can conclude that all the eigenvalues of the of the characteristic equation (9) contains negative real part. Hence the system (2) is asymptotically stable around the equilibrium point.
For delay (τ = 0), equation 8 is given by Suppose the equation 10 has imaginary root λ = iω. For τ > 0, if iω(ω > 0) is the root for equation (10) so: by separating the real and imaginary part on the equation (11), we then obtained Equation (12a) and (12b) each squared and then summed up and we retrieved (13) can be written by f (σ ) := σ 3 + Jσ 2 + Hσ + I = 0 If iω(ω > 0) is pure imaginary from equation (10), so the equation (14) must have positive real roots σ = ω 2 . Equation (14) is a third-degree polynomial equation. Therefore, if I < 0 then the possibility is equation (14) has one or three positive real root. Thus equation (10) has pure imaginary root and as a result there comes bifurcation or the stability changing condition from stable to unstable. Furthermore for I > 0 and because the polynomial (14) has three degrees, then the polynomial is guaranteed to have a real negative root. There is only one possibility of polynomial (14) has real positive roots, which are both real positive roots. In other words, all the roots are real numbers. Therefore there is a condition which led equation (14) to have three real roots which were analyzed using Sturm Sequences approach [11].
. The Forde theorem implicitly states that there is a stability changing of the equilibrium point B 2 (x * , y * , P * ) along with increasing τ(τ > 0). In other words, there is a value τ critical (τ c ) so that equation (11) has a pure imaginary root. For τ > τ c then there are roots to the right of the imaginary axis which causes the equilibrium point B 2 unstable. Furthermore, we will set the value of τ c that led to changes in the stability of the equilibrium point B 2 i.e. multiply equation (12a) with h 1 ω 2 − h 3 and equation (12b) with h 2 ω obtaining τ k = 1 with k = 0, ±1, ±2, . . .. The smallest value of τ is a value that causes a stability change of the equilibrium point B 2 . So, if it meets the qualification of the Theorem 1, then point equilibrium B 2 stable for 0 ≤ τ < τ 0 , when τ = τ 0 bifurcation occurs and for τ > τ 0 equilibrium point B 2 unstable. It is interpreted that the stability of the system at the equilibrium point B 2 is influence by the existence of time delay. This means that the displacement of prey in the present which is influenced by the amount of food (prey) before units of time can be change the nature or behavior of the system.

E. Hopf Bifurcation
Changes in stability from stable to unstable is called bifurcation. Bifurcation occurs at the equilibrium point B 2 (x,ȳ,P) when τ = τ 0 . Next, it will be investigated that the type of bifurcation which occurs in equilibrium point B 2 is Hopf bifurcation in order to prove equilibrium point B 2 (x,ȳ,P) experience by having transversal condition from λ that is If equation (8) is rewritten as then the total derivative of H(λ , τ) = 0 is taken real part so we obtained Equation (12a) and (12b) substituted into equation (17) to obtain It fulfills the transversal requirement and the Hopf bifurcation is occured during τ = τ 0 .
For the case without delay (τ = 0) the result is B 2 (7.39422, 1.94329, 1.39579). From Figure 1 dan 2, the number of predators seems to move up from time 0-10, otherwise the number of preys in the predator area is decreased, although it eventually move towards a stable equilibrium point B 2 (7.39422, 1.94329, 1.39579). The results obtained in  Figure 3 and 4, we can see that when the time is between 0-100, the number of preys and predators is directly proportional, i.e. when the number of prey populations in the predators area increasing, the number of predators population also increases. When the prey population in the predator region begins to decrease, then the number of predators population also decreases, which in turn It means although the predator takes time τ = 3.5 in seeking other prey when the prey population runs out, but the condition of the system is still stable.  Figure 5 and 6, we can see that for τ = 10 and γ = 0.2 system solution is getting farther from equilibrium point B 2 (7.39422, 1.94329, 1.39579). It can be observed from the prey population graphic in predator region oscillating higher and radius from 6 is greater. It means that by time the predator needs in seeking other prey populations when the prey population runs out by the existence of competition coefficient between predators of 0.2, making the condition of the system unstable. This means that the given time delay and the competition among predators greatly affects prey populations in the predator region.  Figure 7 and 8, for τ = 10, a decreasing competition among predators by reducing intra-specific coefficient among predators from 0.2 to 0.18, the solution of x(t), y(t) and P(t) around equilibrium point B 2 (7.36807, 1.808, 1.41864). It is interpreted that with a time delay τ = 10 > τ 0 by reducing competition coefficient among predators make the system stable. This means that the competition among predators is very influential in the ecosystem (2).

IV. CONCLUSIONS
Based on the analysis of the results, three equilibrium points were obtained, i,e B 0 = (0, 0, 0), B 1 = (x 1 , y 1 , 0), and B 2 = (x,ȳ,P). The stability of the system at the equilibrium point B 0 = (0, 0, 0) and B 1 = (x 1 , y 1 , 0) is not affected by the delay so that these two equilibrium points are stable for each τ ≥ 0. It is because of the delay time that is only used by predators to look for other prey population when the existing prey population decreases or runs out. While at the two equilibrium points is free of predators. The time delay affects the stability of the system at the equilibrium point B 2 = (x,ȳ,P), which is a stable equilibrium point for τ ∈ [0, τ 0 ) and occurs at the time of bifurcation τ = τ 0 .
In this study, only a local stability analysis was carried out at the equilibrium point, for this reason in future studies to better analyze the global stability.