Mathematical Model of HIV-1 Dynamics with T Cell Activation

Obias Mulenga Chimbola^{1,2}

Show more

1. Introduction

The most recent global health observatory (GHO, 2014) data of HIV/AIDS published by WHO [1] proclaims that around 78 million people have been infected with the HIV virus and about 39 million people have died of HIV virus. In the current scenario about 35 million people are living with HIV globally and the majority of the cases are in Sub-Saharan Africa where 70% of the total cases are found. The number of people dying from AIDS-related causes like AIDS-related Kaposi’s sarcoma has been reduced with the increased access to antiretroviral therapy in low and middle income countries. This accentuates the importance of treatment of HIV through antiretroviral therapy.

The human immunodeficiency virus-1 (HIV-1) is a retrovirus containing two single strands of RNA. It is an enveloped virus that carries its own reverse transcriptase, an enzyme, responsible for transcribing viral RNA. HIV-1 is transmitted via three primary routes: unprotected sexual intercourse, mother to child (vertical transmission), and intravenously. HIV targets immune cells that carry a CD4 receptor. These cells include the CD4+ T cells and monocytes/macrophages. The CD4 receptor binds with the protein gp120 (found on the surface of the virus) in initiating cellular infection. This is a necessary step but not sufficient. A co-receptor is needed for proper infection. The virus uses a cytokine receptor on the surface of the host cell near the CD4 receptor to finish the lock and key entry process. It has been shown that individuals lacking the co-receptor have a natural immunity to infection by HIV-1 [2]. Once inside the host cell, the viral DNA can be inserted into the host cell’s genome only upon cellular activation, and new viral particles are created and released.

Many HIV-1 in-host models that seek to elucidate the role of cellular immune response in the pathogenesis of HIV-1 infection have been developed. Some of these models could be reviewed in [3] [4] [5] [6]. Most of these models ignore the activation aspect of the naive T cells. T cells must first get activated to acquire the ability to perform their cytotoxic functions. For instance, CD4+ T cells (helper T cells) are activated upon recognition of their cognate antigen that is bound on major histocompatibility complex class II molecules presented on the surface of an activated antigen presentation cell such as macrophage, dendritic cell etc. For a complete activation process, a second non-specific signal is needed. Once a CD4+ T cell, through a T cell receptor, has recognised its cognate antigen, the next step is to transmit a signal from the surface of the T cell, where the recognition takes place, to the nucleus of the T cell. For the CD4+ T cell to switch from naive state to a state of activation, gene expression must be altered inside the nucleus.

Through the help of activated CD4+ T cells, other immune cells are recruited to the sight of an infection. This is normally done through secretion of chemical messengers (cytokines) by the helper T cells that act in an autocrine or paracrine manner. The activated CD8+ T cells, on the contrary, are only able to recognise their cognate antigen that is loaded on the major histocompatibility complex class I of the activated antigen presenting cells. And the recognition is facilitated through their T cell receptors.

With this background, we are motivated to develop an HIV-1 model that includes the activation of both the naive CD4+ and CD8+ T cells. We assume that the activation process of the CD8+ T cells is enhanced by the action of the activated CD4+ T cells. We further assume that both naive and activated CD4+ T cells are susceptible to HIV-1 infection.

In the current paper, we endeavour to address the following questions.

1) Under what conditions does the infection persistent equilibrium point exist?

2) What conditions must be satisfied for the model equilibria to remain stable?

3) What are the different outcomes in the HIV-1 infection?

Our article is organised as follows: In Section 2, a mathematical model is developed which includes naive T cells, activated T cells, productively infected T cells, and HIV-1 virions. We prove positivity and boundedness of solutions. The equilibria are determined and the basic reproduction number derived. The existence theorem regarding the infection persistent equilibrium is formulated and proved. In Section 3, we establish the local and global stability of the equilibria. Appropriate Lyapunov functionals are formulated in the proof of global stability. In Section 4, we carry out the numerical simulations to validate our mathematical analyses. And lastly, in Section 5 we end with the discussions and conclusion.

2. Model Formulation

In this study, we present a mathematical model of the progression of HIV-1 in the presence of the cellular immune response. We focus on the role of CD4 T cell activation in the pathogenesis of the disease by including the interactions of healthy T cells, activated T cells, HIV-1 infected T cells that are assumed to be productively infected and naive CD8 T cells, activated CD8 T cells or effector cells such as the cytotoxic T lymphocytes (CTLs).

${\stackrel{\dot{}}{T}}_{n}={\pi}_{{T}_{n}}-\beta {T}_{n}{V}_{1}-{d}_{{T}_{n}}{T}_{n}-{\alpha}_{1}{T}_{n}{T}_{p},$ (1)

${\stackrel{\dot{}}{T}}_{a}={\alpha}_{1}{T}_{n}{T}_{p}-\beta {T}_{a}{V}_{1}-{d}_{{T}_{a}}{T}_{a},$ (2)

${\stackrel{\dot{}}{T}}_{p}=\beta \left({T}_{n}+{T}_{a}\right){V}_{1}-{d}_{{T}_{p}}{T}_{p}-\rho {E}_{a}{T}_{p},$ (3)

${\stackrel{\dot{}}{E}}_{n}={\pi}_{{E}_{n}}-{d}_{{E}_{n}}{E}_{n}-{\alpha}_{2}{E}_{n}{T}_{a}{T}_{p},$ (4)

${\stackrel{\dot{}}{E}}_{a}={\alpha}_{2}{E}_{n}{T}_{a}{T}_{p}-{d}_{{E}_{a}}{E}_{a},$ (5)

${\stackrel{\dot{}}{V}}_{1}={N}_{{V}_{1}}{d}_{{T}_{p}}{T}_{p}-{d}_{{v}_{1}}{V}_{1}.$ (6)

Equation (1) describes the dynamics of naive T-cells subjected to HIV-1 infection. In this study we assume that both the naive and activated CD4 T cells have the same infectivity rate though it is now known that it mostly the activated CD4 T cells are susceptible to HIV-1 infection. The first term represents the constant source for these cells from the thymus. The second term represents the rate at which the naive CD4 T cells are getting infected by HIV-1. The third term represents the per capita death rate of these cells and the last term represents the removal rate of activated naive CD4 T cells after presented with their cognate antigens by the antigen presenting cells.

Equation (2) represents a class of activated CD4 T cells that are reduced due to HIV-1 infection. This reduction is represented by the second term. And the last term represents the per capita death rate constant.

Equation (3) represents a class of productively infected CD4 T cells. The first term represents the production of these cells from the infection of naive and activated CD4 T cells by HIV-1. The second term is the bursting rate of these cells as they release the HIV-1 virions. And the last term is the rate at which activated effector cells kill the productively infected CD4 T cells. This is the lytic part of the cellular immune response. We could also consider the nonlytic part in which the cytotoxic T lymphocyte inhibits the viral replication through the release of soluble factors.

Equation (4) represents the naive CD8 T cells. The first term represents the constant source of the naive CD8 T cells from the thymus. The second term is the per capita death of the naive CD8 T cells. The third term represents the activation rate of these cells upon exposure to the cognate antigen with the help of the activated CD4 T cells.

Equation (5) represents a class of activated CD8 T cells. The first term is a variable source from the activated naive CD8 T cells and the last term is the per capita death rate constant.

Equation (6) represents the concentration of HIV-1. The first term is the source of HIV-1 virions arising from the bursting of productively infected CD4 T cells and the last term is the per capita death rate constant.

The model parameters above with their definitions are summarised in Table 1.

3.1. Positivity of Solutions and Boundedness of Solutions

3.1.1. Positivity of Solutions

For the model system above to be immunological meaningful, it is important to prove that all its state variables are non-negative for all times. In other words, solutions of the model above with positive initial data remain positive for all $t>0$.

Table 1. Model parameters and their definitions.

We thus show that our model is well posed in the sense that the populations do not become negative, and the populations are bounded. Also, the non-negative orthant is positively invariant, in other words, any trajectory that begins in the non-negative orthant remains there for all time.

For simplicity of notations, we let

$\begin{array}{l}{x}_{1}\left(t\right):={T}_{n}\left(t\right),\mathrm{}{x}_{2}\left(t\right):={T}_{a}\left(t\right),\mathrm{}{x}_{3}\left(t\right):={T}_{p}\left(t\right),\mathrm{}{x}_{4}\left(t\right):={E}_{n}\left(t\right),\mathrm{}{x}_{5}\left(t\right):={E}_{a}\left(t\right),\\ {x}_{6}\left(t\right):={V}_{1}\left(t\right),\mathrm{}{\pi}_{{d}_{{x}_{1}}}:={\pi}_{{T}_{n}},\mathrm{}{d}_{{x}_{1}}:={d}_{{T}_{n}},\mathrm{}{d}_{{x}_{2}}:={d}_{{T}_{a}},\mathrm{}{d}_{{x}_{3}}:={d}_{{T}_{p}},\mathrm{}{\pi}_{{x}_{4}}:={\pi}_{{E}_{n}},\\ {d}_{{x}_{4}}:={d}_{{E}_{n}},\mathrm{}{d}_{{x}_{5}}:={d}_{{E}_{a}},\mathrm{}{d}_{{x}_{6}}:={d}_{{V}_{1}},\mathrm{}{N}_{{x}_{6}}:={N}_{{V}_{1}}\end{array}$

where := means equal by definition.

Denote by ${\mathbb{R}}_{+}^{6}$ the set of points

${x}_{t}=\left({x}_{1}\mathrm{,}{x}_{2}\mathrm{,}{x}_{3}\mathrm{,}{x}_{4}\mathrm{,}{x}_{5}\mathrm{,}{x}_{6}\right)\in {\mathbb{R}}^{6}$ (7)

with positive coordinates and consider the system with initial values

${x}^{0}=\left({x}_{0}^{1}\mathrm{,}{x}_{0}^{2}\mathrm{,}{x}_{0}^{3}\mathrm{,}{x}_{0}^{4}\mathrm{,}{x}_{0}^{5}\mathrm{,}{x}_{0}^{6}\right)\in {\mathbb{R}}_{+}^{6}$ (8)

Lemma 1. *Let*
${x}_{0}^{i}\ge 0,i=1,2,3,4,5,6$. *Then* *the* *solution*
${x}_{t}\in {\mathbb{R}}_{+}^{6}$ *for* *all*
$t>0$ *in* *the* *region*

${\Gamma}^{6}=\left\{\left({x}_{1}\mathrm{,}{x}_{2}\mathrm{,}{x}_{3}\mathrm{,}{x}_{4}\mathrm{,}{x}_{5}\mathrm{,}{x}_{6}\right)|{x}_{i}\in \mathbb{R}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\forall i\mathrm{,}i=\mathrm{1,2,3,4},\mathrm{5,6}\right\}\subseteq {\mathbb{R}}^{6}\mathrm{.}$

The systems (1)-(6) can be rewritten as follows:

${\stackrel{\dot{}}{x}}_{1}={\pi}_{{x}_{1}}-\beta {x}_{1}{x}_{6}-{d}_{{x}_{1}}{x}_{1}-{\alpha}_{1}{x}_{1}{x}_{3},$ (9)

${\stackrel{\dot{}}{x}}_{2}={\alpha}_{1}{x}_{1}{x}_{3}-\beta {x}_{2}{x}_{6}-{d}_{{x}_{2}}{x}_{2},$ (10)

${\stackrel{\dot{}}{x}}_{3}=\beta \left({x}_{1}+{x}_{2}\right){x}_{6}-{d}_{{x}_{3}}{x}_{3}-\rho {x}_{3}{x}_{5},$ (11)

${\stackrel{\dot{}}{x}}_{4}={\pi}_{{x}_{4}}-{d}_{{x}_{4}}{x}_{4}-{\alpha}_{2}{x}_{2}{x}_{3}{x}_{4},$ (12)

${\stackrel{\dot{}}{x}}_{5}={\alpha}_{2}{x}_{2}{x}_{3}{x}_{4}-{d}_{{x}_{5}}{x}_{5},$ (13)

${\stackrel{\dot{}}{x}}_{6}={N}_{{x}_{6}}{d}_{{x}_{3}}{x}_{3}-{d}_{{x}_{6}}{x}_{6}.$ (14)

*Proof*. Under the given initial conditions, it is easy to prove that the solutions of the system (9)-(14) are positive; if not, we assume a contradiction: that there exists a first time
${t}_{i},\mathrm{}i=1,2,\cdots ,6$ such that

${x}_{i}\left({t}_{i}\right)=0,\mathrm{}{\stackrel{\dot{}}{x}}_{i}\left({t}_{i}\right)<0,\mathrm{}{x}_{j}\left(t\right)\ge 0,0\le t\le {t}_{j},\mathrm{}i\ne j.$ (15)

In the first case, we have

$\stackrel{\dot{}}{x}\left({t}_{1}\right)={\pi}_{{x}_{1}}>0,$ (16)

which is a contradiction. It means that ${x}_{1}\left(t\right)\ge 0,\forall t\ge 0$.

In the second case, we have

$\stackrel{\dot{}}{x}\left({t}_{2}\right)={\alpha}_{1}{x}_{1}\left({t}_{2}\right){x}_{3}\left({t}_{2}\right)\ge 0.$ (17)

which is a contradiction. Hence, ${x}_{2}\left(t\right)\ge 0$ for all $t\ge 0$. In the third case,

$\stackrel{\dot{}}{x}\left({t}_{3}\right)=\beta {x}_{1}\left({t}_{3}\right){x}_{6}\left({t}_{3}\right)+\beta {x}_{2}\left({t}_{3}\right){x}_{6}\left({t}_{3}\right)\ge 0.$ (18)

which is a contradiction. Hence, ${x}_{3}\left(t\right)\ge 0$ for all $t\ge 0$. Analogously, we can show that the rest of the variables are non-negative.

3.1.2. Boundedness of Solutions

We already have a lower bound given by the above lemma 1 for the solutions of model (9)-(14) with nonnegative initial values. We now show that the solutions are bounded from above.

We denote by ${C}_{b}\left(I\right)$ the set of all continuous and bounded functions defined on the interval I taking values in ${\mathbb{R}}^{6}$.

Lemma 2. *Let*
$\phi \mathrm{:}\left[\mathrm{0,}\infty \right)\to {\mathbb{R}}^{6}$ *be* *a* *solution* *of* *system* (9)-(14). *If*
$\phi \left(0\right)\in {\mathbb{R}}_{+}^{6}$, *then*
$\phi \in {C}_{b}\left[\mathrm{0,}\infty \right)$.

