EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6473 ISSN 1307-5543 – ejpam.com Published by New York Business Global Optimal Control Strategies in a Nonlinear Smoking Cessation Model Sameera Bano1, Shifa Siddiqui1, Anwar Zeb1,∗, Thoraya N. Alharthi2, Ilyas Khan3,4,5, Osama Oqilat6, Wei Sin Koh7 1 Department of Mathematics, COMSATS University Islamabad, Abbottabad Campus, KPK, Pakistan 2 Department of Mathematics, College of Science, University of Bisha, P.O. Box 551, Bisha 61922, Saudi Arabia 3 Department of Mathematical Sciences, Saveetha School of Engineering, SIMATS, Chennai, Tamil Nadu, India 4 Hourani Center for Applied Scientific Research, Al-Ahliyya Amman University, Amman, Jordan 5 Department of Mathematics, College of Science Al-Zulfi, Majmaah University, Al-Majmaah 11952, Saudi Arabia 6 Hourani Center for Applied Scientific Research, Department of Basic Sciences, Faculty of Arts and Science, Al-Ahliyya Amman University, Amman, Jordan 7 INTI International University, Persiaran Perdana BBN Putra Nilai, 71800 Nilai, Negeri Sembilan, Malaysia Abstract. This study introduces a novel nonlinear smoking cessation model with optimal control strategies and time delays. The model incorporates a convex incidence rate and divides the total population into four classes. We first establish positivity, compute equilibrium points, and deter- mine the basic reproduction number using the next generation method. Global stability is analyzed via the Castillo–Chavez principle for the smoking-free equilibrium and the geometric approach for the endemic equilibrium. Sensitivity analysis is performed to identify influential parameters. The model is extended to include optimal control, where the existence and uniqueness are established, and the impact of time delays is examined. A Hopf bifurcation is shown to occur at a critical delay, indicating oscillatory behavior. Numerical simulations validate the theoretical results and demon- strate improved performance over existing methods in reducing smoking prevalence and associated costs. 2020 Mathematics Subject Classifications: 34H05, 92D30, 93C15 Key Words and Phrases: Convex incidence rate, Castillo-Chavez principle, geometric approach, sensitivity analysis, optimal control, time-delay ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6473 Email addresses: sameera@cuiatd.edu.pk (S. Bano), swadah827@gmail.com (S. Siddiqui), anwar@cuiatd.edu.pk (A. Zeb), talhrthe@ub.edu.sa (T. N. Alharthi), i.said@mu.edu.sa (I. Khan) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 2 of 38 1. Introduction Research in mathematical biology was first started in the 19th century; Initially, it has a slow progress, but later at the end of the 20th century it became more attractive. John led the development of mathematical biology in 1909 and in 1912, he gave the basic law of epidemic spread. Kermack and Mckendrick proposed the first dynamic model in 1927 [1]. Later, his work was further extended by Andrea, who applied a time delay on the model [2]. Many models proposed based on the SIR model, like SEIR [3] presented by NH Shah, other SEIR [4] which include delay information and feedback mechanisms influencing susceptible individuals’ behavior, this model analysis the global stability at both equilibrium points. SIS model [5] investigate model within the time scales framework, analyzing its dynamics and incorporating births and deaths. In mathematical biology models incidence rate plays very significant role. Yorke and London present a model [6] on the outbreak of measles, chickenpox, and mumps. In this model, they use an incidence rate of the form fIS = IS(1 − cI). Liu and his co-workers present a model in which they use the contact rate of the form βIaS (1+αIc) , where β, α, c, and a > 0. A saturated incidence rate model of cholera was proposed by Capasso and Serio, which also highlights the importance of the non-linear contact rate for the spread of the disease. Some other models related to nonlinear incidence rates like nonmonotone incidence rate [7], Harmonic mean type incidence rate [8], Saturated incidence rate [9] and some other model with different incidence rate [10] are given in the literature review. Another interesting contact rate is the convex contact rate, which is expressed as g(s) = βSI(1 + αI). This incidence rate shows the increase rate of infection due to double exposure. βSI is the contact rate produced due to single interaction, while βαSI2 is due to the newly formed infection due to double exposure. A. khan developed a model of COVID-19 [11] and hepatitis B [12] with a convex incidence rate. Other models with convex incidence rate are [13], [14], [15]. Recent studies have used mathematical models to explore how infectious diseases like rabies affect the brain, focusing on transmission dynamics and control strategies. This reflects the expanding role of mathematical biology in tackling neurological health issues [16]. Researcher often use incidence functions to model the spread of infectious diseases [17– 19], but they can also be used to model the spread of behaviors like smoking. Smoking has a knack for influencing the non-smokers and quit smokers to start smoking, making smoking contagious like infectious disease. In 500 BC in America, it was used in rituals and ceremonies. It was brought to Europe in the 15th and 16th centuries. In 1880 James A. Bonsack invented the first smoking machine that increases the smoking rate. It causes 5 million deaths around the world, making it the largest cause of death. The harmful chemical enter our bodies by smoking and affect our hearts, nervous systems, respira- tory system, stomach, etc. Among teenagers, as they grow and develop, smoking leads to headaches, dizziness, and memory loss. Due to the strong addiction to smoking only 20-22 percent of people can quit smoking. In 1997 Castillo Garsow first established a smoking model in which he divides the total community into 3 classes P,L and S that are potential, regular, and quit smokers [20]. Later, Castillo’s work was expanded by Zamman [21], in S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 3 of 38 which he applies two control measures. Okhyar proposed a model in which he divides the population of quit smokers into temporary and permanent quit smokers [22]. Some other smoking-related model includes the works of J.Singh, who developed a fractional model [23], O sharmoni, who presented a model in which quit smokers are divided into tempo- rary and permanent quit smokers and he also added the control strategies to analyze most effective intervention to reduce the smoking impact [24], MI. Chuhan, who developed a smoking model in which homotopic analysis method is used [25]. Smoking remains a leading cause of preventable deaths globally. Current work is related to a new smoking model with intervention-based time delay and controls not previously explored to fill the gap. We offer theoretical insights and numerical validations, with sensitivity analysis and optimization-based techniques to validate efficacy. The main con- tributions of this work are as follows: Formulation of a novel delay-free optimal control model. Stability analysis using Lyapunov and sensitivity theory. Derivation and verification of effective control strategies through numerical simulations. The rest of the paper is organized as follows: Section 2 develops the new model formula- tion. Section 3 provides a theoretical analysis. Section 4 discusses the control problem. Section 5 is related the delay model with stability analysis, Section 6 shows the numerical simulations and Section 7 provides the discussion of model, control and delay applied on system and finally in Section 8, we present the conclusion of this work. 2. Model formation In the proposed model, the total community was divided into four distinct classes, which are PLSQ. Here, P stands for potential smokers, L stands for light smokers, S stands for regular smokers, and Q stands for quit smokers. N is used to denote the total population as N = P + L + S + Q. The smoking cessation model with convex incidence rate is given as: dP dt = Λ− (b+ µ)P − βPS(1 + νS), dL dt = βPS(1 + νS)− (b+ µ+ ζ)L, dS dt = ζL− (b+ µ+ δ)S, dQ dt = δS − (b+ µ)Q. (1) The initial conditions are P0 ≥ 0, L0 ≥ 0, S0 ≥ 0, Q0 ≥ 0. Here β and ν are the smoking transmission coefficient, b represents the death rate due to smoking-related complications. While it is generally associated with active smokers, we also include it in the P (potential smokers) class to capture individuals who may already be affected by early-stage health issues (such as heart or lung problems) or second-hand smoke exposure. This reflects real-world situations where people are at risk due to environmental S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 4 of 38 Figure 1: The graph of the proposed model is given above. exposure or the smoking behavior of friends and family, even before they themselves begin smoking. Λ The continuous influx of the potential smokers (Vulnerable smokers) ν Constant rate β Constant rate ζ The transition rate from light smokers to full smoker δ The smoking cessation rate µ Natural mortality rate b Smoking related mortality rate Table 1: Description of the parameter S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 5 of 38 As the total population is equal to the sum of all classes, this implies that N = P + L+ S +Q, dN dt = Λ− (b+ µ)N. By using the integrating factor e(b+µ)t and after some computation, we get N = Λ (b+µ) as t → ∞. This implies that all the system trajectory is bounded within a feasible region Ω. Ω = {(P,L, S,Q) : 0 ≤ N ≤ Λ (b+µ)}. (2) 3. Positivity and Boundedness Since positivity can be assumed due to non-negative initial conditions and biological feasibility. 3.1. Existence of equilibrium point In this subsection, we discuss the equilibrium points of the PLSQ model. There are two equilibrium points of the model (1) that are smoking-free and endemic equilibrium points. 3.2. Smoking-free equilibrium point To determine the equilibrium points of system (1), we set the right-hand sides of the differential equations (1) equal to zero. This corresponds to the state where all time derivatives vanish, indicating no change in the population compartments over time. In case of a smoking-free equilibrium point, it shows that there is no case of smoking that is L=S=0 and after using the value of S we get Q = 0 and P having value which is Λ b+µ . Smoking-free equilibrium point of the proposed model is E0 = (P0, L0, S0, Q0) = ( Λ (b+ µ) , 0, 0, 0 ) . 3.3. Smoking endemic equilibrium point The model (1) endemic equilibrium point entries are given below: L∗ = (b+µ+δ)S∗ ζ , Q∗ = δS∗ (b+µ) , P ∗ = (b+µ+ζ)(b+µ+δ) βζ(1+νS∗) , S∗ = √ e2−4fg−e 2f . Here, f = νβζ(b+ µ+ δ)(b+ µ+ ζ), e = ζβ(b+ µ+ δ)(b+ µ+ ζ)− νβΛζ2, g = ζ(b+ µ)(b+ µ+ δ)(b+ µ+ ζ)− Λβζ2. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 6 of 38 3.4. Reproductive number In this subsection, R0 is computed with the help of the next-generation method. Ac- cording to this method, the reduced system is written in the form of dY dt = F − V. (3) Here, Y consists of three classes; potential smokers, light smokers, and regular smokers. F contains the transmission rate ( the rate at which the individual starts smoking) and V includes the transition rate ( the rate at which the individual moves between the smoking stages and exists from these compartments) and these techniques are used to compute R0. Here, F = −βPS(1 + νS) βPS(1 + νS) 0  , V =  −Λ + (b+ µ)P (b+ µ+ ζ)L −ζL+ (b+ µ+ δ)S  . The Jacobian of F and V at E0 are: F = 0 0 − βΛ (b+µ) 0 0 βΛ (b+µ) 0 0 0  , V = (b+ µ) 0 0 0 (b+ µ+ ζ) 0 0 −ζ (b+ µ+ δ)  , FV −1 = 0 − βΛζ (µ+b)(µ+b+ζ)(µ+b+δ) − βΛ (µ+b)(µ+b+δ) 0 βΛζ (µ+b)(µ+b+ζ)(µ+b+δ) βΛ (µ+b)(µ+b+δ) 0 0 0  . As R0 is the dominant eigenvalue of FV −1, first 2 eigen values are zero and the third eigen value is βΛζ r1r2r3 . So, R0 = βΛζ r1r2r3 . (4) Here, r1 = (µ+b), r2 = (µ+b+ζ), r3 = (µ+b+δ). F and V represent matrices in R3×3; the vectors used in the next-generation matrix approach are Fi, Vi. We include below a concrete example using assumed parameters to demonstrate that R0 is the dominant eigenvalue of |FV −1 − λI|. Entries (1,2) and (2,2) show this, while (1,3) and (2,3) do not contribute dominantly. Now the smoking endemic equilibrium point E∗ = (P ∗, L∗, S∗, Q∗) in term of the reproductive number are as given: P ∗ = Λ R0r1(1 + νS∗) , L∗ = r3S ∗ ζ , Q∗ = δS∗ r1 , S∗ = √ e2 − 4fg − e 2f . Here, f = νβ2ζ2Λ R0(µ+b) , e = Λβζ2 [ β R0(µ+b) − ν ] , g = Λβζ2 [ 1 R0 − 1 ] . S∗ is positive for R0 > 1. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 7 of 38 3.5. Stability In this section, local as well as global stability of the equilibrium points are evaluated. 3.5.1. Local stability For local stability of the equilibrium points, we will first linearized the model by using Taylor series method. Now let y1 = P − P ∗, y2 = L − L∗, y3 = S − S∗, y4 = Q − Q∗ as the system(1) first three equation are independent of the Q, so we use only first three equation of proposed model. Linearized system is given as: ẏ1 = [ − βS∗(1 + νS∗)− (b+ µ) ] y1 + [ − βP ∗(1 + 2νS∗)]y3, ẏ2 = [ βS∗(1 + νS∗) ] y1 − [ (b+ µ+ ζ) ] y2 + [ βP ∗(1 + 2νS∗) ] y3, ẏ3 = ζy2 − (b+ µ+ δ)y3. Jacobian matrix corresponding to the linearized system is J = −βS∗(1 + νS∗)− r1 0 −βP ∗(1 + 2νS∗) βS∗(1 + νS∗) −r2 βP ∗(1 + 2νS∗) 0 ζ −r3  . (5) Theorem 1. The smoking-free equilibrium point is locally stable for R0 < 1 and unstable if R0 > 1 . Proof. Jacobian matrix at E0 is JE0 = −r1 0 −βΛ r1 0 −r2 βΛ r1 0 ζ −r3  , The first eigen value of the JE0 is λ1 = −r1 and for second and third eigen value we have an equation that is λ2 + (r2 + r3)λ + (r2r3 − βΛζ r1 ) = 0, which is solve by using Routh Horwitz criteria. Clearly, (r2 + r3) > 0 and (r2r3 − βΛζ r1 ) > 0 if R0 < 1. It is proven that the smoking-free equilibrium points are locally asymptotically stable for R0 < 1. Theorem 2. The Smoking endemic equilibrium point are locally stable for R0 > 1 and unstable for R0 < 1. Proof. Jacobian matrix at smoking endemic equilibrium point is JE∗ = −βS∗(1 + νS∗)− r1 0 −βP ∗(1 + 2νS∗) βS∗(1 + νS∗) −r2 βP ∗(1 + 2νS∗) 0 ζ −r3  . (6) After some row operation applied on JE∗ , we get the two eigen values that are λ1 = −r2 < 0, λ2 = −ζr1 < 0. For third eigen value, we have λ3 = βζP ∗(1 + 2νS∗)r1 − r1r2r3 − r2r3βS ∗(1 + νS∗), S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 8 of 38 after putting the value of P ∗ λ3 = r2r3S ∗[νr1 − β(1 + νS∗)2 (1 + νS∗) ] . =⇒ r1ν < β(1 + νS∗)2, where r1 is sum of parameters and these prameter values are between 0 and 1 and β(1 + νS∗)2 contain the square of S∗ which is positive value for R0 > 1 above inequality hold which proves that third eigen value is negative. =⇒ λ3 < 0. As all eigenvalues are negative, it is proven that the endemic equilibrium points are locally asymptotically stable for R0 > 1. 3.5.2. Global stability To proof that the smoking-free and endemic equilibrium points are globally asymptotically stable, we will use the Castillo Chavez principle and the geometric approach [26]-[27], re- spectively. Global stability at smoking free equilibrium point To apply the Castillo Chavez principle, we divide the system (1) into dK dt = E(K,W ), dW dt = F (K,W ). Here, K represents the nonsmokers population (P,Q) and W represents the smoker population (L,S). E0 = (K0, 0̄) is the smoking-free equilibrium point. According to the Castillo Chavez principle if the following two conditions are fulfilled, then smoking free equilibrium point is globally stable at R0 < 1. N1 dK dt = E(K, 0), K0 is globally stable . N2 The formula for F=F(K,W) is F (K,W ) = YW − F̄ (K,W ), where F̄r(K,W ) ≥ 0∀ (K,W ) in feasible region Ω for r=1,2. Y is M-matrix, which is defined as Y = DWF (K0, 0̄), F̄r(K,W ) is the individual component of F̄ (K,W ), DW is the partial derivative of the sub matrix F (K,W ) with respect to smoker class (L,S) and Y = DWF (K0, 0̄) is the Jacobian of the sub matrix F of smoker classes at smoking free equilibrium point. Lemma 1. Smoking-free equilibrium point E0 = (K0, 0̄) is globally asymptotically stable, provided that above conditions are satisfied. Theorem 3. At smoking-free equilibrium point, proposed model (1) is globally asymptot- ically stable if R0 < 1. Proof. As the system (1) is divided into two sub systems E and F which are defined as; dK dt = E(K,W ) = ( Λ− βPS(1 + νS)− (b+ µ)P δS − (b+ µ)Q ) , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 9 of 38 dW dt = F (K,W ) = ( βPS(1 + νS)− (b+ µ+ ζ)L ζL− (b+ µ+ δ)S ) , The smoking-free equilibrium point is E0 = (K0, 0̄), where K0 = ( Λ (b+µ) , 0) and 0̄ = (0, 0). For P = P0 and Q = Q0, dK dt = ( Λ− (b+ µ) Λ (b+µ) 0 ) = 0, (7) dK dt = Λ− (b+ µ)K, Now by using integrating factor e(b+µ)t, we get K → K0 as t → ∞. This proves the condition (N1) that is K0 is globally stable for subsystem dK dt = E(K, 0). For condition (N2) we first defined Y matrix for (K,W ) = (K0, 0̄) Y = [ −(b+ µ+ ζ) βP0(1 + 2νS0) ζ −(b+ µ+ δ) ] F̄ (K,W ) = [ −(b+ µ+ ζ) βP0(1 + 2νS0) ζ −(b+ µ+ δ) ] · [ L S ] − [ −(b+ µ+ ζ) βP (1 + νS) ζ −(b+ µ+ δ) ] · [ L S ] = [ βP0S(1 + 2νS0)− βPS(1 + νS) 0 ] . Now F̄1 ≥ 0, if βP0S(1 + 2νS0)− βPS(1 + νS) ≥ 0. This implies that βP0S(1 + 2νS0) ≥ βPS(1 + νS). And F̄2 = 0. F(K,W) is the sub matrix which contain the smoker population and F̄ (K,W ) is obtain after rearranging the formula in the condition N2 that is F̄ (K,W ) = YW − F (K,W ). Here, W is the smoker class (L,S) and F̄1(K,W ) and F̄2(K,W ) are the entries of the F̄ (K,W ). Hence, both conditions (N1, N2) are satisfied, this proves the global stability of the smoking-free equilibrium point. Global stability at smoking present equilibrium point Now, to check the global stability at the endemic equilibrium point, we apply geometric approach, which was introduced by Li and Muldowney. The lemma and condition for this approach is given as S1 System has a unique equilibrium . S2 Compact absorbing set exists . S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 10 of 38 Lemma 2. If the conditions S1, S2 are met and Bendixson criteria are met then equilirium point is globally asymptotically stable. For the proof of Bendixson criteria, we first suppose the matrix-valued function P such that its inverse exists. P is defined as: P = ( n 2 ) × ( n 2 ) , Bendixson criteria is given as q2 = ∫ T 0 (µ(B)dt) < 0, where the B matrix is defined as B = PfP 1 + PJ2P−1, µ(B) is the Lozinski measure of matrix B and J2 is the second additive compound matrix of Jacobian. Lozinski measure of B matrix with respect to the vector norm ||.|| is defined as: µ(B) = lim h→0 |I +Bh| − 1 h . For system Lozinski measure can also be expressed as µ(B) ≤ sup{si, i = 1, 2}, si = ∥Bij∥+ µ1(Bii), where (i ̸= j, i, j = 1, 2). Above explanation is the general introduction of the Geometric approach. Lemma 3. Let us assume Q1 is a simply connect ste and both conditions S1 and S2 hold; if q2 < 0, then the endemic equilibrium points are globally asymptotically stable. Clarification of above Lemma: Ω = {(P,L, S,Q) such that 0 ≤ P + L+ S +Q ≤ Λ b+µ} and Q1 = {(P,L, S) such that 0 ≤ P + L + S ≤ Λ b+µ} is subset of Ω it is compact absorbing set, simply connected and E∗ is the unique equilibrium point. For the Bendixson criteria we will prove the below theorem 4. Theorem 4. Endemic equilibrium point is globally stable if R0 > 1 and βS(1 + νS) > βPS(1+2νS) L . Proof. The second additive compound matrix for Jacobian is ( considering only first 3 classes PLS) J2 = A1 βP (1 + 2νS) βP (1 + 2νS) ζ A2 0 0 βS(1 + νS) −(b+ µ+ ζ)− (b+ µ+ δ)  , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 11 of 38 where A1 = −βS(1+νS)− (b+µ)− (b+µ+ζ), A2 = −βS(1+νS)− (b+µ)− (b+µ+δ). The P function is defined as P = diag(1, LS , L S ) , now PfP −1 = diag ( 0, L̇L − Ṡ S , L̇ L − Ṡ S ) . To find matrix B, we have the formula B = PfP −1 + PJ2P−1. B = ( B11 B12 B21 B22 ) , where B11 = A1, B12 = ( SβP (1+2νS) L SβP (1+2νS) L ) , B21 = ( ζL S 0 ) , B22 = ( w11 w12 w21 w22 ) . Here, w11 = L̇ L − Ṡ S +A2, w12 = 0, w21 = βS(1 + νS), w22 = L̇ L − Ṡ S − (b+ µ+ ζ)− (b+ µ+ δ). The norm is defined as; ∥(n1, n2, n3)∥ = max(∥n1∥, ∥n2∥+ ∥n3∥). µ(B) ≤ sup{si, i = 1, 2}, (8) si = ∥Bij∥+ µ1(Bii) where (i ̸= j, i, j = 1, 2), µ1(B11) = A1, µ1(B11) = −βS(1 + νS)− (b+ µ)− (b+ µ+ ζ), ∥B12∥ = SβP (1 + 2νS) L , ∥B21∥ = ζL S . µ1(B22) = max { L̇ L − Ṡ S +A2 + βS(1 + νS), L̇ L − Ṡ S − (b+ µ+ ζ)− (b+ µ+ δ) } , after using the value of A2, we get: µ1(B22) = L̇ L − Ṡ S − (b+ µ+ δ)−min{(b+ µ), (b+ µ+ ζ)}, µ1(B22) = L̇ L − Ṡ S − (b+ µ+ δ)− (b+ µ). using all above values in si, where i=1,2 s1 = −βS(1 + νS)− (b+ µ)− (b+ µ+ ζ) + SβP (1 + 2νS) L , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 12 of 38 Using equation 2 of system (1) in s1 L̇ L = βPS(1 + νS) L − (b+ µ+ ζ), =⇒ s1 ≤ −βS(1 + νS)− (b+ µ) + L̇ L + SβP (1 + 2νS) L . s2 = L̇ L − Ṡ S − (b+ µ+ δ)− (b+ µ) + ζL S , use equation 3 of system (1) in above equation s2 Ṡ S = ζL S − (b+ µ+ δ), =⇒ s2 = L̇ L − (b+ µ), put values in equation (8) µ(B) ≤ L̇ L − (b+ µ) + sup{−βS(1 + νS) + SβP (1 + 2νS) L , 0}, we can deduce that if βS(1 + νS) > SβP (1 + 2νS) L , (9) then µ(B) ≤ L̇ L − (b+ µ), q2 = lim t→∞ sup sup 1 t ∫ t 0 l(B) dL < −(b+ µ). (10) The Bendixson criteria is satisfied as q2 < 0. this implies that endemic equilibrium point (P ∗, L∗, S∗) is globally asymptotically stable. From equation 4 of system (1) dQ dt = δS − (b+ µ)Q, limit equation is dQ dt = δS∗ − (b+ µ)Q, from integrating factor e(b+µ)t and let t → ∞ Q(t) → Q∗. Hence, the endemic equilibrium point E∗ is globally asymptotically stable. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 13 of 38 3.6. Sensitivity In this part of the paper, the sensitivity analysis is conducted. Sensitivity analysis is an essential part of infectious disease modeling because it facilitates the researcher to identify which parameter has the greatest influence on the propagation and control of the infection. The forward index method is used to check sensitivity, which involves computing the partial derivative of R0 with respect to each parameter. The sensitivity is checked by using the formula that is xR0 p = ∂R0 ∂p × p R0 . Here p ∈ (Λ, β, ζ, δ, µ, b), xR0 Λ = 1 > 0, xR0 β = 1 > 0, xR0 ζ = (b+ µ) (b+ µ+ ζ) > 0, xR0 δ = −δ (b+ µ+ δ) < 0, xR0 µ = −µ(3µ2 + 3b2 + 2bζ + 2bδ + 6µb+ 2µδ + 2µζ + δζ) (b+ µ)(b+ µ+ ζ)(b+ µ+ δ) < 0, xR0 b = −b(3d2 + 6µb+ 2bδ + 2bζ + 2ζµ+ 2µδ + ζδ + 3µ2) (b+ µ)(b+ µ+ ζ)(b+ µ+ δ) < 0. (11) From the above computation, it is shown that β, ζ, Λ has positive influence on R0 and the index of b, δ, µ show that increasing their value decreases the value of R0. Parameter Sensitive index formula Index values Base value used Source Λ xR0 Λ 1 1.025 assumed β xR0 β 1 0.00038 [28] ζ xR0 ζ 0.3805 0.21 [29] δ xR0 δ -0.5686 0.017 assumed b xR0 b -0.2669 0.0019 [29] µ xR0 µ -1.5451 0.011 [29] Table 2: Sensitive index of the R0 Figure 2 of the sensitivity analysis is constructed by using Matlab. Figure 2 shows that µ, b, δ has a negative impact on disease propagation. By increasing these parameter values will reduce the spread of smoking. The parameter values that are used in the sensitivity analysis are given below. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 14 of 38 Figure 2: Sensitivity index graph. 4. Optimal control We have extended the dynamics to include new compartments and control functions, making it a novel contribution. In this part of the paper, we employ the method from the optimal control theory to forecast the decrease and analyze the eradication of smoking habits in the population. To achieve this, we modified the model (1) with optimal control by using three control measures that are the media campaign for the potential smokers, antismoking gum for light smokers, and medicine for regular smokers. In [30] and [31] researchers found that the media campaign is very helpful in reducing the use of tobacco among adults and young generations and decreases the smoking rate. The modified model of the optimal control problem is given below. dP dt = Λ− (b+ µ)P − βPS(1 + νS)− a0u1P + π1u2L+ π2u3S, dL dt = βPS(1 + νS)− (b+ µ+ ζ + u2)L, dS dt = ζL− (b+ µ+ δ + u3)S, dQ dt = δS − (b+ µ)Q+ a0u1P + (1− π1)u2L+ (1− π2)u3S. , (12) with initial condition P0 ≥ 0, L0 ≥ 0, S0 ≥ 0, Q0 ≥ 0 The goal of the control strategies is to decrease the population of light smokers and regular smokers and enhancing the population of quit smokers and potential smokers, which can be accomplished by reducing the contact between smokers and nonsmokers. The cost (objective function) of the above control problem is given bellow J = ∫ T 0 (D1L+D2S −D3P −D4Q+ 1 2(B1u 2 1 +B2u 2 2 +B3u 2 3))dt. In the above objective function (D1, D2, D3, D4) are positive constant for the classes P,L,S,Q respectively. The second term includes (B1, B2, B3) are the positive weight param- eters associated with the control variables (u1, u2, u3). The objective of the cost function S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 15 of 38 is to decrease smoking and increase quit smokers. U is the admissible control set and is defined as U = {u = (u1, u2, u3)|uj(t), 0 ≤ uj(t) ≤ uj(t) max ≤ 1, t ∈ [0, T ], j = 1, 2, 3} In the above set U, uj(t) are the control variables for j=1,2,3 and set U is the Lebesgue measurable function. In above set U, (umax 1 , umax 2 , and umax 3 ) denotes the maximum achievable values of u1, u2, and u3 respectively, and T is the terminal time. The model (12) all solutions lie in the feasible region Ω, i-e. Ω = { (P,L, S,Q) ∈ R4, 0 ≤ N ≤ Λ b+ µ } . 4.1. Optimal Control Existence Theorem 5. There exists an optimal control (u∗1, u ∗ 2, u ∗ 3) that minimizes the objective function J. Proof. To prove the existence of optimal control, we must satisfied the following condition: (i) Both the control set and the state variable are non empty. (ii) The control set is closed as well as convex. (iii) Objective function integrand is convex. (iv) r(t, x, u) = A(t, x) +B(t, x)u. (v) There exist Y1 > 0, positive numbers Y2, Y3, the integrand of the objective function fulfill. L ≥ Y2(|u1|2 + |u2|2 + |u3|2) Y1 2 − Y3, where Y1 = 2, Y2 = 1 2min(B1, B2, B3), Y3 = B3u3 2 . (i) To proof the above conditions, we use the result by Luke [28] to establish the ex- istence solution of system with bounded coefficient, which verified the Condition 1. The system can be expressed as R(x) = Y x+O(x), where state variables are x(t) = (P (t), L(t), S(t), Q(t))T and matrix Y and O(x) nonlinear function are defined as below: Y =  −(b+ µ+ a0u1) π1u2 π2u3 0 0 −(b+ µ+ ζ + u2) 0 0 0 ζ −(b+ µ+ δ + u3) 0 a0u1 (1− π1)u2 (1− π2)u3 + δ −(b+ µ) , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 16 of 38 O(x) =  Λ− βPS(1 + νS) βPS(1 + νS) 0 0  , let set x1 = (P1, L1, S1, Q1), x2 = (P2, L2, S2, Q2) done in similar manner as it is done in [32] |R(x1)−R(x2)| ≤ ||Y |||x1 − x2|+ |O(x1)−O(x2)|, O(x1)−O(x2) = ( −βP1S1(1 + νS1) + βP2S2(1 + νS2) βP1S1(1 + νS1)− βP2S2(1 + νS2) ) , Now, |O(x1)−O(x2)| = | − βP1S1(1 + νS1) + βP2S2(1 + νS2)| +|βP1S1(1 + νS1)− βP2S2(1 + νS2)|, |O(x1)−O(x2)| ≤ (2β|S1|+ 2νβ|S2 1 |)|P1 − P2| +(2β|P2|+ 2νβ|P2||S1 + S2|)|S1 − S2|, |O(x1)−O(x2)| ≤ 2βΛ (b+µ) ( 1 + νΛ (b+µ) )( |P1 − P2|+ |S1 − S2| ) , |R(x1)−R(x2)| ≤ H|x1 − x2|. Here, H = max{( 2βΛ (b+µ) [1 + νΛ (b+µ) ], ||Y ||} < ∞. This show that function R(x) is uniformly Lipchitz continuous (satisfy the assumption of condition 1). U control set is nonempty because it contain the measurable function and bounded by 0 ≤ uj(t) ≤ umax j where j = 1, 2, 3 for example constant zero control u(t)=(0,0,0). (ii) As control set U is closed by definition and each component of U lies in [0,1] and for any u1 = (u11, u 1 2, u 1 3) ∈ U and u2 = (u21, u 2 2, u 2 3) ∈ U and q ∈ [0, 1] and u = ru1 + (1 − r)u2 for each j, 0 ≤ u1j ≤ umax j ≤ 1, 0 ≤ u2j ≤ umax j ≤ 1 this shows the convex combination defined by 0 ≤ qu1j + (1 − q)u2j ≤ umax j ≤ 1, where j=1,2,3 remain in [0,1], this means u ∈ U . This implies that U is convex. This proves the condition (2) (iii) To prove the integrand function is convex, we let λ ∈ [0, 1], s ∈ (s1, s2, s3) and a ∈ (a1, a2, a3) , x ∈ (P,L, S,Q) and we need to show that L(x, (1− λ)s+ λa)) ≤ (1− λ)L(x, s) + λL(x, a) L(x, (1− λ)s+ λa))− (1− λ)L(x, s)− λL(x, a) S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 17 of 38 = D1L+D2S −D3P −D4Q+ B1 2 ( (1− λ)s1 + λa1 )2 + B2 2 ( (1− λ)s2 + λa2 )2 + B3 2 ( (1− λ)s3 + λa3 )2 − ( 1− λ )( D1L+D2S −D3P −D4Q ) − ( 1− λ )B1 2 s21 − ( 1− λ )B2 2 s22 − ( 1− λ )B3 2 s23 − λ ( D1L+D2S −D3P −D4Q ) −λ B1 2 a21 − λ B2 2 a22 − λ B3 2 a23, = B1 2 ( (1− λ)2s21 + λ2a21 + 2(1− λ)λs1a1 − (1− λ)a21 − λa21 ) + B2 2 ( (1− λ)2s22 + λ2a22 + 2(1− λ)λs2a2 − (1− λ)s22 − λa22 ) + B3 2 ( (1− λ)2s23 + λ2a23 + 2(1− λ)λs3a3 − (1− λ)s23 − λa23 ) , = (λ− 1)(λ) ( B1 2 (s1 − a1) 2 + B2 2 (s2 − a2) 2 + B3 2 (s3 − a3) 2 ) ≤ 0. =⇒ L(x, (1− λ)s+ λa)) ≤ (1− λ)L(x, s) + λL(x, a). Condition (3) holds. (iv) Now x = (P,L, S,Q)T , u = (u1, u2, u3) T . Now for the proof condition 4, we separate the system (12) into sub matrix A(t,x) and B(t,x), where A(t,x) is the matrix of dynamic without control, B(t,x) is the state dependent control coefficient matrix. A(t, x) =  Λ− βPS(1 + νS)− (b+ µ)P βPS(1 + νS)− (b+ µ+ ζ)L ζL− (b+ µ+ δ)S δS − (b+ µ)Q  , A(t, x) is the matrix that represent the part of system that depend only on state vari- able and independent of control and it reflect natural behavior of system, including its inherent nonlinear behavior. B(t, x) =  −a0P π1L π2S 0 −L 0 0 0 −S a0P (1− π1)L (1− π2)S  , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 18 of 38 B(t, x) shows the influence of control on system. The system (12) is linear in control, and can be rewritten as: r(t,x,u)=A(t,x)+B(x,t)u dx dt = r(t, x, u) =  Λ− βPS(1 + νS)− (b+ µ)P βPS(1 + νS)− (b+ µ+ ζ)L ζL− (b+ µ+ δ)S δS − (b+ µ)Q  +  −a0P π1L π2S 0 −L 0 0 0 −S a0P (1− π1)L (1− π2)S  · u1u2 u3  . Condition (4) holds. (v) As integrand of the objective function is L = D1L+D2S −D3P −D4Q+ 1 2 (B1u 2 1 +B2u 2 2 +B3u 2 3) For the proof of condition 5, we use (P,L, S,Q) ≥ 0 and we assume that weighted sum of the smoker class(L,S) are greater than that of nonsmoker class(P,Q) that is D1L+D2S > D3P +D4Q, where Di > 0, i = 1, 2, 3, 4 L > 1 2 (B1u 2 1 +B2u 2 2 +B3u 2 3). As B3u3 2 ≤ B3 2 L ≥ Y2(|u1|2 + |u2|2 + |u3|2) Y1 2 − Y3. Y1 = 2, Y2 = 1 2min(B1, B2, B3), Y3 = B3u3 2 . This prove the condition (5). Hence all the five condition are proof, this implies that optimal control exists that minimizes the objective function. 4.2. Optimal control Characteristic of the PLSQ system The Hamiltonian function is given as: H = L+ ∑4 n=1 λnhn. Here hn = (dPdt , dL dt , dS dt , dQ dt ), where n = 1, 2, 3, 4. and x = [P,L, S,Q]T , L = D1L+D2S −D3P −D4Q+ 1 2(B1u 2 1 +B2u 2 2 +B3u 2 3), Pontraygin’s maximum principle conditions are λ̇n = −∂H ∂x , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 19 of 38 ∂H ∂u = 0, λn(T ) = 0, where x = (P,L, S,Q)T and n = 1, 2, 3, 4. Theorem 6. If u∗1, u∗2, u∗3 are the control pair and P ∗, L∗, S∗, Q∗ are the control path and λn are the adjoint variables then the following conditions should be satisfied λ̇1 = D3 − λ1(−βS(1 + νS)− (b+ µ)− a0u1)− λ2(βS(1 + νS))− λ4(a0u1), λ̇2 = −D1 − λ1(π1u2)− λ2(−(b+ µ+ ζ + u2))− λ3ζ − λ4((1− π1)u2), λ̇3 = −D2 − λ1(−βP (1 + 2νS) + π2u3)− λ2(βP (1 + 2νS)) − λ3(−(b+ µ+ δ + u3))− λ4(δ + (1− π2)u3), λ̇4 = D4 − λ4(−(b+ µ)), (13) along with the transversality conditions: λn(T ) = 0 where n = 1, 2, 3, 4. Furthermore, the optimal pair are given as u∗1 = max ( min ( (λ1 − λ4)a0P ∗ B1 , umax 1 ) , 0 ) , u∗2 = max ( min ( (λ4 − λ1)π1L ∗ + (λ2 − λ4)L ∗ B2 , umax 2 ) , 0 ) , u∗3 = max ( min ( (λ4 − λ1)π2S ∗ + (λ3 − λ4)S ∗ B3 , umax 3 ) , 0 ) . (14) Proof. By applying the first condition of Pontryagin Maximum principle, we get λ̇1 = −∂H ∂P , λ̇2 = −∂H ∂L , λ̇3 = −∂H ∂S , λ̇4 = −∂H ∂Q , and transversality condition of the Pontryagin Maximum principle λn(T ) = 0, which gives λ1(T ) = 0, λ2(T ) = 0, λ3(T ) = 0, λ4(T ) = 0. From the 2 condition of Pontryagin’s maximum principle ∂H ∂u1 = 0, at u1 = u∗1, =⇒ u∗1 = (λ1 − λ4)a0P ∗ B1 , ∂H ∂u2 = 0, at u2 = u∗2 =⇒ u∗2 = (λ4 − λ1)π1L ∗ + (λ2 − λ4)L ∗ B2 , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 20 of 38 ∂H ∂u3 = 0, at u3 = u∗3 =⇒ u∗3 = (λ4 − λ1)π2S ∗ + (λ3 − λ4)S ∗ B3 . Using the control set property, we get: u∗1 =  0, if (λ1−λ4)a0P ∗ B1 ≤ 0, (λ1−λ4)a0P ∗ B1 , if 0 < (λ1−λ4)a0P ∗ B1 < umax 1 , umax 1 , if (λ1−λ4)a0P ∗ B1 ≥ umax 1 . u∗2 =  0, if (λ4−λ1)π1L∗+(λ2−λ4)L∗ B2 ≤ 0, (λ4−λ1)π1L∗+(λ2−λ4)L∗ B2 , if 0 < (λ4−λ1)π1L∗+(λ2−λ4)L∗ B2 < umax 2 , umax 2 , if (λ4−λ1)π1L∗+(λ2−λ4)L∗ B2 ≥ umax 2 . u∗3 =  0, if (λ4−λ1)π2S∗+(λ3−λ4)S∗ B3 ≤ 0, (λ4−λ1)π2S∗+(λ3−λ4)S∗ B3 , if 0 < (λ4−λ1)π2S∗+(λ3−λ4)S∗ B3 < umax 3 , umax 3 , if (λ4−λ1)π2S∗+(λ3−λ4)S∗ B3 ≥ umax 3 . This completes the proof. 4.3. Uniqueness of optimal control system The optimal system consists of the state system (12) with initial condition, adjoint system (13) with transversality condition, and optimality condition (14). Thus, optimality system is given as dP dt = Λ− βPS(1 + νS)− (b+ µ)P − a0u1P + π1u2L+ π2u3S, dL dt = βPS(1 + νS)− (b+ µ+ ζ + u2)L, dS dt = ζL− (b+ µ+ δ + u3)S, dQ dt = δS − (b+ µ)Q+ a0u1P + (1− π1)u2L+ (1− π2)u3S, λ̇1 = D3 − λ1(−βS(1 + νS)− (b+ µ)− a0u1) −λ2(βS(1 + νS))− λ4(a0u1), λ̇2 = −D1 − λ1(π1u2)− λ2(−(b+ µ+ ζ + u2))− λ3ζ −λ4((1− π1)u2), λ̇3 = −D2 − λ1(−βP (1 + 2νS) + π2u3)− λ2(βP (1 + 2νS)) −λ3(−(b+ µ+ δ + u3))− λ4(δ + (1− π2)u3), λ̇4 = D4 − λ4(−(b+ µ)). , (15) where P (0), L(0), S(0), Q(0) ≥ 0 and λi(T ) = 0, for i = 1, 2, 3, 4. u1 = max ( min ( (λ1 − λ4)a0P B1 , umax 1 ) , 0 ) , S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 21 of 38 u2 = max ( min ( (λ4 − λ1)π1L+ (λ2 − λ4)L B2 , umax 2 ) , 0 ) , u3 = max ( min ( (λ4 − λ1)π2S + (λ3 − λ4)S B3 , umax 3 ) , 0 ) . Figure 3: The flowchart of the optimality system is given above. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 22 of 38 The adjoint system has a linear bounded coefficient in each of the adjoint variables because the state system is bounded. Thus, the upper bound of the adjoint system is finite. In order to prove the uniqueness of the optimality system for small time, we state the following statement Lemma 4. The function v∗(r) = max{min{r, c}, d} is Lipchitz continuous in r with c < d for some non-negative constant. The uniqueness of the above model is evaluated in the similar way it is found in [33], [34]. Theorem 7. For a sufficiently short time interval T, the solution of the optimal system (15) is unique. Proof. Contradiction is used to prove this theorem. Contradiction is obtained by taking into account optimum control function Lipschitz continuity, state system, and adjoint system boundedness, and some few basic inequalities. Let the two solutions of system (15) be (P,L, S,Q) and (P̄ , L̄, S̄, Q̄). It is convenient to modify the variable to demonstrate that two solutions are equivalent. Now, let P = eλta1, L = eλta2, S = eλta3, Q = eλta4, P̄ = eλtā1, L̄ = eλtā2, S̄ = eλtā3, Q̄ = eλtā4, λ1 = e−λtb1, λ2 = e−λtb2, λ3 = e−λtb3, λ4 = e−λtb4, λ̄1 = e−λtb̄1, λ̄2 = e−λtb̄2, λ̄3 = e−λtb̄3, λ̄4 = e−λtb̄4. λ is the constant, whose value is chosen at the end. u1 = max ( min ( (b1−b4)a0a1 B1 , umax 1 ) , 0 ) , u2 = max ( min ( (b4−b1)π1a2+(b2−b4)a2 B2 , umax 2 ) , 0 ) , u3 = max ( min ( (b4−b1)π2a3+(b3−b4)a3 B3 , umax 3 ) , 0 ) , ū1 = max ( min ( (b̄1−b̄4)a0ā1 B1 , umax 1 ) , 0 ) , ū2 = max ( min ( (b̄4−b̄1)π1ā2+(b̄2−b̄4)ā2 B2 , umax 2 ) , 0 ) , ū3 = max ( min ( (b̄4−b̄1)π2ā3+(b̄3−b̄4)ā3 B3 , umax 3 ) , 0 ) . Some inequality conditions that we will use later are taken from [33] (w + y)2 ≤ 2(w2 + y2), (w − w̄)(y − ȳ) ≤ (w − w̄)2 + (y − ȳ)2). (ab− āb̄)2 = (ab− āb+ āb− āb̄)2, ≤ max{2b2, 2ā2}[(a− ā)2 + (b− b̄)2], ≤ D0[(a− ā)2 + (b− b̄)2]. (c− c̄)(ab− āb̄) = (c− c̄)(ab− āb+ āb− āb̄), S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 23 of 38 = (c− c̄)(b(a− ā) + ā(b− b̄)), ≤ D1((a− ā)2 + (b− b̄)2 + (c− c̄)2). From first equation of system(15), we get λa1e λt + eλta′1 = λ− βe2λta1a3(1 + νeλta3)− (b+ µ)eλta1 −a0u1a1e λt + π1u2a2e λt + π2u3a3e λt, =⇒ (λ+ b+ µ)a1 + a′1 = λe−λt − βa1a3e λt − νβa1a 2 3e 2λt −a0u1a1 + π1a2u2 + π2u3a3. Similarly, =⇒ (λ+ b+ µ)ā1 + ā1 ′ = λe−λt − βā1ā3e λt − νβā1ā3 2e2λt −a0u1ā1 + π1ā2ū2 + π2ū3ā3, by subtracting and integrating the above two equations, we get =⇒ (λ+ b+ µ) ∫ T 0 (a1 − ā1) 2 dt+ 1 2 (a1(T )− ā1(T )) 2 = −β ∫ T 0 ( eλt(a1a3 − ā1ā3)(a1 − ā1) ) dt − νβ ∫ T 0 ( e2λt(a1a 2 3 − ā1ā3 2)(a1 − ā1) ) dt − a0 ∫ T 0 ((u1a1 − ū1ā1)(a1 − ā1)) dt + π1 ∫ T 0 ((u2a2 − ā2ū2)(a1 − ā1)) dt + π2 ∫ T 0 ((u3a3 − ā3ū3)(a1 − ā1)) dt. (16) Now (a1a3 − ā1ā3)(a1 − ā1) ≤ B1[(a1 − ā1) 2 + (a3 − ā3) 2]. (a1a 2 3 − ā1ā3 2)(a1 − ā1) ≤ B2[(a1 − ā1) 2 + (a3 − ā3) 2]. (a1u1 − ā1ū1)(a1 − ā1) ≤ B3[(a1 − ā1) 2 + (u1 − ū1) 2]. (u1 − ū1) 2 = ( a0 B1 )2 [(b1 − b4)a1 − (b̄1 − b̄4)ā1] 2, ≤ B4[(a1 − ā1) 2 + (b1 − b̄1) 2 + (b4 − b̄4) 2]. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 24 of 38 Use the value of (u1 − ū1) in the above equation (a1u1 − ā1ū1)(a1 − ā1) ≤ B5[(a1 − ā1) 2 + (b1 − b̄1) 2 + (b4 − b̄4) 2]. (a2u2 − ā2ū2)(a1 − ā1) ≤ B6[(u2 − ū2) 2 + (a1 − ā1) 2 + (a2 − ā2) 2]. (u2 − ū2) 2 = ( 1 B2 )2 [π1a2(b4 − b1) + a2(b2 − b4)− π1ā2(b̄4 − b̄1) − ā2(b̄2 − b̄4)] 2, ≤ B7[(a2 − ā2) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2 + (b4 − b̄4) 2]. Use the value of (u2 − ū2) 2,we get (a2u2 − ā2ū2)(a1 − ā1) ≤ B8[(a1 − ā1) 2 + (a2 − ā2) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2 + (b4 − b̄4) 2]. (a3u3 − ā3ū3)(a1 − ā1) ≤ B9[(u3 − ū3) 2 + (a1 − ā1) 2 + (a3 − ā3) 2], (u3 − ū3) 2 = ( 1 B3 )2 [π2a3(b4 − b1) + a3(b3 − b4)− π2ā3(b̄4 − b̄1) − ā3(b̄3 − b̄4)] 2, ≤ B10[(a3 − ā3) 2 + (b1 − b̄1) 2 + (b3 − b̄3) 2 + (b4 − b̄4) 2]. Using the value of (u3 − ū3) 2, we get (a3u3 − ā3ū3)(a1 − ā1) ≤ B11[(a1 − ā1) 2 + (a3 − ā3) 2 + (b1 − b̄1) 2 + (b3 − b̄3) 2 + (b4 − b̄4) 2]. Bi where i=1,2,3,...11 depend on the bound of the state variable. Now substitute the vales in equation (16) =⇒ (λ+ b+ µ) ∫ T 0 (a1 − ā1) 2 dt+ 1 2 (a1(T )− ā1(T )) 2 ≤ M1e λT ∫ T 0 ( (a1 − ā1) 2 + (a3 − ā3) 2 ) dt +K1e 2λT ∫ T 0 ( (a1 − ā1) 2 + (a3 − ā3) 2 ) dt + P1 ∫ T 0 ( (a1 − ā1) 2 + (a2 − ā2) 2 + (a3 − ā3) 2 + (b1 − b̄1 2 ) +(b2 − b̄2) 2 + (b3 − b̄3) 2 + (b4 − b̄4) 2 ) dt. (17) Similarly =⇒ (λ+ b+ µ+ ζ) ∫ T 0 (a2 − ā2) 2 dt+ 1 2 ( a2(T )− ¯a2(T ) )2 S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 25 of 38 ≤ M2e λT ∫ T 0 ( (a1 − ā1) 2 + (a2 − ā2) 2 + (a3 − ā3) 2 ) dt +K2e 2λT ∫ T 0 ( (a1 − ā1) 2 + (a2 − ā2) 2 + (a3 − ā3) 2 ) dt + P2 ∫ T 0 ( (a2 − ā2) 2 + (b1 − b̄1) + (b2 − b̄2) + (b4 − b̄4) 2 ) dt. (18) =⇒ (λ+ b+ µ+ δ) ∫ T 0 (a3 − ā3) 2 dt+ 1 2 (a3(T )− ā3(T )) 2 ≤ P3 ∫ T 0 ( (a2 − ā2) 2 + (a3 − ā3) 2 + (b1 − b̄1) 2 +(b3 − b̄3) 2 + (b4 − b̄4) 2 ) dt. (19) =⇒ (λ+ b+ µ) ∫ T 0 (a4 − ā4) 2 dt+ 1 2 (a4(T )− ā4(T )) 2 ≤ P4 ∫ T 0 ( (a1 − ā1) 2 + (a2 − ā2) 2 + (a3 − ā3) 2 +(a4 − ā4) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2 +(b3 − b̄3) 2 + (b4 − b̄4) 2 ) dt. (20) =⇒ (λ+ b+ µ) ∫ T 0 (b1 − b̄1) 2 dt+ 1 2 ( b1(0)− b̄1(0) )2 ≤ M5e λT ∫ T 0 ( (a3 − ā3) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2) ) dt +K5e 2λT ∫ T 0 ( (a3 − ā3) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2) ) dt + P5 ∫ T 0 ( (a1 − ā1) 2 + (b1 − b̄1) 2 + (b4 − b̄4) 2 ) dt. (21) =⇒ (λ+ b+ µ+ ζ) ∫ T 0 (b2 − b̄2) 2 dt+ 1 2 ( b2(0)− b̄2(0) )2 ≤ P6 ∫ T 0 ( (a2 − ā2) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2 +(b3 − b̄3) 2 + (b4 − b̄4) 2 ) dt. (22) =⇒ (λ+ b+ µ+ δ) ∫ T 0 (b3 − b̄3) 2 dt+ 1 2 ( b3(0)− b̄3(0) )2 S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 26 of 38 ≤ M7e λT ∫ T 0 ( (a1 − ā1) 2 + (b1 − b̄1) 2 + (b2 − b̄2) 2 +(b3 − b̄3) 2 ) dt +K7e 2λT ∫ T 0 ( (a1 − ā1) 2 + (a3 − ā3) 2 + (b1 − b̄1) 2 +(b2 − b̄2) 2 + (b3 − b̄3) 2 ) dt + P7 ∫ T 0 ( (a3 − ā3) 2 + (b1 − b̄1) 2 + (b3 − b̄3) 2 +(b4 − b̄4) 2 ) dt, (23) (λ+ b+ µ) ∫ T 0 (b4 − b̄4) 2 dt+ 1 2 ( b4(0)− b̄4(0) )2 + ā = 0. (24) Here, (Mi,Kj , Pr), where (i = 1, 2, , 5, 7), (j = 1, 2, 3, ..7) and r = 1, 2, 5, 7 depend on the bounds and coefficients of the state variable and the costate variable. From equation (17)-(24)[ (λ+ b+ µ)− (M1 +M2 +M7)e λT − (K1 +K2 +K7)e 2λT −(P1 + P4 + P5)] ∫ T 0 ( (a1 − ā1) 2 ) dt + [ (λ+ b+ µ+ ζ)−M2e λT −K2e 2λT −(P1 + P2 + P3 + P4 + P6)] ∫ T 0 ( (a2 − ā2) 2 ) dt + [ (λ+ b+ µ+ δ)− (M1 +M2 +M5)e λT − (K1 +K2 +K5 +K7)e 2λT −(P1 + P3 + P4 + P7)] ∫ T 0 ( (a3 − ā3) 2 ) dt + [(λ+ b+ µ)− P4] ∫ T 0 ( (a4 − ā4) 2 ) dt + [ (λ+ b+ µ)− (M5 +M7)e λT − (K5 +K7)e 2λT −(P1 + P2 + P3 + P4 + P5 + P6 + P7)] ∫ T 0 ( (b1 − b̄1) 2 ) dt + [ (λ+ b+ µ+ ζ)− (M5 +M7)e λT − (K5 +K7)e 2λT −(P1 + P2 + P4 + P6)] ∫ T 0 ( (b2 − b̄2) 2 ) dt + [ (λ+ b+ µ+ δ)−M7e λT −K7e 2λT S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 27 of 38 −(P1 + P3 + P4 + P6 + P7)] ∫ T 0 ( (b3 − b̄3) 2 ) dt + [(λ+ b+ µ)− (P1 + P2 + P3 +P4 + P5 + P6 + P7)] ∫ T 0 ( (b4 − b̄4) 2 ) dt ≤ 0. (25) The coefficient of each integral of equation (25) are non negative if we choose sufficiently small T and sufficiently large λ. For example, if we fix λ > M1 + M2 + M7 + K1 + K2 + K7 + P1 + P4 + P5 − b − µ and T < 1 2λ ln ((λ+b+µ)−(P1+P4+P5)) M1+M2+M7+K1+K2+K7 . Then the coefficient of integral ∫ T 0 (a1 − ā1) 2 dt is nonnegative. The same method is applied to the remaining integral terms. We obtain the others λs, Ts. Take the maximum of the all λs and minimum of the Ts, then used it as λ and T, then all the integral coefficients will be non-negative. This implies that ā1 = a1, ā2 = a2, ā3 = a3, ā4 = a4, b̄1 = b1, b̄2 = b2, b̄3 = b3, b̄4 = b4, P̄ = P, L̄ = L, S̄ = S, Q̄ = Q, λ̄1 = λ1, λ̄2 = λ2, λ̄3 = λ3, λ̄4 = λ4. Hence, the solution of equation (15) is unique. 5. Delay We modify the classical PLSQ model to incorporate intervention-based dynamics and time delay. Unlike previous models, our structure captures the relapse and partial immu- nity. The modified model is given below dP dt = Λ− (b+ µ)P − βPS(1 + νS), dL dt = βPS(1 + νS)− (µ+ b)L− ζL(t− τ), dS dt = ζL(t− τ)− (µ+ b+ δ)S, dQ dt = δS − (µ+ b)Q. (26) Here, τ is the delay parameter. Now we linearize the delay system by choosing the variable as a = P − P ∗, b = L−L∗, c = S − S∗. As the above system is independent of the Q, so the linearized system is given as ȧ = [−βS∗(1 + νS∗)− (b+ µ)]a+ [−βP ∗(1 + 2νS∗)]c, ḃ = [ βS∗(1 + νS∗) ] a− [ b+ µ+ ζe−λτ ] b+ [ βP ∗(1 + 2νS∗) ] c, ċ = [ ζe−λτ ] b− [ b+ µ+ δ ] c. The Jacobian of the linearized system is J = −βS∗(1 + νS∗)− (b+ µ) 0 −βP ∗(1 + 2νS∗) βS∗(1 + νS∗) −(b+ µ+ ζe−λτ ) βP ∗(1 + 2νS∗) 0 ζe−λτ −(b+ µ+ δ) . (27) S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 28 of 38 The characteristic equation from above equation is λ3 + λ2(d1 + f1e −λτ ) + λ(d2 + f2e −λτ ) + (d3 + f3e −λτ ) = 0. (28) The values of the coefficients of the characteristic equation are as given d1 = −k1 − k4 − k8, d2 = k1k4 + k4k8 + k1k8, d3 = −k1k4k8, f1 = −k5, f2 = k1k5 + k5k8 − k6k7, f3 = −k1k5k8 + k1k6k7 − k2k3k7, k1 = −βS∗(1 + νS∗)− (b+ µ), k2 = −βP ∗(1 + 2νS∗), k3 = βS∗(1 + νS∗), k4 = −(b+ µ), k5 = −ζ, k6 = βP ∗(1 + 2νS∗), k7 = ζ, k8 = −(b+ µ+ δ). 5.1. Stability analysis The stability of the delay system of the smoking cessation model with convex incidence rate is computed in the same way as already done in [35]. We will check the two cases that are (1) when delay parameter τ = 0, (2) when delay parameter τ > 0. 5.1.1. In the absence of delay In this case, we put τ = 0 in equation (28) =⇒ λ3 + λ2(d1 + f1) + λ(d2 + f2) + d3 + f3 = 0, let, D2 = d1 + f1, D1 = d2 + f2, D0 = d3 + f3. Then equation becomes =⇒ λ3 + λ2D2 + λD1 +D0 = 0. According to the Routh Hurwitz criteria, if the below conditions hold, then the equilibrium points are locally stable. det1 = D2 > 0, det2 = ( D2 D0 1 D1 ) > 0, det3 = D2 D0 0 1 D1 0 0 D2 D0  > 0. (29) 5.1.2. Presence of delay In this case we let τ > 0 and λ = ιω where ω > 0. After using the value λ equation (28) becomes ((ιω)3 + (ιω)2d1 + (ιω)d2 + d3) + ((ιω)2f1 + (ιω)f2 + f3)e −ιωτ = 0. As e−ιθ = cos θ − ι sin θ, after some simple calculation we get cosωτ = ω4Y1 + ω2Y2 + Y3 ω4Y4 + ω2Y5 + Y6 . (30) sinωτ = ω5Y7 + ω3Y8 + ωY9 ω4Y4 + ω2Y5 + Y6 . (31) S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 29 of 38 Here, Y1 = f2 − d1f1, Y2 = d1f3 + d3f1 − d2f2, Y3 = −d3f3, Y4 = f2 1 , Y5 = f2 2 −2f1f3, Y6 = f2 3 , Y7 = f1, Y8 = d1f2−d2f1−f3, Y9 = d2f3−d3f2. After utilizing cos2 ωτ + sin2 ωτ = 1, we get ω10U1 + ω8U2 + ω6U3 + ω4U4 + ω2U5 + ωU6 = 0, where, U1 = Y 2 7 , U2 = Y 2 1 + 2Y7Y8 − Y 2 4 , U3 = 2Y1Y2 + Y 2 8 + 2Y7Y9 − 2Y4Y5, U4 = Y 2 2 + 2Y1Y3 + 2Y8Y9 − Y 2 5 − 2Y4Y6, U5 = 2Y2Y3 + Y 2 9 − 2Y5Y6, U6 = Y 2 3 − Y 2 6 . Let suppose ω2 = t t5U1 + t4U2 + t3R3 + t2U4 + tU5 + U6 = 0. (32) One of the positive roots of the equation(32) is ω0 and for ω0 there exists τn0 = 1 ω0 (cos−1 g1(ω0) + 2nπ), where, τ0 = min{τn0 , n = 0, 2, 3..}. g1(ω0) = ω4 0(f2 − d1f1) + ω2 0(d1f3 + d3f1 − d2f2)− d3f3 ω4 0d 2 1 + ω2 0(f 2 2 − 2f3f1) + f2 3 . After differentiating characteristic equation[ ∂λ ∂τ ]−1 = (3λ2 + 2d1λ+ d2) + e−λτ (2f1λ+ f2) λe−λτ (f1λ2 + f2λ+ f3) − τ λ . Now after putting λ = ιω0, we get[ ∂λ ∂τ ]−1 |τ=τ0 = E + Fι G+Hι + ιτ ω , ℜ [ ∂λ ∂τ ]−1 |τ=τ0 = EG+ FH G2 +H2 ̸= 0. Here, E = d2−3ω2 0+f2 cosω0τ+2f1ω0 sinω0τ, F = 2ω0d1+2ω0f1 cosω0τ−f2 sinω0τ, G = −ω2 0f2 cosω0τ−ω3 0f1 sinω0τ+ω0f3 sinω0τ, H = −ω3 0f1 cosω0τ+ω0f3 cosω0τ+ω2 0f2 sinω0τ. Hence, Hopf bifurcation exists at τ = τ0. 6. Numerical simulation In this part of the paper, we use a numerical method to solve the optimal control and delay problem. For numerical simulation fourth order Runge Kutta method is employed. All the numerical simulations are done in MATLAB. A fixed step size of 0.1 was used throughout the simulations. For the delay model, a constant delay value of 10 days was implemented using discrete delay approximation based on past state values. The implementation covers three scenarios: the basic model, a model with delay in the light S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 30 of 38 smoker class, and a model with control interventions applied to potential, light, and regular smokers. The parameter values that are used in the simulation are presented in the table 4. The baseline parameters (β, ζ, δ, µ, ν, b) are kept constant across all three scenarios. The only difference in the controlled model is the inclusion of three intervention controls, which influence the system dynamics. Therefore, the change in behavior between models arises from the control functions, not from changes in the original parameter values. Parameter values source Λ 8 assumed β 0.0038 [28] ζ 0.021 [29] δ 0.017 assumed b 0.0019 [29] µ 0.011 [29] ν 0.0009 [13] u1 0.01 [29] u2 0.01 [29] u3 0.05 [29] a0 0.0051 [29] π1 0.085 [29] π2 0.095 [29] Table 3: Parameter values for proposed model S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 31 of 38 The evolution of each population class over the 300-day simulation period is shown in Fig. (4). Fig. 4(a) demonstrates the natural progression of each class without the effect of external intervention like delay and control. The light smokers population increases rapidly as compared to other classes due to the fast transition into light smokers, while regular and quit smokers gradually rises before stabilizing. Figure 4(b) introduces the three controls in the simulation, namely education campaign, the anti-smoking gum, and the medicine whose values are given in the above table. The control measures decreasing the spread of smoking, lowering the potential, light and regular smokers population and increasing the quit smokers population, which suggest that more people stop smoking. Figure 4(c) shows the impact of the delay (τ = 10) in the simulation. Delay affects the magnitude and timing of the transition between classes, such as delay increases light and potential smokers population as compared to the case of without delay. The rapid increase in light smokers suggests that delays prolong the time smokers spend as light smokers and postpone the transition to regular and quit smokers. Control measures are very effective in stabilizing the system and decreasing the light, regular smokers and increasing quit smokers population compared to with and without delay cases. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 32 of 38 (a) PLSQ model without control and delay (b) Control PLSQ model (c) Delay PLSQ model Figure 4: (a) Plot of proposed model without control and delay. (b) Plot of proposed model with control and without delay. (c) Plot of the proposed model without control and with delay. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 33 of 38 Fig. 5 shows the effect of control and delay on each population separately. In the potential smokers graph, without delay and control curves show the slow transition of potential smokers to other classes, while the delay curve shows that the transition from potential to other classes stabilizes faster than the simple curve, but with the introduction of controls, rapidly potential smokers move to the other classes and reach the stability more quickly than the delay and simple curve. In the case of the light smokers, with delay, the population increases rapidly and then declines as it moves to another class and stabilizes at the same level as the simple curve. This shows that delay causes a temporary rise in light smokers due to the slow transition to quit or regular smokers. while control measures (like anti smoking gums) accelerate the reduction in the light smoker population. The population of regular smokers, without the effect of delay and control population, increases rapidly before stabilization, which shows a swifter transition to the regular smoking population, but with delay, this transition slows down. A control measure (like medicine), further decreases this transition and reaches a lower peak, which shows the effectiveness of the control. Finally, for the quit smoker, there is a gradual increase in the population in all three cases, which indicates a slow growth rate, but the delay curve is slightly lower than the simple curve, showing a minor reduction in growth. On the other hand, control measures speed up the process of quitting smoking and result in the increase in the quit smokers population. In general, by implementing delay (τ = 10), transition becomes slow and prolongs the quitting process. However, control measures increase the rate of quitting and reduce the smoking rate. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 34 of 38 (a) Potential smokers (b) Regular smokers (c) Light smokers (d) Quit smokers Figure 5: (a) The potential smokers population who are vulnerable to smoking with and without control and delay. (b) The population of smokers who smoke occasionally with and without control and delay. (c) The population of smokers who smoke regularly with and without control and delay. (d) The population of individuals who quit smoking permanently with and without control and delay. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 35 of 38 7. Discussion The present study examines a smoking behavior model that incorporates both optimal control strategies and a discrete time delay in the transmission dynamics. By integrating these elements, the model captures more realistic behavioral patterns, such as the gradual transition between smoking states and the delayed impact of intervention measures. This framework extends previous models that neglected either the influence of time delays or the role of systematically designed control strategies. Numerical simulations were performed using the classical fourth-order Runge–Kutta (RK4) method implemented in MATLAB, enabling accurate and efficient approximation of the system’s trajectories over time. The computational results demonstrate that the introduction of a time delay can significantly alter system behavior, potentially inducing oscillations or changing the stability of equilibria, thereby influencing the effectiveness of intervention strategies. Moreover, the application of optimal control provides a systematic means to minimize the prevalence of smokers while balancing the associated implementa- tion costs. Compared with earlier formulations in the literature, our model delivers a more com- prehensive and practical representation of smoking behavior and its control. The combina- tion of time delay and optimal control allows for more precise predictions of intervention outcomes and facilitates the identification of critical time thresholds for policy actions. Such insights are highly valuable for public health authorities, as they inform the timing and intensity of anti-smoking campaigns, resource allocation, and awareness programs. Overall, this work enhances existing methodological approaches by integrating behav- ioral delay effects with mathematically optimized intervention strategies. The results not only strengthen the theoretical understanding of smoking dynamics but also provide ac- tionable guidance for evidence-based policymaking aimed at reducing smoking prevalence and its associated health and economic burdens. 8. Conclusion In this article, we examine the smoking cessation model considering the important feature of the convex incidence rate with and without control and delay. This model is constructed using the current literature in this field. We examine the boundedness and positive of the solution once the model was constructed, provided that the initial conditions are satisfied. Two equilibrium points of the model are determined, namely smoking-free and endemic equilibrium points, through simple computation. Then the basic reproduc- tive number is computed, which is very important in understanding the propagation of smoking. With the help of the reproductive number, sufficient conditions for local and global stability are determined. Sensitivity analysis is performed, which shows that the dynamic behavior of the system is affected by changing the values of the parameters. After that, control problem are established by utilizing three control measures. Here pontryagin maximum principle is utilize to establish the best control strategies. The model is modified with delay, where delay is treated as a bifurcation parameter. In delay model stability is S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 36 of 38 checked with help of the Routh-Hurwitz criteria. Finally, to verify the theoretical results, numerical computation is performed. The effectiveness of control and delay are shown by comparing the curves with and without control and delay. References [1] W. O. Kermack and A. G. McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115:700–721, 1927. [2] Andrea. A s.i.r vector disease model with delay. Mathematical Modelling, 7:793–802, 1986. [3] N. H. Shah and J. Gupta. Seir model and simulation for vector borne diseases. Scientific Research Publisher, 2013. [4] G. Zaman, Y. H. Kang, and I. H. Jung. Stability analysis and optimal vaccination of sir epidemic model. BioSystems, 93(3):240–249, 2008. [5] M. Bohner and S. H. Streipert. The sis model on time scales. Pliska Studia Mathe- matica Bulgarica, 26:11–28, 2016. [6] Wayne P. London and James A. Yorke. Recurrent outbreaks of measles, chickenpox and mumps: I. seasonal variation in contact rates. American Journal of Epidemiology, 98(6):453–468, 1973. [7] D. Xiao and S. Ruan. Global analysis of an epidemic model with non-monotone incidence rate. Mathematical Biosciences, 208(2):419–429, 2007. [8] G. ur Rahman, R. P. Agarwal, and Q. Din. Mathematical analysis of giving up smoking model via harmonic mean type incidence rate. Applied Mathematics and Computation, 354:128–148, 2019. [9] A. Kaddar. Stability analysis in a delayed sir epidemic model with saturated incidence rate. Nonlinear Analysis: Modelling and Control, 15(3):299–306, 2010. [10] Y. Jin, W. Wang, and S. Xiao. An sirs model with a nonlinear incidence rate. Chaos, Solitons Fractals, 34(5):1482–1497, 2007. [11] A. Khan, R. Zarin, G. Hussain, N. A. Ahmad, M. Hafiz, and A. Yusuf. Stability analy- sis and optimal control of covid-19 with convex incidence rate in khyber pakhtunkhwa (pakistan). Results in Physics, 20:103703, 2021. [12] A. Khan, R. Zarin, G. Hussain, A. H. Usman, U. W. Humphries, and J. F. Gomez- Aguilar. Modelling and sensitivity analysis of hbv epidemic model with convex inci- dence rate. Results in Physics, 22:103836, 2021. [13] O. Babasola, O. Kayobe, O. J. Peter, F. C. Onwuegbuche, and F. A. Oguntolu. Time delayed modelling of covid-19 dynamics with a convex incidence rate. Informatics in Medicine Unlocked, 33:101124, 2022. [14] R. Zarin, A. Raouf, A. Khan, and U. W. Humphries. Modelling hepatitis b infection dynamics with a novel mathematical model incorporating convex incidence rate and real data. The European Physical Journal Plus, 138(11):1056, 2023. [15] B. Buonomo and D. Lacitignola. On the dynamics of an seir epidemic model with convex incidence rate. Ricerche di Matematica, 57:261–281, 2008. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 37 of 38 [16] Kamel Guedri et al. Rabies-related brain disorders: Transmission dynamics and epidemic management via educational campaigns and application of nanotechnology. The European Physical Journal Plus, 139(1):1–16, 2024. [17] V. Ambalarajan, A. R. Mallela, V. Sivakumar, P. B. Dhandapani, V. Leiva, C. Martin- Barreiro, and C. Castro. A six-compartment model for covid-19 with transmission dynamics and public health strategies. Scientific Reports, 14(1):22226, 2024. [18] T. Mathivanan, R. Bheeman, P. B. Dhandapani, I. Alraddadi, H. Ahmad, and T. Rad- wan. Monkeypox: A new mathematical model using the caputo-fabrizio fractional derivative. Fractals, page 2540128, 2025. [19] P. R. Murugadoss, V. Ambalarajan, V. Sivakumar, P. B. Dhandapani, and D. Baleanu. Analysis of dengue transmission dynamic model by stability and hopf bifurcation with two-time delays. Frontiers in Bioscience-Landmark, 28(6):117, 2023. [20] C. Castillo-Garsow, G. Jordan-Salivia, and A. Rodriguez-Herrera. Mathematical models for the dynamics of tobacco use, recovery and relapse. 1997. [21] G. Zaman. Qualitative behavior of giving up smoking models. Bulletin of the Malaysian Mathematical Sciences Society. Second Series, 34(2):403–415, 2011. [22] O. Khyar, J. Danane, and K. Allali. Mathematical analysis and optimal con- trol of giving up smoking model. International Journal of Differential Equations, 2021(1):8673020, 2021. [23] J. Singh, D. Kumar, M. A. Qurashi, and D. Baleanu. A fractional model for giving up smoking dynamics. Advances in Difference Equations, 2017:1–16, 2017. [24] O. Sharomi and A. B. Gumel. Curtailing smoking dynamics: A mathematical mod- eling approach. Applied Mathematics and Computation, 195(2):475–499, 2008. [25] A. Zeb, M. I. Chohan, and G. Zaman. The homotopic analysis method for approxi- mating of giving up smoking model in fractional order. Scientific Research Publishing, 2012. [26] D. P. Gao, N. J. Huang, S. M. Kang, and C. Zhang. Global stability analysis of an sveir epidemic model with general incidence rate. Boundary Value Problems, 2018(1):42, 2018. [27] M. Y. Li and J. S. Muldowney. A geometric approach to global-stability problems. SIAM Journal on Mathematical Analysis, 27(4):1070–1083, 1996. [28] D. L. Lukes. Differential Equations: Classical to Controlled, volume 162 of Mathe- matics in Science and Engineering. Academic Press, New York, NY, USA, 1982. [29] S. S. Alzaid and B. S. T. Alkahtani. Asymptotic analysis of a giving up smoking model with relapse and harmonic mean type incidence rate. Results in Physics, 28:104437, 2021. [30] Sarah Durkin, Emily Brennan, and Melanie Wakefield. Mass media campaigns to promote smoking cessation among adults: an integrative review. Tobacco Control, 21(2):127–138, 2012. [31] US Department of Health and Human Services. Preventing tobacco use among youth and young adults: a report of the surgeon general. Technical report, 2012. [32] Nur Ilmayasinta and Heri Purnawan. Optimal control in a mathematical model of smoking. Journal of Mathematical and Fundamental Sciences, 53(3), 2021. S. Bano et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6473 38 of 38 [33] Da-peng Gao and Nan-jing Huang. Optimal control analysis of a tuberculosis model. Applied Mathematical Modelling, 58:46–64, 2018. [34] Yinggao Zhou, Yiting Liang, and Jianhong Wu. An optimal strategy for hiv multi- therapy. Journal of Computational and Applied Mathematics, 263:326–337, 2014. [35] Haojie Yang, Yougang Wang, Soumen Kundu, Zhiqiang Song, and Zizhen Zhang. Dynamics of an sir epidemic model incorporating time delay and convex incidence rate. Results in Physics, 32:105025, 2022. Appendix A: Pseudocode for Numerical Scheme Implementing Optimal Control Strategy Algorithm 1 Forward-Backward Sweep Method for Optimal Control 1: Input: Initial state values S(0), E(0), I(0), R(0), . . . 2: Time step size ∆t 3: Total time T 4: Initial guess for controls u1(t), u2(t), . . . , uk(t) 5: Maximum number of iterations Nmax 6: Convergence tolerance ε 7: Output: Optimal state trajectories S(t), E(t), I(t), R(t), . . . 8: Optimal control functions u∗1(t), u ∗ 2(t), . . . , u ∗ k(t) 9: Step 1: Initialize time grid t ∈ [0, T ] with step size ∆t 10: Step 2: Initialize control functions u1(t), u2(t), . . . , uk(t) 11: Step 3: Repeat for n = 1 to Nmax a) Forward Sweep: – Integrate the state system using initial conditions – Use current control functions u1(t), . . . , uk(t) – Solve using a numerical method (e.g., 4th order Runge-Kutta) b) Backward Sweep: – Initialize adjoint variables λi(T ) = 0 – Integrate the adjoint system backward in time – Use state trajectories from the forward sweep c) Control Update: – Update controls u1(t), . . . , uk(t) using optimality conditions – Apply projection if control bounds exist (e.g., 0 ≤ u ≤ 1) d) Check Convergence: – If change in control functions < ε, then break loop 12: Step 4: Return final state and control trajectories Introduction Model formation Positivity and Boundedness Existence of equilibrium point Smoking-free equilibrium point Smoking endemic equilibrium point Reproductive number Stability Local stability Global stability Sensitivity Optimal control Optimal Control Existence Optimal control Characteristic of the PLSQ system Uniqueness of optimal control system Delay Stability analysis In the absence of delay Presence of delay Numerical simulation Discussion Conclusion