EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6447 ISSN 1307-5543 – ejpam.com Published by New York Business Global Analysis of a Smoking Dynamics Model with Age-Dependent Incidence Function Samia Bushnaq1, Sameera Bano2, Anwar Zeb2,∗ 1 Department of Basic Sciences, Princess Sumaya University for Technology, Amman 11941, Jordan 2 Department of Mathematics COMSATS University Islamabad, Abbottabad Campus, KPK, Pakistan Abstract. Smoking remains a major public health concern, influenced significantly by age structure and household environments. In this study, we develop and analyze an age-structured mathematical model that incorporates age-dependent incidence and household effects to better understand smoking dynamics. We establish the model’s well-posedness by proving the existence and uniqueness of solutions. The global stability of both smoking-free and endemic equilibria is analyzed using a suitably constructed Lyapunov function. For numerical simulations, we apply the nonstandard finite difference (NSFD) method and compare the results with those obtained via the classical fourth-order Runge–Kutta (RK4) scheme. The findings underscore the importance of age-specific interventions and offer a framework for future studies incorporating relapse, control strategies, or stochastic effects. 2020 Mathematics Subject Classifications: 92D30, 34D20, 65L12, 37N25 Key Words and Phrases: Age-Structured Smoking Model, Equilibria, Stability Analysis, Lyapunov Function, Non-Standard Finite Difference Technique 1. Introduction Mathematical modeling has become an indispensable tool in modern epidemiology and pub- lic health research. These models have been extensively employed to study a wide range of diseases, including infectious diseases such as influenza, COVID-19, HIV/AIDS, and tubercu- losis, as well as non-communicable diseases like cancer, diabetes, and cardiovascular conditions [1–4]. By incorporating biological, environmental, social, and behavioral variables, mathemat- ical models help predict disease dynamics, assess the impact of interventions, and formulate effective public health policies [5–7]. For instance, compartmental models such as SIR and SEIR frameworks have been used to analyze the spread of infectious agents, while more com- plex models with delays, spatial heterogeneity, and control strategies have been applied to cancer progression, substance abuse, and mental health disorders [8, 9]. In recent years, the scope of mathematical modeling has expanded to behavioral health problems, including tobacco use and smoking. Smoking remains a persistent global health chal- lenge, contributing significantly to mortality and morbidity rates across both developed and developing nations. It is directly linked to multiple life-threatening conditions, including respi- ratory diseases, stroke, cardiovascular complications, and various forms of cancer. Statistically, a substantial proportion of the adult American population engages in smoking, with preva- lence exceeding that of the adult female population [10–12]. Epidemiological studies reveal that ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6447 Email addresses: S.Bushnaq@psut.edu.jo (S. Bushnaq), sameera@cuiatd.edu.pk (S. Bano), anwar@cuiatd.edu.pk (A. Zeb) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 2 of 12 smokers face a fourteen-fold increase in the mortality risk from cancers of the respiratory and oral systems compared to nonsmokers [13–15]. While healthcare professionals and pharmaceutical researchers are actively striving to com- bat smoking through cessation programs and drug therapies, mathematicians also play a cru- cial role by developing analytical models that improve understanding of smoking behaviors and support evidence-based medical interventions. These models help analyze smoking initi- ation, continuation, cessation, and relapse by integrating behavioral, psychological, and social components into the mathematical framework [16–20]. Traditionally, smoking behavior has been modeled using classical integer-order differential equations [21, 22]. However, recent approaches have employed fractional-order differential equa- tions due to their capability to capture memory and hereditary effects, which are especially rel- evant for addictive behaviors like smoking [23, 24]. Stochastic models also contribute to a more realistic representation by incorporating randomness and uncertainty, reflecting population-level variability in smoking tendencies. Additionally, the inclusion of age-structured modeling has gained prominence in the study of both infectious and non-infectious diseases. Age-specific modeling allows researchers to explore differences in susceptibility, behavior, and exposure across age groups, a feature particularly important in smoking dynamics where younger populations exhibit higher vulnerability to ad- diction. Motivated by these observations, we propose and study an age-structured smoking model. The model’s development and analysis of the endemic and smoking-free equilibrium points come first. The existence and uniqueness of the model solutions are then established. To investigate long-term dynamics, we perform a global stability analysis of the equilibrium states using an appropriate Lyapunov function. Furthermore, we fix the age variable and implement a nonstan- dard finite difference (NSFD) approach to numerically approximate the model’s solutions, which are then compared graphically with those obtained using the classical fourth-order Runge–Kutta (RK4) method via MATLAB simulations. The remainder of the paper is structured as follows. Section 2 introduces the formulation of the age-structured model and discusses the existence and uniqueness of solutions. Section 3 analyzes the equilibrium points and their local stability. Section 4 presents the global stability results using the Lyapunov method. In Section 5, we employ the NSFD technique for numerical simulations. Finally, Section 6 concludes the paper with key findings. 2. Model Formulation and Mathematical Well-Posedness We introduce an age-structured smoking model that incorporates household effects, such as school vacations, that influence smoking initiation rates. The model consists of four compart- ments: susceptible (S), exposed (E), smokers (I), and recovered (R). The system of ordinary differential equations (ODEs) is given by: dS dt = Λ− β(a, t)SI − µS, dE dt = β(a, t)SI − (σ + µ)E, dI dt = σE − (γ + µ)I, dR dt = γI − µR. (1) Here, Λ is the recruitment rate, µ is the natural death rate, σ is the transition rate from the exposed class to the infected class, γ is the transition rate from the infected class to the recovered class, β(a, t) represents the transmission rate, which varies seasonally according to the effect of the household model: S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 3 of 12 β(a, t) = β0 ( 1 + ϵ cos ( 2πt T )) , where ϵ captures the strength of seasonal effects (e.g., school vacations), and T represents the yearly period. By adding equations from 1-4 of the model (1), we get Figure 1: Flow Chart for model (1) dN dt = Λ− µN dN dt + µN = Λ =⇒ N = Λ µ + ( No − Λ µ ) exp(−µt) as t → ∞, we get N = Λ µ , which implies that the given system has a unique solution. 3. Existence of Equilibrium Points At equilibrium, the time derivatives are zero. Setting dS dt = 0, dE dt = 0, dI dt = 0, dR dt = 0, the equilibrium equations become Λ− β(a, t)SI − µS = 0, β(a, t)SI − (σ + µ)E = 0, σE − (γ + µ)I = 0, γI − µR = 0. Solving these equations, we get a smoking-free equilibrium point D0(S0, E0, I0, R0), where S0 = Λ µ . S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 4 of 12 For an endemic equilibrium point, we have I > 0. The given system has an existence of unique positive equilibrium D∗(S ∗, E∗, I∗, R∗), where S∗ = (σ + µ)(γ + µ) βσ , E∗ = µ(γ + µ)(R0 − 1) σβ , I∗ = µ(R0 − 1) β , R∗ = γ(R0 − 1) β . So there exists an endemic equilibrium point for R0 − 1 > 0. 3.1. Stability and Basic Reproduction Number The basic reproductive number, R0, is calculated as the spectral radius of the next-generation matrix: R0 = βσλ µ(σ + µ)(γ + µ) . 3.2. Local Stability Analysis Theorem 1. Smoking-free equilibriums are locally asymptotically stable when R0 < 1. Proof. The smoking-free equilibrium (SFE) is given by (S0, E0, I0, R0) = (Λµ , 0, 0, 0). The Jacobian matrix evaluated at SFE is as follows: J =  −µ 0 −βS0 0 0 −(σ + µ) βS0 0 0 σ −(γ + µ) 0 0 0 γ −µ  . (2) For eigenvalues, we solve |J − λ′I| = 0, |J − λ′I| = ∣∣∣∣∣∣∣∣ −µ− λ′ 0 −βS0 0 0 −(σ + µ)− λ′ βS0 0 0 σ −(γ + µ)− λ′ 0 0 0 γ −µ− λ′ ∣∣∣∣∣∣∣∣ = 0. Expanding by C1 (−µ− λ′) ∣∣∣∣∣∣ −(σ + µ)− λ′ βS0 0 σ −(γ + µ)− λ′ 0 0 γ −µ− λ′ ∣∣∣∣∣∣ = 0, again expand with C1 =⇒ ( − µ− λ′)[(− (σ + µ)− λ′)(((γ + µ) + λ′)(µ+ λ′) ) − σ ( − βS0(µ+ λ′) )] = 0, =⇒ ( − µ− λ′)[(− (σ + µ)− λ′)(λ′2 + µλ′ + µ(γ + µ) + λ′(γ + µ) ) + σβS0µ+ σβS0λ′ ] = 0, S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 5 of 12 =⇒ ( − µ− λ′)(− λ′3 − (γ + 2µ)λ′2 − µ(γ + µ)λ′ − (σ + µ)λ′2 − (σ + µ)(γ + 2µ)λ′ − (σ + µ)µ(γ + µ) + σβS0µ+ σβS0λ′ ) = 0. Thus, it follows the following equation, =⇒ ( − µ− λ′)[λ′3 + (γ + 3µ+ σ)λ′2 + ( µ(γ + µ) + (σ + µ)(γ + 2µ)− βS0σ ) λ′ + (σ + µ)µ(γ + µ)− βS0σµ ] = 0, From the following inequalities: =⇒ βS0σ < (σ + µ)(γ + µ), =⇒ βS0σ < (σ + µ)(γ + 2µ), it follows that R0 < 1, which implies that the characteristic equation associated with the linearized system at the smoking-free equilibrium has negative roots. Moreover, the corresponding characteristic polynomial satisfies the Routh-Hurwitz stability criteria under the condition R0 < 1. Therefore, all eigenvalues of the Jacobian matrix have negative real parts, confirming that the smoking-free equilibrium is locally asymptotically stable. Theorem 2. The positive (endemic) equilibrium is locally asymptotically stable if R0 > 1. Proof. The endemic equilibrium is given by (S∗, E∗, I∗, R∗) = ( (σ + µ)(γ + µ) βσ , µ(γ + µ)(R0 − 1) σβ , µ(R0 − 1) β , γ(R0 − 1) β ) . The Jacobian matrix evaluated at the endemic equilibrium point is J =  −βI∗ − µ 0 −βS∗ 0 βI∗ −(σ + µ) βS∗ 0 0 σ −(γ + µ) 0 0 0 γ −µ  . To determine the eigenvalues, we solve the characteristic equation |J − λ′I| = 0, which becomes |J − λ′I| =  −βI∗ − µ− λ′ 0 −βS∗ 0 βI∗ −(σ + µ)− λ′ βS∗ 0 0 σ −(γ + µ)− λ′ 0 0 0 γ −µ− λ′  = 0. Expanding along the fourth column C4, we obtain (−µ− λ′) ∣∣∣∣∣∣ −βI∗ − µ− λ′ 0 −βS∗ βI∗ −(σ + µ)− λ′ βS∗ 0 σ −(γ + µ)− λ′ ∣∣∣∣∣∣ = 0. This yields one eigenvalue λ′ = −µ. Next, we expand the remaining 3× 3 determinant along the first column C1, leading to:( − βI∗ − µ− λ′)(((σ + µ) + λ′)((γ + µ) + λ′)− βσS∗ ) − β2σI∗S∗ = 0, S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 6 of 12 ( − βI∗ − µ− λ′)(λ′2 + ( (σ + µ) + (γ + µ) ) λ′ + (σ + µ)(γ + µ)− βσS∗ ) − β2σI∗S∗ = 0, λ′3 + ( (σ + µ) + (γ + µ) + βI∗ + µ ) λ′2 + (( (σ + µ)(γ + µ)− βσS∗)+ (βI∗ + µ) ( (σ + µ) + (γ + µ) )) λ′ +(βI∗ + µ) ( (σ + µ)(γ + µ)− βσS∗)+ β2σI∗S∗ = 0. Given that S∗ = (σ + µ)(γ + µ) βσ , the above equation simplifies to the reduced form: λ′3 + ( σ + 3µ+ γ + βI∗ ) λ′2 + (βI∗ + µ) ( σ + 2µ+ γ ) λ′ + β2σI∗S∗ = 0, where I∗ = µ(R0 − 1) β , for R0 > 1. Since all the coefficients of the above cubic polynomial are positive forR0 > 1, the Routh–Hurwitz criteria are satisfied. Therefore, all roots of the characteristic equation have negative real parts. This confirms that the endemic equilibrium is locally asymptotically stable. 4. Global Stability Analysis Using Lyapunov Function 4.1. Smoking-Free Equilibrium Point Theorem 3. Smoking-free equilibrium point is globally asymptotically stable for R0 < 1. Proof. To show global stability at smoking-free equilibrium point, we build a Lyapunov function: V = ∫ A 0 ( E2 2 + I2 2 )da. (3) This function is 1) positive definite, i-e V (t) > 0, when E = I = 0. 2) Zero only at the equilibrium point E = I = 0. 3) Differentiable in time. Taking the time derivative: dV dt = ∫ A 0 (E dE dT + I dI dt )da. By substituting the model (1) equations, we get dV dt = ∫ A 0 ( E ( βSI − (σ + µ)E ) + I ( σE − (γ + µ)I )) da, dV dt = ∫ A 0 ( β(a, t)SEI − (σ + µ)E2 + σEI − (γ + µ)I2 ) da, dV dt = ∫ A 0 ( λ(a, t)SE + σEI − (σ + µ)E2 − (γ + µ)I2 ) da, where λ = β(a, t)I dV dt = − ∫ A 0 ((σ + µ)E2 + (γ + µ)I2)da+ ∫ A 0 (λ(a, t)SE + σEIda). S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 7 of 12 Now, we analyze the sign of dV dt . Using Young’s inequality ab ≤ a2 2 + b2 2 for cross terms. λSE ≤ λ2S2 2ϵ1 + ϵ1E 2 2 , and σEI ≤ σ2E2 2ϵ2 + ϵ2I 2 2 . So dV dt ≤ − ∫ A 0 ((σ + µ)E2 + (γ + µ)I2)da+ ∫ A 0 ( λ2S2 2ϵ1 + ϵ1E 2 2 + σ2E2 2ϵ2 + ϵ2I 2 2 ) da. By rearranging dV dt ≤ − ∫ A 0 (( σ + µ− ϵ1 2 − σ2 2ϵ2 ) E2 − ( γ + µ− ϵ2 2 ) I2 ) da+ ∫ A 0 λ2S2 2ϵ1 da. If smoking prevalence is small, i-e λ is small, i-e λ(t) → 0 as I(t) → 0, then λ2S2 is small near the equilibrium point, so the last term of the above equation is negligible. Choose ϵ1, ϵ2 such that σ − ϵ1 2 − σ2 2ϵ2 > 0, and γ − ϵ2 2 > 0, then dV dt ≤ 0. The equality holds only when E = I = 0. According to Lasalle’s principle of invariance, the smoking-free equilibrium point is globally asymptotically stable. 4.2. Positive Equilibrium point Theorem 4. Positive equilibrium point is globally asymptotically stable. Proof. To show the global stability of the endemic equilibrium point, we construct a Lya- punov function V (t) = 1 2 (E − E∗)2 + 1 2 (I − I∗)2, =⇒ dV dt = (E − E∗) dE dt + (I − I∗) dI dt . From the model (1) equations dE dt = βSI − (σ + µ)E, and dI dt = σE − (γ + µ)I, =⇒ dV dt = (E − E∗) ( βSI − (σ + µ)E ) + (I − I∗) ( σE − (γ + µ)I ) , =⇒ dV dt = (E − E∗) ( (σ + µ)E∗ − (σ + µ)E ) + (E − E∗)(βSI − βS∗I∗) S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 8 of 12 +(I − I∗)(σE − σE∗) + (I − I∗) ( (γ + µ)I∗ − (γ + µ)I ) , =⇒ dV dt = −(σ + µ)(E − E∗)2 + (E − E∗)(βSI − βS∗I∗) + σ(I − I∗)(E − E∗)− (I − I∗)2(γ + µ), =⇒ dV dt = −(σ + µ)(E − E∗)2 − (I − I∗)2(γ + µ) + (E − E∗)(βSI − βS∗I∗) + σ(I − I∗)(E − E∗). Use Young’s inequality on cross terms, we have (E − E∗)(I − I∗) ≤ 1 2 (E − E∗)2 + 1 2 (I − I∗)2, and (E − E∗)(λS − λ∗S) ≤ 1 2ϵ (λS − λ∗S∗)2 + ϵ 2 (E − E∗)2. Thus we can bound dV dt by dV dt ≤ −(σ + µ)(E − E∗)2 − (I − I∗)2(γ + µ) + 1 2ϵ(λS − λ∗S∗)2 +σ 2 (E − E∗)2 + σ 2 (I − I∗)2, =⇒ dV dt ≤ −( σ 2 + µ)(E − E∗)2 − (γ + µ− σ 2 )(I − I∗)2 + 1 2ϵ (λS − λ∗S∗)2. If (λS−λ∗S∗) is small near the equilibrium and if the parameters satisfy γ+µ > σ 2 , then dV dt < 0 near the equilibrium. This implies global asymptotic stability of the endemic equilibrium. S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 9 of 12 5. Numerical Simulations We solve the system numerically using non-standard MATLAB, employing the Euler method and varying β(a, t) to simulate household effects. The results illustrate the impact of sea- sonal variations on the prevalence. The initial conditions are E(0)=50, I(0)=30, R(0)=20 and P(0)=N-E(0)-I(0)-R(0). Parameters values are given in below table; Parameter values Λ 10.25 β0 0.005 σ 0.1 γ 0.03 ϵ 0,05 N 250 Table 1: Parameter values for the proposed model (a) Comparison of Susceptible popu- lation with RK4, Euler, Nonstandard finite difference method (b) Comparison of exposed population with RK4, Euler, Nonstandard finite difference method (c) Comparison of infected population with RK4, Euler, Nonstandard finite difference method (d) Comparison od recovered popula- tion with RK4,Euler, Nonstandard finite difference method (e) population graph with RK4 (f) Population graph with Eu- ler (g) Population graph with NSFD Figure 2: (a) plot of susceptible class in which RK4, Euler and Nonstandard finite difference. (b) plot of infected class in which RK4, Euler and Nonstandard finite difference. (c) plot of exposed class in which RK4, Euler and Nonstandard finite difference. (d) plot of recovered class in which RK4, Euler and Nonstandard finite difference (e) plot of population with RK4 method. (e) plot of population with Euler method. (f) plot of population with nonstandard finite difference method. For different value of b(a, t) the graph of all population by using RK4, Euler, and nonstan- dard finite difference method are given below. S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 10 of 12 (a) Susceptible population using the RK4 technique with various numbers of β(a, t) (b) Infected pop- ulation using the RK4 technique with various numbers of β(a, t) (c) Exposed pop- ulation using the RK4 technique with various numbers of β(a, t) (d) Recovered pop- ulation using the RK4 technique with various numbers of β(a, t) (e) Susceptible pop- ulation using the Eu- ler technique with various numbers of β(a, t) (f) Exposed popula- tion using the Eu- ler technique with various numbers of β(a, t) (g) Infected popula- tion using the Eu- ler technique with various numbers of β(a, t) (h) Recovered popu- lation using the Eu- ler technique with various numbers of β(a, t) (i) Susceptible pop- ulation using the NSFD technique with various num- bers of β(a, t) (j) Exposed popula- tion using the NSFD technique with various numbers of β(a, t) (k) Infected popula- tion using the NSFD technique with various numbers of β(a, t) (l) Recovered pop- ulation using the NSFD technique with various num- bers of β(a, t) Figure 3: (a) plot of susceptible class in which RK4 method is used for different value of β(a, t). (b)plot of exposed class in which RK4 method is used for different value of β(a, t). (c) plot of infected class in which RK4 method is used for different value of β(a, t) (d) plot of recovered class in which RK4 method is used for different value of β(a, t). (e) plot of susceptible class in which Euler method is used for different value of β(a, t) . (f) plot of exposed class in which Euler method is used for different value of β(a, t). (g) plot of infected class in which Euler method is used for different value of β(a, t). (h) plot of recovered class in which Euler method is used for different value of β(a, t). (i) plot of susceptible class in which nonstandard finite difference method is used for different value of β(a, t).(j) plot of exposed class in which nonstandard finite difference method is used for different value of β(a, t). (k) plot of infected class in which nonstandard finite difference method is used for different value of β(a, t).(l) plot of recovered class in which nonstandard finite difference method is used for different value of β(a, t). 6. Conclusions and Discussions The current paper investigated the global asymptotic dynamics of an age-structured smoking model. Our initial step involved formulating the model and verifying that it admits a unique solution. We then analyzed the equilibrium points, including the smoking-free and endemic equilibria. To establish the global asymptotic stability of these points, we constructed an ap- propriate Lyapunov function. The analytical results were supported by numerical simulations, S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 11 of 12 which validated the model behavior for different parameter values. The impact of the incre- mental rate was examined through graphical results, highlighting the sensitivity of the model to behavioral transitions across age groups. Compared to previous models in the literature, our age-structured formulation adds a signif- icant layer of realism, especially in modeling differential susceptibility and recovery rates among distinct age cohorts. This provides a refined framework for analyzing the spread and persistence of smoking behaviors in populations. In line with the reviewer’s recommendation, this work opens several promising avenues for future research. One natural extension is to incorporate **relapse dynamics**, allowing for- mer smokers to return to smoking under certain behavioral or social conditions. Moreover, integrating **intervention strategies**, such as public health campaigns or medical treatments, could enhance the practical relevance of the model. Additionally, the inclusion of **stochastic effects** would enable the modeling of randomness and uncertainty in smoking behaviors and policy outcomes, which is particularly important for small populations or individual-level mod- eling. These enhancements can lead to more comprehensive insights and provide policymakers with robust tools to design targeted anti-smoking strategies. References [1] R. M. Anderson and R. M. May. Infectious Diseases of Humans: Dynamics and Control. Oxford University Press, 1991. [2] H. W. Hethcote. The mathematics of infectious diseases. SIAM Review, 42(4):599–653, 2000. [3] F. Brauer, C. Castillo-Chavez, and Z. Feng. Mathematical Models in Epidemiology. Springer, 2019. [4] O. Diekmann, J. A. P. Heesterbeek, and T. Britton. Mathematical Tools for Understanding Infectious Disease Dynamics. Princeton University Press, 2013. [5] L. G. de Pillis and A. E. Radunskaya. A mathematical tumor model with immune resistance and drug therapy: an optimal control approach. Journal of Theoretical Biology, 238(4):841– 862, 2006. [6] M. Gümüs and Ö. Öcalan. Global asymptotic stability of a nonautonomous difference equation. Journal of Applied Mathematics, 2014:Article ID 395954, 5 pages, 2014. [7] M. Y. Li and J. S. Muldowney. Global stability for the seir model in epidemiology. Math- ematical Biosciences, 125(2):155–164, 1995. [8] J. Wang, M. Zhang, and W. Yang. Dynamics and asymptotic profiles of a local-nonlocal dispersal sir epidemic model with spatial heterogeneity. Journal of Differential Equations, 313:1–34, 2022. [9] C. Liu and J. Shi. Traveling wave solution for a delayed reaction-diffusion two-group sir epidemic model with a generalized nonlinear incidence function. Nonlinear Analysis: Real World Applications, 65:103434, 2022. [10] E. Alzahrani and A. Zeb. Stability analysis and prevention strategies of tobacco smoking model. Boundary Value Problems, 2020:1–13, 2020. [11] A. Zeb, S. Kumar, A. Tesfay, and A. Kumar. A stability analysis on a smoking model with stochastic perturbation. International Journal of Numerical Methods for Heat and Fluid Flow, 32(3):915–930, 2022. [12] S. M. Lestari, Y. Yulida, A. S. Lestia, and M. A. Karim. Stability analysis of a smoking behavior model. Journal of Physics: Conference Series, 2106(1):012025, 2021. [13] A. M. Pulecio-Montoya, L. E. López-Montenegro, and L. M. Benavides. Analysis of a mathematical model of smoking. Contemporary Engineering Sciences, 12(3):117–129, 2019. [14] Z. Alkhudari, S. Al-Sheikh, and S. Al-Tuwairqi. The effect of occasional smokers on the dynamics of a smoking model. International Mathematical Forum, 9(25):1207–1222, 2014. S. Bushnaq, S. Bano, A. Zeb / Eur. J. Pure Appl. Math, 18 (3) (2025), 6447 12 of 12 [15] G. Zaman and S. Gul. Qualitative behavior of giving up smoking models. Bulletin of the Malaysian Mathematical Sciences Society. Second Series, 34(2):403–415, 2011. [16] Z. Alkhudhari, S. Al-Sheikh, and S. Al-Tuwairqi. Global dynamics of a mathematical model on smoking. International Scholarly Research Notices, 2014(1):847075, 2014. [17] A. Zeb, G. Zaman, S. Momani, and S. Shaher. Square-root dynamics of a giving up smoking model. Applied Mathematical Modelling, 37(7):5326–5334, 2013. [18] A. Zeb, A. Bano, E. Alzahrani, and G. Zaman. Dynamical analysis of cigarette smoking model with a saturated incidence rate. AIP Advances, 8(4), 2018. [19] Y. S. Tan, X. X. Li, J. Yang, and R. A. Cheke. Global stability of an age-structured model of smoking and its treatment. International Journal of Biomathematics, 16(01):2250063, 2023. [20] A. Zeb, S. Hussain, O. J. Algahtani, and G. Zaman. Global aspects of age-structured cigarette smoking model. Mediterranean Journal of Mathematics, 15:1–11, 2018. [21] E. Addai, A. Adeniji, O. J. Peter, J. O. Agbaje, and K. Oshinubi. Dynamics of age- structure smoking models with government intervention coverage under fractal-fractional order derivatives. Fractal and Fractional, 7(5):370, 2023. [22] S. Djilali, S. Bentout, A. Zeb, S. Saeed, and T. Tareq. Global stability of hybrid smoking model with nonlocal diffusion. Fractals, 30(08):2240224, 2022. [23] A. Lahrouz, L. Omari, D. Kiouach, and A. Belmaâti. Deterministic and stochastic stability of a mathematical model of smoking. Statistics & Probability Letters, 81(8):1276–1284, 2011. [24] A. Zeb, G. Zaman, I. H. Jung, and M. Khan. Optimal campaign strategies in fractional- order smoking dynamics. Zeitschrift für Naturforschung A, 69(5-6):225–231, 2014.