Optimal harvesting policy of a prey–predator model with Crowley–Martin-type functional response and stage structure in the predator

Abstract. In this paper, a three-dimensional dynamical model consisting of a prey, a mature predator and an immature predator is proposed and analysed. The interaction between prey and mature predator is assumed to be of the Crowley–Martin type, and both the prey and mature predator are harvested according to catch-per-unit-effort (CPUE) hypothesis. Steady state of the system is obtained, stability analysis (local and global both) are discussed to explore the longtime behaviour of the system. The optimal harvesting policy is also discussed with the help of Pontryagin’s maximum principle. The harvesting effort is taken as an effective control instrument to preserve prey and predator and to maintain them at an optimal level.


Introduction
The Lotka-Volterra system of equations was established and analyzed long time back, about 100 years ago. These equations are a mathematical and dynamical model representing the relationship between two or more species. Several attempts have been made to generalize, modify and extend these equations. However, due to complex nature of the biological species, their complete dynamics is still not known and needs to be investigated with care. It has now been established that age plays an important role in deciding the dynamics and evolution of various species. The rates of reproduction and survival largely depend upon age or the developmental stage, and hence it could be remarked that the life history of several species is composed of at least two stages, immature or juvenile and mature or adult, with significantly different biological, physiological and morphological characteristics.
The analysis of stage-structured predator-prey system has attracted good amount of attention recently, as a way to eliminate the shortcomings of classical Lotka-Volterra models [1,2,11,22,28,36]. In these models, a time delay represents the age of maturity of the species. In fishery system, cannibalism has been observed, and various types of cannibalism models have been discussed [6,16,26]. Recently, Bosch and Gabriel [4] and Kar [19] studied the stage and age structure of species without or with time delays.
One of the major aims of ecologists is to gain insight into predator-prey relationship, and one vital aspect of predator-prey relationship is the rate of predation by an average consumer (this is known as the functional response or the "trophic function"). The functional response takes into account of both the predator and prey biological and physiological processes. The functional response largely controls the stability of the system, and they are of several types: Holling I-III, ratio-dependent, Beddington-DeAngelis, Crowley-Martin, Leslie-Gower [3,17,24,25,30,34,37]. There are a few literatures available on predator-prey model with Crowley-Martin-type (CM-type) functional response [27,33,34]. The CM functional response involves the interference among individual of predators engaged in handling or searching the prey, and it is given by where b, a 1 and b 1 are positive parameters that are used for effects of capture rate, handling time and magnitude of interference among predators, respectively. The CM-type functional response reduces to classical and Holling-type II functional response under the following constraints on the parameters: (i) a 1 = 0, b 1 = 0 implies linear mass action (classical Lotka-Volterra) functional response. (ii) a 1 > 0, b 1 = 0 implies Holling-type II (Michaelis-Menten) functional response. (iii) a 1 = 0, b 1 > 0 implies Holling-type II (saturation with respect to predator) functional response.
The effect of intra-specific interference among predators has been studied in a preypredator model with Holling-type II functional response in [13], with Holling-type III functional response in [12] and with Beddington-DeAngelis-type functional response in [15]. In these three studies, spatiotemporal dynamics of the system are also investigated. Guin et al. [14] have also studied spatiotemporal pattern in a prey-predator model with prey refuge and Beddington-DeAngelis-type functional response.
The optimal management and utilization of renewable and natural resources, which is directly related to sustainable development, has been studied extensively by many authors [7-10, 18-20, 23]. Recently, Maiti et al. [27] investigated the dynamics of a preypredator model with CM-type functional response with refuge for the prey species. To the best of authors' knowledge, optimal harvesting of prey-predator with CM-type functional response and with stage structure in predator population has not been studied. Keeping these in view, we propose a three-dimensional model consisting of prey and predator in which predator is divided into two categories: mature and immature. The prey and mature predator are harvested as CPUE hypothesis. The rest of the paper is organized as follows: In Section 2, we formulate the mathematical model and its qualitative properties. Section 3 deals with the existence of all feasible equilibria, and stability analysis is presented in Section 4. Optimal harvesting policy is discussed in Section 5, and numerical simulations are presented in Section 6. Finally, concluding results are given in Section 7 followed by references.

Mathematical model and its qualitative properties
We consider a habitat consisting of a prey and predator system. We assume that the density of prey population or the renewable resource under consideration, represented by x(t) at any time t 0, can be mathematically and dynamically modelled by a logistic equation when the predator is absent. We assume that the predators are classified into two stage groups -mature or adult and immature or juveniles, and their densities are denoted by y(t) and z(t), respectively, at any time t 0. Here we are assuming the fact that only mature predators are capable of attacking the prey and have reproductive ability, while the immature predator does not attack the prey and has no reproductive ability. A good example of such a situation is the case of the Chinese alligator, which can be regarded as a stage structured species since the mature is more than 10 years old, and can be regarded as a predator because almost all the aquatic animals are the chief nutritional source for the alligator. The interference between prey and adult predator is assumed to be of the CM-type. One of the novel features about our model is to account for the universally prevalent intra-specific competition in the consumer growth dynamics [21]. This intraspecific competition is assumed to induce an additional increased death rate, which is proportional to the square of the adult population [12,13,15]. We assume that prey and adult predator are harvested as CPUE hypothesis, and juvenile predators are not harvested. With these assumptions in mind, we propose the following stage-structured prey-predator interaction model: Here r is the specific growth rate of the prey, and K is the carrying capacity. The predator functional response incorporated is the CM type, where α, a and b are positive parameters that are used for the effects of capture rate, handling time and magnitude of interference among predators, respectively. The parameter c is the conversion factor, q 1 is the catchability coefficient of the prey, q 2 is the catchability coefficient of the mature predator or adult species, E is the harvesting effort, δ 0 is the death rate of the matured predator and δ 1 is intra-specific interference coefficient of the adult predator. The parameter β denotes the birth rate of immature predator, β 0 denotes the death rate of the immature predator, and β 1 denotes the proportionality constant of transformation of immature to mature predators. Remark 1. If E = 0, δ 1 = 0 and β 1 = 0, then dynamics of model (1) is well studied in [31].
Next, we present some qualitative properties of our proposed model to show that the model is biologically well behaved.
From the first equation of the model

Now suppose
then we have Hence, it follows that dW dt We also note that if x K and W (t) rK/(4δ m ), then dx/dt 0, dW /dt 0. This shows that all solutions of system (1) starting in Ω remain in Ω for all t > 0.
Theorem 2. Let the following inequalities are satisfied: Then the model system (1) is uniformly persistence, where x m is defined in the proof.
Proof. Permanence or uniform persistence of a system implies that all species will be present in future and none of them will become extinct if they are initially present. System (1) is said to be uniformly persistence if there are positive constants M 1 and M 2 such that each positive solution X(t) = (x(t), y(t), z(t)) of the system with positive initial conditions satisfies Keeping the above in view, if we define This also shows that for any sufficiently small > 0, there exists a T > 0 such that for all t T , the following holds: Now from the first equation of model system (1), for all t T , we can write Hence, it follows that which is true for every > 0, thus Similarly, the third equation of model system (1) yields Taking M 1 = min{x m , y m , z m }, the theorem follows. https://www.mii.vu.lt/NA

