Investigation of fractal-fractional HIV infection by evaluating the drug therapy effect in the Atangana-Baleanu sense

1.   Introduction

Presently, human immunodeficiency virus (HIV) is still a serious health issue since it is a fatal, incurable disease that has claimed the lives of millions of people. HIV is the cause of acquired immunodeficiency syndrome (AIDS), which reduces the body's ability to fight disease and makes it vulnerable to other infections ^[1]. When HIV gets into the body of healthy people, it targets CD$4^{+}$ T-cells, which are white blood cells that function as an essential part of the human immune system and quickly reproduce to damage the CD$4^{+}$ T-cells. The deterioration of CD$4^{+}$ T-cells can have a wide-ranging impact on the protective capacity of a healthy immune system. However, due to medical evolution and effective HIV treatments, people living with HIV can live long and healthy lives. Antiretroviral therapy is the use of HIV drugs to treat HIV infection and protect the immune system by preventing the virus from replicating at various phases of the HIV life cycle ^[2], such as the phase of fusion inhibitors, reverse transcription inhibitors (RTIs), integrase inhibitors and protease, see; ^[3,4,5].

Mathematical modeling of dynamical systems plays a significant role in many fields of applied and natural sciences in term of understanding the dynamic behavior of problems in the real world, such as COVID-$19$ ^[6], the DateJimbo-Kashiwara-Miwa equation ^[7], magnetohydrodynamics ^[8,9], Maxwell materials ^[10], non-Newtonian fluids ^[11] and the problems described in the references cited therein. Many mathematical models have been that are created related to the HIV epidemic. It is beneficial to use a mathematical framework to determine the significance of the interaction between HIV infection and CD$4^{+}$ T-cells. Numerous models have been designed and expanded to represent the transmission phenomena of CD$4^{+}$ T-cells with HIV and describe the infection of HIV and its interaction with the immune system. A large number of publications can be found for examples, e.g., Wang and Li ^[12] created a model of HIV infection with CD$4^{+}$ T-cells and investigated the global dynamics of the model in $2006$. After that, Srivastava et al. ^[13] demonstrated a model under an RTI, which occurs before the virus is produced by the infected cell in $2009$. They divided the class of infected CD$4^{+}$ T-cells into two classes, namely pre-reverse transcription (RT) and post-RT classes which denote the stage of the infected cells in which the reverse transcription is not finished or completed, respectively. For the results, they established stability analysis at equilibrium points (EPs) as well as performed numerical simulations. See more examples in ^{[14,15,16,17,18,19]}.

Fractional calculus is a huge body of knowledge that has piqued the curiosity of numerous researchers. Modeling with fractional calculus tools provides additional advantages in term of describing the dynamics of real-world problems with memory because a fractional-order derivative depends on local conditions and historical events. Some studies have found that fractional-order models provide more accurate and dependable information on dynamical behavior than ordinary differential integer-order models, such as those previously developed for HIV ^[20,21], Rift Valley fever ^[22], dengue ^[23], chemical kinetics ^[24], COVID-$19$ ^[25] and Zika ^[26]. Various fractional operators have been extensively utilized in mathematical models to solve a wide range of real-world problems. There exist different types of operators based on the kernels utilized, such as the Liouville-Caputo operator which involves the singularity of the kernel, the Caputo-Fabrizio (CF) fractional derivative which concerns the exponential kernel, and the Atangana-Baleanu (AB) fractional derivative operator which deals with the Mittag-Leffler (ML) kernel. There is also the newly established fractal-fractional (FF) operator which combines fractional and fractal derivatives was introduced by Atangana ^[27]. The power law, exponential law, and generalized ML law were convoluted with fractal derivatives to create this new operator. The three different kernels depend on two criteria: the fractional order and fractal dimension. The memory and fractal dimension of existing fractional-order derivatives are represented by these derivatives. According to Atangana's study and investigation, the FF order operator is a good alternative to looking at the mathematical model for real-world problems.

However, the memory characteristic of fractional-order systems enables a more accurate prediction and translation of the models by allowing for the incorporation of more historical data. Additionally, the hereditary characteristic describes the genetic profile as well as the age and immune system condition. Furthermore, dealing with fractional-order systems enables a more accurate understanding of the HIV systems by taking into account their memory and hereditary properties, which are the intricate behavioral patterns of biological systems. Using fractional derivatives and applying them to various models, several different sorts of research investigations have been established for example, in $2018$, Bushnaq et al. ^[28] studied the existence theory and stability results of an HIV-AIDS infection model with a CF fractional derivative. In $2020$, a fractional mathematical model of the rotavirus epidemic with the impacts of breastfeeding and vaccination was investigated using the AB derivative; it is discussed in ^[29]. In the same year, Naik et al. ^[30] applied a fractional-order HIV epidemic model in the Caputo sense to study the effects of prostitution in the population on disease transmission. In $2021$, Shah et al. ^[31] incorporated CF into the HIV model proposed in ^[12] involving a source term for the provision of novel CD$4^+$ T-cells depending on the viral load. In $2021$, Kongson et al. ^[32], the existence theory and stability for the HIV CD$4^+$ T-cells under the condition of treatment was investigated using a generalized fractional derivative of the Caputo type. After that, the FF operators were applied to analyze the dynamics of the dengue disease model with the hospitalization class of infected cases ^[33]. In $2022$, a nonlinear fractional order system with a Caputo fractional derivative was proposed for an SEIR epidemic model of HIV transmission; additionally, an approximate solution was established by using the homotopy analysis method ^[34]. For more details and interesting analyses, see the examples of models of the human liver ^[35], cancer ^[36], HIV/AIDS ^{[37,38,39,40,41]} and an epidemic ^[42].

