Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js
Research article

Global stability of a diffusive humoral immunity viral infection model with time delays and two modes of transmission

  • Received: 22 April 2025 Revised: 07 June 2025 Accepted: 10 June 2025 Published: 19 June 2025
  • MSC : 34D40, 35K57, 92D30

  • In this paper, the dynamical behaviors for a diffusive and delayed viral infection model with two modes of transmission were investigated. The uninfected cells dynamics, two infection modes for both virus-to-cell infection and cell-to-cell transmission, and concentration dependence for the latently infected cells, infected cells, viruses, and B cells were modeled by seven general nonlinear functions along with some assumptions. The basic reproduction number was calculated and demonstrated the global properties of the virus model. The theoretical results were illustrated by numerical simulations.

    Citation: Hui Miao. Global stability of a diffusive humoral immunity viral infection model with time delays and two modes of transmission[J]. AIMS Mathematics, 2025, 10(6): 14122-14139. doi: 10.3934/math.2025636

    Related Papers:

    [1] Liang Hong, Jie Li, Libin Rong, Xia Wang . Global dynamics of a delayed model with cytokine-enhanced viral infection and cell-to-cell transmission. AIMS Mathematics, 2024, 9(6): 16280-16296. doi: 10.3934/math.2024788
    [2] Xin Jiang . Threshold dynamics of a general delayed HIV model with double transmission modes and latent viral infection. AIMS Mathematics, 2022, 7(2): 2456-2478. doi: 10.3934/math.2022138
    [3] Ahmed M. Elaiw, Ghadeer S. Alsaadi, Aatef D. Hobiny . Global co-dynamics of viral infections with saturated incidence. AIMS Mathematics, 2024, 9(6): 13770-13818. doi: 10.3934/math.2024671
    [4] Ru Meng, Yantao Luo, Tingting Zheng . Stability analysis for a HIV model with cell-to-cell transmission, two immune responses and induced apoptosis. AIMS Mathematics, 2024, 9(6): 14786-14806. doi: 10.3934/math.2024719
    [5] Mohammed H. Alharbi . HIV dynamics in a periodic environment with general transmission rates. AIMS Mathematics, 2024, 9(11): 31393-31413. doi: 10.3934/math.20241512
    [6] B. S. Alofi, S. A. Azoz . Stability of general pathogen dynamic models with two types of infectious transmission with immune impairment. AIMS Mathematics, 2021, 6(1): 114-140. doi: 10.3934/math.2021009
    [7] A. M. Elaiw, N. H. AlShamrani, A. D. Hobiny . Mathematical modeling of HIV/HTLV co-infection with CTL-mediated immunity. AIMS Mathematics, 2021, 6(2): 1634-1676. doi: 10.3934/math.2021098
    [8] Ahmed M. Elaiw, Amani S. Alsulami, Aatef D. Hobiny . Global properties of delayed models for SARS-CoV-2 infection mediated by ACE2 receptor with humoral immunity. AIMS Mathematics, 2024, 9(1): 1046-1087. doi: 10.3934/math.2024052
    [9] A. M. Elaiw, E. A. Almohaimeed, A. D. Hobiny . Stability of HHV-8 and HIV-1 co-infection model with latent reservoirs and multiple distributed delays. AIMS Mathematics, 2024, 9(7): 19195-19239. doi: 10.3934/math.2024936
    [10] E. A. Almohaimeed, A. M. Elaiw, A. D. Hobiny . Modeling HTLV-1 and HTLV-2 co-infection dynamics. AIMS Mathematics, 2025, 10(3): 5696-5730. doi: 10.3934/math.2025263
  • In this paper, the dynamical behaviors for a diffusive and delayed viral infection model with two modes of transmission were investigated. The uninfected cells dynamics, two infection modes for both virus-to-cell infection and cell-to-cell transmission, and concentration dependence for the latently infected cells, infected cells, viruses, and B cells were modeled by seven general nonlinear functions along with some assumptions. The basic reproduction number was calculated and demonstrated the global properties of the virus model. The theoretical results were illustrated by numerical simulations.



    During the process of viral infection, the interactions among uninfected cells, infected cells, virus, and immune responses play a crucial role in controlling the virus propagation and antiviral defence. Establishing a virus model and analyzing it can effectively predict disease development trends[1,2,3]. Wang et al. [4], Zhang et al. [5], Georgescu et al. [6], Yuan et al. [7], and Hattaf [8] proposed different nonlinear incidence rates describing the infection process in detail in order to comprehensively characterize biological systems, further explaining different biological phenomena in depth.

    It is mentioned in [9,10,11] that different modes of infection have varying impacts on the infection process, such as exhaustion of the immune system, organ damage, and increased antibiotic resistance. Komarova et al. [12], Sigal et al. [13], and Iwami et al. [14] studied the spread of HIV models with two transmission modes. As is known to all, latently infected cells are one of the main reasons why AIDS cannot be completely eradicated. Meanwhile, latently infected cells are not only unaffected by drugs, but can also be activated by antigens. Wang et al. [15] proposed an HIV latent infection model with cell-to-cell transmission, but the humoral immune response has been ignored. Shu et al. [16], Lai et al. [17], and Yang et al. [18] incorporated two modes of viral models without the latently infected cells and humoral immune response. Viral models including logistic growth, multi-stages, and cell-to-cell transmission were also analyzed to exhibit complex dynamic behavior [19,20].

    The immune system protects us from various virus infections. Mathematical modeling of virus infection dynamics is critical to the understanding of complex interaction between immune response and viral infection. In [21], Elaiw et al. considered humoral immunity virus models including latently infected cells, without involving cell-to-cell infection and diffusion. Meanwhile, Lin et al. [22] studied the global dynamics of an HIV infection model which incorporated the cell-to-cell transmission and adaptive immunity. The model presented in [22] has neglected the latently infected cells and diffusion.

    Spatial diffusion can be a specific drug for preventing and treating certain diseases, providing precise guidance on drug carriers. Wang et al. [23] proposed a delayed and diffusive model with linear incidence. Thus, cell mobility plays an important role in different virus infections. Many diffusion viral infection models were studied in [24,25,26]. However, to our knowledge, there are few works that simultaneously consider factors such as latently infected cells, time delays, diffusion, and humoral immune response.

    Given the above discussion, the diffusive and delayed latent virus infection model with humoral immunity is described by the following nonlinear system:

    T(x,t)t=d1ΔT(x,t)+n(T(x,t))π1(T(x,t),V(x,t))π2(T(x,t),G(x,t)),L(x,t)t=d2ΔL(x,t)+(1η)ea1τ1[π1(T(x,tτ1),V(x,tτ1))+π2(T(x,tτ1),G(x,tτ1))](μ+α)h1(L(x,t)),G(x,t)t=d3ΔG(x,t)+ηea1τ2[π1(T(x,tτ2),V(x,tτ2))+π2(T(x,tτ2),G(x,tτ2))]+αh1(L(x,t))σh2(G(x,t)),V(x,t)t=d4ΔV(x,t)+γh2(G(x,t))kh3(V(x,t))ρh3(V(x,t))h4(Z(x,t)),Z(x,t)t=d5ΔZ(x,t)+χ+δh3(V(x,t))h4(Z(x,t))βh4(Z(x,t)), (1.1)

    with initial conditions

    T(x,θ)=ϕ1(x,θ)0,L(x,θ)=ϕ2(x,θ)0,G(x,θ)=ϕ3(x,θ)0,V(x,θ)=ϕ4(x,θ)0,Z(x,θ)=ϕ5(x,θ)0,xˉΩ,θ[τ,0],τ=max{τ1,τ2}, (1.2)

    and homogeneous Neumann boundary conditions

    Tn=Ln=Gn=Vn=Zn=0,t>0,xΩ, (1.3)

    where Ω is a bounded domain in Rn with smooth boundary Ω, and n denotes the outward normal derivative on Ω. Δ is the Laplacian operator. di(i=1,2,3,4,5) are the diffusion coefficients. T(x,t), L(x,t), G(x,t), V(x,t), and Z(x,t) denote the concentration of uninfected cells, latently infected cells, infected cells, viruses, and B cells at position x and time t, respectively. η(0<η<1) is the probability that the uninfected cell will turn into an infected cell. α is the conversion rate. n(T) denotes the growth of the uninfected cells. μh1(L), σh2(G), kh3(V), and βh4(Z) are the death rates of the latently infected cells, infected cells, viruses, and B cells, which only depend on its concentration. γ is the production rate. π1(T,V) and π2(T,G) are the virus-to-cell and cell-to-cell incidence rates, respectively. Let ρh3(V)h4(Z) and δh3(V)h4(Z) be the neutralization rates of viruses and activation rate of B cells, respectively. ea1τ1 and ea1τ2 represent the probability of an infected cell surviving to the stage of τ1 and τ2, respectively. χ is the generation rate of B cells.

    Define

    π11(T)=limV0π1(T,V)V=π1(T,0)V,π21(T)=limG0π2(T,G)G=π2(T,0)G.

    In this paper, we first introduce the following assumptions:

    (A1) n(T) is continuously differentiable, and there exists T0>0 such that n(T0)=0 and n(T0)<0.

    (A2) πi(T,θ) is continuously differentiable; πi(T,θ)>0 for T(0,),θ(0,); πi(T,θ)=0 if and only if T=0 or θ=0. πi(T,θ)T>0 and πi(T,θ)θ>0, for all T>0 and θ>0,i=1,2. πi1(T)>0 and πi1(T)>0 for all T>0,i=1,2.

    (A3) hi is strictly increasing on [0,+), hi(0)=0, hi(0)=1, limθ0hi(θ)=, and there exists ϱi>0 such that hi(θ)ϱiθ for any θ0,i=1,2,3,4.

    (A4) π1(T,V)h3(V) is non-increasing with respect to V(0,+) and π2(T,G)h2(G) is non-increasing with respect to G(0,+).

    In this paper, the purpose is to investigate the dynamical properties of model (1.1). The organization of our paper is as follows: In Section 2, the basic properties of solutions and the existence of equilibria are discussed. In Section 3, the global stability is stated. In Section 4, the numerical simulations are presented to further illustrate the dynamical behavior of the model. Finally, we will give a conclusion.

    Let Y=C(¯Ω,R5) be the Banach space with the supremum norm. For τ0, define C=C([τ,0],Y), which is a Banach space of continuous functions from [τ,0] into Y with the norm and let \mathcal{C}_{+} = \mathcal{C}([-\tau, 0], \mathbb{Y}_{+}) with \mathbb{Y}_{+} = \mathcal{C}(\overline{\Omega}, \mathbb{R}^5_{+}) . We will say that \Phi\in\mathcal{C} if \Phi is a function from \overline{\Omega}\times[-\tau, 0] to \mathbb{R}^5 and is defined by \Phi(x, s) = \Phi(s)(x) . Also, for \zeta > 0 , a function \nu(\cdot): [-\tau, \zeta)\rightarrow \mathbb{Y} induces functions \nu_t\in C for t\in[0, \zeta) , defined by \nu_t(\kappa) = \nu(t+\kappa), \; \kappa\in[-\tau, 0] .

    Theorem 2.1. For any given initial condition \psi\in C satisfying (1.2), there exists a unique non-negative solution of model (1.1)(1.3) defined on \bar{\Omega}\times[0, +\infty) and this solution remains bounded for all t\geq0 .

    Proof: For any \psi = (\psi_1, \psi_2, \psi_3, \psi_4, \psi_5)^T\in C and x\in \overline{\Omega} , we define \mathbb{H} = (\mathbb{H}_1, \mathbb{H}_2, \mathbb{H}_3, \mathbb{H}_4, \mathbb{H}_5):C\rightarrow \mathbb{Y} by

    \begin{align*} \mathbb{H}_1(\psi)(x) = \;& n(\psi_1(x,0))-\pi_1(\psi_1(x,0),\psi_4(x,0))-\pi_2(\psi_1(x,0),\psi_3(x,0)),\\[5pt] \mathbb{H}_2(\psi)(x) = \;& (1-\eta) e ^{-a_1\tau_1}[\pi_1(\psi_1(x,-\tau_1),\psi_4(x,-\tau_1))\\[6pt] &+\pi_2(\psi_1(x,-\tau_1),\psi_3(x,-\tau_1))]-(\mu +\alpha) h_1(\psi_2(x,0)),\\[5pt] \mathbb{H}_3(\psi)(x) = \;& \eta e ^{-a_1\tau_2}[\pi_1(\psi_1(x,-\tau_2),\psi_4(x,-\tau_2)) \\[6pt]& +\pi_2(\psi_1(x,-\tau_2),\psi_3(x,-\tau_2))] +\alpha h_1(\psi_2(x,0))-\sigma h_2(\psi_3(x,0)),\\[5pt] \mathbb{H}_4(\psi)(x) = \;& \gamma h_2(\psi_3(x,0))-kh_3(\psi_4(x,0))-\rho h_3(\psi_4(x,0))h_4(\psi_5(x,0)),\\[5pt] \mathbb{H}_5(\psi)(x) = \;& \chi+\delta h_3(\psi_4(x,0))h_4(\psi_5(x,0))-\beta h_4(\psi_5(x,0)). \end{align*}

    After that, model (1.1)–(1.3) can be written as the following abstract functional differential equation:

    \begin{equation} \begin{array}{rl} W'(t) = & \mathbb{B}W+\mathbb{H}(W_t),\;t > 0,\\[5pt] W(0) = & \psi\in \mathbb{X}, \end{array} \end{equation} (2.1)

    where W = (T, L, G, V, Z)^T , \psi = (\psi_1, \psi_2, \psi_3, \psi_4, \psi_5)^T , and \mathbb{B}W = (d_1\Delta T, d_2\Delta L, d_3\Delta G, d_4\Delta V, d_5\Delta Z)^T . Obviously, \mathbb{H} is locally Lipschitz in \mathbb{Y} . From [27,28,29,30,31], we deduce that model (2.1) has a unique local solution on [0, T_{max}) , where T_{max} is the maximal existence time for the solution of model (2.1).

    It is obvious that a lower-solution of the model (1.1)–(1.3) is 0 = (0, 0, 0, 0, 0) . So, we have T(x, t)\geq0, \; L(x, t)\geq 0, \; G(x, t)\geq 0, \; V(x, t)\geq 0 , and Z(x, t)\geq 0 .

    From the first equation of model (1.1), we have T(t)\leq n(T(t))\leq m-\bar{m}T(t) , which gives \lim\limits_{t\rightarrow +\infty}\sup T(t)\leq { \frac{{m}}{{\bar{m}}}} . Let

    \begin{align*} \mathbb{G}_1(x,t) = \;&(1-\eta) e ^{-a_1\tau_1}T(x,t-\tau_1)+\eta e^{-a_1\tau_2}T(x,t-\tau_2) +L(x,t)+G(x,t) \\[5pt] & +{ \frac{{\sigma}}{{2\gamma}}}V(x,t)+{ \frac{{\sigma\rho}}{{2\gamma\delta}}}Z(x,t), \end{align*}

    and then, it can be obtained that

    \begin{align*} { \frac{{\partial \mathbb{G}_1(x,t)}}{{\partial t}}} \leq \;& (1-\eta) e ^{-a_1\tau_1}d_1\Delta T(x,t-\tau_1) +\eta e ^{-a_1\tau_2}d_1\Delta T(x,t-\tau_2) +d_2\Delta L(x,t) \\[5pt] & +d_3\Delta G(x,t)+{ \frac{{\sigma}}{{2\gamma}}}d_4\Delta V(x,t)+{ \frac{{\sigma\rho}}{{2\gamma\delta}}}d_5\Delta Z(x,t) \\[5pt] & +{ \frac{{\sigma\rho\chi}}{{2\gamma\delta}}} +\delta[(1-\eta) e ^{-a_1\tau_1}+\eta e ^{-a_1\tau_2}] -m_1\mathbb{G}_1(x,t) \\[5pt] = \;& (1-\eta) e ^{-a_1\tau_1}d_1\Delta T(x,t-\tau_1) +\eta e ^{-a_1\tau_2}d_1\Delta T(x,t-\tau_2) +d_2\Delta L(x,t) \\[5pt] & +d_3\Delta G(x,t)+{ \frac{{\sigma}}{{2\gamma}}}d_4\Delta V(x,t)+{ \frac{{\sigma\rho}}{{2\gamma\delta}}}d_5\Delta Z(x,t) \\[5pt] & +A-m_1\mathbb{G}_1(x,t), \end{align*}

    where

    \begin{align*} m_1& = \min\{\bar{m},\mu,{ \frac{{\sigma}}{{2}}},k,\beta\},\\[3pt] A& = { \frac{{\sigma\rho\chi}}{{2\gamma\delta}}} +\delta[(1-\eta) e ^{-a_1\tau_1}+\eta e ^{-a_1\tau_2}]. \end{align*}

    Therefore, \mathbb{G}_1(x, t)\leq\max\Big\{\frac{A}{m_1}, B\Big\}, where

    \begin{align*} B = &\max\limits_{x\in\overline{\Omega}} \Big\{(1-\eta) e ^{-a_1\tau_1}\psi_1(x,-\tau_1)+\eta e^{-a_1\tau_2}\psi_1(x,-\tau_2) +\psi_2(x,0) \\[3pt]& +\psi_3(x,0) +{ \frac{{\sigma}}{{2\gamma}}}\psi_4(x,0)+{ \frac{{\sigma\rho}}{{2\gamma\delta}}}\psi_5(x,0)\Big\}, \end{align*}

    for \forall (x, t)\in\overline{\Omega}\times[0, T_{max}) . Thus, (T(x, t), L(x, t), G(x, t), V(x, t), Z(x, t)) are bounded on \overline{\Omega}\times[0, T_{max}) . Therefore, by the standard theory for semilinear parabolic systems [32], we have T_{max} = +\infty .

    Next, we discuss the existence of equilibria of model (1.1). Inspired by the method in [33,34], we consider the infection and viral production, and define matrices \mathbb{F} and \mathbb{V} as

    \mathbb{F} = \left(\begin{array}{ccc} 0& (1-\eta) e^{-a_1\tau_1}\cdot \frac{\partial \pi_2(T_0,0)}{\partial G}& (1-\eta) e^{-a_1\tau_1}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}\\[5pt] 0 & \eta e^{-a_1\tau_2}\cdot \frac{\partial \pi_2(T_0,0)}{\partial G} & \eta e^{-a_1\tau_2}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}\\[5pt] 0 & 0 & 0 \end{array}\right),

    and

    \mathbb{V} = \left(\begin{array}{ccc} \mu+\alpha& 0& 0\\[5pt] -\alpha & \sigma & 0\\[5pt] 0 & -\gamma & k \end{array}\right).

    Thus, the basic reproductive number, R_0 , can be defined as the spectral radius of the next generation operator \mathbb{F}\mathbb{V}^{-1} , where

    \mathbb{FV}^{-1} = \left(\begin{array}{ccc} a_{11} & a_{12} & a_{13}\\[5pt] a_{21} & a_{22} & a_{23}\\[5pt] 0 & 0 & 0 \end{array}\right),

    where

    \begin{align*} a_{11}& = (1-\eta) e^{-a_1\tau_1}\frac{\alpha}{\mu+\alpha}\Big(\frac{1}{\sigma}\cdot\frac{\partial \pi_2(T_0,0)}{\partial G}+\frac{\gamma}{k\sigma}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}\Big),\\[5pt] a_{12}& = (1-\eta) e^{-a_1\tau_1}\Big(\frac{1}{\sigma}\cdot\frac{\partial \pi_2(T_0,0)}{\partial G}+\frac{\gamma}{k\sigma}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}\Big),\\[5pt] a_{13}& = \frac{(1-\eta) e^{-a_1\tau_1}}{k}\cdot\frac{\partial \pi_1(T_0,0)}{\partial V},\\[5pt] a_{21}& = \eta e^{-a_1\tau_2}\frac{\alpha}{\alpha+\mu}\Big(\frac{1}{\sigma}\cdot\frac{\partial \pi_2(T_0,0)}{\partial G}+\frac{\gamma}{k\sigma}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}\Big),\\[5pt] a_{22}& = \eta e^{-a_1\tau_2}\Big(\frac{1}{\sigma}\cdot\frac{\partial \pi_2(T_0,0)}{\partial G}+\frac{\gamma}{k\sigma}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}\Big),\\[5pt] a_{23}& = \frac{\eta e^{-a_1\tau_2}}{k}\cdot\frac{\partial \pi_1(T_0,0)}{\partial V}. \end{align*}

    Therefore,

    \begin{array}{rl} R_0 = & \Big[\eta e^{-a_1\tau_2}+(1-\eta ) e^{-a_1\tau_1}\frac{\alpha}{\alpha+\mu}\Big]\Big(\frac{\gamma}{k\sigma}\cdot \frac{\partial \pi_1(T_0,0)}{\partial V}+\frac{1}{\sigma}\cdot\frac{\partial \pi_2(T_0,0)}{\partial G}\Big), \end{array}

    which biologically describes the average number of secondary infections produced by one infected cell at the beginning of infection. In the above expression of R_0 , divided into parts as R_0 = R_{01}+R_{02} , where R_{01} = [\eta e^{-a_1\tau_2}+(1-\eta) e^{-a_1\tau_1}\frac{\alpha}{\alpha+\mu}]\cdot\frac{\gamma}{k\sigma}\cdot \frac{\partial \pi_1(T_0, 0)}{\partial V} is the basic reproduction number via the virus-to-cell infection and R_{02} = [\eta e^{-a_1\tau_2}+(1-\eta) e^{-a_1\tau_1}\frac{\alpha}{\alpha+\mu}]\cdot\frac{1}{\sigma}\cdot\frac{\partial \pi_2(T_0, 0)}{\partial G} is the basic reproduction number via the cell-to-cell transmission, respectively.

    To find the equilibria of model (1.1), we need to solve

    \begin{equation} \begin{array}{r l} & n(T)-\pi_1(T,V)-\pi_2(T,G) = 0,\\[6pt] & (1-\eta) e ^{-a_1\tau_1}[\pi_1(T,V) +\pi_2(T,G)]-(\mu +\alpha) h_1(L) = 0,\\[6pt] & \eta e ^{-a_1\tau_2}[\pi_1(T,V) +\pi_2(T,G)] +\alpha h_1(L)-\sigma h_2(G) = 0,\\[6pt] & \gamma h_2(G)-kh_3(V)-\rho h_3(V)h_4(Z) = 0,\\[6pt] & \chi+\delta h_3(V)h_4(Z)-\beta h_4(Z) = 0. \end{array} \end{equation} (2.2)

    When V = 0 , the second and fourth equations of (2.2) lead to G = 0 and L = 0 . From the first equation of (2.2), we obtain n(T) = 0\Rightarrow T = T_0. Solving Z from (2.2) yields \chi-\beta h_4(Z) = 0\Rightarrow Z_0 = h_4^{-1}({ \frac{{\chi}}{{\beta}}}). It always has an infection-free equilibrium E_0 = (T_0, 0, 0, 0, h_4^{-1}({ \frac{{\chi}}{{\beta}}})) .

    Now, we assume that there exists V_1\in\Big(0, h_3^{-1}({ \frac{{\beta}}{{\delta}}})\Big) , the fifth equation of (2.2) leads to Z_1 = h_4^{-1}\Big({ \frac{{\chi}}{{\beta-\delta h_3(V_1)}}}\Big) , and the fourth equation of (2.2) leads to G_1 = h_2^{-1}\Big({ \frac{{k}}{{\gamma}}}h_3(V_1)+{ \frac{{\rho}}{{\gamma}}}h_3(V_1)h_4(Z_1)\Big) .

    Define

    F(T) = n(T)-\pi_1(T,V_1)-\pi_2(T,G_1),

    and then, F(0) = n(0) > 0 and F(T_0) = n(T_0)-\pi_1(T_0, V_1)-\pi_2(T_0, G_1) = -\pi_1(T_0, V_1)-\pi_2(T_0, G_1) < 0 . According to (A_1) and (A_2) , F(T) is a strictly decreasing function of T , which implies that there exists a unique T_1\in(0, T_0) such that F(T_1) = 0 . From the second equation, we obtain L_1 = h_1^{-1}\Big({ \frac{{(1-\eta)e ^{-a_1\tau_1}(\pi_1(T_1, V_1) +\pi_2(T_1, G_1))}}{{\mu+\alpha}}}\Big) . Hence, model (1.1) has unique endemic equilibrium E_1 = (T_1, L_1, G_1, V_1, Z_1), where

    \begin{align*} &T_1\in(0,T_0),\; L_1 = h_1^{-1}\Big({ \frac{{(1-\eta)e ^{-a_1\tau_1}(\pi_1(T_1,V_1) +\pi_2(T_1,G_1))}}{{\mu+\alpha}}}\Big),\\[5pt] &G_1 = h_2^{-1}\Big({ \frac{{k}}{{\gamma}}}h_3(V_1)+{ \frac{{\rho}}{{\gamma}}}h_3(V_1)h_4(Z_1)\Big),\\[5pt] &V_1\in\Big(0,h_3^{-1}({ \frac{{\beta}}{{\delta}}})\Big),\; Z_1 = h_4^{-1}\Big({ \frac{{\chi}}{{\beta-\delta h_3(V_1)}}}\Big). \end{align*}

    In this section, the global stability of the equilibria E_0 and E_1 of model (1.1)-(1.3) will be investigated. Let H(\xi) = \xi-1-\ln \xi, \; \xi\in(0, +\infty) , and it is observed that H(\xi)\geq 0 , \xi > 0 . H(\xi) = 0\Leftrightarrow \xi = 1 . For convenience, for any solution (T(x, t), L(x, t), G(x, t), V(x, t), Z(x, t)) of model (1.1), we set

    \begin{align*} &T(x,t) = T,\; T(x,t-\tau_1) = T_{\tau_1},\;T(x,t-\tau_2) = T_{\tau_2},\; L(x,t) = L,\\[5pt] &G(x,t) = G,\; G(x,t-\tau_1) = G_{\tau_1},\;G(x,t-\tau_2) = G_{\tau_2},\;V(x,t) = V,\\[5pt] &V(x,t-\tau_1) = V_{\tau_1},\;V(x,t-\tau_2) = V_{\tau_2},\;Z(x,t) = Z. \end{align*}

    To state the global stability on E_0 , we need an additional assumption:

    (A_5) \lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\cdot \lim\limits_{G\rightarrow 0}\frac{\pi_2(T,G)}{G} \leq \frac{\pi_{21}(T)/\pi_{11}(T)} {\pi_{21}(T_0)/\pi_{11}(T_0)}\cdot\pi_{21}(T_0)\leq\pi_{21}(T_0)\,\mbox{for}\,T\leq T_0.

    Theorem 3.1. Assume that (A_1) (A_5) hold. If R_0 \leq 1 , then infection-free equilibrium E_0 is globally asymptotically stable.

    Proof: Define a Lyapunov functional

    \begin{align*} U_1(t) = & \int_{\Omega}^{}\bigg\{\frac{\alpha(1-\eta)e^{-a_1\tau_1} +\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(T-T_0-\int_{T_0}^{T(t)} \lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(\theta,V)}\,\mathrm{d}\theta \Big) \\[5pt] & +\frac{\alpha}{\mu+\alpha}L(t)+G(t)+\frac{\sigma(1-R_{02})}{\gamma}V(t) \\[5pt] & +\frac{\sigma\rho(1-R_{02})}{\gamma\delta} \Big(Z-Z_0-\int_{Z_0}^{Z(t)} \frac{h_4(Z_0)} {h_4(\theta)}\,\mathrm{d}\theta \Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha} \int_{-\tau_1}^{0} [\pi_1(T(t+\theta),V(t+\theta))+\pi_2(T(t+\theta),G(t+\theta))]\,\mathrm{d}\theta \\[5pt] & +\eta e^{-a_1\tau_2}\int_{-\tau_2}^{0} [\pi_1(T(t+\theta),V(t+\theta))+\pi_2(T(t+\theta),G(t+\theta))]\,\mathrm{d}\theta\bigg\} \,\mathrm{d}x. \end{align*}

    Calculating the derivative of U_1(t) along the positive solution of model (1.1), we obtain

    \begin{align*} { \frac{{dU_1(t)}}{{dt}}} = & \int_{\Omega}^{}\bigg\{\frac{\alpha(1-\eta)e^{-a_1\tau_1} +\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha}\Big(1-\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\Big)\frac{\partial T} {\partial t} \\[5pt] & +\frac{\alpha}{\mu+\alpha}\frac{\partial L} {\partial t}+\frac{\partial G} {\partial t}+\frac{\sigma(1-R_{02})}{\gamma}\frac{\partial V} {\partial t} +\frac{\sigma\rho(1-R_{02})}{\gamma\delta}\Big(1-\frac{h_4(Z_0)} {h_4(Z)}\Big)\frac{\partial Z} {\partial t} \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha} [\pi_1(T,V)+\pi_2(T,G)-\pi_1(T_{\tau_1},V_{\tau_1})-\pi_2(T_{\tau_1},G_{\tau_1})] \\[5pt] & +\eta e^{-a_1\tau_2} [\pi_1(T,V)+\pi_2(T,G) -\pi_1(T_{\tau_2},V_{\tau_2}) -\pi_2(T_{\tau_2},G_{\tau_2})]\bigg\}\,\mathrm{d}x, \\[5pt] = & \int_{\Omega}^{}\bigg\{\frac{\alpha(1-\eta)e^{-a_1\tau_1} +\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(1-\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\Big)(n(T)-n(T_0)) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} (\pi_1(T,V)+\pi_2(T,G))\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)} \\[5pt] & -R_{02}\sigma h_2(G) -\frac{\sigma(1-R_{02})}{\gamma}h_3(V)(k+\rho h_4(Z_0)) \\[5pt] & +\frac{\sigma\rho(1-R_{02})}{\gamma\delta}\Big(1-\frac{h_4(Z_0)}{h_4(Z)}\Big) (\beta h_4(Z_0)- \beta h_4(Z)) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(1-\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\Big)d_1\Delta T \\[5pt] & +\frac{\alpha}{\mu+\alpha}d_2\Delta L +d_3\Delta G+\frac{\sigma(1-R_{02})}{\gamma}d_4\Delta V \\[5pt] & +\frac{\sigma\rho(1-R_{02})}{\gamma\delta}\Big(1 -\frac{h_4(Z_0)}{h_4(Z)}\Big)d_5\Delta Z\bigg\} \,\mathrm{d}x. \end{align*}

    From assumptions (A_3) and (A_4) , we have

    \begin{align*} & \frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_1(T,V)\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)} \\[5pt] & -\frac{\sigma(1-R_{02})}{\gamma}(k+\rho h_4(Z_0))h_3(V) \\[5pt] = \;& \frac{\sigma(k+\rho h_4(Z_0))}{\gamma}(R_{0}-1)h_3(V), \end{align*}

    and

    \begin{align*} & \frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_2(T,G)\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}-R_{02}\sigma h_2(G) \\[5pt] = \;& \frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \bigg(\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\cdot \lim\limits_{G\rightarrow 0}\frac{\pi_2(T,G)}{G} \\[5pt] & -\frac{\partial\pi_2(T_0,0)}{\partial G}\bigg)G \\[5pt] \leq \;& 0. \end{align*}

    From assumptions (A_1) and (A_2) , we have

    \Big(1-\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\Big)(n(T)-n(T_0)) < 0.

    Moreover, by utilizing assumption (A_5) , we obtain

    \begin{align*} & \lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\cdot \lim\limits_{G\rightarrow 0}\frac{\pi_2(T,G)}{G} \leq \frac{\pi_{21}(T)/\pi_{11}(T)} {\pi_{21}(T_0)/\pi_{11}(T_0)}\cdot\pi_{21}(T_0)\leq\pi_{21}(T_0),\;\mbox{for}\; T\leq T_0. \end{align*}

    Therefore, we obtain

    \begin{align*} { \frac{{dU_1(t)}}{{dt}}} \leq& \int_{\Omega}^{}\bigg\{\frac{\sigma(k+\rho h_4(Z_0))}{\gamma}(R_{0}-1)h_3(V) -\frac{\sigma\rho\beta(1-R_{02})}{\gamma\delta}\frac{(h_4(Z)-h_4(Z_0))^2}{h_4(Z)} \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(1-\lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)} {\pi_1(T,V)}\Big)d_1\Delta T \\[5pt] & +\frac{\alpha}{\mu+\alpha}d_2\Delta L +d_3\Delta G+\frac{\sigma(1-R_{02})}{\gamma}d_4\Delta V \\[5pt] & +\frac{\sigma\rho(1-R_{02})}{\gamma\delta}\Big(1-\frac{h_4(Z_0)}{h_4(Z)}\Big)d_5\Delta Z\bigg\} \,\mathrm{d}x. \end{align*}

    Using the divergence theorem and the homogeneous Neumann boundary conditions, we get

    \begin{align*} & \int_{\Omega} \Delta T\,\mathrm{d}x = \int_{\partial\Omega} \frac{\partial T}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\frac{\Delta T}{\pi_1(T,V_1)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla T\parallel^2}{\pi_1^2(T,V_1)}\,\mathrm{d}x, \\[5pt] & \int_{\Omega} \Delta L\,\mathrm{d}x = \int_{\partial\Omega} \frac{\partial L}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\Delta G\,\mathrm{d}x = \int_{\partial\Omega} \frac{\partial G}{\partial\vec{n}}\,\mathrm{d}x = 0, \\[5pt] & \int_{\Omega} \Delta V\,\mathrm{d}x = \int_{\partial\Omega} \frac{\partial V}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\Delta Z\,\mathrm{d}x = \int_{\partial\Omega} \frac{\partial Z}{\partial\vec{n}}\,\mathrm{d}x = 0, \\[5pt] & \int_{\Omega}\frac{\Delta Z}{h_4(Z)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla Z\parallel^2}{h_4^2(Z)}\,\mathrm{d}x. \end{align*}

    Thus, we obtain

    \begin{align*} { \frac{{dU_1(t)}}{{dt}}} \leq& \int_{\Omega}^{}\bigg\{ \frac{\sigma(k+\rho h_4(Z_0))}{\gamma}(R_{0}-1)h_3(V) -\frac{\sigma\rho\beta(1-R_{02})}{\gamma\delta}\frac{(h_4(Z)-h_4(Z_0))^2}{h_4(Z)}\bigg\} \,\mathrm{d}x \\[5pt] & -\frac{(\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2})d_1}{\mu+\alpha} \lim\limits_{V\rightarrow 0}\frac{\pi_1(T_0,V)}{\pi_1(T,V)} \int_{\Omega} \frac{\parallel\nabla T\parallel^2}{\pi_1^2(T,V_1)}\,\mathrm{d}x \\[5pt] & -\frac{\sigma\rho(1-R_{02})d_5 h_4(Z_0)}{\gamma\delta}\int_{\Omega} \frac{\parallel\nabla Z\parallel^2}{h_4^2(Z)}\,\mathrm{d}x. \end{align*}

    Therefore, { \frac{{dU_1(t)}}{{dt}}}\leq0 . { \frac{{dU_1(t)}}{{dt}}} = 0\Leftrightarrow T = T_0 , L = 0 , G = 0 , V = 0 , and Z = Z_0 . By LaSalle's invariance principle [31], E_0 is globally asymptotically stable when R_0\leq1 .

    Assume that \pi_1(T, V) , \pi_2(T, G) , and h_3(V) satisfy

    \begin{array}{rl} (A_6) & \Big(\frac{\pi_1(T,V)}{\pi_1(T,V_1)}-\frac{h_3(V)}{h_3(V_1)}\Big) \Big(1-\frac{\pi_1(T,V_1)}{\pi_1(T,V)} \Big)\leq0,\\[8pt] & \Big(\frac{\pi_2(T,G)\pi_1(T_1,V_1)}{\pi_2(T_1,G_1)\pi_1(T,V_1)} -\frac{h_3(V)}{h_3(V_1)}\Big) \Big(1-\frac{\pi_1(T,V_1)\pi_2(T_1,G_1)}{\pi_1(T_1,V_1)\pi_2(T,G)} \Big)\leq0. \end{array}

    Theorem 3.2. If R_0 > 1 , and (A_1) - (A_6) hold, then the endemic equilibrium E_1 is globally asymptotically stable.

    Proof: Define a Lyapunov functional

    \begin{align*} U_2(t) = & \int_{\Omega}^{}\bigg\{ \frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(T-T_1-\int_{T_1}^{T(t)} \frac{\pi_1(T_1,V_1)} {\pi_1(\theta,V_1)}\,\mathrm{d}\theta \Big) \\[5pt] & +\frac{\alpha}{\mu+\alpha} \Big(L-L_1-\int_{L_1}^{L(t)}\frac{h_1(L_1)}{h_1(\theta)}\,\mathrm{d}\theta \Big) +\Big(G-G_1-\int_{G_1}^{G(t)}\frac{h_2(G_1)}{h_2(\theta)}\,\mathrm{d}\theta \Big) \\[5pt] & +\frac{\sigma}{\gamma}\Big(V-V_1-\int_{V_1}^{V(t)}\frac{h_3(V_1)}{h_3(\theta)}\,\mathrm{d}\theta \Big)+\frac{\sigma\rho}{\gamma\delta} \Big(Z-Z_1-\int_{Z_1}^{Z(t)}\frac{h_4(Z_1)}{h_4(\theta)}\,\mathrm{d}\theta \Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\pi_1(T_1,V_1) \int_{-\tau_1}^{0}H\Big(\frac{\pi_1(T(t+\theta),V(t+\theta))} {\pi_1(T_1,V_1)}\Big)\,\mathrm{d}\theta \\[5pt] & +\eta e^{-a_1\tau_2}\pi_1(T_1,V_1) \int_{-\tau_2}^{0}H\Big(\frac{\pi_1(T(t+\theta),V(t+\theta))} {\pi_1(T_1,V_1)}\Big)\,\mathrm{d}\theta \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\pi_2(T_1,G_1) \int_{-\tau_1}^{0}H\Big(\frac{\pi_2(T(t+\theta),G(t+\theta))} {\pi_2(T_1,G_1)}\Big)\,\mathrm{d}\theta \\[5pt] & +\eta e^{-a_1\tau_2}\pi_2(T_1,G_1) \int_{-\tau_2}^{0}H\Big(\frac{\pi_2(T(t+\theta),G(t+\theta))} {\pi_2(T_1,G_1)}\Big)\,\mathrm{d}\theta\bigg\} \,\mathrm{d}x. \end{align*}

    Calculating the derivative of U_2(t) along the positive solution of model (1.1), it follows that

    \begin{align*} { \frac{{dU_2(t)}}{{dt}}} = & \int_{\Omega}^{}\bigg\{ \frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(1-\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \Big)d_1\Delta T \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(1-\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \Big)\Big(n(T)-n(T_1)\Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_1(T_1,V_1) \Big(1-\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_2(T_1,G_1)\Big(1-\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_1(T,V)\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_2(T,G)\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \\[5pt] & +\frac{\alpha}{\mu+\alpha}\Big(1-\frac{h_1(L_1)}{h_1(L)} \Big)d_2\Delta L +\alpha h_1(L_1) +\Big(1-\frac{h_2(G_1)}{h_2(G)} \Big)d_3\Delta G \\[5pt] & -\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\frac{h_1(L_1)}{h_1(L)} (\pi_1(T_{\tau_1},V_{\tau_1})+\pi_2(T_{\tau_1},G_{\tau_1})) \\[5pt] & -\alpha\frac{h_2(G_1)}{h_2(G)}h_1(L)+\sigma h_2(G_1) +\frac{\sigma}{\gamma}\Big(1-\frac{h_3(V_1)}{h_3(V)} \Big)d_4\Delta V \\[5pt] & -\eta e^{-a_1\tau_2}\frac{h_2(G_1)}{h_2(G)} (\pi_1(T_{\tau_2},V_{\tau_2})+\pi_2(T_{\tau_2},G_{\tau_2})) -\frac{\sigma k}{\gamma}h_3(V) \\[5pt] & -\sigma\frac{h_2(G)h_3(V_1)}{h_3(V)} +\frac{\sigma k}{\gamma}h_3(V_1)+\frac{\sigma \rho}{\gamma}h_3(V_1)h_4(Z) \\[5pt] & +\frac{\sigma \rho}{\gamma\delta}\Big(1-\frac{h_4(Z_1)}{h_4(Z)} \Big)d_5\Delta Z +\frac{\sigma \rho}{\gamma\delta}\Big(1-\frac{h_4(Z_1)}{h_4(Z)} \Big) \Big(\beta h_4(Z_1)-\beta h_4(Z)\Big) \\[5pt] & -\frac{\sigma \rho}{\gamma}h_3(V_1)h_4(Z_1)-\frac{\sigma \rho}{\gamma}h_3(V)h_4(Z_1) +\frac{\sigma \rho}{\gamma}h_3(V_1)h_4(Z_1)\frac{h_4(Z_1)}{h_4(Z)} \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\pi_1(T_1,V_1) \ln\Big(\frac{\pi_1(T_{\tau_1},V_{\tau_1})}{\pi_1(T,V)} \Big) +\eta e^{-a_1\tau_2}\pi_1(T_1,V_1) \ln\Big(\frac{\pi_1(T_{\tau_2},V_{\tau_2})}{\pi_1(T,V)} \Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\pi_2(T_1,G_1) \ln\Big(\frac{\pi_2(T_{\tau_1},G_{\tau_1})}{\pi_2(T,G)} \Big) +\eta e^{-a_1\tau_2}\pi_2(T_1,G_1) \ln\Big(\frac{\pi_2(T_{\tau_2},G_{\tau_2})}{\pi_2(T,G)} \Big)\bigg\} \,\mathrm{d}x. \end{align*}

    By using

    \begin{align*} &(1-\eta)e^{-a_1\tau_1}(\pi_1(T_1,V_1)+\pi_2(T_1,G_1)) = (\mu+\alpha)h_1(L_1),\\[5pt] & \eta e^{-a_1\tau_2}(\pi_1(T_1,V_1)+\pi_2(T_1,G_1))+\alpha h_1(L_1) = \sigma h_2(G_1) ,\\[5pt] &\gamma h_2(G_1) = k h_3(V_1)+\rho h_3(V_1)h_4(Z_1), \end{align*}

    and by the divergence theorem,

    \begin{align*} \int_{\Omega}\Delta T\,\mathrm{d}x = & \int_{\partial\Omega} \frac{\partial T}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\frac{\Delta T}{\pi_1(T,V_1)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla T\parallel^2}{\pi_1^2(T,V_1)}\,\mathrm{d}x,\\[5pt] \int_{\Omega}\Delta L\,\mathrm{d}x = & \int_{\partial\Omega} \frac{\partial L}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\frac{\Delta L}{h_1(L)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla L\parallel^2}{h_1^2(L)}\,\mathrm{d}x,\\[5pt] \int_{\Omega}\Delta G\,\mathrm{d}x = & \int_{\partial\Omega} \frac{\partial G}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\frac{\Delta G}{h_2(G)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla G\parallel^2}{h_2^2(G)}\,\mathrm{d}x,\\[5pt] \int_{\Omega}\Delta V\,\mathrm{d}x = & \int_{\partial\Omega} \frac{\partial V}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\frac{\Delta C}{h_3(V)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla V\parallel^2}{h_3^2(V)}\,\mathrm{d}x,\\[5pt] \int_{\Omega}\Delta Z\,\mathrm{d}x = & \int_{\partial\Omega} \frac{\partial Z}{\partial\vec{n}}\,\mathrm{d}x = 0,\; \int_{\Omega}\frac{\Delta Z}{h_4(Z)}\,\mathrm{d}x = \int_{\Omega} \frac{\parallel\nabla Z\parallel^2}{h_4^2(Z)}\,\mathrm{d}x. \end{align*}

    Thus, we have

    \begin{align*} { \frac{{dU_2(t)}}{{dt}}} = & \int_{\Omega}^{}\bigg\{ \frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \Big(1-\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)} \Big)\Big(n(T)-n(T_1)\Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_1(T_1,V_1) \Big(\frac{\pi_1(T,V)}{\pi_1(T,V_1)}-\frac{h_3(V)}{h_3(V_1)}\Big) \Big(1-\frac{\pi_1(T,V_1)}{\pi_1(T,V)} \Big) \\[5pt] & +\frac{\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2}}{\mu+\alpha} \pi_2(T_1,G_1) \\[5pt] & \times \Big(\frac{\pi_2(T,G)\pi_1(T_1,V_1)}{\pi_2(T_1,G_1)\pi_1(T,V_1)} -\frac{h_3(V)}{h_3(V_1)}\Big) \Big(1-\frac{\pi_1(T,V_1)\pi_2(T_1,G_1)}{\pi_1(T_1,V_1)\pi_2(T,G)} \Big) \\[5pt] & -\frac{\sigma\rho\chi}{\gamma\delta}\frac{(h_4(Z)-h_4(Z_1))^2}{h_4(Z)h_4(Z_1)} -\frac{\sigma\rho}{\gamma}h_3(V_1)h_4(Z_1)\Big[2-\frac{h_4(Z)}{h_4(Z_1)} -\frac{h_4(Z_1)}{h_4(Z)} \Big] \\[5pt] & -\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\pi_1(T_1,V_1) \Big[H\Big(\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)}\Big) +H\Big(\frac{h_1(L)h_2(G_1)}{h_1(L_1)h_2(G)}\Big) \\[5pt] & +H\Big(\frac{h_2(G)h_3(V_1)}{h_2(G_1)h_3(V)}\Big) +H\Big(\frac{\pi_1(T,V_1)h_3(V)}{\pi_1(T,V)h_3(V_1)}\Big) +H\Big(\frac{h_1(L_1)\pi_1(T_{\tau_1},V_{\tau_1})} {h_1(L)\pi_1(T_1,V_1)}\Big) \Big] \\[5pt] & -\eta e^{-a_1\tau_2}\pi_1(T_1,V_1) \Big[H\Big(\frac{\pi_1(T,V)}{\pi_1(T,V_1)}\Big) +H\Big(\frac{h_2(G)h_3(V_1)}{h_2(G_1)h_3(V)}\Big) +H\Big(\frac{h_2(G_1)\pi_1(T_{\tau_2},V_{\tau_2})}{h_2(G)\pi_1(T_1,V_1)}\Big) \\[5pt] & +H\Big(\frac{\pi_1(T,V_1)h_3(V)}{\pi_1(T,V)h_3(V_1)}\Big)\Big] -\frac{\alpha(1-\eta)e^{-a_1\tau_1}}{\mu+\alpha}\pi_2(T_1,G_1) \Big[H\Big(\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)}\Big) +H\Big(\frac{h_1(L)h_2(G_1)}{h_1(L_1)h_2(G)}\Big) \\[5pt] & +H\Big(\frac{h_1(L_1)\pi_2(T_{\tau_1},G_{\tau_1})}{h_1(L)\pi_2(T_1,G_1)}\Big) +H\Big(\frac{h_2(G)h_3(V_1)}{h_2(G_1)h_3(V)}\Big) +H\Big(\frac{h_3(V)\pi_1(T,V_1)\pi_2(T_1,G_1)} {h_3(V_1)\pi_1(T_1,V_1)\pi_2(T,G)}\Big)\Big] \\[5pt] & -\eta e^{-a_1\tau_2}\pi_2(T_1,G_1) \Big[H\Big(\frac{\pi_1(T_1,V_1)}{\pi_1(T,V_1)}\Big) +H\Big(\frac{h_3(V)\pi_1(T,V_1)\pi_2(T_1,G_1)} {h_3(V_1)\pi_1(T_1,V_1)\pi_2(T,G)}\Big) \\[5pt] & +H\Big(\frac{h_2(G_1)\pi_2(T_{\tau_2},G_{\tau_2})}{h_2(G)\pi_2(T_1,G_1)}\Big) +H\Big(\frac{h_2(G)h_3(V_1)}{h_2(G_1)h_3(V)}\Big) \Big]\bigg\} \,\mathrm{d}x \\[5pt] & -\frac{(\alpha(1-\eta)e^{-a_1\tau_1}+\eta(\mu+\alpha)e^{-a_1\tau_2})d_1\pi_1(T_1,V_1)} {\mu+\alpha}\int_{\Omega} \frac{\parallel\nabla T\parallel^2}{\pi_1^2(T,V_1)}\,\mathrm{d}x \\[5pt] & -\frac{\alpha d_2 h_1(L_1)}{\mu+\alpha}\int_{\Omega} \frac{\parallel\nabla L\parallel^2}{h_1^2(L)}\,\mathrm{d}x -d_3 h_2(G_1)\int_{\Omega} \frac{\parallel\nabla G\parallel^2}{h_2^2(G)}\,\mathrm{d}x \\[5pt] & -\frac{\sigma d_4 h_3(V_1)}{\gamma}\int_{\Omega} \frac{\parallel\nabla V\parallel^2}{h_3^2(V)}\,\mathrm{d}x -\frac{\sigma \rho d_5 h_4(Z_1)}{\gamma\delta}\int_{\Omega} \frac{\parallel\nabla Z\parallel^2}{h_4^2(Z)}\,\mathrm{d}x. \end{align*}

    Hence, { \frac{{dU_2(t)}}{{dt}}}\leq 0 . { \frac{{dU_2(t)}}{{dt}}} = 0\Leftrightarrow T(t) = T_1, L(t) = L_1, G(t) = G_1 , V(t) = V_1 , and Z(t) = Z_1. From LaSalle's invariance principle [31], we have that E_1 is globally asymptotically stable when R_0 > 1 .

    In this section, we present several numerical examples to illustrate the results obtained in Section 3. We will use the finite difference scheme which is proposed in [35,36] for the delayed reaction-diffusion epidemic models. For convenience, we consider model (1.1) under the one-dimensional spatial domain \Omega = [0, 1] . The homogeneous Neumann boundary conditions and the initial conditions are given by

    \begin{equation} \frac{\partial T}{\partial \vec{n}} = \frac{\partial L}{\partial \vec{n}} = \frac{\partial G}{\partial \vec{n}} = \frac{\partial V}{\partial \vec{n}} = \frac{\partial Z}{\partial \vec{n}} = 0,\; t > 0,\;x = 0,\;1, \end{equation} (4.1)

    and

    T(x,\theta) = 80,\;L(x,\theta) = 0.1,\;G(x,\theta) = 0.1,\;V(x,\theta) = 0.01,\;Z(x,\theta) = 0.001,

    for 0\leq x\leq 1 , -\tau\leq\theta\leq0, \tau = \max\{\tau_1, \tau_2\} .

    In model (1.1), we choose n(T(t)) = s-dT(t)+rT(t)\Big(1-\frac{T(t)}{K}\Big) , \pi_1(T(t), V(t)) = \frac{\beta_1 T(t)V(t)}{(1+\eta_1T(t))(1+\eta_2 V(t))} , \pi_2(T(t), G(t)) = \frac{\beta_2 T(t)G(t)}{1+\alpha_1G(t)} , and h_i(\xi) = \xi . We can easily verify that (A_1) - (A_4) hold. For simulations, we take \eta_1 = 0.01, \; \eta_2 = 0.01, \; \alpha_1 = 0.01, \; d_1 = 0.1, \; d_2 = 0.1, \; d_3 = 0.1, \; d_4 = 0.1, \; d_5 = 0.1, \tau_1 = 10, \; \tau_2 = 5, and choose \beta_1 and \beta_2 as free parameters. The values of the other parameters are summarized in Table 1.

    Table 1.  List of parameters.
    Parameter Definition Value Source
    s production rate of uninfected cells 10 \mu l^{-1}day^{-1} [37]
    d death rate of uninfected cells 0.01 day^{-1} [37]
    \eta probability that the uninfected cell
    will become an infected cell 0.49 Assumed
    r logistic growth rate 0.01 day^{-1} [38]
    K carrying capacity 1000 \mu l day^{-1} [38]
    \mu death rate for latently infected cells 0.1 day^{-1} [38]
    \alpha conversion rate 0.05 Assumed
    \sigma death rate for infected cells 0.1 day^{-1} [6]
    \gamma production rate 1 cell^{-1} day^{-1} [6]
    k clearance rate for free virus 3 day^{-1} [37]
    \rho neutralizing rate of viruses 0.5 \mu l day^{-1} [7]
    \chi generation rate of B cells 1.4 \mu l day^{-1} [37]
    \delta proliferation rate of B cells 1.2 \mu l day^{-1} [6,7]
    \beta death rate for B cells 1.2 day^{-1} [6]
    a_1 death rate for latently infected cells during [t-\tau_1, t] 0.01 Assumed
    a_2 death rate for infected cells during [t-\tau_2, t] 0.01 Assumed

     | Show Table
    DownLoad: CSV

    In the following Figures 1 and 2, (a), (b), (c), (d), and (e) are denoted time-series figures of T(t) , L(t) , G(t) , V(t) and Z(t) .

    Figure 1.  Taking \beta_1 = 0.0001 and \beta_2 = 0.0001, we have R_0 = 0.8352 < 1 , and the infection-free equilibrium E_0 = (1000, 0, 0, 0, 1.1667) is globally asymptotically stable.
    Figure 2.  Taking \beta_1 = 0.01 and \beta_2 = 0.01, we have R_0 = 83.52 > 1 , and the endemic equilibrium E_1 = (18.16, 29.62,113.9, 0.9948,223) is globally asymptotically stable.

    In this paper, a diffusive and delayed viral dynamics model with two modes of transmission has been analyzed. Some assumptions about nonlinear functions for n(T) , \pi_1(T(t), V(t)) , \pi_2(T(t), G(t)) , h_1(L(t)) , h_2(G(t)) , h_3(V(t)) , and h_4(Z(t)) are made and the global stabilities of model (1.1) are proved. The contribution is to construct suitable Lyapunov functionals for the diffusive virus model considering the humoral cells, cell-to-cell transmission, two delays, and latently infected cells, and we can extend this method to more complicated models. Furthermore, the formula of the basic reproduction number R_0 is independent of the diffusion coefficient. Without considering either the virus-to-cell infection or cell-to-cell transmission, R_0 could be under-evaluated and the transmission and spread trends of diseases need to be studied.

    Based on the obtained results of this paper, we can directly propose the following questions that need further research. On the one hand, in addition to spatial diffusion, humoral response and delays should be considered, determining whether the results obtained in this paper can be extended to a spatially heterogeneous model with immune response delay, random perturbation effect, and memory effect. On the other hand, the globally asymptotic stability of some classes of multiple infection dynamics models will be a very valuable and significative subject. We leave these problems as possible future works.

    The author declares she has not used Artificial Intelligence (AI) tools in the creation of this article.

    The author would like to express the deepest gratitude to the anonymous referees for their careful reading of the manuscript, several valuable comments, and suggestions for its improvement. This work was supported by the NSFC (Nos. 12371504, 12471471), and the General Research Fund for Shanxi Basic Research Project (Nos. 202403021221218, 202403021221214, 202103021224291).

    The author declares that there are no conflicts of interest.



    [1] X. Yang, Y. M. Su, X. J. Zhuo, T. H. Gao, Global analysis for a delayed HCV model with saturation incidence and two target cells, Chaos Soliton. Fract., 166 (2023), 112950. https://doi.org/10.1016/j.chaos.2022.112950 doi: 10.1016/j.chaos.2022.112950
    [2] A. M. Elaiw, N. H. AlShamrani, Stability of an adaptive immunity pathogen dynamics model with latency and multiple delays, Math. Method. Appl. Sci., 41 (2018), 6645–6672. https://doi.org/10.1002/mma.5182 doi: 10.1002/mma.5182
    [3] S. F. Wang, D. Y. Zou, Global stability of in host viral models with humoral immunity and intracellular delays, Appl. Math. Model., 36 (2012), 1313–1322. https://doi.org/10.1016/j.apm.2011.07.086 doi: 10.1016/j.apm.2011.07.086
    [4] J. L. Wang, S. Q. Liu, The stability analysis of a general viral infection model with distributed delays and multi-staged infected progression, Commun. Nonlinear Sci., 20 (2015), 263–272. https://doi.org/10.1016/j.cnsns.2014.04.027 doi: 10.1016/j.cnsns.2014.04.027
    [5] T. Q. Zhang, X. N. Xu, X. Z. Wang, Dynamic analysis of a cytokine-enhanced viral infection model with time delays and CTL immune response, Chaos Soliton. Fract., 170 (2023), 113357. https://doi.org/10.1016/j.chaos.2023.113357 doi: 10.1016/j.chaos.2023.113357
    [6] P. Georgescu, Y.-H. Hsieh, Global stability for a virus dynamics model with nonlinear incidence of infection and removal, SIAM J. Appl. Math., 67 (2006), 337–353. https://doi.org/10.1137/060654876 doi: 10.1137/060654876
    [7] Z. H. Yuan, X. F. Zou, Global threshold dynamics in an HIV virus model with nonlinear infection rate and distributed invasion and production delays, Math. Biosci. Eng., 10 (2013), 483–498. https://doi.org/10.3934/mbe.2013.10.483 doi: 10.3934/mbe.2013.10.483
    [8] K. Hattaf, Spatiotemporal dynamics of a generalized viral infection model with distributed delays and CTL immune response, Computation, 7 (2019), 21. https://doi.org/10.3390/computation7020021 doi: 10.3390/computation7020021
    [9] K. A. Pawelek, S. Q. Liu, F. Pahlevani, L. B. Rong, A model of HIV-1 infection with two time delays: mathematical analysis and comparison with patient data, Math. Biosci., 235 (2012), 98–109. https://doi.org/10.1016/j.mbs.2011.11.002 doi: 10.1016/j.mbs.2011.11.002
    [10] P. Balasubramaniam, P. Tamilalagan, M. Prakash, Bifurcation analysis of HIV infection model with antibody and cytotoxic T-lymphocyte immune responses and Beddington-DeAngelis functional response, Math. Method. Appl. Sci., 38 (2015), 1330–1341. https://doi.org/10.1002/mma.3148 doi: 10.1002/mma.3148
    [11] J. L. Wang, M. Guo, X. N. Liu, Z. T. 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. https://doi.org/10.1016/j.amc.2016.06.032 doi: 10.1016/j.amc.2016.06.032
    [12] N. L. Komarova, D. Wodarz, Virus dynamics in the presence of synaptic transmission, Math. Biosci., 242 (2013), 161–171. https://doi.org/10.1016/j.mbs.2013.01.003 doi: 10.1016/j.mbs.2013.01.003
    [13] A. Sigal, J. T. Kim, A. B. Balazs, E. Dekel, A. Mayo, R. Milo, et al., Cell-to-cell spread of HIV permits ongoing replication despite antiretroviral therapy, Nature, 477 (2011), 95–98. https://doi.org/10.1038/nature10347 doi: 10.1038/nature10347
    [14] S. Iwami, J. S. Takeuchi, S. Nakaoka, F. Mammano, F. Clavel, H. Inaba, et al., Cell-to-cell infection by HIV contributes over half of virus infection, eLife, 4 (2015), e08150. https://doi.org/10.7554/eLife.08150 doi: 10.7554/eLife.08150
    [15] X. Wang, S. Y. Tang, X. Y. Song, L. B. Rong, Mathematical analysis of an HIV latent infection model including both virus-to-cell infection and cell-to-cell transmission, J. Biol. Dynam., 11 (2017), 455–483. https://doi.org/10.1080/17513758.2016.1242784 doi: 10.1080/17513758.2016.1242784
    [16] H. Y. Shu, Y. M. Chen, L. Wang, Impacts of the cell-free and cell-to-cell infection modes on viral dynamics, J. Dyn. Diff. Equat., 30 (2018), 1817–1836. https://doi.org/10.1007/s10884-017-9622-2 doi: 10.1007/s10884-017-9622-2
    [17] X. L. Lai, X. F. Zou, Modelling HIV-1 virus dynamics with both virus-to-cell infection and cell-to-cell transmission, SIAM J. Appl. Math., 74 (2014), 898–917. https://doi.org/10.1137/130930145 doi: 10.1137/130930145
    [18] Y. Yang, L. Zou, S. G. Ruan, Global dynamics of a delayed within-host viral infection model with both virus-to-cell and cell-to-cell transmissions, Math. Biosci., 270 (2015), 183–191. https://doi.org/10.1016/j.mbs.2015.05.001 doi: 10.1016/j.mbs.2015.05.001
    [19] F. Li, J. L. Wang, Analysis of an HIV infection model with logistic target-cell growth and cell-to-cell transmission, Chaos Soliton. Fract., 81 (2015), 136–145. https://doi.org/10.1016/j.chaos.2015.09.003 doi: 10.1016/j.chaos.2015.09.003
    [20] A. M. Elaiw, N. H. AlShamrani, Global stability of a delayed adaptive immunity viral infection with two routes of infection and multi-stages of infected cells, Commun. Nonlinear Sci., 86 (2020), 105259. https://doi.org/10.1016/j.cnsns.2020.105259 doi: 10.1016/j.cnsns.2020.105259
    [21] A. M. Elaiw, N. H. AlShamrani, Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal, Nonlinear Anal.-Real, 26 (2015), 161–190. https://doi.org/10.1016/j.nonrwa.2015.05.007 doi: 10.1016/j.nonrwa.2015.05.007
    [22] J. Z. Lin, R. Xu, X. H. 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. https://doi.org/10.1016/j.amc.2017.08.004 doi: 10.1016/j.amc.2017.08.004
    [23] K. F. Wang, W. D. Wang, S. P. Song, Dynamics of an HBV model with diffusion and delay, J. Theor. Biol., 253 (2008), 36–44. https://doi.org/10.1016/j.jtbi.2007.11.007 doi: 10.1016/j.jtbi.2007.11.007
    [24] Y. Gao, J. L. Wang, Threshold dynamics of a delayed nonlocal reaction-diffusion HIV infection model with both cell-free and cell-to-cell transmissions, J. Math. Anal. Appl., 488 (2020), 124047. https://doi.org/10.1016/j.jmaa.2020.124047 doi: 10.1016/j.jmaa.2020.124047
    [25] Y. Yang, Y. C. Xu, Global stability of a diffusive and delayed virus dynamics model with Beddington-DeAngelies function and CTL immune response, Comput. Math. Appl., 71 (2016), 922–930. https://doi.org/10.1016/j.camwa.2016.01.009 doi: 10.1016/j.camwa.2016.01.009
    [26] H. Q. Sun, J. L. Wang, Dynamics of a diffusive virus model with general incidence function, cell-to-cell transmission and time delay, Comput. Math. Appl., 77 (2019), 284–301. https://doi.org/10.1016/j.camwa.2018.09.032 doi: 10.1016/j.camwa.2018.09.032
    [27] C. C. Travis, G. F. Webb, Existence and stability for partial functional differential equations, Trans. Amer. Math. Soc., 200 (1974), 395–418. https://doi.org/10.2307/1997265 doi: 10.2307/1997265
    [28] W. E. Fitzgibbon, Semilinear functional differential equations in Banach space, J. Differ. Equations, 29 (1978), 1–14. https://doi.org/10.1016/0022-0396(78)90037-2 doi: 10.1016/0022-0396(78)90037-2
    [29] R. H. Martin, H. L. Smith, Abstract functional differential equations and reaction-diffusion systems, Trans. Amer. Math. Soc., 321 (1990), 1–44. https://doi.org/10.1090/S0002-9947-1990-0967316-X doi: 10.1090/S0002-9947-1990-0967316-X
    [30] R. H. Martin, H. L. Smith, Reaction-diffusion systems with time delays: Monotonicity, invariance, comparison and convergence, J. Reine Angew. Math., 413 (1991), 1–35. https://doi.org/10.1515/crll.1991.413.1 doi: 10.1515/crll.1991.413.1
    [31] J. H. Wu, Theory and applications of partial functional differential equations, New York: Springer, 1996. https://doi.org/10.1007/978-1-4612-4050-1
    [32] D. Henry, Geometric theory of semilinear parabolic equations, Berlin: Springer, 1981. https://doi.org/10.1007/BFb0089647
    [33] O. Diekmann, J. A. P. Heesterbeek, J. A. J. Metz, On the definition and the computation of the basic reproduction ratio R_0 in models for infectious diseases in heterogeneous populations, J. Math. Biol., 28 (1990), 365–382. https://doi.org/10.1007/BF00178324 doi: 10.1007/BF00178324
    [34] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6 doi: 10.1016/S0025-5564(02)00108-6
    [35] H. Q. Sun, J. Li, A numerical method for a diffusive virus model with general incidence function, cell-to-cell transmission and time delay, Physica A, 545 (2020), 123477. https://doi.org/10.1016/j.physa.2019.123477 doi: 10.1016/j.physa.2019.123477
    [36] O. Nave, Modification of semi-analytical method applied system of ODE, Modern Applied Science, 14 (2020), 75–81. https://doi.org/10.5539/mas.v14n6p75 doi: 10.5539/mas.v14n6p75
    [37] Y. Wang, Y. C. Zhou, F. Brauer, J. M. Heffernan, Viral dynamics model with CTL immune response incorporating antiretroviral therapy, J. Math. Biol., 67 (2013), 901–934. https://doi.org/10.1007/s00285-012-0580-3 doi: 10.1007/s00285-012-0580-3
    [38] D. Wodarz, Hepatitis C virus dynamics and pathology: the role of CTL and antibody responses, J. Gen. Virol., 84 (2003), 1743–1750. https://doi.org/10.1099/vir.0.19118-0 doi: 10.1099/vir.0.19118-0
  • Reader Comments
  • © 2025 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(183) PDF downloads(14) Cited by(0)

Figures and Tables

Figures(2)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog