Dynamic analysis of a cytokine-enhanced viral infection model with infection age

1.   Introduction

In recent years, much attention have been paid on the mathematical modeling for HIV infection ^{[1,2,3,4,5,6,7,8,9,10]}. However, most of the above mentioned models assumed that the death rate and virus production rate of infected cells are constants. In reality, the results in ^[11,12] shown that the death rate of infected cells should depend on the infection age of infected cells, i.e., the time since the infection of the cell. When taking the effect of infection age into a model, Nelson et al. ^[13] proposed and studied the following age-structured viral infection model

with boundary condition $i(t, 0) = \beta TV$. Here, $T$ and $V$ denote the concentration of uninfected target cells and free virions, respectively. $i(t, a)$ denotes the density of infected cells of infection age $a$ at time $t$. $\delta(a)$ is the age-dependent death rate of infected cells, $p(a)$ is the viral production rate of an infected cell with age $a$. The global dynamics of model (1.1) has been investigated by constructing Lyapunov functions in ^[14]. Then, many age-structured viral infection models have been studied by researchers ^{[15,16,17,18,19,20,21,22,23,24,25,26]} and references cited in.

Virus-to-cell infection has been regarded as the main infection mode for HIV infection for a long time. However, literature reveals that cell-to-cell infection is a more potent and efficient way of virus propagation rather than virus-to-cell infection ^{[27,28,29,30,31]}. Motivated by this fact that many models have been proposed to study the virus dynamics for a model with cell-to-cell infection ^{[32,33,34,35,36,37,38,39,40,41,42,43,44,45]}, and references cited in. For example, Xu et al. ^[45] studied a viral infection model by taking cell-to-cell infection into consideration in model (1.1) and the global dynamics have been investigated.

As in the above just mentioned references, most of the existing work just considering the death of CD4$^+$ T cells caused by apoptosis. However, it is reported recently that CD4$^+$ T cells death caused by pyroptosis has a larger percentage than apoptosis ^{[46,47,48,49]}. Pyroptosis is a kind of programmed cell death triggered during the procedure of infection, in which the cytoplasmic content of infected cells and pro-inflammatory cytokines are released. When virus enters the CD4$^+$T cells that are unlicensed to viral infection, then the caspase-1 pathway will be activated to induce pyroptosis, which can secrete inflammatory cytokines, and then these inflammatory cytokines establish a chronic inflammation state and can attract more CD4$^+$ T cells to the inflammatory state resulting in more infection and cell death ^{[46,47,48,49]}. In recent years, though some researchers have paid attention on the cytokine-enhanced infection ^{[49,50,51,52]}. However, to the best of our knowledge the above mentioned work have not taken the infection age into consideration. Therefore, motivated by ^{[13,49,50,51,52]}, we propose a novel age-structured viral infection model incorporating virus-to-cell, cell-to-cell and cytokine-enhanced viral infection modes. Namely, we consider

with initial condition

where $M(t)$ denotes the concentration of inflammatory cytokines. Here, $\beta_1TV$ and $\beta_2TM$ are the free-virus infection and cytokine-enhanced viral infection, respectively. $\Lambda$ is the uninfected target cell production rate. $d_1$, $d_3$ and $d_4$ are the natural death rates of uninfected cells, inflammatory cytokines and virus, respectively. $\beta_3(a)$ is the infection-age specific transmission rate of productively infected cells. $d_2(a)$ is the natural death rate of infected cells with age $a$. $\alpha_1(a)$ is the death rate of infected cells which are caused by pyroptotis with age $a$. $\alpha_2(a)$ and $p(a)$ are the production rates of inflammatory cytokines released from infected cells and virus production rate with age $a$. All parameters in system (1.2) are positive, and $\alpha_1(a)$, $\beta_3(a)$, $d_2(a)$, $\alpha_2(a)$ and $p(a)$ are all Lipschitz continuous and belong to $L_+^{\infty}(0, \infty)$. Assume $i(t, a)$ is expected to be small for large $a$ which reflecting the fact that it is essentially zero for large age. Actually, it can be shown that $I(t)$ is bounded. Thus, assume that there exists a maximum age $0 < a_{\dagger} < +\infty$ such that $i(t, a) = 0$ for all $a\geq a_{\dagger}$ implying that no cells can live forever. It is worth mentioning that model (1.2) includes the existing models in ^[13,36,45]. The aim of this paper is to investigate the global dynamics of model (1.2).

The organization of this paper is as follows. Some preliminaries results of the system (1.2) are presented in Section 2. The existence and stability of steady states are analyzed in Section 3. Some numerical simulations are carried out for evaluating the results in Section 4. A brief discussion is presented at the end of the paper.

2.   Preliminaries

Denote $\bar{\delta} = ess\sup\limits_{a\in\mathbb{R}_+}\delta(a)$ and $\underline{\delta} = ess\inf\limits_{a\in\mathbb{R}_+}\delta(a) > 0$, where $\delta\in\{\alpha_1(a), \beta_3(a), d_1(a), \alpha_2(a), p(a)\}$. Let $\mathbb{X} = \mathbb{R}\times L^1(\mathbb{R}_+, \mathbb{R})\times\mathbb{R}\times\mathbb{R}$, $\mathbb{X}_0 = \mathbb{R}\times\{0\}\times L^1(\mathbb{R}_+, \mathbb{R})\times\mathbb{R}\times\mathbb{R}$, $\mathbb{X}_+ = \mathbb{R}_+\times L^1(\mathbb{R}_+, \mathbb{R})\times\mathbb{R}_+\times\mathbb{R}_+$, $\mathbb{X}_{0+} = \mathbb{X}_+\cap\mathbb{X}_0$ with the norm $\|\psi_1, \varphi(\cdot), \psi_2, \psi_3\|_{\mathbb{X}} = |\psi_1|+\int_0^{\infty}|\varphi(a)|da+|\psi_2|+|\psi_3|$. Define a linear operator $\mathcal{B}: Dom(\mathcal{B})\subset\mathbb{X}\rightarrow \mathbb{X}$ with the form

with $Dom(\mathcal{B}) = \mathbb{R}\times\{0\}\times W^{1, 1}(0, \infty)\times\mathbb{R}\times\mathbb{R}$, where $W^{1, 1}(0, \infty)$ is a Sobolev space. The nonlinear operator $\mathcal{F}: Dom(\mathcal{B})\subset\mathbb{X}\rightarrow \mathbb{X}$ is given by

and $\mathcal{F}$ is Lipschitz continuous on bounded sets. Let

where $\mathrm{T}$ represents transposition of a vector. Then, we can reformulate model (1.2) as the following abstract Cauchy problem:

Denote $\rho(\mathcal{B})$ as the resolvent set of $\mathcal{B}$. We will show $\mathcal{B}$ is a Hille-Yosida operator.

Definition 2.1. (^[53]) A linear operator $\mathcal{B}: Dom(\mathcal{B})\subset\mathbb{X}\rightarrow \mathbb{X}$ on a Banach space $(\mathbb{X}, \|\cdot\|)$ is called Hille-Yosida operator if there exist real constants $M_1\geq1$ and $\omega\in\mathbb{R}$ such that $(\omega, +\infty)\subseteq\rho(\mathcal{B})$ and

Lemma 2.1. The operator $\mathcal{B}$ is a Hille-Yosida operator.

Proof. It follows from the definition of $\mathcal{B}$ that

Then, we have

Let $\xi = \left(\tilde{\psi}_1, \left(\begin{array}{r} \tilde{\varphi}_0\\ \tilde{\varphi}(a) \end{array} \right), \tilde{\psi}_2, \tilde{\psi}_3 \right)^{\mathrm{T}}$. Then, we have

where $\mu_0 = \min\{d_1, \underline{\alpha}_1+\underline{d}_2, d_3, d_4\}$. Hence, it follows from the Definition 2.1 that $\mathcal{B}$ is a Hille-Yosida operator.

Let $\chi_0 = \left(T_0, \left(\begin{array}{r} 0\\ i_0 \end{array} \right), M_0, V_0 \right)^{\mathrm{T}}\in\mathbb{X}_{0+}$. Then it follows from ^[53] that the following result holds.

Theorem 2.1. There exists a uniquely determined semi-flow $\{U(t)\}_{t\geq0}$ on $\mathbb{X}_{0+}$ such that for each $\chi_0$, there exists a unique continuous map $U\in C([0, \infty), \mathbb{X}_{0+})$ which is an integrated solution of Cauchy problem (2.1), that is

Let $\mathcal{D} = \{(T(t), i(t, a), M(t), V(t))\in \mathbb{X}_{0+}|T(t)\leq\frac{\Lambda}{d_1}, T(t)+\int_0^{\infty}i(t, a)da\leq\frac{\Lambda}{d_0}, M(t)\leq\frac{\Lambda\bar{\alpha}_2}{d_0d_3}, V(t) \leq\frac{\Lambda\bar{p}}{d_0d_4}\}$, where $d_0 = \min\{d_1, \underline{\alpha_1}+\underline{d_2}\}$. Then, it can be shown that $\mathcal{D}$ is a positively invariant set under semi-flow $\{U(t)\}_{t\geq0}$.

Integrating the second equation of model (1.2) along the characteristic line yields

where $\Gamma(a) = e^{-\int_0^a(\alpha_1(\tau)+d_2(\tau))d\tau}$. Consequently, we have $I(t) = \int_0^{t}i(t-a, 0)\Gamma(a)da +\int_t^{\infty}i_0(a-t)\frac{\Gamma(a)}{\Gamma(a-t)}da$. It follows from model (1.2) that $\frac{dT(t)}{dt}\leq\Lambda -d_1T$, which implies $\limsup\limits_{t\rightarrow \infty}T(t)\leq\frac{\Lambda}{d_1}$. From the first two equations of model (1.2), then we have

which yields

where $d_0 = \min\{d_1, \underline{\alpha_1}+\underline{d_2}\}$. Therefore, $U(t)\chi_0\in\mathcal{D}$ for $\chi_0\in\mathcal{D}$, which implies $\mathcal{D}$ is a positively invariant set. Moreover, the semi-flow ${U(t)}_{t\geq0}$ is point dissipative and $\mathcal{D}$ attracts all positive solutions of model (1.2) in $\mathbb{X}_{0+}$. Thus, we have the following result.

Theorem 2.2. $\mathcal{D}$ is positively invariant set under the semi-flow $\{U(t)\}_{t\geq0}$. Moreover, the semi-flow $\{U(t)\}_{t\geq0}$ is point dissipative and attracts all positive solutions of model (1.2).

Furthermore, we can show that the semi-flow $\{U(t)\}_{t\geq0}$ is asymptotically smooth. In order to give the proof, we rewrite $U = \Phi+\Psi$ where

with

Theorem 2.3. ${U(t)\chi_0: t\geq0}$ has compact closure in $\mathbb{X}$ for $\forall \chi_0\in\mathcal{D}$ if the following two conditions hold: (i) There exists a function $\Delta: \mathbb{R}_+\times\mathbb{R}_+\rightarrow \mathbb{R}_+$ such that $\lim\limits_{t\rightarrow \infty}\Delta(t, r) = 0, \forall r > 0$, and if $\chi_0\in\mathcal{D}$ with $\|\chi_0\|_{\mathbb{X}}\leq r$, then $\|\Phi(t)\chi_0\|_{\mathbb{X}}\leq\Delta(t, r)$ for $t\geq0$; (ii) For $t\geq0$, $\Psi(t)\chi_0$ maps any bounded sets of $\mathcal{D}$ into sets with compact closure in $\mathbb{X}$.

