Mathematical model on HIV and nutrition

ABSTRACT

HIV continues to be a major global health issue, having claimed millions of lives in the last few decades. While several empirical studies support the fact that proper nutrition is useful in the fight against HIV, very few studies have focused on developing and using mathematical modelling approaches to assess the association between HIV, human immune response to the disease, and nutrition. We develop a within-host model for HIV that captures the dynamic interactions between HIV, the immune system and nutrition. We find that increased viral activity leads to increased serum protein levels. We also show that the viral production rate is positively correlated with HIV viral loads, as is the enhancement rate of protein by virus. Although our numerical simulations indicate a direct correlation between dietary protein intake and serum protein levels in HIV-infected individuals, further modelling and clinical studies are necessary to gain comprehensive understanding of the relationship.

KEYWORDS:

• HIV
• nutrition
• reproduction number
• within-host mathematical model
• backward bifurcation
• protein

MATHEMATICS SUBJECT CLASSIFICATION:

• 92D25

1. Introduction

Human Immunodeficiency Virus (HIV) is a virus that attacks cells of the human immune system (i.e. cells that help the human body to fight against various infections), making a person more vulnerable to other infections. HIV, which causes acquired immune deficiency syndrome (AIDS), has been a major public health challenge since its first reported case in 1981 [1]. In 2021, approximately 36.7 (1.7) million adults (children) across the world were living with HIV and about 680,000 HIV/AIDS-related deaths were reported [2,3]. Overall, HIV/AIDS has killed approximately 40.1 million humans globally from the time it was first identified (i.e. from 1981) to 2021 [2].

Dynamic interactions between HIV and the human immune system are complex. These interactions have been studied extensively through rigorous mathematical modelling approaches. Most of the research papers on the subject have focused on data- or biologically-driven mathematical models [4,5], while some have been statistical [6,7] or focused on the epidemiology of HIV [8]. A review paper by Jessica M. Conway and Ruy M. Ribeiro summarizes most of the work on the immunology of HIV and offers a comprehensive basis for current work on the topic [9]. Additional sources of information on viral dynamics and immune response to HIV can be found in [10,11].

The fundamental principle behind most of these models is simple: when viral particles meet susceptible cells, they infect the susceptible cells. Some of these infected cells are killed by immune cells, while the other infected cells produce more viruses, leading to a self-sustaining cycle that generates additional infected cells. Although replication of the virus can be modelled with one density-dependent logistic growth equation [11], more elaborate, yet simple models, that are able to capture more biological aspects of the virus require more variables (see, for example, [10,12,13]). These simple mathematical frameworks account for the dynamics of susceptible and infected cells, as well as the free virus, and have been useful in understanding certain aspects of the disease even though they do not capture full biological details of the dynamic interplay between HIV and the immune system. For example, some models in [10,12,13] account for preferential attacks against activated T cells. The simplest models describing interactions between HIV and the immune system, particularly CD8 T cells incorporate an additional compartment for CD8 T cells [12,14]. In addition to these simple models, more complex models have been developed and used to answer more subtle and biologically significant questions on HIV and the human immune system. These models have broad theoretical and practical applications, ranging from CD8 T cell induced pathology [15,16] through modelling of CD8 T cell escape (i.e. the process by which HIV escape recognition by immune effector cells) [17–19] and HIV latency [20–22]. Unfortunately, none of these mathematical models focus on interactions between nutrition, HIV and the human immune system.

Dynamic interplay between nutrition, HIV and the human immune system have been reported. Specifically, HIV weakens the human immune system thereby reducing the ability of the human immune system to fight against diseases significantly [23,24]. On the other hand, the human immune system requires a good nutritional balance to fight such diseases. The formation of immune cells requires various micro and macro nutrients [25–28]. In the absence of proper nutrition, this immune cell formation process can be weakened, rendering a person who is infected with HIV more vulnerable to other opportunistic infections. Although, several studies have focused on the relationship between malnutrition and HIV [29–33], there are still unanswered questions on the impact of dynamic interactions between HIV, immune response to HIV, and nutrition.

Unlike most infectious diseases in which carbohydrates and fats are broken down to supplement the increased nutritional needs of the immune system, HIV induces a special metabolic effect that initiates a preferential loss of proteins over fats and other macro nutrients [34–39]. In particular, HIV is well known for loss of muscle mass and the ‘wasting syndrome’ (i.e. progressive involuntary weight loss of $10\mathrm{\%}$ of baseline body weight in the setting of a chronic infection and/or chronic diarrhea [40]) [41,42]. In addition, studies have found some abnormalities in serum protein levels of HIV infected individuals. Serum protein consists of approximately $60\mathrm{\%}$ albumin, which is considered to be a major indicator of malnutrition [43], and approximately $40\mathrm{\%}$ globulin. Studies have found higher levels of serum protein and globulin, and lower levels of albumin in HIV infected individuals compared to their uninfected counterparts [44–47]. A summary of some of the results from various studies that examined the total protein, albumin, and globulin levels in individuals with and without HIV infection is presented in Table 1. Despite these interesting studies, little to no mathematical modelling effort has been invested in accounting for the relationship between nutrition and human immune response to HIV and the implications of this relationship to HIV dynamics.

Table 1. Total protein, albumin, and globulin levels in individuals with and without HIV infection.

	HIV Positive		HIV Negative
	Men	Women	Men	Women	Reference
Total protein g/L	$76.0\pm 1.3$	$79.3\pm 1.5$	$72.3\pm 0.1$	$71.2\pm 0.1$
Albumin g/L	$41.7\pm 0.6$	$38.6\pm 0.8$	$44.6\pm 0.1$	$42.0\pm 0.1$	[44]
Total Protein g/dL	$7.72(7.45-8)$		$7.8(7.57-8.02)$
Albumin g/dL	$4.4(4.27-4.54)$		$4.59(4.25-4.93)$		[46]
IgG g/dL	N/A		N/A
Total protein g/L, mean (SD)	$85.4\left(7.9\right)$		$75.5\left(6.1\right)$
Albumin g/L, mean (SD)	$40\left(3.7\right)$		$43\left(4.3\right)$		[47]
IgG g/L, mean (SD)	$27\left(8.7\right)$		$16\left(3.9\right)$
IgG mg/mL	1998 (655–4411)		1458 (710–3122)
IgA mg/dL	294 (47–1200)		232 (110–593)		[48]
IgM mg/dL	206 (62–661)		284 (47–709)

Here, we develop a within-host model for the dynamics of HIV that accounts for interactions between HIV, human immune response to HIV, and nutrition. While this coupled HIV-immune response-nutrition model is new, it is based primarily on well-known within-host HIV models that account for human immune response to HIV [12,49,50]. We incorporate nutrition into the model via protein, which is arguably the most important nutrition factor with regards to HIV. For mathematical tractability, we do not account for separate components of serum protein (such as albumin and globulin). The rest of the paper is organized as follows: The coupled HIV-immune response-nutrition model is developed in Section 2. The immunological reproduction number of the model is computed in Section 3 and used to establish the existence and stability of equilibria to the model system, as well as the existence of a backward bifurcation in Section 3. Numerical simulations to assess the impact of important nutritional parameters on the total protein and HIV viral load are presented in Section 4, while a discussion and concluding remarks are presented in Section 5.

2. The within-host model

2.1. Formation of the within-host model

In this section, we develop a within-host mathematical model for the dynamics of HIV, which accounts for the interactions between HIV, the human immune system, and nutrition. This entails introducing nutrition (specifically, protein) to the target-cell limited model developed in [12] and immune control as in [50]. Since, the primary target of HIV is the CD4 T cells [12,51], we consider CD4 T cells as the target cells cells in our model. These target cells are denoted by T. Infected cells are denoted by $T_i$, while the virus is denoted by V. Since CD8 T cells constitute the dominant defense mechanism against HIV [12], the cellular immune response is assumed to be dominated by CD8 T cells, which are denoted by Z in our model. Protein is denoted by η. Schematics of the model are depicted in Figure 1.

Figure 1. Schematic diagram of the within-host model. CD4 T cells are denoted by T, infected cells are denoted by $T_i$, the virus is denoted by V , CD8 T cells are denoted by Z and nutrition (protein) is denoted by η.

The effect of malnutrition on the immune system is well-documented [29,30]. In particular, studies have shown that malnutrition leads to a decrease in the total number of lymphocytes and CD4 T lymphocytes, while CD8 T lymphocytes are relatively maintained [29,30,52]. Taking this into account, we model the production of target (i.e. CD4 T) cells and the production of immune (i.e. CD8 T) cells as functions of protein, η. Specifically, target (immune) cells are produced at rate r (${\lambda }_z$). Following the target cell limited model in [12], we assume that infected cells ($T_i$) are produced at rate $\rho \left(\eta \right)V$, where V is the population density of the free virus and $\rho \left(\eta \right)$–a function of protein is given by $\rho \left(\eta \right)=\frac{{\rho }_0}{1+A\eta }$. Here, ${\rho }_0$ is the rate at which target cells are infected per free virus per day and A is a non-negative constant. Observe that the ρ is a decreasing function of η. Target cells die at per capita rate d per day, and the loss of infected cells due to viral cytopathicity is at per capita rate, δ per day. The immune-mediated cytotoxic effect on infected cells by CD8 T cells is modelled with the term $\psi \left(\eta \right)T_{i}Z$, where $\psi \left(\eta \right)=\frac{{\psi }_0\eta }{\mathrm{\Psi }+\eta }$ is the rate at which infected cells are killed by activated CD8 T cells, ${\psi }_0$ is the killing rate of infected cells per CD8 T cell per day and Ψ is the half-saturation constant. It should be noted that Ψ allows the killing rate of infected cells to be independent of nutrition, while A in $\rho \left(\eta \right)$ allows us to remove the dependence of the infectivity rate of target cells on nutrition. Virions (V) are produced from infected cells ($T_i$) at rate of p and cleared naturally at per capita rate c. The clearance of virus by immunoglobulins is modelled by the term ${\mu }_V\eta V$, where ${\mu }_V$ is the virus clearance rate by immunoglobulins per gram of protein per day. In the presence of higher protein levels, the immunoglobulin levels are also higher and in this case, we assume that the virus is cleared at a faster rate. Following [12], the activation of CD8 T cells is modelled by the term $b{T}_{i}Z$ that is dependent on the density of infected cells, where b is the antigen activation rate per infected cell per day. The per capita death rate of CD8 T cells ${\mu }_Z$ per day. The generation rate of protein is ${\lambda }_{\eta }$ grams per day, the per capita clearance rate of protein is ${\mu }_{\eta }$ per day, while virus-driven increase in total protein levels is modelled by the term ${\gamma }_{\eta }\eta V$, where ${\gamma }_{\eta }$ is the enhancement rate of protein per virus per day. It is worth noting that the positive term ${\gamma }_{\eta }\eta V$ coincides with the observed high levels of serum protein in HIV infected individuals. While this term may seem somewhat counter-intuitive, it aligns with the biological data at hand [44,46–48], previous studies on the subject[52] and our numerical simulations. We hypothesize that an increased presence of viral activity in body due to HIV infection leads to an increased antibodies (immunoglobulins) and this in turn leads to higher levels of globulin and serum protein. The cause and effect of this significant phenomena is further discussed in Section 5 where all our results are compared with already existing biological data. Brief descriptions of the model parameters and their units, as well as numerical values are presented in Table 2.

