Dengue fever is a critical feverish illness caused by mosquito-borne dengue fever viruses (DENV) [1], which are Flaviviridae flaviviruses. When an infected Aedes aegypti mosquito infects a person, many viruses are transmitted. The transmission of the dengue virus is mostly carried out by humans. A virus-infected mosquito feeds on the blood of a virus-infected person. The virus spreads to tissues such as the duct gland from the mosquito’s gut in 8–10 days. The virus appears to have no negative influence on the mosquito. When a DENV-carrying mosquito bites a human, both the virus and the mosquito’s secretions are injected into the skin. As it circulates throughout the body, it clings to and enters white blood cells, multiplying in them. White blood cells produce a number of signal proteins in response, including interferon, which can lead to symptoms including fever, flu-like symptoms, and excruciating pains.

A variety of organs, including the liver and bone marrow, are regularly damaged by serious infections. They also increase the body’s production of viruses and frequently cause fluid to leak from the circulation into internal organ cavities via the membranes of tiny blood vessels. As a consequence, blood pressure drops to a level where many organs are unable to get enough blood due to decreased blood flow in the blood vessels. The other major effect of dengue fever is bleeding, which was made more likely by bone marrow diseases that decreased the quantity of platelets needed for effective blood coagulation [2].

We provide an associative model that supports the susceptible-infective-recovery-susceptible (SIRS) in the human population and the susceptible-infective (SI) in the vector (mosquito) population while avoiding delays in a trial to examine the dynamics of infection across considerable time periods when susceptible people are born and immunity is lost. Because of vector dynamics, the vector population is typically regarded as being in equilibrium with the human population [3, 4]. The analysis of the delay differential equation in the epidemic model was studied in [5, 6, 7]. The transmission process between humans and mosquitoes takes a long time because of the different parasites’ incubation times [2, 8]. The dynamics of dengue virus transmission with a delayed SIVA model are studied in [9]. Sefidgar et al. [10] discussed the nonlinear system of fractional differential equations that appear in a model of HIV infection of CD4+T cells and proposed the LAM for solving the system. The results showed the effectiveness and efficiency of the method. Wei et al. [11] developed a system for a vector-borne disease with a direct method of transmission, but they also showed how the addition of a latency inside the host-to-vector transmitter term will make the system unstable and cause periodic solutions through Hopf bifurcation.

Omame et al. [12] formulated a mathematical model for the co-infection, COVID-19 and dengue transmission dynamics with optimal control and cost-effectiveness analysis. The results showed that the strategy for implementing control against incident dengue infection is the most cost-effective in controlling dengue and COVID-19 co-infection. Omame et al. [13] design a non-integer ordered model for SARS-CoV-2, dengue, and HIV co-dynamics to assess the impact of SARS-CoV-2 infection on the dynamics of dengue and HIV through fractional derivatives. They showed, using numerical simulations, that keeping the spread of SARS-CoV-2 low would have a significant impact on reducing the co-infections of SARS-CoV-2 and dengue or SARS-CoV-2 and HIV. Based on SIRS-SI in the human and vector populations, Wan and Cui [14] presented a model that utilized two latencies for communication between humans and vector populations with the standard incidence rate, loss of immunity, and rate of recovery from infectiousness. To explore the possibility of equilibrium stability and dynamic behavior, Xu and Zhou [15] projected the dynamics of delayed vector-borne transmission with reinfection. An epidemic model with vaccination and numerous time delays is taken into account in [16] along with the stability and Hopf bifurcation analysis. Baleanu and Babak [17] studied a terminal value problem for nonlinear systems of generalized fractional differential equations and formulated a classical operator and a related weighted space with a generalized fractional operator. The results showed the effects of various choices of weight function on modeling with a TVP. Nowadays, many authors do their research on epidemic models [18, 19, 20] with various strategies. To take into consideration the amount of time required for a viral infection to spread to host and vector populations, Yanxia et al. [21] design an associated upgraded vector-borne epidemic model with two latency periods and reinfection. We created a new model using the information in this article and included three new parameters, including loss of immunity, recovery from infectiousness, and partial protection of the human population.

In this paper, the basic reproduction number for the developed model was determined, and the existence of equilibrium was also examined. The aim of the paper is to explore the stability and Hopf bifurcation of a dynamic model of dengue transmission that incorporates two delays. The numerical simulations were described, and the main conceptual outcomes were exhibited. Fig. 1 provides a schematic overview of this model.