Proof. (i) Let $\Delta(t, r) = re^{-(\underline{\alpha_1}+\underline{d_2})t}$, then we have $\lim\limits_{t\rightarrow \infty}\Delta(t, r) = 0$. For $\chi_0\in\mathcal{D}$ satisfying $\|\chi_0\|_{\mathbb{X}}\leq r$, we have

This completes the proof of (i).

(ii) In order to show (ii) is true. We just need to show that the following conditions hold ^[54].

(a) The supremum of $\int_0^{\infty}w_2(t, a)da$ with respect to $\chi_0\in\mathcal{D}$ is finite;

(b) $\lim\limits_{h\rightarrow \infty}\int_h^{\infty}w_2(t, a)da = 0$ uniformly with respect to $\chi_0\in\mathcal{D}$;

(c) $\lim\limits_{h\rightarrow 0^{+}}\int_0^{\infty}(w_2(t, a+h)-w_2(t, a))da = 0$ uniformly with respect to $\chi_0\in\mathcal{D}$;

(d) $\lim\limits_{h\rightarrow 0^{+}}\int_h^{\infty}w_2(t, a)da = 0$ uniformly with respect to $\chi_0\in\mathcal{D}$.

It follows from the definition of $\mathcal{D}$ that (a), (b) and $(d)$ hold. Thus, we only to show the condition $(c)$ holds. For convenience, denote $K(t) = \int_0^{\infty}\beta_3(a)i(t, a)da$. For sufficiently small $h\in(0, t)$, we have

where

We note that $\Gamma(a)$ is a decreasing function and satisfies $0\leq\Gamma(a)\leq1$. Thus,

Then, it follows from (2.2) and (2.3) that

It follows from ^[38,55] that $TV$, $TM$ and $KT$ are Lipschitz on $\mathbb{R}_+$. Assume $L_{TV}$, $L_{TM}$ and $L_{KT}$ be the Lipschitz coefficients of $TV$, $TM$ and $KT$, respectively. Similar techniques as ^[56], then we have

Thus,

which converges to 0 as $h\rightarrow 0^+$. Therefore, condition (iii) holds. This completes the proof.

3.   Stability of steady states

Define the basic reproduction number of model (1.2) as follows:

Clearly, model (1.2) always has an infection-free steady state $E_0 = (T^0, 0, 0, 0)$ with $T^0 = \frac{\Lambda}{d_1}$. If $\mathcal{R}_0\neq1$, then there exists a unique infection steady state $E_* = (T^*, i^*(a), V^*, M^*)$, which satisfies

Theorem 3.1. If $\mathcal{R}_0 < 1$, then the infection-free steady state $E^0$ is locally stable.

Proof. Denote $\tilde{T}(t) = T(t)-T^0$, $\tilde{i}(t, a) = i(t, a)$, $\tilde{M}(t) = M(t)$ and $\tilde{V}(t) = V(t)$. Linearizing model (1.2) at $E^0$, then we have

In order to analyze the stability of $E_0$, we look for solutions of the form $\tilde{T} = \tilde{T}_0e^{\mu t}$, $\tilde{i}(t, a) = \tilde{i}_0(a)e^{\mu t}$, $\tilde{M} = \tilde{M}_0e^{\mu t}$ and $\tilde{V} = \tilde{V}_0e^{\mu t}$. Substituting the solutions into the above linearized model yields

Then we have $\tilde{i}_0(a) = \tilde{i}_0(0)\Gamma(a)e^{-\mu a}$, $\tilde{M}_0 = \frac{1}{\mu+d_3}\int_0^{\infty}\alpha_2(a)\tilde{i}_0(a)da$, $\tilde{V}_0 = \frac{1}{\mu+d_4}\int_0^{\infty}p(a)\tilde{i}_0(a)da$, substituting $\tilde{i}_0(a)$, $\tilde{M}_0$ and $\tilde{V}_0$ into the last equation of (3.2), we have

Obviously, $\lim\limits_{\mu\rightarrow \infty}G(\mu) = 0$, $G(0) = \mathcal{R}_0$, and a simple computation shows that $G(\mu)$ is a decreasing function with respect to $\mu$. Therefore, if $\mathcal{R}_0 < 1$, then any real root of (3.3) is negative. Thus, $E_0$ is unstable for $\mathcal{R}_0 > 1$. Moreover, we claim that (3.3) has no complex roots with nonnegative real part if $\mathcal{R}_0 < 1$. In fact, if there exists a root $\mu = \xi+\eta i$ with $\xi\geq0$. Then,

Thus, the above arguments imply that every root of (3.3) must have a negative real part, which implies that $E_0$ is locally stable for $\mathcal{R}_0 < 1$.

Theorem 3.2. If $\mathcal{R}_0 > 1$, then $E_*$ is locally stable.

Proof. With the same technique of Theorem 3.1. Let $x(t) = T-T^*$, $y(t, a) = i(t, a)-i^*(a)$, $z = M(t)-M^*$ and $v(t) = V(t)-V^*$. Linearizing model (1.2) at $E_*$ and looking for solutions of the form $x = \tilde{x}e^{\mu t}$, $y(t, a) = \tilde{y}(a)e^{\mu t}$, $z = \tilde{z}e^{\mu t}$ and $v = \tilde{v}e^{\mu t}$ leads to

A simple calculation leads to $\tilde{y}(a) = \tilde{y}(0)e^{-\mu a}\Gamma(a)$, $\tilde{x} = -\frac{\tilde{y}(0)}{\mu+d_1}$, $\tilde{z} = \frac{\tilde{y}(0)}{\mu+d_3}\int_0^{\infty}\alpha_2(a)e^{-\mu a}\Gamma(a)da$, and $\tilde{v} = \frac{\tilde{v}(0)}{\mu+d_4}\int_0^{\infty}p(a)e^{-\mu a}\Gamma(a)da$. Then we have

Obviously, for $Re{\mu}\geq0$ then

and

It follows from (3.4)–(3.6) that there are no characteristic roots with non-negative real part. Thus, $E_*$ is locally stable.

Before we discuss the global stability of model (1.2) by constructing Lyapunov functionals. We present the following result by a procedure similar to ^[22,23], so we omit the proof.

Theorem 3.3. If $\mathcal{R}_0 > 1$, then for each $\chi_0\in\mathbb{X}$ there exists a constant $\rho > 0$ such that

Now, we are in position to investigate the global stability of steady states.

Theorem 3.4. If $\mathcal{R}_0 < 1$, then the infection-free steady state $E_0$ is globally asymptotically stable.

Proof. Define

where

It is easy to show that $\Phi(0) = \mathcal{R}_0$, and

Furthermore, we can obtain

Thus, by using (3.7) and (3.8) and noting that $T^0 = \frac{\Lambda}{d_1}$, we have

Therefore, if $\mathcal{R}_0 < 1$, we have $H'(t)\leq0$ and $H'(t) = 0$ implies that $T = T^0$, $i(t, a) = 0$, $M(t) = 0$ and $V(t) = 0$. Hence, the largest invariant subset of $\{H'(t) = 0\}$ is a singleton $\{E_0\}$, which means $E_0$ is global asymptotically stable for $\mathcal{R}_0 < 1$ by Lyapunov-LaSalle theorem ^[57]. This completes the proof.

Theorem 3.5. If $\mathcal{R}_0 > 1$, then the infection steady state $E_*$ is globally asymptotically stable.

Proof. Define

where

and satisfies

Differentiating and using the steady state identities (3.1) and (3.9), then we have

Since,

Then we have

where $\phi(x) = 1+\ln x-x$ satisfies $\phi(x)\leq 0$ for $x > 0$ and $\phi(x) = 0$ if and only if $x = 1$. Thus, $H'_2\leq0$. Furthermore, it can be shown that the largest compact invariant set of $H'_2 = 0$ is the singleton $\{E_*\}$, which implies $E_*$ is globally asymptotically stable.

4.   Numerical simulations

In this part, we carry out some numerical simulations to illustrate the above obtained theoretical results. Most of the parameters are from ^[13], and the functions $p(a)$, $\alpha_2(a)$, $d_2(a)$ are given with the following forms:

For simulation, we assume $\alpha_1(a) = 0.1$ and $\beta_3(a) = 0.0000075$, $\Lambda = 100$, $\beta_1 = 0.0000046$, $\beta_2 = 0.0000065$, $d_1 = 0.1$, $d_3 = 6.6$, $d_4 = 2.4$, $a_1 = 0.2$, $a_2 = 0.5$, $a_{max} = 15$, $\gamma_1 = 10$, $\gamma_2 = 5$, $\delta_0 = 0.05$, $\delta_m = 0.35$, $\alpha_{20} = 1000$, $\gamma = 1$, $p_{max} = 850$, computation yields $\mathcal{R}_0 = 0.4557 < 1$, which implies that the infection-free steady state is globally asymptotically stable as shown in Figure 1. When $\Lambda = 10$, $\beta_1 = 0.0000046$, $\beta_2 = 0.0000065$, $d_1 = 0.01$, $d_3 = 6.6$, $d_4 = 2.4$, $a_1 = 0.2$, $a_2 = 0.5$, $a_{max} = 15$, $\gamma_1 = 10$, $\gamma_2 = 5$, $\delta_0 = 0.05$, $\delta_m = 0.35$, $\alpha_{20} = 1000$, $\gamma = 1$, $p_{max} = 1880$, computation yields $\mathcal{R}_0 = 14.7759 > 1$, which implies that the infection steady state is globally asymptotically stable as shown in Figure 2. Moreover, it follows from Figure 3 that the existence of cytokine-enhanced effect can lead to a higher peak of viral load. Also, the formula of the basic reproduction number implies that it may be under-evaluated without considering cytokine-enhanced viral infection.

5.   Summary and discussion

In this paper, an age-structured virus infection model in which the cytokine-enhanced viral infection have been taken into consideration. By constructing Lyapunov functionals, we show that the global properties of the model are completely determined by the basic reproduction numbers $\mathcal{R}_0$: if $\mathcal{R}_0 < 1$, then the infection-free equilibrium is globally asymptotically stable and the infection dies out; if $\mathcal{R}_0 > 1$, there exists a unique infection steady state which is globally asymptotically stable. Recall that

The first term of $\mathcal{R}_0$ is induced by viral infection, the second term is induced by cytokine-enhanced viral infection, and the third term corresponds to the cell-to-cell infection mode. Thus, the basic reproduction number will be under-evaluated without considering cytokine-enhanced viral infection, and it may not be enough to eliminate the infection by decreasing the basic reproduction number just for virus-to-cell or cell-to-cell infection.

In most virus infection process, cytotoxic T lymphocytes (CTLs) play a critical role in antiviral defense by attacking virus-infected cells. Thus, it would be very interesting to improve the current work by considering both CTL responses and antibody response in an age-structured viral infection model. Besides, the motion of the virus should also be taken into consideration. This will result in a reaction-diffusion model with age-structure. Whether the improved models can preserve these global results is an interesting problem and we leave this as a future work.

Acknowledgments

This work was supported by National Natural Science Foundation of China (#11701445, #11971379), by Natural Science Basic Research Programm in Shaanxi Province of China (2022JM-042, 2022JM-038, 2022JQ-033, 2020JQ-831).

Conflict of interest

The authors declare there is no conflict of interest.