Using the schematics in Figure 1 together with the variable and parameter descriptions above, we obtain the following coupled HIV-Immune response-nutrition model: (1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1)

Table 2. Brief descriptions of the parameters of the within-host model together with the numerical values and units of the parameters and initial conditions.

Param.	Description	Value	Units
$T\left(0\right)$	Initial value for target cells	10	cells/μL
$T_i\left(0\right)$	Initial value for infected cells	0	cells/μL
$V\left(0\right)$	Initial value for the virus	${10}^{-9}$	virions/μL
$Z\left(0\right)$	Initial value for immune cells	6	cells/μL
$\eta \left(0\right)$	Initial value for total protein	5	gr
$r_0$	Target cell production	10	1/μL.day
r	Target cell production per gram total protein	10/350	1/μL.day.gr
β	Infection rate of target cells	$6.5\times {10}^{-4}$	μL/virion.day
${\rho }_0$	Infection rate of target cells times gram	$6.5\times {10}^{-4}\times 350$	μL.gr/virion.day
Ψ	Half saturation constant for nutrition	175	gr
A	Proportionality constant for nutrition	0.1	1/gr
d	Death rate of uninfected cells	0.01	1/ day
δ	Death rate of infected cells	1	1/day
p	Virus production rate	2000	virions/cell.day
c	Clearance rate of virus	23	1/day
${\lambda }_0$	CD8 T cell recruitment	1	cell/μL.day
${\lambda }_z$	CD8 T cell recruitment per gram total protein	${10}^{-3}/350$	cells/μL.day.gr
${\psi }_0$	Strength of CD8 T response	0.42	$\text{μ}\mathrm{L}$/cell.day
${\mu }_z$	Death rate of CD8 T cells	2	1/day
${\lambda }_{\eta }$	Rate of ingestion of total protein	65	gr/day
${\text{μ}}_{\eta }$	Clearance of total protein in the body	1.1857	1/day
${\gamma }_{\eta }$	Enhancement rate of protein by virus	0.0001	1/day
b	Antigen driven activation rate of CD8 T	1	1/day
${\mu }_V$	Viral clearance by immunoglobulins	0.1057	1/gr.day

2.2. Existence, uniqueness, positiveness and boundedness of solutions

The existence of a unique, positive solution locally in time can easily be verified for the system given by (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ). Let $F:{\mathbb{R}}_{+}^6⟶{\mathbb{R}}^5$ be defined by $F(\tau ,x)=\left(F_1\right(\tau ,x),F_2(\tau ,x),F_3(\tau ,x),F_4(\tau ,x),F_5(\tau ,x\left)\right)$, where $x=(T,T_i,V,Z,\eta )$, and $F_i$ denote the five differential equations in the system. Since F and $\frac{\mathrm{\partial }F}{\mathrm{\partial }x}$ are continuous on ${\mathbb{R}}_{+}^6$, F is locally Lipschitz. Now, the existence of a unique, positive solution defined on some interval $[0,{\tau }_{max})$ with ${\tau }_{max}>0$ follows by Theorem A.4 in [53].

The boundedness of solutions of the system given by (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) cannot be established without further imposing conditions on model parameters. In fact, in Appendix 1, we establish the boundedness of solutions under certain assumptions on parameter values.

2.3. Limitations of the model

Model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) is based on some simplifying assumptions, which might limit its applicability. However, relaxing some of these assumptions will render the model more complex and mathematically intractable. For example, a logistic-like growth function can be used for the growth of CD4 T cells to account for the data measuring the number of dividing helper cells in HIV infected individuals with different CD4 T cell counts [12]. But this will make the within-host model more complicated since we must also take the nutritional aspects into consideration. Also, the immune control component of the model, given by the fourth equation of (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ), can be made more realistic by using an activation function that saturates [12,50], e.g. (2) $\begin{array}{r}\frac{\mathrm{d}Z}{\mathrm{d}\tau }={\lambda }_z\eta +\frac{b_Z{T}_{i}Z}{K_B+T_i}-{\mu }_{z}Z.\end{array}$(2) Here, $b_Z$ is the maximal rate of CD8 T activation. Furthermore, we can improve on the way cellular immunity is modelled. For example, immune impairment due to high viral loads, and hence high infected cell densities, can be taken can be modelled through another saturation function [50], which will lead to the following modified version of the fourth equation of (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ): (3) $\begin{array}{r}\frac{\mathrm{d}Z}{\mathrm{d}\tau }={\lambda }_z\eta +\frac{b_Z{T}_{i}Z}{K_B+T_i}-\frac{{\mathrm{d}}_Z{T}_{i}Z}{K_D+T_i}-{\mu }_{z}Z.\end{array}$(3) Additionally, considering abnormalities related to albumin and globulin concentrations in HIV-infected individuals, albumin and globulin can be modelled separately, instead of the total protein (see, for example, [52]).

The four linear functions of η, namely, $r\eta ,{\mu }_v\eta ,{\lambda }_z\eta$ and ${\gamma }_{\eta }\eta$ were assumed to be linear for the sake of mathematical and analytical simplicity. These four functions, however, are unbounded and increasing to ∞ as $\eta ⟶\mathrm{\infty }$. Due to nutrition's limited role in fighting the virus, it would be more realistic to use saturated functions in place of these linear functions.

3. Analytical results

3.1. The infection-free equilibrium $\left(E_0\right)$

Equilibria are important in determining the long-term behaviour of solutions to autonomous systems of differential equations that cannot be solved analytically [54]. To identify the equilibria of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ), we set the left-hand sides of the system to zero and solve the resulting system of algebraic equations for five state variables $T^{\ast },{T_i}^{\ast },V^{\ast },Z^{\ast }$, and ${\eta }^{\ast }$ (where the ‘*’ denotes ‘the equilibrium value of ’) simultaneously. Specifically, any equilibrium point $(T^{\ast },{T_i}^{\ast },V^{\ast },Z^{\ast },{\eta }^{\ast })$, of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) must satisfy the following system of equations: (4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) Solving for $T^{\ast },{T_i}^{\ast },Z^{\ast }$, and ${\eta }^{\ast }$ in terms of $V^{\ast }$ from the first, third, fourth, and fifth equations of system (4(4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) ) leads to (5) $\begin{array}{rl}{T}^{\ast }& =\frac{r{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(d+\rho ({\eta }^{\ast }\left)V^{\ast }\right)},\end{array}$(5) (6) $\begin{array}{rl}{T_i}^{\ast }& =\frac{(c+\frac{{\lambda }_{\eta }{\mu }_V}{{\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }})V^{\ast }}{p},\end{array}$(6) (7) $\begin{array}{rl}{Z}^{\ast }& =\frac{p{\lambda }_Z{\lambda }_{\eta }}{p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })-b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]},\end{array}$(7) (8) $\begin{array}{rl}{\eta }^{\ast }& =\frac{{\lambda }_{\eta }}{{\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }}.\end{array}$(8) It should be noted that for defined and non-negative equilibria, $0\le V^{\ast }<\frac{{\mu }_{\eta }}{{\gamma }_{\eta }}$. From the second equation of the system (4(4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) ), i.e. (9) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-{T_i}^{\ast }[\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right]=0,\end{array}$(9) and (5(5) $\begin{array}{rl}{T}^{\ast }& =\frac{r{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(d+\rho ({\eta }^{\ast }\left)V^{\ast }\right)},\end{array}$(5) )–(9(9) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-{T_i}^{\ast }[\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right]=0,\end{array}$(9) ), it is clear that $V^{\ast }=0$ is a solution of (4(4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) ). This is the infection-free equilibrium value of $V^{\ast }$, which we will denote by $V_0$. To obtain the infection-free equilibrium values of the other state variables (i.e. $T^{\ast },{T_i}^{\ast },Z^{\ast }$, and ${\eta }^{\ast }$) denoted by $T_0,{T_i}_0,Z_0$, and ${\eta }_0$, respectively, we substitute $V^{\ast }=0$ in Equations (5(5) $\begin{array}{rl}{T}^{\ast }& =\frac{r{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(d+\rho ({\eta }^{\ast }\left)V^{\ast }\right)},\end{array}$(5) )–(9(9) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-{T_i}^{\ast }[\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right]=0,\end{array}$(9) ). This gives (10) $\begin{array}{r}{E}_0:=(T_0,{T_i}_0,V_0,Z_0,{\eta }_0)=(\frac{r{\lambda }_{\eta }}{\mathrm{d}{\mu }_{\eta }},0,0,\frac{{\lambda }_Z{\lambda }_{\eta }}{{\mu }_Z{\mu }_{\eta }},\frac{{\lambda }_{\eta }}{{\mu }_{\eta }}).\end{array}$(10)

3.2. Immunological reproduction number (${\mathfrak{R}}_0$)

The immunological reproduction number, ${\mathfrak{R}}_0$, represents the number of secondary virus particles that one viral particle will produce in an entirely susceptible target cell population through out the period within which the viral particle is capable of producing secondary virus particles. In this Section, we compute the reproduction number of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) using two approaches–the next generation operator approach proposed by van den Driessche and Watmough [55] and an intuitive approach with a direct biological interpretation.

3.2.1. Reproduction number through the next generation operator approach

In this subsection, we use the next generation operator technique proposed by van den Driessche and Watmough [55] to calculate the immunological reproduction number of the within-host model given by system (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ). Following this method and using the same notation as in [55], we define the matrices of new infections (F) and transitions V as follows: $\begin{array}{rl}F:& =\left[\frac{\mathrm{\partial }{\mathcal{F}}_m(0,y_0)}{\mathrm{\partial }{x}_m}\right]=\left[\begin{array}{cc}0& \rho \left({\eta }_0\right)T_0\\ p& 0\end{array}\right]\mathrm{a}\mathrm{n}\mathrm{d}\\ V:& =\left[\frac{\mathrm{\partial }{\mathcal{V}}_m(0,y_0)}{\mathrm{\partial }{x}_m}\right]=\left[\begin{array}{cc}\delta +\psi \left({\eta }_0\right)Z_0& 0\\ 0& c+{\mu }_V{\eta }_0\end{array}\right],\end{array}$where ${\mathcal{F}}_m,m=1,2$ is the vector of new infections, ${\mathcal{V}}_m,m=1,2$ is the vector of new transition terms, $x_1=T_i$, $x_2=V$, and ${ϵ}_0=(0,y_0)$, with $y_0=(T_0,Z_0,{\eta }_0)=(\frac{r{\lambda }_{\eta }}{\mathrm{d}{\mu }_{\eta }},\frac{{\lambda }_Z{\lambda }_{\eta }}{{\mu }_Z{\mu }_{\eta }},\frac{{\lambda }_{\eta }}{{\mu }_{\eta }})$ is the unique unique equilibrium of the infection-free system $y^{\prime }=g(0,y)$. The next-generation matrix is then defined by $F{V}^{-1}=\left[\begin{array}{cc}0& {\displaystyle \frac{\rho \left({\eta }_0\right)T_0}{c+{\mu }_V{\eta }_0}}\\ {\displaystyle \frac{p}{\delta +\psi \left({\eta }_0\right)Z_0}}& 0\end{array}\right]=\left[\begin{array}{cc}0& {\mathfrak{R}}_1\\ {\mathfrak{R}}_2& 0\end{array}\right].$The reproduction number of the model system (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) is then defined to be the spectral radius of the next-generation matrix, $F{V}^{-1}$. The eigenvalues of the $2\times 2$ matrix ($F{V}^{-1}$) are given by the characteristic equation: $det(F{V}^{-1}-\lambda I)=0$, which simplifies to ${\lambda }^2-{\mathfrak{R}}_1{\mathfrak{R}}_2=0.$Thus, the reproduction number from the next-generation operator approach is (11) $\begin{array}{r}{{\mathfrak{R}}_0}^{NG}=\sqrt{{\mathfrak{R}}_1{\mathfrak{R}}_2}=\sqrt{\frac{\rho \left({\eta }_0\right)T_0}{c+{\mu }_V{\eta }_0}\times \frac{p}{\delta +\psi \left({\eta }_0\right)Z_0}}=\sqrt{\frac{p\rho \left({\eta }_0\right)T_0}{(c+{\mu }_V{\eta }_0)(\delta +\psi ({\eta }_0\left)Z_0\right)}}.\end{array}$(11) Since $\rho \left({\eta }_0\right)=\frac{{\rho }_0}{1+A{\eta }_0}$ and $\psi \left({\eta }_0\right)=\frac{{\psi }_0{\eta }_0}{\mathrm{\Psi }+{\eta }_0}$, where ${\eta }_0=\frac{{\lambda }_{\eta }}{{\mu }_{\eta }}$, and since $Z_0=\frac{{\lambda }_Z{\lambda }_{\eta }}{{\mu }_Z{\mu }_{\eta }}$, the expression for the reproduction number in (11(11) $\begin{array}{r}{{\mathfrak{R}}_0}^{NG}=\sqrt{{\mathfrak{R}}_1{\mathfrak{R}}_2}=\sqrt{\frac{\rho \left({\eta }_0\right)T_0}{c+{\mu }_V{\eta }_0}\times \frac{p}{\delta +\psi \left({\eta }_0\right)Z_0}}=\sqrt{\frac{p\rho \left({\eta }_0\right)T_0}{(c+{\mu }_V{\eta }_0)(\delta +\psi ({\eta }_0\left)Z_0\right)}}.\end{array}$(11) )) can be simplified to (12) $\begin{array}{r}{{\mathfrak{R}}_0}^{NG}=\sqrt{\frac{pr{\rho }_0{\lambda }_{\eta }{\mu }_Z{{\mu }_{\eta }}^2(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })}{\mathrm{d}({\mu }_{\eta }+A{\lambda }_{\eta })(c{\mu }_{\eta }+{\mu }_V{\lambda }_{\eta })\left[\delta {\mu }_Z{\mu }_{\eta }\right(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })+{\psi }_0{\lambda }_Z{{\lambda }_{\eta }}^2]}}.\end{array}$(12)

3.2.2. Reproduction number through the intuitive approach

While the immunological reproduction number can be calculated using the next generation operator method, the biological interpretation of this important threshold might not be obvious from the final expression (especially due to the square-root in this expression). In this section, we derive the reproduction number through an intuitive approach in which terms and expressions and associated with meaningful biological interpretations.

Each infected cell produces virus at rate p. Also, infected cells die at rate δ due to viral cytopathogenic effects, or at rate $\psi \left({\eta }_0\right)Z_0$, where $Z_0$ and ${\eta }_0$ (given in (10(10) $\begin{array}{r}{E}_0:=(T_0,{T_i}_0,V_0,Z_0,{\eta }_0)=(\frac{r{\lambda }_{\eta }}{\mathrm{d}{\mu }_{\eta }},0,0,\frac{{\lambda }_Z{\lambda }_{\eta }}{{\mu }_Z{\mu }_{\eta }},\frac{{\lambda }_{\eta }}{{\mu }_{\eta }}).\end{array}$(10) )) are the infection-free equilibrium values of CD8 T cells and protein, respectively, due to immune-mediated cytotoxic effects. Because infected cells are cleared at a total rate of $\delta +\psi \left({\eta }_0\right)Z_0$, they live, on average, for $\frac{1}{\delta +\psi \left({\eta }_0\right)Z_0}$ days. Therefore, on average, each infected cell produces ${\mathfrak{R}}_2=\frac{p}{\delta +\psi \left({\eta }_0\right)Z_0}$ virions during its lifetime. On the other hand, Now, virions are cleared naturally at a per capita rate c, or by immunoglobulins at rate of ${\mu }_V{\eta }_0$. Thus, virions are cleared at a total rate of $c+{\mu }_V{\eta }_0$. Hence, on average, each virion lives for $\frac{1}{c+{\mu }_V{\eta }_0}$ days. During this time, each virion infects ${\mathfrak{R}}_1=\frac{\rho \left({\eta }_0\right)T_0}{c+{\mu }_V{\eta }_0}$ susceptible cells on average. Thus, the total number of cells infected by $\frac{p}{\delta +\psi \left({\eta }_0\right)Z_0}$ virions released from one infected cell is given by ${\mathfrak{R}}_1{\mathfrak{R}}_2=\frac{p\rho \left({\eta }_0\right)T_0}{(c+{\mu }_V{\eta }_0)(\delta +\psi ({\eta }_0\left)Z_0\right)}$. Thus, ${\mathfrak{R}}_0={\mathfrak{R}}_1{\mathfrak{R}}_2=\frac{p\rho \left({\eta }_0\right)T_0}{(c+{\mu }_V{\eta }_0)(\delta +\psi ({\eta }_0\left)Z_0\right)}.$

3.2.3. Comparing the reproduction numbers from the next generation operator and intuitive approaches

When the reproduction number is greater than 1, the reproduction number obtained using the intuitive approach $\left({\mathfrak{R}}_0\right)$ produces a higher value compared to the reproduction number obtained using next generation operator approach $\left({{\mathfrak{R}}_0}^{NG}\right)$. On the other hand, if the reproduction number is less than 1, the relationship between the two reproduction numbers is reversed. Specifically, the reproduction number from the intuitive approach is the square of that from the next generation operator approach, i.e. ${\mathfrak{R}}_0={\left({{\mathfrak{R}}_0}^{NG}\right)}^2.$However, the intuitive derivation process is more informative as it provides meaningful biological interpretations for the two non-zero entries (${\mathfrak{R}}_1$ and ${\mathfrak{R}}_2$) in the next generation operator approach. Irrespective of the difference, the epidemiological/immunological implications of the reproduction number that the disease can be cleared when ${\mathfrak{R}}_0<1$, or can persist when ${\mathfrak{R}}_0>1$, remains the same since ${\mathfrak{R}}_0<1\iff {\mathfrak{R}}_0^2<1$ and ${\mathfrak{R}}_0>1\iff {\mathfrak{R}}_0^2>1$.

3.3. Stability of the infection-free equilibrium

Local stability of the infection-free equilibrium ($E_0=(T_0,{T_i}_0,V_0,Z_0,{\eta }_0)$) of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) will be establish by considering the signs of the eigenvalues of the Jacobian of the model system computed at $E_0$. This Jacobian is given by $J_0:=J\left(E_0\right)=\left[\begin{array}{ccccc}-d& 0& -\rho \left({\eta }_0\right)T_0& 0& r\\ 0& -\delta -\psi \left({\eta }_0\right)Z_0& \rho \left({\eta }_0\right)T_0& 0& 0\\ 0& p& -c-{\mu }_V{\eta }_0& 0& 0\\ 0& b{Z}_0& 0& -{\mu }_Z& {\lambda }_Z\\ 0& 0& {\gamma }_{\eta }{\eta }_0& 0& -{\mu }_{\eta }\end{array}\right].$From the Jacobian ($J_0$), it is clear that three of the eigenvalues are ${\lambda }_1=-d<0,{\lambda }_2=-{\mu }_Z<0$, and ${\lambda }_3=-{\mu }_{\eta }<0$ and that the remaining two eigenvalues are eigenvalues can be determined from the $2\times 2$ sub-matrix $J_2=\left[\begin{array}{cc}-\delta -\psi \left({\eta }_0\right)Z_0& \rho \left({\eta }_0\right)T_0\\ p& -c-{\mu }_V{\eta }_0\end{array}\right].$Since the first three eigenvalues are negative, all we need to do to establish the stability of the infection-free equilibrium is to determine the signs of the two remaining eigenvalues. To this effect, we use Theorem 3.2 from [54], which states that a $2\times 2$ matrix is stable if the trace of the matrix is negative and the determinant is positive. Hence, we need to determine whether $Tr\left(J\right)<0$ and $Det\left(J\right)>0$. The trace of $J_2$ is given by $Tr\left(J_2\right)=-\delta -\psi \left({\eta }_0\right)Z_0-c-{\mu }_V{\eta }_0$, which is negative, while the determinant of $J_2$ is given by $Det\left(J_2\right)=(\delta +\psi ({\eta }_0\left)Z_0\right)(c+{\mu }_V{\eta }_0)-p\rho \left({\eta }_0\right)T_0$ $=(\delta +\psi ({\eta }_0\left)Z_0\right)(c+{\mu }_V{\eta }_0)(1-\frac{p\rho \left({\eta }_0\right)T_0}{(\delta +\psi ({\eta }_0\left)Z_0\right)(c+{\mu }_V{\eta }_0)})$ $=(\delta +\psi ({\eta }_0\left)Z_0\right)(c+{\mu }_V{\eta }_0)(1-{\mathfrak{R}}_0)$ $>0$, if ${\mathfrak{R}}_0<1$. Hence, we have the following standard result:

