1 Introduction

Ecosystem is a system formed by the ecological community and its environment as a unit. Basically, Ecology is a branch of biology which deals with the relationships among organisms and the relationship between them and their surroundings. In order to survive, all living organisms live mutually and affect each other in several ways. Hence, the fundamental concerns in ecology is the stability of ecological systems and the persistence of species within them.

The elementary unit of ecosystem is the predator–prey models in which biomass are grown from their resource masses. The dynamical relationship between predator and their prey is from since long and will remain a significant theme in ecology and mathematical ecology because of its importance and universal existence. The simplest Predator–prey dynamic model is Lotka–Volterra Predator–prey model [1]. The model is modified by several mathematician and ecologist after its original formulation in 1920s. Many researchers have done study in different fields depending on how a species is affected by the other. Therefore, the dynamical interaction between predator and prey plays an vital role in the ecosystem. Predation is a biological interaction in which a predator organism feeds on another living or other organisms which are known as prey. The functional response of predator to prey density refers to the change in prey density per unit when prey is attacked by the predator. Based on this, Holling et al. [2] proposed three type of functional response of the predator which are usually called as Holling type-I, Holling type-II and Holling type-III. The mathematicians and ecologists have a great interest in the Holling type predator–prey models. A number of response functions also exist in the literature which are known as Beddington DeAngelis type response function, Hassel Verley type response function etc.

In the linear Holling type-I functional response, the time taken by a predator in search of prey is ignored. The more realistic functional response is defined as

$$\begin{aligned} f(x)=\frac{mx}{a+x} \end{aligned}$$
(1)

where the positive parameter m and a represented as maximum growth rate of species and half saturation constant. The function (1) is called Holling type-II functional response. Another functional response is Holling type-III, which is defined as

$$\begin{aligned} f(x)= \dfrac{mx^{2}}{a+x^2} \end{aligned}$$

Andrews [3] suggested a function

$$\begin{aligned} f(x)=\dfrac{mx}{a + bx + x^2} \end{aligned}$$

which is called the Monod–Haldane function for low concentrations but includes the inhibitory effect at high concentrations as well. Sugui and Howell [4] suggested Monod–Haldane functional response of the form

$$\begin{aligned} f(x) = \dfrac{x}{bx^2 + m} \end{aligned}$$

where \(b=\frac{1}{i}\) is taken as the inverse measure of inhibitory effect. The above functional response describes the group defense phenomenon which illustrates that when the prey population is significantly large and their defending ability is increased then predation always reduced. Tener et al. [5] suggests one example of group phenomenon. Wolves can attack Lone musk ox successfully. Small herds of musk ox are attacked but with less chance of success. Success rate of attack in large herds is very low. The Monod–Haldane functional response have rich meaningful dynamics [6, 7] due to which researchers are attracted towards it for further study of Monod–Haldane type Models. Lui and Tan [8] proposed Monod–Haldane functional response with impulsive harvesting and stocking. Zhang et al. [9] described non-monotonic functional response for studying group defense capability. Lui et al. [10] studied the existence of periodic cycles in a predator–prey system equipped with group defense mechanism of prey species. The toxic phytoplankton and zooplankton population interaction has been studied by Rimpi et al. [11] using simplified Monod–Haldane type functional response.

The study of predator–prey system has been investigated by many researchers [12,13,14,15,16]. In the recent research, chaotic dynamics is also found by the researchers. Chaotic dynamical behaviour in tri-tophic food chain model has been studied by Upadhyay and Raw [14]. Hu and Cao [15] studied the dynamics of Predator–prey system using nonlinear Michaelis-Menten type harvesting where as Shen et al. [16] investigated canard limit cycles. Mena-Lorca et al. [17] investigated the dynamical behavior of Leslie–Gower predator–prey model with proportional harvesting. Lin and Ho [18], proposed a modified Leslie–Gower type predator–prey model using Holling type-II functional response with time delay. They studied the local and global dynamics of this system. Li and Xiao [19] studied Leslie–Gower predator–prey model using Holling type- III functional response and observed a rich and complex dynamics of the system. Song and Li [20] analyzed and proposed the periodic Modified Leslie–Gower predator–prey model with Holling type-II scheme and they have investigated the impulsive effect of this system. Zhu and Lan [21] have investigated a Leslie–Gower predator prey model with constant harvesting rate in prey. They observed that the interior state can be saddle, stable and unstable, saddle-node, under certain parametric conditions. Zhang et al. [22] investigated proportional harvesting of Leslie–Gower predator–prey model. The dynamics of Leslie–Gower model with proportional harvesting subjected to Allee effect has been investigated by Rojas-Palma and Gonzalez [23]. The prey is considered growing logistically and predator is assumed to follow modified Leslie–Gower type predation. The system exhibits complex dynamical behavior including several local and global bifurcations.

It is very crucial to interpret the effect of one species on the other, in the presence of some process as harvesting of some specific species [24,25,26,27,28,29,30]. In case of interactions, the analysis of harvesting of any species becomes important. Clark [24, 25] took the initiative to study the impact of harvesting. Song and Chen [28] proposed the optimal harvesting of two species system. Kar and Chaudhary [29] studied the harvesting of one or more species in a dynamical system of three species. Different kind of harvesting functions such as constant, linear, proportional and non-linear harvesting are described and studied in literature. Observation indicates that non linear harvesting [27] provides better accuracy in comparison to other possible harvesting function. Linear harvesting function indicating that these function are not suitable for large population and the results indicates that these functions are not adapted for a large population [24,25,26]. The harvesting function discussed by the author [13] is not linear in control variable. Feng et al. [30] reported the quadratic harvesting function to prevent these two unusual aspect. Motivated by the above literature survey, Monod–Haldane type functional response is used to study proposed predator–prey model.

In the present work, the model discussed by the author [31] is modified by taking quadrating harvesting of prey. Therefore, we have considered a predator–prey model with Monod–Haldane functional response, where prey is subjected to quadratic harvesting. In this system, modified Leslie–Gower type predation is assumed for predation. The structure of manuscript is as follows:

In Sect. 2, model formulation is presented. Section 3 illustrates the positivity, boundedness and permanence of the system. Section 4 describes the existence of all possible equilibria of model system. Section 5 represents the local stability and local bifurcation analysis of model system. The optimal harvesting of system is derived in the Sect. 7. In the Sect. 8, numerical simulations are investigated for a suitable set of data to validate the analytical results. The Sect. 9 gives the description of the result obtained and analysis of the study.

2 The model formulation

A predator–prey model with modified Monod–Haldane functional response is considered. Let N(T) and P(T) be the population density of prey and predator at time T. In proposed model, prey population N(T) is considered logistically growing in the absence of predator P(T) with intrinsic growth rate r and carrying capacity K such that

$$\begin{aligned} \dfrac{dN}{dT}=rN\bigg (1-\frac{N}{K}\bigg ) \end{aligned}$$
(2)

The interaction between prey and predator is assumed to follow modified Monod–Haldane functional response. The density of prey takes the following form in the presence of predator:

$$\begin{aligned} \dfrac{dN}{dT}=rN\bigg (1-\frac{N}{K}\bigg )-\dfrac{m_1NP}{a_1+N^2} \end{aligned}$$
(3)

where the parameters \(m_1\) and \(a_1\) represent the catch rate and half saturation rate of the prey w.r.t. the predator.

The density of predator population is assumed to follow the modified Leslie–Gower predation as follows:

$$\begin{aligned} \dfrac{dP}{dT}=sP\bigg (1-\dfrac{m_2P}{a_2+N}\bigg ) \end{aligned}$$
(4)

where parameters \(m_2\) and \(a_2\) are used for the maximum growth rate and half saturation value of the predator, respectively.

In the proposed model, prey is subjected to follow quadratic harvesting with harvesting function \(h_1EN^2\). The constants \(h_1\) and E are known as catch-ability coefficient and total effort, respectively. Combining the Eqs. (3) and (4) with harvesting of prey, the following set of dynamical differential equations:

$$\begin{aligned} \dfrac{dN}{dT}= & {} rN\bigg (1-\frac{N}{K}\bigg )-\dfrac{m_1NP}{a_1+N^2}-h_1EN^2 \nonumber \\ \dfrac{dP}{dT}= & {} sP\bigg (1-\dfrac{m_2P}{a_2+N}\bigg ) \end{aligned}$$
(5)

The model can be non-dimensionalized with the transformations \(x = \dfrac{N}{K}, y = \dfrac{m_1P}{rK^2}, t = rT\), using following dimensionless parameters:

$$\begin{aligned} A_1=\dfrac{a_1}{K^2}, \delta = \dfrac{r}{s}, \beta =\dfrac{sm_2K}{m_1}, A_2=\dfrac{a_2}{K}, h= \dfrac{h_1EK}{r} \end{aligned}$$

The system will take the form as:

$$\begin{aligned} \frac{dx}{dt}= & {} x(1-x)-\dfrac{xy}{A_1 +x^2}-hx^2=x f(x,y)\nonumber \\ \frac{dy}{dt}= & {} y\bigg (\delta -\dfrac{\beta y}{A_2+x} \bigg )= y g(x,y) \end{aligned}$$
(6)

with initial conditions \(x(0)=x_0>0, y(0)=y_0>0\). All the constants \( A_1,h,\delta ,\beta ,A_2 \) of the system (6) are positive. It is easily observed that the functions xf(x, y) and yg(x, y) of the system (6) are continuous in positive octant of \({\mathbb {R}}=\{(x,y):x(0)>0, y(0)>0\}\).

It can be seen that the functions xf(x, y) and yg(x, y) are Lipschitz functions on \({R_2}^+\). Therefore, the solution of the system exists which is unique.

3 Positive invariance and boundedness

3.1 Positive invariance

For \(X(0)=X_0\in {\mathbb {R}}_+^2\), where \(X=(x,y)\in {\mathbb {R}}_+^2\), the system (6) can be formulated in a matrix form as \(\dot{X}=F(X)\) given by

$$\begin{aligned} F(X) = \left[ \begin{array}{cc} x(1-x)-\dfrac{ xy}{A_1 +x^2}-hx^2\\ y\bigg (\delta -\dfrac{\beta y}{A_2+x} \bigg )\end{array} \right] \end{aligned}$$

where \(F:C_+\rightarrow {\mathbb {R}}^2\) and \(F\in C^\infty ( {\mathbb {R}}^2)\).

It is observed that \(F_i(X)\mid _{{X_i}=0}\ge 0\) (for \(i=1,2\)) where \(X(0)\in {\mathbb {R}}_+^2\) such as \(X_i=0\). Therefore, the system (6) is positively invariant in \({\mathbb {R}}_+^2\) \(\forall \) \(t> 0\) (Nagumo [32]).

3.2 Boundedness

Theorem 3.1

All the solution of the system (6) are uniformly bounded.

Proof

Using system (6), the following condition is obtained

$$\begin{aligned} \frac{dx}{dt}<x(1-x) \end{aligned}$$
(7)

having initial condition \(x(0)= x_0>0\).

This shows that the solutions of defined system (6) must satisfy the condition \(x(t)\le 1, \quad \forall t>0\).

Consider a function \(\phi (t)\) such that

$$\begin{aligned} \phi (t)= & {} x(t)+y(t) \\ \dfrac{d\phi (t)}{dt}= & {} x^{\prime }(t)+y^{\prime }(t)\\= & {} x(1-x)- \dfrac{x y}{A_1+x^2}-hx^2+y \bigg ( \delta -\dfrac{\beta y}{A_2+x}\bigg )\\< & {} x(1-x)-\dfrac{x y}{A_1+x^2}+y \bigg ( \delta -\dfrac{\beta y}{A_2+1}\bigg )\\< & {} x(1-x)+\delta y- \bigg ( \dfrac{\beta y^2}{A_2+1}\bigg ) \end{aligned}$$

Assuming a positive constant K and rewriting the above equations as follows:

$$\begin{aligned} \dfrac{d\phi (t)}{dt}+K\phi (t)\le & {} - ( x^2-x )- \frac{\beta }{A_2+1} \bigg ( y^2- \frac{\delta y(A_2+1)}{\beta } \bigg ) +Kx+Ky\\\le & {} - ( x^2-(1+K)x)- \frac{\beta }{A_2+1} \bigg ( y^2- \frac{(\delta +K)(A_2+1) y}{\beta } \bigg ) \end{aligned}$$

and the further simplification gives,

$$\begin{aligned} \dfrac{d\psi (t)}{dt}+K\psi (t)\le & {} -\bigg (x-\frac{1+K}{2}\bigg )^2-\dfrac{\beta }{A_2+1}\bigg (y-\frac{(\delta +K)(A_2+1)}{2\beta }\bigg )^2\\&+\dfrac{(A_2+1)^2(\delta +K)^2}{4\beta ^2}+\dfrac{(1+K)^2}{4}\\\le & {} \dfrac{(A_2+1)^2(\delta +K)^2}{4\beta ^2}+\dfrac{(1+K)^2}{4} \\\le & {} N; \quad \quad N= \dfrac{(A_2+1)^2(\delta +K)^2}{4\beta ^2}+\dfrac{(1+K)^2}{4} \end{aligned}$$

The following inequality is obtained after solving the above differential,

$$\begin{aligned} \psi (t)\le & {} \frac{N}{K}\bigg (1-e^{-Kt}\bigg )+\psi (0)e^{-Kt} \\ 0< & {} \limsup _{t \rightarrow \infty }\psi (t)\le \dfrac{N}{K} \end{aligned}$$

Which gives

$$\begin{aligned} S=\bigg \{(x,y)\in {\mathbf {R}}; 0<x+y \le \dfrac{N}{K}+\psi \text { for any }\psi >0 \bigg \} \end{aligned}$$

This proves that all solutions are bounded. \(\square \)

3.3 Permanence

Definition 3.1

A system (6) is said to be permanent if there is a possibility of getting some positive constants \(k_1\), \(k_2\) and \(K_1\), \(K_2\) such that each positive solution of system (6) with some initial condition \((x_0, y_0)\in ({\mathbb {R}}_+^2)\) satisfies:

$$\begin{aligned}&k_1\le \liminf _{t \rightarrow \infty }x(t,x_0,y_0)\le \limsup _{t \rightarrow \infty }x (t,x_0,y_0)\le K_1\\&k_2\le \liminf _{t \rightarrow \infty }y(t,x_0,y_0)\le \limsup _{t \rightarrow \infty } y(t,x_0,y_0)\le K_2 \end{aligned}$$

Theorem 3.2

The system (6) with initial condition \((x_0, y_0)\) is permanent iff

$$\begin{aligned} h<\dfrac{A_1-M_2}{M_1 A_1} \end{aligned}$$
(8)

Proof

From the first equation of (6), we obtained:

$$\begin{aligned} \frac{dx}{dt}&\le x(1-x)\nonumber \\ \limsup _{t \rightarrow \infty }x(t)&\le max\{x(0), 1\}=1\equiv K_1 \end{aligned}$$
(9)

Similarly, the second equation of (6) provides:

$$\begin{aligned} \dfrac{dy}{dt}&\le y\bigg (\delta -\frac{\beta y}{A_2+L_1}\bigg )\nonumber \\&\le \beta y(A_2+L_1)\bigg (\dfrac{\delta (A_2+L_1)}{\beta }-y\bigg ) \nonumber \\ \limsup _{t \rightarrow \infty }y(t)&\le max\bigg \{y(0), \dfrac{\delta (A_2+L_1)}{\beta }\bigg \}\equiv K_2 \end{aligned}$$
(10)

Again, from first equation, we observe:

$$\begin{aligned} \dfrac{dx}{dt}&=x\bigg (\big (1-x \big )-\frac{y}{A_1+x^2}-h x \bigg ) \\&\ge x\bigg (\big (1-x \big )- \frac{y}{A_1}-h x\bigg ) \\&\ge x\bigg (\bigg ( 1-h M_1- \frac{M_2}{A_1} \bigg )-x\bigg ) \\&\ge x\big (L_1-x\big ) \end{aligned}$$

The condition \( L_1=1-M_1h-\dfrac{M_2}{A_1}>0\) provides \(h<\dfrac{A_1-M_2}{M_1 A_1}\). Therefore, the following can be concluded,

$$\begin{aligned} \liminf _{t \rightarrow \infty }x(t)\ge min\{x(0), L_1 \}\equiv k_1 \end{aligned}$$
(11)

Again, from the second equation, we observe:

$$\begin{aligned} \dfrac{dy}{dt}&=y\bigg (\delta -\frac{\beta y}{A_2+x}\bigg ) \nonumber \\&\ge y\bigg (\delta -\frac{\beta y}{A_2}\bigg ) \nonumber \\&\ge \frac{\beta y}{A_2} \bigg ( \frac{\delta A_2}{\beta }-y\bigg ) \\&\ge \frac{\beta y}{A_2} \bigg ( L_2-y\bigg ); \quad L_2=\dfrac{\delta A_2}{\beta }>0 \nonumber \\ \liminf _{t \rightarrow \infty }y(t)\ge min\{y(0), L_2 \}\equiv k_2\nonumber \end{aligned}$$
(12)

\(\square \)