The model considers a uniform mix of human and mosquito populations, ensuring that every mosquito bite provides the same possibility of spreading the virus (or transmitting the virus from an infected human). Because mosquitoes cannot recover from infection, their infection period ends when they die, owing to their incredibly short lifespan. As a result, the mosquito population has a relatively low immunity class, and mortality rates are comparable across all categories.

$H_N$(t) represents the whole human population at time t, which is divided into three compartments: susceptible humans $S_H$(t), infectious humans $I_H$(t), and recovered humans $R_H$(t). So, the entire human population is $H_N$(t) = $S_H$(t) + $I_H$(t) + $R_H$(t).

Similarly, $V_N$(t) represents the whole mosquito population at time t, which is divided into two compartments: susceptible vector $S_V$(t) and infectious vector $I_V$(t). So, the entire vector population is $V_N$(t) = $S_V$(t) + $I_V$(t).

Assumptions of the model:

(a) The human and vector total population sizes are considered to be constant. It is expected that new humans will enter the human population at any time at ${\mathrm{\Omega }}_H$ rate through birth or immigration, and those susceptible mosquitoes will be recruited at a constant rate ${\mathrm{\Omega }}_V$.

(b) Depending on the sickness, people shift from one class to another as their health improves. This infusion is not contagious since there is no vertical transmission or immigration of affected persons. When the Aedes aegypti mosquito bites the host, all people get infected, and dengue development begins.

(c) Natural death occurs at a rate of $d_h$ and $d_v$ [21] (according to their limited life span) for all humans and mosquitoes respectively, regardless of condition.

(d) Individuals who have recovered in the human population acquire partial immunity ($\sigma$) or loss of immunity ($\rho$) [21].

(e) Illness mortality rates for humans and mosquitoes are ${\mu }_h$ and ${\mu }_v$, respectively.

(f) The term $e^{-d_h{\tau }_1}$ and $e^{-d_v{\tau }_2}$ is the human and mosquito survival rate [22].

(g) Mosquitoes do not die or become infected by infection.

The system of non-linear differential equations for the dengue model is

(1)${\displaystyle \frac{d{S}_H}{dt}}={\mathrm{\Omega }}_H-\alpha {\beta }_h{S}_H\left(t-{\tau }_1\right)I_V\left(t-{\tau }_1\right)e^{-d_h{\tau }_1}$ $-d_h{S}_H+\theta I_H+\rho R_H$

(2)${\displaystyle \frac{d{I}_H}{dt}}=\alpha {\beta }_h{S}_H\left(t-{\tau }_1\right)I_V\left(t-{\tau }_1\right)e^{-d_h{\tau }_1}+\sigma \alpha {\beta }_h{R}_H{I}_V$ $-\left(\delta +d_h+{\mu }_h+\theta \right)I_H$

(3)$$\frac{d{R}_H}{dt}=\delta I_H-\sigma \alpha {\beta }_h{R}_H{I}_V-\left(d_h+\rho \right)R_H$$

(4)$$\frac{d{S}_V}{dt}={\mathrm{\Omega }}_V-\alpha {\beta }_V{S}_V\left(t-{\tau }_2\right)I_H\left(t-{\tau }_2\right)e^{-d_v{\tau }_2}-d_v{S}_V$$

(5)$$\frac{d{I}_V}{dt}=\alpha {\beta }_V{S}_V\left(t-{\tau }_2\right)I_H\left(t-{\tau }_2\right)e^{-d_v{\tau }_2}-\left(d_v+{\mu }_v\right)I_V$$

with the following initial conditions:

(6)$S_H(\psi )=S_{0H}(\psi ),I_H(\psi )=I_{0H}(\psi ),R_H(\psi )=R_{0H}(\psi )$ $S_V(\psi )=S_{0V}(\psi ),I_V(\psi )=I_{0V}(\psi )$ $H_N(\psi )=H_{0N}(\psi ),V_N(\psi )=V_{0N}(\psi )$ $\psi \in [-\tau ,0],\tau =\max \{{\tau }_1,{\tau }_2\}$

The systems Eqn. 1 to Eqn. 5 is reformulated as

