Fractional calculus used in modeling complex dynamical systems where memory and hereditary effects are involved. Fractional calculus is a powerful extension of classical differential equations by incorporating memory effects through derivatives of non-integer order. The fundamental theory of fractional differential equations was developed by Podlubny and Kilbas et al. , while numerical and analytical aspects were extensively studied by Diethelm et al. . In particular, the Caputo fractional derivative widely used in biological modeling, as it allows the use of classical initial conditions while capturing memory-dependent dynamics. As a result, fractional-order models have attracted considerable attention of the researchers working across various scientific fields. Human Immunodeficiency Virus type 1 (HIV-1) continues to be a major health concern due to its complex interaction with the human immune system and its ability to develop a long-term infection.

In recent years, many researchers successfully applied fractional-order models to various biological and epidemiological systems. Ahmed et al. demonstrated that fractional predator–prey and rabies models exhibit more dynamical behavior compared to classical models. Ameen and Novati investigated fractional epidemic models using implicit Adam’s methods and shown the improved numerical performance of fractional approaches. Similar advancements have been made in fractional SEIR epidemic modeling and eco-epidemiological systems confirming that fractional frameworks can better represent complex biological interactions. Arafa et al. introduced a fractional model to describe the interaction between HIV and CD4+ T-cells during primary infection, highlighting the importance of fractional dynamics in representing biological memory effects. Extensions of classical HIV models incorporating delay, immune response, and nonlinear incidence rates have also been investigated. Guo et al. analyzed a delayed HIV-1 infection model with saturation incidence rate and cytotoxic T lymphocyte (CTL) immune response. Further developments include fractional time-delay models that capture the proliferation of CD4+ T-cells under antiretroviral therapy, as studied by Liu et al. . Lichae et al. proposed a fractional differential model incorporating antiviral drug treatment to examine the dynamics of HIV infection more accurately. Similarly, Naik et al. explored the mechanics of viral kinetics under immune control using a fractional-order framework, demonstrating the importance of fractional derivatives in describing treatment effects during primary infection. Different types of fractional operators have also been used to study HIV dynamics. Moore et al. proposed a Caputo–Fabrizio fractional model for HIV/AIDS that includes a treatment compartment, offering a nonsingular kernel formulation. Stability analysis of HIV models using incommensurate fractional-order systems was investigated by Dasbasi . Günerhan et al. further analyzed fractional HIV models using Caputo and proportional Caputo operators, while Nazir et al. examined HIV dynamics using nonsingular kernel fractional derivatives.

Recent research has focused on developing more advanced computational and analytical methods for fractional HIV models. Khater et al. investigated stable computational solutions of an Atangana–Baleanu fractional HIV infection model involving CD4+ T-cells. Kongson et al. analyzed a generalized Caputo fractional derivative model for HIV dynamics with treatment. Additionally, Pitchaimani and Saranya Devi studied analytical properties of delayed fractional-order HIV models. Wang et al. examined the global stability of switched HIV/AIDS models involving fractional derivatives and drug treatment strategies. Further contributions include the dynamical analysis of post-treatment HIV infection models by Pradeesh et al. , which provide insights into long-term viral dynamics. Almoneef et al. recently investigated fractional HIV models under proportional Hadamard–Caputo operators, demonstrating the flexibility of different fractional derivative definitions in modeling biological processes.

To analyze such fractional systems, several analytical and numerical techniques have been developed. El-Sayed et al. presented analytical and numerical methods for multi-term nonlinear fractional differential equations. The fractional differential transform method has been applied to solve fractional differential–algebraic equations, as shown by Ibis et al. . Numerical solution techniques for nonlinear fractional systems have also been studied by Ullah et al. and Ziada , including decomposition-based methods. Most recently, optimal control strategies combined with fractional-order HIV/AIDS models have been explored to improve treatment efficiency. Hussein and Mebrate developed a fractional-order HIV/AIDS model incorporating treatment and optimal control using Caputo derivatives, highlighting the potential of fractional calculus in designing effective intervention strategies.

Many HIV fractional models are existing but do not include CTL immune saturation. Some studies include treatment but ignored fractional memory effects. Few studies analyze existence, stability. Therefore, there is a need to develop a fractional-order HIV model that simultaneously incorporates antiretroviral treatment, CTL immune response with saturation, and rigorous stability analysis.

The present work proposes a fractional-order HIV-1 model formulated in the Caputo sense. The model incorporates treatment efficacy, CTL immune response, and a saturated immune stimulation term avoiding unrealistic immune overactivation. A rigorous theoretical analysis is carried out, including positivity and boundedness of solutions, existence and uniqueness results based on fixed-point theory, derivation of the basic reproduction number $R_0$, and local and global stability analysis of equilibrium points.

In addition to the theoretical investigation, a detailed comparative study of semi-analytical methods is performed moreover the results of semi analytical methods are compared with a fractional predictor–corrector numerical scheme. The comparison is conducted in terms of convergence, accuracy, and computational performance. Numerical experiments further illustrate the impact of fractional order and treatment efficacy on disease progression.

The findings of this study indicate that fractional-order modeling provides a more flexible and realistic description of HIV-1 dynamics. Moreover, the comparative framework developed here offers practical guidance for selecting appropriate solution techniques for nonlinear fractional epidemiological systems.

This section deals with basic definitions of Caputo fractional derivative and integration as given in .

Let $\left.\mathrm{\alpha }\in [\mathrm{n}-1,\mathrm{n}\right)$, where $\mathrm{n}\in \mathbb{N}$, and let $\mathrm{f }$be an $\mathrm{n}$-times continuously differentiable function on $[a,\mathit{ t}]$. The Caputo fractional derivative of order $\mathrm{\alpha }$ of $f\left(t\right)$is defined by

where $\mathrm{\Gamma }\left(\cdot \right)$is the Gamma function.

Let $\alpha >0$ and let $f$ be a function defined on $[a,\mathit{ t}]$. The fractional integral of order $\alpha$ of $f\left(t\right)$ is defined by

where $\mathrm{\Gamma }\left(\cdot \right)$is the Gamma function.

In this section, we formulate a novel fractional-order mathematical model describing the dynamics of HIV-1 infection by incorporating memory effects, antiretroviral treatment efficacy, and immune response mechanisms.

Let the state variables be defined as follows:

$x\left(t\right)$: concentration of HIV-infected ${\mathrm{CD}4}^{+}$T cells at time $t$;

$y\left(t\right)$: concentration of free HIV-1 virions in plasma at time $t$;

$z\left(t\right)$: population of cytotoxic T lymphocyte (CTL) immune cells at time $t$.

All parameters are assumed to be positive constants and are described as follows:

$a$: infection rate of ${\mathrm{CD}4}^{+}$ T cells by free virus;

$b$: natural death rate of infected cells;

$k$: rate of virus production by infected cells;

$c$: clearance rate of free virus particles;

$p$: rate at which CTLs eliminate infected cells;

$r$: activation rate of CTLs due to infected cells;

$q$: saturation constant representing immune exhaustion;

$d$: natural death rate of CTLs;

$\epsilon$($0\le \epsilon \le 1$): efficacy of antiretroviral therapy;

$\alpha$($0<\alpha \le 1$): fractional order describing memory effects in the system.

To account for memory and hereditary effects inherent in biological processes, the classical integer-order derivatives are replaced by Caputo fractional derivatives of order $\alpha$. In recent years, several fractional-order HIV models have been proposed to study the interaction between infected ${\mathrm{CD}4}^{+}$ T cells, viral particles, and immune response . The proposed fractional-order HIV-1 model is given by

with initial conditions:

where $0<\alpha \le 1$ and the term $a(1-\epsilon )y\left(t\right)$ represents new infections reduced by antiretroviral treatment.

The term $\mathit{px}\left(t\right)z\left(t\right)$ models the clearance of infected cells by CTL immune response. The nonlinear term $\frac{\mathit{rx}\left(t\right)}{1+\mathit{qx}\left(t\right)}$ represents immune stimulation with saturation, preventing unrealistic immune overactivation. The fractional derivative introduces memory effects, meaning the current infection rate depends on past states.

When $\alpha =1$, the model reduces to a classical integer-order HIV-1 system.

When $\epsilon =0$, the model represents untreated HIV infection.

When $z\left(t\right)=0$, the model reduces to a Perelson-type viral dynamics model.

These cases confirm the consistency and generality of the proposed formulation.

The proposed model captures essential biological mechanisms of HIV-1 infection, including viral replication, immune response, and therapeutic intervention, while incorporating memory effects through fractional-order derivatives. This model provides a more realistic representation of HIV-1 dynamics and serves as a foundation for analytical and numerical investigations.

In this section, we analyze the qualitative behavior of the proposed fractional-order HIV-1 model .

Theorem 1:

If the initial conditions satisfy $x\left(0\right)\ge 0,\mathit{ y}\left(0\right)\ge 0,\mathit{ z}\left(0\right)\ge 0,$ then the solutions of the fractional HIV-1 system remain non-negative for all $t>0$.

Proof:

Consider the first equation

If $x\left(t\right)=0$, then $D_t^{\alpha }x\left(t\right)=a\left(1-\epsilon \right)y\left(t\right)\ge 0.$

Thus, the trajectory cannot cross the boundary $x=0$.

Similarly, for the virus equation, $D_t^{\alpha }y\left(t\right)=\mathit{kx}\left(t\right)-\mathit{cy}\left(t\right).$

If $y\left(t\right)=0$, then $D_t^{\alpha }y\left(t\right)=\mathit{kx}\left(t\right)\ge 0.$

For the CTL population, $D_t^{\alpha }z\left(t\right)=\frac{\mathit{rx}\left(t\right)}{1+\mathit{qx}\left(t\right)}-\mathit{dz}\left(t\right).$

If $z\left(t\right)=0$, $D_t^{\alpha }z\left(t\right)=\frac{\mathit{rx}\left(t\right)}{1+\mathit{qx}\left(t\right)}\ge 0.$

Therefore, solutions starting in the non-negative region remain non-negative.

Hence, the feasible region $\mathrm{\Omega }=\left\{\right(x,y,z)\in {\mathbb{R}}_{+}^3\}$ is positively invariant.

The system admits a unique disease-free equilibrium given by $E_0=(0,0,0).$

Proof: At equilibrium, set the right-hand sides of the system to zero:

From the second equation $y=\frac{k}{c}x.$

If $x=0$, then $y=0$. Substituting $x=0$ into the third equation gives $z=0.$

Hence the only equilibrium with zero infection is $E_0=(0,0,0).$

Theorem 3:

The basic reproduction number of the system is $R_0=\frac{\mathit{ak}(1-\epsilon )}{\mathit{bc}}$

Proof: The reproduction number measures the expected number of secondary infections produced by one infected cell in a fully susceptible environment. Near the disease-free equilibrium, the infection dynamics are governed by

Using the next-generation matrix approach, the infection matrix $F$ and transition matrix $V$ are

The spectral radius of $F{V}^{-1}$ gives

Theorem 4:

The disease-free equilibrium $E_0$ is locally asymptotically stable if $R_0<1$ and unstable if $R_0>1.$

Proof:

The Jacobian matrix of the system is

Evaluating at $E_0$,

One eigenvalue is ${\lambda }_1=-d<0.$

The remaining eigenvalues satisfy ${\lambda }^2+\left(b+c\right)\lambda +\left(\mathit{bc}-\mathit{ak}\left(1-\epsilon \right)\right)=0.$

For stability , all eigenvalues must satisfy the fractional stability condition $\left|\mathit{arg}\left({\lambda }_i\right)\right|>\frac{\mathit{\alpha \pi }}{2}$

If $\mathit{bc}-\mathit{ak}\left(1-\epsilon \right)>0$ or equivalently $R_0<1$, all eigenvalues lie in the stable region and the disease-free equilibrium is locally asymptotically stable.

Theorem 5:

If $R_0>1$, then the system admits a positive endemic equilibrium $E^{*}=(x^{*},y^{*},z^{*}).$

Proof:

From the equilibrium conditions $k{x}^{*}=c{y}^{*}$

we obtain $y^{*}=\frac{k}{c}{x}^{*}.$

Similarly, $d{z}^{*}=\frac{r{x}^{*}}{1+q{x}^{*}}$

which gives $z^{*}=\frac{r{x}^{*}}{d\left(1+q{x}^{*}\right)}.$

Substituting these into the first equilibrium equation yields

A positive solution $x^{*}>0$ exists when $R_0>1$.

Consequently $y^{*}>0$ and $z^{*}>0$, establishing the endemic equilibrium.

Theorem 6:

The endemic equilibrium $E^{*}$ is locally asymptotically stable if all eigenvalues of the Jacobian evaluated at $E^{*}$ satisfy $\left|\mathit{arg}\left({\lambda }_i\right)\right|>\frac{\mathit{\alpha \pi }}{2}$

Proof:

Linearizing the system around $E^{*}$ gives the Jacobian matrix $J\left(E^{*}\right)$. The corresponding characteristic polynomial is ${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0$

where $A_1,A_2,A_3>0$ for biologically realistic parameters.

If these coefficients satisfy the Routh–Hurwitz conditions and the fractional stability criterion, then all eigenvalues lie in the stable region of the complex plane. Therefore, the endemic equilibrium is locally asymptotically stable.

Theorem 7:

The disease-free equilibrium $E_0$ is globally asymptotically stable if $R_0<1$.

Proof:

Consider the Lyapunov function

Where $A=\frac{k}{c},\mathit{ B}=1.$

This function is positive definite in the feasible region. Taking the Caputo derivative along system trajectories gives

If $R_0<1$, then $\mathit{bc}-\mathit{ak}\left(1-\epsilon \right)>0$

which implies $D_t^{\alpha }V\le 0.$

Thus, the Lyapunov function decreases along system trajectories, ensuring global stability of $E_0$.

Theorem 8:

The fractional HIV-1 system is Ulam–Hyers stable.

Proof:

Let $\tilde{X}\left(t\right)$ be an approximate solution satisfying

with $\mid ϵ\left(t\right)\mid \le \delta$.

Since the function $F\left(X\right)$ is Lipschitz continuous, the fractional Grönwall inequality implies

for some constant $C>0$. Therefore, small perturbations in the system produce bounded deviations from the exact solution, establishing Ulam–Hyers stability.

In this section semi analytical methods and a fractional Adams-Bashforth-Moulton predictor corrector scheme are employed for numerical validation and comparison of the obtained solutions.

The fractional differential transform of a function $f\left(t\right)$ is defined as

and the inverse transform is given by

Let $\mathrm{X}\left(\mathrm{k}\right)\leftrightarrow \mathrm{x}\left(\mathrm{t}\right),\mathrm{ Y}\left(\mathrm{k}\right)\leftrightarrow \mathrm{y}\left(\mathrm{t}\right),\mathrm{ Z}\left(\mathrm{k}\right)\leftrightarrow \mathrm{z}\left(\mathrm{t}\right)\mathbf{.}$

Using the properties of fractional DTM, the transformed system of equations becomes

Infected cells $x\left(t\right)$

Viral load $y\left(t\right)$

Immune response $z\left(t\right)$

For the nonlinear immune term, a series expansion is used:

Thus, the recurrence relation is

From the initial conditions $X\left(0\right)=x_0,\mathit{ Y}\left(0\right)=y_0\mathit{ and\; Z}\left(0\right)=z_0$

For $k\mathit{=}0$,

For $k\mathit{=}1$,

For $k\mathit{=}2$

For $k\mathit{=}3$

By continuing the same procedure, higher-order components $X\left(5\right),\mathit{ Y}\left(5\right)\mathit{ and\; X}\left(5\right)$ can be obtained recursively.

The approximate analytical solution of the system is given by

The Adomian Decomposition Method (ADM) is a technique used to solve linear and nonlinear differential equations without linearization or discretization.

Consider a general nonlinear differential equation

where $L$ is the highest-order linear operator, $R$ is the remaining linear operator, $\mathit{Nu}$ contains the nonlinear terms.

Operating the inverse operator $L^{-1}$ on both sides of equation (8) gives

where $\phi$ is integration constant satisfying the condition $\mathit{L\phi }=0$.

The solution u can be represented as an infinite series of the form

Also assume that the nonlinear term $\mathit{Nu}$ can be written as an infinite series using Adomian polynomials $A_n$ of the form

where the Adomian polynomials $A_{\mathit{n }}$are evaluated using the formula

Then using equations (10), (11) and (12) into equation (9) we get

In which each term is given by the recurrence relation

Applying the fractional integral operator $D^{-\alpha }$ to equation (1) and using (2), the iterations are as follows

Initial Iteration $x_0\left(t\right)=x_0,\mathit{ }{y}_0\left(t\right)=y_0,\mathit{ }{z}_0\left(t\right)=z_0$

Adomian polynomials $A_0=x_0{z}_0,\mathit{ }{\mathit{ B}}_0=x_0-q{x}_0^2$

for $n=1$

Then

for $n=2$

for $n=3$

By continuing the same procedure, higher-order components $x_5\left(t\right),\mathit{ }{y}_5\left(t\right)\mathit{ and }{z}_5\left(t\right)$, can be obtained recursively.

The approximate analytical solution of the system is given by

For a nonlinear problem, consider

The operator $A$ can be decomposed as

$A\left(y\right)=L\left(y\right)-N\left(y\right)$, where, $L$ is linear part and $N$ is nonlinear part. a homotopy $H(y,\xi )$ is constructed as

where $\xi \in [0,1]$ is an embedding parameter and $y_0$ is an initial approximation.

The solution is assumed in the form

Setting $\xi \to 1$ yields the approximate analytical solution. HPM is computationally simple, avoids small-parameter assumptions, and provides accurate solutions for nonlinear fractional differential equations .

Introduce an embedding parameter $\xi \in [0,1]$ and define the homotopy

Assume solutions in power series of $\xi$

equating powers of $\xi$ we get,

Let $A_1=a(1-\epsilon )y_0-b{x}_0-p{x}_0{z}_0$,

Since RHS are constants terms.

Therefore $x_1\left(t\right)=\frac{A_1}{\mathrm{\Gamma }\left(\alpha +1\right)}{t}^{\alpha },\mathrm{ }$

Similarly,

Substitute known expressions

where, $A_2=a(1-\epsilon )B_1-b{A}_1-p(x_0{C}_1+z_0{A}_1)$

Similarly

For $\frac{\mathit{rx}}{1+\mathit{qx}}$ with $x=x_0+\xi x_1$.

Using Taylor expansion about $x_0$,

Thus $C_2=\frac{r{A}_1}{\left(1+q{x}_0{)}^2\right.}-d{C}_1\mathrm{ and }{z}_2\left(t\right)=\frac{C_2}{\mathrm{\Gamma }\left(2\alpha +1\right)}{t}^{2\alpha }$

By continuing the same procedure, higher-order components $x_3\left(t\right),\mathit{ }{y}_3\left(t\right)\mathit{ and }{z}_3\left(t\right)$, can be obtained recursively.

Consequently, the solution of the system can be expressed as a convergent series in terms of the embedding parameter $\xi$. Finally, by setting $\xi =1$, the homotopy series yields the approximate analytical solution of the original fractional-order system as