Therefore, the system with initial condition \((x_0, y_0)\) is permanent for the condition \(h<\dfrac{A_1-M_2}{M_1 A_1}\). This shows that the permanence of the system depends on the harvesting rate of prey. If the harvesting rate exceeds the critical level i.e., \( h>\dfrac{A_1-M_2}{M_1 A_1}\) the system will not remain permanent.

4 Existence of equilibrium points

In this section, the existence of various equilibrium points for the system (6) are calculated which are given as follows:

4.1 Trivial equilibrium point

The existence of the trivial equilibrium point \(P_0(0, 0)\) at origin can be easily observed.

4.2 Boundary equilibrium points

The boundary equilibrium points on x-axis and on y-axis are calculated as follows:

$$\begin{aligned} P_{1}\bigg (\dfrac{1}{1+h}, 0\bigg ) \quad \text {and} \quad P_2\bigg (0, \dfrac{\delta A_2}{\beta }\bigg ) \end{aligned}$$

The points \(P_{1} (\bar{x},0)\) is known as a predator free and \(P_{2}(0, \hat{y})\) as a prey free equilibrium point, respectively.

4.3 Interior equilibrium point

The nonzero positive interior equilibrium point of this system (6) can be calculated after solving the following equations:

$$\begin{aligned}&y^*=\dfrac{\delta (A_2+x^*)}{\beta } \end{aligned}$$
(13)
$$\begin{aligned}&Ax^{*3}-Bx^{*2}+Cx^*-D=0 \end{aligned}$$
(14)

where,

$$\begin{aligned} A= & {} \beta (1+h)\nonumber \\ B= & {} \beta \nonumber \\ C= & {} \beta A_1(1+h)+\delta \nonumber \\ D= & {} (\beta A_1-\delta A_2) \end{aligned}$$
(15)

For the existence of positive interior equilibrium points, some conditions are derived in the following Lemmas 1 and 2 as follows:

Lemma 1

For \(0<\delta <\dfrac{\beta A_1}{A_2}\), system (6) can have either three or one positive real roots. This implies that we can obtain three or one positive interior equilibrium points for the condition.

Lemma 2

For \(\delta >\dfrac{\beta A_1}{A_2}\), system (6) can have two or zero positive real roots. This implies that we can obtain two or zero interior equilibrium states.

From the above lemmas, this can be observed that the system (6) can have at most three positive interior equilibrium points say \(E_1, E_2\) and \(E_3\). For different conditions on parameter \(\delta \), the numbers of equilibrium states of the system (6) may also vary. Therefore, the existence of interior equilibrium points depend on the parameter \(\delta \). However, the existence of equilibrium points for different conditions on parameter \(\delta \) are explained well in the Figs. 1 and 2. These figures depicts that how existence of various equilibrium points depend upon values of parameter \(\delta \).

Fig. 1
figure 1

For the set of parameters (37), a three interior equilibrium points for \(\delta =0.75\). b The value of \(\delta \) increases to \(\delta =0.78315914\), two equilibrium points collide and we get only two interior equilibrium points, c gives only one interior equilibrium point for \(\delta =0.8\) and in d, interior equilibrium point for \(\delta = 1.05\) does not exist

Full size image
Fig. 2
figure 2

For the set of parameters (37), a describes three interior equilibrium points for \(\delta =0.75\), b indicates that as the value of \(\delta \) decreases to \(\delta =0.63271628\), two equilibrium points collide and we get only two interior equilibrium points, c gives only one interior equilibrium point for \(\delta =0.6\) and in d, interior equilibrium point for \(\delta =0.011\) does not exist

Full size image

5 Stability analysis of equilibrium points

In the present section, the local and global asymptotic stability of various equilibrium points is explored. The stability of any equilibrium point of a nonlinear system can be related to the stability of their corresponding variational (or Jacobian) matrix. Therefore, to examine the local stability of various equilibrium points, we shall linearize the system around every equilibrium point. Further, Routh–Hurwitz criteria has been used for sufficient conditions of local stability of these points.

The variational matrix at any point (x, y) of the model (6) is obtained as follows:

$$\begin{aligned} V(x, y)= \left[ \begin{array}{cc} x\bigg (-1+\dfrac{2xy}{(A_1+x^2)^2}-h\bigg )+f &{} -\dfrac{x}{A_1+x^2} \\ \dfrac{\beta y^2}{(A_2+x)^2} &{} -\dfrac{\beta y}{(A_2+x)}+g \end{array} \right] \end{aligned}$$

The stability of various equilibrium points is investigated using above variational matrix.

Theorem 5.1

The trivial point \(P_0(0,0)\) is always an unstable node.

Proof

The variational matrix of the model (6) around (0, 0) is calculated as:

$$\begin{aligned} V(0,0)= \left[ \begin{array}{cc} 1&{} 0\\ 0 &{} \delta \end{array} \right] \end{aligned}$$

This can be seen from the above matrix that point (0, 0) is always an unstable node because the eigenvalues are \(\lambda _1=1>0\) and \(\lambda _2=\delta >0\). \(\square \)

Theorem 5.2

The boundary equilibrium point \(P_1(\bar{x},0)\) is always a saddle point.

Proof

The variational matrix of the model (6) at \( (\bar{x}, 0)\) is calculated as

$$\begin{aligned} V(\bar{x},0)= \left[ \begin{array}{cc} -\,1&{} -\,\dfrac{(1+h)}{1+A_1(1+h)^2} \\ 0 &{} \delta \end{array} \right] \end{aligned}$$

Since the eigenvalues of \( V(\bar{x}, 0)\) are \(\lambda _1=-1<0\) and \(\lambda _2=\delta >0\). This implies that the boundary equilibrium point \(P_1\) is always a saddle point. \(\square \)

Theorem 5.3

  1. 1.

    The boundary equilibrium point \( P_2(0, \hat{y})\) is locally asymptotically stable provided \(\delta >\dfrac{\beta A_1}{A_2+\beta A_1}\).

  2. 2.

    The boundary equilibrium point \( P_2(0, \hat{y})\) is always a saddle point for the condition

    $$\begin{aligned} \dfrac{\beta A_1}{A_2+\beta A_1}<\delta <\dfrac{\beta A_1}{A_2} \end{aligned}$$
  3. 3.

    The boundary equilibrium point \( P_2(0, \hat{y})\) exhibits transcritical bifurcation for \( \delta =\dfrac{\beta A_1}{A_2}\).

Proof

The variational matrix of the model (6) at \( (0, \hat{y})\) is given by

$$\begin{aligned} V(0, \hat{y})= \left[ \begin{array}{cc} 1-\dfrac{\delta A_2}{\beta A_1} &{} 0 \\ \dfrac{\delta ^2 }{\beta } &{} -\delta \end{array} \right] \end{aligned}$$

The variational matrix at boundary point \((0, \hat{y})\) gives the characteristics equation as follows:

$$\begin{aligned} \lambda ^2-\bigg (1-\dfrac{\delta A_2}{\beta A_1}-\delta \bigg )\lambda +\delta \bigg ( \dfrac{\delta A_2}{\beta A_1}-1\bigg )=0 \end{aligned}$$

The trace and determinant of above variational matrix at \((0, \hat{y})\) are calculated as follows:

$$\begin{aligned} tr V =1-\dfrac{\delta A_2}{\beta A_1}-\delta \quad and \quad det V=\delta \dfrac{\delta A_2}{\beta A_1}-1 \end{aligned}$$

The point \((0, \hat{y})\) is locally asymptotically stable for \(1-\dfrac{\delta A_2}{\beta A_1}-\delta <0\) and \(\dfrac{\delta A_2}{\beta A_1}-1>0\) which gives the following condition:

$$\begin{aligned} \delta >min \bigg [\dfrac{\beta A_1}{A_2}, \dfrac{\beta A_1}{A_2+\beta A_1}\bigg ]=\dfrac{\beta A_1}{A_2+\beta A_1} \end{aligned}$$
(16)

Further, the point \((0, \hat{y})\) is saddle for \(tr V<0\) and \(det V<0\). This gives the following condition

$$\begin{aligned} \dfrac{\beta A_1}{A_2+\beta A_1}<\delta <\dfrac{\beta A_1}{A_2} \end{aligned}$$
(17)

There is chance of occurrence of transcritical bifurcation around \((0, \hat{y})\) when \(det V=0\). This implies that the trascritical bifurcation occurs for

$$\begin{aligned} \delta =\dfrac{\beta A_1}{A_2} \end{aligned}$$
(18)

\(\square \)

Theorem 5.4

  1. 1.

    The nonzero interior equilibrium point \((x^*, y^*)\) is locally asymptotically stable provided

    $$\begin{aligned} tr V(x^*, y^*)<0 \quad and \quad det V(x^*, y^*)>0 \end{aligned}$$
  2. 2.

    The equilibrium point \((x^*, y^*)\) exhibits saddle-node bifurcation when

    $$\begin{aligned} tr V(x^*, y^*)<0 \quad and \quad det V(x^*, y^*)=0 \end{aligned}$$

Proof

The variational matrix of the model (6) at \( (x^*, y^*)\) is given by

$$\begin{aligned} V(x^*, y^*)= \left[ \begin{array}{cc} x^*\bigg (-1+\dfrac{2x^*y^*}{(A_1+{x^*}^2)^2}-h\bigg ) &{} -\dfrac{x^*}{A_1+{x^*}^2} \\ \dfrac{\beta {y^*}^2}{(A_2+x^*)^2} &{} \dfrac{-\beta y^*}{A_2+x^*} \end{array} \right] \end{aligned}$$

the characteristic equation and determinant of variational matrix \(V(x^*,y^*)\) at interior equilibrium points is given by

$$\begin{aligned} \lambda ^2- tr V(x^*, y^*)\lambda +det V(x^*, y^*)=0 \end{aligned}$$
(19)

where the trace and determinant are given as follows:

$$\begin{aligned} tr V(x^*, y^*)= & {} \dfrac{2{x^*}^2y^*}{(A_1+x^{*2})^2}-\bigg (\dfrac{\beta y^*}{A_2+x^*}+x^*+hx^*\bigg ),\\ det V (x^*, y^*)= & {} \dfrac{\beta x^* {y^*}^2}{(A_1+x^{*2})(A_2+x^*)}+\dfrac{\beta y^*}{(A_2+x^*)}\bigg ( x^*+hx^*-\dfrac{2{x^*}^2y^*}{(A_1+x^{*2})^2}\bigg ) \end{aligned}$$

Therefore, the point \((x^*, y^*)\) is locally asymptotically stable provided

$$\begin{aligned}&\dfrac{2{x^*}^2y^*}{(A_1+x^{*2})^2}-\bigg (\dfrac{\beta y^*}{A_2+x^*}+x^*+hx^*\bigg )<0 \end{aligned}$$
(20)
$$\begin{aligned}&\dfrac{\beta x^* {y^*}^2}{(A_1+x^{*2})(A_2+x^*)}+\dfrac{\beta y^*}{(A_2+x^*)}\bigg ( x^*+hx^*-\dfrac{2{x^*}^2y^*}{(A_1+x^{*2})^2}\bigg )>0 \end{aligned}$$
(21)

There is possibility of occurrence of saddle-node bifurcation provided

$$\begin{aligned}&\dfrac{2{x^*}^2y^*}{(A_1+x^{*2})^2}-\bigg (\dfrac{\beta y^*}{A_2+x^*}+x^*+hx^*\bigg )<0 \end{aligned}$$
(22)
$$\begin{aligned}&\dfrac{\beta x^* {y^*}^2}{(A_1+x^{*2})(A_2+x^*)}+\dfrac{\beta y^*}{(A_2+x^*)}\bigg ( x^*+hx^*-\dfrac{2{x^*}^2y^*}{(A_1+x^{*2})^2}\bigg )=0 \end{aligned}$$
(23)

\(\square \)

6 Local bifurcations

6.1 Transcritical bifurcation

The variational matrix at \((0,\hat{y})\) for \(\delta =\dfrac{\beta A_1}{A_2}\) gives a zero eigenvalue. Hence, the point \((0, \hat{y})\) becomes a non-hyperbolic equilibrium point. This shows the occurrence of local bifurcation around this point.

Theorem 6.1

The model (6) undergoes a transcritical bifurcation around the point \(P_2(0,\hat{y})\) for the condition

$$\begin{aligned} \delta =\delta ^{tc}=\dfrac{\beta A_1}{A_2} \end{aligned}$$

Proof

Corresponding to zero eigenvalue, the eigenvectors of \(V(0,\hat{y})\) are calculated:

$$\begin{aligned} U=\bigg (1,\dfrac{\delta }{\beta }\bigg ) \quad \text {and} \quad W=(1,0) \end{aligned}$$

Further, we will calculate \(\Delta _1,\Delta _2 \) and \(\Delta _3 \) as follows:

$$\begin{aligned}&\Delta _1=W^T F_\delta (P_2,\delta ^{tc})=0,\quad \quad where \quad \quad F=(F_1,F_2)^T=(xf,yg)^T \\&\Delta _2=W^T(D F_\delta (P_2,\delta ^{tc})U)\ne 0, \end{aligned}$$

where

$$\begin{aligned} DF_{\delta }(P_1,\delta ^{tc})= \left[ \begin{array}{cc} 0 &{} 0 \\ \dfrac{\delta ^2 }{\beta } &{} -\delta \end{array} \right] \end{aligned}$$

Also

$$\begin{aligned} \Delta _3=W^T[D F_\delta ^2 (P_2,\delta ^{tc})(U,U)] =-2-h-\frac{2\delta }{\beta }\ne 0. \end{aligned}$$

As, the value of \(\Delta _1=0\), It can be concluded that the occurrence of saddle-node bifurcation around the point \((0,\hat{y})\) is not possible. Thus, by using Sotomayor’s Theorem [33], the system (6) undergoes a transcritical bifurcation around the point \((0,\hat{y})\). \(\square \)

6.2 Saddle-node bifurcation

For \(tr V(x^*, y^*)< 0\) and \(det V(x^*, y^*)=0\), there is a possibility of saddle-node bifurcation around \((x^*, y^*)\).

Theorem 6.2

The system (6) undergoes a saddle-node bifurcation around the equilibrium point \((x^*, y^*)\) for the condition

$$\begin{aligned} \delta =\delta ^{SN} \end{aligned}$$
(24)

Proof

Corresponding to zero eigenvalue, the eigenvectors of \(V(x^*, y^*)\) and \((V(x^*, y^*))^T\) are obtained as

$$\begin{aligned} U= & {} \bigg (1,(A_1+{x^*}^2)\bigg (-1+\frac{2x^*y^*}{(A_1+(x^*)^2)^2}-h\bigg )\bigg )^T\\= & {} (u_1, u_2)^T \quad and \quad W=(0,1)^T,\quad respectively. \end{aligned}$$

Compute \( \Delta _1\) is as follows:

$$\begin{aligned} \Delta _1 = W^T F_{\delta } (P_3,\delta ^{tc})= y^* \ne 0; \quad F=(F^1, F^2)^T \end{aligned}$$

and

$$\begin{aligned} \Delta _3 =&W^T\big [D^2F(P_3,\delta ^{tc})(U,U)\big ]= \dfrac{-2 \beta {y^*}^2 {u_1}^2}{(A_2+x^*)^3}+ \dfrac{2 \beta y^* u_1 u_2}{(A_2+x^*)^2}- \dfrac{2 \beta {u_2}^2}{A_2+x^*} \ne 0 \end{aligned}$$

Therefore, using Sotomayor’s theorem, the system undergoes a saddle-node bifurcation around \((x^*, y^*)\). \(\square \)

7 Optimal harvesting policy

In this section, we investigate the optimal harvesting policy for the system (5) which maximizes the total discounted net revenue from harvesting. Let c and \(\delta \) be the harvesting cost per unit effort and the instantaneous annual rate of discount fixed by harvesting agencies, respectively. Here E is an important parameter which helps in determining the optimal harvesting policy. Therefore, we consider E as control variable subject to the condition \(E \in [0, E_{max}]\).

The present value J of a continuous time stream of revenues given by [24] is:

$$\begin{aligned} J= \int _0^\infty \!\! e^{-\delta t}{(ph_1N^2-c)E}\,dt \end{aligned}$$
(25)

The main objective is to establish an optimal harvesting policy which provides maximum profit using (25) subjected to the state equations (5).

Pontryagin’s Maximum Principle provides the optimal level of the solution of the problem (25). The associated Hamiltonian function is given by

$$\begin{aligned} \mathcal{H} (t,N,P,E)= & {} e^{-\delta t}(ph_1N^2-c)E+\lambda _1\bigg (rN\bigg (1-\frac{N}{k}\bigg )-\dfrac{m_1 NP}{a_1+N^2}-h_1EN^2\bigg )\nonumber \\&+\lambda _2sP\bigg ( 1-\dfrac{m_2P}{a_2+N}\bigg ) \end{aligned}$$
(26)

where \(\lambda _1(t)\) and \( \lambda _2(t)\) are known as adjoint variables that corresponds to the variables N and P w.r.t. the time t, respectively. Here \(\psi (t)=e^{-\delta t}(ph_1N^2-c)-\lambda _1 h_1N^2 \) is called the switching function.

