A delayed eco-epidemiological model with nonlinear incidence rate and Crowley–Martin functional response for infected prey and predator

Abstract

A prey–predator model with disease in prey population is considered here. The prey population is classified into two categories: susceptible prey and infected prey. The nonlinear incidence rate, at which the susceptible prey becomes infected, is incorporated. Here, predator consumes infected prey only according to Crowely–Martin functional response. Time lag is introduced to convert the susceptible prey to infected one, and the present model becomes delayed eco-epidemic model. The positivity, boundedness and permanence of the solutions of both models are guaranteed here. The existence of equilibrium points is assured. Local as well as global stability analyses of both the models around the equilibrium points are performed and finally summarized in a table, comparing the conditions between two models. The occurrence of bifurcation (saddle-node, transcritical, Hopf and Bogdanov–Takens) in the non-delay model is explored here. The nature and direction of Hopf bifurcation are searched using normal form and centre manifold theorems for delay model. Lastly, the numerical simulation is carried out to verify the analytical findings. An effort is made to search the parameters, by paying special attention to newly introduced parameters, which control the dynamics of both models. Accordingly, stable/unstable zones on $x\tau _\mathrm{c}$-planes, where x is any considered parameter and $\tau _\mathrm{c}$ is the value of critical delay time, are proposed. The stranger attractor is searched in the delay model. Bifurcated diagrams of population with respect to newly introduced parameters are exemplified and discussed. When the population (S, I and P) starts oscillating, their time-averaged values are presented to show their variation due to aforesaid parameters. Moreover, pattern of all the curves $S=S(x)$, $I=I(x)$ and $P=P(x)$, where x is any aforesaid parameters, are certified using their respective gradient.

Appendices

Appendix A

Consider the following positive definite function about the equilibrium point $E^*$

$$\begin{aligned} V(S,I,P)&=\left[ S-S^*-S^*\ln \frac{S}{S^*}\right] \\&\quad +C_1\left[ I-I^*-\ln \frac{I}{I^*}\right] \\&\quad +C_2\left[ P-P^*-P^*\ln \frac{P}{P^*}\right] . \end{aligned}$$

where $C_1$ and $C_2$ are constants. Differentiating V w.r.t. time t along the solution of the system and while choosing $C_1= \frac{1}{1+\alpha I^*}$ and $C_2=\frac{m(1+aI^*)}{m_1(1+bP^*)(1+\alpha I^*)}$, then we have

$$\begin{aligned} \frac{dV}{\mathrm{d}t}= & {} -\left[ b_{11}(S-S^*)^2+b_{22}(I-I^*)^2\right. \\&\left. +\,b_{33}(P-P^*)^2+b_{12}(S-S^*)(I-I^*)\right] . \end{aligned}$$

where $b_{11}=\frac{r}{k},\ b_{22}=C_1\left[ \frac{\alpha \beta S^*}{(1+\alpha I)(1+\alpha I^*)}- \frac{amP^*}{(1+aI)(1+aI^*)(1+bp^*)}\right] ,\ b_{33}= \left[ \frac{m_1bI^*}{(1+aI^*)(1+bP)(1+bP^*)}+e\right] C_2$ and $\ b_{12}=\frac{r}{k}$. Clearly, $\frac{dV}{\mathrm{d}t}=0$ if and only if $S=S^*$, $I=I^*$ and $P=P^*$. The sufficient condition for $\frac{dV}{\mathrm{d}t}$ to be negative definite around $E_3(S^*,I^*,P^*)$ is that $\ b_{22}>0$ and $\ 4b_{11}b_{22}>b_{12}^2$. These conditions are satisfied under condition of Theorem 4.5. Therefore, $E_3(S^*,I^*,P^*)$ is globally asymptotically stable under condition of Theorem 4.5.

Appendix B

$$\begin{aligned} {\bar{L}}_1&=\left[ \frac{2\alpha \beta ^2 K^*I^*}{r(1+\alpha I^*)^3}+2am{\bar{L}}_4\left\{ P^*+\frac{m_1I^*}{{\bar{L}}_3} \right\} \right. \\&\quad +\frac{2bmm_1^2I^*{\bar{L}}_4}{{\bar{L}}_3^2}\left\{ 1+\frac{m_1I^*}{P^*{{\bar{L}}_3}^3}\right\} \\&\quad \left. + \frac{2emm_1^2I^*(1+bP^*)^2}{P^*{\bar{L}}_3^3(1+aI^*)^2} \right] ,\\ {\bar{L}}_2&=\left[ \frac{2\beta K^*}{r(1+\alpha I^*)^2}\left\{ \frac{r}{K^*}+\frac{\beta }{(1+\alpha I^*)^2} \right\} \right. \\&\quad \left. +\frac{2\alpha \beta S^*}{(1+\alpha I^*)^3}+\frac{2mm_1{\bar{L}}_4}{{\bar{L}}_3}\left\{ 1+\frac{m_1 I^*}{P^*{\bar{L}}_3} \right\} \right] ,\\ {\bar{L}}_3&=m_1bI^*+e(1+aI^*)(1+bP^*)^2 \ \text {and} \\ {\bar{L}}_4&=\left\{ (1+aI^*)^3(1+bP^*) \right\} ^{-1}. \end{aligned}$$

Appendix C

The characteristic equation corresponding to $E_2$ is

$$\begin{aligned}&\left[ \lambda -\left( \frac{m_1\hat{I}}{1+a\hat{I}}-d \right) \right] \left[ \left\{ {\lambda }^2 +\lambda \left( c+\frac{r\hat{s}}{k}\right) \right. \right. \\&\quad \left. \left. +\frac{rc}{k}\left( \hat{S}+\hat{I}\right) \right\} + e^{-\lambda \tau }\left\{ -\frac{c^2}{\beta \hat{S}}\lambda \right. \right. \\&\quad \left. \left. +\left( \frac{c^3\hat{I}}{\beta \hat{S}^2}-\frac{c^2r}{\beta k}\right) \right\} \right] =0. \end{aligned}$$

Clearly, one eigenvalue is $\frac{m_1\hat{I}}{1+a\hat{I}}-d<0$ if $\hat{I}<\frac{d}{m_1-ad}$ (under this condition non-delay system was asymptotically stable) and other two eigenvalues are the root of the equation

$$\begin{aligned} \lambda ^2+l_1\lambda +l_2=e^{-\lambda \tau }(l_3\lambda +l_4) \end{aligned}$$

(6.1)

where $l_1=c+\frac{r\hat{s}}{k}$, $l_2=\frac{rc}{k}\left( \hat{S}+\hat{I}\right)$, $l_3=-\frac{c^2}{\beta \hat{S}}$ and $l_4=\frac{c^3\hat{I}}{\beta \hat{S}^2}-\frac{c^2r}{\beta k}$.

Now we shall look for the condition under which the system is locally asymptotically stable and undergoes a Hopf bifurcation near $E_2$. Putting $\lambda =\mu +i\omega$ with $\mu =0$, we have

$$\begin{aligned} l_3\omega \sin {\omega \tau }+l_4\cos {\omega \tau }+\left( \omega ^2-l_2\right) =0 \end{aligned}$$

(6.2)

and

$$\begin{aligned} l_4\sin {\omega \tau }-l_3\omega \cos {\omega \tau }+l_1\omega =0. \end{aligned}$$

(6.3)

With the help of Eqs. (6.2) and (6.3), we have

$$\begin{aligned} \rho _1^2+d_4\rho _1+d_5=0. \end{aligned}$$

(6.4)

where $\rho _1={\hat{\omega }}^2$, $d_4=l_1^2-2l_2-l_3^2$ and $d_5=l_2^2-l_4^2$. Now if $d_4<0,d_5<0$ or $d_4>0,d_5>0$, then this quadratic Eq. (6.4) will have at least one positive root and consequently a positive root for $\omega$, say ${\hat{\omega }}$, and corresponding critical value of $\tau$, say $\tau _\mathrm{c}$, is given by

$$\begin{aligned} \tau _\mathrm{c}^{(n)}&=\frac{1}{{\hat{\omega }}}\arcsin \frac{1}{\Delta _1}\left( l_2l_3{\hat{\omega }}-l_1l_4{\hat{\omega }}-l_3{\hat{\omega }}^3\right) \nonumber \\&\quad +\frac{2n\pi }{{\hat{\omega }}} \qquad \qquad \ \ n=0,1,2,3.... \end{aligned}$$

(6.5)

where $\Delta _1=l_3^2{\hat{\omega }}^2+l_4^2$.

If $E_2$ is locally asymptotically for $\tau =0$, then by G.J. Butler’s lemma, $E_2$ remains stable for $\tau <\tau _\mathrm{c}$. The transversality condition for the Hopf bifurcation at $\tau _\mathrm{c}$ is $Re\left( \frac{d\lambda }{d\tau }\right) _{\tau =\tau _\mathrm{c}}>0$. So differentiating (6.1) with respect to $\tau$ and putting $\tau =\tau _\mathrm{c}$, $\omega ={\hat{\omega }}$, $\mu =0$, we get

$$\begin{aligned} \text {Re}\left( \frac{d\lambda }{d\tau }\right) _{\tau =\tau _\mathrm{c}}=\frac{P_1R_1-Q_1S_1}{P_1^2+Q_1^2}>0 \end{aligned}$$

(6.6)

provided $P_1R_1-Q_1S_1>0$ where $P_1=l_1-l_3\cos {{\hat{\omega }}\tau _\mathrm{c}}+l_3{{\hat{\omega }}\tau _\mathrm{c}}\sin {{\hat{\omega }}\tau _\mathrm{c}}+l_4\tau _\mathrm{c}\cos {{\hat{\omega }}\tau _\mathrm{c}}$, $Q_1=-2{\hat{\omega }}-l_3{\hat{\omega }}\tau _\mathrm{c}\cos {{\hat{\omega }}\tau _\mathrm{c}}-l_3\sin {{\hat{\omega }}\tau _\mathrm{c}}+l_4\tau _\mathrm{c}\sin {{\hat{\omega }}\tau _\mathrm{c}}$, $R_1=l_3{\hat{\omega }}^2\cos {{\hat{\omega }}\tau _\mathrm{c}}-l_4{\hat{\omega }}\sin {{\hat{\omega }}\tau _\mathrm{c}}$, and $S_1=-l_3{\hat{\omega }}^2\sin {{\hat{\omega }}\tau _\mathrm{c}}-l_4{\hat{\omega }}\cos {{\hat{\omega }}\tau _\mathrm{c}}$. From (6.6), it is clear that the condition for the occurrence of Hopf bifurcation at $\tau =\tau _\mathrm{c}$ is established.

Appendix D

Let $\tau =\tau _\mathrm{c}+\mu$ where $\mu =0$ is a Hopf bifurcation value of the system. Define the space of continuous real valued function as $C=C([-1,0],R^3)$. Using the transformation $u_1(t)=S-S^*$, $u_2(t)=I-I^*$ and $u_3(t)=P-P^*$, the system (4.1–4.4) is transformed to a functional differential equation in C as

$$\begin{aligned} \frac{\mathrm{d}U}{\mathrm{d}t}=L_{\mu }(U_t)+H(\mu ,U_t) \end{aligned}$$

(6.7)

where $U=(u_1(t),u_2(t),u_3(t))^T\in R^3$, $U_t(\Theta )=U(t+\Theta ) \ \ \text {for} \ \Theta \in [-1,0]$.

With $\phi (t)=(\phi _1(t),\phi _2(t),\phi _3(t))^T\in C([-\tau ,0],R^3),$$L_\mu :C\rightarrow R^3$ and $H:R\times C\rightarrow R^3$ are given as follows :

$$\begin{aligned} L_\mu (\phi )= & {} (\tau _\mathrm{c}+\mu )[M\phi (0)+N\phi (-1)],\end{aligned}$$

(6.8)

$$\begin{aligned} H(\mu ,\phi )= & {} (\tau _\mathrm{c}+\mu )\nonumber \\&\times \left[ \begin{array}{c} -\frac{r}{k}\phi _1^2(0)-\frac{r}{k}\phi _1(0)\phi _2(0)-\beta F_1\\ \beta F_1-mF_2 \\ m_1F_2-e\phi _3^2(0) \end{array}\right] \nonumber \\ \end{aligned}$$

(6.9)

where

$$\begin{aligned} F_1= & {} \sum \limits _{i+j\ge 2}\frac{1}{i!j!}\left[ \frac{\partial ^{i+j}}{\partial S^i \partial I^j}\left( \frac{IS}{1+\alpha I} \right) \right] _{(S,I)=(S^*,I^*)}\nonumber \\&\phi _1^i(0)\phi _2^j(-1),\end{aligned}$$

(6.10)

$$\begin{aligned} F_2= & {} \sum \limits _{i+j\ge 2}\frac{1}{i!j!}\left[ \frac{\partial ^{i+j}}{\partial I^i \partial P^j}\left( \frac{IP}{(1+aI)(1+bP)} \right) \right] _{(I,P)=(I^*,P^*)}\nonumber \\&\phi _2^i(0)\phi _3^j(0). \end{aligned}$$

(6.11)

Using the Riesz representation theorem, there exists a matrix function $\eta (\theta ,\mu )$ whose components are of bounded variation for $\theta \in [-1,0]$ such that

$$\begin{aligned} L_\mu (\phi )=\int _{-1}^{0}d\eta (\theta ,\mu )\phi (\theta ). \end{aligned}$$

(6.12)

Now we choose

$$\begin{aligned} \eta (\theta ,\mu )=(\tau _\mathrm{c}+\mu )[M\delta (\theta )-N\delta (\theta +1)] \end{aligned}$$

(6.13)

where $\delta$ is the Dirac delta function. For $\phi \in C([-1,0],R^3)$, we define

