key: cord-0781171-jgmg66sr
authors: Tang, Haoxiang; Li, Mingtao; Yan, Xiangyu; Lu, Zuhong; Jia, Zhongwei
title: Modeling the Dynamics of Drug Spreading in China
date: 2021-01-02
journal: Int J Environ Res Public Health
DOI: 10.3390/ijerph18010288
sha: a5bed2f2767af0c19b9645f3f502fc86f79e6eda
doc_id: 781171
cord_uid: jgmg66sr

Drug abuse remains one of the major public health issues at the global level. In this article, we propose a drug epidemic model with a complete addiction–rehabilitation–recovery process, which allows the initiation of new users under the influence of drug addicts undergoing treatment and hidden drug addicts. We first conduct qualitative analyses of the dynamical behaviors of the model, including the existence and positivity of the solutions, the basic reproduction number, global asymptotic stabilities of both the drug-free and the drug-persistent equilibria, as well as sensitivity analysis. Then we use the model to predict the drug epidemic in China during 2020–2030. Finally, we numerically simulate the potential impact of intervention strategies on different drug users. The results show that the drug epidemic will decrease significantly during 2020−2030, and the most effective intervention strategy to eliminate drug epidemics is to strengthen the investigation and rehabilitation admission of hidden drug users.

The phenomenon of drug abuse, which involves the consumption of illicit drugs and nonmedical use of prescription drugs, has become one of the global health issues threatening the safety and sustainability of human society in the 21st century. According to 2019 World Drug Report released by the United Nations Office on Drugs and Crime (UNODC), approximately 271 million people, which constituted 5.5% of the global population aged 15−64, had used drugs in 2016 [1] . From a historical perspective, the world has witnessed a 30% increase in the drug-using population ever since 2009 [1, 2] . In terms of drug type, opioids remained the most lethal group, which resulted in around 66% of overdose-related deaths worldwide in 2017 [1] . The level of manufacture and trafficking of conventional drugs such as cocaine and cannabis remained high, and that of synthetic drugs such as methamphetamine and 3,4-methylenedioxy-n-methylamphetamine (MDMA) even soared in recent years. Some 35 million people suffered from drug use disorders and required treatment service around the world, and the death toll attributed to drug use totaled 585,000 in 2017 [1] .

Despite its amelioration for two consecutive years, the drug situation in China remains a serious issue. According to the 2019 Report of Drug Situation in China, by the end of 2019, the number of drug users in China totaled 2.14 million (excluding the dead, those who went abroad, or those who remained abstinent for at least 3 years), which accounted for 0.16% of the total population [3] . In the same year, law enforcement agencies of China settled 83,000 drug-related crimes and seized 65.1 tons of drugs. Globally, around 43%

To study the dynamics of drug abuse in China, we formulate this drug epidemic model, in which the total population is divided into five mutually exclusive compartments: the susceptible population (S), light drug users (I 1 ), drug addicts undergoing treatment (I 2 ), hidden drug addicts (I 3 ), and recovered individuals (R).

The population of interest is civilians aged between 15 and 64 years, while those outside of this age bracket are supposedly either too young to get in touch with drugs or old enough to mature out [18, 58] . The susceptible compartment, which is defined as civilians with no history of drug use, receives a constant population inflow at an annual rate λ. This inflow takes into account young people reaching 15 years old as well as immigrants, and it is the only approach of population replenishment from outside the system. Light drug users (I 1 ) are defined as individuals who tried drugs but have not reached addiction level, the criteria of which could be referred to Methods for the Identification of Drug Addiction issued by the Ministry of Public Security of China [59] . Drug addicts undergoing treatment (I 2 ) are defined as drug users who have reached addiction level and are currently receiving detoxification rehabilitation. The drug rehabilitation system in China comprises compulsory-isolated detoxification centers, voluntary detoxification facilities, community detoxification and methadone maintenance treatment clinics (MMTs) [60] [61] [62] [63] , among which compulsory-isolated detoxification centers adopt in-patient bases and require isolation for two years, and other forms of detoxification facilities allow certain degrees of free movement. As a consequence, the I 2 compartment is actually a mixture of patients receiving various forms of detoxification rehabilitation. Hidden drug addicts (I 3 ) are defined as drug users who have reached addiction level and remained hidden to the law enforcement system, and the recovered compartment (R) is defined as individuals who have reached abstinence through detoxification treatment.

We assume homogeneous mixing and that the spread of drug-using behavior in the public can be modeled similar to an infectious disease. Susceptible individuals may initiate drug-using behavior and be converted to light drug users following interactions with drug addicts undergoing treatment (I 2 ) or hidden drug addicts (I 3 ), in which process we adopt a bilinear incidence rate. The compartment of drug addicts undergoing treatment (I 2 ), which is a mixture of patients receiving various forms of rehabilitation, when regarded as a whole, manifests a lower mobility and an effective contact rate when compared with their hidden counterparts. Due to the illicit nature of drug-using behaviors in China, hidden drug addicts (I 3 ) would try to evade law enforcement departments and are able to impose a larger impact on the transmission of drug addiction. We assume no fast progressor, which means a susceptible person must first become a light drug user before entering the compartments of addicts. While progressing to a drug addict, an individual would either be discovered and admitted to treatment facilities or remain hidden, corresponding to transfer rates k 1 and k 2 , respectively. Hidden drug addicts (I 3 ) could also be transferred to the I 2 compartment at admission rate α, and those who finish rehabilitation will leave I 2 and enter the recovered group R at recovery rate r. Without direct scientific proof of self-abstinence, we assume no other recovery approach without treatment in this model. Each compartment bears a natural death rate µ, and hidden drug addicts (I 3 ) suffer from an additional death rate µ d resulting from drug abuse [39, 44] . The definitions of all compartment variables and parameters are listed in Table 1 . 

The susceptible individuals at time t I 1 (t)

Light drug users at time t I 2 (t) Drug addicts undergoing treatment at time t I 3 (t)

Hidden drug addicts at time t R(t) Recovered individuals at time t λ Inflow rate into the susceptible individuals µ Natural death rate µ d Additional death rate resulting from drug abuse β 1 Effective contact rate between drug addicts in treatment and susceptibles β 2 Effective contact rate between hidden drug addicts and susceptibles k 1 Progression rate from light drug users to drug addicts in treatment k 2 Progression rate from light drug users to hidden drug addicts α Discovery and admission rate from hidden addicts to addicts in treatment r

Recovery rate of drug addicts undergoing treatment

The total population at time t is given by N(t) = S(t) + I 1 (t) + I 2 (t) + I 3 (t) + R(t). Based on the model assumptions above, we can obtain the model diagram ( Figure 1 ) and the following set of nonlinear differential equations: 

Lemma 1.  with the topology of uniform convergence. System (1) with non-negative initial conditions can be considered as the following initial-value problem:

Let Ω be an open subset in C ×  

Lemma 1. The model system admits a unique solution for non-negative initial conditions

Proof Suppose C [a, b], R 5 is the Banach space of continuous functions mapping the interval [a, b] into R 5 with the topology of uniform convergence. System (1) with non-negative initial conditions can be considered as the following initial-value problem:

. Let Ω be an open subset in R × C and f : Ω → R 5 is the mapping function of System (1). It is obvious that f (t, x) is continuous and Lipschitzian in x in each compact set in Ω. Since the initial point (t 0 , x 0 ) ∈ Ω, according to Cauchy-Lipschitz Theorem, there is a unique solution of System (1) passing through (t 0 , x 0 ) [64] .

The feasible region is an interval where System (1) will be analyzed, and it should be forward invariant for biological reasons. Thus, the following lemma will state the positive invariance and attractiveness of the system's feasible region.

Lemma 2. The feasible region of the model system is defined by

Γ is a positively invariant set, and it is attracting with regard to system (1) for all t > 0.

Proof Adding all the equations of System (1), we obtain

Solving this differential inequality, we have 0 ≤ N(t) ≤ λ µ + N(0) − λ µ e −µt , where N(0) represents the sum of the initial values of all variables. Taking the limit as t → ∞ , we have that 0 ≤ N ≤ λ µ . Thus, the state variables will remain biologically meaningful in the feasible region Γ for all positive initial conditions, and Γ is positively invariant and attractive with respect to system (1). Hence, System (1) is well-posed mathematically and biologically in Γ, and it would suffice to study the dynamics of the system in Γ.

According to their biological meanings, all state variables and parameters are supposed to remain positive during the modeling period. In the following lemma, we will show that with positive initial conditions, each variable would remain non-negative for all t > 0.

Lemma 3. Given the initial conditions S(0) ≥ 0, I 1 (0) ≥ 0, I 2 (0) ≥ 0, I 3 (0) ≥ 0, R(0) ≥ 0, the solutions S(t), I 1 (t), I 2 (t), I 3 (t) and R(t) will remain non-negative for all t > 0.

Proof Assume that t = sup{t > 0 :

Thus, t > 0, and from the first equation of System (1), we have dS

From the second equation of System (1), we obtain dI 1 dt ≥ −(k 1 + k 2 + µ)I 1 .

Solving the differential inequality, we have:

Similarly, it can be shown that I 2 (t) > 0, I 3 (t) > 0 and R(t) > 0 for all t > 0.

In this model, the basic reproduction number is defined as the number of secondary drug users converted from susceptible individuals by a single drug addict (either hidden or undergoing treatment) introduced into a totally susceptible population during his entire addiction period [8, 54] . The basic reproduction number plays a vital role in determination of the existence of the drug epidemic, as well as the analysis of dynamics of the model. In this section, we will obtain the formula of the basic reproduction number with the next-generation matrix method [65] .

First, we rearrange the equations of System (1) according to the order of I 1 , I 2 , I 3 , S, R and let x = (I 1 , I 2 , I 3 , S, R) be the solution to the system; then, System (1) can be rewritten as:

where F represents the new initiate terms, and V corresponds to the interior transfer terms. The corresponding linearized matrices of F and V evaluated at the drug-free equilibrium

The next-generation matrix is given by

Thus, the basic reproduction number R 0 is defined as the spectral radius of the nextgeneration matrix, which is the largest absolute value of its eigenvalues:

The epidemiological meaning of this formula can be understood in the following way. The first term of the expression above represents the role of drug addicts undergoing treatment, i.e., when a drug addict in treatment is introduced into a wholly susceptible population with a size of λ µ , the number of light drug users converted from susceptible individuals under his/her influence in unit time is β 1 λ µ , the proportion of which who progress to I 2 compartment is k 1 k 1 +k 2 +µ , whose average addiction period is 1 r+µ ; hence, the number of I 2 directly generated during this period is k 1 β 1 λ µ(k 1 +k 2 +µ)(r+µ) . On the other hand, Int. J. Environ. Res. Public Health 2021, 18, 288 7 of 25 the proportion of I 1 who progress to I 3 is k 2 k 1 +k 2 +µ , and the proportion of I 3 who progress to I 2 is α α+µ+µ d ; hence, the number of I 2 generated through I 3 during this period is obtained as αk 2 β 1 λ µ(k 1 +k 2 +µ)(r+µ)(α+µ+µ d ) . Adding these two parts gives the first term of the formula. The latter term can be explained in a similar way, which represents the role of hidden drug addicts. When a hidden drug addict is introduced into this wholly susceptible population, the number of light drug users converted from susceptible individuals under his/her influence in unit time is β 2 λ µ , the proportion of which who progress to I 2 compartment is k 2 k 1 +k 2 +µ , and the average addiction period of I 3 is 1 α+µ+µ d ; hence, the number of I 2 generated during this period is obtained as

To sum up, the formula of basic reproduction number represents the overlay of the influence of drug addicts undergoing treatment and hidden drug addicts.

In order to investigate the existence of the drug-free equilibrium and the drugpersistent equilibrium, we set the left side of the equations of System (1) at zero:

Through direct calculations, we obtain:

1. When I 1 = 0, it is easy to acquire I 2 = I 3 = R = 0, S = λ µ . At this point, N = S = λ µ , and dN dt = λ − µN − µ d I 3 = 0 is also satisfied. We refer to this point as the drug-free equilibrium E 0 = λ µ , 0, 0, 0, 0 .

When

, and it can be easily verified that dN dt = λ − µN − µ d I 3 = 0 also holds. We refer to this point as the unique drugpersistent equilibrium E * = S * , I * 1 , I * 2 , I * 3 , R * . For biological reasons, it requires that S * < λ µ , which corresponds to the condition that R 0 > 1.

Theorem 1. The drug-free equilibrium E 0 is globally asymptotically stable when R 0 ≤ 1.

Proof See Appendix A for a detailed proof of Theorem 1.

Theorem 2. The drug-persistent equilibrium E * is globally asymptotically stable when R 0 > 1.

Proof See Appendix B for a detailed proof of Theorem 2.

Due to the critical role of the basic reproduction number R 0 in determination of the persistence of the drug epidemic, it is of vital importance to identify the most effective approach to bring R 0 down to below one. In this section, we calculate the normalized forward sensitivity index (NFSI) of R 0 with respect to each parameter following Arriola and Hyman [66] . Considering their biological and epidemiological meanings, the changes of some parameters are either impractical or unethical (demographic parameters or inherent progression rates, etc.). Hence, we have identified four parameters of interest and calculated their NFSIs:

The signs of the numerators of A k 1 and A α depend on the final values of the parameters, and it is easy to prove that all NFSIs are lower than one. According to their biological meanings, we have β 2 β 1 for most cases; hence, it seems obvious that A β 2 > A β 1 . That is to say, R 0 is more sensitive to changes in β 2 than those in β 1 . The relative sizes of A β 2 , A k 1 , and A α are yet to be determined, depending on the exact value of each parameter.

With the aim of illustrating the sensitivity of the basic reproduction number R 0 with respect to the parameters of interest, we conduct a series of numerical simulations based on artificial training parameters. The parameter ranges are chosen to accommodate reasonable biological meanings, and extra attention is paid when the basic reproduction number R 0 crosses one, which acts as the threshold for the persistence of the drug epidemic.

(1) Comparison between β 1 and β 2 Fix the parameters at λ = 100, µ = 0.007, µ d = 0.025, k 1 = 0.05, k 2 = 0.2, α = 0.05, and r = 0.5. Set β 1 ∈ [1e − 7, 1e − 5] and β 2 ∈ [1e − 7, 1e − 5], and draw the phase plane as well as the contour plot of R 0 with respect to β 1 and β 2 ( Figure 2 ). According to the phase plane (Figure 2a ), the value of R 0 decreases sharply as β 2 decreases, but the change of β 1 has a significantly smaller impact on R 0 . The contour plot ( Figure 2b ) demonstrates similar results, which indicates that R 0 is far more sensitive to β 2 than to β 1 in this case.

(2) Comparison between β 2 and k 1 Fix the parameters at λ = 100, µ = 0.007, µ d = 0.025, k 2 = 0.2, α = 0.05, r = 0.5, and β 1 = 1e − 6. Set β 2 ∈ [1e − 7, 1e − 5] and k 1 ∈ [0.005, 0.5], and draw the phase plane as well as the contour plot of R 0 with respect to β 2 and k 1 (Figure 3 ). Both the phase plane ( Figure 3a ) and the contour plot ( Figure 3b ) indicate that β 2 is more effective than k 1 in controlling R 0 . In this case, increasing k 1 also lowers R 0 to some extent; note that R 0 is negatively correlated with k 1 .

(3) Comparison between β 2 and α Fix the parameters at λ = 100,

and α ∈ [0.005, 0.5], and draw the phase plane as well as the contour plot of R 0 with respect to β 2 and α ( Figure 4 ). As demonstrated in the phase plane ( Figure 4a ) and the contour plot (Figure 4b ), the relative effectiveness of β 2 and α in controlling R 0 seems debatable. In the top-left corner of the contour plot, which corresponds to lower values of α and higher values of β 2 , α turns out more efficient in adjusting the basic reproduction number R 0 .

phase plane (Figure 2a) , the value of 0 R decreases sharply as 2 β decreases, but the change of 1 β has a significantly smaller impact on 0 R . The contour plot ( Figure 2b ) demonstrates similar results, which indicates that 0 R is far more sensitive to 2 β than to 1 β in this case.

. Figure 2 . Phase plane (a) and contour plot (b) of the basic reproduction number 0 R with respect to 1 β (effective contact rate between drug addicts undergoing treatment and the susceptible) and 2 β (effective contact rate between hidden drug addicts and the susceptible). , and draw the phase plane as well as the contour plot of 0 R with respect to 2 β and 1 k (Figure 3 ). Both the phase plane ( Figure 3a ) and the contour plot ( Figure 3b ) indicate that 2 β is more effective than 1 k in controlling 0 R . In this case, increasing 1 k also lowers 0 R to some extent; note that 0 R is negatively correlated with 1 k .

Phase plane (a) and contour plot (b) of the basic reproduction number R 0 with respect to β 1 (effective contact rate between drug addicts undergoing treatment and the susceptible) and β 2 (effective contact rate between hidden drug addicts and the susceptible).

. Figure 2 . Phase plane (a) and contour plot (b) of the basic reproduction number 0 R with respect to 1 β (effective contact rate between drug addicts undergoing treatment and the susceptible) and 2 β (effective contact rate between hidden drug addicts and the susceptible). β is more effective than 1 k in controlling 0 R . In this case, increasing 1 k also lowers 0 R to some extent; note that 0 R is negatively correlated with 1 k . Figure 3 . Phase plane (a) and contour plot (b) of the basic reproduction number R 0 with respect to β 2 (effective contact rate between hidden drug addicts and the susceptible) and k 1 (transfer rate from light drug users to drug addicts undergoing treatment).

With the aid of numerical simulation tools, we are also able to illustrate the global stability of the model equilibria. For this purpose, we fix the parameters and obtain five distinct solutions for five different sets of initial values.

(1) Drug-free equilibrium Fix the parameters at λ = 400, µ = 0.007, µ d = 0.025, k 1 = 0.05, k 2 = 0.2, r = 0.5, α = 0.05, β 1 = 1e − 7, and β 2 = 1e − 6. The solutions of five sets of initial values were demonstrated in Figure 5 , including 3D trajectories in the I 1 − I 2 − I 3 space and long-term time-series plot of each compartment. In this case, R 0 = 0.5498 and the model equilibrium manifests as a drug-free equilibrium. The statement is supported by Figure 5 , where the number of susceptible individuals stabilizes at somewhere above zero, and all other compartments fade away with time. . Phase plane (a) and contour plot (b) of the basic reproduction number 0 R with respect to 2 β (effective contact rate between hidden drug addicts and the susceptible) and α (transfer rate from hidden drug addicts to drug addicts undergoing treatment). (2) Drug-persistent equilibrium Fix the parameters at and the model equilibrium manifests as a drug-persistent equilibrium. This statement is supported by Figure 6 , where all compartments persist and stabilize somewhere above zero. (2) Drug-persistent equilibrium Fix the parameters at λ = 400, µ = 0.007, µ d = 0.025, k 1 = 0.05, k 2 = 0.2, r = 0.5, α = 0.05, β 1 = 1e − 6, and β 2 = 1e − 5. The solutions of five sets of initial values were demonstrated in Figure 6 . In this case, R 0 = 5.4985 and the model equilibrium manifests as a drug-persistent equilibrium. This statement is supported by Figure 6 , where all compartments persist and stabilize somewhere above zero.

(2) Drug-persistent equilibrium Fix the parameters at Figure 6 . In this case, 0 5.4985 R = and the model equilibrium manifests as a drug-persistent equilibrium. This statement is supported by Figure 6 , where all compartments persist and stabilize somewhere above zero. 

Based on the qualitative results of dynamical behaviors acquired above, we now apply the model to the drug epidemic in real world. Since our model possess a complete 

Based on the qualitative results of dynamical behaviors acquired above, we now apply the model to the drug epidemic in real world. Since our model possess a complete addiction-rehabilitation-recovery process, and its basic assumptions are more compatible with the current situation in mainland China, we fit our model to the historical data of drug users in China and aim to predict the scale and trend of drug users in the near future as well as analyze the potential effect of various intervention strategies.

The most comprehensive and publicly available official sources are the Annual Report of Drug Situation in China and Annual Report on Drug Control in China released by the National Narcotics Control Committee (NNCC) [67] . Since 2015, the reports have changed their scopes of statistics and provided the numbers of all existing drug users (excluding those who had died, emigrated, or remained abstinent for three or more years) instead of the numbers of cumulative totals previously reported [3, 68] . According to the Law of Drug Control and the Drug Rehabilitation Ordinance, all drug users identified will receive their corresponding type of rehabilitation treatment; thus, we can safely conclude that the numbers of all existing drug users released by the NNCC correspond to the compartment of drug addicts undergoing treatment (I 2 ) in our model. The numbers of former drug users who had remained abstinent for three or more years reported by the NNCC approximately correspond to the compartment of individuals who have reached abstinence through detoxification treatment (R). Without a direct data source, the initial value of hidden drug addicts (I 3 ) can be ascertained through the initial value of I 2 compartment and the explicit-to-implicit ratio of 1:4 reported by the NNCC and the existing literature [62, [68] [69] [70] . Likewise, the initial value of light drug users is estimated based on dynamical relationships with neighboring compartments. Finally, the initial value of the susceptible individuals (S) is obtained through subtracting all other compartments from the total population aged 15−64 recorded by the National Bureau of Statistics [71] . The historical numbers of population aged 15−64, existing drug users, and former drug users who had remained abstinent for three or more years are listed in Table 2 . 

We obtain demographic parameters mainly from official data released by the authorities, among which the natural death rate µ was acquired from the website of The National Bureau of Statistics, and inflow rate into the susceptible compartment λ was ascertained through the equation that Net population growth = Inflow-Death [71]. Other parameters with explicit sources include the additional death rate of hidden drug addicts µ d and the progression rate from light drug users to hidden drug addicts k 2 [23, 25, 39, 44, 72] . Without direct data, the recovery rate was ascertained according to the rehabilitation durations of the detoxification facilities [60] . Other parameters involve implicit processes that are hard to study directly and will need to be ascertained through fitting procedures. The potential ranges of the parameters were chosen based on dynamical relationships between neighboring compartments and parameters with similar roles in other modeling studies as well [23, 25, 39, 44, 60, 72] . For instance, with around 1 billion susceptible individuals and 10 million hidden drug addicts at the starting point, the incidence rate of new initiates under the influence of hidden drug addicts shall not exceed 1 million/year (which approximates the number of light drug users at the initial phase). According to the bilinear incidence rate adopted, we obtain that an effective contact rate of β 2 < 1e − 6 /10 thousand people*year from this inequality. We list the parameter units, value ranges, final values used, and their sources in Table 3 . 

In this section, we fit our model to historical data of drug abuse in China, where re- Figure 7 that the numbers of light drug users, drug addicts undergoing treatment, and hidden drug addicts experience decreases of 81.22%, 71.98%, and 84.69% respectively in 15 years, so long as the model dynamics remain unperturbed. The sum of all compartments experiences a mild decrease of 4.44%, which is in consistent with the tendencies of historical data. The final values of the parameters acquired through curve fitting are listed in Table 3 , and the corresponding basic reproduction number was calculated as R 0 = 0.087256, which is quite low and accords with the rapid decreases observed in the number of drug addicts undergoing treatment and hidden drug addicts. 

Despite the delightful tendency of decrease generated by model simulation in the previous section, it is still the objective of policy makers and public health workers to further shrink the scale of drug epidemic. In this section, we present hypothesized results of several potential intervention strategies through numerical simulation and compare their effects on drug control.

(1) Intervention 1: anti-drug education and propaganda Previous studies have shown that a misconception of drugs (86−90%) and peer influence (13−44%) are the most common reasons for initiating drug use among Chinese drug users [73, 74] . School-based anti-drug education or preventive propaganda through media publicity has been implemented in China in the past decades and proved itself as an effective tool against drug epidemic [75] . Strengthening anti-drug propaganda could correct misunderstandings of drugs among the susceptible individuals and lower the risk of first exposure to illicit drugs. We assume that the effective contact rates 1 β and 2 β are inversely proportional to the intensity of anti-drug education and propaganda, and that doubling the frequency and budgets of such activities could lower 1 β and 2 β by 50%. Suppose this strategy starts to take effect since the end of 2020, and plot the new curves of the model solution in parallel with the original ones ( Figure 8) . The results show that the implementation of this intervention will be able to lower the number of light drug users (−49.27%) and drug addicts undergoing treatment (−12.79%) by 2030, but its effect on hidden addicts (−9.58%) is relatively smaller. 

Despite the delightful tendency of decrease generated by model simulation in the previous section, it is still the objective of policy makers and public health workers to further shrink the scale of drug epidemic. In this section, we present hypothesized results of several potential intervention strategies through numerical simulation and compare their effects on drug control.

(1) Intervention 1: anti-drug education and propaganda Previous studies have shown that a misconception of drugs (86−90%) and peer influence (13−44%) are the most common reasons for initiating drug use among Chinese drug users [73, 74] . School-based anti-drug education or preventive propaganda through media publicity has been implemented in China in the past decades and proved itself as an effective tool against drug epidemic [75] . Strengthening anti-drug propaganda could correct misunderstandings of drugs among the susceptible individuals and lower the risk of first exposure to illicit drugs. We assume that the effective contact rates β 1 and β 2 are inversely proportional to the intensity of anti-drug education and propaganda, and that doubling the frequency and budgets of such activities could lower β 1 and β 2 by 50%. Suppose this strategy starts to take effect since the end of 2020, and plot the new curves of the model solution in parallel with the original ones ( Figure 8) . The results show that the implementation of this intervention will be able to lower the number of light drug users (−49.27%) and drug addicts undergoing treatment (−12.79%) by 2030, but its effect on hidden addicts (−9.58%) is relatively smaller. (2) Intervention 2: moderate anti-drug propaganda Accounting for the limited resources available for anti-drug education and propaganda, as well as the huge number of susceptible individuals, we might consider an alternative approach instead of an amplification of 100% in the intensities of such activities. Hence, in this intervention strategy, the frequency and budgets of anti-drug education and propaganda are increased by 50%, which correspond to a 33% decrease in effective contact rates 1 β and 2 β . We repeat the rest of the procedures and obtain simulation results in Figure 9 . The results showed that the number of light drug users will be lowered by 33.20% compared to the original curve, and the impact on hidden drug addicts is minimal. (2) Intervention 2: moderate anti-drug propaganda Accounting for the limited resources available for anti-drug education and propaganda, as well as the huge number of susceptible individuals, we might consider an alternative approach instead of an amplification of 100% in the intensities of such activities. Hence, in this intervention strategy, the frequency and budgets of anti-drug education and propaganda are increased by 50%, which correspond to a 33% decrease in effective contact rates β 1 and β 2 . We repeat the rest of the procedures and obtain simulation results in Figure 9 . The results showed that the number of light drug users will be lowered by 33.20% compared to the original curve, and the impact on hidden drug addicts is minimal.

(3) Intervention 3: investigation and admission In contrast to preventive strategies, we now consider improving the intensity of investigation of hidden drug users and their admission into drug rehabilitation facilities. We assume that the transfer rates k 1 and α are proportional to the intensity of investigation and admission, and that doubling the frequency and budgets of such activities, as well as the number of rehabilitation facilities, could increase k 1 and α by 100%. Supposing this change of parameters takes place since the end of 2020, we obtained the simulation results shown in Figure 10 . It could be observed that the number of light drug users (−71.18%) and hidden drug addicts (−71.44%) undergoes sharp decreases compared to the original curves, and that of drug addicts in treatment I 2 , though experiencing a temporary increase, ends up lower than the original curve (−16.82%). (3) Intervention 3: investigation and admission In contrast to preventive strategies, we now consider improving the intensity of investigation of hidden drug users and their admission into drug rehabilitation facilities. We assume that the transfer rates 1 k and α are proportional to the intensity of investigation and admission, and that doubling the frequency and budgets of such activities, as well as the number of rehabilitation facilities, could increase 1 k and α by 100%. Supposing this change of parameters takes place since the end of 2020, we obtained the simulation results shown in Figure 10 . It could be observed that the number of light drug users (−71.18%) and hidden drug addicts (−71.44%) undergoes sharp decreases compared to the original curves, and that of drug addicts in treatment 2 I , though experiencing a temporary increase, ends up lower than the original curve (−16.82%). (4) Intervention 4: moderate investigation On account of the available police strengths and limited rehabilitation capacities of the detoxification facilities, we might consider an alternative approach instead of a thorough search of hidden drug users for the time being [67] . Hence, in this intervention strategy, the frequency and budgets of investigation activities, as well as the number of rehabilitation facilities are increased by 50%, corresponding to increases of 50% in k 1 and α. Repeat the rest of the procedures, and we obtain the results in Figure 11 . A considerable drop is still observed in the I 3 compartment (−46.66%), and the final number of I 2 is also smaller than its original counterpart (−4.78%).

(5) Comparison of the interventions In order to visually compare the effects of different intervention strategies on drug control, we tabulate the final numbers of all compartments in 2030 and their relative growth compared to the original curves in Table 4 . The results show the high efficiency of Interventions 3 and 4 in reducing the number of drug users, especially hidden drug addicts. For instance, though possessing proximate basic reproduction numbers, Interventions 1 and 3 generated totally diverse outcomes in that the declining percentage of hidden drug addicts in Intervention 3 (−71.44%) was around 7 times larger than that in Intervention 1 (9.58%). Similar situations are also observed for Interventions 2 and 4. (4) Intervention 4: moderate investigation On account of the available police strengths and limited rehabilitation capacities of the detoxification facilities, we might consider an alternative approach instead of a thorough search of hidden drug users for the time being [67] . Hence, in this intervention strategy, the frequency and budgets of investigation activities, as well as the number of rehabilitation facilities are increased by 50%, corresponding to increases of 50% in 1 k and α . Repeat the rest of the procedures, and we obtain the results in Figure 11 . A considerable drop is still observed in the 3 I compartment (−46.66%), and the final number of 2 I is also smaller than its original counterpart (−4.78%). Table 4 . The results show the high efficiency of Interventions 3 and 4 in reducing the number of drug users, especially hidden drug ad- Figure 11 . The effect of Intervention 4 (regional investigation, red curves) in comparison with the original drug situation (blue curves). Subplots include time-series plots of light drug users (a), drug addicts undergoing treatment (b), hidden drug addicts (c), and recovered individuals (d). 

In this article, we propose a drug epidemic model with a complete addictionrehabilitation-recovery process, which assumes the conversion of susceptible individuals into light drug users under the influence of drug addicts undergoing treatment and hidden drug addicts. Unlike many previous studies, we discard unrealistic assumptions such as self-detoxification without treatment or permanent immunity to drugs granted by "vaccines" or education [49, 57] . We have acquired qualitative results of the dynamical behaviors of the model, including the feasible region, basic reproduction number, global asymptotic stabilities of the drug-free equilibrium and the drug-persistent equilibrium, as well as sensitivity analysis realized through normalized forward sensitivity indices. Subsequently, we applied the model to the drug epidemic in China and obtained the numerical simulation results via curve fitting and projections. The results show significant decreases in the numbers of all groups of drug users, including light drug users, drug addicts undergoing treatment, and hidden drug addicts. Should the model dynamics remain undisturbed, the predicted drug shrink in the following decade will be a positive signal to the accumulative anti-drug efforts by the Chinese government and public health workers, and it is in accordance with the 2030 Agenda for Sustainable Development by the United Nations [1].

One of the most interesting results could be observed in Section 5.4, where Interventions 3 and 4, which correspond to increasing transfer rates k 1 and α, turned out to be more efficient in reducing the numbers of existing drug users (including light drug users, drug addicts undergoing treatment, and hidden drug addicts) as well as the basic reproduction number than Interventions 1 and 2 (corresponding to lowering the effective contact rates β 1 and β 2 ). This conclusion is in accordance with the sensitivity analyses in Section 4, in that the parameters β 2 = 3.86e − 7, α = 0.124 acquired through curve fitting correspond to the bottom-left corner of Figure 4b , where the basic reproduction number R 0 is more sensitive to changes in α than in β 2 . Likewise, it can be shown that R 0 is far more sensitive to β 2 than to k 1 or β 1 . As a consequence, it would not be surprising that a combination of α and k 1 (Interventions 3 and 4) is slightly more efficient in reducing R 0 than the combination of β 1 and β 2 (Interventions 1, 2). The conclusions above seem to contradict one of the most frequently stated conclusions that "prevention is better than cure" made by several previous studies [18, 20, 21] at first glimpse, but they are actually different aspects of the drug epidemic. Basic reproduction number R 0 focuses on the capability of new addicts to convert susceptible individuals into new initiates, whose effect on the existing drug addicts are realized through complex dynamical processes. On the contrary, transfer rates α and k 1 directly act on the existing number of hidden drug addicts and light drug users, and they proved to be efficient in controlling the scale of the drug epidemic. A direct proof of this phenomenon is the slight effect of Intervention 1 or 2 on the number of hidden drug addicts compared to that of Intervention 3 or 4, despite their similar basic reproduction numbers.

Another fact worth noticing is the tiny value of the basic reproduction number (R 0 = 0.087256 in the baseline). The calculation of this value was based on the formulation obtained in Section 3.2 and the parameter values acquired through curve fitting. Since we used the historical data of 2015−2018 to ascertain the parameter values and the data of 2019 to verify the projected trend, the result showed an overall well goodness of fit, and we have reasons to believe that this value of R 0 is an authentic manifestation of the drug situation in China for the modeling period. Given the illicit nature of drug-taking behaviors in China and the low reported number of existing drug users in official data, the scale of drug epidemic is by no means comparable to those of infectious diseases, and it seems reasonable that the R 0 of drug-using behaviors is smaller than those of infectious diseases by an order of magnitude [11] . Based on this R 0 and the already-decreasing number of drug users observed since 2018, it would not be difficult to understand the larger impact of Intervention 3 or 4 compared with Intervention 1 or 2 on the already-shrinking drug epidemic. In addition, it should also be noted that the number of drug addicts undergoing treatment in Interventions 3 and 4 even ended up lower than the original curves, and we owe this fact to the relatively higher value of β 2 compared with β 1 , as well as the relevant dynamical processes. We believe the discussion above partially explains the observed results, and in the meantime, it offers new insights into formulations of anti-drug strategies and policies. The basic reproduction number R 0 is still an important threshold value for determination of the persistence of the drug epidemic, but cautions should be taken when choosing the appropriate strategy to further eliminate drug spread.

Our study possesses the following novelties: (1) a complete addiction-rehabilitationrecovery process, without unrealistic assumptions such as self-detoxification or permanent immunity; (2) application to the historical data of drug users in China and projection to the future; (3) novel insights in discussions of intervention strategies to accelerate the reduction of existing drug users. Similar to all existing modeling research studies, we acknowledge that our present study bears certain limitations. Above all, the assumption of permanent abstinence and absence of a relapse process is adopted for ease of analysis and application of the model, which is a major contradiction with the real situation around the world [76] . Secondly, the scarcity of historical data occurred due to former scopes of statistics adopted by the NNCC, which only provided accumulative numbers of drug users, and there was no other direct source available to secure the data of interest. Drug-use patterns and demographic characteristics may differ greatly among different subgroups divided according to age, gender, drug type consumed, etc., which requires advanced mathematical tools such as stratified models or partial differential equations [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] . Despite these limitations, our model still offers a universally applicable tool for prediction and analysis of the drug situation, and the complex issues listed above will be considered in our future research.

In this study, we have formulated a drug epidemic model with a complete addictionrehabilitation-recovery process, which allows the generation of new initiates under the influence of drug addicts undergoing treatment and hidden drug addicts. We have established the basic properties of the model system, including the existence, uniqueness, and positivity of the solution, the forward-invariance of the feasible region, and the formulation of the basic reproduction number. We have shown that the drug-free equilibrium is globally asymptotically stable when R 0 ≤ 1, and the drug-persistent equilibrium is globally asymptotically stable when R 0 > 1. We have also carried out sensitivity analysis based on normalized forward sensitivity indices and numerical simulations, and we found that the relative efficiencies of parameters in adjusting R 0 can be debatable. Based on the established qualitative results, we use the model to simulate the drug epidemic in China, generate projections of the future, and provide in-depth discussions of intervention strategies. The simulation results show that the drug epidemic undergoes significant decreases in the following decade, and it would be more efficient to strengthen the investigation and admission of implicit drug users in order to pace up the elimination of drug spread. However, it should also be noticed that our model is a simplification of the real situation and can be further enhanced in its mathematical form. Further research could take into account the heterogeneous nature of the human society and incorporate complexities (e.g., delayed diffusive equations, multi-layer stochastic equations, and co-transmission models of various types of drugs) into their models. Moreover, regional-level modeling studies could also be carried out based on the availability of historical data, which arouses the needs for co-operation among epidemiologists, public health specialists, and governmental authorities. 

The authors declare that they have no conflict of interest.

Proof of Theorem 1. To investigate the local asymptotic stability of the drug-free equilibrium E 0 = λ µ , 0, 0, 0, 0 , we obtain the Jacobian matrix associated with System (1) evaluated at E 0 :

Through direct calculations, we see that the eigenvalues of J E 0 are −µ (multiplicity 2) and the solution to the cubic equation a 3 x 3 + a 2 x 2 + a 1 x + a 0 = 0, where

It is obvious that a 3 > 0, a 2 > 0, and from

we can easily obtain that a 0 > 0. It can be deduced that k 1 β 1 λ µ < (k 1 + k 2 + µ)(r + µ) and k 2 β 2 λ µ < (k 1 + k 2 + µ)(α + µ + µ d ).

Hence, we have a 1 > (r + µ)(α + µ + µ d ) > 0. Moreover, it can also be obtained that a 2 a 1 − a 3 a 0 > (α + µ + µ d )(r + µ)(2k 1 + 2k 2 + 4µ + r + α + µ d ) + k 2 αβ 1 λ µ > 0.

According to Routh-Hurwitz criterion, all eigenvalues of J E 0 have negative real parts when a 0 , a 1 , a 2 , a 3 > 0 and a 2 a 1 − a 3 a 0 > 0, which leads to the conclusion that E 0 is locally asymptotically stable when R 0 < 1 [64] .

To prove the global asymptotic stability of E 0 , we need to construct a suitable Lyapunov function. By applying the method introduced by Li et al. [77] , we define b ≥ 0 as the left eigenvector of the non-negative matrix V −1 F with respect to the eigenvalue R 0 . We can easily obtain that b = (0, β 1 , β 2 ) is a suitable eigenvector and further establish the following Lyapunov function:

Compute the time-derivative of L along the trajectory of System (1): dL dt = µR 0 λ · dI 1 dt + β 1 r+µ · dI 2 dt + αβ 1 +(r+µ)β 2 (r+µ)(α+µ+µ d ) · dI 3 dt = µR 0 λ [β 1 SI 2 + β 2 SI 3 − (k 1 + k 2 + µ)I 1 ] + β 1 r+µ [k 1 I 1 + αI 3 − (r + µ)I 2 ] + αβ 1 +(r+µ)β 2 (r+µ)(α+µ+µ d ) [k 2 I 1 − (α + µ + µ d )I 3 ].

Through tedious algebraic manipulations, we acquire the coefficients of I 1 ∼ I 3 . Coefficient of I 1 :

when R 0 ≤ 1, Coefficient of I 3 : µR 0 λ β 2 S − β 2 ≤ β 2 (R 0 − 1) ≤ 0 when R 0 ≤ 1. The equalities above can be satisfied only when R 0 = 1, which corresponds to the drug-free equilibrium E 0 = λ µ , 0, 0, 0, 0 . By LaSalle's invariance principle, the largest invariant set in Γ = (S, I 1 , I 2 , I 3 , R) ∈ R 5 + ; 0 ≤ S + I 1 + I 2 + I 3 + R ≤ λ µ is reduced to the singleton E 0 in this case. Hence, the unique drug-free equilibrium is globally asymptotically stable; thus, it completes the proof.

Proof of Theorem 2. The Jacobian matrix method seems impractical for the analysis of local stability of the drug-persistent equilibrium E * , since determination of the eigenvalues requires solving a quartic equation, which is difficult to handle even for Routh-Hurwitz criterion. Alternatively, we choose the bifurcation method based on center manifold theory by Castillo-Chavez and Song [78] .

Consider System (1) as an ODE system with bifurcation parameter ϕ:

The choice of the bifurcation parameter ϕ should satisfy f (E 0 , φ) ≡ 0 for ∀φ. In this case, β 1 is obviously an eligible parameter for φ. To proceed, we substitute the variables as S = x 1 , I 1 = x 2 , I 2 = x 3 , I 3 = x 4 , R = x 5 , and System (1) changes to:

Set R 0 = [(α+µ+µ d )k 1 +αk 2 ]β 1 +k 2 (r+µ)β 2 (k 1 +k 2 +µ)(r+µ)(α+µ+µ d ) · λ µ = 1, and based on previous results of the Jacobian matrix J E 0 , it is not difficult to find that zero is an eigenvalue of J E 0 , and all other eigenvalues have negative real parts. Hence, the criterion for the application of the bifurcation method is met. Subsequently, we compute the right and left eigenvectors corresponding to zero eigenvalue. The right eigenvector is given by → w = (w 1 , w 2 , w 3 , w 4 , w 5 ) T , where w 1 = − k 1 +k 2 +µ µ · w 2 , w 3 = µ(k 1 +k 2 +µ)(α+µ+µ d )−k 2 β 2 λ (α+µ+µ d )β 1 λ · w 2 ,

The left eigenvalue is given by

Let f k be the kth equation of System (5) , and according to the bifurcation method,

. The associated non-zero second-order partial derivatives around (E 0 , 0) are: ∂ 2 f 2 ∂x 1 ∂x 3 = ∂ 2 f 2 ∂x 3 ∂x 1 = β 1 , ∂ 2 f 2 ∂x 1 ∂x 4 = ∂ 2 f 2 ∂x 4 ∂x 1 = β 2 , ∂ 2 f 2 ∂x 3 ∂φ = λ µ . The rest of the second-order partial derivatives were all calculated to be 0 and are hence omitted. Take v 2 = w 2 = 1, and it is easy to obtain that:

w i w j ∂ 2 f 2 ∂x i ∂x j (0, 0) = 2(w 1 w 3 φ + w 1 w 4 β 2 ).

It is obvious that w 1 < 0, w 4 > 0, and given the condition of R 0 = 1, it is not hard to prove w 3 > 0. Hence, we have a < 0 proved. On the other hand, we have:

In conclusion, since a < 0 and b > 0, according to the theory of backward bifurcation, the model system experiences forward bifurcation, and the drug-persistent equilibrium is locally asymptotically stable for R 0 > 1 but close to one [40] .

To prove the global stability of the drug-persistent equilibrium E * , we now construct a Volterra-type Lyapunov function following the method described in [77] :

among which: V 1 = S − S * − S * ln S S * + I 1 − I * 1 − I * 1 ln I 1

National Narcotics Control Committee. Drug Situation in China

An Epidemiologic Assessment of Heroin Use

Epidemiology Abuse: Epidemiological and Psychosocial Models of Drug Abuse

A mathematical model of a heroin epidemic: Implications for control policies

A contribution to the mathematical theory of epidemics

Mathematical Models and Dynamics of Infectious Disease

Uncertainty and sensitivity analysis of the basic reproduction number of a vaccinated epidemic model of influenza

Effects of behavioral changes in a smallpox attack model

Modeling and analysis of COVID-19 epidemics with treatment in fractional derivatives using real data from Pakistan

Chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells

The role of prostitution on HIV transmission with memory: A modeling approach

Analysis of a drinking epidemic model

Dynamical analysis of cigarette smoking model with a saturated incidence rate

Epidemics of computer viruses: A complex-network approach

An epidemic model of rumor diffusion in online social networks

Heroin epidemics, treatment and ODE modelling

A note on heroin epidemics

Dynamic behaviour for a nonautonomous heroin epidemic model with time delay

Global behaviour of a heroin epidemic model with distributed delays

A note on global stability for a heroin epidemic model with distributed delay

Global stability for a heroin model with two distributed delays

Global asymptotic properties of a heroin epidemic model with treat-age

Global stability for a heroin model with age-dependent susceptibility

Global dynamics of a heroin epidemic model with age structure and nonlinear incidence

A Heroin Epidemic Model: Very General Non Linear Incidence, Treat-Age, and Global Stability

Mathematical Analysis for an Age-Structured Heroin Epidemic Model

Qualitative analysis on a diffusive age-structured heroin transmission model

Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate

Global dynamical analysis of a heroin epidemic model on complex networks

Threshold dynamics of a delayed multi-group heroin epidemic model in heterogeneous populations

Analysis of an age-structured multi-group heroin epidemic model

Dynamics of stochastic heroin epidemic model with lévy jumps

Dynamics of a stochastic heroin epidemic model

Dynamics of a stochastic heroin epidemic model with bilinear incidence and varying population size

Dynamics of the stochastically perturbed Heroin epidemic model under non-degenerate noises

Numerical treatment of stochastic heroin epidemic model

From heroin epidemics to methamphetamine epidemics: Modelling substance abuse in a South African province

A theoretical model for substance abuse in the presence of treatment

Modelling the Dynamics of Crystal Meth ('Tik') Abuse in the Presence of Drug-Supply Chains in South Africa

Modelling the trends of inpatient and outpatient rehabilitation for methamphetamine in the Western Cape province of South Africa

Modelling Drug Abuse Epidemics in the Presence of Limited Rehabilitation Capacity

On the Role of Imitation on Adolescence Methamphetamine Abuse Dynamics

Estimates of the national trend of drugs use during 2000-2030 in China: A population-based mathematical model

Complete global analysis of an SIRS epidemic model with graded cure and incomplete recovery rates

Global Dynamics of a Discretized Heroin Epidemic Model with Time Delay

Analysis of a Heroin Epidemic Model with Saturated Treatment Function

Dynamics of an age-structured heroin transmission model with vaccination and treatment

Global dynamic of a heroin epidemic model

Global dynamics in a heroin epidemic model with different conscious stages and two distributed delays

Dynamics of synthetic drugs transmission model with psychological addicts and general incidence rate

The Role of Family on the Transmission Model of Methamphetamine

Synthetic drugs transmission: Stability analysis and optimal control

Modelling and stability of a synthetic drugs transmission model with relapse and treatment

Hopf Bifurcation Analysis of a Synthetic Drug Transmission Model with Time Delays

The analysis of a drug transmission model with family education and public health education

European Monitoring Centre for Drugs and Drug Addiction (EMCDDA)

Methods for the Identification of Drug Addiction

Pilot trial of a recovery management intervention for heroin addicts released from compulsory rehabilitation in China

Evaluation of an implementation of methadone maintenance treatment in China

Needle and syringe programs in Yunnan, China yield health and financial return

Tracking the evolution of drug abuse in China

Introduction to Functional Differential Equations

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Lecture notes, forward and adjoint sensitivity analysis: With applications in Dynamical Systems. In Linear Algebra and Optimisation Mathematical and Theoretical Biology Institute

Analysis and considerations of the current detoxification method in China

On enforcement of drug trafficking prohibition between the western China and neighboring countries

Dynamics of a Heroin Epidemic Model with Very Population

Illicit drug initiation among institutionalized drug users in China

Reported Reasons for Initiating Drug Use among Drug-Dependent Adolescents and Youths in Yunnan

A review of the history of drug control publicity and education in China

From Abstinence to Relapse: A Preliminary Qualitative Study of Drug Users in a Compulsory Drug Rehabilitation Center in Changsha, China

Modeling direct and indirect disease transmission using multi-group model

Dynamical Models of Tuberculosis and Their Applications

The time-derivative of V 1 is given by:Based on System (3) at the drug-persistent equilibrium E * , we obtain λ = β 1 S * I * 2 + β 2 S * I * 3 + µS * and k 1 + k 2 + µ =, and the equation above turns into: [78] . It is easy to prove Φ(a) = 1 − a + ln a ≤ 0 for ∀a > 0, and the equality can be reached only when a = 1. Thus, the expressions above can be rewritten as:Hence, it can be deduced that:Similarly, the time-derivative of V 2 is given by:Based on System (3) at the drug-persistent equilibrium E * , we obtain r + µ = k 1 I * 1 +αI * 3 I * 2 , and the equation above turns into:Hence, it is easily obtained that dV 2 dt ≤ k 1 I * 1 [D(I 2 ) − D(I 1 )] + αI * 3 [D(I 2 ) − D(I 3 )]. By similar algebraic manipulations, we can also obtain that: As a consequence, the time-derivative of the Lyapunov function becomes:Hence, the right-hand side of the inequality is identically vanishing, and dV dt ≤ 0 holds whenever R 0 > 1. The equalities above can be satisfied only when S = S * , I 1 = I * 1 , I 2 = I * 2 and I 3 = I * 3 , which corresponds to the drug-persistent equilibrium E * = S * , I * 1 , I * 2 , I * 3 , R * . By LaSalle's invariance principle, the largest invariant set in Γ = (S, I 1 , I 2 , I 3 , R) ∈ R 5 + ; 0 ≤ S + I 1 + I 2 + I 3 + R ≤ λ µ is reduced to the singleton E * in this case. Hence, the drug-persistent equilibrium is globally asymptotically stable, thus completing the proof.