Stability and Bifurcation Analysis of a Discrete Singular Bioeconomic System

The main concern of this paper is to discuss stability and bifurcation analysis for a class of discrete predator-prey interaction with Holling type II functional response and harvesting effort. Firstly, we establish a discrete singular bioeconomic system, which is based on the discretization of a system of differential algebraic equations. It is shown that the discretized system exhibits much richer dynamical behaviors than its corresponding continuous counterpart. Our investigation reveals that, in the discretized system, two types of bifurcations (i.e., period-doubling and Neimark–Sacker bifurcations) can be studied; however, the dynamics of the continuous model includes only Hopf bifurcation. Moreover, the state delayed feedback control method is implemented for controlling the chaotic behavior of the bioeconomic model. Numerical simulations are presented to illustrate the theoretical analysis. The maximal Lyapunov exponents (MLE) are computed numerically to ensure further dynamical behaviors and complexity of the model.


Introduction
Bioeconomics is linked closely to the early development of ideas in fisheries economics due to the pioneering work of Canadian economists Gordon [1] and Anthony Scott (in 1955). eir basic theories used recent developments in modeling of biological fisheries, initially the contributions made by Schaefer in 1954 and 1957 on introducing a systematic connection between fishing mechanism and growth of biological type through the implementation of mathematical modeling verified by experimental studies, and also associated itself to resource protection, ecology, and the environment [2]. ese concepts were developed from the multifishing science environment in Canada at that time. Modeling and fishing science developed rapidly during an innovative and productive period, especially among Canadian fishing researchers of various disciplines. Fishing mortality and population modeling were launched for economists, and novel interdisciplinary methods of modeling became obtainable for the economists, which made it feasible to measure the economical and biological impacts of various fisheries management decisions and fishing activities. Modern bioeconomics related to fisheries science can furnish perception into developing methods to deal with the overexploitation and complexities of overcapacity in marine fisheries where most are affected by lack of solid governance, changing coastal ecosystem dynamics, and natural fluctuations [3].
Moreover, Gordon in [1] suggested economic theory keeping in view the common property of resource, which was based on the effect of the harvest effort on an ecosystem by taking into account an economic perspective, assuming that x(t) and e(t) denote the density of harvested population and the harvest effort in an ecosystem, respectively; then the total cost is equal to ce(t), and the total revenue is equal to pe(t)x(t), where c denotes the cost of harvest effort, and p is used for the unit price of harvested population. en, the economic interest μ for the harvest effort by the harvested population is given by μ � e(t)(px(t) − c). (1) Taking into account Gordon [1] theory, Zhang et al. [4] studied a class of bioeconomic system with implementation of theory for singular systems. eir study was consisted of bifurcation analysis and chaos control for the proposed bioeconomic model. Later on, Liu et al. [5] reported stability, bifurcation analysis, and state feedback control for a class of predator-prey interaction with harvest effort on predator and stage structure for prey. Chakraborty et al. [6] studied a bioeconomic system with implementation of theory of differential algebraic equations. ey investigated stability, Hopf bifurcation, and state feedback control for a class of predator-prey interaction with time-delayed effect. Zhang et al. [7] investigated a singular bioeconomic model for preypredator interaction with diffusion and time delay. Zhang et al. [8] carried out comprehensive study related to theory, applications, complexity, and control of singular bioeconomic systems. Zhang et al. [9] explored the Hopf bifurcation for a predator-prey type bioeconomic system with two delays and predator harvesting. Meng and Zhang [10] discussed the qualitative behavior of a delayed singular bioeconomic predator-prey model without and with stochastic fluctuation. Liu et al. [11] analyzed the local dynamics and Hopf bifurcation for a biological economic model with Holling type II functional response and harvesting effort on prey. Liu et al. [12] proposed a singular fishery model for a class of prey-predator interaction with gestation delay for predator and maturation delay for prey. Liu et al. [13] formulated and discussed a singular predatorprey model by implementing commercial harvesting on predator with gestation delay for predator and maturation delay for prey. In [14], Li et al. studied a singular bioeconomic predator-prey system with Holling type II functional response and nonlinear harvesting on prey. Meng and Wu [15] discussed a singular prey-predator system with two delays, nonlinear predator harvesting and Beddington-DeAngelis functional response. Babaei and Shafiee [16] reported stability analysis, bifurcation, chaotic behavior, and control for a singular bioeconomic model of prey-predator interaction governed by an algebraic equation and 3-dimensional differential equations.
In case of mathematical modeling of predator-prey interaction, the research concerning interspecific interactions has been numerously based on continuous predator-prey systems of two variables. On the other hand, particular species, covering several classes of insects and seasonal plants, have nonoverlapping generations successively, and consequently, their population undergoes in discrete time steps. Populations with nonoverlapping generations can be modeled suitably with difference equations, otherwise known as iterative maps or discrete dynamical systems. Several authors have shown that nonoverlapping generations governed by iterative maps reveal complex and chaotic behavior, and the dynamics in such cases may yield a much richer set of patterns than those examined in continuoustime systems (cf. [17][18][19][20][21][22][23][24]). Furthermore, discrete-time models also have been used for rich dynamics of bioeconomic systems for some classes of predator-prey interaction. For example, in [25], the authors analyzed complex dynamics of a discrete-time bioeconomic system for predator-prey interaction with the implementation of Euler approximations. Wu and Chen [26] implemented the Poincaré scheme for discretization of a singular bioeconomic model and they analyzed period-doubling bifurcation, Neimark-Sacker bifurcation, and stability behavior. Liu et al. [27] studied the chaotic behavior of a discrete singular system related to the bioeconomic model of the prey-predator type.
Taking into account predator interaction with logistic growth and Holling type II functional response for prey population, we have the following system [28]: where x � x(t) and y � y(t) denote state variables for the densities of prey and predator at time t, respectively. Moreover, d is the intrinsic growth rate of prey, r represents the natural death rate of the predator, a is used for half capturing saturation constant, and b represents the maximal growth rate of the predator. Furthermore, d/k denotes the environmental carrying capacity for the prey population. Keeping in view (1) and (2), we obtain the following predator-prey biological economic model with Holling type II functional response with harvest effort: Applying the forward Euler scheme to system (3), we obtain the discrete-time predator-prey biological economic model with Holling type II functional response as follows: x n+1 � x n + hx n d − kx n − y n a + x n − e n , y n+1 � y n + hy n bx n a + x n − r , μ � e n px n − c , where h is the integral step size for the Euler approximation.
In this paper, we discuss some dynamical aspects of the discrete singular model (4). For this, the first existence of 2 Discrete Dynamics in Nature and Society interior (positive) fixed point and local dynamics of system (4) about biologically feasible equilibrium are carried out. Secondly, it is proved that system (4) undergoes perioddoubling bifurcation and Neimark-Sacker bifurcation by varying the economic profit μ as the bifurcation parameter. irdly, a state delayed feedback control strategy is applied to avoid bifurcating and chaotic behavior of bioeconomic model (4). At the end, numerical examples are presented for verification and illustration of our theoretical discussion.

Fixed Points and Stability Analysis
In order to study the qualitative behavior of the solutions of the nonlinear model (4), we study the existence of fixed points and their stability properties. From system (4), we can see that there exists a fixed point X 0 ≔ (x 0 , y 0 , e 0 ) in R + 3 if and only if X 0 is a solution of the following equations: rough a simple calculation, we obtain For biological considerations, we focus on the dynamics of the positive fixed point of the system (4). us, throughout the paper, we assume the conditions for the existence of a unique positive fixed point of system (4) as follows: In μk-plane, the existence of a unique positive fixed point of system (4) is depicted in Figure 1. e generalized Jacobian matrix J(x 0 , y 0 , e 0 ) of system (4) about interior (positive) fixed point (x 0 , y 0 , e 0 ) is computed as follows: en, it is easy to see that the generalized characteristic equation of the Jacobian matrix J(x 0 , y 0 , e 0 ) can be written as which on simplification yields where where H � Let F(λ) � λ 2 + Pλ + Q, then where In order to discuss the stability of the fixed point of (x 0 , y 0 , e 0 ), we need the following lemma.
Keeping in view Lemma 1, the following theorem is presented for local dynamics of system (4) about its positive fixed point. (4) satisfying the following conditions: is nonhyperbolic if one of the following conditions hold true: Next, for a � 3.9, b � 14, c � 3.7, r � 3.4, p � 5.8, d � 21.1, and h � 0.25, the dynamical classification of interior equilibrium of system (4) is depicted in μk-plane (see Figure 3).
Taking into account part (iv.1) of eorem 1, it is easily to observe that the eigenvalues about the equilibrium point (x 0 , y 0 , e 0 ) are given by λ 1 � − 1 and λ 2 � 1 − P with |λ 2 | ≠ − 1. On the other hand, if part (iv.2) of eorem 1 holds true, then one can obtain that the eigenvalues for the equilibrium point (x 0 , y 0 , e 0 ) are conjugate complex numbers with modulus one.
Next, we consider the following set: It is easy to see that the steady state (x 0 , y 0 , e 0 ) undergoes period-doubling (flip) bifurcation whenever parameters vary in a small neighborhood of F B .
Next, we consider the following curve: On the other hand, (x 0 , y 0 , e 0 ) undergoes the Neimark-Sacker (Hopf ) bifurcation whenever parameters vary in a small neighborhood of H B . In Section 3, we discuss the occurrence of period-doubling bifurcation around the