$$\begin{aligned} A(\mu )\phi&=\left\{ \begin{array}{ll} \frac{d\phi (\theta )}{d\theta }, &{}\quad \theta \in [-1,0)\\ \int _{-1}^{0}d\eta (\mu ,s)\phi (s), &{} \quad \theta =0, \end{array}\right. \end{aligned}$$

(6.14)

$$\begin{aligned} R(\mu )\phi&=\left\{ \begin{array}{ll} 0&{}\qquad \theta \in [-1,0)\\ H(\mu ,\phi ), &{}\qquad \theta =0. \end{array}\right. \end{aligned}$$

(6.15)

Then Eq. (6.7) can be transformed into following operator differential equation

$$\begin{aligned} \frac{\mathrm{d}U}{\mathrm{d}t}=A(\mu )U_t+R(\mu )U_t. \end{aligned}$$

(6.16)

For $\psi \in C^1([-1,0],(R^3)^*)$, the adjoint operator $A^*$ of A is defined as

$$\begin{aligned} A^*\psi (s)=\left\{ \begin{array}{ll} -\frac{d\psi (s)}{ds} &{}\qquad for\ s\in (0,1]\\ \int _{-1}^{0}d\eta ^T(t,0)\psi (-t) &{}\qquad for \ s=0 \end{array}\right. \nonumber \\ \end{aligned}$$

(6.17)

and a bilinear inner product is defined as

$$\begin{aligned}&<\psi (s),\phi (\theta )>\nonumber \\&\quad ={\bar{\psi }}(0)\phi (0)\nonumber \\&\qquad -\int _{-1}^{0} \int _{\xi =0}^{\theta }{\bar{\psi }}(\xi -\theta )d\eta (\theta )\phi (\xi )d\xi \end{aligned}$$

(6.18)

where $\eta (\theta )=\eta (\theta ,0)$. We know that $\pm i\omega _0\tau _\mathrm{c}$ are eigenvalues of A(0) and consequently that are of $A^*$ also. Now we need to find eigenvectors of A and $A^*$ corresponding to $+i\omega _0\tau _\mathrm{c}$ and $-i\omega _0\tau _\mathrm{c}$, respectively. Let us assume $q(\theta )=(1,\alpha _1,\beta _1)^Te^{i\omega _0\tau _\mathrm{c}\theta }$ be the eigenvector of A(0) corresponding to $i\omega _0\tau _\mathrm{c}$ where $\alpha _1=\frac{i\omega _0-m_{11}}{m_{12}+n_{12}e^{-i\omega _0\tau _\mathrm{c}}}$ and $\beta _1=\frac{m_{32}}{i\omega _0-m_{33}}\alpha _1$. Similarly assuming $q^*(s)=D(1,\alpha ^*_1,\beta ^*_1)^Te^{i\omega _0\tau _\mathrm{c}s}$ as the eigenvector of $A^*$ corresponding to $-i\omega _0\tau _\mathrm{c}$, we have $\alpha ^*_1=-\frac{i\omega _0+m_{11}}{m_{21}}$ and $\beta ^*_1=-\frac{1}{m_{32}}[(m_{12}+n_{12}e^{i\omega _0\tau _\mathrm{c}})+(m_{22}+n_{22}e^{i\omega _0\tau _\mathrm{c}}+i\omega _0 )\alpha ^*_1]$. Now we have to calculate the value of D from $<q^*(s),q(\theta )>=1$. Using Eq. (6.18),

$$\begin{aligned}&<q^*(s),q(\theta )>\\&\quad ={\bar{D}}(1,\bar{\alpha ^*_1},\bar{\beta ^*_1}) (1,\alpha _1,\beta _1)^T\\&\qquad -\int _{-1}^{0}\int _{\xi =0}^{\theta }{\bar{D}} (1,\bar{\alpha ^*_1},\bar{\beta ^*_1})e^{-i\omega _0\tau _\mathrm{c} (\xi -\theta )}d\eta (\theta )\\&\qquad (1,\alpha _1,\beta _1)^Te^{i\omega _0\tau _\mathrm{c}\xi }d\xi \\&\quad ={\bar{D}}\left[ 1+\alpha _1\bar{\alpha ^*_1}+\beta _1\bar{\beta ^*_1}\right. \\&\qquad \left. +\frac{\beta S^*}{(1+\alpha I^*)^2}\alpha _1\tau _\mathrm{c}(\bar{\alpha ^*_1}-1) e^{-i\omega _0\tau _\mathrm{c}}\right] \end{aligned}$$

where ${\bar{D}}=\left[ 1+\alpha _1\bar{\alpha ^*_1}+\beta _1\bar{\beta ^*_1}+\frac{\beta S^*}{(1+\alpha I^*)^2}\alpha _1\tau _\mathrm{c}(\bar{\alpha ^*_1}-1)e^{-i\omega _0\tau _\mathrm{c}}\right] ^{-1}$.

Now we will calculate co-ordinates to describe the centre manifold $C_0$ at $\mu =0$. For each $U\in D(A)$, we may choose the pair (y, V) where

$$\begin{aligned} y(t)=\,<q^*,U_t>, \ V(t,\theta )= U_t(\theta )-2\text {Re}\left\{ y(t)q(\theta )\right\} .\nonumber \\ \end{aligned}$$

(6.19)

On the centre manifold $C_0$, we have $V(t,\theta )=V(y(t),{\bar{y}}(t),\theta )$ where

$$\begin{aligned} V(y,{\bar{y}},\theta )= & {} V_{20}(\theta )\frac{y^2}{2}+V_{11}(\theta )y{\bar{y}} +V_{02}(\theta )\frac{{\bar{y}}^2}{2}\nonumber \\&+V_{30}(\theta )\frac{y^3}{6}+\cdots \end{aligned}$$

(6.20)

where y and ${\bar{y}}$ are local co-ordinates for centre manifold $C_0$ in the direction of $q^*$ and $\bar{q^*}$ . Clearly, V is real if $U_t$ is real. We consider real solutions only. For solution $U_t\in C_0$ of Eq. (6.16), since $\mu =0$, we have

$$\begin{aligned} \dot{y}&=\,<q^*,\dot{U_t}>\, = i\omega _0\tau _\mathrm{c}y(t)+g(y,{\bar{y}}) \end{aligned}$$

where

$$\begin{aligned} g(y,{\bar{y}})= & {} \bar{q^*}(0)H_0(y,{\bar{y}})=g_{20}\frac{y^2}{2} +g_{11}y{\bar{y}}\nonumber \\&+g_{02}\frac{{\bar{y}}^2}{2}+g_{21}\frac{y^2{\bar{y}}}{2}+\cdots \end{aligned}$$

(6.21)

From Eq. (6.19), we have

$$\begin{aligned} U_t(\theta )&=V(t,\theta )+2\text {Re}\left\{ y(t)q(\theta )\right\} \\&= V(t,\theta )+yq(\theta )+{\bar{y}}{\bar{q}}(\theta ). \end{aligned}$$

For $U_t(\theta )=(U_{1t}(\theta ),U_{2t}(\theta ),U_{3t}(\theta ))$, it follows from Eq. (6.9) that

$$\begin{aligned} g(y,{\bar{y}})&=\bar{q^*}(0)H_0(y,{\bar{y}})=\bar{q^*}(0)H(0,U_t)\nonumber \\&={\bar{D}}(1,\bar{\alpha ^*_1},\bar{\beta ^*_1)}\tau _\mathrm{c}\nonumber \\&\quad \left[ \begin{array}{c} -\frac{r}{k}U_{1t}^2(0)-\frac{r}{k}U_{1t}(0)U_{2t}(0)-\beta F_3\\ \beta F_3-mF_4 \\ m_1F_4-eU_{3t}^2(0) \end{array}\right] \nonumber \\&={\bar{D}}\tau _\mathrm{c}(r_1y^2+r_2y{\bar{y}}+r_3{\bar{y}}^2+r_4y^2{\bar{y}})+\cdots \end{aligned}$$

(6.22)

where

$$\begin{aligned} F_3= & {} \sum \limits _{i+j\ge 2}\frac{1}{i!j!}\left[ \frac{\partial ^{i+j}}{\partial S^i \partial I^j}\left( \frac{IS}{1+\alpha I} \right) \right] _{(S,I)=(S^*,I^*)}\nonumber \\&U_{1t}^i(0)U_{2t}^j(-1),\end{aligned}$$

(6.23)

$$\begin{aligned} F_4= & {} \sum \limits _{i+j\ge 2}\frac{1}{i!j!}\left[ \frac{\partial ^{i+j}}{\partial I^i \partial P^j}\left( \frac{IP}{(1+aI)(1+bP)} \right) \right] _{(I,P)=(I^*,P^*)}\nonumber \\&U_{2t}^i(0)U_{3t}^j(0). \end{aligned}$$

(6.24)

Comparing Eqs. (6.21) and (6.22), we have $g_{20}=2{\bar{D}}\tau _\mathrm{c}r_1,\ g_{11}={\bar{D}}\tau _\mathrm{c}r_2,\ g_{02}=2{\bar{D}}\tau _\mathrm{c}r_3,\ g_{21}=2{\bar{D}}\tau _\mathrm{c}r_4$. Now we calculate $V_{20}(\theta )$ and $V_{11}(\theta )$ as they are present in $g_{21}$. From Eqs. (6.16) and (6.19), we have

$$\begin{aligned} \dot{V}&=\dot{U_t}-\dot{y}q-\dot{{\bar{y}}}{\bar{q}}\nonumber \\&= \left\{ \begin{array}{ll} AV-2\text {Re}[\bar{q^*}(0)H_0q(\theta )] &{} \qquad -1\le \theta <0\\ AV-2\text {Re}[\bar{q^*}(0)H_0q(\theta )]+H_0 &{}\qquad \theta =0 \end{array}\right. \nonumber \\&=AV+G(y,{\bar{y}},\theta ) \end{aligned}$$

(6.25)

where

$$\begin{aligned} G(y,{\bar{y}},\theta )= & {} G_{20}(\theta )\frac{y^2}{2}+G_{11}(\theta )y{\bar{y}}\nonumber \\&+\,G_{02}(\theta )\frac{{\bar{y}}^2}{2}+\cdots \end{aligned}$$

(6.26)

Again near the origin

$$\begin{aligned} \dot{V}&=V_y\dot{y}+V_{{\bar{y}}}\dot{{\bar{y}}}\nonumber \\&=(V_{20}(\theta )y+V_{11}(\theta ){\bar{y}}+\cdots )(i\omega _0\tau _\mathrm{c}y(t))\nonumber \\&\quad +g(y,{\bar{y}}))+(V_{11}(\theta )y+V_{02}(\theta ){\bar{y}}+\cdots )\nonumber \\&\qquad (-i\omega _0\tau _\mathrm{c}{\bar{y}}(t))+{\bar{g}}(y,{\bar{y}})). \end{aligned}$$

(6.27)

Comparing Eqs. (6.26) and (6.27), we get

$$\begin{aligned}&(A(0)-2i\omega _0\tau _\mathrm{c}I_3)V_{20}(\theta )=-G_{20}(\theta ) \ \ \ \text {and} \nonumber \\&\quad A(0)V_{11}(\theta )=-G_{11}(\theta ).\nonumber \\ \end{aligned}$$

(6.28)

Case 1 :$\theta \in [-1,0)$.

From Eqs. (6.19), (6.20), (6.25) and (6.26) for $\theta \in [-1,0)$, we have

$$\begin{aligned} G(y,{\bar{y}},\theta )&=-g(y,{\bar{y}})q(\theta )-{\bar{g}}(y,{\bar{y}}) {\bar{q}}(\theta )\nonumber \\&=\left( g_{20}(\theta )\frac{y^2}{2}+g_{11}(\theta )y{\bar{y}}\right. \nonumber \\&\left. \quad +\,g_{02}(\theta ) \frac{{\bar{y}}^2}{2}+\cdots \right) \nonumber \\&\quad -\,\left( {\bar{g}}_{20}(\theta )\frac{{\bar{y}}^2}{2} +{\bar{g}}_{11}(\theta )y{\bar{y}}+{\bar{g}}_{02}(\theta )\frac{y^2}{2}+\cdots \right) . \end{aligned}$$

(6.29)

Comparing the coefficient of $y^2$ and $y{\bar{y}}$ between Eqs. (6.26) and (6.28),

$$\begin{aligned} G_{20}(\theta )&=-g_{20}q(\theta )-{\bar{g}}_{20}{\bar{q}}(\theta ) \ \ \ \text {and} \nonumber \\ G_{11}(\theta )&=-g_{11}q(\theta )-{\bar{g}}_{11}{\bar{q}}(\theta ). \end{aligned}$$

(6.30)

From Eqs. (6.28) and (6.8) using the definition of A,

$$\begin{aligned} V_{20}(\theta )&=\frac{ig_{20}}{\omega _0\tau _\mathrm{c}}q(0)e^{i\omega _0\tau _\mathrm{c}\theta }\nonumber \\&\quad +\frac{i{\bar{g}}_{02}}{3\omega _0\tau _\mathrm{c}}{\bar{q}}(0)e^{-i\omega _0\tau _\mathrm{c}\theta } +E_1e^{2i\omega _0\tau _\mathrm{c}\theta }. \end{aligned}$$

(6.31)

Similarly,

$$\begin{aligned} V_{11}(\theta )&=-\frac{ig_{11}}{\omega _0\tau _\mathrm{c}}q(0)e^{i\omega _0\tau _\mathrm{c}\theta }\nonumber \\&\quad +\frac{i{\bar{g}}_{11}}{\omega _0\tau _\mathrm{c}}{\bar{q}}(0)e^{-i\omega _0\tau _\mathrm{c}\theta }+E_2 \end{aligned}$$

(6.32)

where $E_1=(E_1^{(1)},E_1^{(2)},E_1^{(3)})$ and $E_2=(E_2^{(1)},E_2^{(2)},E_2^{(3)})$ are constant vector in $R^{3}$.

Case 2 :$\theta =0$.

From Eqs. (6.25) and (6.26),

$$\begin{aligned} G_{20}(0)&=-g_{20}q(0)-{\bar{g}}_{02}{\bar{q}}(0)+2\tau _\mathrm{c}(c_1,c_2,c_3)^T, \end{aligned}$$

(6.33)

$$\begin{aligned} G_{11}(0)&=-g_{11}q(0)-{\bar{g}}_{11}{\bar{q}}(0)+\tau _\mathrm{c}(c_4,c_5,c_6)^T \end{aligned}$$

(6.34)

where $M_1=(c_1,c_2,c_3)^T$ and $M_2=(c_4,c_5,c_6)^T$ are the coefficient of $y^2$ and $y{\bar{y}}$, respectively, in $H_0(y,{\bar{y}})$ where

$$\begin{aligned} c_1&=-\frac{r}{k}(1+\alpha _1)-\frac{\beta \alpha _1}{(1+\alpha I^*)^{2}} e^{-i\omega _0\tau _\mathrm{c}}\\&\quad +\frac{\beta S^*}{(1+\alpha I^*)^{3}} \alpha _1^2e^{-2i\omega _0\tau _\mathrm{c}},\\ c_2&=\frac{\beta \alpha _1}{(1+\alpha I^*)^{2}}e^{-i\omega _0\tau _\mathrm{c}}\\&\quad -\frac{\beta S^*}{(1+\alpha I^*)^{3}}\alpha _1^2e^{-2i\omega _0\tau _\mathrm{c}} -\frac{m\alpha _1\beta _1}{(1+a I^*)^{2}(1+b P^*)^{2}}\\&\quad +\frac{mP^*}{(1+a I^*)^{3}(1+b P^*)}\alpha _1^2+\frac{mI^*}{(1+a I^*)(1+b P^*)^3} \beta _1^2,\\ c_3&=\frac{m_1\alpha _1\beta _1}{(1+a I^*)^{2}(1+b P^*)^{2}}-\frac{m_1P^*}{(1+a I^*)^{3} (1+b P^*)}\alpha _1^2\\&\quad -\left\{ \frac{m_1I^*}{(1+a I^*)(1+b P^*)^3}+e\right\} \beta _1^2,\\ c_4&=-\frac{r}{k}(2+\alpha _1+\bar{\alpha _1})\\&\quad -\frac{\beta }{(1+\alpha I^*)^{2}} \left\{ \bar{\alpha _1}e^{i\omega _0\tau _\mathrm{c}} +\alpha _1e^{-i\omega _0\tau _\mathrm{c}}\right\} \\&\quad +\frac{2\alpha _1\bar{\alpha _1}\beta S^*}{(1+\alpha I^*)^{3}}\\ c_5&=\frac{\beta }{(1+\alpha I^*)^{2}}\left\{ \bar{\alpha _1}e^{i\omega _0\tau _\mathrm{c}} +\alpha _1e^{-i\omega _0\tau _\mathrm{c}}\right\} \\&\quad -\frac{2\alpha _1\bar{\alpha _1}\beta S^*}{(1+\alpha I^*)^{3}}-\frac{m(\alpha _1\bar{\beta _1}+\bar{\alpha _1}\beta _1)}{(1+a I^*)^{2}(1+b P^*)^{2}}\\&\quad +\frac{2\alpha _1\bar{\alpha _1}mP^*}{(1+a I^*)^{3}(1+b P^*)}\\&\quad +\frac{2\beta _1\bar{\beta _1}mI^*}{(1+a I^*)(1+b P^*)^3},\\ c_6&=\frac{m_1(\alpha _1\bar{\beta _1}+\bar{\alpha _1}\beta _1)}{(1+a I^*)^{2} (1+b P^*)^{2}}-\frac{2\alpha _1\bar{\alpha _1}m_1P^*}{(1+a I^*)^{3}(1+b P^*)}\\&\quad -\left\{ \frac{m_1I^*}{(1+a I^*)(1+b P^*)^3}+e\right\} 2\beta _1\bar{\beta _1}. \end{aligned}$$

Now using Eqs. (6.31) and (6.33),

$$\begin{aligned} \left\{ 2i\omega _0\tau _\mathrm{c}I_3-\int _{-1}^{0}e^{2i\omega _0\tau _\mathrm{c}\theta }d\eta (\theta )\right\} E_1=2\tau _\mathrm{c}M_1. \end{aligned}$$

This gives $M_3E_1=2M_1$ where

$$\begin{aligned} M_3= \left[ \begin{array}{ccc} 2i\omega _0-m_{11}&{} -m_{12}-n_{12}e^{-2i\omega _0\tau _\mathrm{c}} &{} 0 \\ -m_{21} &{} 2i\omega _0-m_{22}-n_{22}e^{-2i\omega _0\tau _\mathrm{c}} &{} -m_{23} \\ 0 &{} -m_{32} &{} 2i\omega _0-m_{33}\\ \end{array}\right] . \end{aligned}$$

By simple calculation, we have $E_1^{(i)}=\frac{2\Delta _1^{(i)}}{\Delta _1}$ where $\Delta _1=\det {(M_3)}$, and $\Delta _1^{(i)}$ is the determinant value of $M_3^{(i)}$ obtained by replacing $i \ th \ (i=1,2,3)$ column vector of $M_3$ by another column vector $(c_1,c_2,c_3)^T$. Similarly we have $M_4E_2=2M_2$ where

$$\begin{aligned} M_4= \left[ \begin{array}{ccc} -m_{11}&{} -m_{12}-n_{12}&{} 0 \\ -m_{21} &{} -m_{22}-n_{22} &{} -m_{23} \\ 0 &{} -m_{32} &{}-m_{33}\\ \end{array}\right] . \end{aligned}$$

This gives $E_2^{(j)}=\frac{2\Delta _2^{(j)}}{\Delta _2}$, where $\Delta _2=\det {(M_4)}$ and $\Delta _2^{(j)}$ is the determinant value of $M_4^{(j)}$ obtained by replacing $j \ th \ (j=1,2,3)$ column vector of $M_4$ by another column vector $(c_4,c_5,c_6)^T$. Thus we can compute $V_{20}(\theta )$ and $V_{11}(\theta )$ from Eqs. (6.31) and (6.32) and consequently $g_{21}$ also. Moreover, the following values can be computed as

$$\begin{aligned} C_1(0)&=\frac{i}{2\omega _0\tau _\mathrm{c}}\left( g_{20}g_{11}-2|g_{11}|^2 -\frac{|g_{02}|^2}{3} \right) +\frac{g_{21}}{2}\nonumber ,\\ \mu _2&=-\frac{\text {Re}\left\{ C_1(0)\right\} }{\text {Re} \left\{ \frac{d\lambda (\tau _\mathrm{c})}{d\tau } \right\} }\nonumber ,\\ \beta _2&=2\text {Re}\left\{ C_1(0)\right\} \nonumber ,\\ T_2&=-\frac{\text {Im}\left\{ C_1(0)\right\} +\mu _2 \text {Im} \left\{ \frac{d\lambda (\tau _\mathrm{c})}{d\tau } \right\} }{\omega _0\tau _\mathrm{c}}. \end{aligned}$$

(6.35)

These quantities determine the direction of Hopf bifurcation at $\tau =\tau _\mathrm{c}$.

Appendix E

The derivatives of S, I and P with respect to $b, \alpha$ and $\beta$ are obtained from following equations :

$$\begin{aligned}&rS\left( 1-\frac{S+I}{k}\right) -\frac{\beta IS}{1+\alpha I}-\gamma S=0,\\&\frac{\beta IS}{1+\alpha I}-\frac{mIP}{(1+aI)(1+bP)}-cI=0,\\&\text {and} \ \ \ \frac{m_1IP}{(1+aI)(1+bP)}-\mathrm{d}P-eP^2=0. \end{aligned}$$

b	$\beta$	$\alpha$
$\frac{dS}{db}=-g_1 \frac{dI}{db}$	$\frac{dS}{d\beta }=-g_1\frac{dI}{d\beta }-\frac{kI}{r(1+\alpha I)}$	$\frac{dS}{d\alpha }=-g_1\frac{dI}{d\alpha }+\frac{\beta kI^2}{r(1+\alpha I)^2}$
$\frac{dI}{db}=\frac{(d+eP)(1+aI)^2P+m_1g_2P^2}{m_1(1-g_2g_3)}$	$\frac{dI}{d\beta }=g_2\frac{\mathrm{d}P}{d\beta }$	$\frac{dI}{d\alpha }=g_2\frac{\mathrm{d}P}{d\alpha }$
$\frac{\mathrm{d}P}{db}=g_3 \frac{dI}{db}+P^2$	$\frac{\mathrm{d}P}{d\beta }=g_4\left\{ \frac{S}{1+\alpha I}-\frac{\beta kI}{r(1+\alpha I)^2} \right\}$	$\frac{\mathrm{d}P}{d\alpha }=\frac{(1+aI)(1+bP)^2}{m(1-g_2g_3)}\left\{ \frac{k\beta ^2 I^2}{r(1+\alpha I)^3}-\frac{\beta IS}{(1+\alpha I)^2}\right\}$

where $g_1=1+\frac{\beta k}{r(1+\alpha I)^2}$, $g_2=\frac{(e+bd+2beP)(1+aI)^2}{m_1}$, $g_3=\left[ \frac{amP}{(1+aI)^2(1+bP)}-\frac{\alpha \beta S}{(1+\alpha I)^2}-\frac{\beta g_1}{1+\alpha I}\right] \frac{(1+aI)(1+bP)^2}{m}$ and $g_4=\left[ \frac{m}{(1+aI)(1+bP)^2}+\frac{\alpha \beta g_2S}{(1+\alpha I)^2}+\frac{\beta g_1g_2}{1+\alpha I}-\frac{amg_2P}{(1+aI)^2(1+bP)} \right] ^{-1}$.