(7)${\displaystyle \frac{d{S}_H}{dt}}={\mathrm{\Omega }}_H-\alpha {\beta }_h{S}_H\left(t-{\tau }_1\right)I_V\left(t-{\tau }_1\right)e^{-d_h{\tau }_1}-d_h{S}_H$ $+\theta I_H+\rho \left({\displaystyle \frac{{\mathrm{\Omega }}_H}{{\mathrm{d}}_h}}-S_H-I_H\right)$ ${\displaystyle \frac{d{I}_H}{dt}}=\alpha {\beta }_h{S}_H\left(t-{\tau }_1\right)I_V\left(t-{\tau }_1\right)e^{-d_h{\tau }_1}+$ $\sigma \alpha {\beta }_h\left({\displaystyle \frac{{\mathrm{\Omega }}_H}{{\mathrm{d}}_h}}-S_H-I_H\right)I_V-\left(\delta +d_h+{\mu }_h+\theta \right)I_H$ ${\displaystyle \frac{d{I}_V}{dt}}=\alpha {\beta }_V\left({\displaystyle \frac{{\mathrm{\Omega }}_V}{{\mathrm{d}}_v}}-I_V\left(t-{\tau }_2\right)\right)I_H\left(t-{\tau }_2\right)e^{-d_v{\tau }_2}$ $-\left(d_v+{\mu }_v\right)I_V$

The parameters of the above system Eqn. 7 is described in Table 1.

Theorem 2

For ${\tau }_1$, ${\tau }_2$ $\ge$0, illness-free equilibrium $E^0$ is locally asymptotically stable if $R_0$1 and is unstable if $R_0$$>$1.

Proof

The Jacobian matrix of the system Eqn. 7 at $E^0$ is

The characteristic equation of Eqn. 15 is given by

One negative eigen value is $\lambda =-(d_h+\rho )$ and the remaining eigen values of the characteristic equation are

(17)${\lambda }^2+\left[\left(\delta +d_h+{\mu }_h+\theta \right)+\left(d_v+{\mu }_v\right)\right]\lambda +\left(\delta +d_h+{\mu }_h+\theta \right)$ $\left(d_v+{\mu }_v\right)\left[1-R_0{e}^{-\lambda \left({\tau }_1+{\tau }_2\right)}\left(1+{\displaystyle \frac{\sigma \rho }{d_h}}{e}^{d_h{\tau }_1}{e}^{\lambda {\tau }_1}\right)\right]=0$

$\Rightarrow {\lambda }^2+X\lambda +Y=0$

Where $\mathrm{X}=\left(\delta +d_h+{\mu }_h+\theta \right)+\left(d_v+{\mu }_v\right)$

$Y=\left(\delta +d_h+{\mu }_h+\theta \right)\left(d_v+{\mu }_v\right)$ $\left[1-R_0{e}^{-\lambda \left({\tau }_1+{\tau }_2\right)}\left(1+{\displaystyle \frac{\sigma \rho }{d_h}}{e}^{d_h{\tau }_1}{e}^{\lambda {\tau }_1}\right)\right]$

The quadratic equation is the same in the ODE case for ${\tau }_1$ = ${\tau }_2$ = 0. In that case, the real component of every eigenvalue of the characteristic equation is already negative. The Hurwitz criterion states that when ${\tau }_1$ = ${\tau }_2$ = 0, the illness-free equilibrium $E^0$ at $R_0$1 is locally asymptotically stable, however, when $R_0$$>$1, it is unstable.

Let us first consider $R_0$$>$1. It is simple to demonstrate that Eqn. 17 has a real positive root. Rearranging Eqn. 17 in the form

Suppose $\lambda$ is real. From Eqn. 18, the LHS and RHS are denoted by $A(\lambda )$ and $B(\lambda )$ respectively. We have $A(0)$ = 0 and $\underset{\lambda \to \mathrm{\infty }}{lim}A(\lambda )$ = $\mathrm{\infty }$.

Here, $B(\lambda )$ is a decreasing function of $\lambda$

As a result, Eqn. 17 has a non-negative real root for any ${\tau }_1$$\ge$0 and ${\tau }_2$$\ge$0, and the illness-free equilibrium is destabilizing.

Hence $E^0$ is locally asymptotically stable in system Eqn. 7 if $R_0$1 and is unstable if $R_0$$>$1.

Theorem 3

For ${\tau }_1$=${\tau }_2$ = 0, the endemic equilibrium $E^{*}$ is locally asymptotically stable if $R_0$$>$1 and is unstable if $R_0$1.

Proof

The Jacobian matrix of the system Eqn. 7 of $E^{*}$ is

The characteristic equation of Eqn. 19 is given by

(20)${\lambda }^3+k_2{\lambda }^2+k_1\lambda +k_0+\left(l_2{\lambda }^2+l_1\lambda +l_0\right)e^{-\lambda {\tau }_1}+$ $\left(q_2{\lambda }^2+q_1\lambda +q_0\right)e^{-\lambda {\tau }_2}+\left(r_1\lambda +r_0\right)e^{-\lambda \left({\tau }_1+{\tau }_2\right)}=0$

where

When ${\tau }_1$ = 0 and ${\tau }_2$ = 0

The characteristic Eqn. 20 becomes

(21)$${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0$$

where, $A_1=k_2+l_2+q_2$, $A_2=k_1+l_1+q_1+r_1$, $A_3=k_0+l_0+q_0+r_0$

According to the Routh-Hurwitz criterion, the roots of the characteristic Eqn. 21 does not have a positive real part if and only if the coefficients of $A_i$ are non negative and matrix $H_i$$>$0, for $i=0,\mathrm{\hspace{0.25em}1},\mathrm{\hspace{0.25em}2},\mathrm{\hspace{0.25em}3}$. From this $A_1{A}_2-A_3$$>$0.

Hence the endemic equilibrium $E^{*}$ is locally asymptotically stable when ${\tau }_1$ = ${\tau }_2$ = 0.

Theorem 4

For ${\tau }_1$$>$0 and ${\tau }_2$ = 0, the endemic equilibrium $E^{*}$ is locally asymptotically stable if $R_0$$>$1 and is unstable if $R_0$1 [26, 27].

Proof

When ${\tau }_1$$>$0 and ${\tau }_2=0$, the characteristic Eqn. 20 becomes

(22)$${\lambda }^3+B_1{\lambda }^2+B_2\lambda +B_3+\left(C_1{\lambda }^2+C_2\lambda +C_3\right)e^{-\lambda {\tau }_1}=0$$

where $B_1=k_2+q_2$, $B_2=k_1+q_1$, $B_3=k_0+q_0$, $C_1=l_2$, $C_2=l_1+r_1$ and $C_3=l_0+r_0$

Suppose that $\lambda =i{\omega }_0$, ${\omega }_0$$>$0 then the Eqn. 22 becomes

(23)$B_1{\omega }_0^2-B_3=\left(C_3-C_1{\omega }_0^2\right)\cos {\omega }_0{\tau }_1+C_2{\omega }_0\sin {\omega }_0{\tau }_1$ ${\omega }_0^3-B_2{\omega }_0=-\left(C_3-C_1{\omega }_0^2\right)\sin {\omega }_0{\tau }_1+C_2{\omega }_0\cos {\omega }_0{\tau }_1$

It follows that

(24)${\omega }_0^6+\left(B_1^2-2B_2-C_1^2\right){\omega }_0^4+\left(B_2^2-C_2^2-2B_3{B}_1+2C_3{C}_1\right)$ ${\omega }_0^2+B_3^2-C_3^2=0$

Denote in $z=\omega _0{}^2$ in Eqn. 24

(25)$g(z)=z^3+\left(B_1^2-2B_2-C_1^2\right)z^2+$ $(B_2{}^2-C_2^2-2B_3{B}_1+2C_3{C}_1)z+\left(B_3{}^2-C_3^2\right)=0$

(26)$$g(z)=z^3+h_3{z}^2+h_{2}z+h_1=0$$

where $h_3=B_1{}^2-2B_2-C_1{}^2$, $h_2=B_2{}^2-C_2{}^2-2B_3{B}_1+2C_3{C}_1$, $h_1=B_3{}^2-C_3{}^2$

It is easy to get that $g^{{}^{\prime }}(z)=\mathrm{\hspace{0.25em}3}{z}^2+2h_{3}z+h_2=0$ has two roots $z_1=\frac{(-h_3+\sqrt{h_3{}^2-3h_2})}{3}$ and $z_2=\frac{(-h_3-\sqrt{h_3{}^2-3h_2})}{3}$

Clearly, if $h_1$$\ge$0, $h_2$$\ge$0 and $h_3$$\ge$0, then Eqn. 26 has a negative real root. As a result, Eqn. 22 has negative real parts and has no pure imaginary roots.

Suppose that Eqn. 26 has positive root ${\omega }_0$. A pair of imaginary roots ($\pm i{\omega }_0$), appear in the characteristic Eqn. 22. Let $\lambda ({\tau }_1)$ = $\eta ({\tau }_1)$+$i{\omega }_0({\tau }_1)$ be the eigenvalues Eqn. 22 for a certain starting value of the bifurcation parameter ${\tau }_1$ then we have $\eta \left({\tau }_{1k}^n\right)=0$, ${\omega }_0\left({\tau }_{1k}^n\right)={\omega }_0$.

Denote ${\tau }_{1k}^n$ = $\frac{{\theta }_1+2n\pi }{{\omega }_{0k}}$, n = 0, 1, 2, …. Where ${\theta }_1\in [0,\mathrm{\hspace{0.25em}2}\pi ]$

From Eqn. 23, we get

sin (${\theta }_1$) = $\frac{{\omega }_0^5{C}_1+{\omega }_0^3\left(B_1{C}_2-C_3-C_1{B}_2\right)+{\omega }_0(C_3{B}_2-B_3{C}_2)}{{(C_3-C_1{\omega }_0^2)}^2+C_2^2{\omega }_0^2}$

cos (${\theta }_1$) = $-\frac{{\omega }_0^4(B_1{C}_1-C_2)+{\omega }_0^2\left(C_2{B}_2-B_1{C}_3-B_3{C}_1\right)+B_3{C}_3}{{(C_3-C_1{\omega }_0^2)}^2+C_2^2{\omega }_0^2}$

The following transversality condition also verified.

${\frac{d}{d{\tau }_1}\mathrm{Re}(\lambda (\tau ))\parallel }_{\tau ={\tau }_1}={\frac{d}{d{\tau }_1}\eta (\tau )|}_{\tau ={\tau }_1}>0$ holds.

By continuity, when ${\tau }_1$$>$${\tau }_1^{*}$ and the real part of $\lambda ({\tau }_1)$ turns positive, the steady state is unstable. Furthermore, the Hopf bifurcation occurs when ${\tau }_1$ reaches the crucial value in ${\tau }_1^{*}$ [23]. When ${\tau }_2$ = 0, ${\tau }_1$$>$${\tau }_2$, the equilibrium $E^{*}$ is asymptotically stable. However, if ${\tau }_1^{*}$ remains in a certain right neighbourhood of ${\tau }_1$, the equilibrium $E^{*}$ becomes unstable.

Hence, the Hopf bifurcation occurs when ${\tau }_1^{*}={\tau }_1$.

Theorem 5

For ${\tau }_1$ = 0 and ${\tau }_2$$>$0, the endemic equilibrium $E^{*}$ is locally asymptotically stable if $R_0$$>$1 and is unstable if $R_0$1. Hopf bifurcation occurs when ${\tau }_2={\tau }_2^{*}$, at the equilibrium $E^{*}$ of the system Eqn. 7.

Proof

When ${\tau }_1=0$ and ${\tau }_2$$>$0, the characteristic Eqn. 20 becomes

(27)$${\lambda }^3+D_1{\lambda }^2+D_2\lambda +D_3+\left(E_1{\lambda }^2+E_2\lambda +E_3\right)e^{-\lambda {\tau }_2}=0$$

where $D_1=k_2+l_2$, $D_2=k_1+l_1$, $D_3=k_0+l_0$, $E_1=q_2$, $E_2=q_1+r_1$ and $E_3=q_0+r_0$

Suppose that $\lambda =i{\omega }_1$, ${\omega }_1$$>$0 then the Eqn. 27 becomes

(28)$D_1{\omega }_1^2-D_3=\left(E_3-E_1{\omega }_1^2\right)\cos {\omega }_1{\tau }_2+E_2{\omega }_1\sin {\omega }_1{\tau }_2$ ${\omega }_1^3-D_2{\omega }_1=-\left(E_3-E_1{\omega }_1^2\right)\sin {\omega }_1{\tau }_2+E_2{\omega }_1\cos {\omega }_1{\tau }_2$

It follows that

(29)${\omega }_1^6+\left(D_1^2-2D_2-E_1^2\right){\omega }_1^4+\left(D_2^2-E_2^2-2D_3{D}_1+2E_3{E}_1\right){\omega }_1^2$ $+D_3^2-E_3^2=0$

Denote $z=\omega _1{}^2$ in Eqn. 29

(30)$f(z)=z^3+\left(D_1^2-2D_2-E_1{}^2\right)z^2+$ $\left(D_2^2-E_2^2-2D_3{D}_1+2E_3{E}_1\right)z+\left(D_3^2-E_3^2\right)=0$

Let $j_3=D_1^2-2D_2-E_1^2,j_2=D_2{}^2-E_2^2-2D_3{D}_1+2E_3{E}_1,j_1=\left(D_3{}^2-E_3^2\right)$ then Eqn. 30 is

(31)$$f(z)=z^3+j_3{z}^2+j_{2}z+j_1=0$$

It is easy to get that $f^{{}^{\prime }}(z)=\mathrm{\hspace{0.25em}3}{z}^2+2j_{3}z+j_2=0$ has two roots $z_1$ = $\frac{(-j_3+\sqrt{j_3{}^2-3j_2})}{3}$ and $z_2$ = $\frac{(-j_3-\sqrt{j_3{}^2-3j_2})}{3}$

Clearly, if $j_1$$\ge$0, $j_2$$\ge$0 and $j_3$$\ge$0, then Eqn. 31 has negative real roots. Hence, Eqn. 27 has negative real parts and does not have any pure imaginary roots.

If we assume that Eqn. 31, $z_1$, $z_2$ and $z_3$ are three positive roots, then ${\omega }_1$ = $\sqrt{z_1}$, ${\omega }_2$ = $\sqrt{z_2}$ and ${\omega }_3$ = $\sqrt{z_3}$ are three positive numbers.

From Eqn. 28

cos ${\omega }_1{\tau }_2$ = $\frac{\left(D_1\omega _1{}^2-D_3\right)\left(E_3-E_1\omega _1{}^2\right)+E_2{\omega }_1(\omega _1{}^3-D_2{\omega }_1)}{{\left(E_3-E_1\omega _1{}^2\right)}^2+E_2{}^2\omega _1{}^2}$

Denote

${\tau }_{2k}^n=\frac{1}{{\omega }_{1k}}\arccos \left\{\frac{(D_1{\omega }_1^2-D_3)(E_3-E_1{\omega }_1^2)+E_2{\omega }_1\left({\omega }_1{}^3-D_2{\omega }_1\right)}{{(E_3-E_1{\omega }_1^2)}^2+E_2{}^2\omega _1{}^2}\right\}+\frac{2n\pi }{{\omega }_{1k}}$

and ${\tau }_2^{*}$ = $\underset{1\le k\le m}{min}{\tau }_{2k}^{(0)}$, ${\omega }_1={\omega }_i$, $i=1,\mathrm{\hspace{0.25em}2},\mathrm{\hspace{0.25em}3},\mathrm{\dots }.m$

where ${\omega }_1$ is corresponds to ${\tau }_2^{*}$.

Let $\lambda ({\tau }_2)$ = $\eta ({\tau }_2)$+$i{\omega }_1({\tau }_2)$ be the eigenvalues of the characteristic Eqn. 27 satisfying $\eta \left({\tau }_{2k}^n\right)=0$, ${\omega }_1\left({\tau }_{2k}^n\right)={\omega }_1$.

From Eqn. 27, we get ${(\frac{d\lambda }{d{\tau }_2})}^{-1}$ = $\frac{P_1+i{P}_2}{P_3+i{P}_4}$

where,

Therefore, Re ${{(\frac{d\lambda }{d{\tau }_2})}^{-1}|}_{\lambda =i{\omega }_0}$ = $\frac{P_1{P}_3+P_2{P}_4}{P_3^2+P_4^2}$

If $P_1{P}_3+P_2{P}_4\ne 0$ then Re ${{(\frac{d\lambda }{d{\tau }_2})}^{-1}|}_{\lambda =i{\omega }_0}\ne 0$ hold.

If ${\tau }_1=0$ and $P_1{P}_3+P_2{P}_4\ne 0$ the $E^{*}$ is asymptotically stable when ${\tau }_2\in [0,{\tau }_2^{*})$.

Hence, the Hopf bifurcation occurs when ${\tau }_2^{*}={\tau }_2$.