Bifurcation Analysis
Keeping in view the analysis of Section 3, we discuss the period-doubling bifurcation and Neimark-Sacker bifurcation about the positive fixed point (x 0 , y 0 , e 0 ) in this section. For this, choose parameter μ as bifurcation parameter for investigating the period-doubling bifurcation and Neimark-Sacker bifurcation about (x 0 , y 0 , e 0 ) by implementing the novel normal form theory of discrete singular systems, the center manifold theorem, and the bifurcation theory of discrete systems [29][30][31][32][33][34][35].

Period-Doubling
Bifurcation. First, we start our investigation related to period-doubling bifurcation for system (4) about its equilibrium (x 0 , y 0 , e 0 ) with a variation of parameters in a small neighborhood of F B . For this, we In this case, system (4) is described by the following 2-dimensional map: It is easy to see that system (17) has a unique positive fixed point (x 0 , y 0 , e 0 ) such that the eigenvalues are given by Next, we take μ * as a bifurcation parameter, considering a perturbation for system (17) as follows: where μ * ≪ 1 is taken as a small perturbation parameter. en, system (17) can be described in the following way: g(x, y, e) � 0, en, it is easy to see that Dg(x 0 , y 0 , e) � (e 0 p, 0, px 0 − c) such that the rank of Dg(x 0 , y 0 , e 0 ) � 1. On the other hand, a local parameterization ψ for the 2-dimensional smooth manifold defined by M g � (x, y, e) ∈ R 3 ; g(x, y, e) � 0, and rank Dg(x, y, e) � 1} about (x, y, e) ∈ B(x 0 , y 0 , e 0 ) ⊂ M g is given as follows: for any where X � (x, y, e), X 0 � (x 0 , y 0 , e 0 ), Z � (z 1 , z 2 ), and H: R 2 ⟶ R is a smooth mapping. For further details, the interested reader is referred to [36]. Taking into account the definition of ψ, we obtain the following: for arbitrarily chosen (x, y, e) ∈ B(x 0 , y o , e 0 ). en, system (19) is written as follows: Taking into account (24), we have 6 Discrete Dynamics in Nature and Society In order to compute the coefficients related the normal form, we need some further computation. From the aforementioned computation, one can easily get and from this, it follows that Assume that f iz i (ψ(Z)) represents the derivative of f i (ψ(Z)) with respect to z i and taking f iz i (X) � f iz i (ψ(Z)). In a similar way, one can adopt notations for f iz i z j (ψ(Z)) and f iz i z j z k (ψ(Z)). Next, it is easy to see that Putting X 0 into (28), one has the following: Discrete Dynamics in Nature and Society 7 en, from (28), it is easy to see that Consequently, one has the following: 8 Discrete Dynamics in Nature and Society Putting the value of X 0 in (31), one has the following: en, from (31), it follows that Discrete Dynamics in Nature and Society en, we get Putting the value of X 0 into (34), we obtain that Consequently, it follows from (29), (32), and (35) that Discrete Dynamics in Nature and Society Next, we consider the following nonsingular matrix: implementing the following translation: en, it is easy to see that system (24) can be written as follows: where In order to determine the center manifold W c (0, 0, 0) of (39) about the equilibrium point (0, 0) in a small neighborhood of μ * , we implement an application of the center manifold theorem [29], and it is easy to see that there exists a center manifold of the following form: for u and μ * sufficiently small. Furthermore, we assume that h u, μ * � c 0 u 2 + c 1 uμ * + c 2 μ * 2 + O |u| + μ * 3 . (42) en, it is easy to see that the center manifold satisfies the following:

12
Discrete Dynamics in Nature and Society Taking into account (42) and (43) and then comparing coefficients for (43), one can easily obtain that a 1 a 2 a 11 − a 12 1 + a 1 − a 2 Furthermore, taking into account (39), it is restricted to the center manifold W c (0, 0, 0) as follows: where In order for system (45) to undergo a flip bifurcation, we require that two discriminatory quantities α 1 and α 2 are not zero, where Keeping in view the aforementioned computation and bifurcation theory given in [36,37], we have the following theorem. that (a, b, c, d, k, p, r, μ 1 ) ∈ F B . If α 2 ≠ 0, then system (17) undergoes a flip bifurcation at the fixed point (x 0 , y 0 , e 0 ) when the parameter μ varies in a small neighborhood of μ 1 . Moreover, if α 2 > 0 (resp., α 2 < 0), then the period-2 orbits that bifurcate from (x 0 , y 0 , e 0 ) are stable (resp., unstable).

Discrete Dynamics in Nature and Society 13
In the last Section 5 related to numerical simulation, we choose some parametric values for system (4) such that it undergoes period-doubling bifurcation about positive equilibrium as μ varies in the suitable interval (see Figure 4).

Neimark-Sacker
Bifurcation. Finally, we discuss the Neimark-Sacker bifurcation about equilibrium (x 0 , y 0 , e 0 ) with variation of parameters (a, b, c, d, k, p, r, μ 2 , h) in a small neighborhood of H B . For this, we choose parameters (a, b, c, d, k, p, r, μ 2 , h) arbitrarily from H B , taking into account system (4) with (a, b, c, d, k, p, r, μ 2 , h) ∈ H B . In this case, system (4) is described by the following 2-dimensional map: en, it is easy to observe that map (48) has a unique positive fixed point (x 0 , y 0 , e 0 ). Next, we take μ * as a bifurcation parameter and consider a perturbation corresponding to map (48) given as follows: where Keeping in mind the aforementioned computation, one has the following result.

Chaos Control
Bifurcating behavior, chaos, and unstable fluctuations have always been considered as adverse criteria in biology; therefore, these are damaging for the reproduction of the biological population. Certainly, we require to take action to stabilize the biological population. For further details and applications of chaos control methods, we refer to [37][38][39][40][41][42][43][44][45][46][47][48][49][50] and the references therein. 16 Discrete Dynamics in Nature and Society In this section, we propose the following state delayed feedback control method for system (4): x n+1 � x n + hx n d − kx n − y n a + x n − e n + δ x n − x n− 1 , y n+1 � y n + hy n bx n a + x n − r , μ � e n px n − c , where δ is the feedback gain for the controlled system (63). Next, introducing u n ≔ x n − x n− 1 , we obtain the following controlled system equivalent to system (63): x n+1 � x n + hx n d − kx n − y n a + x n − e n + δu n , y n+1 � y n + hy n bx n a + x n − r , μ � e n px n − c .

Numerical Simulation and Discussion
In this section, our main purpose is to validate theoretical findings with numerical simulation. For this, first of all, the existence of the period-doubling bifurcation is illustrated through particular choice of biologically feasible parametric values. Choose a � 5. 8 (3) is given as follows: Moreover, the roots of P(λ) are − 1 and 0.942995. Consequently, (a, b, c, d, k, p, r, μ 1 , h) ∈ F B , and it follows the correctness of eorem 2. On the other hand, bifurcation diagrams for singular system (3) and corresponding largest Lyapunov exponents (LLE) are depicted in Figure 4. Secondly Moreover, the roots of P(λ) are 0.603567 − 0.797313i and 0.603567 + 0.797313i with |0.603567 ± 0.797313i| � 1. Consequently, (a, b, c, d, k, p, r, μ 2 , h) ∈ H B , and it follows the correctness of eorem 3. On the other hand, the bifurcation diagrams for singular system (3) and corresponding largest Lyapunov exponents (LLE) are depicted in Figure 5. Moreover, in the chaotic region, that is, for μ ∈ [0.678396, 1.1], some phase portraits of system (4) are depicted in Figure 6.
On the other hand, the characteristic polynomial of the Jacobian matrix of system (73) is given as follows: Taking into account the conditions of Lemma 2, we have that the positive fixed point of system (73) Figure 7.   Discrete Dynamics in Nature and Society

Conclusion
We discuss the dynamical behavior of a discrete-time singular bioeconomic model. Moreover, Euler's approximation is implemented to a bioeconomic model governed by the differential-algebraic system proposed in [11]. e Neimark-Sacker bifurcation and period-doubling bifurcation are studied for the discrete bioeconomic model with the implementation of normal form theory, bifurcation theory, and the center manifold theorem. On the other hand, we select μ (the economic profit parameter) as a bifurcation parameter. Our investigations show richer dynamical behaviors for the discrete-time bioeconomic model compared with its continuous counterpart studied in [11]. Numerical computation for maximum Lyapunov exponents ensures further dynamical behavior and complexity of the model. Such type of complex phenomena might result from economic profit [25]. With the variation of the bifurcation parameter μ, the biologically feasible fixed point resembles the stability of economic profit, and later on, the system may sacrifice its stability by undergoing period-doubling or Neimark-Sacker bifurcation, and consequently trajectories tend to a period-doubling cascade or an invariant circle. Our theoretical discussion reveals that if the economic revenue μ increases beyond a certain threshold value μ 1 (respectively, μ 2 ), a phenomenon of period-doubling bifurcation (respectively, Neimark-Sacker bifurcation) occurs. From eorem 2 (respectively, eorem 3), if the economic revenue μ is equal to or larger than the bifurcation value μ 1 (respectively, μ 2 ), the predator population, the prey population, and the harvesting effort will not stay at steady states, which will result in serious biological economic environmental imbalance. us, it is sensible for fishermen to keep the economic revenue within a certain range for the purpose of maintaining the sustainable development of biological resources. On the other hand, the suggested state feedback control method of the delayed type vanishes such fluctuating phenomena and derives the bifurcating behavior of singular discrete-time economic prey-predator system towards a stable situation. Computation reveals that this delayed-type control strategy can be applied only by modifying the effort of harvesting for the prey population density by taking into account the present prey population density and the previous one. Consequently, economists can formulate some strategies to restrain or encourage the harvesting attempts in practical applications, for example, adjusting market price and revenue, abating pollution, making an allowance to fishermen, so that the biological populations can remain at their stable states.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Discrete Dynamics in Nature and Society 21