# Stability analysis of a novel delay differential equation model of HIV
infection of CD4+ T-cells

Hoang Anh NGO
LIX, École Polytechnique
Institut Polytechnique de Paris
91120 Palaiseau, France
&Hung Dang NGUYEN
Vietnamese - German Faculty of Medicine (VGFM)
Pham Ngoc Thach University of Medicine
Ho Chi Minh City 710000, Vietnam›
&Mehmet Dik
Department of Mathematics and Computer Science
Beloit College
WI 53511, United States
Corresponding author, Email address<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate a novel 3-compartment model of HIV infection of
CD4+ T-cells with a mass action term by including two versions: one baseline
ODE model and one delay-differential equation (DDE) model with a constant
discrete time delay.

Similar to various endemic models, the dynamics within the ODE model is fully
determined by the basic reproduction term $R_{0}$. If $R_{0}<1$, the disease-
free (zero) equilibrium will be asymptotically stable and the disease
gradually dies out. On the other hand, if $R_{0}>1$, there exists a positive
equilibrium that is globally/orbitally asymptotically stable within the
interior of a predefined region.

To present the incubation time of the virus, a constant delay term $\tau$ is
added, forming a DDE model. In this model, this time delay (of the
transmission between virus and healthy cells) can destabilize the system,
arising periodic solutions through Hopf bifurcation.

Finally, numerical simulations are conducted to illustrate and verify the
results.

_K_ eywords HIV $\cdot$ Globally asymptotical stability $\cdot$ Periodic
solution $\cdot$ Delay term $\cdot$ Steady state

##  1 Inroduction

In the field of epidemiology, although our knowledge of viral dynamics and
virus-specific immmune responses has not fully developed, numerous
mathematical models have been developed an investigated to describe the
immunological response to HIV infection (for example, [5, 14, 8, 28, 54, 43]
and references therein). The models have been used to explain different
phenomena within the host body, and by directly applying the models to real
clinical data, they can also predicts estimates of many measures, including
the death rate of productively infected cells, the rate of viral clearance or
the viral production rate.

These simple HIV models have played an essential role in providing a better
understanding in the dynamics of this infectious diseases, while providing
very important biological meanings for the (combined) drug therapies used
against it. For more references and detailed meta mathematical analysis on
these models in general, we can refer to survey papers written by Kirschner,
1996 [29] or Perelson and Nelson, 1999 [58]

The simplest HIV model, only considering the dynamics of the virus
concentration, is

$\frac{dV}{dt}=P-cV$ (1)

where

  * •

$P$ is an unknown function representing the rate of production of the virus,

  * •

$V$ is the virus concentration.

The dynamics of the population of target cells (CD4+ T-cells for HIV or
hepatic cells for HBV and HCV) is still not fully understood. Nevertheless, a
reasonable, simple model for this population of cells, which can be extended
further in various models, is

$\frac{dT}{dt}=s-dT+aT\left(1-\frac{T}{T_{\max}}\right)$ (2)

with

  * •

$s$ representing the rate at which new T-cells are created from sources within
the body, such as the thymus, or from the proliferation of existing T-cells,

  * •

$d$ being the death rate per T-cells,

  * •

$a$ is the maximum proliferation rate of target T-cells, when the
proliferation is represented by a logistic function, and

  * •

$T_{\max}$ is the population density of T-cells at which proliferation shuts
off.

Human immunodeficiency virus, or HIV, is a virus belonging to the genus
Lentivirus, part of the family Retroviridae [71]. It has an outer envelope of
lipid and viral proteins, which encloses its core. The virion core contains
two positive-sense single-stranded RNA and the enzyme reverse transcriptase,
an RNA-dependent DNA polymerase.

HIV, like most viruses, cannot reproduce by itself. Therefore, they require a
host cell and its materials to replicate. For HIV, it infects a variety of
immune cells, including helper T cells, lymphocytes, monocytes, and dendritic
cells by attaching to a specific receptor called the CD4 receptor contained in
the cell membrane. Along with a chemokine coreceptor, the virus is granted
entry into the cell. Inside the host cell, the viral RNA is transcribed into
DNA by the enzyme reverse transcriptase. However, the enzyme has no
proofreading capacity, so errors often occur during this process, giving rise
to 1 to 3 mutations per newly synthesized virus particle. The DNA provirus is
then transported into the nucleus and inserts itself into the host cell DNA
with the aid of viral integrase. Thus, the viral genetic code becomes a stable
part of the cell genome, which is then transcribed into a full-length mRNA by
the host cell RNA polymerase. The full-length mRNA would be

  1. (a)

the genomes of progeny virus, which would be transported to the cytoplasm for
assembly,

  2. (b)

translated to produce the viral proteins, including reverse transcriptase and
integrase, and

  3. (c)

spliced, creating new translatable sequences

The nonstructural genes on the virus also encode regulatory proteins that have
diverse effects on the host cell, including down-regulating host cell
receptors like CD4 and major histocompatibility complex class I molecules,
aiding in synthesizing full-length HIV RNAs and enabling transportation of the
viral mRNAs out of the nucleus without being spliced by the host cell.
Altogether, these effects enable viral mRNAs to be correctly translated into
polypeptides and packaged into virions. These components are then transported
to the plasma membrane and assembled into the mature virion, exiting the cell.

A person can contract the virus through one of four routes: sexual contact,
either homo- or heterosexual; transfusions with whole blood, plasma, clotting
factors and cellular fractions of blood; contaminated needles; perinatal
transmission. The virus causes tissue destruction, immunodeficiency and can
progress to acquired immunodeficiency syndrome (AIDS), completely breaking
down the human body’s defense mechanisms. These patients are now more
susceptible to infections that should be harmless to a normal person, such as
P.jiroveci pneumonia or tuberculosis, and the conditions are worse as well. So
far, treatments for the disease mainly target reverse transcriptase, viral
proteases, and viral integration and fusion, dealing with the virus infection
before it progresses to AIDS. Currently, one treatment for HIV is highly
active antiretroviral therapy (HAART), which includes a combination of drugs
including nucleoside/nucleotide analog reverse transcriptase inhibitors,
nonnucleoside reverse transcriptase inhibitors, protease inhibitors, fusion
inhibitors, integrase inhibitors, and coreceptor blockers. These drugs are
administered based on individualized criteria such as tolerability, drug-drug
interactions, convenience/adherence, and possible baseline resistance.
Although HAART can lower the viral load, the virus reemerges if the treatment
is stopped. Therefore, HIV infection is currently both chronic and incurable.
[20]

Whenever the population reaches $T_{\max}$, it will decrease, allowing us to
impose an upper constrain $dT_{\max}<s$. With this constrain, the equation (2)
has a unique equilibrium at

$\hat{T}=\frac{T_{\max}}{2a}\left[a-d+\sqrt{(a-d)^{2}+\frac{4as}{T_{\max}}}\right]$
(3)

In 1989, Perelson [56] proposed a general model for the interaction between
the human immune system and HIV; in the same paper, he also simplified that
general model into a simpler model with four compartments, whose dynamics are
described by a system of four ODEs:

  * •

Concentration of cells that are uninfected ($T$),

  * •

Concentration of cells that are latently infected ($T^{*}$),

  * •

Concentration of cells that are actively infected ($T^{**}$), and

  * •

Concentration of free infectious virus particles ($v$).

Later, he extended his own model in Perelson et al. (1993) [55] by proving
various mathematical properties of the model, choosing parameter values from a
restricted set that give rise to the long incubation period characteristic of
HIV infection, and presenting some numerical solutions. He also observed that
his model exhibits many clinical symptoms of AIDS, including:

  * •

Long latency period,

  * •

Low levels of free virus in the environment, and

  * •

Depletion of CD4+ cells.

Time delay, of one type or another, have been incoporated into biological
models in various research papers (for example, [56]); particularly, by the
similar theoretical analysis to dynamical population system (in [50]), they
also play an important role in the dynamical properties of the HIV infection
models. Generally speaking, systems of delay-differential equations (DDEs)
have much more complicated dynamics than that of ordinary differential
equations (ODEs), as the time delay can cause a stable equilibrium of the ODE
system to become unstable, leading to the fluctuation of popolations. In
studying the viral clearance rate, Perelson et al. (1996) [60] stated that
there are two different types of delay that can occur in an HIV infection
model:

  * •

Pharmacological delay: This delay occurs between the ingestion of drug and its
appearance within cells,

  * •

Intracellular delay: This delay happens between the initial HIV infection of a
cell and the release of virions within the environment.

There has also been various attempts by different authors, trying to come up
with the most realistic model by implementing these delays, in one form or
another (constant delay, discrete delay, continuous delay, etc.). For example,

  * •

Herz et al. (1996) [23], who implemented a discrete delay to represent the
intracellular delay in the HIV model. He showed that the incorporation of the
delay would significantly shorten the estimate for the half-life of free virus
particles.

  * •

Mittler et al. (1998) [47] stated that a $\gamma$ \- distribution delay would
be more realistic to describe the intracellular delay, then implemented it to
the original model proposed by Perelson et al. (1996) [60].

  * •

Mittler et al. [46] et al. and Tam et al. (1999) [70] also derived an analytic
expression for the rate of decline of virus following drug treatment by
assuming the drug to be completely effcacious.

  * •

Song and Neumann (2007) [66] proposed a saturated mass-action term into the
simplified model from Perelson et al. (1996) [60], and later investigated the
drug effectiveness under this saturation infection.

The paper will be organized as follows: First, we will investigate a
simplified ODE model from Perelson et al. (1993) [55] by considering three
main components: the uninfected CD4+ T-cells ($T$), the infected CD4+ T-cells
($I$), and the free virus ($V$) with. This model is also assumed to have a
saturation response of the infection rate. Next, the existence and stability
of he infected steady state are considered. Then, we incorporate a discrete
constant delay into the model to indicate the time range between the infection
of a CD4+ T-cell and the emission of viral particles at the cellular level,
resulting in a system of three delay-differential equations (DDEs). To
understand the dynamics of this delay model and obtain sufficient conditions
for local/global asymptotic stability of the equilibria of all time delay, we
carried out a complete analysis on the transcendental characteristic equations
of the linearized system at both the viral-free equilibrium and the infected
(positive) equilibrium. Finally, numerical simulations are carried out, using
Julia, to confirm the obtained results, before some remarks are included in
the conclusion.

##  2 The proposed ODE model

Simplifying the model proposed in Perelson et al. (1993) [55] by reducing the
number of dimensions and assuming that all of the infected cells have the
ability of producing virus at an equal rate, we propose the following epidemic
model of HIV infection of CD4+ T-cells as follows:

$\displaystyle\frac{dT}{dt}$
$\displaystyle=s-dT+aT\left(1-\frac{T}{T_{\max}}\right)-\frac{\beta
TV}{1+\alpha V}+\rho I$ (4) $\displaystyle\frac{dI}{dt}$
$\displaystyle=\frac{\beta TV}{1+\alpha V}-(\delta+\rho)I$
$\displaystyle\frac{dV}{dt}$ $\displaystyle=qI-cV-k_{1}VT$

where

  * •

$T(t)$ is the concentration of healthy CD4+ T-cells at time $t$ (target
cells),

  * •

$I(t)$ is the concentration of infected CD4+ T-cells at time $t$, and

  * •

$V(t)$ is the viral load of the virions (concentration of free HIV at time
$t$).

In infection modelling, it is very common to augment (4) with a "mass-action"
term in which the rate of infection is given by $\beta TV$. This type of term
is sensible, since the virus must interact with T-cells in order to infect and
the probability of virus encountering a T-cell at a low concentration
environment (where infected cells and viral load’s motions are regarded as
independent) can be assumed to be proportional to the product of the density,
which is called linear infection rate. As a result, it follows that the
classical models can assume that T-cells are infected at rate $-\beta TV$ and
are generated at rate $\beta TV$.

With that simple mass-action infection term, the rates of change of uninfected
cells, $T$, productively infected cells $I$, and free virus $V$, would be

$\displaystyle\frac{dT}{dt}$
$\displaystyle=s-dT+aT\left(1-\frac{T}{T_{\max}}\right)-\beta TV$ (5)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=\beta TV-\delta I$
$\displaystyle\frac{dV}{dt}$ $\displaystyle=qI-cV$

Moreover, although the rate of infection in most HIV models is bilinear for
the virus $V$ and the uninfected target cells $T$, the actual incidence rates
are probably not strictly linear for each variable in over the whole valid
range. For example, a non-linear or less-than-linear response in $V$ could
occur due to the saturation at a high enough viral concentration, where the
infectious fraction is significant for exposure to happen very likely. Thus,
is it reasonable to assume that the infection rate of HIV modelling in
saturated mass action is

$\frac{\beta TV^{x}}{1+\alpha V^{y}},\quad x,y,\alpha>0$ (6)

In this paper, we will investigate the viral model with saturation response of
the infection rate where $x=y=1$, for the sake of simplicity. With that being
said, we will proceed to explain the parameters within the model, with

  * •

$s$ is the rate at which new T-cells are created from source from precursors,

  * •

$d$ is the natural death rate of the CD4+ T-cells,

  * •

$a$ is the maximum proliferation rate (growth rate) of T-cells (this means
that $a>d$ in general),

  * •

$T_{\max}$ is the T-cells population density at which proliferation shuts off
(their carrying capacity),

  * •

$\beta$ is the rate constant of infection of T-cells with free virus,

  * •

$\rho$ is the "cure" rate, or the non-cytolytic loss of infected cells,

  * •

$\delta$ is the death rate of the infected cells,

  * •

$q$ is the reproduction rate of the infected cells, and

  * •

$c$ is the clearance rate constant (loss rate) of the virions.

From the explanations above, we can say that

  * •

$\delta+\rho$ is the total rate of disappearance of infected cells from the
environment,

  * •

$1/\delta$ is the average lifespan of a productively infected cell

  * •

$q/\delta$ is the total number of virions produced by an actively infected
cell during its lifespan, and

  * •

$q$ is the average rate of virus released by each cell.

Under the absence of virus (i.e $I(t)=V(t)=0\quad\forall t>0$), the T-cell
population has a steady state value of

$\displaystyle T_{0}$
$\displaystyle=\frac{T_{\max}}{2a}\left[(a-d)+\sqrt{(a-d)^{2}+\frac{4a}{T_{\max}}}\right]$
(7)

The system (4) needs to be initialized with the following initial conditions

$T(0)>0,\quad I(0)>0,\quad V(0)>0,$ (8)

which leads us to denote that

$R_{+}^{3}=\\{(T,I,V)\in\mathbb{R}^{3}\|T\geq 0,I\geq 0,V\geq 0\\}$ (9)

###  2.1 Equilibria and local stability

The system (4) has two steady states: the uninfected steady state
$E_{0}=\left(T_{0},0,0\right)$ and the (positive) infected steady state
$\bar{E}=\left(\bar{T},\bar{I},\bar{V}\right)$, where:

$\displaystyle\bar{T}$
$\displaystyle=\frac{T_{\max}}{2a}\left[a-d-\delta\frac{q\beta-(\delta+\rho)}{q\alpha(\delta+\rho)}+\sqrt{\left(a-d-\delta\frac{q\beta-(\delta+\rho)}{q\alpha(\delta+\rho)}\right)^{2}-\frac{4a}{T_{\max}}\left(\frac{\delta
c}{q\alpha}-s\right)}\right]$ (10) $\displaystyle\bar{I}$
$\displaystyle=\frac{[q\beta-(\delta+\rho)k_{1}]\bar{T}-(\delta+\rho)c}{q\alpha(\delta+\rho)}$
$\displaystyle\bar{V}$
$\displaystyle=\frac{1}{\alpha}\left[\frac{q\beta\bar{T}}{\alpha(\delta+\rho)(c_{1}+k_{1}T}-1\right]$

Now, we will proceed to analyse the stability of the equilibria of system (4).

Since $T_{0}$ and $\bar{T}$ satisfy

$\displaystyle s-dT_{0}+aT_{0}\left(1-\frac{T_{0}}{T_{\max}}\right)$
$\displaystyle=0$ (11) $\displaystyle
s-d\bar{T}+a\bar{T}\left(1-\frac{\bar{T}}{T_{\max}}\right)$
$\displaystyle=\delta\bar{I}=\frac{\delta}{q\alpha(\delta+\rho)}\left[(q\beta-(\delta+\rho))T-(\delta+\rho)c\right]$

we get that

$\bar{T}>\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)k_{1}}\quad\Rightarrow\quad
s-d\bar{T}+a\bar{T}\left(1-\frac{\bar{T}}{T_{\max}}\right)>0\quad\Rightarrow\quad
T_{0}>\bar{T}$ (12)

and

$\bar{T}<\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)k_{1}}\quad\Rightarrow\quad
s-d\bar{T}+a\bar{T}\left(1-\frac{\bar{T}}{T_{\max}}\right)<0\quad\Rightarrow\quad
T_{0}<\bar{T}$ (13)

Hence,

  * •

If $\bar{T}>\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)k_{1}}$, then
$T_{0}>\bar{T}>\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)k_{1}}$, which means
that $E_{0}(T_{0},0,0)$ is unstable, while the positive equilibrium
$\bar{E}(\bar{T},\bar{I},\bar{V})$ exists.

  * •

If $\bar{T}<\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)k_{1}}$, then
$T_{0}<\bar{T}<\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)k_{1}}$, which means
that $E_{0}(T_{0},0,0)$ is locally asymptotically stable, while the positive
equilibrium $\bar{E}(\bar{T},\bar{I},\bar{V})$ is not feasible, as
$\bar{I}<0,\bar{V}<0$.

Let

$R_{0}=\left(\frac{q\beta-(\delta+\rho)k_{1}}{c(\delta+\rho)}\right)\bar{T}$
(14)

We can see that $R_{0}$ is the bifurcation parameter. When $R_{0}<1$, the
uninfected steady state $E_{0}$ is stable and the infected steady state
$\bar{E}$ does not exist (unphysical). When $R_{0}>1$, $E_{0}$ becomes
unstable and $\bar{E}$ exists.

For system (5), it is known that the basic reproductive ratio is given by:

$R_{01}=\left(\frac{q\beta-(\delta+\rho)k_{1}}{c(\delta+\rho)}\right)T_{0}$
(15)

Once again, we emphasize the large difference of the basic reproduction ratio
between the linear infection rate and the saturation infection rate.

  * •

If $\alpha\to 0$, then
$\bar{T}\to\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)},\quad R_{0}\to 1$;

  * •

If $\alpha\rightarrow+\infty$, then $\bar{T}\to T_{0},R_{0}\to R_{01}$.

The Jacobian matrix of system (4) is:

$\left(\begin{matrix}(a-d)-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}&\rho&-\frac{\beta T}{(1+\alpha V)^{2}}\\\ \frac{\beta V}{1+\alpha
V}&-(\delta+\rho)&\frac{\beta T}{(1+\alpha V)^{2}}\\\
-k_{1}V&q&-c-k_{1}T\end{matrix}\right)$ (16)

Let $E^{*}(T^{*},I^{*},V^{*})$ be any arbitrary equilibrium. Then, the
characteristic equation about $E^{*}$ is:

$\left|\begin{matrix}\lambda+\left((d-a)+\frac{2aT^{*}}{T_{\max}}+\frac{\beta
V^{*}}{1+\alpha V^{*}}\right)&-\rho&\frac{\beta T^{*}}{(1+\alpha
V^{*})^{2}}\\\ -\frac{\beta V^{*}}{1+\alpha
V^{*}}&\lambda+(\delta+\rho)&-\frac{\beta T^{*}}{(1+\alpha V^{*})^{2}}\\\
k_{1}V^{*}&-q&\lambda+(c+k_{1}T^{*})\end{matrix}\right|=0$ (17)

For equilibrium $E_{0}=(T_{0},0,0)$, (17) reduces to

$\left(\lambda-a+d+\frac{2aT_{0}}{T_{\max}}\right)\left[\lambda^{2}+(c+\delta+\rho)\lambda+c(\delta+\rho)-q\beta
T_{0}\right]=0$ (18)

Hence, we can see that $E_{0}=(T_{0},0,0)$ is locally asymptotically stable if
$R_{0}<1$, and it is a saddle point if $\dim W^{s}(E_{0})=2$, or if $\dim
W^{s}(E_{0})=1$ while $R_{0}>1$. As a result, we have the following theorems

######  Theorem 1.

If $R_{0}<1$, $E_{0}=(T_{0},0,0)$ is locally asymptotically stable; else, if
$R_{0}>1$, $E_{0}=E_{0}=(T_{0},0,0)$ is unstable.

######  Theorem 2.

There exists $M>0,M\in\mathbb{R}$ such that for any positive solution
$(T(t),I(t),V(t))$ of system (4),

$T(t)\leq M,I(t)\leq M,V(t)\leq M$ (19)

for all large enough $t$.

###### Proof.

Let $L(t)=T(t)+I(t)$ and assume that $L(0)=T(0)+I(0)=\text{const}=c$.
Calculating the derivative of $L(t)$ using the equations in system (4), we
have:

$\displaystyle\frac{dL(t)}{dt}$
$\displaystyle=\frac{dT(t)}{dt}+\frac{dI(t)}{dt}$ (20)
$\displaystyle=s-dT+aT\left(1-\frac{T}{T_{\max}}\right)-\delta I$
$\displaystyle=-dt-\delta
I-\frac{a}{T_{\max}}\left(T-\frac{T_{\max}}{2a}\right)^{2}+\frac{4s+aT_{\max}}{4}$
$\displaystyle\leq-(T+I)\min{(d,\delta)}-\frac{a}{T_{\max}}\left(T-\frac{T_{\max}}{2a}\right)^{2}+\frac{4s+aT_{\max}}{4}$
$\displaystyle=-hL(t)-M_{0}\left(h=\min{(d,\delta)},M_{0}=\frac{4s+aT_{\max}}{4}\right)$

Let $U(t)=L(t)-\frac{M_{0}}{h}$. This means that

$\displaystyle U(0)$ $\displaystyle=L(0)-\frac{M_{0}}{h}=c-\frac{M_{0}}{h}$
(21) $\displaystyle\frac{dU(t)}{dt}$ $\displaystyle=\frac{dL(t)}{dt}$

The inequality (20) can be rewritten as

$\frac{dU(t)}{dt}\leq(-h)U(t)$ (22)

which yields, according to Gronwall’s inequality,

$\displaystyle U(t)$ $\displaystyle\leq
U(0)\exp\left(\int_{0}^{t}(-h)ds\right)$ (23)
$\displaystyle=\left(c-\frac{M_{0}}{h}\right)\exp\left(\left[-hs\right]_{0}^{t}\right)$
$\displaystyle=\left(c-\frac{M_{0}}{h}\right)\exp(-ht)$ $\displaystyle\leq
c-\frac{M_{0}}{h}$

or

$T(t)+I(t)=L(t)=U(t)+\frac{M_{0}}{h}=c-\frac{M_{0}}{h}+\frac{M_{0}}{h}=c$ (24)

As $T(t)>0,I(t)>0\ \forall i\in\mathbb{Z}^{+}$, we can say that

$V(t)\leq c,I(t)\leq c$ (25)

Moreover, we also know that

$\frac{dV}{dt}=qI-cV-k_{1}VT\leq qI-cV\leq qc-cV=-c(V-q)$ (26)

Setting $V(0)=\text{const}=c_{V}$, using the exact same procedure with
Gronwall’s inequality, we obtain

$V(t)\leq c_{V}\ \forall t\in\mathbb{Z}^{+}$ (27)

With $M=\max{(c,c_{V})}$, we would then conclude that

$T(t)\leq M,I(t)\leq M,V(t)\leq M\ \forall t\in\mathbb{Z}^{+}$ (28)

We can easily see that this set is convex. As a consequence, the system (4) is
dissipative.

The proof is complete. ∎

From this theorem, we define

$D=\left\\{(T,I,V)\in\mathbb{R}^{3},0\leq T,I,V\leq M\right\\}.$ (29)

Denote

$M=d-a+\frac{2a\bar{T}}{T_{\max}},\quad
N=\frac{\beta\bar{V}}{1+\alpha\bar{V}},\quad
P=\frac{\beta\bar{T}}{(1+\alpha\bar{V})^{2}}.$ (30)

Then, the characteristic equation of the system around the equilibrium
$\bar{E}(\bar{T},\bar{I},\bar{V})$ reduces to:

$\lambda^{3}+a_{1}\lambda^{2}+(a_{2}+a_{4})\lambda+(a_{3}+a_{5})=0$ (31)

where

$\displaystyle a_{1}$ $\displaystyle=M+(\delta+\rho+c_{1}+k_{1}\bar{T})$ (32)
$\displaystyle a_{2}$
$\displaystyle=(\delta+\rho)(c_{1}+k_{1}T)+M(\delta+\rho+c_{1}+k_{1}\bar{T})+(-k_{1}\bar{V}P)$
$\displaystyle a_{3}$
$\displaystyle=\rho\left[-N(c_{1}+k_{1}\bar{T})+Pk_{1}\bar{V}\right]+PNq$
$\displaystyle a_{4}$ $\displaystyle=-NP$ $\displaystyle a_{5}$
$\displaystyle=M(\delta+\rho)(c_{1}+k_{1}\bar{T})-P(\delta+\rho)k_{1}\bar{V}$

By the Routh-Hurwitz criterion [34], it follows that all eigenvalues of
equation (31) have negative real parts if and only if

$a_{1}>0,\quad a_{3}+a_{5}>0,\quad a_{1}(a_{2}+a_{4})-(a_{3}+a_{5})>0$ (33)

This leads us to the following theorem

######  Theorem 3.

Suppose that

  1. 1.

$R_{0}>1$,

  2. 2.

$a_{1}>0,\quad a_{3}+a_{5}>0,\quad a_{1}(a_{2}+a_{4})-(a_{3}+a_{5})>0$.

Then, the positive equilibrium $\bar{E}(\bar{T},\bar{I},\bar{V})$ is
asymptotically stable.

######  Theorem 4.

If $R_{0}<1$, then $E_{0}(T_{0},0,0)$ is globally asymptotically stable.

###### Proof.

First of all, as $R_{0}<1$, we would have

$T_{0}<\bar{T}<\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)}$ (34)

which means that

$p<\frac{(c+k_{1}T)(\delta+\rho)}{\beta T}$ (35)

From the system (4), we would have

$\displaystyle\frac{dI}{dt}$ $\displaystyle\leq\beta TV-(\delta+\rho)I,$ (36)
$\displaystyle\frac{dV}{dt}$ $\displaystyle=qI-cV-k_{1}VT.$

Now, we would consider the following comparative system

$\displaystyle\frac{dz_{1}}{dt}$ $\displaystyle=\beta
Tz_{2}-(\delta+\rho)z_{1}$ (37) $\displaystyle\frac{dz_{2}}{dt}$
$\displaystyle=pz_{1}-cz_{2}-k_{1}z_{2}T$

We will consider the following form of Lyapunov function:

$L(\mathbf{X})=V(z_{1},z_{2})=\frac{\delta+\rho}{(\beta
T)^{2}}z_{1}^{2}+\frac{1}{c+k_{1}T}z_{2}^{2}$ (38)

The derivative of the function can be calculated as follows

$\displaystyle\frac{dL}{dt}$ $\displaystyle=\frac{\partial L}{\partial
z_{1}}\frac{dz_{1}}{dt}+\frac{\partial L}{\partial z_{2}}\frac{dz_{2}}{dt}$
(39) $\displaystyle=2\frac{\delta+\rho}{(\beta T)^{2}}z_{1}\left(\beta
Tz_{2}-(\delta+\rho)z_{1}\right)+2\frac{1}{c+k_{1}T}z_{2}\left(qz_{1}-cz_{2}-k_{1}Tz_{2}\right)$
$\displaystyle=-2\left[\left(\frac{\delta+\rho}{\beta
T}z_{1}\right)^{2}+z_{2}^{2}-\left(\frac{\delta+\rho}{\beta
T}z_{1}z_{2}+\frac{q}{c+k_{1}T}\right)z_{1}z_{2}\right]$
$\displaystyle\leq-2\left[\left(\frac{\delta+\rho}{\beta
T}z_{1}\right)^{2}+z_{2}^{2}-\left(\frac{\delta+\rho}{\beta
T}+\frac{\beta+\rho}{\beta T}\right)z_{1}z_{2}\right]$
$\displaystyle=-2\left[\frac{\delta+\rho}{\beta T}z_{1}-z_{2}\right]^{2}\leq
0\quad\forall z_{1},z_{2}$

We can see that the derivative is negative definite everywhere except at
$(0,0)$. This means that $(z_{1},z_{2})=(0,0)$ is globally asymptotically
stable.

As we can also see that

$0\leq I(0)\leq z_{1}(0),\quad 0\leq V(0)\leq z_{2}(0)$ (40)

which means that, if the system (37) admits the initial values
$(z_{1}(0),z_{2}(0))$, we have that

$I(t)\leq z_{1}(t),\quad V(t)\leq z_{2}(t)\quad\forall t>t_{1}$ (41)

or, in other words,

$\lim_{t\to+\infty}I(t)=\lim_{t\to+\infty}V(t)=0$ (42)

From this, using the first equation of the system (4), for an $\epsilon\
in(0,1)$ infinitesimal,

$s+(a-d-\delta\epsilon)T-\frac{aT^{2}}{T_{\max}}\leq\frac{dT(t)}{dt}\leq
s+(a-d)T-\frac{aT^{2}}{T_{\max}}\quad\forall t>t_{2}$ (43)

which shows that

$\lim_{t\to+\infty}T(t)=T_{0}.$ (44)

From (42) and (44), we conclude that the system is globally asymptotically
stable. The proof is complete. ∎

######  Theorem 5.

If $R_{0}>1$, then the system (4) is permanent.

###### Proof.

If $R_{0}>1$, we would have

$(q\beta-(\delta+\rho)k_{1})T_{0}>(q\beta-(\delta+\rho)k_{1})\bar{T}>c(\delta+\rho)$
(45)

We will proceed to prove the weak permanence of this system using
contradiction.

Assume that the system is not weakly permanent, from Theorem 4, there exists a
positive orbit $(T(t),I(t),V(t))$ such that

$\lim_{t\to+\infty}T(t)=T_{0},\quad\lim_{t\to+\infty}I(t)=\lim_{t\to+\infty}V(t)=0$
(46)

Since $T_{0}>\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)}$, combining with
(46), we choose an arbitrary infinitesimal $\epsilon>0$ such that there exists
a $t_{0}>0$, for all $t>t_{0}$,

$\displaystyle\frac{T_{0}-\epsilon}{1+\alpha\epsilon}$
$\displaystyle>\frac{c(\delta+\rho)}{q\beta-(\delta+\rho)}$ (47)
$\displaystyle T(t)$ $\displaystyle>T_{0}-\epsilon,$ $\displaystyle V(t)$
$\displaystyle<\epsilon$

Under these conditions, the system (4) becomes

$\displaystyle\frac{dI}{dt}$ $\displaystyle=\frac{\beta TV}{1+\alpha
V}-(\delta+\rho)I\geq\frac{\beta(T_{0}-\epsilon)}{1+\alpha\epsilon}V-(\delta+\rho)I(t)$
(48) $\displaystyle\frac{dV}{dt}$ $\displaystyle=qI-(c_{1}+k_{1}T)\approx qI-
cV-k_{1}T_{0}$

Consider the following Jacobian matrix

$J_{\epsilon}=\left(\begin{matrix}-(\delta+\rho)&\frac{\beta(T_{0}-\epsilon)}{1+\alpha\epsilon}\\\
q&-(c+k_{1}T_{0})\end{matrix}\right)$ (49)

Since $J_{\epsilon}$ has positive off-diagonal element, according to the
Perron - Frobenius theorem, for the maximum positive eigenvalue $j_{1}$ of
$J_{\epsilon}$, there is an associated positive eigenvector
$v=\left(\begin{matrix}v_{1}\\\ v_{2}\end{matrix}\right)$.

Next, we consider a system associated with the Jacobian matrix $J_{\epsilon}$

$\displaystyle\frac{dz_{1}}{dt}$
$\displaystyle=\frac{\beta(T_{0}-\epsilon)}{1+\alpha\epsilon}z_{2}-(\delta+\rho)z_{1}$
(50) $\displaystyle\frac{dz_{2}}{dt}$
$\displaystyle=qz_{1}-(c+k_{1}T_{0})z_{2}$

Let $z(t)=(z_{1}(t),z_{2}(t))$ be a solution of (50) through $(lv_{1},lv_{2})$
at $t=t_{0}$, where $l>0$ satisfies that

$lv_{1}<I(t_{0}),\quad lv_{2}<V(t_{0})$ (51)

As we know that the semi-flow of (50) is monotone and $J_{\epsilon}v=v>0$,
$z_{i}(t)(t=1,2)$ is strictly increasing, meaning $\lim_{t\
to+\infty}z_{i}(t)=+\infty$. This contradicts the Theorem (2), saying that the
positive solution of (4) is bounded from above. This contradiction says that
there exists no positive orbit of (4) tends to $(T_{0},0,0)$ and
$t\to+\infty$. Combining this and a result provided in [11], we conclude that
the system (4) is permanent.

The proof is complete.

∎

######  Theorem 6.

Assume that $D$ is convex and bounded. Suppose that the system

$\frac{dX}{dt}=F(X),\quad X\in D$ (52)

is competitive, permanent and has the property of stability of periodic
orbits. If $\bar{X_{0}}$ is the only equilibrium point in $\text{int}D$ and if
it is locally asymptotically stable, then it is globally asymptotically stable
in $\text{int}{D}$.

###### Proof.

This matrix can easily be proven by considering the Jacobian matrix and choose
the matrix $H$ as

$H=\left(\begin{matrix}1&0&0\\\ 0&-1&0\\\ 0&0&1\end{matrix}\right)$ (53)

By simple calculation, we obtain that

$H\frac{\partial f}{\partial
x}H=\left(\begin{matrix}(a-d)-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}&-\rho&-\frac{\beta T}{(1+\alpha V)^{2}}\\\ -\frac{\beta V}{1+\alpha
V}&-(\delta+\rho)&-\frac{\beta T}{(1+\alpha V)^{2}}\\\
-k_{1}V&-q&-c-k_{1}T\end{matrix}\right)$ (54)

This means that the system (4) is competitive in $D$, with respect to the
partial order defined by the orthant

$K=\left\\{(T,I,V)\in\mathbb{R}^{3}\|T\leq 0,I\geq 0,V\geq 0\right\\}$ (55)

∎

######  Remark 1.

As $D$ is convex and the system (4) is competitive in $D$. we can say that the
system (4) satisfies the Poincare - Bendixson property. This has been proven
by Hirsch (1990) [24], Zhu and Smith (1994) [76] and Smith and Thieme (1991)
[64] that any three-dimensional competitive system that lie in convex sets
would have the Poincaré - Bendixson property; in other words, any non-empty
compact omega limit set that contains no equilibria must be a closed orbit.

######  Theorem 7.

Let $c=I(0)+T(0)$ and suppose that

  1. 1.

$R_{0}>1$,

  2. 2.

$a_{1}>0,a_{3}+a_{5}>0,a_{1}(a_{2}+a_{4})-(a_{3}+a_{5})>0$.

Then, the positive equilibrium $\bar{E}(\bar{T},\bar{I},\bar{V})$ of system
(4) is globally asymptotically stable provided that one of the following two
assumptions hold

  1. 1.

$T_{\max}\frac{a-d+k_{1}c}{2a}<m<T_{0}<T_{\max}\frac{a-d+\delta+k_{1}c}{2a}$,

  2. 2.

$m>T_{\max}\frac{a-d+\delta+k_{1}c}{2a}$.

As we have already known that the system (4) is competitive and permanent
(from Theorem 5 and Theorem 6), while $\bar{E}(\bar{T},\bar{I},\bar{V})$ is
locally asymptotically stable if the two properties (i) and (ii) of Theorem 7
holds. As a result, in accordance with Theorem 6 (choosing $D=\Omega$),
Theorem 7 if we can prove that the system (4) has the stability of periodic
orbits. We will proceed to prove this under the following proposition.

######  Proposition 1.

Assume condition (iii) or (iv) of 7 hold true. Then, system (4) has the
property of stability of periodic orbits.

###### Proof.

Let $P(t)=((T(t),I(t),V(t))$ be a periodic solution whose orbit $\Gamma$ is
contained in $\text{int}\Omega$. In accordance with the criterion given by
Muldowney in [48], for the asymptotic orbital stability of a periodic orbit of
a general autonomous system, it is sufficient to prove that the linear non-
autonomous system

$\frac{dW(t)}{dt}=\left(DF^{[2]}\left(P(t)\right)\right)W(t)$ (56)

is asymptotically stable, where $DF^{[2]}$ is the second additive compound
matrix of the Jacobian $DF$ [66].

The Jacobian matrix of the system (4) is given by

$J=\left(\begin{matrix}(a-d)-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}&\rho&-\frac{\beta T}{(1+\alpha V)^{2}}\\\ \frac{\beta V}{1+\alpha
V}&-(\delta+\rho)&\frac{\beta T}{(1+\alpha V)^{2}}\\\
-k_{1}V&q&-(c+k_{1}T)\end{matrix}\right)$ (57)

For the solution $P(t)$, the equation (56) becomes

$\displaystyle\frac{dW_{1}}{dt}$
$\displaystyle=-\left(\delta+\rho-(a-d)+\frac{2aT}{T_{\max}}+\frac{\beta
V}{1+\alpha V}\right)W_{1}+\frac{\beta T}{(1+\alpha V)^{2}}(W_{2}+W_{3}),$
(58) $\displaystyle\frac{dW_{2}}{dt}$
$\displaystyle=qW_{1}+\left(a-d-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}-(c+k_{1}T)\right)W_{2}+\rho W_{3},$ $\displaystyle\frac{dW_{3}}{dt}$
$\displaystyle=k_{1}VW_{1}+\frac{\beta V}{1+\alpha
V}W_{2}-(\delta+\rho+c+k_{1}T)W_{3}.$

To prove that the system (58) is asymptotically stable, we shall use the
following Lyapunov function, which is similar to the one found in [38] for the
SEIR model:

$L(W_{1}(t),W_{2}(t),W_{3}(t),T(t),I(t),V(t))=\left\lVert\left(W_{1}(t),\frac{I(t)}{V(t)}W_{2}(t),\frac{I(t)}{V(t)}W_{3}(t)\right)\right\rVert$
(59)

where $\left\lVert\cdot\right\rVert$ is the norm in $\mathbb{R}^{3}$ defined
by

$\left\lVert(W_{1},W_{2},W_{3})\right\rVert=\sup\\{|W_{1}|,|W_{2}+W_{3}|\\}$
(60)

From Theorem 5, we obtain that the orbit of $P(t)$ remains at a positive
distance from the boundary of $\Omega$. Therefore,

$I(t)\geq\eta,\quad
V(t)\geq\eta,\quad\eta=\min\\{\underline{I},\underline{V}\\}\quad\forall
t\to+\infty$ (61)

Hence, the function $L(t)$ is well defined along $P(t)$ and

$L(W_{1},W_{2},W_{3};T,I,V)\geq\frac{\eta}{M}\left\lVert(W_{1},W_{2},W_{3})\right\rVert$
(62)

Along a solution $(W_{1},W_{2},W_{3})$ of the system (58), $L(t)$ becomes

$L(t)=\sup\left\\{\left|W_{1}(t)\right|,\frac{I(t)}{V(t)}\left(\left|W_{2}(t)\right|+\left|W_{3}(t)\right|\right)\right\\}$
(63)

Then, we would have the following inequalities

$\displaystyle D_{+}|W_{1}(t)|$
$\displaystyle\leq-\left(\delta+\rho-(a-d)+\frac{2aT}{T_{\max}}+\frac{\beta
V}{1+\alpha V}\right)|W_{1}|+\frac{\beta T}{(1+\alpha
V)^{2}}\left(|W_{2}(t)|+|W_{3}(t)|\right)$ (64) $\displaystyle
D_{+}|W_{2}(t)|$ $\displaystyle\leq
q|W_{1}(t)|+\left(a-d-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}-(c+k_{1}T)\right)|W_{2}(t)|+\rho|W_{3}(t)|$ $\displaystyle D_{+}|W_{3}(t)|$
$\displaystyle\leq k_{1}V|W_{1}(t)|+\frac{\beta V}{1+\alpha
V}|W_{2}(t)|-(\delta+\rho+c+k_{1}T)|W_{3}(t)|$

From this, we get

$\displaystyle D_{+}\frac{I}{V}(|W_{2}|+|W_{3}|)$
$\displaystyle=\left(\frac{dI/dt}{V}-\frac{IdV/dt}{V^{2}}\right)(|W_{2}|+|W_{3}|)+\frac{I}{V}D_{+}(|W_{2}|+|W_{3}|)$
(65)
$\displaystyle\leq\left(\frac{dI/dt}{I}-\frac{dV/dt}{V}\right)\frac{I}{V}(|W_{2}|+|W_{3}|)+\left(\frac{qI}{V}+k_{1}I\right)|W_{1}|$
$\displaystyle-\left(-a+d+\frac{2aT}{T_{\max}}+(c+k_{1}T)\right)\frac{I}{V}|W_{2}(t)|-(\delta+c+k_{1}T)\frac{I}{V}|W_{3}(t)|$

Thus, we can obtain

$D_{+}L(t)\leq\sup\\{g_{1}(t),g_{2}(t)\\}L(t),$ (66)

where

$\displaystyle g_{1}(t)$
$\displaystyle=-\delta-\rho+a-d-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}+\frac{\beta TV}{I(1+\alpha V)^{2}}$ (67) $\displaystyle g_{2}(t)$
$\displaystyle=\frac{qI}{V}+k_{1}I+\frac{dI/dt}{I}-\frac{dV/dt}{V}-G_{1}$
$\displaystyle G_{1}$
$\displaystyle=\min\left\\{-a+d+\frac{2aT}{T_{\max}}+(c+k_{1}T),\delta+c+k_{1}T\right\\}$

From the second equation of the system (4), we obtain

$\displaystyle g_{1}(t)$
$\displaystyle=-\delta-\rho+a-d-\frac{2aT}{T_{\max}}-\frac{\beta V}{1+\alpha
V}+\frac{\beta TV}{I(1+\alpha V)^{2}}$ (68)
$\displaystyle\leq-\delta-\rho+a-d-\frac{2aT}{T_{\max}}-\frac{\beta
V}{1+\alpha V}+\frac{\beta tV}{I(1+\alpha V)}$
$\displaystyle=a-d-\frac{2aT}{T_{\max}}-\frac{\beta T}{1+\alpha
V}+\frac{dI/dt}{I}$

Here, we consider two different cases.

  * •

Case 1: If Point 3 of Theorem 7 holds, then

$-\delta<a-d-\frac{2aT}{T_{\max}}<0,$ (69)

that is

$G_{1}=-a+d+\frac{2aT}{T_{m}ax}+(c+k_{1}T)$ (70)

Then, we would obtain

$g_{2}(t)=a-d-\frac{2aT}{T_{\max}}+k_{1}I+\frac{dI/dt}{I}=g_{1}(t)+k_{1}I+\frac{\beta
V}{1+\alpha V}>g_{1}(t)$ (71)

Hence,

$\sup\\{g_{1}(t),g_{2}(t)\\}\leq
a-d-\frac{2aT}{T_{\max}}+k_{1}I+\frac{dI/dt}{I}\leq-\mu_{1}+\frac{dI/dt}{I}$
(72)

where

$\mu_{1}>0,\quad a-d-\frac{2aT}{T_{\max}}+k_{1}I\leq-\mu_{1}<0$ (73)

with the assumption that $k_{1}I$ is negligible compare to the term
$a-\frac{2aT}{T_{\max}}$. This assumption would be verified in the examples of
the simulation part below.

  * •

Case 2: If Point 4 of Theorem 7 holds, then

$-a+d+\frac{2aT}{T_{\max}}\leq\delta,$ (74)

which means that $G_{1}=\delta+c+k_{1}T$. Then, we obtain that

$\mu_{2}<0,\quad
g_{1}(t)<g_{2}(t)=k_{1}T-\delta+\frac{dI/dt}{I}\leq-\mu_{2}+\frac{dI/dt}{I}$
(75)

with the same assumption that $k_{1}T<\sigma$ in reasonably practical
scenarios. Hence,

$\sup\\{g_{1}(t),g_{2}(t)\\}\leq-\mu+\frac{dI/dt}{I}$ (76)

Let $\mu=\min\\{\mu_{1},\mu_{2}\\}$. Then, form (72) and (75), we have

$\sup\\{g_{1}(t),g_{2}(t)\\}\leq-\mu+\frac{dI/dt}{I},$ (77)

or

$D_{+}L(t)\leq\left(-\mu+\frac{dI/dt}{I}\right)L(t).$ (78)

According to Gronwall’s inequality, we would have

$\displaystyle L(t)$ $\displaystyle\leq
L(0)\exp\left(\int_{0}^{t}\left[-\mu+\frac{dI/dt}{I}\right]ds\right)$ (79)
$\displaystyle=L(0)\exp\left(\left[-\mu s+\ln(I(s))\right]_{0}^{t}\right)$
$\displaystyle=L(0)\exp(-\mu t)\exp\left(\ln(I(t))-\ln(I(0))\right)$
$\displaystyle=L(0)\exp(-\mu t)\frac{I(t)}{I(0)}$
$\displaystyle\leq\frac{ML(0)}{I(0)}\exp(-\mu t)\to 0\quad\text{as}\quad
t\to+\infty$

From (62), we conclude that

$\left(W_{1}(t),W_{2}(t),W_{3}(t)\right)\to 0\quad\text{as}\quad t\to+\infty$
(80)

This implies that the linear system equation (58) is asymptotically stable,
and therefore the periodic solution is asymptotically orbitally stable. This
proves proposition 1

∎

######  Theorem 8.

Suppose that

  1. 1.

$R_{0}>1$,

  2. 2.

$a_{1}>0,a_{3}+a_{5}>0,a_{1}(a_{2}+a_{4})-(a_{3}+a_{5})>0$.

Then, system (4) has an orbitally asymptotically stable periodic solution.

###### Proof.

First, we perform a change of variables as follows:

$z_{1}(t)=-T(t),\quad z_{2}(t)=I(t),\quad z_{3}(t)=-V(t)$ (81)

Applying this transformation to the system (4), we obtain

$\displaystyle\frac{dz_{1}(t)}{dt}$
$\displaystyle=-s-dz_{1}+az_{1}\left(1+\frac{z_{1}}{T_{\max}}\right)+\frac{\beta
z_{1}z_{3}}{1-\alpha z_{3}}+\rho z_{2}$ (82)
$\displaystyle\frac{dz_{2}(t)}{dt}$ $\displaystyle=\frac{\beta
z_{1}z_{3}}{1-\alpha z_{3}}-(\delta+\rho)z_{2}$
$\displaystyle\frac{dz_{3}(t)}{dt}$
$\displaystyle=-qz_{2}-cz_{3}+k_{1}z_{1}z_{3}$

The Jacobian matrix of the system (82) is then given by

$J(z)=\left(\begin{matrix}a-d+\frac{2az_{1}}{T_{\max}}+\frac{\beta
z_{3}}{1-\alpha z_{3}}&\rho&\frac{\beta z_{1}}{(1+\alpha z_{3})^{2}}\\\
\frac{\beta z_{3}}{1-\alpha z_{3}}&-(\delta+\rho)&\frac{\beta z_{1}}{(1+\alpha
z_{3})^{2}}\\\ k_{1}z_{3}&-q&-c+k_{1}z_{1}\end{matrix}\right)$ (83)

Similar to the definition of the set $D$ at 29, we define set $E$ as:

$E=\left\\{(z_{1},z_{2},z_{3}):z_{1}\leq 0,z_{2}\geq 0,z_{3}\leq 0\right\\}$
(84)

Since $J(z)$ has non-positive off diagonal elements at each point of $E$, (82)
is competitive at $E$. Set $z^{*}=(-T^{*},I^{*},V^{*})$. It is easy to see
that $z^{*}$ is unstable and $\det J(z^{*})<0$. Furthermore, it follows from
Theorem 5 that there exists a compact set $B$ in the interior of $E$ such that
for any $z_{0}\in\text{int}E$, there exists $T(z_{0})>0$ such that
$z(t,z_{0})\in B$ for all $t>T(z_{0})$. Consequently, by Theorem 1.2 in Zhu
and Smith (1994) [76] for the class of three-dimensional competitive systems,
it has an orbitally asymptotically stable periodic solution.

The proof is complete. ∎

##  3 The delay differential equation (DDE) model

In this section, we introduce a time delay into the model (5) to represent the
incubation time that the vectors need to become infectious. The model is
rewritten as follows

$\displaystyle\frac{dT}{dt}$
$\displaystyle=s-dT+aT\left(1-\frac{T}{T_{\max}}\right)-\frac{\beta
TV}{1+\alpha V}+\rho I$ (85) $\displaystyle\frac{dI}{dt}$
$\displaystyle=\frac{\beta T(t-\tau)V(t-\tau)}{1+\alpha
V(t-\tau)}-(\delta+\rho)I$ $\displaystyle\frac{dV}{dt}$ $\displaystyle=qI-cV-
k_{1}VT$

under the initial values

$T(\theta)=T_{0},\quad I(\theta)=I_{0},\quad
V(\theta)=V_{0}\quad\forall\theta\in[-\tau,0]$ (86)

All parameters of this delay model are the same as those of the system (4),
except that the additional positive constant $\tau$ represents the length of
the delay, in days.

This time delay parameter can be explained as follows: At time $t$, only
healthy cells that have been infected by the virus $\tau$ days ago (i.e at
time $t-\tau$ are infectious, provided that they have survived the incubation
period of $\tau$ days and were alive at the time $t-\tau$ when they infect the
healthy cells. As a result, the incidence term of healthy cells in the
derivative of infected cells with respect to time is modified from $\beta
T(t)V(t)$ to $\beta T(t-\tau)V(t-\tau)$.

The reproduction of this delay differential equation can be given the same as
the original ODE model, which is

$R_{01}=\frac{q\beta-(\delta+\rho)k_{1}}{c(\delta+\rho)}T_{0}$ (87)

Its biological meaning is that, if one virus is introduced in the population
of uninfected cells, the total number of secondary infected cells during the
infectious period would be $\frac{q\beta-(\delta+\rho)k_{1}}{c(\delta+\rho)}$.

###  3.1 Local and Global Stability of the Disease-free Equilibrium

Within this section, we will study the local and global stability of the
disease-free equilibrium $E_{0}$ of the delay model in two cases: when
$R_{0}>1$ and when $R_{0}<1$.

######  Theorem 9.

The disease-free equilibrium of the system (85) is locally asymptotically
stable if $R_{0}<1$, and is unstable if $R_{0}>1$.

###### Proof.

Linearizing the system (85) around $E_{0}=(T_{0},0,0)$, we obtain one negative
characteristic root

$\lambda_{1}=a-d-\frac{2aT_{0}}{T_{\max}}$ (88)

and the following transcendental characteristic equation whose roots are the
remaining eigenvalues

$\lambda^{2}+(\delta+\rho+c+k_{1}T_{0})+(c+k_{1}T_{0})(\delta+\rho)-q\beta
T_{0}e^{-\lambda\tau}=0$ (89)

For $\tau=0$, we obtain the exact same quadratic equation as the original ODE
system. In this case, we have proven previously that all eigenvalues of the
characteristic equation (89) have negative real parts. According to the Routh
- Hurwitz criterion, the disease free equilibirum $E_{0}$ will be locally
asymptotically stable when $R_{0}<1$ and is unstable when $R_{0}>1$.

As a result, we now only need to prove that the statement holds true for all
$\tau\neq 0$.

  * •

Case 1: $R_{0}>1$. In this case, we expect that (89) has one positive root and
the disease-free equilibrium is unstable. Indeed, we arrange the
characteristic equation into the form of

$\lambda^{2}+(\delta+\rho+c+k_{1}T_{0})\lambda=q\beta
T_{0}e^{-\lambda\tau}-(c+k_{1}T_{0})(\delta+\rho)$ (90)

Now, suppose that $\delta\in\mathbb{R}$ and denote

$\displaystyle F(\lambda)$
$\displaystyle=\lambda^{2}+(\delta+\rho+c+k_{1}T_{0})\lambda$ (91)
$\displaystyle G(\lambda)$ $\displaystyle=q\beta
T_{0}e^{-\lambda\tau}-(c+k_{1}T_{0})(\delta+\rho)$

We would then have that

$F(0)=0,\quad\lim_{\lambda\to+\infty}F(\lambda)=+\infty$ (92)

while

$\displaystyle G(0)=q\beta
T_{0}-(c+k_{1}T_{0})(\delta+\rho)=c(\delta+\rho)(R_{0}-1)>0,$ (93)
$\displaystyle G^{\prime}(\lambda)<0\quad\forall\lambda>0$

As a result, the two functions must intersect at a point $\lambda>0$, which
means that the equation (89) admits a positive real root, which means that the
disease-free equilibrium is unstable.

  * •

Case 2: $R_{0}<1$. First, we can notice that (89) can not have any non-
negative roots since , while

$\displaystyle F(0)$ $\displaystyle=0,\quad
F^{\prime}(\lambda)>0\quad\forall\lambda\geq 0$ (94) $\displaystyle
G(\lambda)$ $\displaystyle<0,\quad G^{\prime}(\lambda)<0\quad\forall\lambda>0$

As a result, if (89) has roots with non-negative real parts, they must be
complex and should be obtained from a pair of complex conjugate with cross the
imaginary axis. This means that (89) must have a pair of purely imaginary
roots for $\tau>0$.

As a result, we assume that $\lambda=i\omega$, and without loss of generality,
we assume that $\omega>0$ is a root of (89), meaning that

$-\omega^{2}+i\omega(\delta+\rho+c+k_{1}T_{0})+(c+k_{1}T_{0})(\delta+\rho)-q\beta
T_{0}(\cos(\omega\tau)+i\sin(\omega\tau))=0$ (95)

Separating the real and imaginary part, we would have

$\displaystyle-\omega^{2}+(c+k_{1}T_{0})(\delta+\rho)$ $\displaystyle=q\beta
T_{0}\cos(\omega\tau)$ (96) $\displaystyle(\delta+\rho+c+k_{1}T_{0})\omega$
$\displaystyle=-q\beta T_{0}\sin(\omega\tau)$

Squaring and adding up both sides of the two equations above, we obtain the
following fourth-order equation for $\omega$ as

$\omega^{4}+\omega^{2}\left[(\delta+\rho+c+k_{1}T_{0})^{2}-2(c+k_{1}T_{0})(\delta+\rho)\right]+[(c+k_{1}T_{0})(\delta+\rho)]^{2}-(q\beta
T_{0})^{2}=0$ (97)

To reduce this fourth-order equation into a quadratic equation, let
$z=\omega^{2}$ and denote the coefficients as

$\displaystyle a_{1}$
$\displaystyle=(\delta+\rho+c+k_{1}T_{0})^{2}-2(c+k_{1}T_{0})(\delta+\rho)$
(98) $\displaystyle a_{2}$
$\displaystyle=(c+k_{1}T_{0})^{2}(\delta+\rho)^{2}-(q\beta T_{0})^{2}$

The equation (97) can be rewritten as

$z^{2}+a_{1}z+a_{2}=0$ (99)

Since $R_{0}<1$, we would have that

$\displaystyle a_{1}$ $\displaystyle=(\delta+\rho)^{2}+(c+k_{1}T_{0})^{2}>0$
(100) $\displaystyle a_{2}$ $\displaystyle=[(c+k_{1}T_{0})(\delta+\rho)+q\beta
T_{0}]c(\delta+\rho)(1-R_{0})>0$

As a result, this means that the two roots of (99) have positive product,
which means that they have the same sign, regardless of being real or complex.
As these two roots also have negative real products, they would be either
negative real numbers, or complex conjugate with negative real parts. As a
result, the equation (99) can not have any positive real roots, leading to the
fact that there would be no $\omega$ such that $i\omega$ is a root of (89).
Using Rouche’s theorem, we conclude that the real parts of all eigenvalues of
the characteristic equation of the disease-free equilibrium (89) are all
negative for all delay values $\tau>0$.

In conclusion, if $R_{0}<1$, the disease-free equilibrium $E_{0}$ is locally
asymptotically stable.

The proof is complete. ∎

###  3.2 Local and Global Stability of the Positive Equilibrium

To study the stability of the steady states $\bar{E}$, we define

$x(t)=T(t)-\bar{T},\quad y(t)=I(t)-\bar{I},\quad z(t)=V(t)-\bar{V}.$ (101)

Then, the linearized system of (85) at $\bar{E}$ is given by

$\displaystyle\frac{dx(t)}{dt}$
$\displaystyle=\left[-d+a-\frac{2a\bar{T}}{T_{\max}}-\frac{\beta\bar{V}}{(1+\alpha\bar{V})^{2}}\right]x(t)+\rho
y(t)-\frac{\beta\bar{T}}{(1+\alpha\bar{V})^{2}}z(t)$ (102)
$\displaystyle\frac{dy(t)}{dt}$
$\displaystyle=\frac{\beta\bar{V}}{1+\alpha\bar{V}}x(t-\tau)-(\delta+\rho)y(t)+\frac{\beta\bar{V}}{(1+\alpha\bar{V})^{2}}z(t-\tau)$
$\displaystyle\frac{dz(t)}{dt}$
$\displaystyle=-k_{1}\bar{V}x(t)+qy(t)-(c+k_{1}\bar{T})z(t)$

The system (102) can be expressed in matrix form as follows

$\frac{d}{dt}\left(\begin{matrix}x(t)\\\ y(t)\\\
z(t)\end{matrix}\right)=A_{1}\left(\begin{matrix}x(t)\\\ y(t)\\\
z(t)\end{matrix}\right)+A_{2}\left(\begin{matrix}x(t-\tau)\\\ y(t-\tau)\\\
z(t-\tau)\end{matrix}\right),$ (103)

where $A_{1}$ and $A_{2}$ are $3\times 3$ matrices given by

$A_{1}=\left(\begin{matrix}a-d-\frac{2a\bar{T}}{T_{\max}}&\rho&-\frac{\beta\bar{T}}{(1+\alpha\bar{V})^{2}}\\\
0&-(\delta+\rho)&0\\\
-k_{1}\bar{V}&q&-(c+k_{1}\bar{T})\end{matrix}\right),\quad
A_{2}=\left(\begin{matrix}0&0&0\\\
\frac{\beta\bar{V}}{1+\alpha\bar{V}}&0&\frac{\beta\bar{T}}{(1+\alpha\bar{V})^{2}}\\\
0&0&0\end{matrix}\right).$ (104)

The characteristic equation of system (102) is given by

$\Delta(\lambda)=\left|\lambda I-A_{1}-e^{-\lambda\tau}A_{2}\right|=0,$ (105)

that is,

$\lambda^{3}+a_{1}\lambda^{2}+a_{2}\lambda+a_{5}=-e^{-\lambda\tau}(a_{3}+a_{4}\lambda),$
(106)

with $a_{i}(i=1,...,5)$ previously defined in (32).

Next, we shall study the distribution of the roots of the transcendental
equation (106) with respect to $0$ analytically. Based on the point (or
assumption) that the positive steady state of the original ODE model (4) is
stable, we will derive further conditions on the parameters to ensure that the
steady state of the delay model is still stable.

First, we will consider the base case when $\tau=0$. Then, the characeristic
equation (106) will become (89). Now, we will assume that all roots of his
equation, in this case, has all negative real parts, which is equivalent to
the fact that the conditions in Theorem 8 are satisfied. As the delay term
$\tau$ is considered to be continuous on $\mathbb{R}$, from Rouché’s Theorem
[16], the transcendental equaion (106) can only have roots with negative real
parts if and only if it has purely imaginary roots. We will investigate
whether (106) can admit any purely imaginary roots; from which, we will be
able to determine the conditions under which all eigenvalues would have
negative real parts.

We assume that $\lambda=\eta(\tau)+i\omega(\tau)\quad(\omega>0)$ is the
eigenvalue of the characteristic equation (106), where $\eta(\tau)$ and
$\omega(\tau)$ are functions depending on the delay term $\tau$. As the
positive equilibrium $\bar{E}$ of the model (4) is stable, we can say that
$\eta(0)<0$ at $\tau=0$.

If $\eta(\tau_{0})=0$ for some certain values of of $\tau_{0}>0$ (which means
that $\lambda=i\omega(\tau_{0})$ are purely imaginary roots of the
characteristic equation (106), the steady state $\bar{E}$ would lose is
stability and become unstable whenever $\eta(\tau_{0})$ is greater han $0$. In
other words, if there exists no $\omega(\tau_{0})$ such that the condition
above happens, or, if the characteristic equation (106) does not have any
purely imaginary roots for all values of $\tau$, the positive equilibrium
$\bar{E}$ is always stable. In the following part, we will prove that this
statement is indeed correct for equation (106).

Clearly $i\omega\ (\omega>0)$ is a root of equation (106) if and only if

$-i\omega^{3}-a_{1}\omega^{2}+ia_{2}\omega+a_{5}=-a_{3}(\cos(\omega\tau)-i\sin(\omega\tau))-a_{4}\omega(\sin(\omega\tau)+icos(\omega\tau))$
(107)

Separating the real and imaginary parts, we would have

$\displaystyle a_{1}\omega^{2}-a_{5}$
$\displaystyle=a_{3}\cos(\omega\tau)+a_{4}\omega\sin(\omega\tau)$ (108)
$\displaystyle\omega^{3}-a_{2}\omega$
$\displaystyle=-a_{3}\sin(\omega\tau)+a_{4}\omega\cos(\omega\tau)$

Squaring both sides of each equation above and adding up, we obtain the
following sixth-degree equation for $\omega$:

$\omega^{6}+(a_{1}^{2}-2a_{2})\omega^{4}+(a_{2}^{2}-2a_{1}a_{5}-a_{4}^{2})\omega^{2}+(a_{5}^{2}-a_{3}^{2})=0$
(109)

Since this equation contains only even powers of $\omega$, we can reduce the
order by letting once again $z=\omega^{2}$ and

$\displaystyle m_{1}$ $\displaystyle=a_{1}^{2}-2a_{2},$ (110) $\displaystyle
m_{2}$ $\displaystyle=a_{2}^{2}-2a_{1}a_{5}-a_{4}^{2},$ $\displaystyle m_{3}$
$\displaystyle=a_{5}^{2}-a_{3}^{2},$

the equation (109) becomes

$h(z)=z^{3}+m_{1}z^{2}+m_{2}z+m_{3}=0$ (111)

In order to show that the positive equilibrium $\bar{E}$ is locally stable, we
have to prove that the equation (111) does not have any positive real root
which associates to the square of $\omega$; that is, (106) can not have any
purely imaginary roots. The Theorem below provides us with necessary
conditions satisfying the result.

######  Theorem 10.

If $m_{3}\geq 0$ and $m_{2}>0$, the equation (109) has no positive real roots.

###### Proof.

We will proceed to prove the lemma above using contradiction.

Assume that there exists at least one positive real roots for the equation
$h(z)=0$.

Notice that $h(0)=m_{3}\geq 0$. This means that in order for the equation
(109) to have a positive real roots, there exists $z_{0}\geq 0$ such that

$\frac{dh(z_{0})}{dz_{0}}\leq 0.$ (112)

This is equivalent to

$\frac{dh(z_{0})}{dz_{0}}=3z_{0}^{2}+2m_{1}z_{0}+m_{2}\leq 0$ (113)

or

$\frac{-m_{1}-\sqrt{m_{1}^{2}-3m_{2}}}{2}\leq
z\leq\frac{-m_{1}+\sqrt{m_{1}^{2}-3m_{2}}}{2}<0$ (114)

This contradicts our original assumption that $z_{0}\geq 0$, which means that
there does not exist any $z_{0}\geq 0$ such that $\frac{dh(z_{0})}{dz_{0}}>0$,
or the equation does not have any positive real roots.

The proof is complete. ∎

This theorem has implied that there exists no $\omega$ such that $i\omega$ is
an eigenvalue of the characteristic function (106). As a result, from Rouche’s
theorem [16], the real parts of the eigenvalues of (106) are negative for all
$\tau\geq 0$. Summarizing all the above analysis, we have the following
theorem

######  Theorem 11.

Suppose that

  1. 1.

$a_{1}>0,a_{3}+a_{5}>0,a_{1}(a_{2}+a_{4})-(a_{3}+a_{5})>0$;

  2. 2.

$m_{3}\geq 0$ and $m_{2}>0$.

Then, the infected steady state $\bar{E}$ of the delay model (85) is
absolutely stable; that is, $\bar{E}$ is asymptotically stable for all
$\tau\geq 0$.

######  Remark 2.

The Theorem 11 indicates that if the parameters satisfy both of the
conditions, the equilibrium $\bar{E}$ of (85) is asymptotically stable
regardless of the value of the delay (independent asymptotic stability).
However, we also need to note that if any of the conditions in Theorem 11 is
violated (particularly the inequalities in Point 2), the stability of the
equilibrium will then depend on the delay value; and when the delay value
varies, the equilibrium can lose stability, leading to oscillations

For example, if

  1. 1.

If $m_{3}<0$: From equation (111), we would have that

$h(0)<0,\quad\lim_{z\to+\infty}=+\infty,$ (115)

which means that (111) has at least one positive real root, denoted by
$\omega_{0}$.

  2. 2.

If $m_{2}<0$, we would have that

$\frac{-m_{1}+\sqrt{m_{1}^{2}-3m_{2}}}{2}=\frac{-3m_{2}}{m_{1}+\sqrt{m_{1}^{2}+3m_{2}}}>0$
(116)

which means that the equation (111) has one positive real root $\omega_{0}$.

These two cases implies that the characteristic equation (106) has a pair of
purely imaginary roots $\pm i\omega_{0}$.

Next, we would focus on the bifurcation analysis, using the delay term $\tau$
as the bifurcation parameter, in light that the solutions of (106) as function
of this parameter.

Let $\lambda(\tau)=\mu(\tau)+i\omega(\tau)$ be the eigenvalue of (108) such
that for some initial values of the bifurcation parameter $\tau_{0}$, we would
have $\mu(\tau_{0})=0,\quad\omega(\tau_{0})=\omega_{0}$. From the system
(108), we would have:

$\tau_{j}=\frac{1}{\omega_{0}}\arccos\left(\frac{a_{4}\omega_{0}^{4}+(a_{1}a_{3}-a_{2}a_{4})\omega_{0}^{2}-a_{3}a_{5}}{a_{3}^{2}+a_{4}^{2}\omega_{0}^{2}}\right)+\frac{2j\pi}{\omega_{0}}$
(117)

Moreover, we can verify the following transversal condition:

$\frac{d}{d\tau}\Re(\lambda(\tau))\|_{\tau=\tau_{0}}=\frac{d}{d\tau}\mu(\tau)\|_{\tau=\tau_{0}}>0$
(118)

holds. By continuity, the real part of $\lambda(\tau)$ becomes positive when
$\tau>\tau_{0}$ and the steady state becomes unstable. Moreover, a Hopf
bifurcation occurs when $\tau$ passes through the critical value $\tau_{0}$
(see [21]).

To apply the Hopf bifurcation theorem stated in Marsden and McCracken [41], we
state the following theorem

######  Theorem 12.

Suppose that $\omega_{0}$ is the largest positive simple root of (109). Then,
$i\omega(\tau_{0})=i\omega_{0}$ is a simple root of (109), and
$\eta(\tau)+i\omega(\tau)$ is differentiable with respect to $\tau$ in a
neighborhood of $\tau=\tau_{0}$

After previous reasoning, we admit that $i\omega_{0}$ is a simple root of
(109), which is an analytic equation; as a result, using the analytic version
of the implicit function theorem mentioned in Chow and Hale (1982) [12], we
have that $\eta(\tau)+i\omega(\tau)$ is well-defined and analytic in
neighborhood of $\tau=\tau_{0}$

To establish the Hopf bifurcation at $\tau=\tau_{0}$, we need to show that

$\frac{d\Re(\lambda(\tau))}{d\tau}\|_{\tau=\tau_{0}}>0.$ (119)

In order to prove this inequality, we first start with a lemma and its
respective proof.

######  Lemma 1.

Suppose that $z_{1},z_{2},z_{3}$ are the roots of
$h(z)=z^{3}+m_{1}z^{2}+m_{2}z+m_{3}=0\quad(m_{2}<0)$, and
$z_{3}\in\mathbb{R}^{+}$ is the largest positive simple root, then

$\frac{dh(z)}{dz}\|_{z=z_{3}}>0.$ (120)

###### Proof.

We will proceed to prove the lemma above with contradiction.

Assume that the largest positive simple root $z_{3}$ of the equation $h(z)=0$
and

$\frac{dh(z)}{dz}\|_{z=z_{3}}\leq 0.$ (121)

As a result, there exists a $\bar{z}\in\mathbb{R}^{+},\bar{z}>z_{3}$ such that
$h(\bar{z})<0$.

According to the Intermediate Value Theorem, for any
$h\in\mathbb{R},h\in[h(\bar{z},+\infty)$, there always exists
$z_{4}\in[\bar{z},+\infty)$. Taking $h=0$, we would have

$\begin{cases}h(z_{4})=0\\\ z_{4}>\bar{z}>z_{3}\end{cases},$ (122)

which means that $z_{4}$ is the highest positive simple root of $h(z)=0$, not
$z_{3}$. This contradicts our original assumption.

In conclusion, if $z_{3}$ is the largest positive simple root,

$\frac{dh(z)}{dz}\|_{z=z_{3}}>0.$ (123)

The proof is complete. ∎

From the equation (106), derivating both sides with respect to $\tau$, we
obtain

$\displaystyle\left(3\lambda^{2}+2a_{1}\lambda+a_{2}\right)\frac{d\lambda}{d\tau}=$
$\displaystyle\left[-\tau\exp(-\lambda\tau)(-a_{3}-a_{4}\lambda)+\exp(-\lambda\tau)(-a_{4})\right]\frac{d\lambda}{d\tau}$
(124) $\displaystyle-\lambda\exp(-\lambda\tau)(-a_{3}-a_{4}\lambda).$

This gives us

$\displaystyle\left(\frac{d\lambda}{d\tau}\right)^{-1}$
$\displaystyle=\frac{3\lambda^{2}+2a_{1}\lambda+a_{2}+\tau\exp(-\lambda\tau)(-a_{3}-a_{4}\lambda)-\exp(-\lambda\tau)(-a_{4})}{-\lambda\exp(-\lambda\tau)(-a_{3}-a_{4}\lambda)}$
(125)
$\displaystyle=\frac{3\lambda^{2}+2a_{1}\lambda+a_{2}}{-\lambda\exp(-\lambda\tau)(-a_{3}-a_{4}\lambda)}+\frac{a_{4}}{\lambda(a_{3}+a_{4}\lambda)}-\frac{\tau}{\lambda}$
$\displaystyle=\frac{2\lambda^{3}+a_{1}\lambda^{2}-a_{5}}{-\lambda^{2}(\lambda^{3}+a_{1}\lambda^{2}+a_{2}\lambda+a_{5})}+\frac{-a_{3}}{\lambda^{2}(a_{3}+a_{4}\lambda)}-\frac{\tau}{\lambda}$

Thus,

$\displaystyle\text{Sign}\left\\{\frac{d(\Re(\lambda))}{d\tau}\right\\}$
$\displaystyle=\text{Sign}\left\\{\Re\left(\frac{d\lambda}{d\tau}\right)^{-1}\right\\}$
(126)
$\displaystyle=\text{Sign}\left\\{\Re\left[\frac{2\lambda^{3}+a_{1}\lambda^{2}-a_{5}}{-\lambda^{2}(\lambda^{3}+a_{1}\lambda^{2}+a_{2}\lambda+a_{5})}\right]_{\lambda=i\omega_{0}}+\Re\left[\frac{-a_{3}}{\lambda^{2}(a_{3}+a_{4}\lambda)}\right]_{\lambda=i\omega_{0}}\right\\}$
$\displaystyle=\text{Sign}\left\\{\Re\left[\frac{-2\omega_{0}^{3}i-a_{1}\omega_{0}^{2}-a_{5}}{\omega_{0}^{2}(-\omega_{0}^{3}i-a_{1}\omega_{0}^{2}+a_{2}\omega_{0}i+a_{5})}\right]+\Re\left[\frac{-a_{3}}{-\omega_{0}^{2}(a_{3}+a_{4}\omega_{0}i)}\right]\right\\}$
$\displaystyle=\text{Sign}\left\\{\frac{2\omega_{0}^{6}+(a_{1}^{2}-2a_{2})\omega_{0}^{4}-a_{5}^{2}}{\omega_{0}^{2}[(a_{2}\omega_{0}-\omega_{0}^{3})^{2}+(a_{5}-a_{1}\omega_{0}^{2})^{2}]}+\frac{a_{3}^{2}}{\omega_{0}^{2}(a_{4}^{2}\omega_{0}^{2}+a_{3}^{2})}\right\\}$
$\displaystyle=\text{Sign}\left\\{\frac{3\omega_{0}^{4}+2(a_{1}^{2}-2a_{2})\omega_{0}^{2}+(a_{2}^{2}-2a_{1}a_{5}-a_{4}^{2})}{(a_{2}\omega_{0}-\omega_{0}^{3})^{2}+(a_{5}-a_{1}\omega_{0}^{2})^{2}}\right\\}$

Since

$h(z)=z^{3}+m_{1}z^{2}+m_{2}z+m_{3},$ (127)

we would have

$\frac{dh(z)}{dz}=3z^{2}+2m_{1}z+m_{2}=3z^{2}+2(a_{1}^{2}-2a_{2})z+(a_{2}^{2}-2a_{1}a_{5}-a_{4}^{2})$
(128)

As we have assumed that $\omega_{0}$ is the largest positive simple root of
the equation (109), from Lemma 1, we get

$\frac{dh(z)}{dz}\|_{z=\omega_{0}^{2}}>0.$ (129)

Hence,

$\frac{3\omega_{0}^{4}+2(a_{1}^{2}-2a_{2})\omega_{0}^{2}+(a_{2}^{2}-2a_{1}a_{5}-a_{4}^{2})}{(a_{2}\omega_{0}-\omega_{0}^{3})^{2}+(a_{5}-a_{1}\omega_{0}^{2})^{2}}=\frac{\frac{dh(z)}{dz}\|_{z=\omega_{0}^{2}}}{(a_{2}\omega_{0}-\omega_{0}^{3})^{2}+(a_{5}-a_{1}\omega_{0}^{2})^{2}}>0$
(130)

or

$\text{Sign}\left\\{\frac{d(\Re(\lambda))}{d\tau}\right\\}=\text{Sign}\left\\{\frac{3\omega_{0}^{4}+2(a_{1}^{2}-2a_{2})\omega_{0}^{2}+(a_{2}^{2}-2a_{1}a_{5}-a_{4}^{2})}{(a_{2}\omega_{0}-\omega_{0}^{3})^{2}+(a_{5}-a_{1}\omega_{0}^{2})^{2}}\right\\}=1$
(131)

i.e

$\frac{d(\Re(\lambda))}{d\tau}>0$ (132)

The Hopf bifurcation analysis above can be summarized in the following
theorem.

######  Theorem 13.

Suppose that

$a_{1}>0,\quad a_{3}+a_{5}>0,\quad a_{1}(a_{2}+a_{4})-(a_{3}+a_{5})>0$ (133)

and

$R_{0}>1$ (134)

If

$m_{3}<0\quad\lor\quad m_{3}\geq 0,\quad m_{2}<0,$ (135)

the infected steady state $\bar{E}$ of the delay model (85) is asymptotically
stable when $\tau<\tau_{0}$ and unstable when $\tau>\tau_{0}$, where

$\tau_{0}=\frac{1}{\omega_{0}}\arccos\left(\frac{a_{4}\omega_{0}^{4}+(a_{1}a_{3}-a_{2}a_{4})\omega_{0}^{2}-a_{3}a_{5}}{a_{3}^{2}+a_{4}^{2}\omega_{0}^{2}}\right)$
(136)

When $\tau=\tau_{0}$, a Hopf bifurcation occurs; that is, a family of periodic
solutions bifurcates from $\bar{E}$ as $\tau$ passes through the critical
value $\tau_{0}$.

##  4 Numerical simulation

After providing all the analytical tools and qualitatively analysing the
system for patterns on its dynamics, in this section ,we will perform some
numerical analysis on the model to verify the previous results.

###  4.1 Simulation tools

The numerical simulation is conducted on the programming language Julia
through the package DifferentialEquation.jl, A Performant and Feature-Rich
Ecosystem for Solving Differential Equations in Julia by Rackauckas and Nie
(2017) [61].

In order to avoid any stiffness in the ODE/DDE models, the algorithm for the
Method of Steps in Julia is set to Rosenbrock23, which is the same as the
classic ODE solver ode23s in MATLAB.

For the complete version of the Julia notebooks for simulation, please refer
to the Github repository at https://github.com/hoanganhngo610/DDE-HIV-
NGO2020etal.

###  4.2 Simulation results

####  4.2.1 Simulation of the ODE model

Parameters and Variables |  | Values
---|---|---
Dependent variables |  |
$T$ | Uninfected CD4$+$ T-cell population size |  $250$ mm$-3$
$I$ | Infected CD4$+$ T-cell density |  $50$ mm$-3$
$V$ | Initial density of HIV RNA |  $160$ mm$-3$
Parameters and Constants |  |
$s$ | Source term for uninfected CD4$+$ T-cells |  $5$ day$-1$ mm$-3$
$d$ | Natural death rate of CD4$+$ T-cells | 0.01 day-1
$a$ | Growth rate of CD4$+$ T-cell population |  $0.8$ day$-1$
$T_{\max}$ | Maximal population level of CD4$+$ T-cells |  $1500$ mm$-3$
$\beta$ | Rate CD4$+$ T-cells became infected with virus |  $2.4\times 10^{-4}$ mm$-3$
$\alpha$ | Saturated mass-action term | $0.001$
$\rho$ | Rate of cure |  $0.01$ day$-1$
$\delta$ | Blanket death rate of infected CD4$+$ T-cells |  $0.3$ day$-1$
$q$ | Reproduction rate of the infected CD4$+$ T-cells |  $500$ mm$-3$ day$-1$
$c$ | Death rate of free virus |  $8$ day$-1$
Table 1:  Preliminary values of variables and parameters for viral spread. Parameters | Original scenario | Scenario #2 | Scenario #3 | Scenario #4
---|---|---|---|---
$s$ | $5$ | $-$ | $-$ | $-$
$d$ | $0.01$ | $-$ | $-$ | $-$
$a$ | $0.8$ | $8$ | $-$ | $-$
$T_{\max}$ | $1500$ | $-$ | $-$ | $-$
$\beta$ | $2.4\times 10^{-4}$ | $-$ | $0.0024$ | $0.0024$
$\alpha$ | $0.001$ | $0.0001$ | $0.000001$ | $0.000001$
$\rho$ | $0.01$ | $0.01$ | $-$ | $-$
$\delta$ | $0.3$ | $5$ | $-$ | $-$
$q$ | $500$ | $-$ | $2.5$ | $2.5$
$c$ | $8$ | $1.3$ | $3$ | $1.3$
Table 2:  Values of parameters for viral spread in different scenarios.

Within the range of parameters that are proven to be realistic in medical
research, we investigate the behavior of the model within 4 different
scenarios.

  * •

The original scenario: In this scenario, the condition $1$, $2$ and $3$ in
Theorem 7 are satisfied. This means that, the positive equilibrium of the
system (4) is globally asymptotically stable.

Figure 1: The ODE model is locally asymptotically stable with parameters in
the original scenario

  * •

Scenario #2: In this scenario, the conditions $1$, $2$ and $4$ in Theorem 7
are satisfied. This means that, the positive equilibrium of the system (4) is
also globally asymptotically stable.

Figure 2: The ODE model is locally asymptotically stable with parameters in
Scenario #2

  * •

Scenario #3: In this scenario, the conditions $1$ and $2$ of Theorem 3 is
satisfied. This means that, the positive equilibrium of the system (4) is
locally asymptotically stable.

Figure 3: The ODE model is locally asymptotically stable with parameters in
Scenario #3

  * •

Scenario #4: In this scenario, the conditions $1$ and $2$ of Theorem 8 is
satisfied. This means that, the positive equilibrium of the system (4) is
orbitally asymptotically stable.

Figure 4: The ODE model is orbitally asymptotically stable with parameters in
Scenario #4

####  4.2.2 Simulation of the DDE model

Parameters (DDE) | Original scenario (ODE) | Scenario #1 | Scenario #2 | Scenario #3 | Scenario #4
---|---|---|---|---|---
$s$ | $5$ | $-$ | $-$ | $-$ | $-$
$d$ | $0.01$ | $-$ | $-$ | $-$ | $-$
$a$ | $0.8$ | $-$ | $5$ | $5$ | $5$
$T_{\max}$ | $1500$ | $-$ | $-$ | $-$ | $-$
$\beta$ | $2.4\times 10^{-4}$ | $-$ | $-$ | $-$ | $-$
$\alpha$ | $0.001$ | $0.000001$ | $0.000001$ | $0.000001$ | $0.000001$
$\rho$ | $0.01$ | $-$ | $-$ | $-$ | $0.3$
$\delta$ | $0.3$ | $-$ | $-$ | $-$ | $-$
$q$ | $500$ | $-$ | $-$ | $-$ | $-$
$c$ | $8$ | $-$ | $-$ | $-$ | $-$
$\tau$ | N/A | $0.4$ | $10$ | $5$ | $5$
Table 3:  Values of parameters for viral spread in different scenarios.

First of all, instead of keeping $\alpha=0.001$, we modify this parameter into
$\alpha=0.000001$ so that we can observe different behaviors while modifying
other parameters.

  * •

Scenario #1 (DDE): When $\tau=0.4$ and all other variables are kept the same
as the original scenario in the ODE setting (apart from $\alpha$), $T(t)$,
$I(t)$ and $V(t)$ all converges to their positive equilibrium. We say that, in
this setting, the positive equilibrium $\bar{E}$ is globally asymptotically
stable.

Figure 5: The DDE model is globally asymptotically stable with parameters in
Scenario #1, with the delay term $\tau=0.4$.

  * •

Scenario #2 (DDE): When we modify $a=5$ and $\tau=10$, the parameters, in this
setting, satisfy the conditions of Theorem 13. This means that, the positive
equilibrium $\bar{E}$ is orbitally asymptotically stable, or in other words,
there exists a positive periodic solution for all the components of the
system.

Figure 6: The DDE model is globally asymptotically stable with parameters in
Scenario #2, with $a=5$ and the delay term $\tau=10$.

  * •

From the last two scenarios, Scenario #3 (DDE) and Scenario #4 (DDE), we can
draw a conclusion that the solution of the system would return to stability
when the cure rate $\rho$ is increased. For example, if we select $\rho=0.3$
instead of $\rho=0.01$ with all other parameters kept identical, the system
would admit a lower global asymptotic stability. We can conclude that $\rho$
is an important parameter in the sense that increasing it helps us control the
disease.

Figure 7: The DDE model is globally asymptotically stable with parameters in
Scenario #3 and Scenario #4, with $\rho=0.01$ and $\rho=0.3$, respectively.
These graphs show that the cure rate is an important parameter in controlling
the disease.

##  Appendix A List of macros for formatting text, figures and tables

######  Theorem 14 (Gronwall, 1919).

Let $I$ denote an interval of the real line of the form $[a,\inf)$ or
$[a,b]$lr $[a,b)$ with $a<b$ Let $\beta$ and $u$ be real-valued continuous
functions defined on $I$. If $u$ is a differentiable in the interior $I^{0}$
of $I$ (the interval $I$ without the end points $a$ and possibly $b$) and
satisfies the differential inequality

$u^{\prime}(t)\leq\beta(t)u(t),t\in I^{0}$ (137)

then $u$ is bounded by the solution of the corresponding differential equation
$\nu^{\prime}(t)=\beta(t)\nu(t)$:

$u(t)\leq u(a)\exp\left(\int_{a}^{t}\beta(s)ds\right)$ (138)

######  Theorem 15 (Lyapunov’s stability).

Let a function $V(\mathbf{X})$ be continuously differentiable in a
neighbourhood $U$ of the origin. The function $V(\mathbf{X})$ is called the
Lyapunov function for an autonomous system

$\mathbf{X}^{\prime}=f(\mathbf{X})$ (139)

if the following conditions are met:

  1. 1.

$V(\mathbf{X})>0$ for all $\mathbf{X}\in U\setminus\\{0\\}$;

  2. 2.

V(0) = 0;

  3. 3.

$\frac{dV}{dt}\leq 0$ for all $\mathbf{X}\in U$.

Then, if in a neighborhood $U$ of the zero solution $\mathbf{X}=0$ of an
autonomous system there is a Lyapunov function $V(\mathbf{X})$ with a negative
definite derivative $\frac{dV}{dt}$ for all $\mathbf{X}\in U\setminus\\{0\\}$,
then the equilibrium point $\mathbf{X}=0$ of the system is asymptotically
stable.

######  Theorem 16 (Perron - Frobenius theorem).

[17] Let $A$ be a irreducible Metzler matrix (A Metzler matrix is a matrix
whose all of its off-diagonal elements are non-negative). Then, $\lambda_{M}$,
the eigenvalue of $A$ of largest real part is real, and the elements of its
associated eigenvector $v_{M}$ are positive. Moreover, any eigenvector of $A$
with non-negative elements belongs the the span of $v_{M}$.

######  Theorem 17 (Implicit Function Theorem (Chow and Hale, 1982)).

Suppose that

  * •

$X,Y,Z$ are Banach spaces,

  * •

$F:U\times V\to Z$ is continuously differentiable,

  * •

$F(x_{0},y_{0})=0$ and $D_{x}F(x_{0},y_{0})$ has a bounded inverse.

Then, there exists a neighborhood $U_{1}\times V_{1}\in U\times V$ of
$(x_{0},y_{0})$ and a function $f:V_{1}\to U_{1},f(y_{0})=x_{0}$ such that

$F(x,y)=0\text{ for }(x,y)\in U_{1}\times V_{1}\text{ iff }x=f(y)$ (140)

If $F\in C^{k}(U\times V,Z),k\geq 1$ or analytic in a neighborhood
$(x_{0},y_{0})$, then $f\in C^{k}(V_{1},X)$ or is analytic in a neighborhood
of $y_{0}$.

######  Theorem 18 (Poincaré - Bendixson theorem).

[73]

Given a differentiable real dynamical system defined on an open subset of the
plane, every non-empty compact $\omega$-limit set of an orbit, which contains
only finitely many fixed points, is either

  * •

a fixed point,

  * •

a periodic orbit, or

  * •

a connected set composed of a finite number of fixed points together with
homoclinic and heteroclinic orbits connecting these.

Moreover, there is at most one orbit connecting different fixed points in the
same direction. However, there could be countably many homoclinic orbits
connecting one fixed point.

Next, we will give the definition of an additive compound matrix and consider
the particular case when it’s a square matrix [66]. A survey of properties of
additive compound matrices, along with their connections to differential
equations have been investigated in [48, 38]

We will start with the definition of the $k$-th exterior power (or
multiplicative compound) of an $n\times m$ matrix.

######  Definition 1 (Multiplicative compound of a matrix).

Let $A$ be an $n\times m$ matrix of real or complex numbers. Let
$a_{i_{1},i_{2},...,i_{k},j_{1},j_{2},...,j_{k}}$ be the minor of $A$
determined by the rows $(i_{1},...,i_{k})$ and the columns
$(j_{1},...,j_{k})$, $1\leq i_{1}<i_{2}<...<i_{k}\leq n,1\leq
j_{1}<j_{2}<...<j_{k}\leq m$. The $k$-th multiplicative compound matrix
$A^{(k)}$ of $A$ is the $\binom{n}{k}\times\binom{m}{k}$ matrix whose entries,
written in lexicographic order, are $a_{i_{1},...,i_{k},j_{1},...,j_{k}}$.

In particular, when $A$ is an $n\times k$ matrix with columns
$a_{1},a_{2},...,a_{k}$, $A^{(k)}$ is the exterior product $a_{1}\lor
a_{2}\lor...\lor a_{k}$.

In the case $m=n$, the additive compound matrices are defined as follows.

######  Definition 2.

Let $A$ be an $n\times n$ matrix. The $k$-th additive compound $A^{[k]}$ of
$A$ is the $\binom{n}{k}\times\binom{n}{k}$ matrix given by

$A{[k]}=D(I+hA)\|_{h=0}$ (141)

If $B=A^{[k]}$, the following formula for $b_{i,j}$ can be deduced from the
equation (141), For any integer $i=1,...,\binom{n}{k}$, let
$(i)=(i_{1},i_{2},...,i_{k})$ be the $i$-th member in the lexicographic
ordering of all $k$-tuples of integers such that $1\leq
i_{1}<i_{2}<...<i_{k}\leq n$. Then,

$b_{i,j}=\begin{cases}\begin{aligned}
a_{i_{1},i_{1}}+...+a_{i_{k},i_{k}}\quad&\text{ if }(i)=(j)\\\
(-1)^{r+s}a_{i_{s},j_{r}}\quad&\text{ if exactly one entry $i_{s}$ in $(i)$
does not occur in $(j)$ and $j_{r}$ does not occur in $(i)$, }\\\
0\quad&\text{ if $(i)$ differs from $(j)$ in two or more
entries.}\end{aligned}\end{cases}$ (142)

In the extreme cases when $k=1$ and $k=n$, we would have that $A^{[1]}=A$ and
$A^{[n]}=\text{tr}(A)$. For $n=3$ ,which is the case that we are considering
in this paper, we would have the matrices $A^{[k]},k=0,1,2$ as follows:

$A^{[1]}=A,\quad A^{[2]}=\left(\begin{matrix}a_{11}+a_{22}&a_{23}&-a_{13}\\\
a_{32}&a_{11}+a_{33}&a_{12}\\\
-a_{31}&a_{21}&a_{22}+a_{33},\end{matrix}\right),\quad
A^{[3]}=a_{11}+a_{22}+a_{33}$ (143)

## Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the
publication of this paper.

## Acknowledgement

The authors would like to thank Nguyen Tran Hai Yen, undergraduate student at
the Faculty of Biology and Biotechnology, Ho Chi Minh University of Science,
VNU - HCM, Class of 2022 for providing valuable biological insights and ideas
to support this research.

## References

  * [1] A.T. al.  “Quantitative image analysis of HIV-1 infection in lymphoid tissue”  In _Science_ 274, 1996, pp. 985
  * [2] R.. al.  “Constant mean viral copy number per infected cell in tissues regardless of high, low, or undetectable plasma HIV RNA”  In _Journal of Experimental Medicine_ 189, 1999, pp. 1545
  * [3] W. al.  “Kinetics of response in lymphoid tissues to antiretroviral therapy of HIV-1 infection”  In _Science_ 276, 1997, pp. 960
  * [4] R.. Anderson  “Mathematical and statistical studies of the epidemiology of HIV”  In _AIDS_ 4, 1990, pp. 107
  * [5] R.. Anderson and R.. May  “Complex dynamical behavior in the interaction between HIV and the immune system”  In _Cell to Cell Signalling: From Experiments to Theoretical Models_ New York: Academic Press, 1989, pp. 335
  * [6] J.. Bailey, J.. Fletcher, E.. Chuck and R.. Shrager  “A kinetic model of CD4+ lymphocytes with the human immunodeficiency virus (HIV)”  In _BioSystems_ 26, 1992, pp. 177
  * [7] R. Bellman and K.. Cooke  “Differential-Difference Equations”  New York: Academic Press, 1993
  * [8] R.. Boer and A.. Perelson  “Target Cell Limited and Immune Control Models of HIV Infection: A Comparison”  In _Journal of Theoretical Biology_ 190 Elsevier, 1998, pp. 201–214  DOI: 10.1006/jtbi.1997.0548
  * [9] S. Bonhoeffer, R.. May, G.. Shaw and M.. Nowak  “Virus dynamics and drug therapy”  In _Proceedings of the National Academy of Sciences of the United States of America_ 94, 1997, pp. 6971
  * [10] S. Busenberg and K. Cooke  “Vertically Transmitted Diseases”  Berlin: Springer, 1993
  * [11] G. Butler, H.. Freedman and P. Waltma  “Uniform persistence system”  In _Proceedings of the American Mathematical Society_ 96, 1986, pp. 425–430
  * [12] S.. Chow and J.. Hale  “Methods of Bifurcation Theory”  In _Grundlehren der Mathematischen Wissenschaften_ 251 New York, NY, USA: Springer, 1982
  * [13] R.. Culshaw  “Mathematical Models of Cell-to-Cell and Cell-Free Viral Spread of HIV Infection”, 1997
  * [14] Rebecca V. Culshaw and Shigui Ruan  “A delay - differential equation model of HIV infection of CD4+ T-cells”  In _Mathematical Biosciences_ 165 Elsevier, 2000, pp. 27–39
  * [15] J.. Cushing  “Integrodifferential Equations and Delay Models in Population Dynamics”  Heidelberg: Springer, 1977
  * [16] J. Dieudonné  “Foundations of Modern Analysis”  New York: Academic Press, 1960
  * [17] F.. Gantmacher  “The Theory of Matrices”  New York: Chelsea Publishing Company, 1959
  * [18] K. Gopalsamy  “Stability and Oscillations in Delay-Differential Equations of Population Dynamics”  Dordrecht: Kluwer, 1992
  * [19] J.. Hale and P. Waltman  “Persistence in infinite-dimensional systems”  In _SIAM Journal on Mathematical Analysis_ 20, 1989, pp. 388–396
  * [20] R. Harvey  “Microbiology”  Philadelphia, PA: Lippincott Williams & Wilkins, 2012, pp. 295–306
  * [21] B.. Hassard, N.. Kazarinoff and Y.. Wan  “Theory and Applications of Hopf Bifurcation”  Cambridge: Cambridge University, 1981
  * [22] B.. Hassard, N.. Kazarinoff and Y.. Wan  “Theory and Applications of Hopf Bifurcation”  In _London Mathematical Society Lecture Note Series_ 41 Cambridge, UK: Cambridge University Press
  * [23] A… Herz et al.  “Viral dynamics in vivo: limitations on estimates of intracellular delay and virus decay”  In _Proceedings of the National Academy of Sciences of the United States of America_ 93, 1996, pp. 7247
  * [24] M.. Hirsch  “System of differential equations which are competitive or cooperative, IV”  In _SIAM Journal on Mathematical Analysis_ 21, 1990, pp. 1225–1234
  * [25] D. Ho et al.  “Rapid turnover of plasma virions and CD4+ lymphocytes in HIV-1 infection”  In _Nature_ 373, 1995, pp. 123–126
  * [26] T. Hraba, J. Dolezal and S. Celikovsky  “Model-based analysis of CD4+ lymphocyte dynamics in HIV infected individuals”  In _Immunobiology_ 181, 1990, pp. 108
  * [27] N. Intrator, G.. Deocampo and L.. Cooper  “Analysis of immune system retrovirus equations”  In _Theoretical Immunology II_ , 1988, pp. 85
  * [28] T.. Kepler and A.S. Perelson  “Cyclic re-entry of germinal center B cells and the efficiency of affinity maturation”  In _Immunology Today_ 14, 1993, pp. 412–415
  * [29] D.. Kirschner  “Using mathematics to understand HIV immune dynamics”  In _Notices Of The American Mathematical Society_ 43, 1996, pp. 191
  * [30] D.. Kirschner, S. Lenhart and S. Serbin  “Optimal control of the chemotherapy of HIV”  In _Journal of Mathematical Biology_ 35, 1997, pp. 775
  * [31] D.. Kirschner and A.. Perelson  “A model for the immune system response to HIV: AZT treatment studies”  In _Mathematical Population Dynamics: Analysis of Heterogeneity, vol. 1, Theory of Epidemics_ , 1995, pp. 295
  * [32] D.. Kirschner and G.. Webb  “A model for the treatment strategy in the chemotherapy of AIDS”  In _Bulletin of Mathematical Biology_ 58, 1996, pp. 367
  * [33] D.. Kirschner and G.. Webb  “Understanding drug resistance for monotherapy treatment of HIV infection”  In _Bulletin of Mathematical Biology_ 59, 1997, pp. 763
  * [34] Y. Kuang  “Delay-Differential Equations with Applications in Population Dynamics”  New York: Academic Press, 1993
  * [35] J.. LaSalle  “The Stability of Dynamical Systems”  In _SIAM_ , 1976
  * [36] Dan Li and Wanbiao Ma  “Asymptotic properties of a HIV-1 infection model with time delay”  In _Journal of Mathematical Analysis and Applications_ 335 (1), 2007, pp. 683–691
  * [37] M.. Li and L. Wang  “Global stability in some SEIR models”  In _IMA Journal of Applied Mathematics_ 126, 2002, pp. 259–311
  * [38] Y. Li and J.. Muldowney  “Global stability for the SEIR model in epidemiology”, 1995, pp. 155–164
  * [39] A.. Lloyd  “The dependence of viral parameter estimates on the assumed viral life cycle: Limitations of studies of viral load data”  In _Proceedings of the Royal Society of London. Series B_ 268, 2001, pp. 847–854
  * [40] N. MacDonald  “Time Delays in Biological Models”  Heidelberg: Springer, 1978
  * [41] J.. Marsden and M. McCracken  “The Hopf Bifurcation and Its Applications”  In _Applied Mathematical Sciences_ 1 New York, NY, USA: Springer, 1976
  * [42] T.. McKeithan  “Kinetic proofreading in $T$ \- cell receptor signal transduction”  In _Proceedings of the National Academy of Sciences of the United States of America_ 92, 1995, pp. 5042–4046
  * [43] A.. McLean and T… Kirkwood  “A model of human immunodeficiency virus infection in T-helper cell clones”  In _Journal of Theoretical Biology_ 147, 1990, pp. 177
  * [44] A.. McLean and M.. Nowak  “Models of interaction between HIV and other pathogens”  In _Journal of Theoretical Biology_ 155, 1992, pp. 69
  * [45] A.. McLean et al.  “Resource competition as a mechanism for $B$ cell homeostasis”  In _Proceedings of the National Academy of Sciences of the United States of America_ 94, 1997, pp. 5792–5797
  * [46] J.. Mittler, M. Markowitz, D.. Ho and A.. Perelson  “Improved estimates for HIV-1 clearance rate and intracellular delay”  In _AIDS_ 13, 1999, pp. 1415
  * [47] J.. Mittler, B. Sulzer, A.. Neumann and A.S. Perelson  “Influence of delayed viral production on viral dynamics in HIV-1 infected patients”  In _Mathematical Biosciences_ 152 Elsevier, 1998, pp. 143
  * [48] J.. Muldowney  “Compound matrices and ordinary differential equations”  In _Rocky Mountain Journal of Mathematics_ 20, 1990, pp. 857–872
  * [49] A. Neumann et al.  “Hepatitis C viral dynamics in vivo and antiviral efficacy of the interferon-$\alpha$ therapy”  In _Science_ 282, 1998, pp. 103–107
  * [50] M.. Nowak and R.. Bangham  “Population dynamics of immune responses to persistent viruses”  In _Science_ 272, 1996, pp. 74
  * [51] M.. Nowak, A.. Lloyd and G.. al.  “Viral dynamics of primary viremia and antitroviral therapy in simian immunodeficiency virus infection”  In _Journal of Virology_ 71, 1997, pp. 7518–7525
  * [52] M.. Nowak and R.. May  “Mathematical biology of HIV infection: antigenic variation and diversity threshold”  In _Mathematical Biosciences_ 106, 1991, pp. 1
  * [53] M.. Nowak et al.  “Viral dynamics in hepatitis B virus infection”  In _Proceedings of the National Academy of Sciences of the United States of America_ 93, 1996, pp. 4398–4402
  * [54] J.. Percus, O.. Percus and A.. Perelson  “Predicting the size of the T -cell receptor and antibody combining region from consideration of efficient self–nonself discrimination”  In _Proceedings of the National Academy of Sciences of the United States of America_ 90, 1993, pp. 2691–1695
  * [55] A.. Perelson  “Dynamics of HIV Infection of CD4+ T-cells”  In _Mathematical Biosciences_ 114 Elsevier, 1993, pp. 81
  * [56] A.. Perelson  “Modelling the interaction of the immune system with HIV”  In _Mathematical and Statistical Approaches to AIDS Epidemiology_ Berlin: Springer, 1989, pp. 350
  * [57] A.. Perelson, P. Essunger and D.. Ho  “Dynamics of HIV-1 and CD4+ lymphocytes in vivo”  In _AIDS_ 11 (Suppl. A), 1997, pp. S17–S24
  * [58] A.. Perelson and P.. Nelson  “Mathematical analysis of HIV-1 dynamics in vivo”  In _SIAM Review_ 41, 1999, pp. 3
  * [59] A.. Perelson et al.  “Decay characteristics of HIV-1-infected compartments during combination therapy”  In _Nature_ 387, 1997, pp. 188–191
  * [60] A.. Perelson et al.  “HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time”  In _Science_ 271, 1996, pp. 1582
  * [61] Christopher Rackauckas and Qing Nie  “DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia”  In _Journal of Open Research Software_ 5 (1), 2017, pp. 15  DOI: 10.5334/jors.151
  * [62] R.. Root-Bernstein and S.. Merrill  “The necessity of cofactors in the pathogenesis of AIDS: a mathematical model”  In _Journal of Theoretical Biology_ 187, 1997, pp. 135
  * [63] H.. Smith  “Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems”  Providence, RI: American Mathematical Society, 1995
  * [64] H.. Smith and H. Thieme  “Convergence for strongly ordered preserving semiflows”  In _SIAM Journal on Mathematical Analysis_ 22, 1991, pp. 1081–1101
  * [65] X.. Song and S.. Cheng  “A delay-differential equation model of HIV infection of CD4+ $T$ -cells”  In _Journal of the Korean Mathematical Society_ 42 (5), 2005, pp. 1071–1086
  * [66] Xinyu Song and Avidan U. Neumann  “Global stability and periodic solution of the viral dynamics”  In _Journal of Mathematical Analysis and Applications_ 329 Elsevier, 2007, pp. 281–297  DOI: 10.1016/j.jmaa.2006.06.064
  * [67] J.. Spouge, R.. Shrager and D.. Dimitrov  “HIV-1 infection kinetics in tissue culture”  In _Mathematical Biosciences_ 138, 1996, pp. 1
  * [68] G. Stépán  “Retarded Dynamical Systems: Stability and Characteristic Functions”  UK: Longman, 1989
  * [69] N.. Stilianakis, K. Dietz and D. Schenzle  “Analysis of a model for the pathogenesis of AIDS”  In _Mathematical Biosciences_ 145, 1997, pp. 27
  * [70] J. Tam  “Delay effect in a model for virus replication”  In _IMA Journal of Mathematics Applied in Medicine and Biology_ 16, 1999, pp. 29
  * [71] International Committee Taxonomy of Viruses  “Taxonomy” Updated July 2019, Accessed December 13, 2020, National Institutes of Health  URL: https://talk.ictvonline.org/taxonomy/
  * [72] X. Wei et al.  “Viral dynamics in human immunodeficiency virus type 1 infection”  In _Nature_ 373, 1995, pp. 117
  * [73] Wikipedia contributors  “Poincaré–Bendixson theorem — Wikipedia, The Free Encyclopedia” [Online; accessed 2-December-2020], 2020  URL: https://en.wikipedia.org/w/index.php?title=Poincar%5C%C3%5C%A9%5C%E2%5C%80%5C%93Bendixson_theorem&oldid=970525369
  * [74] Junyuan Yang, Xiaoyan Wang and Fengqin Zhang  “A Differential Equation Model of HIV Infection of CD4+ T-Cells with Delay” Article ID 903678, 16 pages  In _Discrete Dynamics in Nature and Society_ Hindawi, 2008  DOI: 10.1155/2008/903678
  * [75] Xueyong Zhou, Xinyu Song and Xiangyun Shi  “A differential equation model of HIV infection of CD4+ $T$ \- cells with cure rate”  In _Journal of Mathematical Analysis and Applications_ 342 Elsevier, 2008, pp. 1342–1355  DOI: 10.1016/j.jmaa.2008.01.008
  * [76] H.. Zhu, H.. Smith and M.. Hirsch  “Stable periodic orbits for a class of three dimensional competitive systems”  In _Journal of Differential Equations_ 110, 1994, pp. 143–156