Next Article in Journal
A Comparison of the Gluco-Regulatory Responses to High-Intensity Interval Exercise and Resistance Exercise
Previous Article in Journal
Blood Pressure and Tooth Loss: A Large Cross-Sectional Study with Age Mediation Analysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling the Dynamics of Drug Spreading in China

1
Department of Biomedical Engineering, College of Engineering, Peking University, Beijing 100871, China
2
School of Mathematics, Taiyuan University of Technology, Taiyuan 030024, China
3
School of Public Health, Peking University, Beijing 100191, China
4
State Key Lab of Bioelectronics, School of Biological Science and Medical Engineering, Southeast University, Nanjing 210096, China
5
Center for Intelligent Public Health, Institute for Artificial Intelligence, Peking University, Beijing 100191, China
*
Author to whom correspondence should be addressed.
Int. J. Environ. Res. Public Health 2021, 18(1), 288; https://doi.org/10.3390/ijerph18010288
Submission received: 3 November 2020 / Revised: 6 December 2020 / Accepted: 30 December 2020 / Published: 2 January 2021

Abstract

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

1. 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% of the People Who Inject Drugs (PWIDs) reside in three countries: China, Russia, and the United States [1]. All the evidence above aroused the need for global attention to the problem of drug abuse as well as targeted intervention to address this issue.
The mathematical model has long been an effective tool in the area of public health, and extensive research has been conducted in terms of drug users and drug-using behaviors. In the 1970s, the analogy between heroin use and communicable disease was proposed, which verified the validity and utility of the epidemiologic approach to studying heroin use [4,5,6]. This paved the way for the utilization of compartmental dynamic models, which could date back to 1926, when Kermack and McKendrick formulated the famous Susceptible–Infected–Recovered (SIR) compartmental model while studying the Great Plague of London during 1665−1666 and the 1906 plague in Bombay [7,8]. Compartmental dynamic modeling has undergone extensive development in the past decades and has proved itself as an effective tool in the research of infectious diseases, including influenza, smallpox, coronavirus, HIV, etc. [9,10,11,12,13]. In recent years, the application of compartmental dynamic models has been extended to other research fields, for instance, dynamics of the spread of alcoholism, cigarette smoking, internet virus, rumors, or drug-using behavior [14,15,16,17].
After several explorative research studies since the 1980s, a classic three-compartment dynamical model was proposed by White and Comiskey in 2006, which paved the way for drug epidemic models to come [18]. The authors utilized ordinary differential equations (ODE) systems and made calculations on several key aspects of the model, including basic reproduction number (R0), drug-free equilibrium, and drug-persistent equilibrium. Later on, with the aid of development in theories of nonlinear dynamic systems and computer-assisted simulation tools, plenty of research studies have sprung up worldwide, adding to this field of the drug epidemic model. A secondary analysis of the White-Comiskey model was conducted by Mulone et al. in 2009, who loosened the assumption of constant inflow rate [19]. Other modification studies, which were mostly heroin epidemic models, took advantage of various mathematical tools to account for practical factors. For example, delayed differential equations were used to simulate processes with known durations [20,21,22,23], partial differential equations were utilized to incorporate the effect of age or treatment duration [24,25,26,27,28,29], multi-layered models were proposed when population heterogeneity was involved [30,31,32,33], and stochastic differential equations were formulated to reflect unexpected fluctuations in reality [34,35,36,37,38]. In addition to theoretical analyses, some synthetic drug epidemic models were applied to real settings and fitted to historical data, most of which were based on methamphetamine epidemics in South Africa [39,40,41,42,43,44].
Despite these modeling efforts, few of them have investigated the drug situation in China [45]. Hence, it is our objective to model the scale and trend of the drug epidemic in China and thoroughly discuss specified intervention strategies. In this article, we formulate a drug epidemic model with a complete addiction–rehabilitation–recovery process. Unlike many other studies, we do not allow self-abstinence without treatment or any other unrealistic assumptions incompatible with the social background in China [46,47,48,49,50,51,52,53,54,55,56,57]. After qualitative analyses of the theoretical behaviors of the model system, we fit the model to historical data of drug abuse in China and make projections of the future. Efficiencies of various interventions are discussed. Our study is innovative from two aspects. First, a new version of the drug epidemic model is applied to the drug abuse epidemic in China. Second, new insights into intervention strategies to curtail the drug epidemic are provided.
This article is arranged as follows: In Section 2, we formulate the model and establish its basic properties. The existence and stability of the model equilibria was discussed in Section 3. In Section 4, sensitivity analysis and some numerical results were provided. Section 5 includes parameter estimation, model fitting, and projections, as well as simulations of intervention efficiencies. In Section 6, we conclude this paper with detailed discussions.

2. Model Assumptions, Formulations, and Basic Properties

2.1. 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 (I1), drug addicts undergoing treatment (I2), hidden drug addicts (I3), 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 (I1) 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 (I2) 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 I2 compartment is actually a mixture of patients receiving various forms of detoxification rehabilitation. Hidden drug addicts (I3) 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 (I2) or hidden drug addicts (I3), in which process we adopt a bilinear incidence rate. The compartment of drug addicts undergoing treatment (I2), 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 (I3) 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 (I3) could also be transferred to the I2 compartment at admission rate α , and those who finish rehabilitation will leave I2 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 (I3) 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.

2.2. 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:
d S d t = λ β 1 S I 2 β 2 S I 3 μ S d I 1 d t = β 1 S I 2 + β 2 S I 3 ( k 1 + k 2 + μ ) I 1 d I 2 d t = k 1 I 1 + α I 3 ( r + μ ) I 2 d I 3 d t = k 2 I 1 ( α + μ + μ d ) I 3 d R d t = r I 2 μ R
with initial conditions S ( 0 ) 0 , I 1 ( 0 ) 0 , I 2 ( 0 ) 0 , I 3 ( 0 ) 0 , R ( 0 ) 0 .

2.3. Existence and Uniqueness of the Solution

Lemma 1.
The model system admits a unique solution for non-negative initial conditions S ( 0 ) 0 , I 1 ( 0 ) 0 , I 2 ( 0 ) 0 , I 3 ( 0 ) 0 , R ( 0 ) 0 .
Proof Suppose C ( [ a , b ] , 5 ) is the Banach space of continuous functions mapping the interval [ a , b ] into 5 with the topology of uniform convergence. System (1) with non-negative initial conditions can be considered as the following initial-value problem:
{ x ˙ ( t ) = f ( t , x ) x ( t 0 ) = x 0 .
Let Ω be an open subset in × C and f : Ω 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].

