EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6711 ISSN 1307-5543 – ejpam.com Published by New York Business Global Modeling the Spread and Control of Drug Addiction: A Mathematical and Sensitivity Analysis Approach Lubaba Yaseen1,∗, Syeda Alishwa Zanib1, Nadeem Abbas2, Wasfi Shatanawi2,3 1 Department of Mathematics, Riphah International University, Main Satyana Road, Faisalabad 44000, Pakistan 2 Department of Mathematics and Sciences, College of Humanities and Sciences, Prince Sultan University, Riyadh, 11586, Saudi Arabia 3 Department of Mathematics, Faculty of Science, The Hashemite University, P.O Box 330127, Zarqa 13133, Jordan Abstract. This paper presents a comprehensive mathematical model that captures the dynamics of drug addiction by dividing the population into distinct groups: individuals at high and low risk of addiction (Si, Sd), those prone to addiction (P ), active drug users (ID), and people receiving treatment either in hospital (TH) or through outpatient programs (Tp), as well as those who have recovered (R). The model effectively describes how individuals move between these groups, reflecting the impact of behavioral, social, and medical factors on the progression of addiction. Key analytical results include the calculation of the basic reproduction number R0 and a sensitivity analysis that identifies the most significant factors that influence the spread and control of drug abuse. These include the transmission rate (β1), the rate at which individuals shift between risk levels (∆), treatment effectiveness (η), and especially the rate at which treatment is initiated (σ). Through numerical simulations, the study shows that increasing the treatment initiation rate as a control strategy could significantly decrease the number of susceptible individuals and lower overall drug addiction rates. Visual phase portraits further clarify how different groups interact over time, and additional parameter analysis underscores the critical roles of transmission, risk changes, and particularly treatment efforts in shaping addiction trends. These results highlight the urgent requirement for quick and effective treatment-based approaches, offering insightful information to inform public health policies targeted at lowering drug addiction and promoting long-term recovery. 2020 Mathematics Subject Classifications: 26A33, 34A08, 03C65 Key Words and Phrases: Drug addiction, mathematical modeling, high-risk and low-risk pop- ulations, treatment interventions, basic reproduction number (R0), sensitivity analysis ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6711 Email addresses: lubabayaseen123@gmail.com (L. Yaseen), 19907@riphahfsd.edu.pk (S.A. Zanib), nabbas@psu.edu.sa (N. Abbas), wshatanawi@psu.edu.sa (W. Shatanawi) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 2 of 29 1. Introduction Drug addiction is a chronic disorder in which individuals compulsively use substances de- spite harmful consequences. It alters brain function and behavior, creating both physical and psychological dependence. The causes of addiction are multifactorial, encompassing genetic, environmental, and psychological influences. Some individuals may turn to drugs as a coping mechanism for stress, trauma, or mental health challenges, while others are influenced by peer pressure or unfavorable social circumstances [1, 2]. Beyond its individ- ual impact, drug addiction imposes heavy social and economic costs, including increased crime rates, strain on healthcare systems, family disruption, and reduced workforce pro- ductivity [3]. Collectively, these outcomes diminish community well-being and public safety. Mathematical models are powerful tools for understanding such complex biological and social processes. They provide structured frameworks to represent interactions, pre- dict outcomes, and reveal patterns within data [4–6]. By quantifying key variables, models highlight how changes in one factor influence system-wide behavior. Applications of math- ematical modeling span diverse fields, including infectious disease dynamics, fluid mechan- ics, and neuroscience. Importantly, these models allow researchers to simulate scenarios that may be impractical, costly, or unethical to test experimentally [7]. This predictive capacity enhances scientific understanding and guides evidence-based decision-making in medicine, neuroscience, and public health. In recent years, fractional mathematical models have gained importance because they can more accurately capture real-world dynamics, particularly memory and hereditary effects that classical integer-order models often over- look [8–11]. Fractional calculus offers greater flexibility for analyzing irregular behaviors such as anomalous diffusion [12], enabling the development of more realistic and effective solutions. The literature on drug addiction reflects its complexity. Research consistently highlights the interplay of biological, psychological, and social factors in shaping addiction outcomes. Studies emphasize the role of genetic predisposition, neurochemical alterations, and environmental stressors, as well as the challenges of relapse and chronicity. For ex- ample, De Angelis et al. (2020) [13] demonstrated the detrimental effects of smoking, alcohol, and drug addiction on female fertility, while Vandaele and Ahmed (2021) [14] explored the transition from voluntary behavior to compulsive use, focusing on underlying brain circuits. Ceceli, Bradberry, and Goldstein (2022) [15] investigated prefrontal cortex dysfunction, showing how addiction impairs decision-making and impulse control. From a modeling perspective, Zanib et al. (2024) [16] proposed a compartmental framework (SD, ED, HD, LD, RD, CD) that distinguishes between heavy and light addiction alongside rehabilitation. Their simulations using the RK4 method in Maple highlighted the impor- tance of early detection and treatment in controlling addiction. Bunaciu et al. (2024) [17] systematically reviewed recovery capital, providing tools to evaluate the strengths that support long-term recovery. Other studies, such as Mamo et al. (2024) [18], examined the dynamic relationship between crime and drug abuse, while Muli (2025) [19] extended this framework to include policing and rehabilitation, drawing parallels with infectious disease models. Alharbi et al. (2025) [20] presented a four-compartment model incorpo- rating social factors, demonstrating through stability and optimal control analyses that L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 3 of 29 targeted interventions can significantly reduce addiction rates and social costs. New ap- proaches have also integrated artificial intelligence (AI). Kim et al. (2025) [21] reviewed how AI-mediated communication enhances treatment, prevention, and control strategies, showing its potential to improve engagement and health outcomes. After reviewing the ex- isting literature, we identified several research gaps, and to address them, we developed a mathematical model that incorporates the important concept of proneness, distinguishing between individuals at high and low risk of addiction. This concept, which reflects the sus- ceptibility of individuals or populations to specific outcomes, is combined with treatment as a control strategy to mitigate drug addiction. Proneness is widely applied in assessing vulnerabilities arising from environmental and social conditions, including disease out- breaks and ecological changes [22–24]. Building on this foundation, the present paper proposes a comprehensive mathematical model of drug addiction dynamics. The model classifies populations according to risk level and treatment status, derives key threshold parameters, and examines both stability and sensitivity. The basic reproduction number is calculated to assess the transmission potential of addiction, while numerical simulations highlight the critical factors influencing addiction spread and control. The findings under- score the importance of timely treatment interventions and provide valuable insights for designing effective public health strategies to reduce substance abuse. 2. Model Formulation In this section, we develop a mathematical model to describe the dynamics of drug abuse, incorporating different treatment types and risk levels. The total human population is divided into seven compartments based on their relationship to drug addiction: • High-risk susceptible individuals, denoted by Si, • Low-risk susceptible individuals, denoted by Sd, • Prone individuals, denoted by P , who have a tendency towards addiction, • Active drug abusers, denoted by ID, • Individuals in treatment but not hospitalized, denoted by Tp, • Individuals receiving treatment in a hospital setting, denoted by TH , • Recovered individuals, denoted by R. We emphasize the distinction between low-risk and high-risk susceptible populations to capture the structure of risk, which is critical for understanding behavioral patterns. Fac- tors such as community environment, personal values, aspirations, and social support systems influence these risk levels. The parameters used in the model are summarized in Table 1. Where, ϵ1 = β1P N , ϵ2 = β1ηP N , L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 4 of 29 And based on these assumptions, the transmission diagram shown in Figure 1 is con- structed. From figure 1, the governing differential equations of the model are given by: dSd dt = ρ∆+ δ2Si − (ϵ1 + µ+ δ1)Sd, dSi dt = ρ(1−∆) + δ1Sd − (ϵ2 + µ+ δ2)Si, dP dt = ϵ2Si + ϵ1Sd − (φ+ µ)P, dID dt = φP + α3Tp + α2R− ( (1− γ)σ + γσ + µ+ ω1 + ϕ ) ID, dTp dt = (1− γ)σID + λ1TH − (α3 + Λ1 + ϱ1 + µ)Tp, dTH dt = Λ1Tp + γσID − (λ1 + ϱ2 + µ)TH , dR dt = ϕID + ϱ1Tp + ϱ2TH − (α2 + µ)R. (2.1) The model is non-negative initial conditions: Sd(0), Si(0), P (0), ID(0), Tp(0), TH(0), R(0) ≥ 0. (2.2) P 𝑺𝒊 𝑺𝒅 𝑰𝑫 R 𝑻𝑯 𝑻𝒑 𝜇𝑃 𝜇𝑆𝑑 𝜖2𝑆𝑖 𝜖1𝑆𝑑 𝛿2𝑆𝑖 𝛿1𝑆𝑑 𝜇𝑆𝑖 𝜑 𝛼3𝑇𝑃 𝜇𝑇𝐻 𝜇𝑅𝜇 + 𝜔1 𝐼𝐷 𝜙𝐼𝐷 𝛾𝜎𝐼𝐷 1 − 𝛾 𝜎 𝐼 𝐷 𝜚1𝑇𝑝 𝜚2𝑇𝐻 Λ1𝑇𝑝 𝜆1𝑇𝐻 𝜌(1 − Δ) 𝜌Δ 𝛼2𝑅 𝜇𝑇𝑃 Figure 1: Graphical representation of drug addiction dynamics. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 5 of 29 Parameter Description ∆ Proportion of new recruits in the population ρ Recruitment rate into the susceptible population γ Proportion of drug abusers hospitalized in rehabilitation centers φ Rate at which prone individuals become drug abusers α3 Transmission rate from interaction between non-hospitalized treated in- dividuals and drug abusers α2 Transmission rate from interaction between recovered individuals and drug abusers ϵ2 Transmission rate from low-risk susceptible to prone individuals ϵ1 Transmission rate from high-risk susceptible to prone individuals δ1 Transition rate from high-risk to low-risk susceptibility δ2 Transition rate from low-risk to high-risk susceptibility λ1 Transfer rate from hospital to non-hospitalized treatment Λ1 Transfer rate from non-hospitalized to hospital treatment ϱ1 Recovery rate of drug abusers under non-hospitalized care ϱ2 Recovery rate of drug abusers under hospital care σ Rate of initiation of treatment from drug abusers ω1 Death rate of drug abusers µ Natural death rate ϕ Rate at which drug abusers convert to recovered individuals Table 1: Description of model parameters. 3. Basic Properties of the Model 3.1. Positively Invariant Region Let the total population at time t be N(t) = Si(t) + Sd(t) + P (t) + ID(t) + Tp(t) + TH(t) +R(t). (3.3) Differentiating both sides of (3.3) with respect to time and substituting from the model equations yields dN dt = ∆− µN, (3.4) where ∆ is the recruitment rate and µ is the natural death rate. The solution to (3.4) is given by N(t) = ( N(0)− ∆ µ ) e−µt + ∆ µ . (3.5) As t → ∞, the exponential term vanishes and the total population approaches the steady- state value N∗ = ∆ µ . (3.6) L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 6 of 29 Therefore, the feasible region for the system, which is positively invariant and biologically meaningful, can be defined as Γ = { (Sd, Si, P, ID, Tp, TH , R) ∈ R7 + : Sd + Si + P + ID + Tp + TH +R ≤ ∆ µ } . (3.7) This ensures that all model solutions with non-negative initial conditions remain bounded and non-negative for all future time, preserving the biological relevance of the model. 3.2. Positivity and Boundedness Theorem 1 (Positivity and Boundedness). For non-negative initial conditions (Sd(0), Si(0), P (0), ID(0), Tp(0), TH(0), R(0)) ∈ R7 +, the solutions of system (2.1) satisfy: (i) Sd(t), Si(t), P (t), ID(t), Tp(t), TH(t), R(t) ≥ 0 ∀t ≥ 0. (ii) lim supt→∞ND(t) ≤ ∆ µ , where ND(t) = ∑ X(t). Proof. Consider the system (2.1). For any compartment, X(t) ∈ {Sd, Si, P, ID, Tp, TH , R}. The right-hand side satisfies: dX dt ∣∣∣∣ X=0 ≥ 0. This quasi-positivity property ensures solutions remain non-negative by Nagumo’s theo- rem. From the total population equation: dND dt = ∆− µND − ω1ID, ≤ ∆− µND. (since ω1ID ≥ 0) Solving the differential inequality using an integrating factor eµt: ND(t) ≤ ( ND(0)− ∆ µ ) e−µt + ∆ µ . Taking the limit superior as t → ∞: lim sup t→∞ ND(t) ≤ ∆ µ . Thus, all solutions eventually enter and remain in the biologically feasible region: Γ = { (Sd, Si, P, ID, Tp, TH , R) ∈ R7 + : ∑ X ≤ ∆ µ } . L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 7 of 29 4. Drug-Free Equilibrium Point The drug-free equilibrium (DFE) [25] of system (2.1) represents the state in which there are no active cases of drug abuse within the population, that is, P = ID = Tp = TH = R = 0. Theorem 2. System (2.1) admits a unique drug-free equilibrium given by, E0 = (S∗ d , S ∗ i , 0, 0, 0, 0, 0) , Where the equilibrium susceptible populations are, S∗ d = ρ(∆µ+ δ2) µ(µ+ δ1 + δ2) , S∗ i = ρ[δ1 + (1−∆)µ] µ(µ+ δ1 + δ2) . Proof. At equilibrium, the time derivatives of all compartments vanish: dSd dt = dSi dt = dP dt = dID dt = dTp dt = dTH dt = dR dt = 0. Substituting the DFE conditions P = ID = Tp = TH = R = 0 into system (2.1), the equations governing the susceptible populations reduce to, 0 = ρ∆+ δ2S ∗ i − (µ+ δ1)S ∗ d , 0 = ρ(1−∆) + δ1S ∗ d − (µ+ δ2)S ∗ i . This system can be expressed in matrix form as,( −(µ+ δ1) δ2 δ1 −(µ+ δ2) )( S∗ d S∗ i ) = ( −ρ∆ −ρ(1−∆) ) . From the first equation, solve for S∗ d as, S∗ d = ρ∆+ δ2S ∗ i µ+ δ1 . Substituting into the second equation yields, ρ(1−∆) + δ1 ( ρ∆+ δ2S ∗ i µ+ δ1 ) = (µ+ δ2)S ∗ i , Substituting the solution for S∗ i back into the equation for S∗ d yields the explicit formulas as stated. The determinant of the coefficient matrix is, det = (µ+ δ1)(µ+ δ2)− δ1δ2 = µ(µ+ δ1 + δ2) > 0, ensuring that the system admits a unique positive solution. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 8 of 29 5. Basic Reproduction Number The basic reproduction number R0 represents the expected number of secondary cases produced by one infected individual in a fully susceptible population. It serves as a critical threshold for epidemic control: • R0 < 1: Infection decrease. • R0 > 1: Infection increase. Theorem 3. Using the next generation matrix method [26], the basic reproduction number for system (2.1) is given by: R0 = β1 [((1−∆)η +∆)µ+ ηδ1 + δ2] (µ+ δ1 + δ2)(φ+ µ) . (5.8) Proof. Consider the infected compartments x = (P, ID, Tp, TH , R)T with: F = New infection matrix, V = Transition matrix, The Jacobian matrices evaluated at DFE E0 are: F =  β1[((1−∆)η+∆)µ+ηδ1+δ2] µ+δ1+δ2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  , (5.9) V =  φ+ µ 0 0 0 0 −φ σ + µ+ ω1 + ϕ −α3 0 −α2 0 −(1− γ)σ α3 + Λ1 + ϱ1 + µ −λ1 0 0 −γσ −Λ1 λ1 + ϱ2 + µ 0 0 −ϕ −ϱ1 −ϱ2 α2 + µ  . (5.10) The basic reproduction number is the spectral radius of FV−1: R0 = ρ(FV−1) = β1[((1−∆)η +∆)µ+ ηδ1 + δ2] (µ+ δ1 + δ2)(φ+ µ) . (5.11) This follows from block matrix inversion properties, as F has rank 1 and V is upper triangular. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 9 of 29 0.2 0.4 0.6 0.8 1.0 1 0.2 0.4 0.6 0.8 1.0 1.50 3.00 4.50 6.00 7.50 0 1 2 3 4 5 6 7 8 9 0 (a) 0.02 0.04 0.06 0.08 0.10 1 0.02 0.04 0.06 0.08 0.10 2 0.75 1.00 1.25 1.50 1.7 5 2.0 02. 252. 50 0.6 0.9 1.2 1.5 1.8 2.1 2.4 0 (b) 0.02 0.04 0.06 0.08 0.10 1 0.2 0.4 0.6 0.8 1.0 1.00 1.25 1.50 1.75 2.00 2.25 2.50 0.9 1.1 1.3 1.5 1.7 1.9 2.1 2.3 2.5 2.7 0 (c) 0.02 0.04 0.06 0.08 0.10 2 0.2 0.4 0.6 0.8 1.0 0.80 1.20 1.60 2.00 2.40 0.15 0.45 0.75 1.05 1.35 1.65 1.95 2.25 2.55 2.85 0 (d) 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 1.40 1.60 1.80 2.00 2.20 2.40 2.60 1.20 1.44 1.68 1.92 2.16 2.40 2.64 0 (e) 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 0.25 0.50 0.75 1.00 1.25 1.50 0.00 0.24 0.48 0.72 0.96 1.20 1.44 1.68 0 (f) Figure 2: Contour plots of R0 for different pairs of model sensitive parameters L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 10 of 29 The six subfigures presented in Figure 2 provide a comprehensive visualization of how the basic reproduction number, R0, responds to simultaneous changes in key pairs of epidemi- ological parameters. Each contour plot displays the combined effect of two parameters on R0, with color gradients indicating the magnitude of R0 across the parameter space. In general, the lighter regions correspond to higher values of R0, signifying increased poten- tial for disease transmission, while darker regions indicate lower values, suggesting more effective control. For instance, in subfigure (2a), increasing either β1 or η leads to a marked rise in R0, highlighting the importance of controlling both transmission rate and relative infectivity. The subfigures (2b) and (2c) show that higher values of δ1 and φ or η tend to reduce R0, demonstrating the benefit of enhancing transitions to lower risk or faster progression to treatment. Conversely, subfigures (2d) and (2e) reveal that increases in δ2 or η can elevate R0, emphasizing the risk associated with transitions to higher suscep- tibility and increased infectivity. Finally, subfigure (2f) illustrates that higher values of ∆ and η are associated with greater R0, underlining the combined impact of recruitment and infectivity on outbreak potential. Collectively, these plots enable clear identification of parameter regions where interventions would be most effective in reducing R0 below the epidemic threshold, thereby guiding optimal public health strategies. 5.1. Local Stability Analysis Theorem 4. The drug-free equilibrium E0 of system (2.1) is locally asymptotically stable if R0 < 1 and unstable if R0 > 1 [27]. Proof. The Jacobian matrix evaluated at E0 of the system of non-linear equations (2.1) and substitute the drug-free equilibrium points. J0 =  −µ− δ1 δ2 − (∆µ+δ2)β1 µ+δ1+δ2 0 0 0 0 δ1 −µ− δ2 β1((∆−1)µ−δ1)η µ+δ1+δ2 0 0 0 0 0 0 −µ2+(((1−∆)η+∆)β1−φ−δ1−δ2)µ+(η δ1+δ2)β1−φ (δ1+δ2) µ+δ1+δ2 0 0 0 0 0 0 φ −σ − µ− ω1 − ϕ α3 0 α2 0 0 0 (1− γ)σ −α3 − Λ1 − ϱ1 − µ λ1 0 0 0 0 γ σ Λ1 −λ1 − ϱ2 − µ 0 0 0 0 ϕ ϱ1 ϱ2 −α2 − µ  , (5.12) The local stability analysis of the drug-free equilibrium (DFE) is conducted by evaluating the Jacobian matrix at the DFE and examining its eigenvalues. The DFE is given by, (S∗ d , S ∗ i , 0, 0, 0, 0, 0) = ( ρ(∆µ+ δ2) µ(µ+ δ1 + δ2) , ρ(δ1 + (1−∆)µ) µ(µ+ δ1 + δ2) , 0, 0, 0, 0, 0 ) . The Jacobian matrix evaluated at this equilibrium is structured into blocks corresponding to non-infected (Sd, Si) and infected (P, ID, Tp, TH , R) compartments. The non-infected subsystem yields eigenvalues with negative real parts due to the trace and determinant L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 11 of 29 conditions of its submatrix: ( −(δ1 + µ) δ2 δ1 −(δ2 + µ) ) , where the trace −(δ1+ δ2+2µ) is negative, and the determinant µ(µ+ δ1+ δ2) is positive. For the infected subsystem, the critical eigenvalue is derived from the entry governing P , which simplifies to λP = (R0 − 1)(φ+ µ), where R0 = (((1−∆)η +∆)µ+ η δ1 + δ2)β1 (µ+ δ1 + δ2) (φ+ µ) . The remaining eigenvalues of the infected subsystem are determined by the submatrix: −(µ+ ω1 + ϕ+ σ) α3 0 α2 σ(1− γ) −(Λ1 + α3 + µ+ ϱ1) λ1 0 γσ Λ1 −(λ1 + µ+ ϱ2) 0 ϕ ϱ1 ϱ2 −(α2 + µ)  . This submatrix has negative diagonal entries and non-positive off-diagonal terms, ensuring eigenvalues with negative real parts if R0 < 1. Thus, the DFE is locally asymptotically stable when R0 < 1, as all eigenvalues of the Jacobian have negative real parts. Con- versely, if R0 > 1, the eigenvalue λP becomes positive, non-stable the DFE and allowing disease persistence. This stability criterion aligns with the biological interpretation that the infection dies out when each infected individual transmits to fewer than one susceptible person on average. 5.2. Global Stability Analysis Theorem 5 (Global Stability of DFE). If R0 < 1, the disease-free equilibrium E0 of system (2.1) is globally asymptotically stable (GAS) in the feasible region Γ. Proof. We employ the Castillo-Chavez stability criteria [28] by partitioning the system into: dX dt = F (X, 0) (Non-infected subsystem) (5.13) dY dt = G(X,Y ) (Infected subsystem) (5.14) where X = (Sd, Si) and Y = (P, ID, Tp, TH , R). Condition 1: The non-infected subsystem must be GAS at X∗ = (S∗ d , S ∗ i ). From (5.13): dSd dt = ρ∆+ δ2Si − (µ+ δ1)Sd dSi dt = ρ(1−∆) + δ1Sd − (µ+ δ2)Si L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 12 of 29 The Jacobian at DFE: JX = ( −(µ+ δ1) δ2 δ1 −(µ+ δ2) ) has eigenvalues with negative real parts since: • Trace: −2µ− δ1 − δ2 < 0 • Determinant: µ(µ+ δ1 + δ2) > 0 Condition 2: The infected subsystem must satisfy Yn(Sd, Si) = AUn − Yn( ˆSd, Si). Yn =  η β1Si N + β1SD N − φ− µ 0 0 0 0 φ − (1− γ)σ − γ σ − µ− ω1 − ϕ α3 0 α2 0 (1− γ)σ −α3 − Λ1 − ϱ1 − µ λ1 0 0 γ σ Λ1 −λ1 − ϱ2 − µ 0 0 ϕ ϱ1 ϱ2 −α2 − µ  (5.15) Replace Sd = S∗ d and Si = S∗ i Yn( ˆSd, Si) =  η β1S∗ i N + β1S∗ d N − φ− µ 0 0 0 0 φ − (1− γ)σ − γ σ − µ− ω1 − ϕ α3 0 α2 0 (1− γ)σ −α3 − Λ1 − ϱ1 − µ λ1 0 0 γ σ Λ1 −λ1 − ϱ2 − µ 0 0 ϕ ϱ1 ϱ2 −α2 − µ  (5.16) Ḡ(Xn, Yn) =  η β1PSi N + β1PSd N − (φ+ µ)P φP + α3ID Tp N + α2ID R N − ((1− γ)σ + γ σ + µ+ ω1 + ϕ) ID (1− γ)σ ID + λ1TH − (α3 + Λ1 + ϱ1 + µ)Tp γ σ ID + Λ1Tp − (λ1 + ϱ2 + µ)TH ϕ ID + ϱ2TH + ϱ1Tp − (α2 + µ)R  (5.17) AUn − Yn( ˆSd, Si) =  ((S∗ i −Si)η−(S∗ d−Sd ))β1P N (Rα2+Tp α3)(N−ID ) N 0 0 0  (5.18) L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 13 of 29 this shows that, Yn( ˆSd, Si), where Un represent an M matrix, it contains a non-negative off-diagonal element. Therefore, the conditions 1 and 2 are satisfied, so the drug free- equilibrium (DFE) globally asymptotically stable if R0 < 1. Here is the complete proof. 5.3. Sensitivity Analysis We conduct a sensitivity analysis to quantify how variations in model parameters influence the basic reproduction number R0 [27]. Using the normalized forward sensitivity index: Γθ = ∂R0 ∂θ × θ R0 , we assess the relative impact of each parameter θ on R0. Positive indices indicate param- eters that increase R0 when raised, while negative indices correspond to reducing effects. Transmission rate (β1): Γβ1 = 1 > 0. Recruitment proportion (∆): Γ∆ = (η − 1)µ∆ [(∆− 1)µ− δ1] η −∆µ− δ2 > 0. Relative infectivity (η): Γη = [(∆− 1)µ− δ1] η [(∆− 1)η −∆]µ− ηδ1 − δ2 > 0. Natural mortality (µ): Γµ = − µ (N1η −N2) (µ+ δ1 + δ2)(φ+ µ)D < 0. whereN1 andN2 represent numerator terms from the original expression. Higher mortality reduces R0 through population turnover. Progression rate (φ): Γφ = − φ φ+ µ < 0. Transition rates (δ1, δ2): Γδ1 = −δ1(η − 1)(∆µ+ δ2) (µ+ δ1 + δ2)D < 0, Γδ2 = −δ2[(∆− 1)µ− δ1](η − 1) (µ+ δ1 + δ2)D < 0. where D = [(∆− 1)η −∆]µ− ηδ1 − δ2. Both rates show context-dependent effects based on η values. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 14 of 29 1 1 2 Model Parameters 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Se ns iti vi ty In de x ( ) 1.00 0.05 0.17 -0.30 -0.91 -0.30 -0.30 Figure 3: Sensitivity indices of R0 for key model parameters, showing their relative impact on disease transmission This figure 3 presents the normalized sensitivity indices of the basic reproduction number R0 with respect to key model parameters. The height and color of each bar indicate the magnitude and direction of each parameter’s influence on R0. The parameter β1 has the highest positive sensitivity index (1.00), showing that increasing it will most strongly raise R0. Parameters ∆ and η also have positive but smaller effects, suggesting that their increase slightly raises R0. In contrast, µ, φ, δ1, and δ2 have negative sensitivity indices, meaning that increasing these parameters will reduce R0. The most substantial negative effect is observed for φ (−0.91), indicating its critical role in lowering disease transmission. The similar negative values for µ, δ1, and δ2 (all −0.30) highlight their moderate but important reducing effects. This analysis helps identify which parameters are most influential for targeted interventions to control the spread of drug abuse or infection. The figure visually guides policymakers on where to focus resources for maximum impact on reducing R0. 6. Numerical and Simulation Results We performed numerical simulations using realistic parameter values and initial conditions to investigate the dynamic behavior of model (2.1). Figures 4-10 illustrate the time evo- lution of each compartment in the population. These results provide insight into how the various subpopulations interact and change over time in response to drug addiction and intervention strategies. A key observation from the simulations is that as the number of in- dividuals receiving treatment increases, both in hospital (TH) and without hospitalization (Tp), the size of the susceptible populations, particularly the high-risk (Si) and low-risk (Sd) groups decreases. This trend is evident in Figures 8 and 9, which show a steady rise in the treatment compartments, while Figures 4 and 5 demonstrate a corresponding decline in susceptible individuals. The reduction in susceptibles is a direct consequence of effective L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 15 of 29 treatment, which removes individuals from the pool at risk of progressing to drug abuse or relapse. Furthermore, the number of prone individuals (P ) and untreated drug abusers (ID) initially increases but subsequently declines as more individuals enter treatment and recovery (Figures 6 and 7). The recovered compartment (R) shows a consistent increase over time (Figure 10), reflecting the cumulative effect of successful interventions. 0 20 40 60 80 100 Time 10 20 30 40 50 Po pu la tio n Sd (Low risk susceptible) Figure 4: Time evolution of the low-risk susceptible population (Sd) in model (2.1). 0 20 40 60 80 100 Time 20 40 60 80 100 Po pu la tio n Si (High risk susceptible) Figure 5: Time evolution of the high-risk susceptible population (Si) in model (2.1). L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 16 of 29 0 20 40 60 80 100 Time 2 4 6 8 10 12 14 Po pu la tio n P (prone) Figure 6: Time evolution of the prone individuals (P ) in model (2.1). 0 20 40 60 80 100 Time 3 4 5 6 7 Po pu la tio n ID (Drug abusers untreated) Figure 7: Time evolution of untreated drug abusers (ID) in model (2.1). 0 20 40 60 80 100 Time 2 3 4 5 6 7 8 9 Po pu la tio n TH (Treatment hospitalization) Figure 8: Time evolution of individuals receiving hospital treatment (TH) in model (2.1). L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 17 of 29 0 20 40 60 80 100 Time 4 6 8 10 Po pu la tio n TP (Treatment no hospitalization) Figure 9: Time evolution of individuals receiving non-hospital treatment (Tp) in model (2.1). 0 20 40 60 80 100 Time 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Po pu la tio n R (Recovering) Figure 10: Time evolution of recovered individuals (R) in model (2.1). The phase plots shown in Figure 11 offer valuable insight into the dynamic relationships between key compartments in the drug addiction model (2.1). Each subfigure illustrates the temporal interaction between two populations, revealing characteristic trajectories and dependencies. For instance, the plots of ID versus TH and ID versus Tp show that as the number of untreated drug abusers (ID) declines, the populations undergoing treatment (both hospitalized and non-hospitalized) increase, highlighting the effectiveness of treat- ment interventions in reducing active drug abuse. Similarly, the phase plots of Sd versus P and Si versus P demonstrate how decreases in susceptible populations coincide with rises and subsequent declines in prone individuals, reflecting transitions from susceptibility to risk and eventual treatment or recovery. Furthermore, phase portraits involving treat- ment and recovery compartments (TH vs. R, Tp vs. R) reveal that increases in treated individuals are followed by growth in the recovered population, confirming the success of intervention strategies. The Tp versus TH plot captures the interaction between differ- ent treatment modalities, while plots such as Sd versus Si and Si versus P illustrate the shifting balance between risk groups as the epidemic progresses. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 18 of 29 2 3 4 5 6 7 ID 2 3 4 5 6 7 8 9 T H Start End 0 20 40 60 80 100 120 140 Ti m e (a) ID vs TH 2 3 4 5 6 7 ID 4 6 8 10 T P Start End 0 20 40 60 80 100 120 140 Ti m e (b) ID vs Tp 10 20 30 40 50 Sd 0 2 4 6 8 10 12 14 P Start End 0 20 40 60 80 100 120 140 Ti m e (c) Sd vs P 10 20 30 40 50 Sd 20 40 60 80 100 S i Start End 0 20 40 60 80 100 120 140 Ti m e (d) Sd vs Si 20 40 60 80 100 Si 0 2 4 6 8 10 12 14 P Start End 0 20 40 60 80 100 120 140 Ti m e (e) Si vs P 2 3 4 5 6 7 8 9 TH 2.5 5.0 7.5 10.0 12.5 15.0 17.5 R Start End 0 20 40 60 80 100 120 140 Ti m e (f) TH vs R 4 6 8 10 TP 2.5 5.0 7.5 10.0 12.5 15.0 17.5 R Start End 0 20 40 60 80 100 120 140 Ti m e (g) Tp vs R 4 6 8 10 TP 2 3 4 5 6 7 8 9 T H Start End 0 20 40 60 80 100 120 140 Ti m e (h) Tp vs TH Figure 11: Phase plots of key compartment pairs illustrating the dynamic interactions in the drug abuse model (2.1). L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 19 of 29 6.1. Parameter Variation Effects on Model (2.1) Figures 12–16 illustrate the impact of varying key model parameters on the dynamics of all compartments in model (2.1). Each figure shows the time evolution of the model populations under different values of the respective parameter, providing direct insight into the system’s sensitivity and potential intervention points. The figure 12 show that increasing the infectivity modification parameter η would be expected to amplify the spread of drug abuse, raising the peaks of P , ID, and treatment compartments, and accelerating the decline of susceptible classes, consistent with the trends observed for β1. 0 20 40 60 80 100 Time (days) 20 40 60 80 100 Po pu la tio n (a) Si 1=0.2 1=0.3 1=0.4 1=0.5 0 20 40 60 80 100 Time (days) 10 20 30 40 50 Po pu la tio n (b) Sd 1=0.2 1=0.3 1=0.4 1=0.5 0 20 40 60 80 100 Time (days) 0 5 10 15 20 25 30 Po pu la tio n (c) P 1=0.2 1=0.3 1=0.4 1=0.5 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 14 Po pu la tio n (d) ID 1=0.2 1=0.3 1=0.4 1=0.5 0 20 40 60 80 100 Time (days) 4 6 8 10 12 14 16 Po pu la tio n (e) TP 1=0.2 1=0.3 1=0.4 1=0.5 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 14 Po pu la tio n (f) TH 1=0.2 1=0.3 1=0.4 1=0.5 0 20 40 60 80 100 Time (days) 0 5 10 15 20 25 Po pu la tio n (g) R 1=0.2 1=0.3 1=0.4 1=0.5 Figure 12: Effect of β1 on the model (2.1) Figure 13 demonstrates that increasing the transition rate from high-risk to low-risk sus- ceptibility (δ1) results in a marked decrease in the high-risk susceptible (Si) and untreated abuser (ID) populations, while the low-risk susceptible (Sd) and recovered (R) popula- tions increase. This indicates that interventions promoting transition to lower risk can effectively reduce the epidemic burden. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 20 of 29 0 20 40 60 80 100 Time (days) 0 20 40 60 80 100 Po pu la tio n (a) Si 1=0.01 1=0.03 1=0.05 1=0.07 0 20 40 60 80 100 Time (days) 10 20 30 40 50 60 Po pu la tio n (b) Sd 1=0.01 1=0.03 1=0.05 1=0.07 0 20 40 60 80 100 Time (days) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Po pu la tio n (c) P 1=0.01 1=0.03 1=0.05 1=0.07 0 20 40 60 80 100 Time (days) 3 4 5 6 7 8 9 10 Po pu la tio n (d) ID 1=0.01 1=0.03 1=0.05 1=0.07 0 20 40 60 80 100 Time (days) 4 6 8 10 12 14 Po pu la tio n (e) TP 1=0.01 1=0.03 1=0.05 1=0.07 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 Po pu la tio n (f) TH 1=0.01 1=0.03 1=0.05 1=0.07 0 20 40 60 80 100 Time (days) 0 5 10 15 20 25 Po pu la tio n (g) R 1=0.01 1=0.03 1=0.05 1=0.07 Figure 13: Effect of δ1 on the model (2.1) As depicted in Figure 14, higher values of the transition rate from low-risk to high-risk susceptibility (δ2) increase the high-risk susceptible (Si), prone (P ), and untreated abuser (ID) populations, while reducing the low-risk susceptible (Sd) and recovered (R) classes. This highlights the importance of minimizing factors that drive individuals into higher risk categories. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 21 of 29 0 20 40 60 80 100 Time (days) 20 40 60 80 100 Po pu la tio n (a) Si 2=0.01 2=0.03 2=0.05 2=0.07 0 20 40 60 80 100 Time (days) 10 20 30 40 50 60 Po pu la tio n (b) Sd 2=0.01 2=0.03 2=0.05 2=0.07 0 20 40 60 80 100 Time (days) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Po pu la tio n (c) P 2=0.01 2=0.03 2=0.05 2=0.07 0 20 40 60 80 100 Time (days) 2 3 4 5 6 7 8 Po pu la tio n (d) ID 2=0.01 2=0.03 2=0.05 2=0.07 0 20 40 60 80 100 Time (days) 4 6 8 10 12 Po pu la tio n (e) TP 2=0.01 2=0.03 2=0.05 2=0.07 0 20 40 60 80 100 Time (days) 2 4 6 8 10 Po pu la tio n (f) TH 2=0.01 2=0.03 2=0.05 2=0.07 0 20 40 60 80 100 Time (days) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Po pu la tio n (g) R 2=0.01 2=0.03 2=0.05 2=0.07 Figure 14: Effect of δ2 on the model (2.1) The figure 15 show that increasing the infectivity modification parameter η would be expected to amplify the spread of drug abuse, raising the peaks of P , ID, and treatment compartments, and accelerating the decline of susceptible classes, consistent with the trends observed for β1. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 22 of 29 0 20 40 60 80 100 Time (days) 20 40 60 80 100 Po pu la tio n (a) Si =0.1 =0.2 =0.3 =0.4 0 20 40 60 80 100 Time (days) 10 20 30 40 50 Po pu la tio n (b) Sd =0.1 =0.2 =0.3 =0.4 0 20 40 60 80 100 Time (days) 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Po pu la tio n (c) P =0.1 =0.2 =0.3 =0.4 0 20 40 60 80 100 Time (days) 3 4 5 6 7 8 9 10 Po pu la tio n (d) ID =0.1 =0.2 =0.3 =0.4 0 20 40 60 80 100 Time (days) 4 6 8 10 12 14 Po pu la tio n (e) TP =0.1 =0.2 =0.3 =0.4 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 Po pu la tio n (f) TH =0.1 =0.2 =0.3 =0.4 0 20 40 60 80 100 Time (days) 0 5 10 15 20 Po pu la tio n (g) R =0.1 =0.2 =0.3 =0.4 Figure 15: Effect of η on the model (2.1) Figure 16 show that increasing the progression rate to treatment or recovery (φ) leads to a faster decline in the untreated abuser population and a more rapid increase in the recovered class, underscoring the value of timely intervention and treatment access. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 23 of 29 0 20 40 60 80 100 Time (days) 0 20 40 60 80 100 Po pu la tio n (a) Si =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 10 20 30 40 50 Po pu la tio n (b) Sd =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 0 5 10 15 20 25 30 35 Po pu la tio n (c) P =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2 4 6 8 10 Po pu la tio n (d) ID =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 4 6 8 10 12 14 16 Po pu la tio n (e) TP =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 14 Po pu la tio n (f) TH =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Po pu la tio n (g) R =0.05 =0.1 =0.15 =0.2 Figure 16: Effect of φ on the model (2.1) Figure 17 and figure 18 demonstrate that increasing the treatment transition rate σ leads to a significant reduction in the population of untreated drug abusers (ID). As more individuals move into treatment, the number of recovered individuals (R) rises steadily. This shift also results in a lower influx of new high-risk and low-risk susceptibles, as fewer untreated abusers are present to influence them. Overall, a higher σ enhances recovery outcomes and curtails the cycle of addiction within the population. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 24 of 29 0 20 40 60 80 100 Time (days) 20 40 60 80 100 Po pu la tio n (a) Si =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 10 20 30 40 50 Po pu la tio n (b) Sd =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 14 Po pu la tio n (c) P =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2 3 4 5 6 7 8 9 Po pu la tio n (d) ID =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 4 6 8 10 12 14 16 Po pu la tio n (e) TP =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2 4 6 8 10 12 14 Po pu la tio n (f) TH =0.05 =0.1 =0.15 =0.2 0 20 40 60 80 100 Time (days) 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Po pu la tio n (g) R =0.05 =0.1 =0.15 =0.2 Figure 17: Effect of the treatment transition rate σ on the model. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 25 of 29 = 0.05 = 0.1 = 0.15 = 0.2 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Active Drug Abusers (ID) Figure 18: Spatial distribution of active drug abusers (ID) at final time for varying treatment rates (σ) across the population. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 26 of 29 Parameter Value Source β1 0.3 [18] η 0.2 assumed ρ 0.1 assumed ∆ 0.4 assumed δ1 0.05 [18] δ2 0.05 assumed ϕ 0.1 assumed µ 0.01 assumed α1 0.02 [18] α2 0.03 assumed α3 0.01 assumed γ 0.5 [18] σ 0.1 assumed ω1 0.05 [18] λ1 0.02 assumed ϱ1 0.01 [18] ϱ2 0.01 assumed Λ1 0.01 assumed Table 2: Parameter values used in the model. Initial Condition Value Si(0) 100 Sd(0) 50 P (0) 10 ID(0) 5 Tp(0) 3 TH(0) 2 R(0) 1 Table 3: Initial conditions of model compart- ments. 7. Conclusion This paper has developed and analyzed a comprehensive compartmental model to describe the dynamics of drug addiction, incorporating multiple population groups differentiated by risk status and treatment phases. The model captures key transitions governing addiction spread and recovery, with analytical results revealing the basic reproduction number R0. The sensitivity analysis reveals that the parameters η (treatment effectiveness) and β1 (transmission rate) are the most sensitive, with increases in these values leading to a significant rise in the spread of addiction within the population. Numerical simulations and phase portraits further illustrate the dynamic interplay among compartments and demonstrate the substantial impact of timely treatment interventions in reducing addiction prevalence. These findings underscore the importance of enhancing treatment accessibility and effectiveness as central strategies for public health. Future work should extend the model to include relapse, spatial heterogeneity, and socio-economic factors to better inform targeted prevention and intervention policies aimed at mitigating the societal burden of drug addiction. L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 27 of 29 Acknowledgements Authors (Dr. Nadeem Abbas and Prof. Dr. Wasfi Shatanawi) would like to thank Prince Sultan University for their support through the TAS research lab. Data Availability No data availability Conflict of Interest The authors declare that there is no conflict of interest. References [1] Louis Sanford Goodman et al. Goodman and gilman’s the pharmacological basis of therapeutics. 1549, 1996. [2] Roy A Wise and George F Koob. The development and maintenance of drug addic- tion. Neuropsychopharmacology, 39(2):254–262, 2014. [3] Pablo Ruisoto and Israel Contador. The role of stress in drug addiction. an integrative review. Physiology & behavior, 202:62–68, 2019. [4] Nadeem Abbas, Wasfi Shatanawi, and Syeda Alishwa Zanib. A spatiotemporal frame- work for malaria control integrating exposure heterogeneity and optimal strategy deployment. European Journal of Pure and Applied Mathematics, 18(3):6423–6423, 2025. [5] Syeda Alishwa Zanib and Muzamil Abbas Shah. A piecewise nonlinear fractional- order analysis of tumor dynamics: Estrogen effects and sensitivity. Modeling Earth Systems and Environment, 10(5):6155–6172, 2024. [6] Tamour Zubair, Syeda Alishwa Zanib, and Muhammad Imran Asjad. A novel defini- tion of the caputo fractional finite difference approach for maxwell fluid. Computa- tional and Applied Mathematics, 43(4):238, 2024. [7] Lore Bellaert, Amine Zerrouk, Deborah Louise Sinclair, Thomas F Martinelli, David Best, Freya Vander Laenen, Dike van de Mheen, and Wouter Vanderplasschen. Corre- lates and stability of recovery capital among persons in long-term recovery from drug addiction. In Handbook of Addiction, Recovery and Quality of Life: Cross-cutting Perspectives from Around the Globe, pages 421–434. Springer, 2025. [8] Adeyemi Olukayode Binuyo. Mathematical modelling of the addiction of drug sub- stances among students in tertiary institutions in nigeria, 2021. [9] Indah Nurun Izzati and Cicik Alfiniyah. Fractional-order model of the drug user transmission. Barekeng, 19(1):511–524, 2025. [10] Kamal Shah, Muhammad Sarwar, Thabet Abdeljawad, et al. On mathematical model L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 28 of 29 of infectious disease by using fractals fractional analysis. Discrete and Continuous Dynamical Systems-S, 17(10):3064–3085, 2024. [11] Kamal Shah, Bahaaeldin Abdalla, Thabet Abdeljawad, and Manar A Alqudah. A fractal-fractional order model to study multiple sclerosis: a chronic disease. Fractals, 32(02):2440010, 2024. [12] EM Moumine, O Balatif, and M Rachik. Modeling and mathematical analysis of drug addiction with the study of the effect of psychological and biological treatment. Mathematical Modeling and Computing, 10(3):935–943, 2023. [13] Cristina de Angelis, Antonio Nardone, Francesco Garifalos, Claudia Pivonello, Andrea Sansone, Alessandro Conforti, Carla Di Dato, Felice Sirico, Carlo Alviggi, Andrea Isidori, et al. Smoke, alcohol and drug addiction and female fertility. Reproductive biology and endocrinology, 18(1):21, 2020. [14] Y Vandaele and SH Ahmed. Habit, choice, and addiction. Neuropsychopharmacology, 46(4):689–698, 2021. [15] Ahmet O Ceceli, Charles W Bradberry, and Rita Z Goldstein. The neurobiology of drug addiction: cross-species insights into the dysfunction and recovery of the prefrontal cortex. Neuropsychopharmacology, 47(1):276–291, 2022. [16] Syeda Alishwa Zanib, Sehrish Ramzan, Nadeem Abbas, Aqsa Nazir, and Wasfi Shatanawi. A mathematical approach of drug addiction and rehabilitation control dynamic. Modeling Earth Systems and Environment, 10(2):2995–3002, 2024. [17] Adela Bunaciu, Ana-Maria Bliuc, David Best, Emily A Hennessy, Matthew J Be- langer, and Christopher SY Benwell. Measuring recovery capital for people recovering from alcohol and drug addiction: A systematic review. Addiction Research & Theory, 32(3):225–236, 2024. [18] Dejen Ketema Mamo, Mathew Ngugi Kinyanjui, Shewafera Wondimagegnhu Teklu, and Gizachew Kefelew Hailu. Mathematical modeling and analysis of the co-dynamics of crime and drug abuse. Scientific Reports, 14(1):26461, 2024. [19] Francis Musili Muli. Mathematical analysis of drug and substance abuse model with treatment and policing. Journal of Applied Mathematics, 2025(1):6666148, 2025. [20] Abdulaziz H Alharbi, MSJ Alzahrani, Fadhel Jday, and Aref Alsehaimi. Mathemat- ical analysis for a dynamic model of illicit drug consumption. AIMS Mathematics, 10(6):14784–14803, 2025. [21] Sunny Jung Kim, Viktor Clark, Jeff T Hancock, Reza Rawassizadeh, Hongfang Liu, Emmanuel A Taylor, and Vanessa B Sheppard. Leveraging artificial intelligence- mediated communication for cancer prevention and control and drug addiction: A systematic review. Translational behavioral medicine, 15(1):ibaf007, 2025. [22] Richard Rosse, Stephen Deutsch, and Melissa Chilton. Cocaine addicts prone to cocaine-induced psychosis have lower body mass index than cocaine addicts resistant to cocaine-induced psychosis-implications for the cocaine model of psychosis prone- ness. Israel Journal of Psychiatry and Related Sciences, 42(1):45, 2005. [23] Diana Sketriene, Damien Battista, Laddawan Lalert, Natcharee Kraiwattanapirom, Han Ngoc Thai, Tanawan Leeboonngam, Lori A Knackstedt, Jess Nithiananthara- jah, Priya Sumithran, Andrew J Lawrence, et al. Compulsive-like eating of high-fat L. Yaseen et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6711 29 of 29 high-sugar food is associated with ‘addiction-like’glutamatergic dysfunction in obesity prone rats. Addiction biology, 27(5):e13206, 2022. [24] Anna L Horton, Erin J Campbell, Timothy D Aumann, Katrina R O’Brien, Andrew J Lawrence, and Robyn M Brown. Addiction-like behaviour towards high-fat high-sugar food predicts relapse propensity in both obesity prone and obesity resistant c57bl/6 j mice. Progress in Neuro-Psychopharmacology and Biological Psychiatry, 121:110654, 2023. [25] Kamal Shah, Hafsa Naz, Thabet Abdeljawad, and Bahaaeldin Abdalla. Study of fractional order dynamical system of viral infection disease under piecewise derivative. CMES-Computer Modeling in Engineering & Sciences, 136(1), 2023. [26] Pauline Van den Driessche and James Watmough. Reproduction numbers and sub- threshold endemic equilibria for compartmental models of disease transmission. Math- ematical biosciences, 180(1-2):29–48, 2002. [27] AS Alqahtani, Sehrish Ramzan, Syeda Alishwa Zanib, Aqsa Nazir, Khalid Masood, and MY Malik. Mathematical modeling and simulation for malaria disease transmis- sion using the cf fractional derivative. Alexandria Engineering Journal, 101:193–204, 2024. [28] Carlos Castillo-Chavez, Zhilan Feng, Wenzhang Huang, et al. On the computation of r (o) and its role on global stability. 2001. Introduction Model Formulation Basic Properties of the Model Positively Invariant Region Positivity and Boundedness Drug-Free Equilibrium Point Basic Reproduction Number Local Stability Analysis Global Stability Analysis Sensitivity Analysis Numerical and Simulation Results Parameter Variation Effects on Model (2.1) Conclusion Acknowledgements