Theorem 3.1

If ${\mathfrak{R}}_0<1$, then the infection-free equilibrium ($E_0$) is locally asymptotically stable. If ${\mathfrak{R}}_0>1$, $E_0$ is unstable.

3.4. Existence of infection equilibria and backward bifurcation

In this Section, we show that the model system (4(4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) ) has a unique infection equilibrium when the reproduction number is greater than unity and that there is a parameter regime within which the model exhibits a backward bifurcation when the reproduction number is less than unity. The existence of a backward bifurcation implies that measures aimed at clearing the virus must be sustained even when the reproduction number is slightly less than one until the reproduction number falls below a smaller threshold value.

3.4.1. Existence of unique infection equilibrium

All equilibria of the model system (4(4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) ) must be non-negative. Hence, as indicated in Section 3.1, we must have $0<V^{\ast }<\frac{{\mu }_{\eta }}{{\gamma }_{\eta }}$. From (8(8) $\begin{array}{rl}{\eta }^{\ast }& =\frac{{\lambda }_{\eta }}{{\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }}.\end{array}$(8) ), (13) $\begin{array}{r}{Z}^{\ast }=\frac{p{\lambda }_Z{\lambda }_{\eta }}{p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })-b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]},\end{array}$(13) where $0<V^{\ast }<\frac{{\mu }_{\eta }}{{\gamma }_{\eta }}$ and $V^{\ast }$ can be determined from (9(9) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-{T_i}^{\ast }[\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right]=0,\end{array}$(9) ). For $Z^{\ast }$ to be positive, the denominator of the right-hand side must be positive. That is, we must have (14) $\begin{array}{rl}& p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })>b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]\\ & \iff b{\lambda }_{\eta }{\mu }_V{V}^{\ast }<({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(p{\mu }_Z-bc{V}^{\ast }).\end{array}$(14)

The left-hand side ($b{\lambda }_{\eta }{\mu }_V{V}^{\ast }$) of the inequality (14(14) $\begin{array}{rl}& p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })>b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]\\ & \iff b{\lambda }_{\eta }{\mu }_V{V}^{\ast }<({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(p{\mu }_Z-bc{V}^{\ast }).\end{array}$(14) ) is a straight line that goes through the origin, while the right-hand side ($({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(p{\mu }_Z-bc{V}^{\ast })$) is a parabola with roots $\widehat{V_1}=\frac{{\mu }_{\eta }}{{\gamma }_{\eta }}$ and $\widehat{V_2}=\frac{p{\mu }_Z}{bc}$ Figure 2. The parabola opens up (since the coefficient of $V^{\ast }$ is positive) and meets the straight line $y=b{\lambda }_G{\mu }_V{V}^{\ast }$ at two distinct points ($V_1$ and $V_2$), with $V_1<V_2$. Then $V_1<\mathrm{m}\mathrm{i}\mathrm{n}\{\widehat{V_1},\widehat{V_2}\}$ and hence $V_1<\frac{{\mu }_{\eta }}{{\gamma }_{\eta }}$. Since this satisfies our requirement for positive infection equilibrium values, $Z^{\ast }>0$.

Next, considering (9(9) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-{T_i}^{\ast }[\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right]=0,\end{array}$(9) ) with $T^{\ast },{T_i}^{\ast },Z^{\ast }$ and ${\eta }^{\ast }$ given by Equations (5(5) $\begin{array}{rl}{T}^{\ast }& =\frac{r{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(d+\rho ({\eta }^{\ast }\left)V^{\ast }\right)},\end{array}$(5) )–(8(8) $\begin{array}{rl}{\eta }^{\ast }& =\frac{{\lambda }_{\eta }}{{\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }}.\end{array}$(8) ), we have (15) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }=\frac{(c+{\mu }_V{\eta }^{\ast })(\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right)}{p}.\end{array}$(15) Let $y_1\left(V\right)=\rho \left(\eta \right)T$ and $y_2\left(V\right)=\frac{(c+{\mu }_V\eta )(\delta +\psi (\eta \left)Z\right)}{p}$. Here, $y_1$ is decreasing on $(0,\widehat{V_1})$, while $y_2$ is increasing on $(0,\widehat{V_1})$ (see Appendix 1). Thus, if a solution exists, it must be unique. We prove the existence of a solution as follows. Since $0<V^{\ast }<\widehat{V_1}$, we set V =0 in $y_1\left(V\right)$ and $y_2\left(V\right)$ to obtain $\begin{array}{rl}{y}_1\left(0\right)& =\frac{r{\rho }_0{\lambda }_{\eta }{\mu }_{\eta }}{\mathrm{d}({\mu }_{\eta }+A{\lambda }_{\eta })}\mathrm{a}\mathrm{n}\mathrm{d}\\ y_2\left(0\right)& =\frac{(c{\mu }_{\eta }+{\mu }_V{\lambda }_{\eta })\left[\delta {\mu }_Z{\mu }_{\eta }\right(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })+{\psi }_0{\lambda }_Z{{\lambda }_{\eta }}^2]}{p{\mu }_Z{\mu }_{\eta }(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })}.\end{array}$Now, (16) $\begin{array}{rl}{\mathfrak{R}}_0& =\frac{pr{\rho }_0{\lambda }_{\eta }{\mu }_Z{{\mu }_{\eta }}^2(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })}{\mathrm{d}({\mu }_{\eta }+A{\lambda }_{\eta })(c{\mu }_{\eta }+{\mu }_V{\lambda }_{\eta })\left[\delta {\mu }_Z{\mu }_{\eta }\right(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })+{\psi }_0{\lambda }_Z{{\lambda }_{\eta }}^2]}>1\\ & ⟹\frac{r{\rho }_0{\lambda }_{\eta }{\mu }_{\eta }}{\mathrm{d}({\mu }_{\eta }+A{\lambda }_{\eta })}>\frac{(c{\mu }_{\eta }+{\mu }_V{\lambda }_{\eta })\left[\delta {\mu }_Z{\mu }_{\eta }\right(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })+{\psi }_0{\lambda }_Z{{\lambda }_{\eta }}^2]}{p{\mu }_Z{\mu }_{\eta }(\mathrm{\Psi }{\mu }_{\eta }+{\lambda }_{\eta })}.\end{array}$(16) That is, $y_1\left(0\right)>y_2\left(0\right)$. Similarly, we set $V=\widehat{V_1}$ to obtain $y_1\left(\widehat{V_1}\right)=0\mathrm{a}\mathrm{n}\mathrm{d}\underset{V⟶\widehat{V_1}}{lim}{y}_2\left(V\right)=\mathrm{\infty }.$Thus, $y_1\left(\widehat{V_1}\right)<y_2\left(\widehat{V_1}\right)$. Since $y_1$ is decreasing on $(0,\widehat{V_1})$, $y_2$ is increasing on $(0,\widehat{V_1})$, $y_1\left(0\right)>y_2\left(0\right)$ and $y_2\left(\widehat{V_1}\right)>y_1\left(\widehat{V_1}\right)$, $y_1\left(V\right)$ and $y_2\left(V\right)$ meet at one point, say, $V^{\ast }$. When ${\mathfrak{R}}_0>1$, Equation (15(15) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }=\frac{(c+{\mu }_V{\eta }^{\ast })(\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right)}{p}.\end{array}$(15) ) has a unique root on $(0,\widehat{V_1})$. Hence, Equation (15(15) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }=\frac{(c+{\mu }_V{\eta }^{\ast })(\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right)}{p}.\end{array}$(15) ) has a unique positive solution $V^{\ast }$. Consequently, the model system (4(4) $\begin{array}{rl}0& =r{\eta }^{\ast }-\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\mathrm{d}{T}^{\ast },\\ 0& =\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-\delta {T_i}^{\ast }-\psi \left({\eta }^{\ast }\right){T_i}^{\ast }{Z}^{\ast },\\ 0& =p{T_i}^{\ast }-c{V}^{\ast }-{\mu }_V{\eta }^{\ast }{V}^{\ast },\\ 0& ={\lambda }_z{\eta }^{\ast }+b{T_i}^{\ast }{Z}^{\ast }-{\mu }_z{Z}^{\ast },\\ 0& ={\lambda }_{\eta }-{\mu }_{\eta }{\eta }^{\ast }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }\end{array}$(4) ) has a unique infection equilibrium, $E_1=(T^{\ast },{T_i}^{\ast },V^{\ast },Z^{\ast },{\eta }^{\ast })$, with $T^{\ast },{T_i}^{\ast },Z^{\ast }$ and ${\eta }^{\ast }$ given by Equations (5(5) $\begin{array}{rl}{T}^{\ast }& =\frac{r{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(d+\rho ({\eta }^{\ast }\left)V^{\ast }\right)},\end{array}$(5) )–(8(8) $\begin{array}{rl}{\eta }^{\ast }& =\frac{{\lambda }_{\eta }}{{\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }}.\end{array}$(8) ) and $V^{\ast }>0$ given by (9(9) $\begin{array}{r}\rho \left({\eta }^{\ast }\right)T^{\ast }{V}^{\ast }-{T_i}^{\ast }[\delta +\psi ({\eta }^{\ast }\left)Z^{\ast }\right]=0,\end{array}$(9) ). This proves the following theorem:

Figure 2. Graphical representation of inequality (14(14) $\begin{array}{rl}& p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })>b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]\\ & \iff b{\lambda }_{\eta }{\mu }_V{V}^{\ast }<({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(p{\mu }_Z-bc{V}^{\ast }).\end{array}$(14) ).

Theorem 3.2

There exists a unique infection equilibrium $E_1=(T^{\ast },{T_i}^{\ast },V^{\ast },Z^{\ast },{\eta }^{\ast })$ of the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) when ${\mathfrak{R}}_0>1$.

