\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 209, pp. 1--27.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/209\hfil Dynamics of a prey-predator system] {Dynamics of a prey-predator system with infection in prey} \author[S, Kant, V. Kumar \hfil EJDE-2017/209\hfilneg] {Shashi Kant, Vivek Kumar} \address{Shashi Kant \newline Department of Applied Mathematics, Delhi Technological University, Delhi 110042, India} \email{onlineskmishra@gmail.com} \address{Vivek Kumar \newline Department of Applied Mathematics, Delhi Technological University, Delhi 110042, India} \email{vivek\_ag@iitkalumni.org} \dedicatory{Communicated by Zhaosheng Feng} \thanks{Submitted May 23, 2016. Published September 8, 2017.} \subjclass[2010]{34C23, 34C25, 34C28} \keywords{Predator-prey model; linear Functional Response; \hfill\break\indent Hopf-bifurcation; stability analysis; time delay} \begin{abstract} This article concerns a prey-predator model with linear functional response. The mathematical model has a system of three nonlinear coupled ordinary differential equations to describe the interaction among the healthy prey, infected prey and predator populations. Model is analyzed in terms of stability. By considering the delay as a bifurcation parameter, the stability of the interior equilibrium point and occurrence of Hopf-bifurcation is studied. By using normal form method, Riesz representation theorem and center manifold theorem, direction of Hopf bifurcation and stability of bifurcated periodic solutions are also obtained. As the real parameters are not available (because it is not a case study). To validate the theoretical formulation, a numerical example is also considered and few simulations are also given. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} The study of prey-predator systems has been a burning topic of research for several years. The pioneer work of Kermack and Mckendrick on Susceptible-Infective-Recovered-Susceptible (SIRS) models \cite{WS} have been widely accepted among researchers and scientific community. After the work of Kermack and Mckendrick \cite{WS} many mathematical models have been published (\cite{ MH, ST, RM, SJB,YX} etc. and references therein). M. Haque et al \cite{MH} proposed and analyzed a predator prey model using standard disease incidence. They observed that the disease in the prey may avoid extinction of predators and its presence can destabilize an otherwise stable configuration of species. In \cite{RM}, Naji and Mustafa investigated the dynamical nature of an eco-epidemiological model by applying nonlinear disease incidence rate among living species of the ecosystem. They proposed and investigated with regards to local and global dynamical nature of Holling type-II model with Susceptible-Infective type of disease in prey \cite{RM}. Jang and Baglama \cite{SJB} proposed a deterministic continues time ecological model with the effect of parasites, where it is assumed that intermediate host for the parasites are the prey species and observed the dynamics of it. They conclude that parasites are in position to affect the dynamics of the predator prey interaction due to infection. Jang and Baglama \cite{SJB} have also proposed a stochastic version of the model and simulated the model numerically to verify the theoretical results. They performed asymptotic dynamics and compared the deterministic and stochastic models \cite{SJB}. Jana and Kar \cite{ST} proposed and analyzed a three dimensional eco-epidemiological model consisting of susceptible prey, infected prey and predator. They introduced time delay in the model for considering the time delay as the time taken by a susceptible prey to become infected. Mathematically, they analyzed the dynamics of the model in terms of existence of non-negative equilibria, boundedness, local and global stability of the interior equilibrium point. They also studied Hopf bifurcation and by using central manifold reduction they investigated the direction of Hopf bifurcation and stability of limit cycles. Many mathematical models have been proposed to understand the evolution of diseases and provided valuable information for control strategies (\cite{KB, LZ, ML, FK} and references therein). Hilker and Schmitz \cite{FK} proved that predator infection counteracts the paradox of enrichment. They discussed the implication for the biological control and resource management on more than one trophic level. Ecology and epidemiology are two different major and important research areas. The basic work of Lotka \cite{AL} and Volterra \cite{VV} on predator-prey models in the form of coupled system of non-linear differential equations may be considered as the first break through in the modern mathematical ecology. Further, overlapping study of ecology and epidemiology termed as eco-epidemiology. In eco-epidemiology, we study prey-predator models with disease dynamics. Thus, eco-epidemiology may be considered as the study of interacting species in which disease spreads. Eco-epidemiology has very important ecological significance. Population growth models with disease spreading often provide complex non-linear mathematical dynamics. In these models the main concern is to study equilibrium points, their stability analysis, periodic solutions, bifurcations, chaotic nature etc. A large number of mathematical and statistical techniques are available to analyze the eco-epidemiological models. While formulating a prey-predator model, it is a basic assumption that reproduction of predator species after the event of predation will not be instantaneous, but it will be mediated by some discrete time lag (delay) essential for the gestation of predator population \cite{GX}. To study mathematical models in ecology more scientifically, peoples coined a new word `time delay'. Time delay has been used in large number of papers e.g. (\cite{BB,FW} are few name to). Mukhopadhyaya and Bhattacharyya \cite{BB} studied the effect of delay on a prey predator model with disease in prey. They have considered Holling type II functional response. Fengyan Wang et al \cite{FW} studied a predator prey model by assuming stages viz. mature and immature with both discrete and distributed delays. They considered delay as length of immature stage. For detailed study of delay differential equations we can refer reader to \cite{YK}. Chattopadhyay and Arino \cite{JO} proposed the following eco-epidemiological model with disease in prey \begin{equation}\label{1.1} \begin{gathered} \frac{dS}{dt}=r(S+I)(1-\frac{S+I}{K})-\beta SI-\eta \gamma_{1}(S)Y,\\ \frac{dI}{dt}=\beta SI-\gamma(I)Y-CI,\\ \frac{dY}{dt}=(\varepsilon\gamma(I)+\eta\varepsilon\gamma_{1}(S)-d)Y, \end{gathered} \end{equation} where, $S$ is the number of sound prey, $I$ is the number of infected prey population, $Y$ is the number of predator population, $\gamma(I)$ and $\eta\gamma_{1}(S)$ are predator functional response functions. They analyzed the model \eqref{1.1} in terms of positivity, uniqueness, boundedness and the study the existence of the Hopf bifurcation. Model \eqref{1.1} may be re-written in simplified form as \begin{equation}\label{1.2} \begin{gathered} \frac{dS}{dt}=rS(1-\frac{S+I}{K})-\beta SI,\\ \frac{dI}{dt}=-cI+\beta SI -p I Y,\\ \frac{dY}{dt}=-dY+ p q IY. \end{gathered} \end{equation} Motivated by model \eqref{1.1}, Samanta\cite{S} proposed a diseased nonautonomous predator-prey system with time delay, which is given as \begin{equation}\label{1.3} \begin{gathered} \frac{dx_{1}(t)}{dt}=x_{1}(t)[r(t)-k_{1}(t)(x_{1}(t)+x_{2}(t))-a_{1}(t)x_{3}(t) -\beta(t)x_{2}(t)],\\ \frac{dx_{2}(t)}{dt}=x_{2}(t)[r(t)-k_{2}(t)(x_{1}(t)+x_{2}(t))-a_{2}(t)x_{3}(t) +\beta(t)x_{1}(t)],\\ \begin{aligned} \frac{dx_{3}(t)}{dt}&=-d(t)x_{3}(t)-b(t)x^{2}_{3}(t)+c_{1}(t)x_{3} (t-\tau)x_{1}(t-\tau)\\ &\quad+c_{2}(t)x_{3}(t-\tau)x_{2}(t-\tau), \end{aligned} \end{gathered} \end{equation} where $x_{1}(t),x_{2}(t)$ and $ x_{3}(t)$ are susceptible, infected and predator population respectively and the corresponding parameters has the meaning as defined in \cite{S}. Time delay is considered as gestation period and disease can be transmitted by contact and spreads among prey species only. Author established some sufficient conditions for the permanence of the system by applying the method of inequality analytical techniques. By the well known method of Lyapunov functional,global asymptotic stability of model \eqref{1.3} has been derived in \cite{S}. Author concluded that the time delay has no effect on the permanence of the system but it has an effect on the global asymptotic stability of model \eqref{1.3}. Model \eqref{1.2} was modified by Hu and Li (2012)\cite{GX} and proposed an autonomous model similar to \eqref{1.3}, their model takes the form \begin{equation}\label{1.4} \begin{gathered} \frac{dS}{dt}=rS({1-\frac{S+I}{k}})-SI\beta-p_{1}SY,\\ \frac{dI}{dt}=-cI+SI\beta-p_{2}IY,\\ \frac{dY}{dt}=-dY+qp_{1}S(t-\tau)Y(t-\tau)+qp_{2}I(t-\tau)Y(t-\tau), \end{gathered} \end{equation} where $S(t),I(t)$ and $ Y(t)$ are susceptible, infected and predator population respectively and parameters used has the meaning as defined in \cite{GX}. They derived stability results and investigate Hopf-bifurcation analysis. They performed stability analysis by using Routh-Hurwitz criteria. The effect of delay on model \eqref{1.4} is considered as a bifurcation parameter for the purpose of the stability of the positive equilibrium. They investigated the Hopf bifurcation. By applying the normal form theory and the center manifold reduction method, the direction of Hopf bifurcations and the stability of bifurcated periodic solutions has been determined in \cite{GX}. The main motivation of the present study is to modify the models \eqref{1.2} and \eqref{1.4} by introducing suitable ecological and biological assumptions. We study the role of time delay as bifurcation parameter by using the normal form theory, Riesz representation theorem and central manifold theorem. The parameters are time independent as considered in \cite{GX,JO}. We have also analyzed the model with and without delay. Detailed ecological and biological assumptions for model formulation are listed in the next section. Rest of the paper is organized as follows. Section 2 deals with mathematical model formation with help of some ecological and biological assumptions. In Section 3 we determine the stability of different equilibrium points for mathematical model without delay. In Section 4 we determine the stability of different equilibrium points for mathematical model with delay. In Section 5 we investigate Hopf-bifurcation and direction of the Hopf-bifurcation including stability of bifurcated periodic solutions. To verify the theoretical frame work, in Section 6 some numerical computation has been performed by considering suitable parameters and initial conditions followed by discussion and future directions in the last Section 7. \section{The model} For mathematical simplicity we impose the following ecological and biological assumptions: \begin{itemize} \item[(A1)] We consider linear functional response as described in \cite{VV}. \item[(A2)] In the absence of disease and predation, prey population follow the logistic rule with the growth rate $r$ $(r>0)$ and carrying capacity $k$ $(k>0)$ \cite{GX}. \item[(A3)] In the presence of disease, prey population is divided into two parts: susceptible (S) and infective (I). Hence, total biomass of the prey population is $S(t)+I(t)$. \item[(A4)] It is considered that by means of contact, disease spreads among the prey species only. \item[(A5)] Only the susceptible prey is assumed to be reproducing offsprings with logistic law i.e. only S has growth rate. However, infected prey population contributes to the carrying capacity. \item[(A6)] Prey population may have possible source of infection (external source) viz. viruses and other seasonal effects. After infection they converted into infected prey (I). The disease dynamics has been omitted. \item[(A7)] Prey population (susceptible(S) and infective (I)) and predator population remains in the same environmental conditions and in the same terrestrial area and zone i.e in same ecosystem. In other words migration (in and out both) has been omitted here. Detail classification of an ecosystem has been ignored. \item[(A8)] It is also assumed that infected prey has high probability of being predated (eaten) by the predator as compare to susceptible prey population. One of the reason of this may be that healthy prey population is more active than infected one. \item[(A9)] It is also assumed that the coefficient of conversing of both the prey to predator are different. One is S-prey to predator and other one from I-prey to predator. \item[(A10)] It is assumed that all the three species susceptive prey, infected prey and predator have their natural death rates. \item[(A11)] Infected Prey has no growth i.e. they are declining only. \item[(A12)] Motivated by studies in \cite{ML, KB, MC, DG,ME}, that linear mass-action incidence is more appropriate than a proportional mixing one in case of direct transmission, we assume that the infection follow the simple law of mass action of the form $\beta SI$ where $\beta$ is the force of infection. \item[(A13)] Initially there may not be infected Prey. It is also assumed that infected prey neither recover nor immune. \item[(A14)] Time delay ($\tau$) is the gestation period of predator. \end{itemize} In view of above assumptions, model \eqref{1.4} takes the form \begin{equation}\label{2.1} \begin{gathered} \frac{dS}{dt}=rS(1-\frac{S+I}{k})-SI\beta-p_{1}SY,\\ \frac{dI}{dt}=SI\beta-p_{2}IY-(d_{2}+c)I,\\ \frac{dY}{dt}=-d_{3}Y+q_{1}p_{1}S(t-\tau)Y(t-\tau)+q_{2}p_{2}I(t-\tau)Y(t-\tau). \end{gathered} \end{equation} We summarize the various nomenclature in Table \ref{table1}. \begin{table}[htbp] \caption{Biological/ecological meaning of the symbols} \label{table1} \begin{center} \renewcommand{\arraystretch}{1.3} \begin{tabular}{|c|l|} \hline $S(t)$& Susceptible(healthy) prey population\\ $I(t)$& Infected prey population\\ $Y(t)$& Predator population\\ $\beta$ & Disease contact rate (force of infection) \\ $p_{1},p_{2}$ & Predation coefficients of susceptible (S) and infected prey (I)\\ $r$&Intrinsic growth rate\\ $k$& Carrying capacity\\ $\tau$ & Gestation period(time delay)\\ $c$ & Death rate of infected prey due to disease \\ $d_{2}$ & Natural death rate of infected prey \\ $d_{3}$ & Natural death rate of predator \\ $q_{1}$ & Coefficient of conversing susceptible prey into predator \\ $q_{2}$ & Coefficient of conversing infected prey into predator\\ \hline \end{tabular} \end{center} \end{table} On the basis of ecological and biological assumption that healthy prey are more active as compare to infected one, the relationship between $q_{1}$ and $q_{2}$ is established as under: \begin{gather*} q_{2} \neq q_{1} \text{ and }0 q_{1} \text{ and }00$, $I(0)=\phi_{2}\geq 0$, $Y(0)=\phi_{3}>0$, where \begin{equation*} \big\{\phi\in\mathbb{C_{+}}: \phi=(\phi_{1},\phi_{2},\phi_{3})\big\}, \end{equation*} where $\mathbb{C_{+}}$ is the Banach space of positive continuous functions $\phi:[-\tau,0]\to{R}_{+}^{3}$ with norm \begin{gather*} \sup_{ [-\tau,0]} \{|\phi_{1}|,|\phi_{2}|,|\phi_{3}|\}, \\ \mathbb{R}_{+}^{3}=\big\{\phi\in\mathbb{C}_{+}: \phi_{i}\geq0,\phi=(\phi_{1},\phi_{2},\phi_{3}),\, i=1,2,3 \big\}. \end{gather*} Model \eqref{1.3} is different with the proposed model \eqref{2.1} in the sense that parameters in \eqref{1.3} are time dependent as contrary to those in \eqref{2.1}. \section{ Model without delay} In absence of time delay $\tau$, model \eqref{2.1} takes the form \begin{equation}\label{3.1} \begin{gathered} \frac{dS}{dt}=rS(1-\frac{S+I}{k})-SI\beta-p_{1}SY,\\ \frac{dI}{dt}=SI\beta-p_{2}IY-(d_{2}+c)I,\\ \frac{dY}{dt}=-d_{3}Y+q_{1}p_{1}S(t)Y(t)+q_{2}p_{2}I(t)Y(t). \end{gathered} \end{equation} \subsection{Equilibria and their feasibility} Model \eqref{3.1} has the following equilibrium points \begin{enumerate} \item $E_{1}=(0,0,0)$, which is trivial equilibrium. \item $E_{2}=(k,0,0)$, this provides the case where prey is infection free and predator is absent. This is called boundary equilibrium. \item $E_{3}=(\widehat{S},0,\widehat{Y})$, where $\widehat{S}$ and $\widehat{Y}$ are given by \begin{equation}\label{3.2} \begin{gathered} \widehat{S}=\frac{d_{3}}{q_{1}p_{1}},\\ \widehat{Y}=\frac{r}{p_{1}}(1-\frac{d_{3}}{kq_{1}p_{1}}), \end{gathered} \end{equation} this provides the case where prey is infection free. \item $E_{4}=(\overline{S},\overline{I},0)$, where $\overline{S}$ and $\overline{I}$ are given by: \begin{equation}\label{3.3}\begin{gathered} \overline{S}=\frac{c+d_{2}}{\beta},\\ \overline{I}=\frac{r(k\beta-c-d_{2})}{\beta(r+k\beta)} \end{gathered} \end{equation} this provides the case where predator is absent. \item $E_{5}=(\widetilde{S},\widetilde{I},\widetilde{Y})$, where $\widetilde{S},\widetilde{I},\widetilde{Y}$ are given by \begin{equation}\label{3.4} \begin{gathered} \widetilde{I}=\frac{d_{3}-q_{1}p_{1}\widetilde{S}}{q_{2}p_{2}},\\ \widetilde{Y}=\frac{\beta \widetilde{S}-c-d_{2}}{p_{2}}, \end{gathered} \end{equation} and \begin{equation}\label{3.5} \begin{gathered} \widetilde{S}A+B=0,\\ A=\big[-\frac{r}{K}+(\frac{r}{K}+\beta)\frac{q_{1}p_{1}}{q_{2}p_{2}} -\frac{p_{1}\beta}{p_{2}}\big],\\ B=\big[r-(\frac{r}{K}+\beta)\frac{d_{3}}{q_{2}p_{2}} +\frac{p_{1}(c+d_{2})}{p_{2}}\big]. \end{gathered} \end{equation} Set $E_{5}$ provides the coexistence of all the three species. Existence (feasibility) conditions of equilibrium points of model ref{3.1} are listed in Table \ref{table2}. \end{enumerate} \begin{table}[htbp] \caption{Existence conditions of equilibrium points of model \eqref{3.1}} \label{table2} \begin{center} \renewcommand{\arraystretch}{1.3} \begin{tabular}{|c|l|} \hline Equilibrium Point & Existence Condition\\ \hline $E_{1}=(0,0,0)$ & always\\ $E_{2}=(k,0,0)$ & always\\ $E_{3}=(\widehat{S},0,\widehat{Y})$ & $kq_{1}p_{1}>d_{3}$\\ $E_{4}=(\overline{S},\overline{I},0)$ & $k\beta>(c+d_{2})$\\ $E_{5}=(\widetilde{S},\widetilde{I},\widetilde{Y})$ & $(d_{3}-q_{1}p_{1}\widetilde{S})>0$, $(\beta\widetilde{S}-c-d_{2})>0$, \\ & either $A<0$ or $B<0$ but not both.\\ \hline \end{tabular} \end{center} \end{table} \begin{remark} \label{rmk1}\rm From Table \ref{table2} it is observed that: \begin{itemize} \item[(i)] The existence of equilibrium points $E_{1}$, $E_{2}$ and $E_{4}$ is independent of parameters $q_{1}$ and $q_{2}$. \item[(ii)] The existence of equilibrium point $E_{3}$ is dependent on $q_{1}$. \item[(iii)] The existence of equilibrium point $E_{5}$ is dependent on $q_{1}$ and $q_{2}$ both. \end{itemize} \end{remark} \subsection{Stability analysis} The variational matrix is given by \[ J = \begin{bmatrix} (r-\frac{2rS}{k}-\frac{rI}{k}-\beta I-p_{1}Y) & (-\frac{rS}{k}-\beta S) & (-p_{1}S)\\ (\beta I) & (\beta S-p_{2} Y-c-d_{2}) & (-p_{2}I) \\ (q_{1}p_{1}Y) & (q_{2}p_{2}Y) & (q_{1} p_{1}S+q_{2}p_{2}I-d_{3}) \end{bmatrix}. \] It is very easy to prove from the above equality that the equilibrium points $E_{1}$, $E_{2}$, $E_{3}$ and $E_{4}$ are unstable. Now, the jacobian matrix at $E_{5}$ is \begin{align*} &J(E_{5}) \\ &= \begin{bmatrix} (r(1-\frac{2\widetilde{S}}{k})-\widetilde{I}(\frac{r}{k}+\beta)-p_{1}\widetilde{Y}) & (-\widetilde{S}(\frac{r}{k}+\beta)) & (-p_{1}\widetilde{S})\\ (\beta \widetilde{I})& (\beta\widetilde{S}-c-d_{2}-p_{2}\widetilde{Y}) & (-p_{2}\widetilde{I}) \\ q_{1}p_{1}\widetilde{Y}& q_{2}p_{2}\widetilde{Y} & (q_{1}p_{1}\widetilde{S}+q_{2}p_{2}\widetilde{I}-d_{3}) \end{bmatrix} \end{align*} and the characteristics equation of $J(E_{5})$ is \begin{equation*} \lambda^{3}+C_{1}\lambda^{2}+C_{2}\lambda+C_{3}=0. \end{equation*} By the Routh-Hurwitz criteria, we can conclude that equilibrium $E_{5}$ is locally stable provided the following conditions are satisfied \begin{equation}\label{3.7} \begin{gathered} C_{i}>0, i=1, 2, 3\\ C_{1}C_{2}-C_{3}>0. \end{gathered} \end{equation} The values of $C_{i}>0, i=1, 2, 3$ are listed at Appendix 1. \begin{remark} \label{rmk2} \rm Stability of the non zero equilibrium point $E_{5}$ depends on $C_{i}>0$, $i=1, 2, 3$ (from Eq. \eqref{3.7}). Since $C_{i}>0$, $i=1, 2, 3$ involves $q_{1}$ and $q_{2}$ both (from Appendix 1). Hence, stability of $E_{5}$ depends on $q_{1}$ and $q_{2}$ both. \end{remark} \section{Model with time delay} Ecologically it is a fact that reproduction of predator after predation is not instantaneous but it will mediated by some time lag, so it may call as gestation period. We record this gestation period as time delay ($\tau$) in our proposed model \eqref{2.1}. It is also clear that since delay is gestation period of predator, hence delay term ($\tau$) appears only in last equation of model \eqref{2.1}. Time delay played a crucial role in analysis. In this section, we will observe the role of time delay. \subsection{Equilibria and their feasibility} It is remarkable that the two models \eqref{2.1} and \eqref{3.1} have the same equilibrium points ecologically. Because of the mathematical point of view we denote them differently. In $R_{+}^{3}$ the system \eqref{2.1} has several possible stationary states (equilibrium points) and they are summarized in Table \ref{table3}. Table \ref{table3} provides a brief ecological meaning of equilibrium points and their applications to real ecosystems. \begin{table}[htbp] \caption{Possible equilibrium points of model \eqref{2.1}} \label{table3} \begin{center} \renewcommand{\arraystretch}{1.3} \begin{tabular}{|l|l|l|} \hline Equilibr. point & Name & Ecological meaning\\ \hline $E_{10}=(0,0,0)$& Trivial & Species will die out. Ecologically not important\\ $E_{11}=(k,0,0)$ & Boundary& Only sound prey survive.\\ && Infection fee. Predator will die out \\ $E_{12}=(\overline{S},\overline{I},0)$ & Boundary & Predator will die out. Infection exists\\ $E_{13}=(\widetilde{S},0,\widetilde{Y})$ & Boundary& No infection. Co-existence of sound prey\\ && and predator \\ $E_{*}=(S_{*},I_{*},Y_{*}) $ & Non zero & Co-existence of all species. \\ & (interior)& Ecologically very important \\ \hline \end{tabular} \end{center} \end{table} \subsection{Stability Analysis} $E_{13}$ and $E_{*}$ play an important role in the controlling of epidemic. The variational matrix of system \eqref{2.1} is written as \begin{align*} J &= \begin{bmatrix} (r(1-\frac{S+I}{k})-\frac{rS}{k}-\beta I-p_{1}Y) & (-\frac{rS}{k}-\beta S) & (-p_{1}S)\\ (\beta I) & (\beta S-p_{2} Y-c-d_{2}) & (-p_{2}I) \\ 0 & 0 & (-d_{3}) \end{bmatrix} \\ &\quad +\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & 0 \\ (q_{1}p_{1}Y) & (q_{2}p_{2}Y) & (q_{1} p_{1}S+q_{2}p_{2}I) \end{bmatrix} e^{-\lambda\tau}, \end{align*} where $\lambda$ being a complex number. In simplified form the above equation may be written as \[ J = \left[\begin{smallmatrix} (r(1-\frac{S+I}{k})-\frac{rS}{k}-\beta I-p_{1}Y) & (-\frac{rS}{k}-\beta S) & (-p_{1}S)\\ (\beta I) & (\beta S-p_{2} Y-c-d_{2}) & (-p_{2}I) \\ (q_{1}p_{1}Y) (e^{-\lambda\tau})& (q_{2}p_{2}Y)(e^{-\lambda\tau}) & (-d_{3})+(q_{1} p_{1}S+q_{2}p_{2}I)(e^{-\lambda\tau}) \end{smallmatrix}\right]. \] We will use the following lemma due to Hu and Li \cite{GX}. \begin{lemma} \label{lem1} Let $A>0, B>0$. Then \begin{itemize} \item If $AB$, all roots of the equation $\lambda+A-B e^{-\lambda\tau}=0$ have negative real parts for any $\tau$. \end{itemize} \end{lemma} Now we study the dynamical behavior of system \eqref{2.1} about different equilibrium points with the help of variational matrix. \subsubsection{Trivial equilibrium point ($E_{10}$) } Proceeding as in Sub-section 3.2.1, it is concluded that $E_{10}$ is unstable. \subsubsection{Boundary equilibrium point ($E_{11}$)} Variational matrix evaluated at $E_{11}$ takes the form \begin{equation*} J(E_{11})=\begin{bmatrix} -r & -(r+k\beta) & (-p_{1}k)\\ 0 & (\beta k-c-d_{2}) & 0\\ 0 & 0 & (q_{1} p_{1}k)e^{-\lambda\tau}-d_{3} \end{bmatrix}, \end{equation*} with corresponding eigenvalues $-r$, $(\beta k-c-d_{2})$ and $(q_{1} p_{1}k)e^{-\lambda\tau}-d_{3}$. Hence, the stability of $E_{11}$ depends on $(\beta k-c-d_{2})$ and $(q_{1} p_{1}k)e^{-\lambda\tau}-d_{3}$. Now, if the following condition is satisfied \begin{equation}\label{4.2} (q_{1} p_{1}k)0$, then $E_{11}$ is unstable always. \subsubsection{Predator free equilibrium ($E_{12}$)} Now as in previous section, we have $\overline{S}=\frac{(c+d_{2})}{\beta}$ and $ \overline{I}=\frac{r(k\beta-c-d_{2})}{\beta(r+k\beta)}$. The variational matrix at $E_{12}$ takes the form \begin{equation*} J(E_{12})=\begin{bmatrix} (r(1-\frac{\overline{S}+\overline{I}}{k})-\frac{r \overline{S}}{k}-\beta \overline{I}) & (-\frac{r \overline{S}}{k}-\beta \overline{S}) & (-p_{1}\overline{S})\\ (\beta \overline{I}) & (\beta \overline{S}-c-d_{2}) & (-p_{2} \overline{I})\\ 0 & 0 & (q_{1} p_{1}\overline{S}+q_{2} p_{2}\overline{I})e^{-\lambda\tau}-d_{3} \end{bmatrix}, \end{equation*} and the characteristics equation corresponding to $J(E_{12})$ is \begin{align*} &[\lambda+ d_{3}-((q_{1} p_{1}\overline{S}+q_{2} p_{2}\overline{I}) e^{-\lambda\tau})][\lambda^{2}+\lambda \bigl(c+d_{2}+\frac{r(c+d_{2})}{k \beta} +\frac{\beta}{(c+d_{2})}-\frac{(c+d_{2})}{k}\bigr) \\ &-\frac{(c+d_{2})^{2}r}{k \beta}-\beta-r(c+d_{2})]=0. \end{align*} If the following condition is satisfied \begin{equation*} d_{3}>(q_{1} p_{1}\overline{S}+ q_{2} p_{2}\overline{I}), \end{equation*} then by Lemma \ref{lem1}, the root of the functions \begin{equation*} \lambda+d_{3}-((q_{1} p_{1}\overline{S}+q_{2} p_{2} \overline{I})e^{-\lambda\tau})), \end{equation*} will have negative real part for any value of $\tau$ and for the equation \begin{align*} &\lambda^{2}+\lambda \Big((c+d_{2})+\frac{r (c+d_{2})}{K\beta} +\frac{\beta}{(c+d_{2})}-\frac{(c+d_{2})}{k}\Big)\\ &+\Big(-\frac{r(c+d_{2})^{2}}{k\beta}-\beta-r(c+d_{2})\Big)=0, \end{align*} the Routh-Hurwitz criteria, may be used for proving the fact that this equation will have roots with negative real parts. Hence, if the equilibrium $E_{12}$ is asymptotically stable, it will mean that predator population will be die-out from system so considered. \subsubsection{Infection free equilibrium $(E_{13})$} $\widetilde{S}$ and $\widetilde{Y}$ are given by $\widetilde{S}=\frac{d_{3}}{q_{1}p_{1}}, \widetilde{Y} =[\frac{r}{p_{1}}(1-\frac{d_{3}}{kq_{1}p_{1}})]$. The variational matrix at $E_{13}$ takes the form \begin{equation*} J(E_{13})=\begin{bmatrix} (r(1-\frac{\widetilde{S}}{k})-\frac{r \overline{S}}{k}-p_{1} \widetilde{Y}) & (-\frac{r \widetilde{S}}{k}-\beta \widetilde {S}) & (-p_{1}\widetilde{S})\\ 0 & (\beta \widetilde{S}-p_{2} \widetilde{Y}-c-d_{2}) & 0\\ (q_{1} p_{1} \widetilde{Y}e^{-\lambda\tau}) & (q_{2} p_{2} \widetilde{Y}e^{-\lambda\tau}) & (q_{1} p_{1}\widetilde{S}e^{-\lambda\tau}-d_{3}) \end{bmatrix}, \end{equation*} One of the eigenvalue is $(\beta\widetilde{S}-p_{2}\widetilde{Y}-c-d_{2})$ and two other eigenvalues are the roots of the expression \begin{align*} &\Big[\lambda^{2}-\lambda\bigl(r(1-\frac{2\widetilde{S}}{k})-p_{1} \widetilde{Y}-d_{3}+(q_{1} p_{1} \widetilde{Y}e^{-\lambda\tau})\bigr) +\bigl(-r d_{3}(1-\frac{2 \widetilde{S}}{k})-d_{3}p_{1}\widetilde{Y}\\ &+d_{3} +r q_{1} p_{1}\widetilde{S}(1-\frac{2\widetilde{S}}{k}) e^{-\lambda\tau}\bigr)\Big]. \end{align*} If the condition \begin{equation*} \widetilde{S}>\frac{(c+d_{2}+p_{2}\widetilde{Y})}{\beta} \end{equation*} is satisfied. Then one eigenvalue $(\beta \widetilde{S}-p_{2} \widetilde{Y}-c-d_{2})$ corresponding to $J(E_{13})$ is positive. Hence, in this case $E_{13}$ is unstable. Let us put, $\lambda=(u+iv)$ with condition for $u$ as $u\geq0$ in the expression \begin{align*} &\Big[\lambda^{2}-\lambda\bigl(r(1-\frac{2\widetilde{S}}{k})-p_{1} \widetilde{Y}-d_{3}+(q_{1} p_{1} \widetilde{Y}e^{-\lambda\tau})\bigr) +\bigl(-r d_{3}(1-\frac{2 \widetilde{S}}{k})\\ &-d_{3}p_{1}\widetilde{Y}+d_{3} +r q_{1} p_{1}\widetilde{S}(1-\frac{2\widetilde{S}}{k}) e^{-\lambda\tau}\bigr)\Big], \end{align*} and on separating real and imaginary parts, we obtained \begin{align*} \text{Real part} &= (u^{2}-v^{2})-\bigl(r(1-\frac{2\widetilde{S}}{k})-p_{1} \widetilde{Y}-d_{3}+(q_{1} p_{1} \widetilde{S}e^{-\lambda\tau} \cos v\tau\bigr) \\ &\quad +\bigl(r q_{1} p_{1}\widetilde{S}(1-\frac{2\widetilde{S}}{K})e^{-\lambda\tau} \cos v\tau \bigr), \end{align*} \begin{align*} \text{Imaginary part} &=2uv+\Big(p_{1} \widetilde{Y}+d_{3}-vr (1-\frac{2\widetilde{S}}{k}) -q_{1} p_{1} \widetilde{S}e^{-\lambda\tau} \cos v\\ &\quad -uq_{1} p_{1} \widetilde{S}e^{-\lambda\tau} \sin v\tau-r q_{1} p_{1}\widetilde{S}(1-\frac{2\widetilde{S}}{k})e^{-\lambda\tau} \sin v\tau\Big). \end{align*} Now, if the following condition is satisfied \begin{equation*} \widetilde{S}<\frac{(c+d_{2}+p_{2}\widetilde{Y})}{\beta}, \end{equation*} and the real part is negative, we can conclude that equilibrium state $E_{13}$ is stable. \subsubsection{Non zero equilibrium $(E_{*}({S}_{*},{I}_{*},{Y}_{*}))$} $S_{*}, I_{*},Y_{*}$ are given by \begin{gather*} S_{*}=\frac{-B}{A}, I_{*}=\frac{d_{3}-q_{1}p_{1}S_{*}}{q_{2}p_{2}}, Y_{*}=\frac{\beta S_{*}-c-d_{2}}{p_{2}},\\ A=\Big[-\frac{r}{k}+(\frac{r}{k}+\beta)\frac{q_{1}p_{1}}{q_{2}p_{2}} -\frac{p_{1}\beta}{p_{2}}\Big] B=\Big[(r) -(\frac{r}{K}+\beta)\frac{d_{3}}{q_{2}p_{2}}+\frac{p_{1}(c+d_{2})}{p_{2}}\Big]. \end{gather*} The variational matrix at $E_{*}$ takes the form \[ J(E_{*}) =\left[\begin{smallmatrix} (r(1-\frac{S_{*}+I_{*}}{k})-\frac{r S_{*}}{k}-p_{1}Y_{*}-\beta I_{*}) & (-\frac{r S_{*}}{k}-\beta S_{*}) & (-p_{1}S_{*})\\ \beta I_{*} & (\beta S_{*}-p_{2} Y_{*}-c-d_{2}) & p_{2}I_{*}\\ (q_{1} p_{1} Y_{*}e^{-\lambda\tau}) & (q_{2} p_{2}Y_{*}e^{-\lambda\tau}) & (q_{1} p_{1}S_{*}+q_{2} p_{2}I_{*})e^{-\lambda\tau}-d_{3}) \end{smallmatrix}\right], \] and the characteristics equation corresponding to $J(E_{*})$ is \begin{equation}\label{4.3} (\lambda^{3}+m_{2}\lambda^{2}+m_{1}\lambda+m_{0}) +(n_{2}\lambda^{2}+n_{1}\lambda+n_{0})e^{-\lambda\tau}=0, \end{equation} where $m_{i}, n_{j}$, $i=0,1,2$; $j=0,1,2$ are listed in Appendix 2. Now we put $\lambda=i\omega$ $(\omega>0)$ in the above equation and separating real and imaginary parts, we obtain \begin{equation}\label{4.4} \begin{gathered} \text{Real part}=\{n_{2}\omega^{2}+n_{0}\}\cos{\omega\tau} +\{n_{1}\omega\sin\omega\tau-m_{2}\omega^{2}+m_{0}\},\\ \text{Imaginary part}=n_{1}\omega\cos\omega\tau-(-n_{2}\omega^{2} +n_{0})\sin\omega\tau+m_{1}\omega-\omega^{3},\\ (\text{Real part})^{2}+(\text{Imaginary part})^{2} =\omega^{6}+p_{0}\omega^{4}+q_{0}\omega^{2}+r_{0}, \end{gathered} \end{equation} where \begin{gather*} p_{0}=(m_{2}^{2}-2 m_{1}-n_{2}^{2}),\\ q_{0}=(m_{1}^{2}-2 m_{2}m_{0}+2 n_{2}n_{0}-n_{1}^{2}),\\ r_{0}=(m_{0}^{2}-n_{0}^{2}). \end{gather*} We refer the following lemma due to \cite{GX, YX} \begin{lemma} \label{lem2} For the polynomial \begin{equation}\label{4.5} h(z)=z^{3}+p_{0}z^{2}+q_{0}z+r_{0}=0, \end{equation} \begin{itemize} \item[(i)] If $r_{0}<0$, then this equation has at least one positive root; \item[(ii)] If $r_{0}\geq0$ and $\triangle=(p_{0}^{2}-3q_{0})\leq 0$, then this equation has no positive roots; \item[(iii)] If $r_{0}\geq0$ and $\triangle=(p_{0}^{2}-3q_{0})>0$, then this equation has positive roots if and only if $z_{1}^{*}=\frac{-p_{0}+\sqrt{\triangle}}{3}$ and $h(z_{1}^{*})\leq0$. \end{itemize} \end{lemma} If we put $z=\omega^{2}$ in $\omega^{6}+p_{0}\omega^{4}+q_{0}\omega^{2}+r_{0}=0$, then we have the equation $z^{3}+p_{0}z^{2}+q_{0}z+r_{0}=0$. If $m_{0}^{2}\geq n_{0}^{2}$, then we will have $r_{0}\geq0$, we have two situations for $\Delta$: \begin{itemize} \item[(i)] $\Delta=(p_{0}^{2}-3q_{0})\leq 0,$ \item[(ii)] $\Delta=(p_{0}^{2}-3q_{0})>0$. \end{itemize} In situation (i), we have to say that $E_{*}$ is stable thus $E_{*}$ is absolutely stable if $r_{0}\geq0$ and $\Delta=(p_{0}^{2}-3q_{0})\leq 0$ and also if we have and $r_{0}\geq0$ and $\Delta=(p_{0}^{2}-3q_{0})>0$ then equation has negative roots if and only if $h(z_{1}^{*})>0$ where $z_{1}^{*}=\frac{-p_{0}+\sqrt{\triangle}}{3}$, thus we have the following theorem for the stability of $E_{*}$ \begin{theorem} \label{thm1} Equilibrium $E_{*}(S_{*},I_{*},Y_{*})$ is absolutely stable if one of the following three conditions holds \begin{itemize} \item[(i)] $\Delta=(p_{0}^{2}-3q_{0})\leq 0$; \item[(ii)] $\Delta=(p_{0}^{2}-3q_{0})>0$ and $z_{1}^{*}=\frac{-p_{0}+\sqrt{\triangle}}{3}<0$; \item[(iii)] $\Delta=(p_{0}^{2}-3q_{0})>0$, $z_{1}^{*}=\frac{-p_{0}+\sqrt{\triangle}}{3}>0$ and $h(z_{1}^{*})>0$ provided $r_{0}\geq 0$. \end{itemize} \end{theorem} Next if we consider the case when $r_{0}<0$ or $\{r_{0}\geq0,\Delta=(p_{0}^{2}-3q_{0})>0,z_{1}^{*}>0 ,h(z_{1}^{*})<0\}$. Then according to lemma \ref{lem2}, \eqref{4.4} will have one positive root say $\omega_{0}$ that is the characteristic equation has a pair of purely imaginary roots say $\pm i \omega_{0}$. Now assume that $\ i \omega_{0}$, $\omega_{0}>0$ is a root of $h(z)$. Solving the eq. \eqref{4.4} for $\tau$, we have (by eliminating $\sin\omega\tau$, we obtain \begin{equation}\label{4.6} \tau=\frac{1}{\omega_{0}}\cos^{-1}(\frac{n_{1}\omega_{0}^{2} \{\omega_{0}-m_{1}\}-\{m_{2}\omega_{0}^{2}-m_{0}\}\{n_{2} \omega_{0}^{2}-n_{0}\}}{n_{1}^{2}\omega_{0}^{2}+n_{2}\omega_{0}^{2}-n_{0}}) +\frac{2k\pi}{\omega_{0}}, \end{equation} for $k=0,1,2,\dots$. We call it as a `critical value' and is denoted by \begin{equation}\label{4.7} \tau_{k}=\frac{1}{\omega_{0}}\cos^{-1}(\frac{n_{1}\omega_{0}^{2}\{\omega_{0} -m_{1}\}-\{m_{2}\omega_{0}^{2}-m_{0}\}\{n_{2} \omega_{0}^{2}-n_{0}\}}{n_{1}^{2}\omega_{0}^{2} +n_{2}\omega_{0}^{2}-n_{0}})+\frac{2k\pi}{\omega_{0}}, \end{equation} for $k=0,1,2,\dots$. This corresponds to the characteristic equation that has purely imaginary roots $\pm i \omega_{0}$. Which is a result similar to that is discussed in \cite{GX}. Transversality condition may also obtained as discussed in \cite{GX}. As discussed in \cite[Theorem 2.4]{GX}, the equilibrium point $E_{*}$ of the system \eqref{2.1} is asymptotically stable when $\tau>\tau_{0}$. $\tau=\tau_{k}$ $(k=0,1,2,3,\dots)$ are Hopf-bifurcation values for the system \eqref{2.1} and $\tau_{k}$ is used as a point for direction of Hopf Bifurcation in next section. \begin{remark} \label{rmk3} \rm From \eqref{4.6} and \eqref{4.7}, it is observed that the delay term $\tau$ (here bifurcation parameter) depends on the values of $m_{i}, n_{j}$, $i=0,1,2$; $j=0,1,2$. Since $m_{i}, n_{j}$, $i=0,1,2$; $j=0,1,2$ depends on $q_{1}$ and $q_{2}$ (see Appendix 2), hence bifurcation parameter ($\tau$) depends on $q_{1}$ and $q_{2}$ both. A little variation in the values of $q_{1}$ and $q_{2}$ may change the bifurcation parameter ($\tau$). Hence, a little variation in the values of $q_{1}$ and $q_{2}$ may change the dynamics of the delayed model \eqref{2.1}. \end{remark} \section{Direction and stability of the Hopf Bifurcation} With the symbols used in \cite{GX} and the procedure explained in \cite{CC}. System \eqref{2.1}, can be translated to the following functional differential equation (FDE) system \begin{equation}\label{5.1} \dot{u}(t)=L_{\mu}(\mu_{t})+F(\mu,u_{t}), \end{equation} where $u_{t}=u(t)\in R^{3}$ and $L_{\mu}:\mathbb{R}\times\mathbb{C}\to\mathbb{R}^{3}$ and $F:\mathbb{R}\times\mathbb{C}\to\mathbb{R}^{3}$ are given by, \begin{gather*} L_{\mu}\phi=(\tau_{k}+\mu)(M_{1}\phi(0)+M_{2}\phi(-1)), \\ F(\mu,\theta)= \begin{bmatrix} -\frac{r}{k}\phi_{1}^{2}(0) -(\frac{r}{k}+\beta)\phi_{1}(0)\phi_{2}(0)-p_{1}\phi_{1}(0)\phi_{3}(0) \\ \beta \phi_{1}(0)\phi_{2}(0)- p_{2}\phi_{2}(0)\phi_{3}(0)\\ q_{1} p_{1}\phi_{1}(-1)\phi_{3}(-1)+ q_{2} p_{2}\phi_{1}(-1)\phi_{2}(-1) \end{bmatrix}, \end{gather*} where \begin{gather*} M_{1}=\begin{bmatrix} (r-\frac{2rS_{*}}{k}-(\frac{r}{k}+\beta)I_{*}-p_{1}Y_{*}) & -(\frac{r}{k}+\beta)S_{*} & (-p_{1}S_{*})\\ \beta I_{*} & (\beta S_{*}-p_{2} Y_{*}-c-d_{2}) & -p_{2}I_{*}\\ 0& 0 & -d_{3}) \end{bmatrix},\\ M_{2}=\begin{bmatrix} 0 & 0& 0\\ 0& 0 & 0\\ q_{1} p_{1}Y_{*}& q_{2} p_{2}Y_{*} & q_{1} p_{1}S_{*}+q_{2} p_{2}I_{*} \end{bmatrix},\\ \phi(0)=(\phi_{1}(0),\phi_{2}(0),\phi_{3}(0))^{T}\in\mathbb{C},\\ \phi(-1)=(\phi_{1}(-1),\phi_{2}(-1),\phi_{3}(-1))^{T}\in\mathbb{C}. \end{gather*} We have considered, $\tau=(\tau_{k}+\mu)$, $\mu=0$ which gives the Hopf bifurcation value for the mathematical model with delay as promised in previous section. Normalizing delay $\tau$ by the time scaling $t\to\frac{t}{\tau}$ then the model is written in the Banach Space $\mathbb{C}\equiv\mathbb{C}([-1,0],\mathbb{R}^{3})$. By the Riesz representation theorem, we found that there exists a matrix function whose components are bounded variation function $\eta(\theta,\mu)$ in $\theta\in[-1,0]$, such that $L_{\mu}\phi=\int_{\Omega}d\eta(\theta,\mu)\phi(\theta)$, $\phi\in\mathbb{C}$, $\Omega\in[-1,0)$. We can choose \begin{equation*} \eta(\theta,\mu)=(\tau_{k}+\mu)M_{1}\delta(\theta)- (\tau_{k}+\mu) M_{2}\delta(\theta+1), \end{equation*} where $\delta(\theta)$ denotes the dirac delta function viz. \begin{equation*} \delta(\theta)= \begin{cases} 0, & \theta\neq0,\\ 1, & \theta=0. \end{cases} \end{equation*} For $\phi\in\mathbb{C}^{1}([-1,0],\mathbb{R}^{3})$, we define \begin{align*} A(\mu)\phi(\theta) &=\begin{cases} \frac{d\phi(\theta)}{d\theta}, & \theta\in[-1,0),\\ \int_{-1}^{0}d\eta(\theta,\mu)\phi(\theta), & \theta=0, \end{cases}\\ &= \begin{cases} \frac{d\phi(\theta)}{d\theta}, & -1\leq\theta<0,\\ \int_{-1}^{0}d\eta(\theta,\mu)\phi(\theta), & \theta=0. \end{cases} \end{align*} and \begin{align*} \mathbb{R}(\mu)\phi(\theta) &= \begin{cases} 0, & \theta\in[-1,0),\\ F(\mu,\phi), & \theta=0, \end{cases}\\ &=\begin{cases} 0, & -1\leq\theta<0,\\ F(\mu,\phi), & \theta=0, \end{cases} \end{align*} with these symbols, FDE system \eqref{5.1} may be written in the form \begin{equation*} \dot{u}(t)=A({\mu})(\mu_{t})+\mathbb{R}(\mu)\mu_{t}, \end{equation*} which is an abstract differential equation where $u_{t}(\theta)=u(t+\theta),-1\leq\theta<0$. Now we come to operator theory, for $\psi\in\mathbb{C}^{1}\big([0,1],(\mathbb{R}^{3})^{*}\bigr)$ we define $A^{*}$, the adjoint operator of $A$, by \begin{equation*} A^{*}\psi(S)= \begin{cases} -\frac{d \psi(S)}{dS},& S\in(0,1],\\ \int_{-1}^{0}d\eta^{T}(S,\mu)\psi(-S), & S =0, \end{cases} \end{equation*} and a bilinear product \begin{equation}\label{5.2} \langle \psi(S),\phi(\theta)\rangle =\overline{\psi}(0)\phi(0)-\int_{1}^{0} \int_{\xi=0}^{\theta}\overline{\psi}^{T}(\xi-\theta)d\eta(\theta)\phi(\xi)d\xi, \end{equation} where $\eta(\theta)=\eta(\theta,0)$. Then $A(0)$ and $A_{*}$ are adjoint operators. From previous section, it is noted that $\pm i \omega_{0}\tau_{k}$ are eigenvalues of $A(0)$. Hence they are also the eigenvalues of $A^{*}$. To determine the poincare normal form of the operator $A$, we first need to evaluate the eigenvectors of $A(0)$ and $A^{*}$ corresponding to $i\omega_{0}\tau_{k}$ and $-i \omega_{0}\tau_{k}$ respectively. Suppose that $q(\theta)=(1,\alpha_{1},\alpha_{2})^{T}\exp(i\omega_{0} \tau_{k}\theta)$ is the eigenvector of $A(0)$ corresponding to $i\omega_{0}\tau_{k}$, then we have $A(0)q(\theta)=i\omega_{0}q(\theta)$ from the definition of $A(0)$, we have \begin{equation*} \Big[M_{1}+M_{2}\exp(i\omega_{0}\tau_{k})\Big] \begin{bmatrix} 1 \\ \alpha_{1} \\ \alpha_{2} \end{bmatrix} =i\omega_{0}\begin{bmatrix} 1 \\ \alpha_{1} \\ \alpha_{2} \end{bmatrix}, \end{equation*} or \begin{equation*} M_{1}\begin{bmatrix} 1 \\ \alpha_{1} \\ \alpha_{2} \end{bmatrix}+ M_{2}\begin{bmatrix} \exp(-i\omega_{0}\tau_{k})\\ \alpha_{1}\exp(-i\omega_{0}\tau_{k} )\\ \alpha_{2}\exp(-i\omega_{0}\tau_{k}) \end{bmatrix}=i\omega_{0}\begin{bmatrix} 1 \\ \alpha_{1} \\ \alpha_{2} \end{bmatrix}. \end{equation*} By simple calculation, we obtain \begin{gather*} \alpha_{1}=\frac{-p_{2}I_{*}(i \omega_{0}-(r-\frac{2rS_{*}}{K} -(\frac{r}{k}+\beta)I_{*}-p_{1}Y_{*}))-p_{1}\beta S_{*}I_{*} }{p_{2} (\frac{r}{k}+\beta) S_{*}I_{*}-p_{1}S_{*}(i \omega_{0}-\beta S_{*} +c+d_{2}+p_{2}Y_{*})},\\ \alpha_{2}=\frac{q_{1}p_{1}Y_{*}\exp(-i\omega_{0}\tau_{k} +q_{2}p_{2 }Y_{*}\exp(-i\omega_{0}\tau_{k}}{i \omega_{0}+d_{3} -q_{1}p_{1}S_{*}+q_{2}p_{2}I_{*}}. \end{gather*} Next, suppose that $q_{*}(s)=B(1,\alpha_{1}^{*},\alpha_{2}^{*})\exp(i\omega_{0}\tau_{k}s)$ is the eigenvector of $A^{*}$ corresponding to $-i \omega_{0}\tau{k}$. Analogously, we have \begin{gather*} \alpha_{1}^{*}=\frac{-p_{1}(\frac{r}{k}+\beta)S_{*}-p_{2}(i\omega_{0} -(r-\frac{2rS_{*}}{k}-(\frac{r}{k}+\beta)I_{*}-p_{1}Y_{*}))}{p_{2} \beta I_{*}-p_{1}(i\omega_{0}+\beta S_{*}-c-d_{2}-p_{2}Y_{*})},\\ \alpha_{2}^{*}=\frac{-p_{1}S_{*}-p_{2}I_{*}\alpha_{1}^{*}}{-i\omega_{0}+d_{3} -(q_{1} p_{1}S_{*}+q_{2} p_{2}I_{*})\exp(-i\omega_{0}\tau_{k}}), \end{gather*} where $B$ has to be calculated. We have the two conditions: \begin{equation*} \langle q^{*},q(\theta)\rangle =1, \quad \langle q^{*},\overline{q(\theta)}\rangle=0, \end{equation*} which may be verified. By equation \eqref{5.2}, we have \begin{align*} &\langle q^{*},q(\theta)\rangle \\ &=\overline{q^{*}}(0)q(0)-\int_{-1}^{0}\int_{\xi=0}^{\infty} \overline{q^{*}}^{T}(\xi-\theta)d\eta(\theta)q(\xi)d\xi\\ &=\overline{B}(1,\overline{\alpha_{1}}^{*},\overline{\alpha_{2}}^{*}) (1,\alpha_{1},\alpha_{2})^{T}-\int_{-1}^{0}\int_{\xi=0}^{\infty} \overline{B}(1,\overline{\alpha_{1}}^{*},\overline{\alpha_{2}}^{*}) \exp({-i\omega_{0}\tau_{k}(\xi-\theta))d\eta(\theta)} \\ &\quad \times(1,\alpha_{1},\alpha_{2})^{T}\exp(i\omega_{0}\tau_{k}\xi)d\xi \\ &=\overline{B}\{1+\alpha_{1}\overline{\alpha_{1}}^{*} +\alpha_{2}\overline{\alpha_{2}}^{*}-\int_{-1}^{0} (1,\overline{\alpha_{1}}^{*},\overline{\alpha_{2}}^{*})\exp(i\omega_{0}\tau_{k}) d \eta(\theta)(1,\alpha_{1},\alpha_{2})^{T}\} \\ &=\overline{B}\{1+\alpha_{1}\overline{\alpha_{1}}^{*} +\alpha_{2}\overline{\alpha_{2}}^{*}+\tau_{k}[q_{2}p_{2} \overline{\alpha_{2}}^{*}Y_{*}+q_{2}p_{2}\alpha_{1}\overline{\alpha_{2}}^{*}Y_{*} +(q_{1} p_{1}S_{*}+q_{2} p_{2}I_{*})\alpha_{2}\overline{\alpha_{2}}^{*}]\\ &\times \exp(-i\omega_{0}\tau_{k})\}, \end{align*} and \begin{align*} 1&=\overline{B}\{1+\alpha_{1}\overline{\alpha_{1}}^{*}+\alpha_{2} \overline{\alpha_{2}}^{*}+\tau_{k}\big[q_{2}p_{2} \overline{\alpha_{2}}^{*}Y_{*}+q_{2}p_{2}\alpha_{1} \overline{\alpha_{2}}^{*}Y_{*}\\ &\quad +(q_{1} p_{1}S_{*}+q_{2} p_{2}I_{*})\alpha_{2} \overline{\alpha_{2}}^{*}\big] \exp(-i\omega_{0}\tau_{k})\}, \end{align*} which gives \begin{align*} \overline{B}&=\Big(1+\alpha_{1}\overline{\alpha_{1}}^{*}+\alpha_{2} \overline{\alpha_{2}}^{*}+\tau_{k}[q_{2}p_{2}\overline{\alpha_{2}}^{*}Y_{*} +q_{2}p_{2}\alpha_{1}\overline{\alpha_{2}}^{*}Y_{*}\\ &\quad +(q_{1} p_{1}S_{*}+q_{2} p_{2}I_{*})\alpha_{2}\overline{\alpha_{2}}^{*}] \exp(-i\omega_{0}\tau_{k})\Big)^{-1}. \end{align*} \subsection{Stability of bifurcated periodic solutions} Firstly, we will investigate the coordinates of the center manifold $\mathbb{C}_{0}$ at $\mu=0$. Let $u_{t}$ be the solution of \eqref{5.1} and define, $z(t)=\langle q^{*},u_{t}\rangle $, $q^{*}$ being the eigenvalue of $A^{*}$ and $W(t,\theta)=u_{t}(\theta)-2Re\{z(t)q(\theta)\}$ on the Center Manifold $\mathbb{C}_{0}$, we have \begin{equation*} W(t,\theta)=W(z(t),\overline{z(t)},\theta), \end{equation*} where \begin{equation}\label{5.3} W(z,\overline{z},\theta)=W_{20}(\theta)\frac{z^{2}}{2}+W_{02}(\theta) \frac{\overline{z}^{2}}{2}+W_{11}(\theta)z\overline{z}+W_{30} \frac{z^{3}}{\lfloor{3}}+\dots\, . \end{equation} In fact, $z$ and $\overline{z}$ are local coordinates for the center manifold $\mathbb{C}_{0}$ in the direction of $q^{*}$ and $\overline{q}^{*}$ respectively. The existence of $\mathbb{C}_{0}$ will provide an opportunity to reduce the system $\dot{u}(t)=L_{\mu}(u_{t})+F(\mu,\mu_{t})$ into an ordinary differential equation ODE( in a single complex variable z) on $\mathbb{C}_{0}$ which is very interesting. $u_{t}$ is the solution of system under consideration. $u_{t}\in \mathbb{C}_{0}$, we have \begin{align*} \dot{z}(t)&=\langle q^{*},\dot{u_{t}}\rangle \\ &=\langle q^{*},A( u_{t})+R (u_{t})\rangle \\ &=\langle q^{*},A( u_{t})>++ $, $W(t,\theta)={u}_{t}(\theta)-2Re\{z(t)q(\theta)\}$ gives us \[ \dot{W} =\dot{u_{t}}-zq-\dot{\overline{{z}}}\overline{q}\\ =\begin{cases} AW-2Re{\overline{q}^{*}(0)F_{0}q(\theta)}, & -1\leq\theta<0, \\ AW-2Re{\overline{q}^{*}(0)F_{0}q(\theta)}+F_{0}, & \theta=0. \end{cases} \] Rewriting the above equation, we obtain \begin{equation}\label{5.4} \dot{W}=AW+H(z,\overline{z},\theta), \end{equation} where \begin{equation}\label{5.5} H(z,\overline{z},\theta)=H_{20}(\theta)\frac{z^{2}}{2} +H_{11}(\theta)z\overline{z}+H_{02}(\theta) \frac{\overline{z^{2}}}{2}+H_{21}(\theta)\frac{z^{2}\overline{z}}{2}+\dots\,. \end{equation} Near the origin on $\mathbb{C}_{0}$, we have \begin{equation*} \dot{W}=W_{z}\dot{z}+W_{\overline{z}}\overline{z}, \end{equation*} by \eqref{5.4} and \eqref{5.5}, we have \begin{gather*} (A-2i\omega_{0}\tau_{k})W_{20}{(\theta)}=-H_{20}(\theta), \\ AW_{11}(\theta)=-H_{11}(\theta), \end{gather*} hence for $-1\leq\theta<0$ we have \begin{equation*} H(z,\overline{z},\theta)=-2Re(\overline{q}^{*}(0)F_{0}q(\theta)) =-g(z,\overline{z})q(\theta)-\overline{g}(z,\overline{z})\overline{q}(\theta). \end{equation*} By comparing the coefficients of $z$ in $H(z,\overline{z},\theta)$, we have \begin{gather*} H_{20}(\theta)=-g_{20}q(\theta)-\overline{g}_{02}\overline{q}(\theta),\\ H_{11}(\theta)=-g_{11}q(\theta)-\overline{g}_{11}\overline{q}(\theta). \end{gather*} Now, we have \begin{gather*} \dot{W}_{20}(\theta)=2 i \omega_{0}\tau_{k}{W}_{20}(\theta)+g_{20}q(\theta) +\overline{g}_{02}\overline{q}(\theta),\\ \dot{W}_{11}(\theta)=g_{11}q(\theta)+\overline{g}_{11}\overline{q}(\theta). \end{gather*} On integrating, we have \begin{gather*} {W}_{20}(\theta)=\frac{i g_{20}}{\omega_{0}\tau_{k}}q(0) \exp(i\omega_{0}\tau_{k}\theta)+\frac{i \overline{g}_{20} \overline{q}(0)}{3\omega_{0}\tau_{k}}\exp(-i\omega_{0}\tau_{k}\theta) +E_{1}\exp(2i\omega_{0}\tau_{k}\theta),\\ {W}_{11}(\theta)=\frac{ g_{21}}{i \omega_{0}\tau_{k}}q(0) \exp(i\omega_{0}\tau_{k}\theta)+\frac{i \overline{g}_{11} \overline{q}(0)}{\omega_{0}\tau_{k}}\exp(-i\omega_{0}\tau_{k}\theta)+E_{2}, \end{gather*} where $E_{1}$ and $E_{2}$ are to be determined. From definitions of $A$, the equation $(A-2i \omega_{0}\tau_{k}){W}_{20}(\theta)=-H_{20}(\theta)$ gives us \begin{equation*} \int_{-1}^{0}d\eta(\theta){W}_{20}(\theta)=2i \omega_{0}\tau_{k}{W}_{20}(0)-H_{20}(0), \end{equation*} which gives us \begin{equation*} H_{20}(0)=-g_{20}q(0)-\overline{g}_{02}\overline{q}(0)+2\tau_{k} \begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)\alpha_{1}-p_{1}\alpha_{2}\\ \beta \alpha_{1}-p_{2}\alpha_{1}\alpha_{2}\\ ( q_{1}p_{1}\alpha_{2}+q_{2}p_{2}\alpha_{1})\exp(-2i\omega_{0}\tau_{k}) \end{bmatrix}. \end{equation*} Now, \begin{gather*} \Big(i\omega_{0}\tau_{k}I-\int_{-1}^{0}\exp(i\omega_{0}\tau_{k}\theta) d\eta(\theta)\Big)q(0)=0, \\ \Big(-i\omega_{0}\tau_{k}I-\int_{-1}^{0} \exp(-i\omega_{0}\tau_{k}\theta)d\eta(\theta)\Big)\overline{q(0)}=0. \end{gather*} We also have \begin{equation*} \Big(2i\omega_{0}\tau_{k}I-\int_{-1}^{0}\exp(i\omega_{0} \tau_{k}\theta)d\eta(\theta)\Big)E_{1}=2\tau_{k}\begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)\alpha_{1}-p_{1}\alpha_{2}\\ \beta \alpha_{1}-p_{2}\alpha_{1}\alpha_{2}\\ ( q_{1}p_{1}\alpha_{2}+q_{2}p_{2}\alpha_{1})\exp(-2i\omega_{0}\tau_{k}) \end{bmatrix},\end{equation*} which leads to \begin{align*} &\begin{bmatrix} a_{11}& S_{*}(\frac{r}{k}+\beta) & p_{1}S_{*}\\ -\beta I_{*} & a_{22} & p_{2}I_{*} \\ -q_{1}p_{1} Y_{*}\exp(-2i\omega_{0}\tau_{k})& -q_{2}p_{2} Y_{*} \exp(-2i\omega_{0}\tau_{k}) & a_{33} \end{bmatrix}E_{1} \\ &=2\begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)\alpha_{1}-p_{1}\alpha_{2}\\ \beta \alpha_{1}-p_{2}\alpha_{1}\alpha_{2}\\ ( q_{1}p_{1}\alpha_{2}+q_{2}p_{2}\alpha_{1})\exp(-2i\omega_{0}\tau_{k}) \end{bmatrix}, \end{align*} where \begin{gather*} a_{11}= 2i\omega_{0}-(r-\frac{2rS_{*}}{k} -(\frac{r}{k}+\beta)I_{*}-p_{1}Y_{*}),\quad a_{22}=2i\omega_{0}-\beta S_{*}+c+d_{2}+p_{2}Y_{*},\\ a_{33}=2i\omega_{0}+d_{3}-(q_{1}p_{1} S_{*}+q_{2}p_{2} I_{*}) \exp(-2i\omega_{0}\tau_{k}), \end{gather*} therefore, \begin{equation}\label{5.6} \begin{aligned} E_{1}&= 2\begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)\alpha_{1}-p_{1}\alpha_{2}\\ \beta \alpha_{1}-p_{2}\alpha_{1}\alpha_{2}\\ ( q_{1}p_{1}\alpha_{2}+q_{2}p_{2}\alpha_{1})\exp(-2i\omega_{0}\tau_{k}) \end{bmatrix}\\ &\quad\times \begin{bmatrix} a_{11}& S_{*}(\frac{r}{k}+\beta) & p_{1}S_{*}\\ -\beta I_{*} & a_{22} & p_{2}I_{*} \\ -q_{1}p_{1} Y_{*}\exp(-2i\omega_{0}\tau_{k})& -q_{2}p_{2} Y_{*}\exp(-2i\omega_{0}\tau_{k}) & a_{33} \end{bmatrix}^{-1}. \end{aligned} \end{equation} provided \[ \begin{bmatrix} a_{11}& S_{*}(\frac{r}{k}+\beta) & p_{1}S_{*}\\ -\beta I_{*} & a_{22} & p_{2}I_{*} \\ -q_{1}p_{1} Y_{*}\exp(-2i\omega_{0}\tau_{k})& -q_{2}p_{2} Y_{*}\exp(-2i\omega_{0}\tau_{k}) & a_{33} \end{bmatrix} \] is invertible. Now, $\int_{-1}^{0}d\eta(\theta){W}_{11}(\theta)=-H_{11}(0)$ and \begin{equation*} H_{11}(0)=-g_{11}q(0)-\overline{g}_{11}\overline{q}(0)+2\tau_{k}\begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)Re(\alpha_{1})-p_{1}Re(\alpha_{2})\\ \beta Re(\alpha_{1})-p_{2}Re(\alpha_{1}\alpha_{2})\\ ( q_{1}p_{1}Re(\alpha_{2})+q_{2}p_{2}Re(\alpha_{1}) \end{bmatrix}, \end{equation*} and hence this leads to the equation \[ \begin{bmatrix} b_{11} & S_{*}(\frac{r}{k}+\beta) & p_{1}S_{*}\\ -\beta I_{*} & b_{22} & p_{2}I_{*} \\ -q_{1}p_{1} Y_{*}& -q_{2}p_{2} Y_{*} & b_{33}\end{bmatrix} E_{2}=2\begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)Re(\alpha_{1})-p_{1}Re(\alpha_{2})\\ \beta Re(\alpha_{1})-p_{2}Re(\alpha_{1}\alpha_{2})\\ ( q_{1}p_{1}Re(\alpha_{2})+q_{2}p_{2}Re(\alpha_{1}) \end{bmatrix}, \] where \begin{gather*} b_{11}=(r-\frac{2rS_{*}}{k}-(\frac{r}{k}+\beta)I_{*}-p_{1}Y_{*}),\quad b_{22}=-\beta S_{*}+c+d_{2}+p_{2}Y_{*},\\ b_{33}=d_{3}-(q_{1}p_{1} S_{*}+q_{2}p_{2} I_{*}). \end{gather*} Then $E_{2}$ can be obtained as \[ E_{2}=2\begin{bmatrix} -\frac{r}{k}-(\frac{r}{k}+\beta)Re(\alpha_{1})-p_{1}Re(\alpha_{2})\\ \beta Re(\alpha_{1})-p_{2}Re(\alpha_{1}\alpha_{2})\\ ( q_{1}p_{1}Re(\alpha_{2})+q_{2}p_{2}Re(\alpha_{1}) \end{bmatrix} \begin{bmatrix} b_{11} & S_{*}(\frac{r}{k}+\beta) & p_{1}S_{*}\\ -\beta I_{*} & b_{22} & p_{2}I_{*} \\ -q_{1}p_{1} Y_{*}& -q_{2}p_{2} Y_{*} & b_{33} \end{bmatrix}^{-1}. \] provided \[ \begin{bmatrix} b_{11} & S_{*}(\frac{r}{k}+\beta) & p_{1}S_{*}\\ -\beta I_{*} & b_{22} & p_{2}I_{*} \\ -q_{1}p_{1} Y_{*}& -q_{2}p_{2} Y_{*} & b_{33} \end{bmatrix} \] is invertible. By putting values of $E_{1}$ and $E_{2}$ we can obtain $W_{20}(\theta)$ and $W_{11}(\theta)$ and hence $g_{21}$ is completely determined. Hence as stated in \cite{CC, GX}, we can obtain the following values: \begin{equation}\label{5.8} \begin{gathered} c_{1}(0)=\frac{i}{2 \omega_{0} \tau_{k}} (g_{11}g_{20}-2\mid g_{11}\mid^{2}-\frac{\mid g_{02}\mid ^{2}}{3})+\frac{g_{21}}{2},\\ \mu_{2}=-\frac{Re(c_{1}(0))}{Re (\lambda^{'}(\tau_{k}))},\\ \beta_{2}=2 Re (c_{1}(0)),\\ T_{2}=-\frac{1}{\omega _{0}\tau_{k}}[\operatorname{Im}({c_{1}(0)}) +\mu _{2} \operatorname{Im}(\lambda^{'}(\tau_{k}))], \end{gathered} \end{equation} which determines the direction and stability of the model with delay at the critical value $\tau_{k}$. Now, we state the following main theorem of this section due to \cite{GX, CC, ST} \begin{theorem} \label{tgm2} \begin{itemize} \item[(i)] The sign of $\mu_{2}$ determined the direction of Hopf bifurcation. If $\mu_{2}>0 (\mu_{2}<0)$ , then the Hopf bifurcation is supercritical (sub critical). \item[(ii)] The stability of bifurcated periodic solutions is determined by $\beta_{2}$. The periodic solutions are stable if $\beta_{2}<0$ and unstable if $\beta_{2}>0$. \item[(iii)] The period of bifurcated periodic solutions is determined by $T_{2}$. The period increases if $T_{2}>0$ and decreases if $T_{2}<0$. \end{itemize} \end{theorem} \section{Numerical example} In this section, we consider a numerical example and generate some numerical simulations to verify our theoretical calculations. As an example, we choose the set of parameters in Table \ref{table4}. The initial values are taken as $S(0)=0.9$, $I(0)=0.9$, $Y(0)=0.2$. \begin{table}[htbp] \caption{Parameter values}\label{table4} \begin{center} \renewcommand{\arraystretch}{1.3} \begin{tabular}{|c|c|c|} \hline \text{Parameter }& \text{Numerical Value} & \text{Source}\\ \hline $r$ & $1/2$& Hu and Li \cite{GX}\\ $k$ & $1$ & Hu and Li \cite{GX}\\ $\beta$ &$ 1 $& Hu and Li \cite{GX}\\ $p_{1}$ & $1/8$& Hu and Li \cite{GX}\\ $p_{2}$& $6$&Hu and Li \cite{GX}\\ $d_{2}$&$1/4$& assumed\\ $d_{3}$&$1/2$& assumed\\ $c$&$1/4$& Hu and Li \cite{GX}\\ $q_{1}$&$1/2$& assumed\\ $q_{2}$&$3/4$& assumed\\ \hline \end{tabular} \end{center} \end{table} It is observed that the system has the equilibrium points $(0,0,0)=E_{10}$, $(1 ,0,0)=E_{11}$, $(\frac{1}{2},\frac{1}{6},0)=E_{12}$ and $E_{5}=(\frac{31}{48},\frac{353}{1728},\frac{7}{6})=E_{*}$ and all are locally stable. The equilibrium point $(\frac{4}{3},0,\frac{-4}{3})=E_{13}$ does not exist. By calculation, it is observed that $\omega_{0}=0.3694$ and $\tau_{0}=1.5326$. Thus positive equilibrium $E_{*}$ is asymptotically stable when $0<\tau<\tau_{0}=1.5326$. The system undergoes a Hopf bifurcation when system crosses through $\tau_{0}$. Few numerical simulations are represented by using MATLAB. It is observed from numerical simulations that the phase diagram of the system changes with slight changes in the initial values. First two figures (\ref{fig1} and \ref{fig2}) are drawn for the value of $\tau=\tau_{0}=1.5326$ and $\tau=3.6249>\tau_{0}=1.5326$ respectively. From theoretical foundation it is observed that system is stable up to the critical value of time delay $\tau$. Critical value is calculated in \eqref{4.7}. Beyond the critical value of $\tau$, Hopf bifurcation occurs. It is also observed in simulations. It is easy to see that system undergoes Hopf bifurcation (See Figures \ref{fig1} and \ref{fig2}). The simulations are taken for the time interval $[0, 1000]$ . For the simulation interval $[0, 1000]$. From Figure \ref{fig1}, it is found that all the three populations are unstable with respect to time. Similar explanation is concluded from Figure \ref{fig3}. Figure \ref{fig3} shows the stable nature of positive equilibrium which is drawn for $\tau=1.1026<1.5326=\tau_{0}$, which is again consistent with the theoretical formulation. Therefore, numerical simulations are consistent with the theoretical formulation. If we choose the different set of values of the parameters $q_{1}$ and $q_{2}$ with the same other parameters and initial values, we may have different figures (see Remark \ref{rmk2}, \ref{rmk3}). \begin{figure}[htb] \begin{center} \includegraphics[width=0.45\textwidth]{fig1} \includegraphics[width=0.45\textwidth]{fig2}\\ \includegraphics[width=0.45\textwidth]{fig3} \includegraphics[width=0.45\textwidth]{fig4}\\ \includegraphics[width=0.45\textwidth]{fig5} \includegraphics[width=0.45\textwidth]{fig6}\\ \includegraphics[width=0.45\textwidth]{fig7} \includegraphics[width=0.45\textwidth]{fig8} \end{center} \caption{Phase portrait of model with parameter values in Table \ref{table4} with initial values $S(0)=0.9$, $I(0)=0.9$, $Y(0)=0.2$} \label{fig1} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.45\textwidth]{fig9} \includegraphics[width=0.45\textwidth]{fig10}\\ \includegraphics[width=0.45\textwidth]{fig11} \includegraphics[width=0.45\textwidth]{fig12}\\ \includegraphics[width=0.45\textwidth]{fig13} \includegraphics[width=0.45\textwidth]{fig14}\\ \includegraphics[width=0.45\textwidth]{fig15} \includegraphics[width=0.45\textwidth]{fig16} \end{center} \caption{Phase portrait of model with parameter values in Table \ref{table4} with initial values $S(0)=0.9$, $I(0)=0.9$, $Y(0)=0.2$} \label{fig2} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.45\textwidth]{fig17} \includegraphics[width=0.45\textwidth]{fig18}\\ \includegraphics[width=0.45\textwidth]{fig19} \includegraphics[width=0.45\textwidth]{fig20}\\ \includegraphics[width=0.45\textwidth]{fig21} \includegraphics[width=0.45\textwidth]{fig22}\\ \includegraphics[width=0.45\textwidth]{fig23} \includegraphics[width=0.45\textwidth]{fig24} \end{center} \caption{Phase portrait of model with parameter values in Table \ref{table4} with initial values $S(0)=0.9$, $I(0)=0.9$, $Y(0)=0.2$} \label{fig3} \end{figure} \section{Discussion and future directions} In the present study, we have considered an eco-epidemiological model with disease in the prey population only. We have considered both cases with delay and without delay for analysis purpose. Local stability of all equilibrium points has been discussed. We have found the time delay as a game changer. This has been observed that time delay $\tau$ may change the stability and even bifurcation may occur. Hopf-bifurcation analysis is presented, by the application of famous normal form theory, Riesz representation theorem and central limit theorem, stability and direction of bifurcated periodic solutions have been investigated. Our analytical results has been compared with those in \cite{GX}, when $d_{2}=0$ and $q_{1}=q_{2}=q$. Stability of the non zero (positive) equilibrium $E_{5}$ indicates the existence and survival of all the three species in the ecosystem. Ecologically this equilibrium is very important, because it provides actual interaction among all the three living components of the ecosystem. Actual balance is maintained under this situation. Ecologists are interested to observe the stability of non zero equilibrium point. It is also remarkable that we classify our analysis in two parts (i) without delay (ii) with delay. Ecologically, corresponding equilibrium points listed for our models \eqref{2.1} and \eqref{3.1} are same. For example, ecologically, $E_{5}$ and $E_{*}$ are same but mathematically both are different. Theoretically, we can consider few examples of different ecosystems. If we consider Desert, Tundra, Savana geographical areas and consider an ecosystem from this area. Then positive equilibrium is more important, because in these areas many predators are going to decline or decay because their prey are also facing natural problems like climate changes, lake of water, low quality of oxygen etc. and they also decaying, therefore in such geographical areas the predator population is not in position to catch their prey as they are limited and therefore they remain hungry for most of the times, finally predator population starts decaying. For the ecological imbalance in such ecosystems, the main reason is the change in climate. Besides cutting of trees and global warming are also important environmental issues. Thus positive equilibrium can never exists, however if it exists, then not stable. If we consider hilly areas and consider an ecosystem from this area, then the cutting of trees and farming are two issues which are responsible for disturbance of climate in this area. In near past, one natural disaster in Uttarakhand (India) occurred, possibly because of cutting of trees and other developmental projects. This changed the climate in hilly areas as well. This occurs the ecological imbalance in these areas and disturbs the natural ecosystems. Prey and predator both decays simultaneously due to unexpected mortality in the system. Thus if the positive equilibrium exists and is stable then the stability losses due to disaster. Similarly, if we choose the aquatic ecosystem of any river in India, then the water pollution is a major factor in this ecosystem. Pollution disturbs the stability of the ecosystem and disturbs the ecological balance. Positive equilibrium may loose stability in such ecosystems as well. If we consider an ecosystem from green forest where pollution, weather conditions etc. are favorable for living species then actual prey-predation interaction will be occurred. In this situation, if the positive equilibrium exists and is stable this will provide the coexistence of all the three species. Thus equilibrium points $E_{5}$ and $E_{*}$, if exists, are stable for an ecosystem from a forest area. The parameters considered for numerical computation are purely imaginary and are suitably selected and are quite similar to result discussed in \cite{GX}. An attempt may be done to estimate real parameters and these parameters may fit with the present mathematical model. This task may be achieved by real data collection from field. Further, to study the model more scientifically, control strategies may also be investigated. This study may have applications with real ecosystems and biomass available in the real world. \section{Appendix} \subsection*{Appendix 1: $C_{i}$, $i=1,2,3$ of \eqref{3.7}} \[ C_{1}=-\Big((r-d_{2}-d_{3}-c)+\widetilde{S}[-\frac{2r}{k}+\beta + q_{1}p_{1}]+\widetilde{I}[-\frac{r}{k}+q_{2}p_{2}-\beta] -\widetilde{Y}(p_{1}+p_{2})\Big), \] \begin{align*} C_{2} &=\widetilde{S}^{2}[\beta q_{1}p_{1}-\frac{2r q_{1}p_{1}}{k} -\frac{2r \beta}{k}]+\widetilde{Y}^{2}[p_{1}p_{2}] +\widetilde{I}^{2}[(\frac{r}{k}+\beta)q_{2}p_{2}]\\ &\quad +\widetilde{S}\widetilde{I}[\beta q_{2}p_{2}-\frac{2r q_{2}p_{2}}{k} -(\frac{r}{k}+\beta)q_{1}p_{1}]+\widetilde{S}\widetilde{Y} [-q_{1}p_{1}p_{2}-q_{1}p_{1}^{2}+\frac{2r p_{2}}{k}-p_{1}\beta]\\ &\quad +\widetilde{I}\widetilde{Y}[-q_{2}p_{1}p_{2}+(\frac{r}{k}+\beta)p_{2}] +\widetilde{S}\Big[-\beta d_{3}-(c+d_{2})q_{1}p_{1}+rq_{1}p_{1} +\frac{2rd_{3}}{k}\\ &\quad +r\beta+\frac{2r(c+d_{2})}{k}\Big]\\ &\quad +\widetilde{I}\Big[-(c+d_{2})q_{2}p_{2}+rq_{2}p_{2} +(\frac{r}{k}+\beta)d_{3}+(\frac{r}{k}+\beta)(c+d_{2})\Big] \\ &\quad +\widetilde{Y}\Big[d_{3}p_{1}-rp_{2}+p_{1}(c+d_{2})\Big] +\Big[d_{3}(c+d_{2})-rd_{3}-r(c+d_{2})\Big], \end{align*} \begin{align*} C_{3}&=-\Big(\widetilde{S}^{3}[-\frac{2r q_{1}p_{1}\beta}{k}] +\widetilde{S}^{2}\widetilde{Y}[\frac{2r q_{1}p_{1}p_{2}}{k}] +\widetilde{S}^{2}\widetilde{I}[\frac{2r q_{1}p_{1}p_{2}}{k}] +\widetilde{I}^{2}\widetilde{Y}[(\frac{r}{k}+\beta)q_{2}p_{2}^{2}]\\ &\quad +\widetilde{Y}^{2}\widetilde{I}[q_{2}p_{2}^{2}p_{1}] +\widetilde{S}\widetilde{I}\widetilde{Y} \Big[-\frac{2r q_{2}p_{2}^{2}\beta}{k}+2(\frac{r}{k} +\beta)q_{1}p_{1}p_{2}-2\beta q_{2}p_{1}p_{2}\Big]\\ &\quad +\widetilde{S}^{2}\Big[ \beta rq_{1}p_{1} +\frac{2r d_{3}\beta}{k}+\frac{2r q_{1}p_{1}(c+d_{2})}{k}\Big]\\ &\quad +\widetilde{I}^{2}[(\frac{r}{k}+\beta)q_{2}p_{2} (c+d_{2})]+\widetilde{Y}^{2}[d_{3}p_{1}p_{2}] \\ &\quad +\widetilde{S}\widetilde{I}\Big[ \beta rq_{2}p_{2} +\frac{2r q_{2}p_{2}(c+d_{2})}{k}+(\frac{r}{k}+\beta) (c+d_{2})q_{1}p_{1}\Big]\\ &\quad +\widetilde{S}\widetilde{Y}[-rq_{1}p_{1}p_{2} +\frac{2rp_{2}d_{3}}{k}+d_{3}\beta p_{1}]\\ &\quad +\widetilde{I}\widetilde{Y}\Big[-rq_{2}p_{2}^{2} +(\frac{r}{k}+\beta)d_{3}p_{2}+(c+d_{2})q_{2}p_{1}p_{2}\Big]\\ &\quad +\widetilde{S}\Big[-d_{3}r\beta-r(c+d_{2})q_{1}p_{1} -\frac{2r(c+d_{2})d_{3}}{k}\Big] \\ &\quad +\widetilde{I}\Big[-r(c+d_{2})q_{2}p_{2}-(\frac{r}{k}+\beta)(c+d_{2})d_{3}\Big]\\ &\quad +\widetilde{Y}\Big[-d_{3}p_{2}r-(c+d_{2})d_{3}p_{1}\Big]+r(c+d_{2})d_{3}\Big). \end{align*} \subsection*{Appendix 2: $m_{i}, n_{j}$, $i=0,1,2$; $j=0,1,2$ of \eqref{4.3}} \[ m_{2}=\Big(\frac{2r}{k}+\frac{p_{1}\beta}{p_{2}} -\frac{q_{1} p_{1}\beta}{q_{2} p_{2}}S_{*}\Big) +\Big(\frac{rd_{3}}{q_{2} p_{2}k}+d_{3}-\frac{p_{1}(c+d_{2})}{p_{2}}-r\Big), \] \begin{align*} m_{1}&=\Big(\{d_{3}(c+d_{2})-r(c+d_{2})\}+S_{*}\{-d_{3}\beta+\frac{2r}{k}d_{3} +r\beta+\frac{2r}{k}(c+d_{2})\} \\ &\quad +I_{*}\{(\frac{r}{k}+\beta)(d_{3}+c+d_{2})\}+Y_{*}\{d_{3}p_{2}+d_{3}p_{1} +(c+d_{2})p_{1}+(r)p_{2}\} \\ &\quad +S_{*}^{2}(-\frac{2r\beta}{k})+Y_{*}^{2}(p_{1}p_{2}) +S_{*}Y_{*}(-\beta p_{1}+\frac{2r}{K}p_{2})+I_{*}Y_{*}(\frac{r}{k}+\beta)p_{2}\Big), \end{align*} \begin{align*} m_{0}&=\Big(\{-r(c+d_{2})d_{3}+S_{*}\{d_{3}\beta r+\frac{2r}{k}d_{3}\}(c+d_{2})\} +I_{*}\{d_{3}(c+d_{2})(\frac{r}{k}+\beta)\}\\ &\quad +Y_{*}\{d_{3}(c+d_{2})p_{1}-rd_{3}p_{2}\} +S_{*}^{2}(-\frac{2r}{k}d_{3}\beta)+Y_{*}^{2}(d_{3}p_{2}p_{1})\\ &\quad +S_{*}Y_{*}\{-d_{3}\beta p_{1}-\frac{2r}{k}d_{3}p_{2}\} +I_{*}Y_{*}\{(\frac{r}{k}+\beta)p_{2}d_{3}\}\Big), \end{align*} $n_{2}=-d_{3}$, \begin{align*} n_{1}&=\Big(S_{*}^{2}\{q_{1}p_{1}\beta-\frac{-2r}{k}q_{1}p_{1}\} +I_{*}^{2}\{-(\frac{r}{k}+\beta)q_{2}p_{2}\}\\ &\quad +S_{*}Y_{*}\{-q_{1}p_{2}p_{1}\} +S_{*}I_{*}\{q_{2}p_{2}\beta+\frac{2r}{k}q_{2}p_{2}\} +Y_{*}I_{*}\{-q_{2}p_{2}p_{1}\} \\ &\quad +S_{*}\{\{-(c+d_{2})-r\}q_{1}p_{1}\} +I_{*}\{-(c+d_{2})q_{2}p_{2}+rq_{2}p_{2}\}\Big), \end{align*} \begin{align*} n_{0}&=-\Big(S_{*}^{2}\{q_{1} p_{1}\beta(r-\frac{2r}{k})+\frac{2r}{k} (c+d_{2})q_{1} p_{1}-q_{1} p_{1}^{2}\beta\} +I_{*}^{2} \{(\frac{r}{k}+\beta)(c+d_{2})q_{2} p_{2}\\ &\quad +S_{*}^{2}Y_{*} \{\frac{2r}{k}q_{1} p_{1}p_{2}+q_{1} p_{1}^{2}\beta\} +S_{*}^{2}I_{*}\{-\frac{2r}{k}q_{2} p_{2}\beta\} +S_{*}I_{*}Y_{*}\{(\frac{r}{k}+\beta)q_{1} p_{1}p_{2} \\ &\quad -2\beta q_{2} p_{1}p_{2}\} +S_{*}Y_{*}\{\{-r+(\frac{r}{k}+\beta)\}q_{1} p_{1} p_{2}\} \\ &\quad +S_{*}I_{*}\{\beta rq_{2} p_{2}+\frac{2r}{k}(c+d_{2})q_{2} p_{2}+(c+d_{2}) (\frac{r}{k}+\beta)q_{1} p_{1}\}\\ &\quad +Y_{*}I_{*}\{(c+d_{2})q_{2} p_{2} p_{1}\} +S_{*}\{(c+d_{2})rq_{1} p_{1}\}+I_{*}\{-(c+d_{2})rq_{2} p_{2}\}\Big). \end{align*} \subsection*{Appendix 3} We have \[ g(z,\overline{z})=g_{20}(\theta)\frac{z^{2}}{2}+g_{02}(\theta) \frac{\overline{z}^{2}}{2}+g_{11}(\theta)z\overline{z}+g_{21} \frac{\overline{z}z^{2}}{\lfloor{3}}+\dots\,, \] \begin{align*} g(z,\overline{z})&=(\overline{q}^{*})^{T}F(z,\overline{z}) \\ &=\tau_{k}\overline{B}(1,\alpha_{1}^{*},\alpha_{2}^{*}) \begin{pmatrix} -\frac{r}{k}u_{1}^{2}(t)-(\frac{r}{k}+\beta)u_{1}(t)u_{2}(t)-p_{1}u_{1}(t)u_{3}(t)\\ \beta u_{1}(t)u_{2}(t)- p_{2}u_{2}(t)u_{3}(t)\\ p_{1}q_{1}u_{1}(t-1)(t)u_{3}(t-1)+p_{2}q_{2}u_{1}(t-1)u_{2}(t-1), \end{pmatrix} \end{align*} further, it is noticed that \begin{gather*} u(t+\theta)=W(t,\theta)+z(t)q(\theta)+\overline{z}(t)\overline{q}(\theta), \\ u_{1}(t)=z+\overline{z}+W^{(1)}(t,0),\\ u_{2}(t)=\alpha_{1}z+\overline{\alpha_{1}}\overline{z}+W^{(2)}(t,0),\\ u_{3}(t)=\alpha_{2}z+\overline{\alpha_{2}}\overline{z}+W^{(3)}(t,0),\\ u_{1}(t-1)=z \exp(-i\omega_{0}\tau_{k})+\overline{z}\exp(i\omega_{0}\tau_{k})+W^{(1)}(t,-1),\\ u_{2}(t-1)=\alpha_{1} z \exp(-i\omega_{0}\tau_{k})+\overline{\alpha_{1}}\overline{z}\exp(i\omega_{0}\tau_{k})+W^{(2)}(t,-1),\\ u_{3}(t-1)=\alpha_{2} z \exp(-i\omega_{0}\tau_{k})+\overline{\alpha_{2}}\overline{z}\exp(i\omega_{0}\tau_{k})+W^{(3)}(t,-1), \end{gather*} hence \begin{align*} g(z,\overline{z}) &=\tau_{k}\overline{B}[-\frac{r}{k}u_{1}^{2}(t) -(\frac{r}{k}+\beta)u_{1}(t)u_{2}(t)-p_{1}u_{1}(t)u_{3}(t) \\ &\quad +\overline{\alpha_{1}}^{*}\{\beta u_{1}(t)u_{2}(t)- p_{2}u_{2}(t)u_{3}(t)\}\\ &\quad +\overline{\alpha_{2}}^{*}\{p_{1}q_{1}u_{1}(t-1)(t)u_{3}(t-1) +p_{2}q_{2}u_{1}(t-1)u_{2}(t-1)\}], \end{align*} putting the values of $u_{1}$, $u_{2}$, $u_{3}$, $u_{1}(t-1)$, $u_{2}(t-1)$, $ u_{3}(t-1)$ etc. in $g(z,\overline{z})$, we obtain \begin{align*} g(z,\overline{z})&=\tau_{k}\overline{B} \Big(-\frac{r}{k}[z+\overline{z}+W^{(1)}(t,0)]^{2} -(\frac{r}{k}+\beta)[z+\overline{z}+W^{(1)}(t,0)] [\alpha_{1}z\\ &\quad +\overline{\alpha_{1}}\overline{z}+W^{(2)}(t,0)] -p_{1}[z+\overline{z}+W^{(1)}(t,0)][\alpha_{2}z +\overline{\alpha_{2}}\overline{z}+W^{(3)}(t,0)] \\ &\quad +\overline{\alpha_{1}}^{*}(\beta[z+\overline{z}+W^{(1)}(t,0)] [\alpha_{1}z+\overline{\alpha_{1}}\overline{z}+W^{(2)}(t,0)]\\ &\quad -p_{2} [\alpha_{1}z+\overline{\alpha_{1}}\overline{z}+W^{(2)}(t,0)] [\alpha_{2}z+\overline{\alpha_{2}}\overline{z}+W^{(3)}(t,0)])\\ &\quad +\overline{\alpha_{2}}^{*}(p_{1}q_{1}[z \exp(-i\omega_{0}\tau_{k}) +\overline{z}\exp(i\omega_{0}\tau_{k})+W^{(1)}(t,-1)] \\ &\quad\times [\alpha_{2} z \exp(-i\omega_{0}\tau_{k}) +\overline{\alpha_{2}}\overline{z} \exp(i\omega_{0}\tau_{k})+W^{(3)}(t,-1)]\\ &\quad +p_{2}q_{2}[z \exp(-i\omega_{0}\tau_{k}) +\overline{z}\exp(i\omega_{0}\tau_{k})+W^{(1)}(t,-1)] \\ &\quad \times[\alpha_{1} z \exp(-i\omega_{0}\tau_{k})+\overline{\alpha_{1}}\overline{z} \exp(i\omega_{0}\tau_{k})+W^{(2)}(t,-1)])\Big), \end{align*} from this equation we can find the values of the coefficients $g_{20}(\theta)$, $g_{02}(\theta)$, $g_{11}(\theta)$, $g_{21}(\theta)$ etc. by comparing the same powers of $z$, we have \begin{align*} g_{20}&=2 \tau_{k}\overline{B}\{-\frac{r}{k}-(\frac{r}{k}+\beta)\alpha_{1} -p_{1}\alpha_{2}+\beta\overline{\alpha_{1}}^{*}\alpha_{1}- \overline{\alpha_{1}}^{*}\alpha_{1}\alpha_{2}p_{2} \\ &\quad +\overline{\alpha_{2}}^{*}(p_{1}q_{1}\alpha_{2} +p_{2}q_{2}\alpha_{1})\exp(-2 i\omega_{0}\tau_{k})\}, \end{align*} \begin{align*} g_{11}&=\tau_{k}\overline{B}(-\frac{2r}{k} +(\overline{\alpha_{1}}^{*})\beta+\overline{\alpha_{2}}^{*}p_{2}q_{2} -\frac{r}{k}+\beta)(\overline{\alpha_{1}}+\alpha_{1})\\ &\quad +(\overline{\alpha_{2}}^{*}p_{1}q_{1}-p_{1})(\overline{\alpha_{2}} +\alpha_{2}) - \overline{\alpha_{1}}^{*}p_{2}(\overline{\alpha_{2}}\alpha_{1} +\overline{\alpha_{1}}\alpha_{2})), \end{align*} \begin{align*} g_{02}&=2\tau_{k}\overline{B}\{-\frac{r}{k}-(\frac{r}{k}+\beta) \overline{\alpha_{1}}-p_{1}\overline{\alpha_{2}} +\beta\overline{\alpha_{1}}^{*}\overline{\alpha_{1}}- \overline{\alpha_{1}}^{*}\overline{\alpha_{1}}\overline{\alpha_{2}}p_{2}\\ &\quad +\overline{\alpha_{2}}^{*}(p_{1}q_{1} \overline{\alpha_{2}}+p_{2}q_{2}\overline{\alpha_{1}}) \exp(2 i\omega_{0}\tau_{k})\}, \end{align*} \begin{align*} g_{21}&=2\tau_{k}\overline{B}\Big(-\frac{r}{k}(2W_{11}^{(1)}(0)+W_{20}^{(1)}(0)) -(\frac{r}{k}+\beta)[W_{11}^{(2)}(0)+\alpha_{1}W_{11}^{(1)}(0) \\ &\quad +\frac{1}{2}\overline{\alpha_{1}}W_{20}^{(1)}(0)+\frac{1}{2}W_{20}^{(1)}(0)] -p_{1}\big[W_{11}^{(3)}(0)+\frac{1}{2}W_{20}^{(3)}(0) +\alpha_{2}W_{11}^{(1)}(0)\\ &\quad +\frac{1}{2}W_{20}^{(1)}(0)\big] +\overline{\alpha_{1}}^{*}\beta[W_{11}^{(2)}(0)+ \frac{1}{2}W_{20}^{(2)}(0)+\alpha_{1}W_{11}^{(1)}(0) +\frac{1}{2}\overline{\alpha_{1}}W_{20}^{(1)}(0)] \\ &\quad -\overline{\alpha_{1}}^{*}p_{2}[\alpha_{1}W_{11}^{(3)}(0) +\alpha_{2}W_{11}^{(2)}(0)+\frac{1}{2}\overline{\alpha_{2}}W_{20}^{(2)}(0) +\frac{1}{2}\overline{\alpha_{1}}W_{20}^{(3)}(0)]\\ &\quad +\overline{\alpha_{2}}^{*}p_{1}q_{1}[W_{11}^{(3)}{(-1)} e^{-i\omega_{0}\tau_{k}} +\alpha_{2}W_{11}^{(1)}{(-1)}e^{-i\omega_{0}\tau_{k}} +\frac{1}{2}W_{20}^{(3)}(-1)e^{i\omega_{0}\tau_{k}}\\ &\quad + \frac{1}{2}W_{20}^{(1)}(-1)e^{i\omega_{0}\tau_{k}}] +\overline{\alpha_{2}}^{*}p_{2}q_{2}[W_{11}^{(2)}(-1)e^{-i\omega_{0}\tau_{k}} +\alpha_{1}W_{11}^{(-1)}e^{-i\omega_{0}\tau_{k}}\\ &\quad +\frac{1}{2}W_{20}^{(2)}(-1)e^{i\omega_{0}\tau_{k}}+\frac{1}{2}W_{20}^{(1)} (-1)\overline{\alpha_{1}}e^{i\omega_{0}\tau_{k}}]\Big). \end{align*} \begin{thebibliography}{99} \bibitem{KB} Berthier, K.; Langlais, M.; Auger, P.; Pontier, D.; Dynamics of a feline virus with two transmission medels with exponentially growing host populations. \newblock \emph{Proc. Roy. Soc. Lond.} 3267 (2000): 2049--2056 {http: //www.ncbi.nlm.nih.gov/pubmed/11416908} \bibitem{JO} Chattopadhyay, J.; Arino, O.; A predator-prey model with disease in prey. \newblock \emph{Non Linear Analysis} 36 (1999): 747--766 {http: //repository.ias.ac.in/7560/1/411.pdf} \bibitem{MH} Haque, M; Zhen, J; Ventrino, E; An eco-epidemiological predator-prey model with standard disease incidence. \newblock \emph{Mathematical Methods in the Applied Sciences} http: //dx.doi.org /10.1002/mma.1071, 2008. \bibitem{FK} Hilker, F. M.; Schmitz, K.; Disease-induced stabilization of predator prey oscillations. \newblock \emph{Journal of Theoretical Biology } 255 (2008): 299--306 {http: //dx.doi.org/10.1016/j.jtbi.2008.08.018} \bibitem{GX} Hu, G.-P.; Li, X.-L.; Stability and Hopf bifurcation for a delayed predator-prey model with disease in the prey. \newblock \emph{Chaos, Solitons and Fractals } 45 (2012): 229--237 {http: //dx.doi.org/10.1016 /j.chaos.2011.11.011} \bibitem{ST} Jana, S.; Kar, T. K.; Modeling and Analysis of a prey-predator system with disease in the prey. \newblock \emph{Chaos, Solitons and Fractals} 47 (2013): 42--53 {http: //dx.doi.org/10.1016 /j.chaos.2012.12.002} \bibitem{MC} De Jong, M. C. M.; et al; How does infection depend on the population size? in D. Mollison(Ed.), Epidemic Models, their structure and Relation in data. \newblock \emph{Cambridge University Press}, 1994. {http: //igitur-archive.library.uu.nl/vet} \bibitem{DG} Greenhalgh, D.; Haque, M.; A predator-prey model with disease in the prey species only. \newblock \emph{Math. Methods Appl. Sci.} 30 (2007): 911--929. {http: //dx.doi.org/10.1002/mma.815} \bibitem{WS} Kermack, W. O.; Mckendrick, A. G.; Contributons to the mathematical Theory of Epidimics, Part 1. \newblock \emph{Proc. Roy. Soc. A} 115 (1927): 700--721 {http: //rspa.royalsocietypublishing.org /content/115/772/700.full.pdf+html} \bibitem{ME} Haque, M.; Venturino, E.; Increase of the prey may decrease the healthy prtedator population in presence of a disease in predator. \newblock \emph{HERMIS} 7 (2006): 38-59 {http: //www.aueb.gr /pympe/hermis/hermis-volume-7} \bibitem{LZ} Litao, H.; Zhien, M. A.; Four Predator Prey Models with Infectious Diseases. \newblock \emph{Mathematical and Computer Modeling} 34 (2001): 849--858 {http: //homepage.math.uiowa.edu/~hethcote/PDFs/2001MathCompMod.pdf} \bibitem{CC} Karaaslanli, C. C.; Bifurcation Analysis and Its Applications, Chapter 1 in: Mykhaylo Andriychuk (Ed.),Numerical Simulation: From Theory to Industry. (2012) 1-34 {http: //dx.doi.org/10.5772/50075} \bibitem{AL} Lotka, A.; Elements of Physical Biology. \newblock \emph{Williams and Wilkins, Baltimore}. \bibitem{ML} Maoxing, L.; Zhen, J.; Haque, M.; An Impulsive predator-prey model with communicable disease in the prey species only. \newblock \emph{Non Linear analysis: Real World Applications} 10 (2009): 3098--3111 {http: //dx.doi.org/10.1016/j.nonrwa.2008.10.010} \bibitem{BB} Mukhopadhyaya, B.; Bhattacharyya, R.; Dynamics of a delay-diffusion prey-predator Model with disease in the prey. \newblock \emph{J. Appl. Math and Computing } 17 (2005): 361--377 {http: //dx.doi.org/10.1007/BF02936062} \bibitem{RM} Naji, R. K.; Mustafa, A. N.; The Dynamics of an Eco-Epidemiological Model with Nonlinear Incidence Rate. \newblock \emph{Journal of Applied Mathematics} {http: //dx.doi.org/10.1155/2012/852631}, 2012. \bibitem{S} Samanta, G. P.; Analysis of a delay nonautonomous predator-prey system with disease in the prey. \newblock \emph{Non Linear Analysis: Modeling and Control} 15 (2010): 97--108. \bibitem{SJB} Sophia Jang, R.-J.; Baglama, J.; Continous-time predator-prey models with parasites. \newblock \emph{J. Bio. Dyn} 3 (2009): 87--98 {http: //dx.doi.org/10.1080/1751375080228325} \bibitem{VV} Volterra, V.; Variazioni e fluttauazionidel numero d individui in species animals conviventii. \newblock \emph{Mem Acd. Linciei} (1926) 231--33. \bibitem{FW} Wang, F.; Kuang, Y.; Ding, C.; Zhang, S.; Stability and bifurcation of a stage-structured predator-prey model with both discete and distributed delays, \newblock \emph{Chaos, Solitons and Fractals } 46 ( 2013): 19--27 {http: //dx.doi.org/10.1016/j.chaos.2012.10.003} \bibitem{YX} Yanni, X.; Chen, L.; Modeling and analysis of a predator-prey model with disease in the prey. \newblock\emph{Mathematical Biosciences} 171 (2001): 59--82 {http: /dx.doi.org/10.1016/S0025-5564(01)00049-9} \bibitem{YK} Yang, K.; Delay differential equations with applications in population dynamics. \newblock\emph{Academic Press, INC.}, 1993. \end{thebibliography} \end{document}