Dynamics of a Delayed Model for the Transmission of Malicious Objects in Computer Network

Abstract

An SEIQRS model for the transmission of malicious objects in computer network with two delays is investigated in this paper. We show that possible combination of the two delays can affect the stability of the model and make the model bifurcate periodic solutions under some certain conditions. For further investigation, properties of the periodic solutions are studied by using the normal form method and center manifold theory. Finally, some numerical simulations are given to justify the theoretical results.

1. Introduction

Computer viruses in network have posed a major threat to our work and life with the rapid popularization of the Internet. Many virus propagation models [1–4] have been proposed to understand the way that computer viruses propagate after Kephart and White [5] proposed the first epidemiological model of computer viruses. In [1], Thommes and Coates proposed a modified version of the SEI model to predict the virus propagation in a network. In [3], Wen and Zhong studied an SIR model on bipartite networks and they proved the existence and the asymptotic stability of the endemic equilibrium by applying the theory of the multigroup model. In [4], Mishra and Jha proposed the following SEIQRS model to describe the transmission of malicious objects in computer network by introducing a new compartment quarantine into the SEIRS model proposed in [2]: 𝑑𝑆(𝑡)(𝑡)𝐼(𝑡)(𝑡)(𝑡),(𝑡)𝑑𝑡=𝐴−𝛽𝑆−𝑑𝑆+𝜂𝑅𝑑𝐸(𝑡)𝐼(𝑡)−𝐸(𝑡),(𝑡)𝑑𝑡=𝛽𝑆𝑑+𝜇𝑑𝐼(𝑡)−𝐼(𝑡),(𝑡)𝑑𝑡=𝜇𝐸𝑑+𝛼+𝛾+𝛿𝑑𝑄(𝑡)−()𝑄(𝑡),(𝑡)𝑑𝑡=𝛿𝐼𝑑+𝛼+𝜀𝑑𝑅(𝑡)(𝑡)−𝑅(𝑡),𝑑𝑡=𝛾𝐼+𝜀𝑄𝑑+𝜂 (1) where 𝑆(𝑡), 𝐸(𝑡), 𝐼(𝑡), 𝑄(𝑡), and 𝑅(𝑡) denote the sizes of nodes in the states susceptible, exposed, infectious, quarantined, and recovered at time 𝑡, respectively. 𝐴 is the rate at which new computers are attached to the network. 𝑑 is the rate at which computers are disconnected to the network. 𝛼 is the crashing rate of computers due to the attack of malicious objects. 𝛽 is the transmission rate. 𝜇, 𝛾, 𝛿, 𝜀, and 𝜂 are the state transition rates.

As is known, an infected computer becomes a recovered one by using antimalicious software and the recovered computer has a temporary immunity, and computer virus models with delay have been studied by many scholars [6–12]. In [6], Ren et al. investigated local and global stability of a delayed viral infection model in computer virus propagation model. In [8], Dong et al. proposed a delayed SEIR computer virus model and studied the problem of Hopf bifurcation of the model by regarding the delay as a bifurcating parameter. Motivated by the work above, Liu [12] incorporated the time delay due to the temporary immunity period into system (1) and proposed the following SEIQRS model with time delay: 𝑑𝑆(𝑡)(𝑡)𝐼(𝑡)(𝑡)(),(𝑡)𝑑𝑡=𝐴−𝛽𝑆−𝑑𝑆+𝜂𝑅𝑡−𝜏𝑑𝐸(𝑡)𝐼(𝑡)−𝐸(𝑡),(𝑡)𝑑𝑡=𝛽𝑆𝑑+𝜇𝑑𝐼(𝑡)−𝐼(𝑡),(𝑡)𝑑𝑡=𝜇𝐸𝑑+𝛼+𝛾+𝛿𝑑𝑄(𝑡)−()𝑄(𝑡),(𝑡)𝑑𝑡=𝛿𝐼𝑑+𝛼+𝜀𝑑𝑅(𝑡)(𝑡)(𝑡)(),𝑑𝑡=𝛾𝐼+𝜀𝑄−𝑑𝑅−𝜂𝑅𝑡−𝜏 (2) where 𝜏≥0 is the time delay due to the temporary immunity period. However, we know that an infected computer needs a period to clean viruses by antivirus software and then becomes a recovered one. Therefore, there is a time delay before the infected computers develop themselves into the recovered ones. And there have been some papers that deal with the research of Hopf bifurcation of dynamical system with multiple delays [13–18]. In [13], Xu and He considered a two-neuron network with resonant bilinear terms and two delays. They studied the problem of Hopf bifurcation by regarding the sum of the two delays as a bifurcation parameter. In [16], Meng et al. studied the Hopf bifurcation of a three-species system with two delays by regarding possible combination of the two delays as a bifurcation parameter. Motivated by the work above, we consider the following SEIQRS computer virus model with two delays in the present paper: 𝑑𝑆(𝑡)(𝑡)𝐼(𝑡)(𝑡)𝑑𝑡=𝐴−𝛽𝑆−𝑑𝑆+𝜂𝑅𝑡−𝜏₁,𝑑𝐸(𝑡)(𝑡)𝐼(𝑡)−𝐸(𝑡),(𝑡)𝑑𝑡=𝛽𝑆𝑑+𝜇𝑑𝐼(𝑡)−()𝐼(𝑡)𝑑𝑡=𝜇𝐸𝑑+𝛼+𝛿−𝛾𝐼𝑡−𝜏₂,𝑑𝑄(𝑡)(𝑡)−()𝑄(𝑡)𝑑𝑡=𝛿𝐼𝑑+𝛼−𝜀𝑄𝑡−𝜏₂,𝑑𝑅(𝑡)𝑑𝑡=𝛾𝐼𝑡−𝜏₂+𝜀𝑄𝑡−𝜏₂−𝑑𝑅(𝑡)−𝜂𝑅𝑡−𝜏₁, (3) where 𝜏₁ is the time delay due to the temporary immunity period and 𝜏₂ is the time delay due to the period that the infected computer uses to clean viruses by antivirus software.

The main purpose of this paper is to investigate the effects of the two delays on system (3) and the remainder of this paper is organized as follows. Sufficient conditions for local stability and existence of local Hopf bifurcation are obtained by analyzing the distribution of the roots of the associated characteristic equation in Section 2. Properties of the Hopf bifurcation are further investigated by using the normal form method and center manifold theory in Section 3. In Section 4, we give a numerical example to support the theoretical results in the paper.

2. Local Stability and Existence of Local Hopf Bifurcation

