1.
Introduction
The human immunodeficiency virus (HIV) continues to pose a serious global health challenge. Despite extensive research, no effective vaccine has been developed for this chronic, lifelong infection. In Taiwan, there is a particular gap in mathematical modeling tools capable of accurately estimating HIV transmission dynamics, especially among men who have sex with men (MSM), a population that bears the highest burden of new infections.
To address this, we incorporate the modified Atangana-Baleanu-Caputo (ABC) fractional operator into a new mathematical model that builds upon foundational work in [1]. This extension allows for deeper insights and more accurate approximations compared to existing models. Currently, there are no established fractional-order models for capturing the dynamics of HIV/AIDS among the Taiwanese MSM population. Therefore, we propose a modified ABC fractional-order model that categorizes individuals by their stage of infection and level of infectivity. This model aims to address critical gaps in HIV/AIDS estimation and provide a robust foundation for targeted intervention strategies.
HIV attacks the immune system, progressively leading to acquired immunodeficiency syndrome (AIDS) [2]. Infected individuals remain contagious for life, and without treatment, the infection is incurable [3]. While the overall incidence of HIV has declined since 2008, sexual transmission remains the primary mode of infection [4]. MSM face a significantly elevated risk of contracting HIV, making it crucial to analyze and predict trends in this group to guide public health responses [5]. Disease transmission within a closed population can be represented by compartmental transitions over time [6]. A widely used framework for modeling infectious disease dynamics is the susceptible-infected-removed (SIR) model [7,8], and numerous extensions have been proposed to accommodate different epidemiological features [9,10,11,12].
Hsieh and Lin [13] used a mathematical modeling approach to assess trends and project future outcomes for the HIV epidemic in Taiwan. While a molecular epidemiological study has also examined HIV transmission among MSM in Taiwan, sample selection bias limited its generalizability [14]. According to a 2012 report by the United Nations Joint Programme on HIV/AIDS (UNAIDS), 35.3 million individuals globally were living with HIV [15]. Recent assessments show rising HIV incidence among MSM, particularly in East Asia, Africa, and Russia [16]. Ko et al. [17] investigated HIV trends and risky sexual behaviors among gay bathhouse attendees in Taiwan from 2004 to 2008. Weber et al. [18] identified key behavioral and demographic risk factors linked to HIV infection among young gay and bisexual men in Canada. HIV transmission in this group is closely connected to several high-risk factors, including unprotected receptive or insertive anal sex [19,20], histories of sexually transmitted infections (STIs) [21], multiple concurrent or sequential partners [22,23], and sexual activity under the influence of alcohol or drugs [24,25].
Over the past decade, numerous mathematical models have been developed to study immune system responses to HIV infection [26,27,28], enhancing our understanding of the disease's biological and epidemiological behavior. The modified ABC fractional derivative model for HIV/AIDS dynamics was introduced in [29]. Recent advancements by Araz and Atangana [30] extended this framework further, transforming it into a piecewise modified ABC operator developed by Baleanu and Al-Refai [31]. This approach blends classical and fractional derivatives, enabling the analysis of crossover behavior in the transition between dynamic regimes over time. The ability to model such behavior significantly improves realism in nonlinear systems.
Beyond filling existing gaps in predictive modeling, the proposed hybrid fractional model provides a more accurate and comprehensive representation of HIV transmission dynamics. This research aims to guide data-driven policy decisions and contribute to the development of effective public health strategies aimed at reducing HIV incidence in high-risk communities. Specifically, it introduces a fractional-order hybrid model designed to enhance prediction accuracy and support more responsive interventions for the growing HIV epidemic among Taiwanese MSM.
This paper is organized as follows: Section 1 presents a detailed introduction to the proposed framework. Section 2 introduces the clinical context and discusses basic fractional-order derivatives relevant to the model. Section 3 outlines the generalized hypothesis and describes the modeling approach. Section 4 analyzes the well-posedness of the proposed system. Section 5 addresses core properties of the MABC-MSM model, including sensitivity analysis, positivity, existence, the Lipschitz condition, and uniqueness. Section 6 discusses the stability of the model, including Ulam-Hyers stability. Section 7 details the numerical strategy based on Lagrange interpolation. Finally, Sections 8 and 9 present the simulation results and concluding discussion, respectively.
2.
Essential principles of the fractional operator
This section presents the fundamental concepts of the modified Atangana-Baleanu-Caputo (MABC) fractional calculus, which form the mathematical foundation for the results developed in this study.
Definition 2.1. [32,33] Let f(t) be a continuous function on the interval (0,1). The fractal-fractional derivative of f(t) of order β and fractal dimension α, with a kernel of the modified ABC type, is defined as
It is straightforward to verify that MABCDβ,α0,tC=0 for any constant C, based on the definition provided in [31].
The next definition introduces the corresponding fractional integral associated with the MABC operator.
Definition 2.2. [32,33] Let f(t) be continuous on the interval (0,1). The fractal-fractional integral of f(t) of order β and dimension α, with a modified AB-type kernel, is defined as
where RLJ denotes the Riemann-Liouville fractional integral.
Lemma 2.1. [31] Let f(n)∈L1(0,∞), with n−1<β<n, n∈N, and β=α−n+1. Then, the following identity holds:
In the special case where 0<β<1, the identity simplifies to
3.
Mathematical fractal-fractional formulation of the MSM model
We propose a fractional-order SIAJB model, developed as an extension of the classical SIAJB framework, to describe the dynamics of HIV transmission among men who have sex with men (MSM) in Taiwan. By incorporating the modified ABC (MABC) fractional operator, the model captures more complex temporal behavior and provides a more accurate representation of HIV dynamics in high-risk populations.
Given that individuals exhibit varying levels of infectivity at different stages of infection, some eventually cease to be infectious. The progression of HIV within the MSM population under treatment interventions is described by the following system of fractal-fractional differential equations:
with initial conditions
3.1. Model description
The entire population N(t) is divided into five compartments at any given time t:
S(t): Susceptible individuals (not infected).
I(t): HIV-infected individuals who are contagious but have not developed AIDS.
J(t): Formerly infectious individuals who no longer engage in risky behavior.
A(t): Individuals who have progressed to AIDS and remain infectious.
B(t): Individuals with AIDS who have ceased high-risk behavior.
The model reflects treatment and behavioral transitions, particularly relevant for MSM undergoing therapy. Table 1 defines the parameters used in the model.
A schematic of the model dynamics is illustrated in Figure 1.
The total population is defined as
4.
Well-posedness of the framework
This section examines the conditions under which the model admits a well-defined and biologically meaningful solution. We consider the system within a finite time interval and a positively invariant region. From the model equations and their sum, we compute the dynamics of the total population N(t). Summing the equations in (4.1) yields
Given that all compartments are non-negative and δ>0, it follows that
by using a standard comparison principle. This inequality defines the threshold level for the total population. Hence, the solution is restricted to the positively invariant region
Here, R5+ denotes the positive cone in five dimensions, including boundary faces. For practical and epidemiological relevance, we exclude the case where MABCDβ,α0,tN(t)≥0 indefinitely, as this could imply an unrealistic scenario of unbounded population growth.
5.
Essential properties of the MABC-MSM system
To deepen our understanding of the mathematical model describing HIV/AIDS transmission among men who have sex with men (MSM), we now conduct a qualitative analysis of the system of nonlinear differential equations defined in (4.1). This analysis focuses on identifying the essential characteristics of the model and determining which parameters most significantly influence the dynamics of disease spread.
5.1. Sensitivity analysis
The behavior of the basic reproduction number, denoted R0, is central to understanding the transmission potential of the infection. As detailed in [2], Table 2 and Figures 2 and 3 for sensitivity analysis help identify which parameters most strongly affect R0, thereby highlighting the key drivers of disease dynamics. In the context of the MABC-MSM model, sensitivity analysis also offers insights into the role of memory and hereditary effects introduced via the fractional-order parameter.
The inclusion of fractional-order derivatives emphasizes the system's sensitivity to prior states, resulting in a slower decline in infection levels and greater responsiveness to intervention strategies. This validates the use of fractional modeling in representing long-term, history-dependent processes typical of chronic diseases such as HIV/AIDS.
The expression for the basic reproduction number is given by
To assess parameter influence, we compute the partial derivatives of R0 with respect to each model parameter:
Positive partial derivatives indicate that increasing the corresponding parameter will increase R0, while negative values suggest a reduction in transmission potential.
In conclusion, the sensitivity analysis indicates that controlling disease transmission is most effective through the management of key parameters, such as β1, β2, and g, which have the most significant influence on R0. The analysis also reinforces the utility of fractional-order models in capturing memory-dependent dynamics, offering a powerful tool for informing prevention strategies and optimizing intervention outcomes.
5.2. Positively invariant set Ω of the MABC-MSM system
Theorem 5.1. The region Ω associated with the MABC-MSM system (4.1) is positively invariant.
Proof. Given appropriate initial conditions, the solutions to the MABC-MSM system (4.1) remain non-negative for all t≥0. Therefore, once the solution enters the positive orthant of R5+, it remains there for all future time. This ensures that the population variables maintain biologically feasible (non-negative) values throughout the system's evolution.
We verify the invariance by evaluating the system at the boundary of each compartment:
Each equation confirms that the derivative is non-negative on the boundary, ensuring that the solution cannot exit the positive orthant. Thus, the region Ω⊂R5+ is positively invariant. □
5.3. Existence conditions for solutions
The existence conditions for the MABC-MSM system lend theoretical rigor to the model's predictions and ensure well-posedness. Using Lemma 2.1 and Definition 2.2, system (4.1) can be rewritten in its equivalent integral form as follows:
For S(t):
For I(t):
For J(t):
For A(t):
For B(t):
The functions Oi(t,x(t)), for i=1,…,5, are defined as
These expressions represent the dynamics of each compartment and are essential for analyzing the existence, uniqueness, and stability conditions of the MABC-MSM model.
5.4. Lipschitz condition
Lemma 5.1. Assume that the functions
are continuous. Then, there exist five positive constants yi, for 1≤i≤5 such that
Under these bounds, the nonlinear functions Oi defined in (5.8) satisfy Lipschitz continuity, with a global Lipschitz constant
where
Proof. Consider the function O1(t,S(t)). For S(t),˜S(t)∈L1[0,1], we have
Letting MO1=μy1+β1y2+β2y4, we rewrite as follows:
Similarly, for the remaining functions with upper bounds,
Since MO=max1≤i≤5MOi<∞, each Oi satisfies a Lipschitz condition on L1[0,1]. □
Theorem 5.2. Assume that the functions Oi satisfy the Lipschitz condition as defined in Lemma 5.1, with MO<1. Then, the MABC-MSM system (4.1) admits a unique solution.
Proof. Define the iterative difference functions
Taking norms and applying the Lipschitz condition, for S(t), we get
Since MOi<1, the above expression implies that ‖Sn+1(t)−S(t)‖→0 as n→∞.
Similarly, we obtain
where C=(1−βΦ(β)+TβΦ(β)Γ(β)).
Thus, all sequences {Sn(t)},{In(t)},{Jn(t)},{An(t)},{Bn(t)} converge, and system (4.1) admits a unique solution. □
5.5. Uniqueness of the solution
Theorem 5.3. Suppose the functions Oi, for i=1,2,…,5, satisfy the Lipschitz condition stated in Lemma 5.1 with a Lipschitz constant
If the inequality
holds, then model (4.1) admits a unique solution.
Proof. Let
be two potential solutions of system (4.1). Define the error terms:
From the integral representation of the model and the Lipschitz continuity of the functions Oi, we have
Rewriting these inequalities,
If
then each ‖Ξi‖=0, implying S(t)=˜S(t), I(t)=˜I(t). In such a case, each ‖Ξi‖=0, which implies that S(t)=˜S(t), I(t)=˜I(t), and so forth. Therefore, the solution is unique. □
6.
Stability analysis
Analyzing the stability properties of the MABC-MSM system offers critical insights into the long-term dynamics of HIV/AIDS transmission and informs the development of effective intervention strategies. In this section, we rigorously establish the Ulam-Hyers (UH) stability of the proposed model.
6.1. Ulam-Hyers stability
Definition 6.1. The MABC-MSM system (4.1) is said to be UH stable if there exist positive constants Ni>0 such that, for any κi>0 and perturbed functions ˜S(t),˜I(t),˜J(t),˜A(t),˜B(t)∈L1[0,1] satisfying
then there exist exact solutions S(t),I(t),J(t),A(t),B(t) to system (4.1) satisfying:
Remark 6.1. Solutions to Eq (5.21) for the set
exist if and only if there exist small perturbations ψ1,ψ2,ψ3,ψ4,ψ5 such that, for any t∈L, the following inequalities hold:
and
Lemma 6.1. Suppose ˜S(t),˜I(t),˜J(t),˜A(t),˜B(t) are solutions to Eq (5.21) with initial conditions as specified in (3.2). Then, for each κi>0, i=1,2,…,5, the following integral inequalities hold:
Proof. Assume that the first inequality in (6.1) is satisfied by ˜S(t). Then, according to Remark 6.1, we have
Integrating both sides and applying the structure of the MABC operator, the solution ˜S(t) can be expressed as
To estimate the deviation from the corresponding exact solution, we compute
Since |ψ1(t)|≤κ1 for all t∈L, it follows that
which is precisely the bound given in inequality (6.5).
Therefore, inequality (6.5) is satisfied. The remaining inequalities for ˜I(t),˜J(t),˜A(t),˜B(t) can be proven using an analogous procedure. □
Theorem 6.1. The modified ABC-MSM system (4.1) is Ulam-Hyers (UH) stable provided that the following condition holds. Assume that the kernels Oi satisfy the Lipschitz condition with a constant
as defined in Lemma 5.1. Then, the UH stability is ensured if
Proof. Let ˜S(t)∈L1[0,1] be an approximate solution to system (4.1), and suppose κ1>0 such that S(t)∈L1[0,1] satisfies the first inequality in (6.4). Then, the exact solution S(t) is given by
By Lemma 6.1 and the triangle inequality, we estimate the deviation
Applying the bound on the perturbation and the Lipschitz condition, we get
Rearranging inequality (6.12), we obtain:
Letting
we conclude
Applying the same reasoning for the remaining state variables yields
where
Therefore, under condition (6.11), the modified ABC-MSM system (4.1) is Ulam-Hyers stable. □
6.2. Existence of a solution via the Leray-Schauder fixed point theorem
The primary focus of this subsection is to establish the existence of a solution (root) to system (4.1) using the fixed point approach. Let ω=Ψ5 be the product space equipped with an appropriate norm, and define L=[0,T]⊂R+. The Banach space Ψ=G(L,R+) is specified by
assuming that \{S(t), I(t), J(t), A(t), B(t)\} \in \omega .
We address the compact initial value problem, formulated as
where
We define the operator \Delta: \omega \to \omega as
Theorem 6.2. Let O: \mathcal{L} \times \omega \to \mathbb{R} be a continuous function satisfying the growth condition
Then, a solution to system (4.1) exists provided that
where
Proof. Consider the operator \Delta defined in (6.20). To show the existence of a solution to system (4.1), it suffices to prove that \Delta has a fixed point.
To this end, define the closed ball V_r \subset \omega by
where
The proof continues by verifying that \Delta maps V_r into itself and is compact and continuous, thereby satisfying the conditions of the Leray-Schauder fixed point theorem. We now divide the proof into the following steps:
Step Ⅰ: Show that \Delta H(t) \in V_r .
Using the bounds on O(t, H(t)) , we estimate
Therefore,
This demonstrates that \Delta H(t) \in V_r .
Step Ⅱ: The operator \Delta is relatively compact, meaning it is equicontinuous, continuous, and uniformly bounded.
The continuity of \Delta follows directly from the continuity of the function O(t, H(t)) . Assume H(t) \in V_r . Then, using inequality (6.21), we obtain
Thus, \Delta is uniformly bounded on V_r .
Next, assume 0 < t_1 < t_2 < T . Then, we consider the difference
and estimate
Using the bound from (6.21), we further obtain
It follows that
Therefore,
By the Arzelà-Ascoli theorem, \Delta is equicontinuous and compact on V_r . Thus, by the Leray-Schauder fixed point theorem, there exists at least one solution to model (4.1). □
7.
Numerical strategy using Lagrange interpolation
In this section, we develop a numerical scheme to approximate the solution of the structural system defined by Eq (4.1). To facilitate this, we employ the compact form of the initial value problem as given in Eq (6.18). Accordingly, we consider the integral form
Discretizing Eq (7.1) at t = t_{p+1} = h(p + 1) , where h is the time step size, yields
To approximate the integrand, we apply the Lagrange interpolation formula. For simplicity, we use the two-point interpolation between t_{q-1} and t_q :
This can be written as
Substituting the interpolation approximation (7.3) into the discretized form (7.2), we obtain
Now, by evaluating the two integrals in (7.4) independently, we obtain the following numerical expression
The numerical solutions of the MABC-MSM system (4.1), derived using the scheme in (7.5) along with the structural decomposition in (6.18), are given as follows
8.
Numerical simulations and discussion
This section presents numerical simulations to analyze the dynamics of the non-linear demographic framework for MSM (men who have sex with men) with therapeutic intervention. The model incorporates fractional-order derivatives to assess how fractional parameters influence disease progression and treatment efficacy. Simulations are conducted using the compact form of the model described in Eq (4.1), and they aim to illustrate the role of fractional dynamics in the stability and behavior of the system.
Several numerical experiments are performed using varying values of the fractional order parameter \beta , demonstrating its effect on the system's compartments. The results, depicted in Figures 1–8, show the approximate solutions for different values of \beta . The simulations reveal that as \beta decreases (i.e., moving away from the classical derivative), the susceptible population S(t) increases, while the infected populations I(t) , J(t) , A(t) , and B(t) decrease. These trends suggest that lower fractional orders lead to quicker stabilization and potentially improved outcomes under treatment dynamics.
As the fractional order approaches lower values, the solution converges more quickly to a steady state, indicating that the system behaves more effectively with sub-classical dynamics. This supports the notion that the fractional-order model provides a more flexible and accurate representation of real-world epidemiological processes than the classical integer-order approach.
The simulation results emphasize the influence of parameter variations on the system's overall behavior, particularly under different fractional dynamics. Additionally, the results vividly illustrate the temporal evolution of HIV/AIDS among MSM populations, making this study highly relevant for decision-making and policy implementation in public health contexts.
To reflect the real-world HIV transmission scenario in Taiwan, the model stratifies the infected population into multiple compartments, representing different stages of infection. This stratification enables a more nuanced understanding of HIV transmission risks, and supports the development of targeted interventions and public health strategies. By integrating early-stage monitoring data and employing a novel data optimization technique, we enhance the precision of model predictions while maintaining general applicability.
The role of the fractional parameter \beta is crucial: lower values of \beta correspond to more rapid stabilization of the system, while higher values approach classical behavior. As \beta increases, the dynamics of each compartment begin to resemble those of the classical integer-order system. This finding highlights the increased complexity and representational power of the fractional-order framework.
The simulation further demonstrates how an individual's disease state evolves over time under varying conditions. The dual-parameter visualization captures the dynamic interplay between the fractional order and fractal dimensions, projecting the system's response to parameter shifts. The use of a fractal-fractional operator thus provides a richer understanding of the intricate behaviors that govern disease transmission and progression.
To ensure epidemiological relevance, the model uses actual surveillance data from Taiwan's Centers for Disease Control (TCDC). Monthly time-series data from 2005 to 2006 (during a period of rapid HIV spread among MSM) was selected to ensure sufficient sampling. Additionally, undiagnosed infection rates between 2005 and 2019 were estimated using data from another study [34], revealing a decline from 42.6% in 2005 to 12.1% in 2019. This trend highlights the effectiveness of Taiwan's healthcare system in controlling HIV.
Taiwan's approach to HIV/AIDS prevention and treatment has been highly effective, encompassing broad access to antiretroviral therapy, harm-reduction programs, and targeted educational campaigns. These efforts have contributed significantly to the reduction in new HIV cases among MSM. The healthcare system, overseen by the Taiwan Centers for Disease Control, serves as a model for global disease prevention and management. The initial conditions used in the simulation were estimated based on epidemiological data from TCDC [35], ensuring the model's alignment with real-world dynamics.
8.1. Limitations
Although the proposed model provides a good fit and offers insightful predictions, it represents an idealized scenario of HIV/AIDS transmission rather than a fully realistic depiction. The dynamics of HIV/AIDS outbreaks depend significantly on high-risk behaviors within subsets of the infected population, which are not explicitly modeled.
The small sample size used in parameter estimation presents another limitation; however, the robustness of the model can be partially verified through sensitivity analysis and adjustments to the estimation period. While this study captures key demographic characteristics of the target population, it does not account for age-specific differences in HIV/AIDS progression and transmission. Future work should aim to disaggregate the population by age group to more accurately reflect the heterogeneity in risk and disease dynamics.
9.
Conclusions
This study applied the modified Atangana-Baleanu-Caputo (MABC) fractional operator to investigate the dynamics of HIV/AIDS among men who have sex with men (MSM) using a nonlinear compartmental modeling approach. The model divides the population into five mutually exclusive compartments, capturing the complexity of HIV progression and treatment effects within the MSM population.
The major contributions of this work are summarized as follows:
● A novel mathematical framework is proposed for modeling HIV/AIDS transmission among MSM, incorporating the MABC fractional operator to capture memory and hereditary effects.
● Conditions for HIV/AIDS eradication or persistence (both with and without immune response) are derived using comparison principles and the method of compact invariant set localization.
● Theoretical analysis based on recursive sequences and fixed point theory establishes solution existence, uniqueness, global stability, persistence, and attractivity.
● Numerical simulations demonstrate enhanced accuracy over traditional methods, validating the utility of the MABC operator in capturing real-world epidemic behavior.
● Model calibration is based on limited epidemiological data from the Taiwan Centers for Disease Control (TCDC), reflecting practical constraints in data availability.
● Parameter estimates suggest the epidemic may peak before reaching a dynamic equilibrium, aligning with observed trends in Taiwan's MSM population.
● Simulation results highlight the importance of focusing on individuals in the infectious stages to curb transmission effectively.
● The influence of fractional–order values is analyzed, demonstrating that the MABC fractional derivative offers significant advantages in modeling complex biological systems such as HIV/AIDS. These findings advocate for broader application of MABC derivatives in mathematical biology and infectious disease modeling.
● Despite its promise, the model also underscores challenges, including data scarcity, difficulty in parameter estimation, computational complexity, and the lack of prior studies utilizing MABC operators in this specific epidemiological context. These issues must be addressed to ensure model robustness and reliability.
● A notable trend of increasing HIV infections among younger MSM in Taiwan was identified, emphasizing the urgent need for targeted public health strategies.
Future research should aim to refine the model by incorporating richer, more granular epidemiological data to improve predictive accuracy. Time-varying parameters, behavioral heterogeneity, and stratification by treatment adherence or co-infection status should be integrated to enhance realism. Additionally, network-based modeling approaches can be employed to capture the structure of interpersonal interactions and further inform public health decision-making in Taiwan and similar contexts.
Author contributions
Sana Ullah Saqib: Writing–original draft, software, investigation, formal analysis, data curation; Ali Hasan: Writing–original draft, software, investigation, review and editing, methodology, visualization; Yin-Tzer Shih: Conceptualization, writing–review and editing, supervision. All authors have read and approved the final version of the manuscript for publication.
Use of Generative-AI tools declaration
The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
The authors would like to thank the anonymous referees for their kind comments and valuable suggestions. The authors gratefully acknowledge the support of the National Science and Technology Council (NSTC), Taiwan, through the NSTC project under Grant No. NSTC 113-2115-M-005-004.
Conflict of interest
The authors declare that they have no conflict of interest regarding the publication of this paper.