Thus, the optimal control E(t) that maximizes \(\mathcal H\) must satisfy the following conditions:

$$\begin{aligned} \overline{E}=\left\{ \begin{array}{ll} E_{max} &{} when \quad \psi (t)>0, i.e., \lambda _1 e^{\delta t}< p-\dfrac{c}{h_1N^2}, \\ 0 &{} when\quad \psi (t)<0, i.e., \lambda _1 e^{\delta t} > p-\dfrac{c}{h_1N^2}. \end{array}\right. \end{aligned}$$

Notice that Hamiltonian \(\mathcal H\) is linear in the control variable E(t) which is a combination of extreme controls and the singular control.

The condition \(E= E_{max}\) when \(\psi (t)>0,\) shows that it is beneficial to harvest the resources up to extent if the profit is positive even after paying all the expenses. The condition \(E=0\) for \(\psi (t)<0\) expresses that the term \(\lambda _1 e^{\delta t}\) exceeds the fisherman’s net economic revenue. This shows that it is no more profitable to harvest. Hence, the fisherman will not exert any profit.

The condition \(\psi (t)=0\) gives \(\mathcal H\) as independent of E(t), i.e., \(\dfrac{d \mathcal H}{dE}=0\). This provides the necessary optimal condition for \(E^*(t)\) on control set \(0\le E(t)\le E_{max}\). Accordingly, the optimal harvesting policy is:

$$\begin{aligned} \overline{E}=\left\{ \begin{array}{ll} E_{max} &{}\quad for \quad \psi (t)>0, \\ 0 &{}\quad for \quad \psi (t)<0, \\ E^*&{}\quad for \quad \psi (t)=0. \end{array}\right. \end{aligned}$$

In order to attain singular control, Pontryagin’s Maximum Principle [34] is utilized. The adjoint variables must satisfy the adjoint equations as follows:

$$\begin{aligned} \dfrac{d\lambda _1}{dt}=-\,\dfrac{\partial \mathcal H}{\partial N}, \quad \dfrac{d\lambda _2}{dt}=-\,\dfrac{\partial \mathcal H}{\partial P} \end{aligned}$$
(27)

Now, the adjoint equations are

$$\begin{aligned} \dfrac{d\lambda _1}{dt}= & {} -\,\dfrac{\partial \mathcal H}{\partial N} =2ph_1NE e^{-\delta t}+\lambda _1\bigg (r-\frac{2rN}{K}+2h_1EN-\frac{m_1P(a_1-N^2)}{(a_1+N^2)^2}\bigg ) \nonumber \\&+\,\lambda _2\bigg (\dfrac{m_2sP^2}{(a_2+N)^2}\bigg ) \end{aligned}$$
(28)
$$\begin{aligned} \dfrac{d\lambda _2}{dt}= & {} -\,\dfrac{\partial \mathcal H}{\partial P}=-\,\bigg ( h_1EN^2 e^{-\delta t}+\lambda _1\bigg (\dfrac{-m_1N}{a_1+N^2} \bigg )+\lambda _2 \bigg (s-\dfrac{2m_2sP}{a_2+N}\bigg )\bigg ) \end{aligned}$$
(29)

The singular solution on \([0, E_{max}]\) can be obtained as

$$\begin{aligned} \dfrac{\partial \mathcal H}{\partial E}=e^{-\delta t}(ph_1N^2-c)+\lambda _1h_1N^2=0 \end{aligned}$$

which gives

$$\begin{aligned} \lambda _1=\dfrac{(c-ph_1N^2)e^{-\delta t}}{h_1N^2} \end{aligned}$$
(30)

Solving the Eq. (29) using Eq. (30), a linear differential equation in \(\lambda _2\) and \((N^*, P^*)\) is obtained s.t.

$$\begin{aligned}&\dfrac{d\lambda _2}{dt}-M_1 \lambda _2= -\,e^{-\delta t}M_2; \nonumber \\&M_1= -s+\dfrac{2 m_2sP^*}{a_2+N^*}\quad \text {and}\quad M_2 = h_1{N^*}^2E-\dfrac{m_1N^*(c_1-ph_1{N^*}^2)}{(a_1+{N^*}^2)(h_1{N^*}^2)} \end{aligned}$$
(31)

Solving Eq. (31), the following is obtained,

$$\begin{aligned} \lambda _2(t)=\dfrac{M_2 e^{-\delta t}}{M_1+\delta } \end{aligned}$$
(32)

The substitution of \(\lambda _2(t)\) in (32) gives the solution of (28)

$$\begin{aligned} \dfrac{d\lambda _1}{dt}-N_1 \lambda _1= -\,e^{-\delta t}N_2 \end{aligned}$$
(33)

where

$$\begin{aligned}&N_1= -\,r+\dfrac{2rN^*}{k}-2h_1EN^*+ m_1P^*\dfrac{a_1-{N^*}^2}{(a_1+{N^*}^2)^2}\\&\quad and\quad N_2 = 2ph_1EN^*+\frac{(m_2)^2s{P^*}^2}{(\delta +m_1)(a_2+N^*)^2} \end{aligned}$$

Solving Eq. (33), we get,

$$\begin{aligned} \lambda _1(t)=-\dfrac{N_2 e^{-\delta t}}{N_1+\delta } \end{aligned}$$
(34)

Substituting the value of \(\lambda _i (i=1, 2)\) with \((N^*,P^*)\), the desired singular path is obtained as follows:

$$\begin{aligned} \dfrac{N_2h_1}{N_1+\delta }{N^*}^2=(ph_1{N^*}^2-c) \end{aligned}$$
(35)

It is observed that \(\lambda _i(t)e^{\delta t}\) \( (i=1, 2)\) is free from time t satisfying the condition at \(\infty \) and remain bounded at \(\infty \). Also,

$$\begin{aligned} \dfrac{N_2h_1}{N_1+\delta }{N^*}^2=(ph_1{N^*}^2-c)\rightarrow 0\quad as \quad \delta \rightarrow \infty \end{aligned}$$
(36)

From the above Eq. (36) we can easily observe that as the discount rate reaches to infinity then economic revenue tends to zero.

Fig. 3
figure 3

Bifurcation diagram of the system w.r.t. parameter \(\delta \) for the set of parameters (37).

Full size image
Fig. 4
figure 4

Dynamical behaviour of system around various equilibrium points for \(\delta =0.75\) for the set of parameters (37)

Full size image
Fig. 5
figure 5

Dynamical behaviour of system around various equilibrium points for \(\delta =0.78315914\) for the set of parameters (37)

Full size image
Fig. 6
figure 6

Dynamical behaviour of system around various equilibrium points for \(\delta =0.8\) for the set of parameters (37)

Full size image
Fig. 7
figure 7

Dynamical behaviour of system around various equilibrium points for \(\delta =1.05\) for the set of parameters (37)

Full size image
Fig. 8
figure 8

Dynamical behaviour of system around various equilibrium points for \(\delta =0.63271628\) for the set of parameters (37)

Full size image
Fig. 9
figure 9

Dynamical behaviour of system around various equilibrium points for \(\delta =0.6\) for the set of parameters (37)

Full size image

8 Numerical simulations

In this section, numerical simulations for the system (6) have been performed to validate the analytical results and to study the dynamical behaviour for different set of parameters. The following set of parametric values is considered:

$$\begin{aligned} A_1=0.037, h=0.0002, \beta =3.38, \delta = 0.75, A_2=0.12 \end{aligned}$$
(37)

Using above data, Figs. 1 and 2 are drawn to describe the various equilibrium points for different values of parameter \(\delta \), keeping other parameter as fixed.

Figures 1 and 2 depict that as the values of \(\delta \) increases or decreases, the number of equilibrium points also varies. In the Fig. 1, we get three different interior equilibrium points say \(E_1,E_2,E_3 \) for \(\delta =0.75\) and it is shown in the Fig. 1a. As the value of \(\delta \) is slightly increased from the value of \(\delta = 0.75\) to \(\delta = 0.78315914 \), the equilibrium points \(E_2\) and \(E_3\) collide at \(\delta = 0.78315914\) and give two interior points only (see the Fig. 1b). Therefore, when \(\delta \) exceeds 0.78315914, then one equilibrium point disappeared and only one equilibrium point has been obtained. It is clear from Fig. 1c that only one equilibrium point \(E_1\) has been obtained for \( \delta =0.8\). However, it is observed that for \( \delta =1.05\), the interior equilibrium point does not exist (see Fig. 1d).

Proceeding in the same way, as \(\delta \) is decreased from \(\delta =0.75\) to \(\delta =0.63271628\), the equilibrium points \(E_1,E_2\) collide and only two interior equilibrium points are obtained, which are shown in Fig. 2b. In Fig. 2c, only one interior equilibrium point \(E_3\) is obtained for \(\delta =0.6\). However, in the Fig. 2d, the value of \( \delta =0.011\) gives no interior equilibrium point.

The dynamical study of different equilibrium points \(E_1,E_2,E_3\) for different choices of \(\delta \) is shown in Figs. 4, 5, 6, 7, 8 and 9. From the Table 1, it is observed that \(E_2\) and \(E_3\) collide at the point \(\delta =\delta _{SN}=0.783159 \). This depicts that the saddle-node bifurcation occurs around equilibrium point \((x^*, y^*)=(0.4001, 0.13902)\) for \(\delta =\delta _{SN}=0.783159\). Its normal form coefficient is \(1.641641e+000\), which has been calculated using the software package MATCONT [35,36,37,38]. This saddle node bifurcation (LP point) is detected in the bifurcation diagram 3 using MATCONT.

From the Table 2, it is observed that as value of \(\delta \) decreases from \(\delta =0.75\), the equilibrium points \(E_1\) and \(E_2\) collide at the point \(\delta =\delta _{SN}=0.632716\). This shows that saddle-node bifurcation occurs around equilibrium point \((x^*, y^*)=(0.142616, 0.049160)\) for \(\delta =\delta _{SN}=0.632716\) and its normal form coefficient has been obtained as \(5.294928e+000\). This saddle node bifurcation (LP point) is detected in the bifurcation diagram 3 using MATCONT.

Table 1 Stability analysis of various equilibrium points for increasing values of \(\delta \) from \(\delta =0.75\)
Full size table
Table 2 Stability analysis of various equilibrium points for decreasing values of \(\delta \) from \(\delta =0.75\)
Full size table

9 Conclusion

In the present work, we have analyzed modified Leslie–Gower predator–prey dynamical system with Monod–Haldane functional response, where prey is subjected to quadratic harvesting. We observed that Monod–Haldane functional response and quadratic harvesting play an important role in the dynamics of the proposed model. From the positivity and boundedness of the system it is clear that the all solutions of the system are always positive and bounded. The system has been found permanent for critical value of harvesting rate of prey. It is observed that the system will not be remain permanent if harvesting exceeds the critical value of harvesting rate. The existence of all equilibrium states give sufficient conditions on parameter \(\delta \). It is found that the system exhibits maximum three nontrivial interior equilibrium states for maximum value of mortality rate of predator. The stability analysis of all steady states have been investigated under some suitable conditions. Routh–Hurwitz criterion is used to investigate the dynamics of each equilibrium points. It is also noticed that the system undergoes local bifurcations like transcritical bifurcation around boundary equilibrium point\((0,\hat{y})\) and saddle-node bifurcation around interior equilibrium point which is investigated using Sotomayor’s theorem. The optimal harvesting policy has been carried out using Pontryagin’s maximum principal to show that at what extent harvesting has to be done. The numerical results are presented in the different figures and tables, which reflect the dynamical behavior of the system effectively.