The role of long-lived plasma cells in viral clearance

ABSTRACT

The adaptive immune system has two types of plasma cells (PC), long-lived plasma cells (LLPC) and short-lived plasma cells (SLPC), that differ in their lifespan. In this paper, we propose that LLPC is crucial to the clearance of viral particles in addition to reducing the viral basic reproduction number in secondary infections. We use a sequence of within-host mathematical models to show that, CD8 T cells, SLPC and memory B cells cannot achieve full viral clearance, and the viral load will reach a low positive equilibrium level because of a continuous replenishment of target cells. However, the presence of LLPC is crucial for viral clearance.

KEYWORDS:

• Viral clearance
• adaptive immune system
• basic reproduction number
• within-host mathematical model
• secondary infection

MATHEMATICS SUBJECT CLASSIFICATION (2020):

• 92C30

1. Introduction

In human immune responses to viral infections, plasma cells (PC) play an important role, as they produce antibodies that neutralize the virus [1]. PC are generally divided into long-lived plasma cells (LLPC) and short-lived plasma cells (SLPC) based on their lifespan. SLPC typically live for a few days then die [1] while LLPC may have a lifespan of years or decades [2–4]. It is generally believed that SLPC produce a large number of antibodies to clear the virus during an infection [5], while LLPC is generally thought to maintain constant antibody levels to prevent future infections [1,6,7]. In this paper, we use mathematical models to demonstrate that viruses reach a positive equilibrium in the absence of LLPC due to the short lifespan of SLPC and the target cells supplementation, thus virus clearance can only be achieved by LLPC.

Virus can trigger both innate and adaptive immune responses. Innate immune responses occur in the early stages of infection and produce cytokines that are nonspecific to the virus to regulate cell functions (e.g. putting target cells into quiescence or signal natural killer cells) [8,9]. Adaptive immune responses are virus-specific, and include a wide variety of immune responses that mainly produce two functions: CD8 T cells are activated, proliferated and differentiate into effector CD8 T cells to kill infected cells, and B cells are activated, proliferate and differentiate into SLPC and LLPC to produce antibodies, with the help of helper T cells which are activated, proliferate and differentiate from CD4 T cells [10]. During the process, memory immune cells (e.g. memory B cells and memory helper T cells) are produced, which reactive to secondary infections much quicker than naïve immune cells [11,12]. We will use multiple mathematical models to show that CD8 T cells, SLPC or memory B cells cannot successfully clear viruses without LLPC.

Viral dynamics models have been extensively applied to the investigation of the antiviral mechanism of immunity against a range of pathogens such as HIV and influenza. Simple within-host models include the uninfected target cells, infected target cells and free virus [13]. More complex modelling that incorporated antibodies [14] or CD8 T cells [15] were also considered without including full regulatory relationships of the adaptive immune system. Lee et al. [10] proposed a more complete within-host model containing dendritic cells, CD4 T cells, B cells, and PC et al. that quantifies the relations between viral replication and adaptive immunity. Ada et al. [16] investigated a within-host model that considers the memory of the adaptive immune responses. We adapt these models to study the effect of the immune responses of B cells and PC, naïve CD8 T cells and effector CD8 T cells, and memory B cells, respectively.

In Section 2, we build a simple host-virus model that includes healthy target cells, infected target cells, virus, B cells and PC, and show that the model cannot achieve viral clearance unless the B cells or PC have a very long lifespan. In Sections 3 and 4, we extend the simple model to include memory B cells and LLPC, respectively, and show that only LLPC eliminate the virus. In Section 5, we extend the simple model to consider CD8 T cells, and show that CD8 T cells cannot eliminate the virus without LLPC. Concluding remarks are given in Section 6.

2. A within-host model with PC

In this section, we develop a target-cell model that considers the dynamics of B cells, PC, and antiviral antibodies, and study its dynamics. We show that this model has a positive equilibrium viral load.

2.1. Model formulation

Our model describes the interaction between healthy target cells (H), infected target cells (I), virus (V), naïve B cells (B), activated B cells ($B_A$), PC (P) and antiviral antibody (A). This model is a simplification of the Lee et al. model [10]. In our model, healthy target cells (H) are recruited at a rate ${\lambda }_H$. Each healthy target cell dies at a rate ${\delta }_H$, and becomes infected by free viral particles (V) at a rate τ. Each infected target cell dies at a rate ${\delta }_I$ and releases n free viral particles upon death. A free viral particle is neutralized by antibody (A) at a rate $k_V$ and is cleared at a rate c. Naïve B cells are recruited at a rate ${\lambda }_B$. A naïve B cell dies at a rate ${\delta }_B$, and is activated by antigen (here represented by the viral load V) and differentiates into activated B cells ($B_A$) at a rate s. $B_A$ then proliferate and differentiate into PC at a rate μ. A plasma cell produces antibodies at a rate ${\pi }_P$ and dies at a rate ${\delta }_P$. An antibody particle is cleared at a rate ${\delta }_A$. Our model is given by the following system of differential equations: (1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) (1b) $\begin{array}{rl}{I}^{\mathrm{\prime }}& =\mathrm{\tau HV}-{\delta }_{I}I,\end{array}$(1b) (1c) $\begin{array}{rl}{V}^{\mathrm{\prime }}& =n{\delta }_{I}I-\mathrm{cV}-k_V\mathrm{AV},\end{array}$(1c) (1d) $\begin{array}{rl}{B}^{\mathrm{\prime }}& ={\lambda }_B-\mathrm{sVB}-{\delta }_{B}B,\end{array}$(1d) (1e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(1e) (1f) $\begin{array}{rl}{P}^{\mathrm{\prime }}& =\mu B_A-{\delta }_{P}P,\end{array}$(1f) (1g) $\begin{array}{rl}{A}^{\mathrm{\prime }}& ={\pi }_{P}P-{\delta }_{A}A.\end{array}$(1g) The parameters are listed in Table 1. We assume that all parameters are positive. Infected target cell dies at a rate ${\delta }_I>{\delta }_H$ due to both virus-induced and immune-system-induced toxicity. The initial values of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) are non-negative. A diagram of this model is shown in Figure 1.

Figure 1. The flow chart for the within-host virus Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ).

Table 1. Parameter value.

Variable	Definitions	Initial value (ref [10])
H	Healthy target cells (EID${}_{50}$/ml)	$5\times {10}^5$
I	Infected target cells (EID${}_{50}$/ml)	0
V	Virus (EID${}_{50}$/ml)	10
B	Naïve B cells (EID${}_{50}$/ml)	${10}^3$
$B_A$	Activated B cells (EID${}_{50}$/ml)	0
$B_m$	Memory B cells (EID${}_{50}$/ml)	0
P	Plasma cells (PC) (EID${}_{50}$/ml)	0
$P_s$	Short-lived plasma cells (SLPC) (EID${}_{50}$/ml)	0
$P_{\ell }$	Long-lived plasma cells (LLPC) (EID${}_{50}$/ml)	0
A	Antiviral antibody (EID${}_{50}$/ml)	0
T	Naïve CD8 T cells (EID${}_{50}$/ml)	${10}^3$
$T_c$	Effector CD8 T cells (EID${}_{50}$/ml)	0
Parameter	Definitions	Value
τ	Infection rate of healthy target cells [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	$7\times {10}^{-5}$
${\lambda }_H$	Recruitment rate of healthy target cells [day${}^{-1}$(EID${}_{50}$/ml)]	$5\times {10}^2$
${\delta }_H$	Death rate of healthy target cells (day${}^{-1}$)	${10}^{-3}$
m	Killing rate of infected target cells by effector CD8 T cells [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	$1.8\times {10}^{-4}$
${\delta }_I$	Death rate of infected target cells (day${}^{-1}$)	1.2
n	Number of viral particles released per infected cell	1.58
c	Viral clearance rate (day${}^{-1}$)	1
$k_V$	Antibody neutralization rate [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	$4\times {10}^{-3}$
${\lambda }_B$	Recruitment of naïve B cells [day${}^{-1}$(EID${}_{50}$/ml)]	2
s	Activation rate of naïve B cells [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	$3\times {10}^{-4}$
μ	Differentiation rate of activated B cells into SLPC (day${}^{-1}$)	${10}^{-3}$
${\delta }_B$	Death rate of naïve B cells (day${}^{-1}$)	$2\times {10}^{-3}$
${\delta }_P$	Death rate of PC (day${}^{-1}$)	0.1
${\pi }_P$	Antibody secretion rate of a plasma cell (day${}^{-1}$)	0.06
${\delta }_A$	Clearance rate of antibodies (day${}^{-1}$)	0.04
f	Fraction of B cells becoming plasma cells in Model (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) )	0.9, 0.1 (Estimated)
$s_M$	Activation rate of memory B cells in Model (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) ) [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	$1.5\times {10}^{-3}$ (Estimated)
φ	The fraction of proliferated B cells becoming SLPC in Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) )	0.99, 0.01
${\delta }_s$	SLPC death rate in Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) (day${}^{-1}$)	0.1
${\delta }_{\ell }$	LLPC death rate in Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) (day${}^{-1}$)	$3\times {10}^{-2}$
${\pi }_s$	Antibody secretion rate by a short-lived plasma cell in Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) (day${}^{-1}$)	0.06
${\pi }_{\ell }$	Antibody secretion rate by a long-lived plasma cell in Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) (day${}^{-1}$)	0.8
${\lambda }_T$	Naïve CD8 T cells production rate in Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	2
σ	Activation rate of naïve CD8 T cells in Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) [day${}^{-1}$(EID${}_{50}$/ml)${}^{-1}$]	$3\times {10}^{-2}$
${\delta }_T$	Death rate of naïve CD8 T cells in Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) (day${}^{-1}$)	$2\times {10}^{-3}$
${\delta }_{T_c}$	Death rate of effector CD8 T cells in Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) (day${}^{-1}$)	0.75

To understand how well our model approximates the behaviour of a more realistic model, we compare the numerical solutions of our Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) with those of the Lee et al. [10] model, which includes other immune system components such as antigen-presenting cells, T helper cells and CD8 T cells with Michaelis–Menten interaction terms. We use their parameter values and initial conditions, and vary ${\delta }_p$ and $k_V$ to study the effect of short-lived plasma cells and antibody neutralization for both models. The results are shown in Figure 2. The figure shows that the two models demonstrate similar qualitative behaviour. In fact, using the same values ${\delta }_p=0.1$ and $k_V=0.004$ as in [10], the solutions of the two models are almost identical. We thus base our study on our simpler Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ).

Figure 2. Comparison of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) and Lee et al. model [10]. The values of ${\delta }_p$ and $k_V$ are shown in column and row titles. The other parameter values and the initial conditions are taken from [10], also listed in Table 1.

2.2. Basic reproduction number

In this section, we compute the basic reproduction number ${\mathcal{R}}_0$ of the viral infection and study the equilibria for Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). First, note that the set ${\mathrm{\Omega }}_1:=\left\{\right(H,I,V,B,B_A,P,A)\in {\mathrm{\Re }}_{+}^7:H+I\le \frac{{\lambda }_H}{{\delta }_H},B+B_A\le \frac{{\lambda }_B}{min\{{\delta }_B,\mu \}}\}$is positively invariant for Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). We study the behaviour of the model in ${\mathrm{\Omega }}_1$.

We use the next generation matrix approach [17] to compute ${\mathcal{R}}_0$ of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). The disease-free equilibrium (DFE) is $(\frac{{\lambda }_H}{{\delta }_H},0,0,\frac{{\lambda }_B}{{\delta }_B},0,0,0)$. After linearizing (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) about the DFE, we use the infection-related terms to find the new-infection matrix F, and use the transition-related terms to find the transition matrix V $F=\left[\begin{array}{ccccc}0& \tau \frac{{\lambda }_H}{{\delta }_H}& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right],V=\left[\begin{array}{ccccc}{\delta }_I& 0& 0& 0& 0\\ -n{\delta }_I& c& 0& 0& 0\\ 0& -\frac{s{\lambda }_B}{{\delta }_B}& \mu & 0& 0\\ 0& 0& -\mu & {\delta }_P& 0\\ 0& 0& 0& -{\pi }_P& {\delta }_A\end{array}\right].$The next generation matrix is thus $F{V}^{-1}=\left[\begin{array}{ccccc}{\displaystyle \frac{\mathrm{n\tau }{\lambda }_H}{c{\delta }_H}}& {\displaystyle \frac{\tau {\lambda }_H}{c{\delta }_H}}& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right].$The basic reproduction number ${\mathcal{R}}_0$ is the spectral radius of the next generation matrix, i.e. (2) ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\mathrm{n\tau }{\lambda }_H}{c{\delta }_H}.$(2) It means a viral particle produces ${\mathcal{R}}_0$ new viral particles on average when all target cells are healthy. If ${\mathcal{R}}_0>1$, then the DFE is unstable, leading to viral infection. If ${\mathcal{R}}_0<1$, then the DFE is locally asymptotically stable, and thus a small amount of viral particles cannot lead to an infection.

2.3. Positive equilibrium

Theorem 2.1

If ${\mathcal{R}}_0>1$, then the system has a positive equilibrium value $V^{\ast }$, which is an increasing function of the PC death rate ${\delta }_P$, and ${\displaystyle \underset{{\delta }_P\to 0}{lim}{V}^{\ast }=0}$. In addition, ${\displaystyle \underset{{\delta }_B\to 0}{lim}{V}^{\ast }=0}$.

Proof.

We investigate the existence of positive equilibrium $(H^{\ast },I^{\ast },V^{\ast },B^{\ast },B_A^{\ast },P^{\ast },A^{\ast })$, which satisfies (3a) $\begin{array}{rl}{\lambda }_H-\tau H^{\ast }{V}^{\ast }-{\delta }_H{H}^{\ast }& =0,\end{array}$(3a) (3b) $\begin{array}{rl}\tau H^{\ast }{V}^{\ast }-{\delta }_I{I}^{\ast }& =0,\end{array}$(3b) (3c) $\begin{array}{rl}n{\delta }_I{I}^{\ast }-c{V}^{\ast }-k_V{A}^{\ast }{V}^{\ast }& =0,\end{array}$(3c) (3d) $\begin{array}{rl}{\lambda }_B-s{V}^{\ast }{B}^{\ast }-{\delta }_B{B}^{\ast }& =0,\end{array}$(3d) (3e) $\begin{array}{rl}s{V}^{\ast }{B}^{\ast }-\mu B_A^{\ast }& =0,\end{array}$(3e) (3f) $\begin{array}{rl}\mu B_A^{\ast }-{\delta }_P{P}^{\ast }& =0,\end{array}$(3f) (3g) $\begin{array}{rl}{\pi }_P{P}^{\ast }-{\delta }_A{A}^{\ast }& =0.\end{array}$(3g) Note that (3a(3a) $\begin{array}{rl}{\lambda }_H-\tau H^{\ast }{V}^{\ast }-{\delta }_H{H}^{\ast }& =0,\end{array}$(3a) ) gives $H^{\ast }=\frac{{\lambda }_H}{\tau V^{\ast }+{\delta }_H},$and (3b(3b) $\begin{array}{rl}\tau H^{\ast }{V}^{\ast }-{\delta }_I{I}^{\ast }& =0,\end{array}$(3b) ) gives $I^{\ast }=\frac{\tau H^{\ast }{V}^{\ast }}{{\delta }_I}=\frac{\tau {\lambda }_H{V}^{\ast }}{{\delta }_I(\tau V^{\ast }+{\delta }_H)}.$Similarly, from (3d(3d) $\begin{array}{rl}{\lambda }_B-s{V}^{\ast }{B}^{\ast }-{\delta }_B{B}^{\ast }& =0,\end{array}$(3d) ), (3e(3e) $\begin{array}{rl}s{V}^{\ast }{B}^{\ast }-\mu B_A^{\ast }& =0,\end{array}$(3e) ), (3f(3f) $\begin{array}{rl}\mu B_A^{\ast }-{\delta }_P{P}^{\ast }& =0,\end{array}$(3f) ) and (3g(3g) $\begin{array}{rl}{\pi }_P{P}^{\ast }-{\delta }_A{A}^{\ast }& =0.\end{array}$(3g) ), we have $\begin{array}{rl}{B}^{\ast }& =\frac{{\lambda }_B}{s{V}^{\ast }+{\delta }_B},\\ B_A^{\ast }& =\frac{s{\lambda }_B{V}^{\ast }}{\mu (s{V}^{\ast }+{\delta }_B)},\\ P^{\ast }& =\frac{s{\lambda }_B{V}^{\ast }}{{\delta }_P(s{V}^{\ast }+{\delta }_B)},\\ A^{\ast }& =\frac{{\pi }_{P}s{\lambda }_{B}V}{{\delta }_A{\delta }_P(s{V}^{\ast }+{\delta }_B)}.\end{array}$Substituting these relationships into (3c(3c) $\begin{array}{rl}n{\delta }_I{I}^{\ast }-c{V}^{\ast }-k_V{A}^{\ast }{V}^{\ast }& =0,\end{array}$(3c) ) gives (4) $A_1{V}^{\ast 2}+B_1{V}^{\ast }+C_1=0,$(4) where $\begin{array}{rl}{A}_1& =\mathrm{\tau s}(c{\delta }_A{\delta }_P+{\pi }_P{k}_V{\lambda }_B)>0,\\ B_1& =c{\delta }_P{\delta }_A\left[s{\delta }_H\right(1-{\mathcal{R}}_0)+{\delta }_B\tau ]+s{k}_V{\pi }_P{\delta }_H{\lambda }_B,\\ C_1& =c{\delta }_H{\delta }_A{\delta }_B{\delta }_P(1-{\mathcal{R}}_0).\end{array}$If ${\mathcal{R}}_0>1$, then $C_1<0$, and thus Equation (4(4) $A_1{V}^{\ast 2}+B_1{V}^{\ast }+C_1=0,$(4) ) has a single positive root, and thus, the system has a unique positive equilibrium.