It should be noted that (8(8) $\begin{array}{rl}{\eta }^{\ast }& =\frac{{\lambda }_{\eta }}{{\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }}.\end{array}$(8) ) shows that the amount of protein at equilibrium (${\eta }^{\ast }$) is an increasing function of the equilibrium viral load ($V^{\ast }$), which is consistent with results from other studies that have found higher total protein concentrations in people with HIV compared to people without HIV [44,47].

3.4.2. Backward bifurcation

In this subsection, we explore the possibility of backward bifurcation for the within-host model given by (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ). Although we will not discuss the theoretical background and biological implications of the phenomenon, we will derive a necessary and sufficient condition for its existence.

We apply the Center Manifold Method proposed in [56,57] to establish the existence of a backward bifurcation in the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ). The first step involves choosing a bifurcation parameter (ϕ) such that ${\mathfrak{R}}_0\left(\varphi \right)$ is an increasing function of ϕ. In this case, we use p as the bifurcation parameter. It can be verified that ${\mathfrak{R}}_0\left(p\right)$ is an increasing function of p. Let $p\ast$ be the value of p for which ${\mathfrak{R}}_0=1$. That is, (17) $\begin{array}{r}{p}^{\ast }=\frac{(c+{\mu }_V{\eta }_0)(\delta +\psi ({\eta }_0\left)Z_0\right)}{\rho \left({\eta }_0\right)T_0}.\end{array}$(17) When $p=p^{\ast }$, the Jacobian of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) evaluated at the infection-free equilibrium ($J_0=J\left(E_0\right)$) has a simple eigenvalue of zero and all other eigenvalues are negative. Specifically, when $p=p^{\ast }$, the eigenvalues of $J\left(E_0\right)$ are $0,-d,-{\mu }_Z,-{\mu }_{\eta }$, and $-(\delta +\psi ({\eta }_0)Z_0+c+{\mu }_V{\eta }_0)$. The right eigenvector for the simple eigenvalue zero is given by $w=(w_1{w}_2{w}_3{w}_4{w}_5{)}^T$, where $\begin{array}{rl}& w_1=\frac{r}{d}-\frac{\rho \left({\eta }_0\right){\mu }_{\eta }{T}_0}{\mathrm{d}{\gamma }_{\eta }{\eta }_0},w_2=\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0},w_3=\frac{{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0},\\ & w_4=\frac{b{\mu }_{\eta }(c+{\mu }_V{\eta }_0)Z_0}{p^{\ast }{\mu }_Z{\gamma }_{\eta }{\eta }_0}+\frac{{\lambda }_Z}{{\mu }_Z},w_5=1.\end{array}$Observe that, $w_1$ may not be non-negative. However the method can still be applied since $w_1$ corresponds to a positive entry in the infection-free equilibrium [57]. Now, if $v=(v_1{v}_2{v}_3{v}_4{v}_5{)}^T$ is the left eigenvector for the simple eigenvalue zero, then $v_1=0,v_2=\frac{p^{\ast }}{\delta +Z_0\psi \left({\eta }_0\right)},v_3=1,v_4=0,v_5=0.$The left eigenvector has all zero components except $v_2$ and $v_3$. Define, $f_i$ for $i=1,\dots ,5$ by (18) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T=:f_1,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,=:f_2,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V=:f_3,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z=:f_1,\\ \frac{\mathrm{d}\eta }{\mathrm{d}t}& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V=:f_5.\end{array}$(18) First, we calculate $b=\sum _{k,m=1}^t{v}_k{w}_m\frac{{\mathrm{\partial }}^2{f}_k}{\mathrm{\partial }{x}_m\mathrm{\partial }\varphi }(0,0)$, where $x_1=T$, $x_2=T_i$, $x_3=V$, $x_4=Z$ and $x_5=\eta$. Plugging in $v_1=v_4=v_5=0$, $\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_m\mathrm{\partial }p}\equiv 0$ and $\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_m\mathrm{\partial }p}\equiv 1$, for m=2, and $\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_m\mathrm{\partial }p}\equiv 0$, for $m\ne 2$ the sum b becomes (19) $\begin{array}{rl}& b=v_2\sum _{m=1}^5{w}_m\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_m\mathrm{\partial }p}(0,0)+v_3\sum _{m=1}^5{w}_m\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_m\mathrm{\partial }p}(0,0)\\ & =v_3{w}_2\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_2\mathrm{\partial }p}(0,0)=w_2=\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0},\end{array}$(19) which is clearly positive. Next, we calculate $a=\sum _{k,m,n=1}^t{v}_k{w}_m{w}_n\frac{{\mathrm{\partial }}^2{f}_k}{\mathrm{\partial }{x}_m\mathrm{\partial }{x}_n}(0,0)$. The non-zero partial derivatives of $f_2$ and $f_3$ at $(0,0)$ are given below: (20) $\begin{array}{rl}\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_3}& =\rho \left({\eta }_0\right),\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_4}=-\psi \left({\eta }_0\right),\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_5}=-\frac{\mathrm{\Psi }{\psi }_0{Z}_0}{(\mathrm{\Psi }+{\eta }_0{)}^2},\\ \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_3\mathrm{\partial }{x}_5}& =-\frac{{\rho }_{0}A{T}_0}{(1+A{\eta }_0{)}^2},\text{and}\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_3\mathrm{\partial }{x}_5}=-{\mu }_V.\end{array}$(20) Since $v_1=v_4=v_5=0$, the formula for a becomes: (21) $\begin{array}{r}a=v_2\sum _{m,n=1}^t{w}_m{w}_n\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_m\mathrm{\partial }{x}_n}+v_3\sum _{m,n=1}^t{w}_m{w}_n\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_m\mathrm{\partial }{x}_n}.\end{array}$(21) We substitute (20(20) $\begin{array}{rl}\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_1\mathrm{\partial }{x}_3}& =\rho \left({\eta }_0\right),\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_4}=-\psi \left({\eta }_0\right),\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_2\mathrm{\partial }{x}_5}=-\frac{\mathrm{\Psi }{\psi }_0{Z}_0}{(\mathrm{\Psi }+{\eta }_0{)}^2},\\ \frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_3\mathrm{\partial }{x}_5}& =-\frac{{\rho }_{0}A{T}_0}{(1+A{\eta }_0{)}^2},\text{and}\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_3\mathrm{\partial }{x}_5}=-{\mu }_V.\end{array}$(20) ) in (21(21) $\begin{array}{r}a=v_2\sum _{m,n=1}^t{w}_m{w}_n\frac{{\mathrm{\partial }}^2{f}_2}{\mathrm{\partial }{x}_m\mathrm{\partial }{x}_n}+v_3\sum _{m,n=1}^t{w}_m{w}_n\frac{{\mathrm{\partial }}^2{f}_3}{\mathrm{\partial }{x}_m\mathrm{\partial }{x}_n}.\end{array}$(21) ) and simplify the above expression to get: (22) $\begin{array}{rl}a& =-\frac{2{\mu }_V{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}+\left(\frac{2p^{\ast }}{\delta +Z_0\psi \left({\eta }_0\right)}\right)\left\{\frac{r{\mu }_{\eta }}{\mathrm{d}{\gamma }_{\eta }{\eta }_0}\rho \right({\eta }_0)-\frac{{{\mu }_{\eta }}^2{\rho \left({\eta }_0\right)}^2{T}_0}{\mathrm{d}{{\gamma }_{\eta }}^2{{\eta }_0}^2}\\ & -\frac{{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}\left[\frac{{\rho }_{0}A{T}_0}{(1+A{\eta }_0{)}^2}\right]-\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right]\left[\frac{\mathrm{\Psi }{\psi }_0{Z}_0}{(\mathrm{\Psi }+{\eta }_0{)}^2}\right]\\ & -\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right][\frac{b{\mu }_{\eta }(c+{\mu }_V{\eta }_0)Z_0}{p^{\ast }{\mu }_Z{\gamma }_{\eta }{\eta }_0}+\frac{{\lambda }_Z}{{\mu }_Z}]\psi \left({\eta }_0\right)\phantom{\frac{(c^2+{\mu }_V{\eta }_0)}{{\mu }_{\eta }^2}}\}.\end{array}$(22) This leads to the following result:

Theorem 3.3

A backward bifurcation occurs in the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ), if and only if $a>0$, where a is given by (22(22) $\begin{array}{rl}a& =-\frac{2{\mu }_V{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}+\left(\frac{2p^{\ast }}{\delta +Z_0\psi \left({\eta }_0\right)}\right)\left\{\frac{r{\mu }_{\eta }}{\mathrm{d}{\gamma }_{\eta }{\eta }_0}\rho \right({\eta }_0)-\frac{{{\mu }_{\eta }}^2{\rho \left({\eta }_0\right)}^2{T}_0}{\mathrm{d}{{\gamma }_{\eta }}^2{{\eta }_0}^2}\\ & -\frac{{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}\left[\frac{{\rho }_{0}A{T}_0}{(1+A{\eta }_0{)}^2}\right]-\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right]\left[\frac{\mathrm{\Psi }{\psi }_0{Z}_0}{(\mathrm{\Psi }+{\eta }_0{)}^2}\right]\\ & -\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right][\frac{b{\mu }_{\eta }(c+{\mu }_V{\eta }_0)Z_0}{p^{\ast }{\mu }_Z{\gamma }_{\eta }{\eta }_0}+\frac{{\lambda }_Z}{{\mu }_Z}]\psi \left({\eta }_0\right)\phantom{\frac{(c^2+{\mu }_V{\eta }_0)}{{\mu }_{\eta }^2}}\}.\end{array}$(22) ).