By a simple computation, it is easy to get that if 𝑅₀=𝐴𝜇𝛽/𝑑(𝑑+𝜇)(𝛼+𝛿+𝛾+𝑑)>1, then system (3) has a unique positive equilibrium 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗), where 𝑆^∗=𝑑+𝜇𝑑+𝛼+𝛿+𝛾𝜇𝛽,𝐸^∗=𝑑+𝛼+𝛿+𝛾𝐼^∗𝜇,𝐼^∗=𝑑+𝜂(𝑑+𝛼+𝜀)𝐴𝜇𝛽−𝑑𝑑+𝜇𝑑+𝛼+𝛿+𝛾𝛽,𝑄^∗=𝛿𝐼^∗𝑑+𝛼+𝜀,𝑅^∗=𝛾𝐼^∗+𝑑+𝜂𝜀𝛿𝐼^∗(),𝑑+𝜂𝑑+𝛼+𝜀 (4) and 𝑅₀ is the basic reproduction number. It is easy to get the linearization of system (3) at 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗): 𝑑𝑆(𝑡)𝑑𝑡=𝑎₁₁𝑆(𝑡)+𝑎₁₃𝐼(𝑡)+𝑏₁₅𝑅𝑡−𝜏₁,𝑑𝐸(𝑡)𝑑𝑡=𝑎₂₁𝑆(𝑡)+𝑎₂₂𝐸(𝑡)+𝑎₂₃𝐼(𝑡),𝑑𝐼(𝑡)𝑑𝑡=𝑎₃₂𝐸(𝑡)+𝑎₃₃𝐼(𝑡)+𝑐₃₃𝐼𝑡−𝜏₂,𝑑𝑄(𝑡)𝑑𝑡=𝑎₄₃𝐼(𝑡)+𝑎₄₄𝑄(𝑡)+𝑐₄₄𝑄𝑡−𝜏₂,𝑑𝑅(𝑡)𝑑𝑡=𝑎₅₅𝑅(𝑡)+𝑏₅₅𝑅𝑡−𝜏₁+𝑐₅₃𝐼𝑡−𝜏₂+𝑐₅₄𝑄𝑡−𝜏₂, (5) where 𝑎₁₁=−𝛽𝐼^∗+𝑑,𝑎₁₃=−𝛽𝑆^∗,𝑎₂₁=𝛽𝐼^∗,𝑎₂₂=−𝑑+𝜇,𝑎₂₃=𝛽𝑆^∗,𝑎₃₂𝑎=𝜇,₃₃=−(𝑑+𝛼+𝛿),𝑎₄₃=𝛿,𝑎₄₄(),𝑎=−𝑑+𝛼₅₅=−𝑑,𝑏₁₅=𝜂,𝑏₅₅=−𝜂,𝑐₃₃𝑐=−𝛾,₄₄=−𝜀,𝑐₅₃=𝛾,𝑐₅₄=𝜀. (6) Thus, the characteristic equation of system (5) is 𝜆⁵+𝑚₄𝜆⁴+𝑚₃𝜆³+𝑚₂𝜆²+𝑚₁𝜆+𝑚₀+𝑛₄𝜆⁴+𝑛₃𝜆³+𝑛₂𝜆²+𝑛₁𝜆+𝑛₀𝑒^−𝜆𝜏1+𝑝₄𝜆⁴+𝑝₃𝜆³+𝑝₂𝜆²+𝑝₁𝜆+𝑝₀𝑒^−𝜆𝜏2+𝑞₃𝜆³+𝑞₂𝜆²+𝑞₁𝜆+𝑞₀𝑒^{−2𝜆𝜏2}+𝑟₃𝜆³+𝑟₂𝜆²+𝑟₁𝜆+𝑟₀𝑒^{−𝜆(𝜏1+𝜏2)}+𝑠₂𝜆²+𝑠₁𝜆+𝑠₀𝑒^{−𝜆(𝜏1+2𝜏2)}=0, (7) where 𝑚₀=𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₃𝑎₂₁𝑎₃₂𝑎₄₄𝑎₅₅,𝑚₁=𝑎₄₄+𝑎₅₅𝑎₁₃𝑎₂₁𝑎₃₂−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₃𝑎₃₂𝑎₄₄𝑎₅₅+𝑎₁₁𝑎₂₂𝑎₃₃𝑎₄₄+𝑎₁₁𝑎₂₂𝑎₅₅𝑎₃₃+𝑎₄₄+𝑎₃₃𝑎₄₄𝑎₅₅𝑎₁₁+𝑎₂₂,𝑚₂=𝑎₂₃𝑎₃₂𝑎₁₁+𝑎₄₄+𝑎₅₅−𝑎₅₅𝑎₁₁𝑎₂₂+𝑎₃₃𝑎₄₄−𝑎₁₃𝑎₂₁𝑎₃₂−𝑎₁₁𝑎₂₂𝑎₃₃+𝑎₄₄−𝑎₃₃𝑎₄₄𝑎₁₁+𝑎₂₂−𝑎₅₅𝑎₁₁+𝑎₂₂𝑎₃₃+𝑎₄₄,𝑚₃=𝑎₁₁𝑎₂₂+𝑎₃₃𝑎₄₄+𝑎₁₁+𝑎₂₂𝑎₃₃+𝑎₄₄−𝑎₂₃𝑎₃₂+𝑎₅₅𝑎₁₁+𝑎₂₂+𝑎₃₃+𝑎₄₄,𝑚₄=−𝑎₁₁+𝑎₂₂+𝑎₃₃+𝑎₄₄+𝑎₅₅,𝑛₀=𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₃𝑎₂₁𝑎₃₂𝑎₄₄𝑏₅₅,𝑛₁=𝑎₁₁𝑎₂₂𝑏₅₅𝑎₃₃+𝑎₄₄+𝑎₃₃𝑎₄₄𝑏₅₅𝑎₁₁+𝑎₂₂+𝑎₁₃𝑎₂₁𝑎₃₂𝑏₅₅−𝑎₂₃𝑎₃₂𝑏₅₅𝑎₁₁+𝑎₄₄,𝑛₂=𝑎₂₃𝑎₃₂−𝑎₁₁𝑎₂₂−𝑎₃₃𝑎₄₄𝑏₅₅−𝑎₁₁+𝑎₂₂𝑎₃₃+𝑎₄₄𝑏₅₅,𝑛₃=𝑏₅₅𝑎₁₁+𝑎₂₂+𝑎₃₃+𝑎₄₄,𝑛₄=−𝑏₅₅,𝑝₀=𝑎₁₁𝑎₂₃−𝑎₁₃𝑎₂₁𝑎₃₂𝑎₅₅𝑐₄₄−𝑎₄₄𝑐₃₃+𝑎₃₃𝑐₄₄𝑎₁₁𝑎₂₂𝑎₅₅,𝑝₁=𝑎₁₁𝑎₂₂𝑐₃₃𝑎₄₄+𝑎₅₅+𝑎₄₄𝑎₅₅𝑐₃₃𝑎₁₁+𝑎₂₂+𝑎₁₁𝑎₂₂𝑐₄₄𝑎₃₃+𝑎₅₅−𝑎₂₃𝑎₃₂𝑐₄₄𝑎₁₁+𝑎₅₅+𝑎₃₃𝑎₅₅𝑐₄₄𝑎₁₁+𝑎₂₂−𝑎₁₃𝑎₂₁𝑎₃₂𝑐₄₄,𝑝₂=𝑎₂₃𝑎₃₂𝑐₄₄−𝑐₃₃𝑎₁₁+𝑎₂₂𝑎₄₄+𝑎₅₅−𝑐₃₃𝑎₁₁𝑎₂₂+𝑎₄₄𝑎₅₅−𝑐₄₄𝑎₁₁+𝑎₂₂𝑎₃₃+𝑎₅₅−𝑐₄₄𝑎₁₁𝑎₂₂+𝑎₃₃𝑎₅₅,𝑝₃=𝑐₃₃𝑎₁₁+𝑎₂₂+𝑎₄₄+𝑎₅₅+𝑐₄₄𝑎₁₁+𝑎₂₂+𝑎₃₃+𝑎₅₅,𝑝₄=−𝑐₃₃+𝑐₄₄,𝑞₀=𝑎₁₁𝑎₂₂𝑎₅₅𝑐₃₃𝑐₄₄,𝑞₃=𝑐₃₃𝑐₄₄,𝑞₁=𝑎₁₁𝑎₂₂+𝑎₁₁𝑎₅₅+𝑎₂₂𝑎₅₅𝑐₃₃𝑐₄₄,𝑞₂=−𝑎₁₁+𝑎₂₂+𝑎₅₅𝑐₃₃𝑐₄₄,𝑟₀=𝑎₂₁𝑎₃₂𝑏₁₅𝑎₄₄𝑐₅₃−𝑎₄₃𝑐₅₄−𝑎₁₁𝑎₂₂𝑏₅₅𝑎₃₃𝑐₄₄+𝑎₄₄𝑐₃₃−𝑎₃₂𝑏₅₅𝑐₄₄𝑎₁₁𝑎₂₃+𝑎₁₃𝑎₂₁,𝑟₁=𝑎₁₁𝑎₂₂+𝑎₁₁𝑎₄₄+𝑎₂₂𝑎₄₄𝑏₅₅𝑐₃₃−𝑎₂₁𝑐₃₂𝑏₁₅𝑐₅₃+𝑎₁₁𝑎₂₂+𝑎₁₁𝑎₃₃+𝑎₂₂𝑎₃₃𝑏₅₅𝑐₄₄+𝑎₂₃𝑎₃₂𝑏₅₅𝑐₄₄,𝑟₂=−𝑎₁₁+𝑎₂₂+𝑎₃₃𝑏₅₅𝑐₄₄−𝑎₁₁+𝑎₂₂+𝑎₄₄𝑏₅₅𝑐₃₃,𝑟₃=𝑐₃₃+𝑐₄₄𝑏₅₅,𝑠₀=𝑎₂₁𝑎₃₂𝑏₁₅𝑐₄₄𝑐₅₃−𝑎₁₁𝑎₂₂𝑏₅₅𝑐₃₃𝑐₄₄,𝑠₁=𝑎₁₁+𝑎₂₂𝑏₅₅𝑐₃₃𝑐₄₄,𝑠₂=−𝑏₅₅𝑐₃₃𝑐₄₄. (8)

