Dynamics of an HIV model with cytotoxic T-lymphocyte memory

We consider a four-dimensional HIV model that includes healthy cells, infected cells, primary cytotoxic T-lymphocyte response (CTLp), and secondary cytotoxic T-lymphocyte response (CTLe). The CTL memory generation depends on CD4+ T-cell help, and infection of CD4+ T cells results in impaired T-cell help. We show that the system has up to five equilibria. By the Routh–Hurwitz theorem and central manifold theorem we obtain some sufficient conditions for the local stability, globally stability of the equilibria, and the bifurcations. We still discover the bistability case where in the system there may coexist two stable equilibria or a stable equilibrium together with a stable limit cycle. Several numerical analyses are carried out to illustrate the validity of our theoretical results.


Introduction
Human immunodeficiency virus (HIV) and acquired immunodeficiency syndrome (AIDS) have become global health problems; there were estimated 38.0 million people living with HIV at the end of 2019 [1]. HIV can be transmitted via the exchange of a variety of body fluids from infected people, such as blood, breast milk, or semen and vaginal secretions, but it is known that current antiretroviral drugs cannot enucleate HIV from the body. Mathematical models have been formulated for various epidemiological diseases like novel coronavirus   [2][3][4] or malaria [5]. Models of HIV as a chronic infectious disease have been investigated in many papers. Some models focus on the size of infection, starting from the population models [6][7][8], whereas other focus on cell infection, starting from the virus models, which have become an important tool in both understanding HIV-1 infection in host and providing valuable insight into HIV pathogenesis. A well-known model for HIV infection is the following system [9,10]: where variables x, y, and v represent the densities of the healthy cells, the infected cells, and the virus at time t, respectively. Here a mass action infection mechanism is adopted. The parameters β and λ stand for the infection and constant growth rates of the healthy cell, respectively. System (1.1) and its variations have been investigated in many papers [11][12][13][14][15][16]. Stephen et al. [17] recently added the term ρTV C+V incorporating the homeostatic proliferation of T-cells, which leads to interesting dynamic results, such as bistability and Hopf bifurcation.
It is well-known that it takes a long period for an HIV to become an AIDS, and in the medical literature [18], it was pointed out that the latent reservoir (i.e., latent infection) was the main obstacle to eradicate the virus. Therefore the four-dimensional mathematical model including the latent infection seems more reasonable [19,20]. In recent years, latent cells were considered in many models, such as Beddington-DeAngelis function response with delay [20], Crowley-Martin function response [21], and general infection function with CTC and VTC transmission [22].
To recover from a viral infection, the cytotoxic T lymphocyte (CTL), which can clear away the infected cells to prevent further viral replications, plays a particularly important role. In 1996, Nowak and Bangham [9] proposed the well-known model with immune response: where the variable z represents the concentration of CTLs. Many authors have studied the infective models with different immune responses, such as lytic and nonlytic immune responses [23,24], cell-mediated immune mechanism or humoral immune mechanism [25][26][27], delayed immune response with drug therapies [28], and general CTL immune response with silent infected cell-to-cell spread [29].
However, after a viral infection, the CTLs that are responsible for clearing away the infected cells become cytotoxic T-lymphocyte precursors (CTLp) and have receptors for detecting the virus from the previous infection [30]. Upon contacting with the virus during a subsequent infection, the precursors differentiate and become cytotoxic T-lymphocyte effectors (CTLe), and these cells are again responsible for clearing away the invading virus. Considering this infective mechanism, Wodarz et al. [31,32] provided the following model with CTL response: where the healthy cells x and the infected cells y are described similarly as in system (1.1). Instead of just one class of CTL response, the CTLp and CTLe are introduced. The CTLp and CTLe are represented by w and z. These precursors emerge at rate cyw and may become effectors at rate cqyw or cleared away naturally at rate bw. Similarly, the effectors are created at rate cqyw and cleared at rate hz.
In model (1.2), there is no virus term whose population is assumed at a quasi-steady state, which is proportional to infected cells. Model (1.2) is completely analyzed by Bernard et al. [33], who have found that the system transforms from one equilibrium to the next as the basic reproductive number R 0 increases. When R 0 increases further, they show that periodic solutions may arise from the third equilibrium via Hopf bifurcation.
In fact, another model with cytotoxic T-lymphocyte memory proposed in [31,32] is given by This model assumes that the target cells are CD4 + T cells; moreover, it includes the additional feature that expansion of the CTLp population is proportional to both antigen and the number of uninfected CD4 + T-cells capable of delivering T-cell help. The memory generation depends on CD4 + T-cell help, and infection of CD4 + T-cells results in impaired T-cell help. We also assume that differentiation into effector functions is independent of CD4 + T-cell help [31]. A detailed explanation of the model can be found in [31]. All the parameters are positive.
Dynamics of system (1.3) is numerically analyzed in [31,32]. In this paper, we provide a rigorous analytical method of system (1.3), and the basic framework is as follows. In Sect. 2, we establish the well-posedness of the model including nonnegativity and boundedness of the solutions, the existence of equilibria, and local stability of the boundary equilibria. The local stability analysis of the positive equilibria and their bifurcations are presented in Sect. 3. Numerical illustrations are given in Sect. 4. Finally, we discuss both mathematical and biological perspectives of the findings in Sect. 5.

The equilibrium and stability of boundary equilibrium
For mathematical simplicity, we do some rescallings in system (1.3).
After changing back to the origin variables x, y, w, z, t, the scaled system is given by where the horizontal lines on the heads of these parameters are removed, and the parameters d, a, c, q, b, h are replaced by d, a, c, q, b, h. Obviously, all the parameters are positive. The basic reproductive number of model (1.3) is R 00 = λβ ad , and for system (2.1), R 00 becomes R 0 = 1 ad .
Proof By variation of constants we find the following solutions of (2.1): which proves the nonnegativity of solutions of system (2.1). Note that the first equation of (2.1) implies dx dt ≤ 1dx. The solution is given by Adding the first two equations of (2.1), we obtain where d 1 = min{d, a}. It has the solution x + y ≤ (x(0) + y(0))e -d 1 t + 1 d 1 , which implies lim sup t→∞ (x(t) + y(t)) ≤ 1 d 1 , and thus x(t) and y(t) are bounded. Supposing that z is unbounded, by the second equation of (2.1) we have lim t→∞ y(t) = 0, which implies lim t→∞ w(t) = 0 from the third equation of (2.1). Then we get lim t→∞ z(t) = 0 from the fourth equation of (2.1), which contradicts with the unboundedness of z. Thus z must be bounded. Lastly, assume that w is unbounded. Based on the boundedness of z and the fourth equation of (2.1), we obtain lim t→∞ y(t) = 0, and from the third equation of (2.1) it follows that lim t→∞ w(t) = 0, which causes a contradiction. Hence w is bounded. The proof is complete. Theorem 2.1 shows that there exists a bounded positive invariant region ⊂ R 4 + for the system. Thus we concentrate on to discuss the dynamics. In fact, the infectionfree equilibrium E 0 = ( 1 d , 0, 0, 0) always exists, and there exists an infectious equilibrium without CTL, E 1 = (a, d(R 0 -1), 0, 0) if R 0 > 1. To find the infectious equilibrium with CTL, it suffices to solve the system If x 1 , x 2 are solutions of (2.4), then x 1 x 2 = q cd > 0 and x 1 + x 2 = -b-c-dq cd . Hence, to ensure the existence of a positive solution, we have b < c + dq.
The determinant of (2.4) is Note that c + dq -2 cdq ≥ 0. Combining this with the condition b < c + dq, we get that ≥ 0 only if 0 < b ≤ c + dq -2 cdq, and y > 0 and z > 0 imply x > q c and x > a by (2.3). Set The positive equilibria of model (2.1) are classified by the sign of . Let us consider three cases.

