A novel hybrid fractional approach to nonlinear dynamics of HIV transmission among men who have sex with men in Taiwan

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 $\beta$ and fractal dimension $\alpha$, with a kernel of the modified ABC type, is defined as

It is straightforward to verify that $^{MABC}\mathbb{D}_{0, t}^{\beta, \alpha}C = 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 $\beta$ and dimension $\alpha$, with a modified AB-type kernel, is defined as

where $^{RL}\mathbb{J}$ denotes the Riemann-Liouville fractional integral.

Lemma 2.1. ^[31] Let $f^{(n)} \in L^1(0, \infty)$, with $n - 1 < \beta < n$, $n \in \mathbb{N}$, and $\beta = \alpha - n + 1$. Then, the following identity holds:

In the special case where $0 < \beta < 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 $\mathbb{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 $\mathbb{N}(t)$. Summing the equations in (4.1) yields

Given that all compartments are non-negative and $\delta > 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, $\mathbb{R}_+^5$ denotes the positive cone in five dimensions, including boundary faces. For practical and epidemiological relevance, we exclude the case where $^{MABC}\mathbb{D}_{0, t}^{\beta, \alpha} \mathbb{N}(t) \geq 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 $\mathcal{R}_0$, 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 $\mathcal{R}_0$, 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 $\mathcal{R}_0$ with respect to each model parameter:

Positive partial derivatives indicate that increasing the corresponding parameter will increase $\mathcal{R}_0$, 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 $\beta_1$, $\beta_2$, and $g$, which have the most significant influence on $\mathcal{R}_0$. 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 $\Omega$ of the MABC-MSM system

Theorem 5.1. The region $\Omega$ 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 \geq 0$. Therefore, once the solution enters the positive orthant of $\mathbb{R}_+^5$, 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 $\Omega \subset \mathbb{R}_+^5$ 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 $O_i(t, x(t))$, for $i = 1, \dots, 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 $y_i$, for $1 \le i \le 5$ such that

Under these bounds, the nonlinear functions $O_i$ defined in (5.8) satisfy Lipschitz continuity, with a global Lipschitz constant

where

Proof. Consider the function $O_1(t, S(t))$. For $S(t), \widetilde{S}(t) \in L^{1}[0, 1]$, we have

Letting $\mathcal{M}_{O_1} = \mu y_1 + \beta_1 y_2 + \beta_2 y_4$, we rewrite as follows:

Similarly, for the remaining functions with upper bounds,

Since $\mathcal{M}_O = \max_{1 \le i \le 5} \mathcal{M}_{O_i} < \infty$, each $O_i$ satisfies a Lipschitz condition on $L^1[0, 1]$. □

Theorem 5.2. Assume that the functions $O_i$ satisfy the Lipschitz condition as defined in Lemma 5.1, with $\mathcal{M}_O < 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 $\mathcal{M}_{O_i} < 1$, the above expression implies that $\|S^{n+1}(t) - S(t)\| \to 0$ as $n \to \infty$.

Similarly, we obtain

where $C = \left(\frac{1 - \beta}{\Phi(\beta)} + \frac{T^{\beta}}{\Phi(\beta) \Gamma(\beta)} \right)$.

Thus, all sequences $\{S^n(t)\}, \{I^n(t)\}, \{J^n(t)\}, \{A^n(t)\}, \{B^n(t)\}$ converge, and system (4.1) admits a unique solution. □

5.5.   Uniqueness of the solution

Theorem 5.3. Suppose the functions $O_i$, for $i = 1, 2, \ldots, 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 $O_i$, we have

Rewriting these inequalities,

If

then each $\|\Xi_i\| = 0$, implying $S(t) = \widetilde{S}(t)$, $I(t) = \widetilde{I}(t)$. In such a case, each $\|\Xi_i\| = 0$, which implies that $S(t) = \widetilde{S}(t)$, $I(t) = \widetilde{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 $N_i > 0$ such that, for any $\kappa_i > 0$ and perturbed functions $\widetilde{S}(t), \widetilde{I}(t), \widetilde{J}(t), \widetilde{A}(t), \widetilde{B}(t) \in L^{1}[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 $\psi_1, \psi_2, \psi_3, \psi_4, \psi_5$ such that, for any $t \in \mathcal{L}$, the following inequalities hold:

and

Lemma 6.1. Suppose $\widetilde{S}(t), \widetilde{I}(t), \widetilde{J}(t), \widetilde{A}(t), \widetilde{B}(t)$ are solutions to Eq (5.21) with initial conditions as specified in (3.2). Then, for each $\kappa_i > 0$, $i = 1, 2, \ldots, 5$, the following integral inequalities hold:

Proof. Assume that the first inequality in (6.1) is satisfied by $\widetilde{S}(t)$. Then, according to Remark 6.1, we have

Integrating both sides and applying the structure of the MABC operator, the solution $\widetilde{S}(t)$ can be expressed as

To estimate the deviation from the corresponding exact solution, we compute

Since $|\psi_1(t)| \leq \kappa_1$ for all $t \in \mathcal{L}$, it follows that

which is precisely the bound given in inequality (6.5).

Therefore, inequality (6.5) is satisfied. The remaining inequalities for $\widetilde{I}(t), \widetilde{J}(t), \widetilde{A}(t), \widetilde{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 $O_i$ satisfy the Lipschitz condition with a constant

as defined in Lemma 5.1. Then, the UH stability is ensured if

Proof. Let $\widetilde{S}(t) \in L^1[0, 1]$ be an approximate solution to system (4.1), and suppose $\kappa_1 > 0$ such that $S(t) \in L^1[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 $\omega = \Psi^5$ be the product space equipped with an appropriate norm, and define $\mathcal{L} = [0, T] \subset \mathbb{R}^+$. The Banach space $\Psi = G(\mathcal{L}, \mathbb{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.