The Existence of a Limit-Cycle of a Discrete-Time Lotka-Volterra Model with Fear Effect and Linear Harvesting

. Modeling the interaction between prey and predator plays an important role in maintaining the balance of the ecological system. In this paper, a discrete-time mathematical model is constructed via a forward Euler scheme, and then studied the dynamics of the model analytically and numerically. The analytical results show that the model has two ﬁxed points, namely the origin and the interior points. The possible dynamical behaviors are shown analytically and demonstrated numerically using some phase portraits. We show numerically that the model has limit-cycles on its interior. This guarantees that there exists a condition where both prey and predator maintain their existence periodically.


Introduction
The journey of mathematical modeling has an important role in the wide area of ecological problems. Several researchers are interested in studying the interaction between two or more populations due to its contribution to describing the existence of populations in their ecosystem [1][2][3][4]. The famous ones are given by Lotka (1925) and Volterra (1926), which is defined by [5,6] dx dt = ax − bxy, where x and y denote the population density of prey and predator. The parameter a, b, c, and d define the intrinsic growth rate, the predation rate on prey, the conversion of predation rate to the birth rate of the predator, and the rate of natural death of predator, respectively. The system (1) becomes the base model for researchers to develop a new model based on the reality that has been found in nature. For example, the Gause-type (GT) predator-prey model which includes the logistic growth rate on prey [7], the Rosenzweig-MacArthur (RM) model, which assumes that the prey is logistically growth and the predation process has a saturation point [8], and the Bazykin-type (BT) model which assumes that there exists intraspecific competition on predator [9]. In this paper, we assume that the existence of the predator affects the birth rate of prey. When the population density of predator increases, the birth rate of prey decreases due to the fear of prey on predator. This phenomenon is called the fear effect and successfully applied to GT, RM, and BT models. See [10][11][12] for detailed mathematical models. By applying where k is the level of fear of prey due to the appearance of the predator. One of the biological components which almost always included in the ecosystem is the appearance of harvesting. This condition involves humans, which need other organisms for their food. Harvesting generally occurs in forestry and fishery [13][14][15]. By assuming both prey and predator are proportionally harvested, we acquire where q 1 and q 2 are the prey and predator harvesting, respectively.
If we consider the operator to develop a model with a deterministic approach, we have continuous-time and discrete-time models. The continuous-time model is given by first-order and fractional-order derivatives as the operator. In this paper, we use the discrete-time model as the operator. This decision is taken by considering the mathematical aspect rather than its biological meanings. The discrete-time model has some advantages rather than its continuous ones, such as complexity in dynamical behaviors and suitability for approximating the real data since data has discrete formulae [16][17][18]. By using the forward Euler method as in [19], we obtain the following discrete-  time model.
According to our literature review, although eq. (4) is a simple predator-prey model with fear effect and harvesting, this model is still relevant and has never been studied by other researchers. Therefore, we trust that our results have a contribution to ecological sciences and application in nature. We organized this paper as follows: All analytical results consist of all feasible fixed points, and their dynamical behaviors are given in Section 2. We explore more dynamical behaviors by giving numerical results on Section 3. We end this paper in Section 4 by presenting the conclusion of our works.

Analytical Results
In this section, we present the dynamical behaviors of the model analytically. We start by investigating the existence of fixed points. The fixed point of model (4) is obtained by solving the following equations.
The first fixed point is given by the origin E 0 = (0, 0), which states the condition when both prey and predator will be extinct. The second fixed point is the interior point given by E 1 = (x * , y * ) where x * = d+q 2 c and y * is the positive solution with respect to y from the following quadratic equation.
This fixed point describes the condition where both prey and predator maintain their existence in the ecosystem. Using Decartes' rules of sign, we confirm that this fixed point exists and is unique only when a > q 1 or the intrinsic growth rate is greater than its harvesting rate. Therefore, we have To study the dynamical behaviors around these fixed points, we apply linearization on the model (4). We obtain the Jacobian matrix as follows.
(iv) non-hyperbolic if h = 2 q 1 −a or 2 d+q 2 . Proof. By substituting E 1 = (x * , y * ) to the Jacobian matrix (6), we get  which has two eigenvalues as follows.
We confirm that [20], all of the given statements on Theorem 1 are analytically proven.
Theorem 2. The interior point E 1 = (x * , y * ) is always a source.
Proof. The Jacobian matrix (6) at E 1 = (x * , y * ) is which gives a polynomial characteristic as follows.

Numerical Results
In this section, we demonstrate some numerical simulations to support the analytical findings in the previous section. Since the model talks about a general predator-prey model without specifically studying an ecological case, we set the parameters experimentally using hypothetical values. We start the simulation by setting the parameter values as follows. a = 4, b = 3, c = 2, d = 2, q 1 = 6, q 2 = 2, k = 0.5, and h = 0.05. (10) According to analytical results, E 1 does not exist as the impact of q 1 = 6 > 4 = a. Moreover, according to Theorem 1, we have a sink E 0 . We obtain the phase portrait and time-series in Figure 1. Four initial values are given around the origin. All of the solutions given by those initial values converge to E 0 = (0, 0). See the phase portrait given by Figure 1(a). From Figure 1(b), although all solutions converge to the origin E 0 , we also confirm that the population density of the predator goes to its peak point before finally tending to E 0 . Now, we replace the value of q 1 with 0.2 and portray the phase portrait and time-series on Figure 2. we find that E 0 losses its stability and a fixed point E 1 ≈ (2, 0.4305) occurs in interior. From Figure 2(a), the solution of model (4), which starts closed enough to E 1 stays away and seems to go to infinity for n → ∞. See Figure 2(b) for the time-series. This condition support Theorem 2.
Finally, by using similar parameter values as in eq. (10), replacing q 1 = 2, and h = 10 −4 , we portray the numerical simulation given by Figure 3. Although E 1 is still unstable, an interesting phenomenon occurs. For each initial value, the solution is still in a similar orbit and isolating the interior point E 1 as in Figure 3(a). This means that we have some limit-cycles on the interior. We also give Figure 2(b) to show the solution is oscillating in constant frequency and amplitude. This end with the numerical simulations.

Conclusion
The dynamics of a discrete-time Lotka-Volterra model with fear effect and linear harvesting have been studied. Two types of fixed points have been founded. The first  fixed point is called the origin, which is asymptotically stable (sink) when the harvesting rate on prey is greater than the intrinsic growth rate of prey and the step-size less enough to some biological conditions depending on the harvesting rate, intrinsic growth rate, and death rate. When the harvesting rate on prey is less than the intrinsic growth rate, the origin becomes unstable, and the second fixed point occurs, namely the interior point. This interior point is always unstable. Furthermore, the detailed dynamics of each fixed point have been studied numerically. We have confirmed that there exist two types of instability of interior point depending on the value of stepsize. For some large values of step-size, the solution oscillates and increases unbounded. When the step-size is decreased and becomes small enough, the instability of the interior point becomes bounded where for each initial condition, the solution goes periodically and isolates the interior point, which is called the limit-cycle. The biological meaning of this condition state that both prey and predator will never be extinct. The density of both populations is maintained periodically.