THE DETERMINISTIC EVOLUTION OF ILLICIT DRUG CONSUMPTION WITHIN A GIVEN POPULATION

We study the NERA model that describes the dynamic evolution of illicit drug usage in a population. The model consists of nonusers (N) and three categories of drug users: the experimental (E) category, the recreational (R) category and the addict (A) category. Two epidemic threshold terms known as the reproduction numbers, R0 and μ are defined and derived. Sensitivity analysis of R0 on the parameters are performed in order to determine their relative importance to illicit drug prevalence. The local and global stability of the equilibrium states are also analysed. We also prove that a transcritical bifurcation occurs at R0 = 1. It is shown that an effective campaign of prevention can help to fight against the prevalence of illicit drug consumption. We demonstrate persistence when R0 > 1 and conditions for the extinction of drug consumption are also established. Numerical simulations are performed to verify our model. Our results show that the NERA model can assist policy makers in targeting prevention for maximum effectiveness and can be used to adopt evidence-based policies to better monitor and quantify drug use trends.


Introduction
Illicit drug use is an insidious and global problem, affecting valuable human lives which has led to inestimable harm on public health as well as both social and legal issues.The United Nations Office on Drug and Crime (UNODC) defines illicit drugs as drugs which are under international control which may or may not have licit medical purposes but which are produced, trafficked or consumed illicitly [17].These drugs can either be naturally occurring, semi synthetic (chemical manipulations of substances extracted from natural materials) or synthetic (created entirely by laboratory manipulation) such as the amphetamine-type including MDMA (ecstasy) or cocaine stimulants, hallucinogens like cannabis, and the opioids or narcotics like heroin or morphine and the sedative hypnotics (benzodiazepines, barbiturates) [17].Today, drugs are a mass phenomenon that affect all social strata, constituting a new critical issue in the context of the alarming spread of HIV/AIDS and other blood-borne infectious diseases such as Hepatitis B and C. Numerous countries have been challenged with a war against the distressing growth of illicit drug use since it possesses morbidity, mortality and substantial socio economic impact through its neuro-psychological effects on quality of life of "victims".Globally, it is estimated that 1 in 20 adults, aged 15-64 years used at least one drug in 2014.Drug-related deaths remain unacceptably high and in 2014, an estimated 207400 deaths for people aged 15-64 were reported worldwide (World Drug Report, 2016).Consequently, illicit drug use represents a complex social and health problem

Some Existing Dynamic Models of Illicit Drug Consumption
Illicit drug consumption is a global menace and presents challenging problems for policy makers.There is lack of evidence upon which to base policies, results from studies done are not being used to the maximum in drug prevention strategies.There seems to be an absence of adequate approaches or models to help policymakers in right decision taking.Thus, drug policy is a highly complicated and politicized arena.For decision takers, it is essential to quantify the trend in the evolution of the different categories of drug users.They also need to be able to identify how each category responds to drug control interventions in order to adopt evidence-based policies.
According to [9], dynamic drug models are an interesting approach to understand this phenomenon via scenario analysis, thereby providing a predictive tool for how drug takers behave and as such becoming a useful device to aid specialist teams in devising appropriate intervention measures and predicting drug use trends.They also provide a means of integrating data from different sources, describing a process to increase understanding and simulating experiments that are practically and ethically not possible in real life.As stated in [1], epidemiological data can just describe the phenomenon retrospectively whereas dynamic modeling is used to forecast possible future developments.For this reason, several mathematical models have been developed to explore the dynamics of drug consumption like an infectious disease that spreads through social contact or peer pressure.The One-Dimensional Drug Model by [8] is the simplest of all possible drug prevalence models.This model is developed with only one state variable, which represents the total number of drug users at a particular time t.It also consists of two control variables: u(t), the expenditures for treatment at time t and v(t), the law enforcement measures.However, this one-state model does not categories the different types of drug users.The LH model of [8] and [2] was based on three groups, namely, heavy users, denoted by H, light users, denoted by L and the group of nonusers.In this model, the number of nonusers is assumed to be very large and constant so that it does not need to be modeled explicitly.Moreover, this work assumes initiation in the transition of nonusers to light users, the rate of escalation of light users to the heavy users category and the rate of desistance of heavy users.A significant limitation in this study is that the model assumes the flow of heavy users to light users, which is not applicable in consumption of many illicit drugs.Furthermore, other authors revisited this study and other models have been drawn such as the initiation model, the model of U.S cocaine consumption and the model of controlled drug demand.In [10], the authors proposed an epidemiological model to explore the dynamics of drinking behaviour.As the spread of epidemic infections and that of drinking have much in common, a simple SDR model was described which consists of three compartments, namely, occasional and moderate drinkers (S), problem drinkers (D), and temporarily recovered (R).Recently, [11] developed a new mathematical model of alcohol abuse with four groups of people.Unlike the work of [10], the authors have included another group, namely, drinkers in treatment.In addition, both models considered moderate and occasional drinkers in the same compartment.However, they could have categorized the two types of drinkers separately so that the flow of drinkers from occasional drinkers to moderate drinkers category could have been analysed.In the present study, we consider the dynamical evolution of three categories of drug users, namely, experimental, recreational and addicts, and the nonusers of age 15-30.This age group has been chosen because addiction to drugs is more common among people in their mid teens to their mid-20s [18].According to the Centers for Disease Control and Prevention, those Americans aged from 15 to 34 are at higher risk of developing major depression [19,21].We have therefore concluded that it is more probable for someone aged between 15 − 30 to fall prey to illicit drug use, consistent with the majority of data from different countries worldwide.
We revisit the NERA model presented by [1] and [7].We present a complete, simple and realistic qualitative analysis of the model.In the next section, we describe the system of ordinary differential equations for the deterministic NERA model as further described.A mathematical analysis of the model is done in section 3. Numerical experiments illustrating the results of the NERA model are performed in section 4 and eventually, we provide the concluding remarks.

Model Formulation
The population is partitioned into four categories, namely, nonuser, experimental, recreational and addict.The proportion of the four categories at time t are denoted by N (t) + E(t) + R(t) + A(t) = 1.The schematic representation of the model is shown in figure 1.The different arrows represent the transition of an individual from one category to another.Figure 1 gives a realistic point of view as nonusers represent a larger part of the population, while addicts make up the smaller portion of the population.We note that movement of the individuals from one class to another is irreversible within the model.
Moreover, in our model, we consider the movement of individuals into the nonuser category only, whereas movement out will be from nonuser category as well as drug user category.However, we ignore the fact that experimental users can become addict directly because to be more realistic, there are not many individuals who become dependent on drugs at an early age of drug taking [22].We consider movement into and out of the population β.This movement is a result of individuals becoming too old to still form part of the population (thus moving out) or has become mature enough and introduced in the population (thus moving in).We note that movement out of population due to death is not regarded in our model since death rate is negligible as compared to the aging rate in this population [23].We assume that the rate at which individuals are moving into and out of the population are the same, that is, we assume that the population of individuals of age 15 − 30 does not experience serious fluctuations, thus keeping β constant.In our model, we consider the mutual influence that drug users can have on nonusers and on each other, thereby introducing nonlinear terms N (t)R(t), N (t)E(t) and R(t)E(t) in the model [1,7].The following system of ordinary differential equations describes the dynamics of N (t), E(t), R(t) and A(t).
The physical meaning of each parameter used in the model is given in Table 1.Rate at which recreational users change to addicts r 5 Rate at which recreational users quit drugs r 6 Rate at which addicts quit drugs β Rate of moving in or out of the population due to ageing

Evaluation of the two reproduction numbers, R 0 and µ
The first reproduction number, denoted by R 0 , is the expected number of secondary drug users produced by a typical drug user in the population.Let F i (X) represent the arrival of new drug users in class i and , where V + i (X) and V − i (X) be the rate of transfer of individuals into and out of class i respectively.Let X = (N, E, R, A) T and we suppose x = (E, R) T , where x denotes the set of infective classes.
Then model ( 1) can be written as where The Jacobian related to F(x) and V(x) respectively are given as Linearising matrices F and V at the drug free equilibrium, P 0 = (1, 0, 0, 0), we get, Therefore, where F V −1 is the next generation matrix of model (1).
According to Theorem 2 in [12], it follows that the spectral radius of F V −1 gives R 0 .Hence, the first reproduction number of model ( 1) is given by For this model, the basic reproduction number R 0 can be defined as the average number of individuals from the nonuser category, who are influenced by a member of the drug user population to become experimental users and recreational users.In [7], it was noted that the drug-free equilibrium is asymptotically stable as all eigenvalues are negative for R 0 < 1.It is further mentioned that the critical point (1, 0, 0, 0) changes from an asymptotically stable node at γ < 1 to an unstable saddle point at γ > 1, where γ = (r1−r3) β .We indeed prove the existence of a transcritical bifurcation at R 0 = 1, later on in the present work.
We next derive the second reproduction number, denoted by µ. µ is defined as the expected number of secondary drug users (recreational users) produced by a typical drug user of the experimental category.
Correspondingly, the matrices F and V are then given as At recreational-addict free equilibrium, , 0, 0 , F and V are given by Thus, We note that µ is the proportion of recreational users that exert an influence on nonusers via the experimental users before leaving the category [7].

Sensitivity analysis of the basic reproduction number, R 0
We study the sensitivity of R 0 with respect to each of its parameters using the normalized forward index approach.We calculate the sensitivity indices of the reproduction number and these indices tell us how crucial each parameter is to drug prevalence.Sensitivity analysis is commonly used to determine the robustness of model predictions to parameter values.Following [5], the sensitivity indices of R 0 are calculated.
From the above calculations, we observe that the basic reproduction number R 0 is directly proportional to the parameter r 1 , that is, R 0 is most sensitive to changes in r 1 .In other words, if r 1 increases, R 0 will also increase with same proportion and if r 1 decreases, R 0 will also decrease with same proportion.We note that the parameters β and r 3 have an inverse relationship with respect to R 0 , so an increase in any of them results in a decrease in R 0 .
β is the rate of moving in or out of the population due to ageing and thus, it is clear that an increase in this rate is neither ethical nor practical.Therefore, our aim is to focus on one of the following two parameters: either r 1 , the influence rate of experimental and recreational on nonusers or r 3 , the rate at which experimental users quit drugs.
As R 0 is more sensitive to changes in r 1 than r 3 , it appears practical to focus efforts on the reduction of r 1 .This sensitivity analysis tells us that prevention is better than treatment.Efforts to increase prevention are more effective in controlling the transition of individuals from nonuser to experimental users than efforts to increase the number of experimental users to quit drugs, that is, the transition of individuals from experimental to nonusers.Our mathematical analysis suggest that the spread of illicit drug use should be controlled through educational campaigns at all social levels to reduce the value of r 1 .These results are provided with a view to inform and assist policy makers in targeting prevention for maximum effectiveness.

Positively invariant
The number of individuals in each category, that is N (t), E(t), R(t) and A(t) must all be non-negative for all t > 0. Thus, model (1) will be analyzed in the following feasible region We demonstrate that the closed set Ω is positively invariant by the following equations Hence, due to lemma 2 in [14] , the closed set Ω is positively invariant for model (1).In other words, every trajectory starting in Ω stays in Ω for all t > 0.

The existence of equilibria
In this section, we investigate the existence of equilibria of system (1).We write the system of equations in (1) as We solve the right hand side of ( 7) by equating it to zero for the critical points using Maple software.Three valid critical points are obtained: (i) P 0 = (1, 0, 0, 0), the drug free equilibrium, DFE (the absence of all drug users), (ii) P * = (N * , E * , R * , A * ) where, P * is known as the recreational-addict free (RAF) equilibrium because the population consists only of nonusers and experimental users.
(iii) P * * = (N * * , E * * , R * * , A * * ).This third critical point is quite complicated and we will consider that for future work.The third equilibrium consists of all the four categories of drug users and is therefore known as the drug endemic equilibrium, DEE.In the present work, the DEE is not the current scope of this study but, we have proved it numerically in Section 4.
In the sections that follows, we conduct the stability analysis of both the drug free and the recreational-addict free equilibria.

Stability analysis of the drug free equilibrium
In this section, we investigate the local and global stability of the drug free equilibrium and we analyse the effect of the basic reproduction number, R 0 on the drug free equilibrium.

Local stability
We study the local stability of the DFE by evaluating the Jacobian matrix, J, of the system (7) which is given as follows: Evaluating the Jacobian matrix at the DFE, P 0 = (1, 0, 0, 0) gives The stability of the drug free equilibrium point is analysed by studying the behaviour of J(P 0 ).Solving J(P 0 ), we obtain the following eigenvalues Clearly, λ 1 , λ 3 and λ 4 are negative and if R 0 < 1, λ 2 is also negative.Our study reveals that under the condition R 0 < 1 , the DFE is locally asymptotically stable as the eigenvalues have all negative real parts.That is, we will reach a situation of almost zero drug consumption in the population.Considering R 0 > 1, we note that the DFE is no longer stable as some eigenvalues are negative while others are positive giving rise to an unstable saddle point.Similar conclusions have been reached in [7].
We use the Theorem 4.1 in [4].
Proof.By setting R 0 = 1, we obtain the following bifurcation parameter Evaluating the Jacobian matrix at the DFE P 0 = (1, 0, 0, 0) when r 1 = r * 1 gives The eigenvalues of the matrix J(P 0 , r * 1 ) are Since there is the presence of a simple zero eigenvalue, ξ 2 = 0 and the other three are all real and negative, the DFE, P 0 , is a non-hyperbolic equilibrium.Hence, we can use the center manifold theory to analyse the dynamics of system (7).

It follows that
Solving (11), we get We next calculate the left eigenvector corresponding to the zero eigenvalue.
Solving (12) gives The coefficients of a and b for system (7) are defined as follows: (P 0 , r * 1 ).
Considering only the non-zero components of the left eigenvector υ, The rest of the second derivatives in the formula for a and b are all zero.
Clearly, we see that the coefficient of a is negative and that of b is positive.Thus, we conclude that our model demonstrates a forward bifurcation at R 0 = 1, ie , when r 1 = β + r 3 = r * 1 .

Global stability
In this section, we discuss the global asymptotically stability of the drug free equilibrium.
Proof.Let v = (N ) and w = (E, R, A), with v ∈ R denoting the proportion of nonusers and w ∈ R 3 denoting the proportion of experimental, recreational and addicts.Let P 0 = (v 0 , 0) where v 0 = 1, represents the drug free equilibrium of our model.We write system (7) as When N = N 0 = 1, P (v, 0) = 0 and dv/dt = P (v, 0) = β − β.v.Now, as t → ∞, v → v 0 = 1, thus dv dt = P (v, 0) = 0. Hence, v 0 is globally asymptotically stable.Now, where Thus, We note that Q(v, w) ≥ 0 and since the off diagonals of A are non-negative, A is clearly a Metzler matrix, and thus meeting all the required conditions.Hence, we conclude that the drug free equilibrium point P 0 is globally asymptotically stable when R 0 < 1.

Stability analysis of the recreational-addict free equilibrium
In this section, we analyse the local stability of the RAF equilibrium.

Local stability
Lemma 3.3.If the determinant of the Jacobian matrix is positive and its trace is negative, then these two criteria are enough to prove the local stability of the recreational-addict free equilibrium.
Remark 3.5.A sufficient condition for the RAF equilibrium to be locally asymptotically stable is Proof.Evaluating the Jacobian matrix at P * , we get The trace and determinant of J(P * ) is obtained using Maple software.

T race[J(P
We note that when Hence, these two conditions are adequate to prove that the RAF equilibrium is locally asymptotically stable in Ω when R 0 > 1.

Analysis of R 0 and µ on the two equilibria
In this section, we analyse the effect of both R 0 and µ on the two equilibra of model (7).
We note that when R 0 < 1, the drug free equilibrium (DFE) is stable.As R 0 > 1, the DFE is no longer stable and another equilibrium emerges, the drug endemic equilibrium which is also known as the recreational-addict free equilibrium.Considering R 0 > 1 and µ < 1, the drug free equilibrium is unstable, while the recreational-addict free equilibrium is asymptotically stable.In other words, the population is recreational-addict free.Further, when R 0 > 1 and µ > 1, both the drug free equilibrium and the recreational-addict free equilibrium are unstable and the population consists of all the four categories.More details of this section can be found in [7].Hence, it is observed that to obtain a drug free population, R 0 should be less than one, else persistence will occur.

Extinction of drug consumption
The results from sensitivity analysis suggest that educational and prevention programs can help to fight against the prevalence of illicit drug consumption.We note that the efficiency of a campaign of prevention relies on the fact that we must adjust r 1 and r 3 appropriately in order to reduce the reproduction number, R 0 .That is, we must reduce r 1 and increase r 3 simultaneously so that R 0 can be reduced accordingly.
obtain a population free of drug users. Proof.
Dividing the numerator and denominator of the RHS of ( 15) by (β + r 3 ), we get Simplifying ( 16), we get We note that in an endemic equilibrium, an effective campaign of prevention will amount to one in which more emphasis must be laid on the transition from N → E rather than taking E → N .We thus confirm the saying "prevention is better than cure".

Persistence of drug consumption
Theorem 3.7.If R 0 > 1, then model ( 7) is uniformly persistent, that is, there exists a constant δ > 0 such that where δ > 0 is independent of the initial data in Ω.
Proof.We use the persistence theory of dynamical systems in [15] (Section 1.3 in Chapter 1) to prove the theorem.
We define the sets We require to show that system (7) is uniformly persistent with respect to (Ω, Ω 0 ).From Section 3.3, we see that both Ω and Ω 0 are positively invariant.

Numerical Experiment
In this section, we use Maple software to carry out numerical simulations in order to verify the analytic results of the deterministic N ERA model.We use the data obtained from [1] to perform our numerical experiments.With slight modification, we consider the parameter values given in Table 2 (page 270 in [1]) to evaluate the consumption of Marijuana in Colorado.
The phase portrait diagrams below show the trajectories for the extinction of drug consumption among the population using the set of parameter values given by When R 0 > 1, we see that consumption of drugs will automatically persist in the population.Clearly, the phase portrait diagrams in Figure 3 demonstrate that as time t increases, all the drug classes of the N ERA model approach the drug endemic equilibrium point, P * * = (N * * , E * * , R * * , A * * ).
We next illustrate the effect of the basic reproduction number, R 0 on the four categories of drug classes.
Using the same set of parameter values in (25), the figure below shows the evolution of N (t), E(t), R(t) and A(t) at R 0 < 1 .From Figure 4, we see that as time t increases, the number of experimental, recreational and addict users decreases to zero while the nonuser sums to 1.In other words, we obtained a drug free population thereby confirming that the drug free equilibrium is stable when R 0 < 1.
Figure 5 shows persistence of drug consumption in the population.We observe that all the drug classes of the N ERA model arrive at the drug endemic equilibrium point, P * * = (N * * , E * * , R * * , A * * ) when R 0 > 1.
From Theorem 3.6, we note that an efficient prevention campaign is one in which r 1 is reduced and r 3 is increased simultaneously so as to lower R 0 .Thus, we use the data in Table 2 and we adjust r 1 and r 3 accordingly.Figures 6 and 7 support Theorem 3.6 that prove the effectiveness of a prevention campaign where drug consumption will be wiped out from a given population.
Figure 6 clearly indicates that when we reach 50 time steps, we have the same number of recreational users and addicts.As we march in time, we see that the proportion of recreational and addict users decreases to zero while the nonusers and the experimental users reach a constant value, that is, they approach the recreational-addict free equilibrium point.We notice that a slight decrease and an increase in the value of r 1 and r 3 respectively, result in thwarting the prevalence of drug consumption.
Figure 7 illustrates the effect of increasing the parameter r 1 further and decreasing the parameter r 3 .We observe that the proportion of experimental users, recreational users and addicts drops to zero.Hence, this is an indication that as more experimental users quit, there will be a decrease in the prevalence of drug use since nonusers become less likely to join the experimental category.In other words, we note that reducing the parameter r 1 and increasing the parameter r 3 will lower the drug epidemic.that when R 0 < 1, the drug free equilibrium P 0 , is stable and P 0 loses its stability, giving rise to another stable equilibrium, P * , when R 0 > 1 at the bifurcation value r * 1 = β + r 3 .
In order to reach a situation of almost zero drug consumption in a given society, the influence rate, r 1 , of experimental users on nonusers should be less than the bifurcation value r * 1 , so that E * and P * will be close to zero, otherwise, persistence will occur.

Conclusion
In this work, we considered a deterministic NERA model which studies the dynamics of illicit drug consumption in a given population.Two reproduction numbers, R 0 and µ have been calculated and were used to analyse the stability of the system.Sensitivity analysis has been conducted on the parameters of R 0 and the calculated sensitivity indices were used to analyse the model.We have demonstrated the local stability of both the drug free and the recreational-addict free equilibrium.Global stability of the drug free equilibrium has also been proved.Moreover, the model was shown to exhibit a forward bifurcation at R 0 = 1.We also showed the prevalence of persistence in drug consumption when the basic reproduction number, R 0 > 1. Conditions for the extinction of drug consumption were also established.Finally, numerical experiments were conducted to confirm our analytical results.The results from sensitivity analysis suggest that effective campaigns of prevention can help to control the spread of illicit drug use.Thus, our results show that the NERA model can be used as a policy control mechanism in tackling illicit drug consumption in a given population.
The presence of a randomness in the influence that the experimental users exert on the nonusers [1] and in the effect of dopamine in the transition of recreational to addict category suggests that a stochastic system taking those two features of drug consumption will be a better representation of the reality.Formulation of such a stochastic NERA model is currently underway.

Figure 1 .
Figure 1.Schematic representation of N ERA model

{r 1 Figure 2 .
Figure 2 demonstrates that when R 0 < 1, we obtained a drug free population.The three plots N (t) − E(t),

Table 1 .
Interpretation of the parameters in N ERA model

Table 2 .
Parameter values for the consumption of Marijuana in Colorado