2.4. 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
Γ = { ( S , I 1 , I 2 , I 3 , R ) + 5 ; 0 S + I 1 + I 2 + I 3 + R λ μ }
with initial conditions S ( 0 ) 0 , I 1 ( 0 ) 0 , I 2 ( 0 ) 0 , I 3 ( 0 ) 0 , R ( 0 ) 0 , Γ 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
d N d t = λ μ N μ d I 3 λ μ N .
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 Γ .

2.5. 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 .
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 : S > 0 , I 1 > 0 , I 2 > 0 , I 3 > 0 , R > 0 } [ 0 , t ] . Thus, t ¯ > 0 , and from the first equation of System (1), we have d S d t = λ β 1 S I 2 β 2 S I 3 μ S .
Based on the Leibniz integral rule,
d d t 0 t ( β 1 I 2 + β 2 I 3 ) d s = β 1 I 2 + β 2 I 3 + 0 t t ( β 1 I 2 + β 2 I 3 ) d s β 1 I 2 + β 2 I 3 d d t [ e μ t + 0 t ( β 1 I 2 + β 2 I 3 ) d s ] ( μ + β 1 I 2 + β 2 I 3 ) e μ t + 0 t ( β 1 I 2 + β 2 I 3 ) d s d d t [ S ( t ) e μ t + 0 t ( β 1 I 2 + β 2 I 3 ) d s ] [ d S ( t ) d t + S ( t ) ( μ + β 1 I 2 + β 2 I 3 ) ] e μ t + 0 t ( β 1 I 2 + β 2 I 3 ) d s = λ e μ t + 0 t ( β 1 I 2 + β 2 I 3 ) d s S ( t ¯ ) e μ t ¯ + 0 t ¯ ( β 1 I 2 + β 2 I 3 ) d s S ( 0 ) 0 t ¯ λ e μ t ¯ + 0 t ¯ ( β 1 I 2 + β 2 I 3 ) d s d t ¯ , S ( t ¯ ) e [ μ t ¯ + 0 t ¯ ( β 1 I 2 + β 2 I 3 ) d s ] [ S ( 0 ) + 0 t ¯ λ e μ t ¯ + 0 t ¯ ( β 1 I 2 + β 2 I 3 ) d s d t ¯ ] > 0 .
From the second equation of System (1), we obtain
d I 1 d t ( k 1 + k 2 + μ ) I 1 .
Solving the differential inequality, we have:
I 1 ( t ) I 1 ( 0 ) e ( k 1 + k 2 + μ ) t > 0 .
Similarly, it can be shown that I 2 ( t ) > 0 , I 3 ( t ) > 0 and R ( t ) > 0 for all t > 0 .

3. Model Equilibria Analysis

3.1. 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:
d x d t = V = [ β 1 S I 2 + β 2 S I 3 0 0 0 0 ] [ ( k 1 + k 2 + μ ) I 1 ( r + μ ) I 2 k 1 I 1 α I 3 ( α + μ + μ d ) I 3 k 2 I 1 β 1 S I 2 + β 2 S I 3 + μ S λ μ R r I 2 ]
where represents the new initiate terms, and V corresponds to the interior transfer terms. The corresponding linearized matrices of and V evaluated at the drug-free equilibrium E 0 = ( λ μ , 0 , 0 , 0 , 0 ) are
F = [ 0 β 1 S 0 β 2 S 0 0 0 0 0 0 0 ] ,   V = [ k 1 + k 2 + μ 0 0 k 1 r + μ α k 2 0 α + μ + μ d ] .
The next-generation matrix is given by F V 1 = [ A 11 A 12 A 13 0 0 0 0 0 0 ] , where
A 11 = [ ( α + μ + μ d ) k 1 + α k 2 ( k 1 + k 2 + μ ) ( r + μ ) ( α + μ + μ d ) β 1 + k 2 ( k 1 + k 2 + μ ) ( α + μ + μ d ) β 2 ] λ μ ,
A 12 = β 1 λ μ ( r + μ ) , A 13 = [ α β 1 ( r + μ ) ( α + μ + μ d ) + β 2 α + μ + μ d ] λ μ .
Thus, the basic reproduction number R 0 is defined as the spectral radius of the next-generation matrix, which is the largest absolute value of its eigenvalues:
R 0 = ρ ( F V 1 ) = [ ( α + μ + μ d ) k 1 + α k 2 ( k 1 + k 2 + μ ) ( r + μ ) ( α + μ + μ d ) β 1 + k 2 ( k 1 + k 2 + μ ) ( α + μ + μ d ) β 2 ] λ μ .
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, 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 k 2 β 2 λ μ ( k 1 + k 2 + μ ) ( α + μ + μ d ) . To sum up, the formula of basic reproduction number represents the overlay of the influence of drug addicts undergoing treatment and hidden drug addicts.

3.2. Existence of the Equilibria