Case I
> 0, which is equivalent to 0 < b < c + dq -2 cdq. Here we consider three cases by the sign of f (a).
1. f (a) < 0. In this case, > 0 is obviously satisfied, and a > q c is equivalent to R 1 > 1, and thus R 1 < R 0 by (2.6). Combining this with R 1 > 1, we get that model (2.1) admits a unique positive equilibrium 2. f (a) = 0. In this case, R 0 = R 1 , which indicates that a is a solution of f (x) = 0, and the other root is x * = q acd . To get a positive equilibrium of model (2.1), x * = q acd > a and a > q c are required, which yields q > a 2 cd and ac > q. Combining this with > 0 (if and only if 0 3. f (a) > 0. Let us consider two cases.
(1) q c ≥ a, which is equivalent to R 1 ≤ 1. To ensure the positive equilibrium of model (2.1), the following conditions are required: (2.7) Note that f ( q c ) > 0 always holds, and > 0 if and only if 0 Combining this with c > dq, we have ad < 1 to ensure that the intersection is nonempty, and hence model (2) q c < a, which is equivalent to R 1 > 1. Model (2.1) admits two positive equilibria if the following conditions are satisfied: (2.8) if and only if q < a 2 cd, and hence model (2.1) admits two positive equilibria if one of the following conditions holds:

Case II
= 0 if and only if b = c + dq -2 cdq, and model (2.1) admits a unique positive equilibrium if Combining these two conditions with = 0 which is equivalent to b = c + dq -2 cdq, the following results are obtained: Case III If < 0, then there is no positive equilibrium, because f (x) = 0 has no real root. The following theorem summarizes all positive equilibria of system (2.1).
Next, we discuss the local stability and global stability of the boundary equilibria. The stability of the equilibria is based on the Jacobian matrix of model (2.1): The characteristic roots at the equilibrium E 0 are given by and Hence the local stability of equilibrium E 0 is given by the following theorem.

Theorem 2.3 The infection-free equilibrium E
Proof To show the global stability of equilibrium E 0 , we use the method of fluctuation lemma employed by Hirsch et al. [34][35][36]. We introduce some notations. For a continuous bounded function g : [0, ∞] → R, let By Theorem 2.1 the solutions x(t), y(t), w(t), and z(t) are always nonnegative and bounded for any nonnegative initial conditions, and the limits lim sup t→∞ and lim inf t→∞ always exist for each of these solutions. By the fluctuation lemma there exists a sequence t n such that (2.11) By a similar argument, for the remaining three equations in model (2.1), we get (2.14) We claim that y ∞ = 0. Otherwise, y ∞ > 0; then it follows from (2.11) and (2.12) that Therefore R 0 ≥ 1, which contradicts the condition R 0 < 1, and thus y ∞ = 0. And from equations (2.13) and (2.14) we get that w ∞ = 0 and z ∞ = 0. The conditions y ∞ = 0, w ∞ = 0, and z ∞ = 0 imply that y(t) → 0, w(t) → 0, and z(t) → 0 as t → ∞, respectively. Thus from the first equation of model (2.1) we have the asymptotic differential equationẋ = 1dx. By simple calculation we get the solution d , which clearly shows that the solution x(t) → 1 d as t → ∞ by the theory of asymptotically autonomous systems. The proof is complete.
Letting (s + d + y 1 )(sx 1 + a) + x 1 y 1 = 0, we have s 2 + (d + y 1 )s + ay 1 = 0. Note that s 1 + s 2 = -(d + y 1 ) < 0 and s 1 s 2 = ay 1 > 0, which imply that s 1 and s 2 have negative real parts. Obviously, s 3 = -h has a negative real part. In order for all the roots of equation (2.15) to have negative real parts, it is required that Theorem 2.5 Equilibrium E 1 is locally asymptotically stable for 1 < R 0 < R 1 or R 0 > 1 > R 1 and is unstable for Remark Note that f (a) = ab (R 1 -1) (R 1 -R 0 ) = -as 4 , from which it follows that E 1 is locally asymptotically stable, whereas system (2.1) may have positive equilibrium under certain conditions by Theorem 2.2 (case I.3). Moreover, system (2.1) may have both stable positive equilibrium and stable boundary equilibrium E 1 .

Stability of positive equilibria and their bifurcations
By an easy calculation the characteristic equation of the positive equilibrium follows: It has a characteristic root s 1 = 0, and s 2 , s 3 , s 4 are determined by the equation where the horizontal lines on the heads of these letters are removed, and we still denote x, y, w, z by x, y, w, z. The third equation of model (3.7) has no linear term, and then the center manifold is a curve tangent to the w-axis.

Theorem 3.2 The infectious equilibrium E 2-is unstable once it exists.
Proof Since the determinant det(E 2-) < 0, there is at least one characteristic root that has no negative real part. Therefore E 2-is unstable.
Next, we discuss the stability of equilibria E 2+ . By (3.2) the characteristic equation at E 2+ is given by Note that α 1 > 0, α 2 > 0, and α 3 > 0, and according to the proof of Propositiony 3.1, we know that α 4 > 0. The relevant Routh-Hurwitz determinants are (3.12) Note that 1 > 0 and the sign of 4 is the same as 3 . By the formulas of x 2+ , y 2+ , and z 2+ , 2 and 3 can be written more explicitly as Next, we give the following lemma to show that if both 2 and 3 can become zero, then 3 will cross zero before 2 does. Remark The proof is similar to that in [33], so here we omit it.
Thus, to consider the stability of the positive equilibrium E 2+ , we only need to consider the possibility of 3 = 0 (see [33]). Note that x, y, z do not contain h and 3 = 0 is just a quadratic equation in terms of h. Let (3.14) Obviously, h * 1 < h * < h * 2 , and thus we have the following theorem.

, +∞) and unstable when h ∈
Remarks 1. The proof is straightforward by considering the sign of the quadratic polynomial A 3 h 2 + B 3 h + C 3 , and the detailed proof is contained in Theorem 3.4, so we omit it.
2. If we use the computing method of paper [33], then 2 and 3 can be written more explicitly as ).

a Hopf bifurcation occurs at
Remark Since the expressions of x, y, w, z are very complicated, it is too complicate to directly discuss the sign of 4 = B 2 3 -4A 3 C 3 . In the next section, by numerical calculations we will demonstrate that the three cases 4 > 0, 4 = 0, and 4 < 0 are possible.

Numerical illustrations
In this section, we demonstrate the theoretical results by numerical simulations. For convenience, we will work on the scaled model (2.1) instead of the original model (1.3).
The values of λ used for simulations in [32] are λ = 1 and λ = 10, and as a bifurcation in [33], this indicates that the values of λ have a fairly large variation; therefore we take λ = 0.5657. Note that β = 0.5 and β = 0.001 were used in [32], whereas β = 3 400 was used in [33], and hence choosing β = 0.0707 is reasonable. The parameter q ∈ [0, 1] represents the decomposition rate; here we choose q = 0.7071. The same parameter values as in [32] are taken, and appropriate value for h is chosen: According to the relationship between the parameters of the original model (1.3) and the scaled model (2.1), the parameter values in the scaled model (2.1) are chosen as follows: Choose a as the bifurcation parameter. With these parameter values, we obtain To get R 0 < 1, a > 2 is required, which in turn yields that the infection-free equilibrium E 0 = ( 1 d , 0, 0, 0) = (2, 0, 0, 0) is globally stable by Theorem 2.4. For a = 3, system (2.1) has a globally stable infection-free equilibrium E 0 , which is shown in Fig. 1.   As shown in Fig. 1, the healthy and infected cells decrease and converge to the equilibrium E 0 directly, whereas CTLe and CTLp firstly increase, then decrease, and finally converge to zero, which revels that the infected cells CTLp and CTLe die out directly after a brief fluctuation.
Next, we will prove Theorem 3.   it is locally stable, as shown in Fig. 3.
As can be seen in Figs. 4(a), 4(b), and 4(c), a stable limit cycle appears, and stable periodic solutions bifurcate from it.
As shown in Figs. 5(a) and 5(b), there are two stable equilibria E 1 and E 2+ , that is, bistability occurs, and the infected cells converge to one of them depending on the initial values.
Besides, when equilibrium E 2+ is unstable, then a stable limit cycle occurs, and E 1 is still stable. To display this phenomenon directly, letting a = 0.1 and h = 0.3, we draw the diagrams with different initial values.
As shown in Figs. 6(a) and 6(b), a stable equilibrium E 1 , an unstable equilibrium E 2+ , and a stable limit cycle, which is bifurcated from E 2+ , appear, and the infected cells converge to one of them depending on the initial values.
For the case of C 3 < 0, in Theorems 3.3 and 3.4, we give another set of parameters d = 3.6 = b, q = 1, c = 28.8, h = 0.8. We do not discuss it now.

Discussion
This paper studies an HIV model proposed by Wodarz et al. [31,32] to describe the interaction between healthy cells and infected cells as well as primary and secondary immune   response. Compared with [33], in this model, we assume that the production of primary immune response is not only connected with infected cells but also with healthy cells. We also assume that virus at its steady state is proportional to infected cells. The structure of equilibria is analyzed in [31,32]. But for a higher-dimensional system, stability and bifurcation analysis is important and complex for the full range of possibilities. Because of adding the healthy cells to the produced CTLp, the model shows rich dynamic behavior on stability and bifurcations.  It is interesting that this model displays the bistability phenomenon, that is, two stable equilibria E 1 and E 2+ or a stable equilibrium E 1 and a stable limit cycle, which is bifurcated from the unstable equilibrium E 2+ . Which one is stable not only depends on relationship of parameters but also depends on the initial values of cells. As shown in Figs. 5(a) and 5(b), high initial virus load close to E 1 leads to the convergence to E 1 , which means that CTL memory fails to establish. Low initial virus load close to E 2+ leads to the convergence to E 2+ , CTL memory successfully establishes, and the virus load first increases and then decreases to stay at a low level. Therefore we must increase dosage to inhibit the virus replication in a brief period and help the immune response establish. As shown in Fig. 6(a), if the initial healthy cells and virus load are the same, and high initial CTL account means that the CTL response establishes, virus load may decrease at primary process, which was considered as the CTL clear away some virus, but over time, it oscillates and cannot be completely eradicated. This phenomenon can be viewed as an individual having a chronic disease that may flare up from time to time, and it is a long struggle between virus and immune response. Initial values with little CTL lead to high virus load. As can been in Fig. 6(c), high or low initial healthy cells do not change the development of disease. However, high initial healthy cells first lead to decreasing the virus load to a very low state, which is almost clear away the virus, which means that the healthy cells play an important role in clearing away the virus under certain circumstances.
The interaction of virus and the host immune system is a complicated and long process as HIV has a long latent period, and the disease cannot cure completely. We have shown rich dynamic patterns, but the model considered here is just a simple one. It is easy to improve and expand the model. For example, we can add the virus equation in model (1.3) (see a more detailed description in [32]), we may consider delay in this model as in [39], drug treatment [13], or latent cells [40] mentioned in the introduction. Such modifications should more precisely react the reality and give us more advice in understanding the infection process, which leads to a more challenging mathematical analysis.