Since the fractional differential operators can be either local or nonlocal which are effective tools for explaining real-world issues, they are needed to develop mathematical models. In addition, the most recent operators are the FF operators, and there has not been much research on their use and application in the literature. In light of the foregoing reasoning, to make our work distinguishable from that of others, we developed the HIV model presented in ^[13] by utilizing the FF derivative in the AB sense because it gives a much more realistic result with a greater degree of freedom than integer order. In terms of the FF derivatives, the considered model explains the memory effect, as well as fractal qualities such as the fractal dimension $\beta$ and fractional order $\alpha$; these qualities are necessary for describing real-world phenomena. The existence and uniqueness of the solutions for an HIV model in the context of FF operators which also utilzes a variety of fixed-point theories (e.g., Banach and Leray-Schauder types), are investigated. We present the stability of EPs using the basic reproduction number (BRN) and analyze the stability of the solutions via various Ulam stability-based functions, such as Ulam-Hyers (UH), generalized UH (GUH), UH-Rassias (UHR) and GUH-Rassias (GUHR) functions. Moreover, we use the Newton polynomial approach, the Adams-Bashforth method and the predictor-corrector technique to find the approximated solutions and achieve the numerical form for the ML law of the proposed model for different fractal dimensions and fractional orders. We give a comparison of the three numerical algorithms and also detect some dynamic behaviors of the solutions using our simulations.

The rest of this paper is structured as follows. Section $2$ gives some foundational information on FF operators in the AB sense, formulations of fixed-point theorems and the fractional HIV model construction. In Section $3$, we demonstrate the solution's existence and positivity, as well as the EPs and BRN, and conduct a stability study on the proposed model. The unique result for the FF-HIV model is established by using Banach's fixed-point theory, while the existence result is investigated by using Leray-Schauder's nonlinear fixed-point theory in Section $4$. Various categories of stability results, known as Ulam-type stability, are extensively verified in Section $5$. In Section $6$, the discussion of numerical simulations utilizing the Newton polynomial, the Adams-Bashforth and predictor-corrector methods are provided to analyze the behaviors of the considered model. Finally, we give the conclusion of our paper in the last section.

2.   Preliminaries

2.1.   Background materials

In this subsection, we survey the associated findings using FF calculus. We look at the fundamentals of FF calculus.

Definition 2.1. (^[27]). Assume $f(t)$ is a differentiable function in $(a, b)$. If $f$ is a fractal dimension on $(a, b)$ with the order $\beta \in (0, 1]$, then the Riemann-Liouville-type FF derivative of $f$ of the order $\alpha \in (0, 1]$ with the generalized ML kernel is defined by

where

is the fractal derivative, $\mathbb{AB}(\alpha) = 1 - \alpha + \alpha/\Gamma(\alpha)$ with $\mathbb{AB}(0) = \mathbb{AB}(1) = 1$ and the ML function is defined by

with $\mathbb{C}$ being the set of the complex numbers.

Definition 2.2. (^[27]). If $f(t)$ is a continuous function in $(a, b)$, then the FF integral of $f$ with the order $\alpha \in (0, 1]$ and fractal dimension $\beta \in (0, 1]$ is defined by

Definition 2.3. (^[27]). Let $f$ be continuous on an open interval $(a, b)$; the fractal Laplace transform of the order $\alpha$ is given by

2.2.   Construction of FF-HIV model (2.5)

Consider the model of drug therapy for HIV infection that is presented in ^[13], which studied the antiretroviral therapy in the phase of RTI. The population was separated into four unknown variables ($\mathcal{T}(t)$, $\mathcal{I}(t)$, $\mathcal{V}(t)$, $\mathcal{L}(t)$), as shown in Figure 1.

They constructed the model as below.

with the initial conditions $\mathcal{T}(0) = \mathcal{T}_{0} \geq 0$, $\mathcal{I}(0) = \mathcal{I}_{0} \geq 0$, $\mathcal{V}(0) = \mathcal{V}_{0} \geq 0$ and $\mathcal{L}(0) = \mathcal{L}_{0} \geq 0$; descriptions of the unknown variables and parameters are described in Tables 1 and 2.

This work further develops the model in ^[13] using a fractional framework, which allows for more precise and realistic predictions of the dynamic behavior over time. Therefore, we are going to establish the HIV model by using the FF operator in the AB sense (FF-HIV model). The proposed model is described as follows:

which is subject to

and, $\alpha$ and $\beta$ represent the fractional and fractal orders, respectively. ${_{t}^{\mathbb{FFM}}}\mathfrak{D}_{a}^{\alpha, \beta}$ is the FF derivative with the fractional order $\alpha \in (0, 1]$ and the fractal order $\beta \in (0, 1]$ with the generalized ML kernel.

3.   Model analysis

3.1.   Existence and positivity of the solution

This subsection analyzes the existence and positivity of the $\mathbb{FF}$-$\mathbb{HIV}$ model (2.5).

Theorem 3.1. There is a unique positive solution for the FF-HIV model $(2.5)$ that remains in $\mathbb{R}_{+}^{4}$. Moreover, the solution is non-negative.

Proof. In the FF-HIV model (2.5), we obtain its existence and uniqueness on the time interval $(0, \infty)$. Next, we are going to show that the non-negative region $\mathbb{R}_{+}^{4}$ is a positivity invariant region. From the FF-HIV model (2.5), we get

If $(\mathcal{T}(0), \mathcal{I}(0), \mathcal{V}(0), \mathcal{L}(0)) \in \mathbb{R}_{+}^{4}$, then, according to (3.1), the solutions cannot escape from the hyperplane.

Since

where $\phi_{\mu. \sigma} = \min\left\{ \mu, \sigma \right\}$, solving (3.2) by using the fractal Laplace transform (2.4), we obtain

Without loss of generality, we can suppose that

Then, we have the following biologically feasible region of the FF-HIV model (2.5):

with respect to the FF-HIV model (2.5), where

The proof is completed.

3.2.   EPs and BRN

This subsection describes how to find the EPs of the FF-HIV model (2.5). We are going to determine the EPs and the model's BRN ($\mathfrak{R}_{0}$). In the FF-HIV model (2.5), there are two types of possible EPs. The first point is the point at which there is no disease in the group, which is called the disease-free EP ($\mathfrak{E}_{0}^{*}$). In the method of obtaining $\mathfrak{E}_{0}^{*}$, we will be determining the right hand side of the FF-HIV model (2.5) when equal to zero, with $\mathcal{I}(t) = 0$, $\mathcal{V}(t) = 0$ and $\mathcal{L}(t) = 0$. Then,

To analyze the stability of $\mathfrak{E}_{0}^{*}$, we will focus on $\mathfrak{R}_{0}$, which can be obtained via the next-generation matrix (NGM) method ^[43,44] for the FF-HIV model (2.5). To solve $\mathfrak{R}_{0}$, we focus solely on the infectious classes of the FF-HIV model (2.5), i.e., $\mathcal{I}$, $\mathcal{V}$ and $\mathcal{L}$. The transmission matrix $\mathcal{F}$ and the transition matrix $\mathcal{V}$ are acquired as follows:

Hence, the NGM is provided by

Then, the spectral radius of the NGM (3.5) provides $\mathfrak{R}_{0} = \rho(\mathcal{F}\mathcal{V}^{-1})$, where

and $\rho$ represents the spectral radius.

Theorem 3.2. The FF-HIV model $(2.5)$ at the disease-free EP ($\mathfrak{E}_{0}$) is locally asymptotically stable (LAS) whenever $\mathfrak{R}_{0} < 1$, with the necessary and sufficient conditions:

Proof. In order to investigate the sufficient conditions for the FF-HIV model (2.5), we must obtain the eigenvalues $\theta_{i}$ of the Jacobian matrix ($\mathfrak{J}(\mathfrak{E}_{0})$) by computing the FF-HIV model (2.5) at $\mathfrak{E}_{0}$, which gives

Then, we obtain the characteristic equation as follows:

where

The coefficients given by $\omega_{i}$ ($i = 0, 1, 2, 3$) are positive. Then, we have the eigenvalue $\theta = - \mu$, which provides negative real parts; the others can be obtained by solving

where

Obviously, the coefficients $\epsilon_1$ and $\epsilon_2$ are positive, while $\epsilon_0$ is positive whenever $\mathfrak{R}_{0} < 1$, which corresponds to the assumption. Since

we have that

Thus, by the Routh-Hurwitz conditions given in ^[45], if $\epsilon_2 > 0, \epsilon_0 > 0$ and $\epsilon_1 \epsilon_2 > \epsilon_0$, we can conclude that all roots of (3.9) have negative real parts. This assures the assumption of (3.7) for all $\alpha \in (0, 1]$. Therefore, the model (2.5) at the steady-state $\mathfrak{E}_{0}$ is LAS if $\mathfrak{R}_{0} < 1$.

As we know, the value of $\mathfrak{R}_{0}$ provides the knowledge needed to estimate the transmission potential of an infectious disease over time. The FF-HIV model (2.5) has an endemic EP ($\mathfrak{E}_{1}^{*}$) when $\mathfrak{R}_{0} > 1$. To solve $\mathfrak{E}_{1}^{*}$, we will be applying the fact that all unknown variables $\mathcal{T}(t)$, $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$ of the FF-HIV model (2.5) are non-negative. It is possible to compute it by equating each equation of the FF-HIV model (2.5) to zero, that is,

which yields that $\mathfrak{E}_{1}^{*} = (\mathcal{T}_{1}^{*}, \mathcal{I}_{1}^{*}, \mathcal{V}_{1}^{*}, \mathcal{L}_{1}^{*})$, where the components of $\mathfrak{E}_{1}^{*}$ are given by

Theorem 3.3. The FF-HIV model $(2.5)$ at the endemic EP ($\mathfrak{E}_{1}^{*}$) is LAS whenever $\mathfrak{R}_{0} > 1$, with the necessary and sufficient conditions given by $(3.7)$ and

where

Proof. The Jacobian matrix computed at $\mathfrak{E}_{1}^{*}$ is provided as follows:

Then, the characteristic equation of $\mathfrak{J}(\mathfrak{E}_{1}^{*})$ is

We can see that $\overline\omega_{1}$ and $\overline\omega_{3}$ are positive, while $\overline\omega_{0}$ is positive if $\mathfrak{R}_{0} > 1$. We utilize the proposition given in ^[46], which applies the Routh-Hurwitz criterion that, if the given condition (3.12) is satisfied, then all real eigenvalues and all real parts of complex conjugate eigenvalues of (3.13) are negative, which ensures the condition (3.7) for all $\alpha \in (0, 1]$. Therefore, the point $\mathfrak{E}_{1}^{*}$ is LAS.

4.   Existence properties for the FF-HIV model (2.5)

In this section, we establish the uniqueness, existence and stability properties with fixed-point theorems. To demonstrate the existence and unique results of solutions for the FF-HIV model (2.5), we define a Banach space $\mathcal{B} = \mathcal{C}(\mathcal{J} \times \mathbb{R}^{4}, \mathbb{R})$ equipped with the norm

Since the integral is differentiable, the FF-HIV model (2.5) can be rewritten as

where

For the sake of discussion, we can rewrite the system (4.1) to have the following form:

where

and

By replacing ${_{t}^{\mathbb{ABR}}}\mathfrak{D}_{a}^{\alpha, \beta}$ with ${_{t}^{\mathbb{ABC}}}\mathfrak{D}_{a}^{\alpha, \beta}$ and applying the AB-fractional integral, we obtain

Define an operator $\mathcal{Q}: \mathcal{B} \to \mathcal{B}$ as

Clearly, the problem (4.3) has the solution if and only if $\mathcal{Q}$ has the fixed points.

4.1.   Uniqueness property

Theorem 4.1. Suppose that $\mathbb{H} \in \mathcal{C}(\mathcal{J} \times \mathbb{R}^{4}, \mathbb{R})$ satisfies the following condition:

($\mathcal{A}_{1}$) There is a number $\mathcal{L}_{\mathbb{H}} > 0$ such that

for all $\mathbb{U}_{1}$, $\mathbb{U}_{2} \in \mathcal{B}$ and $t \in \mathcal{J}$. If

then the problem $(4.3)$ has one solution, which signifies that the FF-HIV model $(2.5)$ has a unique solution.

Proof. First, we will tranform (4.3), corresponding to (2.5), into $\mathbb{U} = \mathcal{Q} \mathbb{U}$ (fixed-point problem).

Let $\mathbb{H}^{*}$ be a non-negative number so that $\sup_{t\in\mathcal{J}}\vert \mathbb{H}(t, 0) \vert = \mathbb{H}^{*} < +\infty$. Set a bounded, closed and convex subset $\mathcal{D}_{r_{1}} = \left\{ \mathbb{U} \in \mathcal{B} : \Vert \mathbb{U} \Vert \leq r_{1} \right\}$, where $r_{\mathbb{U}}$ is chosen so that

Next, the proof is separated into two steps.

Step 1. We present $\mathcal{Q}\mathcal{D}_{r_{1}} \subset \mathcal{D}_{r_{1}}$.

For any $\mathbb{U} \in \mathcal{D}_{r_{1}}$, we obtain

which yields that $\mathcal{Q} \mathcal{D}_{r_{1}} \subset \mathcal{D}_{r_{1}}$.

Step 2. We present that $\mathcal{Q}$ is a contraction.

For each $\mathbb{U}_{1}$, $\mathbb{U}_{2} \in \mathcal{D}_{r_{1}}$ and $t \in \mathcal{J}$, we obtain

where $t_{\min} = \min\{ t \in \mathcal{J} \}$, which yields that

Since $([(1-\alpha)\beta t_{\min}^{\beta-1}]/\mathbb{AB}(\alpha) + [\alpha\Gamma(\beta+1) T^{\alpha-\beta+1}]/[\mathbb{AB}(\alpha) \Gamma(\alpha + \beta)]) \mathcal{L}_{\mathbb{H}} < 1$, by the summary of the Banach contraction principle (defined as in ^[47]), then $\mathcal{Q}$ is called a contraction. Therefore, $\mathcal{Q}$ has the unique fixed-point, which is the unique solution of the FF-HIV model (2.5).

4.2.   Existence property

Theorem 4.2. Suppose that $\mathbb{H} \in \mathcal{C}(\mathcal{J}\times\mathbb{R}^{4}, \mathbb{R})$ satisfies the following condition:

($\mathcal{A}_{2}$) There are non-decreasing functions $\mathcal{G}\in\mathcal{C}(\mathbb{R}^+, \mathbb{R}^+)$ and $q \in \mathcal{C}(\mathcal{J}, \mathbb{R}^{+})$, such that

with $q^{*} = \sup_{t\in \mathcal{J}}\{q(t)\}$.

$(\mathcal{A}_{3})$ There is a number $\mathcal{M}^* > 0$ such that

Then, the problem $(4.3)$, which is consistent with the FF-HIV model $(2.5)$, has at least one solution $\mathbb{U} \in \mathcal{B}$.

Proof. Suppose that $\mathcal{Q}$ is defined by (4.7). In the first procedure, we show that $\mathcal{Q}$ maps bounded sets (balls) into bounded sets in $\mathcal{B}$. For any number $r_2 > 0$, set $\mathcal{D}_{r_2} : = \{\mathbb{U} \in \mathcal{B} : \Vert \mathbb{U} \Vert \leq r_2\}$ is a bounded set (ball) in $\mathcal{B}$. From ($\mathcal{A}_{2}$), for $t \in \mathcal{J}$, we have

where $t_{\min} = \min\{ t \in \mathcal{J}\}$. It follows that

Now, we present that $\mathcal{Q}$ maps bounded sets into equicontinuous sets of $\mathcal{B}$. Let $t_{1}$, $t_{2} \in \mathcal{J}$ given $t_{1} < t_{2}$, and for any $\mathbb{U} \in \mathcal{D}_{r_{2}}$. Hence, we obtain that

Obviously, the above result is independent of $\mathbb{U} \in \mathcal{D}_{r_{2}}$ when $t_{2} \to t_{1}$; then, $\left\vert (\mathcal{Q}\mathbb{U})(t_{2}) - (\mathcal{Q}\mathbb{U})(t_{1}) \right\vert \to 0$. Consequently, $\left\Vert (\mathcal{Q}\mathbb{U})(t_{2}) - (\mathcal{Q}\mathbb{U})(t_{1}) \right\Vert \to 0$ as $t_{2} \to t_{1}$. Thus, $\mathcal{Q}$ is equicontinuous. So, by the Arzelá-Ascoli theorem (defined as in ^[48]), $\mathcal{Q}$ is completely continuous.

Finally, we present that there exists an open set $\mathbb{D} \subset \mathcal{B}$ with $\mathbb{U} \neq \upsilon \mathcal{Q} \mathbb{U}$ for $0 < \upsilon < 1$ and $\mathbb{U} \in \partial \mathbb{D}$. Assume that $\mathbb{U} \in \mathcal{B}$ is a solution of $\mathbb{U} = \upsilon \mathcal{Q} \mathbb{U}$ for each $0 < \upsilon < 1$. Then, for any $t \in \mathcal{J}$, we will present that $\mathcal{Q}$ is bounded; thus, we have that

Taking the norm of both sides of (4.9) for any $t \in \mathcal{J}$, yields that

In view of ($\mathcal{A}_{3}$), there is $\mathcal{M}^{*} > 0$ such that $\Vert \mathbb{U} \Vert \neq \mathcal{M}^{*}$. Define $\mathbb{D} : = \{ \mathbb{U} \in \mathcal{B} : \Vert \mathbb{U} \Vert < \mathcal{M}^{*} + 1 \}$ and $\mathcal{E} = \mathbb{D} \cup \mathcal{D}_{r_{2}}$. Clearly, $\mathcal{Q}: \overline{\mathcal{E}} \to \mathcal{B}$ is continuous and completely continuous. By the choice of $\mathbb{D}$, there is no $\mathbb{U} \in \partial \mathbb{D}$ such that $\mathbb{U} = \upsilon \mathcal{Q} \mathbb{U}$ for some $0 < \upsilon < 1$.

Hence, by using a nonlinear alternative of the Leray-Schauder type (defined as in ^[47]), we can conclude that $\mathcal{Q}$ has the fixed-point $\mathbb{U} \in \overline{\mathbb{D}}$, which suggests that the FF-HIV model (2.5) has at least one solution $\mathbb{U} \in \mathcal{B}$. The proof is completed.

5.   Ulam stability for the FF-HIV model (2.5)

In this section, we are going to prove different types of Ulam's stability for the FF-HIV model (2.5). First, we will provide the definitions of Ulam's stability that will be used in this section. Assume that we have a real number $\epsilon \in \mathbb{R}^{+}$ and a function $\mathcal{K}_{\mathbb{H}} \in \mathcal{C}(\mathcal{J}, \mathbb{R}^{+})$.

Definition 5.1. The FF-HIV model $(2.5)$ is called UH-stable if there is a constant $\mathcal{R}_{\mathbb{H}} > 0$ so that, for every $\epsilon > 0$ and any solution $\mathbb{Z} \in \mathcal{B}$ of

there exists the solution $\mathbb{U} \in \mathcal{B}$ of the FF-HIV model $(2.5)$ with

where $\epsilon = \max(\epsilon_{i})^{T}$ and $\mathcal{R}_{\mathbb{H}} = \max(\mathcal{R}_{\mathbb{H}_{i}})^{T}$, $i = 1, 2, 3, 4$.

Definition 5.2. The FF-HIV model $(2.5)$ is called GUH-stable if there is a function $\mathcal{K}_{\mathbb{H}} \in \mathcal{C}(\mathbb{R}^{+}, \mathbb{R}^{+})$ given $\mathcal{K}_{\mathbb{H}}(0) = 0$ so that, for every solution $\mathbb{Z} \in \mathcal{B}$ of

there exists the solution $\mathbb{U} \in \mathcal{B}$ of the FF-HIV model $(2.5)$ so that

where $\epsilon = \max(\epsilon_{i})^{T}$ and $\mathcal{R}_{\mathbb{H}} = \max(\mathcal{K}_{\mathbb{H}_{i}})^{T}$, $i = 1, 2, 3, 4$.

Definition 5.3. The FF-HIV model $(2.5)$ is called RUH-stable with respect to $\mathcal{K}_{\mathbb{H}}\in\mathcal{C}(\mathcal{J}, \mathbb{R}^{+})$ if there is a constant $\mathcal{R}_{\mathcal{K}_{\mathbb{H}}} > 0$ so that, for every $\epsilon > 0$ and any solution $\mathbb{Z} \in \mathcal{B}$ of $(5.3)$, there is a solution $\mathbb{U} \in \mathcal{B}$ of the FF-HIV model $(2.5)$ with

where $\epsilon = \max(\epsilon_{i})^{T}$, $\mathcal{R}_{\mathcal{K}_{\mathbb{H}}} = \max(\mathcal{R}_{\mathcal{K}_{\mathbb{H}_{i}}})^{T}$ and $\mathcal{K}_{\mathbb{H}} = \max(\mathcal{K}_{\mathbb{H}_{i}})^{T}$, $i = 1, 2, 3, 4$.

Definition 5.4. The FF-HIV model $(2.5)$ is called GUHR-stable with respect to $\mathcal{K}_{\mathbb{H}}(\mathcal{J}, \mathbb{R}^{+})$ if there is a constant $\mathcal{R}_{\mathcal{K}_{\mathbb{H}}} > 0$ so that, for every solution $\mathbb{Z} \in \mathcal{B}$ of

there is the solution $\mathbb{U} \in \mathcal{B}$ of the FF-HIV model $(2.5)$ with

where $\mathcal{R}_{\mathcal{K}_{\mathbb{H}}} = \max(\mathcal{R}_{\mathcal{K}_{\mathbb{H}_{i}}})^{T}$ and $\mathcal{K}_{\mathbb{H}} = \max(\mathcal{K}_{\mathbb{H}_{i}})^{T}$, $i = 1, 2, 3, 4$.

Remark 5.5. It is obvious that

(i) Definition 5.1 implies Definition 5.2.

(ii) Definition 5.3 implies Definition 5.4.

(iii) Definition 5.3 with $\mathcal{K}_{\mathbb{H}}(t) = 1$ implies Definition 5.1.

Remark 5.6. $\mathbb{Z} \in \mathcal{B}$ is the solution of $(5.1)$ if and only if there is $\mathbb{W} \in \mathcal{B}$ (which depends on $\mathbb{Z}$) so that the following properties apply:

($\mathcal{H}_{1}$) $\left\vert \mathbb{W}(t) \right\vert \leq \epsilon$ with $\mathbb{W} = \max(\mathbb{W}_{i})^{T}$ for $i = 1, 2, 3, 4$ and $t \in \mathcal{J}$.

($\mathcal{H}_{2}$) ${_{t}^{\mathbb{ABR}}}\mathfrak{D}_{a}^{\alpha, \beta} \mathbb{Z}(t) = \mathbb{H}(t, \mathbb{Z}) + \mathbb{W}(t)$ for all $t \in \mathcal{J}$.

Remark 5.7. $\mathbb{Z} \in \mathcal{B}$ is the solution of $(5.3)$ if and only if there is $\mathbb{V} \in \mathcal{B}$ (which depends on $\mathbb{Z}$) so that the following properties apply:

($\mathcal{H}_{3}$) $\left\vert \mathbb{V}(t) \right\vert \leq \epsilon \mathcal{K}_{\mathbb{H}}(t)$ with $\mathbb{W} = \max(\mathbb{W}_{i})^{T}$ and $\mathcal{K}_{\mathbb{H}} = (\mathcal{K}_{\mathbb{H}_{i}})^{T}$ for $i = 1, 2, 3, 4$ and $t \in \mathcal{J}$.

($\mathcal{H}_{4}$) ${_{t}^{\mathbb{ABR}}}\mathfrak{D}_{a}^{\alpha, \beta} \mathbb{Z}(t) = \mathbb{H}(t, \mathbb{Z}) + \mathbb{V}(t)$ for all $t \in \mathcal{J}$.

5.1.   UH and GUH stability properties

Lemma 5.8. Assume that $\alpha \in (0, 1]$ and $\beta \in (0, 1]$. If $\mathbb{Z} \in \mathcal{B}$ is the solution of $(5.1)$ then, $\mathbb{Z}$ is the solution of

where $\mathcal{G}_{\mathbb{Z}}(t) = \mathbb{Z}(0) + \frac{(1-\alpha)\beta}{\mathbb{AB}(\alpha)} t^{\beta-1} \mathbb{H}(t, \mathbb{Z}(t))$.

