Dynamic analysis and optimal control of a fractional order HIV/HTLV co-infection model with HIV-specific antibody immune response

1.   Introduction

Acquired Immune Deficiency Syndrome (AIDS) was firstly identified in the United States in 1980s, and the corresponding virus was named Human Immunodeficiency Virus (HIV) in 1986 ^[1]. HIV is a chronic virus infecting the human immune system, which can cause immunodeficiency but not cancer, and it is one of the retroviruses. According to the report of World Health Organization, by 2022, over 37.7 million individuals have been infected with HIV ^[2]. The persistence of this disease has attracted the attention of many scholars in the field of biomathematics. From a mathematical point of view, they constructed some dynamic models of HIV to simulate the progression of the disease, and to understand the control effects of drugs and preventive measures. In 1989, Perelson put forward a mathematical model to analyze the interaction between HIV and the immune system ^[3]. The author in ^[4] introduced the idea of fractional order into HIV modeling, and they established a fractional order HIV model with drug control measures in 2009. In recent years, scholars have used different methods to conduct multi-scale research on HIV, including ordinary differential equations ^[5,6,7,8], stochastic differential equations ^{[9,10,11,12,13,14]}, fractional differential equations ^{[15,16,17,18]}, etc.

Human T-lymphocytic virus (HTLV), which is also a kind of retrovirus, has not spread as widely as HIV. It is mainly prevalent in Caribbean, Central Africa, Japan and South America. People who die from complications caused by HILV infection each year is up to 2 million ^[19]. HIV was once considered as HTLV-Ⅲ. It can be seen that the mechanism of action and mode of transmission of both in patients are similar. Research has shown that when HTLV enters the human body through blood infections, sexual transmission, mother to child transmission, and other means, it can attack $CD4^{+} T$ cells in the human immune system ^[20], causing irreversible damage to the body, such as being more prone to lymphoma, tropical spasmodic paraplegia, and large Karla cell leukemia ^[21].

In recent years, scholars have established many models of HTLV to study its dynamic behaviors ^{[22,23,24,25]}, but the study of HIV and HTLV co-infection is not so much. In 1990, Kobayashi et al. found that the production of tumor necrosis factors by human T cell lines infected with HTLV-1 may lead to their high susceptibility to HIV infection ^[26], which indicates that the study on HIV/HTLV co-infection is a significant work. In 2019, Mendoza et al. studied the proportion of people co-infected with HIV and HTLV in areas with high incidence of HIV in Spain ^[27]. In 2021, Alshakh et al. established an HIV/HTLV co-infection model with effective HIV-specific antibody immune response ^[28] and conducted stability analysis on the various equilibriums of the model. They divided the state variables into seven compartments, namely susceptible $CD4^{+}T$ cells $U(t)$, silent HIV-infected cells $H(t)$, active HIV-infected cells $X(t)$, silent HTLV-infected cells $N(t)$, Tax-expressing HTLV-infected cells $W(t)$, free HIV particles $V(t)$ and HIV-specific antibody $A(t)$. In the same year, Elaiw et al. established an eight-dimensional HIV/HTLV co-infection system ^[29]. Compared with the work in ^[28], they considered the impact of HIV-specific CTLs and HTLV-specific CTLs. In 2023, Elaiw et al. also conducted stability analysis on the host HTLV-I/HIV-1 co-infection model in the presence of macrophages^[30]. In the context of the frequent occurrence of various infectious diseases in the world, HIV patients, as people with low immunity, are much more likely to be infected with other diseases than ordinary people. In fact, there are many works concerning different co-infections, such as HIV/TB co-infection ^{[31,32,33,34,35]}, and HIV/COVID-19 co-infection ^[36,37,38].

Because HIV and HTLV have similar transmission route and action mechanism, many co-infected patients tend to ignore the interference of HTLV virus in the process of medical treatment. After more than 30 years of efforts, though there is still no cure for HIV virus, people have found combined treatment which can effectively alleviate and inhibit the process to AIDS. Combination treatment of synthetase inhibitors and reverse transcriptase inhibitors, can well delay the process of AIDS, so that the patients can live a longer time with high quality. However, up to now, there is still no effective vaccines or drugs for HTLV.

Fractional order calculus has been favored by many scholars in recent years because it extends integer differentiation and integration to any order, and it has more advantages in memory than integer order systems ^[39]. A variety of definitions of fractional order calculus were given in ^[40], such as, Riemann-Liouville (RL) definition, Caputo definition, Grunwald-Letnikov (GL) definition, etc. The initial value of the fractional order system is not only difficult to find, but also has no clear physical meaning. However, as a particular case of the new Hattaf mixed fractional (HMF) derivative ^[41,42], the Caputo derivative defines the initial value conditions of fractional differential equations that have the same meaning as integer order differential equations. This advantage makes the fractional order defined by Caputo widely used in engineering, physics and other fields ^{[43,44,45,46]}. Thus, the Caputo derivative is adopted in this paper.

In summary, motivated by ^[4,28], this paper will consider establishing a fractional order HIV/HTLV co-infection model with optimized control measures. Compared with integer order models, fractional order models have not been widely used in the practical application of epidemic models, and fractional order models can better explain the memory function of the immune system. Moreover, compared to constant control, optimal control not only saves costs but also is more in line with actual treatment situations. Therefore, the research method of this article is worth exploring. The specific model of this article is as follows

where $\alpha$ ($0 < \alpha < 1$) is the factional order derivative in the Caputo sense. Two control functions $u_{1}(t)$ and $u_{2}(t)$ are considered in system $(1.1)$. $u_{1}(t)$ represents the therapeutic effect of synthetase inhibitor, and $u_{2}(t)$ represents the therapeutic effect of reverse transcriptase inhibitor. Two combined drug therapies will be included in this model to analyze their impact on the dynamic behaviors. The detailed biological meanings of state variables and parameters are shown in Table 1.

It can be found that system $(1.1)$ has some defects because it is obtained by replacing the integer derivative by a fractional derivative, which can lead to asymmetry in the system's temporal dimension. The left side of system $(1.1)$ has dimensional (time)$^{-\alpha}$, while the dimension of the right side is (time)$^{-1}$. In order to unify the dimensions, we can modify it by the method in ^[47,48]. After modification, the dimensionless value on the right side of system $(1.1)$ does not need to change, and the dimension is still (time)$^{-1}$, and the dimension of other values is modified to (time)$^{-\alpha}$. The correct form of the modified system $(1.1)$ will become to the following system $(1.2)$.

with initial conditions

2.   Preliminaries

In this section, we list some basic definitions and necessary lemmas of fractional order calculus which are useful.

Definition 2.1. ^[40] The Caputo fractional derivative of order $\alpha > 0$ for a function $f:R^{+}\rightarrow R$ is defined by

where $\alpha\in(n-1, n)$, $n\in\mathrm{N}$.

Definition 2.2. ^[40] The two-parameter Mittag-Leffler function is defined by

When $\beta = 1$, the two-parameter Mittag-Leffler function becomes to the one-parameter Mittag- Leffler function, i.e.,

Lemma 2.3. ^[49] Suppose that $f(t)\in C[a, b]$ and $D^{\alpha}f(t)\in C[a, b]$ for $0 < \alpha\leq1$, then we have

(i) If $D^{\alpha}f(t)\geq0$, $\forall t\in [a, b]$, then $f(t)$ is non-decreasing for each $t\in [a, b]$.

(ii) If $D^{\alpha}f(t)\leq0$, $\forall t\in [a, b]$, then $f(t)$ is non-increasing for each $t\in [a, b]$.

Lemma 2.4. ^[50] The equilibrium of the fractional order differential system

is local asymptotically stable if both of the eigenvalues of the Jacobian matrix

evaluated at the equilibrium satisfies the following condition

3.   Qualitatively analysis of system $(1.2)$ with constant control

In this section, we will discuss a simple case for system $(1.2)$, where the control measures are constant, i.e., $u_{1}(t) \equiv u_{1}$ and $u_{2}(t) \equiv u_{2}$, $\forall t\geq 0$. Then system $(1.2)$ will become to the following system $(3.1)$.

with initial conditions

In the next of this section, we will analyze the dynamics of system $(3.1)$. Firstly, the existence and uniqueness of positive solution is proved. Then the basic reproduction number and several other thresholds will be obtained. After that, the sufficient conditions for the existence and stability of five equilibriums are derived.

3.1.   The existence and uniqueness of positive solution for system $(3.1)$

Denote

Theorem 3.1. System $(3.1)$ with any positive initial value has a unique solution and $\Omega$ is positively invariant for system $(3.1)$.

Proof. Firstly, we will prove that the solution of system $(3.1)$ with any positive initial value is always nonnegative. Based on system $(3.1)$, we have

Observing the second equation above, combined with Lemma $2.3$, it can be seen that when the initial value of H is 0, for any $U, V\geq0$, H is non decreasing, that is, $H\geq0$. Similarly, it can be inferred that as long as the initial values of all state variables are non negative, then each state variable with an positive initial value is nondecreasing, which means $\mathbf{Y}(t)\geq0$ for any $t\geq0$. As a result, the solution $\mathbf{Y}(t)$ will remain in $\mathbf{R}^7_+$.

Secondly, we need to prove that the solution of system $(3.1)$ is bounded above. Three steps are needed to achieve this goal.

Step 1.

which implies that

Since $E_\alpha (-\Phi_{1} t^{\alpha}) \geq 0$ for any $t\geq 0$, then we have

provided that $U(0)+H(0)+X(0)+ \frac{1}{\varepsilon}N(0)\leq \frac{\xi^{\alpha}}{\Phi_{1}}$.

Step 2.

Similarly to the above step, we have

which implies that

Since $E_\alpha (-\Phi_{2}t^{\alpha}) \geq 0$ for any $t\geq 0$, then we have

provided that $\frac{1}{\psi^{\alpha}}V(0) + \frac{1}{\varpi^{\alpha}}A(0)\leq N_{2}$.

Step 3.

From the fifth equation of system $(3.1)$, we can get

and by similar method, we can get

provided that $W(0)\leq N_{1}$.

To sum up, it can be seen that $\Omega$ is positively invariant for system $(3.1)$, so we only need to consider this system within $\Omega$ in the rest of this section.

Thirdly, we will show that system $(3.1)$ with any positive initial value has a unique solution.

Define the right side of system $(3.1)$ as a vector function $f(t, Y(t)):\mathrm{R^{+}}\times\mathrm{R^{7}}\rightarrow \mathrm{R^{7}}$. By using the Theorem $3.1$ and Remark $3.2$ in ^[51], we can prove the existence and uniqueness of the solution for system $(3.1)$. According to the conclusion in ^[51], system $(3.1)$ exists a unique positive solution if the following conditions are satisfied.

(ⅰ) The function $f(t, Y(t))$ is Lebesgue measurable with respect to $t\in\mathrm{R^{+}}.$

(ⅱ) The function $f(t, Y(t))$ is continuous with respect to $Y(t)$ on $\mathrm{R^{7}}.$

(ⅲ) $\frac{\partial f(t, Y(t))}{\partial Y}$ is continuous with respect to $Y(t)$ on $\mathrm{R^{7}}.$

(ⅳ) $\|f(t, Y(t))\|\le \zeta_{2} \|Y\|+k_{2},$ for almost every $t\in R^{+}$ and all $Y\in R^{7}$. Here, $\zeta_{2}$, $k_{2}$ are two positive constants.

It is obvious that the above conditions (ⅰ)-(ⅲ) are satisfied for system $(3.1)$. Next, we only need to prove that the condition (ⅳ) is satisfied for system $(3.1)$.

Denote

where $y_{1}(t) = U(t)$, $y_{2}(t) = H(t)$, $y_{3}(t) = X(t)$, $y_{4}(t) = N(t)$, $y_{5}(t) = W(t)$, $y_{6}(t) = V(t)$, $y_{7}(t) = A(t)$.

Thus, system $(3.1)$ can be rewritten as bellow.

and

Therefore, the above fourth condition is also satisfied. Combining the above arguments we get the desired result.

This completes the proof of this theorem. □

3.2.   Thresholds and the existence of five equilibriums

In this subsection, we will firstly get the four thresholds of system (3.1). Then, the sufficient conditions for the existence of five equilibriums of system (3.1) are obtained.

By using the method of the next generation matrix ^[52], we will get some thresholds. If HIV specific antibodies are ineffective, then the basic reproduction number for HIV mono-infection $R_{01}$, and HTLV mono-infection $R_{02}$, will be obtained as follows

where $\rho(A)$ represents the spectral radius of matrix $A$, and

Define $R_{0} = \max\{R_{01}, R_{02}\}$.

We also denote the following two thresholds,

which are useful for next argument.

Remark 3.2. (i) $R_{01}$ represents the reproduction number for HIV mono-infection with ineffective HIV-specific antibodies.

(ii) $R_{02}$ represents the reproduction number for HTLV mono-infection.

(iii) $R_{03}$ represents the reproduction number for HIV mono-infection with HIV specific antibody immune response.

(iv) $R_{04}$ represents the reproduction number for HTLV and HIV co-infection with HIV specific antibody immune response.

In order to obtain the equilibriums of system $(3.1)$, let the right side of Eq $(3.1)$ equal to zero, we will get the following algebraic equations

After a simple calculation, we will obtain the theorem as follows.

Theorem 3.3. (i) System (3.1) always exists a disease-free equilibrium $E_{0} = (U_{0}, 0, 0, 0, 0, 0, 0)$, where $U_0 = \frac{\xi^{\alpha}}{\eta^{\alpha}}$.

(ii) When $R_{01}>1$, there is an HIV mono-infection without antibody immune response equilibrium $E_{1} = (U_{1}, H_{1}, X_{1}, 0, 0, V_{1}, 0)$, where

(iii) When $R_{03}>1$, there is an HIV mono-infection with antibody immune response equilibrium $E_{2} = (U_{2}, H_{2}, X_{2}, 0, 0, V_{2}, A_{2})$, where

(iv) When $R_{02} > 1$, there is an HTLV mono-infection equilibrium $E_{3} = (U_{3}, 0, 0, N_{3}, W_{3}, 0, 0)$, where

(v) When $R_{04}>1$, $\frac{R_{01}}{R_{02}}>1$, there is an HIV/HTLV coexisting with antibody immune response equilibrium $E_{4} = (U_{4}, H_{4}, X_{4}, N_{4}, W_{4}, V_{4}, A_{4})$, where

3.3.   Stability analysis of the equilibriums

The Jacobian matrix for system $(3.1)$ at an arbitrary equilibrium $E^*$ is

where

It is easy to know that two eigenvalues of $J(E_{0})$ is $\lambda_1 = -\eta^{\alpha} < 0$, $\lambda_2 = -\tau^{\alpha} < 0$, and the remaining eigenvalues are determined by the following equation

where

When $R_{0} < 1$, the coefficients of Eq $(3.5)$ satisfies $a_{1} > 0, a_{2} > 0, a_{3} > 0, a_{4} > 0, a_{5} > 0.$ According to the generalized Routh-Hurwitz criteria for fractional order system ^[53], if the coefficients of Eq $(3.5)$ satisfy the following conditions, then the disease-free equilibrium $E_{0}$ is also locally asymptotically stable.

Remark 3.4. $J (E_{1}),$ $J (E_{2}),$ $J (E_{3})$ and $J (E_{4})$ represent the Jacobian matrices at the corresponding equilibriums $E_{1},$ $E_{2},$ $E_{3}$ and $E_{4}$, respectively.

According to Lemma $2.4$, we can obtain the following stability results.

Theorem 3.5. (i) If $R_{01} > 1$, then the equilibrium $E_{1}$ exists within $\Omega$, and it is locally asymptotically stable if all eigenvalues $\lambda_{i}$ of $J(E_{1})$ satisfy

(ii) If $R_{03} > 1$, then the equilibrium $E_{2}$ exists within $\Omega$, and it is locally asymptotically stable if all eigenvalues $\lambda_{i}$ of $J(E_{2})$ satisfy

(iii) If $R_{02} > 1$, then the equilibrium $E_{3}$ exists within $\Omega$, it is locally asymptotically stable if all eigenvalues $\lambda_{i}$ of $J(E_{3})$ satisfy

(iv) If $R_{04} > 1$ and $\frac{R_{01}}{R_{02}} > 1,$ then the equilibrium $E_{4}$ exists within $\Omega$, and it is locally asymptotically stable if all eigenvalues $\lambda_{i}$ of $J(E_{4})$ satisfy

Next, we will investigate the global stability for different equilibriums.

Theorem 3.6. The disease-free equilibrium $E_{0}$ is globally asymptotically stable if $R_{0} \leq 1$.

Proof. By similar method in ^[28], we will take the following Lyapunov function

Observing the above equation, it can be seen that $L_{0}(U, H, X, N, W, V, A) > 0$ for all U, H, X, N, W, V, $A > 0$, and $L_{0}(U_{0}, 0, 0, 0, 0, 0, 0) = 0.$ We calculate $\mathrm{D}^{\alpha}L_{0}$ along the solutions of system (3.1) as follows

If $R_{0}\leq1$, then $1-R_{01}\geq0$ and $1-R_{02}\geq0$, which means $D^{\alpha}L_{0}|_{(3.1)}\leq 0$. $D^{\alpha}L_{0}|_{(3.1)} = 0$ if and only if $U = U_{0}$, $W = 0$, $V = 0$. Then, from the fifth and sixth equations of system (3.1) we have N = X = 0. In addition, we know that the maximum invariant set of system $(3.1)$ on the set $\{(U, H, X, N, W, V, A)\in \Omega: \mathrm{D}^{\alpha}L_{0}|_{(3.1)} = 0\}$ is the singleton $\{E_{0}\}$. According to the LaSalle's invariance principle ^[54], we know that $E_{0}$ is global asymptotically stable.

This completes the proof of this theorem. □

We also have the following global stability result about other equilibriums.

Theorem 3.7. (i) If $R_{01} > 1$, $\frac{R_{02}}{R_{01}}\leq1$, and $R_{03}\leq1$, then equilibrium $E_{1}$ is globally asymptotically stable.

(ii) If $R_{03} > 1$, and $R_{04}\leq1$, then equilibrium $E_{2}$ is globally asymptotically stable.

(iii) If $R_{02} > 1$, and $\frac{R_{01}}{R_{02}}\leq1$, then equilibrium $E_{3}$ is globally asymptotically stable.

(iv) If $R_{04} > 1$, and $\frac{R_{01}}{R_{02}} > 1$, then equilibrium $E_{4}$ is globally asymptotically stable.

Proof. The proof of this theorem is similar to the method in ^[28], so we omit it. If the order of fractional derivative equals to one (i.e., $\alpha = 1$), and the control parameters are zero ($u_1 = 0, u_2 = 0$), then the result will be exactly the same as in ^[28]. □

4.   The fractional optimal control problem (FOCP)

In this section, we will analyze the fractional-order optimal control system (1.2) where the control parameters are not constant. Similar to the method in Section 3.1, it is easy to prove that system (1.2) with any positive initial value will have a positive solution, and it will remain within $\Omega$.

Our goal is to reduce the number of infected cells while minimizing the cost of medical treatment and the harm caused by treatment. Therefore, we define the following objective function

where $A_1$, $A_2$ are positive constants to keep a balance in the size of $H(t)$ and $X(t)$. $B_1$, $B_2$ are positive weight parameters which are associated with the control variables $u_1(t)$ and $u_2(t)$.

Define the optimal control set as follows

4.1.   Existence of an optimal control pair

Finding an optimal control solution $(u_{1}^{*}, u_{2}^{*})$ to minimize the objective function $J(u_{1}, u_{2})$ is the key to optimal control problems. According to the method in ^[4,55], we can obtain the existence of the optimal control solution as follows.

Theorem 4.1. If the optimal system meets the following conditions

(i) The sets of control variables and corresponding state variables are non-empty.

(ii) The control set $\Gamma$ is closed and convex.

(iii) The right side of system (1.2) is bounded by a linear function of the state and control variables.

(iv) The integrand of the objective function is convex on the control set $\Gamma$.

(v) There exists constant $c_{1}, c_{2} > 0$ and $c_{3} > 1$ such that the integrand $L(H, X, u_1, u_2)$ of the objective functional satisfies

then there exists an optimal control pair $(u_{1}^{*}(t), u_{2}^{*}(t))\in \Gamma$ such that

Proof. By using the existence of solutions for the bounded coefficient system (3.1.1) in ^[56], it is easy to verify that condition (ⅰ) satisfies. According to the definition of the control set $\Gamma$, $\Gamma$ is closed and convex, and condition (ⅱ) is satisfied. The proof of condition (ⅲ) is as follows.

Denote $\vec{z} = (U(t), H(t), X(t), N(t), W(t), V(t), A(t))^T$, then system (1.2) can be rewritten as

where

Denote

then we can get

Therefore

where

The above formula indicates that $G(\vec{z})$ is uniformly Lipschitz continuous, therefore, the right side of system (1.2) is constrained by a linear function of the state and control variables, and the condition (ⅲ) is proved.

Condition (ⅳ) can be proven by definition. Additionally, note that the integrand of the objective function is convex and the properties (ⅴ) is satisfied

where we take $c_{1} = \frac{1}{2}\min\{B_{1}, B_{2}\}, \quad c_{2} = 1, \quad c_{3} = 2.$

Combining the above arguments, we will get the desired result. □

4.2.   The optimal solution for control system $(1.2)$

Define the Lagrangian and Hamiltonian as follows

Lagrangian

Hamiltonian

where $\vartheta_{i}(t), i = 1, 2, \cdots, 7$ represent adjoint variables.

Theorem 4.2. Given the optimal controls $(u_{1}^{*}(t), u_{2}^{*}(t))$ and the corresponding optimal solutions $(U^{*}(t), H^{*}(t), X^{*}(t), N^{*}(t), W^{*}(t), V^{*}(t), A^{*}(t))$, there exists adjoint variables $\vartheta_{i}(t), i = 1, 2, \cdots, 7$ satisfying

with the transversal conditions

Furthermore, for $t\in[0, T],$ the optimal controls $u_{1}^{*}(t), u_{2}^{*}(t)$ are given by

Proof. Using the Pontryagin's maximum principle ^[57,58], the adjoint equations and the transversal conditions for the optimization system can be obtained as follows

with the transversal conditions

By solving the optimality conditions

we will get the optimal control pair $u_{1}^{*}(t)$ and $u_{2}^*(t)$ as follows.

Taking into account the boundedness of the control variables, the control pair $u_{1}^{*}(t), u_{2}^{*}(t)$ will be derived as

Substituting $u_{1}^{*}(t), u_{2}^{*}(t)$ into systems (1.2) and (4.1), we obtain the following optimality system

with initial conditions

and transversal conditions

□

5.   Examples and numerical simulations

In the previous sections, we have made some theoretical predictions about the dynamical behavior of systems $(1.2)$ and $(3.1)$. In this section, we will perform some numerical simulations to verify the theoretical results. In addition, the sensitive analysis of some parameters is also taken for systems $(1.2)$ and $(3.1)$. The simulation results are all based on the Adama-Bashforth-Moulton prediction correction method.

5.1.   Examples and numerical simulations for system $(3.1)$

Example 5.1. Fix the following parameter values: $\eta = 0.01$, $\rho_{2} = 0.005$, $\varpi = 0.1$, $u_{1} = 0.09$, $u_{2} = 0.02$.

Figure 1 shows the variation between threshold $R_{0}$ (i.e., $R_{01}$, $R_{02}$) and parameter $\alpha$, with different values of $\rho_{1}$ ($\rho_{1} = 0.0007, 0.0008, 0.001, 0.003$).

Example 5.2. Fix the following parameter values: $\eta = 0.1$, $\rho_{1} = 0.001$, $\rho_{2} = 0.005$, $\varpi = 0.5$, $u_{1} = 0.3$, $u_{2} = 0.2$.

(ⅰ) In Figure 2a–g, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and $\alpha$ have different values ($\alpha = 0.73, 0.85, 0.90, 1$), In this case $R_{0}\in[0.2634, 0.3125]$.

(ⅱ) In Figure 2h, the value of $\alpha$ is fixed to 0.9, and different initial values are taken with $Y_{0} = [6, 0.5, 0.5, 1, 2, 0.5, 0.5], [10, 0.9, 1, 3, 4, 0.7, 0.8], [15, 1, 3, 1.5, 3, 1, 1.5].$ In this case, we get $R_{0} = 0.2796 < 1$.

Figure 2a–g show that if $R_{0} \leq 1$, then the disease-free equilibrium $E_{0}$ is always asymptotically stable for all $\alpha\in[0.73, 1]$. Figure 2h indicates that different initial values doesn't affect the stability of equilibrium $E_{0}$.

Example 5.3. Fix the following parameter values: $\eta = 0.01$, $\rho_{1} = 0.001$, $\rho_{2} = 0.0003$, $\varpi = 0.001$, $u_{1} = 0.3$, $u_{2} = 0.2$.

(ⅰ) In Figure 3a–g, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and $\alpha$ have different values ($\alpha = 0.73, 0.85, 0.90, 1$). In this case, we get $R_{01} \in[1.4604, 2.4080]$, $\frac{R_{02}}{R_{01}} \in[0.0779, 0.1096]$, $R_{03}\in[0.3648, 0.3692]$.

(ⅱ) In Figure 3h, the value of $\alpha$ is fixed to 0.9, and different initial values are taken with $Y_{0} = [6, 0.5, 0.5, 1, 2, 0.5, 0.5], [10, 0.9, 1, 3, 4, 0.7, 0.8], [15, 1, 3, 1.5, 3, 1, 1.5].$ In this case, we get $R_{01} = 1.9957 > 1$, $\frac{R_{02}}{R_{01}} = 0.0882 < 1$, $R_{03} = 0.3689 < 1$.

Figure 3a–g show that if $R_{01} > 1$, $\frac{R_{02}}{R_{01}} < 1$ and $R_{03} < 1$, then equilibrium $E_{1}$ is always asymptotically stable for all $\alpha\in[0.73, 1]$. Figure 3h indicates that different initial values doesn't affect the stability of equilibrium $E_{1}$.

Example 5.4. Fix the following parameter values: $\eta = 0.01$, $\rho_{1} = 0.001$, $\rho_{2} = 0.0003$, $\varpi = 0.01$, $u_{1} = 0.3$, $u_{2} = 0.2$.

(ⅰ) In Figure 4a–g, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and $\alpha$ have different values ($\alpha = 0.80, 0.85, 0.90, 1$). In this case, we get $R_{03} \in[1.0747, 1.5436]$, $R_{04} \in[0.1064, 0.1202]$.

(ⅱ) In Figure 4h, the value of $\alpha$ is fixed to 0.9, and different initial values are taken with $Y_{0} = [6, 0.5, 0.5, 1, 2, 0.5, 0.5], [10, 0.9, 1, 3, 4, 0.7, 0.8], [15, 1, 3, 1.5, 3, 1, 1.5]$. In this case, we get $R_{03} = 1.2807 > 1$, $R_{04} = 0.1130 < 1$.

(ⅲ) In Figure 5, the value of $\alpha$ is fixed to 0.9, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and different values of $u_{1}, u_{2}$ are taken with $u_{1} = 0, u_{2} = 0; u_{1} = 0.17, u_{2} = 0.19; u_{1} = 0.23, u_{2} = 0.2; u_{1} = 0.3, u_{2} = 0.2; u_{1} = 0.8, u_{2} = 0.6$.

Figure 4a–g show that if $R_{03} > 1$ and $R_{04} < 1$ then equilibrium $E_{2}$ is always asymptotically stable for all $\alpha\in[0.80, 1]$. Figure 4h indicates that different initial values doesn't affect the stability of the equilibrium $E_{2}$. Figure 5 indicates that the parameters $u_{1}$ and $u_{2}$ have significant effect on the control of the disease.

Example 5.5. Fix the following parameter values: $\eta = 0.01$, $\rho_{1} = 0.0001$, $\rho_{2} = 0.003$, $\varpi = 0.001$, $u_{1} = 0.3$, $u_{2} = 0.2$.

(ⅰ) In Figure 6a–g, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and $\alpha$ have different values ($\alpha = 0.82, 0.85, 0.90, 1$). In this case, we get $R_{02} \in[1.1071, 1.8750]$, $\frac{R_{01}}{R_{02}} \in[0.128, 0.2365]$.

(ⅱ) In Figure 6h, the value of $\alpha$ is fixed to 0.9, and different initial values are taken with $Y_{0} = [6, 0.5, 0.5, 1, 2, 0.5, 0.5], [10, 0.9, 1, 3, 4, 0.7, 0.8], [15, 1, 3, 1.5, 3, 1, 1.5]$. In this case, we get $R_{02} = 1.3976 > 1$, $\frac{R_{01}}{R_{02}} \approx 0.1788 < 1$.

Figure 6a–g show that if $R_{02} > 1$ and $\frac{R_{01}}{R_{02}} < 1$, then equilibrium $E_{3}$ is always asymptotically stable for all $\alpha\in[0.82, 1]$. Figure 6h indicates that different initial values doesn't affect the stability of equilibrium $E_{3}$.

Example 5.6. Fix the following parameter values: $\eta = 0.01$, $\rho_{1} = 0.001$, $\rho_{2} = 0.005$, $\varpi = 0.1$, $u_{1} = 0.3$, $u_{2} = 0.2$.

(ⅰ) In Figure 7a–g, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and $\alpha$ have different values ($\alpha = 0.73, 0.85, 0.90, 1$). In this case, we get $R_{04} \in[1.0927, 2.9151]$, $\frac{R_{01}}{R_{02}} \in[1.01132, 1.5369]$.

(ⅱ) In Figure 7h, the value of $\alpha$ is fixed to 0.9, and different initial values are taken with $Y_{0} = [6, 0.5, 0.5, 1, 2, 0.5, 0.5], [10, 0.9, 1, 3, 4, 0.7, 0.8], [15, 1, 3, 1.5, 3, 1, 1.5].$ In this case, we get $R_{04} = 2.0198 > 1$, $\frac{R_{01}}{R_{02}} \approx 1.165 > 1$.

(ⅲ) In Figure 8, the value of $\alpha$ is fixed to 0.9, the initial value is $Y_{0} = [6, \quad 0.5, \quad 0.5, \quad 1, \quad 2, \quad 0.5, \quad 0.5],$ and different values of $u_{1}, u_{2}$ are taken with $u_{1} = 0, u_{2} = 0;\ u_{1} = 0.17, u_{2} = 0.19;\ u_{1} = 0.23, u_{2} = 0.2;\ u_{1} = 0.3, u_{2} = 0.2;\ u_{1} = 0.8, u_{2} = 0.6$.

Figure 7a–g show that if $R_{04} > 1$ and $\frac{R_{01}}{R_{02}} > 1$, then equilibrium $E_{4}$ is always asymptotically stable for all $\alpha\in[0.73, 1]$. Figure 7h indicates that different initial values doesn't affect the stability of equilibrium $E_{4}$. From Figure 8 we can see that the parameters $u_{1}$ and $u_{2}$ have significant effect on HIV/HTLV co-infection, but have no significant effect on $CD4^{+}T$ cells.

Remark 5.1. Figure 1 shows the relationship between $R_{0}$ (i.e., $R_{01}$ and $R_{02}$) and the parameter $\alpha$ under different values of $\rho_{1}$.

Strict observation indicates that when the value of $\rho_{1}$ is relatively small, the value of $R_{0}$ will less than 1, or great than 1, as the value of $\alpha$ is changing. Correspondingly, the disease-free equilibrium $E_{0}$ may be stable or unstable. However, if the value of $\rho_{1}$ is relatively large, then $R_{0}$ will always great than 1, which indicates that the disease-free equilibrium $E_{0}$ is unstable. That is to say, the parameter $\rho_{1}$ is sensitive to the dynamics of system $(1.2)$.

Remark 5.2. (i) Figure 2a–g indicate that the value of $\alpha$ will have effect on the value of the coordinate of equilibrium $E_{0}$, and also have effect on the speed towards the equilibrium $E_{0}$. When $R_{0} \leq1$, the disease-free equilibrium $E_{0}$ is always stable, which is in accordance with the result of Theorem 3.6.

(ii) Figure 2, 3, 4, 6, 7h indicate that the initial values will have no effect on the stability, which is in accordance with Theorem 3.1.

Remark 5.3. (i) Figure 3a–g indicate that the value of $\alpha$ will have effect on the value of the coordinate of equilibrium $E_{1}$, and also have effect on the speed towards the equilibrium $E_{1}$. When $R_{01} > 1$, $\frac{R_{02}}{R_{01}} < 1$ and $R_{03} < 1$, the equilibrium $E_{1}$ is stable, which is in accordance with the result of Theorem 3.7.

(ii) Figure 4a–g indicate that the value of $\alpha$ will have effect on the value of the coordinate of equilibrium $E_{2}$, and also have effect on the speed towards the equilibrium $E_{2}$. When $R_{03} > 1$ and $R_{04} < 1$, the equilibrium $E_{2}$ is stable, which is in accordance with the result of Theorem 3.7.

(iii) Figure 6a–g indicate that the value of $\alpha$ will affect the value of the coordinate of equilibrium $E_{3}$, and also affect the speed towards the equilibrium $E_{3}$. When $R_{02} > 1$ and $\frac{R_{01}}{R_{02}} < 1$, the equilibrium $E_{3}$ is stable, which is in accordance with the result of Theorem 3.7.

(iv) Figure 7a–g indicate that the value of $\alpha$ will affect the value of the coordinate of equilibrium $E_{4}$, and also affect the speed towards the equilibrium $E_{4}$. When $R_{04} > 1$ and $\frac{R_{01}}{R_{02}} > 1$, the equilibrium $E_{4}$ is stable, which is in accordance with the reslut of Theorem 3.7.

Remark 5.4. (i) Figure 5 shows how the values of $u_{1}$ and $u_{2}$ affect the stability of equilibrium $E_{2}$. When the values of $u_{1}$ and $u_{2}$ are relatively small, the equilibrium $E_{2}$ is stable. As the values of $u_{1}$ and $u_{2}$ increasing, the number of HIV active and silent infected cells will decrease. If the values of $u_{1}$ and $u_{2}$ are large enough, then the equilibrium $E_{2}$ will lose its stability and the dynamics will towards to the disease-free equilibrium $E_{0}$.

(ii) Figure 5 also shows that for the special case of HIV mono-infection under constant control, HIV drug control measures have a significant impact on the number of $CD4^{+}T$ cells, the content of HIV particles, and antibody immune response.

Remark 5.5. (i) Figure 8 shows that for the case of HIV/HTLV co-infection under constant control, when the values of $u_{1}$ and $u_{2}$ increase, both HIV active infected cells and silent infected cells will decrease, while HTLV infected cells will increase. That is to say, these two types of infected cells have a competitive relationship. If the values of $u_{1}$ and $u_{2}$ are large enough, then the equilibrium $E_{4}$ will lose its stability, and the dynamics will towards to the equilibrium $E_{3}$.

(ii) Figure 8 also shows that HIV drug control measures have little impact on the number of $CD4^{+}T$ cells in the co-infection model, but have a significant impact on the content of HIV virus particles.

(iii) Comparing Figures 5 and 8, we will find that for the HIV/HTLV co-infection model, the number of $CD4^{+}T$ cells can no longer be used as an authoritative basis for patients to check the treatment effect. On the contrary, HIV virus particles can be used as an important reference index, which is in accordance with the conclusion in ^[59].

5.2.   Examples and numerical simulations for system $(1.2)$ with optimal control

Example 5.6. Fix the following parameter values: $\eta = 0.01$, $\rho_{1} = 0.001$, $\rho_{2} = 0.005$, $\varpi = 0.18$, $A_{1} = 0.2$, $A_{2} = 0.8$, $B_{1} = 0.3$, $B_{2} = 0.7$.

(ⅰ) In Figure 9, two sets of values of $\alpha$ (0.73, 0.9) and initial values $Y_{01} = [6, 0.5, 0.5, 1, 2, 0.5, 0.5 ]$, $Y_{02} = [15, 1, 3, 1.5, 3, 1, 1.5 ]$ are taken. In this case, we get $\frac{R_{01}}{R_{02}} \in[1.296, 16.27]$, $R_{04} \in[1.317, 2.7263]$.

(ⅱ) Figure 10 shows the dynamic behavior of system (1.2) without control ($u_1 = 0, u_2 = 0$) or with optimal control ($u_{1} = u_{1}^{*}(t), u_{2} = u_{2}^{*}(t)$).

Remark 5.6. (i) Figure 9 shows that the values of $\alpha$ will affect the amplitude and stable state of $(1.2)$.

(ii) Figure 9 also shows that the initial values do not change the stability of the optimal control system $(1.2)$.

Remark 5.7. (i) Figure 10 shows that the optimal control can significantly delay and reduce the peak of HIV infected cells and reduce the number of infected cells.

(ii) Figure 10 shows that when optimal control is present, the number of HIV cells decrease in both active and silent infection, while the number of HTLV cells increase in both active and silent infection. In other words, there is competition between the two types of infected cells.

(iii) We also observe that in the co-infection case, the existence of optimal control has little influence on the content of $CD4^{+}T$ cells, but has a significant influence on the content of HIV virus particles. This also indicates that $CD4^{+}T$ cells are no longer the authoritative parameter for detecting disease development in co-infected patients when they seek medical treatment. On the contrary, the number of HIV virus particles can be used as an important reference index, which is in accordance with the conclusion in ^[59].

Remark 5.8. (i) Figure 11a shows that when $\alpha = 0.85$, the duration of synthetase drug treatment is relatively short, as shown by the blue line; when $\alpha = 1$, the treatment time of synthetase drugs is relatively longer, as shown by the red line.

(ii) Figure 11b shows that when $\alpha = 0.85$, the duration of reverse transcriptase therapy is relatively short, as shown by the blue lines; when $\alpha = 1$, the duration of reverse transcriptase therapy is relatively longer, as shown by the red line.

(iii) Figure 11 shows that when optimal control is applied, the fractional order model requires less time than the integer order model for both synthetase drug therapy and reverse transcriptase therapy, which means that less cost of treatment are needed for the fractional order system. Since drug dosage is not fixed in all course of disease development, compared with constant control, optimal control is more practical and can reduce the harm caused by over-treatment.

6.   Discussion and conclusions

In ^[28], the authors proposed an HIV/HTLV co-infection model and analyzed its dynamics. On the basis of ^[28], this article introduces two drug control methods and adopts a fractional order system with Caputo definition to better describe the memory characteristics of the immune system. Finally a fractional order HIV/HTLV co-infection model with optimal control is established. Compared with integer order models, this model has more diverse dynamic behaviors. From the perspective of clinical medicine, the article ^[59] concluded that the number of free HIV virus particles can more accurately reflect the disease progression of patients in the context of HIV/HTLV co-infection than the concentration of $CD4^{+}T$ cells. This article verified this conclusion from a mathematical perspective by analyzing model dynamics and numerical simulations. The specific conclusions of this article are as follows.

For the constant control case, we get the following main results.

$\diamondsuit$ The existence and uniqueness of the positive solution is proved.

$\diamondsuit$ Some critical thresholds ($R_{01}$–$R_{04}$) are derived.

$\diamondsuit$ The sufficient conditions for the existence and stability of five equilibriums are obtained.

For the optimal control case, the expression of the optimal solutions are obtained by using the Pontryagin maximum principle.

In addition, through numerical simulations we get the following results and conclusion.

$\diamondsuit$ From Figure 1, we will find that the value of $\alpha$ has a significant impact on the threshold $R_{01}$ and $R_{02}$, which in turn affects the stability of the disease-free equilibrium.

$\diamondsuit$ Figure 1 also shows that the value of $\rho_{1}$ is very sensitive to the dynamics of the system.

$\diamondsuit$ Figures 2–7 shows that the initial value does not affect the stability of the system. However, the value of $\alpha$ is sensitive to the dynamics of the system, and if affects the rate towards to the equilibriums.

$\diamondsuit$ Comparing Figures 8 and 10, we can see that optimal control can delay and reduce the peak of HIV-infected cells more efficiently than constant control.

$\diamondsuit$ Figures 8 and 10 indicate that in the co-infection case, drug therapy can effectively reduce the number of HIV-infected cells. However, the number of HTLV-infected cells will increase. This show that these two types of infected cells have a competitive relationship.

$\diamondsuit$ Combining Figure 5, Figure 8 and Figure 10, we will see that in the case of HIV/HTLV co-infection, drug control do not cause too much fluctuation in the content of $CD4^{+}T$ cells, while for the case of HIV infection alone, the effect of drug control measures can be well demonstrated on $CD4^{+}T$ cells. This indicates that when co-infected patients seek medical attention, the content of $CD4^{+}T$ cells can no longer be used as an effective reference data for understanding the development process of the disease. On the contrary, the number of HIV virus particles can be used as an important indicator, which is in accordance with the conclusion in ^[59].

$\diamondsuit$ Figure 11 shows that the optimal control is Bang-Bang type, and no singular solution is found.

$\diamondsuit$ Figure 11 also shows that under optimal control strategy, the fractional order model requires less cost than the integer order model, which means that fractional order system is better than the integer system.

There are also some meaningful topics to study in this article, such as

(ⅰ) It is not an instant for an individual to be infected and become infectious, so it is more practical to consider adding time delay in the system. In future, a fractional order HIV/HTLV co-infection model with time delay can be established.

(ⅱ) The state variables in this system all have the same fractional order derivatives. If different fractional order derivatives are used for different state variables, will we get more diverse dynamic behaviors? This is also a direction that can be attempted in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Ruiqing Shi is supported by "Research Project Supported by Shanxi Scholarship Council of China (No. 2021-091)". Yihong Zhang is supported by "2022 Science and Technology Innovation Project for Shanxi Normal University (No. 2022XSY013)". The authors would like to thank the anonymous reviewers for their helpful comments and suggestions, which improved the quality of this paper greatly.

Conflict of interest

The authors declare that they have no competing interests.