Modeling the Dynamics of Drug Spreading in China

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.


Introduction
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%

Basic Assumptions
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.

S(t)
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

Model Formulations
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:

Existence and Uniqueness of the Solution
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 × 

Existence and Uniqueness of the Solution
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].

Feasible Region
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 Γ.

Positivity of Solutions
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.
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.

The Basic Reproduction Number
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.

Existence of the Equilibria
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 .

2.
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.

Global Stability of the Drug-Free Equilibrium
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.

Global Stability of the Drug-Persistent Equilibrium
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.

Sensitivity Analysis
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.

Simulation of Sensitivity
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.
(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 .

Figure 2.
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).

Simulation of Model Equilibria
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.

Application of the Model
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

Application of the Model
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.

Data Source and Variables
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.

Parameter Estimation
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.

Model Fitting and Projections
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.

Evaluation of Intervention Strategies
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.

Evaluation of Intervention Strategies
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).

Discussion
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.

Conclusions
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.

Conflicts of Interest:
The authors declare that they have no conflict of interest.
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].

Appendix B
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].
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