![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
ABSTRACT
Time delays can affect the dynamics of HIV infection predicted by mathematical models. In this paper, we studied two mathematical models each with two time delays. In the first model with HIV latency, one delay is the time between viral entry into a cell and the establishment of HIV latency, and the other delay is the time between cell infection and viral production. We defined the basic reproductive number and showed the local and global stability of the steady states. Numerical simulations were performed to evaluate the influence of time delays on the dynamics. In the second model with HIV immune response, one delay is the time between cell infection and viral production, and the other delay is the time needed for the adaptive immune response to emerge to control viral replication. With two positive delays, we obtained the stability crossing curves for the model, which were shown to be a series of open-ended curves.
1. Introduction
Human immunodeficiency virus type 1 (HIV-1) infects CD4+ T cells, a substantial component of the immune system, and the infection goes through three distinct stages (see review in [Citation31]): primary infection, asymptomatic stage, and acquired immunodeficiency syndrome (AIDS). The primary infection stage lasts a few weeks during which the viral load in blood increases rapidly to the peak level, followed by decline to a steady state. CD4+ T cells decrease initially because of the infection and then undergo a minor increase as the primary infection ends. A high level of viral infection stimulates the development of adaptive immune responses, including CD8+ T cells, which can kill infected cells [Citation55]. The asymptomatic infection stage can take several years, in which there is usually a very slow increase in the viral load and a slow decline in CD4+ T cells. As CD4+ T cells decline to a low level, for example, 200 cells per microlitre, the AIDS stage begins [Citation10]. Infected individuals have a rapid increase in plasma viral load and are also accompanied with various opportunistic infections because of the collapse of the immune system. Current treatment consisting of several antiretroviral drugs can suppress viral load to below the detection limit of standard assays but cannot eradicate the virus. An important reason is that HIV can persist as provirus in resting memory CD4+ T cells [Citation7, Citation13, Citation56]. These latently infected CD4+ T cells live relatively long, are not affected by antiretroviral therapy, but can be activated by their relevant antigens to produce new virions [Citation4].
Mathematical models have been developed to study the dynamics of virions and CD4+ T cells during infection, the influence of antiretroviral therapy, the emergence of drug resistance, immune responses, and low viral load persistence in patients receiving prolonged antiretroviral therapy (see reviews in [Citation34–36, Citation43]). Many models assumed that cells produce virus immediately after they are infected. However, the life cycle of the virus contains a number of intracellular processes; for example, viral entry into a target cell, reverse transcription from HIV RNA to DNA, integration of viral DNA into the host cell's DNA, transcription and translation, assembly and virus budding from the infected cell [Citation34]. A time delay representing the time between cell infection and viral production has been incorporated into models to study its influence on virus dynamics. Assuming that the target CD4+ T cell level is constant and that the treatment is 100% effective in the delay model, Herz et al. obtained the analytic expression of the viral load decline under treatment and used it to analyse the viral load decline data in patients [Citation18]. Nelson et al. extended the analysis and showed that the time delay can affect the estimate of the death rate of infected cells when treatment is not 100% effective [Citation29]. A distributed continuous time delay was used to account for the different time between cell infection and viral production [Citation15, Citation27, Citation28]. More models including the intracellular delay can be found in refs. [Citation3, Citation11, Citation12, Citation20–23, Citation25, Citation49, Citation51–54, Citation57] and references cited therein.
The adaptive immune responses need time to emerge to control viral replication [Citation55]. Ciupe et al. [Citation8] included a time delay representing the time needed to activate the CD8+ T cell response in a model. Pawelek et al. [Citation33] incorporated two delays, one the time needed for infected cells to produce virus and the other the time needed for CD8+ T cells to emerge to kill infected cells. The model was fit to the viral load data from patients during primary infection. Some special cases of the model with two delays have been mathematically analysed [Citation33]. However, the stability switching curves for the infected steady state when both delays are positive have not been investigated.
In this paper, we will study two HIV models each with two time delays. In the first model, we incorporate two time delays into an HIV latent infection model which was introduced to evaluate the contribution of latently infected cell activation to HIV production [Citation38]. One delay is the time between viral entry into a cell and the establishment of HIV latency (i.e. the viral DNA has been integrated into the host cell's DNA), and the other delay is the time between initial cell infection and viral production. We will investigate the local and global stability of the steady states. Numerical simulations will be performed to evaluate the influence of time delays on the dynamics. The second model of HIV immune response with two delays is the same as that in ref. [Citation33]. Some special cases of the model have been studied before [Citation33]. Here we will use the method developed by Gu et al. [Citation16] to determine the stability crossing curves for the infected steady state of the model when both delays are positive.
2. An HIV latent infection model with two delays
In this section, we study an HIV latent infection model with two time delays. The model without delays was introduced to analyse the biphasic viral load decline observed in patients receiving combination therapy of two reverse transcriptase inhibitors and one protease inhibitor [Citation38]. It was mathematically analysed in ref. [Citation32]. The model assumes that once infected, the cell immediately becomes either a productively infected cell that can produce virus, or becomes a latently infected cell that can be activated to produce virus. However, it takes time for the latent infection to be established and also for productively infected cells to produce virus. We introduce two time delays into the model. One is the time between viral entry and latent infection (i.e. the integration of viral DNA into cell's DNA has finished) and the other is the time between viral entry and viral production. The first delay is denoted by and the second delay is denoted by
. Obviously, we have
from the viral life cycle. The model can be described by the following system of delayed differential equations.
(1)
(1)
In the model, is the concentration of uninfected CD4+ T cells at time t,
denotes the concentration of latently infected T cells at time t,
is the concentration of productively infected cells, and
represents the concentration of virions in plasma at t. The parameter s is the generate rate of uninfected CD4+ T cells,
is the per capita death rate of uninfected cells, and β is the infection rate of target cell by virus. A small fraction (f) of infected cells are assumed to result in latency and the remaining become productively infected cells. Latently infected cells die at a rate constant
per cell and productively infected cells die at a rate constant δ per cell. N is the viral burst size, which is the total number of virions released by one infected cell in its lifespan, and c is the viral clearance rate. Latently infected cells can be activated by their relevant antigens to become productively infected cells at a rate constant α. The parameter
represents the death rate of infected cells that have not started to produce virus. Thus,
and
represent the probability of an infected cells surviving the time
and
, respectively, following viral entry. Parameters and values used in simulations are summarized in Table .
Table 1. Parameters and values used in the models and simulations.
We derive the basic reproductive number using the next-generation method. From the infection and viral production term in the model, we define matrices and
as follows
where
is the steady-state level of target cells before infection. The basic reproductive number,
, can be derived as the spectral radius of the next generation operator
. Straightforward calculation yields
Thus, is given by
(2)
(2) Notice that
is the basic reproductive number of the basic model (i.e. model (Equation1
(1)
(1) ) without latently infected cells and without time delays, see ref. [Citation44]),
is the contribution to the basic reproductive number from productively infected cells, and
represents the contribution from activation of latently infected cells. Therefore,
is the basic reproductive number of model (Equation1
(1)
(1) ).
The model has two equilibrium points. One is the infection-free equilibrium and the other is the infected equilibrium
, where
(3)
(3) It is clear that the infected steady state exists if and only if
.
3. Stability analysis
3.1. Local stability
In this section, we study the local stability of the infection-free equilibrium and the infected equilibrium points. The following result suggests that the basic reproductive number provides a threshold value determining if the infection is predicted to die out or persist.
Theorem 3.1.
The infection-free equilibrium of model (Equation1
(1)
(1) ) is locally asymptotically stable when
and unstable when
.
Proof.
We linearize the system (Equation1(1)
(1) ), obtain the Jacobian matrix J at the infection-free steady state
, and calculate the matrix
as follows, where λ is the eigenvalue and I is the identity matrix with the same dimension as
.
Thus, the characteristic equation at the equilibrium
is given by
(4)
(4) When
, the infection-free equilibrium
of Equation (Equation1
(1)
(1) ) was shown to be locally asymptotically stable when
and unstable when
according to the Routh–Hurwitz stability criterion [Citation32]. Next, we show that the conclusion holds for any
and
.
It is clear that Equation (Equation4(4)
(4) ) has a negative root
. All the other roots are determined by the following equation
(5)
(5)
For simplicity, we introduce a new notation
(6)
(6)
Obviously, we have 0<g<1. Using the notation g and the basic reproductive number , the characteristic equation (Equation5
(5)
(5) ) can be simplified as
(7)
(7)
Suppose that the real part of the eigenvalue λ is non-negative. Taking the modulus of the left-hand side of Equation (Equation7(7)
(7) ) and also using Cauchy–Schwartz inequality, we obtain
(8)
(8)
On the other hand, the modulus of the right-hand side of Equation (Equation7(7)
(7) ) satisfies
(9)
(9) That leads to a contradiction. It follows that all the roots of the characteristic equation (Equation4
(4)
(4) ) have negative real parts. Thus, the infection-free steady state is locally asymptotically stable when
.
Next, we turn to the case of with delays
and
. The characteristic equation (Equation4
(4)
(4) ) can be written as the following quasipolynomial
(10)
(10) where
When
, we define a function
. It can be seen that
as
and
. Thus, there is at least one positive real root. This shows that the inflection-free steady-state
is unstable when
.
From the expressions in Equation (Equation3(3)
(3) ), we know that the infected steady state exists if and only if
. The following result shows that the chronic infection is established when
.
Theorem 3.2.
The infected equilibrium of Equation (Equation1
(1)
(1) ) is locally asymptotically stable when it exists
i.e. when
.
Proof.
Similar to the infection-free steady state, we obtain the Jacobian matrix at the infected steady-state and calculate the matrix
, given by
where λ is the eigenvalue,
and
are given in Equation (Equation3
(3)
(3) ). Using the basic reproductive number
, the characteristic equation at the infected steady state
is
(11)
(11)
When the infected equilibrium
of model (Equation1
(1)
(1) ) was shown to be locally asymptotically stable when
(see [Citation32]). We consider the case of
with
and
. Suppose that the eigenvalue λ has a non-negative real part. The characteristic equation (Equation11
(11)
(11) ) can be written as
(12)
(12) Using the notation g defined in Equation (Equation6
(6)
(6) ), the above equation can be rewritten as
(13)
(13)
Using similar arguments as in Theorem 3.1, it can be shown that the modulus of the left-hand side of Equation (Equation13(13)
(13) ) satisfies
(14)
(14) On the other hand, the modulus of the right-hand side of Equation (Equation13
(13)
(13) ) satisfies
(15)
(15) Equation (Equation13
(13)
(13) ) and inequalities (Equation14
(14)
(14) ) and (Equation15
(15)
(15) ) lead to a contradiction. Thus, the infected steady-state
is locally asymptotically stable whenever it exists, that is,
.
3.2. Global stability
In this section, we use the Lyapunov method to study the global stability of steady states of model (Equation1(1)
(1) ). The following result suggests that when
the infection is predicted to die out regardless of the initial condition.
Theorem 3.3.
The infection-free equilibrium of model (Equation1
(1)
(1) ) is globally asymptotically stable when
.
Proof.
We define a Lyapunov functional
where
and
The Lyapunov functional W is non-negative. At the infection-free equilibrium
, we have
.
Calculating the time derivative of along the solution of model (Equation1
(1)
(1) ), we obtain
In view of , we have
It is clear that
when
, and that
if and only if
and
. The largest invariant set in
is
. Therefore, by Lyapunov–LaSalle Invariance Principle [Citation17] and Theorem 3.1, the infection-free equilibrium
is globally asymptotically stable.
The following result shows that when the chronic infection is always established regardless of the initial condition.
Theorem 3.4.
The infected equilibrium of model (Equation1
(1)
(1) ) is globally asymptotically stable when
.
Proof.
We define a different Lyapunov functional
where
and
The Lyapunov functional
is non-negative. At the infected steady-state
, we have
.
The time derivative of along the solution of system (Equation1
(1)
(1) ) yields
From the infected steady state, we have
Using the above equations, we obtain
which is equal to
Let
, where x>0. The function
is non-negative for any x>0, and
if and only if x=1. Using this property, we have
The largest invariant set in
is
. Thus, by Lyapunov–LaSalle Invariance Principle and Theorem 3.2, the infected steady-state
is globally asymptotically stable whenever it exists.
4. Numerical results of model (1)
In this section, we present some numerical results for model (Equation1(1)
(1) ). For modelling simulation, we chose parameter values based on existing experimental data and our previous modelling studies [Citation41, Citation44, Citation45]. In a healthy individual, the CD4+ T cell level is approximately
cells per μ l [Citation5]. We changed the unit to cells/ml and assumed the initial CD4+ T cell level to be
cells/ml before infection. The death rate of uninfected CD4+ T cells was assumed to be
day−1 [Citation37]. Thus, the generation rate of CD4+ T cells can be calculated as
cells ml−1 day−1. The viral infection rate β was chosen to be
ml virion−1 day−1 [Citation44]. A very small fraction of new infection results in latent infection. We assumed that the fraction of latency is f=0.001, as used in the ref. [Citation41]. The death rate of productively infected T cells is
day−1 [Citation26]. The total number of virions produced by one infected cell during its lifespan was assumed to be N=2000 virions cell−1 [Citation44]. The viral clearance rate was assumed to be c=23 day−1 [Citation39]. The death rate of latently infected cells was assumed to be
day−1, which corresponds to a half life of about 6 months, consistent with the half life of memory T cells [Citation43]. The activation rate was assumed to be
day−1 [Citation41]. The death rate (
) of infected cells that have not started to produce virus should be between
and δ, and was assumed to be 0.05 day−1. The time delay
between viral entry and viral DNA integration was chosen to be 0.25 days and the delay
between viral entry and viral production was chosen to be 0.5 days [Citation45]. Parameter values are summarized in Table .
We chose the initial CD4+ T cell count cells per ml, initial viral load
virions per ml,
, and
. Numerical simulation of model (Equation1
(1)
(1) ) shows that the solution approaches an infected steady state after some oscillations (Figure ). The steady states of uninfected CD4+ T cells, latently infected cells, productively infected cells, and viral load are
cells per ml, 358.69 cells per ml,
cells per ml, and
viral RNA copies per ml, respectively. The basic reproductive number using the above parameter values is 2.03. According to the stability results, the infected steady state is globally asymptotically stable, which is illustrated by the numerical simulation shown in Figure .
Figure 1. Numerical simulation of model (Equation1(1)
(1) ) without treatment. The time delay
was fixed to be 0.25 days and
was fixed to be 0.5 days [Citation45]. All the other parameter values are the same as those listed in Table .
![Figure 1. Numerical simulation of model (Equation1(1) dTdt=s−dTT−βVT,dLdt=fβV(t−τ1)T(t−τ1)e−δ1τ1−δLL−αL,dIdt=(1−f)βV(t−τ2)T(t−τ2)e−δ1τ2−δI+αL,dVdt=NδI−cV.(1) ) without treatment. The time delay τ1 was fixed to be 0.25 days and τ2 was fixed to be 0.5 days [Citation45]. All the other parameter values are the same as those listed in Table 1.](/cms/asset/8d816355-1498-4da7-9d3c-690e1d1b5411/tjbd_a_1148202_f0001_c.jpg)
In Figure , we showed the simulation of model (Equation1(1)
(1) ) under drug treatment. The overall drug efficacy was assumed to be 0.9, which means that the infection rate β is reduced by 90% [Citation42]. Using the above steady states as the initial conditions of the model under treatment, the simulation shows that the CD4+ T cell level rebounds to the pre-infection level, and that all other infected cells and viral load decline to 0. In this case, the infection is predicted to die out. The effective reproductive number under treatment is 0.2. The simulation confirms the stability of the infection-free steady state when
.
Figure 2. Numerical simulation of model (Equation1(1)
(1) ) under treatment. The infection rate β was assumed to be reduced by 90%. The initial conditions were assumed to be the steady states of the model before treatment (see Figure ). All the other parameter values are the same as those listed in Table .
![Figure 2. Numerical simulation of model (Equation1(1) dTdt=s−dTT−βVT,dLdt=fβV(t−τ1)T(t−τ1)e−δ1τ1−δLL−αL,dIdt=(1−f)βV(t−τ2)T(t−τ2)e−δ1τ2−δI+αL,dVdt=NδI−cV.(1) ) under treatment. The infection rate β was assumed to be reduced by 90%. The initial conditions were assumed to be the steady states of the model before treatment (see Figure 1). All the other parameter values are the same as those listed in Table 1.](/cms/asset/851b0dd4-51bc-43b1-98dc-ba33c983370e/tjbd_a_1148202_f0002_c.jpg)
We also studied the effect of different time delays on the dynamics of model (Equation1(1)
(1) ) in Figures and . In Figure , the delay
was fixed to be 0.5 days and
varied from 0.1 to 0.4 days. We found there is almost no difference in the dynamics of cells and viral load. In Figure , the delay
was fixed to be 0.25 days and
varied from 0.5 to 0.8 days. As the delay
increases, there is a delay in the peak of target cells, infected cells, and viral load. However, there are very small differences in the steady states. This can be observed from the steady states shown in Equation (Equation3
(3)
(3) ) and the basic reproductive number shown in Equation (Equation2
(2)
(2) ).
Figure 3. Effect of different time delays () on the dynamics of model (Equation1
(1)
(1) ). The delay
was chosen to be 0.1, 0.2, 0.3, or 0.4 days, and the delay
was fixed at 0.5 days. All the other parameter values are the same as those in Figure . The lines of different cases overlap.
![Figure 3. Effect of different time delays (τ1) on the dynamics of model (Equation1(1) dTdt=s−dTT−βVT,dLdt=fβV(t−τ1)T(t−τ1)e−δ1τ1−δLL−αL,dIdt=(1−f)βV(t−τ2)T(t−τ2)e−δ1τ2−δI+αL,dVdt=NδI−cV.(1) ). The delay τ1 was chosen to be 0.1, 0.2, 0.3, or 0.4 days, and the delay τ2 was fixed at 0.5 days. All the other parameter values are the same as those in Figure 1. The lines of different cases overlap.](/cms/asset/dac91493-a225-4266-b36f-13cef90c2792/tjbd_a_1148202_f0003_c.jpg)
Figure 4. Effect of different time delays () on the dynamics of model (Equation1
(1)
(1) ). The delay
was chosen to be 0.5, 0.6, 0.7, or 0.8 days, and the delay
was fixed at 0.25 days. All the other parameter values are the same as those in Figure .
![Figure 4. Effect of different time delays (τ2) on the dynamics of model (Equation1(1) dTdt=s−dTT−βVT,dLdt=fβV(t−τ1)T(t−τ1)e−δ1τ1−δLL−αL,dIdt=(1−f)βV(t−τ2)T(t−τ2)e−δ1τ2−δI+αL,dVdt=NδI−cV.(1) ). The delay τ2 was chosen to be 0.5, 0.6, 0.7, or 0.8 days, and the delay τ1 was fixed at 0.25 days. All the other parameter values are the same as those in Figure 1.](/cms/asset/d462bfcc-e0f2-4158-8d06-e946d1414747/tjbd_a_1148202_f0004_c.jpg)
5. An HIV immune model with two delays
Pawelek et al. [Citation33] studied a model of HIV infection with two time delays. One is the time () needed for infected cells to produce virus after viral entry and the other is the time (
) needed for the adaptive immune response to emerge to control viral replication. The model was fit to the viral load data from 10 patients during primary HIV infection. The model is described by the following system of differential equations.
(16)
(16) In the model,
is the concentration of effector cells at time t, which can kill productively infected cells. The variables T, I and V are the same as those in model (Equation1
(1)
(1) ). The parameter β is the infection rate and
. The parameter
is the killing rate of infected cells by effector cells, p is the generation rate of effector cells stimulated by productively infected cells, and
is the death rate of effector cells. The other parameters are the same as those in model (Equation1
(1)
(1) ). Parameters and values are also summarized in Table .
The model has two steady states: the infection-freesteady-state with
and the infected steady-state
where
(17)
(17)
From , we know that the infected steady state exists if and only if
, which is equivalent to
or
. Note that
is the basic reproductive number of the basic model without the immune response [Citation44].
It has been shown that for the infection-free steady state is locally asymptotically stable when the basic reproductive number
[Citation33]. The infected steady state is locally asymptotically stable when
in the case of
. When
, the analysis of the stability of the infected steady state is challenging. Pawelek et al. [Citation33] showed that even for a simplified model (assuming that the target cell level is constant
in Equation (Equation16
(16)
(16) )), a positive immune delay
is able to destabilize the infected steady state. Specifically, in the case of
, it was shown that the infected steady state is locally asymptotically stable for
and a Hopf bifurcation occurs when
, where
is a threshold value of the immune delay [Citation33].
In this paper, we will investigate the case of and
. Because
is relatively small (i.e.
), we assume
and consider the following simplified model. Stability crossing curves of the infected steady state will be explored for this model.
(18)
(18)
At the infected steady state, the characteristic equation of model (Equation18(18)
(18) ) is
(19)
(19) where λ is the eigenvalue and
. Expanding the product, the above equation can be written as the following equation
where
(20)
(20) Thus, the stability of the infected steady state is determined by the zeros of the characteristic quasipolynomial
The main objective of the remaining part of this section is to identify the stability crossing curves of
in
(i.e. the set of two-dimensional vectors with non-negative components).
We start with the definition of stability crossing curve, which was used in Gu et al. [Citation16].
Definition 5.1.
The stability crossing curve is the set of all crossing points
in
where
is the point such that
has at least one zero on the imaginary axis. Thus,
can be written as
.
From Equation (Equation20(20)
(20) ), it is easy to show that
does not have zeros on the imaginary axis. Let
. We have
, where
and
. Thus, for given
and
,
and
have the same zeros in a neighbourhood of the imaginary axis. Instead of using
, we will study the crossing points from the solutions of the following equation
(21)
(21)
From Equation (Equation20(20)
(20) ), we have that for model (Equation18
(18)
(18) )
(22)
(22) and
(23)
(23)
Next, we recall a few propositions that can be used to characterize the stability crossing curves of a general linear scalar system with two time delays. The proof of these propositions can be found in Gu et al. [Citation16].
Proposition 5.2.
For each
can be a solution of
for some
if and only if
(24)
(24)
(25)
(25) For
satisfying
can be zero of
for some
if and only if
For our model (Equation18(18)
(18) ),
cannot be 0 for any
. Thus,
is a solution of
for some
if and only if
and
satisfy inequalities (Equation24
(24)
(24) ) and (Equation25
(25)
(25) ).
We define the crossing set Ω as the set of all that satisfy Equations (Equation24
(24)
(24) ) and (Equation25
(25)
(25) ). From Proposition 5.2, we have
. Furthermore, all the pairs of
satisfying
Equation (Equation21
(21)
(21) ) can be calculated as follows [Citation16]:
(26)
(26)
(27)
(27) where ∠ represents the argument of a complex number,
are given by
(28)
(28)
(29)
(29) and
,
,
,
are the smallest possible integers such that the corresponding
calculated are non-negative. Note that
and
are mirror images with respect to the real axis.
The following result regarding the crossing set Ω was obtained in ref. [Citation16].
Proposition 5.3.
The crossing set Ω consists of a finite number of intervals of finite length.
Let these intervals be ,
, which are arranged in an order such that the left-end point of
increases with increasing k. Then
.
Let
The stability crossing curves
can be written as the union of
,
, that is,
.
Let ,
, where
and
are the left and right end points, respectively. The end points of the intervals
,
, can be classified into the following three types (Type 1, Type 2, and Type 3). Type 0 is an additional type if the left-end point is equal to 0 [Citation16].
Type 1. is satisfied. In this case,
, and
is connected with
at this end.
Type 2. is satisfied. In this case,
, and
is connected with
at this end.
Type 3. is satisfied. In this case,
, and
is connected with
at this end.
Type 0. is equal to 0. In this case,
satisfies Equations (Equation24
(24)
(24) ) and (Equation25
(25)
(25) ). As
,
and
approach ∞ with asymptotes passing through the points
with slopes of
(30)
(30) where
and
are evaluated by Equations (Equation28
(28)
(28) ) and (Equation29
(29)
(29) ) using
and
, respectively, and
An interval of the crossing set Ω is said to be of type lr if the left-end point of
is of type l and the right-end point is of type r.
It follows from Equations (Equation22(22)
(22) ) and (Equation23
(23)
(23) ) that
and
. Because of
, we have
and
. Thus,
satisfies Equations (Equation24
(24)
(24) ) and (Equation25
(25)
(25) ), which means that 0 is the left-end point of the interval
, belonging to type 0. From Equations (Equation22
(22)
(22) ) and (Equation23
(23)
(23) ), it is also easy to show that
Thus, there is at least one solution ω of the equation
.
In the upper panel of Figure , we plotted (solid curve) and
(dashed curve), respectively, as a function of ω. The basic reproductive number was chosen to be
, the death rate of infected cells was set to
day−1, the viral clearance rate was chosen to be c=23 day−1, and the death rate of effector cells was chosen to be
day−1. These parameter values were chosen based on the fitting of the two-delay model (Equation Equation16
(16)
(16) ) to the data from 10 patients during primary infection [Citation33]. We found there is only one intersection point
between
and y=1. Therefore, there is only one interval
in the crossing set Ω, which is of type 03. Corresponding to
, there are only one type of stability crossing curves
. Furthermore, for type 03,
and
are connected at the left-end point
. The other ends of
and
extend to infinity with slopes
and
(given in Equation (Equation30
(30)
(30) )), respectively. Therefore,
consists of a series of open-ended curves with both ends approaching ∞. These open-ended curves in
were plotted in the lower panel of Figure for the infected steady state of model (Equation18
(18)
(18) ).
Figure 5. Upper panel: plots of (see Equations (Equation22
(22)
(22) ) and (Equation23
(23)
(23) )) as a function of ω. The intersection point of
and 1 is 0.69. Lower panel: stability crossing curves corresponding to the crossing interval (0, 0.69] for the two-delay model (Equation18
(18)
(18) ). The basic reproductive number was chosen to be
. The other parameters were chosen based on the estimates from data fitting in ref. [Citation33]:
day−1, c=23 day−1, and
day−1.
![Figure 5. Upper panel: plots of |a1(jω)|±|a2(jω)| (see Equations (Equation22(22) a1(λ)=−R0cδλ−dER0cδλ3+(R0δ+dE+c)λ2+[dER0δ+c(R0δ+dE)]λ+cdER0δ,(22) ) and (Equation23(23) a2(λ)=(R0−1)dEδλ+c(R0−1)dEδλ3+(R0δ+dE+c)λ2+[dER0δ+c(R0δ+dE)]λ+cdER0δ.(23) )) as a function of ω. The intersection point of |a1(jω)|+|a2(jω)| and 1 is 0.69. Lower panel: stability crossing curves corresponding to the crossing interval (0, 0.69] for the two-delay model (Equation18(18) dIdt=βV(t−τ1)T0−δI−dxEI,dVdt=NδI−cV,dEdt=pI(t−τ2)−dEE.(18) ). The basic reproductive number was chosen to be R0=2.6. The other parameters were chosen based on the estimates from data fitting in ref. [Citation33]: δ=0.3 day−1, c=23 day−1, and dE=0.45 day−1.](/cms/asset/a870ad0f-608c-4c49-9159-95909eccb682/tjbd_a_1148202_f0005_c.jpg)
The upper panel of Figure showed the changes of and
when the basic reproductive number was chosen to be
. The other parameter values are the same as those in Figure . There is again only one solution
of the equation
. Thus, there is only one crossing interval
, which is of type 03. The open-ended stability crossing curves corresponding to this interval were shown in the lower panel of Figure .
6. Conclusion and discussion
This paper investigates the asymptotical behaviour of the solutions of delay differential equation systems that are used to study HIV infection dynamics. The first model was used to determine the contribution of latently infected cell activation to HIV production during combination antiretroviral therapy [Citation38]. We incorporated two delays into the model. One is the time between initial cell infection and the establishment of latent infection, and the other is the time between viral entry into cell and viral production in productively infected cells. Although the second delay was included to evaluate its influence on HIV dynamics and parameter estimates in the literature [Citation18, Citation27–29], the first time delay has not been studied before. We derived the basic reproductive number and performed both local and global stability analysis of the steady states. We showed that the infection-free steady state is globally asymptotically stable when the basic reproductive number is less than 1, and the infected steady state is globally asymptotically stable when the basic reproductive number is greater than 1. Thus, incorporation of the latent infection delay does not destabilize the stability of the steady state. This is consistent with the results of models including intracellular delay [Citation22, Citation23]. Numerical results also confirm the stability analysis. Simulation of the model with different delays shows that the delay between viral entry and viral integration has a very minor effect on the dynamics of cells and viral load. The time delay between viral infection and viral production causes a very small change in the steady states of infected cells and viral load, but postpones the time to reach the viral peak.
Figure 6. Upper panel: plots of as a function of ω. The intersection point of
and 1 is 1.12. Lower panel: stability crossing curves corresponding to the crossing interval (0, 1.12] for the two-delay model (Equation18
(18)
(18) ). The basic reproductive number was chosen to be
. The other parameters are the same as those in Figure .
![Figure 6. Upper panel: plots of |a1(jω)|±|a2(jω)| as a function of ω. The intersection point of |a1(jω)|+|a2(jω)| and 1 is 1.12. Lower panel: stability crossing curves corresponding to the crossing interval (0, 1.12] for the two-delay model (Equation18(18) dIdt=βV(t−τ1)T0−δI−dxEI,dVdt=NδI−cV,dEdt=pI(t−τ2)−dEE.(18) ). The basic reproductive number was chosen to be R0=5. The other parameters are the same as those in Figure 5.](/cms/asset/a3505172-32c1-4088-8579-07d9f955d81f/tjbd_a_1148202_f0006_c.jpg)
The model analysis suggests if treatment is very effective then the basic reproductive number will become less than 1, and both the latent reservoir and the virus are predicted to die out. However, current treatment consisting of several antiretroviral drugs cannot eradicate the latent reservoir (and the virus). One reason may be that current treatment regimens are not sufficiently effective to block all the infection and thus there exists residual ongoing viral replication [Citation40]. Another reason is the intrinsic stability of latently infected CD4+ T cells [Citation4]. Some explanations such as asymmetric cell division [Citation42], programmed expansion and contraction [Citation41], and homeostatic proliferation [Citation41] of latently infected cells have been proposed and tested by mathematical models. Homeostatic proliferation has been experimentally confirmed [Citation6]. Some other models studying the dynamics of the latent reservoir and low-level viral load persistence during effective antiretroviral therapy can be found in refs. [Citation2, Citation14, Citation48, Citation50] and references cited therein. Reviews of mathematical models and implications for treatment can be found in refs. [Citation36, Citation43].
The second model also includes two delays. One is the time needed for infected cells to produce virus following cell infection (i.e. , the intracellular delay), and the other is the time needed for the CD8+ T cell response to emerge to control viral replication (i.e.
, the immune delay). Some special cases of the immune model with two delays have been investigated in a previous study [Citation33]. However, the analysis of the model with two positive delays is challenging. We showed that introducing the intracellular delay does not change the stability results but including the immune delay in the model generates rich dynamics. A Hopf bifurcation was shown to take place at a critical threshold of the immune delay for a simplified model assuming that the target cell level is constant [Citation33]. When the two delays are both positive, the stability crossing curves of the HIV immune model have not been characterized. Several papers in the literature have studied the properties of the solutions of models with two delays for different problems [Citation1, Citation9, Citation19, Citation24, Citation30, Citation46, Citation47]. Gu et al. investigated a general linear scalar system with two delays and showed that the stability crossing curves fall into three categories of curves, namely, closed curves, open-ended curves, and spiral-like curves oriented horizontally, vertically, or diagonally [Citation16]. We applied the method to the characteristic quasipolynomial for the HIV immune model, and found that the stability crossing curves consist of a series of open-ended curves in the delay-parameter space. Overall, this analysis helps to better understand the dynamical behaviour of HIV models with two time delays.
Acknowledgments
The authors would like to thank Professor Keqin Gu for providing help in the simulation of the stability crossing curves. The work was finished when X. Wang visited Oakland University in 2015.
Disclosure statement
No potential conflict of interest was reported by the authors.
Additional information
Funding
References
- M. Adimy, F. Crauste, and S. Ruan, Periodic oscillations in leukopoiesis models with two delays, J. Theoret. Biol. 242 (2006), pp. 288–299. doi: 10.1016/j.jtbi.2006.02.020
- N.M. Archin, N.K. Vaidya, J.D. Kuruc, A.L. Liberty, A. Wiegand, M.F. Kearney, M.S. Cohen, J.M. Coffin, R.J. Bosch, C.L. Gay, J.J. Eron, D.M. Margolis, and A.S. Perelson, Immediate antiviral therapy appears to restrict resting CD4+ cell HIV-1 infection without accelerating the decay of latent infection, Proc. Natl. Acad. Sci. USA 109 (2012), pp. 9523–9528. doi: 10.1073/pnas.1120248109
- H.T. Banks, D.M. Bortz, and S.E. Holte, Incorporation of variability into the modeling of viral delays in HIV infection dynamics, Math. Biosci. 183 (2003), pp. 63–91. doi: 10.1016/S0025-5564(02)00218-3
- J.N. Blankson, D. Persaud, and R.F. Siliciano, The challenge of viral reservoirs in HIV-1 infection, Annu. Rev. Med. 53 (2002), pp. 557–593. doi: 10.1146/annurev.med.53.082901.104024
- M. Bofill, G. Janossy, C.A. Lee, D. MacDonald-Burns, A.N. Phillips, C. Sabin, A. Timms, M.A. Johnson, and P.B. Kernoff, Laboratory control values for CD4 and CD8 T lymphocytes: Implications for HIV-1 diagnosis, Clin. Exp. Immunol. 88 (1992), pp. 243–252. doi: 10.1111/j.1365-2249.1992.tb03068.x
- N. Chomont, M. El-Far, P. Ancuta, L. Trautmann, F.A. Procopio, B. Yassine-Diab, G. Boucher, M.-R. Boulassel, G. Ghattas, J.M. Brenchley, T.W. Schacker, B.J. Hill, D.C. Douek, J.-P. Routy, E.K. Haddad, and R.-P. Sékaly, HIV reservoir size and persistence are driven by T cell survival and homeostatic proliferation, Nat. Med. 15 (2009), pp. 893–900. doi: 10.1038/nm.1972
- T. -W. Chun, L. Stuyver, S.B. Mizell, L.A. Ehler, J.A.M. Mican, M. Baseler, A.L. Lloyd, M.A. Nowak, and A.S. Fauci, Presence of an inducible HIV-1 latent reservoir during highly active antiretroviral therapy, Proc. Natl. Acad. Sci. USA 94 (1997), pp. 13193–13197. doi: 10.1073/pnas.94.24.13193
- M.S. Ciupe, B.L. Bivort, D.M. Bortz, and P.W. Nelson, Estimating kinetic parameters from HIV primary infection data through the eyes of three different mathematical models, Math. Biosci. 200 (2006), pp. 1–27. doi: 10.1016/j.mbs.2005.12.006
- K.L. Cooke and P. van den Driessche, Analysis of an SEIRS epidemic model with two delays, J. Math. Biol. 35 (1996), pp. 240–260. doi: 10.1007/s002850050051
- S.M. Crowe, J.B. Carlin, K.I. Stewart, C.R. Lucas, and J.F. Hoy, Predictive value of CD4 lymphocyte numbers for the development of opportunistic infections and malignancies in HIV-infected persons, J. Acquir. Immune Defic. Syndr. 4 (1991), pp. 770–776.
- R.V. Culshaw, S. Ruan, and G. Webb, A mathematical model of cell-to-cell spread of HIV-1 that includes a time delay, J. Math. Biol. 46 (2003), pp. 425–444. doi: 10.1007/s00285-002-0191-5
- N.M. Dixit and A.S. Perelson, Complex patterns of viral load decay under antiretroviral therapy: Influence of pharmacokinetics and intracellular delay, J. Theoret. Biol. 226 (2004), pp. 95–109. doi: 10.1016/j.jtbi.2003.09.002
- D. Finzi, M. Hermankova, T. Pierson, L.M. Carruth, C. Buck, R.E. Chaisson, T.C. Quinn, K. Chadwick, J. Margolick, R. Brookmeyer, J. Gallant, M. Markowitz, D.D. Ho, D.D. Richman, and R.F. Siliciano, Identification of a reservoir for HIV-1 in patients on highly active antiretroviral therapy, Science 278 (1997), pp. 1295–1300. doi: 10.1126/science.278.5341.1295
- J. Forde, J.M. Volpe, and S.M. Ciupe, Latently infected cell activation: A way to reduce the size of the HIV reservoir? Bull. Math. Biol. 74 (2012), pp. 1651–1672. doi: 10.1007/s11538-012-9729-x
- Z. Grossman, M. Polis, M.B. Feinberg, Z. Grossman, I. Levi, S. Jankelevich, R. Yarchoan, J. Boon, F. de Wolf, J.M.A. Lange, J. Goudsmit, D.S. Dimitrov, and W.E. Paul, Ongoing HIV dissemination during HAART, Nat. Med. 5 (1999), pp. 1099–1104. doi: 10.1038/13410
- K. Gu, S. -I. Niculescu, and J. Chen, On stability crossing curves for general systems with two delays, J. Math. Anal. Appl. 311 (2005), pp. 231–253. doi: 10.1016/j.jmaa.2005.02.034
- J.K. Hale and S.M. Verduyn Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993.
- A.V. Herz, S. Bonhoeffer, R.M. Anderson, R.M. May, and M.A. Nowak, Viral dynamics in vivo: Limitations on estimates of intracellular delay and virus decay, Proc. Natl. Acad. Sci. USA. 93 (1996), pp. 7247–7251. doi: 10.1073/pnas.93.14.7247
- J. Li and Y. Kuang, Analysis of a model of the glucose-insulin regulatory system with two delays, SIAM J. Appl. Math. 67 (2007), pp. 757–776. doi: 10.1137/050634001
- B. Li and S. Liu, A delayed HIV-1 model with multiple target cells and general nonlinear incidence rate, J. Biol. Syst. 21 (2013). doi: 10.1142/S0218339013400123.
- D. Li and W. Ma, Asymptotic properties of a HIV-1 infection model with time delay, J. Math. Anal. Appl. 335 (2007), pp. 683–691. doi: 10.1016/j.jmaa.2007.02.006
- M.Y. Li and H. Shu, Global dynamics of an in-host viral model with intracellular delay, Bull. Math. Biol. 72 (2010), pp. 1492–1505. doi: 10.1007/s11538-010-9503-x
- M.Y. Li and H. Shu, Impact of intracellular delays and target-cell dynamics on in vivo viral infections, SIAM J. Appl. Math. 70 (2010), pp. 2434–2448. doi: 10.1137/090779322
- J. Li, Y. Kuang, and C.C. Mason, Modeling the glucose-insulin regulatory system and ultradian insulin secretory oscillations with two time delays, J. Theoret. Biol. 242 (2006), pp. 722–735. doi: 10.1016/j.jtbi.2006.04.002
- S. Liu and L. Wang, Global stability of an HIV-1 model with distributed intracellular delays and a combination therapy, Math. Biosci. Eng. 7 (2010), pp. 675–685. doi: 10.3934/mbe.2010.7.675
- M. Markowitz, M. Louie, A. Hurley, E. Sun, M. Di Mascio, A.S. Perelson, and D.D. Ho, A novel antiviral intervention results in more accurate assessment of human immunodeficiency virus type 1 replication dynamics and T-cell decay in vivo, J. Virol. 77 (2003), pp. 5037–5038. doi: 10.1128/JVI.77.8.5037-5038.2003
- J.E. Mittler, B. Sulzer, A.U. Neumann, and A.S. Perelson, Influence of delayed viral production on viral dynamics in HIV-1 infected patients, Math. Biosci. 152 (1998), pp. 143–163. doi: 10.1016/S0025-5564(98)10027-5
- P.W. Nelson and A.S. Perelson, Mathematical analysis of delay differential equation models of HIV-1 infection, Math. Biosci. 179 (2002), pp. 73–94. doi: 10.1016/S0025-5564(02)00099-8
- P.W. Nelson, J.D. Murray, and A.S. Perelson, A model of HIV-1 pathogenesis that includes an intracellular delay, Math. Biosci. 163 (2000), pp. 201–215. doi: 10.1016/S0025-5564(99)00055-3
- S. -I. Niculescu, P. Kim, K. Gu, P. Lee, and D. Levy, Stability crossing boundaries of delay systems modeling immune dynamics in leukemia, Discrete Contin. Dyn. Syst. Ser. B 13 (2009), pp. 129–156. doi: 10.3934/dcdsb.2010.13.129
- M.A. Nowak and R.M. May, Virus Dynamics: Mathematical Principles of Immunology and Virology, Oxford University Press, Oxford, 2000.
- S. Pankavich, The effects of latent infection on the dynamics of HIV, Differential Equations Dynam. Systems (2015). doi: 10.1007/s12591-014-0234-6.
- K.A. Pawelek, S. Liu, F. Pahlevani, and L. Rong, A model of HIV-1 infection with two time delays: Mathematical analysis and comparison with patient data, Math. Biosci. 235 (2012), pp. 98–109. doi: 10.1016/j.mbs.2011.11.002
- A.S. Perelson, Modelling viral and immune system dynamics, Nat. Rev. Immunol. 2 (2002), pp. 28–36. doi: 10.1038/nri700
- A.S. Perelson and P.W. Nelson, Mathematical analysis of HIV-1 dynamics in vivo, SIAM Rev. 41 (1999), pp. 3–44. doi: 10.1137/S0036144598335107
- A.S. Perelson and R.M. Ribeiro, Modeling the within-host dynamics of HIV infection, BMC Biol. 11 (2013), pp. 1–10. doi: 10.1186/1741-7007-11-96
- A.S. Perelson, D.E. Kirschner, and R.J. De Boer, Dynamics of HIV infection of CD4+ T cells, Math. Biosci. 114 (1993), pp. 81–125. doi: 10.1016/0025-5564(93)90043-A
- A.S. Perelson, P. Essunger, Y. Cao, M. Vesanen, A. Hurley, K. Saksela, M. Markowitz, and D.D. Ho, Decay characteristics of HIV-1-infected compartments during combination therapy, Nature 387 (1997), pp. 188–191. doi: 10.1038/387188a0
- B. Ramratnam, S. Bonhoeffer, J. Binley, A. Hurley, L. Zhang, J.E. Mittler, M. Markowitz, J.P. Moore, A.S. Perelson, and D.D. Ho, Rapid production and clearance of HIV-1 and hepatitis C virus assessed by large volume plasma apheresis, Lancet 354 (1999), pp. 1782–1785. doi: 10.1016/S0140-6736(99)02035-8
- B. Ramratnam, R. Ribeiro, T. He, C. Chung, V. Simon, J. Vanderhoeven, A. Hurley, L. Zhang, A.S. Perelson, D.D. Ho, and M. Markowitz, Intensification of antiretroviral therapy accelerates the decay of the HIV-1 latent reservoir and decreases, but does not eliminate, ongoing virus replication, J. Acquir. Immune Defic. Syndr. 35 (2004), pp. 33–37. doi: 10.1097/00126334-200401010-00004
- L. Rong and A.S. Perelson, Modeling Latently infected cell activation: Viral and latent reservoir persistence, and viral blips in HIV-infected patients on potent therapy, PLoS Comput. Biol. 5 (2009), p. e1000533. doi: 10.1371/journal.pcbi.1000533
- L. Rong and A.S. Perelson, Asymmetric division of activated latently infected cells may explain the decay kinetics of the HIV-1 latent reservoir and intermittent viral blips, Math. Biosci. 217 (2009), pp. 77–87. doi: 10.1016/j.mbs.2008.10.006
- L. Rong and A.S. Perelson, Modeling HIV persistence, the latent reservoir, and viral blips, J. Theoret. Biol. 260 (2009), pp. 308–331. doi: 10.1016/j.jtbi.2009.06.011
- L. Rong, Z. Feng, and A.S. Perelson, Emergence of HIV-1 drug resistance during antiretroviral treatment, Bull. Math. Biol. 69 (2007), pp. 2027–2060. doi: 10.1007/s11538-007-9203-3
- L. Rong, Z. Feng, and A.S. Perelson, Mathematical analysis of age-structured HIV-1 dynamics with combination antiretroviral therapy, SIAM J. Appl. Math. 67 (2007), pp. 731–756. doi: 10.1137/060663945
- S. Ruan and J. Wei, Periodic solutions of planar systems with two delays, Proc. R. Soc. Edinb. Ser. A 129 (1999), pp. 1017–1032. doi: 10.1017/S0308210500031061
- S. Ruan and J. Wei, On the zeros of transcendental functions with applications to stability of delay differential equations with two delays, Dyn. Contin. Discrete Impuls. Syst. 10 (2003), pp. 863–874.
- D. Sanchez-Taltavull and T. Alarcon, Stochastic modelling of viral blips in HIV-1-infected patients: Effects of inhomogeneous density fluctuations, J. Theoret. Biol. 371 (2015), pp. 79–89. doi: 10.1016/j.jtbi.2015.02.001
- X. Song, S. Wang, and J. Dong, Stability properties and Hopf bifurcation of a delayed viral infection model with lytic immune response, J. Math. Anal. Appl. 373 (2011), pp. 345–355. doi: 10.1016/j.jmaa.2010.04.010
- S. Wang and L. Rong, Stochastic population switch may explain the latent reservoir stability and intermittent viral blips in HIV patients on suppressive therapy, J. Theoret. Biol. 360 (2014), pp. 137–148. doi: 10.1016/j.jtbi.2014.06.042
- S. Wang, X. Song, and Z. Ge, Dynamics analysis of a delayed viral infection model with immune impairment, Appl. Math. Model. 35 (2011), pp. 4877–4885. doi: 10.1016/j.apm.2011.03.043
- X. Wang, Y. Tao, and X. Song, A delayed HIV-1 infection model with Beddington-DeAngelis functional response, Nonlinear Dyn. 62 (2010), pp. 67–72. doi: 10.1007/s11071-010-9699-1
- X. Wang, S. Liu, and L. Rong, Permanence and extinction of a non-autonomous HIV-1 model with time delays, Discrete Contin. Dyn. Syst. Ser. B 19 (2014), pp. 1783–1800. doi: 10.3934/dcdsb.2014.19.1783
- Y. Wang, Y. Zhou, J. Wu, and J. Heffernan, Oscillatory viral dynamics in a delayed HIV pathogenesis model, Math. Biosci. 219 (2009), pp. 104–112. doi: 10.1016/j.mbs.2009.03.003
- D. Wodarz, Killer Cell Dynamics: Mathematical and Computational Approaches to Immunology, Springer, New York, 2007.
- J.K. Wong, M. Hezareh, H.F. Günthard, D.V. Havlir, C.C. Ignacio, C.A. Spina, and D.D. Richman, Recovery of replication-competent HIV despite prolonged suppression of plasma viremia, Science 278 (1997), pp. 1291–1295. doi: 10.1126/science.278.5341.1291
- H. Zhu and X. Zou, Dynamics of a HIV-1 infection model with cell-mediated immune response and intracellular delay, Discrete Contin. Dyn. Syst. Ser. B 12 (2009), pp. 511–524. doi: 10.3934/dcdsb.2009.12.511