Analysis of equilibria
It can be inspected that model (1) has four nonnegative equilibria, namely, P 0 (0, 0, 0), P 1 (x, 0, 0), P 2 (0,ȳ,z), P * (x * , y * , z * ). The equilibrium point P 0 exists obviously. We shall show the existence of the other equilibria as follows: Existence of P 1 . Herex is the positive solution of the following equation: Clearly,x > 0 if the following inequality holds: Thus, the equilibrium P 1 exists under condition (2). Existence of P 2 . Hereȳ andz are the positive solutions of the following equations: From Eqs. (3) and (4) we obtain We note that forȳ andz to be positive, we must have Thus, P 2 exists if inequality (5) holds true. Existence of P * . Here x * , y * and z * are the positive solutions of the following algebraic equations: It is easy to note that if we are able to verify the existence of x * and y * , then existence of z * automatically follows from Eq. (8). We perform the following analysis to show the existence of x * and y * .
From Eq. (6) we note the following: We note that y a > 0 if, in addition to (2), the following inequality holds: (ii) When y = 0, It is easy to see that dy/dx < 0 if the following inequality holds: The above analysis shows that isocline (6) is passing through the points (x a , 0) and (0, y a ); and in Eq. (6), y is a decreasing function of x under conditions (2), (9) and (10). Now we note the following from Eq. (7): This shows that isocline (7) is passing through the point (0, y b ) under condition (5), and it has always a positive slope, thus in Eq. (7), y increases as x increases.
From the above analysis we infer that two isoclines (6) and (7) intersect at a unique point ( Now we are in a position to state the following theorem. Theorem 3. The positive equilibrium P * (x * , y * , z * ) exists and it is unique if conditions (2), (5), (9), (10) and (11) hold true.

Stability analysis
The local stability of each equilibria can be studied by computing the corresponding Jacobian matrix. We note the following regarding the linear stability behavior of these equilibria: 1. P 0 is a saddle point. This follows from the following remarks: • The eigenvalue corresponding to the x-direction is r − q 1 E, which is positive from condition (2).
• Since ββ 1 > (δ 0 +q 2 E)(β 0 +β 1 ) from condition (5), the product of eigenvalues corresponding to the yand z-directions is negative. This, in turn, implies that the equilibrium point P 0 is locally stable only in one direction (either yor z-direction) and is unstable in a two-dimensional space.
2. P 1 is also a saddle point. This follows from the following remarks: • The eigenvalue corresponding to the x-direction is equal to −(r − q 1 E), which is negative from condition (2). • The product of the eigenvalues corresponding to the yand z-directions is given by the following expression: This expression is clearly negative under conditions (2) and (5). Therefore, the equilibrium point P 1 is locally stable in a two-dimensional space and is unstable in a single direction, which is either y-direction or z-direction.
3. The following analysis discusses the stability of P 2 : • P 2 is locally stable or unstable in x-direction depending upon the condition whether r − q 1 E < αȳ/(1 + bȳ) or r − q 1 E > αȳ/(1 + bȳ) holds true, respectively. • The product and sum of the eigenvalues corresponding to the yand z-directions, respectively, is given by the following expressions: respectively.
Expression (12) is positive by (5), which implies that the product of the eigenvalues is positive. The sum of the eigenvalues, that is, expression (13) is clearly negative. The above statements imply that both the eigenvalues are negative. Therefore, the equilibrium point P 2 is locally stable in the two-dimensional space spanned by the unit vectors pointing in yand z-directions, respectively. • Hence, the equilibrium point P 2 is locally stable or a saddle point depending upon the condition whether r − q 1 E < αȳ/(1 + bȳ) or r − q 1 E > αȳ/(1 + bȳ) holds true, respectively.
4. We use the Routh-Hurwitz criterion to study the stability behavior of P * . The Jacobian matrix evaluated at P * is given by The characteristic equation corresponding to the above Jacobian matrix is Hence, P * is locally asymptotically stable under conditions (14).

Remark 2.
It has been noted that all inequalities in Eq. (14) are satisfied if holds. Hence, P * is locally asymptotically stable under condition (15).
We will now prove that P * is globally asymptotically stable under certain conditions in the next theorem.
Theorem 4. Let the following inequalities hold: Then P * is globally asymptotically stable in Ω with respect to all the solutions initiating in the interior of the positive octant.
Proof. Consider the following positive definite function about P * : where K 1 and K 2 are positive constants to be chosen suitably in the subsequent steps. Differentiating V with respect to t along the solutions of model (1), a little algebraic manipulation yields Choosing K 1 = (1 + ax * )/(c(1 + by * )) and K 2 = K 1 β 1 /(βy * ), we note that dV /dt is negative definite under conditions (16) and (17). Hence, V is a Liapunov function with respect to all the solutions initiating in the interior of the positive octant, proving the theorem.
The above theorem implies that under parametric conditions (16) and (17), the predator (juvenile and adult both) and prey densities settle down at their interior equilibrium point irrespective of the initial values of their densities at t = 0.
Remark 3. As long as P * exists, it is interesting to note here that condition (16) implies condition (15) because 1 + ax * > 1 holds.

Optimal harvesting policy
The exploitation of biological resources is commonly practiced in fishery, forestry and wildlife management. A management for biological species such as fishery is needed to maintain an ecological balance, which is disrupted due to overexploitation of these renewable resources. Keeping this in mind, we discuss the optimal harvesting policy that is to be adopted by the regulatory agency so as to maximize the total discounted net revenue obtained from harvesting prey and predator species using harvesting effort as the control instrument. We wish to investigate the 3D curve (x, y, z) with the optimal harvesting effort E so that the system remains at an optimal equilibrium level.
The net economic revenue to the society π(x, y, z, E, t) = net economic revenue to the harvester + net economic revenue to the regulatory agency where c is the harvesting cost per unit effort, which in turn is given by c = c 1 + c 2 , where c 1 is the harvesting cost per unit effort corresponding to the prey species, and c 2 is the harvesting cost per unit effort corresponding to the adult predator species. p 1 is the price per unit biomass of x, and p 2 is the price per unit biomass of y. We take p 1 , p 2 and c to be positive constants.
Our problem is to optimize the objective functional We construct the Hamiltonian where λ 1 , λ 2 and λ 3 are the adjoint variables, E is the control variable subject to the constraint: 0 E E max . Here E max denotes a feasible upper limit of E subject to the infrastructural support available to fishing/harvesting. Suppose E is the optimal control, and x, y, z are the responses. By the maximum principle, for t 0, there exist adjoint variables λ 1 , λ 2 and λ 3 such that At the equilibrium point P * , the above equations reduce to This system of linear differential equations can be solved using the operator method by eliminating λ 2 and λ 3 . Then the reduced differential equation in λ 1 can be written as where The solution of Eq. (19) is where A i s (i = 1, 2, 3) are arbitrary constants, and α i 's (i = 1, 2, 3) are the roots of the auxiliary equation a 3 m 3 + a 2 m 2 + a 1 m + a 0 = 0, It is clear from (21) that λ 1 is bounded if and only if α i < 0 (i = 1, 2, 3) or A i (i = 1, 2, 3) are identical to zero. For robust calculations, we ignore the cases where α i < 0 (i = 1, 2, 3) and take A i (i = 1, 2, 3) identical to zero. Then we have Proceeding in a similar fashion, we obtain Thus, the shadow prices e δt λ i (i = 1, 2, 3) remain constant over time in optimal equilibrium when they satisfy the transversality condition at t = ∞, i.e., when they remain bounded as t → ∞. From (18) we note that Hamiltonian H is linear in the control variable E. Hence, optimal control will be a combination of the bang-bang control and singular control. A necessary condition for singular control to be optimal [5] is that which gives Therefore, we may conclude that the total cost of harvest per unit effort (the left-hand side of (22)) must be equal to the discounted value of the marginal profit of the static effort (the right-hand side of (22)) level.
Substituting the values of λ 1 and λ 2 in (22), we obtain The above equation together with the value of E at the interior equilibrium, namely, gives the optimal equilibrium population x = x δ and y = y δ . When δ → ∞, we have M 1 /N → 0, M 2 /N → 0. Then (23) reduces to and hence π(x ∞ , y ∞ , z, E) = 0. This shows that the economic rent is completely dissipated when the discount rate is infinite. The economic rent can be expressed as We note that M 1 is of O(δ 2 ), M 2 and M 3 are O(δ 3 ) and N is of O(δ 3 ). Thus, π is a decreasing function of δ.

Numerical results
In the following section, we present some numerical simulations to verify our theoretical results proved in the previous sections by using MATLAB R2010a.
For the above set of values of the parameters, conditions in Theorem 4 for the existence of the interior equilibrium is satisfied. Thus, the positive equilibrium point P * (x * , y * , z * ) is given by We also note that condition (14) is satisfied for the set of parameters chosen in (25). Thus, the equilibrium point P * (x * , y * , z * ) is locally asymptotically stable. The time series of x, y and z are presented in Fig. 1. This figure shows that the density of the prey species decreases with time whereas densities of predator species (mature and immature both) increase with time and finally settle down at their steady states. It is also observed here that the density of immature predator settles at a larger value than that of the mature predator and prey. It may be pointed out that the values of parameters chosen in (25) satisfy local stability conditions, but they do not satisfy global stability conditions. Since conditions obtained in Theorem 4 are sufficient (not necessary) for the global stability of P * , hence at this stage, we cannot say anything about the global stability of P * . Now we choose following set of values of parameters: r = 7, K = 10, c = 1, α = 0.5, with different initial conditions. These values of the parameters satisfy the global stability conditions of Theorem 4. The trajectories or solution curves of x, y and z with different initial conditions are plotted in Fig. 2. From this figure we note that all the trajectories starting from the different initial conditions converge to the equilibrium point P * (8.568, 0.009763, 0.003665). This shows that P * is globally asymptotically stable. Figure 3 shows the behavior of x, y and z for different values of the parameter α. Here the rest of parameters have the same values as in (25). We note that if α (capture rate) is small, then all the three species grow and finally attain its respective steady states. If α increases beyond a critical value, then prey population decreases, mature and immature populations increase. If α becomes high, then the density of prey species tends to zero. Figure 4 shows the behavior of x, y and z for different values of the parameter c. Here, again, the rest of parameters have the same values as in (25). In this case, it can be noted that as the value of c increases, the densities of y and z increase, but the density of x decreases.     For the optimal harvesting part, we choose the following set of parameters: r = 7, K = 10, c = 0.1, α = 0.5, Solving (23) and (24) simultaneously, we get the optimal values as x δ = 1.001, y δ = 2.3530, z δ = 2.94125, and the optimal value of E is given by E δ = 5.4339. This value of E is optimal in the sense that for such a value, the harvesting agency gets the maximum revenue for the harvest, and all the three species will coexist at an optimal level. From Fig. 5 one can remark that when the value of E is below E δ , the prey and predators (mature and immature both) survive, but when the value of E is above E δ (case of over harvesting), then the population densities of prey and mature predators tend to zero.
In system (1), let b = 0.1 and the rest of parameters are same as that in (28). Then it is easy to note that system (1) has a unique interior equilibrium E * (27.7097, 2.3121, 2.0552), which is globally asymptotically stable as conditions of Theorem 4 are satisfied. Figures 6(a)-6(c) show the time series analysis, and Fig. 6(d) shows the phase portrait analysis of model system (1) for different values of the parameter b. These figures shows that system is stable if b is small, and if b increases beyond a threshold value, the system becomes unstable. Thus, the parameter b induces a Hopf-bifurcation in the system.
In order to consider the importance of parameter δ 1 , let δ 1 = 2, and the rest of parameters are same as that in (28). Then system (1) has unique positive equilibrium point E * (21.7536, 1.3268, 1.1793), and conditions of Theorem 4 are also satisfied. Hence, the (d) Figure 6. Trajectory of x, y, z and limit cycle with respect to the parameter b. Other parameters are as in (28).

Conclusions
The proposed model consists of three nonlinear differential equations, namely, one for mature predator, one for immature predator and one for the prey. Only the mature predator feeds on the prey, immature predator survives via mature predator and some alternative food. The interaction between prey and mature predator has been taken as the Crowley-Martin type, which is more realistic in nature. For ecological balance, it has been modeled in the system that only prey and mature predators are harvested, while the immature predators are not harvested. An interesting aspect in mathematical ecology is permanence/persistence, which ensures the survival of biological entities for all positive initial conditions. If a system exhibits permanence, then the ecological planning on fixed eventual population can be carried out. Analyzing the system, we have obtained some constraints on the intrinsic growth rate of the prey species for the permanence of the solutions of our system. It has been shown that all solutions of the system are positive and bounded if all the species are initially present. Thus, our proposed model is biologically well behaved. The dynamical modeling of the system's behavior shows that the system under consideration is locally stable around positive interior equilibrium. Also, it has been observed that the system around the positive equilibrium is globally asymptotically stable under certain conditions. We have studied the optimal harvesting policy using the Pontryagin's maximum principle. For economic and biological views of renewable resource management, we studied the exploitation of both prey and adult predators. From the point of view of ecological management, in order to plan harvesting strategies and keep sustainable development of ecosystem, we have used the harvesting effort as a control parameter and obtained its optimal level E δ . If applied effort is less that E δ , all the species will coexist at an optimal level, and ecological balance can be maintained. If applied effort is larger than E δ , then it represents over-exploitation and the prey-predator system will be in the danger of extinction. Our numerical simulation results obtained in Figs. 6-8 show that the parameter b (magnitude of interference among predators) and δ 1 (intra-specific interference among adult predators) play an important role in governing the dynamics of the system. We hope that this study will help to understand the dynamics of prey-predator system with harvesting. In this paper, we have considered the harvesting effort as a constant control variable. This effort may be a dynamic variable, and the taxation policy imposed by a regulatory agency can be thought of as a control variable. This work we leave as a future research.