Research article

Threshold dynamics of an HIV-1 model with both viral and cellular infections, cell-mediated and humoral immune responses

  • Received: 14 April 2018 Accepted: 16 August 2018 Published: 13 December 2018
  • Human specific immunity consists of two branches: humoral immunity and cellular immunity. To protect us from pathogens, cell-mediated and humoral immune responses work together to provide the strongest degree of efficacy. In this paper, we propose an HIV-1 model with cell-mediated and humoral immune responses, in which both virus-to-cell infection and cell-to-cell transmission are considered. Five reproduction ratios, namely, immunity-inactivated reproduction ratio, cell-mediated immunity-activated reproduction ratio, humoral immunity-activated reproduction ratio, cell-mediated immunity-competed reproduction ratio and humoral immunity-competed reproduction ratio, are calculated and verified to be sharp thresholds determining the local and global properties of the virus model. Numerical simulations are carried out to illustrate the corresponding theoretical results and reveal the effects of some key parameters on viral dynamics.

    Citation: Jiazhe Lin, Rui Xu, Xiaohong Tian. Threshold dynamics of an HIV-1 model with both viral and cellular infections, cell-mediated and humoral immune responses[J]. Mathematical Biosciences and Engineering, 2019, 16(1): 292-319. doi: 10.3934/mbe.2019015

    Related Papers:

    [1] Ting Guo, Zhipeng Qiu . The effects of CTL immune response on HIV infection model with potent therapy, latently infected cells and cell-to-cell viral transmission. Mathematical Biosciences and Engineering, 2019, 16(6): 6822-6841. doi: 10.3934/mbe.2019341
    [2] Cameron Browne . Immune response in virus model structured by cell infection-age. Mathematical Biosciences and Engineering, 2016, 13(5): 887-909. doi: 10.3934/mbe.2016022
    [3] Tinevimbo Shiri, Winston Garira, Senelani D. Musekwa . A two-strain HIV-1 mathematical model to assess the effects of chemotherapy on disease parameters. Mathematical Biosciences and Engineering, 2005, 2(4): 811-832. doi: 10.3934/mbe.2005.2.811
    [4] Ran Zhang, Shengqiang Liu . Global dynamics of an age-structured within-host viral infection model with cell-to-cell transmission and general humoral immunity response. Mathematical Biosciences and Engineering, 2020, 17(2): 1450-1478. doi: 10.3934/mbe.2020075
    [5] Jiawei Deng, Ping Jiang, Hongying Shu . Viral infection dynamics with mitosis, intracellular delays and immune response. Mathematical Biosciences and Engineering, 2023, 20(2): 2937-2963. doi: 10.3934/mbe.2023139
    [6] A. M. Elaiw, N. H. AlShamrani, A. D. Hobiny . Stability of an adaptive immunity delayed HIV infection model with active and silent cell-to-cell spread. Mathematical Biosciences and Engineering, 2020, 17(6): 6401-6458. doi: 10.3934/mbe.2020337
    [7] Abdessamad Tridane, Yang Kuang . Modeling the interaction of cytotoxic T lymphocytes and influenza virus infected epithelial cells. Mathematical Biosciences and Engineering, 2010, 7(1): 171-185. doi: 10.3934/mbe.2010.7.171
    [8] Yan Wang, Tingting Zhao, Jun Liu . Viral dynamics of an HIV stochastic model with cell-to-cell infection, CTL immune response and distributed delays. Mathematical Biosciences and Engineering, 2019, 16(6): 7126-7154. doi: 10.3934/mbe.2019358
    [9] Chunyang Qin, Yuming Chen, Xia Wang . Global dynamics of a delayed diffusive virus infection model with cell-mediated immunity and cell-to-cell transmission. Mathematical Biosciences and Engineering, 2020, 17(5): 4678-4705. doi: 10.3934/mbe.2020257
    [10] Bing Li, Yuming Chen, Xuejuan Lu, Shengqiang Liu . A delayed HIV-1 model with virus waning term. Mathematical Biosciences and Engineering, 2016, 13(1): 135-157. doi: 10.3934/mbe.2016.13.135
  • Human specific immunity consists of two branches: humoral immunity and cellular immunity. To protect us from pathogens, cell-mediated and humoral immune responses work together to provide the strongest degree of efficacy. In this paper, we propose an HIV-1 model with cell-mediated and humoral immune responses, in which both virus-to-cell infection and cell-to-cell transmission are considered. Five reproduction ratios, namely, immunity-inactivated reproduction ratio, cell-mediated immunity-activated reproduction ratio, humoral immunity-activated reproduction ratio, cell-mediated immunity-competed reproduction ratio and humoral immunity-competed reproduction ratio, are calculated and verified to be sharp thresholds determining the local and global properties of the virus model. Numerical simulations are carried out to illustrate the corresponding theoretical results and reveal the effects of some key parameters on viral dynamics.


    The human immunodeficiency virus (HIV) is a lentivirus that causes HIV infection and over time acquired immunodeficiency syndrome (AIDS). AIDS leads to progressive failure of the immune system, which allows life-threatening opportunistic infections and cancers to thrive. In the past decades, within host virus models have been investigated in some literatures, which helps us understand the biological interactions between viruses and host cells. Nowak et al. [20] designed a mathematical model including uninfected cells x(t), infected cells y(t) and free virus v(t) to describe the viral dynamics in HIV-1 infection:

    ˙x(t)=sdx(t)β1x(t)v(t),˙y(t)=β1x(t)v(t)ay(t),˙v(t)=ky(t)uv(t), (1.1)

    where uninfected cells x(t) are produced at rate s and die at rate d; β1 is the infection rate of virus-to-cell infection; a is the death rate of infected cells; k denotes the number of free virus particles produced by per infected cell; u is the remove rate of virus. System (1.1) has been further investigated by Perelson and Nelson [21] and Cangelosi et al. [1].

    Faced with different virus infections, immunity system protects us against pathogens. Human specific immunity can be classified into cell-mediated immunity, for which the protective function is associated with cells and humoral immunity, where the protective function exists in the humor [2]. As for cell-mediated immunity, activated effector T cells can detect peptide antigens originating from various types of pathogens and remove virus-infected cells. Some HIV-1 infection models have been proposed to describe the virus dynamics with cell-mediated immune response (see, for example, [15,19,24,26,34]). While, in humoral immunity, matured B cells migrate from bone marrow to lymph nodes or other lymphatic organs, where they begin to eliminate pathogens [23]. There have been several works on virus models with humoral immune response (see, for example, [4,14,28,29,30]). In [6], Fouts et al. pointed out that a guiding principle for HIV vaccine design has been that cellular and humoral immunities work together to provide the strongest degree of efficacy. In [33], Yan and Wang considered both cell-mediated and humoral immune responses and put forward an HIV-1 infection model including both T cells and B cells, which only involves virus-to-cell infection mechanism.

    It is mentioned in [17] that cell-to-cell transmission is a more potent and efficient means of virus propagation than the virus-to-cell infection mechanism. Cell-to-cell spread not only facilitates rapid viral dissemination but also reduce the effectiveness of neutralizing antibodies and viral inhibitors by immune evasion. In [25], Sigal et al. proved that cell-to-cell spread of HIV-1 does reduce the efficacy of antiretroviral therapy, since cell-to-cell transmission can cause multiple infections of target cells, which can in turn reduce the sensitivity to the antiretroviral drugs. In view of this, some mathematical analysis of virus models with cell-to-cell transmission has been performed. For instance, Li and Wang [13] dealt with the global dynamics of an HIV infection model which incorporated direct cell-to-cell transmission. Meanwhile, Lai and Zou [11,12] studied the effect of cell-to-cell transfer of HIV-1 on the virus dynamics.

    It was assumed in system (1.1) that the infection process is governed by the mass-action principle, namely, the infection rate per host or per virus is a constant. In [22], Regoes et al. illustrated that the infection rate is often found to be a sigmoidal rather than a linear function of the parasite dose to which it is exposed, and presented a dose-dependent infection rate (v/ID50)κ/[1+(v/ID50)κ], where ID50 denotes the infectious dose at which 50% of the hosts are infected and κ measures the slope of the sigmoidal curve at ID50. In [10], Huang et al. indicated that the bilinear incidence rate is insufficient to describe the infection process in detail and proposed a class of nonlinear incidence. Besides, to place the model on more sound biological grounds, Xu [31] and Elaiw et al. [5] incorporated a saturation incidence β1v(t)/(1+αv(t)) to replace the mass-action infection rate.

    Motivated by the works of Fouts et al. [6], Yan and Wang [33], Sigal et al. [25] and Regoes et al. [22], in the present paper, we are concerned with the effects of cell-to-cell transmission, saturation incidence, both cell-mediated and humoral immune responses on the global dynamics of HIV-1 infection model. To this end, we consider the following delay differential equations:

    ˙x(t)=sdx(t)β1x(t)v(t)1+αv(t)β2x(t)y(t),˙y(t)=β1emτx(tτ)v(tτ)1+αv(tτ)+β2emτx(tτ)y(tτ)ay(t)p1y(t)z(t),˙v(t)=ky(t)uv(t)p2v(t)w(t),˙z(t)=c1y(t)z(t)b1z(t),˙w(t)=c2v(t)w(t)b2w(t), (1.2)

    where x(t), y(t), v(t), z(t), w(t) denote the concentration of uninfected cells, infected cells, virus, T cells and B cells at time t, respectively, and other parameters are described in Table 1. A simple schematic diagram for the virus infection corresponding to system (1.2) is shown in Figure 1.

    Table 1.  Definitions of frequently used symbols.
    SymbolsDescription
    srecruitment rate of uninfected cells
    ddeath rate of uninfected cells
    β1infection rate of virus-to-cell infection
    β2transmission rate of cell-to-cell transmission
    αsaturation infection rate coefficient
    τthe time between viral entry into a cell and the production of new
    free virus or the time between infected cells spreading virus into
    uninfected cells and the production of new free virus [8]
    emτthe probability of surviving the time period from tτ to t
    adeath rate of infected cells
    uremoval rate of virus
    kaverage number of free virus particles produced by per infected cell
    p1kill ratio of infected cells by T cells
    p2kill ratio of virus by B cells
    b1death rate of T cells
    b2death rate of B cells
    c1maturing rate of new T cells from thymocytes in the thymus
    c2production rate of new B cells by antigenic stimulation

     | Show Table
    DownLoad: CSV
    Figure 1.  Simple schematic diagram of the HIV-1 infection model. (a), (b), (c) and (d) depict the process of cell-mediated immunity, humoral immunity, cell-to-cell infection and virus-to-cell transmission, respectively.

    The initial condition for system (1.2) takes the form

    x(θ)=ϕ1(θ), y(θ)=ϕ2(θ), v(θ)=ϕ3(θ), z(θ)=ϕ4(θ), w(θ)=ϕ5(θ), (1.3)

    where it satisfies that ϕi(θ)0, θ[τ,0), ϕi(0)>0, where ϕiC([τ,0],R5+0), i=1,2,3,4,5, the Banach space of continuous functions mapping the interval [τ,0] into R5+0, where R5+0={(x1,x2,x3,x4,x5):xi0,i=1,2,3,4,5}.

    It can be proved by the fundamental theory of functional differential equations [7] that system (1.2) has a unique solution (x(t),y(t),v(t),z(t),w(t)) satisfying the initial condition (1.3). It is easy to show that all solutions of system (1.2) with initial condition (1.3) are defined on [0,+) and remain positive for all t0.

    This paper is organized as follows. In Section 2, we calculate the reproduction ratios to system (1.2) and establish the existence of feasible equilibria. In Section 3, the local asymptotic stability of each of feasible equilibria is studied. In Section 4, we investigate the global asymptotic stability of each of feasible equilibria. In Section 5, we present numerical simulations to illustrate the theoretical results and study the effects of cell-to-cell transmission, viral production rate, death rate of infected cells and viral removal rate on viral dynamics, respectively. Besides, we perform a sensitivity analysis of reproduction ratios. The paper ends with a conclusion in Section 6.

    Clearly, system (1.2) always has an infection-free equilibrium E0(s/d,0,0,0,0). Denote

    R0=(β1k+β2u)semτaud,

    where R0 is called immunity-inactivated reproduction ratio of system (1.2), which represents the number of newly infected cells produced by one infected cell during its lifespan [3]. It is easy to show that if R0>1, system (1.2) has an immunity-inactivated equilibrium E1(x1,y1,v1,0,0), where

    x1=s(u+αky1)(d+β2y1)(u+αky1)+β1ky1,v1=ky1u,

    and

    y1=(β1ak+β2au+αadkαβ2ksemτ)+Δ2αβ2ak,

    in which,

    Δ=(β1ak+β2au+αadkαβ2ksemτ)24αβ2ak(aduβ1ksemτβ2suemτ).

    Denote

    R1=c1semτ[β2(c1u+αb1k)+β1c1k]a[(c1d+β2b1)(c1u+αb1k)+β1b1c1k]=R01+X1R0,

    where

    X1=ab1(β1k+β2u)[c1(β1k+β2u)+αβ2b1k]+αβ1ab1c1dk2c1semτ(β1k+β2u)[c1(β1k+β2u)+αβ2b1k]>0.

    R1 is called cell-mediated immunity-activated reproduction ratio, which denotes the average number of T cells activated by infectious cells when virus infection is successful and humoral immune response has not been established. If R1>1, in addition to E0 and E1, system (1.2) has a cell-mediated immunity-activated equilibrium E2(x2,y2,v2,z2,0), where

    x2=c1s(c1u+αb1k)(c1d+β2b1)(c1u+αb1k)+β1b1c1k,y2=b1c1,v2=b1kc1u,

    and

    z2=c1semτ[β2b1(c1u+αb1k)+β1b1c1k]b1p1[(c1d+β2b1)(c1u+αb1k)+β1b1c1k]ab1b1p1.

    We further denote

    R2=c2ksemτ[β2u(c2+αb2)+β1c2k]au(c2+αb2)(β2b2u+c2dk)+β1ab2c2ku=R01+X2R0,

    in which

    X2=ab2u(β1k+β2u)[β2u(c2+αb2)+β1c2k]+αβ1ab2c2dk2uc2ksemτ(β1k+β2u)[β2u(c2+αb2)+β1c2k].

    R2 is called humoral immunity-activated reproduction ratio, which denotes the average number of B cells activated by viruses when virus infection is successful and cell-mediated response has not been established. When R2>1, system (1.2) has a humoral immunity-activated equilibrium E3(x3,y3,v3,0,w3), where

    x3=c2ks(c2+αb2)(c2+αb2)(β2b2p2w3+β2b2u+c2kd)+β1b2c2k,y3=b2p2c2kw3+b2uc2k,v3=b2c2,

    where w3 is the positive real root of the following quadratic equation:

    w23+(c2+αb2)(2β2ab2u+ac2dkβ2c2ksemτ)+β1ab2c2kβ2ab2p2(c2+αb2)w3+au(c2+αb2)(β2b2u+c2dk)+β1ab2c2kuβ2ab2p22(c2+αb2)(1R2)=0.

    Denote

    R3=b1c2kb2c1u,R4=c1semτ[β2b1(c2+αb2)+β1b2c1]ab1[(c1d+β2b1)(c2+αb2)+β1b2c1],

    where R3 is called humoral immunity-competed reproduction ratio and represents the average number of B cells activated by viruses under the condition that cell-mediated immune response has been established, while, R4 is called cell-mediated immunity-competed reproduction ratio and represents the average number of T cells activated by infectious cells under the condition that humoral immune response has been established. If R3>1 and R4>1, system (1.2) has an immunity-activated equilibrium E(x,y,v,z,w), where

    x=c1s(c2+αb2)(c1d+β2b1)(c2+αb2)+β1b2c1,y=b1c1,v=b2c2,w=b1c2kb2c1ub2c1p2,

    and

    z=c1semτ[β2b1(c2+αb2)+β1b2c1]ab1[(c1d+β2b1)(c2+αb2)+β1b2c1]b1p1[(c1d+β2b1)(c2+αb2)+β1b2c1],

    in which cell-mediated and humoral immune responses take effect simultaneously.

    In this section, we are concerned with the local asymptotic stability of each of feasible equilibria to system (1.2) by analyzing the distribution of roots of corresponding characteristic equations.

    Theorem 3.1. If R0<1, the infection-free equilibrium E0(s/d,0,0,0,0) of system (1.2) is locally asymptotically stable; if R0>1, E0 is unstable.

    Proof. The characteristic equation of system (1.2) at E0 is

    (λ+b1)(λ+b2)(λ+d)(λ+a)(λ+u)sde(λ+m)τ(λ+b1)(λ+b2)(λ+d)(β2λ+β1k+β2u)=0. (3.1)

    It is clear that (3.1) has negative real roots λ=b1, λ=b2, λ=d and other roots are determined by the following equation:

    (λ+a)(λ+u)sd(β2λ+β1k+β2u)e(λ+m)τ=0. (3.2)

    Denote R0=R01+R02, where

    R01=β1ksemτaudandR02=β2semτad.

    Substituting R0 and R02 into (3.2) yields

    (λa+1)(λu+1)=eλτ(λuR02+R0). (3.3)

    Now, we claim that all roots of (3.3) have negative real parts. Otherwise, there exists a root λ1=Reλ1+iImλ1 with Reλ10. In this case, if R0<1, it is easy to see that

    |λ1a+1||eλ1τ|,|λ1u+1|>|λ1uR02+R0|.

    It follows that

    |(λ1a+1)(λ1u+1)|>|eλ1τ(λ1uR02+R0)|,

    which contradicts to (3.3). Therefore, if R0<1, all roots of (3.1) have negative real parts and E0 is locally asymptotically stable. If R0>1, we denote the left side of (3.2) by G(λ):

    G(λ)=(λ+a)(λ+u)e(λ+m)τsd(β2λ+β1k+β2u), (3.4)

    where G(0)=au(1R0)<0 and G(λ) as λ. Noting that G(λ) is a continuous function in respect to λ, if R0>1, Eq. (3.1) has a positive real root, then E0 is unstable.

    Theorem 3.2. If R0>1, R1<1 and R2<1, the immunity-inactivated equilibrium E1(x1,y1,v1,0,0) of system (1.2) is locally asymptotically stable.

    Proof. The characteristic equation of system (1.2) at E1 is

    (λ+a)(λ+u)[λ(c1y1b1)][λ(c2v1b2)](λ+d+β1v11+αv1+β2y1)=e(λ+m)τβ2x1(λ+u)(λ+d)[λ(c1y1b1)][λ(c2v1b2)]+e(λ+m)τ(λ+d)[λ(c1y1b1)][λ(c2v1b2)]β1kx1(1+αv1)2. (3.5)

    Note that

    R1=H1(c1y1b1)+1<1, (3.6)

    in which

    H1=y1(1+αv1)[β2a(c1u+αb1k)+β1ac1k]+αβ1c1dkx1v1emτay1(1+αv1)[(c1d+β2b1)(c1u+αb1k)+β1b1c1k],

    and

    R2=H2(c2v1b2)+1<1, (3.7)

    where

    H2=y1(1+αv1)[β2au2(c2+αb2)+β1ac2ku]+αβ1c2dkux1v1emτy1(1+αv1)[au(c2+αb2)(β2b2u+c2dk)+β1ab2c2ku].

    It is clear that (3.5) has negative real roots λ=c1y1b1 and λ=c2v1b2, and other roots are determined by the following equation:

    (λ+a)(λ+u)(λ+d+β1v11+αv1+β2y1)e(λ+m)τ(λ+d)[β2x1(λ+u)+β1kx1(1+αv1)2]=0. (3.8)

    For the sake of contradiction, let λ2=Reλ2+iImλ2 with Reλ20. In this case, it is easy to see that

    |λ2+d+β1v11+αv1+β2y1|>|eλ2τ(λ2+d)|.

    Direct calculation shows that

    |(λ2+a)(λ2+u)||β2emτx1(λ2+u)+β1emτkx1(1+αv1)2|=λ2[λ2+u+β1emτkx1u(1+αv1)]+β1emτkx11+αv1β1emτkx1(1+αv1)2>λ2[λ2+u+β1emτkx1u(1+αv1)]>0,

    which contradicts to (3.8). Thus, if R0>1, R1<1 and R2<1, all roots of Eq. (3.5) have negative real parts, and E1 is locally asymptotically stable.

    Theorem 3.3. If R1>1 and R3<1, the cell-mediated immunity-activated equilibrium E2(x2,y2,v2, z2,0) of system (1.2) is locally asymptotically stable.

    Proof. The characteristic equation of system (1.2) at E2 is

    (λ+u)[λ(c2v2b2)][λ(λ+a+p1z2)+c1p1y2z2](λ+d+β1v21+αv2+β2y2)=e(λ+m)τ(λ+d)[λ(c2v2b2)][β2x2λ(λ+u)+β1kx2(1+αv2)2λ]. (3.9)

    Note that R3=(c2v2b2)/b2+1<1. It is clear that (3.9) has negative real root λ=c2v2b2, and other roots are determined by the following equation:

    (λ+u)[λ(λ+a+p1z2)+c1p1y2z2](λ+d+β1v21+αv2+β2y2)=e(λ+m)τ(λ+d)[β2x2λ(λ+u)+β1kx2(1+αv2)2λ]. (3.10)

    Similarly, we claim that all roots of (3.10) have negative real parts. Otherwise, there exists a root λ3=Reλ3+iImλ3 with Reλ30. In this case, it is obvious that

    |λ3+d+β1v21+αv2+β2y2|>|eλ3τ(λ3+d)|.

    It follows that

    |(λ3+u)[λ3(λ3+a+p1z2)+c1p1y2z2]||β2emτx2λ3(λ3+u)+β1emτkx2(1+αv2)2λ3|=λ23[λ3+u+β1emτx2v2y2(1+αv2)]+p1c1y2z2(λ3+u)+β1emτkx21+αv2λ3β1emτkx2(1+αv2)2λ3>λ23[λ3+u+β1emτx2v2y2(1+αv2)]+p1c1y2z2(λ3+u)>0,

    which contradicts to (3.10). Hence, if R1>1 and R3<1, all roots of Eq. (3.9) have negative real parts, and E2 is locally asymptotically stable.

    Theorem 3.4. If R2>1 and R4<1, the humoral immunity-activated equilibrium E3(x3,y3,v3,0,w3) of system (1.2) is locally asymptotically stable.

    Proof. The characteristic equation of system (1.2) at E3 is

    (λ+a)[λ(c1y3b1)][λ(λ+u+p2w3)+c2p2v3w3](λ+d+β1v31+αv3+β2y3)=e(λ+m)τx3(λ+d)[λ(c1y3b1)][β2λ(λ+u+p2w3)+β2c2p2v3w3+β1kλ(1+αv3)2]. (3.11)

    Note that

    R4=(c1y3b1)y3[β2ab1(c2+αb2)+β1ab2c1]+β1emτb2c1dx3ab1y3[(c1d+β2b1)(c2+αb2)+β1b2c1]+1<1. (3.12)

    It is obvious that (3.11) has negative real root λ=c1y3b1, and other roots are determined by the following equation:

    (λ+a)[λ(λ+u+p2w3)+c2p2v3w3](λ+d+β1v31+αv3+β2y3)=e(λ+m)τx3(λ+d)[β2λ(λ+u+p2w3)+β2c2p2v3w3+β1kλ(1+αv3)2]. (3.13)

    Similarly, we claim that all roots of (3.13) have negative real parts. If not, there exists a root λ4=Reλ4+iImλ4 with Reλ40. In this case, it is easy to see that

    |λ4+d+β1v31+αv3+β2y3|>|eλ4τ(λ4+d)|.

    Direct calculation yields

    |(λ4+a)[λ4(λ4+u+p2w3)+c2p2v3w3]||emτx3[β2λ4(λ4+u+p2w3)+β2c2p2v3w3+β1kλ4(1+αv3)2]|=λ4[λ4(λ4+u+p2w3)+c2p2v3w3]+β1emτx3v3y3(1+αv3)(λ24+c2p2v3w3)+β1emτx3v3y3(1+αv3)(u+p2w3)λ4β1emτkx3(1+αv3)2λ4>λ4[λ4(λ4+u+p2w3)+c2p2v3w3]+β1emτx3v3y3(1+αv3)(λ24+c2p2v3w3)>0,

    which contradicts to (3.13). Therefore, if R2>1 and R4<1, all roots of Eq. (3.11) have negative real parts, and E3 is locally asymptotically stable.

    Theorem 3.5. If R3>1 and R4>1, the immunity-activated equilibrium E(x,y,v,z,w) of system (1.2) is locally asymptotically stable.

    Proof. The characteristic equation of system (1.2) at E is

    (λ+d+β1v1+αv+β2y)[λ(λ+a+p1z)+c1p1yz][λ(λ+u+p2w)+c2p2vw]=e(λ+m)τx(λ+d){β2λ[λ(λ+u+p2w)+c2p2vw]+β1kλ2(1+αv)2}. (3.14)

    Similarly, we claim that all roots of (3.14) have negative real parts. Otherwise, there exists a root λ5=Reλ5+iImλ5 with Reλ50. In this case, it is clear that

    |λ5+d+β1v1+αv+β2y|>|eλ5τ(λ5+d)|.

    Direct calculation shows that

    |[λ5(λ5+a+p1z)+c1p1yz][λ5(λ5+u+p2w)+c2p2vw]||β2emτxλ5[λ5(λ5+u+p2w)+c2p2wv]+β1emτkxλ25(1+αv)2|=(λ25+c1p1yz)[λ5(λ5+u+p2w)+c2p2vw]+β1emτxvy(1+αv)λ5(λ25+c2p2vw)+β1emτxλ25y(1+αv)(uv+p2vwky1+αv)>(λ25+c1p1yz)[λ5(λ5+u+p2w)+c2p2vw]+β1emτxvy(1+αv)λ5(λ25+c2p2vw)>0,

    which contradicts to (3.14). Therefore, if R3>1 and R4>1, all roots of Eq. (3.14) have negative real parts, and E is locally asymptotically stable.

    In this section, we study the global stability of each of feasible equilibria to system (1.2) by suitable Lyapunov functionals and LaSalle's invariance principle. First, we discuss the boundedness of solutions.

    Lemma 4.1. Any solution of system (1.2) with initial condition (1.3) is bounded for all t0.

    Proof. Let (x(t),y(t),v(t),z(t),w(t)) be any solution of system (1.2) with initial condition (1.3). Denote

    B1(t)=x(tτ)+emτy(t)+p1c1emτz(t),B2(t)=v(t)+p2c2w(t).

    Calculating the derivatives of B1(t) and B2(t) in respect to t yields

    ˙B1(t)=sdx(tτ)aemτy(t)b1p1c1emτz(t)smin{a,b1,d}B1(t),

    and

    ˙B2(t)=y(t)uv(t)b2p2c2w(t)emτsmin{a,b1,d}min{b2,u}B2(t).

    Therefore, the following set is positively invariant set for system (1.2) :

    Ω={(x,y,v,z,w)|x+emτy+p1c1emτzsmin{a,b1,d},v+p2c2wemτsmin{a,b1,d}min{b2,u}}.

    It is easy to see that x(t),y(t),v(t),z(t) and w(t) are bounded in the invariant set Ω.

    Next, define a function g(x)=x1lnx, which will be used in Lyapunov functionals of this section.

    Theorem 4.2. If R0<1, the infection-free equilibrium E0(s/d,0,0,0,0) of system (1.2) is globally asymptotically stable.

    Proof. Let (x(t),y(t),v(t),z(t),w(t)) be any positive solution of system (1.2) with initial condition (1.3). Define

    V1(t)=x0g(x(t)x0)+l11y(t)+l12v(t)+l13z(t)+l14w(t)+ttτ(β1x(s)v(s)1+αv(s)+β2x(s)y(s))ds,

    where x0=s/d, and constants l11, l12, l13, l14 will be determined later. Calculating the derivative of V1(t) along positive solutions of system (1.2) yields

    ˙V1(t)=(1x0x(t))(sdx(t)β1x(t)v(t)1+αv(t)β2x(t)y(t))+l11(β1emτx(tτ)v(tτ)1+αv(tτ)+β2emτx(tτ)y(tτ)ay(t)p1y(t)z(t))+l12(ky(t)uv(t)p2v(t)w(t))+l13(c1y(t)z(t)b1z(t))+l14(c2v(t)w(t)b2w(t))+β1x(t)v(t)1+αv(t)+β2x(t)y(t)β1x(tτ)v(tτ)1+αv(tτ)β2x(tτ)y(tτ).

    Direct calculation yields

    ˙V1(t)=dx0(2x0x(t)x(t)x0)+β1x0v(t)1+αv(t)l12uv(t)l13b1z(t)l14b2w(t)+(l11emτ1)(β1x(tτ)v(tτ)1+αv(tτ)+β2x(tτ)y(tτ))+(β2x0+l12kl11a)y(t)+(l13c1l11p1)y(t)z(t)+(l14c2l12p2)v(t)w(t). (4.1)

    Choose

    l11=emτ,l12=emτaβ2x0k>0,l13=emτp1c1,l14=p2emτaβ2x0c2k>0. (4.2)

    Thus, we obtain from (4.1) and (4.2) that

    ˙V1(t)dx0(2x0x(t)x(t)x0)+(R01)emτaukv(t)l13b1z(t)l14b2w(t).

    It follows that ˙V1(t)0 with equality holding if and only if x=x0,y=v=z=w=0. It can be verified that M0={E0}Ω is the largest invariant subset of {(x(t),y(t),v(t),z(t),w(t)):˙V1(t)=0}. Noting that if R0<1, E0 is locally asymptotically stable, thus we obtain the global asymptotic stability of E0 from LaSalle's invariance principle.

    Theorem 4.3. If R0>1, R1<1 and R2<1, the immunity-inactivated equilibrium E1(x1,y1,v1,0,0) of system (1.2) is globally asymptotically stable.

    Proof. Let (x(t),y(t),v(t),z(t),w(t)) be any positive solution of system (1.2) with initial condition (1.3). Define

    V2(t)=x1g(x(t)x1)+l21y1g(y(t)y1)+l22v1g(v(t)v1)+l23z(t)+l24w(t)+β1x1v11+αv1ttτg(x(s)v(s)(1+αv1)x1v1(1+αv(s)))ds+β2x1y1ttτg(x(s)y(s)x1y1)ds,

    where constants l21, l22, l23 and l24 will be determined later. Calculating the derivative of V2(t) along positive solutions of system (1.2), we have

    ˙V2(t)=(1x1x(t))(sdx(t)β1x(t)v(t)1+αv(t)β2x(t)y(t))+l21(1y1y(t))(β1emτx(tτ)v(tτ)1+αv(tτ)+β2emτx(tτ)y(tτ)ay(t)p1y(t)z(t))+l22(1v1v(t))(ky(t)uv(t)p2v(t)w(t))+l23(c1y(t)z(t)b1z(t))+l24(c2v(t)w(t)b2w(t))+β1x1v11+αv1[g(x(t)v(t)(1+αv1)x1v1(1+αv(t)))g(x(tτ)v(tτ)(1+αv1)x1v1(1+αv(tτ)))]+β2x1y1[g(x(t)y(t)x1y1)g(x(tτ)y(tτ)x1y1)]. (4.3)

    Substituting s=dx1+β1x1v1/(1+αv1)+β2x1y1, β1emτx1v1/(1+αv1)+β2emτx1y1=ay1, ky1=uv1 into (4.3) yields

    ˙V2(t)=dx1(2x1x(t)x(t)x1)+l21ay1+l22uv1l22v1uv1y1y(t)v(t)l22uv(t)+β1x1v11+αv1[1+v(t)(1+αv1)v1(1+αv(t))x1x(t)]l21emτβ1x1v11+αv1x(tτ)v(tτ)(1+αv1)y1x1v1(1+αv(tτ))y(t)+β1x1v11+αv1lnx(tτ)v(tτ)(1+αv(t))x(t)v(t)(1+αv(tτ))+β2x1y1(1x1x(t)l21emτx(tτ)y(tτ)x1y(t))+β2x1y1lnx(tτ)y(tτ)x(t)y(t)+(l21emτ1)(β1x(tτ)v(tτ)1+αv(tτ)+β2x(tτ)y(tτ))+(β2x1+l22kl21a)y(t)+(l21p1y1l23b1)z(t)+(l22p2v1l24b2)w(t)+(l23c1l21p1)y(t)z(t)+(l24c2l22p2)v(t)w(t). (4.4)

    Choose

    l21=emτ,l22=β1x1v1ky1(1+αv1),l23=emτp1c1,l24=β1p2x1v1c2ky1(1+αv1). (4.5)

    From (4.4) and (4.5), we obtain that

    ˙V2(t)=dx1(2x1x(t)x(t)x1)+emτp1c1y1b1c1z(t)+β1p2x1v1ky1(1+αv1)c2v1b2c2w(t)α(v(t)v1)2v1(1+αv1)(1+αv(t))g(1+αv(t)1+αv1)β1x1v11+αv1[g(x1x(t))+g(y(t)v1y1v(t))]β1x1v11+αv1g(x(tτ)v(tτ)(1+αv1)y1x1v1(1+αv(tτ))y(t))β2x1y1[g(x1x(t))+g(x(tτ)y(tτ)x1y(t))]. (4.6)

    From (3.6) and (3.7), we derive that c1y1<b1 and c2v1<b2. Since function g(x)=x1lnx is always positive except for x=1 where g(x)=0. It follows from (4.6) that ˙V2(t)0 with equality holding if and only if x=x1,y=y1,v=v1,z=w=0. It can be proved that M1={E1}Ω is the largest invariant subset of {(x(t),y(t),v(t),z(t),w(t)):˙V2(t)=0}. Noting that if R0>1, R1<1 and R2<1, E1 is locally asymptotically stable, hence we obtain the global asymptotic stability of E1 from LaSalle's invariance principle.

    Theorem 4.4. If R1>1 and R3<1, the cell-mediated immunity-activated equilibrium E2(x2,y2,v2, z2,0) of system (1.2) is globally asymptotically stable.

    Proof. Let (x(t),y(t),v(t),z(t),w(t)) be any positive solution of system (1.2) with initial condition (1.3). Define

    V3(t)=x2g(x(t)x2)+l31y2g(y(t)y2)+l32v2g(v(t)v2)+l33z2g(z(t)z2)+l34w(t)+β1x2v21+αv2ttτg(x(s)v(s)(1+αv2)x2v2(1+αv(s)))ds+β2x2y2ttτg(x(s)y(s)x2y2)ds,

    where constants l31, l32, l33 and l34 will be determined later. Calculating the derivative of V3(t) along positive solutions of system (1.2), we obtain that

    ˙V3(t)=(1x2x(t))(sdx(t)β1x(t)v(t)1+αv(t)β2x(t)y(t))+l31(1y2y(t))(β1emτx(tτ)v(tτ)1+αv(tτ)+β2emτx(tτ)y(tτ)ay(t)p1y(t)z(t))+l32(1v2v(t))(ky(t)uv(t)p2v(t)w(t))+l33(1z2z(t))(c1y(t)z(t)b1z(t))+l34(c2v(t)w(t)b2w(t))+β1x2v21+αv2[g(x(t)v(t)(1+αv2)x2v2(1+αv(t)))g(x(tτ)v(tτ)(1+αv2)x2v2(1+αv(tτ)))]+β2x2y2[g(x(t)y(t)x2y2)g(x(tτ)y(tτ)x2y2)]. (4.7)

    Substituting s=dx2+β1x2v2/(1+αv2)+β2x2y2, β1emτx2v2/(1+αv2)+β2emτx2y2=ay2+p1y2z2, ky2=uv2, c1y2z2=b1z2 into (4.7) yields

    ˙V3(t)=dx2(2x2x(t)x(t)x2)+l31ay2+l32uv2+l33b1z2l32v2uv2y2y(t)v(t)l32uv(t)+β1x2v21+αv2[1+v(t)(1+αv2)v2(1+αv(t))x2x(t)]l31emτβ1x2v21+αv2x(tτ)v(tτ)(1+αv2)y2x2v2(1+αv(tτ))y(t)+β1x2v21+αv2lnx(tτ)v(tτ)(1+αv(t))x(t)v(t)(1+αv(tτ))+β2x2y2(1x2x(t)l31emτx(tτ)y(tτ)x2y(t))+β2x2y2lnx(tτ)y(tτ)x(t)y(t)+(l31emτ1)(β1x(tτ)v(tτ)1+αv(tτ)+β2x(tτ)y(tτ))+(β2x2+l32kl33c1z2l31a)y(t)+(l32p2v2l34b2)w(t)+(l31p1y2l33b1)z(t)+(l33c1l31p1)y(t)z(t)+(l34c2l32p2)v(t)w(t). (4.8)

    Choose

    l31=emτ,l32=β1x2v2ky2(1+αv2),l33=emτp1c1,l34=β1p2x2v2c2ky2(1+αv2). (4.9)

    From (4.8) and (4.9), we obtain that

    ˙V3(t)=dx2(2x2x(t)x(t)x2)+β1p2x2v2ky2(1+αv2)c2v2b2c2w(t)α(v(t)v2)2v2(1+αv2)(1+αv(t))g(1+αv(t)1+αv2)β1x2v21+αv2[g(x2x(t))+g(y(t)v2y2v(t))]β1x2v21+αv2g(x(tτ)v(tτ)(1+αv2)y2x2v2(1+αv(tτ))y(t))β2x2y2[g(x2x(t))+g(x(tτ)y(tτ)x2y(t))]. (4.10)

    Noting that R3=(c2v2b2)/b2+1<1, it is clear that c2v2<b2. It follows from (4.10) that ˙V3(t)0 with equality holding if and only if x=x2,y=y2,v=v2,z=z2,w=0. It can be verified that M3={E2}Ω is the largest invariant subset of {(x(t),y(t),v(t),z(t),w(t)):˙V3(t)=0}. Noting that if R1>1 and R3<1, E2 is locally asymptotically stable, thus we obtain the global asymptotic stability of E2 from LaSalle's invariance principle.

    Theorem 4.5. If R2>1 and R4<1, the humoral immunity-activated equilibrium E3(x3,y3,v3,0,w3) of system (1.2) is globally asymptotically stable.

    Proof. Let (x(t),y(t),v(t),z(t),w(t)) be any positive solution of system (1.2) with initial condition (1.3). Define

    V4(t)=x3g(x(t)x3)+l41y3g(y(t)y3)+l42v3g(v(t)v3)+l43z(t)+l44w3g(w(t)w3)+β1x3v31+αv3ttτg(x(s)v(s)(1+αv3)x3v3(1+αv(s)))ds+β2x3y3ttτg(x(s)y(s)x3y3)ds,

    where constants l41, l42, l43 and l44 will be determined later. Calculating the derivative of V4(t) along positive solutions of system (1.2), we obtain that

    ˙V4(t)=(1x3x(t))(sdx(t)β1x(t)v(t)1+αv(t)β2x(t)y(t))+l41(1y3y(t))(β1emτx(tτ)v(tτ)1+αv(tτ)+β2emτx(tτ)y(tτ)ay(t)p1y(t)z(t))+l42(1v3v(t))(ky(t)uv(t)p2v(t)w(t))+l43(c1y(t)z(t)b1z(t))+l44(1w3w(t))(c2v(t)w(t)b2w(t))+β1x3v31+αv3[g(x(t)v(t)(1+αv3)x3v3(1+αv(t)))g(x(tτ)v(tτ)(1+αv3)x3v3(1+αv(tτ)))]+β2x3y3[g(x(t)y(t)x3y3)g(x(tτ)y(tτ)x3y3)]. (4.11)

    Substituting s=dx3+β1x3v3/(1+αv3)+β2x3y3, β1emτx3v3/(1+αv3)+β2emτx3y3=ay3, ky3=uv3+p2v3w3, c2v3w3=b2w3 into (4.11) yields

    ˙V4(t)=dx3(2x3x(t)x(t)x3)+l41ay3+l42uv3+l44b2w3l42kv3y(t)v(t)(l42u+l44c2w3)v(t)+β1x3v31+αv3[1+v(t)(1+αv3)v3(1+αv(t))x3x(t)]l41emτβ1x3v31+αv3x(tτ)v(tτ)(1+αv3)y3x3(1+αv(tτ))v3y(t)+β1x3v31+αv3lnx(tτ)v(tτ)(1+αv(t))x(t)v(t)(1+αv(tτ))+β2x3y3(1x3x(t)l41emτx(tτ)y(tτ)x3y(t))+β2x3y3lnx(tτ)y(tτ)x(t)y(t)+(l41emτ1)(β1x(tτ)v(tτ)1+αv(tτ)+β2x(tτ)y(tτ))+(β2x3+l42kl41a)y(t)+(l41p1y3l43b1)z(t)+(l42p2v3l44b2)w(t)+(l43c1l41p1)y(t)z(t)+(l44c2l42p2)v(t)w(t). (4.12)

    Choose

    l41=emτ,l42=β1x3v3ky3(1+αv3),l43=emτp1c1,l44=β1p2x3v3c2ky3(1+αv3). (4.13)

    It follows from (4.12) and (4.13) that

    ˙V4(t)=dx3(2x3x(t)x(t)x3)+emτp1c1y3b1c1z(t)α(v(t)v3)2v3(1+αv3)(1+αv(t))g(1+αv(t)1+αv3)β1x3v31+αv3[g(x3x(t))+g(y(t)v3y3v(t))+g(x(tτ)v(tτ)(1+αv3)y3x3v3(1+αv(tτ))y(t))]β2x3y3[g(x3x(t))+g(x(tτ)y(tτ)x3y(t))]. (4.14)

    According to (3.12), it is easy to see that c1y3<b1. It follows from (4.14) that ˙V4(t)0 with equality holding if and only if x=x3,y=y3,v=v3,z=0,w=w3. It can be proved that M4={E3}Ω is the largest invariant subset of {(x(t),y(t),v(t),z(t),w(t)):˙V4(t)=0}. Noting that if R2>1 and R4<1, E3 is locally asymptotically stable, hence we obtain the global asymptotic stability of E3 from LaSalle's invariance principle.

    Theorem 4.6. If R3>1 and R4>1, the immunity-activated equilibrium E(x,y,v,z,w) of system (1.2) is globally asymptotically stable.

    Proof. Let (x(t),y(t),v(t),z(t),w(t)) be any positive solution of system (1.2) with initial condition (1.3). Define

    V5(t)=xg(x(t)x)+l51yg(y(t)y)+l52vg(v(t)v)+l53zg(z(t)z)+l54wg(w(t)w)+β1xv1+αvttτg(x(s)v(s)(1+αv)xv(1+αv(s)))ds+β2xyttτg(x(s)y(s)xy)ds,

    where constants l51, l52, l53 and l54 will be determined later. Calculating the derivative of V5(t) along positive solutions of system (1.2), we have

    ˙V5(t)=(1xx(t))(sdx(t)β1x(t)v(t)1+αv(t)β2x(t)y(t))+l51(1yy(t))(β1emτx(tτ)v(tτ)1+αv(tτ)+β2emτx(tτ)y(tτ)ay(t)p1y(t)z(t))+l52(1vv(t))(ky(t)uv(t)p2v(t)w(t))+l53(1zz(t))(c1y(t)z(t)b1z(t))+l54(1ww(t))(c2v(t)w(t)b2w(t))+β1xv1+αv[g(x(t)v(t)(1+αv)xv(1+αv(t)))g(x(tτ)v(tτ)(1+αv)xv(1+αv(tτ)))]+β2xy[g(x(t)y(t)xy)g(x(tτ)y(tτ)xy)]. (4.15)

    Substituting s=dx+β1xv/(1+αv)+β2xy, β1emτxv/(1+αv)+β2emτxy=ay+p1yz, ky=uv+p2vw, c1yz=b1z, c2vw=b2w into (4.15) yields

    ˙V5(t)=dx(2xx(t)x(t)x)+β1xv1+αv[1+v(t)(1+αv)v(1+αv(t))xx(t)]+β1xv1+αv[lnx(tτ)v(tτ)(1+αv(t))x(t)v(t)(1+αv(tτ))l51emτx(tτ)v(tτ)y(1+αv)xvy(t)(1+αv(tτ))]+β2xy(1xx(t)l51emτx(tτ)y(tτ)xy(t)+lnx(tτ)y(tτ)x(t)y(t))+l51ay+l52uv+l53c1yz+l54c2vwl52vuv+p2vwyy(t)v(t)+(l51emτ1)(β1x(tτ)v(tτ)1+αv(tτ)+β2x(tτ)y(tτ))+(β2x+l52uv+p2vwyl51al53c1z)y(t)(l54c2w+l52u)v(t)+(l51p1l53c1)yz(t)+(l52p2l54c2)vw(t)+(l53c1l51p1)y(t)z(t)+(l54c2l52p2)v(t)w(t). (4.16)

    Choose

    l51=emτ,l52=β1x(1+αv)(u+p2w),l53=p1c1emτ,l54=β1p2xc2(1+αv)(u+p2w). (4.17)

    From (4.16) and (4.17), we can obtain that

    ˙V5(t)=dx(2xx(t)x(t)x)α(v(t)v)2v(1+αv)(1+αv(t))g(1+αv(t)1+αv)β1xv1+αv[g(xx(t))+g(y(t)vyv(t))+g(x(tτ)v(tτ)y(1+αv)xvy(t)(1+αv(tτ)))] (4.18)
    β2xy[g(xx(t))+g(x(tτ)y(tτ)xy(t))].

    It follows from (4.18) that ˙V5(t)0 with equality holding if and only if x=x,y=y,v=v,z=z,w=w. It can be verified that M5={E}Ω is the largest invariant subset of {(x(t),y(t),v(t),z(t),w(t)):˙V5(t)=0}. Noting that if R3>1 and R4>1, E is locally asymptotically stable, we therefore obtain the global asymptotic stability of E from LaSalle's invariance principle.

    In this section, we want to illustrate the theoretical results for system (1.2) by numerical simulations. Besides, we investigate the effects of cell-to-cell transmission, viral production rate, death rate of infected cells and viral remove rate on viral dynamics. Furthermore, sensitivity analysis is used to quantify the range of variables in reproduction ratios and identify the key factors giving rise to reproduction ratios, which can be helpful to design treatment strategies and provide insights on evaluating effective antiviral drug therapies.

    Following [18,26,27,32], we choose appropriate parameters and simulate each of feasible equilibria, respectively.

    Case 1: Corresponding parameters are listed in Case 1 of Table 2. The immunity-inactivated reproduction ratio is calculated as R0=0.5640<1. From Theorem 3.1, we derive that infection-free equilibrium E0 is locally asymptotically stable, which is illustrated in Figure 2.

    Table 2.  List of parameters.
    Parameters (units)Case 1Case 2Case 3Case4Case5Source
    s  (cellsml/day) 50 23 100 23 100Assumed
    d (/day) 0.0046 0.0065 0.0046 0.0046 0.0046 [26]
    β1 (mlvirion/day) 4.8×107 4.8×107 4.8×107 4.8×107 4.8×107 [26]
    β2 (mlvirion/day) 4.7×107 4.7×109 4.7×107 4.7×107 4.7×107 [26]
    α 0.01 0.0001 0.01 0.01 0.01Assumed
    m 1.39 1.39 1.39 1.39 1.39 [26]
    τ (day) 0.5 0.3 0.5 0.5 0.5 [26]
    a (/day) 0.015 0.032 0.008 0.01 0.008Assumed
    p1 (cellsml/day) 0.005 0.005 0.001 0.005 0.001 [27]
    k (cellsvirion/day) 1.1349 7.3 1.1349 11.349 11.349 [26]
    u (/day) 0.5 0.25 0.05 0.05 0.05 [26]
    p2 (μg/day) 0.01 0.01 0.01 0.01 0.01 [27]
    c1 (cellsml/day) 0.002 0.021 0.002 0.002 0.002Assumed
    b1 (/day) 0.12 0.25 0.02 0.12 0.02Assumed
    c2 (cellsvirion/day) 0.0006 0.0013 0.00013 0.0013 0.0013 [27]
    b2 (/day) 0.12 0.46 0.12 0.12 0.12Assumed

     | Show Table
    DownLoad: CSV
    Figure 2.  The temporal solutions of x(t), y(t), v(t), z(t) and w(t) versus t of system (1.2) where R0=0.5640<1.

    Case 2: Corresponding parameters are listed in Case 2 of Table 2. By simple computing, we obtain that R0=1.0217>1, R1=0.9635<1 and R2=0.9625<1. From Theorem 3.2, we derive that immunity-inactivated equilibrium E1 is locally asymptotically stable, which is in accord with Figure 3.

    Figure 3.  The temporal solutions of x(t), y(t), v(t), z(t) and w(t) versus t of system (1.2) where R0=1.0217>1, R1=0.9635<1 and R2=0.9625<1.

    Case 3: Corresponding parameters are listed in Case 3 of Table 2. Similarly, we obtain that R1=5.1140>1 and R3=0.2459<1. From Theorem 3.3, we derive that cell-mediated immunity-activated equilibrium E2 is locally asymptotically stable, which is illustrated in Figure 4.

    Figure 4.  The temporal solutions of x(t), y(t), v(t), z(t) and w(t) versus t of system (1.2) where R1=5.1140>1 and R3=0.2459<1.

    Case 4: Corresponding parameters are listed in Case 4 of Table 2. Likewise, we obtain that R2=14.1830>1 and R4=0.2108<1. From Theorem 3.4, we derive that humoral immunity-activated equilibrium E3 is locally asymptotically stable, which is in keeping with in Figure 5.

    Figure 5.  The temporal solutions of x(t), y(t), v(t), z(t) and w(t) versus t of system (1.2) where R2=14.1830>1 and R4=0.2108<1.

    Case 5: Corresponding parameters are listed in Case 5 of Table 2. By calculation, we obtain that R3=24.5895>1 and R4=3.7395>1. From Theorem 3.5, we derive that immunity-activated equilibrium E is locally asymptotically stable, which is consistent with observation in Figure 6.

    Figure 6.  The temporal solutions of x(t), y(t), v(t), z(t) and w(t) versus t of system (1.2) where R3=24.5895>1 and R4=3.7395>1.

    In order to investigate the effect of cell-to-cell transmission, we carry out some numerical simulations to show the contribution of cell-to-cell transmission during the whole infection. First, we let β2 as zero to compare the virus infection without cell-to-cell transmission with the infection which has both transmissions. Figure 7 (β2=0, β2=4.7×107) shows that cell-to-cell transmission is of benefit to HIV-1 transmission and the time to reach the peak level of virus is shorter. Then, we increase β2 to study the change of the peak level of infected cells and virus, and the time to reach the peak level. Figure 7 (β2=4.7×107, β2=4.7×106, β2=4.7×105) shows that infected cells and virus reach the peak level more quickly as β2 increases, meanwhile, the peak level become larger as β2 increases, too. Therefore, cell-to-cell transmission plays an important role in the whole virus infection.

    Figure 7.  The effect of β2 on the dynamical behavior of system (1.2).

    Viral production rate also has great influence on the dynamical behavior of the model. We set the viral production rate k as 11.349, 34.047, 68.094. In Figure 8, we observe that the time to reach the peak level of infected cells and virus becomes shorter as k increases, which means that larger viral production rate contributes to the viral infection. Meanwhile, T cells and B cells increase more quickly as k increases, especially, larger viral production rate can stimulate more B cells. Hence, the peak level of infected cells and virus decreases as k increases. In terms of the prevention and treatment of HIV, it implies that antiretroviral therapies, such as, reverse transcriptase inhibitors and protease inhibitors are effective methods to decrease k, namely, to inhibit virus reproduction.

    Figure 8.  The effect of k on the dynamical behavior of system (1.2).

    Usually, the death rate of infected cells is larger than the death rate of uninfected cells due to the fact that HIV infection can kill more host cells. We present some numerical simulations to study the effect of death rate of infected cells on the dynamical behavior of the model. We can observe from Figure 9 that, infected cells and virus increase more slowly as a increases, which indicates that increasing the death rate of infected cells can slow down the virus infection. Humoral immunity is mainly used to clear virus in our humor, so the viral remove rate has an effect on viral infection as well. Figure 10 implies that as the viral remove rate increases, infected cells and virus increase more slowly, which has the similar results to a. In the clinic treatment of HIV, promoting body's immune capacity contributes to increasing the death rate of infected cells and viral remove rate.

    Figure 9.  The effect of a on the dynamical behavior of system (1.2).
    Figure 10.  The effect of u on the dynamical behavior of system (1.2).

    Sensitivity analysis is used to quantify the range of variables in reproduction ratios and to identify the key factors giving rise to reproduction ratios. In [9,16], Latin hypercube sampling (LHS) is found to be a more efficient statistical sampling technique which has been introduced to the field of disease modelling.

    LHS allows an un-biased estimate of the reproduction ratios, with the advantage that it requires fewer samples than simple random sampling to achieve the same accuracy. For each parameter of reproduction ratios, a probability density function is defined based on experimental data and stratified into N equiprobable serial intervals. A single value is then selected randomly from every interval and this is done for every parameter. In this way, an parameter value from each sampling interval is used only once in the analysis but the entire parameter space is equitably sampled in an efficient manner. Distributions of the reproduction ratios can then be derived directly by running the model N times with each of the sampled parameter sets.

    In terms of the prevention and treatment of HIV, we pay more attention to antiretroviral therapies, which is directly related to viral production rate and viral remove rate. Figure 11 shows the scatter plots of R0, R1 and R2 in respect to k and u, which implies that k is a positively correlative variable with R0 and R2, while u is a negatively correlative variable. As for R1, we find that the correlation between k and R1 or u and R1 is not clear.

    Figure 11.  Scatter plots of R0, R1 and R2 in respect to k and u.

    In [16], Marino et al. mentioned that Partial Rank Correlation Coefficients (PRCCs) provide a measure of the strength of a linear association between the parameters and the reproduction ratios. Furthermore, PRCCs are useful for identifying the most important parameters. The positive or negative of PRCCs respectively denote the positive or negative correlation with the reproduction ratios, and the sizes of PRCCs measure the strength of the correlation. First, we investigate the immunity-inactivated reproduction ratio R0, as we can see in Figure 12, β1 and k are positively correlative variables with R0 while others are negatively correlative variables. In order of correlative strength, it goes: β1, d, a, k, u and β2. Similarly, we obtain the PRCCs of R1 and R2 (see Figure 13). Specially, we observe that k and u is weakly correlative in respect to R1, which accords with the scatter plots of R1.

    Figure 12.  Tornado plot of PRCCs in regard to R0.
    Figure 13.  Tornado plots of PRCCs in regard to R1 and R2.

    In this paper, we have considered an HIV-1 infection model to describe cell-to-cell transmission, saturation incidence, both cell-mediated and humoral immune responses. By a complete mathematical analysis, the threshold dynamics of the model is established and it can be fully determined by reproduction ratios. If R0<1, the infection-free equilibrium E0 is locally and globally asymptotically stable; if R0>1, R1<1 and R2<1, the immunity-inactivated equilibrium E1 is locally and globally asymptotically stable; if R1>1 and R3<1, the cell-mediated immunity-activated equilibrium E2 is locally and globally asymptotically stable; if R2>1 and R4<1, the humoral immunity-activated equilibrium E3 is locally and globally asymptotically stable; if R3>1 and R4>1, the immunity-activated equilibrium E is locally and globally asymptotically stable.

    Numerical simulations vividly illustrate our main results of stability analysis for system (1.2). Besides, we have investigated the effects of cell-to-cell transmission, viral production rate, death rate of infected cells and viral remove rate on viral dynamics. It is worth mentioning that as the infection rate of cell-to-cell transmission β2 increases, virus load rises quickly and largely, which implies that cell-to-cell transmission facilitates virus spread. Furthermore, we perform a sensitivity analysis of reproduction ratios, which implies some useful consequences on the prevention and treatment of HIV-1.

    It is easy to see that immunity-inactivated reproduction ratios R0 is the sum of the reproduction ratio determined by virus-to-cell infection, R01, and that determined by cell-to-cell transmission, R02. In other words, immunity-inactivated reproduction ratio R0 becomes larger when the model includes cell-to-cell transmission. Meanwhile, we find that our research includes some existing work. When β2=0 and α=0, our virus model is similar to the model in [33] and the immunity-inactivated reproduction ratio R0 reduces to R01. Based on the model in [33], Wang et al. [27] consider nonlinear incidence and continuous intracellular delay, which is similar to our model with β2=0 only. Besides, when we only consider one of the immune responses, our model reduces to the models in [14] and [26].

    This work was supported by the National Natural Science Foundation of China (Nos.11871316, 11801340, 11371368, 11331009), Shanxi Scientific Data Sharing Platform for Animal Diseases (201605D121014), and the Science and Technology Innovation Team of Shanxi Province (201605D131044-06).

    The authors declare there is no conflict of interest.



    [1] R. A. Cangelosi, E. J. Schwartz and D. J. Wollkind, A quasi-steady-state approximation to the basic target-cell-limited viral dynamics model with a non-cytopathic effect, Front. Microbiol., 9 (2018), 54.
    [2] J. Charles, T. Paul and W. Mark, Immunobiology, 5nd edition, Garland Science, New York, 2001.
    [3] P. van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48.
    [4] A. M. Elaiw and N. H. AlShamrani, Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal, Nonlinear Anal. RWA, 26 (2015), 161–190.
    [5] A. M. Elaiw, A. A. Raezah and K. Hattaf, Stability of HIV-1 infection with saturated virus-target and infected-target incidences and CTL immune response, Int. J. Biomath., 10 (2017), 1750070.
    [6] T. R. Fouts, K. Bagley, I. J. Prado, et. al., Balance of cellular and humoral immunity determines the level of protection by HIV vaccines in rhesus macaque models of HIV infection, Proc. Natl. Acad. Sci., 13 (2015), 992–999.
    [7] J. K. Hale and S. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993.
    [8] K. Hattaf and N. Yousfi, A class of delayed viral infection models with general incidence rate and adaptive immune response, Int. J. Dynam. Control, 4 (2016), 254.
    [9] A. Hoare, D. G. Regan and D. P. Wilson, Sampling and sensitivity analyses tools (SaSAT) for computational modelling, Theor. Biol. Med. Model., 5 (2008), 4.
    [10] G. Huang, Y. Takeuchi and W. Ma, Lyapunov functionals for delay differential equations model of viral infections, SIAM J. Appl. Math., 70 (2010), 2693–2708.
    [11] X. Lai and X. Zou, Modeling cell-to-cell spread of HIV-1 with logistic target cell growth, J. Math. Anal. Appl., 426 (2015), 563–584.
    [12] X. Lai and X. Zou, Modeling HIV-1 virus dynamics with both virus-to-cell infection and cell-tocell transmission, SIAM J. Appl. Math., 74 (2014), 898–917.
    [13] F. Li and J. Wang, Analysis of an HIV infection model with logistic target-cell growth and cellto- cell transmission, Chaos Soliton Fract., 81 (2015), 136–145.
    [14] J. Lin, R. Xu and X. Tian, Threshold dynamics of an HIV-1 virus model with both virus-to-cell and cell-to-cell transmissions, intracellular delay, and humoral immunity, Appl. Math. Comput., 315 (2017), 516–530.
    [15] C. Lv, L. Huang and Z. Yuan, Global stability for an HIV-1 infection model with Beddington- DeAngelis incidence rate and CTL immune response, Commun. Nonlinear Sci. Numer. Simulat., 19 (2014), 121–127.
    [16] S. Marino, I. B. Hogue and C. J. Ray, A methodology for performing global uncertainty and sensitivity analysis in systems biology, J. Theor. Biol., 254 (2008), 178–196.
    [17] N. Martin and Q. Sattentau, Cell-to-cell HIV-1 spread and its implications for immune evasion, Curr. Opin. HIV AIDS, 4 (2009), 143–149.
    [18] A. Murase, T. Sasaki and T. Kajiwara, Stability analysis of pathogen-immune interaction dynamics, J. Math. Biol., 51 (2005), 247–267.
    [19] Y. Nakata, Global dynamics of a cell mediated immunity in viral infection models with distributed delays, J. Math. Anal. Appl., 375 (2011), 14–27.
    [20] M. Nowak, S. Bonhoeffer, G. Shaw and R. May, Anti-viral drug treatment: Dynamics of resistance in free virus and infected cell populations, J. Theor. Biol., 184 (1997), 203–217.
    [21] A. S. Perelson and P. W. Nelson, Mathematical Analysis of HIV-1: Dynamics in Vivo, SIAM Review, 41 (1999), 3–44.
    [22] R. R. Regoes, D. Ebert and S. Bonhoeffer, Dose-dependent infection rates of parasites produce the Allee effect in epidemiology, Proc. R. Soc. Lond. Ser. B, 269 (2002), 271–279.
    [23] E. J. Schwartz, N. K. Vaidya, K. S. Dorman, S. Carpenter and R. H. Mealey, Dynamics of lentiviral infection in vivo in the absence of adaptive immune responses, Virology, 513 (2018), 108–113.
    [24] H. Shu, L. Wang and J. Watmough, Global stability of a nonlinear viral infection model with infinitely distributed intracellular delays and CTL immune responses, SIAM J. Appl. Math., 73 (2013), 1280–1302.
    [25] A. Sigal, J. T. Kim, A. B. Balazs, E. Dekel, A. Mayo, R. Milo and D. Baltimore, Cell-to-cell spread of HIV permits ongoing replication despite antiretroviral therapy, Nature, 477 (2011), 95–98.
    [26] J. Wang, M. Guo, X. Liu and Z. Zhao, Threshold dynamics of HIV-1 virus model with cell-to-cell transmission, cell-mediated immune responses and distributed delay, Appl. Math. Comput., 291 (2016), 149–161.
    [27] J. Wang, J. Pang, T.Kuniya and Y. Enatsu, Global threshold dynamics in a five-dimensional virus model with cell-mediated, humoral immune responses and distributed delays, Appl. Math. Comput., 241 (2014), 298–316.
    [28] S. Wang and D. Zou, Global stability of in-host viral models with humoral immunity and intracellular delays, Appl. Math. Model., 36 (2012), 1313–1322.
    [29] T. Wang, Z. Hu, F. Liao and W. Ma, Global stability analysis for delayed virus infection model with general incidence rate and humoral immunity, Math. Comput. Simulat., 89 (2013), 13–22.
    [30] T.Wang, Z. Hu and F. Liao, Stability and Hopf bifurcation for a virus infection model with delayed humoral immunity response, J. Math. Anal. Appl., 411 (2014), 63–74.
    [31] R. Xu, Global stability of an HIV-1 infection model with saturation infection and intracellular delay, J. Math. Anal. Appl., 375 (2011), 75–81.
    [32] J. Xu, Y. Geng and Y. Zhou, Global dynamics for an age-structured HIV virus infection model with cellular infection and antiretroviral therapy, Appl. Math. Comput., 305 (2017), 62–83.
    [33] Y. Yan and W. Wang, Global stability of a five-dimensional model with immune responses and delay, Discrete and Continuous Dynamical Systems - Series B, 17 (2012), 401–416.
    [34] H. Zhu and X. Zou, Dynamics of a HIV-1 infection model with cell-mediated immune response and intracellular delay, Discrete and Continuous Dynamical Systems - Series B, 12 (2009), 511– 524.
  • This article has been cited by:

    1. Yan Wang, Minmin Lu, Daqing Jiang, Viral dynamics of a latent HIV infection model with Beddington-DeAngelis incidence function, B-cell immune response and multiple delays, 2021, 18, 1551-0018, 274, 10.3934/mbe.2021014
    2. Ahmed Elaiw, Noura AlShamrani, Stability of a general adaptive immunity virus dynamics model with multistages of infected cells and two routes of infection, 2020, 43, 0170-4214, 1145, 10.1002/mma.5923
    3. N.H. AlShamrani, Stability of a general adaptive immunity HIV infection model with silent infected cell-to-cell spread, 2021, 09600779, 110422, 10.1016/j.chaos.2020.110422
    4. A. M. Elaiw, N. H. AlShamrani, A. D. Hobiny, I. A. Abbas, Global stability of an adaptive immunity HIV dynamics model with silent and active cell-to-cell transmissions, 2020, 10, 2158-3226, 085216, 10.1063/5.0017214
    5. Ahmed M. Elaiw, Matuka A. Alshaikh, Global stability of discrete pathogen infection model with humoral immunity and cell-to-cell transmission, 2020, 130, 09600779, 109458, 10.1016/j.chaos.2019.109458
    6. A. M. Elaiw, N. H. AlShamrani, STABILITY OF A DELAYED ADAPTIVE IMMUNITY HIV INFECTION MODEL WITH SILENT INFECTED CELLS AND CELLULAR INFECTION, 2021, 11, 2156-907X, 964, 10.11948/20200124
    7. Najmeh Akbari, Rasoul Asheghi, Optimal control of an HIV infection model with logistic growth, celluar and homural immune response, cure rate and cell-to-cell spread, 2022, 2022, 1687-2770, 10.1186/s13661-022-01586-1
    8. A. M. Elaiw, N. H. AlShamrani, E. Dahy, A. A. Abdellatif, Aeshah A. Raezah, Effect of Macrophages and Latent Reservoirs on the Dynamics of HTLV-I and HIV-1 Coinfection, 2023, 11, 2227-7390, 592, 10.3390/math11030592
    9. Najmeh Akbari, Rasoul Asheghi, Maryam Nasirian, Stability and Dynamic of HIV-1 Mathematical Model with Logistic Target Cell Growth, Treatment Rate, Cure Rate and Cell-to-cell Spread, 2022, 26, 1027-5487, 10.11650/tjm/211102
    10. Noura H. AlShamrani, Matuka A. Alshaikh, Ahmed M. Elaiw, Khalid Hattaf, Dynamics of HIV-1/HTLV-I Co-Infection Model with Humoral Immunity and Cellular Infection, 2022, 14, 1999-4915, 1719, 10.3390/v14081719
    11. Ahmed M. Elaiw, Abdualaziz K. Aljahdali, Aatef D. Hobiny, Discretization and Analysis of HIV-1 and HTLV-I Coinfection Model with Latent Reservoirs, 2023, 11, 2079-3197, 54, 10.3390/computation11030054
    12. Lili Yuan, Xiao-Dong Gao, Yufei Xia, Optimising the oil phases of aluminium hydrogel-stabilised emulsions for stable, safe and efficient vaccine adjuvant, 2022, 16, 2095-0179, 973, 10.1007/s11705-021-2123-1
    13. A. M. Elaiw, N. H. Alshamrani, E. Dahy, A. A. Abdellatif, Stability of within host HTLV-I/HIV-1 co-infection in the presence of macrophages, 2023, 16, 1793-5245, 10.1142/S1793524522500668
    14. Rui Xu, Chenwei Song, Dynamics of an HIV infection model with virus diffusion and latently infected cell activation, 2022, 67, 14681218, 103618, 10.1016/j.nonrwa.2022.103618
    15. Ahmed M. Elaiw, Abdulaziz K. Aljahdali, Aatef D. Hobiny, Dynamical Properties of Discrete-Time HTLV-I and HIV-1 within-Host Coinfection Model, 2023, 12, 2075-1680, 201, 10.3390/axioms12020201
    16. Noura H. AlShamrani, Reham H. Halawani, Ahmed M. Elaiw, Analysis of general HIV-1 infection models with weakened adaptive immunity and three transmission modalities, 2024, 106, 11100168, 101, 10.1016/j.aej.2024.06.033
    17. Chong Chen, Yinggao Zhou, Dynamic analysis of HIV model with a general incidence, CTLs immune response and intracellular delays, 2023, 212, 03784754, 159, 10.1016/j.matcom.2023.04.029
    18. Shan Hu, Ruiqi Liu, Sheng Hu, Jianfeng Ye, Sheng Chen, Xiaozhou Ye, Yun Wang, Preparation of OVA‐loaded γ‐PGA/CS nanocomposite delivery vector and its immune response, 2024, 141, 0021-8995, 10.1002/app.54949
    19. Ke Guo, Donghong Zhao, Zhaosheng Feng, Lyapunov functionals for a virus dynamic model with general monotonic incidence, two time delays, CTL and antibody immune responses, 2024, 158, 08939659, 109212, 10.1016/j.aml.2024.109212
    20. Ahmed M. Elaiw, Aeshah A. Raezah, Matuka A. Alshaikh, Global Dynamics of Viral Infection with Two Distinct Populations of Antibodies, 2023, 11, 2227-7390, 3138, 10.3390/math11143138
    21. Noura H. AlShamrani, Reham H. Halawani, Ahmed M. Elaiw, Effect of Impaired B-Cell and CTL Functions on HIV-1 Dynamics, 2023, 11, 2227-7390, 4385, 10.3390/math11204385
    22. A. M. Elaiw, E. A. Almohaimeed, A. D. Hobiny, Modeling the co-infection of HTLV-2 and HIV-1 in vivo, 2024, 32, 2688-1594, 6032, 10.3934/era.2024280
    23. Alberto Vegas Rodriguez, Nieves Velez de Mendizábal, Sandhya Girish, Iñaki F. Trocóniz, Justin S. Feigelman, Modeling the Interplay Between Viral and Immune Dynamics in HIV: A Review and Mrgsolve Implementation and Exploration, 2025, 18, 1752-8054, 10.1111/cts.70160
    24. A. M. Elaiw, E. A. Almohaimeed, Within-host dynamics of HTLV-2 and HIV-1 co-infection with delay, 2025, 19, 1751-3758, 10.1080/17513758.2025.2506536
  • Reader Comments
  • © 2019 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(5042) PDF downloads(741) Cited by(24)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog