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

: In this paper, a fractional order HIV / HTLV co-infection model with HIV-specific antibody immune response is established. Two cases are considered: constant control and optimal control. For the constant control system, the existence and uniqueness of the positive solutions are proved, and then the su ffi cient conditions for the existence and stability of five equilibriums are obtained. For the second case, the Pontryagin’s Maximum Principle is used to analyze the optimal control, and the formula of the optimal solution are derived. After that, some numerical simulations are performed to validate the theoretical prediction. Numerical simulations indicate that in the case of HIV / HTLV co-infection, the concentration of CD 4 + T cells is no longer suitable as an e ff ective reference data for understanding the development process of the disease. On the contrary, the number of HIV virus particles should be used as an important indicator for reference.


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-III.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 coinfected 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 HTLVinfected 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 α ( 0 < α < 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) −α , 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) −α .The correct form of the modified system (1.1) will become to the following system (1.2).The death rates of Hiv specific antibodies 0.1 with initial conditions U(0), H(0), X(0), N(0), W(0), V(0), A(0) ≥ 0.

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 α > 0 for a function f : R + → R is defined by Definition 2.2.[40] The two-parameter Mittag-Leffler function is defined by When β = 1, the two-parameter Mittag-Leffler function becomes to the one-parameter Mittag-Leffler function, i.e., Lemma 2.4.[50] The equilibrium of the fractional order differential system x(0) = x 0 , y(0) = y 0 , is local asymptotically stable if both of the eigenvalues of the Jacobian matrix evaluated at the equilibrium satisfies the following condition arg(eig(A)) > απ 2 .
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) ≡ u 1 and u 2 (t) ≡ u 2 , ∀t ≥ 0. Then system (1.2) will become to the following system (3.1).
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.
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 ≥ 0, H is non decreasing, that is, H ≥ 0. 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 Y(t) ≥ 0 for any t ≥ 0. As a result, the solution Y(t) will remain in 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.
which implies that Since E α (−Φ 1 t α ) ≥ 0 for any t ≥ 0, then we have Step 2. Similarly to the above step, we have Since E α (−Φ 2 t α ) ≥ 0 for any t ≥ 0, then we have From the fifth equation of system (3.1),we can get and by similar method, we can get To sum up, it can be seen that Ω is positively invariant for system (3.1),so we only need to consider this system within Ω 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)) : R + × R 7 → 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 ( It is obvious that the above conditions (i)-(iii) are satisfied for system (3.1).Next, we only need to prove that the condition (iv) is satisfied for system (3.1). Denote Thus, system (3.1) can be rewritten as bellow.
Therefore, the above fourth condition is also satisfied.Combining the above arguments we get the desired result.
This completes the proof of this theorem.□

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 ρ(A) represents the spectral radius of matrix A, and We also denote the following two thresholds, which are useful for next argument.(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.

Stability analysis of the equilibriums
The Jacobian matrix for system (3.1) at an arbitrary equilibrium E * is where 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.
According to Lemma 2.4, we can obtain the following stability results.Next, we will investigate the global stability for different equilibriums.
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 D α L 0 along the solutions of system (3.1) as follows 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) ∈ Ω : D α 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.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., α = 1), and the control parameters are zero (u 1 = 0, u 2 = 0), then the result will be exactly the same as in [28].□

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 Ω.
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

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.(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)) ∈ Γ 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 (i) satisfies.According to the definition of the control set Γ, Γ is closed and convex, and condition (ii) is satisfied.The proof of condition (iii) is as follows.
Denote ⃗ 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 then we can get where The above formula indicates that G(⃗ 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 (iii) is proved.
Figure 4a-g show that if R 03 > 1 and R 04 < 1 then equilibrium E 2 is always asymptotically stable for all α ∈ [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.
Strict observation indicates that when the value of ρ 1 is relatively small, the value of R 0 will less than 1, or great than 1, as the value of α is changing.Correspondingly, the disease-free equilibrium E 0 may be stable or unstable.However, if the value of ρ 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 ρ 1 is sensitive to the dynamics of system (1.2).
Remark 5.2.(i) Figure 2a-g indicate that the value of α 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 ≤ 1, 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 α 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, 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 α 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 α 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 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 α 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 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: η = 0.01, ρ 1 = 0.001, ρ 2 = 0.005, ϖ = 0.18, A 1 = 0.2, A 2 = 0.8, B 1 = 0.3, B 2 = 0.7.
Remark 5.6.(i) Figure 9 shows that the values of α 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 α = 0.85, the duration of synthetase drug treatment is relatively short, as shown by the blue line; when α = 1, the treatment time of synthetase drugs is relatively longer, as shown by the red line.
(ii) Figure 11b shows that when α = 0.85, the duration of reverse transcriptase therapy is relatively short, as shown by the blue lines; when α = 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.with control for =0.9 without control for =0.9

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.♢ The existence and uniqueness of the positive solution is proved.♢ Some critical thresholds (R 01 -R 04 ) are derived.♢ 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.♢ From Figure 1, we will find that the value of α has a significant impact on the threshold R 01 and R 02 , which in turn affects the stability of the disease-free equilibrium.
♢ Figure 1 also shows that the value of ρ 1 is very sensitive to the dynamics of the system.♢ Figures 2-7 shows that the initial value does not affect the stability of the system.However, the value of α is sensitive to the dynamics of the system, and if affects the rate towards to the equilibriums.
♢ 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.
♢ 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.
♢ Combining Figure 5, Figure 8 and Figure 10, we will see that in the case of HIV/HTLV coinfection, drug control do not cause too much fluctuation in the content of CD4 + T cells, while for the

Remark 3 .
2. (i) R 01 represents the reproduction number for HIV mono-infection with ineffective HIVspecific 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.

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 Γ 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 Γ.

Figure 1 .
Figure 1.The relationship between R 01 , R 02 and the parameter α, with different values of ρ 1 .

Figure 2 .Example 5 . 3 .
Figure 2. (a)-(g) are time series of system (3.1) for different values of α.(h) is the phase portrait of U(t), H(t) and N(t) for different initial values.

Figure 3 .Example 5 . 4 .
Figure 3. (a)-(g) are time series of system (3.1) for different values of α.(h) is phase portrait of U(t), H(t) and V(t) for different initial values.

Figure 4 .
Figure 4. (a)-(g) are time series of system (3.1) for different values of α.(h) is phase portrait of U(t), H(t) and A(t) for different initial values.

Figure 5 .Figure 6 .
Figure 5.The influence of different values of parameters u 1 and u 2 on the stability of equilibrium E 2 .

Figure 7 .
Figure 7. (a)-(g) are time series of system (3.1) for different values of α.(h) is the phase portrait of U(t), H(t) and N(t) for different initial values.

Figure 8 .
Figure 8.The influence of different values of parameters u 1 and u 2 on the stability of equilibrium E 4 .

Figure 9 .
Figure 9.Time series of system (1.2) with optimal control for different values of α and initial values.

Figure 10 .
Figure 10.Time series of system (1.2) without control (in red line) or with optimal control (in blue line).

Table 1 .
The biological meanings of the variables and parameters for system (1.1).