To show that $V^{\ast }$ is an increasing function of ${\delta }_P$, we calculate $\mathrm{\partial }{V}^{\ast }/\mathrm{\partial }{\delta }_P$ from (4(4) $A_1{V}^{\ast 2}+B_1{V}^{\ast }+C_1=0,$(4) ), $\mathrm{\tau sc}{\delta }_A{V^{\ast }}^2+c{\delta }_A\left[s{\delta }_H\right(1-{\mathcal{R}}_0)+{\delta }_B\tau ]V^{\ast }+c{\delta }_H{\delta }_A{\delta }_B(1-{\mathcal{R}}_0)+(2A_1{V}^{\ast }+B_1)\frac{\mathrm{\partial }{V}^{\ast }}{\mathrm{\partial }{\delta }_P}=0.$The coefficient of the derivative $2A_1{V}^{\ast }+B_1>0$ (because this is the slope of the tangent line of the quadratic (4(4) $A_1{V}^{\ast 2}+B_1{V}^{\ast }+C_1=0,$(4) ) at the larger root $V^{\ast }$). In addition, using (4(4) $A_1{V}^{\ast 2}+B_1{V}^{\ast }+C_1=0,$(4) ), $\begin{array}{rl}& \mathrm{\tau sc}{\delta }_A{V^{\ast }}^2+c{\delta }_A\left[s{\delta }_H\right(1-{\mathcal{R}}_0)+{\delta }_B\tau ]V^{\ast }+c{\delta }_H{\delta }_A{\delta }_B(1-{\mathcal{R}}_0)\\ & =-\frac{1}{{\delta }_P}(\mathrm{\tau s}{\pi }_P{k}_V{\lambda }_B{V^{\ast }}^2+s{k}_V{\pi }_P{\delta }_H{\lambda }_B{V}^{\ast })<0.\end{array}$Thus $\frac{\mathrm{\partial }{V}^{\ast }}{\mathrm{\partial }{\delta }_P}>0.$We also see that, if ${\delta }_P\to 0$ or ${\delta }_B\to 0$, $C_1\to 0$. So ${\displaystyle \underset{{\delta }_P\to 0}{lim}{V}^{\ast }=0}$ and ${\displaystyle \underset{{\delta }_B\to 0}{lim}{V}^{\ast }=0}$.

2.4. Viral clearance

Naively, we may regard viral clearance as ${\displaystyle \underset{t\to 0}{lim}V\left(t\right)=0}$. In this sense, Theorem 2.1 shows that if the PC or the B cells have a very long lifespan, the virus can be cleared, suggesting that either LLPC or memory B cells may lead to viral clearance. However, $V\left(t\right)$ approaching 0 can be an overly strong condition for viral clearance. We can see in the bottom left panel of Figure 2 that, even though the model solutions remain positive and approach equilibrium, temporarily the viral count can become very small. Intuitively, when the viral load is less than one viral particle (e.g. 0.5 particles), even though $V\left(t\right)$ remains positive, the virus can be considered clear in the host. Like [10], we use ${\mathrm{EID}}_{50}$ as the unit of viral load, which is the dilution leading to $50\mathrm{\%}$ probability infection. Thus, we use $V\left(t\right)=1$ as the viral clearance threshold in this paper, i.e. if $V\left(t\right)<1$, we consider the virus cleared in the host.

In Figure 3, we present the value of ${\delta }_p$ and $k_V$ that makes the positive equilibrium $V^{\ast }=1$. Above the line, $V^{\ast }<1$ i.e. the virus if eventually cleared in the host; below the line, $V^{\ast }>1$, i.e. the virus cannot guarantee to be cleared (but may be temporarily cleared and reinfected). This figure shows that viral clearance needs small ${\delta }_p$ (long life span for plasma cells) or large $k_V$ (high antibody clearance rate). For a fixed antibody clearance rate, this figure emphasizes the vital role of LLPC in virus clearance. These agree with Theorem 2.1.

Figure 3. When positive equilibrium $V^{\ast }=1$, the relationship between ${\delta }_P$ and $k_V$ in Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). The other parameter values are listed in Table 1.

On the other hand, viral clearance may be temporary, because the host may be reinfected if $V^{\prime }\left(t\right)>0$. If $V\left(t\right)$ hits the threshold at $t_1$ (i.e. $V\left(t_1\right)=1$), then reaches a trough at $t_m$ and increases afterward, then $V^{\prime }\left(t\right)>0$ for $t>t_m$, that is, if the host is re-exposed to the virus, then they can be reinfected. This is true if $V\left(t\right)>1$ for some $t>t_m$, and especially if the equilibrium $V^{\ast }>1$. In these cases, the protection is temporary. The duration of primary infection is $[0,t_1]$, and the period of protection is $[t_1,t_m]$.

In Figure 4, we show the effects of ${\delta }_p$, $k_V$ on the period of primary infection, period of temporal protection and positive equilibrium ($V^{\ast }$), respectively. Figure 4(a) suggests that the first infection period is directly proportional to the values of ${\delta }_p$ and is inversely proportional to the values of $k_V$. Figure 4(b) illustrates a positive correlation between the antibody neutralization rate $k_V$ and the temporary viral clearance period. In contrast, there is a negative correlation between the plasma cell death rate ${\delta }_p$ and the temporary viral clearance period. Note that, below the line of the temporal protection period being 0 in panel (b), $V\left(t\right)>1$ for all $t>1$ and the period of primary infection becomes infinity, i.e. the virus cannot be cleared even temporarily in the host. Figure 4(c) reveals that positive equilibrium ($V^{\ast }$) is low when the antibody neutralization rate ($k_V$) is high, and the plasma cell death rate (${\delta }_p$) is less. This figure shows that increasing $k_V$ and decreasing ${\delta }_p$ values enhances the host's ability to control viral replication and contribute to viral clearance.

