A fit of CD4⁺ T cell immune response to an infection by lymphocytic choriomeningitis virus

1.   Introduction

Human immunodeficiency virus (HIV) takes the most important CD4$^{+}$T lymphocytes in the human immune system as the main target, destroys people's CD4$^{+}$T and makes the body lose its immune function ^[1]. HIV is a major problem facing the human race and poses a serious health threat to human society. While there has been remarkable advancement in the development of antiretroviral therapy (ART) and prevention strategies, currently there are still many people living with HIV ^[2]. From the National Health Commission of the People's Republic of China, to the end of 2020, a total of 1.053 million people in China were infected with HIV, 351,000 deaths were reported, and the number of infections is expected to increase to 1.6 million in 2022. Therefore, it is beneficial to study the effective control measures and find the optimal control by using mathematical models.

Early HIV models were primarily devoted to the virus-to-cell infection ^{[3,4,5,6,7,8,9,10]}. With the development of science and the improvement of medical standards, some studies show that virus can also spread by direct cell-to-cell transmission ^{[11,12,13,14,15]}. Thus, many researchers have begun working on an HIV model incorporating virus-to-cell and cell-to-cell transmission ^{[16,17,18,19,20]}. For example, by an HIV model, Wang et al. ^[16] showed the existence, positivity and boundedness of the model solution. Lai et al. ^[19] demonstrated global threshold dynamics by the basic reproduction number. For the control of acquired immune deficiency syndrome (AIDS), the main drug treatment is to prevent new HIV infections by blocking the transformation of viral RNA into DNA in T cells and reducing the number of viral particles ^[21]. With regard to AIDS control, there have already been some results ^[22,23,24]. Liu et al. ^[22], Akbari et al.^[23] and Guo et al. ^[24] proposed an optimal control problem for an HIV infection model with cell-to-cell spread. However, these studies were focused on optimal control given over the entire infection period $T$, which will increase the cost of control measures. If threshold level $N_t$ can be introduced, such that the control measures are not implemented when the total number of infected cells in an infected person is at a relatively low level (threshold level $N_t$), while otherwise effective measures are taken to suppress progression of viruses and infected cells, the time and cost will be reduced. The idea has been widely used in engineering ^[25,26], but little research has been done in epidemic diseases.

In addition, for the case that the number of infected cells in the body exceeds the threshold level ($N_t$), from the perspectives of epidemiology and economics, how to control the spread of virus is a valuable question. As indicated by discussion above, in this paper, we analyze the sliding mode dynamics and optimal control of the HIV model with virus-to-cell and cell-to-cell transmission. Here, because the timing of controlling infected cells and viruses is uncertain, dynamic programming is considered, which demonstrates that not only can infected cells and viruses be controlled in time, but also the goal of minimizing the concentration of infected cells and viruses with a low cost of application control can be achieved. The main novelties are summarized as follows:

● About the HIV model with cell-to-cell transmission, the majority of the existing results only discussed dynamics, this article proposes a piecewise control function concerning threshold policy and discusses sliding mode dynamics.

● Differently from previous works on optimal control over the time period $[0, T]$ for the HIV model, in this paper, optimal control strategies for infected cells and viruses are achieved when the total number of infected cells in the body exceeds the certain tolerance threshold level ($N_t$). Moreover, our results generalize and improve some published results in the literature, such as ^[23,24].

The whole organization of this work is as follows. Section 2 describes the different components of the HIV model, and it then further extends a new three-dimensional Filippov model with two control measures. Section 3 investigates the sliding mode dynamics of the model and shows the existence of a unique positive pseudo-equilibrium. In Section 4, the optimal control problem is discussed. We first give the objective function and prove the uniqueness and existence of the viscosity solution of the HJB equation, and we then obtain the optimal control through the Hamiltonian function. The theoretical results are verified by numerical simulations in Section 5. Finally, conclusions and outlook for further work are given in Section 6.

2.   Model and preliminary knowledge

Inspired by ^[17,19], the classical HIV model with virus-to-cell and cell-to-cell transmission is

where $x(t)$ denotes the concentration of uninfected target cells at time $t$, $y(t)$ is the concentration of infected cells at time $t$, and $v(t)$ denotes the concentration of virus particles at time $t$. $\Lambda$ is the recruitment rate of healthy target cells, $\beta$ is the rate at which an uninfected cell becomes infected by an infectious virus, $\gamma$ represents the infection rate of productively infected cells, $k$ shows the generation rate of virus particles, $a$ expresses the loss rate of infected cells, $b$ represents the natural death rate of uninfected cells, and $c$ indicates the clearance rate of virions.

Based on the existing HIV models, we give the sliding control system to maintain the number of viruses and infected cells below the threshold level. Research ^[27] has found that current drug treatment consisting of five antiretroviral drugs can suppress viral replication to a low level or increase the CD4+T cell, the two main types of HIV resistance: reverse transcriptase inhibitors (RTIs) and protease inhibitors (PIs) ^[28,29]. RTIs prevent new HIV infections by blocking the transformation of viral RNA into DNA in T cells, and PIs reduce the number of viral particles produced by actively infected T cells ^[21]. We represent by $\mu_1$ the RTIs control variable and by $\mu_2$ the PIs control variable. The control system is given as follows:

with

$\mu_1$ is the culling rate of infected cells, $\mu_2$ is the clearance rate of virus particles, and they are constants. The critical level of the total number of infected cells is represented by $N_t$. For convenience, $y(t)-N_t$ is defined as $\omega(M) = y(t)-N_t$ with vector $M = (x(t), y(t), v(t))^{\top}\in \mathbb{R}_+^3$ and $\mathbb{R}_+^3 = \{M = (x, y, v)|x \geq 0, y \geq 0, v \geq 0\}$.

Remark 2.1. When $\varepsilon = 0$, model (2.2) becomes model (2.1) and the control measures are taken when $\varepsilon = 1$.

3.   The sliding mode dynamics of HIV model

In this section, according to system (2.2) with (2.3), we study the dynamics of the system. First, we define $G_1 = \{M\in \mathbb{R}_+^3|\omega(M) < 0\}$, $G_2 = \{M\in \mathbb{R}_+^3|\omega(M) > 0\}$. Furthermore, we describe the manifold $G_s$ as $G_s = \{M\in \mathbb{R}_+^3|\omega(M) = 0\}$ and the normal vector perpendicular to $G_s$ is shown as $n = (0, 1, 0)^{\top}$. Then, we consider the following Filippov system:

where

Based on this, we present the following definitions of various equilibriums and sliding domain.

Definition 3.1. ^[30] If $f_1(M^*) = 0$ with $\omega(M^*) < 0$, or $f_2(M^*) = 0$ with $\omega(M^*) > 0$, the point $M^*$ is called a real equilibrium of system (3.1).

Definition 3.2. ^[30] If $f_1(M^*) = 0$ with $\omega(M^*) > 0$, or $f_2(M^*) = 0$ with $\omega(M^*) < 0$, the point $M^*$ is called a virtual equilibrium of system (3.1).

Definition 3.3. ^[30] If it is an equilibrium of the sliding mode of system (3.1), the point $M^*$ is called a pseudo-equilibrium.

Definition 3.4. ^[30] $S$ is the sliding domain, if $\langle n, f_1\rangle > 0$ and $\langle n, f_2\rangle < 0$ on $S\subset G_s$.

3.1.   Analysis in region $G_1$

In this subsection, we calculate the basic reproduction number and analyze the stability of equilibrium. The dynamics of system (3.1) in region $G_1$ are indicated by

In subsystem (3.2), the disease-free equilibrium $E_1^0 = \left(\frac{\Lambda}{a}, 0, 0\right)$ and the endemic equilibrium $E_1^* = (x_1^*, y_1^*, v_1^*)$ are given, where

Next, the following theorem about the local asymptotic stability of equilibria $E_1^0$ and $E_1^*$ are given.

Theorem 3.5. The disease-free equilibrium $E_1^0$ is locally asymptotically stable if $R_{01} < 1$; the endemic equilibrium $E_1^*$ exists and is locally asymptotically stable if $R_{01} > 1$.

Proof. Using the next generation matrix, the basic reproduction number $R_{01}$ of system (3.2) can be deduced, and system (3.2) can be represented as

where

The Jacobian matrix of $F(M)$ and $V(M)$ at the equilibrium point $E_{1}^{0}$ is

where

Then, the generation matrix of system (3.2) is

So, the spectral radius of $FV^{-1}$ is

where $\rho$ is the spectral radius of a matrix. Therefore, the basic reproduction number of system (3.2) is

For subsystem (3.2), The Jacobian matrix is

If $R_{01} < 1$, $c b-\frac{c\gamma \Lambda+k \beta \Lambda}{a} > 0$, $c+b-\frac{\gamma\Lambda}{a} > 0$ can be obtained, and $(c+b-\frac{\gamma\Lambda}{a})^2-4(c b-\frac{c\gamma \Lambda+k \beta \Lambda}{a}) > 0$. The characteristic equation at $E_1^0$ is

which indicates that

Thus, the eigenvalues of $J(E_1^0)$ are negative when $R_{01} < 1$, so we can obtain that $E_1^0$ is locally asymptotically stable.

If $R_{01} > 1$, then $\frac{\gamma c \Lambda+k\beta \Lambda}{acb} > 1$, and there exists the endemic equilibrium $E_1^*$. Thus, the characteristic equation at $E_1^*$ is

where

Moreover, we can get

According to the Routh-Hurwitz Criterion, all eigenvalues of $J(E_1^*)$ have negative real parts. Hence, this represents that $E_1^*$ is locally asymptotic stable.

3.2.   Analysis in region $G_2$

In region $G_2$, we give the dynamics of system (3.1):

In subsystem (3.3), the disease-free equilibrium $E_2^0 = \left(\frac{\Lambda}{a}, 0, 0\right)$ and the endemic equilibrium $E_2^* = (x_2^*, y_2^*, v_2^*)$ are given, where

Next, the following theorem about the local asymptotic stability of equilibria $E_2^0$ and $E_2^*$ are given.

Theorem 3.6. The disease-free equilibrium $E_2^0$ is locally asymptotically stable if $R_{02} < 1$. The endemic equilibrium $E_2^*$ exists and is locally asymptotically stable if $R_{02} > 1$.

Proof. Using the next generation matrix, we deduce the basic reproduction number $R_{02}$ of system (3.3):

For subsystem (3.3), the Jacobian matrix is

Let $c+\mu_2 = c_1$, $b+\mu_1 = b_1$. Similar to the proof of Theorem 3.5, if $R_{02} < 1$, then all eigenvalues of $J(E_2^0)$ are negative, and thus the local asymptotic stability of $E_2^0$ can also be concluded. Furthermore, all eigenvalues of $J(E_2^*)$ have negative real parts, if $R_{02} > 1$, which shows that $E_2^*$ is locally asymptotically stable.

Remark 3.7. The basic reproduction number $R_{02}$ is related to control variables $\mu_1$ and $\mu_2$, and when $\mu_1 = \mu_2 = 0$, that is, the control measures are not implemented, we have $R_{01} = R_{02}$.

3.3.   Analysis in region $G_s$

In order to study the dynamics of sliding mode of system (3.1) in this subsection, we initially examine the existence of the sliding mode. The manifold $G_s$ is defined as $y(t) = N_t$, and we have

and

which shows that $\sigma_2(x(t), v(t)) < \sigma_1(x(t), v(t))$. Because $\langle n, f_1\rangle > 0$ and $\langle n, f_2\rangle < 0$, we can obtain

and

It is obvious that $H_1(v(t))$ is positive on the interval $(0, N_t)$, and so is $H_2(v(t))$. Then, the sliding domain $S\subset G_s$ is defined as follows:

Moreover, in order to obtain sliding mode equations in the region $S$, we use the Utkin equivalent control method, as used in ^[31]. First, according to $\omega(y(t)) = y(t)-N_t = 0$ and $d\omega/dt = 0$,

By solving Eq (3.7), the function about $\varepsilon$ is obtained:

Substitute Eq (3.8) into system (2.2), and the dynamics of sliding mode in $S$ can be given by differential equations as follows:

Next, we exhibit the following proposition about the existence of equilibrium of sliding mode (3.9).

Proposition 3.8. The existence of a unique positive pseudo-equilibrium $E_s^* = (\frac{\Lambda}{a+\beta v^*+\gamma N_t}, N_t, v^*)$ of sliding mode (3.9) is obtained if $\frac{c-k}{b} < \frac{\mu_2}{\mu_1} < \frac{c}{b}$ and $E_s^*\in S$.

Proof. Let $\frac{dx(t)}{dt} = 0$ and $\frac{dv(t)}{dt} = 0$. It obtains that

and

After substituting Eq (3.10) into (3.11), we get the following equation of $v(t)$:

where

If $\frac{c-k}{b} < \frac{\mu_2}{\mu_1} < \frac{c}{b}$, then $\tau_1 < 0, \tau_3 > 0$ and $\Delta > 0$ can be received. The Vieta Theorem shows

Due to $\Delta > 0$, the existence of the root is proven, and there are two roots, one positive root and one negative root from (3.13). We can conclude that Eq (3.12) has a unique positive root denoted by $v^*$. Thus, the equilibrium point $E_s^*$ can be obtained from (3.10):

Furthermore, $E_s^*$ is a unique pseudo-equilibrium if $E_s^*\in S\subset G_s$ holds.

Similar to the previous theorem, it is worth considering the local asymptotical stability of $E_s^*$, and then the following theorem is given.

Theorem 3.9. Under the same conditions of Proposition 3.8, $E_s^*$ is locally asymptotically stable on the sliding domain $S$.

Proof. In system (3.9), the Jacobian matrix of the first two equations is

From Proposition 3.8, $\frac{c-k}{b} < \frac{\mu_2}{\mu_1} < \frac{c}{b}$, we set

and obtain

Thus, we get

since $J_{11}+J_{22} < 0$, and we can obtain that

Therefore, all eigenvalues of (3.14) have negative real parts. That $E_s^*$ is locally asymptotically stable can be obtained easily.

Remark 3.10. The sliding domain $S$ represents the region between the two lines associated with $x(t), v(t)$ belonging to the plane $y = N_t$, and the local asymptotical stability of pseudo-equilibrium $E_{s}^{*}$ in sliding domain $S$ is related to the control variables, which needs to satisfy the condition of the parameters.

4.   Optimal control problem

In order to reduce the number of infected cells and viruses, while keeping the cost to apply the control at the minimum level at any time, we show an optimal control problem on the basis of system (2.2). In this section, any time refers to when the number of infected cells is greater than threshold level $N_t$, that is, when $\varepsilon = 1$ in system (2.2). If the control variable in system (2.2) is a time dependent variable, how do we find the optimal control? Thus, we represent by $u_1(t)$ the RTIs control variable and by $u_2(t)$ the PIs control variable. First, we have $u(t) = (u_1(t), u_2(t))^\top \in \mathcal{V}[s, T] = \{u(\cdot):[s, T]\rightarrow U|\; u_1(t)$$\; and\; u_2(t)\; are \; measurable: 0 \leq u_1(t) \leq 1, 0 \leq u_2(t)\leq 1\}$, where $U$ is a metric space and convex, $T > 0$, and the control problem of the model is given by

Let $(s, M)\in[0, T)\times\mathbb{R}^3$, and we consider the following control system over $[s, T]$:

where $M(t) = (x(t), y(t), v(t))^{\top}\in\mathbb{R}^3$. Next, we construct the following objective functional:

where $A_{i} \; (i = 1, 2, 3, 4)$ are weights to make the terms of the integrand keep balance. The term $\int_{s}^{T}(A_1 u_1^2(t)+A_2 u_2^2(t))dt$ gives the total cost of using the control strategy, $h(M(T))$ is the penalty function corresponding to the terminal state, and $\int_{s}^{T}(A_3 y(t)+A_4 v(t))dt$ represents the total number of infected cells and viruses over the time period $T$. It is important to find an optimal control pair $(u_1^*(t), u_2^*(t)), t\in[s, T]$ such that

Further, the value function is as follows:

Before further study, we give the following assumption.

Assumption 1. $(U, \tilde{d})$ is a separable metric space.

Then, we would like to study $V(\cdot, \cdot)$ in great detail and present the following results called Bellman's principle of optimality by ^[32].

Theorem 4.1. Let 1 hold, and $U$ is convex. Then, for any $(s, M_0)\in[0, T)\times\mathbb{R}^3$, we have

Proof. Let us define

By (4.6), we obtain

Therefore, we take the infimum over $u(t)\in \mathcal{V}[s, T]$ and get

and there exists a $u_{\varepsilon}(t)\in \mathcal{V}[s, T]$ for $\forall \varepsilon > 0$ such that

where $M_{\varepsilon}(\cdot) = M(\cdot; s, M_0, u_{\varepsilon}(\cdot))$. Combining (4.9) and (4.10), (4.6) is obvious. This completes the proof.

(4.6) is called the dynamic programming equation, and the equation is difficult to handle. Thus, we construct the Hamilton-Jacobi-Bellman (HJB) equation as follows.

Theorem 4.2. Let 1 hold, and $U$ is convex. Suppose $V\in C^{1}([0, T]\times \mathbb{R}^{3})$. Then, $V(s, M_0)$ is a solution to the following terminal value problem of a first-order partial differential equation:

where

We call (4.11) the Hamilton-Jacobi-Bellman (HJB) equation associated with (4.5).

Proof. Fix a $u\in U$. Let $M(t)$ be the state trajectory corresponding to the control $u(t)\equiv u$. According to (4.6) with $\hat{s}\downarrow s$, we get

which results in

For any $\varepsilon > 0, \; 0\leq s\leq \hat{s}\leq T$ with $\hat{s}-s > 0$ small enough, there exists a $u = u_{\varepsilon, \hat{s}}(\cdot)\in\mathcal{V}[s, T]$ such that

so it follows that

Combining (4.14) and (4.16), we complete the proof.

For further study, the definition of viscosity solution of (4.11) by ^[32] is given.

Definition 4.3. A function $V\in C[0, T]\times\mathbb{R}^{3})$ is called a viscosity subsolution (or supersolution) of (4.11) if

and for any $\varphi\in C^{1}([0, T]\times\mathbb{R})$, whenever $V-\varphi$ attains a local maximum(or minimum) at $(t, M)\in[0, T]\times\mathbb{R}^3$, we have

In this case that $V$ is both a viscosity subsolution and supersolution of (4.11), it is a viscosity solution of (4.11).

If we get the value function $V$ by solving the HJB equation, then an optimal pair could be constructed. Thus, the characterization of the value function of the viscosity solution to the HJB equation is given below.

Theorem 4.4. Let 1 hold, and $U$ is convex. Then, the value function $V(\cdot, \cdot)$ satisfies

for some $K > 0$. Moreover, $V$ is the only viscosity solution of (4.11) in the class $C([0, T]\times\mathbb{R}^{3})$.

Proof. We have

where $K = \max\{A_1 u_1^2(t)+A_2 u_2^2(t)+A_3 y(t)+A_4 v(t)\} > 0$. Then, similar to [32, Theorem 2.5], the value function $V$ clearly satisfies the condition and is the only viscosity solution of (4.11).

Next, in order to discuss the existence of the optimal control pair, we define $H(t, M, u(t), p)$ as the Hamiltonian as follows:

Theorem 4.5. There exists an optimal control pair $\left(u_{1}^{*}(t), u_{2}^{*}(t)\right)$ and a corresponding optimal state $\left(x^{*}(t), y^{*}(t), v^{*}(t)\right)$ such that

Proof. Note that both the control variables and state variables are nonnegative, and the objective functional (4.3) is convex with respect to the control variables. The convexity and closure of the control set $U$ are obtained according to the definition of it. Furthermore, the optimal control is bounded. Therefore, we conclude that there exists an optimal control $u_{1}^{*}(t), u_{2}^{*}(t)$ for $J\left(u_{1}^{*}(t), u_{2}^{*}(t)\right) = \min _{u_{1}(t), u_{2}(t)\in U} J\left(u_{1}(t), u_{2}(t)\right)$.

The HJB equation and the value function have been discussed. Further, we figure out the optimal control as follows.

Theorem 4.6. Let $u_{1}^{*}(t), u_{2}^{*}(t)$ be optimal control variables, and $x^{*}(t), y^{*}(t), v^{*}(t)$ are corresponding optimal state variables. There exists adjoint process $p(t) = (p_{1}(t), p_{2}(t), p_{3}(t))^\top$, satisfying the following adjoint equation:

Furthermore, the optimal control is given as follows:

where $B_{1} = \frac{y^{*}(t) p_{2}(t)}{2 A_{1}}, B_{2} = \frac{v^{*}(t) p_{3}(t)}{2 A_{2}}$.

Proof. Because the Hamiltonian function $H$ is given. Moreover, by using the optimal condition, and we obtain $u_{1}^{*}(t)$ and $u_{2}^{*}(t)$

Hence,

So, $B_{1} = \frac{y^{*}(t) p_{2}(t)}{2 A_{1}}, B_{2} = \frac{v^{*}(t) p_{3}(t)}{2 A_{2}}$ are taken, and the optimal control $u_{1}^{*}(t), u_{2}^{*}(t)$ follows:

Thus, the optimal value can be obtained.

Remark 4.7. In practical problems, when the parameters are known, the optimal control $u_{1}^{*}(t)$ and $u_{2}^{*}(t)$ can be calculated by using computer programming for Eqs (4.1), (4.23) and (4.27), that is, the dosage intensity of RTIs and PIs at each time can be calculated. If patients are treated in such a proportion, an optimal control strategy to minimize the costs and the number of viruses and infected cells can be obtained.

5.   Numerical simulations

In this section, several numerical simulations are presented to prove and verify the above theoretical results by MATLAB software.

5.1.   Numerical simulations of sliding mode

In this subsection, numerical simulations are performed to illustrate the theoretical results. We discuss the dynamics behavior of system (3.1). Since the choice of the threshold level $N_t$ is different, system (3.1) will show various dynamics. There are three cases to consider. First, we discretize system (3.1) as follows:

In Table 1, all parameter values are exhibited. Some values are collected from different papers ^[6,17], and others are assumed.

Case 1: $E_1^*$ is a real equilibrium, and $E_2^*$ is a virtual equilibrium.

This situation is established. Suppose that the following conditions are satisfied:

In this case, both equilibriums are located on the same side of the plane $N_t$, namely, the region $G_1$. Thus, when (5.2) is satisfied, we obtain that $E_1^*$ can achieve stability. In Figure 1 $E_1^*\in G_1$ achieves stability with $N_t = 60$, and the possible trajectories of this figure are shown: A trajectory that starts in region $G_1$ will converge to $E_1^*$ as $t\rightarrow +\infty$ with hitting and sliding down on the sliding domain $S\subset G_s$; a trajectory of the initial point in region $G_2$ will cross the manifold $G_s$, then enter region $G_1$ and finally converge to $E_1^*$; a trajectory that begins in region $G_1$ or $G_2$ will hit and slide to the boundary of the sliding domain $S\subset G_s$ before moving towards $E_1^*$.

Case 2: $E_2^*$ is a real equilibrium, and $E_1^*$ is a virtual equilibrium.

This situation is established. Suppose that the following conditions are satisfied:

In this case, both equilibriums are located on the same side of the plane $N_t$, namely, the region $G_2$. Thus, when (5.3) is satisfied, we obtain that $E_2^*$ can achieve stability. In Figure 2 $E_2^*\in G_2$ achieves stability with $N_t = 15$, and the possible trajectories of this figure are shown: A trajectory with initial condition in region $G_2$ will converge to $E_2^*$ as $t\rightarrow +\infty$ with hitting and sliding up on the sliding domain $S\subset G_s$; a trajectory with initial point in region $G_1$ will pass through the manifold $G_s$ from $G_1$ to $G_2$; a trajectory that starts in region $G_1$ or $G_2$ will hit and slide up on the boundary of the sliding domain $S\subset G_s$ before moving towards $E_2^*$.

Case 3: Both $E_1^*$ and $E_2^*$ are virtual equilibria.

This situation is established. Suppose that the following conditions are satisfied:

For this case, we can find that $E_1^*$ and $E_2^*$ belong to regions $G_2$ and $G_1$, respectively. If $E_S^*\in S\subset G_s$, it is a pseudo-equilibrium. Thus, when (5.4) is satisfied, we obtain that $E_s^*$ is stable when it exists. Figure 3 $E_s^*\in S$ achieves stability with $N_t = 25$, and the possible trajectories of this figure are shown: A trajectory with initial condition in region $G_1$ will converge to $E_s^*$ as $t\rightarrow +\infty$; a trajectory that begins in region $G_2$ will cross the manifold $G_s$, then enter region $G_1$ and finally converge to $E_s^*$; a trajectory with initial condition in region $G_1$ or $G_2$ will hit and slide down on $S\subset G_s$ before converging to $E_s^*$.

5.2.   Numerical simulations of the optimal control

In this subsection, we also discretize the optimal control system (4.1) as

The parameter values $\Lambda$, $a$, $b$, $\beta$, $\gamma$, $k$ and $c$ are chosen as shown in Table 1, and $A_1 = 20, A_2 = 50, A_3 = 0.1, A_4 = 0.1$. Then, we compare the effects of different control intensities by the following figures. The cost of control strategies must be considered, and the cost of each measure is different, so we want to know the change under only one control strategy compared with two control strategies. Thus, we make simulations under only $u_1(t)$ (see Figure 4) or $u_2(t)$ (see Figure 5) and optimal control (see Figure 6).

When initial value $(x_0, y_0, v_0) = (50,150,800)$, Figure 4 illustrates that if the control strategy is only applied to infected cells ($u_2(t) = 0$), the concentration of infected cells and viruses will decrease. Meanwhile, the number of infected cells and viruses decreased with the increase of $u_1(t)$ intensity. Figure 5 shows that if the control strategy is only applied to viruses ($u_1(t) = 0$), the density of infected cells and viruses will also decrease, and the number of infected cells and viruses decreased with the increase of $u_2(t)$ intensity.

The expressions of optimal control $u_1^{*}(t)$ and $u_2^{*}(t)$ are obtained through calculation (in Eq (4.27)), and then we obtain the optimal control states of uninfected cells, infected cells and viruses in Figure 6, which shows that the concentration of infected cells and viruses decreases to some extent after the control is applied, and the combination of multi-drug works better than a single-drug approach. Figure 7 shows the values of control variables $u_1(t)$ and $u_2(t)$ in each time. Eventually, the concentration of infected cells and viruses and control intensity gradually decrease and stabilize. Thus, the control strategies in our model have significant influence on the spread of HIV.

6.   Conclusions and discussion

In this paper, we have proposed an HIV model with cell-to-cell transmission and analyzed the sliding mode dynamics of the model and optimal control problem. First, we extended a novel Filippov model (3.1), which indicates that corresponding control measures (i.e., antiretroviral drugs RTIs and PIs) are triggered once the total number of infected cells reaches the threshold level $N_t$, where RTIs prevent new HIV infections by blocking the transformation of viral RNA into DNA in T cells, and PIs reduce the number of viral particles produced by actively infected T cells ^[21]. Further, the sliding domain and sliding mode dynamics of system (2.2) have been examined. In addition, the simulation results show that the model solution is either near a real equilibrium point or near a pseudo-equilibrium point according to the different threshold levels. It is worth mentioning that some parameters can be simulated through actual data, and the basic reproduction number can be calculated to judge the stability. In addition, the number of viruses and infected cells can reduce to a previously desired level when the threshold level is chosen properly. Because different patients have different initial viral loads, an individualized therapy is suggested, which shows that the choice of a treatment strategy for a given patient should depend on HIV viruses and infected cells and proposed threshold level.

Moreover, the cost of treatment is beyond the reach of many infected patients. Therefore, we have introduced an optimal therapy to minimize the cost of treatment and reduce the viral load and the number of infected cells. Then, the efficacies of RTIs and PIs and their combinations have been measured. In addition, we have discussed an efficient numerical method based on optimal control to determine the best treatment strategy for HIV infection. Our results indicate that with the increase of treatment intensity, the number of infected cells and viruses decreases, while the density of CD4$^{+}$T increases. Due to the multiple transmission routes of HIV, the combined use of multiple drugs is better than the use of a single drug. From a biological point of view, it can be concluded that optimal control is adopted when the number of infected cells in the patient is higher than threshold level $N_t$, and at this time, the economic cost can be considered to select the optimal control measures. Control measures are not necessary when the concentration of infected cells is low. Optimal control path is shown in Figure 7.

The results of this paper have practical implications for controlling HIV transmission. However, our work is only a preliminary exploration of the impact of some control measures on HIV transmission, which can be improved in many aspects. The number of viruses in the infected person varies with age, so it is interesting and challenging to consider the sliding mode control of HIV model with age structure. We can also refer to the methods in the literature ^[33,34] to study the fractional-order HIV model. These issues will be the focus of our future research.

Acknowledgments

The research was supported by the Key Research and Development Program of Ningxia (2020BEB04007) and the Fundamental Research Funds for the Central Universities, North Minzu University (2020KYQD17).

Conflict of interest

The authors declare there is no conflict of interest.