Assuming that ${\mu }_V=A=b=\mathrm{\Psi }={\psi }_0=0$, the expression for a (i.e. (22(22) $\begin{array}{rl}a& =-\frac{2{\mu }_V{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}+\left(\frac{2p^{\ast }}{\delta +Z_0\psi \left({\eta }_0\right)}\right)\left\{\frac{r{\mu }_{\eta }}{\mathrm{d}{\gamma }_{\eta }{\eta }_0}\rho \right({\eta }_0)-\frac{{{\mu }_{\eta }}^2{\rho \left({\eta }_0\right)}^2{T}_0}{\mathrm{d}{{\gamma }_{\eta }}^2{{\eta }_0}^2}\\ & -\frac{{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}\left[\frac{{\rho }_{0}A{T}_0}{(1+A{\eta }_0{)}^2}\right]-\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right]\left[\frac{\mathrm{\Psi }{\psi }_0{Z}_0}{(\mathrm{\Psi }+{\eta }_0{)}^2}\right]\\ & -\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right][\frac{b{\mu }_{\eta }(c+{\mu }_V{\eta }_0)Z_0}{p^{\ast }{\mu }_Z{\gamma }_{\eta }{\eta }_0}+\frac{{\lambda }_Z}{{\mu }_Z}]\psi \left({\eta }_0\right)\phantom{\frac{(c^2+{\mu }_V{\eta }_0)}{{\mu }_{\eta }^2}}\}.\end{array}$(22) )) simplifies to $a=\frac{Q{\mu }_{\eta }\rho \left({\eta }_0\right)(r{\gamma }_{\eta }{\eta }_0-T_0{\mu }_{\eta }\rho ({\eta }_0\left)\right)}{\mathrm{d}{{\gamma }_{\eta }}^2{{\eta }_0}^2}$. Furthermore, since all parameters involved in a are positive, $a>0$ if and only if $r{\gamma }_{\eta }{\eta }_0-T_0{\mu }_{\eta }\rho \left({\eta }_0\right)>0$. By substituting ${\eta }_0=\frac{{\lambda }_{\eta }}{{\mu }_{\eta }},T_0=\frac{r{\lambda }_{\eta }}{\mathrm{d}{\mu }_{\eta }}$, and $\rho \left({\eta }_0\right)={\rho }_0$ in (22(22) $\begin{array}{rl}a& =-\frac{2{\mu }_V{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}+\left(\frac{2p^{\ast }}{\delta +Z_0\psi \left({\eta }_0\right)}\right)\left\{\frac{r{\mu }_{\eta }}{\mathrm{d}{\gamma }_{\eta }{\eta }_0}\rho \right({\eta }_0)-\frac{{{\mu }_{\eta }}^2{\rho \left({\eta }_0\right)}^2{T}_0}{\mathrm{d}{{\gamma }_{\eta }}^2{{\eta }_0}^2}\\ & -\frac{{\mu }_{\eta }}{{\gamma }_{\eta }{\eta }_0}\left[\frac{{\rho }_{0}A{T}_0}{(1+A{\eta }_0{)}^2}\right]-\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right]\left[\frac{\mathrm{\Psi }{\psi }_0{Z}_0}{(\mathrm{\Psi }+{\eta }_0{)}^2}\right]\\ & -\left[\frac{(c+{\mu }_V{\eta }_0){\mu }_{\eta }}{p^{\ast }{\gamma }_{\eta }{\eta }_0}\right][\frac{b{\mu }_{\eta }(c+{\mu }_V{\eta }_0)Z_0}{p^{\ast }{\mu }_Z{\gamma }_{\eta }{\eta }_0}+\frac{{\lambda }_Z}{{\mu }_Z}]\psi \left({\eta }_0\right)\phantom{\frac{(c^2+{\mu }_V{\eta }_0)}{{\mu }_{\eta }^2}}\}.\end{array}$(22) ), we obtain the condition $a>0\iff \frac{{\rho }_0{\mu }_{\eta }}{\mathrm{d}{\gamma }_{\eta }}<1$. Hence we have the following corrollary:

Corollary 3.4

Let ${\mu }_V=A=b=\mathrm{\Psi }={\psi }_0=0$ and $\frac{{\rho }_0{\mu }_{\eta }}{\mathrm{d}{\gamma }_{\eta }}<1$. Then the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) exhibits a backward bifurcation.

4. Model parameters and numerical simulation results

4.1. Model parameters

Most of the standard parameters of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) are available in the literature. However, the nutritional parameters remain unclear and hence unknown. This is due to the lack of within-host models that couple the dynamics of HIV and nutrition. Therefore, we use available data to estimate numerical values for the nutritional parameters.

The normal range for total protein in the human body is 60 to 80 grams per litre of blood [58]. Here, we consider an average value of $\frac{60+80}{2}=70$ grams per litre. The amount of blood circulating within a human depends on the human's size and weight. But the average adult has nearly 5 litres of circulating blood [59]. Therefore, an average adult has a total of $70\times 5=350$ grams of protein in their body. Hence, we assume that the equilibrium value of η is ${\eta }^{\ast }=350$ grams of protein.

The estimated lifetime of CD4 T cells is 100 days. Since on average this is the same as $\frac{1}{d}$, d=0.01 per day [50]. We assume that CD4 T cells are produced at a rate of $r_0=10$ cells per μ L of blood per day [50], so that the recruitment rate of CD4 T cells per gram of total protein is $r=\frac{10}{350}$ cells per μ L of blood per day per gram of protein. The infection rate of target cells is taken to be $\text{β}=6.5\times {10}^{-4}\text{μ}\mathrm{L}$ of blood per virion per day [50], so that the infection rate of target cells is ${\rho }_0=6.5\times {10}^{-4}\times 350$ grams of protein μ L of blood per virion per day. The proportionality constant for nutrition, A, is assumed to be 0.1 per gram of protein. Infected cells die at rate $\delta =1$ per day, and the strength of CD8 T cell response is ${\psi }_0=0.42\mu$L of blood per CD8 T cell per day [50]. The half-saturation constant for nutrition, Ψ, represents the density at which half of the equilibrium protein intake is reached and is given by $\mathrm{\Psi }=\frac{1}{2}\times 350=175$ grams of protein. The number of virions produced by an infected CD4 T cell during its lifetime, N, is estimated to be between a few hundred to tens of thousands of virions. Following [50], we take and the average value of N=2000 virions per CD4 cell, so that the virus production rate, p, is given by $p=N\delta =2000$ virions per CD4 cell per day. The viral clearance rate is c=23 per day [50]. The CD8 T cell production rate is given by ${\lambda }_0={10}^{-3}$ CD8 cells per μL of blood per day [50]. Hence, the recruitment rate of CD8 T cells per gram of total protein is given by ${\lambda }_Z=\frac{{10}^{-3}}{350}$ CD8 cells per μL of blood per gram of protein per day. The antigen driven activation rate of CD8 T cells is b=1 per day [50], and the CD8 T cells die at a rate of ${\mu }_Z=2$ per day [50]. The rate of ingestion of total protein, ${\lambda }_{\eta }$, is calculated as the average weight of 180 pounds multiplied by the Dietary Reference Intake (DRI) of 0.36 grams of protein per pound per day [60]. This is approximately 65 grams per day.

The rate of clearance of virus particles by immunoglobulins, ${\mu }_V$, is calculated by solving the third equation of Model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) for ${\mu }_V$. This gives ${\mu }_V=\frac{p{T_i}^{\ast }-c{V}^{\ast }}{{\eta }^{\ast }{V}^{\ast }}$. We assume a typical viral load of $V^{\ast }=100$ virions per μL of blood, and ${T_i}^{\ast }=3$ cells per μL of blood [12]. These estimates, together with p=2000 virions per CD4 cell per day and c=23 per day gives ${\mu }_V=0.1057$ per gram of protein per day. Similarly, solving the fifth equation of (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) at equilibrium we get ${\mu }_{\eta }=\frac{{\lambda }_{\eta }+{\gamma }_{\eta }{\eta }^{\ast }{V}^{\ast }}{{\eta }^{\ast }}$. Assuming that ${\gamma }_{\eta }=0.0001$ per day, we get ${\mu }_{\eta }=1.1857$ per day.

Numerical simulations, as well as predictions of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) depend on the initial conditions. Hence, we estimate the initial conditions that will be used for the numerical simulations here. A normal CD4 T count is from 500 to $\mathrm{1,400}$ cells per cubic millimeter of blood [61]. That is, in $1\mu \mathrm{L}$ of peripheral blood, there are 1000 CD4 T cells. Following [62], we fix the density of target cells before infection (i.e. $T\left(0\right)$), at $1\mathrm{\%}$ of the CD4 T cell density in peripheral blood and assume that there are no infected cells initially (i.e. $T_i\left(0\right)=0$). Based on [62], we take the initial viral load ($V\left(0\right)$), to be ${10}^{-9}$ virions per $\mu \mathrm{L}$ of blood to represent the presence of a small number of virions following primary infection. The normal range for an absolute number of all CD8 T cells, which includes both killer and suppressor T cells, is 150 to 1, 000 cells per cubic millimeter of blood [63]. For simplicity, we assume an average of $\approx 600$ cells per $\mu \mathrm{L}$ of blood. The numerical values of these model parameters and initial conditions are summarized in Table 2.

4.2. Numerical simulation results

In this section the within-host model given by (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) will be simulated using the parameter values presented in Table 2 to assess the impact of some of the important nutritional parameters on the behaviour of protein and viral levels. Primary HIV infection is defined as the first six months after acquiring HIV and is considered to be the most crucial period in terms of disease management [64]. During this period, plasma viremia reaches a peak (usually after 21-28 days) [65]. The intensity and timing of the peak, as well as the equilibrium viral load are important factors that define the trajectory of the disease. Hence, it is important to understand the impact of important model parameters on the viral load, as well as other variables during this period.

4.2.1. Assessing the impact of the viral production rate (p)

The model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) is simulated using the parameter values in Table 2 to explore the impact of the viral production rate (p), on the viral load and total protein. Results of the simulations depicted in Figure 3 show that the equilibrium viral load is strongly dependent on the viral production rate Figure 3(a), which and is consistent with results from previous studies [66]. As expected, increases in the viral production rate will lead to increases in the peak viremia size and total protein, as well as speed up the time for which the viremia and total protein peak. In particular, for the baseline viral production rate of p=2000, the viremia (total protein) peaks on day 31 (34), with a peak size of 2079 virions per μ L and 481 grams of total protein (red curves in Figure 3). Reducing the viral production rate from its baseline value of p=2000 to p=1500 (i.e. by $25\mathrm{\%}$), will lead to a $12.79\mathrm{\%}$ ($3.54\mathrm{\%}$) reduction in peak size of the viral load (total protein) and a 8-day increase in the time that the viremia peaks, as well as a 7-day increase in the time for the total protein to peak (comparing the red and blue curves in Figure 3). At equilibrium, a $29.99\mathrm{\%}$ ($5.33\mathrm{\%}$) reduction in the viral load (total protein) is recorded. However, increasing the viral production rate from its baseline value of p=2000 to p=3000 (i.e. by $50\mathrm{\%}$), will lead to a $14.57\mathrm{\%}$ ($3.95\mathrm{\%}$) increase in the peak size of the viral load (total protein) and a 9-day decrease in the time that the viremia peaks, as well as a 9-day decrease in the time that the total protein peaks (comparing the red and purple curves in Figure 3). For this scenario, a $59.17\mathrm{\%}$ ($15.48\mathrm{\%}$) increase in the equilibrium viral load (total protein) is recorded. In summary, there is a positive correlation between the viral production rate (p) and the viral load and serum protein and a negative correlation between the times at which the viral load and total proteins peak.

