Optimal Control Strategy Applied to Dynamic Model of Drug Abuse Incident for Reducing Its Adverse Effects ========================================================================================================= * Md. Azmir Ibne Islam * Md. Haider Ali Biswas ## Abstract The present world is facing a devastating reality as drug abuse prevails in every comer of a society. The progress of a country is obstructed due to the excessive practice of taking drugs by the young generation. Like other countries, Bangladesh is also facing this dreadful situation. The multiple use of drug substances leads an individual to a sorrowful destination and for this reason, the natural behavior of human mind is disrupted. An addicted individual may regain his normal life by proper monitoring and treatment. The objective of this study is to analyze a mathematical model on the dynamics of drug abuse in the perspective of Bangladesh and reduce the harmful consequences with effective control policies using the idea of optimal control theory. The model has been solved analytically introducing a specific optimal goal. Numerical simulations have also been performed to review the behaviors of analytical findings. The analytical results have been verified with the numerical simulations. The analysis of this paper shows that it is possible to control drug addiction if there is less interaction among general people with the addicted individuals. Family based care, proper medical treatment, awareness and educational programs can be the most effective ways to reduce the adverse effects of drug addiction in a shortest possible time. Keywords * Drug abuse or addiction * Mathematical model * Optimal control * Pontryagin’s maximum principle ## 1. Introduction Drug addiction has become a major problem nowadays. The present world is experiencing the devastating effects of drug abuse day by day. It works as a hindrance to present civilization. It can also restrain us from accomplishing our goals or dreams of life. Generally, the use of substances that are injurious to health is known as drug abuse. It is believed that substance abuse induces sleep or produces excitement or insensibility which has a negative effect on the mental and physical health of an individual. Addiction is not only confined with one single cause but also a number of factors (genetic, brain chemistry, environmental issues and psychological behaviors) are working together behind the practice of consuming drugs. As addiction can trap anyone, people are greatly influenced by the social factors around themselves. The developments in human and national economic growth are about to decline as consuming of illicit drug abuse is increasing rapidly. Drugs have very harmful consequences both for human body and mind. Drug addiction is detrimental to health, even sometimes leads consumers to death. Also the structure and function of brain are greatly hampered due to seeking drugs randomly. As drug addiction can decrease self-control and human ability to make sound decisions, the world is now undergoing a big challenge to meet this crisis. But it can be managed effectively like other chronic and infectious diseases. Fortunately, there are treatments that can help people to counteract with the disruptive effects of addiction and regain control of lives. Research shows that treatment in rehabilitation centers ensures a better life for the patients who are habitual drug addicts [30]. An addict can achieve a sustained recovery in his or her future life without drugs after getting treatment from rehabilitation centers. Drug addiction is significantly gaining attention in whole world and increasing different kinds of social and public health problems nowadays. Drug addiction can be controlled by raising public awareness among the addicts, introducing voluntary activities and enforcing laws against the dealers. Bangladesh has achieved a significant success in socioeconomic development in recent years. Diverse challenges have been taken to increase the life expectancy of people, reduce population growth, increase educational facilities and improve the health of maternal and child. Though the country is advancing with its valuable efforts, another problem has recently arisen which is working as an impediment towards those developments. Such specific problem that becomes the major threat for our country is known as ‘drug addiction’. This human devastation agent has spread its steps to every nook and corner across the country. Bangladesh, though not a drug producing country, is vulnerable for drug abuse because of its geographical location [1]. The tendency of drug addiction can mostly be seen among the young generations who are the leaders of future world. As drug addiction is increasing rapidly throughout the country, the youngsters are getting closure to addiction more and more, ruining their valuable lives and tending the country to a blind future. Many factors such as family problem, social factors, financial and legal issues, school or college environments, physical inability, failure of romance, job stress, poor self-esteem, high level sports competition, loneliness, low cost and easy access to drugs and local trafficking zone etc. are acting behind this worst situation. Mathematical modeling has provided a quantitative insight into diverse fields. Several models have been published in recent years considering various problems regarding practical life. The literature and development of mathematical epidemiology are well documented and can be found in [3, 5, 6, 14, 28, 42, 44, 51]. Mathematical modeling has also launched its steps to explain different social problems and find out the alternative ways for reducing those problems. On the other hand, optimal control theory is another useful part of mathematics that describes the control policies [4, 7, 8, 12, 27, 40, 48, 52, 53]. Mathematical Model and Optimal Control are closely related to each other [40]. Optimal Control deals with the problem of finding a control law for a system so that a certain optimality criterion is achieved by maximizing or minimizing a particular cost function. The formulation of an optimal control problem requires a mathematical model to be controlled. In short, a mathematical model works as a framework for an optimal control problem [40]. Drug addiction has become an important issue in recent years as it ruins the future of young generation. A model on heroin epidemics describes the global behaviour with distributed time delays. Heroin is fairly responsible for the transmission of human immunodeficiency virus (HIV) and hepatitis C virus in China [6, 41]. The impact and adverse effects of an illicit drug use have been analyzed and further extended with two optimal strategies. The objective of the first optimal approach is to reduce the social influence between the susceptible and drug users, while the second aim is to increase the treatment facilities in rehabilitation [48]. A perspective, based on the principles of mathematical epidemiology, shows the consequences of opiate addiction in Europe specifically in Ireland. Individuals, who have a habitual tendency to use opiate drug, are responsible for the destruction of their own physical fitness and mental abilities and also put their families into a great trouble. Later it has been realized that prevention is indeed better than cure [67]. Computational models on the basis of neurobiology of drug abuse have been formulated with the help of mathematical modeling in order to illustrate the potential contribution due to addiction [2]. A model correlated with drug abuse has recently showed that the repeated use of drug substance is merely responsible for the increase of infectious diseases around the world [58] while on the other hand, a tobacco model has been established to obtain the periodic motion from tobacco addiction with complete recovery, relapse rate and addiction stages [3]. In this paper, we have developed a model on the consequences of drug abuse. The desire of this paper is to discover the fact that how the tendency of drug seeking mentality is prevailing among the susceptibles and continuing its chain to covert an individual into a light and heavy addict. We have also established the optimal control policy for minimizing drug addiction so that the young generation may become conscious about the adverse consequences of this deadly drug substance. The valuable efforts of rehabilitation centers can play a vital role in controlling yaba addiction by giving proper treatment to the addicted fellows. The main aim of our present paper is to limit the use of drugs, reduce and minimize the harmful effects of drug abuse as well as take an effective effort to maintain an addiction free environment. The specific objectives of this study are to propose a mathematical model on the dynamics of drug abuse and formulate an optimal control strategy. As drug addiction is increasing rapidly throughout the country, some effective plans are needed for reducing the harmful consequences. In this stage, our model will a play a significant role and can also motivate the general people. Our main purpose is to increase consciousness among young generation about the adverse effects of drug addiction, raise public awareness and minimize the harmful consequences so that this unfavorable situation can be controlled. ## 2. Model Formulation A population size *N*(*t*) is subdivided into five compartments: susceptible compartment *S*(*t*) (individuals who are not drug users, but at a high risk of taking drugs), light drug users *L*(*t*), heavy drug users *H*(*t*), drug users under treatment in rehabilitation R(t) and quitter population Q(t) (individuals who are fully determined not to take drug substances anymore). The schematic diagram of the model is shown in Figure 2.1. ![Figure 2.1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F1.medium.gif) [Figure 2.1.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F1) Figure 2.1. Compartmental representation of drug addiction model Taking the Figure 2.1 into consideration, the model can be represented by the following system of non-linear ordinary differential equations: ![Formula][1] where, *S*(0) > 0,L(0) > 0,*H*(0) > 0,R(0) > 0,Q(0) > 0 and *N = S + L + H + R + Q*. In the above model, the parameter *r* represents the recruitment rate of population, *µ* is the natural mortality rate, *α* is the interaction rate among the susceptibles and light drug users whereas *β* is the effective rate at which light users convert into heavy drug users. The constant *δ* represents the removal rate from addiction without treatment and *γ* is the rate at which heavy addicts are being sent to rehabilitation for treatment. In real life, despite getting treatment from rehabilitation centres, some addicts may again become a heavy drug user. Therefore, it is assumed that *p* is the relapsing rate for which a drug addict under treatment in rehabilitation converts into a heavy user again. Finally, *θ* is the fully recovered rate from addiction. ## 3. Formulation of Optimal Control Model In this section, the drug abuse model has been analyzed newly by introducing an objective function using the well-known optimal control theory. The model is extended to observe the behaviors in the presence of control measures. We have used three control policies in order to minimize the number of drug addicted population. Then the model has been analyzed with several optimality conditions and numerical simulations have also been performed which show an effective outcome that we have expected as our assumption. The schematic diagram of model (2.1) can be reformulated as optimal control form as shown in Figure 3.1. ![Figure 3.1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F2.medium.gif) [Figure 3.1.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F2) Figure 3.1. Compartmental representation of drug addiction model with control Three control variables *u*1, *u*2, *u*3 have been introduced where *u*1 is the awareness and educational programs, *u*2 is the family based care and *u*3 represents the effectiveness of rehabilitation centers. It is to be mentioned that when the effectiveness of rehabilitation centers increases, the relapse case decreases. After considering these three control policies, the model (2.1) can be written in optimal control problem as ![Formula][2] Here the objective functional is ![Formula][3] where, *A*, *B* and *C* are the weight parameters balancing the costs for reducing the addicted individuals. The dynamic constraints can be reformulated as ![Formula][4] with the initial conditions, *S*(0) = *S*, *L*(0) = *L*, *H*(0) = *H*, *R*(0) = *R*, *Q*(0) = *Q* Equation (3.1) is the Lagrange type objective functional, where ![Graphic][5] and *g*(*t*, *x*, *u*) is the right side of the system of differential equations (2.1) updated with control variables i.e. *x*(*t*) = (*S*(*t*), *L*(*t*), *H* (*t*), *R*(*t*), *Q*(*t*)), *u* (*t*) = (*u*1(t), *u*2(t), *u*3(t)) and *g* (*t*, *x*, *u*) = (S′, *L*′, *H*′, *R*′, *Q*′) In other words, the system of differential equations (3.2) along with the objective functional (3.1) describe the controlled model. Now, we have to find a set for optimal control strategies ![Graphic][6] such that ![Formula][7] where, the control set ![Graphic][8] ## 4. The necessary and sufficient conditions of optimal control In this section, we are going to study the sufficient condition for the existence of an optimal control of the system of ordinary differential equations (3.1) [23]. Later on, Pontryagin’s Maximum Principle is used to characterize the optimal control functions and derive the necessary conditions for our control problem [40, 55]. ### 4.1. The existence of an optimal control In this subsection, the existence of an optimal control for the model (3.1) will be examined. In order to establish the existence of an optimal control and the uniqueness of the optimality system, the boundedness of solutions to the system (2.1) for a finite time interval is needed. We begin by examining the priori boundedness of the state solutions of model (2.1). For the model system (2.1) to be epidemiologically meaningful, all the state variables should be nonnegative for all time. We need to show that the total population is bounded for all *t* ≥ 0. From model (2.1), the rate of change of total population is ![Formula][9] In the absence of drug addiction, we obtain ![Formula][10] Now solving this, we have ![Formula][11] From (4.2), it is clear that the total population *N*(*t*) will approach the threshold ![Graphic][12] as *t* → *∞*. This implies that if the initial total population *N* is less than ![Graphic][13] i.e. if ![Graphic][14] then ![Graphic][15]. So, definitely ![Graphic][16] is the upper bound of N. On the other hand, if ![Graphic][17], then *N*(*t*) will decrease to ![Graphic][18] as *t* → *∞*. This means that if ![Graphic][19] then the solution (*S*(*t*), *L*(*t*), *H*(*t*), *R*(*t*), *Q*(*t*)) enters the region Ω or approaches it asymptotically. So it can be concluded that the region Ω is positively invariant under the flow induced by the model (2.1). Therefore, the model is both mathematically and epidemiologically well-posed in the region Ω. It is therefore sufficient to study the dynamics of the model (2.1) in Ω. ### Theorem 4.1. There exists an optimal control ![Graphic][20] to problem (3.1). ### Proof. To prove this theorem, the concept described in [23] has been followed. Considering ![Graphic][21] as the right-hand side of constraint, we would like to satisfy the following conditions: 1. *q* is of class C1 and there exists a constant *c* such that ![Formula][22] 2. The admissible set Q of all solutions to system (3.2) with corresponding control in *U* is nonempty. 3. ![Graphic][23] 4. The control set *U =* [0,1] × [0,1] × [0,1] is closed, convex and compact. 5. The integrand of the objective functional is convex in U. To verify the above conditions, we have ![Formula][24] Therefore ![Graphic][25] is of class C1 and ![Graphic][26]. ![Formula][27] ![Formula][28] Since *S*(*t*), *L*(*t*), *H*(*t*), *R*(*t*), *Q*(*t*) are bounded, there exists a constant *c* such that ![Formula][29] So, condition (1) is satisfied. From condition (1), it is easily visible that there exists a unique solution to the system (3.2) for a constant control, which further implies that condition (2) holds. ![Formula][30] therefore, condition (3) also holds. Condition (4) is obvious from the definition and the proof of this theorem can be completed by verifying condition (5). In order to show the convexity of the integrand in the objective functional ![Graphic][31], we have to prove that ![Formula][32] where, ![Graphic][33] and ![Graphic][34] and ![Graphic][35] are two control vectors with *ψ* ϵ [0,1]. It follows that ![Formula][36] ![Formula][37] Furthermore, we have ![Graphic][38] ![Formula][39] ![Formula][40] Finally, this is the complete proof. ### 4.2. Characterization of an optimal control Since there is an optimal control for the objective function (3.1) subject to the model system (3.2), our target is to derive the necessary conditions for this optimal control [40, 55]. Using Pontryagin’s Maximum Principle, the Hamiltonian function *Hm* w.r.t (*u*1,*u*2,*u*3) is ![Formula][41] ![Formula][42] where, the variables *pi*, *i =* 1, 2, 3, 4 and 5 are adjoint or co-state variables. The optimality system of equations can be obtained by considering the appropriate partial derivatives of Hamiltonian (4.13) with respect to the associated state variables. The following theorem is a consequence of the maximum principle. ### Theorem 4.2. For an optimal control pair ![Graphic][43] and corresponding solutions to the state system *S***, L***, H***, R***, Q** that minimize the objective functional (3.1), there exists adjoint variables *pi*, *i =* 1, 2, 3, 4 and 5, satisfying the following adjoint equations: ![Formula][44] ![Formula][45] Again, the optimal pair may be characterized by the piecewise continuous functions ![Formula][46] ![Formula][47] ![Formula][48] where, ![Formula][49] and *pi*, *i =* 1, 2, 3, 4 and 5 are the solutions of the adjoint system (4.14). ### Proof. Considering ![Graphic][50] be the optimum values (*S*(*t*)*,L*(*t*)*,H*(*t*)*,R*(*t*)*,Q*(*t*)), the adjoint system (4.14) is obtained directly from the application of Pontryagin’s Maximum principle [40, 55]. The differential equations governing the adjoint variables are obtained by differentiating the Hamiltonian function. ![Graphic][51] where, *Yi* =[*S*, *L*, *H*, *R*, *Q*]and *i=*1, 2, 3, 4 and 5. Or, ![Formula][52] and *pi*(*T*) = 0, *i =* 1, 2, 3, 4 and 5. We simply differentiate the Hamiltonian *Hm* with respect to the controls (*u*1,*u*2,*u*3) at the optimal control functions and then equate the resulting expressions to zero, ![Formula][53] Solving for ![Graphic][54] and ![Graphic][55], we derive the optimal controls ![Graphic][56], ![Graphic][57] and ![Graphic][58] as given by the Theorem 4.1. ![Formula][59] Finally, the optimality system consists of the state system (3.2) with the initial conditions, adjoint system (4.14) with transversality conditions (4.15) and optimality condition (4.16) - (4.18). Thus, we have the following optimality system: ![Formula][60] ![Formula][61] ### 4.3. Uniqueness of the optimality system Due to the boundedness for the state system, the adjoint system has bounded coefficients which is linear in each of the adjoint variables. Therefore, the adjoint system has finite upper bounds. The following proposition is needed for the proof of the uniqueness of the optimality system for the small time window. ### Proposition 4.1. The function *u**(m) = min{ max{*m*, *a*},*b* }which is Lipschitz continuous in *m*, with *a < b* for some fixed nonnegative constants [34]. Now, we are about to prove the uniqueness of the optimality system in the same way as described in [25, 34, 71]. ### Theorem 4.3. The solution to the optimality system (4.22) is unique for sufficiently small *T*. ### Proof. The proof can be done by contradiction. That means, the boundedness for the state system and the adjoint system, the Lipschitz continuity of optimal control functions and some elementary inequalities play key roles to reach the contradiction. Considering (*S*, *L*, *H*, *R*, *Q*, *p*1, *p*2, *p*3,*p*4,*p*5)and ![Graphic][62] are the two non-identical solutions to the optimal system (4.22). In order to show that the two solutions are equivalent, it is convenient to make a change of variables. Let, ![Formula][63] where, *ζ* is a positive constant to be chosen later. After introducing the new variables, the optimality conditions can be converted to ![Formula][64] ![Formula][65] ![Formula][66] ![Formula][67] For simplicity of notation, we generally ignore the dependence on time in the following except in the case that a specific time is intended. Now, the difference of the resulting equations for *xi*, ![Graphic][68], *yi* and ![Graphic][69] has been considered. From the first equation of (4.22), we get ![Formula][70] Similarly, ![Formula][71] It follows from the above two equations that ![Formula][72] Multiplying the equation (4.31) by ![Graphic][73] and then integrating both sides from 0 to *T*, we obtain ![Formula][74] In order to simplify the right-hand expressions of (4.32), we need some elementary inequalities. By the elementary inequality (*a* + *b*)2 ≤ 2(*a*2 + *b*2), we have ![Formula][75] where, *c* depends on bounds for *x*1, ![Graphic][76]. Another inequality is also needed to obtain the specific estimation: If 0 ≤ *b* ≤ *a*, then ![Graphic][77]. Using **Proposition 4.1** and the above inequalities, we have the following estimation of ![Graphic][78]. We are dividing our estimation in five cases as follows. Considering the cases: ![Formula][79] We finally can conclude ![Formula][80] where, ![Graphic][81] depends on bounds for ![Graphic][82], *y*1, *y*2. Based on the above arguments, we find ![Formula][83] where, *M*1, *N*1 are the appropriate upper bounds. Similarly, we can obtain the following inequalities for ![Graphic][84] and ![Graphic][85] with *i* = 2,3,4,5 and *j* =1,2,3,4,5: ![Formula][86] ![Formula][87] ![Formula][88] ![Formula][89] ![Formula][90] ![Formula][91] ![Formula][92] ![Formula][93] ![Formula][94] where, *Mi* (*i =* 1,2,…,8)and *Nj* (*j =* 1,2,…,9) depend upon the coefficients and the bounds of the state variables and co-state variables. It can easily be seen that the coefficients of each integral term is nonnegative for sufficiently small *T*. This implies that ![Formula][95] and ![Formula][96] Hence the solution of (4.22) is unique for small time. So the proof is complete. ## 5. Numerical Simulations of Optimal System In this section, a ***sensitivity analysis*** has been performed at first in order to establish which parameters have a significant influence on the reproduction number. For this analysis, the approach described in [48] has been adopted. Sensitivity analysis is the study of how the solutions change with parameters. Results of sensitivity analysis make it easy to investigate any type of dynamic model. The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives [48]. The basic reproduction number for the model (2.1) is ![Formula][97] Neglecting the recruitment rate *r* and natural mortality rate *µ*, we have made a sensitivity analysis on the basic reproduction number *R* with respect to the remaining parameters *α, β, δ* and *γ*. View this table: [Table 5.1.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/T1) Table 5.1. Sensitivity indices of *R* to parameters Parameters whose sensitivity index values are near -1 or +1 suggest that a change in the magnitudes will create a significant impact on either decreasing or increasing the size of basic reproductive number *R*. From **Table 5.1**, it can easily be observed that the interaction rate parameter *α* is the most sensitive parameter to *R* as it has the highest positive sensitivity index value i.e. + 1.00. Therefore, an increase in the values of interaction rate parameter *α* will increase the size of basic reproductive number *R* which means that drug addiction will prevail for a long time. Similarly, the effective rate *β* with a positive sensitivity index value 0.071085 is also responsible for drug addiction. From Figure 5.1, it is observed that an increase (or decrease) in the magnitude of the parameters *α* and *β* will result in an increase (or decrease) in the prevalence of drug addiction. On the other hand, the parameters *γ* and *δ* are shown to have the greatest potential to reduce the harmful situations of drug addiction when the values of these parameters are maximized. From this analysis, it is realized that the interaction rate and effective rate have a large impact on the prevalence of drug addiction. So, we should focus on these two parameters so that the harmful situation of drug addiction can be minimized. That means, we are successful to figure out that the interaction rate parameter *α* and effective rate parameter *β* are mainly liable for the harmful situation. Now we want to minimize those harmful situations by applying our three controls. For this, we have considered four scenarios. In first scenario, we have considered only the awareness and educational programs (*u*1). In second scenario, only the family based care (*u*2) is taken as control. The effectiveness of rehabilitation centers (*u*3) is considered as a control in third scenario. At last, all the controls are considered simultaneously in the fourth scenario. We have an intension to show that social awareness (awareness programs and family based care) is the best way to solve the situation instead of sending the addicted people to rehab centers. By this, the pressure on the rehab centers will be less and the economical cost that is required in rehab centers for treatment will be minimum. To find the answers of all questions, numerical simulations using optimal controls have been performed by using MATLAB with the set of parameter values given in **Table 5.2**. Due to the deficiency of data on drug addiction, precise adjustment is however detrimental to make. Nevertheless, some of the parameters have been assumed in some realistic point of view and the ranges are considered based on the guidance of previous literatures on drug abuse epidemics. Finally, the purpose of this simulation is to state the behaviors of our model and describe the dynamics of drug addiction. ![Figure 5.1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F3.medium.gif) [Figure 5.1.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F3) Figure 5.1. Effects of parameter variations on *R* View this table: [Table 5.2.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/T2) Table 5.2. Description of parameters and respective values ### First Scenario In this scenario, only the awareness and educational programs (*u*1) is considered whereas the other two controls: family based care (*u*2) and the effectiveness of rehabilitation centers (*u*3) are set to zero (i.e. *u*1 ≠ 0, *u*2 = 0, *u*3 = 0). Then the following figures have been obtained. ### Second Scenario In this scenario, only the family based care (*u*2) is taken in absence of other two controls: awareness and educational programs (*u*1) and the effectiveness of rehabilitation centers (*u*3) are set to zero (i.e. *u*2 ≠ 0, *u*1 = 0, *u*3 = 0). Then the following results have been obtained. ### Third Scenario In this scenario, we activate the control *u*3 (effectiveness of rehabilitation centers) whereas the other two controls: awareness and educational programs (*u*1) and family based care (*u*2) are set to zero (i.e. *u*3 ≠ 0, *u*1 = 0, *u*2 = 0). Then we have obtained the following figures. ### Fourth Scenario In this scenario, all the controls *u*1, *u*2, *u*3 are considered simultaneously (i.e. *u*1, *u*2, *u*3 ≠ 0). Using all these controls combinedly, the dynamics of light addicted, heavy addicted, addicted population in rehab and quitter (recovered) population have been shown in following figures. **Figures 5.2 - 5.17** show the dynamics of all scenarios for the optimal controls u1, u2and *u*3. To describe the whole phenomena in an easy manner, the following bar charts have been constructed. Since the main purpose is to minimize the number of light and heavy addicted population, all the scenarios illustrated before have been analyzed. ![Figure 5.2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F4.medium.gif) [Figure 5.2.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F4) Figure 5.2. Dynamics of light addicted population using awareness and educational programs (*u*1) as optimal control ![Figure 5.3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F5.medium.gif) [Figure 5.3.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F5) Figure 5.3. Dynamics of heavy addicted population using awareness and educational programs (*u*1) as optimal control ![Figure 5.4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F6.medium.gif) [Figure 5.4.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F6) Figure 5.4. Dynamics of addicted population in rehabilitation using awareness and educational programs (*u*1) as optimal control ![Figure 5.5.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F7.medium.gif) [Figure 5.5.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F7) Figure 5.5. Dynamics of quitter population using awareness and educational programs (*u*1) as optimal control ![Figure 5.6.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F8.medium.gif) [Figure 5.6.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F8) Figure 5.6. Dynamics of light addicted population using family based care (*u*2) as optimal control ![Figure 5.7.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F9.medium.gif) [Figure 5.7.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F9) Figure 5.7. Dynamics of heavy addicted population using family based care (*u*2) as optimal control ![Figure 5.8.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F10.medium.gif) [Figure 5.8.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F10) Figure 5.8. Dynamics of addicted population in rehabilitation using family based care (*u*2) as optimal control ![Figure 5.9.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F11.medium.gif) [Figure 5.9.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F11) Figure 5.9. Dynamics of quitter population using family based care (*u*2) as optimal control ![Figure 5.10.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F12.medium.gif) [Figure 5.10.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F12) Figure 5.10. Dynamics of light addicted population using effectiveness of rehabilitation centers (*u*3) as optimal control ![Figure 5.11.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F13.medium.gif) [Figure 5.11.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F13) Figure 5.11. Dynamics of heavy addicted population using effectiveness of rehabilitation centers (*u*3) as optimal control ![Figure 5.12.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F14.medium.gif) [Figure 5.12.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F14) Figure 5.12. Dynamics of addicted population in rehabilitation using effectiveness of rehabilitation centers (*u*3) as optimal control ![Figure 5.13.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F15.medium.gif) [Figure 5.13.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F15) Figure 5.13. Dynamics of quitter population using effectiveness of rehabilitation centers (*u*3) as optimal control ![Figure 5.14.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F16.medium.gif) [Figure 5.14.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F16) Figure 5.14. Dynamics of light addicted population using *u*1, *u*2, *u*3 as optimal control ![Figure 5.15.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F17.medium.gif) [Figure 5.15.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F17) Figure 5.15. Dynamics of heavy addicted population using *u*1, *u*2, *u*3 as optimal control ![Figure 5.16.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F18.medium.gif) [Figure 5.16.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F18) Figure 5.16. Dynamics of addicted population in rehabilitation using *u*1, *u*2, *u*3 as optimal control ![Figure 5.17.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F19.medium.gif) [Figure 5.17.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F19) Figure 5.17. Dynamics of quitter population using *u*1, *u*2, *u*3 as optimal control From **Figure 5.18**, it can be observed that the awareness and educational programs (*u*1) is more effective than the other two controls, when we use only one control. Moreover, family based care (*u*2) can also be useful because family is the place where the addicted fellows can get a better circumstance. But when we use the combination of all controls simultaneously, then it shows a highest impact in minimizing the number of light addicted population. **Figure 5.19** shows that the first scenario can be effective at the time of applying only one control i.e. the awareness and educational programs (*u*1) is more effective to diminish the size of heavy addicted population. But if we apply all the controls, the number of heavy addicted population can be minimized within a shortest possible time. ![Figure 5.18.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F20.medium.gif) [Figure 5.18.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F20) Figure 5.18. Approximate peak value of light addicted population for all scenarios ![Figure 5.19.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F21.medium.gif) [Figure 5.19.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F21) Figure 5.19. Approximate peak value of heavy addicted population for all scenarios From **Figure 5.20**, it is observed that although treatments are available in the rehabilitation centers, these treatments (*u*3) are not adequate enough to minimize the number of addicted patients. So, we must have to apply the other two controls: family based care (*u*2), awareness and educational programs (*u*1) at the very early period to control the whole situation. By this, the pressure on the rehabilitation centers will be minimum and we may get a better result from the very beginning which fulfills our main desire. Moreover, it is clear that a significant change or less addicted patients can be obtained if we use all the controls *u*1, *u*2 and *u*3 at the same time. ![Figure 5.20.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F22.medium.gif) [Figure 5.20.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F22) Figure 5.20. Approximate peak value of addicted population in rehab for all scenarios **Figure 5.21** shows the number of fully recovered population after 4 years. We can see that the approximate value is about 2,000 when there is no control. Applying only one control, the size of the population increases to near about 6,000 in first scenario for *u*1, 4,000 in second scenario for *u*2 and 2,200 in third scenario for *u*3. So at the time of using only one control, the awareness and educational programs *u*1 is better than the other two controls. When all the controls *u*1, *u*2 and *u*3 are applied combinedly, then the recovered population increases to 7,000 in fourth scenario which is more effective. ![Figure 5.21.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F23.medium.gif) [Figure 5.21.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F23) Figure 5.21. Approximate value of quitter (recovered) population after 4 years for all scenarios The profiles of the optimal controls have been displayed in **Figure 5.22**. It is noticed that the profiles of the optimal controls *u*1 and u2 are analogous which indicates that both *u*1 and *u*2 require an initial strong start that should be continued for a greater period of time. That means, awareness and educational program with family based care should be maximized to reduce the harmful situation. On the other hand, the control profile u3 demonstrates similar trends. From the figure, it is revealed that supports from the rehabilitation centers should be increased gradually andffcient treatment services should be provided and maintained throughout the treatment period. ![Figure 5.22.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/06/2020.05.02.20088468/F24.medium.gif) [Figure 5.22.](http://medrxiv.org/content/early/2020/05/06/2020.05.02.20088468/F24) Figure 5.22. Profiles of the optimal controls *u*1, *u*2, *u*3 with weight parameters *A* = 0.5, *B* = 0.4 and *C* = 0.2 From the simulation results, it can be understood that the combination of educational campaigns, family care and effective treatment is very effective and plays a vital role in minimizing the harmful consequences of drug abuse in a shortest possible time compared to the case of without controls. ## 6. Conclusions After performing analytical and numerical simulations, we have found that the control, awareness and educational programs, can be the most effective way to minimize the addicted population in a shortest possible time. The family based care is also a useful medium for controlling the adverse situation. By this, the pressure on the rehabilitation centers will be less as such institutions always work for the betterment of social problems. If the controls, family based care, awareness and educational programs, are possible to apply at the very beginning, the economical cost that is required in rehabilitation centers for treatments will be minimum. If we observe the controls in practical life, we must see that all the controls are effective from their own sides to reduce the harmful effects of drug addiction. Normally, family is the best place where the family members, especially the young generation, may find a suitable circumstance. The parents and guardians always give guidance to their lovable sons, daughters or relatives and take care of them about their studies. Such guidelines must help the young generation not to become addicted. This control works till the children are at home and in touch of their families. But every students has to go outside of their family to colleges or universities for their educational purposes. This is the most crucial time when the young generation may become addicted as they are not in the boundaries of their families. In this situation, our assumed control, awareness and educational programs, can be the most constructive way to prevent them from being addicted. Awareness and educational programs always focus on the harmful effects of taking drug substances so that the young generation may become conscious and stay away from the evil deeds. Despite taking all the initiatives, i.e. applying family based care and awareness programs, anyone may become a drug addicted person. In that situation, he or she might get a proper treatment from the rehabilitation centers. But rehabilitation centers are very rare in rural areas. Some families also feel hesitate to send their addicted relatives to the rehabilitation centers for treatment. Moreover, the treatments in rehabilitation centers are quite expensive. This is the main reason why the effects of this control are less and it takes a long time to minimize the addicted population than the other two controls. That means, a minimum number of addicted population can be found in a shortest possible time if family based care, awareness and educational programs are applied. Moreover, the effects of drug addiction can be controlled in an effective manner if we apply all the controls at the same time. In short, we have formulated a general model showing the phenomena of drug addiction and specified the ways about how to control it from a locality in a shortest possible time. Our model will increase consciousness among the young generation and may help the social workers to make plans regarding drug addicted population. If it is possible to apply all the preventive ways that have been shown in this model, anyone can hope for a better country which will be free from drug addiction. ## Data Availability The data used to support the findings of this study are included within the article. ## Conflicts of Interest The authors declare that there are no conflicts of interest. ## Data Availability The data used to support the findings of this study are included within the article. ## Authors’ Contributions MHAB designed the study, performed the conceptualization and methodological analysis and model formulation of the first draft of the manuscript. MAII analyzed the model analytically and wrote the programming codes to perform the computational analysis and also contributed to literature searches and calculated the real data to estimate the parameters. All authors have read and agreed to the publish the final version of the manuscript. ## Acknowledgment The first author was supported by the Ministry of Science and Technology, Government of the People’s Republic of Bangladesh in his MSc program with the NST Fellowship 2017-18, Grant No. 39.00.0000.012.02.009.17-662, Serial: 128, Date-10.01.2018. This partial financial support is greatly acknowledged. ## Footnotes * Mathematics Subject Classification (MSC 2010): 34A12, 34A34, 34H05, 35B35, 49K15, 93A30, 92D30. * Received May 2, 2020. * Revision received May 2, 2020. * Accepted May 6, 2020. * © 2020, Posted by Cold Spring Harbor Laboratory The copyright holder for this pre-print is the author. All rights reserved. The material may not be redistributed, re-used or adapted without the author's permission. ## References 1. [1].Annual drug report, Department of Narcotics Control, Ministry of Home Affairs, Government of the People’s Republic of Bangladesh, 2016. 2. [2]. S. H. Ahmed, M. Graupner and B. Gutkin, Computational approaches to the neurobiology of drug addiction, Pharmacopsychiatry, 42(1): 144–152, 2009. 3. [3]. Y. Bae, Chaotic Dynamics in Tobacco’s Addiction Model, International Journal of Fuzzy Logic and Intelligent Systems, 14(4): 322–331, 2014. 4. [4]. M. H. A. Biswas, Optimal control of Nipah Virus (NIV) infections: a Bangladesh scenario, Journal of Pure and Applied Mathematics, 12(1): 77–104, 2014. 5. [5]. M. H. A. Biswas, L. T. Paiva and M. D. R. de Pinho, A SEIR model for control of infectious diseases with constraints, Mathematical Biosciences and Engineering, 11(4): 761–784, 2014. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3934/mbe&link_type=DOI) 6. [6]. M. H. A. Biswas, AIDS epidemic worldwide and the millennium development strategies: A light for lives, HIV & AIDS Review, 11(4): 87–94, 2012. 7. [7]. M. H. A. Biswas, Necessary conditions for optimal control problems with and without state constraints: a comparative study, WSEAS Transactions on Systems and Control, 6(6): 217–228, 2011. 8. [8]. M. H. A. Biswas, T. Rahman and N. Haque, Modeling the potential impacts of global climate change in Bangladesh: an optimal control approach, Journal of Fundamental and Applied Sciences, 8(1): 1–19, 2016. 9. [9]. K. Blayneh and Y. Cao, H.D. Kwon, Optimal control of vector-borne diseases: treatment and prevention, Dis. Cont. Dyn. Syst. B, 11(3): 587–611, 2009. 10. [10]. S.M. Blower, A.R. McLean and T.C. Porco, P.M. Small, P.C. Hopewell, M.A. Sanchez, A.R. Moss, The intrinsic transmission dynamics of tuberculosis epidemic, Nat. Med., 1: 815–821, 1995. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nm0895-815&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=7585186&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1995RM63500040&link_type=ISI) 11. [11]. S.M. Blower, P.M. Small and P.C. Hopewell, Control strategies for tuberculosis epidemics: new models for old problems, Science, 273: 497–500, 1996. [Abstract](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIyNzMvNTI3NC80OTciO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8wNS8wNi8yMDIwLjA1LjAyLjIwMDg4NDY4LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 12. [12]. E. Bonyah, K. Badu and S. K. A. Addo, Optimal control application to an Ebola model, Asian Pacific Journal of Tropical Biomedicine, 6(4): 283–289, 2016. 13. [13]. S. Bowong and A.M.A. Alaoui, Optimal interventions strategies for tuberculosis, Commun. Nonlinear Sci. Numer. Simul, 18: 1441–1453, 2013. 14. [14]. D. N. Burghes and M.S. Borrie, Modelling with differential equations, Ellis Horwood Limited, New York, 1981. 15. [15]. L. Cai, X. Li, M. Ghosh and B. Guo, Stability analysis of an HIV/AIDS epidemic model with treatment, Journal of ComPutational and APPliedMathematics, 229(1): 313–323, 2009. 16. [16].1. M.A. Horn, 2. G. Simonett, 3. G. Webb C. Castillo-Chavez and Z. Feng, Mathematical models for the disease dynamics of tuberculosis, in: M.A. Horn, G. Simonett, G. Webb (Eds.), Advances in Mathematical Population Dynamics-Molecules, Cells and Man, Vanderbilt University Press, Nashville, 117–128, 1998. 17. [17]. C. Castillo-Chavez and B. Song, Dynamical models of tuberculosis and their applications, Mathematical Biosciences and Engineering, 1: 361–404, 2004. 18. [18]. C.Y. Chiang and L.W. Riley, Exogenous reinfection in tuberculosis, Lancet. Infect. Dis., 5: 629–636, 2005. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1473-3099(05)70240-1&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=16183517&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000232273700018&link_type=ISI) 19. [19]. S. Choi, E. Jung and S.M. Lee, Optimal intervention strategy for prevention tuberculosis using a smoking-tuberculosis model, J. Theor. Biol., 380: 256–270, 2015. 20. [20]. T. Cohen and M. Murray, Modeling epidemics of multidrug-resistant m. tuberculosis of heterogeneous fitness, Nat. Med, 10 (10): 1117–1121, 2004. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nm1110&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15378056&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000224245800043&link_type=ISI) 21. [21]. C. L. Dym, Principles of mathematical modeling, Second Edition, Academia Press, New York, 2004. 22. [22]. C. Dye, G.P. Garnett, K. Sleeman and B.G. Williams, Prospects for worldwide tuberculosis control under the who dots strategy. directly observed short-course therapy, Lancet., 352(9144): 1886–1891, 1998. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0140-6736(98)03199-7&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=9863786&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000077544600009&link_type=ISI) 23. [23]. W.H. Fleming and R.W. Rishel, Deterministic and Stochastic Optimal Control, Springer-Verlag, New York, 1975. 24. [24]. Z. Feng, C. Castillo-Chavez and A.F. Capurro, A model for tuberculosis with exogenous reinfection, Theor. PoPul. Biol., 57: 235–247, 2000. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1006/tpbi.2000.1451&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10828216&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000087411400004&link_type=ISI) 25. [25]. H. Gaff and E. Schaefer, Optimal control applied to vaccination and treatment strategies for various epidemiological models, Math. Biosci. Eng, 6: 469–492, 2009. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3934/mbe.2009.6.469&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19566121&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) 26. [26]. F.R. Gantmacher, Applications of the theory of matrices, Interscience Publisher, New York, 1959. 27. [27]. D. Gao and N. Huang, Optimal control analysis of a tuberculosis model, APPlied Mathematical Modelling, 58: 47–64, 2018. 28. [28]. R. J. Garten and S Lai, J. Zhang, et al., Rapid transmission of hepatitis C virus among young injecting heroin users in Southern China, International Journal of EPidemiology, 33, 182–188, 2004. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/ije/dyh019&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15075167&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000220615000033&link_type=ISI) 29. [29]. M. Gomes, A. Franco, M. Gomes and G. Medley, The reinfection threshold promotes variability in tuberculosis epidemiology and vaccine efficacy, Proc. R. Soc. B, 271(1539): 617–623, 2004. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rspb.2003.2606&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15156920&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000220568100010&link_type=ISI) 30. [30]. M. Hasan and A. S. M. Shahin, Drug rehabilitation center based survey on drug dependence in Dhaka city, Update Dental College Journal, 3(1): 32–36, 2013. 31. [31]. M. A. I. Islam and M. H. A. Biswas, Mathematical Analysis of Dynamic Model of Drug Addiction in Bangladesh, Abstract Proceedings of the International Conference on Advances in Computational Mathematics (ICACM 2017), Department of Mathematics, University of Dhaka, Bangladesh on May 2728, 2017. 32. [32]. E. Jung, S. Lenhart and Z. Feng, Optimal control of treatments in a two-strain tuberculosis model, Dis. Cont. Dyn. Syst. Ser. B, 2(4): 473–482, 2002. 33. [33]. H. K. Khalil, Nonlinear Systems, Third Edition, Prentice Hall, New Jersey, 2001. 34. [34]. S. Khajanchi and D. Ghosh, The combined effects of optimal control in cancer remission, Appl. Math. Comput., 271: 375–388, 2015. 35. [35]. M. R. Khatun and M. H. A. Biswas, Mathematical modeling applied to renewable fishery management, Mathematical Modelling of Engineering problems, 6(1): 121–128, 2019. 36. [36]. M. S. Khatun and M. H. A. Biswas, Optimal control strategies for preventing hepatitis B infection and reducing chronic liver cirrhosis incidence, Infectious Disease Modelling, 5: 91–100, 2020. 37. [37]. A. Kumar and P.K. Srivastava, Vaccination and treatment as control interventions in an infectious disease model with their cost optimization, Commun. Nonlinear Sci. Numer. Simul., 44: 334–343, 2017. 38. [38]. O. Lavi, M. M. Gottesman and D. Levy, The dynamics of drug resistance: A mathematical perspective, Drug Resistance Updates, 15(1-2): 90–97, 2012. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.drup.2012.01.003&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22387162&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) 39. [39]. J. Lee, J. Kim and H.D. Kwon, Optimal control of an influenza model with seasonal forcing and age-dependent transmission rates, J. Theor. Biol., 317: 310–320, 2013. 40. [40]. S. Lenhart and J. T. Workman, Optimal control applied to biological models, Chapman & Hall, New York, 2007. 41. [41]. J.L. Liu and T.L. Zhang, Global behaviour of a heroin epidemic model with distributed delays, Applied Mathematics Letters, 24(1): 1685–1692, 2011. 42. [42]. J.L. Liu and T.L. Zhang, Global stability for a tuberculosis model, Math. Comput. Model, 54: 836–845, 2011. 43. [43]. M. K. Mondal, M. Hanif and M. H. A Biswas, A mathematical analysis for controlling the spread of Nipah virus infection, International Journal of Modelling and Simulation, 37(3): 185–197, 2017. 44. [44]. M. K. Mondal, M. A. Alim, M. F. Rahman and M. H. A. Biswas, Mathematical analysis of financial model on market price with stochastic volatility, Journal of Mathematical Finance, 7(2): 351–365, 2017. 45. [45]. D.P. Moualeu, M. Weiser, R. Ehrig and P. Deuflhard, Optimal control for a tuberculosis model with undected cases in cameroon, Commun. Nonlinear Sci. Numer. Simul., 20: 986–1003, 2015. 46. [46]. G. Mulone and B. Straughan, A note on heroin epidemics, Mathematical Biosciences, 218, 138–141, 2009. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19563739&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) 47. [47]. S. Mushayabasa and C.P. Bhunu, Modeling the impact of early therapy for latent tuberculosis patients and its optimal control analysis, J. Biol. Phys., 39: 723–747, 2013. 48. [48]. S. Mushayabasa and G. Tapedzesa, Modeling illicit drug use dynamics and its optimal control analysis, Computational and Mathematical Methods in Medicine, 2015(2): 1–11, 2015. 49. [49]. J. Mushanyuzy, F. Nyabadza, G. Muchatibayaz and A. G. R Stewartz, The role of family in initiating methamphetamine abuse treatment: insights through a mathematical model, Journal of aPPlied & computational mathematics, 5(1): 1–10, 2016. 50. [50]. F. Nyabadza, J. B. H. Njagarah and R. J. Smith, Modelling the dynamics of crystal meth (‘Tik’) abuse in the presence of drug-supply chains in South Africa, Bulletin of Mathematical Biology, 75(1): 24–48, 2013. 51. [51]. N. Nyerere, L. S. Luboobi and Y. N. Gyekye, Modeling the Effect of Screening and Treatment on the Transmission of Tuberculosis Infections, Mathematical theory and Modeling, 4(7): 51–62, 2014. 52. [52]. F. T. Oduro, G. Apaaboah and J. Baafi, Optimal Control of Ebola Transmission Dynamics with Interventions, British Journal of Mathematics & Computer Science, 19(1): 1–19, 2016. 53. [53]. O. O. Onyejekwe and E. Z. Kebede, Epidemiological Modeling of Measles Infection with Optimal Control of Vaccination and Supportive Treatment, Applied and Computational Mathematics, 4(4): 264–274, 2015. 54. [54]. L. S Pontryagin, V. G Boltyanskii, R.V Gamkrelidze and E. F Mishchenco, The mathematical theory of optical processes. Wiley, New York, 1962. 55. [55]. L.S. Pontryagin, Mathematical Theory of Optimal Processes, CRC Press, 1987. 56. [56]. P. Rodrigues, C.J. Silva and D.F.M. Torres, Cost-effectiveness analysis of optimal control measures for tuberculosis, Bull. Math. Biol., 76(10): 2627–2645, 2014. 57. [57]. S. L. Ross, Differential equations, Third Edition, Jhon Wiley & Sons, U.K., 2004. 58. [58]. C. Rossi, The role of dynamic modelling in drug abuse epidemiology, Bulletin on Narcotics, 54(1): 33–44, 2002. 59. [59]. S. K. Sahani and M. H. A. Biswas, Mathematical modeling applied to understand the dynamical behavior of HIV infection, Open Journal of Modelling and Simulation, 5, 145–157, 2017. 60. [60]. C.J. Silva and D.F.M. Torres, Optimal control applied to tuberculosis models. the IEA-EEF european congress of epidemiology 2012: epidemiology for a fair and healthy society, Eur. J. Epidemiol, 27: S140–141, 2012. 61. [61]. C.J. Silva and D.F.M. Torres, Optimal control for a tuberculosis model with reinfection and postexposure interventions, Math. Biosci., 244(2): 154–164, 2013. 62. [62]. C.J. Silva and D.F.M. Torres, A TB-HIV/AIDS coinfection model and optimal control treatment, Dis. Cont. Dyn. Syst. A, 35(9): 4639–4663, 2015. 63. [63]. S. H. Strogatz, Nonlinear dynamics and chaos, First Edition, Perseus books publishing, New York, 1994. 64. [64]. P. Van den Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences, 180: 29–48, 2002. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0025-5564(02)00108-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12387915&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000179220600004&link_type=ISI) 65. [65]. S. Verver, R.M. Warren, N. Beyers, M. Richardson, G.D. van der Spuy, M.W. Borgdorff, D.A. Enarson, M.A. Behr and P.D. van Helden, Rate of reinfection tuberculosis after successful treatment is higher than rate of new tuberculosis, Am. J. Respir. Crit. Care. Med, 171: 1430–1435, 2005. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1164/rccm.200409-1200OC&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15831840&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000229711200017&link_type=ISI) 66. [66]. E. Vynnycky and P.E. Fine, The natural history of tuberculosis: the implications of age-dependent risks of disease and the role of reinfection, Epidemiol. Infect., 119(2): 183–201, 1997. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1017/S0950268897007917&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=9363017&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1997YD95900010&link_type=ISI) 67. [67]. E. White and C. Comiskey, Heroin epidemics, treatment and ODE modeling, Mathematical Biosciences, 208(1): 312–324, 2007. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17174346&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F06%2F2020.05.02.20088468.atom) 68. [68].World Drug Report, United Nations office on drugs and crime, New York, 2016. 69. [69]. Y.L. Yang, S.Y. Tang, X.H. Ren, H.W. Zhao, C.P. Guo, Global stability and optimal control for a tuberculosis model with vaccination and treatment, Dis. Cont. Dyn. Syst. Ser. B, 21: 1009–1022, 2016. 70. [70]. I. Zeiler, J.P. Caulkins, D. Grass, G. Tragler, Keeping options open:an optimal control model with trajectories that reach a DNSS point in positive time, SIAM J. Control. Optim., 48 (6): 3698–3707, 2010. 71. [71]. Y.G. Zhou, Y.T. Liang, J.H. Wu, An optimal strategy for HIV multitherapy, J. Comput. Appl. Math., 263: 326–337, 2014. 72. [72].[https://docs.google.com/forms/d/e/1FAIpQLSdGoafSkSXptwI9BMJRN8wl6j3WpXkfck9Xkql8CJ0-Qy\_kKA/viewform](https://docs.google.com/forms/d/e/1FAIpQLSdGoafSkSXptwI9BMJRN8wl6j3WpXkfck9Xkql8CJ0-Qy_kKA/viewform). [1]: /embed/graphic-2.gif [2]: /embed/graphic-4.gif [3]: /embed/graphic-5.gif [4]: /embed/graphic-6.gif [5]: /embed/inline-graphic-1.gif [6]: /embed/inline-graphic-2.gif [7]: /embed/graphic-7.gif [8]: /embed/inline-graphic-3.gif [9]: /embed/graphic-8.gif [10]: /embed/graphic-9.gif [11]: /embed/graphic-10.gif [12]: /embed/inline-graphic-4.gif [13]: /embed/inline-graphic-5.gif [14]: /embed/inline-graphic-6.gif [15]: /embed/inline-graphic-7.gif [16]: /embed/inline-graphic-8.gif [17]: /embed/inline-graphic-9.gif [18]: /embed/inline-graphic-10.gif [19]: /embed/inline-graphic-11.gif [20]: /embed/inline-graphic-12.gif [21]: /embed/inline-graphic-13.gif [22]: /embed/graphic-11.gif [23]: /embed/inline-graphic-14.gif [24]: /embed/graphic-12.gif [25]: /embed/inline-graphic-15.gif [26]: /embed/inline-graphic-16.gif [27]: /embed/graphic-13.gif [28]: /embed/graphic-14.gif [29]: /embed/graphic-15.gif [30]: /embed/graphic-16.gif [31]: /embed/inline-graphic-17.gif [32]: /embed/graphic-17.gif [33]: /embed/inline-graphic-18.gif [34]: /embed/inline-graphic-19.gif [35]: /embed/inline-graphic-20.gif [36]: /embed/graphic-18.gif [37]: /embed/graphic-19.gif [38]: /embed/inline-graphic-21.gif [39]: /embed/graphic-20.gif [40]: /embed/graphic-21.gif [41]: /embed/graphic-22.gif [42]: /embed/graphic-23.gif [43]: /embed/inline-graphic-22.gif [44]: /embed/graphic-24.gif [45]: /embed/graphic-25.gif [46]: /embed/graphic-26.gif [47]: /embed/graphic-27.gif [48]: /embed/graphic-28.gif [49]: /embed/graphic-29.gif [50]: /embed/inline-graphic-23.gif [51]: /embed/inline-graphic-24.gif [52]: /embed/graphic-30.gif [53]: /embed/graphic-31.gif [54]: /embed/inline-graphic-25.gif [55]: /embed/inline-graphic-26.gif [56]: /embed/inline-graphic-27.gif [57]: /embed/inline-graphic-28.gif [58]: /embed/inline-graphic-29.gif [59]: /embed/graphic-32.gif [60]: /embed/graphic-33.gif [61]: /embed/graphic-34.gif [62]: /embed/inline-graphic-30.gif [63]: /embed/graphic-35.gif [64]: /embed/graphic-36.gif [65]: /embed/graphic-37.gif [66]: /embed/graphic-38.gif [67]: /embed/graphic-39.gif [68]: /embed/inline-graphic-31.gif [69]: /embed/inline-graphic-32.gif [70]: /embed/graphic-40.gif [71]: /embed/graphic-41.gif [72]: /embed/graphic-42.gif [73]: /embed/inline-graphic-33.gif [74]: /embed/graphic-43.gif [75]: /embed/graphic-44.gif [76]: /embed/inline-graphic-34.gif [77]: /embed/inline-graphic-35.gif [78]: /embed/inline-graphic-36.gif [79]: /embed/graphic-45.gif [80]: /embed/graphic-46.gif [81]: /embed/inline-graphic-37.gif [82]: /embed/inline-graphic-38.gif [83]: /embed/graphic-47.gif [84]: /embed/inline-graphic-39.gif [85]: /embed/inline-graphic-40.gif [86]: /embed/graphic-48.gif [87]: /embed/graphic-49.gif [88]: /embed/graphic-50.gif [89]: /embed/graphic-51.gif [90]: /embed/graphic-52.gif [91]: /embed/graphic-53.gif [92]: /embed/graphic-54.gif [93]: /embed/graphic-55.gif [94]: /embed/graphic-56.gif [95]: /embed/graphic-57.gif [96]: /embed/graphic-58.gif [97]: /embed/graphic-59.gif