Case 1 (𝜏₁=𝜏₂=0). When 𝜏₁=𝜏₂=0, (7) becomes 𝜆⁵+𝐴₁₄𝜆⁴+𝐴₁₃𝜆³+𝐴₁₂𝜆²+𝐴₁₁𝜆+𝐴₁₀𝐴=0,₁₀=𝑚₀+𝑛₀+𝑝₀+𝑞₀+𝑟₀+𝑠₀,𝐴₁₁=𝑚₁+𝑛₁+𝑝₁+𝑞₁+𝑟₁+𝑠₁,𝐴₁₂=𝑚₂+𝑛₂+𝑝₂+𝑞₂+𝑟₂+𝑠₂,𝐴₁₃=𝑚₃+𝑛₃+𝑝₃+𝑞₃+𝑟₃,𝐴₁₄=𝑚₄+𝑛₄+𝑝₄=𝛽𝐼^∗+5𝑑+2𝛼+𝜇+𝛿+𝜀+𝛾+𝜂>0. (9)
Let det₁=𝐴₁₄. Obviously, det₁>0. Therefore, if the condition (𝐻₁): (10) holds, then the positive equilibrium 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗) of system (3) is locally asymptotically stable without delay. Consider det₂=𝐴₁₄1𝐴₁₂𝐴₁₃>0,det₃=𝐴₁₄𝐴10₁₂𝐴₁₃𝐴₁₄0𝐴₁₁𝐴₁₂>0,det₄=𝐴₁₄𝐴100₁₂𝐴₁₃𝐴₁₄1𝐴₁₀𝐴₁₁𝐴₁₂𝐴₁₃00𝐴₁₀𝐴₁₁>0,det₅=𝐴₁₄𝐴1000₁₂𝐴₁₃𝐴₁₄𝐴10₁₀𝐴₁₁𝐴₁₂𝐴₁₃𝐴₁₄00𝐴₁₀𝐴₁₁𝐴₁₂0000𝐴₁₀>0. (10)

Case 2 (𝜏₁>0, 𝜏₂=0). When 𝜏₂=0, (7) becomes the following form: 𝜆⁵+𝐴₂₄𝜆⁴+𝐴₂₃𝜆³+𝐴₂₂𝜆²+𝐴₂₁𝜆+𝐴₂₀+𝐵₂₄𝜆⁴+𝐵₂₃𝜆³+𝐵₂₂𝜆²+𝐵₂₁𝜆+𝐵₂₀𝑒^−𝜆𝜏1=0, (11) where 𝐴₂₀=𝑚₀+𝑝₀+𝑞₀,𝐴₂₁=𝑚₁+𝑝₁+𝑞₁,𝐴₂₂=𝑚₂+𝑝₂+𝑞₂,𝐴₂₃=𝑚₃+𝑝₃+𝑞₃,𝐴₂₄=𝑚₄+𝑝₄,𝐵₂₀=𝑛₀+𝑟₀+𝑠₀,𝐵₂₁=𝑛₁+𝑟₁+𝑠₁,𝐵₂₂=𝑛₂+𝑟₂+𝑠₂,𝐵₂₃=𝑛₃+𝑟₃,𝐵₂₄=𝑛₄. (12) Let 𝜆=𝑖𝜔₁ (𝜔₁>0) be a root of (11). Then, we obtain 𝐵₂₁𝜔₁−𝐵₂₃𝜔³₁sin𝜏₁𝜔₁+𝐵₂₄𝜔⁴₁−𝐵₂₂𝜔²₁+𝐵₂₀cos𝜏₁𝜔₁=𝐴₂₂𝜔²₁−𝐴₂₄𝜔⁴₁−𝐴₂₀,𝐵₂₁𝜔₁−𝐵₂₃𝜔³₁cos𝜏₁𝜔₁−𝐵₂₄𝜔⁴₁−𝐵₂₂𝜔²₁+𝐵₂₀sin𝜏₁𝜔₁=−𝜔⁵₁+𝐴₂₃𝜔³₁−𝐴₂₁𝜔₁. (13) It follows that 𝜔₁¹⁰+𝑐₂₄𝜔⁸₁+𝑐₂₃𝜔⁶₁+𝑐₂₂𝜔⁴₁+𝑐₂₁𝜔²₁+𝑐₂₀=0, (14) with 𝑐₂₀=𝐴²₂₀−𝐵²₂₀,𝑐₂₁=𝐴²₂₁−𝐵²₂₁−2𝐴₂₀𝐴₂₂+2𝐵₂₀𝐵₂₂,𝑐₂₂=𝐴²₂₂−𝐵²₂₂+2𝐴₂₀𝐴₂₄−2𝐴₂₁𝐴₂₃−2𝐵₂₀𝐵₂₄+2𝐵₂₁𝐵₂₃,𝑐₂₃=𝐴²₂₃+2𝐴₂₁−2𝐴₂₂𝐴₂₄−𝐵²₂₃+2𝐵₂₂𝐵₂₄,𝑐₂₄=𝐴²₂₄−𝐵²₂₄−2𝐴₂₃. (15) Let 𝜔²₁=𝑣₁, then (14) becomes 𝑣⁵₁+𝑐₂₄𝑣⁴₁+𝑐₂₃𝑣³₁+𝑐₂₂𝑣²₁+𝑐₂₁𝑣₁+𝑐₂₀=0. (16)
If all the parameters of system (3) are given, one can get all the roots of (16) by the software package Matlab. In order to give the main results in this paper, we make the following assumption.
(𝐻₂₁)  (16) has at least one positive real root.
If the condition (𝐻₂₁) holds, then there exists a 𝑣₁₀ such that (11) has a pair of purely imaginary roots ±𝑖𝜔₁₀=±𝑖√𝑣₁₀. For 𝜔₁₀, the corresponding critical value of time delay is 𝜏₁₀=1𝜔₁₀𝑝arccos²⁸𝜔⁸₁₀+𝑝₂₆𝜔⁶₁₀+𝑝₂₄𝜔⁴₁₀+𝑝₂₂𝜔²₁₀+𝑝₂₀𝑞₂₈𝜔⁸₁₀+𝑞₂₆𝜔⁶₁₀+𝑞₂₄𝜔⁴₁₀+𝑞₂₂𝜔²₁₀+𝑞₂₀, (17) where 𝑝₂₀=−𝐴₂₀𝐵₂₀,𝑝₂₂=𝐴₂₀𝐵₂₂−𝐴₂₁𝐵₂₁+𝐴₂₂𝐵₂₀,𝑝₂₄=𝐴₂₁𝐵₂₃+𝐴₂₃𝐵₂₁−𝐴₂₀𝐵₂₄−𝐴₂₂𝐵₂₂−𝐴₂₄𝐵₂₀,𝑝₂₆=𝐴₂₂𝐵₂₄+𝐴₂₄𝐵₂₂−𝐴₂₃𝐵₂₃−𝐵₂₁,𝑝₂₈=𝐵₂₃−𝐴₂₄𝐵₂₄,𝑞₂₀=𝐵²₂₀,𝑞₂₂=𝐵²₂₁−2𝐵₂₀𝐵₂₂,𝑞₂₄=𝐵²₂₂+2𝐵₂₀𝐵₂₄−2𝐵₂₁𝐵₂₃,𝑞₂₆=𝐵²₂₃−2𝐵₂₂𝐵₂₄,𝑞₂₈=𝐵²₂₄. (18) Taking the derivative of 𝜆 with respect to 𝜏₁ in (11), one can obtain 𝑑𝜆𝑑𝜏₁⁻¹=−5𝜆⁴+4𝐴₂₄𝜆³+3𝐴₂₃𝜆²+2𝐴₂₂𝜆+𝐴₂₁𝜆𝜆⁵+𝐴₂₄𝜆⁴+𝐴₂₃𝜆³+𝐴₂₂𝜆²+𝐴₂₁𝜆+𝐴₂₀+4𝐵₂₄𝜆³+3𝐵₂₃𝜆²+2𝐵₂₂𝜆+𝐵₂₁𝜆𝐵₂₄𝜆⁴+𝐵₂₃𝜆³+𝐵₂₂𝜆²+𝐵₂₁𝜆+𝐵₂₀−𝜏₁𝜆. (19) Thus, Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔10}=𝑓^′₁𝑣_1∗𝐵₂₁𝜔₁₀−𝐵₂₃𝜔³₁₀²+𝐵₂₄𝜔⁴₁₀−𝐵₂₂𝜔²₁₀+𝐵₂₀², (20) where 𝑓₁(𝑣₁)=𝑣⁵₁+𝑐₂₄𝑣⁴₁+𝑐₂₃𝑣³₁+𝑐₂₂𝑣²₁+𝑐₂₁𝑣₁+𝑐₂₀ and 𝑣_1∗=𝜔²₁₀.
Obviously, if the condition (𝐻₂₂)  𝑓^′₁(𝑣_1∗)≠0 holds, then Re[𝑑𝜆/𝑑𝜏₁]⁻¹_{𝜆=𝑖𝜔10}≠0. According to the Hopf bifurcation theorem in [19], we have the following results for system (3).

Theorem 1. For system (3), if the conditions (𝐻₂₁)-(𝐻₂₂) hold, then the positive equilibrium 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗) of system (3) is asymptotically stable for 𝜏₁∈[0,𝜏₁₀) and system (3) undergoes a Hopf bifurcation at the positive equilibrium 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗) when 𝜏₁=𝜏₁₀.

Case 3 (𝜏₁=0, 𝜏₂>0). When 𝜏₁=0, (7) becomes 𝜆⁵+𝐴₃₄𝜆⁴+𝐴₃₃𝜆³+𝐴₃₂𝜆²+𝐴₃₁𝜆+𝐴₃₀+𝐵₃₄𝜆⁴+𝐵₃₃𝜆³+𝐵₃₂𝜆²+𝐵₃₁𝜆+𝐵₂₀𝑒^−𝜆𝜏2+𝐶₃₃𝜆³+𝐶₃₂𝜆²+𝐶₃₀𝑒^{−2𝜆𝜏2}=0, (21) where 𝐴₃₀=𝑚₀+𝑛₀,𝐴₃₁=𝑚₁+𝑛₁,𝐴₃₂=𝑚₂+𝑛₂,𝐴₃₃=𝑚₃+𝑛₃,𝐴₃₄=𝑚₄+𝑛₄,𝐵₃₀=𝑝₀+𝑟₀,𝐵₃₁=𝑝₁+𝑟₁,𝐵₃₂=𝑝₂+𝑟₂,𝐵₃₃=𝑝₃+𝑟₃,𝐵₃₄=𝑝₄,𝐶₃₀=𝑞₀+𝑠₀,𝐶₃₁=𝑞₁+𝑠₁,𝐶₃₂=𝑞₂+𝑠₂,𝐶₃₃=𝑞₃. (22) Multiplying 𝑒^𝜆𝜏2 on both sides of (21), we have 𝐵₃₄𝜆⁴+𝐵₃₃𝜆³+𝐵₃₂𝜆²+𝐵₃₁𝜆+𝐵₂₀+𝜆⁵+𝐴₃₄𝜆⁴+𝐴₃₃𝜆³+𝐴₃₂𝜆²+𝐴₃₁𝜆+𝐴₃₀𝑒^𝜆𝜏2+𝐶₃₃𝜆³+𝐶₃₂𝜆²+𝐶₃₀𝑒^−𝜆𝜏2=0. (23) Let 𝜆=𝑖𝜔₂ (𝜔₂>0) be the root of (23), then we obtain 𝑀₃₁sin𝜏₂𝜔₂+𝑀₃₂cos𝜏₂𝜔₂=𝑀₃₃,𝑀₃₄cos𝜏₂𝜔₂−𝑀₃₅sin𝜏₂𝜔₂=𝑀₃₆, (24) where 𝑀₃₁=𝐴₃₄𝜔⁴₂+𝐶₃₂−𝐴₃₂𝜔²₂+𝐴₃₀−𝐶₃₀,𝑀₃₂=𝜔⁵₂−𝐴₃₃+𝐶₃₃𝜔³₂+𝐴₃₁+𝐶₃₁𝜔₂,𝑀₃₃=𝐵₃₃𝜔³₂−𝐵₃₁𝜔₂,𝑀₃₄=𝐴₃₄𝜔⁴₂−𝐴₃₂+𝐶₃₂𝜔²₂+𝐴₃₀+𝐶₃₀,𝑀₃₅=𝜔⁵₂−𝐴₃₃−𝐶₃₃𝜔³₂+𝐴₃₁−𝐶₃₁𝜔₂,𝑀₃₆=𝐵₃₂𝜔²₂−𝐵₃₄𝜔⁴₂−𝐵₃₀. (25) Then, we obtain cos𝜏₂𝜔₂=𝑝₃₈𝜔⁸₂+𝑝₃₆𝜔⁶₂+𝑝₃₄𝜔⁴₂+𝑝₃₂𝜔²₂+𝑝₃₀𝜔₂¹⁰+𝑞₃₈𝜔⁸₂+𝑞₃₆𝜔⁶₂+𝑞₃₄𝜔⁴₂+𝑞₃₂𝜔²₂+𝑞₃₀,sin𝜏₂𝜔₂=𝑝₃₇𝜔⁷₂+𝑝₃₅𝜔⁵₂+𝑝₃₃𝜔³₂+𝑝₃₁𝜔₂𝜔₂¹⁰+𝑞₃₈𝜔⁸₂+𝑞₃₆𝜔⁶₂+𝑞₃₄𝜔⁴₂+𝑞₃₂𝜔²₂+𝑞₃₀, (26) where 𝑝₃₀=𝐴₃₀𝐶₃₀−𝐵₃₀,𝑝₃₁=𝐴₃₁+𝐶₃₁𝐵₃₀−𝐴₃₀+𝐶₃₀𝐵₃₁,𝑝₃₂=𝐵₃₁𝐶₃₁−𝐴₃₁+𝐴₃₂𝐵₃₀+𝐴₃₀𝐵₃₂−𝐵₃₀𝐶₃₂−𝐵₃₂𝐶₃₀,𝑝₃₃=𝐵₃₁𝐴₃₂+𝐶₃₂+𝐵₃₃𝐴₃₀+𝐶₃₀−𝐵₃₀𝐴₃₃+𝐶₃₃−𝐵₃₂𝐴₃₁+𝐶₃₁,𝑝₃₄=𝐵₃₁𝐴₃₃−𝐶₃₃+𝐵₃₂𝐶₃₂−𝐴₃₂+𝐵₃₃𝐴₃₁−𝐶₃₁+𝐵₃₄𝐶₃₀−𝐴₃₀−𝐴₃₄𝐵₃₀,𝑝₃₅=𝐵₃₂𝐴₃₃+𝐶₃₃+𝐵₃₄𝐴₃₁+𝐶₃₁−𝐵₃₃𝐴₃₂+𝐶₃₂−𝐴₃₄𝐵₃₁,𝑝₃₆=𝐴₃₄𝐵₃₂+𝐴₃₂𝐵₃₄−𝐴₃₃𝐵₃₃−𝐵₃₁+𝐵₃₃𝐶₃₃−𝐵₃₄𝐶₃₂,𝑝₃₇=𝐴₃₄𝐵₃₃−𝐵₃₂−𝐴₃₃+𝐶₃₃𝐵₃₄,𝑝₃₈=𝐵₃₃−𝐴₃₄𝐶₃₄,𝑞₃₀=𝐴²₃₀−𝐶²₃₀,𝑞₃₂=𝐴²₃₁−𝐶²₃₁−2𝐴₃₀𝐴₃₂+2𝐶₃₀𝐶₃₂,𝑞₃₄=𝐴²₃₂−𝐶²₃₂−2𝐴₃₀𝐴₃₄+2𝐶₃₁𝐶₃₃−2𝐴₃₁𝐴₃₃,𝑞₃₆=𝐴²₃₃−𝐶²₃₃+2𝐴₃₁−2𝐴₃₂𝐴₃₄,𝑞₃₈=𝐴²₃₄−2𝐴₃₃. (27) Then, we obtain 𝜔₂²⁰+𝑐₃₉𝜔₂¹⁸+𝑐₃₈𝜔₂¹⁶+𝑐₃₇𝜔₂¹⁴+𝑐₃₆𝜔₂¹²+𝑐₃₅𝜔₂¹⁰+𝑐₃₄𝜔⁸₂+𝑐₃₃𝜔⁶₂+𝑐₃₂𝜔⁴₂+𝑐₃₁𝜔²₂+𝑐₃₀=0, (28) where 𝑐₃₀=𝑞²₃₀−𝑝²₃₀,𝑐₃₁=2𝑞₃₀𝑞₃₂−2𝑝₃₀𝑝₃₂−𝑝²₃₁,𝑐₃₂=𝑞²₃₂−𝑝²₃₂+2𝑞₃₀𝑞₃₄−2𝑝₃₀𝑝₃₄−2𝑝₃₁𝑝₃₃,𝑐₃₃=2𝑞₃₀𝑞₃₆+2𝑞₃₂𝑞₃₄−2𝑝₃₀𝑝₃₆−2𝑝₃₁𝑝₃₅−2𝑝₃₂𝑝₃₄−𝑝²₃₃,𝑐₃₄=𝑞²₃₄−𝑝²₃₄+2𝑞₃₀𝑞₃₈+2𝑞₃₂𝑞₃₆−2𝑝₃₀𝑝₃₈−2𝑝₃₁𝑝₃₇−2𝑝₃₂𝑝₃₆−2𝑝₃₃𝑝₃₅,𝑐₃₅=2𝑞₃₂𝑞₃₈−2𝑝₃₂𝑝₃₈−2𝑝₃₃𝑝₃₇−2𝑝₃₄𝑝₃₆−𝑝²₃₅,𝑐₃₆=𝑞²₃₆−𝑝²₃₆+2𝑞₃₄𝑞₃₈+2𝑞₃₂−2𝑝₃₄𝑝₃₈−2𝑝₃₅𝑝₃₇,𝑐₃₇=2𝑞₃₄+2𝑞₃₆𝑞₃₈−𝑝₃₆𝑝₃₈−𝑝²₃₇,𝑐₃₈=𝑞²₃₈−𝑝²₃₈+2𝑞₃₆,𝑐₃₉=2𝑞₃₈. (29) Let 𝜔²₂=𝑣₂, then (28) becomes 𝑣₂¹⁰+𝑐₃₉𝑣⁹₂+𝑐₃₈𝑣⁸₂+𝑐₃₇𝑣⁷₂+𝑐₃₆𝑣⁶₂+𝑐₃₅𝑣⁵₂+𝑐₃₄𝑣⁴₂+𝑐₃₃𝑣³₂+𝑐₃₂𝑣²₂+𝑐₃₁𝑣₂+𝑐₃₀=0. (30)
Similar as in Case 2, we make the following assumption. (𝐻₃₁) (30) has at least one positive real root. If the condition (𝐻₃₁) holds, then there exists a 𝑣₂₀ such that (23) has a pair of purely imaginary roots ±𝑖𝜔₂₀=±𝑖√𝑣₂₀. For 𝜔₂₀, the corresponding critical value of time delay is 𝜏₂₀=1𝜔₂₀𝑝arccos³⁸𝜔⁸₂₀+𝑝₃₆𝜔⁶₂₀+𝑝₃₄𝜔⁴₂₀+𝑝₃₂𝜔²₂₀+𝑝₃₀𝜔¹⁰₂₀+𝑞₃₈𝜔⁸₂₀+𝑞₃₆𝜔⁶₂₀+𝑞₃₄𝜔⁴₂₀+𝑞₃₂𝜔²₂₀+𝑞₃₀. (31) Differentiating two sides of (23) with respect to 𝜏₂, we have 𝑑𝜆𝑑𝜏₂⁻¹=𝑔₃₁(𝜆)+𝑔₃₂(𝜆)𝑒^𝜆𝜏2+𝑔₃₃(𝜆)𝑒^−𝜆𝜏2𝑔₃₄(𝜆)𝑒^−𝜆𝜏2−𝑔₃₅(𝜆)𝑒^𝜆𝜏2−𝜏₂𝜆, (32) where 𝑔₃₁(𝜆)=4𝐵₃₄𝜆³+3𝐵₃₃𝜆²+2𝐵₃₂𝜆+𝐵₃₁,𝑔₃₂(𝜆)=5𝜆⁴+4𝐴₃₄𝜆³+3𝐴₃₃𝜆²+2𝐴₃₂𝜆+𝐴₃₁,𝑔₃₃(𝜆)=3𝐶₃₃𝜆²+2𝐶₃₂𝜆+𝐶₃₁,𝑔₃₄(𝜆)=𝐶₃₃𝜆⁴+𝐶₃₂𝜆³+𝐶₃₁𝜆²+𝐶₃₀𝑔𝜆,₃₅(𝜆)=𝜆⁶+𝐴₃₄𝜆⁵+𝐴₃₃𝜆⁴+𝐴₃₂𝜆³+𝐴₃₁𝜆²+𝐴₃₀𝜆. (33) Thus, Re𝑑𝜆𝑑𝜏₂⁻¹_{𝜆=𝑖𝜔20}=𝑃_3𝑅𝑄_3𝑅+𝑃_3𝐼𝑄_3𝐼𝑄²_3𝑅+𝑄²_3𝐼, (34) where 𝑃_3𝑅=5𝜔⁴₂₀−3𝐴₃₃+𝐶₃₃5𝜔²₂₀+𝐴₃₁+𝐶₃₁cos𝜏₂₀𝜔₂₀+4𝐴₃₄𝜔³₂₀+2𝐶₃₂−𝐴₃₂𝜔₂₀sin𝜏₂₀𝜔₂₀+𝐵₃₁−3𝐵₃₃𝜔²₂₀,𝑃_3𝐼=5𝜔⁴₂₀−3𝐴₃₃−𝐶₃₃5𝜔²₂₀+𝐴₃₁−𝐶₃₁sin𝜏₂₀𝜔₂₀−4𝐴₃₄𝜔³₂₀−2𝐴₃₂+𝐶₃₂𝜔₂₀cos𝜏₂₀𝜔₂₀+2𝐵₃₂𝜔₂₀−4𝐵₃₄𝜔³₂₀,𝑄_3𝑅=𝐴₃₃+𝐶₃₃𝜔⁴₂₀−𝜔⁶₂₀−𝐴₃₁+𝐶₃₁𝜔²₂₀cos𝜏₂₀𝜔₂₀−𝐴₃₄𝜔⁵₂₀+𝐴₃₂−𝐶₃₂𝜔³₂₀+𝐴₃₀−𝐶₃₀𝜔₂₀×sin𝜏₂₀𝜔₂₀,𝑄_3𝐼=𝐴₃₃−𝐶₃₃𝜔⁴₂₀−𝜔⁶₂₀−𝐴₃₁−𝐶₃₁𝜔²₂₀sin𝜏₂₀𝜔₂₀+𝐴₃₄𝜔⁵₂₀−𝐴₃₂+𝐶₃₂𝜔³₂₀+𝐴₃₀+𝐶₃₀𝜔₂₀×cos𝜏₂₀𝜔₂₀. (35)
Obviously, if the condition (𝐻₃₂)  𝑃_3𝑅𝑄_3𝑅+𝑃_3𝐼𝑄_3𝐼≠0 holds, then Re[𝑑𝜆/𝑑𝜏₂]⁻¹_{𝜆=𝑖𝜔20}≠0. According to the Hopf bifurcation theorem in [19], we have the following results for system (3).