*Proof*. Because of Lemma 1, it only remains to prove the existence of an upper bound to the nonnegative solutions of system (9)-(14). Let
$y\left(t\right)={\displaystyle {\sum}_{i=1}^{5}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{x}_{i}$ be the total concentration of the T cells. Then,

$y\left(t\right)\le y\left(0\right)+\frac{{\pi}_{{x}_{1}}+{\pi}_{{x}_{4}}}{{\Lambda}_{k}}\mathrm{,}\text{\hspace{0.17em}}\text{where}\mathrm{}{\Lambda}_{k}=\underset{j}{min}\left\{{d}_{{x}_{j}}\right\}\mathrm{,}j=\mathrm{1,2,3},\mathrm{4,5.}$ (19)

Hence, we obtain

${x}_{6}\left(t\right)\le \frac{{N}_{{x}_{6}}{d}_{{x}_{3}}}{{d}_{{x}_{6}}}\left(y\left(0\right)+\frac{{\pi}_{{x}_{1}}+{\pi}_{{x}_{4}}}{{\Lambda}_{k}}\right)\mathrm{.}$ (20)

3.2. Equilibria

To find the steady states of the model system (9)-(14), we set each differential equation to solve and solve the resulting system simultaneously from which we obtain two equilibria. These are infection free equilibrium, ${\epsilon}^{0}$ and the infection persistent equilibrium ${\epsilon}^{\mathrm{*}}$ in which are the coordinates are strictly positive.

${\epsilon}^{0}=\left({x}_{1}^{0},{x}_{2}^{0},{x}_{3}^{0},{x}_{4}^{0},{x}_{5}^{0},{x}_{6}^{0}\right),\mathrm{}\text{where}\mathrm{}{x}_{1}^{0}=\frac{{\pi}_{{x}_{1}}}{{d}_{{x}_{1}}},\mathrm{}{x}_{4}^{0}=\frac{{\pi}_{{x}_{4}}}{{d}_{{x}_{4}}},\mathrm{}{x}_{i}^{0}=0,\mathrm{}i=2,3,5,6.$ (21)

Derivation of the Basic Reproduction Number, ${\mathcal{R}}_{0}$

The method of Castillo-Chavez et al. [7] is used to compute the reproductive number for the model (9)-(14) and representing the total number of newly infected CD4 T cells arising out of one infected CD4 T cell.

Heterogeneity is defined using groups defined by fixed characteristics:

$(\begin{array}{l}\frac{\text{d}X}{\text{d}t}=f\left(X,Y,Z\right),\\ \frac{\text{d}Y}{\text{d}t}=g\left(X,Y,Z\right),\\ \frac{\text{d}Z}{\text{d}t}=h\left(X,Y,Z\right).\end{array}$ (22)

where $X\in {\mathbb{R}}^{5}\mathrm{,}Y\in \mathbb{R}\mathrm{,}Z\in \mathbb{R}$, and $h\left(X\mathrm{,0,0}\right)=0$. Assuming that $g\left({X}^{\mathrm{*}}\mathrm{,}Y\mathrm{,}Z\right)=0$ implicitly determines a function $Y=\stackrel{\u02dc}{g}\left({X}^{\mathrm{*}}\mathrm{,}Z\right)$. Let $A={D}_{Z}h\left({X}^{\mathrm{*}}\mathrm{,}\stackrel{\u02dc}{g}\left({X}^{\mathrm{*}}\mathrm{,0}\right)\mathrm{,0}\right)$ and further assume that A can be written in the form $A=M-D$, with $M\ge 0$ (that is ${m}_{ij}\ge 0$ ) and $D>0$, a diagonal matrix. The reproduction numbers are then evaluated from the matrix $M{D}^{-1}$.

The cell population subgroups are divided between X non-infected cells, Y infected but not infectious cells and Z infectious class. Then, $X=\left({x}_{1},{x}_{2},{x}_{4},{x}_{5}\right),\mathrm{}Y={x}_{3}$ and $Z={x}_{6}$.

Let ${U}_{0}=\left({X}^{\mathrm{*}}\mathrm{,0,0}\right)$ denote the virus free equilibrium, that is, $f\left({X}^{*},0,0\right)=g\left({X}^{*},0,0\right)=h\left({X}^{*},0,0\right)=0$ and $Y=\stackrel{\u02dc}{g}\left({X}^{\mathrm{*}}\mathrm{,}Z\right)$, where

$\stackrel{\u02dc}{g}\left({X}^{*},Z\right)=\frac{\beta {\pi}_{{x}_{1}}Z+\beta {d}_{{x}_{1}}{x}_{2}^{0}Z}{{d}_{{x}_{1}}\left({d}_{{x}_{1}+\rho {x}_{5}^{0}}\right)}$ (23)

Computing $A={D}_{Z}h\left({X}^{\mathrm{*}}\mathrm{,}\stackrel{\u02dc}{g}\left({X}^{\mathrm{*}}\mathrm{,0}\right)\mathrm{,0}\right)$ gives

$\frac{{N}_{{x}_{6}}\beta {\pi}_{{x}_{1}}}{{d}_{{x}_{1}}}-{d}_{{x}_{6}}\mathrm{,}$ (24)

Therefore,

${\mathcal{R}}_{0}=\frac{{N}_{{x}_{6}}{\pi}_{{x}_{1}}\beta}{{d}_{{x}_{1}}{d}_{{x}_{6}}}\mathrm{.}$ (25)

Let us now state the infection persistent equilibrium, ${\epsilon}^{\mathrm{*}}$.

${\epsilon}^{\mathrm{*}}=\left({x}_{1}^{\mathrm{*}}\mathrm{,}{x}_{2}^{\mathrm{*}}\mathrm{,}{x}_{3}^{\mathrm{*}}\mathrm{,}{x}_{4}^{\mathrm{*}}\mathrm{,}{x}_{5}^{\mathrm{*}}\mathrm{,}{x}_{6}^{\mathrm{*}}\right)\mathrm{,}$ (26)

where

$\{\begin{array}{l}{x}_{1}^{*}=\frac{{\pi}_{{x}_{1}}}{\beta {x}_{6}^{*}+{d}_{{x}_{1}}+{\alpha}_{1}{x}_{3}^{*}},\\ {x}_{2}^{*}=\frac{{\alpha}_{1}{x}_{3}^{*}}{\beta {x}_{6}^{*}+{d}_{{x}_{2}}}\frac{{\pi}_{{x}_{1}}}{\beta {x}_{6}^{*}+{d}_{{x}_{1}}+{\alpha}_{1}{x}_{3}^{*}},\\ {x}_{4}^{*}=\frac{{\Pi}_{{x}_{4}}\left({d}_{{x}_{2}}+\beta {x}_{6}^{*}\right)\left({d}_{{x}_{1}}+{\alpha}_{1}{x}_{3}^{*}+\beta {x}_{6}^{*}\right)}{{d}_{{x}_{4}}\left({d}_{{x}_{2}}+\beta {x}_{6}^{*}\right)\left({d}_{{x}_{1}}+{\alpha}_{1}{x}_{3}^{*}+\beta {x}_{6}^{*}\right)+{\Pi}_{{x}_{1}}{\alpha}_{1}{\alpha}_{2}{x}_{3}^{*2}},\\ {x}_{5}^{*}=\frac{{\Pi}_{{x}_{1}}{\Pi}_{{x}_{4}}{\alpha}_{1}{\alpha}_{2}{x}_{3}^{*2}}{{d}_{{x}_{5}}\left({d}_{{x}_{4}}\left({d}_{{x}_{2}}+\beta {x}_{6}^{*}\right)\left({d}_{{x}_{1}}+{\alpha}_{1}{x}_{3}^{*}+\beta {x}_{6}^{*}\right)+{\Pi}_{{x}_{1}}{\alpha}_{1}{\alpha}_{2}{x}_{3}^{*2}\right)},\\ {x}_{6}^{*}=\frac{{N}_{{x}_{6}}{d}_{{x}_{3}}{x}_{3}^{*}}{{d}_{{x}_{6}}}.\end{array}$ (27)

Theorem 1. (*Existence* *of* *the* *infection* *persistence* *equilibrium*)

Consider the system (9)-(14). The infection persistent equilibrium, ${\epsilon}^{\mathrm{*}}$, exists if ${\mathcal{R}}_{0}>1$, where ${x}_{3}^{\mathrm{*}}$ is positive solution to the quintic equation

${Q}_{1}\left({x}_{3}\right)={A}_{5}{x}_{3}^{5}+{A}_{4}{x}_{3}^{4}+{A}_{3}{x}_{3}^{3}+{A}_{2}{x}_{3}^{2}+{A}_{1}{x}_{3}=0$ (28)

with

$\{\begin{array}{l}{A}_{5}=\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{3}}^{2}{d}_{{x}_{5}}{\phi}_{4}{\phi}_{5}}{{d}_{{x}_{6}}}+\frac{{N}_{{x}_{6}}{\Pi}_{{x}_{1}}{\pi}_{{x}_{4}}{\alpha}_{1}{\alpha}_{2}\beta {d}_{{x}_{3}}\rho {\phi}_{5}}{{d}_{{x}_{6}}}>0,\\ {A}_{4}={d}_{{x}_{5}}{\phi}_{2}{\phi}_{4}+\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{3}}^{2}{d}_{{x}_{5}}{\phi}_{3}{\phi}_{5}}{{d}_{{x}_{6}}}+{\pi}_{{x}_{1}}{\pi}_{{x}_{4}}{\alpha}_{1}{\alpha}_{2}{d}_{{x}_{2}}\rho {\phi}_{5}+\frac{{N}_{{x}_{6}}{\pi}_{{x}_{1}}{\pi}_{{x}_{4}}{\alpha}_{1}{\alpha}_{2}\beta {d}_{{x}_{1}}{d}_{{x}_{3}}\rho}{{d}_{{x}_{6}}},\\ {A}_{3}={d}_{{x}_{5}}{\phi}_{1}{\phi}_{4}+{d}_{{x}_{5}}{\phi}_{2}{\phi}_{3}+{\pi}_{{x}_{1}}{\pi}_{{x}_{4}}{\alpha}_{1}{\alpha}_{2}{d}_{{x}_{1}}{d}_{{x}_{2}}\rho +\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{1}}{d}_{{x}_{2}}{d}_{{x}_{3}}^{2}{d}_{{x}_{4}}{d}_{{x}_{5}}{\phi}_{5}}{{d}_{{x}_{6}}},\\ {A}_{2}={d}_{{x}_{5}}{\phi}_{1}{\phi}_{3}+{d}_{{x}_{1}}{d}_{{x}_{2}}{d}_{{x}_{4}}{d}_{{x}_{5}}{\phi}_{4},\\ {A}_{1}={d}_{{x}_{1}}{d}_{{x}_{2}}{d}_{{x}_{4}}{d}_{{x}_{5}}{\phi}_{1},\end{array}$ (29)