Figure 3. The impact of the viral production rate (p) on (a) the viral load and (b) total protein. (a) Simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) depicting the impact of the viral production rate (p) on the viral load. The other parameters used for the simulations are presented in Table 2. (b) Simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) depicting the impact of the viral production rate (p) on total protein. The other parameters used for the simulations are presented in Table 2.

4.2.2. Assessing the impact of the protein intake rate $\left({\lambda }_{\eta }\right)$

The model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) is simulated using the parameter values in Table 2 to investigate the relationship between the dietary protein intake (${\lambda }_{\eta }$) and the viral load and total protein. The results show that increasing the dietary protein intake will trigger an increase in the viral load and total protein level, as well as the time for the viral load and protein to peak Figure 4. Specifically, if all parameters are maintained at their baseline values given in Table 2 (including the baseline dietary intake value of ${\lambda }_{\eta }=65$), the viral load (total protein) will attain its first peak on day 33 (35), with a viral peak size of 2164 virions per μL and a total protein peak size of 531 grams (red curves in Figure 4). If the baseline value of the protein intake rate is reduced by approximately $15\mathrm{\%}$ (i.e. setting ${\lambda }_{\eta }=55$), a $8.60\mathrm{\%}$ reduction from the baseline peak size of the viral load will be recorded (comparing the red and blue curves in 4(a)), while a $18.64\mathrm{\%}$ reduction from the baseline peak size of total protein will be recorded (comparing the red and blue curves in 4(b)). For this scenario, a 3-day reduction in the time for both the viral load and total protein to peak will be recorded, while a $5.81\mathrm{\%}$ ($16.51\mathrm{\%}$) reduction in the equilibrium level of the viral load (total protein) will be recorded. On the other hand, an approximate increase of $31\mathrm{\%}$ in the baseline protein intake rate (i.e. setting ${\lambda }_{\eta }=85$) will result in a $6.98\mathrm{\%}$ increase in the equilibrium viral load and a $33.02\mathrm{\%}$ increase in the total protein levels (comparing the red and purple curves in 4). In summary, for the values of ${\lambda }_{\eta }$ considered in this study, the viral load (total protein) attains a peak between $5^{\mathrm{t}\mathrm{h}}$ and $6^{\mathrm{t}\mathrm{h}}$ ($5^{\mathrm{t}\mathrm{h}}$ and $7^{\mathrm{t}\mathrm{h}}$) weeks of primary infection.

Figure 4. The impact of the protein intake (${\lambda }_{\eta }$) on (a) the viral load and (b) total protein. (a) Simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) depicting the impact of the protein intake (${\lambda }_{\eta }$) on the viral load. The other parameters used for the simulations are presented in Table 2. (b) Simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) depicting the impact of the protein intake (${\lambda }_{\eta }$) on total protein. The other parameters used for the simulations are presented in Table 2.

Additional simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) were carried out to investigate the factors that contribute to the observed positive correlation between the dietary protein intake and serum protein Figure 5(a). For this case, we set the virus production rate (p) to 6000 (which is somewhat unrealistic), the strength of CD8 T response (${\psi }_0$) to $200\mu$L per cell per day, vary the protein intake rate from 55 to 85 and maintain the remaining parameters at their baseline values given in Table 2. As expected, increasing the protein intake rate (${\lambda }_{\eta }$), leads to an increase in the viral load initially Figure 5(a). However, this changes after the third week of primary infection and after the beginning of the third week, increasing ${\lambda }_{\eta }$ is associated with a decreasing viral load Figure 5(a). This implies that there is a negative correlation between the intake rate of total protein and viral load when the strength of the CD8 T response is increased to an extremely high level along with higher viral production rates.

The discussion on the amount of virus present during the primary infection aligns with previous observations: a high level of virus in the bloodstream, which is later controlled by the body's immune system [64]. However, there is one critical observation regarding viral loads and serum protein levels that is unique to this study. Specifically, peak viral loads are almost always quickly followed by peak serum protein levels until equilibrium is attained (see, for example, Figure 5(b)).

Figure 5. (a) The impact of the protein intake (${\lambda }_{\eta }$) on the viral load and (b) the dynamics of the viral load (V) and total protein (η). (a) Simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) illustrating the impact of the protein intake rate (${\lambda }_{\eta }$) on viral load. For these simulations, $p=6000,{\psi }_0=200$, and the other parameter values are as given in Table 2. (b) Simulations of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) illustrating the dynamics of the viral load (V) and total protein (η). Parameter values are as fixed at their baseline values given in Table 2.

4.2.3. Assessing the impact of the enhancement rate of protein by the virus $\left({\gamma }_{\eta }\right)$

In this Section, the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) is simulated to assess the impact of the enhancement rate of protein by the virus (${\gamma }_{\eta }$) on both the viral load and total protein. It should be noted that there is little to no experimental data or research on ${\gamma }_{\eta }$ that could be used to evaluate or approximate ${\gamma }_{\eta }$. This is simply due to the lack of mathematical modelling with regards to the role of protein in HIV. Therefore, for numerical simulations we use an assumed value of 0.0001 per day for ${\gamma }_{\eta }$, the enhancement rate of protein by virus. The results obtained and illustrated in Figure 6 show that the peak viral load (respectively, peak total protein) occur around the same time (on day 33 for the viral load and day 35 for total protein) during primary infection (comparing the primary infection peaks in Figure 6(a,b), respectively). This suggests that the time at which the peak viremia or protein occurs is independent of ${\gamma }_{\eta }$. Furthermore, the simulations show that increasing the enhancement rate of protein by the virus leads to increased equilibrium viral loads and total proteins, although the difference in equilibrium viral load is indiscernible. Interestingly, an increasing enhancement rate of protein by the virus will be associated with decreasing peak viral load. In particular, if the enhancement rate of protein by the virus is increased from its baseline value of ${\gamma }_{\eta }=0.0001$ to ${\gamma }_{\eta }=0.0002$, a decrease of $23.98\mathrm{\%}$ in the peak size of the viral load and a slight increase of $7.62\mathrm{\%}$ in the equilibrium viral load will be observed (comparing the red and yellow curves in Figure 6(a)), while an increase of $27.13\mathrm{\%}$ in the peak total protein levels and an increase of $35.11\mathrm{\%}$ in the equilibrium protein level will be observed (comparing the red and yellow curves in Figure 6(b)). Additional increases from the baseline value of the enhancement rate of protein by the virus will lead to more increases in the equilibrium values of both the viral load and total protein (comparing the red curve with the other curves in Figure 6).

Figure 6. The impact of the protein intake (${\gamma }_{\eta }$) on (a) the viral load and (b) total protein. (a) Time series plot of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) showing the effects of the enhancement rate of protein by the virus (${\gamma }_{\eta }$) on the viral load (V). The values of the other parameter used for the simulations are as given in Table 2. (b) Time series plot of the model (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) showing the effects of the enhancement rate of protein by the virus (${\gamma }_{\eta }$) on the total protein level (η). The values of the other parameter used for the simulations are as given in Table 2.

5. Discussion

In this study, we developed a within-host mathematical model for the dynamics of HIV, which accounts for interactions between HIV, the human immune system, and nutrition. The model is based on the target cell limited model developed in [12]. The first phase involved introducing immune control in this basic model based on the mathematical frameworks in [12,50]. The main novelty of the model framework involves the introduction of nutrition in the form of protein primarily due to its importance in the context of HIV. The goal of this study is to clear some of the misconception regarding HIV and the role of malnutrition in HIV infection and to shed some light on related unanswered questions.

Rigorous mathematical analysis of the model including existence, uniqueness and positivity of solutions was carried out. The immunological reproduction number was computed (using two different approaches) and used to establish standard results related to the existence and stability of an infection-free equilibrium and the existence of endemic equilibria. In particular, it was shown that a unique infection equilibrium exists when the immunological reproduction number is greater than one and that when the immunological reproduction number is less than one, there is a parameter regime within which a backward bifurcation occurs. Additionally, necessary and sufficient conditions for the existence of this backward bifurcation were derived using the centre manifold theory method.

According to a recent study which focused on the nutritional state of HIV infected adults in the United States, the total protein levels of HIV infected individuals were found to be higher compared to their non-infected counterparts [44]. According to the study, non-infected women have, on average, $71.2(\pm 0.1)$ grams of protein per litre of blood, whereas HIV infected women have $79.3(\pm 1.5)$ grams per litre, which is almost a $12\mathrm{\%}$ increase compared to non-infected women. Infected men, on average, have $76.0(\pm 1.3)$ grams of total protein per litre of blood, compared to $72.3(\pm 0.1)$ grams per litre of non-infected men, which is a $5\mathrm{\%}$ boost in total protein levels compared to non-infected men. According to another study in which serum protein electrophoresis was performed on 70 HIV-positive and 42 HIV-negative controls, the control group, on average, has 75.5 grams of protein per litre of blood, whereas the HIV infected group, on average, has 85.4 grams per litre, $13.1\mathrm{\%}$ more than the control group [47]. A similar study on serum protein electrophoresis pattern in patients living with HIV in Iran, however, reveals that, unlikely in the United States, a total protein concentrations of HIV infected individuals were lower than that of the control group, i.e. non-infected individuals [46]. The dietary protein intake of individuals involved in the study mentioned above, Thuppal et al. [44], leads us to some interesting hypothesis. According to the study, the daily protein intake of HIV infected women in the United States, on average, is $83.5(\pm 6.1)$ grams, and this is approximately $16\mathrm{\%}$ higher than the average daily protein intake of non-infected women, which is $71.3(\pm 0.5)$ grams. While this increase in intake of protein can be labelled as being cautious or vigilant about their diets after diagnosed with HIV, it does not necessarily imply that there is any significant correlation between this increased dietary protein intake and increased serum protein levels. Interestingly, our numerical simulations also suggest that the dietary protein intake is in direct correlation with total protein when other parameters of the model are kept constant. This however may not be an accurate portrayal of the situation.

Firstly, the usual scientific consensus is that an increased protein intake does not cause high levels of protein in blood [67]. Whether or not this is true in case of HIV infected individuals remains to be understood. Secondly, the same study discussed previously, Thuppal et al. [44], reveals that while HIV infected men in the United States have total protein levels that are $5\mathrm{\%}$ higher than that of non-infected men, unlike women, HIV infected men do not have an increased dietary protein intake compared to non-infected men. In fact, the daily protein intake of HIV infected men, on average, is $104.9(\pm 6.0)$ grams and that of non-infected men is $104.4(\pm 0.7)$ grams. The difference is almost indiscernible. Finally, the results we obtained with respect to the relationship between dietary protein intake and serum protein levels could simply be due to the lack of complexity of the model. The maintenance of blood protein levels in the human body is a complicated process that involves multiple organs (such as lever, kidneys) and mechanisms working hand-in-hand to regulate protein production, removal and homeostasis in the blood. For example, the clearance of albumin alone is a complicated process which is due to $6\mathrm{\%}$ renal activity, $10\mathrm{\%}$ gastrointestinal activity and $84\mathrm{\%}$ catabolic clearances [68]. In our model, the protein ingestion is considered constant and the removal of protein is assumed to be proportional to the amount of protein in blood, which may not accurately explain the situation. Our model, as insightful as it is, may not be able to capture and account for all these subtle biological details. Now, a question that arises naturally is, ‘what other factors, may they be internal or external, would contribute to the observed increment in total protein levels in HIV infected individuals?’. While we still do not possess a full understanding of what factors cause the total protein levels to go up, and more importantly, how these factors contribute to these high levels of blood protein in HIV infected individuals, it is safe to assume that this is caused, or at least catalyzed, by the presence of the virus in the body.

The increased total protein during HIV infection could be explained, at least to a certain extent, by observing the dynamics between HIV, immune system and nutrition. Due to the infection, there is an increased presence of antibodies (immunoglobulins), and this leads to higher levels of globulin, and consequently, to higher total protein levels in blood. This argument is backed by many studies. According to [47], Immunoglobulin G (IgG), the most common antibody, levels in HIV infected individuals are much higher than that of non-infected individuals. The average IgG level in 70 HIV infected individuals is 27.0 grams per litre, a massive $68.75\mathrm{\%}$ higher than an average of 16.0 grams per litre in 42 non-infected individuals (see Table 1). These facts are backed by many other studies [48,69]. In fact, our results suggest that an increased enhancement rate of protein by the virus, ${\gamma }_{\eta }$, leads to increased equilibrium protein levels. Furthermore, our numerical simulations suggest a positive correlation between ${\gamma }_{\eta }$ and equilibrium viral loads. This is to be expected as a higher ${\gamma }_{\eta }$ value is an indicator of high viral activity in the body.

During the course of our study, numerical simulations were carried out to study the effect of the ingestion rate of protein, ${\lambda }_{\eta }$, on HIV viral loads. Our results suggest that with an increasing rate of ingestion of protein, the viral loads also increase, although the difference in equilibrium viral loads is almost indiscernible. While this is a somewhat complex issue and there seems to be no clear consensus in the scientific literature, consistent with our results, a study on the impact of high protein intake on viral load and hematological parameters in HIV infected patients suggests that excess dietary protein and L-lysine can increase the risk of high HIV replication, subsequent acceleration of immunosuppression and the disease progression [70].

The relationship between p, the viral production rate, and serum protein levels is interesting, and somewhat expected. A study published on the PNAS (Proceedings of the National Academy of Sciences) on HIV replication found that when HIV is supplemented with human serum in vitro, the HIV replication rates were increased [5]. While the protein content in human serum is about $7\mathrm{\%}$, this could be due to a positive correlation between the viral production rate and serum protein. This hypothesis is further confirmed by our numerical simulations.

In conclusion, the viral production rate is in positive correlation with HIV viral loads, as is the enhancement rate of protein by virus. The effect of dietary protein intake on HIV viral loads and serum protein should further be studied with the aid of clinical studies coupled with mathematical modelling in order to draw more accurate conclusions.

Disclosure statement

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

Additional information

Funding

M. Martcheva is supported partially through grant National Science Foundation DMS-1951975.

Appendices

Appendix 1

Here, we establish the boundedness of solutions to the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) under certain assumptions on model parameters. Adding the first, second, third, and fifth equations of (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) we get (A1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }+\frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }+\frac{\mathrm{d}V}{\mathrm{d}\tau }+\frac{\mathrm{d}\eta }{\mathrm{d}\tau }& =r\eta -\mathrm{d}T-\delta T_i-\psi \left(\eta \right)T_{i}Z+p{T}_i-cV-{\mu }_V\eta V+{\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V,\\ & \le r\eta -\mathrm{d}T-\delta T_i+p{T}_i-cV+{\lambda }_{\eta }-{\mu }_{\eta }\eta -({\mu }_V-{\gamma }_{\eta })\eta V.\end{array}$(A1) Assuming that ${\mu }_V>{\gamma }_{\eta }$, (A2) $\begin{array}{r}\frac{\mathrm{d}T}{\mathrm{d}\tau }+\frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }+\frac{\mathrm{d}V}{\mathrm{d}\tau }+\frac{\mathrm{d}\eta }{\mathrm{d}\tau }\le r\eta -\mathrm{d}T-\delta T_i+p{T}_i-cV+{\lambda }_{\eta }-{\mu }_{\eta }\eta .\end{array}$(A2) Now, let $c_1=max\{r,p\}$, and $c_2=min\{d,\delta ,c,{\mu }_{\eta }\}$. Then, (A3) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }+\frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }+\frac{\mathrm{d}V}{\mathrm{d}\tau }+\frac{\mathrm{d}\eta }{\mathrm{d}\tau }& \le {\lambda }_{\eta }+c_1(\eta +T_i)-c_2(T+T_i+V+\eta ),\\ & \le {\lambda }_{\eta }+c_1(T+T_i+V+\eta )-c_2(T+T_i+V+\eta ),\\ & ={\lambda }_{\eta }-(c_2-c_1)(T+T_i+V+\eta ).\end{array}$(A3) Assuming $c_2>c_1$, we get (A4) $\begin{array}{r}\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{s}\mathrm{u}\mathrm{p}(T+T_i+V+\eta )\le \frac{{\lambda }_{\eta }}{c_2-c_1}.\end{array}$(A4) Hence, we have the following:

Theorem A.1

Let $c_1=max\{r,p\}$ and $c_2=min\{d,\delta ,c,{\mu }_V\}$. Assume ${\mu }_V>{\gamma }_{\eta }$ and $c_2>c_1$. Then, the solutions of the system (1(1) $\begin{array}{rl}\frac{\mathrm{d}T}{\mathrm{d}\tau }& =r\eta -\rho \left(\eta \right)TV-\mathrm{d}T,\\ \frac{\mathrm{d}{T}_i}{\mathrm{d}\tau }& =\rho \left(\eta \right)TV-\delta T_i-\psi \left(\eta \right)T_{i}Z,\\ \frac{\mathrm{d}V}{\mathrm{d}\tau }& =p{T}_i-cV-{\mu }_V\eta V,\\ \frac{\mathrm{d}Z}{\mathrm{d}\tau }& ={\lambda }_z\eta +b{T}_{i}Z-{\mu }_{z}Z,\\ \frac{\mathrm{d}\eta }{\mathrm{d}\tau }& ={\lambda }_{\eta }-{\mu }_{\eta }\eta +{\gamma }_{\eta }\eta V.\end{array}$(1) ) are bounded.

Appendix 2

Recall that, $y_1\left(V\right)=\rho \left(\eta \right)T$ and $y_2\left(V\right)=\frac{(c+{\mu }_V\eta )(\delta +\psi (\eta \left)Z\right)}{p}$. First we shall prove that $y_1$ is decreasing on $(0,\widehat{V_1})$. First, observe that the derivative of $\rho \left(V\right)$ is given by ${\rho }^{\prime }\left(V\right)=-\frac{{\rho }_0{\gamma }_{\eta }A{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }V+A{\lambda }_{\eta }{)}^2}<0$. Secondly, T is given by Equation (5(5) $\begin{array}{rl}{T}^{\ast }& =\frac{r{\lambda }_{\eta }}{({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(d+\rho ({\eta }^{\ast }\left)V^{\ast }\right)},\end{array}$(5) ) and a direct calculation shows that the derivative of T is given by $T^{\prime }\left(V\right)=-r{\rho }_0{\lambda }_{\eta }\frac{{\gamma }_{\eta }({\mu }_{\eta }-{\gamma }_{\eta }V)+({\mu }_{\eta }-2{\gamma }_{\eta }V)}{{\left\{\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }\left)\left[{\rho }_0{V}^{\ast }\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+d({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast }+A{\lambda }_{\eta })]\right\}}^2},$which is clearly negative on the interval $(0,\widehat{V_1})$ since $(0,\frac{\widehat{V_1}}{2})\subset (0,\widehat{V_1})$. Therefore, both $\rho \left(\eta \right)$ and T are decreasing, and hence $y_1\left(V\right)$ is decreasing on $(0,\widehat{V_1})$.

Now we shall prove that $y_2\left(V\right)$ is increasing on $(0,\widehat{V_1})$. Since η and $\psi \left(\eta \right)$ are clearly increasing on the given interval, it suffices to show that Z is increasing. Referring back to the discussion following (14(14) $\begin{array}{rl}& p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })>b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]\\ & \iff b{\lambda }_{\eta }{\mu }_V{V}^{\ast }<({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })(p{\mu }_Z-bc{V}^{\ast }).\end{array}$(14) ), it is worth noting that the straight line given by $y=b{\lambda }_{\eta }{\mu }_{V}V$ is increasing and the parabola $y=({\mu }_{\eta }-{\gamma }_{\eta }V)(p{\mu }_Z-bcV)$ is decreasing on $(0,\widehat{V_1})$. Therefore, the denominator of Z in (13(13) $\begin{array}{r}{Z}^{\ast }=\frac{p{\lambda }_Z{\lambda }_{\eta }}{p{\mu }_Z({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })-b{V}^{\ast }\left[c\right({\mu }_{\eta }-{\gamma }_{\eta }{V}^{\ast })+{\lambda }_{\eta }{\mu }_V]},\end{array}$(13) ) is decreasing and hence Z is increasing on the given interval. This completes the proof.