EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6142 ISSN 1307-5543 – ejpam.com Published by New York Business Global Understanding Transmission Dynamics of HIV/AIDS with Sensitivity Analysis and Optimal Control Strategies Saiwan M. Fatah1, Sarbaz H. A. Khoshnaw2,∗, Arkan N. Mustafa1 1 Department of Mathematics, College of Education, University of Sulaimani, Sulaymaniah 46001, Iraq 2 Department of Mathematics, College of Basic Education, University of Raparin, Kurdistan Region of Iraq Abstract. In this article, we propose a comprehensive mathematical model that investigates the dynamics of HIV/AIDS transmission, concentrating on two distinct groups: males and females. Within this theoretical framework, seven classes of individuals are considered for each gender. We compute the basic reproduction number R0, which is crucial to determine the risk of disease transmission among these populations. Furthermore, Routh-Hurwitz criterion is used in stability analysis of disease-free and endemic equilibrium points. The stability analysis of the model pro- vides that at the disease-free equilibrium point the system is locally asymptotically stable when R0 < 1, and if R0 > 1 the system is locally asymptotically unstable at the endemic equilibrium point. Additionally, the sensitivity of both reproduction number and the state variables with the key parameters are explored in order to revealing the most influential parameters that af- fect the disease transmission. This study develops an optimal control framework incorporating prevention, screening, and treatment interventions to combat HIV/AIDS transmission in gender- stratified populations. Using fourth-order Runge-Kutta methods implemented in MATLAB, we numerically solve the control system and simulate population dynamics under various intervention scenarios. Our results demonstrate that the combined implementation of all three interventions yields superior infection reduction compared to isolated measures, with significant decreases in prevalence rates observed for both genders. These findings provide quantitative evidence for pub- lic health decision-making, offering both a validated mathematical tool for epidemic modeling and practical insights for optimizing resource allocation in HIV/AIDS control programs. The frame- work advances current methodologies by systematically evaluating intervention effectiveness while accounting for gender-specific transmission dynamics. 2020 Mathematics Subject Classifications: 92C42, 34C60, 92D25, 92D30, 37N25 Key Words and Phrases: HIV/AIDS model, Stability, Sensitivity, Optimal Control, Numerical Simulation ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6142 Email addresses: saiwan.fatah@univsul.edu.iq (Saiwan M. Fatah), sarbaz.hamza@uor.edu.krd (Sarbaz H. A. Khoshnaw),arkan.mustafa@univsul.edu.iq (Arkan N. Mustafa) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 2 of 33 1. Introduction In 1981, for the first time acquired immunodeficiency syndrome (AIDS) was initially recognized as a new infectious disease within the community of homosexual in Los Angeles and San Francisco [1, 2]. Over the past forty years, the human immunodeficiency virus (HIV) has caused the deaths of numerous individuals, including men, women, and children worldwide, for instance in 2003 each day there were 15000 new infected individuals and 8000 died with HIV in the world [1]. The World Health Organization (WHO) reports that 71.3-112.8 million people have contracted HIV, and that 35.7-51.1 million people have lost their lives to the disease since the virus’s inception. There were 39.9 million HIV-positive individuals worldwide by the end of 2023 and in the year 2023 the number of people who died of HIV related illness is 630000 globally [3]. Worldwide, it is estimated that 0.6−0.7% of individuals aged 15 to 49 are infected with HIV, however the impact of the pandemic significantly differs among countries and regions. The World Health Organization reports that the most severely affected is the African Region continues with 3.4% of adults living with HIV which is almost two-thirds of all HIV-positive people globally. There are several ways in which HIV transmitted from an infected to a victim and the most common ones are vertical transmission from mother to child, homosexual contact, blood transfusions, and shared drug injection needles [4–6]. Understanding HIV infection illness and researching the epidemiological dynamics of HIV are made possible in large part by mathematical modeling. During the years 1986-1988 a series of researches on the role of mathematical modelling on HIV transmission and the transmission dynamics of AIDS were published (see for example [7–10]). Researchers interested in HIV dynamics and transmission including both applied mathematicians and biologists. Through the last four decades a huge number of researches on the development of the mathematical models of the disease and dynamics behavior, stability and sensitivity analysis of the model, and the role of basic reproduction number R0 in controlling the spread of AIDS have been published. Some researchers have developed and examined the stability of equilibrium points within HIV/AIDS models. Contingent upon whether the reproduction number R0 exceeds or falls below unity, they conducted an analysis on the stability and examined at both endemic and disease-free equilibrium points [6, 11, 12]. A mathematical model of HIV was created by Doyle et al. [13] to simulate how the virus spreads through sexual contact in a heterosexual population and the stability at the equilibrium points are carried out. They designed their model to explore the implications of a basic constraint in any sexual mixing process heterosexual population. Kaur et al. [14] investigated the impact of screening awareness and counseling on HIV transmission in the endemic region using a mathematical model in nonlinear form. They studied the impact of suitable choice of parameters to be consistent with screening awareness and counselling the disease can be removed in the population in a reasonable period of time. Naresh et ai. In [15] and Waziri et al. in [16] proposed nonlinear mathematical models to for examine the vertical transmission of HIV dynamics. In order to minimize the risk of transmitting the infection to future generations, providing treatment to individuals with HIV before they develop AIDS can significantly decrease the probability of transmission. S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 3 of 33 The obtained results in their study showed that in 1997 in the U.S fewer than 300 infants acquired HIV through vertical transmission. During the last few years, in studying the transmission and development of HIV models within the population and to the reflection of dynamics HIV epidemic, many scholars adopted the idea of compartment modelling, which is the division of the total population and corresponding categories. These models can be used in prediction of HIV incidence in a shorter long-term, they also become an important tool in analyzing the transmission and control of HIV [17]. Sensitivity analysis and optimal control of human immunodeficiency virus have studied by many researchers to show the effect of parameters on the dynamics of the transmission of HIV in a population. Various strategies such as treatment, prevention, and screening are used to control the disease [18–23]. Mathematical model is a very important scientific tool in representing the biological problems to the form of mathematical equations. In this instance, we examine and evaluate an HIV model where people are categorized as male or female. We discuss the role of the basic reproduction number in the study of stability analysis. Moreover, the elasticity of the parameters in which increase or decrease the value of the reproduction number is analyzed. Additionally, different control strategies are used in order to reduce or stop of HIV transmission in human population. The article is arranged as follows: In Section 1, the parameters and variables are introduced and the mathematical model is formulated. In Section 2, the formulated model is analyzed mathematically to find the invariant region and provided that the solutions of the model are positive. In Section 3, the equilibrium points and the basic reproduction number R0 are computed. In Section 4, the local stability of the system of equations at the equilibrium points are studied. In Section 5, local sensitivity analysis is discussed by using three different techniques. Section 6 describes the expansion of the model with the incorporation of optimal control and analyzed. In Section 7, the numerical simulation for different combinations of strategies is illustrated. In Section 8, conclusion is provided. 2. Model Equation In this section, the proposed model equation is formulated. The total number of individuals N(t) at time t is split into seven sub-populations for females and seven sub- populations for males which are fourteen sub-populations in total: susceptible populations Sf (t) and Sm(t), exposed populations Ef (t) and Em(t) which have got the sickness, but are not yet infectious, infected populations If (t) and Im(t), undetected populations Uf (t) and Um(t), detected populations Df (t) and Dm(t), treating populations Tf (t) and Tm(t), and AIDS populations Af (t) and Am(t). Therefore, N(t) = Sf (t)+Sm(t) + Ef (t) + Em(t) + If (t) + Im(t) + Uf (t) + Um(t) +Df (t) +Dm(t) + Tf (t) + Tm(t) +Af (t) +Am(t) At a constant rate Λ, the population enters the susceptible compartment. Infection occurs only through successful contact between susceptible females Sf and infected males Im and S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 4 of 33 susceptible males Sm and infected females If , with probabilities βf or βm respectively. Abstinence avoiding all forms of vaginal, anal, and oral sex are the most efficient methods of preventing the sexual transmission of HIV. For the suggested model, we account for the use of male or female condoms, which can prevent HIV transmission. The parameters ϵ and κ represent the effectiveness of condoms and the percentage of people who use condoms respectively. This leads to a reduction in transmission rates βf and βm by a factor of 1−κ ∈ where 0 < ϵ < 1 and 0 < κ < 1. Higher values of κ and ϵ indicate a better condom a situation in which more individuals use condoms when having sex respectively. In this model, the transition rates from Sf to Ef and Sm to Em are (1 − κϵ)βf and (1 − κϵ)βm respectively. Besides, the transitions from Ef and Em to If and Im are denoted by αf and αm respectively. Additionally, we have transitions from If to both Df and Uf which are symbolized by γf and θf respectively, also from Im to both Dm and Um which are symbolized by γm and θm respectively. Furthermore, there are a transition from Df and Dm to Tf and Tm as a result of an identification of HIV through medical test and denoted by ρf and ρm respectively. Moreover, the transitions from both Uf and Tf to Af are given by λf and ωf respectively and from both Um and Tm to Am are denoted by λm and ωm respectively. The natural death rate of each compartment is symbolized by η and there is an additional death rate denoted by µ due to the disease for all compartments except Sf , Ef , Sm, and Em. The transmission of the disease from an infected female to a healthy male or from an infected male to a healthy female is only through sexual contact. In Figure 1, the male and female subpopulations are shown where coupling represented by the dashed arrows. Figure 1: Flow diagram of the model S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 5 of 33 Based on the above diagram (Figure 1) and the considered assumptions, the following system of differential equations represents the mathematical model dSf dt = Λ− (1− κϵ)βfImSf − ηSf dSm dt = Λ− (1− κϵ)βmIfSm − ηSm dEf dt = (1− κϵ)βfImSf − (αf + η)Ef dEm dt = (1− κϵ)βmIfSm − (αm + η)Em dIf dt = αfEf − (µ+ η + γf + θf ) If dIm dt = αmEm − (µ+ η + γm + θm) Im dDf dt = γfIf − (µ+ η + ρf )Df dDm dt = γmIm − (µ+ η + ρm)Dm dTf dt = ρfDf − (µ+ η + ωf )Tf dTm dt = ρmDm − (µ+ η + ωm)Tm dUf dt = θfIf − (µ+ η + λf )Uf dUm dt = θmIm − (µ+ η + λm)Um dAf dt = λfUf + ωfTf − (µ+ η)Af dAm dt = λmUm + ωmTm − (µ+ η)Am (1) with the initial conditions Sf (0), Sm(0) > 0, and Ef (0), Em(0), If (0), Im(0), Df (0), Dm(0), Tf (0), Tm(0), Uf (0), Um(0), Af (0), Am(0) ≥ 0 In Table 1, the parameter values and initial variables of the model 1 are provided. Table 1: Initial estimated values of model states and param- eters Symbols Biological Meaning Estimated Values Sm(0) Initial susceptible males 19000000 Assumed Sf (0) Initial susceptible females 20000000 Assumed Em(0) Initial exposed males 100000 Assumed Ef (0) Initial exposed females 120000 Assumed Im(0) Initial infectious males 60000 Assumed If (0) Initial infectious females 65000 Assumed Dm(0) Initial detected males 25000 Assumed Df (0) Initial detected females 28000 Assumed Um(0) Initial undetected males 52000 Assumed Uf (0) Initial undetected males 60000 Assumed Tm(0) Initial treated males 9500 Assumed Tf (0) Initial treated females 11000 Assumed Am(0) Initial symptomatic AIDS males 2500 Assumed Af (0) Initial symptomatic AIDS females 2600 Assumed Λ Net flow rate into the susceptible 100,000 [22] S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 6 of 33 Symbols Biological Meaning Estimated Values η Natural mortality rate 0.0222 [24] µ Disease-induced death rate 0.333 [25] ϵ Efficiency rate of condom [0,1] [23] κ Rate of people who use condoms [0,1] [23] βm Male effective contact rate 0.025 [22] βf Female effective contact rate 0.015 [22] αm Rate at which expose males are infected 0.75 Assumed αf Rate at which exposed females are infected 0.8 Assumed γm Rate at which infected males are detected 0.512 [22] γf Rate at which infected female are detected 0.675 [22] ρm Treatment rate of detected male [0,1] Assumed ρf Treatment rate of detected female [0,1] Assumed ωm Progression rate of males under treatment to AIDS 0.004 [22] ωf Progression rate of females under treatment to AIDS 0.004 [22] θm Transmission rate of infected males to undetected 0.12 Assumed θf Transmission rate of infected females to undetected 0.12 Assumed λm Progression rate of undetected males to AIDS 0.34 Assumed λf Progression rate of undetected females to AIDS 0.34 Assumed 3. Model Analysis This section looks at the model solution’s positivity and locates the invariant area. First, we discuss the invariant region in the following theorem. Theorem 1. With a non-negative initial condition in R14 + , all solutions of the HIV/AIDS model in system (1) stays positive for all time t > 0. Proof. From the first equation of model (1) we have that dSf dt = Λ− (1− κϵ)βfImSf − ηSf > − [(1− κϵ)βfIm − η]Sf , so that, dSf dt > − [(1− κϵ)βfIm + η]Sf Solving the above inequality, we obtain Sf (t) > Sf (0)e − ∫ t 0 [(1−κϵ)βf Im(x)+η]dx. It is given that Sf (0) > 0 and clearly e− ∫ t 0 [(1−κϵ)βf Im(x)+η]dx is nonnegative. Therefore, Sf (t) is nonnegative for all time t > 0. By the same manner, we can show that the solutions for the rest of variables of system (1) are nonnegative. Hence, we obtained nonnegative solutions for the model. S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 7 of 33 Now, in the following theorem we find the positive invariant area for the model (1). Theorem 2. The initial conditions are bounded within a bounded region Σ ⊂ R5 +, ensuring that all solutions of the model remain bounded within this specified region. So that, Σ = {Sf (t), Sm(t), Ef (t), Em(t), If (t), Im(t), Uf (t), Um(t), Df (t), Dm(t), Tf (t), Tm(t), Af (t), Am(t) ∈ R14 + : N(t) ≤ 2Λ η } . Proof. N(t) denotes the total population number at time t which is the sum of all populations and given by N(t) = Sf (t)+Sm(t) + Ef (t) + Em(t) + If (t) + Im(t) + Uf (t) + Um(t) +Df (t) +Dm(t) + Tf (t) + Tm(t) +Af (t) +Am(t) Differentiating both sides of the above equation and simplifying, we acquire the fol- lowing result dN dt = 2Λ− ηN − µN + µ (Sf + Sm) dN dt ≤ 2Λ− ηN Applying the comparison theorem to the above inequality we get that, N(t) ≤ N(0)e−ηt + 2Λ η ( 1− e−ηt ) N(t) ≤ 2Λ η + ( N(0)− 2Λ η ) e−ηt It follows that N(t) → 2Λ η as t → ∞. Hence, N(t) ≤ 2Λ η . We conclude that the invariant region Σ represents the feasible solution region of the model. From theorem 1 and theorem 2 , it is sufficient to take into account the dynamics of the system 1 in the region Σ where the positiveness, existence, and uniqueness of the solution hold. 4. Equilibrium Points and Basic Reproduction Number R0 In this section, we compute two key points for system 1: the disease-free equilibrium (DFE) in the absence of the disease, and the endemic equilibrium (EE) during an outbreak. In addition, we utilize the next-generation matrix in calculating the basic reproduction number, R0. In the model 1, we set S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 8 of 33 Γf = αf + η, Γm = αm + η, ∆f = µ+ η + γf + θf , ∆m = µ+ η + γm + θm Ωf = µ+ η + ρf , Ωm = µ+ η + ρm, Ψf = µ+ η + ωf , Ψm = µ+ η + ωm Φf = µ+ η + λf , Φm = µ+ η + λm, ζ = µ+ η The system 1 becomes, dSf dt = Λ− (1− κϵ)βfImSf − ηSf dSm dt = Λ− (1− κϵ)βmIfSm − ηSm dEf dt = (1− κϵ)βfImSf − ΓfEf dEm dt = (1− κϵ)βmIfSm − ΓmEm dIf dt = αfEf −∆fIf dIm dt = αmEm −∆mIm dDf dt = γfIf − ΩfDf dDm dt = γmIm − ΩmDm dTf dt = ρfDf −ΨfTf dTm dt = ρmDm −ΨmTm dUf dt = θfIf − ΦfUf dUm dt = θmIm − ΦmUm dAf dt = λfUf + ωfTf − ζAf dAm dt = λmUm + ωmTm − ζAm (2) 4.1. Disease-Free Equilibrium Point (E0) The equilibrium point without the disease for the system 2 represents a stable state in which the population is uninfected by HIV/AIDS. We set all the variables to be zero except Sf and Sm. Hence, Ef = Em = If = Im = Uf = Um = Df = Dm = Tf = Tm = Af = Am = 0 Then we equal the right-hand side of all equations in system (2) to zero, we obtain that Sf = Sm = Λ η , and Ef = Em = If = Im = Uf = Um = Df = Dm = Tf = Tm = Af = Am = 0. Thus, the disease-free equilibrium is E0 = (Sf , Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um, Am) = ( Λ η , 0, 0, 0, 0, 0, 0, Λ η , 0, 0, 0, 0, 0, 0 ) . 4.2. Basic Reproduction Number (R0) We compute R0 by utilizing the next generation matrix. The matrix serves the purpose of calculating the reproductive number, denoted as R0. This technique defines R0 as the greatest absolute eigenvalue within the next-generation matrix. In creating this matrix, we categorize all terms of the disease classes of the model into two groups, v represents the matrix for transition cases, while f represents the matrix for newly infected cases. S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 9 of 33 The disease classes are dEf dt = (1− κϵ)βfImSf − ΓfEf dEm dt = (1− κϵ)βmIfSm − ΓmEm dIf dt = αfEf −∆fIf dIm dt = αmEm −∆mIm dDf dt = γfIf − ΩfDf dDm dt = γmIm − ΩmDm dTf dt = ρfDf −ΨfTf dTm dt = ρmDm −ΨmTm dUf dt = θfIf − ΦfUf dUm dt = θmIm − ΦmUm dAf dt = λfUf + ωfTf − ζAf dAm dt = λmUm + ωmTm − ζAm (3) In the system 3, the newly infected cases appear in the equations dEf dt and dEm dt which are f1 = (1 − κϵ)βfImSf and f7 = (1 − κϵ)βmIfSm respectively and the last twelve equations are zero f2 = f3 = f4 = f5 = f6 = f8 = f9 = f10 = f11 = f12 = 0. The rest of terms in each equation are v1, v2, . . . , v12. Hence f and v are as follows f =  f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12  =  (1− κϵ)βfImSf 0 0 0 0 0 (1− κϵ)βmIfSm 0 0 0 0 0  and v =  v1 v2 v3 v4 v5 v6 v7 v8 v9 v10 v11 v12  =  Ef −αfEf +∆fIf −γfIf +ΩfDf −ρfDf +ΨfTf −θfIf +ΦfUf −λfUf − ωfTf + ζAf ΓmEm −αmEm +∆mIm −γmIm +ΩmDm −ρmDm +ΨmTm −θmIm +ΦmUm −λmUm − ωmTm + ζAm  The Jacobian of f and v at E0 are symbolized by F and V respectively. The next- generation matrix could be computed from FV −1 which is a 12 × 12 matrix and its absolute eigenvalues are provided below λ∗ 1 = 0, λ∗ 2 = 0, λ∗ 3 = 0, λ∗ 4 = 0, λ∗ 5 = 0, λ∗ 6 = 0, λ∗ 7 = 0, λ∗ 8 = 0, λ∗ 9 = 0, λ∗ 10 = 0, |λ∗ 11| = |λ∗ 12| = Λ(1− ϵκ) √ ∆f∆mΓfΓmαfαmβfβm ∆f∆mΓfΓmη . R0 is the greatest absolute eigenvalue of the next-generation matrix which is |λ∗ 11| or |λ∗ 12|. Therefore, R0 is expressed as R0 = Λ(1− ϵκ) η √ αfαmβfβm ∆f∆mΓfΓm . or S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 10 of 33 R0 = Λ(1− ϵκ) η √ αfαmβfβm (µ+ η + γf + θf ) (µ+ η + γm + θm) (αf + η) (αm + η) The basic reproduction number, R0, in this HIV model is strongly influenced by key epidemiological factors. It increases with higher transmission rates between partners ( βf and βm ), as well as with a greater rate of new individuals entering the susceptible population (Λ). On the other hand, R0 decreases when the natural death rate (η), disease- related death rate (µ), and the rates of treatment and behavior change due to interventions (γf , γm, θf and θm) are higher. Additionally, reducing the effectiveness of intervention parameters (ϵκ) leads to a rise in R0. To bring R0 below 1 an important threshold for controlling the spread of HIV, public health strategies should focus on three main areas: raising γf and γm, and lowering βf , and βm and enhancing prevention efforts (maximizing the impact of ϵκ ). This approach is in line with established HIV control strategies, which emphasize that combining early treatment with strong prevention measures is the most effective way to reduce the spread of the virus. 4.3. Endemic Equilibrium Point (E1) The endemic equilibrium refers to a state of stability when the illness continues to exit among the population. To determine this equilibrium, we solve system 2 for variables in terms of the parameters by equating the expressions on the right side of the system. Therefore, the obtained endemic equilibrium point is E1 = ( S∗ f , E ∗ f , I ∗ f , D ∗ f , T ∗ f , U ∗ f , A ∗ f , S ∗ m, E∗ m, I∗m, D∗ m, T ∗ m, U∗ m, A∗ m ) where , S∗ f = ∆mΓm (Λαfβm(ϵκ− 1)−∆fΓfη) αfβm(ϵκ− 1) (∆mΓmη + Λαmβf (1− ϵκ)) , I∗f = Λ2αfαmβfβm(ϵκ− 1)2 −∆f∆mΓfΓmη2 ∆fΓfβm(ϵκ− 1) (Λαmβf (ϵκ− 1)−∆mΓmη) , E∗ f = ∆f αf I∗f , D∗ f = γf Ωf I∗f , T ∗ f = γfρf ΩfΨf I∗f , U∗ f = θf Φf I∗f , A∗ f = ΩfΨfλfθf +Φfγfωfρf ΩfΦfΨfζ , S∗ m = ∆fΓf (Λαmβf (ϵκ− 1)−∆mΓmη) αmβf (ϵκ− 1) (∆fΓfη + Λαfβm(1− ϵκ)) , I∗m = −∆f∆mΓfΓmη2 + Λ2αfαmβfβm(ϵκ− 1)2 ∆mΓmβf (ϵκ− 1) (Λαfβm(ϵκ− 1)−∆fΓfη) ′ , E∗ m = ∆m αm I∗m, D∗ m = γm Ωm I∗m, T ∗ m = γmρm ΩmΨm I∗m, U∗ m = θm Φm I∗m, A∗ m = ΩmΨmλmθm +Φmγmωmρm ΩmΦmΨmζ I∗m. S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 11 of 33 Since ϵ, κ ∈ (0, 1), hence (ϵκ − 1) is always negative and clearly S∗ f , S ∗ m > 0. Further- more, I∗f , I ∗ m > 0 if and only if R0 > 1. Therefore, all components of E1 are positive if and only if R0 > 1. 5. Stability Analysis of The Model at Equilibrium Points When there is no disease, the model possesses a unique disease-free steady-state de- noted as E0. Conversely, if there is the disease, the model features an endemic equilibrium which is unique and denoted by E1. In the following theorems, we study local stability at both disease-free and endemic equilibrium points. Theorem 3. The disease-free equilibrium point is locally asymptotically stable if and only if R0 < 0. Proof. The Jacobian matrix of the system at the disease-free equilibrium point E0 is a 14× 14 which is J(E0) = −η 0 0 0 0 0 0 0 0 (ϵκ−1)βfΛ η 0 0 0 0 0 −Γf 0 0 0 0 0 0 0 (−ϵκ+1)βfΛ η 0 0 0 0 0 αf −∆f 0 0 0 0 0 0 0 0 0 0 0 0 0 γf −Ωf 0 0 0 0 0 0 0 0 0 0 0 0 0 ρf −Ψf 0 0 0 0 0 0 0 0 0 0 0 θf 0 0 −Φf 0 0 0 0 0 0 0 0 0 0 0 0 ωf λf −ζ 0 0 0 0 0 0 0 0 0 (ϵκ−1)βmΛ η 0 0 0 0 −η 0 0 0 0 0 0 0 0 (−ϵκ+1)βmΛ η 0 0 0 0 0 −Γm 0 0 0 0 0 0 0 0 0 0 0 0 0 αm −∆m 0 0 0 0 0 0 0 0 0 0 0 0 0 γm −Ωm 0 0 0 0 0 0 0 0 0 0 0 0 0 ρm −Ψm 0 0 0 0 0 0 0 0 0 0 0 θm 0 0 −Φm 0 0 0 0 0 0 0 0 0 0 0 0 ωm λm −ζ  The eigenvalues of the Jacobian matrix are provided below γ1,2 = −η, γ3,4 = −ζ, γ5 = −Ωf , γ6 = −Ωm, γ7 = −Ψf , γ8 = −Ψm, γ9 = −Φf , γ10 = −Φm S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 12 of 33 and the last four eigenvalues γ11, γ12, γ13, γ14 are the roots of the following degree four polynomial f(γ) = γ4 + D̃1γ 3 + D̃2γ 2 + D̃3γ + D̃4, where D̃1 = Γf + Γm +∆f +∆m D̃2 = ΓfΓm +∆f∆m +∆mΓm +∆mΓf +∆fΓm +∆fΓf D̃3 = ∆mΓfΓm +∆fΓfΓm +∆f∆mΓm +∆f∆mΓf D̃4 = ∆f∆mΓfΓm − αfαmβfβmΛ2 η2 (1− ϵκ)2 Clearly D̃1 > 0, D̃2 > 0, and D̃3 > 0, and also D̃4 > 0 if and only if R0 < 1 ( D̃1D̃2 − D̃3 ) D̃3 − D̃1D̃4 = ((Γf + Γm +∆f +∆m)(ΓfΓm +∆f∆m +∆mΓm +∆mΓf ∆fΓm +∆fΓf )− (∆mΓfΓm +∆fΓfΓm +∆f∆mΓm +∆f∆mΓf ))(∆mΓfΓm +∆fΓfΓm +∆f∆mΓm +∆f∆mΓf )− (Γf + Γm +∆f +∆m)2(∆f∆mΓfΓm − αfαmβfβmΛ2 η2 (ϵκ− 1)2) By expanding the above expression, we obtain that ( D̃1D̃2 − D̃3 ) D̃3 − D̃1 2D̃4 > 0. Thus, Routh-Hurwitz criterion yields that all roots of the polynomial f(γ) have nega- tive real parts. This follows that all eigenvalues γ11, γ12, γ13, γ14 have negative real parts. Hence, all eigenvalues of the Jacobian matrix J (E0) have negative real parts. It is con- cluded that, the disease-free equilibrium point E0 is locally asymptotically stable if and only if R0 < 1. Theorem 4. The endemic equilibrium point E1 is locally asymptotically stable if and only if R0 > 1. Proof. At the endemic equilibrium point E1, the Jacobian matrix J (E1) has 14 eigen- values which are given below γ∗1 = −Ωf , γ∗2 = −Ψf , γ∗3 = −Φf , γ∗4 = −Ωm, γ∗5 = −Ψm, γ∗6 = −Φm, γ∗7,8 = ζ and the roots of the following degree six polynomial are the last six eigenvalues f (γ∗) = γ∗6 + C̃1γ ∗5 + C̃2γ ∗4 + C̃3γ ∗3 + C̃4γ ∗2 + C̃5γ ∗ + C̃6, where C̃1 =W̃1I ∗ m + W̃2I ∗ f +∆f +∆m + Γf + Γm + 2η C̃2 = W̃1W̃2I ∗ f I ∗ m +∆fW̃1I ∗ m +∆fW̃2I ∗ f +∆mW̃1I ∗ m +∆mW̃2I ∗ f + ΓfW̃1I ∗ m + ΓfW̃2I ∗ f + ΓmW̃1I ∗ m + ΓmW̃2I ∗ f + W̃1ηI ∗ m + W̃2ηI ∗ f +∆f∆m +∆fΓf +∆fΓm + 2η∆f +∆mΓf +∆mΓm + 2η∆m + ΓfΓm + 2ηΓf + 2ηΓm + η2, S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 13 of 33 C̃3 = W̃1W̃2∆fI ∗ f I ∗ m + W̃1W̃2∆mI∗f I ∗ m + W̃1W̃2ΓfI ∗ f I ∗ m + W̃1W̃2ΓmI∗f I ∗ m + W̃1∆f∆mI∗m + W̃2∆f∆mI∗f + W̃1∆fΓfI ∗ m + W̃2∆fΓfI ∗ f + W̃1∆fΓmI∗m + W̃2∆fΓmI∗f + W̃1∆fηI ∗ m + W̃2∆fηI ∗ f + W̃1∆mΓfI ∗ m + W̃2∆mΓfI ∗ f + W̃1∆mΓmI∗m + W̃2∆mΓmI∗f + W̃1∆mηI∗m + W̃2∆mηI∗f + W̃1ΓfΓmI∗m + W̃2ΓfΓmI∗f + W̃1ΓfηI ∗ m + W̃2ΓfηI ∗ f + W̃1ΓmηI∗m + W̃2ΓmηI∗f +∆f∆mΓf +∆f∆mΓm + 2∆f∆mη +∆fΓfΓm + 2∆fΓfη + 2∆fΓmη +∆fη 2 +∆mΓfΓm + 2∆mΓfη + 2∆mΓmη +∆mη2 + 2ΓfΓmη + Γfη 2 + Γmη2, C̃4 = ∆f∆mW̃1W̃2I ∗ f I ∗ m +∆fΓfW̃1W̃2I ∗ f I ∗ m +∆fΓmW̃1W̃2I ∗ f I ∗ m +∆mΓfW̃1W̃2I ∗ f I ∗ m +∆mΓmW̃1W̃2I ∗ f I ∗ m + ΓfΓmW̃1W̃2I ∗ f I ∗ m − W̃1W̃2αfαmS∗ fS ∗ m +∆f∆mΓfW̃1I ∗ m +∆f∆mΓfW̃2I ∗ f +∆f∆mΓmW̃1I ∗ m +∆f∆mΓmW̃2I ∗ f +∆f∆mW̃1ηI ∗ m +∆f∆mW̃2ηI ∗ f +∆fΓfΓmW̃1I ∗ m +∆fΓfΓmW̃2I ∗ f +∆fΓfW̃1ηI ∗ m +∆fΓfW̃2ηI ∗ f +∆fΓmW̃1ηI ∗ m +∆fΓmW̃2ηI ∗ f +∆mΓfΓmW̃1I ∗ m +∆mΓfΓmW̃2I ∗ f +∆mΓfW̃1ηI ∗ m +∆mΓfW̃2ηI ∗ f +∆mΓmW̃1ηI ∗ m +∆mΓmW̃2ηI ∗ f + ΓfΓmW̃1ηI ∗ m + ΓfΓmW̃2ηI ∗ f +∆f∆mΓfΓm + 2∆f∆mΓfη + 2∆f∆mΓmη +∆f∆mη2 + 2∆fΓfΓmη +∆fΓfη 2 +∆fΓmη2 + 2∆mΓfΓmη +∆mΓfη 2 +∆mΓmη2 + ΓfΓmη2, C̃5 = ∆f∆mΓfW̃1W̃2I ∗ f I ∗ m +∆f∆mΓmW̃1W̃2I ∗ f I ∗ m +∆fΓfΓmW̃1W̃2I ∗ f I ∗ m +∆mΓfΓmW̃1W̃2I ∗ f I ∗ m − 2W̃1W̃2αfαmηS∗ fS ∗ m +∆f∆mΓfΓmW̃1I ∗ m +∆f∆mΓfΓmW̃2I ∗ f +∆f∆mΓfW̃1ηI ∗ m +∆f∆mΓfW̃2ηI ∗ f +∆f∆mΓmW̃1ηI ∗ m +∆f∆mΓmW̃2ηI ∗ f +∆fΓfΓmW̃1ηI ∗ m +∆fΓfΓmW̃2ηI ∗ f +∆mΓfΓmW̃1ηI ∗ m +∆mΓfΓmW̃2ηI ∗ f + 2∆f∆mΓfΓmη +∆f∆mΓfη 2 +∆f∆mΓmη2 +∆fΓfΓmη2 +∆mΓfΓmη2, C̃6 = ∆f∆mΓfΓmW̃1W̃2I ∗ f I ∗ m − W̃1W̃2αfαmη2S∗ fS ∗ m +∆f∆mΓfΓmW̃1ηI ∗ m +∆f∆mΓfΓmW̃2ηI ∗ f +∆f∆mΓfΓmη2, here, W̃1 = (1− ϵκ)βf > 0, W̃2 = (1− ϵκ)βm > 0 If R0 > 1, it is easy to show that C̃1, C̃2, C̃3, C̃4, C̃5, C̃6 are all positive as follows: Clearly C̃1 > 0, C̃2 > 0, and C̃3 > 0. The only negative term in C̃4 is W̃1W̃2αfαmS∗ fS ∗ m, and ∆f∆mΓfΓm−W̃1W̃2αfαmS∗ fS ∗ m = 0. Hence C̃4 > 0 . In C̃5 there is only one negative term which are 2W̃1W̃2αfαmηS∗ fS ∗ m, and 2∆f∆mΓfΓmη−2W̃1W̃2αfαmηS∗ fS ∗ m = 0. Thus, S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 14 of 33 C̃5 > 0. For C̃6, since ∆f∆mΓfΓmη2 − W̃1W̃2αfαmη2S∗ fS ∗ m = 0, it follows that C̃6 > 0. The necessary condition for stability of Routh-Hurwitz Criterion holds as C̃l > 0, for i = 1, 2, 3, 4, 5, 6 if and only if R0 > 1. The Routh-Hurwitz Criterion array for f (γ∗) is given by λ∗6 : a0 a2 a4 a6 0 λ∗5 : a1 a3 a5 0 0 λ∗4 : b1 b2 b3 b4 0 λ∗3 : c1 c2 c3 c4 0 λ∗2 : d1 d2 d3 d4 0 λ∗1 : e1 e2 e3 e4 0 λ∗0 : f1 f2 f3 f4 0 where a0 = 1, a1 = c̃1, a2 = C̃2, a3 = C̃3, a4 = C̃4, a5 = C̃5, a6 = c̃6, a7 = 0, b1 = a1a2 − a0a3 a1 , b2 = a1a4 − a0a5 a1 , b3 = a1a6 − a0a7 a1 , b4 = 0, c1 = b1a3 − a1b2 b1 , c2 = b1a5 − a1b3 b1 , c3 = b1a7 − a1b4 b1 , c4 = 0, d1 = c1b2 − b1c2 c1 , d2 = c1b3 − b1c3 c1 , d3 = c1b4 − b1c4 c1 , d4 = 0, e1 = d1c2 − c1d2 d1 , e2 = d1c3 − c1d3 d1 , e3 = d1c4 − c1d4 d1 , e4 = 0, f1 = e1d2 − d1e2 e1 , f2 = e1d3 − d1e3 e1 , f3 = e1d4 − d1e4 e1 , f4 = 0. For R0 > 1 all entries of the first column of Routh-Hurwitz Criterion array are positive. Therefore, all roots of f (γ∗) have negative real parts. Thus, it is concluded that all eigen- values γ∗i for i = 1, 2, . . . 14 have negative real parts. Therefore, the endemic equilibrium point E1 is locally asymptotically stable if and only if R0 > 1. 6. Sensitivity Analysis of The Model and Basic Reproduction Number R0 Sensitivity analysis is used to determine how variables react to change in the parame- ters. The normalized forward sensitivity index of R0 with respect to a parameter p is as follows: εR0 p = ∂R0 ∂p × p R0 We already shown that the expression of R0 is given by R0 = Λ(1− ϵκ) η √ αfαmβfβm (µ+ η + γf + θf ) (µ+ η + γm + θm) (αf + η) (αm + η) . S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 15 of 33 R0 ’s normalized forward sensitivity indices for its parameters are as follows: εR0 Λ = ∂R0 ∂Λ × Λ R0 = 1, εR0 αf = ∂R0 ∂αf × αf R0 = η 2 (αf + η) , εR0 αm = ∂R0 ∂αm × αm R0 = η 2 (αf + η) , εR0 βf = ∂R0 ∂βf × βf R0 = 1 2 , εR0 βm = ∂R0 ∂βm × βm R0 = 1 2 , εR0 κ = ∂R0 ∂κ × κ R0 = ϵκ ϵκ− 1 , εR0 ϵ = ∂R0 ∂ϵ × ϵ R0 = ϵκ ϵκ− 1 , εR0 µ = ∂R0 ∂µ × µ R0 = − (2µ+ 2η + γf + θf + γm + θm)µ 2 (µ+ η + γf + θf ) (µ+ η + γm + θm) , εR0 γf = ∂R0 ∂γf × γf R0 = − γf 2 (µ+ η + γf + θf ) , εR0 γm = ∂R0 ∂γm × γm R0 = − γm 2 (µ+ η + γm + θm) , εR0 θf = ∂R0 ∂θf × θf R0 = − θf 2 (µ+ η + γf + θf ) , ∂R0 ∂θm × θm R0 = − θm 2 (µ+ η + γm + θm) , εR0 η = ∂R0 ∂η × η R0 = − 6η4 + 5η3h1 + 4η2µh2 + 4η2h3 + 3ηµh4 + 3ηαfαmh5 + 3ηh6 + 2αfαmh7 2 (αf + η) (αm + η) (µ+ η + γf + θf ) (µ+ η + γm + θm) where, h1 = 2µ+ αfαm + γf + γm + θf + θm h2 = µ+ 2αf + 2αm + γf + γm + θf + θm h3 = αfαm + (αf + αm) (γf + γm + θf + θm) + γfγm + γfθm + γmθf + θfθm h4 = (αf + αm) (µ+ γf + γm + θf + θm) + 2αfαm h5 = γf + γm + θf + θm h6 = (αf + αm) (γfγm + γfθm + γmθf + θfθm) h7 = µ2 + µγf + µγm + µθf + µθm + γfγm + γfθm + γmθf + θfθm In the below table, the sensitivity index for the estimated values of parameters are given From the information in Table 2, we can observe positive sensitivity indices for the parameters Λ, αf , αm, βf , and βm while the other parameters κ, ϵ, µ, η, γf , γm, θf , and θm exhibit negative sensitivity indices. Additionally, these parameters play a crucial role in determining the value of the reproduction number R0, which aids in analyzing the persistence or extinction of the disease within the population. It’s noteworthy that an S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 16 of 33 Table 2: Sensitivity Index of R0 Parameters Estimated Values Sensitivity Index Sign Sources Λ 100000 1 Positive [22] αf 0.8 0.01350036488 Positive Assumed αm 0.75 0.01437451438 Positive Assumed βf 0.015 0.5 Positive [22] βm 0.025 0.5 Positive [22] κ [0, 1] -0.2500000000 Negative [23] ϵ [0, 1] -0.25000000000 Negative [23] µ 0.333 -0.3134162666 Negative [25] η 0.0222 -1.048769296 Negative [24] γf 0.675 -0.2934272300 Negative [22] γm 0.512 -0.2593192868 Negative [22] θf 0.12 -0.05216484090 Negative Assumed θm 0.12 -0.06077795785 Negative Assumed increase in positive parameter values leads to an elevation in the value of R0, indicating a substantial rise in disease transmission. In contrast, increasing negative parameter values reduces the reproduction number to some amount. We now analyze the local sensitivities. In the model 1, there are fourteen compart- ments and nineteen parameters, we can use three different methods to calculate the local sensitivities which are fullnormalization, half-normalization, and non-normalization. First, the non-normalization’s equation is given by εaiπj = ∂ai ∂πj , where εaiπj denotes the sensitivity coefficient of each variable ai with respect to each parameter πj . The half-normalization sensitivities formula is therefore defined as follows: εaiπj = ∂ai ∂πj × 1 ai Lastly, full-normalization sensitivities equation is defined by εaiπj = ∂ai ∂πj × πj ai We have employed MATLAB for computational simulations using the estimated pa- rameters and initial variables provided in Table 1. The work’s findings offer a significant advancement in our understanding of the dynamics of the model. This aids in the identi- fication of crucial model parameters as well as the ways in which the parameters impact each model state. S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 17 of 33 Figure 2: Local sensitivity analysis using techniques non-normalization (a), half-normalization (b,c), and full- normalization (d) Figure 2 illustrates that the non-normalization technique significantly affects the sen- sitivity of almost all of the model’s variables to the model parameters η, µ, αm, γf , γm, θf , and θm, while they are less sensitive with the rest parameters of the model; see Fig- ure 2(a). In addition, the half-normalization method depicts that the model states Sf and Sm are highly sensitive to the parameters of the model βf and βm respectively, whereas the other system variables are less sensitive to the system parameters; see Fig- ures 2(b)and 2(c). Finally, by utilizing the full-normalization technique demonstrates that almost all model variables are extremely sensitive to the parameters of the model Λ, η, µ, αf , αm, γf , γm, θf , θm, λf , and λm, while they are less sensitive to the other model parameters, as shown in Figure 2(d). 7. Optimal Control In this section, we incorporate three time-dependent control strategies to our model 1 with the goal of determining the best effective method for controlling and reducing the number of people with AIDS within a specific period of time. The optimum control prob- lem introduces the use of the following time-dependent control measures. • Preventive control measure: symbolized by u1 include HIV/AIDS awareness and S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 18 of 33 education campaigns through health centers and social media to safeguard the indi- viduals who are unaware of their susceptibility to HIV. • Screening control measure: denoted by u2 to assist undiagnosed infected persons in screening themselves. • Treatment control measure, represented by u3, is utilized by patients to reduce viral load in the body and prevent virus progression. Integrating the control measures u1, u2 and u3 into model 1 yields the optimal control model that results as follows: dSf dt = Λ− (1− u1) (1− κϵ)βfImSf − ηSf dSm dt = Λ− (1− u1) (1− κϵ)βmIfSm − ηSm dEf dt = (1− u1) (1− κϵ)βfImSf − (αf + η)Ef dEm dt = (1− u1) (1− κϵ)βmIfSm − (αm + η)Em dIf dt = αfEf − (u2 + µ+ η + γf + θf ) If dlm dt = αmEm − (u2 + µ+ η + γm + θm) Im dDf dt = (u2 + γf ) If − (u3 + µ+ η + ρf )Df dDm dt = (u2 + γm) Im − (u3 + µ+ η + ρm)Dm dTf dt = u3Af + (u3 + ρf )Df − (µ+ η + (1− u3)ωf )Tf dTm dt = u3Am + (u3 + ρm)Dm − (µ+ η + (1− u3)ωm)Tm dUf dt = θfIf − (µ+ η + λf )Uf dUm dt = θmIm − (µ+ η + λm)Um dtf dt = λfUf + (1− u3)ωfTf − (u3 + µ+ η)Af dAm dt = λmUm + (1− u3)ωmTm − (u3 + µ+ η)Am (4) where Ef (0), Em(0), If (0), Im(0), Df (0), Dm(0), Tf (0), Tm(0), Uf (0), Um(0), Af (0), Am(0) ≥ 0. S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 19 of 33 To analyze the levels of optimal controls, it is necessary to establish the set U of Lebesgue measurable controls U = {(u1(t), u2(t), u3(t)) : 0 ≤ u1 ≤ 1, 0 ≤ u2 ≤ 1, 0 ≤ u3 ≤ 1, 0 ≤ t ≤ tf} . The objective functional used in this paper extends the HIV control framework [1, 18]. by incorporating gender-stratified compartments (Sf , Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um, Am) and three intervention controls. This work adopts the Pontryagin Maximum Principle discussed in [26] to derive necessary conditions for optimality. The proposed control method applies the Pontryagin Maximum Principle to a gender- stratified HIV model, enabling the derivation of optimal intervention strategies across three control measures: prevention, screening, and treatment. Unlike traditional control approaches that often assume uniform populations or apply single interventions in isola- tion, this method accounts for gender-specific disease dynamics and interaction patterns. This stratification allows for more precise targeting of interventions, improving the ef- fectiveness of resource allocation. Additionally, by formulating the control problem as an optimal control system, the method provides a rigorous mathematical framework for identifying time-dependent strategies that minimize infection levels. Our objective is to determine the optimal controls u∗1, u ∗ 2 and u∗3 along with the optimal solutions S∗ f , E ∗ f , I ∗ f , D ∗ f , T ∗ f , U ∗ f , A ∗ f , S ∗ m, E∗ m, I∗m, D∗ m, T ∗ m, U∗ m and A∗ m by setting the terminal time tf in order to minimize the objective functional J defined by: J (u1, u2, u3) = min u1,u2,u3 ∫ tf 0 ( b1Ef + b2If + b3Df + b4Tf + b5Uf + b6Af + b7Em + b8Im+ b9Dm + b10Tm + b11Um + b12Am + 1 2 ( w1u21 + w2u22 + w3u23 ))dt (5) Here, bi, i = 1, . . . , 12 and w1, w2, and w3 are positive constants. The terms 0.5w1u 2 1, 0.5w2u 2 2 and 0.5w3u 2 3 represent costs associated with the controls. Our goal is to reduce the population sizes of infected compartments by intervening with control measures u1(t), u2(t) and u3(t) as well as the related costs. We have the following for the three optimal controllers u∗1, u ∗ 2 and u∗3 : J (u∗1, u ∗ 2, u ∗ 3) = min {J (u1, u2, u3) , u1, u2, u3 ∈ U} Theorem 5. (Existence of optimal control solution) There is the optimal control (u∗1, u ∗ 2, u ∗ 3) along with their associated state variable solutions S∗ f , E ∗ f , I ∗ f , D ∗ f , T ∗ f , U ∗ f , A ∗ f , S ∗ m, E∗ m, I∗m, D∗ m, T ∗ m, U∗ m and A∗ m for the state initial value problem defined in (4) and (5) that minimizes J (u1, u2, u3) over the set U . Proof. We apply the method used in [27] theorem 4 to proof the existence of optimal control solution. By applying the Pontryagin Maximum Principle (PMP) as outlined in [26], we derive the essential conditions that are met by the optimal pair. Hence, the Hamiltonian function S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 20 of 33 (H) can be obtained by this principle which is defined as H (Sf , Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um, Am, u1, u2, u3) = b1Ef + b2If + b3Df + b4Tf + b5Uf + b6Af + b7Em + b8Im + b9Dm + b10Tm + b11Um + b12Am + 1 2 ( w1u 2 1 + w2u 2 2 + w3u 2 3 ) + λ1 dSf dt + λ2 dEf dt + λ3 dIf dt + λ4 dDf dt + λ5 dTf dt + λ6 dUf dt + λ7 dAf dt + λ8 dSm dt + λ9 dEm dt + λ10 dIm dt + λ11 dDm dt + λ12 dTm dt + λ13 dUm dt + λ14 dAm dt (6) Here, λi, where i = 1, 2, . . . , 14, denotes the adjoint variables associated with the state variables Sf ,Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um and Am, respectively, which must be ascertained utilising Pontryagin’s maximum principle to establish the existence of optimal pairs. Theorem 6. Let Sf , Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um and Am denote op- timum state solutions associated with the optimal control variables u1, u2 and u3 within the optimal control model. Correspondingly, there exist co-state variables λi, where i = 1, 2, . . . , 14 such that: dλi dt = −∂H ∂j Considering transversality or final time conditions, where λi (tf ) = 0 for i = 1, 2, . . . , 14, j = Sf , Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um, Am, and H is the Hamil- tonian function as defined in equation 6. Additionally, the optimal controls u∗1, u ∗ 2 and u∗3 are u∗1 = max { min { 1, (1−kϵ)[βf ImSf (λ2−λ1)+βmIfSm(λ9−λ8)] w1 } , 0 } , u∗2 = max { min { 1, If (λ3−λ4)+Im(λ10−λ11) w2 } , 0 } , u∗3 = max { min { 1, (λ4−λ5)Df+(λ7−λ5)(Af+ωfTf)+(λ11−λ12)Dm+(λ14−λ12)(Am+ωmTm) w3 } , 0 } . Proof. The Hamiltonian H corresponding to the optimum model is expressed as H = b1Ef + b2If + b3Df + b4Tf + b5Uf + b6Af + b7Em + b8Im + b9Dm + b10Tm + b11Um + b12Am + 1 2 ( w1u 2 1 + w2u 2 2 + w3u 2 3 ) + λ1 dSf dt + λ2 dEf dt + λ3 dIf dt + λ4 dDf dt + λ5 dTf dt + λ6 dUf dt + λ7 dAf dt + λ8 dSm dt + λ9 dEm dt + λ10 dIm dt + λ11 dDm dt + λ12 dTm dt + λ13 dUm dt + λ14 dAm dt So, it can be written as S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 21 of 33 H = b1Ef + b2If + b3Df + b4Tf + b5Uf + b6Af + b7Em + b8Im + b9Dm + b10Tm + b11Um + b12Am + 1 2 ( w1u 2 1 + w2u 2 2 + w3u 2 3 ) + λ1 (Λ− (1− u1) (1− κϵ)βfImSf − ηSf ) + λ2 ((1− u1) (1− κϵ)βfImSf − (αf + η)Ef ) + λ3 (αfEf − (u2 + µ+ η + γf + θf ) If ) + λ4 ((u2 + γf ) If − (u3 + µ+ η + ρf )Df ) + λ5 (u3Af + (u3 + ρf )Df − (µ+ η + (1− u3)ωf )Tf ) + λ6 (θfIf − (µ+ η + λf )Uf ) + λ7 (λfUf + (1− u3)ωfTf − (u3 + µ+ η)Af ) + λ8 (Λ− (1− u1) (1− κϵ)βmIfSm − ηSm) + λ9 ((1− u1) (1− κϵ)βmIfSm − (αm + η)Em) + λ10 (αmEm − (u2 + µ+ η + γm + θm) Im) + λ11 ((u2 + γm) Im − (u3 + µ+ η + ρm)Dm) + λ12 (u3Am + (u3 + ρm)Dm − (µ+ η + (1− u3)ωm)Tm) + λ13 (θmIm − (µ+ η + λm)Um) + λ14 (λmUm + (1− u3)ωmTm − (u3 + µ+ η)Am) Based on the PMP’s second condition, adjoint variables λi, where i = 1, 2, . . . , 14 exist such that dλi dt = −∂H ∂j , where i = 1, 2..14, and j = Sf , Ef , If , Df , Tf , Uf , Af , Sm, Em, Im, Dm, Tm, Um, Am Hence, the expressions for the adjoint equations can be formulated as: dλ1 dt = − ∂H ∂Sf = (1− u1) (1− κϵ)βfIm (λ1 − λ2) + ηλ1, dλ2 dt = − ∂H ∂Ef = −b1 + ηλ2 + αf (λ2 − λ3) , dλ3 dt = −∂H ∂If = −b2 + (µ+ η)λ3 + (u2 + γf ) (λ3 − λ4) + θf (λ3 − λ6) + (1− u1) (1− κϵ)βmSm (λ8 − λ9) , dλ4 dt = − ∂H ∂Df = −b3 + (µ+ η)λ4 + (u3 + ρf ) (λ4 − λ5) , dλ5 dt = − ∂H ∂Tf = −b4 + (µ+ η)λ5 + (1− u3)ωf (λ5 − λ7) , dλ6 dt = − ∂H ∂Uf = −b5 + (µ+ η)λ6 + λf (λ6 − λ7) dλ7 dt = − ∂H ∂Af = −b6 + (u3 + µ+ η)λ7 − u3λ5 dλ8 dt = − ∂H ∂Sm = (1− u1) (1− κϵ)βmIf (λ8 − λ9) + ηλ8 dλ9 dt = − ∂H ∂Em = −b7 + ηλ9 + αm (λ9 − λ10) S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 22 of 33 dλ10 dt = − ∂H ∂Im = −b8 + (µ+ η)λ10 + (u2 + γm) (λ10 − λ11) + θm (λ10 − λ13) + (1− u1) (1− κϵ)βfSf (λ1 − λ2) , dλ11 dt = − ∂H ∂Dm = −b9 + (µ+ η)λ11 + (u3 + ρm) (λ11 − λ12) , dλ12 dt = − ∂H ∂Tm = −b10 + (µ+ η)λ12 + (1− u3)ωm (λ12 − λ14) dλ13 dt = − ∂H ∂Um = −b11 + (µ+ η)λ13 + λm (λ13 − λ14) dλ14 dt = − ∂H ∂Am = −b12 + (u3 + µ+ η)λ14 − u3λ12 Further, the optimality conditions give us ∂H ∂ui ∣∣∣∣ ui=u∗ i = 0, where i = 1, 2, 3 Then solving for u∗1, u ∗ 2 and u∗3, we obtain u∗1 = (1−κϵ)[βf ImSf (λ2−λ1)+βmIfSm(λ9−λ8)] w1 , u∗2 = If (λ3−λ4)+Im(λ10−λ11) w2 , and u∗3 = (λ4−λ5)Df+(λ7−λ5)(Af+ωfTf)+(λ11−λ12)Dm+(λ14−λ12)(Am+ωmTm) w3 . Hence, taking into account the control constraints 0 ≤ ui ≤ 1 where i = 1, 2, 3, the optimal control variables are expressed as: u∗1 =  0 if (1−κϵ)[βf ImSf (λ2−λ1)+βmIfSm(λ9−λ8)] w1 < 0 g1 w1 if 0 < (1−κϵ)[βf ImSf (λ2−λ1)+βmIfSm(λ9−λ8)] w1 < 1 1 if 1 < (1−κϵ)[βf ImSf (λ2−λ1)+βmIfSm(λ9−λ8)] w1 u∗2 =  0 if If (λ3−λ4)+Im(λ10−λ11) w2 < 0 g2 w2 if 0 < If (λ3−λ4)+Im(λ10−λ11) w2 < 1 1 if 1 < If (λ3−λ4)+Im(λ10−λ11) w2 u∗3 =  0 if (λ4−λ5)Df+(λ7−λ5)(Af+ωfTf)+(λ11−λ12)Dm+(λ14−λ12)(Am+ωmTm) w3 < 0 g3 w3 if 0 < (λ4−λ5)Df+(λ7−λ5)(Af+ωfTf)+(λ11−λ12)Dm+(λ14−λ12)(Am+ωmTm) w3 < 1 1 if 1 < (λ4−λ5)Df+(λ7−λ5)(Af+ωfTf)+(λ11−λ12)Dm+(λ14−λ12)(Am+ωmTm) w3 where S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 23 of 33 g1 = (1− κϵ) [βfImSf (λ2 − λ1) + βmIfSm (λ9 − λ8)] , g2 = If (λ3 − λ4) + Im (λ10 − λ11) g3 = (λ4 − λ5)Df + (λ7 − λ5) (Af + ωfTf ) + (λ11 − λ12)Dm + (λ14 − λ12) (Am + ωmTm) Concisely, the optimal controls can be expressed in a compact form as: u∗1 = max {min {1, k1} , 0} , u∗2 = max {min {1, k2} , 0} , u∗3 = max {min {1, k3} , 0} , where. k1 = (1−κϵ)[βf ImSf (λ2−λ1)+βmIfSm(λ9−λ8)] w1 k2 = If (λ3−λ4)+Im(λ10−λ11) w2 k3 = (λ4−λ5)Df+(λ7−λ5)(Af+ωfTf)+(λ11−λ12)Dm+(λ14−λ12)(Am+ωmTm) w3 Therefore, we formulate the optimality system by combining the set of equations for the state variables and the set of equations for the adjoint variables with the initial conditions and final time conditions respectively dSf dt = Λ− (1− u1) (1− κϵ)βfImSf − ηSf dEf dt = (1− u1) (1− κϵ)βfImSf − (αf + η)Ef dIf dt = αfEf − (u2 + µ+ η + γf + θf ) If dDf dt = (u2 + γf ) If − (u3 + µ+ η + ρf )Df dTf dt = u3Af + (u3 + ρf )Df − (µ+ η + (1− u3)ωf )Tf dUf dt = θfIf − (µ+ η + λf )Uf dAf dt = λfUf + (1− u3)ωfTf − (u3 + µ+ η)Af dSm dt = Λ− (1− u1) (1− κϵ)βmIfSm − ηSm dEm dt = (1− u1) (1− κϵ)βmIfSm − (αm + η)Em dIm dt = αmEm − (u2 + µ+ η + γm + θm) Im dDm dt = (u2 + γm) Im − (u3 + µ+ η + ρm)Dm S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 24 of 33 dTm dt = u3Am + (u3 + ρm)Dm − (µ+ η + (1− u3)ωm)Tm dUm dt = θmIm − (µ+ η + λm)Um dAm dt = λmUm + (1− u3)ωmTm − (u3 + µ+ η)Am dλ1 dt = (1− u1) (1− κϵ)βfIm (λ1 − λ2) + ηλ1 dλ2 dt = −b1 + ηλ2 + αf (λ2 − λ3) dλ3 dt = −b2 + (µ+ η)λ3 + (u2 + γf ) (λ3 − λ4) + θf (λ3 − λ6) + (1− u1) (1− κϵ)βmsm (λ8 − λ9) dλ4 dt = −b3 + (µ+ η)λ4 + (u3 + ρf ) (λ4 − λ5) dλ5 dt = −b4 + (µ+ η)λ5 + ωf (λ5 − λ7) dλ6 dt = −b5 + (µ+ η)λ6 + λf (λ6 − λ7) dλ7 dt = −b6 + (µ+ η)λ7 dλ8 dt = (1− u1) (1− κϵ)βmIf (λ8 − λ9) + ηλ8 dλ9 dt = −b7 + ηλ9 + αm (λ9 − λ10) dλ10 dt = −b8 + (µ+ η)λ10 + (u2 + γm) (λ10 − λ11) + θm (λ10 − λ13) + (1− u1) (1− κϵ)βfSf (λ1 − λ2) dλ11 dt = −b9 + (µ+ η)λ11 + (u3 + ρm) (λ11 − λ12) dλ12 dt = −b10 + (µ+ η)λ12 + ωm (λ12 − λ14) dλ13 dt = −b11 + (µ+ η)λ13 + λm (λ13 − λ14) dλ14 dt = −b12 + (µ+ η)λ14 With conditions λi (tf ) = 0 where i = 1, 2, . . . , 14, Sf (0), Sm(0) > 0, and Ef (0), Em(0), If (0), Im(0), Df (0), Dm(0), Tf (0), Tm(0), Uf (0), Um(0), Af (0), Am(0) ≥ 0 S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 25 of 33 8. Numerical Simulations This section qualitatively examines the influence of optimum control measures on the transmission of HIV/AIDS within a community. This work utilizes the fourth-order Runge-Kutta technique which enables precise modelling of HIV transmission dynamics and the effects of various interventions across time. This method allows the evaluation of how various combinations of strategies preventative, screening, and treatment might collectively impact the different compartments including Infected and AIDS-stage indi- viduals. Numerical simulations are performed to demonstrate the analytical results from the aforementioned study by using the initial conditions and parameter values provided in Table 1 and setting the terminal time tf = 20. Additionally, the values of parameters in the objective functional are assumed and given in Table 3. Table 3: Values of parameters in the objective functional Cost Parameters Values Cost Parameters Values b1 15 b2 20 b3 20 b4 25 b5 10 b6 18 b7 15 b8 20 b9 20 b10 25 b11 10 b12 18 w1 10 w2 15 w3 20 We used MATLAB software to provide graphical representations of the state system with and without optimum control, as well as the adjoint systems within different scenar- ios. The optimal control profiles for each strategy and state system for both male and female are displayed together in all the graphs. Figure 3 highlights the effect of preventative efforts in the absence of screening and treatment strategies in controlling HIV/AIDS epidemic, showing significant reductions in both infections and AIDS cases over time for both genders. It is shown in Figure 3(a) that the infections peak sharply around year 2 which is about six million for both males and females. Contrastingly, the infected populations reduced dramatically with prevention and staying below 200,000 in the period. On the other hand, as demonstrated in Figure 3(b) the AIDS compartments peak around year 5 and reach about 450 thousand with slightly fewer female cases, whereas with prevention strategy is much lower about 25,000 cases. The AIDS cases with prevention approach zero by the year 20 while without prevention it remains significant. Prevention not only reduces case the number of males and females but also shortens and flattens the epidemic curves for both infections and AIDS. The graphs indicate that females are somewhat less impacted in both the infected and AIDS stages; S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 26 of 33 Figure 3: Simulation of optimal control of Infected and AIDS population for both male and female with prevention however, preventative measures reduce this gender disparity. Figure 4: Simulation of optimal control of Infected and AIDS population for both male and female with screening In Figure 4, the simulation of both infected and AIDS populations for males and females are depicted with and without screening control interventions and excluding prevention and treatment. The Figures 4(a) and 4(b) demonstrate the significant role of screening in decreasing the HIV/AIDS epidemic. While not as dramatically effective as prevention strategies, screening shows substantial benefits in reducing both infection rates and AIDS cases for males and females. Figure 5 illustrates the dynamics of infected and AIDS population’s progression over a 20-year period, comparing scenarios with and without treatment interventions in the absence of prevention and screening controls for both males and females. With treatment, there is minimal visible impact on the infection rates and the curves for treated and untreated scenarios are almost identical this is provided in Figure 5(a). It clearly can be S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 27 of 33 Figure 5: Simulation of optimal control of Infected and AIDS population for both male and female with treatment seen in Figure 5(b) that the peak of AIDS population of both genders with treatment is substantially lowered to about 170,000 cases for males and 160,000 for females and approach zero by year 20. This underscores the importance of treatment in improving quality of life and survival for those who infected with HIV. The model suggests that comprehensive HIV/AIDS management should combine treatment with other strategies like prevention and screening to address both infection rates and disease progression. Figure 6: Simulation of optimal control of Infected and AIDS population for both male and female with prevention and screening Figure 6 shows how the combination of preventive and screening strategies affects both infected and AIDS population’s dynamics over the period. There is a remarkable drop in using both prevention and screening strategies together in infected and AIDS compartments for both male and female populations this can be seen in Figures 6(a) and 6(b). The intervention scenario considerably flattens the epidemic curves, eliminating new infections and bringing AIDS incidence to a minimum by year 20. The significant decrease in both compartments signifies that this integrated strategy address both new infections by prevention and the progression to AIDS through screening and probable early treatment S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 28 of 33 commencement. Figure 7: Simulation of optimal control of Infected and AIDS population for both male and female with prevention and treatment Figure 7 depicts the infection and AIDS populations with and without prevention and treatment strategies and how the interventions affect the reduction of both compartments for males and females. Combining prevention and treatment controls significantly de- creases infection rates and AIDS cases for both sexes, maintaining them at minimal levels over the time period. Integrating preventative and treatment techniques can successfully manage the transmission of infection and progression to AIDS over the period. Figure 8: Simulation of optimal control of Infected and AIDS population for both male and female with screening and treatment Figure 8 demonstrates the substantial influence of screening and treatment strategies on the management of both infection and AIDS populations for men and women over 20 years, even without preventive efforts. Through screening and treatment, as shown in Figure 8(a) the number of infections for both genders reach a lower high of around 4 mil- lion whereas AIDS cases for both sexes peak significantly lower, at around 100,000 which S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 29 of 33 can be seen in Figure 8(b). Screening and treatment interventions substantially decrease infective and AIDS populations for both sexes, but not as markedly as when prevention was incorporated in the prior scenarios. Figure 9: Simulation of optimal control of Infected and AIDS population for both male and female with prevention, screening and treatment Figure 9 presents how the infected and AIDS populations change with and without combining prevention, screening and treatment strategies for both genders over time. In Figure 9(a), it is shown that the number of infected individuals is dramatically reduced with the interventions. This indicates that the combination of prevention, screening and treatment is highly effective in spreading the disease. On the other hand, the interven- tions have a considerable impact on minimizing the number of AIDS cases, similar to the infection graph which is provided in Figure 9(b). This signifies that the therapies are also effective in stopping the progression to AIDS. The treatments appear to exert an imme- diate and sustained impact, preventing the sharp increase in infections and progression to AIDS cases, even two decades later, a significant disparity persists between intervention and non-intervention situations In Figure 10, the impact of different combinations of prevention, screening, and treat- ment interventions on infected and AIDS populations for both genders over 20 years period are illustrated. In Figure 10(a), the green line is the female infected population which peaks at around 6 million about year 2. In addition, there is some reduction of women infected individuals in the use of a single intervention, while the combination of two interventions presents a greater reduction and all three strategies together is the highly effective in min- imizing infections which is drawn by black color. The graphs of female AIDS cases with various scenarios are provided in Figure 10(b). Women AIDS population peaks approxi- mately at 4.3 million about year 5. Moreover, in the case of including a single and double interventions exhibit different level of reduction whereas all three strategies together yield the fewest female AIDS cases. Further, in Figure 10(c) it is shown that the male infected population has analogous patterns to the infected female population with all different S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 30 of 33 Figure 10: Simulation of optimal control of Infected and AIDS population for both males and females with all scenarios combinations of the interventions, although male infections appear to peak slightly higher than female infections in the absence of interventions. Additionally, the illustration of male AIDS cases with different combinations of interventions are given in Figure 10(d) which closely have a similar pattern as female AIDS cases. However, male AIDS individ- uals peak at around 4.5 million at year 5 without interventions which is higher than the peak of female AIDS cases. The results show that the simultaneous use of all three interventions (prevention, screening, and treatment) is the most successful method for lowering infections and AIDS cases in both genders. Furthermore, the patterns for males and females are quite compa- rable, suggesting that interventions are equally successful for both genders. Although all interventions are beneficial, treatment alone appears to be less successful than preventive or screening alone, particularly in decreasing initial infection rates. To sum up, the pre- sented model strongly supports the implementation of comprehensive HIV/AIDS control measures that include prevention, screening, and treatment. It illustrates that although S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 31 of 33 single interventions are advantageous, a comprehensive strategy is significantly more successful in managing the pandemic for both genders over time. 9. Limitations and Future Directions While this study provides valuable insights into gender-specific HIV transmission dy- namics, several limitations should be considered. The deterministic modeling approach cannot account for random transmission events that occur in real populations. Although gender stratification offers important granularity, additional factors like age differences and risk-group variations would further improve the model’s accuracy. Our analysis also as- sumes perfect adherence to control measures, while actual implementation faces behavioral challenges like risk compensation. Nevertheless, despite these simplifications, the current model establishes a useful foundation for assessing gender-targeted HIV interventions and can be progressively refined to better reflect real-world complexity. This framework can be extended by incorporating age-structured transmission, stochas- tic dynamics, and behavioral feedback mechanisms to more accurately represent real-world HIV transmission. Future research should include cost-effectiveness analyses to inform optimal allocation of intervention resources and validate the model using country-specific data. Incorporating machine learning techniques may enhance predictive accuracy, while network-based extensions could improve targeting of high-risk populations. Furthermore, the model is adaptable to other infectious diseases, providing policymakers with a flexible tool for evaluating integrated prevention strategies under resource constraints 10. Conclusion This paper presents a framework for understanding HIV/AIDS transmission dynamics between male and female populations. We analyzed the model’s equilibrium points and determined the basic reproduction number, R0. The stability analysis showed that the disease-free equilibrium is stable when R0 < 1, while the endemic equilibrium is stable when R0 > 1. To identify the most influential factors, we conducted local sensitivity analysis using fullnormalization, half-normalization, and non-normalization methods. We also examined how key parameters affect R0. Optimal control strategies prevention, screening, and treatment were then introduced. The control problem was formulated and solved numerically, revealing the most effi- cient use of resources and timing for interventions. Results showed that combining all three strategies significantly reduces the number of infections and AIDS cases in both sexes. Numerical simulations supported the effectiveness of these interventions in reduc- ing HIV/AIDS prevalence over time. These findings highlight the significance of focused, evidence-driven strategies in public health policy. This study not only improves our understanding of HIV/AIDS dynamics in S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 32 of 33 different populations but also provides the path for effective intervention strategies that can have a dramatic impact on public health outcomes. References [1] Priti Kumar Roy et al. Mathematical models for therapeutic approaches to control HIV disease transmission. Springer, 2015. [2] Xiaodan Sun, Hiroshi Nishiura, and Yanni Xiao. Modeling methods for estimating hiv incidence: a mathematical review. Theoretical Biology and Medical Modelling, 17(1), 2020. [3] World Health Organization. Hiv/aids, 2024. [4] Ram Naresh, Agraj Tripathi, and Dileep Sharma. Modelling and analysis of the spread of aids epidemic with immigration of hiv infectives. Mathematical and computer modelling, 49(5-6):880–892, 2009. [5] Fred Brauer, Carlos Castillo-Chavez, Zhilan Feng, et al. Mathematical models in epidemiology, volume 32. Springer, 2019. [6] DA Munawwaroh, R Heri, SUS Khabibah, HP Anindita, et al. Analysis stability of hiv/aids epidemic model of different infection stage in closed community. In JOUR- NAL of Physics: Conference Series, volume 1524, page 012130. IOP Publishing, 2020. [7] Roy M Anderson. The epidemiology of hiv infection: variable incubation plus infec- tious periods and heterogeneity in sexual activity. JOURNAL of the Royal Statistical Society: Series A (Statistics in Society), 151(1):66–93, 1988. [8] Roy M Anderson. The role of mathematical models in the study of hiv transmission and the epidemiology of aids. JOURNAL of Acquired Immune Deficiency Syndromes, 1(3):241–256, 1988. [9] Robert M May and Roy M Anderson. Commentary transmission dynamics of hiv infection. Nature, 326(137):10–1038, 1987. [10] RM Anderson, GF Medley, RM May, and AM Johnson. A preliminary study of the transmission dynamics of the human immunodeficiency virus (hiv), the causative agent of aids. Mathematical Medicine and Biology: a JOURNAL of the IMA, 3(4):229–263, 1986. [11] Hai-Feng Huo, Rui Chen, and Xun-Yang Wang. Modelling and stability of hiv/aids epidemic model with treatment. Applied Mathematical Modelling, 40(13-14):6550– 6559, 2016. [12] H. R Joshi. Optimal control of an hiv immunology model. Applied Mathematics and Computation, 178(2):262–272, 2006. [13] Murray Doyle, David Greenhalgh, and Steve Blythe. Equilibrium analysis of a math- ematical model for the spread of aids in a two sex population with mixing constraints. JOURNAL of Biological Systems, 6(02):159–185, 1998. [14] Navjot Kaur, Mini Ghosh, and S Bhatia. Modelling the role of awareness and screen- ing of infectives in the transmission dynamics of hiv. World J Modelling and Simul, 12(2):97–111, 2016. [15] Ram Naresh, Agraj Tripathi, and Sandip Omar. Modelling the spread of aids epidemic S. M. Fatah, S. H. A. Khoshnaw, A. N. Mustafa / Eur. J. Pure Appl. Math, 18 (2) (2025), 6142 33 of 33 with vertical transmission. Applied Mathematics and Computation, 178(2):262–272, 2006. [16] Abdallah S Waziri, Estomih S Massawe, and Oluwole Daniel Makinde. Mathematical modelling of hiv/aids dynamics with treatment and vertical transmission. Appl. Math, 2(3):77–89, 2012. [17] Nawei Chen, Shenglong Chen, Xiaoyu Li, and Zhiming Li. Modelling and analysis of the hiv/aids epidemic with fast and slow asymptomatic infections in china from 2008 to 2021. Math Biosci Eng, 20(12):20770–20794, 2023. [18] Kumama Regassa Cheneke, Koya Purnachandra Rao, and Geremew Kenassa Edessa. Bifurcation and stabillity analysis of hiv transmission model with optimal control. JOURNAL of Mathematics, 2021(1):7471290, 2021. [19] DM Bortz and PW Nelson. Sensitivity analysis of a nonlinear lumped parameter model of hiv infection dynamics. Bulletin of mathematical biology, 66:1009–1026, 2004. [20] Marjono Marsudi and Ari Andari. Sensitivity analysis of effect of screening and hiv therapy on the dynamics of spread of hiv. Applied Mathematical Sciences, 8(155):7749–7763, 2014. [21] Honar J Hamad, Sarbaz HA Khoshnaw, and Muhammad Shahzad. Model analy- sis for an hiv infectious disease using elasticity and sensitivity techniques. AIMS Bioengineering, 11(3):281–300, 2024. [22] Augustine S Mbitila and Jean M Tchuenche. Hiv/aids model with early detection and treatment. International Scholarly Research Notices, 2012(1):185939, 2012. [23] Dipo Aldila, Ranandha P Dhanendra, Sarbaz HA Khoshnaw, Juni Wijayanti Puspita, Putri Zahra Kamalia, and Muhammad Shahzad. Understanding hiv/aids dynamics: insights from cd4+ t cells, antiretroviral treatment, and country-specific analysis. Frontiers in Public Health, 12:1324858, 2024. [24] WHO/UNAIDS. UNAIDS report on the global AIDS epidemic 2010. UNAIDS/10.11E JC1958E,2010, 2010. [25] Z Mukandavire and W Garira. Sex-structured hiv/aids model to analyse the effects of condom use with application to zimbabwe. JOURNAL of mathematical biology, 54(5):669–699, 2007. [26] Lev Semenovich Pontryagin. Mathematical theory of optimal processes. Routledge, 2018. [27] Tigabu Kasia Ayele, Emile Franc Doungmo Goufo, and Stella Mugisha. Mathematical modeling of hiv/aids with optimal control: a case study in ethiopia. Results in Physics, 26:104263, 2021.