Synthetic drugs transmission: stability analysis and optimal control

The synthetic drugs that are becoming increasingly popular in drugmarkets are some of the most destructive drugs which have made headlines causing serious social and health care issues in the past few years. In this work, a synthetic drug transmission model with general contact rate and Holling Type-II functional responses among susceptible and drug addicts (both psychological and physiological) is proposed. Sensitivity analysis provides that controlling the contact rate among the drug addicts and the susceptible is better than the treatment. Further, an optimal control problem has been formulated tominimize the cost and drug addiction by choosing the counselling treatment (including awareness programmes) as a control variable. Numerical analysis has indicated that if the control policy is implemented, then it will be economically viable for long term. In the proposed model system, it is observed how counselling can prevent the psychologically addicted persons from taking more drugs. ARTICLE HISTORY Received 17 October 2018 Accepted 21 May 2019


Introduction
Infectious diseases often have serious consequences on population and in fact, it can affect their overall social well-being, health and other related developments. Moreover, the prevention and control of disease transmission are important as it is related to high morbidity and mortality. In recent years, infectious diseases such as influenza, flu, etc. have posed major challenges across the World. These diseases not only increase the disease burden but also add an economic burden on the society which affects the development plans of many countries (The Economic and Social; Russell, 2004;World Health Organization, 2006). In this context, we would like to mention that approximately 1.4 million deaths are caused by various types of hepatitides (Hepatitis A-E) only each year (World Hepatitis Day, 2014). Like these infectious diseases, drug addiction is another threat to the society. In fact, deaths caused by only over-consumption of drugs have increased in these years. Some drugs are natural, means they are derived from natural plants without using chemicals such as opium poppies (heroin, morphine, codeine), coca leaves (cocaine), psilocybin mushrooms (shrooms), marijuana, etc. But except natural drugs, there exist synthetic drugs which are created using chemicals (not natural ingredients). New synthetics appear in the market constantly, therefore, it is not possible to give a compile list. Chemists who choose to evade arrest can simply shift the formula slightly and come up with something new that might not be listed in the text of laws that ban drugs. For illustration, K2 also known as Spice, Ecstasy (Molly) and bath salts are all types of synthetic drugs. Their strength, composition and ingredients are unknown to the consumers as synthetic drugs are created in illegal laboratories to bypass regulations prohibiting controlled substances. Attractive names and colourful appearances can sometimes hide their bad effects. Synthetic drugs are extremely dangerous and can cause addiction, severe health issues and even death. Synthetic drugs appeared first time in the United States around 2009 but their popularity is still increasing for the last nine years (Tricia Escobedo, 2013). They are especially popular among teenagers for their high level of accessibility (source). According to a study of the Center for Substance Abuse Research (CESAR) conducted in 2013 at the University of Maryland at College Park, about 12% of high school students are found to take synthetic drugs regularly (Tomo Drug Testing, 2017). Most of them are not aware of the dangerous effects of these substances.
Synthetic drugs are known as 'new drugs' or 'club drugs' as they mainly appear in entertainment venues. These drugs directly affect the central nervous system (CNS) as a new type of chemical synthetic mental drugs. Methamphetamine, MDMA , triazolam, etc. are the most common synthetic drugs. Synthetic drugs affect our psychology making more addictive towards those rather than the traditional drugs. These drugs have their own symptoms, i.e. when someone takes these drugs, he quickly appears psychotic symptoms such as excitement, mania, depression, hallucination and so on. Then gradually his behaviour becomes out of control resulting in a series of social problems including serious violent crimes. Moreover, many infectious diseases such as hepatitis and AIDS can easily infect drug-users (Kapp, 2008). Statistical data from South African police showed that between 2002 and 2006, only drug-related crimes increased in a huge range, in fact from 621 to more than 3,000 in Cape Town. From the 2005 antenatal survey, HIV infections increased almost 2.6% reaching in 15.7% from 2003 to 2005 in the Western Cape Province of South Africa (Kalula & Nyabadza, 2011).
Medical and healthcare facilities have improved significantly in the past few years but for the complex nature of diseases, the policymakers find a bit difficult to decide the proper detection and prevention. Different tools and techniques are designed to predict the system dynamics and also to suggest suitable control intervention. Mathematical modelling is one of the most useful tools for this purpose (Brauer & Castillo-Chavez, 2001;Gaff & Schaefer, 2009;Lenhart & Workman, 2007). Several research works have been done based on the pharmaceutical as well as non-pharmaceutical control interventions such as vaccination, quarantine, treatment, etc. Gaff and Schaefer (2009), Lenhart and Workman (2007), Joshi, Lenhart, Li, and Wang (2006).
A survey has given a graph that shows us the number of synthetic drug-users is continuously increasing. For instance, from 2014 drug users in China has increased by 7.7% and by the end of 2015, there are almost 2.354 million drug users in China which include 1.34 million synthetic drug abusers which accounted for 57.1% (www.nncc626.com/2016-02/18/c 128731173 2.html.). International practice reveals that the ratio of explicit drug addicts to implicit ones is 1:5. Although, in every country, registered addicts number are far below than actual drug addicts. Also, young people are increasingly falling prey which increases synthetic drugs transmission. Take China for example, from China's drug situations report in 2015 (www.nncc626.com/2016-02/18/c 128731173 2.html.), 62.4% of the existing drug addicts is below 35 years and the main reason of this much addiction is curiousness and excitement. Someone may think that they will not become an addict after taking once but it is nothing but a misconception as literature reveals only using once can lead to addiction as long as enough. As per the information of China National Narcotics Control Commission, there is hardly any mature treatment of drugs addiction around the world. Medications and psychotherapy are the major treatment measures for drug addicts. At the time of admitting to the drug institutions, the addict is first going through a process of detoxification which helps them to recover cerebral neurons and improve sleep and followed by psychotherapy which makes them mentally strong enough to restrain themselves from sucking drugs repeatedly. But for the seized addicts, these methods are not so helpful. By the records of the survey done in 2015 the people whose health condition is deteriorated by synthetic drugs are almost 54.6% of total seized addicts in China (www.nncc626.com/2016-02/18/c 128731173 2.html.). So it is very important to take rational measures to monitor the present rate of spreading of synthetic drugs.
By forming epidemic models, one mainly tries to simulate and reveals the nature of epidemics and provides theoretical rules and results for preventing and controlling diseases (www.nncc626.com/2016-02/18/c 128731173 2.html., www.nncc626.com.) (Cen, Feng, & Zhao, 2014;Fang, Li, Martcheva, & Cai, 2014;Feng, Cen, Zhao, & Velasco-Hernandez, 2015;Huo, Chen, & Wang, 2016;Mulone & Straughan, 2009;Mushanyu, Nyabadza, & Stewart, 2015;Nyabadza, Njagarah, & Smith, 2013;Saha & Samanta, 2019a, 2019bSamanta, 2011;Samanta, Sen, & Maiti, 2016;White & Comiskey, 2007). White and Comiskey (2007) first studied a simple heroin model where the authors divided the total population into susceptible, heroin addicts who were not in treatment and heroin addicts who were under treatment. Mulone and Straughan (2009) studied the qualitative behaviours on White-Comiskey model. Samanta and Sharma (2012) modified heroin epidemic model by dividing the class of drug users without treatment into two classes: light drug users and hard drug users. Samanta (2011) analysed a nonautonomous heroin epidemic model with distributed time delay. Fang et al. (2014) set up a heroin model with two distributed delays and with the help of Lyapunov function, global asymptotic stabilities of steady states have been proved. But in all these models, only traditional drugs have been considered. But compared to the models related to the traditional drug, fewer research works have been done for synthetic drugs (Ma, Liu, Xiang, & Li, 2018;Mushanyu et al., 2015;Nyabadza et al., 2013). Analysis of methamphetamine transmission of Western Cape province of South Africa can be found in the work of Mushanyu et al. (2015). The accidental drug-users may feel guilty after taking drug first time because they may take it by mistake or bewitched by friends. So many of them try to avoid the repetition of contacting with synthetic drugs consciously once again and in general, they cannot be addictive. So, characterizing the effects of synthetic drugs transmission (caused by psychology and physiology) is necessary. Taking into account these factors, we attempt to analyse the transmission of synthetic drugs in this work and also try to control the drug epidemic by implementing some effective measures.
Diseases exhibit a lot of economic burdens which includes productivity loss, health carerelated expenses, losses due to disease-related death and loss of employment, etc. Also, implementation costs in policy choice include the cost for treatment, vaccination, etc. So, policymakers have to implement such a policy that can control the spread of disease as well as can minimize the overall cost incurred during a certain time period. Based on the availability of resources, a single control intervention or multiple control interventions can be implemented. Several research works have already been done in epidemiology using various control interventions including both pharmaceutical and non-pharmaceutical interventions (Behncke, 2000;Castilho, 2006;Gaff, Schaefer, & Lenhart, 2011;Joshi, Lenhart, Hota, & Agusto, 2015;Kassa & Ouhinou, 2015;Lenhart & Workman, 2007;Zeiler, Caulkins, Grass, & Tragler, 2010). For example, Behncke (2000) have incorporated pharmaceutical treatment and health-promotional campaigns as control measures and analysed their effect together with the importance of financial support for an SIR model. Qualitative analysis indicates that control policies reduced the disease level while financial support promoted the campaigns which helped to suppress the disease transmission rate.
In the present work, we have analysed how the awareness campaigns and regular counselling can possibly reduce the addiction towards drugs in a certain population. We have introduced a compartmental model with two stages of addiction (psychological and physiological) and assumed that psychologically addicted individuals can go for regular counselling to restrain themselves from taking drugs. The proposed model takes into consideration that the individuals who are in treatment can move back to the physiologically addicted phase if they stop counselling therapy and pharmaceutical treatments suddenly. Overview of the rest of the article is as follows: in the following section, we have formulated a synthetic drugs model with psychological and physiological addicts. In Section 3, the well-posed behaviour of the basic mathematical model has been shown with the positivity and boundedness of the solution. In Section 4, the basic reproduction number (R 0 ) of the underlying model system is derived using next generation matrix method. In Section 5, we perform equilibrium analysis and their existence criterion depending on R 0 . Our study includes the stability analysis of the equilibrium points of the system which has been presented in Section 6. Section 7 contains the sensitivity analysis of the system which helps to find the effectiveness of the parameters on reproduction number. Numerical simulations for the proposed model system (without control) have been performed in Section 8 to support the related analytical results. In the next part, we have formulated the corresponding optimal control problem in Section 9. In the subsequent section, we have performed numerical simulations for the proposed model system (with control) to validate the analytical findings and in the last section, we end our paper with a brief discussion.

Mathematical model: basic equations
Synthetic drugs or club drugs are chemically created to impersonate another drug such as marijuana, cocaine or morphine and are more harmful than the traditional drugs. These drugs can affect on our CNS causing strong psychological dependence and can spread mainly among youths in a high rate (Dangerous Synthetic Drugs, 2013). Here we have divided the total human population (N(t)) into four classes at time t such as susceptible (S(t)), psychologically addicts (P 1 (t)), physiologically addicts (P 2 (t)) and addicts under treatment (T(t)). Also we have incorporated Michaelis-Menten (Holling type-II) functional response to describe the transmission among (S, P 1 ) and (S, P 2 ). The susceptible primarily become psychological addicts while first coming in contact with a drug addict, but after getting habituated with taking drugs, the individual is likely to become a physiological addict. A psychological or physiological addict will enter into the treatment compartment at the time of taking treatment and rehabilitation. We shall use S(t) = S, P 1 (t) = P 1 , P 2 (t) = P 2 and T(t) = T for sake of calculations. The corresponding model dynamics is presented by the following system of non-linear ODE's: All the model parameters are assumed to be positive constants with following interpretation: : This work focuses on drug users mainly in the age group 15−64 years (Kelly, Carvalho, & Teljeur, 2003). Considering this fact, the rate of entering this age group every year is taken as the recruitment rate (of susceptible) assumed to be a constant. β 1 , β 2 : Probability of transmission from susceptible to physiological and psychological drug addicts respectively. A: Average number of contacts with others per unit time. α: Proportion of psychological addicts who become physiological drug addicts by taking drugs in a regular basis, i.e. escalation rate from psychological to physiological addicts. η: Rate at which some drug-users in treatment may escape and reenter the physiological addict state, i.e. relapse rate. γ , ρ : Per capita pharmaceutical treatment rates for psychological and physiological addicts respectively. δ: Natural death rate.
Generally, a susceptible one can easily be caught into the web of drug when the individual comes in contact with a physiological addict as well as a psychological addict. For the sake of calculations, it is assumed that the probability of the susceptible to become a drug-user after coming in contact with a physiological addict is a multiple of the probability of the susceptible to become a drug-user after coming in contact with a psychological addict. So, we are taking β 2 (N) = β (say) and β 1 (N) = bβ into the underlying model as the effective contact rate where b is a positive constant.
Hence our model comes into the following form: (1) From system (1), by adding all equations, we get So, the total population is N = δ , which is a constant. Let us scale the state variables with respect to the total population: Then system (1) becomes: where a = A N , with initial conditions:

Non-negativity and boundedness
Let us discuss on positivity and boundedness of system (2) to ensure that the model is well-posed or well behaved.

Non-negativity of solutions
Theorem 3.1: All solutions of system (2) that start in R 4 + remain non-negative forever.
Proof: As R.H.S of (2) is completely continuous and locally Lipschitzian on C, (s, p 1 , p 2 , r) of system (2) with the initial conditions (3) exists and is unique on [0, ς ], where 0 < ς < +∞ (Hale, 1977). From the first equation of system (2), we get From the second equation of system (2), we get From the third equation of system (2), we get From the fourth equation of system (2), we get Theorem 3.2: All solutions of system (2) that start in R 4 + are uniformly bounded.

Proof:
The solution N(t) of the above differential equation has the following property: where

Basic reproduction number
Basic reproduction number (Brauer & Castillo-Chavez, 2001) is defined as the number of psychologically addicted individuals who have started to take drugs for the first time after getting touch with a single drug addict individual during his or her effective addiction period when introduced into susceptible population. In the following we will find the basic reproduction number R 0 of system (2) with the help of next generation matrix method formulation (van den Driessche & Watmough, 2002). Let x = (p 1 , p 2 , r, s) and a 0 = α + δ + γ , a 1 = ρ + δ, a 2 = η + δ, then system (2) can be written as Now, FV −1 is the next generation matrix of system (2) and basic reproduction number of the system will be the spectral radius of matrix FV −1 denoted and defined as (van den Driessche & Watmough, 2002): where a 0 = α + δ + γ , a 1 = ρ + δ and a 2 = η + δ.

Local stability of E 0
From (5), Jacobian matrix corresponding to E 0 is So, λ 1 = −δ(< 0) and other three eigenvalues are the roots of the equation: By Routh-Hurwitz criterion E 0 is locally asymptotically stable when Q 1 , Q 3 > 0 and Q 1 Q 2 − Q 3 > 0. Hence we lead to the theorem:

Local stability of E *
Now, the Jacobian matrix corresponding to E * (s * , p * 1 , p * 2 , r * ) is given by The characteristic equation is given by By the Routh-Hurwitz criterion (Kot, 2001), it follows that all eigenvalues of the characteristic equation have negative real part if and only if (i) B i > 0 for i = 1, 2, 3, 4; (ii) So, we arrive to the following result: Theorem 6.2: The synthetic drug addiction equilibrium E * of system (2) is LAS if Routh-Hurwitz criterion is fulfilled which is consistent with R 0 > 1 1+a .

Global stability of E 0
Theorem 6.3: The synthetic drug-free equilibrium E 0 of system (2) is globally asymptotically stable (GAS) if R 0 < 1 1+a .

Proof:
To check the global stability of DFE, let us consider the following positive definite function about E 0 : Here L(p 1 , p 2 , r) is a positive definite function for all (s, p 1 , p 2 , r) other than E 0 . The time derivative of L computed along the solutions of system (2) is as follows: Furthermore dL dt = 0 if and only if p 1 = 0. Therefore, the largest compact invariant set in {(s, p 1 , p 2 , r) ∈ : dL dt = 0}, when R 0 < 1 1+a , is the singleton E 0 . Hence by LaSallea's invariance principle (LaSalle, 1976), E 0 is globally asymptotically stable in when R 0 < 1 1+a .

Global stability of E *
Theorem 6.4: If a(a + s * ) < ab < (a + s * ) 2 holds, then the synthetic drug addiction equilibrium E * (s * , p * 1 , p * 2 , r * ) of system (2) is globally asymptotically stable in the region where Proof: Let us consider the positive solution of the system (2) about E * (s * , p * 1 , p * 2 , r * ) and construct the following Lyapunov function V as: Differentiating V(t) with respect to t along the solution of system (2), we get By the steady-state of the equilibrium point, we have Also, dV dt | E * = 0. Then by Lyapunov LaSallea's theorem (LaSalle, 1976), E * is globally asymptotically stable in the interior of for R 0 > 1 1+a .

Sensitivity analysis
The basic reproduction number (R 0 ) of system (2) depends on seven parameters, namely, per capita contact rate (β), average number of contacts with others per unit time (a), escalation rate from psychological to physiological addicts (α), relapse rate (η), per capita treatment rates for psychological and physiological addicts, respectively (γ , ρ) and natural death rate (δ). Among these parameters, we cannot control the parameters: a, α, η, δ. Therefore, to examine the sensitivity of R 0 to the parameters β, ρ and γ , normalized forward sensitivity index with respect to each of these parameters are computed as follows (Arriola & Hyman, 2005): So, it is clear that the basic reproduction number (R 0 ) is most sensitive to changes in β, probability of transmission from susceptible to drug addicts (both physiological and psychological). If β increases R 0 increases in same proportion and if β decreases R 0 also decreases in same proportion. On the other hand, ρ and γ have an inversely proportional relationship with R 0 , i.e. the size of the increase in any of them causes a decrease in R 0 and a decrease in any of them causes an increase in R 0 . But the increase in ρ, i.e. the treatment rate for physiological addicts can not help as much as the treatment rate for psychological addicts γ can do. So, it is better to focus either on β, the contact rate and γ , treatment rate for psychological addicts. As R 0 is more sensitive to changes in β than γ , so it is sensible to focus on β, per capita contact rate to control the drug abuse in population.

Numerical simulations
In this section, we present computer simulations of different solutions of system (2) using MATLAB. We have taken the total population as N = 54 million in South Africa. Based on the range provided in Table 1, we have made two tables (Table 2 and Table 3) for our numerical simulation and verification. Considering β = 0.0008 with the parametric values of Table 2, we have got R 0 = 0.0019(< 1 1+a = 0.17), while for β = 0.8 we have got R 0 = 1.92(> 0.17). Moreover, β = 0.8 along with the parameter values from Table 2 provides drug addiction equilibrium E * = (0.47, 0.01, 0.27, 0.24). It is observed that for R 0 < 0.17 the susceptible population (s) only persists and psychologically addicted (p 1 ), physiologically addicted (p 2 ) population and those who are under treatment (r) are going to extinct, i.e. the population approaches to the drug-free equilibrium or disease-free equilibrium E 0 (1, 0, 0, 0). We have taken mainly six initial points in the drawing of the stability diagram in three-dimension of E * as:  Table 2.
In the section of sensitivity analysis, we have already showed analytically that the contact rate is more effective to control the drug abuse rather than providing treatment to psychologically and physiologically addicts. The graphs in Figure 8 is consistent with our analytical findings. It is observed that with increasing value of β, R 0 increases monotonically. On the other hand, increasing value of ρ and γ provide decaying graph of R 0 . For the diagram in Figure 8, we have taken β = 0.0008 and the data from Table 2 for (a) and from Table 3 for (b).

Optimal control problem
In this section, we have formulated an optimal control problem corresponding to system (1) considering the effect of 'awareness campaigns, counselling and other preventions' as control policy. We have taken system (1) instead of system (2) to get a better observation when the control policy is implemented. It is aimed to observe how a suitable choice of the control variable can prevent the individuals from taking drugs and also to optimize the cost incurred in the implementation of the control policy. First, we have to describe the control policy and then to determine the cost corresponding to the policy.
Enhancing the response of psychological individuals via awareness programs and regular counselling: Promoting the awareness programs and counselling in a regular interval can induce the behavioural change among psychological individuals. Awareness campaigns not only prevent the population from taking drugs but also make them aware about the consequences. Now the availability of resources related to medical diagnosis, campaigns, treatments, etc. cannot be unlimited, in view of this a saturated treatment rate function uP 1 1+ξ P 1 has been incorporated in system (1). Here represents the treatment rate (via counselling) together with the impact of awareness campaigns and u is the intensity of treatment with saturation constant as ξ −1 . There are different costs involved like diagnosis, medicines and other related costs when counselling is provided. So, u can be used as a possible tool to trigger the responsiveness of psychological individuals with 0 ≤ u ≤ 1. Here 0 represents no improvement through counselling period, whereas 1 is indicating full improvement. Hence the control intensity u fully depends on the effort of the psychological individuals to stop themselves from taking drugs.
In the following, it is aimed to determine the optimal treatment via counselling with minimum cost by implementing the control. From the previous discussions, we have got that the acceptable set for the control variable u(t) as follows: Here T f represents the final time upto which the control policy can be executed. In this case, u(t) is a measurable and bounded function.

Determination of total cost
Let us first determine the total cost which has to be minimized for control intervention in the proposed model system.
Cost involved in counselling and awareness programs: The total cost which is associated with disease burden and counselling policy for psychological patients is given by The cost associated with psychological patients for losing man power and corresponding wealth is represented by w 1 P 1 (t) (Gaff & Schaefer, 2009;Joshi et al., 2006;Kassa & Ouhinou, 2015). It also indicates the opportunity loss, i.e. the lose of productivity due to addiction. The term w 2 u 2 (t) includes both the cost for awareness campaigns and the cost regarding the counselling policy such as medication charges, diagnosis charges, etc. Hence we incorporate a quadratic non-linear term u 2 (t) to present the cost (Gaff & Schaefer, 2009;Joshi et al., 2006;Kassa & Ouhinou, 2015).
The following control problem is considered based on the above-mentioned discussions along with the cost functional J 1 that has to be minimized: subject to the following model: with the initial conditions S 0 > 0, P 1,0 ≥ 0, P 2,0 ≥ 0 and T 0 ≥ 0. Here the functional J 1 represents the total cost that means the sum of the costs as stated. The integrand is which denotes the current value of cost at time t. Parameters w 1 and w 2 are positive weights (assumed as constants) which balance the units of integrand also Gaff Schaefer (2009); Kassa Ouhinou (2015). Let us denote u(t) = u. The existence of optimal control u * in is guaranteed which mainly minimizes the cost functional J 1 .
Proof: Proof has been given in the appendix.
Further, with the help of Pontryagin's Maximum Principle, we characterize the optimal control u * of the system in the following.

Numerical results and discussion
Till now we have analysed the stability conditions (both local and global) of the equilibrium points and also find the optimal control variable for corresponding optimal control problem. The control variable minimizes the total cost that has been considered. Here we have performed numerical simulations to validate our analytical results and also to observe the involvement of control variable in the system dynamics.
We shall numerically solve the control system (8) along with (7) for the set of parameter values provided in Table 4 with the initial size of populations and information as: S(0) = 56.18, P 1 (0) = 5.11, P 2 (0) = 1.43 and T(0) = 0.45. With the help of MATLAB, the graphical scenarios for various cases have been obtained. Forward-backward sweep method has been used to solve the optimal control problem. First we have to solve the optimal state system 'forward in time' and then we need to solve the adjoint state system 'backward in time' . In the next step, these optimal controls are updated using Hamiltonian for optimality of the optimal system and for doing this we need to use the 'steepest descent method' (Kirk, 2012;Wang, 2009). The process will continue until the criterion for convergence is satisfied. Time period of study and application of control interventions is around 1 year, i.e. 12 months. Figure 9 shows the population graphs with respect to time in the absence of the control, i.e. u = 0. For this situation, at T f = 12, the population is (S * , P * 1 , P * 2 , T * ) = (56.3388, 0.1543, 6.0797, 0.8066). It is observed that the growth of psychological addicts rapidly decreases during first 5 months and then slowly decreases after around half a year. It is also noted that there are more addicted population in physiological state than in psychological state. Moreover, there are a significant number of physically addicted population present in the scenario around last few months which will create the economic burden   Table 4.
in terms of productivity loss, morbidity, mortality and in procuring protective measures during the period. Now we shall discuss about the effect of control intervention. The positive weights have been considered as w 1 = 1.6 and w 2 = 10 ( Gaff & Schaefer, 2009;Kassa & Ouhinou, 2015). In Figure 11, the population profiles have been considered by taking the counselling process (u) as control parameter. The optimality of the system has been determined also. For this situation, at T f = 12, the population is (S * , P * 1 , P * 2 , T * ) = (56.3387, 0.1539, 6.0801, 0.8065). Parameters have been taken from Table 4 along with weights as mentioned. The counselling process works better to prevent the urge among population of taking drugs regularly. The addicted population (psychological) are lower than the case when control has not been incorporated. Also the recovery rate from the addiction is quite higher in the early months in this case. The pictures show that there will more physiologically addicted individuals if the addicted stop taking treatment and counselling during their treatment period. It is true that the individuals who are psychologically addicted will become cautious through regular counselling and this will lead to a lower growth rate of physiologically addicted also. The corresponding intensity profile of optimal counselling policy has been shown in Figure 10. It can be observed that in early stage the control intervention work with its higher intensity but later the intensity decreases gradually. It proves that a certain time is needed to convince a psychologically addicted person that up taking drugs in a regular basis is harmful and even can cause physical damages. But once a person start understanding the affects, it becomes easy to make them agree to take medicines and to do the other needful to get rid of this addiction.
We have performed the cost design analysis for optimal control policy as ' cost effectiveness is one of the important characteristics to decide the fitness' (see Figure 12). Cost profile along with the trajectory of psychologically addicted individuals for the optimal control policy can be observed from Figure 12(a ,b). It is trivial when the control policy is used, the optimal cost will be comparatively low than the situation when no control is applied. Smaller number of addicted individuals will reduce the overall total cost via opportunity loss in this case.  Table 4. Figure 11. Profiles of populations with optimal control u * . Parameters are as in Table 4.    Table 4.

Effect of saturated treatment (through counselling) on optimal control variable
This section contains how the counselling resources effect the system dynamics when the control policy is applied with optimal intensity. The saturation on treatment (ξ ) as well as the treatment rate ( ) have been varied in the subsequent figures. Figure 13 depict the profiles of psychologically addicted population and corresponding cost for various values of ξ . Associated optimal control has been drawn in Figure 14. Increasing ξ provides an increased value of psychologically addicted and associated cost also, because higher level of addicted population is related to productivity loss also. From the optimal profile of control variable, it has been concluded that higher saturation in treatment (i.e. smaller value of ξ) will provide lower efforts on the control variable. Therefore, if the saturation in counselling   Table 4. is sufficiently high, then it will be economically viable and needs comparatively lesser efforts to implement the control.
On the other hand, when treatment rate ( ) is varied from 0.1 to 1, the addicted (psychologically) population decreases following a reduction in associated cost also (Figure 15). From Figure 16, it can be observed that with increasing value of , the maximum intensity increases for the control that represents the treatment via counselling. Also, higher treatment rate helps the optimal treatment to work with comparatively higher intensity.

Conclusion
Synthetic drugs are created by man-made chemicals, but not by the natural ingredients. For instance, K2 (Spice), Ecstasy (Molly) and bath salts are synthetic drugs. Synthetic drugs are new types of artificial drugs which are quite acceptable among the drug addicts. Spreading of synthetic drugs has its own characteristics unlike traditional drugs. For example, people believe that the traditional drugs have more toxicity than the synthetic ones, so if someone sucks the synthetic drug a little, he will not become addict. Based on this wrong psychology, many people forget to take the proper precautions and be captured by synthetic drugs gradually.
In this work, we have proposed a model where people is addicted to drugs both psychologically and physiologically. Next generation matrix method helps us to find the basic reproduction number R 0 and this R 0 gives (or consistent with) the local and global stability conditions of the drug-free and drug addiction equilibria. The stability analysis gives us that for R 0 < 1 1+a , the drug-free equilibrium is globally asymptotically stable but for R 0 > 1 1+a , the drug addiction equilibrium is globally asymptotically stable. It means for changing the value of R 0 from more than 1 1+a to less than 1 1+a , the drug addiction could be eliminated. Moreover, analysing the sensitivity of parameters γ , ρ and β about R 0 we have reached to the conclusion that controlling the spread of the synthetic drugs is better than giving treatment to addicts. So, designing the rational measures to control the drug transmission is important.
The next section of this paper is about an optimal control problem relative to the drug abuse epidemic model so as to minimize the drug addiction as well as to minimize the cost of treatment. We have redefined our model by taking the effect of counselling and awareness campaigns as control variable and determined the total cost. The existence of optimal control function has also been proved. Pontryagin's Maximum Principle helps to determine the analytical characterization of the optimal control path. The analytical results and simulations are quite meaningful as the work presented here is dealt with current addiction towards drugs. Through the numerical computations, we can deduce certain observations that have already been discussed.
These days a large number of population, especially the young generation is exposed to world of drugs due to various reasons. These populations can be considered as the focal point for counselling. Because by considering them as susceptible population, it is comparatively easier to assess how best to introduce regular counselling to the psychological addicts in the model. Schools and families should aware the youths about the importance of health education as well as the Government also has a duty to initiate some serious steps to increase the awareness among the people. With the help of campaigns and social programmes, people may realize the danger of synthetic drugs and reduce the curiosity, which could lead to a lower contact rate. The model qualifies the impact of counselling on psychological addicted population because through numerical simulations we have come to the conclusion that the counselling therapy can potentially reduce the addicted population from taking drugs. Moreover, the effect of optimal response due to counselling can minimize the cost burden as well as the number of addicted individuals. The policy can minimize the overall economic load also. So, implementation of a proper control policy will be effective as well as economical.
From the above findings along with the characteristics of control set , we have which is equivalent as (10).