EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6451 ISSN 1307-5543 – ejpam.com Published by New York Business Global Exploring the Effectiveness of an Optimal Control Model for Marijuana Consumption Atta Ullah1,∗, Hamzah Sakidin1, Shehza Gul1, Thoraya N. Alharthi2, Diaa S. Metwally3, Yaman Hamed1, Kamal Shah4, Aceng Sambas5,6 1 Department of Applied Science, Universiti Teknologi PETRONAS, Seri Iskandar 32610, Perak, Malaysia 2 Department of Mathematics, College of Science, University of Bisha, P.O. Box 551, Bisha 61922, Saudi Arabia 3 Department of Accounting, Faculty of Business, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, 11432, Saudi Arabia 4 Department of Mathematics, University of Malakand, Chakdara Dir Lower, KPK, Pakistan 5 Faculty of Informatics and Computing, Universiti Sultan Zainal Abidin, Campus Besut, 22200 Terengganu, Malaysia 6 Department of Mechanical Engineering, Universitas Muhammadiyah Tasikmalaya, Tamansari Gobras 46196 Tasikmalaya, Indonesia Abstract. The primary aims of this research are to reduce marijuana use within the general pop- ulation, addressing its harmful effects and its status as an illegal substance that poses ongoing risks to public health in developing countries. To tackle this issue, we have developed a mathe- matical model that translates real-world marijuana consumption into mathematical terms using first-order nonlinear ordinary differential equations. A key objective of this study is to modify the “Non-smokers, Experimental-smokers, Recreational-smokers, Addicts, Hospitalized individu- als, and Prisoners” as referred to the NERAPH model, by incorporating optimal control measures. This modified model is designed to estimate the initial rate of marijuana transmission within these groups. To implement effective control strategies, we conduct sensitivity analyses to identify the most impactful parameters. These findings enable us to recommend targeted control variables for the most sensitive parameters, leading to the development of an optimal control model aimed at reducing marijuana consumption as efficiently as possible. Finally, we compare the outcomes of the basic and optimal control models through numerical results to evaluate the effectiveness of each approach. This comparison provides a clear measure of how optimal control techniques enhance the model’s ability to reduce marijuana use more effectively than the existing basic model. 2020 Mathematics Subject Classifications: 92D25, 92D30, 93B30, 93C95, 91A80 Key Words and Phrases: Basic mathematical model, basic reproduction number, sensitivity analysis, optimal control model, marijuana consumption ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6451 Email addresses: attaullahmastan@gmail.com (A. Ullah), hamzah.sakidin@utp.edu.my (H. Sakidin), adangul909@gmail.com (S. Gul), talhrthe@ub.edu.sa (T. N. Alharthi), dmetwally@imamu.edu.sa (D. S. Metwally), yaman.hamed@utp.edu.my (Y. Hamed), kamal@uom.edu.pk (K. Shah), acengs@umtas.ac.id (A. Sambas) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 2 of 36 1. Introduction Marijuana abuse has become one of the world’s biggest issues in the beginning of the modern era. Media outlets, both digital and conventional, regularly report on marijuana smoking. Populations, financial stability, and health may suffer because of the increase in marijuana misuse. This illness has a major negative impact on individuals, communities, and the country, especially the younger generation [1]. Marijuana use is a major public health concern [2]. According to recent studies, 18.7 million Americans admitted to con- suming marijuana, and 75% of those who did so report using it regularly [3]. Furthermore, 72 million Americans have smoked at least one marijuana cigarette. A major contributor to the worldwide epidemic of illnesses in the form of “disability- adjusted life years” (DALYs) is addiction to drugs disorders [4]. It has been discovered that consumption of drugs is linked to a higher risk of unintentional injuries, self-harm, infection with “human immunodeficiency virus” (HIV), “acquired immunodeficiency syndrome” (AIDS), and the liver cirrhosis, all which collectively add to the overall burden of disease. The medical cost connected with drug use disorders is further increased by indications of significant overlap amongst drug use problems along with other mental health conditions [5]. The use of marijuana is sometimes viewed as a “gateway” drug for more harmful substances or as a warning sign for other illegal substances. Adolescents who use gateway drugs are more prone to transition to more harmful substances [6]. David et al. noted the fact that there has been much discussion about how risky marijuana use is. Even while marijuana’s effects can seem harmless in comparison to those of other drugs, many people think that it frequently serves as the catalyst for someone to begin experimenting with narcotics [7]. Marijuana use is becoming more common as more states legalize marijuana for recre- ational as well as for medical purposes. According to national surveys, over 2 million indi- viduals in the United States with existing cardiovascular conditions have either used or are currently using marijuana in its various forms, such as inhalation and vaping. Cannabi- noid receptors are present in multiple tissues and cells, including platelets, adipose tissue, and myocytes. Data gathered from observations indicate potential links between mari- juana consumption and a wide range of adverse cardiovascular risks. It is worth noting that the potency of marijuana is increasing, and smoking marijuana shares similar car- diovascular health risks with tobacco smoking [8]. Marijuana can be used as medicine because it contains a complex blend of chemicals and cannabinoids. The positive impact ratio of marijuana has not been adequately characterized because there are more than 200 cannabinoids found in marijuana, some of which have been recognized and some of which are unknown. To extract cannabinoids, cannabis leaves are either vaporized or cooked. They are then smoked, typically in pipes or hand-rolled cigarettes (joints). However, some users combine it with food and brew it like tea [9]. High rates of marijuana use have a number of negative effects on a society, including an increase in crime and the resulting loss of lives and property, an increase in the number of murders and accidents, a significant loss of efficient man hours, a decrease in the financial gain of individuals due to marijuana A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 3 of 36 use, and ongoing, significant costs for substance abuse treatment and avoidance programs [10]. The short-term effects of marijuana use involve the rapid transfer of THC (the ac- tive compound) from the lungs to the bloodstream, which then distributes it to the brain and other organs. When marijuana is consumed orally, the absorption of THC occurs at a slower rate, leading to effects typically manifesting within 30 minutes to an hour. By overstimulating brain regions abundant in receptors, marijuana can have short-term impacts. The long-term effects of marijuana use also influence brain development. If in- dividuals start using marijuana during adolescence, it can hinder cognitive functions like thinking, memory, and learning and impact the establishment of crucial connections be- tween different brain regions involved in these functions. Researchers are still investigating the duration of marijuana’s effects and whether certain changes may be permanent [11]. Marijuana use leads to an elevated heart rate, which can persist for up to three hours after smoking. This increase in heart rate poses a potential risk for heart attacks. Individuals who are older or have pre-existing heart conditions may face a higher risk if they engage in marijuana use [12]. If a pregnant woman consumes marijuana, it can have detrimental effects on specific areas of the developing fetal brain. Infants who receive exposure to marijuana in utero have an increased risk of facing challenges related to memory, concentration, and ability to solve problems in comparison to those who didn’t get exposure. Furthermore, there exists an elevated probability of both neurological and behavioral issues in infants who received exposure to marijuana in utero [13]. Consistent consumption of marijuana may lead to the presence of THC in mother’s milk, which could potentially affect the brain’s growth of the breastfeeding infant. Additional investigation is imperative to attain a more profound understanding of the association between marijuana consumption and childbearing [14]. DeFilippis et al. provided a summary of the cardiovascular implications associated with marijuana use, including potential interactions with medications, and highlighted the need for further research to establish clearer guidelines regarding its safety for the cardiovascular system. They recommend that healthcare providers should actively screen and test for marijuana use, particularly when caring for young patients who present with cardiovascular conditions. This emphasis on screening aims to ensure appropriate management and treatment in clinical settings [8]. The study examined the timing of marijuana use in relation to various cardiovascular events, focusing on young patients who did not have any cardiovascular risk factors apart from recent marijuana consumption. The observed effects encompassed conditions such as atrial fibrillation, acute coronary syndromes, and instances of sudden death [15]. The use of marijuana can result in the emergence of a substance use disorder, which is a medical condition characterized by the inability to cease using the substance despite experiencing negative consequences in one’s health and personal life. Severe cases of substance use disorders are commonly referred to as addiction. Studies indicate that an estimated 9 and 30 percent of marijuana users may develop varying degrees of marijuana use disorder [16]. The impact of marijuana is facilitated by the endocannabinoid system, which involves the distribution of cannabinoid receptors in various tissues and cell types. While cannabi- A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 4 of 36 noid receptors are found in significant concentrations in the central and peripheral nervous systems, they also exist in platelets, adipose tissue, myocytes, the liver, the pancreas, and skeletal muscle. As a result, cannabinoids from external sources can affect multiple systems within the body [17, 18]. The growth and development of many children and adolescents are impacted by prenatal marijuana use. According to research, prenatal alcohol consump- tion is linked to deficiencies in the following areas: general mental growth and academic performance, cognition, visual skills, concentration, and capacity for problem-solving [19]. Both people with fetal alcohol syndrome (FAS) and people who were exposed to lower levels of drinking, have been found to have deficiencies in executive functioning, an um- brella term that includes the capacity to organize, concentrate, figure out solutions, and use focused on objectives behavior’s [20, 21]. A wide range of mental deficiencies, health risks, and medical issues are also linked to marijuana usage. Regular use of marijuana may be linked to small cognitive deficits in concentration as well as complicated executive skills, even if the use of marijuana for a long time does not cause any substantial deficits in mental performance [22]. For in- stance, compared to non-users, marijuana users over a decade showed a significant decline in neurological assessments of response time, speed, and consistency [23]. Additionally, compared to light users, heavy users showed considerably worse abnormalities in focus and executive performance [24]. Marijuana has an impact on young people everywhere. Young individuals have the potential to use marijuana, grow it, or transport it. Young individuals are drawn to marijuana use for a variety of reasons, including those related to their families, schools, and companions, economic circumstances, and the level of their physical surroundings. According to the majority of studies, the important risk phases for beginning to use marijuana are early teenage and late adolescence [25]. A rise in criminal activity and the consequent loss of life and assets, a rise in killings and incidents, a signifi- cant loss of productive individual hours, a decline in the earnings that people make caused by marijuana consumption, and the continuous, substantial expenses of addiction preven- tion and rehabilitation activities are just a few of the adverse impacts that substantial amounts of marijuana consumption have on the community [10]. The “American Association of Pediatrics” (AAP) is a particular academy that is con- scious of the paucity of evidence required to support recommendations made by profes- sionals. In 2018, the American Academy of Pediatrics declared that it is necessary to prevent the use of marijuana during breastfeeding due to the lacking information required for determining the effects of smoking marijuana on newborn behaviors. While smok- ing marijuana deserves to be regarded as a drug that interacts with the nursing field or whether preterm babies should be fed with emitted milk from moms who have acknowl- edged consuming marijuana recreationally were not addressed in the American Association of Pediatrics declaration [26]. Natural system has become increasingly interested in mathematical models due to their applications in several medical and health fields [27–30]. It has been demonstrated that mathematical models use ordinary differential equations can be useful in understanding the changing behaviour of natural phenomena. However, most biological processes exhibit memories or consequences in their behaviour, but traditional integer-order mathemati- A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 5 of 36 cal models neglect these effects. However, because fractional-order systems may reflect phenomena related to memory and influenced by genetic features, they are more suitable than integer-order ones in many disciplines [31–33]. Consequently, there has been a lot of interest lately in creating mathematical models of biology using “fractional differential equations (FDEs)” [34–37]. Compartmental models are regarded as one of the most often used kinds of epidemi- ological models for both observational and experimental research. Individuals in these models exhibit a wide range of characteristics while being in just a few of discrete stages; some of these conditions simply have declares that describe the people’s distinct char- acteristics, which can be epidemiologically significant in a dynamic system including an addictive agent or host. It has been recognized that using numerical techniques to find out “non-linear ordinary differential equations” with local as well as nonlocal variables is a useful mathematical tool [38]. 1.1. Preliminaries A mathematical framework for marijuana consumption in a precise population was proposed by Dauhoo et al. The stated population consists of smokers and non-smokers. Smokers are further divided into three different classes: recreational, experimental, and addicted. The model that includes marijuana use and population growth for consideration is called the “Non-users, Experimental-users, Recreational-users, and Addicted” (NERA) framework [39] is the subsequent system of mathematical equations as described in (1). Additionally, the graphical illustration of the earlier mathematical model is shown in Figure 1. Ṅ = β − (β + r1E + r1R)N + r3E + r5R+ r6A, Ė = −(β + r3 − r1N + r2R)E + r1NR, Ṙ = −(β + r4 + r5 − r2E)R, Ȧ = r4R− (β + r6)A. (1) A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 6 of 36 Figure 1: Graphical illustrations of the existing model [40, 41]. Ullah et al. [42], updated the Dauhoo et al. [39] model by incorporating two new significant classes known as addicts acquiring medical treatment (hospitalized class) and people taken into custody by the police (prisoners) to represent the said problem as shown in Equation (9). There are two primary categories that comprises the entire population: marijuana smokers and non-smokers. The expected outcome is the establishment of five subgroups of users, each of which represents a separate phase or level of substance abuse. However, the existing studies have a notable lack of research that integrates ”Optimal Control Technique” to develop and analyse intervention strategies in a mathematically rigorous way. Optimal control offers a powerful framework for identifying time-sensitive control measures for a faster recovery rate, yet its full potential remains underexplored in this context. The absence of “Optimal Control Technique” represents a significant gap. The study develops and examines a control-based system of nonlinear differential equations governed by “Pontryagin’s Maximum Principle”. Through this, necessary conditions are derived for optimality and construct time-dependent control functions aimed at minimizing marijuana use with a high recovery rate of marijuana users. Finally, the study conducts sensitivity analysis to evaluate the influence of key parameters, ensuring the robustness of both the model and the control strategy. By combining structural expansion with rigorous optimal control application, this research presents a comprehensive and practically relevant framework for understanding and mitigating marijuana consumption and its consequences across diverse populations, aiming to achieve significant recovery in the shortest possible time. The article proposed an optimal control model and is organized as follows. Section 2 presents the modified mathematical model and briefly introduces its associated compart- ments, excluding the involvement of optimal control. Fundamental qualitative properties, A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 7 of 36 including positivity, boundedness, the basic reproduction number, and sensitivity analysis, are also discussed in this section. The key part of this study which provides an optimal control model for marijuana consumption, is covered in Section 3. Section 4 illustrates and discusses the numerical results obtained from the proposed mathematical models. Fi- nally, Section 5 summarized the outcomes and significance of the proposed optimal control model and provides directions for future work. 2. Mathematical Model without Controls Model formulation is an important process where our understanding of a natural system is turned into a mathematical model. It involves two key steps: first we are constructing a mental picture and then converting it into a mathematical framework [42]. The preliminary phase involves identifying the essential elements, commonly referred to as state variables, as well as the fluctuations that designate the dynamics of the issue at hand. By employing the idea of preservation, the formulas of the theoretical frame- work delineate the frequency of fluctuation of the “state variables” as the composite of all arriving flows subtracted by the leaving fluxes in every segment. This could be suc- cessfully demonstrated through an illustration or workflow track, in which the features of the current situation are illustrated by barriers, interconnected by arrows that signify procedural relationships. In the subsequent phase, procedures will be clearly identified and articulated as mathematical equations [42]. According to [43], educating highly skilled individuals with advanced abilities and knowledge in mathematical modeling has become increasingly popular because of the expanding scope and utilization of mathematical representations in an abundance of fields. Formulating a hypothesis for a phenomenon that has been discovered and then testing the hypothesis by predicting the results of multiple tests under specific circumstances is a fundamental scholarly procedure. Typically, predictions and experimental results are compared. While differences indicate that the justification must be reformulated, consistency allows the argument to be recognized as a legitimate concept. To address differences with observational or experimental information, a framework that defines the major elements of the phenomenon, typically expressed mathematically, can be iteratively refined during the revision approach. An improved NERA model that focusses on comprehending the dynamics of a social pandemic will be presented in this part. The study population is categorised into two primary groups: individuals who smoke marijuana and those who do not. Five distinct stages are used to further categorise smokers, each of which represents a unique path towards addiction. The different classes within the population are characterized as follows: susceptible individuals (UN ), experimental users (UE), recreational smokers (UR), addicts (UA), hospitalized individuals (UH), and prisoner’s (UP ). The “entire population (T)” at a specific “time (t)” is represented by Equation (2). T (t) = UN (t) + UE(t) + UR(t) + UA(t) + UH(t) + UP (t). (2) A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 8 of 36 People who are vulnerable to the harmful effects of marijuana smoking are classified as being in the susceptible class (UN ). The mean lifespan of those who were engaged is ap- proximately 14 years old, and the rate at which new members enter class (UN ) is indicated by §. Due to mortality, some individuals naturally depart from this class. In addition, some people are inspired by recreational smokers, which leads them to begin smoking marijuana for fun. These people are then assigned to the experimental group. A subgroup of susceptible individuals begins engaging in recreational smoking because of their inter- actions with addicts, leading them to be categorized as recreational smokers. The size of this subgroup is influenced by the rate at which they interact with addicted individuals. A huge interaction rate classifies a larger influx of people into the experimental (UE) and recreational smokers (UR) classes. The interaction term between the susceptible (UN ) and recreational smokers (UR) classes is represented as ′m1URU ′ N , while the interaction term between the susceptible (UN ) and addicts (UA) classes is denoted as ′m4UAU ′ N . Hence, we can express the change in the susceptible population using the following differential Equation (3), which captures the dynamics of this population change. U̇N = § −m1URUN −m4UAUN − µUN +m3UE +m5UR + x4UH + x5(1− α3)UP . (3) The contact rate ′m′ 1 represents the frequency of interaction between susceptible indi- viduals and members of the recreational smoker’s class, while ′m′ 4 denotes the interaction ratio between susceptible individuals and addicts. On the other hand, ′m′ 3, ′m′ 5, ′x′4, and ′x′5 are the ratios that indicate the proportion of individuals from the experimental (UE), recreational (UR), hospitalized (UH), and prisoner’s (UP ) classes, respectively, who choose to join the susceptible class when they quit smoking due to various reasons. The natural death rate is represented as µUN . The term U̇N shows the change that occurs in the non-smoking class over time. Furthermore, the interaction term ′m1URU ′ N signifies the transition of individuals from the susceptible class (UN ) to the experimental class (UE). Thus, the time-dependent change in the experimental class is represented by the following differential Equation (4). U̇E = m1URUN − (m2 +m3 + µ)UE . (4) The rate of recruitment rate for individuals in the “non-users” class into the exper- imental category (UE) is ′m1URU ′ N . Experimental users transition into the recreational category at the ratio of m2UE . Some smokers of the experimental category stopped smok- ing because of the elders’ advice at a rate of ′m′ 3 and become vulnerable. Natural mortality affects some members of the experimental category at a rate of µUE . The symbol U̇E rep- resents the change over time that takes place in the experimental class. Therefore, the change in (UR) over time is represented by the differential Equation (5) that follows. U̇R = m4UAUN +m2UE − (x1 +m5 + µ)UR. (5) Individuals from the experimental class are recruited at a rate of m2UE , and those from the susceptible class at ′m4UAU ′ N respectively in the recreational smoker class. In A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 9 of 36 this class, the natural death rate is µUR. After completing the necessary time in the recreational class at the pace of x1UR, an individual’s transitions to the category of addicts. Due to their constrained environment, many people in this class develop susceptibility at the ratio of ′m′ 5. The subsequent system of “non-linear differential equation” as shown in (6) is formulated to identify the category of addicted users. U̇A = x1UR + x5α3UP − (x2 + x3 + e+ µ)UA. (6) After finishing the duration of their term in a recreational group, some people join the addicted class at a rate of x1UR, whereas after completing their jail sentence, some people rejoin the addicted compartment with a ratio of x5α3. Police pursue the gang members continuously and they are captured at a rate of ′x′3. Some patients get admitted in the hospital at a rate of ′x′2 for treatment. This class’s members pass away naturally at the pace of µUA, while through police involvement death rate is ′e′ for some of them. The subsequent differential Equation (7) determines when a member of the hospitalized class enters and exits. U̇H = x2UA − (x4 + µ)UH . (7) People are recruited at a rate of ′x′2 from the addicted class into the hospitalized class, and they rejoin the non-user class at a rate of ′x′4 once they have recovered. Some members of this class pass away suddenly at a rate of µUH . Like this, people enter the prisoner’s category from the addicted category with a ratio of ′x′3, and drop out of class at a rate of x5(1 − α3). The natural mortality rate of individuals in this class is denoted by µUP . The prisoner’s compartment is governed by the ordinary differential equation presented in Equation (8) as follows. U̇P = x3UA − (x5 + µ)UP . (8) The simplified version of the updated framework is seen in Figure 2. Moreover, utiliz- ing the provided schematic representation, the subsequent non-linear ordinary differential equations have been formulated, which encapsulate the comprehensive behavior associated with marijuana consumption. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 10 of 36 Figure 2: Graphical illustrations of the proposed modified model [42]. U̇N = § −m1URUN − (m4UA + µ)UN +m3UE +m5UR + x4UH + x5(1− α3)UP , U̇E = m1URUN −m3UE − (m2 + µ)UE , U̇R = m4UAUN +m2UE − x1UR − (m5 + µ)UR, U̇A = x1UR + x5α3UP − (x2 + x3 + e+ µ)UA, U̇H = x2UA − (x4 + µ)UH , U̇P = x3UA − (x5 + µ)UP . (9) A detailed description of the associated parameters is available in Table 1. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 11 of 36 Table 1: Description of parameters and their values Sym. Parameters Explanation Values Ref. § Birth rate of individuals 0.0015875 day−1 [44] m1 Impact of experimental smokers on non-smokers 0.446 day−1 [39] m2 Ratio of experimental users to recreational users 0.5 day−1 [39] m3 Proportion of experimental smokers who stopped smoking after guidance 0.17 day−1 [39] m4 Impact of addicts on non-smokers (susceptible) 0.001201 day−1 [40] m5 Ratio of recreational users who stop due to con- strained environment 0.002 day−1 [39] x1 Frequency of recreational users becoming ad- dicts after transition 0.025 day−1 [40] x2 Rate of hospital admission among heavy mari- juana users 0.22 day−1 [39] x3 Percentage of marijuana addicts facing incarcer- ation 0.0157871 day−1 [45] x4 Addicts’ probability of recovery after treatment 0.2010 day−1 [40] x5 Prisoner rehabilitation rate 0.0331 day−1 [45] α3 Rejoining probability of addicts after jail sen- tence 0.03 day−1 [40] e Police encountering addicts ratio 0.0005 day−1 [40] µ Natural death rate 0.006 day−1 [44] 2.1. Fundamental Qualitative Properties Some qualitative properties of the model (9) solutions are discussed here [46]. 2.1.1. Positivity and Boundedness It is known that the model (9) causes changes in the human population; thus, all param- eters of the given model tabulated in Table 1 are non-negative. Additionally, we need to ensure that every state variable used in the model is always non-negative [46]. Theorem 1. The state variables, UN (t), UE(t), UR(t), UA(t), UH(t), and UP (t), of the model (9) with initial conditions T (t) = T (0) are non-negative for all t > 0 [46]. Proof. Using the model (9) first equation, the result obtained: dUN dt ≥ −(m1UR +m4UA + µ)UN (t), (10) To achieve: d dt ( UN (t) exp ( µt+ ∫ t 0 M(x) dx )) ≥ 0, (11) A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 12 of 36 where M = m1UR +m4UA. The result of Equation (11) is: UN (t) ≥ UN (0) exp ( − ( µt+ ∫ t 0 M(x) dx )) > 0, ∀t > 0. (12) It can also be demonstrated that the other state variables UE(t), UR(t), UA(t), UH(t), and UP (t) are non-negative for all t > 0. Additionally, takes into consideration the feasible region defined by Γ ⊂ R6 +, where Γ = { (UN , UE , UR, UA, UH , UP ) ∈ R6 + : T ≤ § µ } . (13) It is demonstrated that Γ is positively invariant [46]. Theorem 2. The region Γ is positively invariant and attracting with respect to the proposed mathematical model (9) [46]. Proof. The total population’s (T ) rate of change is provided by dT dt = § − µT − eUA, (14) To achieve, dT dt ≤ § − µT, (15) The outcome of Equation (15) is according to the standard comparison theorem [47]. T (t) ≤ T (0)e−µt + § µ ( 1− e−µt ) . (16) T (t) ≤ § µ implies that the region Γ is positively invariant regarding the suggested model (9), as follows from T (0) ≤ § µ . Additionally, if T (0) ≤ § µ then either the outcome of the suggested model (9) reaches the region Γ in limited time, or T (t) converges to § µ while the users approach zero as t → ∞. Hence, the region Γ is attracts individuals. Declared alternatively, any solutions in R6 + eventually comes to Γ. It is sufficient to keep the behaviour of the model (9) in Γ to be considered, because of Theorems 1 and 2. The model can be seen as mathematically well-posed in this region [48]. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 13 of 36 2.2. Basic Reproduction Number The total number of secondary cases created by a single addicted individual in a fully non-smoking population is known as R0, fundamental reproduction number [49]. The “Next Generation Matrix Method” determines the initial rate of transmission [40, 50]. R0 = ρ(FV−1), (17) In the above setting (Equation 17), the “spectral radius” is denoted as ′ρ′, while the “Jacobian matrix” of the function ′f ′ is expressed as F = If . f = f1 f2 f3  = m1URUN m4UAUN 0  , (18) The column in Equation 18 symbolizes those who acquire addiction. F = F11 F12 F13 F21 F22 F23 F31 F32 F33  = 0 m1UN 0 0 0 m4UN 0 0 0  . (19) In similar way, the “Jacobian matrix” of the function v′ is expressed as V ′ = Jv′ ; where V ′ = v′1 v′2 v′3  =  −(m2 +m3 + µ)UE m2UE − (x1 +m5 + µ)UR x1UR − (x2 + x3 + e+ µ)UA  , (20) The individuals coming into or going out the category of addicts, apart from the people coming from the category of vulnerable, are identified in the column of the array V ′ as illustrated in Equation 20. V ′ = V ′ 11 V ′ 12 V ′ 13 V ′ 21 V ′ 22 V ′ 23 V ′ 31 V ′ 32 V ′ 33  , (21) V ′ = −(m2 +m3 + µ) 0 0 m2 −(x1 +m5 + µ) 0 0 x1 −(x2 + x3 + e+ µ)  . (22) The principal eigenvalue of the matrix FV ′−1, and consequently the basic reproduction number R0, is as follows: R0 = √ m1m2§ µ(m2 +m3 + µ)(x1 +m5 + µ) . (23) A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 14 of 36 2.2.1. Initial Transmission Rate within a Biological Context In Equation 23, ′m′ 1 signifies the influence ratio of experimental users on non-users (suscep- tible), while ′m′ 2 identifies the impact rate of recreational users on susceptible. Therefore, the term ′m1m2§′ of the initial transmission rate R0 indicates that specific individuals from the category of non-smokers are going to begin smoking marijuana and then become part of the smoker’s category driven by the effective value of marijuana users on susceptible. This implies that the term ′m1m2§′ indicates the transmission of marijuana from smokers to non-smokers individuals. The additional term (parameters) in R0 solely contribute to the magnitude of the initial transmission rate. 2.3. Sensitivity Analysis Assessing the sensitivity of parameter values is crucial for designing and controlling methods and for orienting future investigations. We can quantify the amount of change in a state variable whenever one of its parameters is changed using native sensitive indexes. To compute the sensitivity assessment, we follow Arriola’s technique [51]. The standardized advance sensitivity rating measures the proportionate variation in a variable in response to the proportionate variation in a parameter. When the variable depends smoothly on the parameter, this index can be alternatively described using partial derivatives. Definition. Variations in certain parameters influence the associated variables, and this proportional change is referred to as the level of sensitivity of the characteristics. The sensitiveness of the given function being used χ in relation to a specific parameter η is established, as indicating in references [51–53], under the condition that the function is differentiable with respect to that parameter [38, 54]. Υχ η = ∂χ ∂η · η χ (24) Table 2 provides the parameter sensitivity indexes, while Figure 3 illustrates the graph- ical results of the sensitive parameter indices. Table 2: Sensitivity indices of parameters Parameters Parameter values Sensitivity indices µ 0.006 −0.5954 x1 0.025 −0.3788 m1 0.446 +0.5000 m2 0.5 +0.1309 m3 0.17 −0.1258 m5 0.002 −0.0303 § 0.0015875 +0.5000 A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 15 of 36 Figure 3: Graphical illustration of sensitivity indices [42]. 3. Mathematical Model with Controls In this crucial section, we integrated optimal control mechanisms into the fundamental model (9) to improve it. Optimal control, one of the currently dynamic topics in control theory, has many real-world applications in economic growth, robotics, chemical engineer- ing, military aircraft, and many branches of physics [55–57]. The optimal control problem and its integration have evolved to perfection for linear time-invariant processes [58, 59]. Optimal control of nonlinear phenomena, which have been the subject of generations of intensive research, is considerably more difficult [60, 61]. Our main goal is to control the spread of addiction, including through controlling the level at which trial individuals affect non-users ′m′ 1. At +0.5000, we discovered that ′m′ 1 had the greatest sensitivity score. This suggests a clear correlation regarding the fundamental reproduction numbers: an increase of 20% in the number of reproductions correlates to a 10 percent increase in the impact rate of trial individuals over susceptible. The attraction rate ′x′1 across enjoyable and addictive users had the second-most sen- sitivity index, at -0.3788. This shows that the impact prevalence of addicts against recre- ational users is inversely correlated with the reproductive rate; for every 10% rise in the impact percentage of drug addicts, the basic reproductive number decreases by 3.7%. ′m′ 2 is the third greatest sensitive parameter, its sensitivity value is +0.1309, which is associated with the impact rate of addicts on people who are not users and the subsequent spread of susceptible to users for recreational purposes. This suggests that there is a 2.6 A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 16 of 36 percent rise in the number of reproductions for every 20 percent rise in the impact rate of addicted over nonsmoking compartment. ′m′ 3 is the fourth most sensitive parameter, the sensitivity index of this parameter is -0.1258, and it concerns the percentage of experimental smokers who become non- smokers after receiving instruction. This suggests that a 4.8 percent reduction in the total number of reproductions occurs for every 50% increase in the direction provided by religious academics and mentors to the trial users. The next critical parameter ′m′ 5, has a sensitivity index of -0.0303 and relates to the percentage of users who go back to non-smoking class during their leisure time because of their restrictive surroundings. Additionally, there exists an inverse relationship between this parameter’s sensitivity and the fundamental reproductive ratio. Remarkably, in this case, there is a direct proportionality between the number of reproductions and the level of interaction between trial users and people who do not smoke Additionally, a decrease in the human death rate is brought about by an increase in the recovering ratio of the population, which in turn lowers the initial rate of marijuana transmission. Here ′m′ 1, ′m′ 3 and ′x′1 are the crucial(targeted) parameters. Therefore, instead of focusing on every parameter, we only address the crucial parameters at a time: the rate at which trial users affect non-users, the rate at which trial smokers come back to non-smokers, and the rate at which recreational smokers quit smoking and comeback to non-smokers compartment. When these important factors are altered, other parameters will follow accordingly, either by increasing or decreasing. To develop the best control strategy, we identify the following control variables such as u1 represents the frequency of interactions between testing smokers and those who do not smoke; u2 represents the rate at which people return from the trial smokers section to the non-smokers section on the counsel of religious academics and older people; and u3 represents the rate at which leisure smokers transmit to the addicted smokers class because members of that class end up in hospitals and prisons, respectively. By introducing these control measures u1, u2 and u3, to get the formulated optimal control model is shown by Equation (25). U̇N = § −m1(1− u1)URUN − (m4UA + µ)UN + (m3 + u2)UE +m5UR + x4UH + x5(1− α3)UP , U̇E = m1(1− u1)URUN − (m3 + u2)UE − (m2 + µ)UE , U̇R = m4UAUN +m2UE − (x1 + u3)UR − (m5 + µ)UR, U̇A = x1UR + x5α3UP − (x2 + x3 + e+ µ)UA, U̇H = x2UA − (x4 + µ)UH , U̇P = x3UA − (x5 + µ)UP . (25) Reducing the influx of non-smokers into the trial marijuana consumer class is the aim of the best control approach. As a result, to accelerate the progress of rehabilitation for smokers who use it for experimental purposes at the first time and those who smoke for recreational purposes, respectively. In addition, minimizing the costs associated with A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 17 of 36 control implementation is accomplished by utilizing the most minimal realistic control variables u1, u2 and u3. Avoiding the costs associated with applying the controls u1, u2 and u3. Additionally, minimising the density of experimental and recreational smokers are the two principal goals of the objective functional J . This dual-use strategy aims to strike an appropriate equilibrium between the advantages of public health and economic factors, making sure that the intervention is both successful in reducing smoking rates and economical in terms of the use of resources [62]. J (u1, u2, u3) = ∫ T 0 ( a1UE + a2UR + a3UA + 1 2 ( d1u 2 1 + d2u 2 2 + d3u 2 3 )) dt. (26) The weight constants in this context are a1 for the significance of decreasing the per- centage of testing smokers, a2 for the impact percentage of recreational users over those who do not smoke and a3 for the attraction rate of smokers in the addicted class and recreational class. In addition, d1 represents the weight constant associated with the in- teraction between experimental smokers and non-smokers. d2 denotes the weight constant related to the interaction between recreational smokers and non-smokers, whereas d3 refers to the weight constant corresponding to the interaction between the addicted class and recreational smokers. The expressions 1 2(d1u 2 1), 1 2(d2u 2 2), and 1 2(d3u 2 3) represent the overall expense related to the execution of anti-marijuana smoking initiatives. This expense arises from efforts directed at establishing an environment that benefits non-smokers, assist- ing them in avoiding interactions with individuals categorized as experimental marijuana smokers, recreational marijuana smokers, and those considered addicted. The expenditure is directly related to the square of the respective control functions. Put differently, our implementation of non-linear control functions adheres to the methodologies employed in optimal control problems, as outlined in the academic literature [63, 64]. The main goal is to identify the optimal controls u∗1, u ∗ 2, and u∗3 that meet the specified criteria described below. J (u∗1, u ∗ 2, u ∗ 3) = min {J (u1, u2, u3) | (u1, u2, u3) ∈ D} , (27) where the collection of acceptable controls designated by D is D = u1, u2, u3 ∣∣∣∣∣∣∣ 0 ≤ u1min ≤ u1(t) ≤ u1max ≤ 1, 0 ≤ u2min ≤ u2(t) ≤ u2max ≤ 1, 0 ≤ u3min ≤ u3(t) ≤ u3max ≤ 1, t ∈ [0, T ]  . (28) 3.1. Solution of the Proposed Optimal Control Model To demonstrate the existence of the control measure, we examine the control structure (25) where t = 0 represents the starting position. A non-negative constrained equilibria exists for the current framework with “bounded Lebesgue measurable” entities and non- negative starting circumstances, as indicated by sources [65]. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 18 of 36 Initially the ‘Lagrangian’ and ‘Hamiltonian’ are established to get the optimal solution for the system. Moreover, ‘Lagrangian’ is provided for controlling the issue. Equation (25) gives L(t) = a1UE(t) + a2UR(t) + a3UA(t) + 1 2 ( d1u 2 1(t) + d2u 2 2(t) + d3u 2 3(t) ) , (29) Subsequently, we calculate the Hamiltonian H to the smallest value of the Lagrangian using the following method: H(t) = a1UE(t) + a2UR(t) + a3UA(t) + 1 2 ( d1u 2 1(t) + d2u 2 2(t) + d3u 2 3(t) ) + λ1 [ § −m1(1− u1)URUN − (m4UA + µ)UN + (m3 + u2)UE +m5UR + x4UH + x5(1− α3)UP ] + λ2 [ m1(1− u1)URUN − (m3 + u2 +m2 + µ)UE ] + λ3 [ m4UAUN +m2UE − (x1 + u3 +m5 + µ)UR ] + λ4 [ x1UR + x5α3UP − (x2 + x3 + e+ µ)UA ] + λ5 [ x2UA − (x4 + µ)UH ] + λ6 [ x3UA − (x5 + µ)UP ] . (30) The adjoint variables related to the state variables UN , UE , UR, UA, UH and UP were designated by the symbols λ1, λ2, λ3, λ4, λ5 and λ6 correspondingly. 3.2. The Existence and Characterization of the Optimal Control Problem We first demonstrate that the framework (25) has a solution, and then we will demon- strate that an optimal control exists [66]. 3.2.1. Existence of an optimal control Theorem 3. The set of an optimal control (u∗1, u ∗ 2, u ∗ 3) minimize J over D, subject to the initial conditions specified at t = 0. Proof. The target function described previously Equation (26) is convex for the control setting (u∗1, u ∗ 2, u ∗ 3), as both the state and control variables are non-negative [67]. Here D = {(u1, u2, u3)}|ui(t) is “Lebesgue measurable” on [0, tf ], 0 ≤ ui(t) ≤ 1, for i = 1, 2, 3. The set D is obviously closed and convex. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 19 of 36 Additionally, there are bounds to the ideal optimal control scheme. As a result, the set of ideal control is minimal. Conversely, the target functional has quadratic values for u1, u2 and u3. Consequently, integrand in the objective functional a1UE(t) + a2UR(t) + a3UA(t) + 1 2 ( d1u 2 1(t) + d2u 2 2(t) + d3u 2 3(t) ) , on control set D is evidently convex. So, we can select a single value Q > 1 along with the positive values s1 > 0, and s2 > 0 such that J (u1, u2, u3) ≥ s1 ( |u1|2 + |u2|2 + |u3|2 )Q/2 − s2 as the state variables are limited. Thus, the existence of the optimal control is proved □. We continue through the system’s control settings in the next outcome, ’Pontryagin’s Maximum Principle’ (PMP) [44]; Assume that T = U∗ N , U∗ E , U ∗ R, U ∗ A, U ∗ H , U∗ P be the column-wise vector of the current state variables associated with the controls (u∗1, u ∗ 2, u ∗ 3). And let (P, u) represent the con- trol system solution, for P ∈ M and u = (u∗1, u ∗ 2, u ∗ 3) then there exist a vector function represented by the symbols λ = (λ1, λ2, . . . , λn) satisfies the following requirements: dP dt = ∂H(P, t, u, λ) ∂λ , 0 = ∂H(P, t, u, λ) ∂u , ∂H(P, t, u, λ) ∂P = λ′(t). (31) By the previously mentioned derivation, which it follows. u∗ =  0, if ∂H ∂u < 0, 0 ≤ u∗ ≤ 1, if ∂H ∂u = 0, 1, if ∂H ∂u > 0. (32) The necessary limitations on the Hamiltonian H of the system (30) are applied. 3.2.2. Characterization of the optimal control Theorem 4. Let U∗ N (t), U∗ E(t), U ∗ R(t), U ∗ A(t), U ∗ H(t), U∗ P (t) be the best possible responses with corresponding optimal controls (u∗1(t), u ∗ 2(t), u ∗ 3 (t) for the ideal controls (25) and (26). Consequently, there exist adjoint or costate variables λi, for i=1,2,3,4,5,6 that satisfy. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 20 of 36 dλ1 dt = (λ1 − λ2)m1(1− u1)U ∗ R + (λ1 − λ3)m4U ∗ A + λ1µ, dλ2 dt = (λ2 − λ1)(m3 + u2) + (λ2 − λ3)m2 + λ2µ− a1, dλ3 dt = (λ1 − λ2)m1(1− u1)U ∗ N + (λ3 − λ4)x1 + (λ3 − λ1)m5 + λ3u3 + λ3µ− a2, dλ4 dt = (λ1 − λ3)m4U ∗ N + (λ4 − λ5)x2 + (λ4 − λ6)x3 + λ4(e+ µ)− a3, dλ5 dt = (λ5 − λ1)x4 + λ5µ, dλ6 dt = (λ6 − λ1)x5 + (λ1 − λ4)x5α3 + λ6µ. (33) Using transversality (or boundary) conditions [68, 69]. λ1(tf ) = 0, λ2(tf ) = 0, λ3(tf ) = 0, λ4(tf ) = 0, λ5(tf ) = 0, and λ6(tf ) = 0. (34) Additionally, for t ∈ [0, tf ] the control setting u∗1, u ∗ 2 and u∗3 are provided by u∗1 = max { min { (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 } , 0 } . (35) u∗2 = max { min { (λ2 − λ1)U ∗ E d2 , 1 } , 0 } . (36) u∗3 = max { min { λ3U ∗ R d3 , 1 } , 0 } . (37) Proof. We differentiate the HamiltonianH with respect to the state variables UN , UE , UR, UA, UH and UP to find the costate equations and the transversality conditions (33). More- over, to get (35)-(37), we differentiate the Hamiltonian H with respect to control measures u1, u2, and u3. Here is how the Hamiltonian H is described: H(t) = a1UE(t)+a2UR(t)+a3UA(t)+ 1 2 ( d1u 2 1(t) + d2u 2 2(t) + d3u 2 3(t) ) + 6∑ i=1 λi(t)Pi(UN , UE , UR, UA, UH , UP ), where P ∈ M. In addition, PMP can be used to get the adjoint variables and transversality criteria for t ∈ [0, tf ] [70–73], that are suitable. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 21 of 36 λ̇1 = − ∂H ∂UN = (λ1 − λ2)m1(1− u1)U ∗ R + (λ1 − λ3)m4U ∗ A + λ1µ, λ̇2 = − ∂H ∂UE = (λ2 − λ1)(m3 + u2) + (λ2 − λ3)m2 + λ2µ− a1, λ̇3 = − ∂H ∂UR = (λ1 − λ2)m1(1− u1)U ∗ N + (λ3 − λ4)x1 + (λ3 − λ1)m5 + λ3u3 + λ3µ− a2, λ̇4 = − ∂H ∂UA = (λ1 − λ3)m4U ∗ N + (λ4 − λ5)x2 + (λ4 − λ6)x3 + λ4(e+ µ)− a3, λ̇5 = − ∂H ∂UH = (λ5 − λ1)x4 + λ5µ, λ̇6 = − ∂H ∂UP = (λ6 − λ1)x5 + (λ1 − λ4)x5α3 + λ6µ. (38) According to the optimality requirements, for t ∈ [0, tf ] we have the following. ∂H ∂u1 = d1u ∗ 1 + λ1m1U ∗ RU ∗ N − λ2m1U ∗ RU ∗ N = 0 at u∗i (t) for i = 1, 2, 3. ⇒ u∗1(t) = (λ2 − λ1)m1U ∗ RU ∗ N d1 , In the same direction, u∗2(t) = (λ2 − λ1)U ∗ E d2 , u∗3(t) = λ3U ∗ R d3 . By employing the controlling space’s assets, we obtained u∗1(t) =  0, if (λ2 − λ1)m1U ∗ RU ∗ N d1 ≤ 0, (λ2 − λ1)m1U ∗ RU ∗ N d1 , if 0 < (λ2 − λ1)m1U ∗ RU ∗ N d1 < 1, 1, if (λ2 − λ1)m1U ∗ RU ∗ N d1 ≥ 1. (39) Follow the same way, to get u∗2(t) =  0, if (λ2 − λ1)U ∗ E d2 ≤ 0, (λ2 − λ1)U ∗ E d2 , if 0 < (λ2 − λ1)U ∗ E d2 < 1, 1, if (λ2 − λ1)U ∗ E d2 ≥ 1. (40) A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 22 of 36 And u∗3(t) =  0, if λ3U ∗ R d3 ≤ 0, λ3U ∗ R d3 , if 0 < λ3U ∗ R d3 < 1, 1, if λ3U ∗ R d3 ≥ 1. (41) In compact form, this can be rewritten as follows: u∗1 = max { min ( (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 ) , 0 } , u∗2 = max { min ( (λ2 − λ1)U ∗ E d2 , 1 ) , 0 } , u∗3 = max { min ( λ3U ∗ R d3 , 1 ) , 0 } . The characterization of the most efficient controls (u∗1, u ∗ 2, u ∗ 3) is denoted by the calcula- tion (35) to (37). The ideal control and state are determined by addressing the optimality arrangement, which consist of the costate structure (33) and (34), the current state frame- work (25) with boundary criteria, and the classification of the ideal control (35) to (37). We employ the very first and transversality requirements in conjunction with the formu- lation of the optimal control (u∗1(t), u ∗ 2(t), u ∗ 3(t)) provided by (35)-(37) to determine the optimality system. Furthermore, with respect to ideal control measures (u1, u2, u3) the Lagrangian’s second derivative is positive. So, the control is optimum(maximum) at the set (u∗1, u ∗ 2, u ∗ 3). By putting the outcomes of u∗1, u ∗ 2, and u∗3 in the model (25) we obtained the subsequent system (42). A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 23 of 36 U̇∗ N = § −m1 ( 1−max { min ( (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 ) , 0 }) U∗ RU ∗ N − (m4U ∗ A + µ)U∗ N + ( m3 +max { min ( (λ2 − λ1)U ∗ E d2 , 1 ) , 0 }) U∗ E +m5U ∗ R + x4U ∗ H + x5(1− α3)U ∗ P , U̇∗ E = m1 ( 1−max { min ( (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 ) , 0 }) U∗ RU ∗ N − ( m3 +max { min ( (λ2 − λ1)U ∗ E d2 , 1 ) , 0 }) U∗ E − (m2 + µ)U∗ E , U̇∗ R = m4U ∗ AU ∗ N +m2U ∗ E − ( x1 +max { min ( λ3U ∗ R d3 , 1 ) , 0 }) U∗ R − (m5 + µ)U∗ R, U̇∗ A = x1U ∗ R + x5α3U ∗ P − (x2 + x3 + e+ µ)U∗ A, U̇∗ H = x2U ∗ A − (x4 + µ)U∗ H , U̇∗ P = x3U ∗ A − (x5 + µ)U∗ P . (42) where the Hamiltonian at H∗ (t, U∗ N , U∗ E , U ∗ R, U ∗ A, U ∗ H , U∗ P , λ1, λ2, λ3, λ4, λ5, λ6, u ∗ 1, u ∗ 2, u ∗ 3) is : A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 24 of 36 H∗(t) = a1U ∗ E(t) + a2U ∗ R(t) + a3U ∗ A(t) + 1 2 [ d1 ( max { min ( (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 ) , 0 })2 + d2 ( max { min ( (λ2 − λ1)U ∗ E d2 , 1 ) , 0 })2 + d3 ( max { min ( λ3U ∗ R d3 , 1 ) , 0 })2 ] + λ1 [ § −m1 ( 1−max { min ( (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 ) , 0 }) U∗ RU ∗ N − (m4U ∗ A + µ)U∗ N + ( m3 +max { min ( (λ2 − λ1)U ∗ E d2 , 1 ) , 0 }) U∗ E +m5U ∗ R + x4U ∗ H + x5(1− α3)U ∗ P ] + λ2 [ m1 ( 1−max { min ( (λ2 − λ1)m1U ∗ RU ∗ N d1 , 1 ) , 0 }) U∗ RU ∗ N − ( m3 +max { min ( (λ2 − λ1)U ∗ E d2 , 1 ) , 0 }) U∗ E − (m2 + µ)U∗ E ] + λ3 [ m4U ∗ AU ∗ N +m2U ∗ E − ( x1 +max { min ( λ3U ∗ R d3 , 1 ) , 0 }) U∗ R − (m5 + µ)U∗ R ] + λ4 [x1U ∗ R + x5α3U ∗ P − (x2 + x3 + e+ µ)U∗ A] + λ5 [x2U ∗ A − (x4 + µ)U∗ H ] + λ6 [x3U ∗ A − (x5 + µ)U∗ P ] (43) We solve the systems (42) and (43) analytically to determine the state structure and the best control. Using an iterative process, we take consideration of the parameter values shown in Table 1 to get the numerical findings regarding the optimal state model. 4. Numerical Results and Discussion This section evaluates the generated mathematical framework’s numerical outcomes with and without control mechanisms in place. The purpose of the simulation is to find out how well control factors work as preventive measures against the general population’s transmission of marijuana use. We employ the ‘fourth-order Runge-Kutta’ (RK4) tech- nique for this. The RK4 method is a well-liked computational technique for examining “ordinary differential equations” (ODEs), and it has been recognized for its many ad- vantages. When it comes to the numerical modeling of a collection of ten “first-order ODEs” with defined transversality conditions, the RK4 method is recognized as the best approach. Its efficacy and reliability across a range of fields are acknowledged, as evidenced by comparative findings [74]. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 25 of 36 Employing suggested control measures across the simulation time frame, we first an- alyze problem (25) towards timescale using the “fourth-order Runge-Kutta procedure.” Next, we exploit the state formulations in the present cycle to handle structure (33) by the application of scenario (34) and the backward approach. Utilizing a convex combina- tion of the control measures from the earlier cycle and the corresponding quantities from systems (35) to (37), we modify our control measures. Until the convergence occurs, this repeated practice is carried out. Within the domain of interest [0, 600], the display time T is specified in days. The values of the given parameters used for the numerical results are presented in Table 1. Additionally, the initial values for the state variables (UN (t0), UE(t0), UR(t0), UA(t0), UH(t0), UP (t0)) = (1000, 100, 100, 100, 50, 50) are chosen to reflect a hypothetical but realistic population structure where non-users are in the majority, and user-related are proportionally smaller. These values are in line with initial population assumptions used in similar studies [50, 75–77] and serve as a baseline for understanding control effects. While the weight constants in the objective functional a1 = 0.181, a2 = 0.010, a3 = 0.014, d1 = 0.081, d2 = 0.010 and d3 = 0.014 are selected based on sensitivity tuning to ensure an effective balance between the intensity of control efforts and their influence on reducing marijuana prevalence. Slightly higher weights are assigned to compartments associated with more sever outcomes such as addicted, hospital- ized individuals and prisoners, reflecting their greater importance in shaping interventions’ priorities. The overall distribution of weights aligns with commonly adopted practices in related optimal control models [50, 75]. When interpreting the graphs, note that individuals without control measures are indi- cated by solid lines, while those with control measures are indicated by dash-dotted lines. In Figure 4, we plot the experimental individuals into two systems, (9) and (25). The solid line represents the population of experimental individuals in system, (9) without control measures, while the dash-dotted line represents the population of experimental individuals in system (25) with control measures. By day 18, the density of experimental smokers can be reduced with the recovery of 100 individuals if the control variables are maintained; otherwise, it will take a longer time to reduce the number of experimental smokers to zero. Similarly, the number of recreational smokers can be reduced within 22 days with the total of 122 recovered people if control measures are implemented; otherwise, it will take much longer to achieve this reduction, as described in Figure 5. Figure 6 describes the popu- lation of addicted individuals under both optimality and non-optimality systems. From this figure, we clearly see that the number of addicted smokers can be reduced within 27 days with the recovery of 100 individuals if control measures are implemented. Without these measures, it would not be possible to achieve the same reduction within the same time-frame. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 26 of 36 Figure 4: The graph compares the control periods for the experimental class UE , under conditions with and without control measures. Figure 5: The graph compares the control periods for the recreational class UR, under conditions with and without control measures. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 27 of 36 Figure 6: The graph compares the control periods for the addicted class UA, under conditions with and without control measures. Figure 7: The graph compares the control periods for the hospitalized class UH , under conditions with and without control measures. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 28 of 36 Figure 8: The graph compares the control periods for the prisoners class UP , under conditions with and without control measures. Figure 9: The plot represents six adjoint variables: λ1, λ2, λ3, λ4, λ5 and λ6. To determined the adjoint equations using a “fourth order backward Runge-Kutta method” due to the transversality conditions. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 29 of 36 Figure 10: The figure illustrates an intervention aimed at reducing the influence of experimental smokers on non-smokers by creating a restricted environment. Figure 11: The figure depicts an intervention involving the creation of a supportive environment and guidance from elders or religious scholars, aimed at helping experimental smokers return to the non-smokers’ group. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 30 of 36 Figure 12: The figure illustrates an intervention that accelerates the transition of recreational smokers into the addicted smokers’ group, as individuals in this category are more likely to end up in hospitals or prisons. In the case of using optimal controls the individual admitted to the hospital for treat- ment can be stabilized within 42 days, with a recovery of 50 individuals. However, under a non-optimal system, it takes a much longer time to achieve the same reduction, as shown in Figure 7. Similarly, for the prison population, if control measures are implemented, the number of individuals can be reduced within 140 days with a recovery of 50 individu- als. Without these control measures, it will take a much longer time to achieve the same output, as seen in Figure 8. The plot in Figure 9 illustrates six adjoint variables λ1, λ2, λ3, λ4, λ5 and λ6 within the optimality system. We solve these adjoint equations using a “fourth order backward Runge-Kutta method” due to the transversality conditions. Figures 10-12 show the effects of the control variables u1, u2, and u3 respectively. u1 represents the frequency of interac- tions between experimental smokers and non-smokers; u2 denotes the rate at which trial smokers revert to non-smokers compartment on the advice of religious leaders and elders; and u3 reflects the rate at which recreational users move to addiction, often resulting in hospitalizations and incarceration. The implementation of optimal control measures significantly accelerates the reduction of smoking populations across various categories, including experimental, recreational, and addicted smokers. The analysis demonstrates that with appropriate interventions, smoking can be effectively controlled within a period of 4 to 5 months. The comparison between optimal and non-optimal systems reveals substantial differences in recovery times, emphasizing the importance of targeted strategies in managing smoking-related health A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 31 of 36 outcomes. These findings underscore the critical role of well-designed interventions in reducing the societal impact of smoking. 5. Conclusions This study presents a mathematical analysis of general population trends in mari- juana consumption, emphasizing optimal time-dependent intervention techniques. The mathematical model comprises a network of six-dimensional first-order non-linear ordi- nary differential equations, with three key control variables (u1, u2, u3) incorporated for effective recovery rate and time management. Specifically, u1 represents the frequency of interactions between experimental smokers and non-smokers; u2 denotes the rate at which trial smokers revert to non-smokers compartment on the advice of religious lead- ers and elders; and u3 reflects the rate at which recreational users progress to addiction, often leading to hospitalizations and incarceration. While the model effectively demon- strates a significant reduction in marijuana consumption and enhanced recovery through optimal control strategies, one notable limitation is the exclusion of an economic cost anal- ysis. The study focuses primarily on minimizing marijuana usage and time to recovery, without incorporating the total cost associated with implementing control interventions a factor that is crucial for real-world policy planning. The most effective control concept has been employed to determine optimal strategy for addressing marijuana consumption. For numerical simulation, the “fourth order forward-backward Runge-Kutta method” is employed to simulate both controlled and uncontrolled scenarios based on “Pontryagin’s maximum principle”. Computational findings demonstrate the effectiveness of the pro- posed control mechanism in reducing overall marijuana use when applied concurrently within the same region. The numerical outcomes of the proposed optimal control tech- niques confirm that the suggested strategy effectively mitigates and prevents the spread of marijuana use, achieving a higher recovery rate within a shorter time frame as compared to both the updated basic model and previous studies. Future research may focus on extend- ing the model by incorporating a fractional-order framework to better capture memory effects and complex behavioural dynamics. Additionally, integrating economic cost func- tions into the optimal control framework could enable the design of more practical and cost-efficient strategies for combating marijuana consumption in the general population. Acknowledgements The authors thank the reviewers for their constructive and valuable comments, which led to the improvement of the entire manuscript. References [1] DP Leong, KK Teo, S Rangarajan, P Lopez-Jaramillo, A Avezum Jr, A Orlandini, et al. World population prospects 2019. department of economic and social affairs population dynamics. new york (ny): United nations; 2019 (https://population. un. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 32 of 36 org/wpp/download/, accessed 20 september 2020). the decade of healthy ageing. geneva: World health organization. World, 73(7):362, 2018. [2] United States. Substance Abuse and Mental Health Services Administration. Office of Applied Studies. Results from the 2005 National Survey on Drug Use and Health: National Findings. Department of Health and Human Services, Substance Abuse and Mental Health . . . , 2006. [3] H. J. Steadman, M. W. Deane, J. P. Morrissey, M. L. Westcott, S. Salasin, and S. Shapiro. A SAMHSA research initiative assessing the effectiveness of jail diversion programs for mentally ill persons. Psychiatric Services, 50(12):1620–1623, 1999. [4] L. Degenhardt et al. The global burden of disease attributable to alcohol and drug use in 195 countries and territories, 1990–2016: a systematic analysis for the global burden of disease study 2016. The Lancet Psychiatry, 5(12):987–1012, 2018. [5] C. Hjorthøj et al. Association between alcohol and substance use disorders and all- cause and cause-specific mortality in schizophrenia, bipolar disorder, and unipolar depression: a nationwide, prospective, register-based study. The Lancet Psychiatry, 2(9):801–808, 2015. [6] A. L. Bretteville-Jensen, H. O. Melberg, and A. M. Jones. Sequential patterns of drug use initiation–can we believe in the gateway theory? The BE Journal of Economic Analysis Policy, 8(2), 2008. [7] D. M. Fergusson and L. J. Horwood. Does cannabis use encourage other forms of illicit drug use? Addiction, 95(4):505–520, 2000. [8] E. M. DeFilippis et al. Marijuana use in patients with cardiovascular disease: Jacc review topic of the week. Journal of the American College of Cardiology, 75(3):320– 332, 2020. [9] T. T. Yusuf. Modelling marijuana smoking epidemics among adults: An optimal control panacea. Journal of Modelling Simulation, Identification, and Control, 2:83– 97, 2014. [10] C. Rossi. The role of dynamic modelling in drug abuse epidemiology. Bulletin on Narcotics, 54(1-2):33–44, 2002. [11] M. H. Meier et al. Persistent cannabis users show neuropsychological decline from childhood to midlife. Proceedings of the National Academy of Sciences, 109(40):E2657–E2664, 2012. [12] K. C. Young-Wolff et al. Trends in self-reported and biochemically tested marijuana use among pregnant females in california from 2009-2016. JAMA, 318(24):2490–2491, 2017. [13] G. A. Richardson, C. Ryan, J. Willford, N. L. Day, and L. Goldschmidt. Prenatal alcohol and marijuana exposure: effects on neuropsychological outcomes at 10 years. Neurotoxicology and Teratology, 24(3):309–320, 2002. [14] J. A. Galli, R. Andari Sawaya, and F. K. Friedenberg. Cannabinoid hyperemesis syndrome. Current Drug Abuse Reviews, 4(4):241–249, 2011. [15] S. Rezkalla and R. A. Kloner. Cardiovascular effects of marijuana. Trends in Cardio- vascular Medicine, 29(7):403–407, 2019. [16] D. S. Hasin et al. Prevalence of marijuana use disorders in the united states between A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 33 of 36 2001-2002 and 2012-2013. JAMA Psychiatry, 72(12):1235–1242, 2015. [17] K. Mackie. Cannabinoid receptors: where they are and what they do. Journal of Neuroendocrinology, 20:10–14, 2008. [18] S. Zou and U. Kumar. Cannabinoid receptors and the endocannabinoid system: sig- naling and function in the central nervous system. International Journal of Molecular Sciences, 19(3):833, 2018. [19] S. N. Mattson, A. M. Goodman, C. Caine, D. C. Delis, and E. P. Riley. Executive functioning in children with heavy prenatal alcohol exposure. Alcoholism: Clinical and Experimental Research, 23(11):1808–1815, 1999. [20] M. D. Lezak, D. B. Howieson, D. W. Loring, and J. S. Fischer. Neuropsychological Assessment. Oxford University Press, USA, 2004. [21] D. Tranel. Development of the concept of ’executive function’ and its relationship to the frontal lobes. In Handbook of Neuropsychology. 1994. [22] W. Hall and N. Solowij. Adverse effects of cannabis. The Lancet, 352(9140):1611– 1616, 1998. [23] S. S. Mendhiratta, V. K. Varma, R. Dang, A. K. Malhotra, K. Das, and R. Nehra. Cannabis and cognitive functions: a re-evaluation study. British Journal of Addiction, 83(7):749–753, 1988. [24] H. G. Pope and D. Yurgelun-Todd. The residual cognitive effects of heavy marijuana use in college students. JAMA, 275(7):521–527, 1996. [25] A. Strashny. Age of substance use initiation among treatment admissions aged 18 to 30. Technical report, Substance Abuse and Mental Health Services Administration, 2016. [26] S. A. Ryan et al. Marijuana use during pregnancy and breastfeeding: implications for neonatal and childhood outcomes. Pediatrics, 142(3), 2018. [27] C. J. Silva et al. Optimal control of the covid-19 pandemic: controlled sanitary deconfinement in portugal. Scientific Reports, 11(1):3451, 2021. [28] I. Area, F. Ndairou, J. J. Nieto, C. J. Silva, and D. F. Torres. Ebola model and optimal control with vaccination constraints. arXiv preprint, 2017. [29] H. M. Srivastava, I. C. Area Carracedo, and J. Nieto. Power-series solution of com- partmental epidemiological models. Mathematical Biosciences and Engineering, 2021. [30] B. Li and Z. Eskandari. Dynamical analysis of a discrete-time sir epidemic model. Journal of the Franklin Institute, 360(12):7989–8007, 2023. [31] D. Baleanu, S. S. Sajjadi, A. Jajarmi, and Ö. Defterli. On a nonlinear dynamical system with both chaotic and nonchaotic behaviors: a new fractional analysis and control. Advances in Difference Equations, 2021(1):234, 2021. [32] F. Ndäırou, I. Area, J. J. Nieto, C. J. Silva, and D. F. Torres. Fractional model of covid-19 applied to galicia, spain and portugal. Chaos, Solitons & Fractals, 144:110652, 2021. [33] A. Ullah, A. Ahmad, U. Khan, and A. Abdussamad. A time-fractional model for brinkman-type nanofluid with variable heat and mass transfer. City University In- ternational Journal of Computational Analysis, 5(1):11–30, 2022. [34] D. Baleanu, S. S. Sajjadi, J. H. Asad, A. Jajarmi, and E. Estiri. Hyperchaotic behav- A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 34 of 36 iors, optimal control, and synchronization of a nonautonomous cardiac conduction system. Advances in Difference Equations, 2021:1–24, 2021. [35] V. Erturk, E. Godwe, D. Baleanu, P. Kumar, J. Asad, and A. Jajarmi. Novel fractional-order lagrangian to describe motion of beam on nanowire. Preprint, 2021. [36] D. Baleanu, S. Zibaei, M. Namjoo, and A. Jajarmi. A nonstandard finite difference scheme for the modeling and nonidentical synchronization of a novel fractional chaotic system. Advances in Difference Equations, 2021(1):308, 2021. [37] A. Jajarmi and F. A. Ghassabzade. Optimal control and general fractional description for a complex biological system. Progress in Fractional Differentiation and Applica- tions, 9(3), 2023. [38] M. Ur Rahman, G. Alhawael, and Y. Karaca. Multicompartmental analysis of middle eastern respiratory syndrome coronavirus model under fractional operator with next- generation matrix methods. Fractals, 31(10):2340093, 2023. [39] M. Z. Dauhoo, B. S. N. Korimboccus, and S. B. Issack. On the dynamics of illicit drug consumption in a given population. IMA Journal of Applied Mathematics, 78(3):432– 448, 2013. [40] A. Ullah et al. Sensitivity analysis-based control strategies of a mathematical model for reducing marijuana smoking. AIMS Bioengineering, 10(4):491–510, 2023. [41] A. Ullah, H. Sakidin, S. Gul, K. Shah, Y. Hamed, and T. Abdeljawad. Mathematical model with sensitivity analysis and control strategies for marijuana consumption. Partial Differential Equations in Applied Mathematics, 10:100657, 2024. [42] A. Ullah, H. Sakidin, K. Shah, Y. Hamed, and T. Abdeljawad. A mathematical model with control strategies for marijuana smoking prevention. Electronic Research Archive, 32(4):2342–2362, 2024. [43] S. M. Moghadas and M. Jaberi-Douraki. Mathematical Modelling: A Graduate Text- book. John Wiley & Sons, 2018. [44] M. Zamir, G. Zaman, and A. S. Alshomrani. Sensitivity analysis and optimal control of anthroponotic cutaneous leishmania. PLOS ONE, 11(8):e0160513, 2016. [45] A. Ullah et al. Validation of the neraph model for cannabis consumption through sen- sitivity analysis. Semarak International Journal of Fundamental and Applied Math- ematics, 5(1):44–60, 2025. [46] J. O. Akanni, S. Olaniyi, and F. O. Akinpelu. Global asymptotic dynamics of a nonlinear illicit drug use system. Journal of Applied Mathematics and Computing, 66:39–60, 2021. [47] H. L. Smith and P. Waltman. The Theory of the Chemostat: Dynamics of Microbial Competition. Cambridge University Press, 1995. [48] H. W. Hethcote. The mathematics of infectious diseases. SIAM Review, 42(4):599– 653, 2000. [49] A. Ullah. Comprehensive validation: Enhancing the modified nerah model via rig- orous sensitivity analysis. Journal of Advanced Research in Applied Sciences and Engineering Technology, 64(4):58–73, 2025. [50] M. Zamir, T. Abdeljawad, F. Nadeem, A. Wahid, and A. Yousef. An optimal control analysis of a covid-19 model. Alexandria Engineering Journal, 60(3):2875–2884, 2021. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 35 of 36 [51] L. Arriola and J. M. Hyman. Sensitivity analysis for uncertainty quantification in mathematical models. In Gerardo Chowell, James M. Hyman, Lúıs M.A. Betten- court, and Carlos Castillo-Chavez, editors, Mathematical and Statistical Estimation Approaches in Epidemiology, pages 195–247. Springer, 2009. [52] X. Zhu, P. Xia, Q. He, Z. Ni, and L. Ni. Ensemble classifier design based on pertur- bation binary salp swarm algorithm for classification. Computer Modeling in Engi- neering & Sciences, 135(1):653–671, 2023. [53] Y. Ramzan, B. M. Fadhl, S. Niazai, A. U. Awan, and K. Guedri. Decoding the transmission and subsequent disability risks of rabineurodeficiency syndrome without recuperation. Scientific Reports, 15(1):17322, 2025. [54] B. Fatima, M. Yavuz, M. U. Rahman, and F. S. Al-Duais. Modeling the epidemic trend of middle eastern respiratory syndrome coronavirus with optimal control. Math- ematical Biosciences and Engineering, 20(7):11847–11874, 2023. [55] L. Tang, L. Zhao, and J. Guo. Research on pricing policies for seasonal goods based on optimal control theory. ICIC Express Letters, 3(4), 2009. [56] V. Manousiouthakis and D. J. Chmielewski. On constrained infinite-time nonlinear optimal control. Chemical Engineering Science, 57(1):105–114, 2002. [57] T. Notsu, M. Konishi, and J. Imai. Optimal water cooling control for plate rolling. In Second International Conference on Innovative Computing, Information and Control (ICICIC 2007), pages 313–313. IEEE, 2007. [58] A. E. Bryson. Applied Linear Optimal Control Hardback with CD-ROM: Examples and Algorithms. Cambridge University Press, 2002. [59] B. Li, T. Zhang, and C. Zhang. Investigation of financial bubble mathematical model under fractal-fractional caputo derivative. Fractals, 31(05):2350050, 2023. [60] A. Jajarmi, N. Pariz, A. V. Kamyad, and S. Effati. A novel modal series representation approach to solve a class of nonlinear optimal control problems. Min. J, 1:2, 2011. [61] R. Padhi and M. Kothari. Model predictive static programming: a computationally efficient technique for suboptimal control design. International Journal of Innovative Computing, Information and Control, 5(2):399–411, 2009. [62] S. Olaniyi, J. Akanni, and O. Adepoju. Optimal control and cost-effectiveness analysis of an illicit drug use population dynamics. Journal of Applied Nonlinear Dynamics, 12(01):133–146, 2023. [63] S. Abimbade, S. Olaniyi, O. Ajala, and M. Ibrahim. Optimal control analysis of a tuberculosis model with exogenous re-infection and incomplete treatment. Optimal Control Applications and Methods, 41(6):2349–2368, 2020. [64] K. O. Okosun. On the dynamics malaria-dysentery co-infection model. Journal of Biological Systems, 28(02):453–474, 2020. [65] N. Chitnis, J. M. Hyman, and J. M. Cushing. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bulletin of Mathematical Biology, 70:1272–1296, 2008. [66] A. Kouidere, L. E. Youssoufi, H. Ferjouchia, O. Balatif, and M. Rachik. Optimal control of mathematical modeling of the spread of the covid-19 pandemic with high- lighting the negative impact of quarantine on diabetics people with cost-effectiveness. A. Ullah et al. / Eur. J. Pure Appl. Math, 18 (3) (2025), 6451 36 of 36 Chaos, Solitons & Fractals, 145:110777, 2021. [67] D. L. Lukes. Differential Equations: Classical to Controlled. Academic Press, 1982. [68] D. Baleanu, A. Jajarmi, S. S. Sajjadi, and D. Mozyrska. A new fractional model and optimal control of a tumor-immune surveillance with non-singular derivative operator. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(8), 2019. [69] A. Jajarmi and M. Hajipour. An efficient finite difference method for the time- delay optimal control problems with time-varying delay. Asian Journal of Control, 19(2):554–563, 2017. [70] O. Balatif, B. Khajji, and M. Rachik. Mathematical modeling, analysis, and optimal control of abstinence behavior of registration on the electoral lists. Discrete Dynamics in Nature and Society, 2020(1):9738934, 2020. [71] L. S. Pontryagin. Mathematical Theory of Optimal Processes. Routledge, 2018. [72] D. Kada, A. Kouidere, O. Balatif, M. Rachik, and E. H. Labriji. Mathematical modeling of the spread of covid-19 among different age groups in morocco: Optimal control approach for intervention strategies. Chaos, Solitons & Fractals, 141:110437, 2020. [73] A. Kouidere, D. Kada, O. Balatif, M. Rachik, and M. Naim. Optimal control approach of a mathematical modeling with multiple delays of the negative impact of delays in applying preventive precautions against the spread of the covid-19 pandemic with a case study of brazil and cost-effectiveness. Chaos, Solitons & Fractals, 142:110438, 2021. [74] N. Ahmad, S. Charan, and V. P. Singh. Study of numerical accuracy of runge- kutta second, third and fourth order method. International Journal of Computer and Mathematical Sciences, 4:111, 2015. [75] N. Ali, M. I. Chohan, and G. Zaman. Optimal control of a time delayed hiv-1 infection model. European Journal of Pure and Applied Mathematics, 12(2):506–518, 2019. [76] N. Ali, G. Zaman, and A. S. Alshomrani. Optimal control strategy of hiv-1 epidemic model for recombinant virus. Cogent Mathematics, 4(1):1293468, 2017. [77] A. A. Lashari, S. Aly, K. Hattaf, G. Zaman, I. H. Jung, and X.-Z. Li. Presentation of malaria epidemics using multiple optimal controls. Journal of Applied Mathematics, 2012(1):946504, 2012.