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

Comparative study between the simulation and experimental results of H2 production from water vapour plasmolysis

  • Received: 17 June 2020 Accepted: 27 August 2020 Published: 03 September 2020
  • In the present study, the kinetics of hydrogen production from water vapour using dielectric barrier discharge (DBD) plasma in a cylindrical reactor were analyzed. The simulation analysis was carried out for both models with and without the dissociative attachment reaction to predict and compare the concentration of produced hydrogen gas from water vapour. The effect of water vapour input temperature ranges of 523–623 K and plasma voltage in the range of 12–14 kV were investigated. It was revealed that the hydrogen concentration increased with the input water vapour temperature and plasma voltage increased in both the simulation models. It was seen that the H2 concentration results of the simulation model with the dissociative attachment reaction (H-) were nearly same as the H2 concentrations of the water vapour plasmolysis experimental results. Moreover, it can be concluded that the dissociative attachment reaction was controlled the H2 generation from water vapour plasmolysis. It was remarkable that the conversion rates of the simulation model included the dissociative attachment reaction has more acceptable results to the experimental data compared to the simulation model deselecting the dissociative attachment reaction (H-). Also, it was seen as the main reason for the difference between simulation and experimental results.

    Citation: Mostafa El-Shafie, Shinji Kambara, Yukio Hayakawa, Fahad Rehman. Comparative study between the simulation and experimental results of H2 production from water vapour plasmolysis[J]. AIMS Energy, 2020, 8(5): 835-858. doi: 10.3934/energy.2020.5.835

    Related Papers:

    [1] 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. Mathematical Biosciences and Engineering, 2019, 16(1): 292-319. doi: 10.3934/mbe.2019015
    [2] 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. Mathematical Biosciences and Engineering, 2021, 18(1): 274-299. doi: 10.3934/mbe.2021014
    [3] 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
    [4] 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
    [5] Jinhu Xu, Yicang Zhou . Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay. Mathematical Biosciences and Engineering, 2016, 13(2): 343-367. doi: 10.3934/mbe.2015006
    [6] 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
    [7] 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
    [8] A. M. Elaiw, A. S. Shflot, A. D. Hobiny . Stability analysis of general delayed HTLV-I dynamics model with mitosis and CTL immunity. Mathematical Biosciences and Engineering, 2022, 19(12): 12693-12729. doi: 10.3934/mbe.2022593
    [9] A. M. Elaiw, N. H. AlShamrani . Analysis of an HTLV/HIV dual infection model with diffusion. Mathematical Biosciences and Engineering, 2021, 18(6): 9430-9473. doi: 10.3934/mbe.2021464
    [10] Farzad Fatehi, Yuliya N. Kyrychko, Konstantin B. Blyuss . Time-delayed model of autoimmune dynamics. Mathematical Biosciences and Engineering, 2019, 16(5): 5613-5639. doi: 10.3934/mbe.2019279
  • In the present study, the kinetics of hydrogen production from water vapour using dielectric barrier discharge (DBD) plasma in a cylindrical reactor were analyzed. The simulation analysis was carried out for both models with and without the dissociative attachment reaction to predict and compare the concentration of produced hydrogen gas from water vapour. The effect of water vapour input temperature ranges of 523–623 K and plasma voltage in the range of 12–14 kV were investigated. It was revealed that the hydrogen concentration increased with the input water vapour temperature and plasma voltage increased in both the simulation models. It was seen that the H2 concentration results of the simulation model with the dissociative attachment reaction (H-) were nearly same as the H2 concentrations of the water vapour plasmolysis experimental results. Moreover, it can be concluded that the dissociative attachment reaction was controlled the H2 generation from water vapour plasmolysis. It was remarkable that the conversion rates of the simulation model included the dissociative attachment reaction has more acceptable results to the experimental data compared to the simulation model deselecting the dissociative attachment reaction (H-). Also, it was seen as the main reason for the difference between simulation and experimental results.


    Viral dynamics is a field of applied mathematics which employs mathematical models to describe the changes over time of infected cells and the viral load. Nowak et al. [1] and Nowak and May [2] proposed the following basic viral dynamic model,

    {du(t)dt=sdu(t)βu(t)v(t),dw(t)dt=βu(t)v(t)δw(t),dv(t)dt=Nδw(t)cv(t), (1.1)

    where u(t), w(t), and v(t) are the numbers of uninfected cells, productively infected cells, and virus particles at time t, respectively. See the references for the biological meanings of the parameters. This basic model has been modified to study different viral infections, which include hepatitis C virus (HCV) [3,4], human immunodeficiency virus (HIV) [5,6,7,8], human T-cell leukemia virus (HTLV) [9,10,11], and so on.

    During the process of viral infection, specific immune response plays an important role. Specific immune response includes cell-mediated immunity (which depends on cytotoxic T lymphocytes response (CTLs)) and humoral immunity (which depends on antibody response). Since the work of Nowak and Bangham [12], much has been done on mathematical models on immune response against infected cells [13,14,15,16].

    Nowadays, time delays have been taken into account in order to better understand viral dynamics. Usually, distributed time delays [17,18,19] and discrete time delays [20,21,22] have been incorporated into viral dynamic models. In particular, based on (1.1), Zhu and Zou [20] proposed the following viral dynamic model with time delay and CTL immune response,

    {du(t)dt=sdu(t)βu(t)v(t),dw(t)dt=βemˆτu(tˆτ)v(tˆτ)δw(t)pw(t)z(t),dv(t)dt=Nδw(t)cv(t),dz(t)dt=qw(t)z(t)bz(t), (1.2)

    where z(t) denotes the density of immune effectors at time t. Here the delay ˆτ represents the time from a virus entering a target cell to the production of new free virus particles. We refer to [20] for the meanings of the other parameters.

    Note that both models (1.1) and (1.2) and most existing ones are described by ordinary differential equations. The cells and free virus particles are assumed to be uniform in location. In other words, the effect of spatial heterogeneity is ignored. For example, the lymphoid tissues are among the primary sites of HIV infection and replication. The lymphoid tissues consist of many lymph nodes with different sizes. The different tissue architecture and composition and biophysical parameters can influence the spread and replication of the virus [23]. To understand the viral pathogenesis better, it is necessary to consider the spatial aspects of the tissues. In [24], Wang and Wang proposed the following mathematical model of HBV infection with spatial dependence,

    {u(x,t)t=sdu(x,t)βu(x,t)v(x,t),w(x,t)t=βu(x,t)v(x,t)δw(x,t),v(x,t)t=DΔv(x,t)+Nδw(x,t)cv(x,t), (1.3)

    where u(x,t), w(x,t), and v(x,t) are the densities of uninfected cells, productively infected cells, and free virus particles at spatial position x and time t, respectively. D is the diffusion coefficient and Δ is the Laplace operator. Using the geometric singular perturbation method, they studied the existence of traveling waves. Since then a lot of works have followed in this direction (see, for example, [25,26,27,28,29]).

    In (1.1), (1.2), and (1.3), only the cell-free transmission (newly released free virus particles infect uninfected cells [2]) is considered. Recent experimental studies [30,31] prove that a healthy cell can be infected when it comes with close contact of an infected cell (cell-to-cell transmission [32,33]). Sigal et al. [34] found that the cell-to-cell spread of HIV can still permit ongoing replication even with an antiretroviral therapy. Consequently, viral dynamic models incorporating both transmission modes have been formulated and studied (to name a few, see [35,36,37,38,39]). We should mention that the incidences in these works are bilinear. Incidence is the number of new infections per unit of time. It depends on the infectivity of viruses and behavior of cells. Thus it is reasonable to be nonlinear in general. For example, the saturated incidence rate βuv1+αv is used in [40] and the Beddington-DeAngelis incidence function is used in [41]. In a recent work, Sun and Wang [42] also used a general incidence f(u,v) in a diffusive viral dynamic model.

    Based on the above discussion, in this paper, we propose and study the following delayed diffusive viral dynamic model with cell-mediated immunity, cell-to-cell transmission, and general incidences,

    {u(x,t)t=sdu(x,t)f(u(x,t),v(x,t))g(u(x,t),w(x,t)),w(x,t)t=emτf(u(x,tτ),v(x,tτ))+emτg(u(x,tτ),w(x,tτ))δw(x,t)pw(x,t)z(x,t),v(x,t)t=D1Δv(x,t)+Nδw(x,t)cv(x,t), xΩ, t>0,z(x,t)t=D2Δz(x,t)+qw(x,t)z(x,t)bz(x,t), xΩ, t>0, (1.4)

    where z(x,t) denotes the densities of immune effectors at spatial position x and time t. Ω is a general open bounded domain in Rn with smooth boundary Ω. We consider model (1.4) with the homogeneous Neumann boundary conditions

    vn=0, zn=0, xΩ, t>0, (1.5)

    where n denotes the outward normal derivative on Ω. We also assume the initial conditions

    u(x,θ)=ϕ1(x,θ)0,w(x,θ)=ϕ2(x,θ)0,v(x,θ)=ϕ3(x,θ)0,z(x,θ)=ϕ4(x,θ)0,(x,θ)¯Ω×[τ,0], (1.6)

    where ϕi's (i=1, 2, 3, 4) are bounded and uniformly continuous functions on ¯Ω×[τ,0].

    In (1.4), intracellular delays for both transmission modes are assumed to be the same. In general, the intracellular delay in the cell-to-cell transmission is less than that in the cell-free infection [30,39]. However, the difference is not large enough. As for some existing studies (for example, [35,38]), for simplicity of presentation, we make the above assumption on the delays.

    In (1.4), the incidences due to the cell-free transmission and the cell-to-cell transmission are given by the nonlinear functions f(u,v) and g(u,w), respectively. As in [26], we always make the following assumption on them in the sequel.

    (A1) The nonlinear incidence functions f and g satisfy the following properties.

    (i) f(u,v)0 and g(u,w)0 for u0, v0, and w0, and the equalities hold if and only if uv=0 and uw=0;

    (ii) There exists η1>0 and η2>0 such that f(u,v)η1u and g(u,w)η2u for u0, v0, and w0;

    (iii) f(u,v)u and g(u,w)u are continuous with f(u,v)u>0 and g(u,w)u>0 for u0, v>0, and w>0;

    (iv) f(u,v)v and g(u,w)w are continuous with f(u,v)v0 and g(u,w)w0 for u0, v0, and w0;

    (v) vf(u,v)vf(u,v)0 and wg(u,w)wg(u,w)0 for u0, v0, and w0.

    Note that, by Assumption (A1), for any u>0,

    (f(u,v)v)v=f(u,v)vvf(u,v)v2<0,

    which implies that f(u,v)v is decreasing on (0,). In particular, for any u>0 and v>0,

    f(u,v)vlimv0+f(u,v)v=f(u,0)v.

    An analog also holds for g. Thus we have

    f(u,v)f(u,0)vv and g(u,w)g(u,0)ww for u0, v0, and w0. (1.7)

    The rest of the paper is organized as follows. In section 2, we consider the existence, uniqueness, positivity, and boundedness of solutions to system (1.4)–(1.6). Then we study the existence of homogeneous steady states in section 3, which depend on the basic reproduction number of infection and the basic reproduction number of immunity. The main part is section 4, where we discuss the local and global dynamics of system (1.4)–(1.6) by analyzing the characteristic equations and constructing suitable Lyapunov functionals. These results are supported with numerical simulations in section 5. The paper ends with a brief conclusion.

    Let X:=C(¯Ω,R4) be the Banach space equipped with the supremum norm X. For τ0, define C=C([τ,0],X), which is a Banach space equipped with the norm ϕ=maxθ[τ,0]ϕ(θ)X. If σ>0 and U:[τ,σ)X, then for t[0,σ), UtC is defined by Ut(θ)=U(t+θ) for θ[τ,0]. Denote X+=C(¯Ω,R4+) and C+=C([τ,0],X+). Then both (X,X+) and (C,C+) are strongly ordered spaces. According to Corollary 4 in [43], we have the following result on the well-posedness. The arguments are standard and hence are omitted here. Interested readers can refer to, for example, a recent paper by Gao and Wang [44].

    Theorem 2.1. For each ϕ=(ϕ1,ϕ2,ϕ3,ϕ4)C+, system (1.4)–(1.6) has a unique solution U(,t,ϕ)=(u(,t,ϕ),w(,t,ϕ),v(,t,ϕ),z(,t,ϕ)) on [0,) with U0(,ϕ)=ϕ. Moreover, Ut(,ϕ)C+ for t0 and U(,t,ϕ) is a classical solution.

    Let Φ(t):C+C+ be the solution semiflow associated with (1.4)–(1.6), that is, Φ(t,ϕ)=Ut(,ϕ), where U(,t,ϕ) is the solution of (1.4)–(1.6) with the initial condition ϕC+.

    The following result gives some properties of solutions.

    Lemma 2.2. For ϕC+, the following statements hold for the solution U(,t,ϕ) of (1.4)–(1.6).

    (i) lim suptu(x,t,ϕ)sd, lim suptw(x,t,ϕ)emτ(η1+η2)sdδ, lim suptv(x,t,ϕ)Nemτ(η1+η2)sdc, and lim suptz(x,t,ϕ)emτ(η1+η2)sdmin{δ,b} uniformly for all xΩ.

    (ii) u(,t,ϕ)>0 for t>0 and lim inftu(x,t,ϕ)sd+η1+η2 uniformly for all xΩ.

    (iii) If w(,t0,ϕ) or v(\cdot, t_0, \phi)\not\equiv0 for some t_0\geq0 , then w(x, t, \phi) > 0 and v(x, t, \phi) > 0 for all x\in\Omega and t > t_0+\tau .

    {\rm(iv)} If z(\cdot, t_0, \phi)\not\equiv0 for some t_0\geq0 , then z(x, t, \phi) > 0 for all x\in\Omega and t > t_0 .

    Proof. For simplicity of notation, in the proof here we omit \phi from the expressions of the solution.

    (ⅰ) First, we have

    \frac{\partial u(x, t)}{\partial t}\le s-du(x, t),

    which implies that \limsup\limits_{t\to\infty}u(x, t)\le \frac{s}{d} uniformly for all x\in \Omega . Next, by Assumption (A1) and the second equation in (1.4), we have

    \frac{\partial w(x, t)}{\partial t}\le e^{-m\tau}(\eta_1+\eta_2) u(x, t)-\delta w(x, t).

    Then \limsup\limits_{t\to\infty}w(x, t)\le \frac{e^{-m\tau}(\eta_1+\eta_2)s}{d\delta} uniformly for x\in \Omega follows easily from this and \limsup\limits_{t\to\infty}u(x, t)\le \frac{s}{d} uniformly for x\in \Omega . Similarly, adding the second and fourth equations of (1.4) yields

    \frac {\partial(w(x, t)+z(x, t))}{\partial t}\le e^{-m\tau}(\eta_1+\eta_2)u(x, t)-\min \{\delta, b\}\big(w(x, t)+z(x, t)\big).

    It follows that \limsup\limits_{t\to\infty}(w(x, t)+z(x, t))\le \frac {e^{-m\tau}(\eta_1+\eta_2)s}{d\min\{\delta, b\}} uniformly for x\in\Omega and hence \limsup\limits_{t\to\infty}z(x, t)\le \frac {e^{-m\tau}(\eta_1+\eta_2)s}{d\min\{\delta, b\}} uniformly for x\in\Omega . Now, \limsup\limits_{t\to\infty}w(x, t)\le \frac{e^{-m\tau}(\eta_1+\eta_2)s}{d\delta} uniformly for x\in \Omega together with the third equation of (1.4) (Lemma 1 in [45]), and comparison theorem, gives \limsup\limits_{t\to\infty}v(x, t)\le \frac{Ne^{-m\tau}(\eta_1+\eta_2)s}{dc} uniformly for x\in\Omega .

    (ⅱ) Noting that \frac{\partial u(x, t)}{\partial t} \geq s-(d+\eta_1+\eta_2)u(x, t) , one can easily get

    u(x, t)\ge e^{-(d+\eta_1+\eta_2)t}u(x, 0)+\frac{s}{d+\eta_1+\eta_2}-\frac{e^{-(d+\eta_1+\eta_2)t}s}{d+\eta_1+\eta_2}

    for x\in \Omega and t\ge 0 . Then (ⅱ) follows immediately.

    (ⅲ) Note that the operator D_1\Delta -c\mathrm {Id} generates a positive semigroup on C(\overline{\Omega}, \mathbb{R}) , where \mathrm {Id} is the identity operator. Thus if w(\cdot, t_0)\not\equiv0 , then from the third equation of (1.4), we see that v(\cdot, t)\not\equiv 0 for t > t_0 . Without loss of generality, we assume that v(\cdot, t_0)\not\equiv 0 . We first show that v(\cdot, t) > 0 for t > t_0 . By Theorem 2.1, v(x, t) satisfies

    \begin{cases} \frac{\partial v(x, t)}{\partial t} \geq D_1\Delta v(x, t)-cv(x, t), \quad t \gt t_0, \ x\in\Omega, \\ \frac{\partial v(x, t)}{\partial \vec{n}} = 0, \quad t \gt t_0, \ x\in\partial\Omega. \end{cases}

    Let \bar{v}(x, t) be the solution of

    \begin{cases} \frac{\partial \bar{v}(x, t)}{\partial t} = D_1\Delta \bar{v}(x, t)-c\bar{v}(x, t), \quad t \gt t_0, \ x\in\Omega, \\ \frac{\partial \bar{v}(x, t)}{\partial \vec{n}} = 0, \quad t \gt t_0, \ x\in\partial\Omega, \\ \bar{v}(x, t_0) = v(x, t_0), \ x\in \overline{\Omega}. \end{cases}

    Then \bar{v}(x, t) > 0 for x\in\Omega and t > t_0 . In fact, suppose, by contradiction, there exist x_0\in\Omega and \hat{t} > t_0 such that \bar{v}(x_0, \hat{t}) = 0 . Then, according to the strong maximum principle [46], \bar{v}(x, t)\equiv0 for each t\geq t_0 , contradicting with \bar{v}(\cdot, t_0)\not\equiv0 . Applying the comparison theorem, we know that v(x, t)\ge\bar{v}(x, t) > 0 for t > t_0 and x\in \Omega . We now prove that w(x, t) > 0 for x\in \Omega and t > t_0+\tau . Otherwise, there exist \bar{x}\in\Omega and \bar{t} > t_0+\tau such that w(\bar{x}, \bar{t}) = 0 . As w(x, t)\ge 0 , we have \frac{\partial w(\bar{x}, \bar{t})}{\partial t} = 0 . This is impossible since

    \frac{\partial w(\bar{x}, \bar{t})}{\partial t} = e^{-m\tau}f(u(\bar{x}, \bar{t}-\tau), v(\bar{x}, \bar{t}-\tau)) +e^{-m\tau}g(u(\bar{x}, \bar{t}-\tau), w(\bar{x}, \bar{t}-\tau)) \gt 0

    by Assumption (A1) (ⅱ) due to u(\bar{x}, \bar{t}-\tau) > 0 and v(\bar{x}, \bar{t}-\tau) > 0 . This proves statement (ⅲ).

    (ⅳ) The proof is similar to that of (ⅲ) on v(x, t) > 0 for x\in \Omega and t > t_0 and hence we omit it here. This completes the proof.

    Lemma 2.2 tells us that \Phi is point dissipative. Then it follows from Theorem 2.1.8 in [47] that \Phi(t) is compact for all t > \tau . This, together with Theorem 3.4.8 in [48], gives the following result.

    Theorem 2.3. The semiflow \Phi has a global compact attractor \mathcal{A} in \mathcal{C}^+ . Moreover, u(x, t, \phi)\le \frac{s}{d} for all x\in \overline{\Omega} , t\ge 0 , and \phi\in \mathcal{A} .

    System (1.4) with (1.5) always has a unique infection-free steady state P_0 = (u_0, 0, 0, 0) , where u_0 = s/d . Applying the result of Wang and Zhao (Theorem 3.4 in [49]), we can obtain the expression of the basic reproduction number of infection, R_0 , which is given by

    \begin{eqnarray*} R_0 = \frac{Ne^{-m\tau}}{c}\cdot\frac{\partial f(\frac{s}{d}, 0)}{\partial v}+\frac{e^{-m\tau}}{\delta}\cdot\frac{\partial g(\frac{s}{d}, 0)}{\partial w}. \end{eqnarray*}

    Denote

    R_{01} = \frac{Ne^{-m\tau}}{c}\cdot\frac{\partial f(\frac{s}{d}, 0)}{\partial v}\quad \text{and} \quad R_{02} = \dfrac{e^{-m\tau}}{\delta}\cdot\dfrac{\partial g(\frac{s}{d}, 0)}{\partial w}.

    Then R_{01} is the number of secondly infected cells through the cell-free transmission and it is referred to as the basic reproduction number from the cell-free transmission; while R_{02} is the number of secondly infected cells through the cell-to-cell transmission and it is referred to as the basic reproduction number from the cell-to-cell transmission [38].

    In the following, we discuss the existence of homogeneous steady states for (1.4) (the stability results in section 4 indicate that they are the only possible steady states). Clearly, a homogeneous steady state P = (u, w, v, z) satisfies

    \begin{align} &s-du-f(u, v)-g(u, w) = 0, \end{align} (3.1a)
    \begin{align} &f(u, v)e^{-m\tau}+g(u, w)e^{-m\tau}-\delta w-pwz = 0, \end{align} (3.1b)
    \begin{align} &N\delta w-cv = 0, \end{align} (3.1c)
    \begin{align} &qwz-bz = 0. \end{align} (3.1d)

    It follows from (3.1d) that z = 0 (which corresponds to the immunity-free infected steady states) or w = \frac{b}{q} (which, when z\neq 0 , corresponds to the infected-immune steady states).

    We firstly consider the case where z = 0 . It follows from (3.1c) that v = \frac{N\delta w}{c} . Multiplying both sides of (3.1b) by e^{m\tau} and then adding up the resultant and (3.1a) to get u = \frac{s-\delta w e^{m\tau}}{d} . It is necessary that w\in(0, \frac{s}{\delta e^{m\tau}}) . Substituting u = \frac{s-\delta w e^{m\tau}}{d} and v = \frac{N\delta w}{c} into (3.1a), we see that w is a positive zero of H_1 , where

    \begin{eqnarray} H_1(w) = f\left(\dfrac{s-\delta w e^{m\tau}}{d}, \dfrac{N\delta w}{c}\right)+g\left(\dfrac{s-\delta w e^{m\tau}}{d}, w\right)-\delta we^{m\tau}. \end{eqnarray} (3.2)

    According to Assumption (A1), we have H_1(0) = 0 , H_1 (\frac{s}{\delta e^{m\tau}})-s < 0 , and

    H_1^{\prime}(0) = \delta e^{m\tau}\left(\frac{Ne^{-m\tau}}{c}\cdot\frac{\partial f(\frac{s}{d}, 0)}{\partial v} +\frac{e^{-m\tau}}{\delta}\cdot\frac{\partial g(\frac{s}{d}, 0)}{\partial w}\right)-\delta e^{m\tau} = \delta e^{m\tau}(R_0-1).

    If R_0 > 1 , then H_1'(0) > 0 . This, together with H_1(0) = 0 , implies that H_1(w) is positive for all sufficiently small w > 0 . By the Intermediate Value Theorem, H_1 has at least one zero in (0, \frac{s}{\delta e^{m\tau}}) and hence (1.4) has at least one immunity-free infected steady state. In fact, there is only one such steady state by the claim that H_1'(w_1) < 0 for any immunity-free infected steady state, which is proved as follows. Note that \delta e^{m\tau} = \frac{f(u_1, v_1)}{w_1}+\frac{g(u_1, w_1)}{w_1} and w_1 = \frac{cv_1}{N\delta} . By Assumption (A1),

    \begin{eqnarray*} H_1^{\prime}(w_1)& = & -\frac{\delta e^{m\tau}}{d}\frac{\partial f(u_1, v_1)}{\partial u}+\frac{N\delta}{c}\frac{\partial f(u_1, v_1)}{\partial v} -\frac{\delta e^{m\tau}}{d}\frac{\partial g(u_1, w_1)}{\partial u}\\ && +\ \frac{\partial g(u_1, w_1)}{\partial w}-\frac{N\delta}{cv_1}f(u_1, v_1)-\frac{1}{w_1}g(u_1, w_1)\\ & = & -\frac{\delta e^{m\tau}}{d}\frac{\partial f(u_1, v_1)}{\partial u}-\frac{\delta e^{m\tau}}{d}\frac{\partial g(u_1, w_1)}{\partial u} +\ \frac{N\delta}{cv_1}\left(v_1\frac{\partial f(u_1, v_1)}{\partial v}-f(u_1, v_1)\right)\\ && +\ \frac{1}{w_1}\left(w_1\frac{\partial g(u_1, w_1)}{\partial w}-g(u_1, w_1)\right)\\ & \lt & 0. \end{eqnarray*}

    This proves the claim. Next, we assume that R_0 < 1 . Then H_1^{\prime}(0) = \delta e^{m\tau}(R_0-1) < 0 , which combined with H_1(0) = 0 implies that H_1(w) < 0 for w > 0 sufficiently small. Using the above claim, we can easily see that there is no immunity-free infected steady state when R_0 < 1 . Moreover, H_1(w) < 0 for w\in (0, \frac{s}{\delta e^{m\tau}}] . Finally, we assume that R_0 = 1 . We use contradictive arguments to show that there is no immunity-free infected steady state in this case. Otherwise, assume that H_1(w) has a positive zero say w^* . Then from the above claim H_1(w) > 0 for w < w^* and closely enough to w^* . Note that H_1(w) depends continuously on the parameters and H_1(w) < 0 for w\in (0, \frac{s}{\delta e^{m\tau}}] when R_0 < 1 . Fix w\in (0, w^*) . Choose a sequence of parameters such that the basic reproduction number R_0 < 1 and tends to 1 . Then H_1(w) tends to H_1(w^*) > 0 , a contradiction to the fact that the limit is less than or equal to 0. This proves that there is no immunity-free infected steady state when R_0 = 1 .

    Now we study the case where w = \frac{b}{q} . This, together with (3.1c), yields v = \frac{N\delta b}{cq} . As before, add up (3.1a) and (3.1b) multiplied by e^{m\tau} to get z = \frac{s-du-\delta e^{m\tau}w}{pe^{m\tau}w} , which necessarily requires u\in(0, \frac{s}{d}-\frac{\delta b}{dq}e^{m\tau}) . Substituting w = \frac{b}{q} and z = \frac{s-du-\delta e^{m\tau}w}{pe^{m\tau}w} into (3.1a), we see that u is a positive zero of H_2 , where

    \begin{eqnarray} H_2(u) = f\left(u, \frac{N\delta b}{cq}\right)+g\left(u, \frac{b}{q}\right)-s+du. \end{eqnarray} (3.3)

    With Assumption (A1), we have H_2(0) = -s < 0 and

    \begin{eqnarray*} H_2^{\prime}(u) = \frac{\partial f(u, \frac{N\delta b}{cq})}{\partial u}+\frac{\partial g(u, \frac{b}{q})}{\partial u}+d \gt 0. \end{eqnarray*}

    Therefore, in order for model (1.4) to have an infected-immune steady state (if exists there is a unique one), it is necessary and sufficient that H_2(\frac{s}{d}-\frac{\delta b}{d q}e^{m\tau}) = H_1(\frac{b}{q}) > 0. Recall that when R_0 \leq 1, \ H_1(w) < 0 for w \in (0, \frac{s}{\delta e^{m\tau}}] ; while when R_0 > 1, \ H_1(w) > 0 for w\in (0, w_1) and H_1(w) < 0 for w\in(w_1, \frac{s}{\delta e^{m\tau}}) . It follows that H_1(\frac{b}{q}) > 0 if and only if R_0 > 1 and \frac{b}{q} < w_1 . Denote,

    R_1 = \dfrac{qw_1}{b}.

    As q is the average number of immune effectors produced from contacting with a productively infected cell and \frac{1}{b} is the average life of an immune effector, it follows that R_1 is the total number of immune effectors produced at the immunity-free infected steady state. Thus R_1 is called the basic reproduction number of immunity.

    Summarizing the above discussion, we have obtained the following result on the existence of homogeneous steady states.

    Theorem 3.1. For model (1.4) with (1.5), the following statements on the existence of homogeneous steady states are true.

    {\rm(i)} If R_0\le1 , then there is only the infection-free steady state P_0 .

    {\rm (ii)} If R_1\leq1 < R_0 , then besides P_0 , there is also a unique immunity-free infected steady state P_1 = (u_1, w_1, v_1, 0) , where w_1 is the only positive zero of H_1 defined by (3.2), u_1 = \frac{s-\delta w_1 e^{m\tau}}{d} and v_1 = \frac{N\delta w_1}{c} .

    {\rm (iii)} If R_1 > 1 (it is necessary that R_0 > 1 ), then in addition to P_0 and P_1 , there is also a unique infected-immune steady state P_2 = (u_2, w_2, v_2, z_2) , where u_2 is the only positive zero of H_2 defined by (3.3), w_2 = \frac{b}{q} , v_2 = \frac{N\delta b}{cq} , and z_2 = \frac{s-du_2-\delta w_2 e^{m\tau}}{pw_2e^{m\tau}} .

    In the main part of this paper, we establish the stability of each steady state obtained in Theorem 3.1.

    Let P^* = (u^*, w^*, v^*, z^*) be an arbitrary homogeneous steady state. The linearization of (1.4) at P^* is

    \begin{equation} \frac{\partial Q}{\partial t} = \mathbb{L}\Delta Q+AQ+BQ_{\tau}, \end{equation} (4.1)

    where

    \begin{eqnarray*} \mathbb{L} & = &\mathrm{diag}(0, 0, D_1, D_2), \\ Q & = & (u(x, t), w(x, t), v(x, t), z(x, t)), \\ Q_{\tau} & = & (u_{\tau}, w_{\tau}, v_{\tau}, z_{\tau}) = \big(u(x, t-\tau), w(x, t-\tau), v(x, t-\tau), z(x, t-\tau)\big), \\ A & = & \begin{pmatrix} -d-\frac{\partial f(u^*, v^*)}{\partial u}-\frac{\partial g(u^*, w^*)}{\partial u} & -\frac{\partial g(u^*, w^*)}{\partial w} & -\frac{\partial f(u^*, v^*)}{\partial v} & 0\\ 0 & -\delta-pz^* & 0 & -pw^*\\ 0& N\delta & -c & 0 \\ 0 & qz^* & 0 & qw^*-b \end{pmatrix}, \\ B & = & \begin{pmatrix} 0 & 0 & 0 & 0\\ \left(\frac{\partial f(u^*, v^*)}{\partial u}+\frac{\partial g(u^*, w^*)}{\partial u}\right)e^{-m\tau} & \frac{\partial g(u^*, w^*)}{\partial w} e^{-m\tau} & \frac{\partial f(u^*, v^*)}{\partial v}e^{-m\tau} & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \end{pmatrix}. \end{eqnarray*}

    Denote 0 = \mu_0 < \mu_1 < \mu_2 < \cdots < \mu_n < \cdots to be all the eigenvalues of the operator -\Delta on \Omega with the homogeneous Neumann boundary condition. Then P^* is locally asymptotically stable if, for any i\in \mathbb{N} = \{0, 1, 2, \dots\} , every solution of the characteristic equation

    \begin{eqnarray} |\lambda E+\mu_i\mathbb{L}-A-Be^{-\lambda\tau}| = 0 \end{eqnarray} (4.2)

    has a negative real part and P^* is unstable if there exists i_0\in \mathbb{N} such that (4.2) has a solution with a positive real part.

    We first study the local stability of P_0 .

    Proposition 4.1. The infection-free steady state P_0 of (1.4) is locally asymptotically stable if R_0 < 1 and unstable if R_0 > 1 .

    Proof. By (4.2), the characteristic equation at P_0 is

    (\lambda+b+\mu_i D_2)(\lambda+d)\left[\left(\lambda+c+\mu_iD_1\right)\left(\lambda+\delta-\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-(\lambda+m)\tau}\right)- N\delta\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-(\lambda+m)\tau}\right] = 0.

    Obviously, the stability of P_0 is determined by

    \begin{equation} (\lambda+c+\mu_iD_1)\left(\lambda+\delta-\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau}e^{-\lambda\tau}\right)-N\delta\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}e^{-\lambda\tau} = 0. \end{equation} (4.3)

    Firstly, suppose that R_0 < 1 . We claim that all roots of (4.3) have negative real parts. Otherwise, there exists i_0\in \mathbb{N} such that (4.3) has a root \lambda_0 with \mathrm{Re}(\lambda_0)\ge 0 . Then

    1 = \frac{1}{\lambda_0+\delta}\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau}e^{-\lambda_0\tau} +\frac{N\delta}{(\lambda_0+\delta)(\lambda_0+c+\mu_{i_0}D_1)}\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}e^{-\lambda_0\tau}.

    It follows that

    \begin{eqnarray*} 1 & = & \left| \frac{1}{\lambda_0+\delta}\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau}e^{-\lambda_0\tau} +\frac{N\delta}{(\lambda_0+\delta)(\lambda_0+c+\mu_{i_0}D_1)}\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}e^{-\lambda_0\tau}\right|\\ &\leq& \left|\frac{1}{\lambda_0+\delta}\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau}e^{-\lambda_0\tau}\right| +\left|\frac{N\delta}{(\lambda_0+\delta)(\lambda_0+c+\mu_{i_0}D_1)}\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}e^{-\lambda_0\tau}\right|\\ &\leq& \frac{1}{\delta}\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau} +\frac{N}{c+\mu_{i_0}D_1}\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}\\ & \le & \frac{1}{\delta}\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau} +\frac{N}{c}\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}\\ & = & R_0, \end{eqnarray*}

    a contradiction to R_0 < 1 . This proves the claim and hence P_0 is locally asymptotically stable when R_0 < 1 .

    Secondly, assume R_0 > 1 . For i\in \mathbb{N} , denote

    F(\lambda, i) = (\lambda+c+\mu_iD_1)\left(\lambda+\delta-\frac{\partial g(\frac{s}{d}, 0)}{\partial w} e^{-m\tau}e^{-\lambda\tau}\right)-N\delta\frac{\partial f(\frac{s}{d}, 0)}{\partial v} e^{-m\tau}e^{-\lambda\tau}.

    Recall that \mu_0 = 0 . We have

    F(0, 0) = c\delta(1-R_0) \lt 0

    and

    F(\lambda, 0) = \lambda^2+(c+\delta)\lambda+c\delta-\left(c\delta R_0+\frac{\partial g(\frac{s}{d}, 0)}{\partial w}e^{m\tau}\right)e^{-\lambda\tau}\to \infty\ \text{as}\ \lambda\to\infty.

    By the Intermediate Value Theorem, F(\lambda, 0) has a positive zero and hence (4.3) has at least one positive zero for i = 0 . This means that P_0 is unstable when R_0 > 1 .

    In fact, P_0 is globally stable if it is locally stable.

    Theorem 4.2. If R_0\leq1 , then the infection-free steady state P_0 of (1.4) is globally attractive. In particular, P_0 is globally asymptotically stable when R_0 < 1 .

    Proof. It suffices to show that P_0 is globally attractive in \mathcal{A} . For this purpose, we consider the Lyapunov functional

    \begin{eqnarray*} W(t) & = & \int_{\Omega}\left(e^{m\tau}w(x, t)+\frac{1}{N}e^{m\tau}v(x, t)+\frac{p}{q}e^{m\tau}z(x, t)\right)\mathrm{d}x\\ && + \int_{\Omega}\left(\int_{t-\tau}^{t}f(u(x, \theta), v(x, \theta))\mathrm{d}\theta +\int_{t-\tau}^{\infty}g(u(x, \theta), w(x, \theta))\text{d}\theta\right)\mathrm{d}x. \end{eqnarray*}

    Calculating the time derivative of W(t) along solutions of model (1.4), we have

    \begin{eqnarray*} \frac{\mathrm{d} W(t)}{\mathrm{d} t} = \int_{\Omega}\left(f(u(x, t), v(x, t))-\frac{c}{N}e^{m\tau}v(x, t)-\frac{pb}{q}e^{m\tau}z(x, t)\right)\mathrm {d}x +\frac{1}{N}e^{m\tau}\int_{\Omega}D_1\Delta v(x, t)\mathrm {d}x. \end{eqnarray*}

    It follows from the homogeneous Neumann boundary condition (1.5) and the Divergence Theorem that

    \int_{\Omega}\Delta v(x, t)\mathrm {d}x = \int_{\partial\Omega}\frac{\partial v(x, t)}{\partial \vec{n}}\mathrm {d}x = 0.

    Moreover, by Theorem 2.3, u(x, t)\le \frac{s}{d} for x\in \Omega and t\ge 0 . With the help of (1.7), for v(x, t) > 0 , we have

    \begin{eqnarray*} f(u(x, t), v(x, t))-\frac{c}{N}e^{m\tau}v(x, t) & = & \frac{c}{N}e^{m\tau}v(x, t)\left(\frac{N}{c}e^{-m\tau}\frac{f(u(x, t), v(x, t))}{v(x, t)}-1\right)\\ &\leq& \frac{c}{N}e^{m\tau}v(x, t)\left(\frac{N}{c}e^{-m\tau}\frac{f(\frac{s}{d}, v(x, t))}{v(x, t)}-1\right)\\ &\le & \frac{c}{N}e^{m\tau}v(x, t)\left(\frac{N}{c}e^{-m\tau}\frac {\partial f(\frac{s}{d}, 0)}{\partial v}-1\right)\\ &\leq& \frac{c}{N}e^{m\tau}v(x, t)\left(R_0-1\right). \end{eqnarray*}

    The above inequality holds automatically for v(x, t) = 0 and also observe that the inequality is strict for u(x, t) < \frac{s}{d} and v(x, t) > 0 . Therefore,

    \begin{eqnarray*} \frac{\mathrm{d} W(t)}{\mathrm{d} t} &\leq& \int_{\Omega}\left(\frac{c}{N}e^{m\tau}v(x, t)\left(R_0-1\right) -\frac{pb}{q}e^{m\tau}z(x, t)\right)\mathrm {d}x\le 0. \end{eqnarray*}

    Moreover, \frac{\mathrm{d} W(t)}{\mathrm{d} t} = 0 if and only if v(x, t) = 0 and z(x, t) = 0 . In fact, if v(x_0, t_0)\neq 0 , then there exists a neighborhood N_{(x_0, t_0)} of (x_0, t_0) such that v(x, t)\neq0 for (x, t)\in N_{(x_0, t_0)} . Then by the observation, u(x, t) = \frac{s}{d} for (x, t)\in N_{(x_0, t_0)} . This, together with the first equation of (1.4) and Assumption (A1), implies that v(x, t) = 0 for (x, t)\in N_{(x_0, t_0)} , a contradiction. Then it is easy to see that the largest invariant subset of \frac{\mathrm{d} W(t)}{\mathrm{d} t} = 0 is \{P_0\} . By LaSalle's Invariance Principle (see Theorem 5.3.1 in [50] or Theorem 3.4.7 in [51]), the infection-free steady state P_0 is globally attractive. In particular, this together with Proposition 4.1, tells us that P_0 is globally asymptotically stable when R_0 < 1 .

    Next we consider the stability of the immunity-free infected steady state P_1 . For convenience of notations, denote

    \begin{array}{lll} f_1 = f(u_1, v_1), \quad & g_1 = g(u_1, w_1), \quad & f_{1u} = \frac{\partial f(u_1, v_1)}{\partial u}, \\ f_{1v} = \frac{\partial f(u_1, v_1)}{\partial v} , & g_{1u} = \frac{\partial g(u_1, w_1)}{\partial u}, &g_{1w} = \frac{\partial g(u_1, w_1)}{\partial w}. \end{array}

    Theorem 4.3. Suppose R_0 > 1 . Then the immunity-free infected steady state P_1 of (1.4) is locally asymptotically stable if R_1 < 1 and unstable if R_1 > 1 .

    Proof. From (4.2), we know that the characteristic equation at P_1 is given by

    (\lambda-qw_1+b+\mu_iD_2)\rho_i(\lambda) = 0, \ i\in \mathbb{N},

    where

    \rho_i(\lambda) = (\lambda+d+f_{1u}+g_{1u})(\lambda+\delta)(\lambda+c+\mu_iD_1)-\big[(\lambda+d)(\lambda+c+\mu_iD_1)g_{1w} +(\lambda+d)N\delta f_{1v}\big]e^{-(\lambda+m)\tau}.

    Clearly, the eigenvalue \lambda = b(R_1-1)-\mu_iD_2 < 0 for i\in \mathbb{N} when R_1 < 1 but when R_1 > 1 , with i = 0 , we have a positive eigenvalue \lambda = b(R_1-1) . Thus P_1 is unstable if R_1 > 1 . Now, we assume that R_1 < 1 . Then the stability of P_1 is determined by the roots of \rho_i(\lambda) = 0 , which is equivalent to

    \begin{eqnarray} 1 = \frac{\lambda+d}{\lambda+d+f_{1u}+g_{1u}}\left(\frac{g_{1w}}{\lambda+\delta}e^{-(\lambda+m)\tau} +\frac{N\delta f_{1v}}{(\lambda+\delta)(\lambda+c+\mu_iD_1)}e^{-(\lambda+m)\tau}\right). \end{eqnarray} (4.4)

    We claim that all solutions of (4.4) have negative real parts. Otherwise, suppose that there exists i_{1}\in \mathbb{N} such that (4.4) has a solution \lambda_1 with \mathrm{Re}(\lambda_{1})\geq0 . Then

    \begin{eqnarray} 1 & = & \left|\frac{\lambda_{1}+d}{\lambda_{1}+d+f_{1u}+g_{1u}}\left(\frac{g_{1w}}{\lambda_{1}+\delta}e^{-m\tau}e^{-\lambda_{1}\tau} +\frac{N\delta f_{1v}}{(\lambda_{1}+\delta)(\lambda_{1}+c+\mu_{i_{1}}D_1)}e^{-m\tau}e^{-\lambda_{1}\tau}\right)\right|\\ & \lt & \left|\frac{g_{1w}}{\lambda_{1}+\delta}e^{-m\tau}e^{-\lambda_{1}\tau}\right| +\left|\frac{N\delta f_{1v}}{(\lambda_{1}+\delta)(\lambda_{1}+c+\mu_{i_{1}}D_1)}e^{-m\tau}e^{-\lambda_{1}\tau}\right|\\ & \lt & \frac{g_{1w}}{\delta}e^{-m\tau}+\frac{N f_{1v}}{c}e^{-m\tau}. \end{eqnarray} (4.5)

    However, from the steady state Eqs (3.1b) and (3.1c), we have

    \frac{g(u_1, w_1)}{\delta w_1}e^{-m\tau}+\frac{Nf(u_1, v_1)}{cv_1}e^{-m\tau} = 1.

    This and Assumption (A1) (v) together give us

    \frac{g_{1w}}{\delta}e^{-m\tau}+\frac{N f_{1v}}{c}e^{-m\tau}\leq\frac{g(u_1, w_1)}{\delta w_1}e^{-m\tau}+\frac{Nf(u_1, v_1)}{cv_1}e^{-m\tau} = 1,

    which is a contradiction with (4.5). This proves the claim and hence P_1 is locally asymptotically stable when R_1 < 1 < R_0 .

    Before studying the global stability of P_1 , we establish the persistence of infection.

    From the linearized system at P_0 (see (4.1)), we have the following cooperative system for (w, v) ,

    \begin{equation} \begin{cases} \frac {\partial w(x, t)}{\partial t} & = e^{-m\tau}\frac {\partial f(u_0, 0)}{\partial v} v(x, t-\tau)+e^{-m\tau}\frac{\partial g(u_0, 0)}{\partial w}w(x, t-\tau)-\delta w(x, t), \\ \frac{\partial v(x, t)}{\partial t} & = D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t). \end{cases} \end{equation} (4.6)

    With similar arguments as those for Lemma 3 and Lemma 4 in Lou and Zhao [45], we can obtain the following results.

    Lemma 4.4. There exists a principal eigenvalue \bar{\lambda}(u_0, \tau)\triangleq \bar{\lambda}(P_0, \tau) of (4.6) associated with a strongly positive eigenvector. Moreover, \bar{\lambda}(u_0, \tau) has the same sign as \lambda(u_0)\triangleq \bar{\lambda}(u_0, 0) .

    Lemma 4.5. R_0-1 and \lambda(u_0) have the same sign.

    Theorem 4.6. Suppose R_1\le 1 < R_0 . Then the infection is persistent, that is, there exists \varepsilon > 0 such that

    \liminf\limits_{t\to\infty}u(x, t, \phi)\geq\varepsilon, \quad \liminf\limits_{t\to\infty}w(x, t, \phi)\geq\varepsilon, \quad \liminf\limits_{t\to\infty}v(x, t, \phi)\geq\varepsilon

    uniformly for all x\in\overline{\Omega} , where \phi\in \mathcal{W}_1: = \left\{\phi\in\mathcal{C}^+ : w(\cdot, 0)\not\equiv0\ \mathit{\text{and}}\ v(\cdot, 0)\not\equiv0 \right\} .

    Proof. Define

    \partial\mathcal{W}_1: = \mathcal{C}^+ \setminus \mathcal{W}_1 = \left\{\phi\in\mathcal{C}^+ : w(\cdot, 0)\equiv0 \mbox{ or $v(\cdot, 0)\equiv0$} \right\}.

    By Lemma 2.2 and the second equation of (1.4), we know that \Phi(t)\mathcal{W}_1\subseteq\mathcal{W}_1 for all t\geq0 . Denote

    \mathcal{M}_{\partial}: = \left\{\phi\in\partial\mathcal{W}_1 : \Phi(t)\phi\in\partial\mathcal{W}_1 \mbox{ for $t\geq0$}\right\}.

    Claim 1. \omega(\phi) = \{(u_0, 0, 0, 0)\} for \phi\in\mathcal{M}_{\partial} , where \omega(\phi) is the omega limit set of the orbit \mathcal{O}^+(\phi): = \{\Phi(t)\phi:t\geq0\} .

    Since \phi\in\mathcal{M}_{\partial} , for all t\ge 0 , either w(x, t, \phi)\equiv0 or v(x, t, \phi)\equiv0 . If w(x, t, \phi)\equiv0 for all t\geq0 , then \lim\limits_{t\to\infty}v(x, t, \phi) = 0 uniformly for x\in\overline{\Omega} from the third equation of (1.4). Now, suppose that w(x, t_1, \phi)\not\equiv0 for some t_1\geq0 . Then by Lemma 2.2, w(x, t, \phi) > 0 for all t\geq t_1+\tau and x\in \Omega . Thus v(x, t, \phi)\equiv0 for all t\geq t_1+\tau . This, combined with the third equation of (1.4), implies that w(x, t, \phi)\equiv 0 for x\in \overline{\Omega} and t\ge t_1+\tau . Then, in either case, \lim\limits_{t\to\infty}v(x, t, \phi) = \lim\limits_{t\to\infty}w(x, t, \phi) = 0 uniformly for x\in \overline {\Omega} . Thus u is asymptotic to

    \frac{\partial u(x, t)}{\partial t} = s-du(x, t).

    By Corollary 4.3 in [52], we get \lim\limits_{t\to\infty}u(x, t, \phi) = u_0 uniformly for x\in\overline{\Omega} . The above discussion tells us that w(x, t, \phi)\equiv 0 for all t large enough. Then we can easily see from the fourth equation of (1.4) that \lim\limits_{t\to\infty}z(x, t, \phi) = 0 uniformly for x\in\overline{\Omega} . This proves \omega(\phi) = \{(u_0, 0, 0, 0)\} .

    Since R_1\leq1 < R_0 , by Lemma 4.4 and Lemma 4.5, there exists a sufficiently small \varepsilon_0 > 0 such that the following linear system

    \left\{ \begin{array}{rcl} \frac {\partial w(x, t)}{\partial t} & = &e^{-m\tau}\left(\frac {\partial f(u_0, 0)}{\partial v}-\varepsilon_0\right) v(x, t-\tau)+e^{-m\tau}\left(\frac{\partial g(u_0, 0)}{\partial w}-\varepsilon_0\right)w(x, t-\tau)-(\delta+p\varepsilon_0) w(x, t), \\ \frac{\partial v(x, t)}{\partial t} & = & D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t) \end{array} \right.

    has a positive principal eigenvalue \bar{\lambda}(u_0-\varepsilon_0) with positive eigenfunction (w_{\varepsilon_0}, v_{\varepsilon_0}) . By the continuity in Assumption (A1), there exists \delta_0\in (0, \varepsilon_0] such that

    \frac{\partial f(u, v)}{\partial v}\ge \frac {\partial f(u_0, 0)}{\partial v}-\varepsilon_0 \quad \mbox{and}\quad \frac{\partial g(u, w)}{\partial w}\ge \frac {\partial g(u_0, 0)}{\partial w}-\varepsilon_0

    for all u_0-\delta_0\le u\le u_0+\delta_0 , 0\le v\le \delta_0 , and 0\le w\le \delta_0 .

    Claim 2. \{(u_0, 0, 0, 0)\} is a uniform weak repeller for \mathcal{W}_1 in the sense that

    \limsup\limits_{t\to\infty}\|\Phi(t)\phi-(u_0, 0, 0, 0)\|\geq\delta_0 \mbox{ for } \phi\in\mathcal{W}_1.

    Suppose, by contradiction, there exists \phi_1^*\in\mathcal{W}_1 such that \limsup\limits_{t\to\infty}\|\Phi(t)\phi_1^*-(u_0, 0, 0, 0)\| < \delta_0 . Then there exists t_2 > 0 such that u(x, t, \phi_1^*) > u_0-\delta_0\ge u_0-\varepsilon_0 , w(x, t, \phi_1^*)\le \delta_0 , and v(x, t, \phi_1^*)\le \delta_0 for t\geq t_2 and x\in\overline{\Omega} . It follows from Assumption (A1) and the choice of \delta_0 that w and v satisfy

    \begin{equation} \left\{ \begin{array}{rcl} \frac{\partial w(x, t)}{\partial t} &\geq & e^{-m\tau}\left(\frac{\partial f(u_0, 0)}{\partial v}-\varepsilon_0\right)v(x, t-\tau) +e^{-m\tau}\left(\frac{\partial g(u_0, 0)}{\partial w}-\varepsilon_0\right) w(x, t-\tau)-(\delta +p\varepsilon_0)w(x, t), \\ & &t\geq t_2, \ x\in\Omega, \\ \frac{\partial v(x, t)}{\partial t} & = & D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t), \quad t\geq t_2, \ x\in\Omega, \\ \frac{\partial v(x, t)}{\partial \vec{n}}& = &0, \quad t\geq t_2, \quad x\in\partial\Omega. \end{array} \right. \end{equation} (4.7)

    Due to w(x, t, \phi_1^*) > 0 and v(x, t, \phi_1^*) > 0 for t > 0 and x\in \Omega , there exists \kappa_1 > 0 such that (w(x, t_2+\theta, \phi_1^*), v(x, t_2+\theta, \phi_1^*))\ge \kappa_1 e^{\bar{\lambda}(u_0-\varepsilon_0)(t_2+\theta)}(w_{\varepsilon_0}(x), v_{\varepsilon_0}(x)) for all x\in \Omega and \theta\in [-\tau, 0] . Then it follows from the comparison principle that w(x, t, \phi_1^*)\ge \kappa_1 e^{\bar{\lambda}(u_0-\varepsilon_0)t}w_{\varepsilon_0}(x) and v(x, t, \phi_1^*)\ge \kappa_1 e^{\bar{\lambda}(u_0-\varepsilon_0)t}v_{\varepsilon_0}(x) for all x\in \Omega and t\ge t_2 , a contradiction to the fact that both w(x, t, \phi_1^*) and v(x, t, \phi_1^*) are bounded. This proves Claim 2.

    Define a continuous function \mathcal{P}_1 : \mathcal{C}^+\rightarrow[0, \infty) by

    \mathcal{P}_1(\phi): = \min\left\{\min\limits_{x\in\overline{\Omega}}\phi_2(x, 0), \ \min\limits_{x\in\overline{\Omega}}\phi_3(x, 0)\right\}\ \mbox{for all $\phi\in\mathcal{C}^+$.}

    Clearly, \mathcal{P}_1^{-1}(0, \infty)\subset\mathcal{W}_1 , and \mathcal{P}_1 has the property that if \mathcal{P}_1(\phi) = 0 and \phi\in\mathcal{W}_1 or \mathcal{P}_1(\phi) > 0 , then \mathcal{P}_1(\Phi(t)\phi) > 0 for all t > 0 . Hence, \mathcal{P}_1 is a generalized distance function for the semiflow \Phi(t) [53]. According to the above discussions, we obtain that any forward orbit of \Phi(t) in \mathcal{M}_{\partial} converges to (u_0, 0, 0, 0) , which is isolated in \mathcal{C}^+ and \mathcal{W}^s(u_0, 0, 0, 0)\cap \mathcal{W}_1 = \varnothing , where \mathcal{W}^s(u_0, 0, 0, 0) is the stable manifold of (u_0, 0, 0, 0) . Moreover, there is no cycle in \partial\mathcal{W}_1 from (u_0, 0, 0, 0) to (u_0, 0, 0, 0) . Applying Theorem 3 in [53], we know that there exists an \bar{\varepsilon} > 0 such that \min\{\mathcal{P}_1(\phi)\} > \bar{\varepsilon} for any \phi\in\mathcal{W}_1 . It follows that

    \liminf\limits_{t\to\infty}w(x, t)\geq\bar{\varepsilon} \ \mbox{and}\ \liminf\limits_{t\to\infty}v(x, t)\geq\bar{\varepsilon} \ \mbox{uniformly for all $x\in\overline{\Omega}$.}

    This combined with Lemma 2.2 finishes the proof with \varepsilon = \min\{\bar{\varepsilon}, \frac{s}{d+\eta_1+\eta_2}\} .

    In order to study the global stability of P_1 , define G:(0, \infty)\ni x\to x-1-\ln x . Obviously, G(x) > 0 for x\in (0, \infty) and G attains its global minimum only at x = 1 . We also need the following assumption.

    \bf(A2) The nonlinear incidence functions f(u, v) and g(u, w) satisfy the following conditions.

    {\rm (i)} For any u > 0 ,

    \begin{cases} \frac{v}{v_1} \leq \frac{u_1f(u, v)}{uf(u_1, v_1)} \lt 1 \quad & \mbox{if }\ 0 \lt v \lt v_1, \\ 1 \leq \frac{u_1f(u, v)}{uf(u_1, v_1)} \lt \frac{v}{v_1} & \mbox{if } \ v_1 \lt v. \end{cases}

    {\rm (ii)} For any u > 0 ,

    \begin{cases} \frac{w}{w_1} \leq \frac{u_1g(u, w)}{ug(u_1, w_1)} \lt 1 \quad & \mbox{if }\ 0 \lt w \lt w_1, \\ 1 \leq \frac{u_1g(u, w)}{ug(u_1, w_1)} \lt \frac{w}{w_1} & \mbox{if }\ w_1 \lt w. \end{cases}

    Theorem 4.7. Suppose that R_1\leq1 < R_0 and Assumption (A2) are satisfied. Then the immunity-free infected steady state P_1 is globally attractive in

    \mathcal{C}_1^+ = \{\phi\in\mathcal{C}^+|\mathit{\mbox{there exists $t_3\in\mathbb{R}^+$ such that $w(\cdot, t_3, \phi)\not\equiv0$ or $v(\cdot, t_3, \phi)\not\equiv0$}}\}.

    In particular, it is globally asymptotically stable in \mathcal{C}_1^+ if further R_1 < 1 .

    Proof. According to Lemma 2.2 and Theorem 4.6, we know that there exists \varepsilon > 0 such that \liminf\limits_{t\to\infty}u(x, t, \phi)\geq\varepsilon , \liminf\limits_{t\to\infty}w(x, t, \phi)\geq\varepsilon , and \liminf\limits_{t\to\infty}v(x, t, \phi)\geq\varepsilon for \phi\in\mathcal{C}_1^+ . Without loss of generality, we define a Lyapunov functional

    L(t) = \int_{\Omega}L(x, t)\mathrm{d}x,

    where

    \begin{eqnarray*} L(x, t) & = & u_1G\left(\frac{u(x, t)}{u_1}\right)+e^{m\tau}w_1G\left(\frac{w(x, t)}{w_1}\right)+ \frac{f(u_1, v_1)}{cv_1}v_1G\left(\frac{v(x, t)}{v_1}\right)+\frac{p}{q}e^{m\tau}z(x, t)\\ & & +f(u_1, v_1)\int_{t-\tau}^{t}G\left(\frac{f(u(x, \theta), v(x, \theta))}{f(u_1, v_1)}\right)\mathrm {d}\theta +g(u_1, w_1)\int_{t-\tau}^{t}G\left(\frac{g(u(x, \theta), w(x, \theta))}{g(u_1, w_1)}\right)\mathrm {d}\theta. \end{eqnarray*}

    Calculating the time derivative of L(x, t) along solutions of (1.4) yields

    \begin{eqnarray*} \frac{\partial L(x, t)}{\partial t}& = & \left(1-\frac{u_1}{u(x, t)}\right)[s-du(x, t)-f(u, v)-g(u, w)]+\frac{p}{q}e^{m\tau}[qw(x, t)z(x, t)-bz(x, t)]\\ && +e^{m\tau}\left(1-\frac{w_1}{w(x, t)}\right)\left[e^{-m\tau}f(u_{\tau}, v_{\tau})+e^{-m\tau}g(u_{\tau}, w_{\tau})-\delta w(x, t)-pw(x, t)z(x, t)\right]\\ && + \frac{f(u_1, v_1)}{cv_1}\left(1-\frac{v_1}{v(x, t)}\right)[D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t)]+\frac{p}{q}e^{m\tau}D_2\Delta z(x, t)\\ && + f(u, v)-f(u_{\tau}, v_{\tau})+f(u_1, v_1)\ln{\frac{f(u_{\tau}, v_{\tau})}{f(u, v)}}\\ && +g(u, w)-g(u_{\tau}, w_{\tau})+g(u_1, w_1)\ln{\frac{g(u_{\tau}, w_{\tau})}{g(u, w)}}. \end{eqnarray*}

    With

    s = du_1+f(u_1, v_1)+g(u_1, w_1), \ \delta = \frac{e^{-m\tau}f(u_1, v_1)+e^{-m\tau}g(u_1, w_1)}{w_1}, \ c = \frac{N\delta w_1}{v_1},

    we have

    \begin{eqnarray*} && \frac{\partial L(x, t)}{\partial t}\\ & = & du_1\left(2-\frac{u_1}{u(x, t)}-\frac{u(x, t)}{u_1}\right)+\frac{p}{q}e^{m\tau}(R_1-1)z(x, t)\\ && + \frac{f(u_1, v_1)}{cv_1}\left(1-\frac{v_1}{v(x, t)}\right)D_1\Delta v(x, t)+\frac{p}{q}e^{m\tau}D_2\Delta z(x, t)-f(u_1, v_1)\frac{w(x, t)v_1}{v(x, t)w_1}\\ && + f(u_1, v_1)\left[3-\frac{u_1}{u(x, t)}-\frac{f(u_{\tau}, v_{\tau})w_1}{f(u_1, v_1)w(x, t)}+ \frac{f(u, v)u_1}{f(u_1, v_1)u(x, t)}-\frac{v(x, t)}{v_1}+\ln{\frac{f(u_{\tau}, v_{\tau})}{f(u, v)}}\right]\\ && + g(u_1, w_1) \left[2-\frac{u_1}{u(x, t)}-\frac{g(u_{\tau}, w_{\tau})w_1}{g(u_1, w_1)w(x, t)} +\frac{g(u, w)u_1}{g(u_1, w_1)u(x, t)}-\frac{w(x, t)}{w_1}+\ln{\frac{g(u_{\tau}, w_{\tau})}{g(u, w)}}\right]\\ & = & du_1\left(2-\frac{u_1}{u(x, t)}-\frac{u(x, t)}{u_1}\right)+\frac{p}{q}e^{m\tau}(R_1-1)z(x, t)\\ && + \frac{f(u_1, v_1)}{cv_1}\left(1-\frac{v_1}{v(x, t)}\right)D_1\Delta v(x, t)+\frac{p}{q}e^{m\tau}D_2\Delta z(x, t)-f(u_1, v_1)G\left(\frac{w(x, t)v_1}{v(x, t)w_1}\right)\\ && + f(u_1, v_1) \left[G\left(\frac{f(u, v)u_1}{f(u_1, v_1)u(x, t)}\right)-G\left(\frac{v(x, t)}{v_1}\right) -G\left(\frac{u_1}{u(x, t)}\right)-G\left(\frac{f(u_{\tau}, v_{\tau})w_1}{f(u_1, v_1)w(x, t)}\right)\right]\\ && + g(u_1, w_1)\left[G\left(\frac{g(u, w)u_1}{g(u_1, w_1)u(x, t)}\right)-G\left(\frac{w(x, t)}{w_1}\right) -G\left(\frac{u_1}{u(x, t)}\right)-G\left(\frac{g(u_{\tau}, w_{\tau})w_1}{g(u_1, w_1)w(x, t)}\right)\right]. \end{eqnarray*}

    Clearly,

    \int_{\Omega}du_1\left(2-\frac{u(x, t)}{u_1}-\frac{u_1}{u(x, t)}\right)\mathrm {d}x\leq0

    and

    \int_{\Omega}\left(\frac{pb}{q}e^{m\tau}(R_1-1)z(x, t)\right)\mathrm {d}x\leq0 \ \mbox{since $R_1\le 1$.}

    Using the Divergence Theorem and the homogeneous Neumann boundary conditions of (1.5), we have

    \int_{\Omega}\Delta v(x, t)\mathrm{d}x = \int_{\partial\Omega}\frac{\partial v(x, t)}{\partial\vec{n}}\mathrm {d}x = 0, \quad \int_{\Omega}\Delta z(x, t)\mathrm{d}x = \int_{\partial\Omega}\frac{\partial z(x, t)}{\partial\vec{n}}\mathrm {d}x = 0,

    and

    \begin{eqnarray*} 0& = & \int_{\partial\Omega}\frac{1}{v(x, t)}\nabla v(x, t)\ \cdot\ \vec{n} \ \mathrm {d}x \\ & = & \int_{\Omega}\nabla\left(\frac{1}{v(x, t)}\nabla v(x, t)\right)\mathrm {d}x \\ & = & \int_{\Omega}\left(\frac{1}{v(x, t)}\Delta v(x, t)-\frac{1}{v^2(x, t)}\|\nabla v(x, t)\|^2 \right)\mathrm {d}x. \end{eqnarray*}

    The latter gives

    \int_{\Omega}\frac{1}{v(x, t)}\Delta v(x, t)\mathrm {d}x = \int_{\Omega}\frac{1}{v^2(x, t)}\|\nabla v(x, t)\|^2\mathrm{d}x\geq0

    and hence

    \int_{\Omega}\frac{f(u_1, v_1)}{cv_1}\left(1-\frac{v_1}{v(x, t)}\right)D_1\Delta v(x, t)\mathrm{d}x\leq 0.

    To summarize, we have obtained

    \begin{eqnarray*} \frac{\mathrm {d} L(t)}{\mathrm{d} t}&\leq& f(u_1, v_1) \int_{\Omega}\left[G\left(\frac{f(u, v)u_1}{f(u_1, v_1)u(x, t)}\right)-G\left(\frac{v(x, t)}{v_1}\right)\right]\mathrm{d}x \\ && + g(u_1, w_1) \int_{\Omega}\left[G\left(\frac{g(u, w)u_1}{g(u_1, w_1)u(x, t)}\right)-G\left(\frac{w(x, t)}{w_1}\right)\right]\mathrm{d}x \\ && - g(u_1, w_1) \int_{\Omega}\left[G\left(\frac{u_1}{u(x, t)}\right) +G\left(\frac{g(u_{\tau}, w_{\tau})w_1}{g(u_1, w_1)w(x, t)}\right)\right]\mathrm{d}x \\ && - f(u_1, v_1) \int_{\Omega}\left[G\left(\frac{u_1}{u(x, t)}\right)+G\left(\frac{f(u_{\tau}, v_{\tau})w_1}{f(u_1, v_1)w(x, t)}\right) +G\left(\frac{w(x, t)v_1}{v(x, t)w_1}\right)\right]\mathrm {d}x. \end{eqnarray*}

    Note that the monotonicity of G(x) on each side of x = 1 and Assumption (A2) give us

    G\left(\frac{f(u, v)u_1}{f(u_1, v_1)u(x, t)}\right) \leq G\left(\frac{v(x, t)}{v_1}\right) \quad \mbox{and} \quad G\left(\frac{g(u, w)u_1}{g(u_1, w_1)u(x, t)}\right) \leq G\left(\frac{w(x, t)}{w_1}\right).

    Thus

    \begin{eqnarray*} \frac{\mathrm{d} L(t)}{\mathrm{d} t} & \leq & -g(u_1, w_1) \int_{\Omega} \left[G\left(\frac{u_1}{u(x, t)}\right) +G\left(\frac{g(u_{\tau}, w_{\tau})w_1}{g(u_1, w_1)w(x, t)}\right)\right] \mathrm{d}x \\ && - f(u_1, v_1)\int_{\Omega}\left[G\left(\frac{u_1}{u(x, t)}\right) +G\left(\frac{f(u_{\tau}, v_{\tau})w_1}{f(u_1, v_1)w(x, t)}\right) +G\left(\frac{w(x, t)v_1}{v(x, t)w_1}\right)\right] \mathrm{d}x \\ & \le & 0. \end{eqnarray*}

    Moreover, \frac{\mathrm {d} L(t)}{\mathrm {d} t} = 0 if and only if u(x, t) = u_1 , w(x, t) = w_1 , v(x, t) = v_1 , and z(x, t) = 0 . Then the largest invariant subset of \frac{\mathrm{d} L(t)}{\mathrm {d} t} = 0 is \{P_1\} . By LaSalle's Invariance Principle (see Theorem 5.3.1 in [50] or Theorem 3.4.7 in [51]), the immunity-free infected steady state P_1 is globally attractive in \mathcal{C}_1^+ when R_1\leq1 < R_0 . This, together with Theorem 4.3, implies that P_1 is globally asymptotically stable in \mathcal{C}_1^+ if further R_1 < 1 .

    For convenience of notations, denote

    \begin{array}{lll} f_2 = f(u_2, v_2), \quad & g_2 = g(u_2, w_2), \quad & f_{2u} = \frac{\partial f(u_2, v_2)}{\partial u}, \\ f_{2v} = \frac{\partial f(u_2, v_2)}{\partial v}, &g_{2u} = \frac{\partial g(u_2, w_2)}{\partial u}, & g_{2w} = \frac{\partial g(u_2, w_2)}{\partial w}. \end{array}

    Theorem 4.8. If R_1 > 1 , then the infected-immune steady state P_2 is locally asymptotically stable.

    Proof. According to (4.2), the characteristic equation at P_2 is

    \begin{eqnarray} \begin{aligned} 0 = & (\lambda+d+f_{2u}+g_{2u})(\lambda+c+\mu_iD_1)(\lambda+b-qw_2+\mu_iD_2)(\lambda+\delta+pz_2) \\ &+qw_2pz_2(\lambda+c+\mu_iD_1)(\lambda+d+f_{2u}+g_{2u}) \\ &-(\lambda+d)(\lambda+c+\mu_iD_1)(\lambda+b-qw_2+\mu_iD_2)g_{2w}e^{-(\lambda+m)\tau} \\ & -(\lambda+d)(\lambda+b-qw_2+\mu_iD_2)N\delta f_{2v}e^{-(\lambda+m)\tau}. \end{aligned} \end{eqnarray} (4.8)

    We claim that all roots of (4.8) have negative real parts. Otherwise, suppose for some i_2\in \mathbb{N} , it has a root \lambda_2 with \mathrm {Re}(\lambda_2)\ge 0 . Since w_2 = \frac {b}{q} , we have

    \begin{eqnarray*} 1 = \frac{(\lambda_2+\mu_{i_2}D_2)(\lambda_2+d)}{(\lambda_2+d+f_{2u}+g_{2u}) \big[(\lambda_2+\mu_{i_2}D_2)(\lambda_2+\delta+pz_2)+pbz_2\big]}\left(g_{2w}e^{-(\lambda_2+m)\tau}+\frac{N\delta f_{2v}e^{-(\lambda_2+m)\tau}}{\lambda_2+c+\mu_{i_2}D_1}\right), \end{eqnarray*}

    which implies

    1 \lt \left|\frac{\lambda_2+\mu_{i_2}D_2}{(\lambda_2+\mu_{i_2}D_2)(\lambda_2+\delta+pz_2)+pbz_2}\right| \times\left(g_{2w}e^{-m\tau}+\frac{N\delta f_{2v}e^{-m\tau}}{c}\right).

    With similar arguments as those in the proof of Theorem 4.3, we can obtain

    g_{2w}e^{-m\tau}+\frac{N\delta f_{2v}e^{-m\tau}}{c}\le \delta+pz_2.

    Thus we have arrived at

    1 \lt \left|\frac{\lambda_2+\mu_{i_2}D_2}{(\lambda_2+\mu_{i_2}D_2)(\lambda_2+\delta+pz_2)+pbz_2}\right|\times(\delta+pz_2),

    which is impossible as one can check that |(\lambda_2+\mu_{i_2}D_2)(\lambda_2+\delta+pz_2)+pbz_2| > |(\lambda_2+\mu_{i_2}D_2)(\delta+pz_2)| .

    This completes the proof.

    To establish the global stability of P_2 , we need the persistence of immunity.

    From the linearized system at P_1 (see (4.1)), we have the following cooperative system for (w, v, z) ,

    \begin{equation} \begin{cases} \frac {\partial w(x, t)}{\partial t} & = e^{-m\tau}\frac {\partial f(u_1, v_1)}{\partial v} v(x, t-\tau)+e^{-m\tau}\frac{\partial g(u_1, w_1)}{\partial w}w(x, t-\tau)-(\delta +pz_1)w(x, t)-pw_1z(x, t), \\ \frac{\partial v(x, t)}{\partial t} & = D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t), \\ \frac{\partial z(x, t)}{\partial t} & = D_2\Delta z(x, t)+qz_1w(x, t)+(qw_1-b)z(x, t). \end{cases} \end{equation} (4.9)

    With similar arguments as those for Lemma 3 and Lemma 4 in Lou and Zhao [45], we can obtain the following results.

    Lemma 4.9. There exists a principal eigenvalue \hat{\lambda}(P_1, \tau) of (4.9) associated with a strongly positive eigenvector. Moreover, \hat{\lambda}(P_1, \tau) has the same sign as \hat{\lambda}(P_1, 0) .

    Lemma 4.10. R_1-1 and \hat{\lambda}(P_1, 0) have the same sign.

    Theorem 4.11. Suppose that R_1 > 1 (it is necessary that R_0 > 1 ) and (A2) holds. Then the immunity is persistent, that is, there exists \epsilon > 0 such that

    \liminf\limits_{t\to\infty}u(x, t, \phi)\geq\epsilon, \quad \liminf\limits_{t\to\infty}w(x, t, \phi)\geq\epsilon, \quad \liminf\limits_{t\to\infty}v(x, t, \phi)\geq\epsilon, \quad \liminf\limits_{t\to\infty}z(x, t, \phi)\geq\epsilon

    uniformly for all x\in\overline{\Omega} , where \phi\in \mathcal{W}_2: = \left\{\phi\in\mathcal{C}^+ : w(\cdot, 0)\not\equiv0, \ v(\cdot, 0)\not\equiv0, \ { and}\ z(\cdot, 0)\not\equiv0 \right\} .

    Proof. The proof is quite similar to that of Theorem 4.6. Denote

    \partial\mathcal{W}_2: = \mathcal{C}^+\setminus \mathcal{W}_2 = \left\{\phi\in\mathcal{C}^+ : w(\cdot, 0)\equiv0, \mbox{ or } v(\cdot, 0)\equiv0, \mbox { or } z(\cdot, 0)\equiv0 \right\}.

    Set M_0 = \{P_0\} and M_1 = \{P_1\} .

    According to Lemma 2.2, we know that w(x, t, \phi) > 0 , v(x, t, \phi) > 0 , and z(x, t, \phi) > 0 for all t > 0 and x\in\Omega , \phi\in\mathcal{W}_2 , which implies that \Phi(t)\mathcal{W}_2\subseteq\mathcal{W}_2 for all t\geq0 . Define

    \mathcal{M}^*_{\partial}: = \left\{\phi\in\partial\mathcal{W}_2 : \Phi(t)\phi\in\partial\mathcal{W}_2 \ \text{for}\ t\geq0\right\}.

    Claim 3. Let \phi\in \mathcal{M}^*_{\partial} . Then \omega(\phi) = M_0 or M_1 .

    Sine \phi\in \mathcal{M}^*_{\partial} , for any t\ge 0 , we have either w(x, t, \phi) \equiv 0 , or v(x, t, \phi)\equiv 0 , or z(x, t, \phi)\equiv 0 . If z(x, t_4, \phi)\not\equiv0 for some t_4\ge 0 , then by Lemma 2.2, z(x, t, \phi) > 0 for t > t_4 and x\in \Omega . Then either w(x, t, \phi)\equiv0 or v(x, t, \phi)\equiv 0 for each t > t_4 . By the proof of Claim 1, we know that \omega(\phi) = M_0 . Now, suppose that z(x, t, \phi)\equiv 0 for all t\ge 0 . If for each t\ge 0 , either w(x, t, \phi)\equiv0 or v(x, t, \phi)\equiv 0 , then by Claim 1, \omega(\phi) = M_0 . If there exists \tilde{t}\ge 0 such that w(x, \tilde{t}, \phi)\not \equiv 0 and v(x, \tilde{t}, \phi)\not\equiv 0 . Then by Theorem 4.6, there exists \xi > 0 such that

    \liminf\limits_{t\to\infty}w(x, t, \phi)\ge \xi\ \mbox{and}\ \liminf\limits_{t\to\infty}v(x, t, \phi)\ge \xi \ \mbox{uniformly in $\overline{\Omega}$.}

    Now consider the reduced system of (1.4) with z = 0 . Modifying the Lyapunov functional L(t) in the proof of Theorem 4.7 by ignoring the term \frac{p}{q}e^{m\tau}z(x, t) in L(x, t) , we can show that the solution of the reduced system converges to (u_1, w_1, v_1) and hence \omega(\phi) = M_1 . This proves Claim 3.

    Claim 4. Both M_0 and M_1 are uniform weak repellers for \mathcal{W}_2 . Since \mathcal{W}_2\subset \mathcal{W}_1 , by Claim 2, M_0 is a uniform repeller for \mathcal{W}_2 . The proof of M_1 being a uniform repeller of \mathcal{W}_2 is similar as that of Claim 2 by using Lemma 4.9 and Lemma 4.10. Therefore, we omit the detail here.

    Define a continuous function \mathcal{P}_2 : \mathcal{C}^+\rightarrow[0, \infty) by

    \mathcal{P}_2(\phi): = \min\left\{\min\limits_{x\in\overline{\Omega}}\phi_2(x, 0), \ \min\limits_{x\in\overline{\Omega}}\phi_3(x, 0), \ \min\limits_{x\in\overline{\Omega}}\phi_4(x, 0)\right\} \ \mbox{for $\phi\in\mathcal{C}^+$.}

    It is easy to see that \mathcal{P}_2^{-1}(0, \infty)\subset\mathcal{W}_2 , and \mathcal{P}_2 has the property that if \mathcal{P}_2(\phi) = 0 and \phi\in\mathcal{W}_2 or \mathcal{P}_2(\phi) > 0 , then \mathcal{P}_2(\Phi(t)\phi) > 0 for all t > 0 . Thus \mathcal{P}_2 is a generalized distance function for the semiflow \Phi(t) . As M_0 and M_1 are repellers, we know that both M_0 and M_1 are isolated, and \mathcal{W}^s(M_i)\cap \mathcal{W}_2 = \emptyset for i = 0 and 1 . Moreover, no subset of \{M_0, M_1\} forms a cycle in \partial \mathcal{W}_2 . By Smith and Zhao [53,Theorem 3], there exists a \bar{\epsilon} > 0 such that \min\{\mathcal{P}_2(\phi)\} > \bar{\epsilon} for any \phi\in\mathcal{W}_2 . Then as for Theorem 4.6, with \varepsilon = \min\{\bar{\epsilon}, \frac {s}{d+\eta_1+\eta_2}\} finishes the proof.

    As for the global stability of P_2 , we make the following assumption to establish the global stability of P_2 .

    \bf(A3) The nonlinear incidence functions f(u, v) and g(u, w) satisfy the following conditions.

    {\rm (i)} For any u > 0 ,

    \left\{ \begin{array}{ll} \frac{v}{v_2} \leq \frac{u_2f(u, v)}{uf(u_2, v_2)} \lt 1 \quad &\mbox{if }\ 0 \lt v \lt v_2, \\ 1 \leq \frac{u_2f(u, v)}{uf(u_2, v_2)} \lt \frac{v}{v_2}\quad & \mbox{if }\ v_2 \lt v. \end{array} \right.

    {\rm (ii)} For any u > 0 ,

    \left\{ \begin{array}{ll} \frac{w}{w_2} \leq \frac{u_2g(u, w)}{ug(u_2, w_2)} \lt 1 \quad & \mbox{if }\ 0 \lt w \lt w_2, \\ 1 \leq \frac{u_2g(u, w)}{ug(u_2, w_2)} \lt \frac{w}{w_2} \quad& \mbox{if }\ w_2 \lt w. \end{array} \right.

    Theorem 4.12. Suppose that R_1 > 1 and Assumptions (A2) and (A3) are satisfied. Then the infected-immune steady state P_2 is globally asymptotically stable in

    \mathcal{C}_2^+ = \left \{\phi\in\mathcal{C}^+\left | \begin{array}{ll} \mathit{\mbox{there exists $t_5\in\mathbb{R}^+$ such that either $w(\cdot, t_5, \phi)\not\equiv0$ or $v(\cdot, t_5, \phi)\not\equiv0$, }}\\ \mathit{\mbox{there exists $t_6\in \mathbb{R}^+$ such that $z(\cdot, t_6, \phi)\not\equiv0$}} \end{array} \right. \right\}.

    Proof. It follows from Lemma 2.2 and Theorem 4.11 that there exists an \varepsilon > 0 such that

    \liminf\limits_{t\to\infty} u(x, t, \phi)\ge \varepsilon, \ \liminf\limits_{t\to\infty}w(x, t, \phi)\ge \varepsilon, \ \liminf\limits_{t\to\infty}v(x, t, \phi)\ge \varepsilon, \ \liminf\limits_{t\to\infty}z(x, t, \phi)\ge \varepsilon

    unfiormly in \overline{\Omega} and \phi\in \mathcal{C}_2^+ . Without loss of generality, we define a Lyapunov functional

    I(t) = \int_{\Omega}I(x, t)\mathrm{d}x,

    where

    \begin{eqnarray*} I(x, t)& = & u_2G\left(\frac{u(x, t)}{u_2}\right)+e^{m\tau}w_2G\left(\frac{w(x, t)}{w_2}\right) +\frac{f(u_2, v_2)}{cv_2}v_2G\left(\frac{v(x, t)}{v_2}\right)+\frac{p}{q}e^{m\tau}z_2G\left(\frac{z(x, t)}{z_2}\right) \\ &&+f(u_2, v_2)\int_{t-\tau}^{t}G\left(\frac{f(u(x, \theta), v(x, \theta))}{f(u_2, v_2)}\right) \mathrm {d}\theta+g(u_2, w_2)\int_{t-\tau}^{t}G\left(\frac{g(u(x, \theta), w(x, \theta))}{g(u_2, w_2)}\right)\mathrm{d}\theta. \end{eqnarray*}

    Calculate the time derivative of I(x, t) along the solutions of (1.4) to get

    \begin{eqnarray*} \frac{\partial I(x, t)}{\partial t}& = &\left(1-\frac{u_2}{u(x, t)}\right)[s-du(x, t)-f(u, v)-g(u, w)] \\ && + \frac{p}{q}e^{m\tau}\left(1-\frac{z_2}{z(x, t)}\right)[qw(x, t)z(x, t)-bz(x, t)] \\ && +e^{m\tau}\left(1-\frac{w_2}{w(x, t)}\right)\left[e^{-m\tau}f(u_{\tau}, v_{\tau})+e^{-m\tau}g(u_{\tau}, w_{\tau})-\delta w(x, t)-pw(x, t)z(x, t)\right] \\ && +\frac{f(u_2, v_2)}{cv_2}\left(1-\frac{v_2}{v(x, t)}\right)[D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t)] \\ && +\frac{p}{q}e^{m\tau}\left(1-\frac{z_2}{z(x, t)}\right)D_2\Delta z(x, t) +f(u, v)-f(u_{\tau}, v_{\tau})+f(u_2, v_2)\ln{\frac{f(u_{\tau}, v_{\tau})}{f(u, v)}} \\ && +g(u, w)-g(u_{\tau}, w_{\tau})+g(u_2, w_2)\ln{\frac{g(u_{\tau}, w_{\tau})}{g(u, w)}}. \end{eqnarray*}

    With the following relations,

    \begin{eqnarray*} s & = & du_2+f(u_2, v_2)+g(u_2, w_2), \\ \delta & = & \frac{e^{-m\tau}f(u_2, v_2)+e^{-m\tau}g(u_2, w_2)-pw_2z_2}{w_2}, \\ c & = & \frac{N\delta w_2}{v_2}, \\ w_2 & = & \frac{b}{q}, \end{eqnarray*}

    we get

    \begin{eqnarray*} \frac{\partial I(x, t)}{\partial t}& = & du_2\left(2-\frac{u_2}{u(x, t)}-\frac{u(x, t)}{u_2}\right)+ \frac{f(u_2, v_2)}{cv_2}\left(1-\frac{v_2}{v(x, t)}\right)D_1\Delta v(x, t) \\ && +\frac{p}{q}e^{m\tau}\left(1-\frac{z_2}{z(x, t)}\right)D_2\Delta z(x, t)-f(u_2, v_2)\frac{w(x, t)v_2}{v(x, t)w_2} \\ && + f(u_2, v_2)\left[3-\frac{u_2}{u(x, t)}-\frac{f(u_{\tau}, v_{\tau})w_2}{f(u_2, v_2)w(x, t)} +\frac{f(u, v)u_2}{f(u_2, v_2)u(x, t)}-\frac{v(x, t)}{v_2}+\ln{\frac{f(u_{\tau}, v_{\tau})}{f(u, v)}}\right] \\ && +g(u_2, w_2)\left[2-\frac{u_2}{u(x, t)}-\frac{g(u_{\tau}, w_{\tau})w_2}{g(u_2, w_2)w(x, t)} +\frac{g(u, w)u_2}{g(u_2, w_2)u(x, t)}-\frac{w(x, t)}{w_2}+\ln{\frac{g(u_{\tau}, w_{\tau})}{g(u, w)}}\right]. \end{eqnarray*}

    Then

    \begin{eqnarray*} \frac {\mathrm{d}I(t)}{\mathrm {d}t} & = & \int_{\Omega}du_2\left(2-\frac{u_2}{u(x, t)}-\frac{u(x, t)}{u_2}\right) \mathrm{d}x +\int_{\Omega}\frac{f(u_2, v_2)}{cv_2}\left(1-\frac{v_2}{v(x, t)}\right)D_1\Delta v(x, t) \mathrm{d}x \\ && +\int_{\Omega}\frac{p}{q}e^{m\tau}\left(1-\frac{z_2}{z(x, t)}\right)D_2\Delta z(x, t) \mathrm{d}x -f(u_2, v_2)\int_{\Omega}G\left(\frac{w(x, t)v_2}{v(x, t)w_2}\right)\mathrm{d}x \\ && + f(u_2, v_2)\int_{\Omega}\left[G\left(\frac{f(u, v)u_2}{f(u_2, v_2)u(x, t)}\right)-G\left(\frac{v(x, t)}{v_2}\right) -G\left(\frac{u_2}{u(x, t)}\right)-G\left(\frac{f(u_{\tau}, v_{\tau})w_2}{f(u_2, v_2)w(x, t)}\right)\right]\mathrm{d}x \\ && +g(u_2, w_2)\int_{\Omega}\left[G\left(\frac{g(u, w)u_2}{g(u_2, w_2)u(x, t)}\right)-G\left(\frac{w(x, t)}{w_2}\right) -G\left(\frac{u_2}{u(x, t)}\right)-G\left(\frac{g(u_{\tau}, w_{\tau})w_2}{g(u_2, w_2)w(x, t)}\right)\right]\mathrm{d}x. \end{eqnarray*}

    Similarly as in the proof of Theorem 4.7, we can show

    \begin{eqnarray*} \int_{\Omega}du_2\left(2-\frac{u(x, t)}{u_2}-\frac{u_2}{u(x, t)}\right)\mathrm{d}x& \leq &0, \\ \int_{\Omega}\frac{f(u_2, v_2)}{cv_2}\left(1-\frac{v_2}{v(x, t)}\right)D_1\Delta v(x, t)\mathrm{d}x & \leq & 0, \\ \int_{\Omega}\frac{p}{q}e^{m\tau}\left(1-\frac{z_2}{z(x, t)}\right)D_2\Delta z(x, t) \mathrm{d}x & \leq & 0, \\ G\left(\frac{f(u, v)u_2}{f(u_2, v_2)u(x, t)}\right) & \leq & G\left(\frac{v(x, t)}{v_2}\right), \\ G\left(\frac{g(u, w)u_2}{g(u_2, w_2)u(x, t)}\right) & \leq & G\left(\frac{w(x, t)}{w_2}\right). \end{eqnarray*}

    Therefore, we have \frac{\mathrm{d}I(t)}{\mathrm {d}t}\le 0 . Moreover, \frac{\mathrm{d} I(t)}{\mathrm{d} t} = 0 if and only if u(x, t) = u_2 , w(x, t) = w_2 , v(x, t) = v_2 , z(x, t) = z_2 . Then the largest invariant subset of \frac{\mathrm{d} I(t)}{\mathrm{d} t} = 0 is \{P_2\} . By LaSalle's Invariance Principle (see Theorem 5.3.1 in [50] or Theorem 3.4.7 in [51]), the infected-immune steady state P_2 is globally attractive in \mathcal{C}_2^+ . This, together with Theorem 4.8, implies the global asymptotic stability of P_2 in \mathcal{C}_2^+ .

    In this section, we perform some numerical simulations to illustrate the results obtained in section 4. Let f(u, v) = \frac{\beta_1 uv}{1+\alpha_1v} and g(u, w) = \frac{\beta_2 uw}{1+\alpha_2w} . One can easily verify that f and g satisfy (A1)–(A3). Then the model (1.4) becomes

    \begin{eqnarray} \begin{cases} \frac{\partial u(x, t)}{\partial t} = s-du(x, t)-\frac{\beta_1 u(x, t)v(x, t)}{1+\alpha_1v(x, t)}-\frac{\beta_2 u(x, t)w(x, t)}{1+\alpha_2w(x, t)}, \ x\in\Omega, \ t \gt 0, \\ \frac{\partial w(x, t)}{\partial t} = e^{-m\tau}\left(\frac{\beta_1 u(x, t-\tau)v(x, t-\tau)}{1+\alpha_1v(x, t-\tau)}+\frac{\beta_2 u(x, t-\tau)w(x, t-\tau)}{1+\alpha_2w(x, t-\tau)}\right) \\ \qquad\qquad\quad -\delta w(x, t)-pw(x, t)z(x, t), \ x\in\Omega, \ t \gt 0, \\ \frac{\partial v(x, t)}{\partial t} = D_1\Delta v(x, t)+N\delta w(x, t)-cv(x, t), \ x\in\Omega, \ t \gt 0, \\ \frac{\partial z(x, t)}{\partial t} = D_2\Delta z(x, t)+qw(x, t)z(x, t)-bz(x, t), \ x\in\Omega, \ t \gt 0, \end{cases} \end{eqnarray} (5.1)

    subject to the homogeneous Neumann boundary conditions

    \begin{eqnarray*} \frac{\partial v}{\partial \vec{n}} = 0, \ \frac{\partial z}{\partial \vec{n}} = 0, \ x\in\partial\Omega, \ t \gt 0. \end{eqnarray*}

    For (5.1), the basic reproduction number of infection is given by

    \begin{eqnarray*} R_0 = \frac{N\beta_1s}{cd}e^{-m\tau}+\frac{\beta_2s}{\delta d}e^{-m\tau} \end{eqnarray*}

    and the basic reproduction number of immunity is given by

    R_1 = \dfrac{qw_1}{b} = \dfrac{q\left(B_2+\sqrt{B_2^2-4B_1B_3}\right)}{2bB_1},

    where

    \begin{array}{ll} B_1 = N\delta^2e^{m\tau}(\beta_1\alpha_2+\beta_2\alpha_1+\alpha_1\alpha_2 d) \gt 0, \\ B_2 = \delta de^{m\tau}(N\delta\alpha_1+c\alpha_2)-N\delta\beta_1\alpha_2 s+N\delta^2\beta_1e^{m\tau}-N\delta\beta_2\alpha_1 s+\delta e^{m\tau}\beta_2c, \\ B_3 = -dc\delta e^{m\tau}(R_0-1) \lt 0\ \mbox{since}\ R_0 \gt 1. \end{array}

    For simulations, we take \alpha_1 = 0.01, \ \alpha_2 = 0.01 , D_1 = 0.0017 , D_2 = 0.0001 , and the values of the other parameters are summarized in Table 1. Moreover, \Omega = [0, 4] and the initial condition used is

    \begin{array}{ll} u(x, \theta) = 23+0.2\cos\frac{\pi x}{2}, \quad & w(x, \theta) = 0.7+0.2\cos\frac{\pi x}{2}, \\ v(x, \theta) = 3+0.2\cos\frac{\pi x}{2}, & z(x, \theta) = 2+0.2\cos\frac{\pi x}{2} \end{array}
    Table 1.  Parameter values for simulation.
    Parameters Ranges value Units References
    s 0\sim10 10 \rm{cells\; ml^{-1}day^{-1}} [54]
    d 0.0001\sim0.2 0.01 {\rm day}^{-1} [22]
    \beta_1 4.6\times10^{-8}\sim0.5 variable \rm{ml^{-1}day^{-1}} [55]
    \beta_2 1\times10^{-5}\sim0.7 2.4\times10^{-5} \rm{ml^{-1}day^{-1}} [39]
    m \alpha\in[d, \delta] 0.05 {\rm day}^{-1} [55]
    \tau 0\sim1.5 0.5 \rm days [55]
    \delta 0.00019\sim1.4 1 {\rm day}^{-1} [22]
    p 0.0001\sim4.048 0.024 \rm{ml^{-1}day^{-1}} [22,55]
    N 6.25\sim23599.9 2000 {\rm viron\; cells}^{-1} [22]
    c 2\sim36 23 {\rm day}^{-1} [39,22]
    q 0.0051\sim3.912 0.15 {\rm day}^{-1} [22]
    b 0.004\sim8.087 0.5 {\rm day}^{-1} [39,22]

     | Show Table
    DownLoad: CSV

    for x\in[0, 4] and \theta\in[-0.5, 0] .

    Firstly, we take \beta_1 = 1\times10^{-5} . Then R_0 = 0.8715 < 1 . By Theorem 4.2, the infection-free steady state P_0 = (1000, 0, 0, 0) is globally asymptotically stable (see Figure 1).

    Figure 1.  When R_0 < 1 , the infection-free steady state P_0 is globally asymptotically stable. Parameter values are given in the text.

    Next, we choose \beta_1 = 2.4\times10^{-5} . Then R_0 = 2.0588 > 1 and R_1 = 0.2989 < 1 . From Theorem 4.7, the immunity-free infected steady state P_1 = (897.8483, 0.9963, 86.6344, 0) is globally asymptotically stable (see Figure 2).

    Figure 2.  When R_0 = 2.0588 > 1 and R_1 = 0.2989 < 1 , the immunity-free infected steady state P_1 is globally asymptotically stable. See the text for the parameter values.

    Finally, with \beta_1 = 8.4\times10^{-5} , we get R_0 = 7.1474 > 1 and R_1 = 1.1594 > 1 . By Theorem 4.12, the infected-immune steady state P_2 = (609.8631, 3.3333,289.8550, 5.8964) is globally asymptotically stable (see Figure 3).

    Figure 3.  When R_1 > 1 , the infected-immune steady state P_2 is globally asymptotically stable. See the text for the parameter values.

    In this paper, we have proposed and studied a reaction-diffusion virus infection model by incorporating time delays, general incidence functions, and cell-to-cell transmission.

    We have proved that the global dynamics of system (1.4)–(1.6) is determined by the basic reproduction number of infection R_0 and the basic reproduction number of immunity R_1 . By analyzing the characteristic equations and constructing Lyapunov functionals, we have obtained the following conclusions: if R_0 < 1 , then the infected-free steady state P_0 is globally asymptotically stable; if R_1\leq1 < R_0 , then the immunity-free infected steady state P_1 is globally asymptotically stable under additional Assumption (A2); if R_1 > 1 , then the infected-immune steady state P_2 is globally asymptotically stable under additional Assumptions (A2) and (A3). We mention that most commonly used incidences satisfy (A1)–(A3). Some examples are the Holling type Ⅱ incidence f(u, v) = \frac{\beta u v}{1+\alpha v} [40], Beddington-DeAnglis incidence [41], and f(u, v) = ku\ln(1+\frac{\beta v}{k}) [56].

    Chen is supported by NSERC of Canada. Wang is supported by the NSFC (No. 11771374), the Nanhu Scholars Program for Young Scholars of Xinyang Normal University.

    The authors declare that there is no conflict of interest.



    [1] Gupta NM (2016) Factors affecting the efficiency of a water splitting photo catalyst: A perspective. Renewable Sustainable Energy Rev 71: 585-601.
    [2] Dubey PK, Tripathi P, Tiwari RS, et al. (2014) Synthesis of reduced graphene oxide-TiO2 nanoparticle composite systems and its application in hydrogen production. Int J Hydrogen Energy 39: 16282-16292. doi: 10.1016/j.ijhydene.2014.03.104
    [3] Dincer I (2002) Technical, environmental and exergetic aspects of hydrogen energy systems. Int J Hydrogen Energy 27: 265-285. doi: 10.1016/S0360-3199(01)00119-7
    [4] Dincer I (2012) Green methods for hydrogen production. Int J Hydrogen Energy 37: 1954-71. doi: 10.1016/j.ijhydene.2011.03.173
    [5] Barreto L, Makihira A, Riahi K (2003) The hydrogen economy in the 21st century: A sustainable development scenario. Int J Hydrogen Energy 28: 267-284. doi: 10.1016/S0360-3199(02)00074-5
    [6] Midilli A, Dincer I (2007) Key strategies of hydrogen energy systems for sustainability. Int J Hydrogen Energy 32: 511-524. doi: 10.1016/j.ijhydene.2006.06.050
    [7] Abe JO, Popoola IPA, Ajenifuja E, et al. (2019) Hydrogen energy, economy and storage: Review and recommendation. Int J Hydrogen Energy 44: 15072-86. doi: 10.1016/j.ijhydene.2019.04.068
    [8] Hoffman P (2019) The Forever Fuel, The Story of Hydrogen, 1st edition, Taylor & Francis group, eBook.
    [9] Holladay DJ, Hu J, King LD, et al. (2009) An overview of hydrogen production technologies. Catal Today 139: 244-260. doi: 10.1016/j.cattod.2008.08.039
    [10] El-Shafie M, Kambara S, Hayakawa Y (2019) Hydrogen production technologies overview. J Power Energy Eng 7: 107-154. doi: 10.4236/jpee.2019.71007
    [11] Simpson AP, Lutz AE (2007) Exergy analysis of hydrogen production via steam methane reforming. Int J Hydrogen Energy 32: 4811-20. doi: 10.1016/j.ijhydene.2007.08.025
    [12] Stefanidis GD, Vlachos DG (2010) Intensification of steam reforming of natural gas: choosing combustible fuel and reforming catalyst. Chem Eng Sci 65: 398-404. doi: 10.1016/j.ces.2009.06.007
    [13] Šingliar M (2007) Solar energy using for hydrogen production. Pet Coal 49: 40-47.
    [14] Guoxin H, Hao H, Yanhong L (2009) Hydrogen-rich gas production from pyrolysis of biomass in an auto generated steam atmosphere. Energy Fuels 23: 1748-53. doi: 10.1021/ef800988r
    [15] Bockris JOM, Dandapani B, Cocke D, et al. (1985) On the splitting of water. Int J Hydrogen Energy 10: 179-201. doi: 10.1016/0360-3199(85)90025-4
    [16] El-Shafie M, Kambara S, Hayakawa Y, et al. (2019) Preliminary results of hydrogen production fromwater vapor decomposition using DBD plasma in a PMCR reactor. Int J Hydrogen Energy 44: 20239-48. doi: 10.1016/j.ijhydene.2019.05.199
    [17] Givotov V, Fridman A, Krotov M, et al. (1981) Plasmochemical methods of hydrogen production. Int J Hydrogen Energy 6: 441-9.
    [18] Fridman AA (2008) Plasma chemistry. Cambridge University Press.
    [19] Maehara T, Toyota H, Kuramoto M, et al. (2006) Radio frequency plasma in water. Jpn J Appl Phys 45: 8864-8. doi: 10.1143/JJAP.45.8864
    [20] Nguyen SVT, Foster JE, Gallimore AD (2009) Operating a radiofrequency plasma source on water vapor. Rev Sci Instrum 80:083503. doi: 10.1063/1.3202250
    [21] Burlica R, Shih KY, Locke B (2010) Formation of H2 and H2O2 in a water-spray gliding arc non thermal plasma reactor. Ind Eng Chem Res.
    [22] Mutaf-Yardimci O, Saveliev A, Fridman A, et al. (1998) Employing plasma as catalyst in hydrogen production. Int J Hydrogen Energy 23: 1109-11. doi: 10.1016/S0360-3199(98)00005-6
    [23] Xi Zhen L, Liu CJ, Eliasson B (2003) Hydrogen production from methanol using corona discharges. Chin Chem Lett 14: 631-3.
    [24] Sarmiento B, Brey JJ, Viera IG, et al. (2007) Hydrogen production by reforming of hydrocarbons and alcohols in a dielectric barrier discharge. J Power Sources 169: 140-3. doi: 10.1016/j.jpowsour.2007.01.059
    [25] Rico VJ, Cotrino J, Gallardo V, et al. (2009) Hybrid catalytic-DBD plasma reactor for the production of hydrogen and preferential CO-oxidation (CO-PROX) at reduced temperatures. Chem Commun 41: 6192-4.
    [26] Wang W, Zhu C, Cao Y (2010) DFT study on pathways of steam reforming of ethanol under cold plasma conditions for hydrogen generation. Int J Hydrogen Energy 35: 1951-6. doi: 10.1016/j.ijhydene.2009.12.170
    [27] Yan Z, Chen L, Wang H (2008) Hydrogen generation by glow discharge plasma electrolysis of ethanol solutions. J Phys D: Appl Phys 41: 155205. doi: 10.1088/0022-3727/41/15/155205
    [28] Rehman F, Lozano-Parada JH, Zimmerman WB (2012) A kinetic model for H2 production by plasmolysis of water vapours at atmospheric pressure in a dielectric barrier discharge micro channel reactor. Int J Hydrogen Energy 37: 17678-90. doi: 10.1016/j.ijhydene.2012.08.113
    [29] Jasinski M, Dors M, Mizeraczyk J (2008) Production of hydrogen via methane reforming using atmospheric pressure microwave plasma. J Power Sour 181: 41-5. doi: 10.1016/j.jpowsour.2007.10.058
    [30] Sobacchi M, Saveliev A, Fridman A, et al. (2002) Experimental assessment of a combined plasma/catalytic system for hydrogen production via partial oxidation of hydrocarbon fuels. Int J Hydrogen Energy 27: 635-42. doi: 10.1016/S0360-3199(01)00179-3
    [31] Avtaeva S, General A, Kel'man V (2010) Kinetic model for low density non-stationary gas discharge in water vapour. J Phys D: Appl Phys 43: 315201. doi: 10.1088/0022-3727/43/31/315201
    [32] Dolan T (1993) Electron and ion collisions with water vapour. J Phys D: Appl Phys 26: 4. doi: 10.1088/0022-3727/26/1/002
    [33] Medodovic S, Locke B (2009) Primary chemical reactions in pulsed electrical discharge channels in water. J Phys D: Appl Phys 42: 049801. doi: 10.1088/0022-3727/42/4/049801
    [34] Lukes P, Clupek M, Babicky V, et al. (2008) Role of solution conductivity in the electron impact dissociation of H2O induced by plasma processes in the pulsed corona discharge in water. HAKONE XI, 11th international symposium on high pressure, low temperature plasma chemistry, contributed papers, Olexron Island.
    [35] Shirafuji T, Morita T, Sakai O, et al. (2009) Enhancement of OH production rate in plasma on water by mixing Ar. Proceedings of the international symposium on plasma on plasma chemistry, Bochum, Germany.
    [36] Locke B, Sato M, Sunka P, et al. (2006) Electrohydraulic discharge and nonthermal plasma for water treatment. Ind Eng Chem Res 45: 882-905. doi: 10.1021/ie050981u
    [37] Zimmerman WBJ (2006) Multiphysics modeling with finite element methods, series on stability, vibration and control of systems, series A. London: World Scientific Publishing Company.
    [38] Lozano-Parada JH, Zimmerman WB (2010) The role of kinetics in the design of plasma microreactors. Chem Eng Sci 65: 4925-30. doi: 10.1016/j.ces.2010.03.056
    [39] Zhang Y, Wen XH, Yang WH (2007) Excitation temperatures of atmospheric argon in dielectric barrier discharges. Plasma Sources Sci Technol 16: 441. doi: 10.1088/0963-0252/16/3/003
    [40] Nehra V, Kumar A, Dwivedi H (2008) Atmospheric non-thermal plasma sources. Int J Eng 2: 53.
    [41] Abdul-Majeed WS, Parada JHL, Zimmerman WB (2011) Optimization of a miniaturized DBD plasma chip for mercury detection in water samples. Anal Bioanal Chem 401: 2713-22. doi: 10.1007/s00216-011-5118-9
    [42] Zito JC, Arnold DP, Durscher RJ, et al. (2010) Investigation of impedance characteristics and power delivery for dielectric barrier discharge plasma actuators. In: 48th AIAA Aerospace Sciences Meeting.
    [43] Kostov K, Honda R, Alves L, et al. (2009) Characteristics of dielectric barrier discharge reactor for material treatment. Braz J Phys 39: 322-5.
    [44] Rehman F, Liu Y, Zimmerman WBJ (2016) The role of chemical kinetics in using O3 generation as proxy for hydrogen production from water vapour plasmolysis. Int J Hydrogen Energy 41: 6180-92. doi: 10.1016/j.ijhydene.2016.02.071
    [45] Shih KY, Locke BR (2010) Optical and electrical diagnostics of the effects of conductivity on liquid phase electrical discharge. IEEE Trans Plasma Sci 99: 1-10.
    [46] El-Shafie M, Kambara S, Hayakawa Y (2020) One-dimensional simulation of hydrogen production kinetic models by water vapor plasmolysis in a DBD plate reactor. J Theor Appl Phys 14: 181-194. doi: 10.1007/s40094-020-00376-3
  • This article has been cited by:

    1. A. M. Elaiw, N. H. AlShamrani, Analysis of an HTLV/HIV dual infection model with diffusion, 2021, 18, 1551-0018, 9430, 10.3934/mbe.2021464
    2. Qing Ge, Xia Wang, Libin Rong, A delayed reaction–diffusion viral infection model with nonlinear incidences and cell-to-cell transmission, 2021, 14, 1793-5245, 10.1142/S179352452150100X
    3. Baolin Li, Fengqin Zhang, Xia Wang, A delayed diffusive HBV model with nonlinear incidence and CTL immune response, 2022, 45, 0170-4214, 11930, 10.1002/mma.8547
    4. Zakaria Yaagoub, Karam Allali, Fractional HBV infection model with both cell-to-cell and virus-to-cell transmissions and adaptive immunity, 2022, 165, 09600779, 112855, 10.1016/j.chaos.2022.112855
    5. Weixin Wu, Zengyun Hu, Long Zhang, Zhidong Teng, Traveling waves for a diffusive virus infection model with humoral immunity, cell‐to‐cell transmission, and nonlinear incidence, 2023, 0170-4214, 10.1002/mma.9291
    6. Yuequn Gao, Ning Li, Fractional order PD control of the Hopf bifurcation of HBV viral systems with multiple time delays, 2023, 83, 11100168, 1, 10.1016/j.aej.2023.10.032
  • Reader Comments
  • © 2020 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(5135) PDF downloads(156) Cited by(3)

Figures and Tables

Figures(18)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog