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

1.   Introduction

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

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

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

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

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

with initial conditions

and homogeneous Neumann boundary conditions

where $\Omega$ is a bounded domain in $\mathbb{R}^n$ with smooth boundary $\partial \Omega$, and $\frac{\partial}{\partial \vec{n}}$ denotes the outward normal derivative on $\partial \Omega$. $\Delta$ is the Laplacian operator. $d_i\, (i = 1, 2, 3, 4, 5)$ are the diffusion coefficients. $T(x, t)$, $L(x, t)$, $G(x, t)$, $V(x, t)$, and $Z(x, t)$ denote the concentration of uninfected cells, latently infected cells, infected cells, viruses, and B cells at position $x$ and time $t$, respectively. $\eta\; (0 < \eta < 1)$ is the probability that the uninfected cell will turn into an infected cell. $\alpha$ is the conversion rate. $n(T)$ denotes the growth of the uninfected cells. $\mu h_1(L)$, $\sigma h_2(G)$, $kh_3(V)$, and $\beta h_4(Z)$ are the death rates of the latently infected cells, infected cells, viruses, and B cells, which only depend on its concentration. $\gamma$ is the production rate. $\pi_1(T, V)$ and $\pi_2(T, G)$ are the virus-to-cell and cell-to-cell incidence rates, respectively. Let $\rho h_3(V)h_4(Z)$ and $\delta h_3(V)h_4(Z)$ be the neutralization rates of viruses and activation rate of B cells, respectively. $e ^{-a_1\tau_1}$ and $e ^{-a_1\tau_2}$ represent the probability of an infected cell surviving to the stage of $\tau_1$ and $\tau_2$, respectively. $\chi$ is the generation rate of B cells.

Define

In this paper, we first introduce the following assumptions:

$(A_1)$ $n(T)$ is continuously differentiable, and there exists $T_0 > 0$ such that $n(T_0) = 0$ and $n'(T_0) < 0$.

$(A_2)$ $\pi_i(T, \theta)$ is continuously differentiable; $\pi_i(T, \theta) > 0$ for $T\in(0, \infty), \; \theta\in(0, \infty);$ $\pi_i(T, \theta) = 0$ if and only if $T = 0$ or $\theta = 0.$ $\frac{\partial \pi_i(T, \theta)}{\partial T} > 0$ and $\frac{\partial \pi_i(T, \theta)}{\partial \theta} > 0,$ for all $T > 0$ and $\theta > 0, i = 1, 2.$ $\pi_{i1}(T) > 0$ and $\pi_{i1}'(T) > 0$ for all $T > 0, i = 1, 2.$

$(A_3)$ $h_i$ is strictly increasing on $[0, +\infty)$, $h_i(0) = 0$, $h_i'(0) = 1$, $\lim_{\theta\rightarrow 0}h_i(\theta) = \infty$, and there exists $\varrho_i > 0$ such that $h_i(\theta)\geq\varrho_i\theta$ for any $\theta\geq0, i = 1, 2, 3, 4$.

$(A_4)$ $\frac{\pi_1(T, V)}{h_3(V)}$ is non-increasing with respect to $V\in(0, +\infty)$ and $\frac{\pi_2(T, G)}{h_2(G)}$ is non-increasing with respect to $G\in(0, +\infty)$.

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

2.   Positivity, boundedness, and equilibrium

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

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

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

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

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

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

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

and then, it can be obtained that

where

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

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

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

and

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

where

Therefore,

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

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

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

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

Define

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

3.   Stability analysis

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

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

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

Proof: Define a Lyapunov functional

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

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

and

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

Moreover, by utilizing assumption $(A_5)$, we obtain

Therefore, we obtain

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

Thus, we obtain

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

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

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

Proof: Define a Lyapunov functional

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

By using

and by the divergence theorem,

Thus, we have

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

4.   Numerical simulations

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

and

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

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

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

5.   Conclusions

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

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

Use of Generative-AI tools declaration

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

Acknowledgments

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

Conflict of interest

The author declares that there are no conflicts of interest.