Review Special Issues

Protein chainmail variants in dsDNA viruses

  • First discovered in bacteriophage HK97, biological chainmail is a highly stable system formed by concatenated protein rings. Each subunit of the ring contains the HK97-like fold, which is characterized by its submarine-like shape with a 5-stranded β sheet in the axial (A) domain, spine helix in the peripheral (P) domain, and an extended (E) loop. HK97 capsid consists of covalently-linked copies of just one HK97-like fold protein and represents the most effective strategy to form highly stable chainmail needed for dsDNA genome encapsidation. Recently, near-atomic resolution structures enabled by cryo electron microscopy (cryoEM) have revealed a range of other, more complex variants of this strategy for constructing dsDNA viruses. The first strategy, exemplified by P22-like phages, is the attachment of an insertional (I) domain to the core 5-stranded β sheet of the HK97-like fold. The atomic models of the Bordetella phage BPP-1 showcases an alternative topology of the classic HK97 topology of the HK97-like fold, as well as the second strategy for constructing stable capsids, where an auxiliary jellyroll protein dimer serves to cement the non-covalent chainmail formed by capsid protein subunits. The third strategy, found in lambda-like phages, uses auxiliary protein trimers to stabilize the underlying non-covalent chainmail near the 3-fold axis. Herpesviruses represent highly complex viruses that use a combination of these strategies, resulting in four-level hierarchical organization including a non-covalent chainmail formed by the HK97-like fold domain found in the floor region. A thorough understanding of these structures should help unlock the enigma of the emergence and evolution of dsDNA viruses and inform bioengineering efforts based on these viruses.

    Citation: Z. Hong Zhou, Joshua Chiou. Protein chainmail variants in dsDNA viruses[J]. AIMS Biophysics, 2015, 2(2): 200-218. doi: 10.3934/biophy.2015.2.200

    Related Papers:

    [1] Cruz Vargas-De-León, Alberto d'Onofrio . Global stability of infectious disease models with contact rate as a function of prevalence index. Mathematical Biosciences and Engineering, 2017, 14(4): 1019-1033. doi: 10.3934/mbe.2017053
    [2] Yoichi Enatsu, Yukihiko Nakata, Yoshiaki Muroya . Global stability for a class of discrete SIR epidemic models. Mathematical Biosciences and Engineering, 2010, 7(2): 347-361. doi: 10.3934/mbe.2010.7.347
    [3] Jinliang Wang, Hongying Shu . Global analysis on a class of multi-group SEIR model with latency and relapse. Mathematical Biosciences and Engineering, 2016, 13(1): 209-225. doi: 10.3934/mbe.2016.13.209
    [4] Rongjian Lv, Hua Li, Qiubai Sun, Bowen Li . Model of strategy control for delayed panic spread in emergencies. Mathematical Biosciences and Engineering, 2024, 21(1): 75-95. doi: 10.3934/mbe.2024004
    [5] Yu Ji . Global stability of a multiple delayed viral infection model with general incidence rate and an application to HIV infection. Mathematical Biosciences and Engineering, 2015, 12(3): 525-536. doi: 10.3934/mbe.2015.12.525
    [6] Abdennasser Chekroun, Mohammed Nor Frioui, Toshikazu Kuniya, Tarik Mohammed Touaoula . Global stability of an age-structured epidemic model with general Lyapunov functional. Mathematical Biosciences and Engineering, 2019, 16(3): 1525-1553. doi: 10.3934/mbe.2019073
    [7] Gonzalo Robledo . Feedback stabilization for a chemostat with delayed output. Mathematical Biosciences and Engineering, 2009, 6(3): 629-647. doi: 10.3934/mbe.2009.6.629
    [8] Pan Yang, Jianwen Feng, Xinchu Fu . Cluster collective behaviors via feedback pinning control induced by epidemic spread in a patchy population with dispersal. Mathematical Biosciences and Engineering, 2020, 17(5): 4718-4746. doi: 10.3934/mbe.2020259
    [9] Dongmei Li, Bing Chai, Weihua Liu, Panpan Wen, Ruixue Zhang . Qualitative analysis of a class of SISM epidemic model influenced by media publicity. Mathematical Biosciences and Engineering, 2020, 17(5): 5727-5751. doi: 10.3934/mbe.2020308
    [10] Jianquan Li, Zhien Ma, Fred Brauer . Global analysis of discrete-time SI and SIS epidemic models. Mathematical Biosciences and Engineering, 2007, 4(4): 699-710. doi: 10.3934/mbe.2007.4.699
  • First discovered in bacteriophage HK97, biological chainmail is a highly stable system formed by concatenated protein rings. Each subunit of the ring contains the HK97-like fold, which is characterized by its submarine-like shape with a 5-stranded β sheet in the axial (A) domain, spine helix in the peripheral (P) domain, and an extended (E) loop. HK97 capsid consists of covalently-linked copies of just one HK97-like fold protein and represents the most effective strategy to form highly stable chainmail needed for dsDNA genome encapsidation. Recently, near-atomic resolution structures enabled by cryo electron microscopy (cryoEM) have revealed a range of other, more complex variants of this strategy for constructing dsDNA viruses. The first strategy, exemplified by P22-like phages, is the attachment of an insertional (I) domain to the core 5-stranded β sheet of the HK97-like fold. The atomic models of the Bordetella phage BPP-1 showcases an alternative topology of the classic HK97 topology of the HK97-like fold, as well as the second strategy for constructing stable capsids, where an auxiliary jellyroll protein dimer serves to cement the non-covalent chainmail formed by capsid protein subunits. The third strategy, found in lambda-like phages, uses auxiliary protein trimers to stabilize the underlying non-covalent chainmail near the 3-fold axis. Herpesviruses represent highly complex viruses that use a combination of these strategies, resulting in four-level hierarchical organization including a non-covalent chainmail formed by the HK97-like fold domain found in the floor region. A thorough understanding of these structures should help unlock the enigma of the emergence and evolution of dsDNA viruses and inform bioengineering efforts based on these viruses.


    Infectious diseases (such as influenza, malaria, cholera, tuberculosis, hepatitis, AIDS, etc) have always seriously threatened humans' life and health. With in-depth understanding of infectious diseases, scientists have been continuing to explore effective methods to prevent and control the outbreaks of various infectious diseases. It is well known that mathematical models have played very important roles in analysis of control strategies for disease transmission [1,2,3,4,5,6,7,8,9,10].

    When studying the long-term evolutionary behavior of an ecological system, as pointed in [11], the equilibrium of biological system may not be the desirable one, and smaller value is required. This can be achieved by introducing suitable feedback control variable. The feedback control mechanism might be implemented through harvesting or culling procedures or certain biological control schemes [12]. In addition, in a control system, the time delay factor generally exists in the signal transmission process. Thus, feedback control with coupled time delay may have better biological significance [12] and has been extensively introduced into some important population ecological systems (see, for example, [13,14,15,16] and the references cited therein).

    In recent years, feedback control has also been successfully applied to some infectious disease dynamical systems. For example, in [17], the authors considered the following SI epidemic model with two feedback control variables:

    $ {˙S(t)=S(t)(raS(t)bI(t)c1u1(t)),˙I(t)=I(t)(bS(t)μfI(t)c2u2(t)),˙u1(t)=e1u1(t)+d1S(t),˙u2(t)=e2u2(t)+d2I(t).
    $
    (1.1)

    In model (1.1), the state variables $ S(t) $ and $ I(t) $ represent the numbers of susceptibles and infectives at time $ t $, respectively; $ u_{1}(t) $ and $ u_{2}(t) $ are feedback control variables. The number of susceptibles grows according to the regulation of a logistic curve with the capacity $ r/a $ ($ r > 0 $, $ a > 0 $) and a constant recruitment rate $ r $; the constant $ b > 0 $ is the transmission rate when susceptibles contact with infectives; the constants $ \mu > 0 $ and $ f > 0 $ are the death rates of the infectives with respect to single and mutiple of infectives, respectively; the constants $ c_1 > 0 $, $ c_2 > 0 $, $ d_1 > 0 $, $ d_2 > 0 $, $ e_1 > 0 $ and $ e_2 > 0 $ are the feedback control parameters. By constructing suitable Lyapunov functions, the authors established threshold dynamics of model (1.1) completely determined by the threshold parameter $ \gamma_{0} = (br-a\mu)e_1/(c_1d_1\mu) $. The results in [17] indicate that, by appropriately choosing feedback control parameters, it can make the disease infection endemic or extinct. In [18], the author considered a two-group SI epidemic model with feedback control only in the susceptible individuals, and showed that the disease outbreaks can be controlled by adjusting feedback control parameters. In addition, in [19], the authors further extend model (1.1) to the case of patchy environment.

    Since the authors of [20] have introduced a nonlinear incidence rate $ g(I)S $ into classic Kermack-McKendrick SIR model, nonlinear incidence rate has been further introduced into more general SIR/SIRS epidemic models with time delays or infection age etc (see, for example, [21,22,23,24,25,26,27] and the references cited therein). Usually, the function $ g(I) $ takes the following two types: (ⅰ) saturated, such as $ g(I) = bI/(1+kI) $, or $ g(I) = bI^2/(1+kI^2) $; (ⅱ) unimodal, such as $ g(I) = bI/(1+kI^2) $, here $ k > 0 $ is constant. In biology, $ bI $ or $ bI^2 $ measures the infection force of the disease, $ 1/(1+kI) $ or $ 1/(1+kI^2) $ measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infective individuals [21,22].

    Recently, in [28], the authors further extended model (1.1) to the following more general case with the saturated incidence rate $ bSI/(1+kI) $ and feedback controls:

    $ {˙S(t)=S(t)(raS(t)bI(t)1+kI(t)c1u1(t)),˙I(t)=I(t)(bS(t)1+kI(t)μfI(t)c2u2(t)),˙u1(t)=e1u1(t)+d1S(t),˙u2(t)=e2u2(t)+d2I(t),
    $
    (1.2)

    and some sufficient conditions for global asymptotic stability of the disease-free equilibrium and the endemic equilibrium of model (1.2) are established by the method of Lyapunov functions. In addition, the authors also considered permanence and existence of almost periodic solutions for a class of non-autonomous system based on model (1.2).

    Motivated by the above works and model (1.2), we further consider the following SI epidemic model with saturated incidence rate, two feedback control variables and four time delays:

    $ {˙S(t)=S(t)(raS(t)bI(t)1+kI(t)c1u1(tτ1)),˙I(t)=I(t)(bS(t)1+kI(t)μfI(t)c2u2(tτ2)),˙u1(t)=e1u1(t)+d1S(tτ3),˙u2(t)=e2u2(t)+d2I(tτ4).
    $
    (1.3)

    The biological significance of all the parameters of model (1.3) are the same as in model (1.2) except time delays $ \tau_{i}\geq0 $ ($ i = 1, 2, 3, 4 $). In model (1.3), $ u_1(t) $ and $ u_2(t) $ are introduced as control variables. Usually, there always exist time delays in the transmission of information. Therefore, $ \tau_1 $ and $ \tau_2 $ can be understood as the result of transmission of information, and while $ \tau_3 $ and $ \tau_4 $ represent usual feedback control delays.

    The purpose in this paper focuses on global dynamics of the equilibria of model (1.3) by constructing appropriate Lyapunov functionals, and our results further extend and improve works in [17,28].

    The rest of the paper is organized as follows. In Section 2, we provide some preliminary results, including the well-posedness and dissipativeness of the solutions of model (1.3), the expression of the basic reproduction number and the classification of the equilibria of model (1.3). In Sections 3 and 4, we establish some sufficient conditions for global asymptotic stability and global attractivity of the disease-free equilibrium and the endemic equilibrium of model (1.3), which are the main results of this paper. In the last section, the conclusions and some numerical simulations are given.

    Let $ C^{+} = C([-\tau, 0], \mathbb{R}_{+}^{4}) $ be the Banach space of continuous functions mapping the interval $ [-\tau, 0] $ into $ \mathbb{R}_{+}^{4} $ equipped with the supremum norm, where $ \tau = \max\{\tau_1, \tau_2, \tau_3, \tau_4\} $. The initial condition of model (1.3) is given as follows,

    $ S(θ)=ϕ1(θ),I(θ)=ϕ2(θ),u1(θ)=ϕ3(θ),u2(θ)=ϕ4(θ),θ[τ,0],
    $
    (2.1)

    where $ \phi = (\phi_{1}(\theta), \phi_{2}(\theta), \phi_{3}(\theta), \phi_{4}(\theta))\in C^{+} $.

    By using the standard theory of delay differential equations (DDEs) (see, for example, [29,30,31]), we can easily establish the following result.

    Theorem 2.1. The solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3) with the initial condition (2.1) is existent, unique and nonnegative on $ [0, \infty) $, and satisfies

    $ lim supt+S(t)ra,lim supt+I(t)rbaf,lim supt+u1(t)rd1ae1,lim supt+u2(t)rbd2afe2.
    $
    (2.2)

    Moreover, the following bounded set

    $ \Omega: = \left\{\phi\in C^{+}: \; \|\phi_{1}\|\leq\frac{r}{a},\; \|\phi_{2}\|\leq\frac{rb}{af},\; \|\phi_{3}\|\leq\frac{rd_{1}}{ae_{1}},\; \|\phi_{4}\|\leq \frac{rbd_{2}}{afe_{2}} \right\} $

    is positively invariant with respect to model (1.3).

    Proof. It is not difficult to show that the solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3) with the initial condition (2.1) is existent, unique and nonnegative on $ [0, \infty) $. Let us consider ultimate boundedness of model (1.3). According to the first equation of model (1.3), we have that for $ t\geq 0 $,

    $ ˙S(t)S(t)(raS(t)),
    $
    (2.3)

    which implies $ \limsup_{t\rightarrow+\infty}S(t)\leq \frac{r}{a} $. For any sufficiently small $ \varepsilon > 0 $, there exists a $ \widehat{t} > 0 $ such that $ S(t) < \frac{r}{a}+\varepsilon $ for $ t\geq \widehat{t} $. Further, according to the second equation of model (1.3), we have for $ t\geq \widehat{t} $,

    $ { \dot{I}(t)\leq I(t)\left[b\left(\frac{r}{a}+\varepsilon\right)-fI(t)\right], } $

    which implies

    $ lim supt+I(t)rbaf+bfε.
    $
    (2.4)

    Since inequality (2.4) holds for arbitrary $ \varepsilon > 0 $, we obtain $ \limsup_{t\rightarrow+\infty}I(t)\leq \frac{rb}{af} $. Similarly, according to the last two equations of model (1.3), we can obtain $ \limsup_{t\rightarrow+\infty}u_{1}(t)\leq\frac{rd_{1}}{ae_{1}} $, $ \limsup_{t\rightarrow+\infty}u_{2}(t)\leq\frac{rbd_{2}}{afe_{2}} $.

    Let $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ be the solution of model (1.3) with the initial function $ \phi = (\phi_{1}, \phi_{2}, \phi_{3}, \phi_{4})\in \Omega $. { For $ t\geq0 $, we have $ \dot{S}(t)\leq S(t)\left(r-aS(t)\right) $, which implies that for $ t\geq0 $,

    $ S(t)raϕ1(0)ϕ1(0)+[raϕ1(0)]ertra,
    $

    where $ \phi_{1}(0)\leq \frac{r}{a} $ is used. Further combining the second equation of model (1.3), for $ t\geq0 $, we have $ \dot{I}(t)\leq I(t)(\frac{rb}{a}-fI(t)) $, which implies that for $ t\geq0 $,

    $ I(t)rbafϕ2(0)ϕ2(0)+[rbafϕ2(0)]erbatrbaf,
    $

    where $ \phi_{2}(0)\leq \frac{rb}{af} $ is used. Thus, for $ t\geq0 $, we have $ \dot{u}_{1}(t)\leq \frac{rd_{1}}{a}-e_{1}u_{1}(t), \; \; \dot{u}_{2}(t) \leq \frac{rbd_{2}}{af}-e_{2}u_{2}(t) $. This implies that for $ t\geq0 $,

    $ u1(t)rd1ae1+[ϕ3(0)rd1ae1]ee1trd1ae1,u2(t)rbd2afe2+[ϕ4(0)rbd2afe2]ee2trbd2afe2,
    $

    where $ \phi_{3}(0)\leq \frac{rd_{1}}{ae_{1}} $ and $ \phi_{4}(0)\leq \frac{rbd_{2}}{afe_{2}} $ are used. Hence, it has that $ \Omega $ is positively invariant with respect to model (1.3).

    The proof is completed.

    Obviously, model (1.3) always has a trivial equilibrium $ \widetilde{E} = (0, 0, 0, 0) $ and a disease-free equilibrium $ E_{0} = (S^{0}, 0, u_{1}^{0}, 0) $, where

    $ S^{0} = \frac{re_{1}}{ae_{1}+d_{1}c_{1}}, \; \; u_{1}^{0} = \frac{rd_{1}}{ae_{1}+d_{1}c_{1}}. $

    Then, by the methods in [32,33], we can derive the expression of the basic reproduction number of model (1.3) as follows. First we define matrices $ \mathbb{F} $ and $ \mathbb{V} $ as

    $ \mathbb{F} = \left( bS0000
    \right),\; \; \; \mathbb{V} = \left( μ0d2e2
    \right). $

    Then the basic reproduction number $ R_0 $ is defined as the spectral radius of $ \mathbb{F}\mathbb{V}^{-1} $. Therefore,

    $ R_{0}: = \rho(\mathbb{F}\mathbb{V}^{-1}) = \frac{bS^{0}}{\mu} = \frac{bre_{1}}{\mu(ae_{1}+d_{1}c_{1})}. $

    Suppose $ (S, I, u_{1}, u_{2}) $ is an endemic equilibrium (positive equilibrium) of model (1.3), where $ S > 0 $, $ I > 0 $, $ u_{1} > 0 $, $ u_{2} > 0 $ satisfy the following equations

    $ raSbI1+kIc1u1=0,bS1+kIμfIc2u2=0,e1u1+d1S=0,e2u2+d2I=0.
    $
    (2.5)

    From Eq (2.5), it is not difficult to obtain the following relationships

    $ u2=d2e2I,u1=d1e1S,S=e1ae1+d1c1(rbI1+kI).
    $
    (2.6)

    Through Eq (2.6) and combining the second equation of Eq (2.5), we can obtain that $ I $ satisfies the following equation,

    $ F(I)\equiv \frac{be_{1}}{(ae_{1}+d_{1}c_{1})(1+kI)} \left(r-\frac{bI}{1+kI}\right)-\mu-\left(f+\frac{d_{2}c_{2}}{e_{2}}\right)I = 0. $

    According to Eq (2.6), in order to ensure that $ S > 0 $, we need to consider the following two cases:

    (ⅰ) $ rk < b $, $ 0 < I < \frac{r}{b-rk}\equiv \widetilde{I} $;

    (ⅱ) $ rk\geq b $, $ I > 0 $.

    Clearly, for both case (ⅰ) and case (ⅱ), we have that

    $ \dot{F}(I) = -\frac{be_{1}k}{(ae_{1}+d_{1}c_{1})(1+kI)^{2}} \left(r-\frac{bI}{1+kI}\right) -\frac{b^{2}e_{1}}{(ae_{1}+d_{1}c_{1})(1+kI)^{3}} -\left(f+\frac{d_{2}c_{2}}{e_{2}}\right) \lt 0. $

    Hence, $ F(I) $ is monotonically decreasing with respect to $ I $ and

    $ \lim\limits_{I\rightarrow 0^{+}}F(I) = \frac{bre_{1}}{ae_{1}+d_{1}c_{1}}-\mu = \mu(R_{0}-1). $

    If $ R_{0}\leq1 $, then $ \lim_{I\rightarrow 0^{+}}F(I)\leq0 $ and $ F(I) = 0 $ has no positive roots. If $ R_{0} > 1 $, then $ \lim_{I\rightarrow 0^{+}}F(I) > 0 $.

    For case (ⅰ), we have that

    $ \lim\limits_{I\rightarrow \widetilde{I}^{-}}F(I) = -\mu-\left(f+\frac{d_{2}c_{2}}{e_{2}}\right)\widetilde{I} \lt 0. $

    For case (ⅱ), we have that

    $ F\left(\frac{rbe_{1}e_{2}}{(ae_{1}+d_{1}c_{1})(fe_{2}+d_{2}c_{2})}\right) \lt \frac{rbe_{1}}{ae_{1}+d_{1}c_{1}} -\left(f+\frac{d_{2}c_{2}}{e_{2}}\right) \frac{rbe_{1}e_{2}}{(ae_{1}+d_{1}c_{1})(fe_{2}+d_{2}c_{2})} = 0. $

    Therefore, for cases (ⅰ) and (ⅱ), $ F(I) = 0 $ has a unique positive root $ I = I^{*} $ if $ R_{0} > 1 $.

    From the above discussions, we have the following result.

    Theorem 2.2. The following statements are true.

    (i) Model (1.3) always has a trivial equilibrium $ \widetilde{E} = (0, 0, 0, 0) $.

    (ii) Model (1.3) always has a disease-free equilibrium $ E_{0} = (S^{0}, 0, u_{1}^{0}, 0) $.

    (iii) Only for $ R_{0} > 1 $, model (1.3) has a unique endemic equilibrium $ E^{*} = (S^{*}, I^{*}, u_{1}^{*}, u_{2}^{*}) $, where

    $ S^{*} = \frac{e_{1}}{ae_{1}+d_{1}c_{1}}\left(r-\frac{bI^{*}}{1+kI^{*}}\right), \; \; u_{1} = \frac{d_{1}}{e_{1}}S^{*}, \; \; u_{2} = \frac{d_{2}}{e_{2}}I^{*}, $

    and $ I^{*} $ is the unique positive root of the equation $ F(I) = 0 $.

    Remak 2.1. In fact, the classifications of the equilibria of models (1.2) and (1.3) are exactly the same for any $ \tau_{i}\geq0 $. Clearly, comparing with the reference [28], our Theorem 2.2 gives more complete classification of the equilibria of model (1.2) and clearer expression of the basic reproduction number $ R_0 $.

    Let $ \overline{E} = (\overline{S}, \overline{I}, \overline{u}_{1}, \overline{u}_{2}) $ be any equilibrium of model (1.3). In order to investigate local stability of the equilibrium $ \overline{E} $, we easily have that the characteristic equation of the corresponding linearized system of model (1.3) at $ \overline{E} $ is given by

    $ |λ(r2a¯Sb¯I1+k¯Ic1¯u1)b¯S(1+k¯I)2c1¯Seλτ10b¯I1+k¯Iλ(b¯S(1+k¯I)2μ2f¯Ic2¯u2)0c2¯Ieλτ2d1eλτ30λ+e100d2eλτ40λ+e2|=0.
    $
    (3.1)

    At the trivial equilibrium $ \widetilde{E} = (0, 0, 0, 0) $, the characteristic Eq (3.1) becomes

    $ (\lambda-r)(\lambda+\mu)(\lambda+e_{1})(\lambda+e_{2}) = 0, $

    which has a positive real root $ \lambda = r $. Hence, $ \widetilde{E} $ is unstable for any $ \tau_{i}\geq0 $ ($ i = 1, 2, 3, 4 $).

    At the disease-free equilibrium $ E_{0} = (S^{0}, 0, u^{0}_{1}, 0) $, the characteristic Eq (3.1) becomes

    $ (λ+μbS0)(λ+e2)[λ2+(aS0+e1)λ+ae1S0+d1c1S0eλ(τ1+τ3)]=0.
    $
    (3.2)

    It is clear that Eq (3.2) has two real roots $ \lambda_{1} = -e_{2} < 0 $ and $ \lambda_{2} = -\mu+bS^{0} = \mu(R_{0}-1) $. Obviously, when $ R_0 > 1 $, Eq (3.2) has a positive real root $ \lambda_{2} > 0 $, and hence, $ E_{0} $ is unstable for any $ \tau_{i}\geq0 $ ($ i = 1, 2, 3, 4 $). When $ R_0 < 1 $, then $ \lambda_{2} < 0 $. When $ R_0 = 1 $, then $ \lambda_{2} = 0 $ is a simple root of Eq (3.2).

    Let

    $ F1(λ,τ1,τ3)λ2+(aS0+e1)λ+ae1S0+d1c1S0eλ(τ1+τ3)=0.
    $
    (3.3)

    The distribution of the roots of Eq (3.3) in the complex plane has been discussed in detail in [30,31,34]. Therefore, we have the following conclusions.

    Lemma 3.1. The following statements are true.

    (i) If $ d_{1}c_{1}\leq ae_{1} $, then all the roots of Eq (3.3) have negative real parts.

    (ii) If $ d_{1}c_{1} > ae_{1} $, then all the roots of Eq (3.3) have negative real parts for $ \tau_{1}+\tau_{3} < \tau_{13}^{0} $, and Eq (3.3) has at least one root which has positive real part for $ \tau_{1}+\tau_{3} > \tau_{13}^{0} $, where

    $ \tau_{13}^{0} = \frac{1}{\omega}\arccos\left[ \frac{\omega^{2}-ae_{1}S^{0}}{d_{1}c_{1}S^{0}}\right], $
    $ \omega = \left(\frac{-\left[(aS^{0})^{2}+e_{1}^{2}\right] +\sqrt{\left[(aS^{0})^{2}+e_{1}^{2}\right]^{2}-4(S^{0})^{2}(ae_{1}+d_{1}c_{1}) (ae_{1}-d_{1}c_{1})}}{2}\right)^{\frac{1}{2}}. $

    According to the discussions above and Lemma 3.1, it follows from stability theory and Hopf bifurcation theorem for DDEs (see, for example, [29,30,31]) that the following results hold.

    Theorem 3.1. The trivial equilibrium $ \widetilde{E} $ of model (1.3) is unstable for any $ \tau_{i}\geq0 \; (i=1,2,3,4

    ) $.

    Theorem 3.2. For any $ \tau_{2}\geq0 $ and $ \tau_{4}\geq0 $, the following statements are true.

    (i) If $ R_{0} > 1 $, then the disease-free equilibrium $ E_{0} $ is unstable for any $ \tau_{1}\geq0 $ and $ \tau_{3}\geq0 $.

    (ii) Assume that $ d_{1}c_{1}\leq ae_{1} $. If $ R_{0} < 1 $, then the disease-free equilibrium $ E_{0} $ is locally asymptotically stable for any $ \tau_{1}\geq0 $ and $ \tau_{3}\geq0 $; If $ R_{0} = 1 $, then the disease-free equilibrium $ E_{0} $ is linearly stable for any $ \tau_{1}\geq0 $ and $ \tau_{3}\geq0 $.

    (iii) Assume that $ d_{1}c_{1} > ae_{1} $. If $ R_{0} < 1 $, then the disease-free equilibrium $ E_{0} $ is locally asymptotically stable for $ \tau_{1}+\tau_{3} < \tau_{13}^{0} $, and is unstable for $ \tau_{1}+\tau_{3} > \tau_{13}^{0} $. Moreover, model (1.3) undergoes a Hopf bifurcation at the disease-free equilibrium $ E_{0} $ when $ \tau_{1}+\tau_{3} = \tau_{13}^{0} $.

    Remak 3.1. Theorem 3.2 indicates that time delays $ \tau_2 $ and $ \tau_4 $ do not affect local asymptotic stability of the disease-free equilibrium $ E_{0} $, and under the condition of $ d_{1}c_{1}\leq ae_{1} $, time delays $ \tau_1 $ and $ \tau_3 $ also do not affect local asymptotic stability of $ E_{0} $. But under the condition $ d_{1}c_{1} > ae_{1} $, for larger time delay $ \tau_1 $ or $ \tau_3 $, stability of $ E_0 $ will be lost.

    In the following discussions, we establish some sufficient conditions for global asymptotic stability of the disease-free equilibrium $ E_{0} $.

    Theorem 3.3. Assume that $ d_{1}c_{1}\leq ae_{1} $. For any $ \tau_{i}\geq0 $ $ (i = 1, 2, 3, 4) $, the following statements are true.

    (i) If $ R_{0} < 1 $, then the disease-free equilibrium $ E_{0} $ is globally asymptotically stable in $ \Omega_{1}: = \{\phi\in \Omega:\phi_{1}(0) > 0\} $.

    (ii) If $ R_{0} = 1 $, then the disease-free equilibrium $ E_{0} $ is globally attractive in $ \Omega_{1} $.

    Proof. First, it is easy to show that the set $ \Omega_{1} $ is positively invariant for model (1.3). If $ R_0 < 1 $, by Theorem 3.2, we only need to show that the disease-free equilibrium $ E_{0} $ is globally attractive.

    Define a Lyapunov functional $ L_{1} $ on $ \Omega_{1} $ as follows,

    $ L_{1} = V_{1} +\frac{a}{2}\int_{-\tau_{3}}^{0}(\phi_1(\xi)-S^{0})^{2}d\xi +\frac{c_1e_{1}}{2d_1}\int_{-\tau_{1}}^{0} (\phi_{3}(\xi)-u_{1}^{0})^{2}d\xi, $

    where

    $ V_{1} = \phi_1(0)-S^{0}-S^{0}\ln\frac{\phi_1(0)}{S^{0}}+\phi_2(0) +\frac{c_1}{2d_1}(\phi_{3}(0)-u_{1}^{0})^{2}. $

    It is clear that $ L_{1} $ is continuous on $ \Omega_{1} $ and satisfies the condition (ii) of Lemma 3.1 in [35] on $ \partial \Omega_{1} = \overline{\Omega}_{1}\setminus \Omega_{1} $.

    Calculating the derivative of $ L_{1} $ along any solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3), it follows that, for $ t\geq0 $,

    $ dL1dt=dV1dt+a2(S(t)S0)2a2(S(tτ3)S0)2+c1e12d1(u1(t)u01)2c1e12d1(u1(tτ1)u01)2,
    $
    (3.4)

    where

    $ dV1dt=(S(t)S0)[a(S0S(t))bI(t)1+kI(t)+c1(u01u1(tτ1))]+I(t)[bS(t)1+kI(t)μfI(t)c2u2(tτ2)]+c1d1(u1(t)u01)[e1(u1(t)u01)+d1(S(tτ3)S0)]=a(S(t)S0)2[μbS01+kI(t)]I(t)fI2(t)c2I(t)u2(tτ2)c1e1d1(u1(t)u01)2+c1(S(t)S0)(u01u1(tτ1))+c1(u1(t)u01)(S(tτ3)S0),
    $
    (3.5)

    here $ r = aS^{0}+c_{1}u_{1}^{0} $ and $ e_{1}u_{1}^{0} = d_{1}S^{0} $ are used. Using the following inequality of arithmetic and geometric means,

    $ Λ1c1(S(t)S0)(u01u1(tτ1))+c1(u1(t)u01)(S(tτ3)S0)d1c1ae1[a2(S(t)S0)2+c1e12d1(u1(tτ1)u01)2+a2(S(tτ3)S0)2+c1e12d1(u1(t)u01)2],
    $

    we further have that

    $ dL1dta2(1d1c1ae1)[(S(t)S0)2+(S(tτ3)S0)2][μbS01+kI(t)]I(t)fI2(t)c2I(t)u2(tτ2)c1e12d1(1d1c1ae1)[(u1(t)u01)2+(u1(tτ1)u01)2].
    $
    (3.6)

    Note that, if $ R_{0}\leq1 $, we have that

    $ [μbS01+kI(t)]I(t)=[μbS0+μkI(t)1+kI(t)]I(t)=[μ(1R0)1+kI(t)+μkI(t)1+kI(t)]I(t)0.
    $
    (3.7)

    Assume that $ d_{1}c_{1} < ae_{1} $ and $ R_{0}\leq1 $. By inequalities (3.6) and (3.7), we can obtain that $ \frac{dL_{1}}{dt}\leq0 $ for $ t\geq0 $. Let $ M $ be the largest invariant set in the following set $ G $:

    $ G: = \left\{\phi\in \overline{\Omega}_{1}:\; L_{1} \lt \infty \; \text{and}\; \frac{dL_{1}}{dt} = 0\right\}. $

    Then, it follows from inequalities (3.6) and (3.7) that

    $ G\subset \left\{\phi\in \Omega_{1}: \; \phi_{1}(0) = S^{0},\; \phi_{2}(0) = 0,\; \phi_{3}(0) = u_{1}^{0}\right\}. $

    We can easily have from model (1.3) and the invariance of $ M $ that $ M = \{E_{0}\} $. Therefore, it follows from Lemma 3.1 in [35] that the disease-free equilibrium $ E_{0} $ is globally attractive.

    Assume that $ d_{1}c_{1} = ae_{1} $ and $ R_{0}\leq1 $. In this case, it is not easy to conclude that the largest invariant set $ M $ is the singleton $ \{E_{0}\} $. Hence, it is necessary to analyze Eq (3.4).

    Note that

    $ a2(S(t)S0)2c1e12d1(u1(tτ1)u01)2+c1(S(t)S0)(u01u1(tτ1))=d1c12e1(S(t)S0+e1d1(u1(tτ1)u01))2,a2(S(tτ3)S0)2c1e12d1(u1(t)u01)2+c1(u1(t)u01)(S(tτ3)S0)=d1c12e1(S(tτ3)S0e1d1(u1(t)u01))2=d1c12e1(S(tτ3)e1d1u1(t))2.
    $

    Hence, Eq (3.4) can be rewritten as

    $ dL1dt=d1c12e1[(S(t)S0+e1d1(u1(tτ1)u01))2+(S(tτ3)e1d1u1(t))2][μ(1R0)1+kI(t)+μkI(t)1+kI(t)]I(t)fI2(t)c2I(t)u2(tτ2).
    $
    (3.8)

    By inequality (3.7) and Eq (3.8), we can obtain that $ \frac{dL_{1}}{dt}\leq0 $ for $ t\geq0 $. Let $ M_{1} $ be the largest invariant set in the following set $ G_1 $:

    $ G_{1}: = \left\{\phi\in \overline{\Omega}_{1}:\; L_{1} \lt \infty \; \text{and}\; \frac{dL_{1}}{dt} = 0\right\}. $

    Then, it follows from Eq (3.8) that

    $ G_{1}\subset \left\{\phi\in \Omega_{1}: \; \phi_{1}(0)+\frac{e_{1}}{d_{1}}\phi_{3}(-\tau_{1}) = S^{0}+\frac{e_{1}}{d_{1}}u_{1}^{0}, \; \phi_{1}(-\tau_{3})-\frac{e_{1}}{d_{1}}\phi_{3}(0) = 0, \; \phi_{2}(0) = 0\right\}. $

    For any $ \varphi = (\varphi_{1}, \varphi_{2}, \varphi_{3}, \varphi_{4})\in M_{1} $, let $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ be the solution of model (1.3) with the initial function $ \varphi $. From the invariance of $ M_{1} $, we have that $ (S_{t}, I_{t}, u_{1t}, u_{2t})\in M_{1}\subset G_{1} $ for any $ t\in R $. Obviously, $ I(t)\equiv0 $ for any $ t\in R $, and then from the fourth equation of model (1.3) and the invariance of $ M_{1} $, it is not difficult to obtain $ u_{2}(t)\equiv0 $ for any $ t\in R $. In addition, according to the first and third equations of model (1.3), we can obtain, for any $ t\in R $,

    $ ˙S(t)=S(t)(raS(t)c1u1(tτ1))=S(t)(rd1c1e1S(t)c1u1(tτ1))=0,˙u1(t)=e1u1(t)+d1S(tτ3)=0.
    $

    Thus, there exist constants $ \delta_1 $ and $ \delta_2 $ such that $ S(t)\equiv \delta_{1} $ and $ u_{1}(t)\equiv \delta_{2} $ for any $ t\in R $. It is not difficult to find that $ \delta_1 $ and $ \delta_2 $ satisfy

    $ \delta_{1}+\frac{e_{1}}{d_{1}}\delta_{2} = S^{0}+\frac{e_{1}}{d_{1}}u_{1}^{0}, \; \; \delta_{1}-\frac{e_{1}}{d_{1}}\delta_{2} = 0, $

    which imply that $ \delta_1 = S^{0} $ and $ \delta_2 = u_{1}^{0} $. Hence, $ S(t)\equiv S^{0} $ and $ u_{1}(t)\equiv u_{1}^{0} $ for any $ t\in R $. This shows that $ M_{1} = \{E_{0}\} $. Then, it follows from Lemma 3.1 in [35] that the disease-free equilibrium $ E_{0} $ is globally attractive.

    The proof is completed.

    Remak 3.2. Note the conclusion (ii) of Theorem 3.2, where we see that Theorem 3.3 gives complete conclusion of the global dynamics of the disease-free equilibrium $ E_{0} $ in the case of $ d_{1}c_{1}\leq ae_{1} $.

    Now, we continue to discuss global dynamics of the disease-free equilibrium $ E_{0} $ in the absence of condition $ d_{1}c_{1}\leq ae_{1} $. The following lemmas will be used.

    Lemma 3.2. $ ($Barbalat's lemma [37,36]$) $ Let $ x(t) $ be a real valued differentiable function defined on some half line $ [a, +\infty) $, $ a\in(-\infty, +\infty) $. If

    (i) $ \lim_{t\rightarrow+\infty}x(t) = \alpha $; $ |\alpha| < \infty $.

    (ii) $ \dot{x}(t) $ is uniformly continuous for $ t > a $.

    Then $ \lim_{t\rightarrow+\infty}\dot{x}(t) = 0 $.

    Lemma 3.3. Let $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ be any solution of model (1.3) with the initial condition (2.1), then the following statements are true.

    (i) If $ rb\leq a\mu $ (which implies $ R_{0} < 1 $), then $ \lim_{t\rightarrow+\infty}I(t) = 0 $, $ \lim_{t\rightarrow+\infty}u_{2}(t) = 0 $.

    (ii) If $ rb > a\mu $, then $ \limsup_{t\rightarrow+\infty}I(t)\leq I_{M} $, where

    $ I_{M} = \left\{ (kμ+f)2+4fk(braμ)(kμ+f)2fk>0,k>0,rbaμaf,k=0.
    \right. $

    Proof. By inequality (2.2), for arbitrary $ \varepsilon > 0 $, there exists a $ T > 0 $ such that $ S(t) < \frac{r}{a}+\varepsilon $ for $ t > T $. Then, it follows from the second equation of model (1.3) that, for $ t > T $,

    $ ˙I(t)I(t)(b(ra+ε)1+kI(t)μfI(t))=I(t)1+kI(t)[fkI2(t)+(μk+f)I(t)b(ra+ε)+μ].
    $
    (3.9)

    If $ rb\leq a\mu $, by inequality (3.9), it follows that, for $ t > T $,

    $ ˙I(t)I(t)1+kI(t)[(μk+f)I(t)bε],
    $

    which implies

    $ lim supt+I(t)bεμk+f.
    $
    (3.10)

    Since $ I(t)\geq0 $ and inequality (3.10) holds for arbitrary $ \varepsilon > 0 $, we obtain that $ \lim_{t\rightarrow+\infty}I(t) = 0 $. Further, according to the last equation of model (1.3), $ \lim_{t\rightarrow+\infty}u_{2}(t) = 0 $ can be easily obtained.

    If $ rb > a\mu $, it follows from inequality (3.9) that, for $ t > T $,

    $ \dot{I}(t)\leq \left\{ fkI(t)1+kI(t)(I(t)I1(ε))(I(t)I2(ε)),k>0,f(I(t)I2(ε)),k=0.
    \right. $

    where

    $ I1(ε)=(μk+f)(μk+f)2+4fk(bra+bεμ)2fk<0,k>0I2(ε)={(μk+f)+(μk+f)2+4fk(bra+bεμ)2fk>0,k>0,b(r+aε)aμaf,k=0.
    $

    Similarly, it follows that

    $ lim supt+I(t)I2(ε),andlim supt+I(t)I2(0)=IM.
    $
    (3.11)

    The proof is completed.

    For convenience, let us give the following conditions $ (H_1) $, $ (H_2) $ and $ (H_3) $:

    $ (H1)c1(e1+2d1)2τ1+c1r2τ3<a.(H2)e12τ1+r2a(a+b+2c1)τ3<e1d1.(H3)rc1b2aτ3<f+μk.
    $

    We can obtain the following result.

    Theorem 3.4. Assume that $ R_{0}\leq1 $. For any $ \tau_{2}\geq0 $ and $ \tau_{4}\geq0 $, the following statements are true.

    (i) If $ rb\leq a\mu $ and conditions $ (H_{1}) $–$ (H_{2}) $ hold, then the disease-free equilibrium $ E_{0} $ is globally attractive in $ X_{1}: = \{\phi\in C^{+}:\; \phi_{1}(0) > 0\} $.

    (ii) If $ rb > a\mu $ and conditions $ (H_{1}) $–$ (H_{3}) $ hold, then the disease-free equilibrium $ E_{0} $ is globally attractive in $ X_{1} $.

    Proof. It is easy to show that $ X_{1} $ is positively invariant for model (1.3). Let $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ be the solution of model (1.3) with any initial function $ \phi\in X_{1} $. By inequality (2.2), for any sufficiently small $ \varepsilon_{0} > 0 $, there exists a $ t_{1} > 0 $ such that $ S(t) < \frac{r}{a}+\varepsilon_{0}\equiv S_{1}(\varepsilon_{0}) $ for $ t > t_{1} $.

    We continue to analyze $ \frac{dV_{1}}{dt} $ given by Eq (3.5). From Eq (3.5) and inequality (3.7), we have that, for $ t\geq 0 $,

    $ dV1dt=a(S(t)S0)2μ(1R0)1+kI(t)I(t)μk1+kI(t)I2(t)fI2(t)c2I(t)u2(tτ2)c1e1d1(u1(t)u01)2+Λ1.
    $
    (3.12)

    Note that, for $ t > \tau_1+\tau_3 $, $ \Lambda_{1} $ can be rewritten as

    $ Λ1=c1(S(t)S0)(u1(t)u1(tτ1))+c1(u1(t)u01)(S(tτ3)S(t)):=Λ2+Λ3,
    $

    where

    $ Λ2=c1(S(t)S0)(u1(t)u1(tτ1))=c1(S(t)S0)ttτ1˙u1(ξ)dξ=c1(S(t)S0)ttτ1(e1(u1(ξ)u01)+d1(S(ξτ3)S0))dξ,Λ3=c1(u1(t)u01)(S(tτ3)S(t))=c1(u1(t)u01)ttτ3˙S(ξ)dξ=c1(u1(t)u01)ttτ3S(ξ)(a(S0S(ξ))bI(ξ)1+kI(ξ)+c1(u01u1(ξτ1)))dξ.
    $

    In addition, for $ t > \tau_1+\tau_3 $, we have

    $ Λ2c1e12ttτ1((S(t)S0)2+(u1(ξ)u01)2)dξ+c1d12ttτ1((S(t)S0)2+(S(ξτ3)S0)2)dξ=c1(e1+d1)2τ1(S(t)S0)2+c1e12ttτ1(u1(ξ)u01)2dξ+c1d12ttτ1(S(ξτ3)S0)2dξ.
    $
    (3.13)

    Similarly, for $ t > t_{1}+\tau_1+\tau_3 $, we have

    $ Λ3c1S(ε0)|u1(t)u01|ttτ3(a|S0S(ξ)|+bI(ξ)1+kI(ξ)+c1|u01u1(ξτ1)|)dξc1aS(ε0)2ttτ3((u1(t)u01)2+(S0S(ξ))2)dξ+c1bS(ε0)2ttτ3((u1(t)u01)2+I2(ξ)(1+kI(ξ))2)dξ+c21S(ε0)2ttτ3((u1(t)u01)2+(u01u1(ξτ1))2)dξ=c1S(ε0)2(a+b+c1)τ3(u1(t)u01)2+c1aS(ε0)2ttτ3(S0S(ξ))2dξ+c1bS(ε0)2ttτ3I2(ξ)(1+kI(ξ))2dξ+c21S(ε0)2ttτ3(u01u1(ξτ1))2dξ,
    $
    (3.14)

    here $ S(t) < S_{1}(\varepsilon_{0}) $ for $ t > t_{1} $ is used. For $ t > t_{1}+\tau_1+\tau_3 $, let us define the following function $ L_2 $,

    $ L_{2} = \sum\limits_{i = 1}^{6}V_{i}, $

    where $ V_1 $ is defined as in the proof of Theorem 3.3,

    $ V2=c1aS(ε0)2ttτ3tθ(S(ξ)S0)2dξdθ,V3=c1d12ttτ1tθ(S(ξτ3)S0)2dξdθ+c1d1τ12ttτ3(S(ξ)S0)2dξ,V4=c1bS(ε0)2ttτ3tθI2(ξ)(1+kI(ξ))2dξdθ,V5=c1e12ttτ1tθ(u1(ξ)u01)2dξdθ,V6=c21S(ε0)2ttτ3tθ(u1(ξτ1)u01)2dξdθ+c21S(ε0)2τ3ttτ1(u1(ξ)u01)2dξ.
    $

    For $ t > t_{1}+\tau_1+\tau_3 $, we calculate the derivatives of $ V_{i} \; (i = 2, \cdots, 6) $ as follows,

    $ dV2dt=c1aS(ε0)2[τ3(S(t)S0)2ttτ3(S(ξ)S0)2dξ],
    $
    (3.15)
    $ dV3dt=c1d12[τ1(S(t)S0)2ttτ1(S(ξτ3)S0)2dξ],
    $
    (3.16)
    $ dV4dt=c1bS(ε0)2[τ3I2(t)(1+kI(t))2ttτ3I2(ξ)(1+kI(ξ))2dξ],
    $
    (3.17)
    $ dV5dt=c1e12[τ1(u1(t)u01)2ttτ1(u1(ξ)u01)2dξ],
    $
    (3.18)
    $ dV6dt=c21S(ε0)2[τ3(u1(t)u01)2ttτ3(u1(ξτ1)u01)2dξ].
    $
    (3.19)

    Combining (3.12)–(3.19), for $ t > t_{1}+\tau_1+\tau_3 $, we finally have

    $ dL2dt=6i=1dVidt[ac1(e1+2d1)2τ1c1aS(ε0)2τ3](S(t)S0)2μ(1R0)1+kI(t)I(t)+c1bS(ε0)2τ3I2(t)(1+kI(t))2μk1+kI(t)I2(t)fI2(t)c2I(t)u2(tτ2)[c1e1d1c1e12τ1c1S(ε0)2(a+b+2c1)τ3](u1(t)u01)2.
    $
    (3.20)

    Let us show global attractivity of the disease-free equilibrium $ E_{0} $ in the following two cases.

    Case (ⅰ) $ rb\leq a\mu $ (which implies $ R_{0} < 1 $) and $ (H_1) $–$ (H_2) $ hold.

    From conditions $ (H_{1}) $–$ (H_{2}) $, it has that, for sufficiently small $ \varepsilon_{0} > 0 $, the inequalities

    $ ac1(e1+2d1)2τ1c1aS(ε0)2τ3>0,e1d1e12τ1S(ε0)2(a+b+2c1)τ3>0
    $
    (3.21)

    hold. Furthermore, from Lemma 3.3, we have that $ \lim_{t\rightarrow\infty}I(t) = 0 $ and $ \lim_{t\rightarrow\infty}u_{2}(t) = 0 $. Thus, there exists a $ t_{2} > 0 $ such that, for $ t > t_{2} $,

    $ c_{1}bS(\varepsilon_{0})\tau_{3}I(t) \lt \frac{\mu d_{1}c_{1}}{ae_{1}+d_{1}c_{1}}. $

    In addition, note that

    $ 1-R_{0} = 1-\frac{bre_{1}}{\mu(ae_{1}+d_{1}c_{1})} = \frac{\mu(ae_{1}+d_{1}c_{1})-bre_{1}}{\mu(ae_{1}+d_{1}c_{1})} \geq\frac{d_{1}c_{1}}{ae_{1}+d_{1}c_{1}}, $

    from which we have that, for $ t > t_{2} $,

    $ Λ4:=μ(1R0)1+kI(t)I(t)+c1bS(ε0)2τ3I2(t)(1+kI(t))2I(t)1+kI(t)[μd1c1ae1+d1c1c1bS(ε0)2τ3I(t)]μd1c12(ae1+d1c1)I(t)1+kI(t)0.
    $
    (3.22)

    Hence, it follows from inequalities (3.20)–(3.22) that $ \frac{dL_{2}}{dt}\leq0 $ for $ t > \widetilde{T}: = \max\{t_{1}+\tau_{1}+\tau_3, t_{2}\} $. This indicates that, for $ t > \widetilde{T} $, the function $ L_{2}(t) $ is monotonically decreasing and bounded. Thus, the limitation $ \lim_{t\rightarrow+\infty}L_{2}(t) $ exists.

    In addition, by Theorem 2.1, it is not difficult to show that, for $ t > \widetilde{T} $, the second derivative $ L''_2(t) $ is also bounded. This implies that $ L'_{2}(t) $ is uniformly continuous for $ t > \widetilde{T} $. Therefore, it follows from Lemma 3.1 that $ \lim_{t\rightarrow\infty}L'_{2}(t) = 0 $. Again from inequalities (3.20)–(3.22), it follows that

    $ \lim\limits_{t\rightarrow\infty}S(t) = S^{0},\; \; \; \; \lim\limits_{t\rightarrow\infty}u_{1}(t) = u_{1}^{0}. $

    This shows that the disease-free equilibrium $ E_{0} $ is globally attractive.

    Case (ⅱ) $ rb > a\mu $ and $ (H_1) $–$ (H_3) $ hold.

    Similarly, from condition $ (H_{3}) $, we have that the inequality

    $ c1bS(ε0)2τ3<f+μk
    $
    (3.23)

    holds for sufficiently small $ \varepsilon_{0} > 0 $. From Lemma 3.3, there exists a $ t_{3} > 0 $ such that $ I(t) < I_M+\varepsilon_{0} $ for $ t > t_{3} $. Then, we have that, for $ t > t_3 $,

    $ Λ5:=c1bS(ε0)2τ3I2(t)(1+kI(t))2μk1+kI(t)I2(t)fI2(t)=I2(t)(1+kI(t))2[f(1+kI(t))2+μk(1+kI(t))c1bS(ε0)2τ3]I2(t)(1+k(IM+ε0))2[f+μkc1bS(ε0)2τ3]0.
    $
    (3.24)

    Hence, by inequalities (3.20), (3.21) and (3.24), we can also obtain $ \frac{dL_{2}}{dt}\leq0 $ for $ t > \hat{T}: = \max\{t_{1}+\tau_{1}+\tau_3, t_{3}\} $. By the same arguments as in Case (i), we can obtain

    $ \lim\limits_{t\rightarrow\infty}S(t) = S^{0}, \; \; \lim\limits_{t\rightarrow\infty}I(t) = 0, \; \; \lim\limits_{t\rightarrow\infty}u_{1}(t) = u_{1}^{0}. $

    Furthermore, from the last equation of model (1.3), we can easily have $ \lim_{t\rightarrow\infty}u_{2}(t) = 0 $. This shows that the disease-free equilibrium $ E_{0} $ is globally attractive.

    The proof is completed.

    From Theorems 3.2 and 3.4, we have the following two corollaries.

    Corollary 3.1. Assume that $ R_{0} < 1 $, $ d_{1}c_{1} > ae_{1} $ and $ \tau_{1}+\tau_{3} < \tau_{13}^{0} $. For any $ \tau_{2}\geq0 $ and $ \tau_{4}\geq0 $, the following statements are true.

    (i) If $ rb\leq a\mu $ and conditions $ (H_{1})-(H_{2}) $ hold, then the disease-free equilibrium $ E_{0} $ is globally asymptotically stable in $ X_{1} $.

    (ii) If $ rb > a\mu $ and conditions $ (H_{1})-(H_{3}) $ hold, then the disease-free equilibrium $ E_{0} $ is globally asymptotically stable in $ X_{1} $.

    Corollary 3.2. Assume that $ R_{0} < 1 $ and $ \tau_{1} = \tau_{3} = 0 $. For any $ \tau_{2}\geq0 $ and $ \tau_{4}\geq0 $, then the disease-free equilibrium $ E_{0} $ is globally asymptotically stable in $ X_{1} $.

    Remak 3.3. If $ \tau_{i} = 0 \; (i=1,2,3,4

    ) $, then model (1.3) reduces into model (1.2). Clearly, Theorems 3.2, 3.3 and 3.4 extend and improve Theorem 1 in [28]. Further, if $ k = 0 $, model (1.2) becomes the model discussed in [17]. Hence, Theorems 3.2, 3.3 and 3.4 also include Theorem 2.1 in [17] as a special case.

    Theoretical analysis of the distribution of the characteristic roots of characteristic equation (3.1) at the endemic equilibrium $ E^{*} = (S^{*}, I^{*}, u_{1}^{*}, u_{2}^{*}) $ usually involves some complicated computations, since there are four time delays $ \tau_{i}\geq0 \; (i=1,2,3,4

    ) $ in model (1.3). However, the numerical simulations in Section 5 show that, each of time delays $ \tau_{i} $ ($ i = 1, 2, 3, 4 $) can destroy stability of the endemic equilibrium $ E^{*} $ of model (1.3) by properly choosing parameters. Hence, it is natural to consider the following two problems.

    (ⅰ) Under what conditions, the time delays $ \tau_{i}\geq0 \; (i=1,2,3,4

    ) $ are harmless for global asymptotic stability of the endemic equilibrium $ E^{*} $ of model (1.3).

    (ⅱ) Under what conditions, the time delays $ \tau_{i}\geq0 \; (i=1,2,3,4

    ) $ may be harmful for global asymptotic stability of the endemic equilibrium $ E^{*} $ of model (1.3).

    In this section, we study the above problems (i)–(ii) by constructing suitable Lyapunov functionals.

    For convenience, let us denote the following conditions,

    $ (H4)τ1=τ2=τ3=τ4=0.(H5)d1c1ae1,τ10,τ30,τ2=τ4=0.(H6)d2c2fe2,τ20,τ40,τ1=τ3=0.(H7)d1c1ae1,d2c2fe2,τ10,τ20,τ30,τ40.
    $

    For problem (ⅰ) above, we have the following result.

    Theorem 4.1. Assume that $ R_{0} > 1 $. If one of conditions $ (H_4) $–$ (H_7) $ holds, then the endemic equilibrium $ E^{*} $ is globally asymptotically stable in $ X_{2}: = \{\phi\in C^{+}: \phi_{i}(0) > 0, \; i = 1, 2\} $.

    Proof. It is easy to show that the set $ X_{2} $ is positively invariant for model (1.3). Since $ E^{*} $ is the equilibrium of model (1.3), the following equalities hold,

    $ raSbI1+kIc1u1=0,bS1+kIμfIc2u2=0,e1u1+d1S=0,e2u2+d2I=0.
    $
    (4.1)

    Define a Lyapunov functional $ W_{1} $ on $ X_{2} $ as follows,

    $ W1=(1+kI)(ϕ1(0)SSlnϕ1(0)S)+ϕ2(0)IIlnϕ2(0)I+c1(1+kI)2d1(ϕ3(0)u1)2+c22d2(ϕ4(0)u2)2.
    $

    It is clear that $ W_{1} $ is continuous on $ X_{2} $ and positive definite with respect to $ E^{*} $.

    Calculating the derivative of $ W_{1} $ along any solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3), it follows that, for $ t\geq0 $,

    $ dW1dt=(1+kI)(S(t)S)[a(SS(t))+bI1+kIbI(t)1+kI(t)+c1(u1u1(tτ1))]+(I(t)I)[bS(t)1+kI(t)bS1+kI(t)+bS1+kI(t)bS1+kI+f(II(t))+c2(u2u2(tτ2))]+c1(1+kI)d1(u1(t)u1)[e1(u1(t)u1)+d1(S(tτ3)S)]+c2d2(u2(t)u2)[e2(u2(t)u2)+d2(I(tτ4)I)],
    $

    here Eq (4.1) is used. Note that

    $ (1+kI^{*})(S(t)-S^{*})\left(\frac{bI^{*}}{1+kI^{*}}-\frac{bI(t)}{1+kI(t)}\right) +(I(t)-I^{*})\left(\frac{bS(t)}{1+kI(t)}-\frac{bS^{*}}{1+kI(t)}\right) = 0, $

    from which we have that, for $ t\geq0 $,

    $ dW1dt=a(1+kI)(S(t)S)2[f+bkS(1+kI(t))(1+kI)](I(t)I)2c1e1(1+kI)d1(u1(t)u1)2c2e2d2(u2(t)u2)2+(1+kI)Υ1+Π1,
    $
    (4.2)

    where

    $ Υ1=c1(S(t)S)(u1u1(tτ1))+c1(S(tτ3)S)(u1(t)u1),Π1=c2(I(t)I)(u2u2(tτ2))+c2(I(tτ4)I)(u2(t)u2).
    $

    We consider global asymptotic stability of the endemic equilibrium $ E^{*} $ in the following four cases.

    If condition $ (H_4) $ holds, it has that $ \Upsilon_{1} = \Pi_{1} = 0 $ and $ \frac{dW_{1}}{dt} $ is negative definite with respect to $ E^{*} $. Hence, the endemic equilibrium $ E^{*} $ is globally asymptotically stable (see, for example, [29,30]).

    If condition $ (H_5) $ holds, we have that $ \Pi_{1} = 0 $. Let us consider another functional as follows,

    $ W_{2} = W_{1}+\frac{a(1+kI^{*})}{2}\int_{-\tau_{3}}^{0}(\phi_1(\xi)-S^{*})^2d\xi +\frac{c_1e_{1}(1+kI^{*})}{2d_1} \int_{-\tau_{1}}^{0}\left(\phi_3(\xi)-u_{1}^{*}\right)^{2}d\xi. $

    $ W_{2} $ is continuous on $ X_{2} $, positive definite with respect to $ E^{*} $, and satisfies condition (ii) of Lemma 3.1 in [35] on $ \partial X_{2} = \overline{X_{2}}\setminus X_{2} $.

    Calculating the derivative of $ W_{2} $ along any solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3), it follows that, for $ t\geq0 $,

    $ dW2dt=dW1dt+a(1+kI)2[(S(t)S)2(S(tτ3)S)2]+c1e1(1+kI)2d1[(u1(t)u1)2(u1(tτ1)u1)2].
    $
    (4.3)

    By Eqs (4.2) and (4.3), and the following inequality of arithmetic and geometric means:

    $ Υ1d1c1ae1[a2(S(t)S)2+c1e12d1(u1(tτ1)u1)2+a2(S(tτ3)S)2+c1e12d1(u1(t)u1)2],
    $
    (4.4)

    we have that, for $ t\geq0 $,

    $ dW2dta(1+kI)2(1d1c1ae1)[(S(t)S)2+(S(tτ3)S)2]c1e1(1+kI)2d1(1d1c1ae1)[(u1(t)u1)2+(u1(tτ1)u1)2][f+bkS(1+kI(t))(1+kI)](I(t)I)2c2e2d2(u2(t)u2)2.
    $

    It follows from condition $ (H_{5}) $ that $ \frac{dW_{2}}{dt}\leq0 $ for $ t\geq0 $. Hence, the endemic equilibrium $ E^{*} $ is stable. Furthermore, $ \frac{dW_{3}}{dt} = 0 $ implies $ I(t) = I^{*} $ and $ u_{2}(t) = u_{2}^{*} $.

    Let $ M_{2} $ be the largest invariant set in the set

    $ \Gamma_{2}: = \left\{\phi\in \overline{X}_{2}:\; W_{2} \lt \infty \; \text{and}\; \frac{dW_{2}}{dt} = 0\right\}. $

    Then, it follows that

    $ \Gamma_{2} \subset \{\phi\in X_{2}: \phi_{2}(0) = I^{*},\; \phi_{4}(0) = u_{2}^{*}\}. $

    From model (1.3) and the invariance of $ M_{2} $, we can easily get that $ M_{2} = \{E^{*}\} $. Thus, it follows from Lemma 3.1 in [35] that the endemic equilibrium $ E^{*} $ is globally asymptotically stable.

    If condition $ (H_6) $ holds, it follows that $ \Upsilon_{1} = 0 $. Let us consider a functional as follows,

    $ W_{3} = W_{1} +\frac{f}{2}\int_{-\tau_{4}}^{0}\left(\phi_2(\xi)-I^{*}\right)^{2}d\xi +\frac{c_2e_{2}}{2d_2}\int_{-\tau_{2}}^{0}\left(\phi_4(\xi)-u_{2}^{*}\right)^{2}d\xi. $

    $ W_{3} $ is continuous on $ X_{2} $, positive definite with respect to $ E^{*} $, and satisfies condition (ii) of Lemma 3.1 in [35] on $ \partial X_{2} = \overline{X_{2}}\setminus X_{2} $.

    Calculating the derivative of $ W_{3} $ along any solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3), it follows that, for $ t\geq0 $,

    $ dW3dt=dW1dt+f2[(I(t)I)2(I(tτ4)I)2]+c2e22d2[(u2(t)u2)2(u2(tτ2)u2)2].
    $
    (4.5)

    By Eqs (4.2) and (4.5), and the following inequality of arithmetic and geometric means:

    $ Π1d2c2fe2[f2(I(t)I)2+c2e22d2(u2(tτ2)u2)2+f2(I(tτ4)I)2+c2e22d2(u2(t)u2)2],
    $
    (4.6)

    we have that, for $ t\geq0 $,

    $ dW3dta(1+kI)(S(t)S)2c1e1(1+kI)d1(u1(t)u1)2f2(1d2c2fe2)[(I(t)I)2+(I(tτ4)I)2]bkS(1+kI(t))(1+kI)(I(t)I)2c2e22d2(1d2c2fe2)[(u2(t)u2)2+(u2(tτ2)u2)2].
    $

    It follows from condition $ (H_{6}) $ that $ \frac{dW_{3}}{dt}\leq0 $ for $ t\geq0 $. Furthermore, $ \frac{dW_{3}}{dt} = 0 $ implies $ S(t) = S^{*} $ and $ u_{1}(t) = u_{1}^{*} $. By the same arguments as in the situation of condition $ (H_5) $, we can also show that the endemic equilibrium $ E^{*} $ is globally asymptotically stable.

    If condition $ (H_7) $ holds, let us consider the following Lyapunov functional,

    $ W4=W1+a(1+kI)20τ3(ϕ1(ξ)S)2dξ+f20τ4(ϕ2(ξ)I)2dξ+c1e1(1+kI)2d10τ1(ϕ3(ξ)u1)2dξ+c2e22d20τ2(ϕ4(ξ)u2)2dξ.
    $

    $ W_{4} $ is continuous on $ X_{2} $, positive definite with respect to $ E^{*} $, and satisfies condition (ii) of Lemma 3.1 in [35] on $ \partial X_{2} = \overline{X_{2}}\setminus X_{2} $.

    Calculating the derivative of $ W_{4} $ along any solution $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ of model (1.3), it follows that, for $ t\geq0 $,

    $ dW4dt=dW1dt+a(1+kI)2[(S(t)S)2(S(tτ3)S)2]+f2[(I(t)I)2(I(tτ4)I)2]+c1e1(1+kI)2d1[(u1(t)u1)2(u1(tτ1)u1)2]+c2e22d2[(u2(t)u2)2(u2(tτ2)u2)2].
    $

    Further, by Eq (4.2) and inequalities (4.4) and (4.6), we have that, for $ t\geq0 $,

    $ dW4dta(1+kI)2(1d1c1ae1)[(S(t)S)2+(S(tτ3)S)2]f2(1d2c2fe2)[(I(t)I)2+(I(tτ4)I)2]bkS(1+kI(t))(1+kI)(I(t)I)2c1e1(1+kI)2d1(1d1c1ae1)[(u1(t)u1)2+(u1(tτ1)u1)2]c2e22d2(1d2c2fe2)[(u2(t)u2)2+(u2(tτ2)u2)2].
    $

    It follows from condition $ (H_{7}) $ that $ \frac{dW_{4}}{dt}\leq0 $ for $ t\geq0 $. Furthermore, if $ d_{1}c_{1} < ae_{1} $, then $ \frac{dW_{4}}{dt} = 0 $ implies that $ S(t) = S^{*} $ and $ u_{1}(t) = u_{1}^{*} $. By the same arguments as in the situation of condition $ (H_5) $, we can show that the endemic equilibrium $ E^{*} $ is globally asymptotically stable. If $ d_{1}c_{1} = ae_{1} $, also by the same arguments as in the situation of $ d_{1}c_{1} = ae_{1} $ in Theorem 3.3, we can show that $ E^{*} $ is globally asymptotically stable.

    The proof is completed.

    Remak 4.1. In the situation of condition $ (H_7) $, Theorem 4.1 indicates that the time delays $ \tau_{i}\geq0 \; (i = 1, 2, 3, 4) $ are harmless for global asymptotic stability of the endemic equilibrium $ E^{*} $.

    Now, let us consider problem (ⅱ). Note that $ R_{0} = \frac{bre_{1}}{\mu(ae_{1}+d_{1}c_{1})} > 1 $ implies $ rb > a\mu $. Let us denote the following conditions,

    $ (H8)c1(e1+2d1)2τ1+c1r2τ3+c2bIM2(1+kI)(1+kIM)τ4<a.(H9)c2(e2+2d2)2τ2+c1br2aτ3+c22[fIM+bkSIM(1+kI)(1+kIM)]τ4<f+bkS(1+kI)(1+kIM).(H10)e12τ1+r2a(a+b1+kI+2c1)τ3<e1d1.(H11)e22τ2+12[bIM1+kIM+fIM+bkSIM(1+kI)(1+kIM)+2c2IM]τ4<e2d2.
    $

    In conditions $ (H_8) $–$ (H_{11}) $ above, the definition of $ I_{M} $ is given in Lemma 3.3 (ⅱ). We have the following result.

    Theorem 4.2. Assume that $ R_{0} > 1 $. If conditions $ (H_{8}) $–$ (H_{11}) $ hold, then the endemic equilibrium $ E^{*} $ is globally attractive in $ X_{2} $.

    Proof. Let $ (S(t), I(t), u_{1}(t), u_{2}(t)) $ be the solution of model (1.3) with any initial function $ \phi\in X_{2} $. From Theorem 2.1 and Lemma 3.3, for sufficient small $ \varepsilon_1 > 0 $, there exists a $ T_{1} > 0 $ such that, for $ t > T_{1} $,

    $ S(t) \lt \frac{r}{a}+\varepsilon_{1}: = S_{1}(\varepsilon_{1}), \; \; I(t) \lt I_{M}+\varepsilon_{1}: = I_{M}(\varepsilon_{1}). $

    Hence, from Eq (4.2), it has that, for $ t > T_{1} $,

    $ dW1dta(1+kI)(S(t)S)2[f+bkS(1+kIM(ε1))(1+kI)](I(t)I)2c1e1(1+kI)d1(u1(t)u1)2c2e2d2(u2(t)u2)2+(1+kI)|Υ1|+|Π1|.
    $
    (4.7)

    For simplicity, denote

    $ A(\varepsilon_{1}): = \left(f+ \frac{bkS^{*}} {(1+kI^{*})(1+kI_{M}(\varepsilon_{1}))}\right)I_{M}(\varepsilon_{1}). $

    Note that $ \Upsilon_{1} $ and $ \Pi_{1} $ can be rewritten as

    $ Υ1=c1(S(t)S)(u1(t)u1(tτ1))+c1(S(tτ3)S(t))(u1(t)u1):=Υ2+Υ3,Π1=c2(I(t)I)(u2(t)u2(tτ2))+c2(I(tτ4)I(t))(u2(t)u2):=Π2+Π3.
    $

    Let us give appropriate estimations on $ |\Upsilon_{2}| $, $ |\Upsilon_{3}| $, $ |\Pi_{2}| $ and $ |\Pi_{3}| $.

    It follows that, $ t > \tau_1+\tau_2+\tau_3+\tau_4+T_1 $,

    $ |Υ2|=|c1(S(t)S)ttτ1˙u1(ξ)dξ|=|c1(S(t)S)ttτ1(e1(u1(ξ)u1)+d1(S(ξτ3)S))dξ|c1e12ttτ1((S(t)S)2+(u1(ξ)u1)2)dξ+c1d12ttτ1((S(t)S)2+(S(ξτ3)S)2)dξ=c1(e1+d1)2τ1(S(t)S)2+c1e12ttτ1(u1(ξ)u1)2dξ+c1d12ttτ1(S(ξτ3)S)2dξ,
    $
    (4.8)
    $ |Υ3|=|c1(u1(t)u1)ttτ3˙S(ξ)dξ|=|c1(u1(t)u1)ttτ3S(ξ)[a(SS(ξ))+bI1+kIbI(ξ)1+kI(ξ)+c1(u1u1(ξτ1))]dξ|=|c1(u1(t)u1)ttτ3S(ξ)[a(SS(ξ))+b(II(ξ))(1+kI)(1+kI(ξ))+c1(u1u1(ξτ1))]dξ|
    $
    $ c1S1(ε1)|u1(t)u1|ttτ3(a|SS(ξ)|+b1+kI|II(ξ)|+c1|u1u1(ξτ1)|)dξc1aS1(ε1)2ttτ3((u1(t)u1)2+(SS(ξ))2)dξ+c1bS1(ε1)2(1+kI)ttτ3((u1(t)u1)2+(II(ξ))2)dξ+c21S1(ε1)2ttτ3((u1(t)u1)2+(u1u1(ξτ1))2)dξ=c1S1(ε1)2[a+b1+kI+c1]τ3(u1(t)u1)2+c1aS1(ε1)2ttτ3(SS(ξ))2dξ+c1bS1(ε1)2(1+kI)ttτ3(II(ξ))2dξ+c21S1(ε1)2ttτ3(u1u1(ξτ1))2dξ,
    $
    (4.9)
    $ |Π2|=|c2(I(t)I)ttτ2˙u2(ξ)dξ|=|c2(I(t)I)ttτ2(e2(u2(ξ)u2)+d2(I(ξτ4)I))dξ|c2e22ttτ2((I(t)I)2+(u2(ξ)u2)2)dξ+c2d22ttτ2((I(t)I)2+(I(ξτ4)I)2)dξc2(e2+d2)2τ2(I(t)I)2+c2e22ttτ2(u2(ξ)u2)2dξ+c2d22ttτ2(I(ξτ4)I)2dξ.
    $
    (4.10)
    $ |Π3|=|c2(u2(t)u2)ttτ4˙I(ξ)dξ|=|c2(u2(t)u2)ttτ4I(ξ)[bS(ξ)1+kI(ξ)bS1+kI+f(II(ξ))+c2(u2u2(ξτ2))]dξ|=|c2(u2(t)u2)ttτ4[bI(ξ)1+kI(ξ)(S(ξ)S)+bkSI(ξ)(1+kI)(1+kI(ξ))(II(ξ))+fI(ξ)(II(ξ))+c2I(ξ)(u2u2(ξτ2))]dξ|c2|u2(t)u2|ttτ4[bIM(ε1)1+kIM(ε1)|S(ξ)S|+A(ε1)|II(ξ)|+c2IM(ε1)|u2u2(ξτ2)|]dξc2bIM(ε1)2(1+kIM(ε1))ttτ4((u2(t)u2)2+(S(ξ)S)2)dξ+c22A(ε1)ttτ4((u2(t)u2)2+(I(ξ)I)2)dξ+c22IM(ε1)2ttτ4((u2(t)u2)2+(u2(ξτ2)u2)2)dξ=c22[bIM(ε1)1+kIM(ε1)+A(ε1)+c2IM(ε1)]τ4(u2(t)u2)2+c2bIM(ε1)2(1+kIM(ε1))ttτ4(S(ξ)S)2dξ+c22A(ε1)ttτ4(I(ξ)I)2dξ+c22IM(ε1)2ttτ4(u2(ξτ2)u2)2dξ.
    $
    (4.11)

    For $ t > \tau_1+\tau_2+\tau_3+\tau_4+T_{1} $, let us define the following function,

    $ U = W_{1}+(1+kI^{*})U_1+U_2, $

    where

    $ U1=c1d12[ttτ1tθ(S(ξτ3)S)2dξdθ+τ1ttτ3(S(ξ)S)2dξ]+c1aS1(ε1)2ttτ3tθ(S(ξ)S)2dξdθ+c1bS1(ε1)2(1+kI)ttτ3tθ(II(ξ))2dξdθ,+c21S1(ε1)2[ttτ3tθ(u1(ξτ1)u1)2dξdθ+τ3ttτ1(u1(ξ)u1)2dξ]+c1e12ttτ1tθ(u1(ξ)u1)2dξdθ,
    $
    $ U2=c2d22[ttτ2tθ(I(ξτ4)I)2dξdθ+τ2ttτ4(I(ξ)I)2dξ]+c22A(ε1)ttτ4tθ(I(ξ)I)2dξdθ+c2bIM(ε1)2(1+kIM(ε1))ttτ4tθ(S(ξ)S)2dξdθ+c22IM(ε1)2[ttτ4tθ(u2(ξτ2)u2)2dξdθ+τ4ttτ2(u2(ξ)u2)2dξ]+c2e22ttτ2tθ(u2(ξ)u2)2dξdθ.
    $

    Computing the derivatives of $ U_1 $ and $ U_2 $, we easily get that, for $ t > \tau_1+\tau_2+\tau_3+\tau_4+T_{1} $,

    $ dU1dt=[c1d12τ1+c1aS1(ε1)2τ3](S(t)S)2c1d12ttτ1(S(ξτ3)S)2dξc1aS1(ε1)2ttτ3(S(ξ)S)2dξ+c1bS1(ε1)2(1+kI)τ3(I(t)I)2c1bS1(ε1)2(1+kI)ttτ3(I(ξ)I)2dξ+[c1e12τ1+c21S1(ε1)2τ3](u1(t)u1)2c21S1(ε1)2ttτ3(u1(ξτ1)u1)2dξc1e12ttτ1(u1(ξ)u1)2dξ,
    $
    (4.12)
    $ dU2dt=c2bIM(ε1)2(1+kIM(ε1))τ4(S(t)S)2c2bIM(ε1)2(1+kIM(ε1))ttτ4(S(ξ)S)2dξ+[c2d22τ2+c22A(ε1)τ4](I(t)I)2c2d22ttτ2(I(ξτ4)I)2dξc22A(ε1)ttτ4(I(ξ)I)2dξ+[c2e22τ2+c22IM(ε1)2τ4](u2(t)u2)2c22IM(ε1)2ttτ4(u2(ξτ2)u2)2dξc2e22ttτ2(u2(ξ)u2)2dξ.
    $
    (4.13)

    In summary, by (4.7)–(4.13), we have that, for $ t > \tau_1+\tau_2+\tau_3+\tau_4+T_{1} $,

    $ dUdt=dW1dt+(1+kI)dU1dt+dU2dt(1+kI)[a(2d1+e1)c12τ1c1aS1(ε1)2τ3c2bIM(ε1)2(1+kI)(1+kIM(ε1))τ4](S(t)S)2[f+bkS(1+kIM(ε1))(1+kI)(2d2+e2)c22τ2c1bS1(ε1)2τ3c22A(ε1)τ4](I(t)I)2(1+kI)c1[e1d1e12τ1S1(ε1)2(a+b1+kI+2c1)τ3](u1(t)u1)2c2[e2d2e22τ212(bIM(ε1)1+kIM(ε1)+A(ε1)+2c2IM(ε1))τ4](u2(t)u2)2:=(1+kI)Q1(ε1)(S(t)S)2Q2(ε1)(I(t)I)2(1+kI)c1Q3(ε1)(u1(t)u1)2c2Q4(ε1)(u2(t)u2)2.
    $
    (4.14)

    Furthermore, from conditions $ (H_8) $–$ (H_{11}) $, we see that, for sufficiently small $ \varepsilon_{1} > 0 $, we have that $ Q_i(\varepsilon_{1}) > 0 \; (i = 1, 2, 3, 4) $. This shows that $ \frac{dU}{dt}\leq0 $ for $ t > \tau_1+\tau_2+\tau_3+\tau_4+T_{1} $. By similar arguments as in the proof of Theorem 3.4, we can have

    $ \lim\limits_{t\rightarrow+\infty}S(t) = S^{*}, \; \; \lim\limits_{t\rightarrow+\infty}I(t) = I^{*}, \; \; \lim\limits_{t\rightarrow+\infty}u_{1}(t) = u_{1}^{*}, \; \; \lim\limits_{t\rightarrow+\infty}u_{2}(t) = u_{2}^{*}. $

    This proves that the endemic equilibrium $ E^{*} $ is globally attractive.

    The proof is completed.

    Remak 4.2. If $ \tau_{i} = 0 \; (i=1,2,3,4

    ) $, then model (1.3) reduces into model (1.2). Clearly, Theorem 4.1 extends and improves Theorem 2 in [28]. Further, if $ k = 0 $, Theorem 4.1 also include Theorem 2.2 in [17] as a special case.

    In this paper, we consider the SI epidemic model (1.3) with two feedback control variables and four time delays. In biology, model (1.3) has more general biological significance. Then, by skillfully constructing appropriate Lyapunov functionals, and combining Lyapunov–LaSalle invariance principle and Barbalat's lemma, some sufficient conditions for global dynamics of the equilibria of model (1.3) are established. In the case of $ d_{1}c_{1}\leq ae_{1} $, Theorem 3.3 gives complete conclusion on global asymptotic stability of the disease-free equilibrium $ E_{0} $. In the case of $ d_{1}c_{1} > ae_{1} $, in Theorem 3.4, global attractivity of the disease-free equilibrium $ E_{0} $ is considered under conditions $ (H_1) $–$ (H_3) $. Note that, in the case, it has from Theorem 3.2 that local asymptotic stability of the disease-free equilibrium $ E_{0} $ also depends on the time delays $ \tau_1 $ and $ \tau_3 $. Hence, The set of conditions $ (H_1) $–$ (H_3) $ has certain rationality. Furthermore, as a special case, Theorems 3.2, 3.3 and 3.4 improve and generalize Theorem 1 in [28] and Theorem 2.1 in [17].

    We also establish some sufficient conditions for global dynamics of the endemic equilibrium $ E^{*} $ of model (1.3). Theorem 4.1 shows that, if condition $ (H_5) $, or $ (H_6) $, or $ (H_7) $ holds, the time delays $ \tau_1 $ and $ \tau_3 $, or $ \tau_2 $ and $ \tau_4 $, or $ \tau_i \; (i = 1, 2, 3, 4) $ are harmless for global asymptotic stability of the endemic equilibrium $ E^{*} $ of model (1.3). If condition $ (H_4) $ holds, i.e., $ \tau_i = 0 \; (i = 1, 2, 3, 4) $, we see that Theorem 4.1 includes Theorem 2 in [28] and Theorem 2.2 in [17] as a special case. Furthermore, in Theorem 4.2, under condition $ R_0 > 1 $, a set of sufficient conditions $ (H_8) $–$ (H_{11}) $ is obtained to ensure global attractivity of the endemic equilibrium $ E^{*} $ of model (1.3). Note that the subsequent numerical simulations imply that any one of time delays $ \tau_i \; (i = 1, 2, 3, 4) $ may destroy local asymptotic stability of the endemic equilibrium $ E^{*} $ of model (1.3) and result in the occurrence of periodic oscillations etc.. Hence, in the set of conditions $ (H_8) $–$ (H_{11}) $, it should be feasible to have some certain limits on the lengths of time delays $ \tau_i \; (i = 1, 2, 3, 4) $.

    In the following, let us give some numerical simulations to summarize the applications of Theorem 3.4, Corollary 3.1 and Theorem 4.2.

    Firstly, let us choose the values of a set of parameters as follows,

    $ r=1.2,a=1,b=0.4,k=1,c1=0.8,μ=0.25,f=1,c2=1,e1=0.6,d1=1.2,e2=1,d2=1.
    $
    (5.1)

    By computations, we have that $ 0.96 = d_1c_1 > ae_1 = 0.6 $, $ 0.48 = rb > a\mu = 0.25 $, $ R_{0}\approx0.738 < 1 $, $ E_{0} = (0.462, 0, 0.923, 0) $ and $ \tau_{13}^{0}\approx4.542 $. Furthermore, conditions $ (H_1) $–$ (H_3) $ become the following inequalities,

    $ (\hat{H})\; \; \; 1.2\tau_1+0.48\tau_3 \lt 1,\; \; \; 0.3\tau_1+1.8\tau_3 \lt 0.5. $

    It has from Theorem 3.2 that the disease-free equilibrium $ E_{0} $ is locally asymptotically stable for $ \tau_{1}+\tau_{3} < \tau_{13}^{0}\approx4.542 $ and unstable for $ \tau_{1}+\tau_{3} > \tau_{13}^{0}\approx4.542 $. If time delays $ \tau_{1} $ and $ \tau_{3} $ satisfy more restrictive condition $ (\hat{H}) $, it has from Theorem 3.4 or Corollary 3.1 that the disease-free equilibrium $ E_{0} $ is also globally asymptotically stable.

    Let us choose $ \tau_1 = 0.75 $ and $ \tau_{3} = 0.15 $. We see that conditions $ \tau_{1}+\tau_{3} < \tau_{13}^{0}\approx4.542 $ and $ (\hat{H}) $ are satisfied. Figure 1(a) shows that the disease-free equilibrium $ E_{0} $ is globally asymptotically stable.

    Figure 1.  The phase trajectories and solution curves of the model (1.3) with $ R_{0}\approx0.738 < 1 $. (a) $ (H_1) $–$ (H_3) $ hold, and $ E_{0} $ is globally asymptotically stable. (b) $ (H_{1}) $–$ (H_{3}) $ do not hold, but $ \tau_1+\tau_3 < \tau_{13}^{0} $. $ E_{0} $ is locally asymptotically stable. (c) $ (H_{1}) $–$ (H_{3}) $ do not hold, but $ \tau_1+\tau_3 > \tau_{13}^{0} $ holds. $ E_{0} $ is unstable and periodic oscillations occur.

    Let us choose $ \tau_1 = 2 $ and $ \tau_{3} = 1.25 $. We see that condition $ (\hat{H}) $ does not hold, but condition $ \tau_{1}+\tau_{3} < \tau_{13}^{0}\approx4.542 $ is still satisfied. Figure 1(b) shows that the disease-free equilibrium $ E_{0} $ is locally asymptotically stable.

    Let us further choose $ \tau_1 = 2.75 $ and $ \tau_{3} = 2 $. We see that condition $ \tau_{1}+\tau_{3} > \tau_{13}^{0}\approx4.542 $ holds. Figure 1(c) shows that the disease-free equilibrium $ E_{0} $ becomes unstable and periodic oscillations occur. In Figure 1(a)–(c), for simplicity, $ \tau_2 $ and $ \tau_4 $ are fixed as $ \tau_2 = 3 $ and $ \tau_4 = 6 $, respectively.

    Secondly, let us give numerical simulations in the situation of the basic reproduction number $ R_{0} > 1 $.

    Let us choose the values of set of parameters as follows,

    $ r=2,a=0.8,b=1,k=1.5,c1=1.1,μ=0.25,f=0.8,c2=0.8,e1=0.8,d1=1,e2=0.5,d2=1.8.
    $
    (5.2)

    By computations, we have that $ R_{0}\approx3.678 > 1 $, $ E^{*}\approx(0.870, 0.130, \; 1.087, 0.467) $. Furthermore, conditions $ (H_8) $–$ (H_{11}) $ approximately become the following inequalities,

    $ (\tilde{H})\; \; \; \left\{ 1.54τ1+1.1τ3+0.132τ4<0.8,1.64τ2+1.375τ3+0.481τ4<1.246,0.4τ1+4.796τ3<0.8,0.25τ2+1.570τ4<0.278.
    \right. $

    If time delays $ \tau_i \; (i = 1, 2, 3, 4) $ satisfy condition $ (\tilde{H}) $, it has from Theorem 4.2 that the endemic equilibrium $ E^{*} $ is globally attractive.

    Let us choose $ \tau_1 = 0.4 $, $ \tau_{2} = 0.35 $, $ \tau_3 = 0.1 $ and $ \tau_{4} = 0.12 $. We see that condition $ (\tilde{H}) $ is satisfied. Figure 2(a) shows that the endemic equilibrium $ E^{*} $ is globally attractive.

    Figure 2.  The phase trajectories and solution curves of the model (1.3) with $ R_{0}\approx3.678 > 1 $. (a) $ (H_{8}) $–$ (H_{11}) $ hold, and $ E^{*} $ is globally attractive. (b) $ (H_{8}) $–$ (H_{11}) $ do not hold, but $ E^{*} $ may be still attractive. (c) $ (H_{8}) $–$ (H_{11}) $ do not hold, and $ E^{*} $ is unstable and periodic oscillations occur.

    Let us choose larger values here, $ \tau_1 = 1.6 $, $ \tau_{2} = 1.4 $, $ \tau_3 = 0.4 $ and $ \tau_{4} = 0.48 $. We see that condition $ (\tilde{H}) $ does not hold. Figure 2(b) show that the endemic equilibrium $ E^{*} $ may be still attractive.

    Let us further choose $ \tau_1 = 2 $, $ \tau_{2} = 1.75 $, $ \tau_3 = 0.6 $ and $ \tau_{4} = 0.8 $. We see that condition $ (\tilde{H}) $ does not hold. Figure 2(c) shows that the endemic equilibrium $ E^{*} $ becomes unstable and periodic oscillations occur.

    Moreover, Figure 3(a)–(d) show that any one of time delays $ \tau_{i} $ ($ i = 1, 2, 3, 4 $) can destroy local asymptotic stability of the endemic equilibrium $ E^{*} $ of model (1.3). Here, the values of the parameters of model (1.3) are the same as Eq (5.2).

    Figure 3.  The phase trajectories and solution curves of the model (1.3) with $ R_{0}\approx3.678 > 1 $.

    At the end of the paper, in view of Figure 1(b) and Figure 2(b) above, we would like to point out that conditions $ (H_1) $–$ (H_3) $ in Theorem 3.3 and conditions $ (H_8) $–$ (H_{11}) $ in Theorem 4.2 are actually conservative and worth of further improving. In addition, for global asymptotic stability of the endemic equilibrium $ E^{*} $ of model (1.3), in the case of $ d_{1}c_{1} > a_{1}e_{1} $ or $ d_{2}c_{2} > f e_{2} $, to give some sufficient conditions which are different from conditions $ (H_{8}) $–$ (H_{11}) $ may be also interesting, since Figure 4(a)–(d) show that model (1.3) may have richer dynamic behaviors. Here, the values of the parameters of model (1.3) are also the same as Eq (5.2). Further, it would be interesting to extend model (1.3) to the case of non-autonomous model and to consider the uniform persistence and existence of almost periodic solutions etc. [38,39].

    Figure 4.  The phase trajectories and solution curves of the model (1.3) with $ R_{0}\approx3.678 > 1 $.

    This paper is supported by Beijing Natural Science Foundation (No.1202019) and National Natural Science Foundation of China (No.11971055). The authors would like to thank the reviewers and the editor for their careful reading, helpful comments and suggestions that greatly improved the paper.

    All authors declare no conflicts of interest in this paper.

    [1] Duda RL (1998) Protein Chainmail: Catenated Protein in Viral Capsids. Cell 94: 55-60.
    [2] Wikoff WR, Liljas L, Duda RL, et al. (2000) Topologically linked protein rings in the bacteriophage HK97 capsid. Science 289: 2129-2133. doi: 10.1126/science.289.5487.2129
    [3] Akita F, Chong KT, Tanaka H, Y et al. (2007) The Crystal Structure of a Virus-like Particle from the Hyperthermophilic Archaeon Pyrococcus furiosus Provides Insight into the Evolution of Viruses. J Mol Biol 368: 1469-1483. doi: 10.1016/j.jmb.2007.02.075
    [4] Sutter M, Boehringer D, Gutmann S, et al. (2008) Structural basis of enzyme encapsulation into a bacterial nanocompartment. Nat Struct Mol Biol 15: 939-947. doi: 10.1038/nsmb.1473
    [5] Gelbart WM, Knobler CM (2009) Virology. Pressurized viruses. Science 323: 1682-1683.
    [6] Zhang X, Guo H, Jin L, et al. (2013) A new topology of the HK97-like fold revealed in Bordetella bacteriophage by cryoEM at 3.5 A resolution. eLife 2: e01299.
    [7] Zhou ZH, Hui WH, Shah S, et al. (2014) Four Levels of Hierarchical Organization, Including Noncovalent Chainmail, Brace the Mature Tumor Herpesvirus Capsid against Pressurization. Struct Lond Engl 1993.
    [8] Lander GC, Evilevitch A, Jeembaeva M, et al. (2008) Bacteriophage lambda stabilization by auxiliary protein gpD: timing, location, and mechanism of attachment determined by cryoEM. Struct Lond Engl 1993 16: 1399-1406.
    [9] Parent KN, Khayat R, Tu LH, et al. (2010) P22 coat protein structures reveal a novel mechanism for capsid maturation: stability without auxiliary proteins or chemical crosslinks. Struct Lond Engl 1993 18: 390-401.
    [10] Baker ML, Jiang W, Rixon FJ, et al. (2005) Common ancestry of herpesviruses and tailed DNA bacteriophages. J Virol 79: 14967-14970. doi: 10.1128/JVI.79.23.14967-14970.2005
    [11] Tso D, Hendrix RW, Duda RL (2014) Transient contacts on the exterior of the HK97 procapsid that are essential for capsid assembly. J Mol Biol 426: 2112-2129. doi: 10.1016/j.jmb.2014.03.009
    [12] Prevelige Jr PE (2008) Send for Reinforcements! Conserved Binding of Capsid Decoration Proteins. Structure 16: 1292-1293. doi: 10.1016/j.str.2008.08.003
    [13] Rizzo AA, Suhanovsky MM, Baker ML, et al. (2014). Multiple functional roles of the accessory I-domain of bacteriophage P22 coat protein revealed by NMR structure and CryoEM modeling. Struct. Lond Engl 1993 22: 830-841.
    [14] Chen D-H, Baker ML, Hryc CF, et al. (2011) Structural basis for scaffolding-mediated assembly and maturation of a dsDNA virus. Proc Natl Acad Sci 108: 1355-1360. doi: 10.1073/pnas.1015739108
    [15] Parent KN, Gilcrease EB, Casjens SR, et al. (2012) Structural evolution of the P22-like phages: Comparison of Sf6 and P22 procapsid and virion architectures. Virology 427: 177-188. doi: 10.1016/j.virol.2012.01.040
    [16] Parent KN, Tang J, Cardone G, et al. (2014). Three-dimensional reconstructions of the bacteriophage CUS-3 virion reveal a conserved coat protein I-domain but a distinct tailspike receptor-binding domain. Virology 464-465: 55-66.
    [17] Guo F, Liu Z, Fang P-A, et al. (2014) Capsid expansion mechanism of bacteriophage T7 revealed by multistate atomic models derived from cryo-EM reconstructions. Proc Natl Acad Sci 111: E4606-E4614. doi: 10.1073/pnas.1407020111
    [18] Fokine A, Leiman PG, Shneider MM, et al. (2005) Structural and functional similarities between the capsid proteins of bacteriophages T4 and HK97 point to a common ancestry. Proc Natl Acad Sci U S A 102: 7163-7168. doi: 10.1073/pnas.0502164102
    [19] Yang F, Forrer P, Dauter Z, et al. (2000) Novel fold and capsid-binding properties of the λ-phage display platform protein gpD. Nat Struct Mol Biol 7: 230-237. doi: 10.1038/73347
    [20] Baker ML, Hryc CF, Zhang Q, et al. (2013) Validated near-atomic resolution structure of bacteriophage epsilon15 derived from cryo-EM and modeling. Proc Natl Acad Sci U S A 110: 12301-12306. doi: 10.1073/pnas.1309947110
    [21] Iwai H, Forrer P, Plückthun A, et al. (2005) NMR solution structure of the monomeric form of the bacteriophage λ capsid stabilizing protein gpD. J Biomol NMR 31: 351-356. doi: 10.1007/s10858-005-0945-7
    [22] Morais MC, Choi KH, Koti JS, et al. (2005) Conservation of the Capsid Structure in Tailed dsDNA Bacteriophages: the Pseudoatomic Structure of ϕ29. Mol Cell 18: 149-159. doi: 10.1016/j.molcel.2005.03.013
    [23] Roizman B, Knipe DM, Whitley RJ (2007) Herpes simplex viruses. In Fields Virology, (Philadelphia: Lippincott-Williams & Wilkins), 2502-1601.
    [24] Mocarski ES, Shenk T, Pass RF (2007) Cytomegaloviruses. In Fields Virology, (Philadelphia: Lippincott-Williams & Wilkins), 2702-2772.
    [25] Chang Y, Cesarman E, Pessin MS, et al. (1994). Identification of herpesvirus-like DNA sequences in AIDS-associated Kaposi's sarcoma. Science 266: 1865-1869. doi: 10.1126/science.7997879
    [26] Ganem D (2007) Kaposi's sarcoma-associated herpesvirus. In Fields Virology, (Philadelphia: Lippincott-Williams & Wilkins), 2847-2888.
    [27] Rickinson AB, Kieff E (2007) Epstein-Barr Virus. In Fields Virology, (Philadelphia: Lippincott-Williams & Wilkins), 2656-2700.
    [28] Dai X, Gong D, Wu T-T, et al. (2014) Organization of capsid-associated tegument components in Kaposi's sarcoma-associated herpesvirus. J Virol 88: 12694-12702. doi: 10.1128/JVI.01509-14
    [29] Hui WH, Tang Q, Liu H, et al. (2013) Protein interactions in the murine cytomegalovirus capsid revealed by cryoEM. Protein Cell 4: 833-845. doi: 10.1007/s13238-013-3060-7
    [30] Forterre P, Krupovic M (2012) The Origin of Virions and Virocells: The Escape Hypothesis Revisited. In Viruses: Essential Agents of Life, G. Witzany, ed. (Springer Netherlands), 43-60.
    [31] Heinemann J, Maaty WS, Gauss GH, et al. (2011) Fossil record of an archaeal HK97-like provirus. Virology 417: 362-368. doi: 10.1016/j.virol.2011.06.019
    [32] Bujnicki JM (2002) Sequence permutations in the molecular evolution of DNA methyltransferases. BMC Evol Biol 2: 3. doi: 10.1186/1471-2148-2-3
    [33] Peisajovich SG, Rockah L, Tawfik DS (2006) Evolution of new protein topologies through multistep gene rearrangements. Nat Genet 38: 168-174. doi: 10.1038/ng1717
    [34] Vogel C, Morea V (2006) Duplication, divergence and formation of novel protein topologies. Bio Essays 28: 973-978.
  • This article has been cited by:

    1. Manuel De la Sen, Asier Ibeas, Santiago Alonso-Quesada, On the Supervision of a Saturated SIR Epidemic Model with Four Joint Control Actions for a Drastic Reduction in the Infection and the Susceptibility through Time, 2022, 19, 1660-4601, 1512, 10.3390/ijerph19031512
    2. Ke Guo, Wanbiao Ma, Global dynamics of a delayed air pollution dynamic model with saturated functional response and backward bifurcation, 2022, 21, 1534-0392, 3831, 10.3934/cpaa.2022124
    3. Meng Li, Ke Guo, Wanbiao Ma, Uniform Persistence and Global Attractivity in a Delayed Virus Dynamic Model with Apoptosis and Both Virus-to-Cell and Cell-to-Cell Infections, 2022, 10, 2227-7390, 975, 10.3390/math10060975
    4. Qing Yang, Xinhong Zhang, Daqing Jiang, Asymptotic behavior of a stochastic SIR model with general incidence rate and nonlinear Lévy jumps, 2022, 107, 0924-090X, 2975, 10.1007/s11071-021-07095-7
    5. Karrar Qahtan Al-Jubouri, Raid Kamel Naji, Arvind Kumar Misra, Dynamics Analysis of a Delayed Crimean‐Congo Hemorrhagic Fever Virus Model in Humans, 2024, 2024, 1110-757X, 10.1155/2024/4818840
    6. Luyao Zhao, Mou Li, Wanbiao Ma, Hopf bifurcation and stability analysis of a delay differential equation model for biodegradation of a class of microcystins, 2024, 9, 2473-6988, 18440, 10.3934/math.2024899
  • Reader Comments
  • © 2015 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(17418) PDF downloads(1324) Cited by(10)

Figures and Tables

Figures(8)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog