Review Special Issues

Reversing the nutrient drain through urban insect farming—opportunities and challenges

  • Received: 05 September 2018 Accepted: 19 November 2018 Published: 28 November 2018
  • Cities consume the majority of proteins produced globally but have unsustainable, linear food systems from production to consumption to disposal, resulting in significant nutrient losses. The industrial rearing of insects is a promising strategy for converting otherwise lost nutrients back into protein-rich animal feed and fertilizer, particularly to supplement local food production. The black soldier fly (BSF), Hermetia illucens, has been identified as a candidate for industrial rearing. BSF has a superior feed conversion ratio and cycle-time compared to other edible insects and can convert and recover nutrients from a vast variety of organic materials to protein, oil and chitin making it an attractive solution for the management of urban organic solid waste. With an increasing awareness of the environmental urgency and interest in the economic potential of the technology, this review discusses the technological factors confounding the upscaling of insect farming in urban and peri-urban contexts using BSF as a case study. These include the challenges of feed homogenisation and pre-treatment, of integrating insect life-cycle factors (e.g. mating) with bioprocess engineering concepts (which complicates automation), of meeting the nutritional requirements of the larvae at different stages of growth in order to maximize bioconversion and product quality, and of elucidating the impact of microbiome on complex behaviours and bioconversion. A multidisciplinary effort is therefore required to lead urban insect farming to full development to ultimately contribute to future food security.

    Citation: Yingyu Law, Leo Wein. Reversing the nutrient drain through urban insect farming—opportunities and challenges[J]. AIMS Bioengineering, 2018, 5(4): 226-237. doi: 10.3934/bioeng.2018.4.226

    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
  • Cities consume the majority of proteins produced globally but have unsustainable, linear food systems from production to consumption to disposal, resulting in significant nutrient losses. The industrial rearing of insects is a promising strategy for converting otherwise lost nutrients back into protein-rich animal feed and fertilizer, particularly to supplement local food production. The black soldier fly (BSF), Hermetia illucens, has been identified as a candidate for industrial rearing. BSF has a superior feed conversion ratio and cycle-time compared to other edible insects and can convert and recover nutrients from a vast variety of organic materials to protein, oil and chitin making it an attractive solution for the management of urban organic solid waste. With an increasing awareness of the environmental urgency and interest in the economic potential of the technology, this review discusses the technological factors confounding the upscaling of insect farming in urban and peri-urban contexts using BSF as a case study. These include the challenges of feed homogenisation and pre-treatment, of integrating insect life-cycle factors (e.g. mating) with bioprocess engineering concepts (which complicates automation), of meeting the nutritional requirements of the larvae at different stages of growth in order to maximize bioconversion and product quality, and of elucidating the impact of microbiome on complex behaviours and bioconversion. A multidisciplinary effort is therefore required to lead urban insect farming to full development to ultimately contribute to future food security.


    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 $ \hat{\tau} $ 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 $ \Delta $ 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 $ \frac{\beta uv}{1+\alpha 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 $. $ \Omega $ is a general open bounded domain in $ \mathbb{R}^n $ with smooth boundary $ \partial\Omega $. We consider model (1.4) with the homogeneous Neumann boundary conditions

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

    where $ \frac{\partial }{\partial \vec{n}} $ denotes the outward normal derivative on $ \partial\Omega $. 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 $ \phi_i $'s ($ i = 1 $, $ 2 $, $ 3 $, $ 4 $) are bounded and uniformly continuous functions on $ \overline{\Omega}\times[-\tau, 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.

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

    $ {\rm (i)} $ $ f(u, v)\ge 0 $ and $ g(u, w)\ge 0 $ for $ u\ge0 $, $ v\ge 0 $, and $ w\ge0 $, and the equalities hold if and only if $ uv = 0 $ and $ uw = 0 $;

    $ {\rm (ii)} $ There exists $ \eta_1 > 0 $ and $ \eta_2 > 0 $ such that $ f(u, v)\leq\eta_1u $ and $ g(u, w)\leq\eta_2u $ for $ u\ge0 $, $ v\ge 0 $, and $ w\geq 0 $;

    $ {\rm (iii)} $ $ \frac{\partial f(u, v)}{\partial u} $ and $ \frac{\partial g(u, w)}{\partial u} $ are continuous with $ \frac{\partial f(u, v)}{\partial u} > 0 $ and $ \frac{\partial g(u, w)}{\partial u} > 0 $ for $ u\ge0 $, $ v > 0 $, and $ w > 0 $;

    $ {\rm (iv)} $ $ \frac{\partial f(u, v)}{\partial v} $ and $ \frac{\partial g(u, w)}{\partial w} $ are continuous with $ \frac{\partial f(u, v)}{\partial v}\ge0 $ and $ \frac{\partial g(u, w)}{\partial w}\geq0 $ for $ u\ge0 $, $ v\ge0 $, and $ w\geq 0 $;

    $ {\rm (v)} $ $ v\frac{\partial f(u, v)}{\partial v}-f(u, v)\leq0 $ and $ w\frac{\partial g(u, w)}{\partial w}-g(u, w)\leq0 $ for $ u\ge0 $, $ v\ge 0 $, and $ w\geq 0 $.

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

    $ \frac {\partial (\frac{f(u, v)}{v})}{\partial v} = \frac {\frac{\partial f(u, v)}{\partial v}v-f(u, v)}{v^2} \lt 0, $

    which implies that $ \frac{f(u, v)}{v} $ is decreasing on $ (0, \infty) $. In particular, for any $ u > 0 $ and $ v > 0 $,

    $ \frac {f(u, v)}{v}\le \lim\limits_{v\to 0^+}\frac {f(u, v)}{v} = \frac {\partial f(u, 0)}{\partial 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 $ \bf{X}: = \mathcal{C}(\overline{\Omega}, \mathbb{R}^4) $ be the Banach space equipped with the supremum norm $ \| \cdot \|_{\bf{X}} $. For $ \tau\geq0 $, define $ \mathcal{C} = \mathcal{C}([-\tau, 0], \bf{X}) $, which is a Banach space equipped with the norm $ \| \phi \| = \max_{\theta\in[-\tau, 0]}\| \phi(\theta) \|_{\bf{X}} $. If $ \sigma > 0 $ and $ U: [-\tau, \sigma)\to \bf{X} $, then for $ t\in [0, \sigma) $, $ U_t\in \mathcal{C} $ is defined by $ U_t(\theta) = U(t+\theta) $ for $ \theta\in [-\tau, 0] $. Denote $ \bf{X}^+ = \mathcal{C}(\overline{\Omega}, \mathbb{R}^4_+) $ and $ \mathcal{C^+} = \mathcal{C}([-\tau, 0], \bf{X}^+) $. Then both $ (\bf{X}, \bf{X}^+) $ and $ (\mathcal{C}, \mathcal{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 $ \phi = (\phi_1, \phi_2, \phi_3, \phi_4)\in\mathcal{C}^+ $, system (1.4)–(1.6) has a unique solution $ U(\cdot, t, \phi) = (u(\cdot, t, \phi), w(\cdot, t, \phi), v(\cdot, t, \phi), z(\cdot, t, \phi)) $ on $ [0, \infty) $ with $ U_0(\cdot, \phi) = \phi $. Moreover, $ U_t(\cdot, \phi)\in\mathcal{C}^+ $ for $ t\ge 0 $ and $ U(\cdot, t, \phi) $ is a classical solution.

    Let $ \Phi(t): \mathcal{C}^+\to \mathcal{C}^+ $ be the solution semiflow associated with (1.4)–(1.6), that is, $ \Phi(t, \phi) = U_t(\cdot, \phi) $, where $ U(\cdot, t, \phi) $ is the solution of (1.4)–(1.6) with the initial condition $ \phi\in \mathcal{C}^+ $.

    The following result gives some properties of solutions.

    Lemma 2.2. For $ \phi\in \mathcal{C}^+ $, the following statements hold for the solution $ U(\cdot, t, \phi) $ of (1.4)–(1.6).

    $ {\rm(i)} $ $ \limsup\limits_{t\to\infty}u(x, t, \phi)\le \frac{s}{d} $, $ \limsup\limits_{t\to\infty}w(x, t, \phi)\le \frac{e^{-m\tau}(\eta_1+\eta_2)s}{d\delta} $, $ \limsup\limits_{t\to\infty}v(x, t, \phi)\le \frac{Ne^{-m\tau}(\eta_1+\eta_2)s}{dc} $, and $ \limsup\limits_{t\to\infty}z(x, t, \phi)\le \frac {e^{-m\tau}(\eta_1+\eta_2)s}{d\min\{\delta, b\}} $ uniformly for all $ x\in \Omega $.

    $ {\rm (ii)} $ $ u(\cdot, t, \phi) > 0 $ for $ t > 0 $ and $ \liminf\limits_{t\to\infty}u(x, t, \phi)\geq\frac{s}{d+\eta_1+\eta_2} $ uniformly for all $ x\in \Omega $.

    $ {\rm(iii)} $ If $ w(\cdot, t_0, \phi)\not\equiv0 $ 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

    $ {v(x,t)tD1Δv(x,t)cv(x,t),t>t0, xΩ,v(x,t)n=0,t>t0, xΩ.
    $

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

    $ {ˉv(x,t)t=D1Δˉv(x,t)cˉv(x,t),t>t0, xΩ,ˉv(x,t)n=0,t>t0, xΩ,ˉv(x,t0)=v(x,t0), x¯Ω.
    $

    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

    $ R0=Nemτcf(sd,0)v+emτδg(sd,0)w.
    $

    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

    $ sduf(u,v)g(u,w)=0,
    $
    (3.1a)
    $ f(u,v)emτ+g(u,w)emτδwpwz=0,
    $
    (3.1b)
    $ Nδwcv=0,
    $
    (3.1c)
    $ qwzbz=0.
    $
    (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

    $ H1(w)=f(sδwemτd,Nδwc)+g(sδwemτd,w)δwemτ.
    $
    (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),

    $ H1(w1)=δemτdf(u1,v1)u+Nδcf(u1,v1)vδemτdg(u1,w1)u+ g(u1,w1)wNδcv1f(u1,v1)1w1g(u1,w1)=δemτdf(u1,v1)uδemτdg(u1,w1)u+ Nδcv1(v1f(u1,v1)vf(u1,v1))+ 1w1(w1g(u1,w1)wg(u1,w1))<0.
    $

    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

    $ H2(u)=f(u,Nδbcq)+g(u,bq)s+du.
    $
    (3.3)

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

    $ H2(u)=f(u,Nδbcq)u+g(u,bq)u+d>0.
    $

    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

    $ Qt=LΔQ+AQ+BQτ,
    $
    (4.1)

    where

    $ L=diag(0,0,D1,D2),Q=(u(x,t),w(x,t),v(x,t),z(x,t)),Qτ=(uτ,wτ,vτ,zτ)=(u(x,tτ),w(x,tτ),v(x,tτ),z(x,tτ)),A=(df(u,v)ug(u,w)ug(u,w)wf(u,v)v00δpz0pw0Nδc00qz0qwb),B=(0000(f(u,v)u+g(u,w)u)emτg(u,w)wemτf(u,v)vemτ000000000).
    $

    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

    $ |λE+μiLABeλτ|=0
    $
    (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

    $ (λ+c+μiD1)(λ+δg(sd,0)wemτeλτ)Nδf(sd,0)vemτeλτ=0.
    $
    (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

    $ 1=|1λ0+δg(sd,0)wemτeλ0τ+Nδ(λ0+δ)(λ0+c+μi0D1)f(sd,0)vemτeλ0τ||1λ0+δg(sd,0)wemτeλ0τ|+|Nδ(λ0+δ)(λ0+c+μi0D1)f(sd,0)vemτeλ0τ|1δg(sd,0)wemτ+Nc+μi0D1f(sd,0)vemτ1δg(sd,0)wemτ+Ncf(sd,0)vemτ=R0,
    $

    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

    $ W(t)=Ω(emτw(x,t)+1Nemτv(x,t)+pqemτz(x,t))dx+Ω(ttτf(u(x,θ),v(x,θ))dθ+tτg(u(x,θ),w(x,θ))dθ)dx.
    $

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

    $ dW(t)dt=Ω(f(u(x,t),v(x,t))cNemτv(x,t)pbqemτz(x,t))dx+1NemτΩD1Δv(x,t)dx.
    $

    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

    $ f(u(x,t),v(x,t))cNemτv(x,t)=cNemτv(x,t)(Ncemτf(u(x,t),v(x,t))v(x,t)1)cNemτv(x,t)(Ncemτf(sd,v(x,t))v(x,t)1)cNemτv(x,t)(Ncemτf(sd,0)v1)cNemτv(x,t)(R01).
    $

    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,

    $ dW(t)dtΩ(cNemτv(x,t)(R01)pbqemτz(x,t))dx0.
    $

    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

    $ f1=f(u1,v1),g1=g(u1,w1),f1u=f(u1,v1)u,f1v=f(u1,v1)v,g1u=g(u1,w1)u,g1w=g(u1,w1)w.
    $

    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

    $ 1=λ+dλ+d+f1u+g1u(g1wλ+δe(λ+m)τ+Nδf1v(λ+δ)(λ+c+μiD1)e(λ+m)τ).
    $
    (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

    $ 1=|λ1+dλ1+d+f1u+g1u(g1wλ1+δemτeλ1τ+Nδf1v(λ1+δ)(λ1+c+μi1D1)emτeλ1τ)|<|g1wλ1+δemτeλ1τ|+|Nδf1v(λ1+δ)(λ1+c+μi1D1)emτeλ1τ|<g1wδemτ+Nf1vcemτ.
    $
    (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) $,

    $ {w(x,t)t=emτf(u0,0)vv(x,tτ)+emτg(u0,0)ww(x,tτ)δw(x,t),v(x,t)t=D1Δv(x,t)+Nδw(x,t)cv(x,t).
    $
    (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\{ w(x,t)t=emτ(f(u0,0)vε0)v(x,tτ)+emτ(g(u0,0)wε0)w(x,tτ)(δ+pε0)w(x,t),v(x,t)t=D1Δv(x,t)+Nδw(x,t)cv(x,t)
    \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

    $ {w(x,t)temτ(f(u0,0)vε0)v(x,tτ)+emτ(g(u0,0)wε0)w(x,tτ)(δ+pε0)w(x,t),tt2, xΩ,v(x,t)t=D1Δv(x,t)+Nδw(x,t)cv(x,t),tt2, xΩ,v(x,t)n=0,tt2,xΩ.
    $
    (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 $,

    $ {vv1u1f(u,v)uf(u1,v1)<1if  0<v<v1,1u1f(u,v)uf(u1,v1)<vv1if  v1<v.
    $

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

    $ {ww1u1g(u,w)ug(u1,w1)<1if  0<w<w1,1u1g(u,w)ug(u1,w1)<ww1if  w1<w.
    $

    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

    $ L(x,t)=u1G(u(x,t)u1)+emτw1G(w(x,t)w1)+f(u1,v1)cv1v1G(v(x,t)v1)+pqemτz(x,t)+f(u1,v1)ttτG(f(u(x,θ),v(x,θ))f(u1,v1))dθ+g(u1,w1)ttτG(g(u(x,θ),w(x,θ))g(u1,w1))dθ.
    $

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

    $ L(x,t)t=(1u1u(x,t))[sdu(x,t)f(u,v)g(u,w)]+pqemτ[qw(x,t)z(x,t)bz(x,t)]+emτ(1w1w(x,t))[emτf(uτ,vτ)+emτg(uτ,wτ)δw(x,t)pw(x,t)z(x,t)]+f(u1,v1)cv1(1v1v(x,t))[D1Δv(x,t)+Nδw(x,t)cv(x,t)]+pqemτD2Δz(x,t)+f(u,v)f(uτ,vτ)+f(u1,v1)lnf(uτ,vτ)f(u,v)+g(u,w)g(uτ,wτ)+g(u1,w1)lng(uτ,wτ)g(u,w).
    $

    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

    $ L(x,t)t=du1(2u1u(x,t)u(x,t)u1)+pqemτ(R11)z(x,t)+f(u1,v1)cv1(1v1v(x,t))D1Δv(x,t)+pqemτD2Δz(x,t)f(u1,v1)w(x,t)v1v(x,t)w1+f(u1,v1)[3u1u(x,t)f(uτ,vτ)w1f(u1,v1)w(x,t)+f(u,v)u1f(u1,v1)u(x,t)v(x,t)v1+lnf(uτ,vτ)f(u,v)]+g(u1,w1)[2u1u(x,t)g(uτ,wτ)w1g(u1,w1)w(x,t)+g(u,w)u1g(u1,w1)u(x,t)w(x,t)w1+lng(uτ,wτ)g(u,w)]=du1(2u1u(x,t)u(x,t)u1)+pqemτ(R11)z(x,t)+f(u1,v1)cv1(1v1v(x,t))D1Δv(x,t)+pqemτD2Δz(x,t)f(u1,v1)G(w(x,t)v1v(x,t)w1)+f(u1,v1)[G(f(u,v)u1f(u1,v1)u(x,t))G(v(x,t)v1)G(u1u(x,t))G(f(uτ,vτ)w1f(u1,v1)w(x,t))]+g(u1,w1)[G(g(u,w)u1g(u1,w1)u(x,t))G(w(x,t)w1)G(u1u(x,t))G(g(uτ,wτ)w1g(u1,w1)w(x,t))].
    $

    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

    $ 0=Ω1v(x,t)v(x,t)  n dx=Ω(1v(x,t)v(x,t))dx=Ω(1v(x,t)Δv(x,t)1v2(x,t)v(x,t)2)dx.
    $

    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

    $ dL(t)dtf(u1,v1)Ω[G(f(u,v)u1f(u1,v1)u(x,t))G(v(x,t)v1)]dx+g(u1,w1)Ω[G(g(u,w)u1g(u1,w1)u(x,t))G(w(x,t)w1)]dxg(u1,w1)Ω[G(u1u(x,t))+G(g(uτ,wτ)w1g(u1,w1)w(x,t))]dxf(u1,v1)Ω[G(u1u(x,t))+G(f(uτ,vτ)w1f(u1,v1)w(x,t))+G(w(x,t)v1v(x,t)w1)]dx.
    $

    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

    $ dL(t)dtg(u1,w1)Ω[G(u1u(x,t))+G(g(uτ,wτ)w1g(u1,w1)w(x,t))]dxf(u1,v1)Ω[G(u1u(x,t))+G(f(uτ,vτ)w1f(u1,v1)w(x,t))+G(w(x,t)v1v(x,t)w1)]dx0.
    $

    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

    $ f2=f(u2,v2),g2=g(u2,w2),f2u=f(u2,v2)u,f2v=f(u2,v2)v,g2u=g(u2,w2)u,g2w=g(u2,w2)w.
    $

    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

    $ 0=(λ+d+f2u+g2u)(λ+c+μiD1)(λ+bqw2+μiD2)(λ+δ+pz2)+qw2pz2(λ+c+μiD1)(λ+d+f2u+g2u)(λ+d)(λ+c+μiD1)(λ+bqw2+μiD2)g2we(λ+m)τ(λ+d)(λ+bqw2+μiD2)Nδf2ve(λ+m)τ.
    $
    (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

    $ 1=(λ2+μi2D2)(λ2+d)(λ2+d+f2u+g2u)[(λ2+μi2D2)(λ2+δ+pz2)+pbz2](g2we(λ2+m)τ+Nδf2ve(λ2+m)τλ2+c+μi2D1),
    $

    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) $,

    $ {w(x,t)t=emτf(u1,v1)vv(x,tτ)+emτg(u1,w1)ww(x,tτ)(δ+pz1)w(x,t)pw1z(x,t),v(x,t)t=D1Δv(x,t)+Nδw(x,t)cv(x,t),z(x,t)t=D2Δz(x,t)+qz1w(x,t)+(qw1b)z(x,t).
    $
    (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\{ vv2u2f(u,v)uf(u2,v2)<1if  0<v<v2,1u2f(u,v)uf(u2,v2)<vv2if  v2<v.
    \right. $

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

    $ \left\{ ww2u2g(u,w)ug(u2,w2)<1if  0<w<w2,1u2g(u,w)ug(u2,w2)<ww2if  w2<w.
    \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 | there exists t5R+ such that either w(,t5,ϕ)0 or v(,t5,ϕ)0there exists t6R+ such that z(,t6,ϕ)0
    \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

    $ I(x,t)=u2G(u(x,t)u2)+emτw2G(w(x,t)w2)+f(u2,v2)cv2v2G(v(x,t)v2)+pqemτz2G(z(x,t)z2)+f(u2,v2)ttτG(f(u(x,θ),v(x,θ))f(u2,v2))dθ+g(u2,w2)ttτG(g(u(x,θ),w(x,θ))g(u2,w2))dθ.
    $

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

    $ I(x,t)t=(1u2u(x,t))[sdu(x,t)f(u,v)g(u,w)]+pqemτ(1z2z(x,t))[qw(x,t)z(x,t)bz(x,t)]+emτ(1w2w(x,t))[emτf(uτ,vτ)+emτg(uτ,wτ)δw(x,t)pw(x,t)z(x,t)]+f(u2,v2)cv2(1v2v(x,t))[D1Δv(x,t)+Nδw(x,t)cv(x,t)]+pqemτ(1z2z(x,t))D2Δz(x,t)+f(u,v)f(uτ,vτ)+f(u2,v2)lnf(uτ,vτ)f(u,v)+g(u,w)g(uτ,wτ)+g(u2,w2)lng(uτ,wτ)g(u,w).
    $

    With the following relations,

    $ s=du2+f(u2,v2)+g(u2,w2),δ=emτf(u2,v2)+emτg(u2,w2)pw2z2w2,c=Nδw2v2,w2=bq,
    $

    we get

    $ I(x,t)t=du2(2u2u(x,t)u(x,t)u2)+f(u2,v2)cv2(1v2v(x,t))D1Δv(x,t)+pqemτ(1z2z(x,t))D2Δz(x,t)f(u2,v2)w(x,t)v2v(x,t)w2+f(u2,v2)[3u2u(x,t)f(uτ,vτ)w2f(u2,v2)w(x,t)+f(u,v)u2f(u2,v2)u(x,t)v(x,t)v2+lnf(uτ,vτ)f(u,v)]+g(u2,w2)[2u2u(x,t)g(uτ,wτ)w2g(u2,w2)w(x,t)+g(u,w)u2g(u2,w2)u(x,t)w(x,t)w2+lng(uτ,wτ)g(u,w)].
    $

    Then

    $ dI(t)dt=Ωdu2(2u2u(x,t)u(x,t)u2)dx+Ωf(u2,v2)cv2(1v2v(x,t))D1Δv(x,t)dx+Ωpqemτ(1z2z(x,t))D2Δz(x,t)dxf(u2,v2)ΩG(w(x,t)v2v(x,t)w2)dx+f(u2,v2)Ω[G(f(u,v)u2f(u2,v2)u(x,t))G(v(x,t)v2)G(u2u(x,t))G(f(uτ,vτ)w2f(u2,v2)w(x,t))]dx+g(u2,w2)Ω[G(g(u,w)u2g(u2,w2)u(x,t))G(w(x,t)w2)G(u2u(x,t))G(g(uτ,wτ)w2g(u2,w2)w(x,t))]dx.
    $

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

    $ Ωdu2(2u(x,t)u2u2u(x,t))dx0,Ωf(u2,v2)cv2(1v2v(x,t))D1Δv(x,t)dx0,Ωpqemτ(1z2z(x,t))D2Δz(x,t)dx0,G(f(u,v)u2f(u2,v2)u(x,t))G(v(x,t)v2),G(g(u,w)u2g(u2,w2)u(x,t))G(w(x,t)w2).
    $

    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

    $ {u(x,t)t=sdu(x,t)β1u(x,t)v(x,t)1+α1v(x,t)β2u(x,t)w(x,t)1+α2w(x,t), xΩ, t>0,w(x,t)t=emτ(β1u(x,tτ)v(x,tτ)1+α1v(x,tτ)+β2u(x,tτ)w(x,tτ)1+α2w(x,tτ))δw(x,t)pw(x,t)z(x,t), xΩ, t>0,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,
    $
    (5.1)

    subject to the homogeneous Neumann boundary conditions

    $ vn=0, zn=0, xΩ, t>0.
    $

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

    $ R0=Nβ1scdemτ+β2sδdemτ
    $

    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

    $ B1=Nδ2emτ(β1α2+β2α1+α1α2d)>0,B2=δdemτ(Nδα1+cα2)Nδβ1α2s+Nδ2β1emτNδβ2α1s+δemτβ2c,B3=dcδemτ(R01)<0 since R0>1.
    $

    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

    $ u(x,θ)=23+0.2cosπx2,w(x,θ)=0.7+0.2cosπx2,v(x,θ)=3+0.2cosπx2,z(x,θ)=2+0.2cosπx2
    $
    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] Alexandratos N, Bruinsma J (2012) World agriculture towards 2030/2050: The 2012 revision. ESA Working Paper No. 12-03, Food and Agriculture Organization of the United Nations (FAO) , Rome.
    [2] Godfray HCJ, Beddington JR, Crute IR, et al. (2010) Food security: The challenge of feeding 9 billion people. Science 327: 812–818. doi: 10.1126/science.1185383
    [3] Oliva-Teles A, Enes P, Peres H (2015) Replacing fishmeal and fish oil in industrial aquafeeds for carnivorous fish, In: Davis DA, Feed and Feeding Practices in Aquaculture, 8 Eds., Oxford: Woodhead Publishing, 203–233.
    [4] Skøt J, Lipper L, Thomas G, et al. (2016) The state of food and agriculture: Climate change, agriculture and food security. Food and Agriculture Organization of the United Nations (FAO), Rome.
    [5] Food and Agriculture Organization of the United Nations (FAO) (2013) Food wastage footprint: Impact on natural resources. Summary report. Rome.
    [6] Tilman D, Fargione J, Wolff B, et al. (2001) Forecasting agriculturally driven global environmental change. Science 292: 281–284. doi: 10.1126/science.1057544
    [7] Zhao C, Liu B, Piao S, et al. (2017) Temperature increase reduces global yields of major crops in four independent estimates. P Natl Acad Sci USA 114: 9326–9331. doi: 10.1073/pnas.1701762114
    [8] Adhikari B, Barrington S, Martinez J (2009) Urban food waste generation: Challenges and opportunities. Int J Environ Waste Manage 3: 4–21. doi: 10.1504/IJEWM.2009.024696
    [9] Hoornweg D, Bhada-Tata P, Kennedy C (2013) Environment: Waste production must peak this century. Nature 502: 615–617. doi: 10.1038/502615a
    [10] Gustavsson J, Cederberg C, Sonesson U, et al. (2011) Global fodd losses and food waste. extent, causes and prevention. Food and Agriculture Organization of the United Nations (FAO), Rome, Italy.
    [11] Diener S, Zurbrügg C, Tockner K (2009) Conversion of organic material by black soldier fly larvae: Establishing optimal feeding rates. Waste Manag Res 27: 603–610. doi: 10.1177/0734242X09103838
    [12] St-Hilaire S, Sheppard C, Tomberlin JK, et al. (2007) Fly Prepupae as a Feedstuff for Rainbow Trout, Oncorhynchus mykiss. J World Aquacul Soc 38: 59–67. doi: 10.1111/j.1749-7345.2006.00073.x
    [13] Newton GL, Booram CV, Barker RW, et al. (1977) Dried Hermetia illucens larvae meal as a supplement for swine. J Anim Sci 44: 395–400. doi: 10.2527/jas1977.443395x
    [14] Lee J, Kim YM, Park YK, et al. (2018) Black soldier fly (Hermetia illucens) larvae enhances immune activities and increases survivability of broiler chicks against experimental infection of Salmonella Gallinarum. J Vet Med Sci 80: 736–740. doi: 10.1292/jvms.17-0236
    [15] Kupferschmidt K (2015) Feature: Why insects could be the ideal animal feed. Science.
    [16] Alltech (2017) 7th Annual Alltech Global Feed Survey 2018. Available from: https://go.alltech.com/alltech-feed-survey.
    [17] Gilbert R (2004) The world animal feed industry. Protein Sources for the Animal Feed Industry. Food and Agriculture Organization of the United Nations (FAO). Rome, Italy.
    [18] SEAFISH (2016) Fishmeal and fish oil facts and figures. Available from: http://www.seafish.org/media/Publications/SeafishFishmealandFishOilFactsandFigures_201612.pdf.
    [19] Food and Agriculture Organization of the United Nations (FAO) (2018)The State of World Fisheries and Aquaculture: Meeting the sustainable development goals. Rome.
    [20] Masuda T, Goldsmith PD (2009) World soybean production: Area harvested, yield, and long-term projections. Int Food Agribus Man 12: 143–163.
    [21] Department of Agriculture (US) (2018) Fish Meal Production by Country in 1000 MT. Available from: https://www.indexmundi.com/agriculture/?commodity=fish-meal: Index Mundi.
    [22] Tacon AGJ, Metian M (2008) Global overview on the use of fish meal and fish oil in industrially compounded aquafeeds: Trends and future prospects. Aquaculture 285: 146–158. doi: 10.1016/j.aquaculture.2008.08.015
    [23] Trostle R (2010) Global Agricultural Supply and Demand: Factors Contributing to the Recent Increase in Food Commodity Prices. Available from: https://www.hsdl.org/?abstract&did=485652.
    [24] Verkerk MC, Tramper J, van Trijp JCM, et al. (2007) Insect cells for human food. Biotechnol Adv 25: 198–202. doi: 10.1016/j.biotechadv.2006.11.004
    [25] Dobermann D, Swift JA, Field LM (2017) Opportunities and hurdles of edible insects for food and feed. Nutr Bull 42: 293–308. doi: 10.1111/nbu.12291
    [26] van Huis A (2013) Potential of insects as food and feed in assuring food security. Annu Rev Entomol 58: 563–583. doi: 10.1146/annurev-ento-120811-153704
    [27] Hilkens W, de Klerk B (2017) The cultivation of insects: A small sector with great opportunities. Netherlands: ABN AMRO and BOM.
    [28] Tomberlin JK, Sheppard DC, Joyce JA (2002) Selected Life-History Traits of Black Soldier Flies (Diptera: Stratiomyidae) Reared on Three Artificial Diets. Ann Entomol Soc Am 95: 379–386. doi: 10.1603/0013-8746(2002)095[0379:SLHTOB]2.0.CO;2
    [29] Sheppard C (1983) House Fly and Lesser Fly Control Utilizing the Black Soldier Fly in Manure Management Systems for Caged Laying Hens. Environ Entomol 12: 1439–1442. doi: 10.1093/ee/12.5.1439
    [30] Čičková H, Newton GL, Lacy RC, et al. (2015) The use of fly larvae for organic waste treatment. Waste Manage 35: 68–80. doi: 10.1016/j.wasman.2014.09.026
    [31] Nguyen TTX, Tomberlin JK, Vanlaerhoven S (2015) Ability of Black Soldier Fly (Diptera: Stratiomyidae) larvae to recycle food waste. Environ Entomol 44: 406–410. doi: 10.1093/ee/nvv002
    [32] Li W, Li Q, Zheng L, et al. (2015) Potential biodiesel and biogas production from corncob by anaerobic fermentation and black soldier fly. Bioresource Technol 194: 276–282. doi: 10.1016/j.biortech.2015.06.112
    [33] Spranghers T, Ottoboni M, Klootwijk C, et al. (2016) Nutritional composition of black soldier fly (Hermetia illucens) prepupae reared on different organic waste substrates. J Sci Food Agr 97: 2594–2600.
    [34] Cammack AJ, Tomberlin KJ (2017) The impact of diet protein and carbohydrate on select life-history traits of the black soldier fly Hermetia illucens (L.) (Diptera: Stratiomyidae). Insects 8: 56.
    [35] Cheng JYK, Chiu SLH, Lo IMC (2017) Effects of moisture content of food waste on residue separation, larval growth and larval survival in black soldier fly bioconversion. Waste Manage 67: 315–323.
    [36] Lardé G (1989) Investigation on some factors affecting larval growth in a coffee-pulp bed. Biol Wastes 30: 11–19. doi: 10.1016/0269-7483(89)90139-0
    [37] Gupta S, Lee JJL, Chen WN (2018) Analysis of improved nutritional composition of potential functional food (Okara) after probiotic solid-state fermentation. J Agr Food Chem 66: 5373–5381. doi: 10.1021/acs.jafc.8b00971
    [38] Jucker C, Erba D, Leonardi MG, et al. (2017) Assessment of vegetable and fruit substrates as potential rearing media for Hermetia illucens (diptera: stratiomyidae) larvae. Environ Entomol 46: 1415–1423. doi: 10.1093/ee/nvx154
    [39] Rehman Ku, Rehman A, Cai M, et al. (2017) Conversion of mixtures of dairy manure and soybean curd residue by black soldier fly larvae (Hermetia illucens L.). J Clean Prod 154: 366–373.
    [40] Tschirner M, Simon A (2015) Influence of different growing substrates and processing on the nutrient composition of black soldier fly larvae destined for animal feed. J Insects Food Feed 1: 249–259. doi: 10.3920/JIFF2014.0008
    [41] Ma J, Lei Y, Rehman KU, et al. (2018) Dynamic effects of initial pH of substrate on biological growth and metamorphosis of black soldier fly (Diptera: Stratiomyidae). Environ Entomol 47: 159–165. doi: 10.1093/ee/nvx186
    [42] Paz ASP, Carrejo NS, Rodríguez CHG (2015) Effects of larval density and feeding rates on the bioconversion of vegetable waste using black soldier fly larvae Hermetia illucens (L.), (Diptera: Stratiomyidae). Waste Biomass Valorization 6: 1059–1065. doi: 10.1007/s12649-015-9418-8
    [43] Harnden LM, Tomberlin JK (2016) Effects of temperature and diet on black soldier fly, Hermetia illucens (L.) (Diptera: Stratiomyidae), development. Forensic Sci Int 266: 109–116. doi: 10.1016/j.forsciint.2016.05.007
    [44] Liu X, Chen X, Wang H, et al. (2017) Dynamic changes of nutrient composition throughout the entire life cycle of black soldier fly. PLoS One 12: e0182601. doi: 10.1371/journal.pone.0182601
    [45] Dortmans B, Diener S, Verstappen B, et al. (2017) Black soldier fly biowaste processing. A step-by step guide.
    [46] Kok R (2017) Insect Production and Facility Design, In: Van Huis A, Tomberlin JK, Insects as food and feed: From production to consumption, 8 Eds., Wageningen: Wageningen Academic Publishers, 143–172.
    [47] Basset Y, Cizek L, Cuénoud P, et al. (2012) Arthropod diversity in a tropical forest. Science 338: 1481–1484. doi: 10.1126/science.1226727
    [48] Engel P, Moran NA (2013) The gut microbiota of insects-diversity in structure and function. FEMS Microbiol Rev 37: 699–735. doi: 10.1111/1574-6976.12025
    [49] The Human Microbiome Project C (2012) Structure, function and diversity of the healthy human microbiome. Nature 486: 207–214. doi: 10.1038/nature11234
    [50] Chaucheyras-Durand F, Durand H (2010) Probiotics in animal nutrition and health. Benefic Microbes 1: 3–9. doi: 10.3920/BM2008.1002
    [51] Zheng L, Crippen TL, Holmes L, et al. (2013) Bacteria mediate oviposition by the black soldier fly, Hermetia illucens (L.), (Diptera: Stratiomyidae). Sci Rep 3: 2563.
    [52] Chandler JA, Lang JM, Bhatnagar S, et al. (2011) Bacterial communities of diverse drosophila species: Ecological context of a host–microbe model system. PLoS Genet 7: e1002272. doi: 10.1371/journal.pgen.1002272
    [53] Klammsteiner T, Walter A, Heussler C, et al. (2018) Hermetia illucens (Diptera: Stratiomyidae) larvae in waste valorization and diet-based shifts in their gut microbiome. Available from: http://uest.ntua.gr/naxos2018/proceedings/pdf/86_NAXOS2018_Klammsteiner_etal.pdf.
    [54] Vogel H, Muller A, Heckel DG, et al. (2018) Nutritional immunology: Diversification and diet-dependent expression of antimicrobial peptides in the black soldier fly Hermetia illucens. Dev Comp Immunol 78: 141–148. doi: 10.1016/j.dci.2017.09.008
    [55] Erickson MC, Islam M, Sheppard C, et al. (2004) Reduction of Escherichia coli O157:H7 and Salmonella enterica serovar Enteritidis in chicken manure by larvae of the black soldier fly. J Food Prot 67: 685–690. doi: 10.4315/0362-028X-67.4.685
    [56] Liu Q, Tomberlin JK, Brady JA, et al. (2008) Black soldier fly (Diptera: Stratiomyidae) larvae reduce Escherichia coli in dairy manure. Environ Entomol 37: 1525–1530. doi: 10.1603/0046-225X-37.6.1525
  • 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
  • © 2018 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(9292) PDF downloads(2516) Cited by(15)

Figures and Tables

Figures(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog