Stability and Direction of Hopf Bifurcation with Optimal Control Analysis of HIV Transmission Dynamics

1. Introduction

AIDS (Acquired Immunodeficiency Syndrome) can develop from HIV (Human Immunodeficiency Virus), a chronic infection that targets the immune system, if treatment is not received. Effective treatment enables people with HIV to live long, healthy lives even though there is no cure [1]. Significant advancements in long-acting treatments and cure-focused research are highlighted by recent HIV research. Clinical trials are showing promising results in both prevention and potential cures, and scientists are moving beyond daily pills to treatments that can be administered only twice a year [2]. Researchers have made great strides in 2025 with long-acting injectables and innovative PrEP options, which lessen the burden of daily medication and enhance adherence for HIV-positive individuals. Significant advancements in long-acting injectables and innovative PrEP options, which lessen the burden of daily medication and increase adherence for individuals living with HIV, have been reported by researchers in 2025. Data on a investigational twice-yearly regimen that combined lenacapavir with broadly neutralizing antibodies (bNAbs) was presented by Gilead Sciences. The regimen met its primary endpoint in Phase 2 trials and was designated as a breakthrough therapy [3]. The first HIV cure clinical trial, which was carried out in South Africa, also showed promising early results, highlighting the global movement toward functional cures. These advancements are significant because they improve quality of life and advance the scientific community’s long-term objective of HIV eradication [4].

HIV transmission dynamics emphasize the critical role of Highly Active Antiretroviral Therapy (HAART) in reducing viral load and thereby lowering the probability of transmission. Mathematical models consistently show that HAART suppresses viral replication, driving the basic reproduction number ($R_0$) below unity and stabilizing the system toward a disease-free equilibrium [5,6]. Interleukin-2 (IL-2), a cytokine that promotes ${\mathrm{CD4}}^{+}$ T-cell proliferation, has been investigated as an adjunct to HAART. Clinical trials demonstrated that IL-2 administration increases ${\mathrm{CD4}}^{+}$ counts, enhancing immune reconstitution, though without significant improvement in long-term clinical outcomes compared to HAART alone [7,8]. Thus, while HAART directly reduces transmission risk by lowering viral load, IL-2 contributes indirectly by bolstering immune recovery, highlighting a complementary but limited role in transmission dynamics.

Mathematical modeling of viral dynamics provides a rigorous framework for designing therapeutic strategies and assessing the efficacy of antiviral interventions. Extensive research has examined the global stability of models that capture the progression of viral infections, particularly in the context of HIV. The host immune response is a critical determinant in limiting viral proliferation, and its modulation significantly influences disease outcomes. Modern antiretroviral agents are capable of suppressing viral replication and reducing immune activation, thereby mitigating apoptosis. Nevertheless, while Highly Active Antiretroviral Therapy (HAART) has proven effective in controlling viral load, it only partially restores immune function. Complete eradication of HIV through HAART alone remains unattainable, as viral rebound typically occurs upon cessation of therapy. Consequently, there is growing interest in complementary therapeutic approaches, such as the administration of interleukin-2 (IL-2), which may enhance immune reconstitution and support long-term viral control [9,10,11].

Hopf bifurcation analysis is an important mathematical characteristic [12]. It is helpfun in understanding the oscillatory behavior of human diseases such as HIV infection dynamics [13]. Such oscillations, captured through Hopf bifurcation, provide insight into how treatment delays, cure rates, or immune saturation effects influence the persistence or suppression of HIV. Recent studies have rigorously examined these dynamics, showing that bifurcation thresholds mark transitions between stable infection states and oscillatory regimes, thereby offering valuable guidance for therapeutic strategies and long-term disease management [14].

Interleukin-2 (IL-2) is a well-established cytokine that functions as a critical growth factor for T cells, regulating both their proliferation and differentiation across the entire T-cell compartment. Following antigenic stimulation, IL-2 is secreted by ${\mathrm{CD4}}^{+}$ and ${\mathrm{CD8}}^{+}$ T-cell subsets within peripheral lymphoid organs such as the spleen and lymph nodes. In the context of HIV infection, intermittent administration of highly active antiretroviral therapy (HAART) in combination with immune modulators like IL-2 has been proposed as a therapeutic strategy to suppress viral replication. Mechanistically, IL-2 enhances ${\mathrm{CD4}}^{+}$ T-cell activation, which can transiently reactivate latent HIV reservoirs. This reactivation allows cytotoxic T lymphocytes (CTLs) specific to HIV to recognize and eliminate infected cells, thereby contributing to viral clearance [12].

To better understand these complex interactions, mathematical models have been developed to simulate the effects of drug therapy in HIV patients. These models provide insights into viral dynamics, the emergence of drug resistance, and the decay kinetics of HIV-infected cellular compartments [15]. Furthermore, they have been applied to study thymic reconstitution in HIV-1 patients undergoing HAART, as well as to design optimal control strategies that integrate both HAART and IL-2 therapy [16]. Collectively, such modeling approaches highlight the potential of combining immunotherapy with antiretroviral treatment to improve long-term outcomes in HIV management.

Optimal control theory provides a rigorous mathematical framework for the design of drug administration strategies in HIV treatment, enabling researchers to balance therapeutic efficacy against the minimization of adverse side effects and economic costs [17]. By embedding drug scheduling protocols into dynamical models of HIV infection, one can formulate control functions that simultaneously suppress viral replication and preserve immune competence [18,19]. Within this context, the Pontryagin Maximum Principle (PMP) supplies necessary conditions for optimality, allowing the derivation of adjoint equations and the construction of the Hamiltonian function that encapsulates system dynamics, control variables, and associated cost functionals [20,21,22]. The Hamiltonian formalism facilitates the identification of optimal trajectories for antiretroviral therapy (ART), guiding both dosage and timing decisions.

Recent studies have extended these models to incorporate compartments for remission, latent reservoirs, and patient adherence, thereby reflecting the complex biological and behavioral dimensions of HIV treatment [23]. Such analyses demonstrate that optimal control strategies, derived through PMP and Hamiltonian methods, can significantly improve long-term outcomes by reducing viral load while maintaining immune system integrity [18,19]. Boukary et al. [24] developed a mathematical model considering vertical transmission and applied optimal control to reduce the basic reproduction number and stabilize disease-free equilibria. Similarly, Rana et al. [25] introduced a deterministic SHATR model (Susceptible–HIV infected–AIDS infected–Treatment–Recovered) and used nonlinear stability analysis with control strategies to evaluate rehabilitation and treatment outcomes. More recently, Kumar et al. [26] emphasized stochastic elements in HIV spread, showing how optimal control can adapt drug interventions under uncertainty. These approaches demonstrate that mathematical optimization not only refines drug application design but also enhances long-term treatment planning, offering insights into personalized therapy schedules and public health strategies.

In previous studies, the focus has primarily been on stability analysis and optimal control, while phenomena such as periodic oscillations, limit cycles, and Hopf bifurcations have not been adequately captured by mathematical models. Although investigations incorporating time delays have addressed these aspects, models without delay have not yet examined Hopf bifurcations, stability properties, or the direction of bifurcating periodic solutions.

In this study, we undertake a comprehensive investigation of the reconstitution dynamics of ${\mathrm{CD4}}^{+}$ T cells and their regulatory influence on cytotoxic T lymphocytes (CTLs) in HIV-infected individuals undergoing Highly Active Antiretroviral Therapy (HAART) in conjunction with interleukin-2 (IL-2) administration. The proposed mathematical modeling framework is designed to elucidate the interactions among distinct T lymphocyte subsets during HIV progression and to determine optimal conditions for immune-based therapeutic strategies associated with HAART. Moreover, the framework integrates principles of optimal control theory to inform the rational scheduling of IL-2 as an adjuvant therapy alongside HAART, with the overarching objective of enhancing immune restoration, suppressing viral persistence, and ultimately advancing toward a functional cure.

3. Basic Properties

We now analyze the qualitative properties of the system (1) with respect to existence, uniqueness, and positivity of solutions.

3.1. Existence and Uniqueness of Solutions

The right-hand side of system (1) is a polynomial in $(T,y,C)$ and hence belongs to $C^{\infty }\left({\mathbb{R}}^3\right)$. Therefore, it is locally Lipschitz continuous on ${\mathbb{R}}^3$. By the Picard–Lindelöf Theorem, for any initial condition $(T_0,y_0,C_0)\in {\mathbb{R}}^3$, there exists a unique local solution. We emphasize that while uniqueness and existence are guaranteed, closed-form explicit solutions are not generally available.

3.2. Positive Invariance

We now show that the nonnegative orthant ${\mathbb{R}}_{+}^3$ is positively invariant under the dynamics.

From the second equation of (1),

$$\dot{y}=\left((1-{\eta }_1){\beta }_{1}T-{\mu }_2-{\gamma }_{1}C\right)y,$$

which is linear in y. Thus, we have

$$y\left(t\right)=y\left(0\right)exp\left({\int }_0^t\left((1-{\eta }_1){\beta }_{1}T\left(s\right)-{\mu }_2-{\gamma }_{1}C\left(s\right)\right)ds\right).$$

If $y\left(0\right)=0$, then $y\left(t\right)\equiv 0$. If $y\left(0\right)>0$, then $y\left(t\right)>0$ for all t in the interval of existence. Hence $y\left(t\right)\ge 0$ for all $t\ge 0$. Similarly, the third equation of (1) is linear in C:

$$\dot{C}=\left({\gamma }_{2}y-{\mu }_3+{\epsilon }_2\right)C,$$

yielding

$$C\left(t\right)=C\left(0\right)exp\left({\int }_0^t\left({\gamma }_{2}y\left(s\right)-{\mu }_3+{\epsilon }_2\right)ds\right).$$

Thus $C\left(t\right)\ge 0$ for all $t\ge 0$ whenever $C\left(0\right)\ge 0$.

For the nonnegativity of T, we consider the first equation of (1),

$$\dot{T}=a+({\epsilon }_1-{\mu }_1)T-\alpha T^2-(1-{\eta }_1){\beta }_{1}Ty.$$

We note that, on the boundary $T=0$, $\dot{T}=a\ge 0$, so the vector field points inward. Also, the quadratic and bilinear terms are non-positive for $T,y\ge 0$, ensuring that $T\left(t\right)$ cannot cross into negative values. Therefore $T\left(t\right)\ge 0$ for all $t\ge 0$ whenever $T\left(0\right)\ge 0$.

3.3. Boundedness of the System

The right-hand side of system (1) consists of smooth polynomial functions in the variables T, y, and C, together with the system parameters. Therefore, local existence, uniqueness, and continuity of solutions are guaranteed. In this subsection, we establish that the solutions of the system are uniformly bounded under suitable parameter conditions.

Theorem 1. 

Let $X\left(t\right)=\left(T\right(t),y(t),C(t\left)\right)$ be a solution of (1) with nonnegative initial data. If ${\mu }_1>{\epsilon }_1$ and ${\mu }_3>{\epsilon }_2$, then $X\left(t\right)$ is uniformly bounded for all $t\ge 0$.

Proof. 

Firstly, we consider the first equation of (1) to obtain the inequality

$${\displaystyle \frac{dT}{dt}}\le a-({\mu }_1-{\epsilon }_1)T.$$

This is a linear differential inequality. By the comparison Theorem [29], $T\left(t\right)$ is bounded above by the solution of

$${\displaystyle \frac{du}{dt}}=a-({\mu }_1-{\epsilon }_1)u,u\left(0\right)=T\left(0\right).$$

Solving this equation, we get

$$u\left(t\right)={\displaystyle \frac{a}{{\mu }_1-{\epsilon }_1}}+\left(T\left(0\right)-{\displaystyle \frac{a}{{\mu }_1-{\epsilon }_1}}\right)e^{-({\mu }_1-{\epsilon }_1)t}.$$

Hence, we can write

$$T\left(t\right)\le max\left\{T\left(0\right),{\displaystyle \frac{a}{{\mu }_1-{\epsilon }_1}}\right\},t\ge 0.$$

For the boundedness of $y\left(t\right)$, we recall

$$\dot{y}=\left((1-{\eta }_1){\beta }_{1}T-{\mu }_2-{\gamma }_{1}C\right)y,$$

which implies

$$y\left(t\right)=y\left(0\right)exp\left({\int }_0^t\left((1-{\eta }_1){\beta }_{1}T\left(s\right)-{\mu }_2-{\gamma }_{1}C\left(s\right)\right)ds\right).$$

Since $T\left(t\right)$ and $C\left(t\right)$ are bounded and continuous, the integral remains finite for all $t>0$, so $y\left(t\right)$ is bounded.

For the boundedness of $C\left(t\right)$, we consider the third equation of (1)

$$\dot{C}=({\gamma }_{2}y-{\mu }_3+{\epsilon }_2)C.$$

Since ${\mu }_3>{\epsilon }_2$, the coefficient of C is bounded above by ${\gamma }_{2}y-({\mu }_3-{\epsilon }_2)$. Since we already proved $y\left(t\right)$ is bounded, the growth rate of C is bounded, and hence $C\left(t\right)$ itself is bounded.

Combining these results, we conclude that $X\left(t\right)=\left(T\right(t),y(t),C(t\left)\right)$ is uniformly bounded for all $t\ge 0$ if ${\mu }_1>{\epsilon }_1$ and ${\mu }_3>{\epsilon }_2$. □

3.4. Equilibrium Analysis

We now determine the equilibrium points of system (1). The model (1) admits three equilibria, namely the followings:

(i)

the disease-free equilibrium $E_{\mathrm{DF}}\left(T_{\mathrm{DF}},0,0\right),$ where

$$T_{\mathrm{DF}}={\displaystyle \frac{({\epsilon }_1-{\mu }_1)+\sqrt{{({\epsilon }_1-{\mu }_1)}^2+4\alpha a}}{2\alpha }},$$

(ii)

the CTL-free endemic equilibrium $E_{\mathrm{CF}}\left(\overline{T},y_{\mathrm{CF}}^{*},0\right)$, where

$$y_{\mathrm{CF}}^{*}={\displaystyle \frac{a+({\epsilon }_1-{\mu }_1)\overline{T}-\alpha {\left(\overline{T}\right)}^2}{(1-{\eta }_1){\beta }_1\overline{T}}},\overline{T}={\displaystyle \frac{{\mu }_2}{(1-{\eta }_1){\beta }_1}}.$$

$E_{\mathrm{CF}}$ is feasible if $y_{\mathrm{CF}}^{*}>0$.

(iii)

the endemic equilibrium $E_{\mathrm{EE}}\left(T_{\mathrm{EE}}^{*},{\displaystyle \frac{{\mu }_3-{\epsilon }_2}{{\gamma }_2}},{\displaystyle \frac{(1-{\eta }_1){\beta }_1{T}_{\mathrm{EE}}^{*}-{\mu }_2}{{\gamma }_1}}\right),$ where

$$T_{\mathrm{EE}}^{*}={\displaystyle \frac{\Delta +\sqrt{{\Delta }^2+4\alpha a}}{2\alpha }},$$

and $\Delta =\left({\epsilon }_1-{\mu }_1-(1-{\eta }_1){\beta }_1\frac{{\mu }_3-{\epsilon }_2}{{\gamma }_2}\right)$. $E_{\mathrm{EE}}$ is biologically feasible if $T_{\mathrm{EE}}^{*}>{\mu }_2/\left((1-{\eta }_1){\beta }_1\right)$ and ${\mu }_3>{\epsilon }_2$.

3.5. The Basic Reproduction Number

In this section, we derive the basic reproduction number ${\mathcal{R}}_0$ using the next-generation matrix (NGM) method. We follow the standard framework in which ${\mathcal{R}}_0$ is defined as the spectral radius of $F{V}^{-1}$, where F collects the rates of appearance of new infections and V collects the transition terms out of infected compartments.

Here, we note that the DFE is $(T_{\mathrm{DF}},0,0)$, where

$$T_{\mathrm{DF}}={\displaystyle \frac{({\epsilon }_1-{\mu }_1)+\sqrt{{({\epsilon }_1-{\mu }_1)}^2+4\alpha a}}{2\alpha }}.$$

Now, we choose the infected compartments. The only epidemiologically infected compartment is the population of infected ${\mathrm{CD4}}^{+}$ T cells, y. CTLs (C) mediate immune response but are not infected states; T represents susceptible/uninfected ${\mathrm{CD4}}^{+}$ T cells. Therefore, the infected subsystem is one-dimensional with state y.

We write the differential equation for y as follows

$$\dot{y}=\underset{F\left(y\right)}{\underbrace{(1-{\eta }_1){\beta }_{1}Ty}}-\underset{V\left(y\right)}{\underbrace{\left({\mu }_2+{\gamma }_{1}C\right)y}}.$$

At the DFE $E_0(T_{\mathrm{DF}},0,0)$, the linearization of F and V with respect to y is the followings:

$$F={\left.{\displaystyle \frac{\partial }{\partial y}}\left((1-{\eta }_1){\beta }_{1}Ty\right)\right|}_{\mathrm{DFE}}=(1-{\eta }_1){\beta }_1{T}_{\mathrm{DF}},V={\left.{\displaystyle \frac{\partial }{\partial y}}\left(({\mu }_2+{\gamma }_{1}C)y\right)\right|}_{\mathrm{DFE}}={\mu }_2.$$

Since the infected subsystem is one-dimensional, the next-generation matrix $F{V}^{-1}$ reduces to the scalar

$${\mathcal{R}}_0={\displaystyle \frac{F}{V}}={\displaystyle \frac{(1-{\eta }_1){\beta }_1{T}_{\mathrm{DF}}}{{\mu }_2}}={\displaystyle \frac{(1-{\eta }_1){\beta }_1}{{\mu }_2}}·{\displaystyle \frac{({\epsilon }_1-{\mu }_1)+\sqrt{{({\epsilon }_1-{\mu }_1)}^2+4\alpha a}}{2\alpha }}.$$

(3)

This expression makes explicit how RTI efficacy (${\eta }_1$) and IL-2 augmentation (${\epsilon }_1$) modulate the infection potential through the susceptible pool at the DFE.

Theorem 2. 

By the NGM theory[30,31], the DFE $(T_{\mathrm{DF}},0,0)$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Remark 1. 

A small introductions of infected ${\mathrm{CD4}}^{+}$ T cells die out, and the system returns to the DFE when ${\mathcal{R}}_0<1$. The immune system successfully eliminates the initial viral challenge, preventing further spread to uninfected cells. When ${\mathcal{R}}_0>1$, the DFE becomes unstable and an endemic equilibrium with $y^{*}>0$ exist (subject to feasibility conditions identified in the equilibrium analysis). In this situation, the immune response is insufficient to eradicate the virus before additional cells become infected, thereby sustaining the infection.

4. Stability of Equilibria

For local stability at a state $(T,y,C)$, linearize the system with the Jacobian $J(T,y,C)$. The equilibrium is locally asymptotically stable if and only if all eigenvalues of J have negative real parts. For a $3\times 3$ system, use the Routh–Hurwitz conditions on the characteristic polynomial.

(i) The Jacobian at disease-free equilibrium, $E_{DF}$ is

$$J(T_{DF}^{*},0,0)=\left(\begin{array}{ccc}-{\mu }_1-2\alpha T_{DF}^{*}+{\epsilon }_1& -(1-{\eta }_1){\beta }_1{T}_{DF}^{*}& 0\\ 0& (1-{\eta }_1){\beta }_1{T}_{DF}^{*}-{\mu }_2& 0\\ 0& 0& -{\mu }_3+{\epsilon }_2\end{array}\right)$$

At $E(T_{DF}^{*},0,0)$, the eigenvalues are

$$\left\{\begin{array}{c}{\lambda }_1=-{\mu }_1-2\alpha T_{DF}^{*}+{\epsilon }_1<0,\hfill \\ {\lambda }_2=(1-{\eta }_1){\beta }_1{T}_{DF}^{*}-{\mu }_2<0,\hfill \\ {\lambda }_3=-{\mu }_3+{\epsilon }_2<0.\hfill \end{array}\right.$$

For the stability of the equilibrium, we require the conditions as:

$${\epsilon }_1<{\mu }_1+2\alpha T_{DF}^{*},(1-{\eta }_1){\beta }_1{T}_{DF}^{*}<{\mu }_2,{\epsilon }_2<{\mu }_3.$$

(ii) At $E_{CF}$ the Jacobian matrix is determined as

$$J=\left(\begin{array}{ccc}-{\mu }_1-2\alpha T_{CF}^{*}-(1-{\eta }_1){\beta }_1{y}_{CF}^{*}+{\epsilon }_1& -(1-{\eta }_1){\beta }_1{T}_{CF}^{*}& 0\\ (1-{\eta }_1){\beta }_1{y}_{CF}^{*}& (1-{\eta }_1){\beta }_1{T}_{CF}^{*}-{\mu }_2& -{\gamma }_1{y}_{CF}^{*}\\ 0& 0& {\gamma }_2{y}_{CF}^{*}-{\mu }_3+{\epsilon }_2\end{array}\right)$$

At $(T_{CF}^{*},y_{CF}^{*},0)$ one eigenvalue is ${\lambda }_3={\gamma }_2{y}_{CF}^{*}-{\mu }_3+{\epsilon }_2<0⟺{\epsilon }_2<{\mu }_3-{\gamma }_2{y}_{CF}^{*}.$ Other two eigenvalues satisfy

$${\lambda }^2-A_1\lambda +A_2=0,$$

(4)

where

$$\begin{array}{ccc}\hfill A_1& =& -{\mu }_1-2\alpha T_{CF}^{*}-(1-{\eta }_1){\beta }_1{y}_{CF}^{*}+{\epsilon }_1+(1-{\eta }_1){\beta }_1{T}_{CF}^{*}-{\mu }_2,\hfill \\ \hfill A_2& =& \left[-{\mu }_1-2\alpha T_{CF}^{*}-(1-{\eta }_1){\beta }_1{y}_{CF}^{*}+{\epsilon }_1\right]\left[(1-{\eta }_1){\beta }_1{T}_{CF}^{*}-{\mu }_2\right]\hfill \\ & & +{(1-{\eta }_1)}^2{\beta }_1^2{T}_{CF}^{*}{y}_{CF}^{*}.\hfill \end{array}$$

Thus the equilibrium is stable if

$$A_1>0,A_2>0.$$

(5)

(iii) Local stability at a endemic equilibrium point $(T,y,C)$: The Jacobian matrix at $E^{*}(T,y,C)$ is

$$J=\left(\begin{array}{ccc}{J}_{11}& J_{12}& 0\\ J_{21}& J_{22}& J_{23}\\ 0& J_{32}& J_{33}\end{array}\right),$$

with

$$\begin{array}{cc}\hfill J_{11}& =-{\mu }_1-2\alpha T-(1-{\eta }_1){\beta }_{1}y+{\epsilon }_1,\hfill \\ \hfill J_{12}& =-(1-{\eta }_1){\beta }_{1}T,\hfill \\ \hfill J_{21}& =(1-{\eta }_1){\beta }_{1}y,\hfill \\ \hfill J_{22}& =(1-{\eta }_1){\beta }_{1}T-{\mu }_2-{\gamma }_{1}C,\hfill \\ \hfill J_{23}& =-{\gamma }_{1}y,\hfill \\ \hfill J_{32}& ={\gamma }_{2}C,\hfill \\ \hfill J_{33}& ={\gamma }_{2}y-{\mu }_3+{\epsilon }_2.\hfill \end{array}$$

The characteristic equation of J at $E_{EE}$ be

$${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0,$$

(6)

where,

$$\begin{array}{ccc}\hfill a_1& =& -(J_{11}+J_{22}+J_{33}),\hfill \\ \hfill a_2& =& (J_{11}{J}_{22}-J_{12}{J}_{21})+(J_{11}{J}_{33}-J_{13}{J}_{31})+(J_{22}{J}_{33}-J_{23}{J}_{32}),\hfill \\ \hfill a_3& =& -(J_{11}{J}_{22}{J}_{33}-J_{11}{J}_{23}{J}_{32}-J_{12}{J}_{21}{J}_{33}).\hfill \end{array}$$

The Routh–Hurwitz conditions for all roots of $p\left(\lambda \right)$ to have negative real parts are

$$\begin{array}{c}\hfill a_1>0,a_2>0,a_3>0,a_1{a}_2>a_3.\end{array}$$

(7)

Therefore, the local asymptotic stability conditions at the point $(T,y,C)$ are

$$\begin{array}{cc}\hfill \left(i\right)& -(J_{11}+J_{22}+J_{33})>0,\hfill \\ \hfill \left(ii\right)& (J_{11}{J}_{22}-J_{12}{J}_{21})+\left(J_{11}{J}_{33}\right)+(J_{22}{J}_{33}-J_{23}{J}_{32})>0,\hfill \\ \hfill \left(iii\right)& -(J_{11}{J}_{22}{J}_{33}-J_{11}{J}_{23}{J}_{32}-J_{12}{J}_{21}{J}_{33})>0,\hfill \\ \hfill \left(iv\right)& \left[\left(\right]-(J_{11}+J_{22}+J_{33})\right]·\left[(J_{11}{J}_{22}-J_{12}{J}_{21})+J_{11}{J}_{33}+(J_{22}{J}_{33}-J_{23}{J}_{32})\right]\hfill \\ \hfill & >-J_{11}{J}_{22}{J}_{33}+J_{11}{J}_{23}{J}_{32}+J_{12}{J}_{21}{J}_{33}.\hfill \end{array}$$

(8)

Now, we study the local Hopf bifurcation of $E^{*}$. Any of the parameters of the model may be a bifurcation parameter. Thus we assume, $\theta$ as the generic bifurcating parameter of the system.

Theorem 3. 

The system (1), undergoes Hopf bifurcation around the endemic equilibrium point $E^{*}$ whenever the critical parameter value $\theta ={\theta }^{*}$ contained in the following domain

$${\Gamma }_{HB}=\left\{\alpha \in R^{+}:a_1\left({\theta }^{*}\right)a_2\left({\theta }^{*}\right)-a_3\left({\theta }^{*}\right)=0,a_2>0,\dot{a_3}-\left(\dot{a_1}{a}_2+a_1\dot{a_2}\right)\ne 0\right\}$$

(9)

Proof. 

Using the condition, $a_1{a}_2-a_3=0$, the characteristic equation (6) becomes

$$\begin{array}{c}\hfill ({\lambda }^2+a_2)(\lambda +a_1)=0,\end{array}$$

(10)

which has three roots namely ${\lambda }_1=+i\sqrt{a_2},{\lambda }_2=-i\sqrt{a_2}$ and ${\lambda }_3=-a_1$. Therefore, a pair of purely imaginary eigenvalues exists for $a_1{a}_2-a_3=0$.

Now we verify the transversality condition. For this, we differentiate the characteristic equation (6) with respect to ${\beta }_1$ to obtain

$$\begin{array}{ccc}& & {\displaystyle \frac{d\lambda }{d{\beta }_1}}={\displaystyle \frac{{\lambda }^2\dot{a_1}+\lambda \dot{a_2}+\dot{a_3}}{3{\lambda }^2+2\lambda a_1+a_2}}{|}_{\lambda =i\sqrt{a_2}}\hfill \\ & & ={\displaystyle \frac{\dot{a_3}-\left(\dot{a_1}{a}_2+a_1\dot{a_2}\right)}{2(a_1^2+a_2)}}+i\left[{\displaystyle \frac{\sqrt{a_2}(a_1\dot{a_3}+a_2\dot{a_2}-a_1\dot{a_1}{a}_2)}{2a_2(a_1^2+a_2)}}\right].\hfill \end{array}$$

Therefore,

$$\begin{array}{c}\hfill {\displaystyle \frac{dRe\lambda }{d{\beta }_1}}{|}_{{\beta }_1={\theta }^{*}}={\displaystyle \frac{\dot{a_3}-(\dot{a_1}{a}_2+a_1\dot{a_2})}{2(a_1^2+a_2)}}\ne 0\iff \dot{a_3}-\left(\dot{a_1}{a}_2+a_1\dot{a_2}\right)\ne 0.\end{array}$$

(11)

Thus Hopf bifurcation occurs at ${\beta }_1={\theta }^{*}.$ □

4.1. Stability and Direction of Hopf bifurcation

We follow the center–manifold and normal–form method of Hassard, Kazarinoff, and Wan [32] to determine the stability and direction of the Hopf bifurcation.

Let us denote $u_1=T-T_{*}$, $u_2=y-y_{*}$, and $u_3=C-C_{*}$, and write the system near $E^{*}=(T_{*},y_{*},C_{*})$ in vector form

$$\dot{u}=Ju+F\left(u\right),u={(u_1,u_2,u_3)}^{\top },$$

(12)

where J is the Jacobian at $E^{*}$ and $F\left(u\right)=\mathcal{O}\left(\right|u{|}^2)$ collects all nonlinear terms. Denote the components of the vector field by $F_1,F_2,F_3$, i.e.,

$${\dot{u}}_i=\sum _{j=1}^3{J}_{ij}{u}_j+F_i\left(u\right),i=1,2,3,$$

(13)

and expand $F_i\left(u\right)$ into a Taylor series up to cubic order:

$$\begin{array}{cc}\hfill F_i\left(u\right)& =\sum _{p+q+r=2}{\displaystyle \frac{1}{p!q!r!}}{\displaystyle \frac{{\partial }^{p+q+r}{F}_i}{\partial u_1^p\partial u_2^q\partial u_3^r}}{|}_{u=0}{u}_1^p{u}_2^q{u}_3^r\hfill \\ \hfill & +\sum _{p+q+r=3}{\displaystyle \frac{1}{p!q!r!}}{\displaystyle \frac{{\partial }^{p+q+r}{F}_i}{\partial u_1^p\partial u_2^q\partial u_3^r}}{|}_{u=0}{u}_1^p{u}_2^q{u}_3^r+{\mathcal{O}\left(\right|u|}^4).\hfill \end{array}$$

Equivalently, using multi-index notation,

$$a_{pqr}^{\left(i\right)}:={\displaystyle \frac{1}{p!q!r!}}{\partial }_{u_1}^p{\partial }_{u_2}^q{\partial }_{u_3}^r{F}_i\left(0\right),i=1,2,3,$$

so that $F_i\left(u\right)={\sum }_{p+q+r\ge 2}{a}_{pqr}^{\left(i\right)}{u}_1^p{u}_2^q{u}_3^r$.

Suppose that at $\theta ={\theta }^{*}$ the characteristic polynomial satisfies $a_1{a}_2-a_3=0$ with $a_2>0$ and $a_1>0$, so that J has a simple pair of purely imaginary eigenvalues ${\lambda }_{1,2}=\pm i{\omega }_0$ with ${\omega }_0=\sqrt{a_2}$ and a third eigenvalue ${\lambda }_3=-a_1<0$. Let $q\in {\mathbb{C}}^3$ and $q^{*}\in {\mathbb{C}}^3$ be the right and left eigenvectors associated with $i{\omega }_0$:

$$Jq=i{\omega }_{0}q,q^{*\top }J=i{\omega }_0{q}^{*\top },$$

normalized by $\langle q^{*},q\rangle =1$, where $\langle ·,·\rangle$ is the standard complex inner product. On the two-dimensional center manifold, the solution can be represented as

$$u=zq+\overline{z}\overline{q}+W(z,\overline{z}),z\in \mathbb{C},$$

where $W=\mathcal{O}\left(\right|z{|}^2)$ is tangent to the stable subspace. Projecting onto $q^{*}$ gives the scalar amplitude $z\left(t\right)=\langle q^{*},u\left(t\right)\rangle$. The reduced dynamics for z has the normal form

$$\dot{z}=i{\omega }_{0}z+{\displaystyle \frac{1}{2}}{g}_{20}{z}^2+g_{11}z\overline{z}+{\displaystyle \frac{1}{2}}{g}_{02}{\overline{z}}^2+{\displaystyle \frac{1}{6}}{g}_{30}{z}^3+{\displaystyle \frac{1}{6}}{g}_{03}{\overline{z}}^3+{\displaystyle \frac{1}{2}}{g}_{21}{z}^2\overline{z}+{\displaystyle \frac{1}{2}}{g}_{12}z{\overline{z}}^2+{\mathcal{O}\left(\right|z|}^4),$$

(14)

where the complex coefficients $g_{jk}$ are computed from the quadratic and cubic terms of F and the eigenvectors $q,q^{*}$.

Let us define the symmetric bilinear form $B:{\mathbb{C}}^3\times {\mathbb{C}}^3\to {\mathbb{C}}^3$ and the symmetric trilinear form $C:{\mathbb{C}}^3\times {\mathbb{C}}^3\times {\mathbb{C}}^3\to {\mathbb{C}}^3$ by

$$B(\xi ,\eta )={\left(B_1(\xi ,\eta ),B_2(\xi ,\eta ),B_3(\xi ,\eta )\right)}^{\top },$$

$$B_i(\xi ,\eta )=\sum _{p+q+r=2}{a}_{pqr}^{\left(i\right)}{\xi }_1^p{\xi }_2^q{\xi }_3^{r-1}{\eta }_3+\mathrm{sym},$$

$$C(\xi ,\eta ,\zeta )={\left(C_1(\xi ,\eta ,\zeta ),C_2(\xi ,\eta ,\zeta ),C_3(\xi ,\eta ,\zeta )\right)}^{\top },$$

$$C_i(\xi ,\eta ,\zeta )=\sum _{p+q+r=3}{a}_{pqr}^{\left(i\right)}{\xi }_1^{p_1}{\eta }_1^{p_2}{\zeta }_1^{p_3}\dots ,$$

where “sym” indicates symmetrization over arguments and the sums are formed from the second- and third-order partial derivatives at $u=0$. In practice one computes B and C directly from the Taylor coefficients $a_{pqr}^{\left(i\right)}$.

The leading normal-form coefficients are (Hassard’s formulas)

$$\begin{array}{cc}\hfill g_{20}& =\langle q^{*},B(q,q)\rangle ,\hfill \\ \hfill g_{11}& =\langle q^{*},B(q,\overline{q})\rangle ,\hfill \\ \hfill g_{02}& =\langle q^{*},B(\overline{q},\overline{q})\rangle ,\hfill \\ \hfill g_{30}& =\langle q^{*},C(q,q,q)+3B\left(q,{(J-i{\omega }_{0}I)}^{-1}B(q,q)\right)\rangle ,\hfill \\ \hfill g_{21}& =\langle q^{*},C(q,q,\overline{q})+2B\left(q,{(J-i{\omega }_{0}I)}^{-1}B(q,\overline{q})\right)+B\left(\overline{q},{(J-i{\omega }_{0}I)}^{-1}B(q,q)\right)\rangle ,\hfill \\ \hfill g_{12}& =\langle q^{*},C(q,\overline{q},\overline{q})+2B\left(\overline{q},{(J-i{\omega }_{0}I)}^{-1}B(q,\overline{q})\right)+B\left(q,{(J+i{\omega }_{0}I)}^{-1}B(\overline{q},\overline{q})\right)\rangle ,\hfill \\ \hfill g_{03}& =\langle q^{*},C(\overline{q},\overline{q},\overline{q})+3B\left(\overline{q},{(J+i{\omega }_{0}I)}^{-1}B(\overline{q},\overline{q})\right)\rangle .\hfill \end{array}$$

Here ${(J\pm i{\omega }_{0}I)}^{-1}$ denote the inverses restricted to the stable subspace (they exist since $\mp i{\omega }_0$ are simple eigenvalues).

We introduce the following quantities:

$$C_1\left(0\right)={\displaystyle \frac{i}{2{\omega }_0}}\left(g_{11}{g}_{20}-2|g_{11}{|}^2-\frac{1}{3}{\left|g_{02}\right|}^2\right)+{\displaystyle \frac{1}{2}}{g}_{21}.$$

(15)

The first Lyapunov coefficient is

$$l_1={\displaystyle \frac{1}{2{\omega }_0}}\mathrm{Re}\left(C_1\left(0\right)\right).$$

(16)

If the Hopf crossing is transversal, then can write

(i)

If $l_1<0$, the Hopf bifurcation is supercritical: a stable limit cycle is born for parameter values on the side where the equilibrium becomes unstable.

(ii)

If $l_1>0$, the Hopf bifurcation is subcritical: an unstable limit cycle exists on the side where the equilibrium is still stable.

For completeness, in the notation often used in Hassard [32], one defines

$${\mu }_2^{HB}=-{\displaystyle \frac{\mathrm{Re}\left(C_1\left(0\right)\right)}{\mathrm{Re}\left({\eta }^{\prime }\left({\theta }^{*}\right)\right)}},{\beta }_2=2\mathrm{Re}\left(C_1\left(0\right)\right),T_2=-{\displaystyle \frac{\mathrm{Im}\left(C_1\left(0\right)\right)+{\mu }_2\mathrm{Im}\left({\eta }^{\prime }\left({\theta }^{*}\right)\right)}{{\omega }_0}},$$

(17)

where $\eta \left(\theta \right)$ denotes the eigenvalue branch with $\eta \left({\theta }^{*}\right)=i{\omega }_0$. These parameters encode, respectively, the direction of the bifurcation (${\mu }_2$), the orbital stability of periodic solutions (${\beta }_2$), and the period variation ($T_2$).

We have the following Theorem.

Theorem 4. 

The coefficients in (17) determine the qualitative behavior:

(i)

If ${\mu }_2^{HB}>0$, the bifurcating periodic orbits exist for $\alpha >{\theta }^{*}$; if ${\mu }_2^{HB}<0$, they exist for $\theta <{\theta }^{*}$.

(ii)

The periodic orbits are orbitally stable if ${\beta }_2<0$ and unstable if ${\beta }_2>0$.

(iii)

The period varies according to $T_2$: it increases if $T_2>0$ and decreases if $T_2<0$.

Equivalently, in terms of the first Lyapunov coefficient, the Hopf bifurcation is supercritical if $l_1<0$ and subcritical if $l_1>0$.

5. System with Optimal Drug Dosing

The primary objective of this section is to design an optimal control strategy that reduces the population of infected ${\mathrm{CD4}}^{+}$ T cells while simultaneously minimizing the overall cost associated with drug administration. In addition, the formulation seeks to maximize the concentration of healthy ${\mathrm{CD4}}^{+}$ T cells. To achieve this, the following two bounded control inputs are introduced:

(i)

$u_1\left(t\right)$: dosage of reverse transcriptase inhibitor (RTI),

(ii)

$u_2\left(t\right)$: dosage of interleukin-2 (IL-2) therapy,

with $0\le u_i\left(t\right)\le 1$ for $i=1,2$ (maximal dosing corresponds to $u_i\left(t\right)=1$, and no treatment to $u_i\left(t\right)=0$[27]). RTI reduces the infection rate by the factor $(1-{\eta }_1{u}_1)$, where ${\eta }_1\in [0,1]$ is the efficacy. IL-2 enriches uninfected T cell and CTL responses; its effect enters linearly via $(1-{\eta }_2{u}_2)$ multiplying the respective endogenous expansion rates ${\epsilon }_1,{\epsilon }_2$.

The reverse transcriptase inhibitor (RTI) is assumed to reduce the infection rate by a factor of $(1-{\eta }_1{u}_1)$, where ${\eta }_1$ denotes the efficacy of the drug and $u_1$ is the control input. Similarly, interleukin-2 (IL-2) therapy enhances the proliferation of uninfected T cells and cytotoxic T lymphocyte (CTL) responses, modeled by ${\eta }_2{u}_2$, where ${\eta }_2$ represents the efficacy of IL-2 and $u_2$ is the corresponding control input. Incorporating these controls, the modified system (1) is expressed as:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dT}{dt}}& =& a-{\mu }_{1}T-\alpha T^2-(1-{\eta }_1{u}_1){\beta }_{1}Ty+(1-{\eta }_2{u}_2){\epsilon }_{1}T,\hfill \\ \hfill {\displaystyle \frac{dy}{dt}}& =& (1-{\eta }_1{u}_1){\beta }_{1}Ty-{\mu }_{2}y-{\gamma }_{1}yC,\hfill \\ \hfill {\displaystyle \frac{dC}{dt}}& =& {\gamma }_{2}yC-{\mu }_{3}C+(1-{\eta }_2{u}_2){\epsilon }_{2}C,\hfill \end{array}$$

(18)

with initial conditions

$$T\left(0\right)=T_0,y\left(0\right)=y_0,C\left(0\right)=C_0.$$

(19)

5.1. Objective Functional

We aim to minimize drug usage while promoting immune health and suppressing infection. The cost functional is

$$\begin{array}{c}\hfill J\left(u\right)={\int }_{t_i}^{t_f}\left[P{u}_1^2\left(t\right)+Q{u}_2^2\left(t\right)-R{T}^2\left(t\right)-S{C}^2\left(t\right)\right]dt,\end{array}$$

(20)

where $P,Q>0$ are weight constants on the benefit of the cost, and $R,S>0$ penalty multipliers. The admissible control set is

$$U=\left\{(u_1,u_2):u_i:[t_i,t_f]\to [0,1]\mathrm{measurable},i=1,2\right\}.$$

The problem here is to find $u^{*}\in U$ such that $J\left(u^{*}\right)={min}_{u\in U}J\left(u\right)$ subject to the state system in (18).

5.2. Existence of the Optimal Control Pair

The existence of an optimal pair $u^{*}=(u_1^{*},u_2^{*})$ is ensure under the following conditions:

(i)

the control set U is nonempty, closed, convex, and bounded;

(ii)

the right-hand side of (18) is jointly continuous in $(x,u)$ and locally Lipschitz in $x=(T,y,C)$ for each pair $u=(u_1,u_2)$;

(iii)

the integrand in (20) is convex in u and bounded below by an $L^1$ function in t.

Here, the dynamics are affine in u, smooth in x, and the integrand is strictly convex in u (due to the presence of $P{u}_1^2+Q{u}_2^2$) and are bounded. Therefore, an optimal control pair $u^{*}$ exists for the system (18).

5.3. The Hamiltonian

We formulate the Hamiltonian as

$$\begin{array}{ccc}\hfill H& =& P{u}_1^2+Q{u}_2^2-R{T}^2-S{C}^2\hfill \\ & & +{\xi }_1\left(a-{\mu }_{1}T-\alpha T^2-(1-{\eta }_1{u}_1){\beta }_{1}Ty+(1-{\eta }_2{u}_2){\epsilon }_{1}T\right)\hfill \\ & & +{\xi }_2\left((1-{\eta }_1{u}_1){\beta }_{1}Ty-{\mu }_{2}y-{\gamma }_{1}yC\right)\hfill \\ & & +{\xi }_3\left({\gamma }_{2}yC-{\mu }_{3}C+(1-{\eta }_2{u}_2){\epsilon }_{2}C\right),\hfill \end{array}$$

(21)

where ${\xi }_1,{\xi }_2,{\xi }_3$ are adjoint variables.

According to maximum principle [33], the necessary conditions for the optimality of the system are:

(i)

The objective functional (20) is minimized subject to the state system (18) with given initial data in (2),

(ii)

Adjoint (costate) equations satisfy: ${\dot{\xi }}_i=-{\displaystyle \frac{\partial H}{\partial x_i}}$ with terminal conditions ${\xi }_i\left(t_f\right)=0$, $i=1,2,3$,

(iii)

Optimal control parameters, ($u_i^{*}\left(t\right),i=1,2$) minimize the Hamiltonian i.e. $H(T,y,C,{\xi }_1,{\xi }_2,{\xi }_3,u^{*}\left(t\right))={min}_{u\in {[0,1]}^2}H(·,u)$ a.e. on $[t_i,t_f]$.

5.4. Adjoint system

Using the second condition mentioned above, we computing the partial derivatives of (21) to get the adjoint system. The adjoint system is obtained as

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{\xi }_1}{dt}}& =& 2RT+{\xi }_1\left({\mu }_1+2\alpha T+(1-{\eta }_1{u}_1){\beta }_{1}y-(1-{\eta }_2{u}_2){\epsilon }_1\right)-{\xi }_2(1-{\eta }_1{u}_1){\beta }_{1}y,\hfill \\ \hfill {\displaystyle \frac{d{\xi }_2}{dt}}& =& {\xi }_1(1-{\eta }_1{u}_1){\beta }_{1}T-{\xi }_2\left((1-{\eta }_1{u}_1){\beta }_{1}T-{\mu }_2-{\gamma }_{1}C\right)-{\xi }_3{\gamma }_{2}C,\hfill \\ \hfill {\displaystyle \frac{d{\xi }_3}{dt}}& =& 2SC+{\xi }_2{\gamma }_{1}y-{\xi }_3\left({\gamma }_{2}y-{\mu }_3+(1-{\eta }_2{u}_2){\epsilon }_2\right),\hfill \end{array}$$

(22)

with transversality ${\xi }_i\left(t_f\right)=0$ for $i=1,2,3$. Solving this system as a boundary value problem , we get the adjoint variables, ${\xi }_i,i=1,2,3$.

5.5. Characterization of the Optimal Control Pair

Using the third condition mentioned above we get the optimal control pair. Thus, we set the gradient of the Hamiltonian with respect to the controls to zero:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial H}{\partial u_1}}& =& 2P{u}_1+{\eta }_1{\beta }_{1}Ty({\xi }_1-{\xi }_2)=0,\hfill \\ \hfill {\displaystyle \frac{\partial H}{\partial u_2}}& =& 2Q{u}_2-{\eta }_2\left({\xi }_1{\epsilon }_{1}T+{\xi }_3{\epsilon }_{2}C\right)=0.\hfill \end{array}$$

The corresponding unconstrained optimal controls are

$$\begin{array}{ccc}\hfill u_1\left(t\right)& =& {\displaystyle \frac{{\eta }_1{\beta }_{1}Ty({\xi }_2-{\xi }_1)}{2P}},u_2\left(t\right)={\displaystyle \frac{{\eta }_2\left({\xi }_1{\epsilon }_{1}T+{\xi }_3{\epsilon }_{2}C\right)}{2Q}}.\hfill \end{array}$$

(23)

Projecting onto the admissible set $[0,1]$ yields the bounded optimal controls:

$$\begin{array}{ccc}\hfill u_1^{*}\left(t\right)& =& max\left\{0,min\left\{1,{\displaystyle \frac{{\eta }_1{\beta }_{1}Ty({\xi }_2-{\xi }_1)}{2P}}\right\}\right\},\hfill \\ \hfill u_2^{*}\left(t\right)& =& max\left\{0,min\left\{1,{\displaystyle \frac{{\eta }_2\left({\xi }_1{\epsilon }_{1}T+{\xi }_3{\epsilon }_{2}C\right)}{2Q}}\right\}\right\}.\hfill \end{array}$$

(24)

From the above analysis, we have the following Theorem:

Theorem 5. 

The objective cost function $J\left(u^{*}\left(t\right)\right)$ over U attains its minimum for the optimal control $u^{*}\left(t\right)$ corresponding to the interior equilibrium $(T^{**},y^{*},C^{*})$. Moreover, there exist adjoint variables ${\xi }_1,{\xi }_2,{\xi }_3$ satisfying the following system of equations:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{\xi }_1}{dt}}& =& 2PT+{\xi }_1\left[d_t+2\alpha T+(1-{\eta }_1{u}_1)y-{\epsilon }_1(1-{\eta }_1{u}_1)\right]-{\xi }_2\left[(1-{\eta }_1{u}_1)y\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_2}{dt}}& =& {\xi }_1\left[(1-{\eta }_1{u}_1)T\right]-{\xi }_2\left[(1-{\eta }_1{u}_1)T-(d_T+{\delta }_T)-{\gamma }_{1}C\right]-{\xi }_3{\gamma }_{2}C,\hfill \\ \hfill {\displaystyle \frac{d{\xi }_3}{dt}}& =& 2QC+{\xi }_2{\gamma }_{1}y-{\xi }_3\left[{\gamma }_{2}y-{\mu }_3+{\epsilon }_2(1-{\eta }_2{u}_2)\right],\hfill \end{array}$$

(25)

together with the transversality condition ${\xi }_i\left(t_f\right)=0$ for $i=1,2,3$. Furthermore, the optimal controls are given by

$$\begin{array}{ccc}\hfill u_1^{*}\left(t\right)& =& max\left(0,min\left(1,{\displaystyle \frac{{\eta }_{1}Ty({\xi }_2-{\xi }_1)}{2A}}\right)\right),\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill u_2^{*}\left(t\right)& =& max\left(0,min\left(1,{\displaystyle \frac{{\xi }_1{\eta }_{1}T+{\xi }_3{\eta }_{2}C}{2B}}\right)\right).\hfill \end{array}$$

(27)

Remark 2. 

The optimality system comprises the coupled state equations (18), adjoint equations (22), and control laws (24), with initial conditions for $(T,y,C)$ and terminal conditions for $({\xi }_1,{\xi }_2,{\xi }_3)$. This defines a two-point boundary value problem.

6. Numerical Results

In this section, we present computer simulation results for model system (1) by using MatLab (R2016a) algorithms. We choose the set of values of parameters from Table 1. Firstly, we simulate the results from the system without optimal control, then we present the simulated results from the optimal control problem.

6.1. Dynamics of the system without control

Time series solution of the system is plotted in Figure 2, taking the values of the model parameters from Table 1. It is observed that the coexistence equilibrium is asymptotically stable. It is seen that the system moves towards the endemic equilibrium. That means the system is stable in nature.

Bifurcation diagram of the system is shown in Figure 3 by varying ${\gamma }_2$. The system bifurcates into periodic solution when gamma passes the critical value, ${\gamma }_2^{*}=0.0515$ (approx.).

Time series solution is plotted for $\beta$, and the bifurcating periodic periodic solution is unstable (see Figure 4). Also, for $\beta =0.0045<{\beta }^{*}$, periodic solution is observed in Figure 5.

Regions of stability of the equilibrium points are plotted in Figure 6 and Figure 7. The figures show that for higher infection, the endemic equilibrium is feasible. Also it is stable in nature for the set of parameters used for these figures. Also, $\alpha$ has a stabilizing role in the dynamics.

6.2. Dynamics of the System with Optimal Control

The forward–backward sweep method is a widely used numerical technique for solving optimal control problems. It works by iteratively integrating the state equations forward in time with an initial guess for the controls, and then integrating the adjoint (costate) equations backward in time using terminal conditions. At each iteration, the controls are updated pointwise according to the optimality conditions derived from the Hamiltonian, typically with bounds enforced. This forward–backward cycle is repeated until the state, adjoint, and control trajectories converge, yielding an approximate solution to the optimal control problem.

From optimal control studies (Figure 8 and Figure 9), several interesting results have been obtained. The CTL effector population is found to decrease in all the cases. Thus, optimal control approach will help in designing an innovative cost-effective safe therapeutic regimen of HAART and IL-2 where the uninfected cell population will be enhanced with simultaneous decrease in the infected cell population. Moreover, successful immune reconstitution can also be achieved with increase in precursor CTL population. Mathematical modeling of viral dynamics thus enables maximization of therapeutic outcome even in case of multiple therapies with specific goal of reversal of immunity impairment.

Figure 8, represents the control pair $u^{*}\left(t\right)$ for RTI drug for the parameter set as given in Table 1. The RTI drug is administered at nearly full level for 20 days approximately and after that it is reduced to zero at 30 days. In Figure 9, we see that during the treatment period, the infected T cells increase almost linearly, and the CTL responses also increases, whereas the virus producing cell linearly decreases.

7. Discussion and Conclusion

In this study, we derived a mathematical framework to elucidate the dynamics of the human immune system in response to HIV infection under the combined influence of Highly Active Antiretroviral Therapy (HAART) and Interleukin-2 (IL-2). Biologically, HIV infection is characterized by the progressive depletion of ${\mathrm{CD4}}^{+}$ T lymphocytes, which compromises immune surveillance and predisposes patients to opportunistic infections. Our model captures these processes by incorporating viral replication, immune activation, and therapeutic intervention, thereby providing a quantitative lens through which the interplay of infection and treatment can be examined.

Analysis of the model system revealed that stability and oscillatory behavior depend critically on the infection rate parameter ${\beta }_1$ and immune response rate, ${\gamma }_2$. For certain values of these parameters, the immune system either stabilizes or destabilizes, reflecting biological scenarios of viral suppression or resurgence. Notably, bifurcating periodic solutions observed at ${\beta }_1=0.0035$ correspond to oscillatory viral loads and T cell counts, a phenomenon consistent with clinical observations of viral blips during therapy. Using normal form theory we see that the bifurcating periodic solution is subcritical type.

From a control-theoretic perspective, we established the existence of optimal therapeutic strategies by applying the Pontryagin Minimum Principle. The derived adjoint equations and control laws formalize the biological trade-off between drug-induced viral suppression and immune stimulation. The resulting two-point boundary value problem (TPBVP) encapsulates the coupled dynamics of viral replication, immune recovery, and drug administration. Computationally, the forward–backward sweep method yielded dosing policies that balance pharmacological costs against immunological benefits, thereby reflecting the clinical objective of maximizing patient health while minimizing toxicity and expense.

Biologically, IL-2 serves as an immune activator that enhances the proliferation and survival of ${\mathrm{CD4}}^{+}$ T cells, while HAART suppresses viral replication by targeting reverse transcriptase and other viral enzymes. Our findings demonstrate that the synergistic application of these therapies not only improves the uninfected T cell population but also reduces the reservoir of infected cells. This dual effect translates into prolonged immune competence and extended patient survival. Importantly, the optimal control strategy was shown to reduce side effects and improve cost-effectiveness, aligning with clinical goals of sustainable long-term therapy.

In conclusion, the integration of mathematical modeling with biological interpretation underscores the potential of optimal control theory in designing effective HIV treatment regimens. The combination of HAART and IL-2, when administered under an optimized dosing policy, enhances immune recovery, suppresses viral persistence, and improves life expectancy in HIV-infected individuals. This framework provides a rigorous foundation for future quantitative investigations into combination therapies and highlights the importance of mathematical biology in guiding clinical decision-making.