1. Introduction

Mathematical modeling has been proven to be valuable in understanding virus infection and immune responses. Models and analysis can provide important insights into the dynamics of viral load in vivo and offer helpful suggestions for clinical treatment. Many of the models are based on a differential equation system, which describes the coupled changes in target cells, infected cells, and free virus particles through time in a single compartment (i.e. the blood) of an infected individual. An example is the application of models to hepatitis B virus (HBV) infection. Persistent infection with HBV is a major health problem worldwide. HBV infection can lead to cirrhosis and primary hepatocellular carcinoma [2,42]. Chronic HBV infection is usually the result of exposure to virus early in life, leading to viral persistence in the absence of strong antibody or cellular immune responses [8]. Treatment of HBV carriers aims to either inhibit viral replication or enhance immunological responses against the virus, or both [26]. Based on the clinical experiment of chronic HBV carriers treated with various doses of lamivudine, Nowak and Bangham [23] used a basic mathematical model to study HBV dynamics. Nowak et al.[24] provided a quantitative understanding of HBV replication dynamics in vivo, estimated the turnover rates of infected cells and virus, and suggested the optimal timing of drug treatment and immunotherapy in chronic HBV infection. Some other within-host virus dynamics models have also been developed to study HBV [33,34,35], HIV [17,18,46], and hepatitis C virus (HCV) infection [27,29].

The age structure of population has been widely investigated in epidemiological models [3,20,37,41]. Because of its flexibility in modeling viral production and mortality of infected cells, age structure of infected cells has also been incorporated into within-host virus infection models [9,22,28,39,40]. For example, Nelson et al. [22] developed and analyzed an age-structured HIV model that allowed the variation in the production rate of virus and the death rate of infected CD4+ T cells. For some special functions, they performed a local stability analysis of the nontrivial equilibrium solution. They used numerical methods to show that the time to reach peak viral levels in the blood depends on both initial conditions and the way in which viral production ramps up. Because the age structure of infection allows the incorporation of different classes of antiretroviral drugs that target different stages of viral lifecycle, Rong et al.[28] used the age-structured model to compare the treatment effectiveness of administrating different drugs [28]. They conducted analysis of the model under treatment for general functions of the death rate of infected cells and viral production rate. Using an age-structured model, Gilchrist et al. [9] also explored how an infected cell's viral production rate can affect the relative fitness of a virus within a host. They performed an invasion analysis to discuss the strategy for achieving the maximum relative viral fitness. Recently Wang et al.[39] analyzed an age-structured HIV model with both virus-to-cell infection and cell-to-cell transmission.

The global stability of age-structured within-host models is the focus of a few recent studies. Huang et al.[12] studied the global asymptotic behavior of an age-structured HIV infection model. Browne et al.[5] studied the within-host viral infection with an explicit age-since-infection structure of infected cells. Browne [4] also considered the global stability of within-host viral infection with multiple strains. Shen et al. analyzed a model that links the between-host and within-host dynamics of HIV infection [30]. More age-structured within-host models can be found in the references [1,16,29,36].

Viruses can infect different populations of target cells. For example, HIV mainly infect CD4+ T cells. However, other cells such as macrophages [13] and dendritic cells [25] are also known to be susceptible to HIV infection. Similar to HIV infection, HCV can also infect different classes of cells. HCV replicates mainly in the hepatocytes of the liver but the virus can also replicate in peripheral blood mononuclear cells [7]. Thus, multi-compartment mathematical models are needed to study virus infection in different populations of cells [38]. In this note, we will study a within-host virus dynamics model including the age of infection and multiple populations of target cells. We will prove that the solution of the system is positive and bounded. Using the methods semigroup and Lyapunov function, we will investigate the global asymptotic property of the infected steady state of the model. Under some special scenarios, we will show that the age of infection model is equivalent to an ordinary differential equation system with or without time delays.

---

2. Model description

We consider a general age-structured within-host virus dynamics model that includes multiple classes of target cells. The population is divided into $2n+1$ classes: uninfected cells, $T_{j},$ infected cells, $i_{j},$ where $j=1,2,...,n,$ and free virus, $V.$ The model is given by the following system:

$\left\{\begin {array}{l l} \frac{d T_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}, \\ \frac{\partial i_{j}(a,t)}{\partial t}+\frac{\partial i_{j}(a,t)}{\partial a}=-\delta_{j} (a)i_{j}(a,t), \\ \frac{d V(t)}{d t}= \sum\limits_{j=1}^{n}\int_{0}^{\infty}p_{j}(a)i_{j}(a,t)da-cV(t), \end {array} \right.$

(2.1)

for $t>0,$ with initial conditions

$\begin {array}{l l}T_{j}(0)=T^{0}_{j}\geq0,~i_{j}(a,0)=i^{0}_{j}(a), j=1, 2,...,n; V( 0)=V^{0}\geq0,\\ \end {array}$

(2.2)

and boundary conditions

$i_{j}(0,t)=\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}.$

(2.3)

For each class (denoted by the subscript $j$, $j$=1, 2, ..., n) of target cells, $T$ represents the population of uninfected cells, $i(a,t)$ denotes the population of infected cells with the infection age $a$ at time $t$, and $V$ is the population of free virus. The parameter $s_{j}$ is the production rate of uninfected cells and $d_{j}$ is the death rate of uninfected cells. Here we use a saturation-dependent functional response ($kVT/(1+\alpha V)$ with $\alpha>0$) to describe the infection of cells by virus. A similar function has been used in other within-host models [31,36,43,44]. The distribution function $i^{0}_{j}(a)\in L^{1}_{+}((0,+\infty),\mathcal {R})$ is the initial condition. The function $\delta_{j}(a)$ is the age-dependent per capita death rate of infected cells and $p_{j}(a)$ is the viral production rate of an infected cell with age $a$. The parameter $c$ is the viral clearance rate.

An example of the death rate of infected cells can be chosen to be the same as that in [22], which is given as follows:

$\delta_{j}(a)=\left\{\begin {array}{l l} \delta^{0}_{j}, a <a_{1},\\ \delta^{0}_{j}+\delta^{m}_{j}(1-e^{-\gamma(a-a_{1})}), a\geq a_{1}, \end {array}\right.$

(2.4)

where $a_1$ is the age at which infected cells express sufficient viral genome on the surface and are susceptible to killing by immune cells. Thus, the death rate of infected cells increases from $\delta_j^0$ at age $a_1$ to the maximum value $\delta_j^0+\delta_j^m$.

Integrating the second equation in (2.1) along the characteristic line $t-a=\mbox{constant},$ we get the following formula

$i_{j}(a,t)=\left\{\begin {array}{l l} B_{j}(t-a)\sigma_{j} (a), ~~for~~~a <t,\\ i_{j}(a-t,0) \frac{\sigma_{j} (a)}{\sigma_{j} (a-t)},~~for~~a\geq t, \end {array}\right.$

(2.5)

where

$B_{j}(t)=\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}$

and

$\sigma_{j} (a)=\exp\Big(-\int_{0}^{a}\delta_{j} (\theta)d\theta\Big).$

Using the above solution, model (2.1) can be written as the following system.

$\left\{\begin {array}{l l} \frac{d T_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}, \\ \frac{d V(t)}{d t}= \sum\limits_{j=1}^{n}\Bigg( \int_{0}^{t}p_{j}(a)\sigma_{j} (a)B_{j}(t-a)da+\widetilde{F_{j}}(t)\Bigg)-cV(t), \end {array} \right.$

(2.6)

where

$\widetilde{F_{j}}(t)=\int_{t}^{\infty}p_{j}(a)i^{0}_{j}(a-t)\frac{\sigma_{j} (a)}{\sigma_{j} (a-t)}da.$

It is clear that $\widetilde{F_{j}}(t)\rightarrow 0$ as $t\rightarrow \infty, i=1,2,...,n.$

For cells with age-dependent viral production and death rates, we define

$N_{j}=\int_{0}^{\infty}p_{j}(a)\sigma_{j} (a)da.$

$N_j$ defines the total number of virions produced by the infected cell of the $j$-th class in its life span, which is called the viral burst size of the $j$-th class. The basic reproductive number is given by

$R_{0}=\sum\limits_{j=1}^{n}\frac{N_{j}s_{j}k_{j}}{cd_{j}},$

which represents the total number of newly infected cells produced by one infected cell during its lifetime in a fully susceptible environment (i.e. assuming all cells are susceptible).

---

3. Integrated semigroup formulation and equilibria

In order to take into account the boundary condition, we expand the state space. Denote

$\mathcal {M}=\mathcal {R}\times L^{1}\Big((0,+\infty),\mathcal {R}\Big), \mathcal {M_{+}}=\mathcal {R_{+}}\times L^{1}\Big((0,+\infty),\mathcal {R}\Big),$

$\mathscr{N}=\mathcal {R}\times \{0\}\times W^{1,1}\Big((0,+\infty),\mathcal {R}\Big),$

$\mathscr{P}=\mathcal {R}\times \{0\}\times L^1\Big((0,+\infty),\mathcal {R}\Big), \mathscr{P_{+}}=\mathcal {R_{+}}\times \{0\}\times L^1\Big((0,+\infty),\mathcal {R}\Big),$

where $W^{1,1}$ is a Sobolev space. Let

$\mathscr{X}=\Big(\prod\limits_{1}^{n}\mathcal {M}\Big)\times\mathcal {R}, \mathscr{X}_{+}=\Big(\prod\limits_{1}^{n}\mathcal {M_{+}}\Big)\times\mathcal {R}_{+}$

and consider the linear operator $\mathscr{A}: Dom(\mathscr{A})\subset\mathscr{X}\rightarrow\mathscr{X}$ defined by

$\mathscr{A}\left( \begin{array}{c} T_{1}\\ \left( \begin{array}{c}0\\ i_{1}\end{array}\right)\\ \vdots\\ \begin{array}{c} T_{n}\\ \left( \begin{array}{c}0\\ i_{n}\end{array}\right)\\ \end{array}\\ V \end{array} \right) =\left( \begin{array}{c} -d_{1}T_{1}\\ \left( \begin{array}{c} i_{1}(0)\\ -i^{'}_{1}-\delta(a)i_{1} \end{array}\right)\\ \vdots\\ -d_{n}T_{n}\\ \left( \begin{array}{c} i_{n}(0)\\ -i^{'}_{n}-\delta(a)i_{n} \end{array}\right)\\ -cV \end{array} \right),$

with

$Dom(\mathscr{A})=\Big(\prod\limits_{1}^{n}\mathcal {P}\Big)\times\mathcal {R}.$

Then $\overline{Dom(\mathscr{A})}=\Big( \prod\limits_{1}^{n}\mathcal { N}\Big)\times\mathcal {R}$ is not dense in $\mathscr{X}.$ We consider a nonlinear map $\mathscr{F}: \overline{Dom(\mathscr{A})}\rightarrow\mathscr{X}$, which is defined by

$\mathscr{F}\left( \begin{array}{c} T_{1}\\ \left( \begin{array}{c}0\\ i_{1}\end{array}\right)\\ \vdots\\ T_{n}\\ \left( \begin{array}{c}0\\ i_{n}\end{array}\right)\\ V \end{array} \right) =\left( \begin{array}{c} s_{1}-\frac{k_{1} T_{1}(t)V(t)}{1+\alpha_{1} V(t)}\\ \left( \begin{array}{c} \psi T_{1}(t)\\ 0_{L^{1}} \end{array}\right)\\ \vdots\\ s_{n}-\frac{k_{ n} T_{n}(t)V(t)}{1+\alpha_{n} V(t)}\\ \left( \begin{array}{c} \psi T_{n}(t)\\ 0_{L^{1}} \end{array}\right)\\ \int_{0}^{\infty}p_{n}(a)i_{n}(a,x,t)da \end{array} \right),$

and let

$u(t)=\left( \begin{array}{c} T_{1}, \left( \begin{array}{c}0\\ i_{1}(\cdot, t)\end{array}\right), \ldots T_{n}, \left( \begin{array}{c}0\\ i_{n}(\cdot, t)\end{array}\right), V \end{array} \right)^{T}.$

Set

$\mathscr{X}_{0}=\overline{Dom(\mathscr{A})}=\Big(\prod\limits_{1}^{n}\mathcal {P}\Big)\times\mathcal {R}$

and

$\mathscr{X}_{0+}=\overline{Dom(\mathscr{A})}\cap\mathscr{X}_{+} =\Big(\prod\limits_{1}^{n}\mathcal {P_{+}}\Big)\times\mathcal {R}_{+}.$

On the basis of the above formulation, system (2.1)$\sim$(2.3) can be rewritten as the following abstract Cauchy problem:

$u(t)=\mathscr{A}u(t)+\mathscr{F}\Big(u(t)\Big), ~\mbox{for}~t\geq0, ~\mbox{with}~u(0)=x\in\mathscr{X}_{0+}.$

(3.1)

By applying the results given in Hale [10], Magal [19], Magal and Thieme [21], Thieme [32], Yang et al. [45] and Wang et al. [40,41] we can show the existence and uniqueness of the semiflow $\{U(t)\}_{t\geq0}$ on $\mathscr{X}_{0+}$ generated by system (3.1) and further have the following result.

Theorem 3.1. System (3.1) generates a unique continuous semiflow $\{U(t)\}_{t\geq0}$ on $\mathscr{X}_{0+}$ that is asymptotically smooth and bounded dissipative. Furthermore, the semiflow $\{U(t)\}_{t\geq0}$ has a compact global attractor $\mathscr{A}\subset\mathscr{X}_{0+}$.

Define

$\mathcal {M}_{0}=\bigg \{(T_{1},i_{1}(a),...,T_{n},i_{n}(a),V)^{T}\in \mathcal{X}_{0+}\big|V+\sum\limits_{j=1}^{n}\int_{0}^{\infty}i_{j}(a)da>0\bigg\}$

and

$\partial\mathcal{M}_{0}=\mathcal{X}_{0+}\setminus \mathcal{M}_{0}.$

Following Theorems 4.1 and 4.2 in [45], we can get the following theorems. Here we omit the proofs.

Theorem 3.2. $\mathcal {M}_{0}$ and $\partial\mathcal {M}_{0}$ are both positively invariant under the semiflow ${U(t)}_{t\geq0}$ generated by system (3.1) on $\mathcal{X}_{0+}.$ Moreover, the infection-free steady state $E_{0}$ of problem (2.1)$\sim$(2.3) is globally asymptotically stable for the semiflow ${U(t)}_{t\geq0}$ restricted to $\partial\mathcal {M}_{0}$.

Theorem 3.3. Assuming $R_0\leq1$, the semiflow ${U(t)}_{t\geq0}$ generated by system (3.1) is uniformly persistent with respect to the pair $(\partial\mathcal {M}_{0}, \mathcal {M}_{0})$; that is, there exists $\epsilon>0,$ such that for each $y\in\mathcal {M}_{0},$

$\liminf\limits_{t\rightarrow+\infty} d(U(t)y, \partial\mathcal {M}_{0})\geq\epsilon.$

Furthermore, there exists a compact subset $\mathcal {A}_{0}\subset\mathcal {M}_{0}$ which is a global attractor for ${U(t)}_{t\geq0}$ in $\mathcal {M}_{0}.$

---

4. The global results of steady states

In this section, we focus on the global asymptotic properties of the steady states of system (2.1)$\sim$(2.3). Generally, it can be challenging to obtain the global properties of a model with saturation response of the infection rate, especially for the model with age structure.

Theorem 4.1. System (2.1) always has an infection-free steady state $E^{0}(T_{10},0,..., T_{n0},0,0)$. When the basic reproductive ratio is greater than 1, system (2.1) has a unique positive infected steady state $E^{*}(T^{*}_{1},i_{1}^{*}(a),...,T^{*}_{n},i_{n}^{*}(a),V^{*}).$

Proof. It is clear that there always exists an infection-free steady state $E^{0}(T_{10},0,..., T_{n0},0,0)$ for system (2.1). To obtain the infected steady state $E^{*}(T^{*}_{1},$ $i_{1}^{*}(a),...,T^{*}_{n}, i_{n}^{*}(a),V^{*})$, we solve the following algebraic equations.

$\left\{\begin {array}{l l}s_{j}-d_{j}T^{*}_{j}- \frac{k_{j} T^{*}_{j}V^{*}}{1+\alpha_{j} V^{*}}=0 \\ i^{*}_{j}(a)= \frac{k_{j} T^{*}_{j}V^{*}}{1+\alpha_{j} V^{*}}\sigma_{j} (a)\\ \sum \limits_{j=1}^{n}\frac{N_{j}k_{j}T^{*}_{j}}{1+\alpha_{j}V^{*}}=c, \end {array}\right.$

(4.1)

From the first equation of (4.1), we have

$T^{*}_{j}=\frac{s_{j}}{d_{j}+\frac{k_{j}V^{*}}{1+\alpha_{j} V^{*}}}.$

Substituting into the last equation of (3.2), we get

$\sum\limits_{j=1}^{n}\frac{N_{j}k_{j}s_{j}}{d_{j}+(k_{j}+d_{j}\alpha_{j})V^{*}}=c.$

Let

$F(V)=\sum\limits_{i=1}^{n}\frac{N_{j}k_{j}s_{j}}{d_{j}+(k_{j}+d_{j}\alpha_{j})V}-c.$

On one hand, $F(V)$ is continuous and monotonically decreasing for $V\in[0, +\infty)$. We also have

$F(0)=\sum\limits_{j=1}^{n}\frac{N_{j}k_{j}s_{j}}{d_{j}}-c=\sum\limits_{j=1}^{n}\frac{N_{j}k_{j}s_{j}}{d_{j}}(1-\frac{1}{R_{0}})>0$

for $R_{0}>1.$ On the other hand,

$F(\infty)=0-c <0.$

Therefore, equation $F(V)=0$ has one positive solution $V^{*}.$ Thus, when $R_{0}>1$, systems (2.1) has a unique positive infected steady state $E^{*}.$

Let $\tilde{E}(\tilde{T_{1}},\tilde{i_{1}}(a),\tilde{T_{2}},\tilde{i_{2}}(a),...,\tilde{T_{n}},\tilde{i_{n}}(a),\tilde{V})$ be any arbitrary steady state of system (2.1)$\sim$(2.3). The linearized system of (2.1)$\sim$(2.3) is

$\left\{\begin {array}{l l} (\lambda+d_{j}+\frac{k_{j}\tilde{V}}{1+\alpha_{j} \tilde{V}})T_{j}+\frac{k_{j}\tilde{T_{j}}}{(1+\alpha_{j} \tilde{V})^{2}}V=0,\\ \frac{d i_{j}(a)}{d a}=-(\lambda+\delta_{j}(a))i_{j}(a),\\ i_{j}(0)=\frac{k_{j}\tilde{V}}{1+\alpha_{j} \tilde{V}}T_{j}+\frac{k_{j}\tilde{T_{j}}}{(1+\alpha_{j} \tilde{V})^{2}}V,\\ (\lambda+c)V-\sum\limits_{j=1}^{n} \int_{0}^{\infty}p_{j}(a)i_{j}(a)da=0. \end {array}\right.$

(4.2)

Solving the second equation of (4.2), we have

$i_{j}(a)=\Big(\frac{k_{j}\tilde{V}}{1+\alpha_{j} \tilde{V}}T_{j}+\frac{k_{j}\tilde{T_{j}}}{(1+\alpha_{j } \tilde{V})^{2}}V\Big)\sigma_{j} (a)e^{-\lambda a}.$

(4.3)

Substituting equation (4.3) into (4.2), we have

$\left\{\begin {array}{l l} (\lambda+d_{1}+\frac{k_{1}\tilde{V}}{1+\alpha_{1} \tilde{V}})T_{1}+\frac{k_{1}\tilde{T_{1}}}{(1+\alpha_{1} \tilde{V})^{2}}V=0,\\ (\lambda+d_{2}+\frac{k_{2}\tilde{V}}{1+\alpha_{2} \tilde{V}})T_{2}+\frac{k_{2}\tilde{T_{2}}}{(1+\alpha_{2} \tilde{V})^{2}}V=0,\\ ...\\ (\lambda+d_{n}+\frac{k_{n}\tilde{V}}{1+\alpha_{n} \tilde{V}})T_{n}+\frac{k_{n}\tilde{T_{n}}}{(1+\alpha_{n} \tilde{V})^{2}}V=0,\\ \Big(\lambda+c-\sum\limits_{j=1}^{n} \frac{k_{j}\tilde{T_{j}}}{(1+\alpha_{j} \tilde{V})^{2}}N_{j}(\lambda)\Big)V-\sum\limits_{j=1}^{n}\frac{k_{j}\tilde{V}}{1+\alpha_{j} \tilde{V}}N_{j}(\lambda)T_{j}=0. \end {array}\right.$

For the infection-free steady state $E_{0}$, the characteristic equation of the linearized system of (2.1)$\sim$(2.3) is

$\Big(\prod\limits_{j=1}^{n}(\lambda+d_{j})\Big)\Big(\lambda+c-\sum\limits_{j=1}^{n} k_{j}T_{j0}N_{j}(\lambda)\Big)=0.$

Let

$F(\lambda)=\lambda+c-\sum\limits_{j=1}^{n} k_{j}T_{j0}N_{j}(\lambda).$

Then

$F'(\lambda)>0$

and

$F(0)=c(1-R_{0})>0$

for $R_{0}\leq1$. Therefore, all the eigenvalues of the infection-free steady state $E_{0}$ are negative.

For $R_{0}>1,$ we have

$F(0)=c(1-R_{0}) <0, \lim\limits_{t\rightarrow+\infty}F(\lambda)=+\infty.$

Thus, $F(\lambda)=0$ has a positive root and the infection-free steady state $E_{0}$ is unstable. We summarize in the following results.

Lemma 4.1. If $R_{0}\leq1$, then the infection-free steady state $E_{0}$ is locally asymptotically stable. While $R_{0}>1$, $E_{0}$ is unstable.

Theorem 4.2. If $R_0\leq1$, then the infection-free steady state $E_{0}$ of model (2.1)$\sim$ (2.3) is globally asymptotically stable.

Proof. From the equations of target cells in system (2.1), we obtain

$\lim\limits_{t\rightarrow+\infty}T_{j}( t)\leq\frac{s_{j}}{d_{j}},~~j=1,\ldots,n.$

One can choose $\varepsilon>0$ small enough such that there exists $t_{1}$ such that $T_{j}(t)\leq \frac{s_{j}}{d_{j}}+\varepsilon$ for all $t\geq t_{1}$.

$\left\{\begin {array}{l l} \frac{dV(t)}{d t}=\sum\limits_{j=1}^{n}\int_{0}^{\infty}p_{j}(a)i_{j}(a,t)da-cV(t), t>0,\\ V( 0)=V^{0}\geq0. \end {array}\right.$

(4.4)

From (4.4), we obtain the following inequality.

$\frac{dV(t)}{d t}\leq\sum\limits_{j=1}^{n} N_{j}k_{j}(\frac{s_{j}}{d_{j}}+\varepsilon)-cV(t), t>0,$

(4.5)

where $\varepsilon$ is chosen as before. Considering $R_0\leq1$, from (4.5), we have

$\lim\limits_{t\rightarrow+\infty}\sup\limits_{}V(t)=0.$

By comparison, it follows that

$\lim\limits_{t\rightarrow+\infty}V( t)=0.$

Therefore, for $\varepsilon>0$ sufficiently small there exists a $t_{2}>0$ such that $0 <V(t) <\varepsilon$ for $t\geq t_{2}$.

Again from the equations of target cells in system (2.1), we can get

$\frac{d T_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)} \geq s_{j}-d_{j}T_{j}(t)-k_{j}\varepsilon T_{j}( t), ~~t>t_{2},$

which, together with the arbitrariness of $\varepsilon>0$, yields that

$\lim\limits_{t\rightarrow+\infty}T_{j}( t)\geq\frac{s_{j}}{d_{j}}.$

Thus,

$\lim\limits_{t\rightarrow+\infty}T_{j}( t)=\frac{s_{j}}{d_{j}}.$

Combining with Lemma 4.1, we can conclude that the infection-free steady state $E_{0}$ of model (2.1)$\sim$(2.3) is a globally asymptotically stable.

Theorem 4.3. When $R_0>1$, the infected steady state $E^{*}$ of model (2.1) $\sim$(2.3) is globally asymptotically stable.

Proof. First, we define a positive function as used in ref.[12]

$\beta_{j}(a)=\int_{a}^{\infty}p_{j}(\epsilon)e^{-\int_{a}^{\epsilon}\delta_{j}(\eta)d\eta}d\epsilon.$

Note that $\beta_{j}(a)>0$ for $0\leq a <+\infty,$ and $\beta_{j}(0)=N_{j}.$ Similar to the discussion in [12], we know that $\beta_{j}(a)$ is bounded and $\beta_{j}'(a)=\delta_{j}(a)\beta_{j}(a)-p_{j}(a).$

We consider the following Lyapunov function

$W(t) = N_{1}\Bigg(T_{1}(t)-T_{1}^{*}-T_{1}^{*}\ln \frac{T_{1}(t)}{T_{1}^{*}}\Bigg)\\ +\int_{0}^{+\infty}\beta_{1}(a)i_{1}^{*}(a)\Bigg(\frac{i_{1}(a,t)}{i_{1}^{*}(a)}-1-\ln \frac{i_{1}(a,t)}{i_{1}^{*}(a)}\Bigg)da\\ +\sum\limits_{j=2,\cdots,n}N_{j}\Bigg(T_{j}(t)-T_{j}^{*}-T_{j}^{*}\ln \frac{T_{j}(t)}{T_{j}^{*}}\Bigg)\\ +\sum\limits_{j=2,\cdots,n}\int_{0}^{\infty}\beta_{j}(a)i_{j}^{*}(a)\Bigg(\frac{i_{j}(a,t)}{i_{j}^{*}(a)}-1-\ln \frac{i_{j}(a,t)}{i_{j}^{*}(a)}\Bigg)da\\ +\Bigg(V(t)-V^{*}-V^{*}\ln\frac{V(t)}{V^{*}}\Bigg).$

It is clear to see that $W$ is nonnegative and the point $E^{*}$ is the global minimum point. Calculating the time derivative of $W$ along the solution of system (2.1)$\sim$(2.3), we obtain

$\frac{dW(t)}{dt} = N_{1}\Bigg(1-\frac{T_{1}^{*}}{T_{1}(t)}\Bigg)\Bigg(s_{1}-d_{1}T_{1}(t)-\frac{k_{1}T_{1}(t)V(t)}{1+\alpha_{1} V(t)}\Bigg) \\ +\int_{0}^{+\infty} \beta_{1}(a)i_{1}^{*}(a)\frac{\partial}{\partial t}\Bigg(\frac{i_{1}(a,t)}{i_{1}^{*}(a)}-1-\ln \frac{i_{1}(a,t)}{i_{1}^{*}(a)}\Bigg)da\\ +\sum\limits_{j=2,\cdots,n} N_{j}\Bigg(1-\frac{T_{j}^{*}}{T_{j}(t)}\Bigg)\Bigg(s_{j}-d_{j}T_{j}(t)-\frac{k_{j}T_{j}(t)V(t)}{1+\alpha_{j} V(t)}\Bigg)\\ +\sum\limits_{j=2,\cdots,n}\int_{0}^{\infty} \beta_{j}(a)i_{j}^{*}(a)\frac{\partial}{\partial t}\Bigg(\frac{i_{j}(a,t)}{i_{j}^{*}(a)}-1-\ln \frac{i_{j}(a,t)}{i_{j}^{*}(a)}\Bigg)da\\ +\Bigg(1-\frac{V^{*}}{V(t)}\Bigg)\Bigg(\sum\limits_{j=1}^{n}\int_{0}^{\infty}p_{j}(a)i_{j}(a,t)da-cV(t)\Bigg).$

It follows from $s_{j}=d_{j}T_{j}^{*}+\frac{k_{j}T_{j}^{*}V^{*}}{1+\alpha_{j} V^{*}}$ and $\sum\limits_{j=1}^{n}\frac{N_{j}k_{j}T_{j}^{*}}{1+\alpha_{j} V^{*}}=c$ that

$\frac{dW(t)}{dt}=term ①+term ②+term ③,$

where

$term ① = N_{1}\Bigg[\Big(1-\frac{T_{1}^{*}}{T_{1}(t)}\Big)d_{1}\Big(T_{1}^{*}-T_{1}(t)\Big)+\frac{k_{1}T_{1}^{*}V^{*}}{1+\alpha_{1} V^{*}} \\ -\frac{k_{1}T_{1}(t)V(t)}{1+\alpha_{1} V(t)}+\frac{k_{1}T_{1}^{*}V(t)}{1+\alpha_{1} V} -\frac{T_{1}^{*}}{T_{1}(t)}\frac{k_{1}T_{1}^{*}V^{*}}{1+\alpha_{1} V^{*}}\Bigg] \\ +\int_{0}^{\infty} \beta_{1}(a)\Big(1-\frac{i_{1}^{*}(a)}{i_{1}(a,t)}\Big)\frac{\partial i_{1}(a,t)}{\partial t}da,$

$term ② = \sum\limits_{j=2,\cdots,n}N_{j}\Bigg[\Big(1-\frac{T_{j}^{*}}{T_{j}(t)}\Big)d_{j}\Big(T_{j}^{*}-T_{j}(t)\Big) +\frac{k_{j}T_{j}^{*}V^{*}}{1+\alpha_{j} V^{*}}\\ -\frac{k_{j}T_{j}(t)V(t)}{1+\alpha_{j} V(t)}+\frac{k_{j}T_{j}^{*}V(t)} {1+\alpha_{j} V}-\frac{T_{j}^{*}}{T_{j}(t)}\frac{k_{j}T_{j}^{*}V^{*}}{1+\alpha_{j} V^{*}}\Bigg]\\ +\sum\limits_{j=2,\cdots,n}\int_{0}^{\infty} \beta_{j}(a)\Big(1-\frac{i_{j}^{*}(a)}{i_{j}(a,t)}\Big)\frac{\partial i_{j}(a,t)}{\partial t}da\\ +\int_{0}^{\infty}p_{1}(a)i_{1}(a,t)da-\frac{V^{*}}{V(t)}\int_{0}^{\infty}p_{1}(a)i_{1}(a,t)da,$

$term ③ = \frac{N_{1}k_{1}T_{1}^{*}V^{*}}{1+\alpha_{1} V^{*}}-\frac{N_{1}k_{1}T_{1}^{*}V(t)}{1+\alpha_{1} V^{*}} +\sum\limits_{j=2,\cdots,n}\int_{0}^{\infty}p_{j}(a)i_{j}(a,t)da\\ -\sum\limits_{j=2,\cdots,n}\frac{V^{*}}{V(t)}\int_{0}^{\infty}p_{j}(a)i_{j}(a,t)da\\ +\sum\limits_{j=2,\cdots,n}\Bigg[\frac{N_{j}}{N_{1}}\frac{k_{j}T_{j}^{*}V^{*}}{1+\alpha_{j} V^{*}} -\frac{N_{j}}{N_{1}}\frac{k_{j}T_{j}^{*}V(t)}{1+\alpha_{j} V^{*}}\Bigg].$

We further obtain that

$\frac{dW(t)}{dt} = -d_{1}N_{1}T_{1}^{*}\Big(2-\frac{T_{1}(t)}{T_{1}^{*}}-\frac{T_{1}^{*}}{T_{1}(t)}\Big)\\ -\beta_{1}(a)i_{1}^{*}(a)\Big(\frac{i_{1}(a,t)}{i_{1}^{*}(a)}-1-\ln \frac{i_{1}(a,t)}{i_{1}^{*}(a)}\Big)|_{a=\infty}\\ -\frac{N_{1}k_{1}T_{1}^{*}V^{*}}{1+\alpha_{1} V^{*}}\Bigg[\frac{T_{1}^{*}}{T_{1}(t)}-1+\frac{V(t)}{V^{*}}-\frac{V(t)}{V^{*}}\frac{1+\alpha_{1} V^{*}}{1+\alpha_{1} V(t)} \\ +\ln \frac{T_{1}(t)V(t)(1+\alpha_{1} V^{*})}{T_{1}^{*}V^{*}(1+\alpha_{1} V(t))}\Bigg]\\ -\int_{0}^{\infty}p_{1}(a)i_{1}^{*}(a)\Big(\frac{V^{*}i_{1}(a,t)}{Vi_{1}^{*}(a)}-1-\ln \frac{i_{1}(a,t)}{i_{1}^{*}(a)}\Big)da\\ -\sum\limits_{j=2,\cdots,n}N_{j}d_{j}T_{j}^{*}\Big(2-\frac{T_{j}(t)}{T_{j}^{*}}-\frac{T_{j}^{*}}{T_{j}(t)}\Big)\\ -\sum\limits_{j=2,\cdots,n}\beta_{j}(a)i_{j}^{*}(a)\Big(\frac{i_{j}(a,t)}{i_{j}^{*}(a)}-1-\ln \frac{i_{j}(a,t)}{i_{j}^{*}(a)}\Big)|_{a=\infty}\\ -\sum\limits_{j=2,\cdots,n}\frac{N_{j}k_{j}T_{j}^{*}V^{*}}{1+\alpha_{1} V^{*}}\Bigg[\frac{T_{j}^{*}}{T_{j}(t)}-1+\frac{V(t)}{V^{*}} \\ -\frac{V(t)}{V^{*}}\frac{1+\alpha_{j} V^{*}}{1+\alpha_{j} V(t)}+\ln \frac{T_{j}(t)V(t)(1+\alpha_{j} V^{*})}{T_{j}^{*}V^{*}(1+\alpha_{j} V(t))}\Bigg]\\ -\sum\limits_{j=2,\cdots,n}\int_{0}^{\infty}p_{j}(a)i_{j}^{*}(a)\Bigg(\frac{V^{*}i_{j}(a,t)}{Vi_{j}^{*}(a)}-1-\ln \frac{i_{j}(a,t)}{i_{j}^{*}(a)}\Bigg)da.$

Because

$1+\frac{V(t)}{V^{*}}-\frac{1+\alpha_{j} V(t)}{1+\alpha_{j} V^{*}}-\frac{V(t)}{V^{*}}\frac{1+\alpha_{j} V^{*}}{1+\alpha_{j} V(t)}=\frac{(V(t)-V^{*})^{2}}{V^{*}(1+\alpha_{j} V^{*})(1+\alpha_{j} V(t))}$

and

$\ln \frac{T_{j}^{*}V^{*}(1+\alpha_{j} V(t))}{T_{j}(t)V(t)(1+\alpha_{j} V^{*})}+\ln \frac{i_{j}(a,t)}{i_{j}^{*}(a)}=\ln\frac{T_{j}^{*}}{T_{j}(t)}+\ln \frac{1+\alpha_{j} V(t)}{1+\alpha_{j} V^{*}}+\ln \frac{i_{j}(a,t)V^{*}}{i_{j}^{*}(a)V(t)},$

we have

$\frac{dW(t)}{dt} \\ = -d_{1}N_{1}T_{1}^{*}\Bigg(2-\frac{T_{1}(t)}{T_{1}^{*}}-\frac{T_{1}^{*}}{T_{1}(t)}\Bigg)\\ -\beta_{1}(a)i_{1}^{*}(a)\Bigg(\frac{i_{1}(a,t)}{i_{1}^{*}(a)}-1-\ln \frac{i_{1}(a,t)}{i_{1}^{*}(a)}\Bigg)|_{a=\infty}\\ -\frac{N_{1}k_{1}T_{1}^{*}V^{*}}{1+\alpha_{1} V^{*}}\Bigg[\Bigg(\frac{T_{1}^{*}}{T_{1}(t)}-1-\ln\frac{T_{1}^{*}}{T_{1}(t)}\Bigg)+\Bigg(\frac{1+\alpha_{1} V(t)}{1+\alpha_{1} V^{*}}-1-\ln \frac{1+\alpha_{1} V(t)}{1+\alpha_{1} V^{*}}\Bigg)\\ +\frac{(V(t)-V^{*})^{2}}{V^{*}(1+\alpha_{1} V^{*})(1+\alpha_{1} V(t))}\Bigg]\\ -\int_{0}^{\infty}p_{1}(a)i_{1}^{*}(a)\Bigg(\frac{V^{*}i_{1}(a,t)}{V(t)i_{1}^{*}(a)}-1-\ln \frac{V^{*}i_{1}(a,t)}{V(t)i_{1}^{*}(a)}\Bigg)da\\ -\sum\limits_{j=2,\cdots,n}N_{j}d_{j}T_{j}^{*}\Bigg(2-\frac{T_{j}(t)}{T_{j}^{*}}-\frac{T_{j}^{*}}{T_{j}(t)}\Bigg)\\ -\sum\limits_{j=2,\cdots,n}\beta_{j}(a)i_{j}^{*}(a)\Bigg(\frac{i_{j}(a,t)}{i_{j}^{*}(a)}-1-\ln \frac{i_{j}(a,t)}{i_{j}^{*}(a)}\Bigg)|_{a=\infty}\\ -\sum\limits_{j=2,\cdots,n}\frac{N_{j}k_{j}T_{j}^{*}V^{*}}{1+\alpha_{1} V^{*}}\Bigg[\Bigg(\frac{T_{j}^{*}}{T_{j}(t)}-1-\ln\frac{T_{j}^{*}}{T_{j}(t)}\Bigg) \\ +\Bigg(\frac{1+\alpha_{j} V(t)}{1+\alpha_{j} V^{*}}-1-\ln \frac{1+\alpha_{j} V(t)}{1+\alpha_{j} V^{*}}\Bigg)\\ +\frac{N_{j}k_{j}T_{j}^{*}}{1+\alpha_{1} V^{*}}\frac{(V(t)-V^{*})^{2}}{(1+\alpha_{j} V^{*})(1+\alpha_{j} V(t)\Big)}\Bigg]\\ -\sum\limits_{j=2,\cdots,n}\int_{0}^{\infty}p_{j}(a)i_{j}^{*}(a)\Bigg(\frac{V^{*}i_{j}(a,t)}{V(t)i_{j}^{*}(a)}-1 -\ln \frac{V^{*}i_{j}(a,t)}{V(t)i_{j}^{*}(a)}\Bigg)da.$

Because the arithmetical mean is greater than or equal to the geometrical mean, we know that $2-\frac{T_{j}(t)}{T_{j}^{*}}-\frac{T_{j}^{*}}{T_{j}(t)}$ is less than or equal to zero. It follows that $\frac{dW}{dt}=0$ if and only if $T_{j}(t)=T_{j}^{*},~i_{j}(a,t)=i_{j}^{*}(a),$ and $V(t)=V^{*}$. Hence, every solution of system (2.1) converges to $E^{*}$, which means that the infected steady state is globally asymptotically stable from the LaSalle's invariance principle.

---

5. Related models

In this section, we show that under some special cases the age of infection model (2.1)$\sim$(2.3) can be reduced to an ordinary differential equation (ODE) system with or without time delays.

Case Ⅰ. Assume that $\delta_{j}(a)=\delta_{j},$ $p_{j}(a)=p_{j},$ where $\delta_{j}$ and $p_{j}$ are positive constants. Model (2.1) can be rewritten as

$\left\{\begin {array}{l l} \frac{d T_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}, \\ \frac{\partial i_{j}(a,t)}{\partial t}+\frac{\partial i_{j}(a,t)}{\partial a}=-\delta_{j}i_{j}(a,t), \\ i_{j}(0,t)=\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}\\ \frac{d V(t)}{d t}= \sum\limits_{j=1}^{n} p_{j} \int_{0}^{\infty}i_{j}(a,t)da-cV(t). \end {array} \right.$

(5.1)

We further assume that in this case $t$ is larger than all possible infection ages, and consequently $i_{j}(a,t)$ is expressed by the first half of (2.5), that is,

$i_{j}(a,t)=i_{j}(0,t-a)\exp\Bigg(- \int_{0}^{a}\delta_{j}(\theta)d\theta\Bigg)=i_{j}(0,t-a)e^{-\delta_{j}a} ~~\mbox{for} ~~a <t.$

Setting

$I_{j}(t)=\int_{0}^{\infty}i_{j}(a,t)da,$

(5.2)

which represents the total number of infected cells in the $j$th class at time $t,$ we have

$\begin {array}{l l} \frac{dI_{j}(t)}{dt}= \int_{0}^{\infty}\frac{\partial i_{j}(a,t)}{\partial t}da= -\int_{0}^{\infty}\Bigg(\frac{\partial i_{j}(a,t)}{\partial a}+\delta_{j}i_{j}(a,t)\Bigg)da \\= -i_{j}(a,t)|^{a=\infty}_{a=0}-\delta_{j}\int_{0}^{\infty}i_{j}(a,t)da. \end {array}$

From the boundary condition, one can see that $i_{j}(0,t)=\frac{k_{j}V(t)T_{j}(t)}{1+\alpha_{j} V(t)}.$ Since $T_{j}(t)$ and $V(t)$ are bounded on $[0,+\infty),$ it can be concluded that

$\lim\limits_{a\rightarrow+\infty}i_{j}(a,t)=\lim\limits_{a\rightarrow+\infty}i_{j}(0,t-a)e^{-\delta_{j}a}= \lim\limits_{a\rightarrow+\infty}e^{-\delta_{j}a}\frac{k_{j}V(t)T_{j}(t)}{1+\alpha_{j} V(t)}=0.$

Hence,

$\frac{dI_{j}(t)}{dt}=\frac{k_{j}V(t)T_{j}(t)}{1+\alpha_{j} V(t)}-\delta_{j}I(t).$

Thus, model (5.1) is equivalent to a standard ODE model with multiple target cell populations, given by

$\left\{\begin {array}{l l} \frac{dT_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}, \\ \frac{d I_{j}(a,t)}{d t}=\frac{k_{j}V(t)T_{j}(t)}{1+\alpha_{j} V(t)}-\delta_{j}I_{j}(a,t), \\ \frac{d{V(t)}}{d t}= \sum\limits_{j=1}^{n} p_{j}I_{j}(t)-cV(t). \end {array} \right.$

(5.3)

Note that in the above ODE model, the viral burst size is $N_{j}=p_{j}/\delta_{j}.$

Case Ⅱ. Assume that it takes time $\tau$ for virus to enter into the target cell, and that there is a time delay $\omega(\omega >\tau)$ between cell infection and viral production. The death rate of infected cells and the viral production rate become

$\delta_{j}(a)=\left\{\begin {array}{l l}\delta_{j}', ~~ a\geq \tau,\\ 0, ~~~0\leq a <\tau, \end {array} \right.$

(5.4)

and

$p_{j}(a)=\left\{\begin {array}{l l}p_{j}', ~~a\geq \omega,\\ 0, ~~~0\leq a <\omega. \end {array} \right.$

(5.5)

Note that $i_{j}(\tau, t) = e^{-\tau\delta_{j}'}i(0, t-\tau).$ We have

$\frac{dI_{j}(t)}{dt}=e^{-\tau\delta_{j}'}\frac{k_{j}V(t-\tau)T_{j}(t-\tau)}{1+\alpha_{j} V(t-\tau)}-\delta_{j}'I_{j}(t).$

Using function (5.5) we have

$\int_{0}^{\infty}p_{j}(a)i_{j}(a,t)da = p_{j}'\int_{\omega}^{\infty}i_{j}(a,t)da \\ = p_{j}'\int_{0}^{\infty}i_{j}(a+\omega,t)da\\ = p_{j}'\int_{0}^{\infty}e^{-\omega\delta_{j}'}i_{j}(a,t-\omega)da \\ = e^{-\omega\delta_j'}p_{j}'I_{j}(t-\omega).$

Hence, model (2.1) with the assumptions (5.4) and (5.5) can be reformulated equivalently as the following delay differential equation (DDE) system:

$\left\{\begin {array}{l l} \frac{d T_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}, \\ \frac{d I_{j}(a,t)}{d t}=e^{-\tau\delta_{j}'}\frac{k_{j}V(t-\tau)T_{j}(t-\tau)}{1+\alpha_{j} V(t-\tau)}-\delta_{j}'I_{1}(a,t), \\ \frac{dV(t)}{d t}= \sum\limits_{j=1}^{n} e^{-\omega\delta_{j}'}p_{j}'I_{j}(t-\omega)-cV(t). \end {array} \right.$

(5.6)

In the above DDE model, $\tau$ and $\omega$ are two intracellular delays describing the time required for virus to enter a target cell and for the infected cell to produce virus, respectively.

A special case of model (5.6) is with the assumption $p_{j}(a)=\bar{p}_{j}=\mbox{constant}$. In this case, system (5.6) becomes

$\left\{\begin {array}{l l} \frac{dT_{j}(t)}{d t}=s_{j}-d_{j}T_{j}(t)-\frac{k_{j} T_{j}(t)V(t)}{1+\alpha_{j} V(t)}, \\ \frac{d I_{j}(a,t)}{d t}=e^{-\tau\delta_{j}'}\frac{k_{j}V(t-\tau)T_{j}(t-\tau)}{1+\alpha_{j} V(t-\tau)}-\delta_{j}'I_{j}(a,t), \\ \frac{dV(t)}{d t}= \sum\limits_{j=1}^{n}\bar{p}_{j}I_{j}(t)-cV(t), \end {array} \right.$

(5.7)

which contains only one time delay.

---

6. Summary and Discussion

Viral infection of different classes of target cells is important in understanding the virus dynamics within infected individuals. HIV infection in macrophages may contribute to the early-stage viral transmission, persistence, and virus dissemination throughout the body [6]. Macrophages are resistant to the cytopathic effect of HIV and can produce virus for a longer period of time [14]. Thus, the production of virus by infected macrophages may explain the viral load explosion in the advanced stage of HIV infection. A recent model included the infection of macrophages to explain the three stages of HIV infection [11]. Damage to monocyte/macrophage lineage cells, although less obvious, provides the information to predict the onset of opportunistic infections and progression to AIDS [15]. Infection of peripheral blood mononuclear cells by HCV may explain the high levels of immunological disorders found in chronically infected HCV patients [7].

In this note, we developed and studied a within-host viral dynamic model including multiple populations of target cells and the age of viral infection. The global asymptotic properties for the model are obtained. When the basic reproductive number is below unity, the infection is predicted to die out. When the basic reproductive number exceeds unity, there exists a unique infected steady state which is globally asymptotically stable. This means that the virus is able to establish the infection within the host. With some assumptions the age-structured model can be reduced to an ODE or DDE system. This model can also be extended to study virus dynamics with multiple compartments or coinfection by multiple types/strains of viruses.

The model with multiple populations of target cells can be used to evaluate the relative contribution of viral production from different compartments. This can improve the understanding of viral evolution and disease progression. However, very limited (spatial) data are available for each cell population or compartment. Thus, it is challenging to estimate parameters and verify models. Another limitation of the model is that it cannot account for the long-term HIV dynamics observed in patients on prolonged antiretroviral therapy. Viral infection is predicted to die out if the basic reproductive number is below 1. However, virus can persist for a prolonged period of time even in patients under long-term antiretroviral therapy. HIV latency can be incorporated into the model to study the long-term virus dynamics under therapy.

---