and

$\{\begin{array}{l}{\phi}_{1}={d}_{{x}_{1}}{d}_{{x}_{2}}{d}_{{x}_{3}}\left(1-{\mathcal{R}}_{0}\right)<0,\\ {\phi}_{2}={d}_{{x}_{2}}{d}_{{x}_{3}}{\phi}_{5}-\frac{{N}_{{x}_{6}}^{2}{\pi}_{{x}_{1}}{\beta}^{2}{d}_{{x}_{3}}^{2}}{{d}_{{x}_{6}}^{2}}+\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{1}}{d}_{{x}_{3}}^{2}}{{d}_{{x}_{6}}}-\frac{{N}_{{x}_{6}}{\pi}_{{x}_{1}}{\alpha}_{1}\beta {d}_{{x}_{3}}}{{d}_{{x}_{6}}},\\ {\phi}_{3}={d}_{{x}_{2}}{d}_{{x}_{4}}{\phi}_{5}+\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{1}}{d}_{{x}_{3}}{d}_{{x}_{4}}}{{d}_{{x}_{6}}},\\ {\phi}_{4}={\pi}_{{x}_{1}}{\alpha}_{1}{\alpha}_{2}+\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{3}}{d}_{{x}_{4}}{\phi}_{5}}{{d}_{{x}_{6}}}>0,\\ {\phi}_{5}={\alpha}_{1}+\frac{{N}_{{x}_{6}}\beta {d}_{{x}_{3}}}{{d}_{{x}_{6}}}>0.\end{array}$ (30)

*Proof*. From (28),
${Q}_{1}\left({x}_{3}\right)=0\Rightarrow {x}_{3}^{*}=0$ (which results in the virus free steady state,
${\epsilon}^{0}$ ), or

${Q}_{2}\left({x}_{3}\right)={A}_{5}{x}_{3}^{4}+{A}_{4}{x}_{3}^{3}+{A}_{3}{x}_{3}^{2}+{A}_{2}{x}_{3}^{2}+{A}_{2}{x}_{3}+{A}_{1}=0,$ (31)

To establish the existence of a positive solution of ${Q}_{2}\left({x}_{3}\right)$, we argue as follows.

Note

${Q}_{2}\left(0\right)={A}_{1}={d}_{{x}_{1}}^{2}{d}_{{x}_{2}}^{2}{d}_{{x}_{3}}{d}_{{x}_{4}}{d}_{{x}_{5}}\left(1-{\mathcal{R}}_{0}\right)<0.$ (32)

Due to continuity of ${Q}_{2}\left({x}_{3}\right)$ we note that

$\begin{array}{l}\underset{{x}_{3}\to \infty}{lim}{Q}_{2}\left({x}_{3}\right)=+\infty ,\mathrm{}\text{there}\text{\hspace{0.17em}}\text{exists}\text{\hspace{0.17em}}\text{a}\text{\hspace{0.17em}}\text{positive}\text{\hspace{0.17em}}\text{value}\text{\hspace{0.17em}}\text{of}\text{\hspace{0.17em}}{x}_{3},\text{say},{x}_{3}^{*}\in \left(0,+\infty \right)\\ \text{such}\text{\hspace{0.17em}}\text{that}\text{\hspace{0.17em}}{Q}_{2}\left({x}_{3}^{*}\right)=0.\end{array}$ (33)

Observation 1

$\underset{{x}_{3}^{*}\to 0}{lim}{x}_{1}^{*}=\frac{{\pi}_{{x}_{1}}}{{d}_{{x}_{1}}},\mathrm{}\underset{{x}_{3}^{*}\to \infty}{lim}{x}_{1}^{*}=0,$ (34)

$\underset{{x}_{3}^{*}\to 0}{lim}{x}_{j}^{*}=0,\mathrm{}\underset{{x}_{3}^{*}\to \infty}{lim}{x}_{j}^{*}=0,\mathrm{}j=2,5$ (35)

$\underset{{x}_{3}^{*}\to 0}{lim}{x}_{4}^{*}=\frac{{\pi}_{{x}_{4}}}{{d}_{{x}_{4}}},$ (36)

$\underset{{x}_{3}^{*}\to \infty}{lim}{x}_{4}^{*}=\frac{{\pi}_{{x}_{4}}}{{d}_{{x}_{4}}+\Psi},\mathrm{}\text{where}\mathrm{}\text{\hspace{0.05em}}\Psi =\frac{{\alpha}_{1}{\alpha}_{2}{d}_{{x}_{6}}^{2}{\pi}_{{x}_{1}}}{\beta {N}_{{x}_{6}}{d}_{{x}_{3}}\left(\beta {N}_{{x}_{6}}{d}_{{x}_{3}}+{\alpha}_{1}{d}_{{x}_{6}}\right)}$ (37)

$\underset{{x}_{3}^{*}\to \infty}{lim}{x}_{5}^{*}=\frac{\Psi {\pi}_{{x}_{4}}}{{d}_{{x}_{4}}{d}_{{x}_{5}}+\Psi}.$ (38)

4. Stability Analysis of Equilibria

4.1. Local Stability of Virus Free Equilibrium Point, ${\epsilon}^{0}$

Theorem 2. *For* *system* (9)-(14), *if*
${\mathcal{R}}_{0}<1$, *then* *the* *infection* *free* *equilibrium*,
${\epsilon}^{0}$, *is* *locally* *asymptotically* *stable* *and* *unstable* *otherwise*.

*Proof*. Consider the Jacobian matrix,
$J\left({\epsilon}^{0}\right)$, evaluated at
${\epsilon}^{0}$, of the system (9)-(14).

$J\left({\epsilon}^{0}\right)=\left(\begin{array}{cccccc}-{d}_{{x}_{1}}& 0& -{\alpha}_{1}{x}_{1}^{0}& 0& 0& -\beta {x}_{1}^{0}\\ 0& -{d}_{{x}_{2}}& {\alpha}_{1}{x}_{1}^{0}& 0& 0& 0\\ 0& 0& -{d}_{{x}_{3}}& 0& 0& \beta {x}_{1}^{0}\\ 0& 0& 0& -{d}_{{x}_{4}}& 0& 0\\ 0& 0& 0& 0& -{d}_{{x}_{5}}& 0\\ 0& 0& {N}_{{x}_{6}}{d}_{{x}_{3}}& 0& 0& -{d}_{{x}_{6}}\end{array}\right)$ (39)

The corresponding characteristic equation has the following eigenvalues

${\lambda}_{1}=-{d}_{{x}_{1}},\mathrm{}{\lambda}_{2}=-{d}_{{x}_{2}},\mathrm{}{\lambda}_{3}=-{d}_{{x}_{4}},\mathrm{}{\lambda}_{4}=-{d}_{{x}_{5}}$

and the other two eigenvalues are obtained from the following reduced matrix, ${J}_{1}$

${J}_{1}=\left(\begin{array}{cc}-{d}_{{x}_{3}}& \beta {x}_{1}^{0}\\ {N}_{{x}_{6}}{d}_{{x}_{3}}& -{d}_{{x}_{6}}\end{array}\right)$ (40)

From (40), we investigate the nature of eigenvalues by noting the following.

$Trace\left({J}_{1}\right)=-\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)<0,$ (41)

$Det\left({J}_{1}\right)={d}_{{x}_{3}}{d}_{{x}_{6}}\left(1-{\mathcal{R}}_{0}\right)>0,\mathrm{}\text{since}\mathrm{}\text{\hspace{0.05em}}{\mathcal{R}}_{0}<1.$ (42)

By (41) and (42), we see that all the eigenvalues are negative and thus the infection free equilibrium, ${\epsilon}^{0}$, is locally asymptotically stable.

If ${\mathcal{R}}_{0}>1$, then the corresponding eigenvalues of the sub-matrix (40) are

${\lambda}^{-}=\frac{-\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)-\sqrt{{\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)}^{2}-4{d}_{{x}_{3}}{d}_{{x}_{6}}\left(1-{\mathcal{R}}_{0}\right)}}{2}<0,$ (43)

${\lambda}^{+}=\frac{-\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)+\sqrt{{\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)}^{2}-4{d}_{{x}_{3}}{d}_{{x}_{6}}\left(1-{\mathcal{R}}_{0}\right)}}{2}>0,$ (44)

where the discriminant, $\Delta $, is

$\Delta ={\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)}^{2}-4{d}_{{x}_{3}}{d}_{{x}_{6}}\left(1-{\mathcal{R}}_{0}\right)={\left({d}_{{x}_{3}}-{d}_{{x}_{6}}\right)}^{2}+4{d}_{{x}_{3}}{d}_{{x}_{6}}{\mathcal{R}}_{0}>0.$ (45)

4.2. Global Stability of Infection Free Equilibrium, ${\epsilon}^{0}$

Theorem 3. *Consider* *the* *system* (9)-(14). *The* *infection* *free* *equilibrium*,
${\epsilon}^{0}$, *is* *globally* *asymptotically* *stable* *if*
${\alpha}_{i}=0,\mathrm{}i=1,2$ *and*
${\mathcal{R}}_{0}<1$.

*Proof*. Consider the Lyapunov function, V, defined by

$V=\left({x}_{1}-{x}_{1}^{*}-{x}_{1}^{*}\mathrm{ln}\frac{{x}_{1}}{{x}_{1}^{*}}\right)+{x}_{2}+{x}_{3}+\left({x}_{4}-{x}_{4}^{*}-{x}_{4}^{*}\mathrm{ln}\frac{{x}_{4}}{{x}_{4}^{*}}\right)+{x}_{5}+\frac{\beta {\pi}_{{x}_{1}}}{{d}_{{x}_{1}}{d}_{{x}_{6}}}{x}_{6}.$ (46)

Differentiating with respect to time, t, we get

$\begin{array}{c}\stackrel{\dot{}}{V}=\left(1-\frac{{x}_{1}^{0}}{{x}_{1}}\right)\left[{\pi}_{{x}_{1}}-\beta {x}_{1}{x}_{6}-{d}_{{x}_{1}}{x}_{1}-{\alpha}_{1}{x}_{1}{x}_{3}\right]+{\alpha}_{1}{x}_{1}{x}_{3}-\beta {x}_{2}{x}_{6}-{d}_{{x}_{2}}{x}_{2}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\beta {x}_{1}{x}_{6}+\beta {x}_{2}{x}_{6}-{d}_{{x}_{3}}{x}_{3}-\rho {x}_{3}{x}_{5}+\left(1-\frac{{x}_{4}^{0}}{{x}_{4}}\right)\left[{\pi}_{{x}_{4}}-{d}_{{x}_{4}}{x}_{4}-{\alpha}_{2}{x}_{2}{x}_{3}{x}_{4}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\alpha}_{2}{x}_{2}{x}_{3}{x}_{4}-{d}_{{x}_{5}}{x}_{5}+\frac{\beta {\pi}_{{x}_{1}}}{{d}_{{x}_{1}}{d}_{{x}_{6}}}\left[{N}_{{x}_{6}}{d}_{{x}_{3}}{x}_{3}-{d}_{{x}_{6}}{x}_{6}\right]\end{array}$ (47)

Using the conditions satisfied by system (9)-(14) at ${\epsilon}^{0}$ and rearranging the terms, we obtain

$\begin{array}{l}\stackrel{\dot{}}{V}={d}_{{x}_{1}}{x}_{1}^{0}\left(2-\frac{{x}_{1}^{0}}{{x}_{1}}-\frac{{x}_{1}}{{x}_{1}^{0}}\right)+{d}_{{x}_{4}}{x}_{4}^{0}\left(2-\frac{{x}_{4}^{0}}{{x}_{4}}-\frac{{x}_{4}}{{x}_{4}^{0}}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}+{d}_{{x}_{3}}\left({\mathcal{R}}_{0}-1\right){x}_{3}-{d}_{{x}_{2}}{x}_{2}-\rho {x}_{3}{x}_{5}-{d}_{{x}_{5}}{x}_{5}\le 0.\end{array}$ (48)

Since the arithmetic mean is greater than or equal to the geometrical mean, the first two terms of (48) are less than or equal to zero. And the global stability of ${\epsilon}^{0}$ follows from LaSalle’s Invariance principle [8].

4.3. Local Stability of Infection Persistent Equilibrium, ${\epsilon}^{\ast}$

Theorem 4. *For* *system* (9)-(14), *if*
${\mathcal{R}}_{0}>1$, *then* *the* *infection* *persistent* *equilibrium*,
${\epsilon}^{\mathrm{*}}$, *is* *locally* *asymptotically* *stable* *and* *unstable* *otherwise*.

*Proof*. Consider the Jacobian matrix,
$J\left({\epsilon}^{0}\right)$, evaluated at
${\epsilon}^{0}$, of the system (9)-(14).

$J\left({\epsilon}^{0}\mathrm{,}{\beta}^{\mathrm{*}}\right)=\left(\begin{array}{cccccc}-{d}_{{x}_{1}}& 0& -{\alpha}_{1}{x}_{1}^{0}& 0& 0& -{\beta}^{\mathrm{*}}{x}_{1}^{0}\\ 0& -{d}_{{x}_{2}}& {\alpha}_{1}{x}_{1}^{0}& 0& 0& 0\\ 0& 0& -{d}_{{x}_{3}}& 0& 0& {\beta}^{\mathrm{*}}{x}_{1}^{0}\\ 0& 0& 0& -{d}_{{x}_{4}}& 0& 0\\ 0& 0& 0& 0& -{d}_{{x}_{5}}& 0\\ 0& 0& {N}_{{x}_{6}}{d}_{{x}_{3}}& 0& 0& -{d}_{{x}_{6}}\end{array}\right)$ (49)

with $\beta $ chosen as our bifurcation parameter where $\beta ={\beta}^{*}=\frac{{d}_{{x}_{1}}{d}_{{x}_{6}}}{{N}_{{x}_{6}}{\pi}_{{x}_{1}}}$.

The corresponding characteristic equation of the matrix, (49), has the following eigenvalues

${\lambda}_{1}=-{d}_{{x}_{1}},\mathrm{}{\lambda}_{2}=-{d}_{{x}_{2}},\mathrm{}{\lambda}_{3}=-{d}_{{x}_{4}},\mathrm{}{\lambda}_{4}=-{d}_{{x}_{5}}$

and the other two eigenvalues are obtained from the following reduced matrix, ${J}_{{1}^{\mathrm{*}}}$

${J}_{1}=\left(\begin{array}{cc}-{d}_{{x}_{3}}& {\beta}^{*}{x}_{1}^{0}\\ {N}_{{x}_{6}}& -{d}_{{x}_{6}}\end{array}\right)$ (50)

From which we get the eigenvalues

${\lambda}_{5}=0\text{and}\mathrm{}{\lambda}_{6}=-\left({d}_{{x}_{3}}+{d}_{{x}_{6}}\right)$ (51)

Note that we have a simple zero eigenvalue and the rest of the eigenvalues are negative. Therefore, we can now apply the center manifold theory to prove the local stability of the infection persistent equilibrium, ${\epsilon}^{\mathrm{*}}$. The left eigenvector associated with the zero eigenvalue of $J\left({\epsilon}^{\mathrm{*}}\mathrm{,}{\beta}^{\mathrm{*}}\right)$ is $v=\left({v}_{1}\mathrm{,}{v}_{2}\mathrm{,}{v}_{3}\mathrm{,}{v}_{4}\mathrm{,}{v}_{5}\mathrm{,}{v}_{6}\right)$ where

${v}_{i}=0,\mathrm{}i=1,2,4,5;\mathrm{}{v}_{3}={N}_{{x}_{6}}{x}_{6},\mathrm{}{v}_{6}={v}_{6}>0.$ (52)

and the corresponding right eigenvector is ${w}^{\text{T}}={\left({w}_{1}\mathrm{,}{w}_{2}\mathrm{,}{w}_{3}\mathrm{,}{w}_{4}\mathrm{,}{w}_{5}\mathrm{,}{w}_{6}\right)}^{\text{T}}$ where

$\{\begin{array}{l}{w}_{1}=\frac{-{d}_{{x}_{6}}\left({\alpha}_{1}{\pi}_{{x}_{1}}+{d}_{{x}_{1}}{d}_{{x}_{3}}\right)}{{N}_{{x}_{6}}{d}_{{x}_{1}}^{2}{d}_{{x}_{3}}}{w}_{6},\\ {w}_{2}=\frac{{\alpha}_{1}{\pi}_{{x}_{1}}}{{N}_{{x}_{6}}{d}_{{x}_{1}}{d}_{{x}_{3}}}{w}_{6},\\ {w}_{3}=\frac{{d}_{{x}_{6}}}{{N}_{{x}_{6}}{d}_{{x}_{3}}}{w}_{6},\\ {w}_{j}=0,\mathrm{}j=4,5,\\ {w}_{6}={w}_{6}>0\end{array}$ (53)

According to coefficients a and b defined in Theorem 4.1 of Castillo-Chavez and Song [9], it follows that

$\begin{array}{c}a={\displaystyle \underset{k,i,j=1}{\overset{6}{\sum}}}{v}_{k}{w}_{i}{w}_{j}\frac{{\partial}^{2}{f}_{k}}{{\partial}_{{x}_{i}}{\partial}_{{x}_{j}}}\left({\epsilon}^{0},{\beta}^{*}\right)\\ ={v}_{3}\left[2{w}_{1}{w}_{6}\frac{{\partial}^{2}{f}_{3}}{{\partial}_{{x}_{1}}{\partial}_{{x}_{6}}}+2{w}_{2}{w}_{6}\frac{{\partial}^{2}{f}_{3}}{{\partial}_{{x}_{2}}{\partial}_{{x}_{6}}}\right]\\ =-\frac{2{N}_{{x}_{6}}{\beta}^{*2}{w}_{6}^{2}{x}_{1}^{0}{v}_{6}}{{d}_{{x}_{1}}}<0.\end{array}$ (54)

and

$\begin{array}{c}b={\displaystyle \underset{k,i=1}{\overset{6}{\sum}}}{v}_{k}{w}_{i}\frac{{\partial}^{2}{f}_{k}}{{\partial}_{{x}_{i}}{\partial}_{{\beta}^{*}}}\\ ={N}_{{x}_{6}}{v}_{6}\left[{w}_{1}\frac{{\partial}^{2}{f}_{3}}{{\partial}_{{x}_{1}}{\partial}_{{\beta}^{*}}}+{w}_{2}\frac{{\partial}^{2}{f}_{3}}{{\partial}_{{x}_{2}}{\partial}_{{\beta}^{*}}}+{w}_{3}\frac{{\partial}^{2}{f}_{3}}{{\partial}_{{x}_{3}}{\partial}_{{\beta}^{*}}}+{w}_{6}\frac{{\partial}^{2}{f}_{3}}{{\partial}_{{x}_{6}}{\partial}_{{\beta}^{*}}}\right]\\ ={N}_{{x}_{6}}\left({x}_{1}^{0}+{x}_{2}^{0}\right){v}_{6}{w}_{6}>0.\end{array}$ (55)

We now use the following theorem whose proof is found in Castillo-Chavez and Song [9].

Theorem 5. *Consider* *the* *following* *general* *system* *of* *ordinary* *differential* *equations* *with* *a* *parameter*
$\varphi $

$\frac{\text{d}x}{\text{d}t}=f\left(x\mathrm{,}\varphi \right)\mathrm{,}f\mathrm{:}{\mathbb{R}}^{n}\times \mathbb{R}\to \mathbb{R}\mathrm{}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{and}\mathrm{}f\in {\u2102}^{2}\left({\mathbb{R}}^{n}\times \mathbb{R}\right)\mathrm{,}$ (56)

where 0 is an equilibrium of the system, that is $f\left(\mathrm{0,}\varphi \right)=0$ for all $\varphi $ and assume

A1: $A={D}_{x}f\left(0,0\right)=\frac{\partial {f}_{i}}{\partial {x}_{j}}\left(0,0\right)$ is the linearization of system (9)-(14) around the equilibrium 0 with $\varphi $ evaluated at 0. Zero is a simple eigenvalue of A and other eigenvalues of A have negative real parts;

A2: Matrix A has a right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue. Let
${f}_{k}$ be the k^{th} component of f and

$\begin{array}{l}a={\displaystyle \underset{k,i,j}{\overset{n}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{v}_{k}{w}_{i}{w}_{j}\frac{{\partial}^{2}{f}_{k}}{\partial {x}_{i}\partial {x}_{j}}\left(0,0\right),\\ b={\displaystyle \underset{k,i}{\overset{n}{\sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{v}_{k}{w}_{i}\frac{{\partial}^{2}{f}_{k}}{\partial {x}_{i}\partial \varphi}\left(0,0\right),\end{array}$ (57)

The local dynamics of system (9)-(14) around 0 are totally governed by a and b.

1) $a>0,b>0$. When $\varphi <0$ with $\left|\varphi \right|\ll 1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when $0<\varphi \ll 1$, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium;

2) $a<0,b<0$. When $\varphi <0$ with $\left|\varphi \right|\ll 1$, 0 is unstable; when $0<\varphi \ll 1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium;

3) $a>0,b<0$. When $\varphi <0$ with $\left|\varphi \right|\ll 1$, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0<\varphi \ll 1$, 0 is stable and a positive unstable equilibrium appears;

4) $a<0,b>0$. When $\varphi $ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

Thus, $a<0$ and $b>0$. Using Theorem 5 item (4), we have established the following result which only holds for ${\mathcal{R}}_{0}>1$ but close to 1:

Theorem 6. *The* *unique* *virus* *present* *equilibrium* *state* *guaranteed* *by* *Theorem* 5 *is* *locally* *asymptotically* *stable* *for*
${\mathcal{R}}_{0}$ *near* 1.

4.4. Global Stability of Infection Persistent Equilibrium, ${\epsilon}^{\ast}$

Theorem 7. *Consider* *the* *system* (9)-(14). *If*
${\alpha}_{1}=\rho =0$ *and*
${\mathcal{R}}_{0}>1$, *then* *the* *infection* *persistent* *equilibrium*,
${\epsilon}^{\mathrm{*}}$, *is* *globally* *asymptotically* *stable* (*GAS*).

Define a Lyapunov function, L, as follows

$\begin{array}{c}L=\left({x}_{1}-{x}_{1}^{*}-{x}_{1}^{*}\mathrm{ln}\left(\frac{{x}_{1}}{{x}_{1}^{*}}\right)\right)+\left({x}_{2}-{x}_{2}^{*}-{x}_{2}^{*}\mathrm{ln}\left(\frac{{x}_{2}}{{x}_{2}^{*}}\right)\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left({x}_{3}-{x}_{3}^{*}-{x}_{3}^{*}\mathrm{ln}\left(\frac{{x}_{3}}{{x}_{3}^{*}}\right)\right)+\frac{1}{{N}_{{x}_{6}}}\left({x}_{6}-{x}_{6}^{*}-{x}_{6}^{*}\mathrm{ln}\left(\frac{{x}_{6}}{{x}_{6}^{*}}\right)\right)\end{array}$ (58)

Differentiating along the solutions of system (9)-(14), we get

$L=\left(1-\frac{{x}_{1}^{*}}{{x}_{1}}\right){\stackrel{\dot{}}{x}}_{1}+\left(1-\frac{{x}_{2}^{*}}{{x}_{2}}\right){\stackrel{\dot{}}{x}}_{2}+\left(1-\frac{{x}_{3}^{*}}{{x}_{3}}\right){\stackrel{\dot{}}{x}}_{3}+\frac{1}{{N}_{{x}_{6}}}\left(1-\frac{{x}_{6}^{*}}{{x}_{6}}\right){\stackrel{\dot{}}{x}}_{6},$ (59)

Using the conditions of the infection persistent equilibrium point, ${\epsilon}^{\mathrm{*}}$ we obtain

$\begin{array}{c}L=\left(1-\frac{{x}_{1}^{*}}{{x}_{1}}\right)\left[\beta {x}_{1}^{*}{x}_{6}^{*}+{d}_{{x}_{1}}{x}_{1}^{*}-\beta {x}_{1}{x}_{6}-{d}_{{x}_{1}{x}_{1}}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(1-\frac{{x}_{2}^{*}}{{x}_{2}}\right)\left[-\beta {x}_{2}{x}_{6}+\beta {x}_{2}^{*}{x}_{6}^{*}\frac{{x}_{2}}{{x}_{2}^{*}}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(1-\frac{{x}_{3}^{*}}{{x}_{3}}\right)\left[\beta {x}_{1}{x}_{6}+\beta {x}_{2}{x}_{6}-\beta {x}_{1}^{*}{x}_{6}^{*}\frac{{x}_{3}}{{x}_{3}^{*}}-\beta {x}_{2}^{*}{x}_{6}^{*}\frac{{x}_{3}}{{x}_{3}^{*}}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{{N}_{{x}_{6}}}\left(1-\frac{{x}_{6}^{*}}{{x}_{6}}\right)\left[{N}_{{x}_{6}}{d}_{{x}_{3}}{x}_{3}-{N}_{{x}_{6}}{d}_{{x}_{3}}{x}_{3}^{*}\frac{{x}_{6}}{{x}_{6}^{*}}\right]\end{array}$ (60)

We now implement the following transformation

${y}_{i}=\frac{{x}_{i}}{{x}_{i}^{*}},\mathrm{}i=1,2,\cdots ,6$ (61)

$\begin{array}{c}L=\beta {x}_{1}^{*}{x}_{6}^{*}\left(1-\frac{1}{{y}_{1}}\right)\left(1-{y}_{1}{y}_{6}\right)+{d}_{{x}_{1}}{x}_{1}^{*}\left(1-\frac{1}{{y}_{1}}\right)\left(1-{y}_{1}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\beta {x}_{2}^{*}{x}_{6}^{*}\left(1-\frac{1}{{y}_{2}}\right)\left({y}_{2}-{y}_{2}{y}_{6}\right)+\beta {x}_{1}^{*}{x}_{6}^{*}\left(1-\frac{1}{{y}_{3}}\right)\left({y}_{1}{y}_{6}-{y}_{3}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\beta {x}_{2}^{*}{x}_{6}^{*}\left(1-\frac{1}{{y}_{3}}\right)\left({y}_{2}{y}_{6}-{y}_{3}\right)+{d}_{{x}_{3}}{x}_{3}^{*}\left(1-\frac{1}{{y}_{6}}\right)\left({y}_{3}-{y}_{6}\right)\\ =\beta {x}_{1}^{*}{x}_{6}^{*}\left(3-\frac{1}{{y}_{1}}-\frac{{y}_{3}}{{y}_{6}}-\frac{{y}_{1}{y}_{6}}{{y}_{3}}\right)+{d}_{{x}_{1}}{x}_{1}^{*}\left(2-\frac{1}{{y}_{1}}-{y}_{1}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\beta {x}_{2}^{*}{x}_{6}^{*}\left(1-\frac{{y}_{3}}{{y}_{6}}-\frac{{y}_{2}{y}_{6}}{{y}_{3}}+{y}_{2}\right)\end{array}$ (62)

Since the arithmetic mean is greater than or equal to the geometrical mean, all the three terms of (62) are less than or equal to zero. And the global stability of ${\epsilon}^{\mathrm{*}}$ follows from LaSalle’s Invariance principle [8].

5. Numerical Simulations

To illustrate the analytic results obtained above, we give some simulations using the parameter values in Table 2 and Table 3. Numerical results are displayed in the following figures.

In the following Figure 1 we demonstrate the existence of the virus free equilibrium point,
${\epsilon}^{0}$. This signifies the case when there is no infection present. In other words, it is a scenario where the pathogen that had initially invaded the host has been cleared either through the action of the immune response or treatment (or both). Of course in our case, there is no treatment. So this clearance is all attributed to the immune response. In Figure 2(a) we simulate the dynamics of naive CD4+ T cells and productively infected cells. The naive CD4+ T cells are assumed to have an initial concentration of 10^{6} cells per mL. We also assume that the productively infected CD4+ T cell population is initially zero. We observe from the figure that the naive CD4+ T cell population declines from its maximum value and attains its minimum value after about one and half years

Table 2. Model parameters and their definitions.

Table 3. Model parameters and their definitions.

(a) (b) (c) (d) (e)

Figure 1. Existence of virus free equilibrium point, ${\epsilon}^{0}$. (a) Naive and productively infected CD4+ T cell dynamics; (b) Productively infected and activated CD4+ T cell dynamics; (c) Naive and activated CD8+ T cell dynamics; (d) Activated CD8+ T and productively infected CD4+ T cell dynamics; (e) Naive CD4+ T cell and HIV-1 dynamics.

before reaching a steady state value in about two and half years. This is quite typical of HIV-1 infection for those patients who do not access highly active antiretroviral therapy (HAART).

In the following Figure 2 we demonstrate the existence of an infection persistent equilibrium point, ${\epsilon}^{\mathrm{*}}$. This represents a scenario when the host is invaded by a pathogen, herein referred to as HIV-1, and the infection is established. We simulate our model with ${\mathcal{R}}_{0}=1.1667$, where the basic reproduction number, ${\mathcal{R}}_{0}$, results from utilising the data in Table 1. In Figure 3, we simulate the long term dynamics of HIV-1 infection. In Figure 3(c), naive CD4+ T cells and HIV-1 dynamics are simulated. It is assumed, in this study, that the initial innoculum is 3 virions. It is interesting to observe that after about five years, the viral load reduces to about 200 cells per microlitre. It is at this level when a patient is said to have acquired immune deficiency syndrome (AIDS). It is a condition when a patient’s immunity is weakened to a point where it can hardly fight off any opportunistic infection. Note that in about eight and half years the HIV-1 virus grows unboundedly and the CD4 count declines to zero.

In Figure 3, we simulate the long term dynamics of HIV-1 infection over a period of 3500 days.

6. Discussion and Conclusion

From the analysis of the infection persistent steady state equilibrium point, ${\epsilon}^{\mathrm{*}}$,

(a) (b) (c) (d) (e)

Figure 2. Existence of infection persistent equilibrium point, ${\epsilon}^{\mathrm{*}}$. (a) Naive and productively infected CD4+ T cell dynamics; (b) Productively infected and activated CD4+ T cell dynamics; (c) Naive and activated CD8+ T cell dynamics; (d) Activated CD8+ T and productively infected CD4+ T cell dynamics; (e) Naive CD4+ T cell and HIV-1 dynamics.

(a) (b) (c) (d) (e)

Figure 3. Evolution of HIV infection over a period of 3500 days. (a) Naive and productively infected CD4+ T cell dynamics; (b) Productively infected and activated CD4+ T cell dynamics; (c) Naive CD4+ T cell and HIV-1 dynamics; (d) Naive and activated CD8+ T cell dynamics; (e) Activated CD8+ T and productively infected CD4+ T cell dynamics.

we observe from (34) that when the population of productively infected CD4+ T cells at steady state approaches zero, the population of the uninfected naive CD4+ T cells, in the long run, rebounds to the non-infected steady state. The numerical simulation attests to this observation in Figure 1(a). And as expected when the productively infected CD4+ T cells increase in number substantially, the population of uninfected CD4+ T cells declines to very low levels before attaining its steady state value. This is illustrated in Figure 2(a).

In Equation (35) we note that when the level of infected CD4+ T cells approaches zero, then the steady state values, ${x}_{j}^{*},\mathrm{}j=2,5$ also approach zero due to low levels of antigenic stimulation. Note that when the number of infected CD4+ T cells increases unboundedly it may result in exhaustion which may ultimately impair the activation mechanism of precursors. From Equation (36), we observe that when the steady state value of infected CD4+ T cells approaches zero, the steady state value of the naive CD8+ T cells approaches the non-infected steady state. This is due to the fact that there are now a small number of activated CD4+ T cells available to aid in the activation of these precursors. Hence, there will be less number of activated CD8+ T cells. From (37), we

see that the non-infected steady state value, $\frac{{\pi}_{{x}_{4}}}{{d}_{{x}_{4}}}$, of the naive CD8+ T cells is

decreased due to the activation mechanism. Also, from (37) and (38) we see that when the activation rates approach zero, the steady state of the naive CD8+ T cells rebounds to its non-infected state and that of activated CD8+ T cells approaches zero. Hence, there will be a small number of cytotoxic T cells to kill the infected CD4+ T cells.

Acknowledgements

The author wishes to extend his gratitude to the Botswana International University of Science and Technology (BIUST), Mulungushi University in Zambia and the Simmon’s Foundation for their support. He is also indebted to the reviewers’ constructive comments which helped to improve the manuscript.

References

[1] WHO (2014) Global Health Observatory (GHO) Data. HIV/AIDS, WHO Press, Geneva.

[2] D’Souza, M.P. and Harden, V.A. (1983) Chemokines and HIV-1 Second Receptors: Confluence of Two Fields Generates Optimism in AIDS Research. Nature Medicine, 2, 1293-1300.

https://doi.org/10.1038/nm1296-1293

[3] Wang, Z. and Xu, R. (2012) Stability and Hopf Bifurcation in a Viral Infection Model with Nonlinear Incidence Rate and Delayed Immune Response. Communications in Nonlinear Science and Numerical Simulation, 17, 964-978.

https://doi.org/10.1016/j.cnsns.2011.06.024

[4] Cao, X., Roy, A.K., et al. (2020) Global Dynamics of HIV Infection with Two Disease Transmission Routes—A Mathematical Model. Communications in Mathematical Biology and Neuroscience, 2020, 8.

[5] Dong, Y. and Ma, W. (2012) Global Properties for a Class of Latent HIV Infection Dynamics Model with CTL Immune Response, International Journal of Wavelets, Multiresolution and Information Processing, 10, Article ID: 1250045.

https://doi.org/10.1142/S0219691312500452

[6] Allali, K., Harroudi, S. and Torres, D.F.M. (2018) Analysis and Optimal Control of an Intracellular Delayed HIV Model with CTL Immune Response. Mathematics in Computer Science, 12, 111-127.

https://doi.org/10.1007/s11786-018-0333-9

[7] Castillo-Chavez, C., Feng, Z. and Huang, W. (2002) On the Computation of and Its Role on Global Stability. In: Castillo-Chavez, C., Blower, S., Driessche, P. van D., Kirschner, D. and Yakubu, A.-A., Eds., Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction (Minneapolis, MN, 1999), Vol. 125 of the IMA Volumes in Mathematics and its Applications, Springer, New York, 229-250.

https://doi.org/10.1007/978-1-4757-3667-0_13

[8] Lasalle, J.P. (1976) The Stability of Dynamical Systems. Regional Conference Series in Applied Mathematics, Vol. 25, Society for Industrial and Applied Mathematics.

[9] Castillo-Chavez, C. and Song, B. (2004) Dynamical Models of Tuberculosis and Their Applications. Mathematical Biosciences and Engineering, 1, 361-404.

https://doi.org/10.3934/mbe.2004.1.361

[10] Perelson, A.S. and Kirschner, D.E. (1993) Dynamics of HIV Infection of CD4+ T Cells. Mathematical Biosciences, 114, 181-125.

https://doi.org/10.1016/0025-5564(93)90043-A

[11] Szomolay, B. and Lungu, E.M. (2014) A Mathematical Model for the Treatment of AIDS-Related Kaposi’s Sarcoma. Journal of Biological Systems, 22, 495-522.

https://doi.org/10.1142/S0218339014500247

[12] Kirschner, D.E. and Perelson, A.S. (1995) A Model for the Immune System Response to HIV: AZT Treatment Studies. In: Arino, O., Axelrod, D., Kimmel, M. and Langlais, M., Eds., Mathematical Population Dynamics: Analysis of Heterogeneity and Theory of Epidemics, Wuerz Publishing, Winnipeg, 295.

[13] Hadjiandreou, M.M., et al. (2009) Long-Term HIV Dynamics Subject to Continuous Therapy and Structured Treatment Interruptions. Chemical Engineering Science, 64, 1600-1617.

https://doi.org/10.1016/j.ces.2008.12.010

[14] Wen, Q. and Lou, J. (2009) The Global Dynamics of a Model about HIV-1 Infection In Vivo. Ricerche di Matematica, 58, 77-90.

https://doi.org/10.1007/s11587-009-0048-y