Figure 4. Contour plot analysis in the Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). (a) The dependence of the primary infection period on ${\delta }_p$ (plasma cell death rate) and $k_V$ (antibody neutralization rate). (b) The dependence of the temporary protection period on ${\delta }_p$ and $k_V$. (c) The dependence of the positive equilibrium ($V^{\ast }$) on ${\delta }_p$ and $k_V$. The other parameter values and the initial conditions are listed in Table 1.

Figure 5 shows the dynamics of viral load with varying parameters, initial values and parameter values are shown in Table 1. These two figures reflect that the positive equilibrium may lose stability and there exist stable positive periodic solutions around the equilibrium, which is possibly from a Hopf-bifurcation. Moreover, although increase $k_V$, the concentration levels of free viruses show oscillatory behaviour with different parameters. This implies that reinfection may still happen without LLPC even if $V^{\ast }<1$, and the virus may keep circulating in the host.

Figure 5. The viral load with different parameters in the Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). Parameter values and initial values are shown in Table 1. Panel (a) shows that, with a large antibody neutralization rate, the positive equilibrium may lose stability and a stable limit cycle emerges. Panel (b) show that, with an even larger $k_V$, the positive equilibrium (approximately the time average of the periodic solution) will be below the viral clearance threshold V = 1 (the red line). However, because of the oscillation, reinfection is still possible. (a) ${\delta }_P=0.1$, $k_V=0.4$. (b) ${\delta }_P=0.1$, $k_V=2$.

In the next section, we will extend Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) to include LLPC, and show that the presence of LLPC indeed leads to viral clearance.

3. A within-host model with LLPC

In this section, we split the PC in Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) into SLPC and LLPC, without considering memory B cells. We will show that, the viral load eventually approaches 0 after an infection.

3.1. Model formulation

Let $P_s$ and $P_{\ell }$ be the SLPC and LLPC counts. Assume that a fraction φ of proliferated $B_A$ cells becomes SLPC, and the fraction $1-\phi$ becomes LLPC. A short-lived plasma cell dies at a constant rate ${\delta }_s$, while the death rate of LLPC is negligible and is set to 0. A short-lived plasma cell or a long-lived plasma cell produces antiviral antibodies at a constant rate ${\pi }_s$ or ${\pi }_{\ell }$, respectively. The extended model is given by the following equations: (5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) (5b) $\begin{array}{rl}{I}^{\mathrm{\prime }}& =\mathrm{\tau HV}-{\delta }_{I}I,\end{array}$(5b) (5c) $\begin{array}{rl}{V}^{\mathrm{\prime }}& =n{\delta }_{I}I-\mathrm{cV}-k_V\mathrm{AV},\end{array}$(5c) (5d) $\begin{array}{rl}{B}^{\mathrm{\prime }}& ={\lambda }_B-\mathrm{sVB}-{\delta }_{B}B,\end{array}$(5d) (5e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(5e) (5f) $\begin{array}{rl}{P}_s^{\mathrm{\prime }}& =\mathrm{\phi \mu }{B}_A-{\delta }_s{P}_s,\end{array}$(5f) (5g) $\begin{array}{rl}{P}_{\ell }^{\mathrm{\prime }}& =(1-\phi )\mu B_A,\end{array}$(5g) (5h) $\begin{array}{rl}{A}^{\mathrm{\prime }}& ={\pi }_s{P}_s+{\pi }_{\ell }{P}_{\ell }-{\delta }_{A}A.\end{array}$(5h) The meaning of the parameters are listed in Table 1. All parameters are positive and the initial values for Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) are nonnegative. The flow chart of this model is given in Figure 6.

Figure 6. The flow chart for the within-host virus Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ).

3.2. Basic reproduction number

We again use the next-generation matrix approach to compute the basic reproduction number ${\mathcal{R}}_0$ of the viral infection. It is easy to prove that ${\mathrm{\Omega }}_3:=\left\{\right(H,I,V,B,B_A,P_s,P_{\ell },A)\in {\mathrm{\Re }}_{+}^8:H+I\le \frac{{\lambda }_H}{{\delta }_H},B+B_A\le \frac{{\lambda }_B}{min\{{\delta }_B,\mu \}}\}$is a positive invariant set of Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ).

The disease-free equilibrium (DFE) is $(\frac{{\lambda }_H}{{\delta }_H},0,0,\frac{{\lambda }_B}{{\delta }_B},0,0,P_{\ell 0},\frac{{\pi }_{\ell }}{{\delta }_A}{P}_{\ell 0})$. Linearize about Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) at the DFE, then pick the infection-related terms to find the new-infection matrix $F=\left[\begin{array}{cccccc}0& \tau \frac{{\lambda }_H}{{\delta }_H}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],$and the transition matrix $V=\left[\begin{array}{cccccc}{\delta }_I& 0& 0& 0& 0& 0\\ -n{\delta }_I& c+k_V\frac{{\pi }_{\ell }}{{\delta }_A}{P}_{\ell 0}& 0& 0& 0& 0\\ 0& -\frac{s{\lambda }_B}{{\delta }_B}& \mu & 0& 0& 0\\ 0& 0& -\mathrm{\phi \mu }& {\delta }_s& 0& 0\\ 0& 0& -(1-\phi )\mu & 0& 0& 0\\ 0& 0& 0& -{\pi }_s& -{\pi }_{\ell }& {\delta }_A\end{array}\right].$Note that V has a single 0 eigenvalue, corresponding to the family of DFE given by arbitrary $P_{\ell 0}$. This also gives a 0 eigenvalue in the Jacobian matrix. We can remove this variable in the calculation of ${\mathcal{R}}_0$. So the basic reproduction number ${\mathcal{R}}_0$ is (6) ${\mathcal{R}}_0=\frac{\mathrm{n\tau }{\lambda }_H}{(c+k_V\frac{{\pi }_{\ell }}{{\delta }_A}{P_{\ell }}_0){\delta }_H}.$(6) From (6(6) ${\mathcal{R}}_0=\frac{\mathrm{n\tau }{\lambda }_H}{(c+k_V\frac{{\pi }_{\ell }}{{\delta }_A}{P_{\ell }}_0){\delta }_H}.$(6) ), it's easy to see that ${\mathcal{R}}_0$ decreases with the initial value of LLPC increases. In the primary infection, we have $P_{\ell 0}=0$. The basic reproduction number is also given by (2(2) ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\mathrm{n\tau }{\lambda }_H}{c{\delta }_H}.$(2) ), i.e. it is the same as the basic reproduction number of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). On the other hand, a large initial number of LLPC may reduce ${\mathcal{R}}_0$ to below unity, preventing a secondary infection.

3.3. Primary and secondary infections

The following theorem shows that the virus is always cleared by using Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ).

Theorem 3.1

Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) has no positive equilibrium, and its solutions starting in the positively invariant set ${\mathrm{\Omega }}_3$ satisfy $I\left(\mathrm{\infty }\right)=0,V\left(\mathrm{\infty }\right)=0,B_A\left(\mathrm{\infty }\right)=0,P_s\left(\mathrm{\infty }\right)=0,P_{\ell }\left(\mathrm{\infty }\right)<\mathrm{\infty },A\left(\mathrm{\infty }\right)=\frac{{\pi }_{\ell }}{{\delta }_A}{P}_{\ell }^{\ast }.$

Proof.

We first prove that Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) has no positive equilibrium. Let $(H^{\ast },I^{\ast },V^{\ast },B^{\ast },B_A^{\ast },{P_s}^{\ast },{P_{\ell }}^{\ast },A^{\ast })$ be an equilibrium of Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ). Assume that $V^{\ast }>0$, then from (5g(5g) $\begin{array}{rl}{P}_{\ell }^{\mathrm{\prime }}& =(1-\phi )\mu B_A,\end{array}$(5g) ), we have $(1-\phi )\mu B_A^{\ast }=0,$and thus $B_A^{\ast }=0$. From (5e(5e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(5e) ), $s{V}^{\ast }{B}^{\ast }-\mu B_A^{\ast }=0,$and thus $B^{\ast }=0$. From (5d(5d) $\begin{array}{rl}{B}^{\mathrm{\prime }}& ={\lambda }_B-\mathrm{sVB}-{\delta }_{B}B,\end{array}$(5d) ), ${\lambda }_B-s{V}^{\ast }{B}^{\ast }-{\delta }_B{B}^{\ast }=0.$This is impossible because ${\lambda }_B>0$. This contradiction originates from the assumption that $V^{\ast }>0$. Thus, $V^{\ast }=0$, i.e. Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) does not have a positive equilibrium.