In order to investigate the existence of the drug-free equilibrium and the drug-persistent equilibrium, we set the left side of the equations of System (1) at zero:
λ β 1 S I 2 β 2 S I 3 μ S = 0 β 1 S I 2 + β 2 S I 3 ( k 1 + k 2 + μ ) I 1 = 0 k 1 I 1 + α I 3 ( r + μ ) I 2 = 0 k 2 I 1 ( α + μ + μ d ) I 3 = 0 r I 2 μ R = 0 .
Through direct calculations, we obtain:
[ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + ( r + μ ) k 2 β 2 ( r + μ ) ( α + μ + μ d ) S I 1 = ( k 1 + k 2 + μ ) I 1 .
  • When I 1 = 0 , it is easy to acquire I 2 = I 3 = R = 0 , S = λ μ . At this point, N = S = λ μ , and d N d t = λ μ 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 I 1 0 , it is obvious that S * = ( r + μ ) ( α + μ + μ d ) ( k 1 + k 2 + μ ) [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + ( r + μ ) k 2 β 2 and I 1 * = λ μ S * k 1 + k 2 + μ , I 2 * = ( α + μ + μ d ) k 1 + α k 2 ( r + μ ) ( α + μ + μ d ) I 1 * , I 3 * = k 2 α + μ + μ d I 1 * , R * = r μ I 2 * . At this point, we obtain N = S * + I 1 * + I 2 * + I 3 * + R * = [ ( α + μ + μ d ) ( k 1 + μ ) + k 2 ( α + μ ) ] λ + k 2 μ μ d S * μ ( k 1 + k 2 + μ ) ( α + μ + μ d ) , and it can be easily verified that d N d t = λ μ N μ d I 3 = 0 also holds. We refer to this point as the unique drug-persistent 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 .

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

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

4. Sensitivity and Numerical Simulations

4.1. 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:
A β 1 = | R 0 β 1 | | β 1 R 0 | = [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + k 2 ( r + μ ) β 2 , A β 2 = | R 0 β 2 | | β 2 R 0 | = k 2 ( r + μ ) β 2 [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + k 2 ( r + μ ) β 2 , A k 1 = | R 0 k 1 | | k 1 R 0 | = k 1 k 1 + k 2 + μ | ( α + μ + μ d ) ( k 1 + μ ) β 1 k 2 [ α β 1 + ( r + μ ) β 2 ] | [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + k 2 ( r + μ ) β 2 , A α = | R 0 α | | α R 0 | = α α + μ + μ d | ( μ + μ d ) k 2 β 1 ( r + μ ) k 2 β 2 | [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + k 2 ( r + μ ) β 2 .
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.

4.2. 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.
(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 [ 1 e 7 , 1 e 5 ] and β 2 [ 1 e 7 , 1 e 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 = 1 e 6 . Set β 2 [ 1 e 7 , 1 e 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 , μ = 0.007 , μ d = 0.025 , k 1 = 0.05 , k 2 = 0.2 , r = 0.5 , and β 1 = 1 e 6 . Set β 2 [ 1 e 7 , 1 e 5 ] 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 .

4.3. 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 = 1 e 7 , and β 2 = 1 e 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.
(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 = 1 e 6 , and β 2 = 1 e 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.

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

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

5.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 < 1 e 6 /10 thousand people*year from this inequality. We list the parameter units, value ranges, final values used, and their sources in Table 3.

5.3. Model Fitting and Projections

In this section, we fit our model to historical data of drug abuse in China, where reported numbers of existing drug users during 2015–2019 were used to model the growth of drug addicts undergoing treatment. Instead of using all five data points in the curve-fitting procedure, we set the point of 2019 aside and used it for verification. Least square curve fitting was realized through the fminsearch function in the Matlab program (Mathworks Corp, Natick, MA, USA), and numerical solutions of the model system were obtained by the Runge–Kutta method of order 4. The iteration procedures of fitting are plotted in Figure 7a, where blue crosses are historical data of existing drug users during 2015–2018. From Figure 7b to Figure 7f, time series of I 2 , N , I 1 , I 3 and R are plotted along with their projections until 2030. Blue crosses in Figure 7b have the same meanings as in Figure 7a, and the red cross marks the corresponding data in 2019. Crosses in Figure 7c represent historical data of population aged 15−64 recorded by The National Bureau of Statistics, which acted as verification with the sum of all compartments. It can be seen from the results in 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.

5.4. 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.
(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%).
(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.

6. Discussion

In this article, we propose a drug epidemic model with a complete addiction–rehabilitation–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.86 e 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–rehabilitation–recovery 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.

7. Conclusions

In this study, we have formulated a drug epidemic model with a complete addiction–rehabilitation–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.

Author Contributions

Z.J. and Z.L. designed the study and established the model, M.L., X.Y. and H.T. analyzed the model and performed simulations, X.Y. and H.T. explained the results, X.Y. and H.T. wrote the initial manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Nature Science Foundation (grant no. U1611264, 91546203 and 11801398), Beijing Advanced Discipline Construction Project (grant no. BMU2019GJJXK005) and National key R&D Program of China (grant no. 2016YFA0501600, 2016YFA0501603).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Publicly available datasets were analyzed in this study. The data can be found in the following link: [http://www.nncc626.com/2020-06/24/c_1210675813.htm] and the extended links therein.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Appendix A

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 :
J | E 0 = [ μ 0 β 1 λ μ β 2 λ μ 0 0 ( k 1 + k 2 + μ ) β 1 λ μ β 2 λ μ 0 0 k 1 ( r + μ ) α 0 0 k 2 0 ( α + μ + μ d ) 0 0 0 r 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
a 3 = 1 , a 2 = α + 3 μ + μ d + k 1 + k 2 + r , a 1 = ( k 1 + k 2 + μ ) ( α + μ + μ d ) + ( r + μ ) ( k 1 + k 2 + 2 μ + α + μ d ) ( k 1 β 1 + k 2 β 2 ) λ μ , a 0 = β 1 λ μ [ k 2 α + k 1 ( α + μ + μ d ) ] + ( r + μ ) [ ( k 1 + k 2 + μ ) ( α + μ + μ d ) k 2 β 2 λ μ ] .
It is obvious that a 3 > 0 , a 2 > 0 , and from
R 0 = [ ( α + μ + μ d ) k 1 + α k 2 ] β 1 + k 2 ( r + μ ) β 2 ( k 1 + k 2 + μ ) ( r + μ ) ( α + μ + μ d ) λ μ < 1 ,
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 + μ ) ( 2 k 1 + 2 k 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:
L = b V 1 x = μ R 0 λ I 1 + β 1 r + μ I 2 + α β 1 + ( r + μ ) β 2 ( r + μ ) ( α + μ + μ d ) I 3 .
Compute the time-derivative of L along the trajectory of System (1):
d L d t = μ R 0 λ d I 1 d t + β 1 r + μ d I 2 d t + α β 1 + ( r + μ ) β 2 ( r + μ ) ( α + μ + μ d ) d I 3 d t = μ R 0 λ [ β 1 S I 2 + β 2 S I 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 :
β 1 k 1 r + μ + α k 2 β 1 + ( r + μ ) k 2 β 2 ( r + μ ) ( α + μ + μ d ) ( k 1 + k 2 + μ ) μ R 0 λ = 0 ,
Coefficient of I 2 :
μ R 0 λ β 1 S β 1 β 1 ( R 0 1 ) 0
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 ) + 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. □

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].
Consider System (1) as an ODE system with bifurcation parameter φ :
d x d t = f ( x , ϕ ) , f : 5 × 5 , f 2 ( 5 × ) .
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:
f 1 = λ β 1 x 1 x 3 β 2 x 1 x 4 μ x 1 f 2 = β 1 x 1 x 3 + β 2 x 1 x 4 ( k 1 + k 2 + μ ) x 2 f 3 = k 1 x 2 + α x 4 ( r + μ ) x 3 f 4 = k 2 x 2 ( α + μ + μ d ) x 4 f 5 = r x 3 μ x 5 .
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 , w 4 = k 2 α + μ + μ d w 2 , w 5 = μ ( k 1 + k 2 + μ ) ( α + μ + μ d ) k 2 β 2 λ ( α + μ + μ d ) β 1 λ μ r w 2 .
The left eigenvalue is given by v = ( 0 , v 2 , v 3 , v 4 , 0 ) , where
v 3 = β 1 λ ( r + μ ) μ v 2 , v 4 = μ ( k 1 + k 2 + μ ) ( r + μ ) k 1 β 1 λ μ k 2 ( r + μ ) v 2
Let f k be the kth equation of System (5), and according to the bifurcation method,
a = k , i , j = 1 5 v k w i w j 2 f k x i x j ( 0 , 0 ) , and b = k , i = 1 5 v k w i 2 f k x i ϕ ( 0 , 0 ) .
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:
a = v 2 i , j = 1 5 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:
b = v 2 i 5 w i 2 f 2 x i ϕ ( 0 , 0 ) = w 3 λ μ > 0 .
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]:
V = μ R 0 ( r + μ ) λ V 1 + β 1 V 2 + β 2 ( r + μ ) + α β 1 α + μ + μ d V 3 ,
among which:
V 1 = S S * S * ln S S * + I 1 I 1 * I 1 * ln I 1 I 1 * V 2 = I 2 I 2 * I 2 * ln I 2 I 2 * ,   V 3 = I 3 I 3 * I 3 * ln I 3 I 3 * .
The time-derivative of V 1 is given by:
d V 1 d t = ( 1 S * S ) d S d t + ( 1 I 1 * I 1 ) d I 1 d t = ( 1 S * S ) ( λ β 1 S I 2 β 2 S I 3 μ S ) + ( 1 I 1 * I 1 ) [ β 1 S I 2 + β 2 S I 3 ( k 1 + k 2 + μ ) I 1 ] .
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 + μ = β 1 S * I 2 * + β 2 S * I 3 * I 1 * , and the equation above turns into:
d V 1 d t = ( 1 S * S ) [ β 1 S * I 2 * ( 1 S I 2 S * I 2 * ) + β 2 S * I 3 * ( 1 S I 3 S * I 3 * ) + μ ( S * S ) ] + ( 1 I 1 * I 1 ) [ β 1 S * I 2 * ( S I 2 S * I 2 * I 1 I 1 * ) + β 2 S * I 3 * ( S I 3 S * I 3 * I 1 I 1 * ) ] Θ 1 = ( 1 S * S ) ( S * S ) 0 d V 1 d t β 1 S * I 2 * ( 2 S * S + I 2 I 2 * I 1 I 1 * I 1 * S I 2 I 1 S * I 2 * ) + β 2 S * I 3 * ( 2 S * S + I 3 I 3 * I 1 I 1 * I 1 * S I 3 I 1 S * I 3 * ) .
Let F = 2 S * S + I 2 I 2 * I 1 I 1 * I 1 * S I 2 I 1 S * I 2 * , G = 2 S * S + I 3 I 3 * I 1 I 1 * I 1 * S I 3 I 1 S * I 3 * , and define D ( x ) = x x * + ln x x * , Φ ( a ) = 1 a + ln a as in [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:
F = Φ ( S * S ) ln ( S * S ) + Φ ( I 1 * S I 2 I 1 S * I 2 * ) ln ( I 1 * S I 2 I 1 S * I 2 * ) + I 2 I 2 * I 1 I 1 * I 2 I 2 * I 1 I 1 * ln ( I 1 * I 2 I 1 I 2 * ) = D ( I 1 ) D ( I 2 ) , G = Φ ( S * S ) ln ( S * S ) + Φ ( I 1 * S I 3 I 1 S * I 3 * ) ln ( I 1 * S I 3 I 1 S * I 3 * ) + I 3 I 3 * I 1 I 1 * I 3 I 3 * I 1 I 1 * ln ( I 1 * I 3 I 1 I 3 * ) = D ( I 1 ) D ( I 3 ) .
Hence, it can be deduced that:
d V 1 d t β 1 S * I 2 * [ D ( I 1 ) D ( I 2 ) ] + β 2 S * I 3 * [ D ( I 1 ) D ( I 3 ) ] .
Similarly, the time-derivative of V 2 is given by:
d V 2 d t = ( 1 I 2 * I 2 ) d I 2 d t = ( 1 I 2 * I 2 ) [ k 1 I 1 + α I 3 ( r + μ ) I 2 ] .
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:
d V 2 d t = ( 1 I 2 * I 2 ) [ k 1 I 1 * ( I 1 I 1 * I 2 I 2 * ) + α I 3 * ( I 3 I 3 * I 2 I 2 * ) ] = k 1 I 1 * ( I 1 I 1 * I 2 I 2 * I 1 I 2 * I 1 * I 2 + 1 ) + α I 3 * ( I 3 I 3 * I 2 I 2 * I 2 * I 3 I 2 I 3 * + 1 ) = k 1 I 1 * [ I 1 I 1 * I 2 I 2 * + Φ ( I 1 I 2 * I 1 * I 2 ) ln ( I 1 I 2 * I 1 * I 2 ) ] + α I 3 * [ I 3 I 3 * I 2 I 2 * + Φ ( I 2 * I 3 I 2 I 3 * ) ln ( I 2 * I 3 I 2 I 3 * ) ] k 1 I 1 * ( I 1 I 1 * I 2 I 2 * ln I 1 I 1 * + ln I 2 I 2 * ) + α I 3 * ( I 3 I 3 * I 2 I 2 * ln I 3 I 3 * + ln I 2 I 2 * ) .
Hence, it is easily obtained that d V 2 d t 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:
d V 3 d t k 2 I 1 * ( I 1 I 1 * I 3 I 3 * ln I 1 I 1 * + ln I 3 I 3 * ) = k 2 I 1 * [ D ( I 3 ) D ( I 1 ) ] .
As a consequence, the time-derivative of the Lyapunov function becomes:
d V d t = μ R 0 ( r + μ ) λ d V 1 d t + β 1 d V 2 d t + β 2 ( r + μ ) + α β 1 α + μ + μ d d V 3 d t , μ R 0 ( r + μ ) λ { β 1 S * I 2 * [ D ( I 1 ) D ( I 2 ) ] + β 2 S * I 3 * [ D ( I 1 ) D ( I 3 ) ] } + β 1 { k 1 I 1 * [ D ( I 2 ) D ( I 1 ) ] + α I 3 * [ D ( I 2 ) D ( I 3 ) ] } + k 2 I 1 * [ β 2 ( r + μ ) + α β 1 ] α + μ + μ d [ D ( I 3 ) D ( I 1 ) ] . 0
Hence, the right-hand side of the inequality is identically vanishing, and d V d t 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 ) + 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. □

References

  1. UNODC. World Drug Report 2019; Sales No. E.19.XI.8; United Nations Publication: New York, NY, USA, 2019. [Google Scholar]
  2. UNODC. World Drug Report 2010; Sales No. E.10.XI.13; United Nations Publication: New York, NY, USA, 2010. [Google Scholar]
  3. National Narcotics Control Committee. Drug Situation in China (2019), Beijing. 2020. Available online: http://www.nncc626.com/2020-06/25/c_1210675877.htm (accessed on 27 October 2020).
  4. Greene, M.H. An Epidemiologic Assessment of Heroin Use. Am. J. Public Health 1974, 64 (Suppl. 12), 1–10. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Jacobs, P.E. Epidemiology Abuse: Epidemiological and Psychosocial Models of Drug Abuse. J. Drug Educ. 1976, 6, 259–271. [Google Scholar] [CrossRef]
  6. Mackintosh, D.R.; Stewart, G.T. A mathematical model of a heroin epidemic: Implications for control policies. J. Epidemiol. Community Health 1979, 33, 299–304. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 1927, 115, 700–721. [Google Scholar]
  8. Ma, Z.E.; Zhou, Y.C.; Wang, W.D.; Jin, Z. Mathematical Models and Dynamics of Infectious Disease, 1st ed.; Science Press: Beijing, China, 2004. [Google Scholar]
  9. Samsuzzoha, M.; Singh, M.; Lucy, D. Uncertainty and sensitivity analysis of the basic reproduction number of a vaccinated epidemic model of influenza. Appl. Math. Model. 2013, 37, 903–915. [Google Scholar] [CrossRef]
  10. Del Valle, S.; Hethcote, H.; Hyman, J.M.; Castillo-Chavez, C. Effects of behavioral changes in a smallpox attack model. Math. Biosci. 2005, 195, 228–251. [Google Scholar] [CrossRef] [PubMed]
  11. Naik, P.A.; Yavuz, M.; Qureshi, S.; Zu, J.; Townley, S. Modeling and analysis of COVID-19 epidemics with treatment in fractional derivatives using real data from Pakistan. Eur. Phys. J. Plus 2020, 135, 1–42. [Google Scholar] [CrossRef]
  12. Naik, P.A.; Owolabi, K.M.; Yavuz, M.; Zu, J. Chaotic dynamics of a fractional order HIV-1 model involving AIDS-related cancer cells. Chaos Solitons Fractals 2020, 140, 110272. [Google Scholar] [CrossRef]
  13. Naik, P.A.; Yavuz, M.; Zu, J. The role of prostitution on HIV transmission with memory: A modeling approach. Alex. Eng. J. 2020, 59, 2513–2531. [Google Scholar] [CrossRef]
  14. Sharma, S.; Samanta, G.P. Analysis of a drinking epidemic model. Int. J. Dyn. Control 2015, 3, 288–305. [Google Scholar] [CrossRef]
  15. Zeb, A.; Bano, A.; Alzahrani, E.; Zaman, G. Dynamical analysis of cigarette smoking model with a saturated incidence rate. AIP Adv. 2018, 8, 045317. [Google Scholar] [CrossRef]
  16. Yang, L.-X.; Yang, X.; Liu, J.; Zhu, Q.; Gan, C. Epidemics of computer viruses: A complex-network approach. Appl. Math. Comput. 2013, 219, 8705–8717. [Google Scholar] [CrossRef]
  17. Cheng, J.-J.; Liu, Y.; Shen, B.; Yuan, W.-G. An epidemic model of rumor diffusion in online social networks. Eur. Phys. J. B 2013, 86, 29. [Google Scholar] [CrossRef]
  18. White, E.; Comiskey, C. Heroin epidemics, treatment and ODE modelling. Math. Biosci. 2007, 208, 312–324. [Google Scholar] [CrossRef] [PubMed]
  19. Mulone, G.; Straughan, B. A note on heroin epidemics. Math. Biosci. 2009, 218, 138–141. [Google Scholar] [CrossRef]
  20. Samanta, G.P. Dynamic behaviour for a nonautonomous heroin epidemic model with time delay. J. Appl. Math. Comput. 2009, 35, 161–178. [Google Scholar] [CrossRef]
  21. Liu, J.; Zhang, T. Global behaviour of a heroin epidemic model with distributed delays. Appl. Math. Lett. 2011, 24, 1685–1692. [Google Scholar] [CrossRef] [Green Version]
  22. Huang, G.; Liu, A.P. A note on global stability for a heroin epidemic model with distributed delay. Appl. Math. Lett. 2013, 26, 687–691. [Google Scholar] [CrossRef]
  23. Fang, B.; Li, X.-Z.; Martcheva, M.; Cai, L.-M. Global stability for a heroin model with two distributed delays. Discret. Contin. Dyn. Syst. B 2014, 19, 715–733. [Google Scholar] [CrossRef]
  24. Fang, B.; Li, X.-Z.; Martcheva, M.; Cai, L.-M. Global asymptotic properties of a heroin epidemic model with treat-age. Appl. Math. Comput. 2015, 263, 315–331. [Google Scholar] [CrossRef]
  25. Fang, B.; Li, X.; Martcheva, M.; Cai, L. Global stability for a heroin model with age-dependent susceptibility. J. Syst. Sci. Complex. 2015, 28, 1243–1257. [Google Scholar] [CrossRef]
  26. Yang, J.Y.; Li, X.X.; Zhang, F.Q. Global dynamics of a heroin epidemic model with age structure and nonlinear incidence. Int. J. Biomath. 2016, 9, 1650033. [Google Scholar] [CrossRef]
  27. Djilali, S.; Touaoula, T.M.; Miri, S.E.-H. A Heroin Epidemic Model: Very General Non Linear Incidence, Treat-Age, and Global Stability. Acta Appl. Math. 2017, 152, 171–194. [Google Scholar] [CrossRef]
  28. Liu, L.; Liu, X. Mathematical Analysis for an Age-Structured Heroin Epidemic Model. Acta Appl. Math. 2019, 164, 193–217. [Google Scholar] [CrossRef]
  29. Duan, X.; Li, X.-Z.; Martcheva, M. Qualitative analysis on a diffusive age-structured heroin transmission model. Nonlinear Anal. Real World Appl. 2020, 54, 103105. [Google Scholar] [CrossRef]
  30. Liu, X.; Wang, J. Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate. J. Nonlinear Sci. Appl. 2016, 9, 2149–2160. [Google Scholar] [CrossRef]
  31. Yang, J.; Wang, L.; Li, X.; Zhang, F. Global dynamical analysis of a heroin epidemic model on complex networks. J. Appl. Anal. Comput. 2016, 6, 429–442. [Google Scholar]
  32. Liu, L.; Liu, X.; Wang, J. Threshold dynamics of a delayed multi-group heroin epidemic model in heterogeneous populations. Discret. Contin. Dyn. Syst. Ser. B 2016, 21, 2615–2630. [Google Scholar] [CrossRef] [Green Version]
  33. Wang, J.L.; Wang, J.; Kuniya, T. Analysis of an age-structured multi-group heroin epidemic model. Appl. Math. Comput. 2019, 347, 78–100. [Google Scholar] [CrossRef]
  34. Li, G.; Yang, Q.; Wei, Y. Dynamics of stochastic heroin epidemic model with lévy jumps. J. Appl. Anal. Comput. 2018, 8, 998–1010. [Google Scholar]
  35. Liu, S.; Zhang, L.; Xing, Y. Dynamics of a stochastic heroin epidemic model. J. Comput. Appl. Math. 2019, 351, 260–269. [Google Scholar] [CrossRef]
  36. Liu, S.; Zhang, L.; Zhang, X.; Li, A. Dynamics of a stochastic heroin epidemic model with bilinear incidence and varying population size. Int. J. Biomath. 2019, 12. [Google Scholar] [CrossRef]
  37. Wei, Y.; Yang, Q.; Li, G. Dynamics of the stochastically perturbed Heroin epidemic model under non-degenerate noises. Phys. A Stat. Mech. Appl. 2019, 526, 120914. [Google Scholar] [CrossRef]
  38. Rafiq, M.; Raza, A.; Iqbal, M.U.; Butt, Z.; Naseem, H.A.; Akram, M.A.; Butt, M.K.; Khaliq, A.; Ain, -Q.-U.; Azam, S. Numerical treatment of stochastic heroin epidemic model. Adv. Differ. Equ. 2019, 2019, 1–17. [Google Scholar] [CrossRef] [Green Version]
  39. Nyabadza, F.; Hove-Musekwa, S.D. From heroin epidemics to methamphetamine epidemics: Modelling substance abuse in a South African province. Math. Biosci. 2010, 225, 132–140. [Google Scholar] [CrossRef] [PubMed]
  40. Kalula, A.S.; Nyabadza, F. A theoretical model for substance abuse in the presence of treatment. S. Afr. J. Sci. 2012, 108, 12. [Google Scholar]
  41. Nyabadza, F.; Njagarah, J.B.H.; Smith, S.R. Modelling the Dynamics of Crystal Meth (‘Tik’) Abuse in the Presence of Drug-Supply Chains in South Africa. Bull. Math. Biol. 2012, 75, 24–48. [Google Scholar] [CrossRef]
  42. Mushanyu, J.; Nyabadza, F.; Stewart, A.G.R. Modelling the trends of inpatient and outpatient rehabilitation for methamphetamine in the Western Cape province of South Africa. BMC Res. Notes 2015, 8, 797. [Google Scholar] [CrossRef] [Green Version]
  43. Mushanyu, J.; Nyabadza, F.; Muchatibaya, G.; Stewart, A.G.R. Modelling Drug Abuse Epidemics in the Presence of Limited Rehabilitation Capacity. Bull. Math. Biol. 2016, 78, 2364–2389. [Google Scholar] [CrossRef]
  44. Mushanyu, J.; Nyabadza, F.; Muchatibaya, G.; Stewart, A.G.R. On the Role of Imitation on Adolescence Methamphetamine Abuse Dynamics. Acta Biotheor. 2017, 65, 37–61. [Google Scholar] [CrossRef]
  45. Su, S.; Fairley, C.K.; Mao, L.; Medland, N.A.; Jing, J.; Cheng, F.; Zhang, L. Estimates of the national trend of drugs use during 2000–2030 in China: A population-based mathematical model. Addict. Behav. 2019, 93, 65–71. [Google Scholar] [CrossRef] [PubMed]
  46. Muroya, Y.; Li, H.; Kuniya, T. Complete global analysis of an SIRS epidemic model with graded cure and incomplete recovery rates. J. Math. Anal. Appl. 2014, 410, 719–732. [Google Scholar] [CrossRef]
  47. Abdurahman, X.; Zhang, L.; Teng, Z. Global Dynamics of a Discretized Heroin Epidemic Model with Time Delay. Abstr. Appl. Anal. 2014, 2014, 1–10. [Google Scholar] [CrossRef] [Green Version]
  48. Wangari, I.M.; Stone, L. Analysis of a Heroin Epidemic Model with Saturated Treatment Function. J. Appl. Math. 2017, 2017, 1–21. [Google Scholar] [CrossRef] [Green Version]
  49. Duan, X.; Li, X.-Z.; Martcheva, M. Dynamics of an age-structured heroin transmission model with vaccination and treatment. Math. Biosci. Eng. 2019, 16, 397–420. [Google Scholar] [CrossRef]
  50. Memarbashi, R.; Pourhossieni, M. Global dynamic of a heroin epidemic model. UPB Sci. Bull. Ser. A 2019, 81, 115–126. [Google Scholar]
  51. Abdurahman, X.; Teng, Z.; Zhang, L. Global dynamics in a heroin epidemic model with different conscious stages and two distributed delays. Int. J. Biomath. 2019, 12. [Google Scholar] [CrossRef]
  52. Ma, M.; Liu, S.; Xiang, H.; Li, J. Dynamics of synthetic drugs transmission model with psychological addicts and general incidence rate. Phys. A Stat. Mech. Appl. 2018, 491, 641–649. [Google Scholar] [CrossRef]
  53. Naowarat, S.; Kumat, N. The Role of Family on the Transmission Model of Methamphetamine. J. Phys. Conf. Ser. 2018, 1039, 012036. [Google Scholar] [CrossRef]
  54. Saha, S.; Samanta, G.P. Synthetic drugs transmission: Stability analysis and optimal control. Lett. Biomath. 2019, 6, 1–31. [Google Scholar] [CrossRef]
  55. Liu, P.; Zhang, L.; Xing, Y. Modelling and stability of a synthetic drugs transmission model with relapse and treatment. J. Appl. Math. Comput. 2019, 60, 465–484. [Google Scholar] [CrossRef]
  56. Zhang, Z.; Yang, F.; Xia, W. Hopf Bifurcation Analysis of a Synthetic Drug Transmission Model with Time Delays. Complexity 2019, 2019, 1–17. [Google Scholar] [CrossRef] [Green Version]
  57. Li, J.; Ma, M.J. The analysis of a drug transmission model with family education and public health education. Infect. Dis. Model. 2018, 3, 74–84. [Google Scholar] [CrossRef] [PubMed]
  58. European Monitoring Centre for Drugs and Drug Addiction (EMCDDA). Population Surveys Methodology. In Methods and Definitions, General Population Surveys, Detailed Notes; EMCDDA Statistical Bulletin; EMCDDA: Lisbon, Portugal, 2004. [Google Scholar]
  59. Ministry of Public Security. Methods for the Identification of Drug Addiction, 2010. Chin. J. Drug Abuse Prev. Treat. 2011, 17, 63–64. [Google Scholar]
  60. Hser, Y.-I.; Fu, L.; Wu, F.; Du, J.; Zhao, M. Pilot trial of a recovery management intervention for heroin addicts released from compulsory rehabilitation in China. J. Subst. Abus. Treat. 2012, 44, 78–83. [Google Scholar] [CrossRef] [Green Version]
  61. Marienfeld, C.; Liu, P.; Wang, X.; Schottenfeld, R.; Zhou, W.; Chawarski, M. Evaluation of an implementation of methadone maintenance treatment in China. Drug Alcohol Depend. 2015, 157, 60–67. [Google Scholar] [CrossRef] [Green Version]
  62. Zhang, L.; Yap, L.; Zhuang, X.; Wu, Z.; Wilson, D.P. Needle and syringe programs in Yunnan, China yield health and financial return. BMC Public Health 2011, 11, 250. [Google Scholar] [CrossRef] [Green Version]
  63. Jia, Z.; Liu, Z.; Chu, P.; McGoogan, J.M.; Cong, M.; Shi, J.; Lu, L. Tracking the evolution of drug abuse in China, 2003–2010: A retrospective, Self-Controlled study. Addiction 2015, 110, 4–10. [Google Scholar] [CrossRef] [Green Version]
  64. Hale, J.K.; Lunel, S.M.V. Introduction to Functional Differential Equations; Springer: New York, NY, USA, 1993. [Google Scholar]
  65. Driessche, P.V.D.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 2002, 180, 29–48. [Google Scholar] [CrossRef]
  66. Arriola, L.; Hyman, J. Lecture notes, forward and adjoint sensitivity analysis: With applications in Dynamical Systems. In Linear Algebra and Optimisation Mathematical and Theoretical Biology Institute; Arizona State University: Tempe, AZ, USA, 2005. [Google Scholar]
  67. National Narcotics Control Committee. Report on Drug Control in China (2014), Beijing. 2014. Available online: http://www.nncc626.com/2014-09/12/c_126979288.htm (accessed on 27 October 2020).
  68. National Narcotics Control Committee. Drug Situation in China (2014), Beijing. 2015. Available online: http://www.nncc626.com/2015-06/24/c_127945747_2.htm (accessed on 27 October 2020).
  69. Du, X. Analysis and considerations of the current detoxification method in China. Chin. J. Drug Depend. 2005, 14, 392–398. [Google Scholar]
  70. Li, B.; Li, C.; Shi, E. On enforcement of drug trafficking prohibition between the western China and neighboring countries. J. Fujian Public Saf. Coll. 2004, 1, 14–17. [Google Scholar]
  71. National Bureau of Statistics. National Data, Beijing. 2020. Available online: https://data.stats.gov.cn/easyquery.htm?cn=C01 (accessed on 27 October 2020).
  72. Wang, X.; Yang, J.; Li, X. Dynamics of a Heroin Epidemic Model with Very Population. Appl. Math. 2011, 2, 732–738. [Google Scholar] [CrossRef] [Green Version]
  73. Li, X.; Zhou, Y.; Stanton, B. Illicit drug initiation among institutionalized drug users in China. Addiction 2002, 97, 575–582. [Google Scholar] [CrossRef] [PubMed]
  74. Yang, L.; Li, J.; Zhang, Y.; Zhang, W.; Dai, F.; Ren, Z.; Maycock, B. Reported Reasons for Initiating Drug Use among Drug-Dependent Adolescents and Youths in Yunnan, China. Am. J. Drug Alcohol Abus. 2009, 35, 445–453. [Google Scholar] [CrossRef] [PubMed]
  75. Yuan, Y.J.; Wang, Y. A review of the history of drug control publicity and education in China. Chin. J. Drug Depend. 2019, 28, 267–275. [Google Scholar]
  76. Yang, M.; Mamy, J.; Gao, P.; Xiao, S. From Abstinence to Relapse: A Preliminary Qualitative Study of Drug Users in a Compulsory Drug Rehabilitation Center in Changsha, China. PLoS ONE 2015, 10, e0130711. [Google Scholar] [CrossRef]
  77. Li, M.-T.; Jin, Z.; Sun, G.-Q.; Zhang, J. Modeling direct and indirect disease transmission using multi-group model. J. Math. Anal. Appl. 2017, 446, 1292–1309. [Google Scholar] [CrossRef]
  78. Castillo-Chávez, C.; Song, B. Dynamical Models of Tuberculosis and Their Applications. Math. Biosci. Eng. 2004, 1, 361–404. [Google Scholar] [CrossRef]
Figure 1. Flow diagram of the drug epidemic model.
Figure 1. Flow diagram of the drug epidemic model.
Ijerph 18 00288 g001
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 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).
Ijerph 18 00288 g002
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).
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).
Ijerph 18 00288 g003
Figure 4. 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 α (transfer rate from hidden drug addicts to drug addicts undergoing treatment).
Figure 4. 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 α (transfer rate from hidden drug addicts to drug addicts undergoing treatment).
Ijerph 18 00288 g004
Figure 5. Numerical simulations with five sets of initial conditions with respect to a drug-free equilibrium, including 3D trajectories in the I 1 I 2 I 3 space (a) and long-term time-series plots of drug addicts undergoing treatment (b), susceptible individuals (c), light drug users (d), hidden drug addicts (e), and the recovered individuals (f).
Figure 5. Numerical simulations with five sets of initial conditions with respect to a drug-free equilibrium, including 3D trajectories in the I 1 I 2 I 3 space (a) and long-term time-series plots of drug addicts undergoing treatment (b), susceptible individuals (c), light drug users (d), hidden drug addicts (e), and the recovered individuals (f).
Ijerph 18 00288 g005
Figure 6. Numerical simulations with five sets of initial conditions with respect to a drug-persistent equilibrium, including 3D trajectories in the I 1 I 2 I 3 space (a) and long-term time-series plots of drug addicts undergoing treatment (b), susceptible individuals (c), light drug users (d), hidden drug addicts (e), and the recovered individuals (f).
Figure 6. Numerical simulations with five sets of initial conditions with respect to a drug-persistent equilibrium, including 3D trajectories in the I 1 I 2 I 3 space (a) and long-term time-series plots of drug addicts undergoing treatment (b), susceptible individuals (c), light drug users (d), hidden drug addicts (e), and the recovered individuals (f).
Ijerph 18 00288 g006
Figure 7. Model fitting and projection results, including iteration procedures of curve fitting (a) and time-series plots of drug addicts undergoing treatment (b), sum of all compartments (c), light drug users (d), hidden drug addicts (e), and recovered individuals (f).
Figure 7. Model fitting and projection results, including iteration procedures of curve fitting (a) and time-series plots of drug addicts undergoing treatment (b), sum of all compartments (c), light drug users (d), hidden drug addicts (e), and recovered individuals (f).
Ijerph 18 00288 g007
Figure 8. The effect of Intervention 1 (anti-drug education and propaganda, 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 the recovered individuals (d).
Figure 8. The effect of Intervention 1 (anti-drug education and propaganda, 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 the recovered individuals (d).
Ijerph 18 00288 g008
Figure 9. The effect of Intervention 2 (regional anti-drug propaganda, 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 the recovered individuals (d).
Figure 9. The effect of Intervention 2 (regional anti-drug propaganda, 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 the recovered individuals (d).
Ijerph 18 00288 g009
Figure 10. The effect of Intervention 3 (investigation and admission, 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).
Figure 10. The effect of Intervention 3 (investigation and admission, 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).
Ijerph 18 00288 g010
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).
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).
Ijerph 18 00288 g011
Table 1. Descriptions of model variables and parameters.
Table 1. Descriptions of model variables and parameters.
Variable/
Parameter
Description
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
Table 2. Historical data of population aged 15−64, number of existing drug users, and former drug users who had remained abstinent for three or more years during 2015−2019.
Table 2. Historical data of population aged 15−64, number of existing drug users, and former drug users who had remained abstinent for three or more years during 2015−2019.
YearPopulation Aged 15−64 *Existing Drug Users *Former Drug Users Abstinent for 3 Years *
2015100,361234.5114.8
2016100,260250.5141.1
201799,829255.3167.9
201899,357240.4207.3
201998,914214.8253.3
* The unit of all numbers is 10,000 people.
Table 3. Parameter units, value ranges, final values used, and sources.
Table 3. Parameter units, value ranges, final values used, and sources.
ParameterUnitRangeValueSource
λ 10 thousand people/year(235, 610)400Estimated from [71]
μ /year(0.0064, 0.00716)0.007[71]
μ d /year(0.021, 0.102)0.025[39,44]
β 1 /10 thousand people*year(1 × 10−9, 1 × 10−6)1.2481 × 10−7Curve fit
β 2 /10 thousand people*year(1 × 10−9, 1 × 10−6)3.8611 × 10−7Curve fit
k 1 /year(0.05, 0.3)0.176Curve fit
k 2 /year(0.05, 0.5)0.2[23,25,72]
α /year(0.05, 0.6)0.124Curve fit
r /year(0.33,0.6)0.45Estimated from [60]
Table 4. Comparison of Interventions 1−4 and their corresponding basic reproduction numbers.
Table 4. Comparison of Interventions 1−4 and their corresponding basic reproduction numbers.
OriginalIntervention 1Intervention 2Intervention 3Intervention 4
2020 *2030 *2030 *Increase (%)2030 *Increase (%)2030 *Increase (%)2030 *Increase (%)
S 97,368.4894,533.2194,592.850.06%94,572.980.04%94,566.470.04%94,552.610.02%
I 1 77.4223.2811.81−49.27%15.55−33.20%6.71−71.18%12.31−47.12%
I 2 206.8665.7057.30−12.79%60.10−8.52%54.65−16.82%62.56−4.78%
I 3 496.84143.57129.81−9.58%134.40−6.39%41.00−71.44%76.58−46.66%
R 639.971140.081116.05−2.11%1124.13−1.40%1261.7110.67%1216.286.68%
R 0 0.0872560.0436280.0584620.0428530.057216
* The unit of all numbers of compartments is 10,000 people.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tang, H.; Li, M.; Yan, X.; Lu, Z.; Jia, Z. Modeling the Dynamics of Drug Spreading in China. Int. J. Environ. Res. Public Health 2021, 18, 288. https://doi.org/10.3390/ijerph18010288

AMA Style

Tang H, Li M, Yan X, Lu Z, Jia Z. Modeling the Dynamics of Drug Spreading in China. International Journal of Environmental Research and Public Health. 2021; 18(1):288. https://doi.org/10.3390/ijerph18010288

Chicago/Turabian Style

Tang, Haoxiang, Mingtao Li, Xiangyu Yan, Zuhong Lu, and Zhongwei Jia. 2021. "Modeling the Dynamics of Drug Spreading in China" International Journal of Environmental Research and Public Health 18, no. 1: 288. https://doi.org/10.3390/ijerph18010288

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop