Coexistence and harvesting optimal policy in three species food chain model with general Holling type functional response

In this paper, we have discussed harvesting of prey and intermediate predator species. Both are subjected to Holling type I–V functional response. Conditions for local and global stability of the nonnegative equilibria are verified. The permanent coexistence criterion of the model system and existence of optimal equilibrium solution of the control problem are demonstrated. Maximum sustainable yield and maximal net present revenue are determined. To confirm analytical results, numerical solution has been carried out using the Matlab™ ODE solver ODE45 and the simulations show the model system reveals complex behavior (such as oscillations), which reflects the real situation.


| INTRODUCTION
The population size and coexistence of species in ecological system are dynamic and complex. One of the important objectives in mathematical ecology is the study of dynamical relationship among the population sizes of predators and the respective prey. Species interactions in an ecosystem is expressed mathematically in a nonlinear differential equation form. In this study we proposed a general nonlinear mathematical model that describes three interacting species with harvesting prey and intermediate predator population. Real life example is plantherbivore-Omnivore. The three interacting species are two predators (a generalist we call from now onwards intermediate predator Y , a specialist we call from onwards top predator Z), and one source prey species X where both are subjected with a general Holling type I-IV functional responses.
A three species predator-prey model with prey switching where, top predator and intermediate predator are competed with each other for same prey with new Holling type functional response. It demonstrates an intermediate predator is the key species in regulating the feeding mechanism between the prey and the top predator (Mortoja et al., 2020).
The dynamical behaviors of three species such as toxin-producing Phytoplankton, Zooplankton, and Fish in a fishery system is studied with Holling type II functional response (Panja & Mondal, 2015). Further they extended predator-prey interaction among Phytoplankton, Zooplankton, and Fish in the presence of Environmental toxin and includes harvesting on the Fish (Panja et al., 2017).
Transition to chaotic behavior is established (Upadhyay & Raw, 2011) via period-doubling bifurcation and some sequences of distinctive period-halving bifurcation leading to limit cycles are observed, furthermore, Holling type IV predator response function is considered. The rich dynamics in a tritrophic food chain mathematical model, consists of three species: prey, intermediate predator, and top predator while both the prey and intermediate predator are being subjected to Holling type I-IV functional response is investigated in Dawed et al. (2020). Some research have been proposed to understand the dynamics in a food chain (one predator-two prey). Coexistence and harvesting control policy in a food chain model with mutual defense of prey is explained well in Tchepmo Djomegni et al. (2019) and in contrast to many approaches, this study considers mutualism (for defense against predators) between the two groups of prey. As it is explained in Tchepmo Djomegni et al. (2019), in a group defense when prey density is sufficiently large, predators can go extinct. To maintain the sustainability, one can harvest the abundant prey population (Chen et al., 2013). Due attention has been given to study the impact of harvesting the dynamics in predator-prey systems . Harvesting has social, economic, and ecological interests. Socially, it can be used for domestic feeding. Economically, it is used for food, clothes, and medicine purpose. From ecological perspective, harvesting can be used for resource renewal and sustainability. Furthermore, in the paper by Bergland et al. (2019), a predatory prey model of livestock prey density with sigmoidal harvesting rate is proposed and studied.
Functional response is an indispensable component of predation models which determines both prey death rate and the predator increase rate. C. S. Holling has suggested three classes of functions to cover the most interesting cases of functional responses (type IV came later). A Holling type functional response is described as a predator's instantaneous per capita feeding rate and is a function of prey abundance. This means that the consumption rate of an individual predator depends on the prey density. The detail about each Holling type response function is found in Dawed et al. (2020).
Dynamical complexities of a three species food web with square root functional response is investigated where the prey species shows a herd behavior (Bera et al., 2016). Furthermore, discrete time-delay mathematical models for tri-trophic food chains have been developed and studied in recent past (Das & Samanta, 2020;Maiti et al., 2008;Maiti & Samanta, 2006;Manna et al., 2017;Pathak et al., 2009). The of effect discrete time delay has been investigated in these papers.
In this paper, we consider a food chain model (a top predator, an intermediate predator, and a source prey) with intermediate predator and prey harvesting. We observe the ecological interest and investigate the dynamics in the food chain and the impact of harvesting on the persistence of the model system. The necessary conditions for the permanent coexistence of the model, and criteria for local and global stability of the nonnegative equilibria are presented. Existence of optimal equilibrium solution of the control problem, maximum sustainable yield, and maximal net present revenue are determined. Numerical simulations using Matlab programming are performed to investigate the dynamic outcome of the system. This model shows complex behavior (such as oscillations), which reflects the real situation. These dynamic outcomes are supported by numerical simulations in the paper.
interaction with general Holling type functional response and harvesting prey and intermediate predator is given as follows: Descriptions of the state variables and parameters are given in Table 1. Here , called fishing effort (Fishing mortality) is introduced.
Proof. We can rewrite Equation (4) as Similarly, y τ ( ) and z τ ( ) are also positive.  ≥ ≥ ≥ , and the solution of the system initiating in the nonnegative octant is bounded. Therefore, the solution of the model system with nonnegative initial condition exists and is unique (Allen, 2007).

| Boundedness
Proposition 3.2. All solutions x τ y τ z τ ( ( ), ( ), ( )) of the model system (4)-(6) with positive initial condition x y z ( , , ) 0 0 0 are bounded in the region Proof. Boundedness argument for x: From the first equation of the model it is true that  Boundedness argument for z: Consider the sum relation m x κy z = + + κ β . After differentiation with respect τ it takes the form Using the relation z m x κy = − − κ β the differential equation reduces to the form as One can easily show that the right hand side function of the differential equation is strict concave with unique equilibrium point ( ) Therefore, the solution of the model system is bounded in the region Ω. □

| Existence of equilibrium points
The axial equilibrium points of the system (4)-(6) are E E h (0, 0, 0), (1 − , 0, 0) and E x y z ( *, *, *) xyz exist with some conditions. The equilibrium point E x y ( *, *, 0) xy exists as the prey and intermediate predator population grow to the difference of their carrying capacities and harvesting in the absence of top predator. To obtain suitable expressions for x* and y* of E xy , let both x and y be considered different from zero while z is zero. The equations now reduce to Considering y 0 ≠ , solve for y from (8), it can be obtained as y* = Substituting this expression for y* into (7) resembles that x* is a positive solution of the equation Therefore, if (9) has a positive solution, say x*, then the equilibrium point To show existence of positive solution for (9) at each Holling type functional response, we use the theorem below.
Theorem 3.3. Consider a polynomial of degree n, with n odd, in the form p x x a x a x a x a ( ) = (−1) The proof of the theorem is found in Dawed et al. (2020). Now we can show existence of positive solution for (9) for each Holling type functional response.
. Substituting this into (9) and solving for x*, . Substituting this into (9) and making arrangements, then x* is the solution of which exists (by Theorem 3.3), provided that ( ) The equilibrium point E xy bespeaks the prospect of prey and intermediate predators to coexist certain. The stability analysis results will provide conditions for such scenario to happen.
The equilibrium point in the yz plane say, E yz , is the solution of the system of equations . The existence of y* for the functional response HT -I-IV is given in Table 2.
We can show existence of positive solution for Equation (10) of each Holling type functional response as follows: Substituting this into Equation (10) and solving for x*, we get The value of y* for the functional response HT I-IV is given in Table 3. This implies that the interior equilibrium point exists when δh ε We also note that no equilibrium points of the form c d ( *, 0, *) and z (0, 0, *) exist. This is true as top predator feed solely on intermediate predator and do not have alternative source of food.

| Dynamical behavior
The dynamical behavior of the equilibrium can be studied by computing the Jacobian matrix corresponding to each equilibrium point. Observing the sign of the eigenvalues of this matrix the stability of an equilibrium point is determined (Allen, 2007). The Jacobian matrix J x y z ( , , ) of the system (4)-(6) is The three eigenvalues of this matrix are The trivial equilibrium point is saddle point which is unstable if h < 1. Biologically, this implies a purely extinction is not possible. However, if h > 1, and ε δh 0 < < hold, E 0 is stable, which implies biologically a pure extinction of the species. ( The three eigenvalues are otherwise it is a saddle point with stable manifold in the xz plane and unstable manifold in y direction. ( The Jacobian matrix at E y takes the form The three eigenvalues are The Jacobian matrix at E xy takes the form Here, we used notations Hence, conditions for local stability of the equilibrium point E xy are and it takes the form In this case the characteristic equation becomes Expanding each term we have The corresponding Jacobian matrix is given by . The eigenvalues say λ λ λ , , and  . We note that a purely extinction of all species is not possible (as E 0 is a saddle point). However, when the prey and intermediate predator are harvested more than their intrinsic grow rates (q E r > X 1 and q E r > Y 2 ) a purely extinction of all the species will happen, that is, E 0 is stable. We also note that prey can survive alone when the intermediate predator is harvested more than it grows, as the works by (Blé et al., 2018;Ghosh et al., 2018;Kumar & Agarwal, 2017). We remark that the coexistence of prey and intermediate predator, intermediate and top predators, and all the three species is possible. In addition, we note that the type and stability of the equilibrium points could change as parameters value vary. Furthermore, oscillators behavior and limit cycles could happen which will be investigated numerically.

| Global stability analysis
In this section we show global stability of the equilibrium points E E , Then, Hence, E x y ( *, *, 0) xy is globally asymptotically stable in the re- Similarly for (3b) E y z (0, *, *) yz is globally asymptotically stable in the region The Lyapunov theory is used to make conclusions about trajectories of a system (especially nonlinear) without finding the trajectories (i.e., solving the differential equation). The approach is based on the work (Upadhyay & Raw, 2011) and others.
Theorem 3.6. The equilibrium point E x y z ( *, *, *) xyz is globally asymptotically stable if Proof. We consider the positive definite Lyapunov function where w 1 and w 2 are positives constants that will choose later. Differentiating V with respect to time τ,

| Persistence
Persistence of a system means no state of the system approaches to zero, that is, there can be no extinction of any of the population that makes up the biological system. As it is explained, E (0, 0, 0) 0 is unstable in the x direction and asymptotically stable in the z direction, E h (1 − , 0, 0) x exists on the positive x axis. It is asymptotically stable in the x direction. An equilibrium point E (0, , 0) y ε δh ε − also exists and it is asymptotically stable in the y direction. No equilibrium point on the z axis.
In addition to the above assumptions the following hold: A 1 : The right hand sides of the system in (4)-(6) are C 1 in x y z ( , , ). A 2 : All solution of the system (4)-(6) with nonnegative initial conditions are bounded in forward time according to Proposition (3.2). A 3 : E x and E y exist and they are hyperbolic saddle points if ε φ h δh . Since E xy and E yz are locally asymptotically stable, and also they are globally asymptotically stable, there are no limit cycles around these points (Freedman & Waltman, 1984). These observations are written in theorem form as follows.
Theorem 3.7. The sufficient conditions for persistence of the model system (4)-(6)

| HARVESTING
In this section we discuss maximum sustainable yield, bionomic equilibrium, and harvesting optimal policy.

| Maximum sustainable yield
Assuming the equilibrium point x y z ( *, *, *) is locally asymptotically stable, the sustainable yield Y h ( ) during the equilibrium situation is given in Kot (2001) as where x* is a positive solution of x h x kφ x y y φ γ β + ( − 1) + ( ) * = 0, and * = .
x y 2 Here, we manipulate maximum sustainable yield for two Holling types functional response combinations: Case I, HT-I/HT-I: φ x αx ( ) = x , and φ y ωy ( ) = y . Substituting these into Equation (12)

| Bionomic equilibrium and harvesting optimal policy
The economic equilibrium means that the total revenue from harvesting of populations equals to the total cost of harvesting efforts together with biological equilibrium point x y z ( , , ). As in Abid et al. (2019), Chakraborty et al. (2012), Rojas-Palma andGonzález-Olivares (2012), Tchepmo Djomegni et al. (2019), we consider the total cost is proportional to the harvesting effort, written as T CE = c . We also assume that the total profit from the harvesting of population is proportional to the harvesting yield, written as T PY E = ( ) r , where P is the unit price of the stock and Y E qEX ( ) = denotes the harvesting yield function. The profit function denoted Π is then given by T T Π = − r c . Let C states harvesting cost per unit effort of population X and Y and P 1 and P 2 denotes the price per unit biomass of population X and Y , respectively. The profit function of harvesting for the intermediate predator and prey populations in the free fishing zone is given by: After scaling Equation (13)  , and c = r C q X 1 . A bionomic equilibrium of the model implies both biological and economic equilibrium. Biological equilibrium occurs when interior equilibrium point of the model system exists and economic equilibrium is achieved when the economic rent is completely degraded, that is, π h ( ) = 0. Thus, the bionomic equilibrium point x y z h ( , , , ) 0 0 0 0 satisfies the following conditions: provided that h 0 exists, h 0 is the critical threshold effort that determines the profitability of the investment in the harvest. Of course, h h = 0 means no profit. If h h > 0 , the total cost of the harvest will exceed the total selling price from the harvest. Then the investment will result in a loss and may be terminated. However, if h h < 0 , the investment is profitable. In this situation, policy must be drawn to protect the sustainability of species while maximum profits are targeted. Thus, we will define an optimal control problem and examine the optimal effort to apply where will produce maximum profit and secure continuity as it is explained in Tchepmo Djomegni et al. (2019).
We assume nonzero profit. Hence, the present profit is defined as J x y h e π x y h dτ ( , , ) = ( , , ) , δ 1 is the discount rate of the net revenue. We want to solve the following optimal control problem J x y h h τ h max ( , , ) subject to (4) − (6) and to 0 ( ) . max ≤ ≤ (21) Since the system Equations (4)-(6) and the objective functional J are linear in the control variable h, the optimal solution of the control problem Equation (21) is either a bang-bang or a singular control (Lenhart & Workman, 2007). But we will show that it cannot be a combination of a bang-bang and singular control which is similar to (Tchepmo Djomegni et al., 2019). For the singular control, the switching function is zero on nontrivial time intervals, where the Hamiltonian H is given in equation (24). Hence the optimal harvest policy (on the effort) will be (Lenhart & Workman, 2007) h τ The corresponding Hamiltonian function of the problem becomes where the functions λ λ , 1 2 and λ 3 are the adjoin variables. The condition Q τ ( ) = 0 implies.
Theorem 4.1. The optimal equilibrium solution x y z h ( *, *, *, *) of the control problem equation Proof. Standard results established the existence of the optimal control h* and the corresponding optimal steady state x y z h ( *, *, *, *) (Fleming & Rishel, 1975). The adjoint equations are Since we are looking for the optimal equilibrium solution, we treat x y z , , , and h as constants. Moreover, we consider the interior equilibrium x y z h ( *, *, *, *) (as the sustainability of species must be preserved), that is, solution of Equations (15)-(17). Hence the adjoint equations are reduced to Solving Equation (25) for λ 1 and substituting the result into Equation (31), we obtain λ a e a λ λ βφ y z where a p x δp y c a p δh εy. Solving Equation (32) for λ 2 and substituting into Equation (33), then λ 3 satisfies the second order ordinary differential equation.
This leads to the solution where ℵ 0 and ℵ 1 are constants of integration, Λ = , Λ = a a βφ y z a a βφ y z 1 We note that Λ + Λ = − Solving λ 2 from Equation (32) and then λ 1 from Equation (25), we obtain λ τ δ a e ( ) = , condition for E E , = (0.7, 0, 0) x x is locally asymptotically stable. Stability of E x is also shown in Figure 2b. , the equilibrium point E xy is stable. In addition, for ε 0.25 ≥ , we observe that the interior equilibrium point E xyz is asymptotically stable. This implies the three species coexist in