We then prove that $P_A\left(\mathrm{\infty }\right)$ exists. Equation (5g(5g) $\begin{array}{rl}{P}_{\ell }^{\mathrm{\prime }}& =(1-\phi )\mu B_A,\end{array}$(5g) ) shows that $P_{\ell }^{\prime }\ge 0$, and thus $P_{\ell }\left(t\right)$ is monotonically increasing. We will show that $P_{\ell }\left(t\right)$ is bounded by reduction to absurdity. If $P_{\ell }\left(t\right)$ is unbounded, that is, $P_{\ell }\left(\mathrm{\infty }\right)=\mathrm{\infty }$, and thus, from (5h(5h) $\begin{array}{rl}{A}^{\mathrm{\prime }}& ={\pi }_s{P}_s+{\pi }_{\ell }{P}_{\ell }-{\delta }_{A}A.\end{array}$(5h) ), $A\left(\mathrm{\infty }\right)=\mathrm{\infty }$. Because $H\left(t\right)$ is bounded in ${\mathrm{\Omega }}_3$, we get (7) $n{I}^{\prime }\left(t\right)+V^{\prime }\left(t\right)=\mathrm{n\tau H}\left(t\right)V\left(t\right)-\mathrm{cV}\left(t\right)-k_{v}V\left(t\right)A\left(t\right)<-\mathrm{CV}\left(t\right)$(7) for some $C>0$ and large enough t. From (7(7) $n{I}^{\prime }\left(t\right)+V^{\prime }\left(t\right)=\mathrm{n\tau H}\left(t\right)V\left(t\right)-\mathrm{cV}\left(t\right)-k_{v}V\left(t\right)A\left(t\right)<-\mathrm{CV}\left(t\right)$(7) ), integrating on both sides gives $0\le C{\int }_0^{+\mathrm{\infty }}V\left(t\right)\mathrm{d}t<\mathrm{\infty }.$Note that $B\left(t\right)$ is bounded in ${\mathrm{\Omega }}_3$, and thus, integrating (5e(5e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(5e) ) gives ${\int }_0^{+\mathrm{\infty }}{B}_A\left(t\right)\mathrm{d}t\le \frac{s{\lambda }_B}{\mu {\delta }_B}{\int }_0^{+\mathrm{\infty }}V\left(t\right)\mathrm{d}t<\mathrm{\infty }$Integrating (5g(5g) $\begin{array}{rl}{P}_{\ell }^{\mathrm{\prime }}& =(1-\phi )\mu B_A,\end{array}$(5g) ) gives $0\le P_{\ell }\left(\mathrm{\infty }\right)-P_{\ell }\left(0\right)={\int }_0^{+\mathrm{\infty }}(1-\phi )\mu B_A\left(t\right)\mathrm{d}t<\mathrm{\infty },$which contradicts the assumption that $\underset{t\to \mathrm{\infty }}{lim}{p}_l\left(t\right)\to \mathrm{\infty }$. Therefore, $P_{\ell }\left(\mathrm{\infty }\right)$ is bounded.

In addition, (8) ${\int }_0^{+\mathrm{\infty }}{B}_A\left(t\right)\mathrm{d}t<\mathrm{\infty }$(8) and thus $B_A\left(\mathrm{\infty }\right)=0$. Substituting $B_A\left(\mathrm{\infty }\right)$ into (5g(5g) $\begin{array}{rl}{P}_{\ell }^{\mathrm{\prime }}& =(1-\phi )\mu B_A,\end{array}$(5g) ), and solving the equation gives $P_s\left(\mathrm{\infty }\right)=0$. Let $P_{\ell }\left(\mathrm{\infty }\right)=P_{\ell }^{\ast }$, substituting $P_s\left(\mathrm{\infty }\right)$ and $P_{\ell }\left(\mathrm{\infty }\right)$ into (5h(5h) $\begin{array}{rl}{A}^{\mathrm{\prime }}& ={\pi }_s{P}_s+{\pi }_{\ell }{P}_{\ell }-{\delta }_{A}A.\end{array}$(5h) ) gives ${\displaystyle A\left(\mathrm{\infty }\right)=\frac{{\pi }_{\ell }}{{\delta }_A}{P}_{\ell }^{\ast }}$.

Also, using the expressions of (8(8) ${\int }_0^{+\mathrm{\infty }}{B}_A\left(t\right)\mathrm{d}t<\mathrm{\infty }$(8) ) and integrating (5e(5e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(5e) ), we get ${\int }_0^{\mathrm{\infty }}V\left(t\right)B\left(t\right)\mathrm{d}t<\mathrm{\infty }.$From (5d(5d) $\begin{array}{rl}{B}^{\mathrm{\prime }}& ={\lambda }_B-\mathrm{sVB}-{\delta }_{B}B,\end{array}$(5d) ) and (5e(5e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(5e) ), $(\mathrm{VB}{)}^{\prime }=V^{\prime }B+B^{\prime }V={\delta }_I\mathrm{IB}+{\lambda }_{B}V-\mathrm{cVB}-k_V\mathrm{AVB}-s{V}^{2}B-{\delta }_B\mathrm{VB}.$Note that H and B are bounded in ${\mathrm{\Omega }}_3$, and thus integrating on both sides gives, ${\delta }_I{\int }_0^{+\mathrm{\infty }}I\left(t\right)B\left(t\right)\mathrm{d}t+{\lambda }_B{\int }_0^{+\mathrm{\infty }}V\left(t\right)\mathrm{d}t<\mathrm{\infty }.$Thus, ${\int }_0^{+\mathrm{\infty }}V\left(t\right)\mathrm{d}t<\mathrm{\infty }.$Therefore, $V\left(\mathrm{\infty }\right)=0$. Again, since H is bounded in ${\mathrm{\Omega }}_3$, solving (5b(5b) $\begin{array}{rl}{I}^{\mathrm{\prime }}& =\mathrm{\tau HV}-{\delta }_{I}I,\end{array}$(5b) ) gives $I\left(\mathrm{\infty }\right)=0$.

Theorem 3.1 guarantees that the viral load eventually reaches 0. Compared with the results in Sections 2 and 3, the result shows that the presence of LLPC indeed leads to viral clearance. But the clearance is a slow process, as shown in Figure 7. This is because the solution approaches another disease-free equilibrium along a centre manifold. The outcome of the infection is dependent on the role of LLPC, which would drive the total viral load to the disease-free equilibrium.

Figure 7. The viral load in the Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ). The red line presents the viral clearance threshold V = 1. Panel (a) shows the viral dynamics with realistic parameter values (see Table 1), and Panel (b) shows the viral dynamics with a much larger antibody neutralization rate $k_V$. The parameter values and initial values are shown in Table 1. (a) $k_V=4\times {10}^{-3}$. (b) $k_V=0.4$.

A secondary infection occurs after the recovery of the primary infection. In a secondary infection, the initial count of LLPC is the final LLPC count in the primary infection, i.e. $P_{\ell 0}=P_{\ell }^{\ast }$. Then, for the secondary infection, the basic reproduction number of Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ) is ${\mathcal{R}}_{02}=\frac{\mathrm{n\tau }{\lambda }_H}{(c+k_V\frac{{\pi }_{\ell }}{{\delta }_A}{P}_{\ell }^{\ast }){\delta }_H}.$Numerical simulation in Figure 8, ${\mathcal{R}}_{02}$ increased in waves with τ, and ${\mathcal{R}}_{02}$ was associated with $P_{\ell }^{\ast }$, shows that ${\mathcal{R}}_{02}\le 1$. Interestingly, ${\mathcal{R}}_{02}$ does not depend monotonically on τ. That may be due to complex interactions between the viral dynamics and immune response. That is, LLPC also protect patients from reinfection by the same virus.

Figure 8. The curves $R_{02}$ are plotted as a function of τ in a secondary infection of Model (5(5a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(5a) ). The initial values are: $P_{\ell }\left(0\right)=P_{\ell }^{\ast }$, $A\left(0\right)=A^{\ast }$, where $P_{\ell }^{\ast }=P_{\ell }\left(\mathrm{\infty }\right)$, $A^{\ast }=A\left(\mathrm{\infty }\right)$ in the primary infection. Parameter values and other initial values are shown in Table 1.

In the next section, we will extend this model to include memory B cells, and show that the positive equilibrium still exists with memory B cells, i.e. memory B cells alone cannot lead to viral clearance.

4. A within-host model with memory B cells

In this section, we incorporate memory B cells into Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ), and show that the viral clearance cannot be achieved without LLPC.

4.1. Model formulation

We assume that a fraction f of the activated B cells $B_A$ becomes PC after proliferation, while the fraction 1−f becomes memory B cells (denoted as $B_M$). A memory B cell is activated again by antigen (here represented by V) and becomes an activated B cell at a rate $s_M$. We assume that $s_M>s$, because the reactivation is a faster process than the activation and proliferation of naïve B cells [18]. The memory B cells do not die. The extended model is given by the following equations: (9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) (9b) $\begin{array}{rl}{I}^{\mathrm{\prime }}& =\mathrm{\tau HV}-{\delta }_{I}I,\end{array}$(9b) (9c) $\begin{array}{rl}{V}^{\mathrm{\prime }}& =n{\delta }_{I}I-\mathrm{cV}-k_V\mathrm{AV},\end{array}$(9c) (9d) $\begin{array}{rl}{B}^{\mathrm{\prime }}& ={\lambda }_B-\mathrm{sVB}-{\delta }_{B}B,\end{array}$(9d) (9e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}+s_{M}V{B}_M-\mu B_A,\end{array}$(9e) (9f) $\begin{array}{rl}{P}^{\mathrm{\prime }}& =\mathrm{f\mu }{B}_A-{\delta }_{P}P,\end{array}$(9f) (9g) $\begin{array}{rl}{B}_M^{\mathrm{\prime }}& =(1-f)\mu B_A-s_{M}V{B}_M,\end{array}$(9g) (9h) $\begin{array}{rl}{A}^{\mathrm{\prime }}& ={\pi }_{P}P-{\delta }_{A}A.\end{array}$(9h) Similar to Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ), we assume that all parameters are positive, the initial values are non-negative and ${\delta }_I>{\delta }_H$. The parameters are listed in Table 1. The flow chart for this model is shown in Figure 9.

Figure 9. The flow chart for the within-host virus Model (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) ).

4.2. Basic reproduction number

Again, we use the next generation matrix approach to compute the basic reproduction number ${\mathcal{R}}_0$. First, note that the set ${\mathrm{\Omega }}_2:=\left\{\right(H,I,V,B,B_A,P,B_M,A)\in {\mathrm{\Re }}_{+}^8:H+I\le \frac{{\lambda }_H}{{\delta }_H},B\le \frac{{\lambda }_B}{{\delta }_B}\}$is positively invariant for Model (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) ).

Model (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) ) has a DFE: $(\frac{{\lambda }_H}{{\delta }_H},0,0,\frac{{\lambda }_B}{{\delta }_B},0,0,B_{M0},0)$, where $B_{M0}>0$ is the initial memory B cell count from previous infections. Linearizing (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) ) about the DFE gives the new infection matrix $F=\left[\begin{array}{cccccc}0& \tau \frac{{\lambda }_H}{{\delta }_H}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$and the transition matrix $V=\left[\begin{array}{cccccc}{\delta }_I& 0& 0& 0& 0& 0\\ -n{\delta }_I& c& 0& 0& 0& 0\\ 0& -\frac{s{\lambda }_B}{{\delta }_B}-s_M{B}_{M0}& \mu & 0& 0& 0\\ 0& 0& -\mathrm{\mu f}& {\delta }_P& 0& 0\\ 0& s_M{B}_{M0}& -\mu (1-f)& 0& 0& 0\\ 0& 0& 0& -{\pi }_P& 0& {\delta }_A\end{array}\right].$Note that V has a single 0 eigenvalue, corresponding to the family of DFE given by arbitrary $B_{M0}$. This also gives a 0 eigenvalue in the Jacobian matrix. We can remove this variable in the calculation of ${\mathcal{R}}_0$. The basic reproduction number ${\mathcal{R}}_0$ is also given by (2(2) ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\mathrm{n\tau }{\lambda }_H}{c{\delta }_H}.$(2) ), i.e. it is the same as the basic reproduction number of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). Thus, memory B cells do not change the basic reproduction number, because their effects on viral suppression come into play only after they are activated by antigens, and is a higher order term near the DFE.

4.3. Positive equilibrium

We investigate the existence of positive equilibrium $(H^{\ast },I^{\ast },V^{\ast },B^{\ast },{B_A}^{\ast },P^{\ast },{B_M}^{\ast },A^{\ast })$ of Model (9(9a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(9a) ), which satisfies (10a) $\begin{array}{rl}{\lambda }_H-\tau H^{\ast }{V}^{\ast }-{\delta }_H{H}^{\ast }& =0,\end{array}$(10a) (10b) $\begin{array}{rl}\tau H^{\ast }{V}^{\ast }-{\delta }_I{I}^{\ast }& =0,\end{array}$(10b) (10c) $\begin{array}{rl}n{\delta }_I{I}^{\ast }-c{V}^{\ast }-k_V{A}^{\ast }{V}^{\ast }& =0,\end{array}$(10c) (10d) $\begin{array}{rl}{\lambda }_B-s{V}^{\ast }{B}^{\ast }-{\delta }_B{B}^{\ast }& =0,\end{array}$(10d) (10e) $\begin{array}{rl}s{V}^{\ast }{B}^{\ast }+s_M{V}^{\ast }{B}_M^{\ast }-\mu B_A^{\ast }& =0,\end{array}$(10e) (10f) $\begin{array}{rl}\mathrm{f\mu }{B}_A^{\ast }-{\delta }_P{P}^{\ast }& =0,\end{array}$(10f) (10g) $\begin{array}{rl}\mathrm{f\mu }{B}_A^{\ast }-{\delta }_P{P}^{\ast }& =0,\end{array}$(10g) (10h) $\begin{array}{rl}(1-f)\mu B_A^{\ast }-s_M{V}^{\ast }{B}_M^{\ast }& =0,\end{array}$(10h) (10i) $\begin{array}{rl}{\pi }_P{P}^{\ast }-{\delta }_A{A}^{\ast }& =0.\end{array}$(10i) Note that (10a(10a) $\begin{array}{rl}{\lambda }_H-\tau H^{\ast }{V}^{\ast }-{\delta }_H{H}^{\ast }& =0,\end{array}$(10a) ) gives $H^{\ast }=\frac{{\lambda }_H}{\tau V^{\ast }+{\delta }_H},$and (10b(10b) $\begin{array}{rl}\tau H^{\ast }{V}^{\ast }-{\delta }_I{I}^{\ast }& =0,\end{array}$(10b) ) gives $I^{\ast }=\frac{\tau H^{\ast }{V}^{\ast }}{{\delta }_I}=\frac{\tau V^{\ast }{\lambda }_H}{{\delta }_I(\tau V^{\ast }+{\delta }_H)}.$Similarly, from (10c(10c) $\begin{array}{rl}n{\delta }_I{I}^{\ast }-c{V}^{\ast }-k_V{A}^{\ast }{V}^{\ast }& =0,\end{array}$(10c) ), (10h(10h) $\begin{array}{rl}(1-f)\mu B_A^{\ast }-s_M{V}^{\ast }{B}_M^{\ast }& =0,\end{array}$(10h) ), (10f(10f) $\begin{array}{rl}\mathrm{f\mu }{B}_A^{\ast }-{\delta }_P{P}^{\ast }& =0,\end{array}$(10f) ) and (10i(10i) $\begin{array}{rl}{\pi }_P{P}^{\ast }-{\delta }_A{A}^{\ast }& =0.\end{array}$(10i) ), we have (11) $\begin{array}{rl}{A}^{\ast }& =\frac{-\mathrm{c\tau }{V}^{\ast }+\mathrm{n\tau }{\lambda }_H-c{\delta }_H}{k_V(\tau V^{\ast }+{\delta }_H)},\end{array}$(11) (12) $\begin{array}{rl}{P}^{\ast }& =\frac{{\delta }_A(-\mathrm{c\tau }{V}^{\ast }+\mathrm{n\tau }{\lambda }_H-c{\delta }_H)}{{\pi }_P{k}_V(\tau V^{\ast }+{\delta }_H)},\end{array}$(12) (13) $\begin{array}{rl}{B}_A^{\ast }& =\frac{{\delta }_P{\delta }_A(-\mathrm{c\tau }{V}^{\ast }+\mathrm{n\tau }{\lambda }_H-c{\delta }_H)}{\mathrm{\mu f}{\pi }_P{k}_V(\tau V^{\ast }+{\delta }_H)},\end{array}$(13) (14) $\begin{array}{rl}{B_m}^{\ast }& =\frac{(1-f){\delta }_P{\delta }_A(-\mathrm{c\tau }{V}^{\ast }+\mathrm{n\tau }{\lambda }_H-c{\delta }_H)}{s_M{V}^{\ast }f{\pi }_P{k}_V(\tau V^{\ast }+{\delta }_H)},\end{array}$(14) (15) $\begin{array}{rl}{B}^{\ast }& =\frac{{\delta }_P{\delta }_A(-\mathrm{c\tau }{V}^{\ast }+\mathrm{n\tau }{\lambda }_H-c{\delta }_H)}{{\pi }_P{k}_V(\tau V^{\ast }+{\delta }_H)s{V}^{\ast }}.\end{array}$(15) Substituting these relationships into (10d(10d) $\begin{array}{rl}{\lambda }_B-s{V}^{\ast }{B}^{\ast }-{\delta }_B{B}^{\ast }& =0,\end{array}$(10d) ) gives (16) $A_2{V^{\ast }}^2+B_2{V}^{\ast }+C_2=0.$(16) Here $\begin{array}{rl}{A}_2& =\mathrm{s\tau }(c{\delta }_A{\delta }_P+{\pi }_P{k}_V{\lambda }_B)>0,\\ B_2& ={\delta }_P{\delta }_A\left[\mathrm{cs}{\delta }_H\right(1-{\mathcal{R}}_0)+\mathrm{c\tau }{\delta }_B]+s{k}_V{\pi }_P{\delta }_H{\lambda }_B,\\ C_2& =c{\delta }_H{\delta }_P{\delta }_A{\delta }_B(1-{\mathcal{R}}_0).\end{array}$Note that, for ${\mathcal{R}}_0>1$, $C_2<0$ and $A_2>0$. Thus, there exists a unique positive root for $V^{\ast }$ as long as ${\mathcal{R}}_0>1$. The simulation results closely resemble those of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). When compared to Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) with the same parameter values, the primary infection period increases, the temporary protection period decreases, and the positive equilibrium remains approximately the same. Consequently, it appears that memory B cells are incapable of achieving successful viral clearance.

However, there are other immune responses such as the CD8 T cell responses. In the next section, we consider the naïve CD8 T cells and effector CD8 T cells, and study whether they may lead to viral clearance without LLPC.

5. A within-host model with CD8 T cells

In this section, we incorporate naïve CD8 T cells and effector CD8 T cells into Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ), and show that the viral load may reach a positive equilibrium without LLPC.

5.1. Model formulation

Our new model contains two more classes: naïve CD8 T cells (T) and effector CD8 T cells ($T_c$). Naïve CD8 T cells are recruited at a constant rate ${\lambda }_T$. A naïve CD8 T cell is activated to become an effector CD8 T cell by antigen (here represented by V) at a rate σ. A naïve CD8 T cell dies at a rate ${\delta }_T$, and an effector CD8 T cell dies at a rate ${\delta }_{T_c}$. In addition, an infected target cell is killed by effector CD8 T cells at a rate m. A diagram of the extended model is shown in Figure 10. The extended model is given by the following system: (17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) (17b) $\begin{array}{rl}{I}^{\mathrm{\prime }}& =\mathrm{\tau HV}-\mathrm{mI}{T}_c-{\delta }_{I}I,\end{array}$(17b) (17c) $\begin{array}{rl}{V}^{\mathrm{\prime }}& =n{\delta }_{I}I-\mathrm{cV}-k_V\mathrm{AV},\end{array}$(17c) (17d) $\begin{array}{rl}{B}^{\mathrm{\prime }}& ={\lambda }_B-\mathrm{sVB}-{\delta }_{B}B,\end{array}$(17d) (17e) $\begin{array}{rl}{B}_A^{\mathrm{\prime }}& =\mathrm{sVB}-\mu B_A,\end{array}$(17e) (17f) $\begin{array}{rl}{P}^{\mathrm{\prime }}& =\mu B_A-{\delta }_{P}P,\end{array}$(17f) (17g) $\begin{array}{rl}{A}^{\mathrm{\prime }}& ={\pi }_{P}P-{\delta }_{A}A,\end{array}$(17g) (17h) $\begin{array}{rl}{T}^{\mathrm{\prime }}& ={\lambda }_T-\mathrm{\sigma VT}-{\delta }_{T}T,\end{array}$(17h) (17i) $\begin{array}{rl}{T}_c^{\mathrm{\prime }}& =\mathrm{\sigma VT}-{\delta }_{T_c}{T}_c.\end{array}$(17i) The meaning of parameters are listed in Table 1. All parameters are positive and all the initial values of Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) are nonnegative.

Figure 10. The flow chart for the within-host virus Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ).

5.2. Basic reproduction number

In this section, we compute the basic reproduction number ${\mathcal{R}}_0$ of the viral infection. First, note that the set $\begin{array}{rl}{\mathrm{\Omega }}_4& :=\left\{\right(H,I,V,B,B_A,P,A,T,T_c)\in {\mathrm{\Re }}_{+}^9:H+I\le \frac{{\lambda }_H}{{\delta }_H},B+B_A\le \frac{{\lambda }_B}{min\{{\delta }_B,\mu \}},\\ & T+T_c\le \frac{{\lambda }_T}{min\{{\delta }_T,{\delta }_{T_c}\}}\}\end{array}$is positively invariant for Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ). The DFE is $(\frac{{\lambda }_H}{{\delta }_H},0,0,\frac{{\lambda }_B}{{\delta }_B},0,0,0,\frac{{\lambda }_T}{{\delta }_T},0)$.

We linearize Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) about the DFE, and pick the infection-related terms to find the new-infection matrix $F=\left[\begin{array}{cccccc}0& \tau \frac{{\lambda }_H}{{\delta }_H}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$and the transition-related terms to find the transition matrix $V=\left[\begin{array}{cccccc}{\delta }_I& 0& 0& 0& 0& 0\\ -n{\delta }_I& c& 0& 0& 0& 0\\ 0& -\frac{s{\lambda }_B}{{\delta }_B}& \mu & 0& 0& 0\\ 0& 0& -\mu & {\delta }_P& 0& 0\\ 0& 0& 0& -{\pi }_P& {\delta }_A& 0\\ 0& -\frac{\sigma {\lambda }_T}{{\delta }_T}& 0& 0& 0& {\delta }_{T_c}\end{array}\right].$The next generation matrix is thus $F{V}^{-1}=\left[\begin{array}{cccccc}\frac{\mathrm{n\tau }{\lambda }_H}{c{\delta }_H}& \frac{\tau {\lambda }_H}{c{\delta }_H}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right].$The basic reproduction number ${\mathcal{R}}_0$ is the spectral radius of the next generation matrix, which is also given by (2(2) ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\mathrm{n\tau }{\lambda }_H}{c{\delta }_H}.$(2) ), i.e. it is the same as the basic reproduction number of Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ).

5.3. Positive equilibrium

Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) has an endemic steady state ($H^{\ast },I^{\ast },B^{\ast },B_A^{\ast },P^{\ast },A^{\ast },T^{\ast },T_c^{\ast }$), where $H^{\ast }>0$, $I^{\ast }>0$, $B^{\ast }>0$, $B_A^{\ast }>0$, $P^{\ast }>0$, $A^{\ast }>0$, $T^{\ast }>0$, and $T_c^{\ast }>0$. Let the right hand side of the Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) equal zero, we obtain $\begin{array}{rl}{H}^{\ast }& =\frac{{\lambda }_H}{\tau V^{\ast }+{\delta }_H},\\ I^{\ast }& =\frac{\tau {\lambda }_H{V}^{\ast }}{(\tau V^{\ast }+{\delta }_H)(\frac{\mathrm{m\sigma }{\lambda }_T{V}^{\ast }}{(\sigma V^{\ast }+{\delta }_T){\delta }_{T_c}}+{\delta }_I)},\\ B^{\ast }& =\frac{{\lambda }_B}{s{V}^{\ast }+{\delta }_B},\\ B_A^{\ast }& =\frac{s{V}^{\ast }{\lambda }_B}{\mu (s{V}^{\ast }+{\delta }_B)},\\ P^{\ast }& =\frac{s{V}^{\ast }{\lambda }_B}{{\delta }_P(s{V}^{\ast }+{\delta }_B)},\\ A^{\ast }& =\frac{{\pi }_{P}s{V}^{\ast }{\lambda }_B}{{\delta }_A{\delta }_P(s{V}^{\ast }+{\delta }_B)},\\ T^{\ast }& =\frac{{\lambda }_T}{\sigma V^{\ast }+{\delta }_T},\\ {T_c}^{\ast }& =\frac{\sigma V^{\ast }{\lambda }_T}{{{\delta }_T}_c(\sigma V^{\ast }+{\delta }_T)}.\end{array}$Then we have the following equation as to $V^{\ast }$: (18) $A_3{V^{\ast }}^3+B_3{V^{\ast }}^2+C_3{V}^{\ast }+D_3=0,$(18) where $\begin{array}{rl}{A}_3& =\mathrm{s\tau \sigma }(m{\lambda }_T+{\delta }_I{\delta }_{T_c})(c{\delta }_P{\delta }_A+{\pi }_P{k}_V{\lambda }_B),\\ B_3& =\left[{\delta }_P\right(\left(\right(1-{\mathcal{R}}_0)\mathrm{cs}{\delta }_H+\mathrm{c\tau }{\delta }_B)\sigma +\mathrm{cs\tau }{\delta }_T){\delta }_A+{\lambda }_{B}s{k}_V{\pi }_P(\sigma {\delta }_H+\tau {\delta }_T\left)\right]{\delta }_I{\delta }_{T_c}\\ & +\left(c{\delta }_P\right(s{\delta }_H+\tau {\delta }_B){\delta }_A+{\delta }_H{\lambda }_{B}s{k}_V{\pi }_P)\mathrm{m\sigma }{\lambda }_T,\\ C_3& =\left({\delta }_P\right(c{\delta }_H{\delta }_B(1-{\mathcal{R}}_0)\sigma +\left(\right(1-{\mathcal{R}}_0)\mathrm{cs}{\delta }_H+\mathrm{c\tau }{\delta }_B){\delta }_T){\delta }_A+{\delta }_H{\delta }_T{\lambda }_{B}s{k}_V{\pi }_P){\delta }_I{\delta }_{T_c}\\ & +\mathrm{delt}{a}_A{\delta }_B{\delta }_H{\delta }_P{\lambda }_T\mathrm{cm\sigma },\\ D_3& =c{\delta }_H{\delta }_A{\delta }_B{\delta }_P{\delta }_T{\delta }_I{\delta }_{T_c}(1-{\mathcal{R}}_0).\end{array}$Note that the coefficients $B_3$, $C_3$ and $D_3$ depend on ${\mathcal{R}}_0$. Since $D_3<0$ when ${\mathcal{R}}_0>1$, Equation (18(18) $A_3{V^{\ast }}^3+B_3{V^{\ast }}^2+C_3{V}^{\ast }+D_3=0,$(18) ) has one or three positive roots by Descartes' rule of sign. Thus, Model (17(17a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(17a) ) at least has a positive solution $V^{\ast }$ when ${\mathcal{R}}_0>1$. This means the virus cannot be cleared when CD8 T cells are considered. The simulation results closely resemble those obtained from Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ). When compared to Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) using identical parameter values, both the primary infection period and the positive equilibrium exhibit a decrease, while the temporary protection period shows an increase. These findings indicate that CD8 T cells play a crucial role in promoting viral clearance.

6. Concluding remarks

We develop a sequence of within-host mathematical models to study viral clearance via the adaptive immune system. We first created a simple Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) that considers the activation, proliferation and neutralization process, keeping track of the uninfected target cells, infected target cells, virus, B cells, PC and antiviral antibody. We find that the virus can be effectively controlled or cleared if the lifespans of PC or B cells are very long. To understand the effect of the long lifespan of B cells, we then considered memory B cells, and show that the positive viral equilibrium exists. On the other hand, after incorporating LLPC into Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ), we show that LLPC eliminate the positive equilibrium, and the viral load eventually approaches 0 without considering memory B cells. To understand if other factors, such as CD8 T cell response, may also contribute to viral clearance, we extend Model (1(1a) $\begin{array}{rl}{H}^{\mathrm{\prime }}& ={\lambda }_H-\mathrm{\tau HV}-{\delta }_{H}H,\end{array}$(1a) ) to include CD8 T cell and effector T cell dynamics without considering LLPC, and show that a positive viral equilibrium still exists. We conclude that LLPC are the key factor that has significant influence on the full clearance of the virus.

The virus may be cleared if the viral load is small enough, not only because the viral load reaching a fraction of a viral particle is meaningless, but also because stochasticity may cause viral clearance when only a small number of viral particles are left. We thus consider the virus being cleared if the viral load is below a threshold (specifically, less than $1{\mathrm{EID}}_{50}$). Using the realistic parameter values in [10], We show that the viral load in the absence of LLPC is always above the threshold, so viral clearance cannot be achieved. However, if the effectiveness of antibodies (or antibody amount) is increased, the virus may indeed be cleared without LLPC. Yet this clearance is only temporary. As the healthy target cell count quickly recovers, re-exposure to the viral causes reinfection only days or a few weeks after the infection. This may account for the observed re-infections during the same wave of the COVID-19 pandemic [19,20]. Even if the equilibrium viral load may be lower than the threshold, it may take a long time for the viral load to drop below the threshold. This may account for the existence of long-COVID cases [21,22].

Our models do not consider the full spectrum of adaptive immune response, such as regulatory T cells, and dendritic cells. But we numerically demonstrated that the behaviour of our model is similar to that of a more realistic model in [10]. Our models also suggest that there may exist temporary protection due to variations in parameter values without LLPC. The virus load may have a positive equilibrium or oscillate behaviour with different parameters. Temporary protection may delay secondary viral infections, but the virus exists eventually. Thus, our models are sufficient to demonstrate the key role of LLPC in viral clearance.

Acknowledgments

We thank the anonymous reviewers for their constructive comments and suggestions.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

We thank the anonymous reviewers for their constructive comments. This work was supported by the National Natural Science Foundation of China (No. 11771075, 12271088) (ML), Natural Science Foundation of Shanghai (No. 21ZR1401000, No. 23ZR1401700) (ML), and two discovery grants of Natural Sciences and Engineering Research Council Canada (RE, JM), two NSERC EIDM grants (OMNI and MfPH) (JM), and the China Scholarship Council (202106630037) (MZ).