Proof. Assume that $\mathbb{Z}$ is the solution of (5.1). By applying ($\mathcal{H}_{2}$) in Remark 5.6, we get

Hence, the solution of (5.9) is given by

By using ($\mathcal{H}_{1}$) in Remark 5.6, we have that

Hence, the property (5.8) is achieved.

Theorem 5.9. Assume $\mathbb{H}\in \mathcal{C}(\mathcal{J} \times \mathbb{R}, \mathbb{R})$ for any $\mathbb{U} \in \mathcal{B}$. If ($\mathcal{A}_{1}$) and $(4.8)$ are satisfied, then the FF-HIV model $(2.5)$ is UH-stable on $\mathcal{J}$.

Proof. Let $\epsilon > 0$ be a real number and $\mathbb{Z} \in \mathcal{B}$ be the solution of (5.1). Suppose that $\mathbb{U} \in \mathcal{B}$ is the unique solution of the problem

where $\mathbb{U}(t)$ is provided by

when $\mathcal{G}_{\mathbb{U}}(t) = \mathbb{U}(0) + \frac{(1-\alpha)\beta}{\mathbb{AB}(\alpha)} t^{\beta-1} \mathbb{H}(t, \mathbb{U}(t)).$

From Lemma 5.8 with ($\mathcal{A}_{1}$) in Theorem 4.1, we have that

Setting

yields that

Therefore, we can conclude that the FF-HIV model (2.5) is UH-stable.

Corollary 5.10. Setting $\mathcal{K}_{\mathbb{H}}(\epsilon) = \mathcal{R}_{\mathbb{H}} \epsilon$ with $\mathcal{K}_{\mathbb{H}}(0) = 0$ in Theorem 5.9, implies that the FF-HIV model $(2.5)$ is GUH-stable.

5.2.   UHR and GUHR stability properties

Next, we will prove the UHR and GUHR stability results. Assume the following:

($\mathcal{H}_{5}$) There exist an increasing function $\mathcal{K}_{\mathbb{H}} \in \mathcal{B}$ and a number $\mathfrak{C}_{\mathcal{K}_{\mathbb{H}}} > 0$ so that

Lemma 5.11. Assume that $\alpha \in (0, 1]$ and $\beta \in (0, 1]$. If $\mathbb{Z} \in \mathcal{B}$ is the solution of (5.3), then $\mathbb{Z}$ is the solution of

where $\mathcal{G}_{\mathbb{Z}}(t) = \mathbb{Z}(0) + \frac{(1-\alpha)\beta}{\mathbb{AB}(\alpha)} t^{\beta-1} \mathbb{H}(t, \mathbb{Z}(t))$.

Proof. Let $\mathbb{Z}$ be a solution of (5.3). By applying ($\mathcal{H}_{4}$) in Remark 5.7, we obtain the following:

Hence, the solution of (5.13) is given by

By using ($\mathcal{H}_{3}$) in Remark 5.7, we have that

Then, the property (5.12) is achieved.

Theorem 5.12. Assume $\mathbb{H}\in \mathcal{C}(\mathcal{J} \times \mathbb{R}, \mathbb{R})$ for every $\mathbb{U} \in \mathcal{B}$. If ($\mathcal{A}_{1}$), ($\mathcal{H}_{5}$) and $(4.8)$ are satisfied, then the FF-HIV model $(2.5)$ is UHR-stable on $\mathcal{J}$.

Proof. Let $\epsilon > 0$ be a real constant and $\mathbb{Z} \in \mathcal{B}$ be the solution of (5.6). Suppose that $\mathbb{U} \in \mathcal{B}$ is the unique solution of

where $\mathbb{U}(t)$ is provided by

when $\mathcal{G}_{\mathbb{U}}(t) = \mathbb{U}(0) + \frac{(1-\alpha)\beta}{\mathbb{AB}(\alpha)} t^{\beta-1} \mathbb{H}(t, \mathbb{U}(t)).$

From Lemma 5.11 with ($\mathcal{A}_{1}$) in Theorem 4.1 and ($\mathcal{H}_{5}$), we obtain that

Setting

yields that

Hence, we can conclude that the FF-HIV model (2.5) is UHR-stable.

Corollary 5.13. Setting $\epsilon = 1$ in $\left\vert \mathbb{Z}(t) - \mathbb{U}(t) \right\vert \leq \mathcal{R}_{\mathcal{K}_{\mathbb{H}}} \epsilon \mathcal{K}_{\mathbb{H}}(t)$, implies that the FF-HIV model $(2.5)$ is GUHR-stable.

6.   Numerical schemes with the generalized ML kernel for the FF-HIV model (2.5)

In this section, we apply three powerful algorithms for the FF-HIV model (2.5) by implementing the Newton polynomial approach, the Adams-Bashforth method and the predictor-corrector method. The numerical simulations are described and presented, as well as a comparison of the three methods, which will be handled later in the subsection on dynamical discussion.

In order to obtain the numerical schemes, we recall (4.3) and (4.6), which was obtained after integrating (4.3) on both sides. We now have

The following result was obtained by using the AB-fractional integral; we do not show $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$ because they contain the same steps.

We rewrite (6.2)–(6.5) at $t = t_{n+1} = (n+1) \Delta t$, which gives

6.1.   Numerical method for the Newton polynomial technique

With an effective idea of an algorithm that has been studied in ^[49], a numerical approximation form for solving the FF-HIV model (4.3) using the Newton approach will be considered.

By using the approximation of the integrals in (6.6), we have that

Applying the Newton polynomial for the term $t^{\beta-1} \mathbb{U}_{1}(t, \mathcal{T}, \mathcal{I}, \mathcal{V}, \mathcal{L})$, yields

Plugging (6.11) into (6.10), we get

Integration by parts to solve the integrals given in (6.12) yields

where the constants $\Lambda_{n, j}$, $\Xi_{n, j}$ and $\Omega_{n, j}$ are given by

Plugging (6.13)–(6.15) into (6.12), we obtain the numerical scheme as follows:

Clearly, we also have that

where $\Lambda_{n, j}$, $\Xi_{n, j}$ and $\Omega_{n, j}$ are provided by (6.16)-(6.18), respectively.

6.2.   Numerical method for the Adams-Bashforth technique

This subsection shows the numerical algorithm for the FF-HIV model (2.5), which was developed by applying the Adams-Bashforth method, which is a well-known technique based on two-step Lagrange polynomials. By using (6.6)–(6.9); then, the approximation of (6.6)–(6.9) are formulated in the form of

Now, we estimate the three functions $s^{\beta-1} \mathbb{U}_{i}(s, \mathcal{T}, \mathcal{I}, \mathcal{V}, \mathcal{L})$, $i = 1, 2, 3$, on $[t_{j}, t_{j+1}]$ by using two-step Lagrange interpolation polynomials with the step size $h = t_{j+1} - t_{j}$. By direct computations, we develop the algorithms that yield the numerical simulations for the model FF-HIV (2.5) as follows:

Repeating all steps for $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$, we also have that

where

and

6.3.   Numerical method for the predictor-corrector technique

This subsection shows the numerical algorithm for the FF-HIV model (2.5), which was developed by applying the predictor-corrector method. By using the approximation of the integrals as in (6.6)–(6.9), we have that

Now, we apply the Adams-type predictor-corrector tool to obtain the numerical approximation of the right-hand sides of (6.27)–(6.30). Assume that the solution belongs to $[0, T]$; then, we set $\Delta t = T/N$, $t_{k} = k\Delta t$ for $k = 0, 1, 2, \ldots, N$. So, the corrector scheme for the FF-HIV model (2.5) are given as follows:

Repeating all steps for $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$, we also have that

where

Further, the predictor terms $\mathcal{T}_{n+1}^{P}$, $\mathcal{I}_{n+1}^{P}$, $\mathcal{V}_{n+1}^{P}$ and $\mathcal{L}_{n+1}^{P}$ are described as

and

6.4.   Numerical simulations and discussion

Employing the numerical schemes in the subsections 6.1–6.3, we conducted the numerical simulations for the FF-HIV model (2.5) by using the real data for the parameter values as assumed in ^[13], as shown in Table 3. Here, the approximate solutions of the model (2.5) have been verified for various fractional-order values of $\alpha = \{ 1.000, 0.995, 0.985, 0.975, 0.965, 0.955, 0.945, 0.935 \}$ and the fractal dimensions $\beta = \{ 1.000, 0.995, 0.985, 0.975, 0.965, 0.955, 0.945, 0.935 \}$, which have been separated into three cases.

Note that we present the pairs of numerical simulations as a comparison of the Newton polynomial, the Adams-Bashforth and the predictor-corrector techniques. The quantity of time $t$ in these graphs denotes the number of days.

Case (i): If we set the fractal dimension $\beta$ and the {fractional-order} $\alpha$ values as $1.000$, $0.995$, $0.985$, $0.975$, $0.965$, $0.955$, $0.945$ and $0.935$ with $(\mathcal{T}(0), \mathcal{I}(0), \mathcal{V}(0), \mathcal{L}(0)) = (200, 25, 25, 25)$. Figure 2 shows the dynamics of the density of the susceptible group of CD$4^+$ T-cells ($\mathcal{T}(t)$), the group of infected CD$4^+$ T-cells before ($\mathcal{I}(t)$) and after the completion of the reverse transcription ($\mathcal{V}(t)$), and the group of the virus as time passes ($\mathcal{L}(t)$), respectively, under the conditions of the Newton polynomial technique in Subsection 6.1 (Figure 2a), the Adams-Bashforth technique in Subsection 6.2 (Figure 2b) and the predictor-corrector technique in Subsection 6.3 (Figure 2c). Figure 3a–c represent that, when the fractional-order $\alpha$ values increase to integer order and the fractal dimension $\beta = 1.000$, the behavior of the susceptible group of CD$4^+$ T-cells rapidly increases from the start day, and then steadily decreases. After that, it increases again and finally converges to $\mathcal{T}_1^* \approx 290.6250$ at the end of the simulation period. Figures 4 and 5 show that, when the fractional-order $\alpha$ values increase from $0.935$ to the integer order and the fractal dimension $\beta = 0.987$, the behavior of the group of infected CD$4+$ T-cells before and after the completion of the reverse transcription rapidly increases to the peak values and steadily decreases, respectively, and also converges to $\mathcal{I}_1^* \approx 40.5357$ and $\mathcal{V}_{1}^{*} \approx 24.9451$ according to their asymptotic stabilities at the end. Figure 6 shows that, when the fractional-order $\alpha$ values increase from $0.935$ to the integer order and the fractal dimension $\beta = 1.000$, the behavior of the group of the virus as time passes rapidly increases and decreases. After that, it rapidly increases to the peak point and steadily decreases to converge to $\mathcal{L}_1^* \approx 2702.3810$ at the end. As shown in Figures 3–6, by using three numerical algorithms, the dynamic behavior of each of the four variables varies in sightly different ways when the fractional order is slightly changed. Their behaviors oscillate decreasing and increasing until, tending toward stabilization.

Case (ii): In this case, we fix $\alpha$ and vary $\beta$ as $1.000$, $0.995$, $0.985$, $0.975$, $0.965$, $0.955$, $0.945$ and $0.935$ with $(\mathcal{T}(0), \mathcal{I}(0), \mathcal{V}(0), \mathcal{L}(0)) = (200, 25, 25, 25)$. Figure 7 shows the dynamics of the {densities} of $\mathcal{T}(t)$, $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$ under the conditions of the Newton polynomial technique (Figure 7a), the Adams-Bashforth technique (Figure 7b) and the predictor-corrector technique in Subsection 6.3 (Figure 7c). Figure 8a-c show that, when the fractional-order $\alpha = 0.987$ and the fractal dimension $\beta$ increases to integer order, the behavior of $\mathcal{T}(t)$ rapidly increases from the start day, and then steadily decreases tending to $\mathcal{T}_1^* \approx 290.6250$. Figures 9 and 10 show, that when the fractional-order $\alpha = 0.987$ and the fractal dimension $\beta$ increases from $0.935$ to the integer order, $\mathcal{I}(t)$ and $\mathcal{V}(t)$ rapidly decrease and increase to the peak values, respectively. Further, they steadily decrease again from the peak values to $\mathcal{I}_1^* \approx 40.5357$ and $\mathcal{V}_{1}^{*} \approx 24.9451$ according to their asymptotic stabilities. Figure 11 indicates that, when the fractional-order $\alpha = 0.987$ and the fractal dimension $\beta$ increases from $0.935$ to the integer order, $\mathcal{L}(t)$ quickly increases and then decreases. After that, it increases again to the peak point and oscillates as it decreases to $\mathcal{L}_1^* \approx 2702.3810$. The dynamics of the four variables, as illustrated in Figure 8 to Figure 11, exhibit behavior that is similar to Case (i). Additionally, after using the three numerical techniques, we discovered that, when the fractional order has tiny variation, the system behaves in a somewhat different manner.

Case (iii). If we vary both the fractional-order $\alpha$ and fractal dimension $\beta$, each as $1.000$, $0.995$, $0.985$, $0.975$, $0.965$, $0.955$, $0.945$ and $0.935$ with $(\mathcal{T}(0), \mathcal{I}(0), \mathcal{V}(0), \mathcal{L}(0)) = (1000, 50, 70, 80)$, which differ from the other two cases. Observing Figures 12–16, one of the noticeable aspects of the behavior of the system is that, even if we start with the other initial condition, the system will converge to the steady state.

In order to compare the results of the three numerical algorithms, i.e., the Newton polynomial, the Adams-Bashforth and the predictor-corrector techniques for all four state functions $\mathcal{T}(t)$, $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$, we illustrate the extra case with the step size $\Delta t = 0.1$ and $\alpha = \beta = 0.945$, as seen in Figures 17a–d. They display the results of the numerical simulations, which compare the numerical solutions of the three algorithms for the concentration of susceptible CD$4^+$ T-cells, infected CD$4^+$ T-cells before and after reverse transcription completion and the virus density. We show some numerical results for the population in four states given $\alpha = \beta = 0.945$ and $\Delta t = 0.1$, which can be seen in Tables 4–7. These graphical simulation and numerical results reveal that these three numerical algorithms produce similar outcomes with minor differences.

Moreover, given that $\alpha$ and $\beta$ varies in each case, we observe that the fractional dimension $\beta$ is independent of the EPs. Therefore, we obtained the values of the parameters $\mathfrak{R}_0 \approx 3.4409$, $\overline{\omega}_0 \approx 0.0071$, $\overline{\omega}_1 \approx 0.0893$, $\overline{\omega}_2 \approx 2.0760$ and $\overline{\omega}_3 \approx 3.6649$ and detected that all of these parameters satisfy the conditions of Theorem 3.3 i.e., that $\mathfrak{R}_0 > 1$ and $\overline\omega_{1}\overline\omega_{2}\overline\omega_{3} > \overline\omega_{1}^2+\overline\omega_{3}^2 \overline\omega_{0}$. This guarantees locally asymptotic stability at the endemic EP $\mathfrak{E}_{1}^{*} \approx (290.6250, 40.5357, 24.9451, 2702.3810)$.

For our discussion, it is clearly visible from all of the images for Case (i)–(iii) that the approximate solutions for the different fractional orders converge to their solution for an integer-order model. The $\mathcal T(t)$, $\mathcal{I}(t)$, $\mathcal{V}(t)$ and $\mathcal{L}(t)$ densities oscillate initially before tending to stabilize at their EPs. Furthermore, the three numerical algorithms gave similar numerical results with a very tiny difference.

7.   Conclusions

As we know, differential equations with a fractional order can be used to help overcome some of the flaws in biological systems that are connected to the idea of memory additionally, the use of differential equations in mathematical modeling helps us to better comprehend the behavior of dynamic biological systems. These are good reasons for employing the FF derivative in the AB sense for the model of the interaction between HIV infection and CD$4^+$ T-cells in the presence of RTI that has been taken before the process of producing the virus was initiated. The significant outcomes and the parade of theoretical results were used to accomplish the requirements of our aim. First, the model was analyzed by deriving the unique positive solution and finding the sufficient conditions for it to be LAS at the EPs by using the BRN and the Routh-Hurwitz criteria. Second, the existence and uniqueness results are proven with the help of famous fixed-point theorems such as Banach's and Leray-Schauder's. Third, in order to emphasize the existence of solutions, various forms of Ulam's stability, such as UH, GUH, UHR and GUHR forms were used to achieve the purpose. Finally, a discussion of the dynamical behavior has been provided to confirm the accuracy of the theoretical results. The numerical simulations were distinguished for three different cases by using the Newton polynomial, the Adams-Basthforth and the predictor-corrector techniques in the numerical scheme. Overall, the fractional operator utilized in the study according to the findings satisfied all of the required theoretical conditions for the considered model, and the parameters were shown to have a substantial influence on the ecological system's stability. A slight variation in the FF orders was found to cause a tiny change in the behavior of the considered model. The three numerical algorithms were nearly identical in terms of the numerical results and system dynamic behavior, with only slight variance between them.

In summary, we found that the significant advantage of our considered FF-order HIV epidemic model is that it is more realistic, effective and efficient than the classical model because it improves the precision of the model by allowing for more flexibility, which helps us to achieve better results. Regarding the future research direction, the FF derivative in the AB sense can be applied to various mathematical models for the investigation and analysis of real-world problems. Alternativelt, there are some interesting fractional-order derivatives, such as CF, AB, Atangana-Koca and Atangana-Gomez derivatives, which can be applied to the HIV integer-order model and are worthy of further study. We continue to hope that our efforts may help researchers in various areas of applied sciences and engineering.

Acknowledgment

J. Kongson, C. Thaiprayoon and A. Neamvonk would like to extend their appreciation to the Faculty of Science at Burapha University in Thailand for the financial support (Grant no. $SC06/2564$). W. Sudsutad was partially supported by Ramkhamhaeng University and J. Alzabut is thankful to Prince Sultan University and OSTİM Technical University for their endless support. This work was financially supported by the Faculty of Science at Burapha University in Thailand (Grant no. $SC06/2564$).

Conflict of interest

The authors declare that there is no conflict of interest.