Theorem 2. For system (3), if the conditions (𝐻₃₁)-(𝐻₃₂) hold, then the positive equilibrium 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗) of system (3) is asymptotically stable for 𝜏₂∈[0,𝜏₂₀) and system (3) undergoes a Hopf bifurcation at the positive equilibrium 𝑃^∗(𝑆^∗,𝐸^∗,𝐼^∗,𝑄^∗,𝑅^∗) when 𝜏₂=𝜏₂₀.

Case 4 (𝜏₁>0, 𝜏₂>0, 𝜏₂∈(0,𝜏₂₀)). We consider system (3) under the condition that 𝜏₂ is in its stable interval and 𝜏₁ is a bifurcation parameter.
Let 𝜆=𝑖𝜔_1∗ (𝜔_1∗>0) be the root of (7), then we obtain 𝑀₄₁sin𝜏₁𝜔_1∗+𝑀₄₂cos𝜏₁𝜔_1∗=𝑀₄₃,𝑀₄₁cos𝜏₁𝜔_1∗−𝑀₄₂sin𝜏₁𝜔_1∗=𝑀₄₄, (36) where 𝑀₄₁=𝑛₁𝜔_1∗−𝑛₃𝜔³_1∗+𝑟₁𝜔_1∗−𝑟₃𝜔³_1∗cos𝜏₂𝜔_1∗−𝑟₀−𝑟₂𝜔²_1∗sin𝜏₂𝜔_1∗+𝑠₁𝜔_1∗cos2𝜏₂𝜔_1∗−𝑠₀−𝑠₂𝜔²_1∗sin2𝜏₂𝜔_1∗,𝑀₄₂=𝑛₄𝜔⁴_1∗−𝑛₂𝜔²_1∗+𝑛₀+𝑟₁𝜔_1∗−𝑟₃𝜔³_1∗sin𝜏₂𝜔_1∗+𝑟₀−𝑟₂𝜔²_1∗cos𝜏₂𝜔_1∗+𝑠₁𝜔_1∗sin2𝜏₂𝜔_1∗+𝑠₀−𝑠₂𝜔²_1∗cos2𝜏₂𝜔_1∗,𝑀₄₃=𝑝₃𝜔³_1∗−𝑝₁𝜔_1∗sin𝜏₂𝜔_1∗−𝑝₄𝜔⁴_1∗−𝑝₂𝜔²_1∗+𝑝₀cos𝜏₂𝜔_1∗+𝑞₃𝜔³_1∗−𝑞₁𝜔_1∗sin2𝜏₂𝜔_1∗+𝑞₂𝜔²_1∗−𝑞₀cos𝜏₂𝜔_1∗+𝑚₂𝜔²_1∗−𝑚₄𝜔⁴_1∗−𝑚₀,𝑀₄₄=𝑝₃𝜔³_1∗−𝑝₁𝜔_1∗cos𝜏₂𝜔_1∗+𝑝₄𝜔⁴_1∗−𝑝₂𝜔²_1∗+𝑝₀sin𝜏₂𝜔_1∗+𝑞₃𝜔³_1∗−𝑞₁𝜔_1∗cos2𝜏₂𝜔_1∗−𝑞₂𝜔²_1∗−𝑞₀sin𝜏₂𝜔_1∗+𝑚₃𝜔³_1∗−𝜔⁵_1∗−𝑚₁𝜔_1∗. (37) Then, we can obtain 𝑓₄₀𝜔_1∗+2𝑓₄₁𝜔_1∗cos𝜏₂𝜔_1∗+2𝑓₄₂𝜔_1∗sin𝜏₂𝜔_1∗+2𝑓₄₃𝜔_1∗cos2𝜏₂𝜔_1∗+2𝑓₄₄𝜔_1∗sin2𝜏₂𝜔_1∗+2𝑓₄₅𝜔_1∗cos𝜏₂𝜔_1∗cos2𝜏₂𝜔_1∗+2𝑓₄₆𝜔_1∗cos𝜏₂𝜔_1∗sin2𝜏₂𝜔_1∗+2𝑓₄₇𝜔_1∗sin𝜏₂𝜔_1∗sin2𝜏₂𝜔_1∗+2𝑓₄₈𝜔_1∗sin𝜏₂𝜔_1∗cos2𝜏₂𝜔_1∗=0, (38) where 𝑓₄₀𝜔_1∗=𝜔¹⁰_1∗+𝑚²₄−𝑛²₄+𝑝²₄−2𝑚₃𝜔⁸_1∗+𝑚²₃−𝑛²₃+𝑝²₃+𝑞²₃−𝑟²₃+2𝑚₁−2𝑚₂𝑚₄+2𝑛₂𝑛₄−2𝑝₂𝑝₄𝜔⁶_1∗+𝑚²₂−𝑛²₂+𝑝²₂+𝑞²₂−𝑟²₂−𝑠²₂+2𝑚₀𝑚₄−2𝑚₁𝑚₃+2𝑛₀𝑛₄+2𝑛₁𝑛₃−2𝑝₁𝑝₃−2𝑞₁𝑞₃+2𝑟₁𝑟₃𝜔⁴_1∗+𝑚²₁−𝑛²₁+𝑝²₁+𝑞²₁−𝑟²₁−𝑠²₁−2𝑚₀𝑚₂+2𝑛₀𝑛₂−2𝑝₀𝑝₂−2𝑞₀𝑞₂+2𝑟₀𝑟₂+2𝑠₀𝑠₂𝜔²_1∗+𝑚²₀−𝑛²₀+𝑝²₀+𝑞²₀−𝑟²₀−𝑠²₀,𝑓₄₁𝜔_1∗=𝑚₄𝑝₄−𝑝₃𝜔⁸_1∗+𝑚₃𝑝₃+𝑝₁−𝑚₂𝑝₄−𝑚₄𝑝₂−𝑛₃𝑟₃+𝑛₄𝑟₂𝜔⁶_1∗+𝑚₀𝑝₄−𝑚₁𝑝₃+𝑚₂𝑝₂−𝑚₃𝑝₁+𝑚₄𝑝₀+𝑛₁𝑟₃−𝑛₂𝑟₂+𝑛₃𝑟₁−𝑛₄𝑟₀𝜔⁴_1∗+𝑚₁𝑝₁−𝑚₀𝑝₂−𝑚₂𝑝₀+𝑛₀𝑟₂−𝑛₁𝑟₁+𝑛₂𝑟₀𝜔²_1∗+𝑚₀𝑝₀−𝑛₀𝑟₀,𝑓₄₂𝜔_1∗=−𝑝₄𝜔⁹_1∗+𝑚₃𝑝₄−𝑚₄𝑝₃+𝑛₄𝑟₃+𝑝₂𝜔⁷_1∗+𝑚₂𝑝₃−𝑚₁𝑝₄−𝑚₃𝑝₂+𝑚₄𝑝₁−𝑛₂𝑟₃+𝑛₃𝑟₂−𝑛₄𝑟₁−𝑝₀𝜔⁵_1∗+𝑚₁𝑝₂−𝑚₀𝑝₃−𝑚₂𝑝₁+𝑚₃𝑝₀+𝑛₀𝑟₃−𝑛₁𝑟₂+𝑛₂𝑟₁−𝑛₃𝑟₀𝜔³_1∗+𝑚₀𝑝₁−𝑚₁𝑝₀−𝑛₀𝑟₁+𝑛₁𝑟₀𝜔_1∗,𝑓₄₃𝜔_1∗=−𝑞₃𝜔⁸_1∗+𝑚₃𝑞₃−𝑚₄𝑞₂+𝑞₁𝜔⁶_1∗+𝑚₂𝑞₂−𝑚₁𝑞₃−𝑚₃𝑞₁+𝑚₄𝑞₀+𝑛₃𝑠₁𝜔⁴_1∗+𝑚₁𝑞₁−𝑚₀𝑞₂−𝑚₂𝑞₀−𝑛₁𝑠₁𝜔²_1∗,𝑓₄₄𝜔_1∗=𝑞₂𝜔⁷_1∗+𝑛₃𝑠₂−𝑚₃𝑞₂−𝑞₀𝜔⁵_1∗+𝑚₁𝑞₂+𝑚₃𝑞₀−𝑛₁𝑠₂−𝑛₃𝑠₀𝜔³_1∗+𝑛₁𝑠₀−𝑚₁𝑞₀𝜔_1∗,𝑓₄₅𝜔_1∗=𝑝₃𝑞₃−𝑝₄𝑞₂𝜔⁶_1∗+𝑝₂𝑞₂+𝑝₄𝑞₀−𝑝₁𝑞₃−𝑝₃𝑞₁+𝑟₂𝑠₂−𝑟₃𝑠₁𝜔⁴_1∗+𝑝₁𝑞₁−𝑝₀𝑞₂−𝑝₂𝑞₀−𝑟₀𝑠₂+𝑟₁𝑠₁−𝑟₂𝑠₀𝜔²_1∗+𝑝₀𝑞₀+𝑟₀𝑠₀,𝑓₄₆𝜔_1∗=−𝑝₄𝑞₃𝜔⁷_1∗+𝑝₂𝑞₃−𝑝₃𝑞₂+𝑝₄𝑞₁+𝑟₃𝑠₂𝜔⁵_1∗+𝑝₁𝑞₂−𝑝₀𝑞₃−𝑝₂𝑞₁−𝑝₃𝑞₀−𝑟₁𝑠₂+𝑟₂𝑠₁−𝑟₃𝑠₀𝜔³_1∗+𝑝₀𝑞₁−𝑝₁𝑞₀−𝑟₁𝑠₀−𝑟₀𝑠₁𝜔_1∗,𝑓₄₇𝜔_1∗=𝑝₃𝑞₃−𝑝₄𝑞₂𝜔⁶_1∗+𝑝₂𝑞₂−𝑝₁𝑞₃−𝑝₃𝑞₁−𝑝₄𝑞₀−𝑟₂𝑠₂+𝑟₃𝑠₁𝜔⁴_1∗+𝑝₁𝑞₁−𝑝₀𝑞₂−𝑝₂𝑞₀+𝑟₀𝑠₂−𝑟₁𝑠₁+𝑟₂𝑠₀𝜔²_1∗+𝑝₀𝑞₀−𝑟₀𝑠₀,𝑓₄₈𝜔_1∗=𝑝₄𝑞₃𝜔⁷_1∗+𝑝₃𝑞₂−𝑝₂𝑞₃−𝑝₄𝑞₁−𝑟₃𝑠₂𝜔⁵_1∗+𝑝₀𝑞₃−𝑝₁𝑞₂+𝑝₂𝑞₁−𝑝₃𝑞₀+𝑟₁𝑠₂−𝑟₂𝑠₁−𝑟₃𝑠₀𝜔³_1∗+𝑝₁𝑞₀+𝑝₀𝑞₁+𝑟₀𝑠₁−𝑟₁𝑠₀𝜔_1∗. (39) In order to give the main results in this paper, we make the following assumption.
(𝐻₄₁) (38) has at least one positive real root. If the conditions (𝐻₄₁) hold, then there exists a 𝜔^∗₁₀ such that (7) has a pair of purely imaginary roots ±𝑖𝜔^∗₁₀. For 𝜔^∗₁₀, the corresponding critical value of time delay is 𝜏^∗₁₀=1𝜔^∗₁₀𝑀arccos⁴¹𝑀₄₄+𝑀₄₂𝑀₄₃𝑀²₄₁+𝑀²_{42𝜔1∗=𝜔∗10}. (40) Differentiating two sides of (7) with respect to 𝜏₁, we have 𝑑𝜆𝑑𝜏₁⁻¹=𝑔₄₁(𝜆)+𝑔₄₂(𝜆)𝑒^−𝜆𝜏1+𝑔₄₃(𝜆)𝑒^−𝜆𝜏2+𝑔₄₄(𝜆)𝑒^{−𝜆(𝜏1+𝜏2)}+𝑔₄₅(𝜆)𝑒^{−2𝜆𝜏2}+𝑔₄₆(𝜆)𝑒^{−𝜆(𝜏1+2𝜏2)}×𝑔₄₇(𝜆)𝑒^−𝜆𝜏1+𝑔₄₈(𝜆)𝑒^{−𝜆(𝜏1+𝜏2)}+𝑔₄₉(𝜆)𝑒^{−𝜆(𝜏1+2𝜏2)−1}−𝜏₁𝜆, (41) where 𝑔₄₁(𝜆)=5𝜆⁴+4𝑚₄𝜆³+3𝑚₃𝜆²+2𝑚₂𝜆+𝑚₁,𝑔₄₂(𝜆)=4𝑛₄𝜆³+3𝑛₃𝜆²+2𝑛₂𝜆+𝑛₁,𝑔₄₃(𝜆)=−𝜏₂𝑝₄𝜆⁴+4𝑝₄−𝜏₂𝑝₃𝜆³+3𝑝₃−𝜏₂𝑝₂𝜆²+2𝑝₂+𝜏₂𝑝₁𝜆+𝑝₁−𝜏₂𝑝₀,𝑔₄₄(𝜆)=−𝜏₂𝑞³₃+3𝑞₃−𝜏₂𝑞₂²+2𝑞₂−𝜏₂𝑝₁𝜆+𝑞₁−𝜏₂𝑞₀,𝑔₄₅(𝜆)=−𝜏₂𝑠₂𝜆²+2𝑠₂−𝜏₂𝑠₁𝜆+𝑠₁−𝜏₂𝑠₀,𝑔₄₆(𝜆)=3𝑟₃𝜆²+2𝑟₂𝜆+𝑟₁,𝑔₄₇(𝜆)=𝑛₄𝜆⁵+𝑛₃𝜆⁴+𝑛₂𝜆³+𝑛₁𝜆²+𝑛₀𝑔𝜆,₄₈(𝜆)=𝑟₃𝜆⁴+𝑟₂𝜆³+𝑟₁𝜆²+𝑟₀𝑔𝜆,₄₉(𝜆)=𝑠₂𝜆³+𝑠₁𝜆²+𝑠₀𝜆. (42) Thus, Re𝑑𝜆𝑑𝜏₁⁻¹_{𝜆=𝑖𝜔∗10}=𝑃_4𝑅𝑄_4𝑅+𝑃_4𝐼𝑄_4𝐼𝑄²_4𝑅+𝑄²_4𝐼, (43) where 𝑃_4𝑅=5𝜔^∗₁₀⁴−3𝑚₃𝜔^∗₁₀²+𝑚₁+2𝑛₂𝜔^∗₁₀−4𝑛₄𝜔^∗₁₀³+2𝑟₂𝜔^∗₁₀cos𝜏₂𝜔^∗₁₀+3𝑟₃𝜔^∗₁₀²−𝑟₁sin𝜏₂𝜔^∗₁₀+2𝑠₂−𝜏₂𝑠₁𝜔^∗₁₀cos2𝜏^∗₂𝜔^∗₁₀−𝜏₂𝑠₂𝜔^∗₁₀²+𝑠₁−𝜏₂𝑠₀sin2𝜏₂𝜔^∗₁₀sin𝜏^∗₁₀𝜔^∗₁₀+𝑛₁−3𝑛₃𝜔^∗₁₀²+2𝑟₂𝜔^∗₁₀sin𝜏₂𝜔^∗₁₀+𝑟₁−3𝑟₃𝜔^∗₁₀²cos𝜏₂𝜔^∗₁₀+2𝑠₂−𝜏₂𝑠₁𝜔^∗₁₀sin2𝜏₂𝜔^∗₁₀+𝜏₂𝑠₂𝜔^∗₁₀²+𝑠₁−𝜏₂𝑠₀cos2𝜏₂𝜔^∗₁₀cos𝜏₂𝜔^∗₁₀+2𝑝₂−𝜏₂𝑝₁𝜔^∗₁₀−4𝑝₄−𝜏₂𝑝₃𝜔^∗₁₀³sin𝜏₂𝜔^∗₁₀+𝜏₂𝑝₂−3𝑝₃𝜔^∗₁₀²−𝜏₂𝑝₄𝜔^∗₁₀⁴+𝑝₁−𝜏₂𝑝₀×cos𝜏₂𝜔^∗₁₀+𝜏₂𝑞₃𝜔^∗₁₀³+2𝑞₂−𝜏₂𝑞₁𝜔^∗₁₀sin2𝜏₂𝜔^∗₁₀+𝜏₂𝑞₂−3𝑞₃𝜔^∗₁₀²+𝑞₁−𝜏₂𝑞₀cos2𝜏₂𝜔^∗₁₀,𝑃_4𝐼=2𝑚₂𝜔^∗₁₀−4𝑚₄𝜔^∗₁₀³+2𝑟₂𝜔^∗₁₀cos𝜏₂𝜔^∗₁₀+3𝑟₃𝜔^∗₁₀²−𝑟₁sin𝜏₂𝜔^∗₁₀+2𝑠₂−𝜏₂𝑠₁𝜔^∗₁₀cos2𝜏₂𝜔^∗₁₀−𝜏₂𝑠₂𝜔^∗₁₀²+𝑠₁−𝜏₂𝑠₀sin2𝜏₂𝜔^∗₁₀cos𝜏₂𝜔^∗₁₀+3𝑛₃𝜔^∗₁₀²−𝑛₁−2𝑟₂𝜔^∗₁₀sin𝜏₂𝜔^∗₁₀−𝑟₁−3𝑟₃𝜔^∗₁₀²cos𝜏₂𝜔^∗₁₀−2𝑠₂−𝜏₂𝑠₁𝜔^∗₁₀sin2𝜏₂𝜔^∗₁₀−𝜏₂𝑠₂𝜔^∗₁₀²+𝑠₁−𝜏₂𝑠₀cos2𝜏₂𝜔^∗₁₀sin𝜏^∗₁₀𝜔^∗₁₀+2𝑝₂−𝜏₂𝑝₁𝜔^∗₁₀−4𝑝₄−𝜏₂𝑝₃𝜔^∗₁₀³cos𝜏₂𝜔^∗₁₀−𝜏₂𝑝₂−3𝑝₃𝜔^∗₁₀²−𝜏₂𝑝₄𝜔^∗₁₀⁴+𝑝₁−𝜏₂𝑝₀×sin𝜏₂𝜔^∗₁₀+𝜏₂𝑞₃𝜔^∗₁₀³+2𝑞₂−𝜏₂𝑞₁𝜔^∗₁₀cos2𝜏₂𝜔^∗₁₀−𝜏₂𝑞₀−3𝑞₃𝜔^∗₁₀²+𝑞₁−𝜏₂𝑞₀sin2𝜏₂𝜔^∗₁₀,𝑄_4𝑅=𝑛₄𝜔^∗₁₀⁵−𝑛₂𝜔^∗₁₀³+𝑛₀𝜔^∗₁₀+𝑟₀𝜔^∗₁₀−𝑟₂𝜔^∗₁₀³cos𝜏₂𝜔^∗₁₀−𝑟₃𝜔^∗₁₀⁴−𝑟₁𝜔^∗₁₀²sin𝜏₂𝜔^∗₁₀+𝑠₀𝜔^∗₁₀−𝑠₂𝜔^∗₁₀³cos𝜏₂𝜔^∗₁₀+𝑠₁𝜔^∗₁₀²sin2𝜏₂𝜔^∗₁₀sin𝜏^∗₁₀𝜔^∗₁₀+𝑛₃𝜔^∗₁₀⁴−𝑛₁𝜔^∗₁₀²+𝑟₀𝜔^∗₁₀−𝑟₂𝜔^∗₁₀³sin𝜏₂𝜔^∗₁₀+𝑟₃𝜔^∗₁₀⁴−𝑟₁𝜔^∗₁₀²cos𝜏₂𝜔^∗₁₀+𝑠₀𝜔^∗₁₀−𝑠₂𝜔^∗₁₀³sin2𝜏₂𝜔^∗₁₀−𝑠₁𝜔^∗₁₀²cos2𝜏₂𝜔^∗₁₀cos𝜏^∗₁₀𝜔^∗₁₀,𝑄_4𝐼=𝑛₄𝜔^∗₁₀⁵−𝑛₂𝜔^∗₁₀