---

Journal of Applied Nonlinear Dynamics 12(3) (2023) 465--483 | DOI:10.5890/JAND.2023.09.004

Md Sarijul Islam, Sahabuddin Sarwardi

Department of Mathematics and Statistics, Aliah University, IIA/27, New Town, Kolkata - 700 160, India

Download Full Text PDF

 

Abstract

In this article, we have studied a mathematical model for one prey and two predator with disease in predator. Here the total population is divided into three classes, namely, prey, sound predator and infected predator. We consider two different nonlinear prey refuge coefficients against the sound predator and infected predator respectively. Also both the sound predator and infected predator are harvested in the proposed dynamical system. The possibility of existence of bionomic equilibrium has been considered. The optimal harvesting policy is studied by using Pontryagin's maximal principle. The optimal harvesting efforts corresponding to the optimal solution have been derived. The positive invariance and boundedness of the solutions of the system are shown. The existence of feasible equilibria and their stability analyses are performed. It is found that the zero equilibrium point is unstable, while other equilibrium points are asymptotically stable under certain conditions. We have investigated that harvesting parameters bear an important role to control the spread of infection. Moreover, the increment of predation rates of sound and infected predator change the stability of prey only equilibrium point to infected predator-free equilibrium point and sound predator-free equilibrium point respectively. Suitable graphical representation with proper discussions are performed in the Numerical simulation section to support the proposed dynamical system.} [\hfill Eco-epidemiological model\par \hfill Nonlinear prey refuge\par \hfill Stability\par \hfill Hopf bifurcation\par \hfill Optimal harvesting policy\par \hfill Numerical simulations][Jorge Duarte][12 December 2021][20 February 2022][1 July 2023][2023 L\&H Scientific Publishing, LLC. All rights reserved.] \maketitle %\thispagestyle{fancy} \thispagestyle{firstpage} \renewcommand{\baselinestretch}{1} \normalsize \section{Introduction} Mathematical biology has become one of the most important subject to the researchers due to the variety of biological process in both the ecology and mathematics. The dynamical relationship between a predator and a prey is one of the preponderant themes in both ecology and mathematical ecology. Mathematical models bear a major role in the fields of study. Mathematical modeling is a systematic methodology by which one can construct model and understand the biological phenomena of the model system. Alfred James Lotka \cite{Lotka} and Vito Volterra \cite{Volterra} have done fundamental work on formulating and analyzing the first predators prey model. After them, Holling \cite{Holling1959,Holling1965} has introduced some more realistic predator prey model considering some types of functional response for different species with the phenomena of predation. Predator functional response \cite{Ghosh2021, Molla2018} describes the number of prey consumed by one predator per unit time for given quantities of the prey and predator. There are many important biological factors like predation, migration, refuge, harvesting, cannibalism, fear effect among prey species that brace a major role to change the dynamical behaviour of a model system \cite{Hassell, McNair1986, Mukherjee2016, Haque2014, Sk2021, Barman2021, Huang2020, Barman@, Roy2020}. Prey refuge is a biological traits by which a prey species survives itself from predation and reduces the chance of extinction due to predation. It has been shown that refuges have both stabilizing \cite{Hassell} and destabilizing effect \cite{McNair1986}. The conventional ways in which the effect of refuge are used by the preys have been included in predator prey models is to consider either a constant number or a constant proportion of the prey population. Many authors have pointed out, the refuges, which protect a constant proportion of prey have a stronger stabilizing effect on the population of dynamics than the refuges, which protect a constant number of preys. For more references, one could refer to \cite{Molla2020,Manarul,Gonzalez2012}. Another important factor is harvesting. In natural ecosystem, many of the species are being harvested. Harvest of populations are widely practiced in fishery and forestry. Many authors have studied \cite{Thota2020, Pal2014, Harekrishna, Bhattacharya, Pal2013, Liu1, Jana2015} the effects of harvesting on the eco-epidemic predator prey model. Some of authors have considered harvesting in their model without disease in the prey or predator population and the others of them have considered harvesting with disease in prey or predator population. In the present article, we have considered disease in predator population with predator harvesting. The bionomic equilibrium and the optimal harvesting policy \cite{Pal2014, Kar9, Kar6} have been also studied in the present investigation. Eco-epidemiology has become one of the most important research subject in it's own right. It is a major branch in mathematical biology, which commonly regards to both the ecological and epidemiological problems. Firstly, Anderson and May \cite{Anderson1981} merged the above two fields and formulated predator prey model where prey species were infected by some disease. In the following time, many researchers have studied different predator-prey models in the presence of disease. During last two decades, many sophisticated work have been done on this field \cite{Mainul1, Sarwardi1, Rahman1, Sieber, Melese, Sarwardi2010, Peng2018, Bezabih2020, Mainul10, Huang2019}. Spread of disease due contact of infected prey or infected predator can be controlled (cf. \cite{Amalia2018, Sowndarrajan2020}). An important issue in such types of models is the choice of the functional response. Absos et al \cite{Absos2017} considered Holling type III response function, Haque et al \cite{Mainul1} considered ratio-dependent functional response in their study. Pusawidjayanti et al \cite{Pusawidjayanti2015} considered Holling type I response function incorporating constant prey refuge with predator harvesting in their model. Our model system is different from them in the sense that we have considered predator harvesting with Holling type I response function incorporating prey refuge which is constant proportion of the prey. The outline of our article is organized as follows: In Section $2$, we have formed an eco-epidemiological model with some basic assumptions. Section $3$ deals with some basic properties of the model with positivity and boundedness of solutions. The existence and feasibility of the equilibria are derived in Section $4$. The system behaviors: the local stability at feasible equilibria, global stability around infected predator free equilibrium point and positive coexistence equilibrium point, different types of bifurcation at the points of equilibria are discussed in Section $5$. The bionomic equilibrium is studied in Section $6$ and the optimal harvesting policy is discussed in Section $7$. To justify the theoretical findings, we have performed numerical computation in Section $8$. Besides, Discussion and conclusion part is provided in Section $9$. % \section{Model formulation} To describe the system mathematically, we have divided the populations into three categories, namely, prey population, sound predator and infected predator. Let us assume at any time $t$, the biomass of prey species be $x(t)$ and that of sound predator and infected predator are $y(t)$ and $z(t)$ respectively. In the absence of predator, the prey population grows logistically with intrinsic growth rate $r$ and environmental carrying capacity $K$. This function is represented by $f(x)=rx(1-\frac{x}{K})$. We have considered that both the sound and infected predator consume prey population with Holling type I functional response at the rate $m$ and $n$ respectively. Holling type I functional response is also known as mass action type functional response. It has the form $px$, where $p$ is the predation rate and $x$ denotes the prey population. In this type of predation, it is considered that the predator do not need any handling time or saturation time to catch the prey. The conversion rate of the biomass from the prey population to the sound and infected predator are taken as $e_1$ and $e_2$ respectively. The nonlinear prey refuge coefficients against the sound and infected predators are $\delta_1$ and $\delta_2$ respectively. The sound predator becomes infected due to their direct contact with the infected predator and let us suppose $\sigma$ as the coefficient of disease transmission factor. Both the predators are harvested. We have taken into account that $E_1$ and $E_2$ are efforts of harvesting of sound predator $y$ and infected predator $z$ respectively. Further, we take $q_1$ and $q_2$ as the catchability coefficients for $y(t)$ and $z(t)$, respectively. Based on the above assumptions, we have the following mathematical model: \begin{equation}\label{Paper1} \left\{ \begin{aligned}\frac{dx}{dt} &= rx(1-\frac{x}{K})-mxy-nxz,\\ \frac{dy}{dt} &=e_1m(1-\delta_1 y)xy-\sigma yz- E_1 q_1 y,\\ \frac{dz}{dt} &=e_2n(1-\delta_2 z)xz+\sigma yz-E_2 q_2 z,\end{aligned}\right. \end{equation} subject to the initial conditions $x(0)\geq 0, ~y(0)\geq 0, ~z(0)\geq 0$. The system parameters $r$, $K$, $m$, $n$, $e_1$, $e_2$, $\sigma$, $\delta_1$, $\delta_2$, $E_1$, $E_2$, $q_1$, $q_2$ associated with the model system are all positive. The biological meaning and unit of these parameters are shown in the Table $2$. \section{Mathematical study of the system} In this section, some basic properties such as positive invariance and boundedness of the solutions of the system (\ref{Paper1}) are illustrated. Using Birkhoff and Rota \cite{Birkhoff1978}, we have also shown that the system is uniformly bounded under some conditions. \subsection{Positive invariance} Biologically, positivity implies the survival of the populations i.e., ecologically well behaved in nature. The system (\ref{Paper1}) can be written as: \begin{eqnarray} \frac {dX}{dt} = f(X,t)=\begin{pmatrix}f_1\\f_2\\f_3\end{pmatrix}=\begin{pmatrix} rx\left(1-\frac{x}{K}\right)-mxy-n xz\\ e_1m(1-\delta_1y) xy-\sigma yz -E_1 q_1 y \\e_2n(1-\delta_2 z)xz+\sigma yz-E_2 q_2 z\end{pmatrix}\end{eqnarray} where, $ X(t)=(x(t), y(t), z(t))^T$. \begin{theorem} The solution with positive initial conditions exists and unique in $[0, \alpha),$ where $0<\alpha<\infty$ (cf. Hale \cite{Hale1989}, Gard and Hallam \cite{GaHa}).\end{theorem} \begin{proof} We consider the system of equation as: \begin{eqnarray}&& \dot x=rx(1-\frac{x}{K})-mxy-n xz=x R_1(x,y,z),\nonumber\\ &&\dot y=e_1m(1-\delta_1y) xy-\sigma yz -E_1 q_1 y=y R_2(x,y,z),\nonumber \\ && \dot z= e_2n(1-\delta_2 z)xz+\sigma yz-E_2 q_2 z=z R_3(x,y,z),\nonumber \end{eqnarray} where $R_i$'s $(i=1,2,3)$ are smooth functions. Therefore, $ X(t) = X(0)e^{\int_{0}^{t} R(x,y,z)dt}\geq 0,$ where $R \equiv (R_1,R_2,R_3)^T$. Thus, the solution with positive initial conditions exists and unique in $[0, \alpha),$ where $0<\alpha<\infty$. \end{proof} \subsection{Boundedness of the system} \begin{theorem} The prey population of the system (\ref{Paper1}) is always bounded from the above. \end{theorem} \begin{proof} From the first equation of the system (\ref{Paper1}), we have \begin{eqnarray}\frac{dx}{dt}\leq {rx(1-\frac{x}{K})}.\nonumber\end{eqnarray} Therefore, by using simple and standard theorem of differential inequality, we have \begin{eqnarray}\lim_{t\to \infty}\sup\,\,x\leq K.\nonumber\end{eqnarray} Hence, the prey population of the system (\ref{Paper1}) is always bounded from above. \end{proof} \begin{theorem} For $e_1, ~e_2>0$ and $\phi<\min\{E_1q_1,E_2q_2\}$, all the solutions of the system starting in $\mathbb{R}^3_+$ are uniformly bounded. \end{theorem} \begin{proof} Considering a function $W=(e_1+e_2)x+y+z$, we have \begin{align} \frac{dW}{dt}+\phi W &=(e_1+e_2)\frac{dx}{dt}+\frac{dy}{dt}+\frac{dz}{dt}+\phi(e_1x+e_2x+y+z)\nonumber\\ &\leq[r(e_1+e_2)(1-\frac{x}{K})+(\phi e_1+\phi e_2)]x+(\phi-E_1q_1)y+(\phi-E_2q_2)z\nonumber\\ &\leq{(e_1+e_2)[rx(1-\frac{x}{K})+\phi x]},\,\,\text{if}\,\, \phi<\min\{E_1q_1,E_2q_2\}\nonumber\\ &\leq\frac{K}{4r}(e_1+e_2)(r+\phi)^2=\rho (say).\nonumber \end{align} So, $\exists$ a positive number $\rho$ can be found such that $\frac{dW}{dt}+\phi W\leq \rho$. Applying the theory of differential inequality \cite{Birkhoff1978}, we have $00\}$. \end{proof} \section{Points of equilibria} The system of equation has the following equilibria: $(i)$ Trivial equilibrium: $P_1(0,0,0)$. $(ii)$ The prey-only equilibrium: $P_2(K,0,0)$ and this equilibrium point is always feasible. $(iii)$ The prey-free equilibrium: $P_3(0,\frac{E_2q_2}{\sigma},\frac{-E_2q_2}{\sigma})$ and this equilibrium point is always infeasible. So it is not worthy to analyze further. $(iv)$ The sound predator-free equilibrium: $P_4(x_4,0,z_4)$ where $z_4=\frac{r(K-x_4)}{Kn}$ and $x_4$ is a positive root of the equation: \begin{eqnarray} \delta_2 e_2rx^2+(Ke_2n-\delta_2 e_2rK)x-KE_2q_2=0. \end{eqnarray} This equilibrium point is feasible if $K>x_4$. $(v)$ The infected predator-free equilibrium : $P_5(x_5,y_5,0)$ where $y_5=\frac{r(K-x_4)}{Km}$ and $x_4$ is a positive root of the equation: \begin{eqnarray} \delta_1e_1rx^2+(Ke_1m-K\delta_1e_1r)x-E_1q_1K=0 \end{eqnarray} This equilibrium point is feasible if $K>x_5$. $(vi)$ The coexistence equilibrium point $P_6(x_6,y_6,z_6)$, where $y_6=\frac{e_1e_2mn\delta_2x_6^2+\sigma e_2nx_6-\sigma E_2q_2-e_2n\delta_2x_6E_1q_1}{\sigma^2+e_1e_2\delta_1\delta_2mnx_6^2}$ and $z_6=\frac{e_1m\sigma x_6+e_1m\delta_1 x_6E_2q_2-e_1m\delta_1 x_6^2e_2n-E_1q_1\sigma}{\sigma^2+e_1e_2\delta_1\delta_2mnx_6^2}$ and $x_6$ is a positive root of the following cubic equation: \begin{eqnarray} \sum_{i=0}^{i=3}B_iz^i=0,\end{eqnarray} where the coefficients are given as follows: \begin{align} B_0&=E_2q_2mK\sigma-n\sigma E_1q_1K-\sigma^2rK,\nonumber\\ B_1&=-E_2q_2mKe_1\delta_1n+n\sigma e_1mK+\sigma^2r-ne_2\delta_2mE_1q_1K-e_2nmK\sigma,\nonumber\\ B_2&=e_2n^2mKe_1\delta_1+ne_2\delta_2m^2e_1K-ne_2\delta_2me_1\delta_1rK,\nonumber \\ B_3&=ne_2\delta_2me_1\delta_1 r.\nonumber \end{align} This equilibrium point is feasible if $\sigma>\max\{\frac{e_2n\delta_2x_6(E_1q_1-x_6e_1m)}{e_2nx_6-E_2q_2},\frac{e_1m\delta_1x_6(E_2q_2-x_6e_2n)}{e_1mx_6-E_1q_1}\}$. \begin{proposition} If the basic reproduction number $R_0=\frac{K(e_2n+\sigma)}{E_2q_2}<1$ then $\frac{dz}{dt}<0$, that means the infection will not spread. This characteristic is related to a similar threshold phenomenon in epidemic theory \cite{Anderson81, Diekmann}. This $R_0$ is called basic reproductive ratio. In epidemic theory, if the value of $R_0>1$ then the epidemics outbreak will occur. This number $R_0$ can be thought of as the number of all individuals who will get the disease in all time following `successful' contact with a sick individual, that is, it is considered as the expected number of secondary cases produced by an infective during its whole infectious period. \end{proposition} \section{Stability analysis} The Jacobian matrix $J$ of the system (\ref{Paper1}) at any arbitrary point $(\bar{x},\bar{y},\bar{z})$ is given by $J(\bar{x},\bar{y},\bar{z})=(a_{ij})_{3\times3},$ where $a_{11}=r-\frac{2rx}{K}-my-nz, ~a_{12}=-mx, a_{13}=-nx, ~a_{21}=e_1my(1-\delta_1 y), ~a_{22}=e_1mx-2e_1m\delta_1xy-\sigma z-E_1q_1, ~a_{23}=-\sigma y, ~a_{31}=e_2nz(1-\delta_2 z), ~a_{32}=\sigma z, ~a_{33}=e_2nx-2e_2n\delta_2 xz+\sigma y-E_2q_2$. We denote $J_s=J$, the Jacobian evaluated at $P_s$, where $s=1,2,3,4,5,6$. The determinant of $J_s=\det(J)$, the trace of $J_s=tr(J)$ and $C_s=tr(J_s)M_s-\det(J_s)$, where $M_s$ denotes the sum of the second order principal minors of $J_s$. \subsection{Local stability analysis of the system} In the following local stability at different equilibrium point are discussed $(i)$ At the trivial equilibrium point $P_1(0,0,0)$, the eigenvalues of $J_1$ are $r, ~-E_1q_1, ~-E_2q_2$. Therefore, it is an unstable hyperbolic saddle critical point in the direction orthogonal to $yz$ plane. $(ii)$ At the axial equilibrium point $P_2(K,0,0)$, the eigenvalues of $J_2$ are $-r, ~e_1mK-E_1q_1$ and $e_2nK-E_2q_2$. If $K<\min\{\frac{E_1q_1}{e_1m}, ~\frac{E_2q_2}{e_2n}\}$, ~$P_2$ is locally asymptotically stable. $(iii)$ Since prey-free equilibrium point $P_3$ is not feasible, we skip this. $(iv)$ For the sound predator-free equilibrium point $P_4(x_4,0,z_4)$, $a_{11}=r-\frac{2rx_4}{K}, ~a_{12}=-mx_4, ~a_{13}=-nx_4, ~a_{21}=0, ~a_{22}=e_1mx_4-\sigma z_4-E_1q_1, ~a_{23}=0, ~a_{31}=e_2nz_4(1-\delta_2 z_4), ~a_{32}=\sigma z_4, ~a_{33}=e_2nx_4-2e_2\delta_2nx_4z_4-E_2q_2$. It is clear that when $K<\frac{2rx_4}{r-nz_4}, ~x_4<\min\{\frac{\sigma z_4+E_1q_1}{e_1m}, ~\frac{E_2q_2}{e_2n(1-2\delta_2 x_4)}\}$ and $\delta_2<\frac{1}{z_4}$, we get $a_{11}<0, ~a_{22}<0, ~a_{33}<0, ~a_{32}>0, ~a_{31}>0$. Then obviously $tr(J_4)<0, ~\det(J_4)=a_{22}(a_{11}a_{33}-a_{31}a_{13})<0, ~M_4=a_{11}(a_{22}+a_{33})+a_{22}a_{33}-a_{13}a_{31}>0$ and $C_4=tr(J_4)M_4-\det(J_4)=a_{11}^2(a_{22}+a_{33})+a_{22}^2(a_{11}+a_{33})+a_{33}^2(a_{11}+a_{22})+2a_{11}a_{22}a_{33}-a_{11}a_{13}a_{31}-a_{13}a_{31}a_{33}<0$. Hence the Routh-Hurwitz condition is satisfied. Therefore, $P_4$ is locally asymptotically stable for $K<\frac{2rx_4}{r-nz_4}, ~x_4<\min\{\frac{\sigma z_4+E_1q_1}{e_1m}, ~\frac{E_2q_2}{e_2n(1-2\delta_2 x_4)}\}$ and $\delta_2<\frac{1}{z_4}$. $(v)$ For the infected predator-free equilibrium point $P_5(x_5,y_5,0)$, $a_{11}=r-\frac{2rx_5}{K}-my_5, ~a_{12}=-mx_5, ~a_{13}=-nx_5, ~a_{21}=e_1my_5(1-\delta_1 y_5), ~a_{22}=e_1mx_5-2e_1m\delta_1x_5y_5-E_1q_1, ~a_{23}=-\sigma y_5, ~a_{31}=0, ~a_{32}=0, ~a_{33}=e_2nx_5+\sigma y_5-E_2q_2.$ It is evident that if $K<\frac{2rx_5}{r-my_5}, ~x_5<\min\{\frac{E_2q_2}{e_2n}, ~\frac{E_1q_1}{e_1m(1-2\delta_1 y_5)}\}$ and $\delta_1<\frac{1}{y_5}$, then $a_{11}<0, ~a_{21}>0, ~a_{22}<0, ~a_{33}<0$. Then $tr(J_5)<0, ~\det(J_5)=a_{33}(a_{11}a_{22}-a_{12}a_{21})<0, ~M_5=a_{33}(a_{11}+a_{22})+a_{11}a_{22}-a_{12}a_{21}>0$ and $C_5=tr(J_5)M_5-\det(J_5)=a_{11}^2(a_{22}+a_{33})+a_{22}^2(a_{11}+a_{33})+a_{33}^2(a_{11}+a_{22})+2a_{11}a_{22}a_{33}-a_{11}a_{12}a_{21}-a_{22}a_{12}a_{21}<0$. Hence by Routh-Hurwitz criteria $P_5$ is locally asymptotically stable for $K<\frac{2rx_5}{r-my_5}, ~x_5<\min\{\frac{E_2q_2}{e_2n}, ~\frac{E_1q_1}{e_1m(1-2\delta_1 y_5)}\}$ and $\delta_1<\frac{1}{y_5}$. \subsection{Local stability analysis for the positive coexistence equilibrium point $P_6(x_6,y_6,z_6)$} \begin{theorem} The system $(1)$ is locally asymptotically stable at $P_6$ if $(i)\sigma0, ~\xi_{31}>0$, also using the condition $(i)$ we get $\det(J_6)<0$. Using both the condition $(i)$ and $(ii)$ we get, $M_6>0$ and $C_6<0$. Thus the conditions of Routh-Hurwitz are satisfied for the jacobian matrix $J_6$, that means, all the characteristic roots of $J_6$ have negative real parts. Hence the claim. \end{proof} \subsection{Global stability analysis around infected predator free equilibrium point $P_5(x_5, y_5, 0)$} \begin{theorem} The system $(1)$ is globally asymptotically stable if $E_2q_2+s_2e_1m\delta_1 nx_5y_5>s_2e_1nx_5+s_2\sigma y_5+e_2nK+\sigma K.$ \end{theorem} \begin{proof} Let us consider $W_1=\{(x, y, z)\in\mathbb{R}^3_+, ~x>0, ~y>0, ~z\geq 0 \}, ~where, ~s_1>0, ~s_2>0$ are two real numbers and a scalar function $L\colon\mathbb{R}^3\to R$ defined by\\ \begin{equation} L=s_1(x-x_5-x_5\ln\frac{x}{x_5})+s_2(y-y_5-y_5\ln\frac{y}{y_5})+z. \end{equation} Then we have, \begin{align} \frac{dL}{dt}&=s_1(\frac{x-x_5}{x})\frac{dx}{dt}+s_2(\frac{y-y_5}{y})\frac{dy}{dt}+\frac{dz}{dt}\nonumber\\ &=-\frac{s_1r}{K}(x-x_5)^2-m(s_1-s_2e_1+y_5e_1\delta_1 s_2)-s_2e_1m\delta_1 x(y-y_5)^2-s_1nzx\nonumber\\&\quad +s_1nzx_5-s_2\sigma yz+s_2\sigma y_5z,\nonumber\\ &\leq-\frac{s_1r}{K}(x-x_5)^2-(s_1-s_2e_1+s_2e_1m\delta_1y_5)(x-x_5)(y-y_5)- s_2e_1m\delta_1x(y-y_5)^2\nonumber\\ &\quad +(s_1nx_5+s_2\sigma y_5+e_2nK+\sigma K-E_2q_2).\nonumber \end{align} Taking $s_1=s_2e_1(1-m\delta_1y_5)$ and $E_2q_2+s_2e_1m\delta_1nx_5y_5>s_2e_1nx_5+s_2\sigma y_5+e_2nK+\sigma K$, we obtain $\frac{dL}{dt}\leq 0.$ \end{proof} \subsection{Global stability analysis around the positive coexistence equilibrium point $P_6(x_6,y_6,z_6)$} \begin{theorem} The given model system is globally asymptotically stable under the conditions: $(i)$ $ M^{[1]}k_2$ and $(iii)$ $ \max\{\delta_1, \delta_2\}<\frac{1}{2x_6}$.\end{theorem} \begin{proof} Let us consider $W_2=\{(x, y, z)\in\mathbb{R}^3_+, x>0, y>0, z>0 \}, ~k_1>0, ~k_2>0,~k_3>0$ are three real numbers and a scalar function $L\colon\mathbb{R}^3\to R$ defined by\\ \begin{equation} L_1=k_1(x-x_6-x_6ln\frac{x}{x_6})+k_2(y-y_6-y_6ln\frac{y}{y_6})+k_3(z-z_6-z_6ln\frac{z}{z_6}). \end{equation} Then, \begin{align} \frac{dL_1}{dt}=&k_1\frac{(x-x_6)}{x}\frac{dx}{dt}+k_2\frac{(y-y_6)}{y}\frac{dy}{dt}+k_3\frac{(z-z_6)}{z}\frac{dz}{dt},\nonumber\\ =&k_1[-\frac{r}{K}(x-x_6)-m(y-y_6)-n(z-z_6)](x-x_6)\nonumber\\ &+k_2[-\sigma (z-z_6)+e_1m(x-x_6)+e_1m\delta_1(x_6y_6-xy)](y-y_6)\nonumber\\ &+k_3[\sigma (y-y_6)+e_2n(x-x_6)+e_2n\delta_2(z_6x_6-zx)](z-z_6),\nonumber\\ =&-\frac{k_1r}{k}(x-x_6)^2-k_2x_6e_1m\delta_1(y-y_6)^2-e_2n\delta_2x_6(z-z_6)^2\nonumber\\ &+[k_2e_1m(1-\delta_1 y)-k_1m](x-x_6)(y-y_6)+\sigma(k_3-k_2)(y-y_6)(z-z_6)\nonumber\\ &+[k_3e_2n(1-\delta_2 z)-k_1n](x-x_6)(z-z_6)\nonumber\\ \leq &-[(\frac{rk_1}{k}-\frac{m(k_2e_1-k_1)}{2}-\frac{n(k_3e_2-k_1)}{2})(x-x_6)^2\nonumber\\ &+(k_2e_1m\delta_1 x_6-\frac{m(k_2e_1-k_1)}{2}-\frac{\sigma(k_3-k_2)}{2})(y-y_6)^2\nonumber\\ &+(k_3e_2n\delta_2 x_6-\frac{\sigma}{2}(k_3-k_2)-\frac{(k_3e_2n-k_1n)}{2})(z-z_6)^2].\nonumber \end{align} Then by the conditions $(i), ~(ii)$ and $(iii)$, we get $\frac{dL_1}{dt}\leq 0$. Thus, $L_1$ is a Lyapunov function and by Lasalle's theorem \cite{Hale1989}, it follows that the system under the conditions $(i), ~(ii)$ and $(iii)$ possesses global stability. \end{proof} \subsection{Hopf bifurcation analysis at the coexistence equilibrium point $P_6(x_6,y_6,z_6)$} \begin{theorem} The system possesses a Hopf bifurcation around $P_6$ when $\sigma$ passes through $\sigma_h$ where $\sigma_h$ satisfy the equation $A_1A_2=A_3$ and $A_i$ satisfy the characteristic equation of the system (\ref{Paper1}) at $P_6(x_6,y_6,z_6)$, that is \begin{eqnarray} \tau^3+A_1\tau^2+A_2\tau+A_3=0 \end{eqnarray} \end{theorem} \begin{proof} The characteristic equation of our model system (\ref{Paper1}) at $P_6(x_6,y_6,z_6)$ is given by $(8)$, where $A_1=-\text{tr}(J_6)$, $A_2=M_6$, $A_3=-\det(J_6)$ and $A_1A_2-A_3=-C_6$. Hopf bifurcation will occur if and only if there exists $\sigma=\sigma_h$ such that $(i)\;A_1(\sigma)A_2(\sigma)=A_3(\sigma)$ with all $A_1(\sigma), ~A_2(\sigma), ~A_3(\sigma)>0$ and $(ii)\;(\frac{d}{d\sigma})(\Re(\tau(\sigma)))|_{\sigma=\sigma_h}\neq 0$. Now when $\sigma=\sigma_h, ~A_1A_2=A_3$ with $A_1, ~A_2, ~A_3>0$ the characteristic equation becomes \begin{equation} (\tau^2+A_2)(\tau+A_1)=0. \end{equation} The Eq. $(9)$ has the roots $\tau_1=i\sqrt{A_2}, ~\tau_2=-i\sqrt{A_2}$ and $\tau_3=-A_1$, so that there exists a strictly negative real eigenvalues and a pair of purely imaginary eigenvalues. For $\sigma$ in a neighbourhood of $\sigma_h$, let the roots be of the form $\tau_1(\sigma)=q_1(\sigma)+iq_2(\sigma), ~\tau_2(\sigma)=q_1(\sigma)-iq_2(\sigma)$ and $\sigma_3(\sigma)=-q_3(\sigma)$, where $q_1(\sigma), ~q_2(\sigma)$ and $q_3(\sigma)$ are real. Next, we shall verify the transversality condition $\frac{d}{d\sigma}(\Re(\tau_j(\sigma)))|_{\sigma=\sigma_h}\neq 0, j=1,2$. Substituting $\tau(\sigma)=q_1(\sigma)+iq_2(\sigma)$ into the characteristic equation and taking the derivative, we have \begin{eqnarray} \Omega(\sigma)\dot{q_1}(\sigma)-\phi(\sigma)\dot{q_2}(\sigma)+\Gamma(\sigma)=0,\\ \phi(\sigma)\dot{q_1}(\sigma)+\Omega(\sigma)\dot{q_2}(\sigma)+\Theta(\sigma)=0, \end{eqnarray} where we consider the following: \[ \begin{split} & \Omega(\sigma)=3(q_1(\sigma))^2+2A_1(\sigma)q_1(\sigma)+A_2(\sigma)-3(q_2(\sigma))^2, \\ & \phi(\sigma)=6q_1(\sigma)q_2(\sigma)+2A_1(\sigma)q_2(\sigma), \\ & \Gamma(\sigma)=(q_1(\sigma))^2\dot{A_1}(\sigma)+q_1(\sigma)\dot{A_2} (\sigma)+\dot{A_3}(\sigma)-\dot{A_1}(\sigma)(q_2(\sigma))^2, \\ & \Theta(\sigma)=2q_1(\sigma)q_2(\sigma)\dot{A_1}(\sigma)+\dot{A_2}(\sigma)q_2(\sigma). \end{split} \] Now, $\frac{d}{d\sigma}(\Re(\tau_j(\sigma)))|_{\sigma=\sigma_h}=-\frac{\Omega\Gamma+\phi\Theta}{\Omega^2+\phi^2}|_{\sigma=\sigma_h}\neq 0$ and $q_3(\sigma_h)=-A_1(\sigma_h)\neq 0.$ This completes the proof. \end{proof} \begin{theorem} The system (\ref{Paper1}) possesses a transcritical bifurcation around the equilibrium point $P_2(K,0,0)$ at $m=m^{[tc]}$, where $m^{[tc]}=\frac{E_1q_1}{e_1K}$. \end{theorem} \begin{proof} One of the eigenvalues of $J_2$ will be zero if and only if $\det(J_2)=0,$ which gives $m=m^{[tc]}$. The other eigenvalues of $J_2$ are $-r$ and $e_2nK-E_2q_2$ evaluated at $m^{[tc]}=\frac{E_1q_1}{e_1K}$. Next we consider $v$ and $w$ as the eigenvectors of the matrix $J_2$ and its transpose $J_2^{T}$ respectively corresponding to the eigenvalue $0$. We see that $v=\begin{pmatrix}\frac{-mK}{r}\\1\\0\end{pmatrix}$ and $w=\begin{pmatrix}0\\1\\0\end{pmatrix}$. Since $w^{T

References

1. [1]	Lotka, A.J. (1956), Elements of Mathematical Biology, Dover: New York.

2. [2]	Volterra, V. and D'Ancona U. (1931), La concorrenza vitale tra le specie dell'ambiente marino, Vlle Congr. Int. acquicult et de p\^{e}che, Paris, 1-14.

3. [3]	Holling, C.S. (1959), The components of predation as revealed by a study of small-mammal predation of the European Pine Sawfly, The Canadian Entomologist, 91(5), 293-320.

4. [4]	Holling, C.S. (1965), The functional response of predators to prey density and its role in mimicry and population regulation, The Memoirs of the Entomological Society of Canada, 97(S45), 5-60.

5. [5]	Ghosh, U., Mondal, B., Rahaman, M.S., and Sarkar, S. (2021), Stability analysis of three species food chain model with linear functional response via imprecise and parametric approach, Journal of Computational Science, 54, 101423.

6. [6]	Molla, H., Rahman, M.S., and Sarwardi, S. (2019), Dynamics of a predator-prey model with Holling type II functional response incorporating a prey refuge depending on both the species, International Journal of Nonlinear Sciences and Numerical Simulation, 20(1), 89-104.

7. [7]	Hassell, M.P. (2020), The Dynamics of Arthopod Predator-Prey Systems, (MPB-13), Volume 13, Princeton University Press.

8. [8]	McNair, J.N. (1986), The effects of refuges on predator-prey interactions: a reconsideration, Theoretical Population Biology, 29(1), 38-63.

9. [9]	Mukherjee, D. (2016), The effect of refuge and immigration in a predator-prey system in the presence of a competitor for the prey, Nonlinear Analysis: Real World Applications, 31, 277-287.

10. [10]	Haque, M., Rahman, M.S., Venturino, E., and Li, B.L. (2014), Effect of a functional response-dependent prey refuge in a predator-prey model, Ecological Complexity, 20, 248-256.

11. [11]	Sk, N., Tiwari, P.K., Kang, Y., and Pal, S. (2021), A nonautonomous model for the interactive effects of fear, refuge and additional food in a prey-predator system, Journal of Biological Systems, 29(01), 107--145.

12. [12]	Barman, D., Roy, J., and Alam, S. (2021), Dynamical behaviour of an infected predator-prey model with fear effect, Iranian Journal of Science and Technology, Transactions A: Science, 45(1), 309-325.

13. [13]	Huang, Y., Zhu, Z., and Li, Z. (2020), Modeling the Allee effect and fear effect in predator-prey system incorporating a prey refuge, Advances in Difference Equations, 2020(1), 1-13.

14. [14]	Barman, D., Roy, J., Alrabaiah, H., Panja, P., Mondal, S.P., and Alam, S. (2021), Impact of predator incited fear and prey refuge in a fractional order prey predator model, Chaos, Solitons $\&$ Fractals, 142, 110420.

15. [15]	Roy, J., Barman, D., and Alam, S. (2020), Role of fear in a predator-prey system with ratio-dependent functional response in deterministic and stochastic environment, Biosystems, 197, 104176.

16. [16]	Molla, H., Rahman, M.S., and Sarwardi, S. (2020), Incorporating prey refuge in a prey-predator model with Beddington-Deangelis type functional response: a comparative study on intra-specific competition, Discontinuity, Nonlinearity, and Complexity, 9(3), 395-419.

17. [17]	Manarul Haque, M. and Sarwardi, S. (2018), Dynamics of a harvested prey-predator model with prey refuge dependent on both species, International Journal of Bifurcation and Chaos, 28(12), 1830040.

18. [18]	González-Olivares, E. and Ramos-Jiliberto, R. (2012), Comments to ``The effect of prey refuge in a simple predator-prey model"[Ecol. Model. 222 (September (18))(2011) 3453-3454], Ecological Modelling, 232, 158-160.

19. [19]	Thota, S. (2020), A mathematical study on a diseased prey-predator model with predator harvesting, Asian Journal of Fuzzy and Applied Mathematics, (ISSN: 2321-564X), 8(2).

20. [20]	Pal, D. and Mahapatra, G.S. (2014), A bioeconomic modeling of two-prey and one-predator fishery model with optimal harvesting policy through hybridization approach, Applied Mathematics and Computation, 242, 748-763.

21. [21]	Das, H., Shaikh, A.A., and Sarwardi, S. (2020), Mathematical analysis of an eco-epidemic model with different functional responses of healthy and infected predators on prey species, Journal of Applied Nonlinear Dynamics, 9(4), 667-684.

22. [22]	Bhattacharyya, R. and Mukhopadhyay, B. (2010), On an eco-epidemiological model with prey harvesting and predator switching: local and global perspectives, Nonlinear Analysis: Real World Applications, 11(5), 3824-3833.

23. [23]	Pal, D., Mahaptra, G.S., and Samanta, G.P. (2013), Optimal harvesting of prey-predator system with interval biological parameters: a bioeconomic model, Mathematical biosciences, 241(2), 181-187.

24. [24]	Liu, X. and Huang, Q. (2020), Analysis of optimal harvesting of a predator-prey model with Holling type IV functional response, Ecological Complexity, 42, 100816.

25. [25]	Jana, S., Guria, S., Das, U., Kar, T.K., and Ghorai, A. (2015), Effect of harvesting and infection on predator in a prey-predator system, Nonlinear Dynamics, 81(1), 917-930.

26. [26]	Kar, T.K., Chattopadhyay, S.K., and Pati, C.K. (2009), A bio-economic model of two-prey one-predator system, Journal of applied mathematics $\&$ informatics, 27(5/6), 1411-1427.

27. [27]	Kar, T.K., Misra, S., and Mukhopadhyay, B. (2006), A bioeconomic model of a ratio-dependent predator-prey system and optimal harvesting, Journal of Applied Mathematics and Computing, 22(1), 387-401.

28. [28]	Anderson, R.M. and May, R.M. (1981), The population dynamics of microparasites and their invertebrate hosts, Philosophical Transactions of the Royal Society of London. B, Biological Sciences, 291(1054), 451-524.

29. [29]	Haque, M. and Venturino, E. (2007), An ecoepidemiological model with disease in predator: the ratio dependent case, Mathematical methods in the Applied Sciences, 30(14), 1791-1809.

30. [30]	Sarwardi, S., Haque, M., and Venturino, E. (2011), A Leslie-Gower Holling-type II ecoepidemic model, Journal of Applied Mathematics and Computing, 35(1), 263-280.

31. [31]	Rahman, M.S. and Chakravarty, S. (2013), A predator-prey model with disease in prey, Nonlinear Analysis: Modelling and Control, 18(2), 191-209.

32. [32]	Sieber, M., Malchow, H., and Hilker, F.M. (2014), Disease-induced modification of prey competition in eco-epidemiological models, Ecological Complexity, 18, 74-82.

33. [33]	Melese, D., Muhye, O., and Sahu, S.K. (2020), Dynamical Behavior of an Eco-epidemiological Model Incorporating Prey Refuge and Prey Harvesting, Applications $\&$ Applied Mathematics, 15(2).

34. [34]	Sarwardi, S., Haque, M., and Venturino, E. (2011), Global stability and persistence in LG-Holling type II diseased predator ecosystems, Journal of Biological Physics, 37(1), 91-106.

35. [35]	Peng, M., Zhang, Z., Lim, C.W., and Wang, X. (2018), Hopf bifurcation and hybrid control of a delayed ecoepidemiological model with nonlinear incidence rate and Holling type II functional response, Mathematical Problems in Engineering, 2018.

36. [36]	Bezabih, A.F., Edessa, G.K., and Koya, P.R. (2020), Mathematical eco-epidemiological model on prey-predator system, Mathematical Modeling and Applications, 5(3), 183-190.

37. [37]	Haque, M., Rahman, M.S., and Venturino, E. (2013), Comparing functional responses in predator-infected eco-epidemics models, BioSystems, 114(2), 98-117.

38. [38]	Huang, C., Zhang, H., Cao, J., and Hu, H. (2019), Stability and Hopf bifurcation of a delayed prey-predator model with disease in the predator, International Journal of Bifurcation and Chaos, 29(07), 1950091.

39. [39]	Amalia, R.D. and Arif, D.K. (2018), Optimal control of predator-prey mathematical model with infection and harvesting on prey, In Journal of Physics: Conference Series, 974(1), p. 012050, IOP Publishing.

40. [40]	Sowndarrajan, P.T., Nyamoradi, N., Shangerganesh, L., and Manimaran, J. (2020), Mathematical analysis of an optimal control problem for the predator-prey model with disease in prey, Optimal Control Applications and Methods, 41(5), 1495-1509.

41. [41]	Shaikh, A.A., Das, H., and Ali, N. (2018), Study of LG-Holling type III predator-prey model with disease in predator, Journal of Applied Mathematics and Computing, 58(1), 235-255.

42. [42]	Pusawidjayanti, K., Suryanto, A., and Wibowo, R.B.E. (2015), Dynamics of a predator-prey model incorporating prey refuge, predator infection and harvesting, Applied Mathematical Sciences, 9(76), 3751-3760.

43. [43]	Birkhoff, G. and Rota, G. (1978), Ordinary differential equations, John Wiley and Sons, New York: 342 QA372 B, 58.

44. [44]	Hale, J. (1989), Ordinary Differential Equation, Klieger Publishing Company, Malabar.

45. [45]	Gard, T.C. and Hallam T.G. (1979), Persistence in Food web-1, Lotka-Volterra food chains, Bulletin of Mathematical Biology, 41(6), 877-891

46. [46]	Anderson, R.M. and May, R.M. (1981), The population dynamics of microparasites and their invertebrate hosts, Philosophical Transactions of the Royal Society of London. B, Biological Sciences, 291(1054), 451-524.

47. [47]	Diekmann, O., Heesterbeek, J.A.P., and Metz, J.A. (1990), On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious diseases in heterogeneous populations, Journal of Mathematical Biology, 28(4), 365-382.

48. [48]	Sotomayor, J. (1973), Generic bifurcations of dynamical systems: Dynamical systems, Academic Press, 549-560.

49. [49]	Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V., and Mishchenko, E.F. (1962), The Mathematical Theory of Optimal Process, Wiley, New York.