Role of imitation and limited rehabilitation capacity on the spread of drug abuse.

OBJECTIVES: We formulate a mathematical model for the spread of drug abuse using non linear ordinary differential equations. The model seeks to investigate both peer influence and limited rehabilitation effects on the dynamics of drug abuse. Peer-influence is modelled through the mechanism of imitation and limited rehabilitation is described using a special treatment function. Center manifold theory is used to show that the model exhibits the phenomenon of backward bifurcation. Matlab has been used to carry out numerical simulations to support theoretical findings.
RESULTS: The model analysis shows that the model has multiple equilibria. It has been shown that the classical [Formula: see text]-threshold is not the key to control drug abuse within a population. In fact drug abuse problems may persist in the population even with subthreshold values of [Formula: see text]. This was shown to result, in particular when, [Formula: see text], [Formula: see text] and [Formula: see text] are high enough such that [Formula: see text], [Formula: see text] and [Formula: see text]. The results suggest the need for comprehensive and accessible substance abuse treatment services to curtail drug abuse.

Introduction

Drug abuse has increased in recent years and is now an epidemic globally. The magnitude of the world drug problem becomes more apparent when considering that more than 1 out of 10 drug users is a problem drug user and the vast majority of these individuals continue to have no access to treatment [1]. There continues to be a large “treatment gap” for substance abuse problems as many countries have a large shortfall in the provision of services. According to the United Nations Office on Drugs and Crime [1], only one out of every six problem drug users in the world has access to treatment. Generally, the number of patients in need of rehabilitation often exceeds the carrying capacities of drug treatment facilities, especially those funded by the state.

Several mathematical models describing the spread of psycho-social ills in a community have been proposed, see for example, drug epidemics [2-9], alcoholism [10-16], smoking [17-19]. The basic assumption in most drug abuse models is that there is a direct proportional relationship between the number of drug users in need of treatment and the available health care resources present. In this paper, we develop a mathematical model that takes into account the possibility of the number of drug abusers in need of rehabilitation exceeding the capacity of rehabilitation centers. Recruitment into rehabilitation (inpatient or outpatient) is denoted by H(U) and defined as follows:where U represents the proportion of individuals abusing drugs, is the maximum rehabilitation uptake per unit of time and measures the extent of the effect of the problem of demand for treatment. Firstly, observe that for small U, , that is, when the number of drug users is not too large, then the rate of entering treatment is proportional to the number of drug users present. Secondly, observe that for large U, , this means that the rate of entering rehabilitation takes a maximum value . Finally, when , we again obtain the result as in the first case, , that is, the function returns to a linear function mostly used in previous drug abuse models. Amongst drug abusers who are seeking help through rehabilitation, we have that a proportion p of these individuals are recruited into inpatient rehabilitation and the complementary proportion are recruited into outpatient rehabilitation. It is also important to note that epidemic models including treatment functions of the form (1) are found in [20-23].

We also include peer influence effects on the spread of drug abuse by assuming that the recruitment process happens through the mechanism of imitation. In this paper, we use the recruitment function given in [11]. Compared to previous drug epidemic models [2-9], a key novelty of our model is the inclusion of both imitation and limited rehabilitation on the dynamics of drug abuse.

The paper is arranged as follows; in “Model formulation” section, we formulate and establish the basic properties of the model. The model is analysed for stability in “Model analysis” section. In “Numerical simulations” section, we carry out some numerical simulations. Parameter estimation is also presented in this section. The paper is concluded in the “Conclusions” section.

Main text

Model formulation

The model divides the population into four classes, S(t), U(t), and . The class S(t) represents the population at risk of being initiated into drug abuse. The class U(t) denotes those abusing drugs, denotes those in rehabilitation as out-patients and denotes those in rehabilitation as in-patients. The total local population is thus given byThe general population enter the susceptible population at a rate , that is, the demographic process of individuals reaching age 15 in the modelling time period. Susceptible individuals become drug users upon contact with individuals in compartments U or . This results from the assumption that those in inpatient rehabilitation do not have contact with the population at risk. The per capita contact rate is a product of the effective number of contacts , between drug users not in rehabilitation and the susceptible population, and the probability , that a contact results into initiation into drug use, that is, . The per capita contact rate is a product of the effective number of contacts , between drug users in outpatient rehabilitation and the susceptible population, and the probability , that a contact results into initiation into drug use, that is, . Individuals under outpatient rehabilitation quit drug abuse permanently at a rate and individuals under inpatient rehabilitation quit drug abuse permanently at a rate . The general population experience natural death at a rate . Drug users undergoing outpatient rehabilitation relapse into drug use at a rate whereas those undergoing inpatient rehabilitation relapse at a rate . The relapse is thus assumed to be a voluntary process, that is not influenced by interaction with users. We allow the transfer from outpatient to inpatient rehabilitation, this happens at a rate . We also allow the transfer from inpatient to outpatient rehabilitation, this rate is represented by . We assume that individuals in each compartment are indistinguishable and there is homogeneous mixing. We have the following general set of nonlinear ordinary differential equations:with the initial conditions:whereHere , with signifying that the chance of initiating drug abuse habit upon contact with an individual in U or is the same, signifying a reduced chance of initiating drug abuse habit upon contact with an individual in as compared to an individual in U, signifies an increased rate of initiating drug abuse habit upon contact with an individual in as compared to an individual in U.

Model analysis

Model properties

Invariant region

It follows from system (2) thatThen, . Thus, the feasible region for system (2) isIt is easy to verify that the region is positively invariant with respect to system (2), see for instance [3-5].

The drug-free equilibrium and the abuse reproduction number

Model system (2) always has a drug-free equilibrium . Denote the abuse reproduction number of model system (2) bywithWe can clearly note that and so . Also, and . Therefore, is non-negative. The abuse reproduction number of the model, is the average number of secondary cases generated by one drug user during his/her duration of drug use in a population of completely potential drug users.

Local stability of the drug-free steady state

Theorem 1

The drug-free equilibrium is locally asymptotically stable when and is unstable when .

Proof

The Jacobian matrix of model system Eq. (2) at is given bywhere and are defined as before and , . The local stability of the drug-free equilibrium is determined by the following submatrix of ,Since all off-diagonal elements are positive, we now consider matrix . We claim that is an M—matrix. Multiplying matrix by the positive matrixwe havewhere is a positive matrix given byThen, it follows from M—matrix theory that all eigenvalues of have negative real parts, which implies the local asymptotic stability of the drug-free equilibrium if . On the other hand, it can be shown that the determinant of is given byThus, if , then matrix has eigenvalues with negative real parts, which implies the stability of the drug-free equilibrium. This completes the proof.

The drug-persistent equilibrium point

The drug-persistent equilibrium always satisfiesFrom the last two equations of (5) we have thatwhereSubstituting expressions (6) into the first equation of (5), we getSubstituting expressions (6) and (7) into the second equation of (5) leads to the following sixth order polynomial equationSolving (8) gives which corresponds to the drug-free equilibrium orwhere the coefficients are in Additional file 1: Appendix S1. We can clearly note that, and . The number of possible positive real roots of the polynomial (9) can be determined using the Descartes Rule of Signs. The number of positive roots are at most five.

Backward bifurcation

Conditions for the existence of backward bifurcation follow from Theorem 4.1 proven in [24]. Let us make the following change of variables:

, so that . We now use the vector notation . System (2) can be written in the form , wherewithLet be the bifurcation parameter, corresponds toThe Jacobian matrix of system (2) at when is given bywhere and are defined as before and , .

System (10), with has a simple eigenvalue, hence the center manifold theory can be used to analyse the dynamics of system (2) near . It can be shown that , has a right eigenvector given by , whereFurther, the left eigenvector of , associated with the zero eigenvalue at is given by , whereThe computations of a and b are necessary in order to apply Theorem 4.1 in Castillo-Chavez and Song [24]. For system (10), the associated non-zero partial derivatives of F at the drug-free equilibrium are in Additional file 1: Appendix S2. It thus follows thatwherewithNote that , and . Also note that if , and then and if , and . Lastly,We thus have the following result

Theorem 2

If , and , then model system (2) has a backward bifurcation at .

Results and discussion

Numerical simulations

Parameter estimation

Since we can rarely enumerate the incidence of drug users, data from treatment centers can be used as proxy for estimating parameters for drug related issues. We use data obtained from previous mathematical models with inpatient and outpatient rehabilitation [4, 5]. Some of the parameter values will be obtained from literature.

Parameter values used for numerical simulations are given in Table 1.

Table 1

Parameter values used in numerical simulations

Parameter	Range	Value	Source
β_1	0.10–0.21	0.105	[7]
β_2	0–0.10	0.063	[6]
ω	0–1	0.62	[5]
α	0–0.05024	0.02827	[4]
p	0–1	0.352	[4]
η_1	0–1	0.24	Assumed
η_2	0–1	0.13	Assumed
δ_1	0.001–1	0.01	[4]
δ_2	0.01–1	0.3142	[4]
ρ_1	0–0.054	0.0382	[4]
ρ_2	0–0.0235	0.0020	[4]
γ_1	0–0.06012	0.02961	[4]
γ_2	0–0.008	0.003	[4]
Λ	0.028–0.080	0.04	[7]
μ	0.019–0.021	0.020	[25]

Numerical results

We carry out detailed numerical simulations using matlab to support our theoretical findings. The initial conditions used are: , , , .

Figures 1 and 2 illustrate the effect of varying parameters and on the prevalence of drug abuse. Figures 1 and 2 demonstrate that increasing and results in an increase in the prevalence of drug abuse. This is a reflection that limited rehabilitation and imitation are of major concern in the fight against drug abuse.

Fig. 1

Effects of varying on the prevalence of drug abuse, starting from 0 up to 1.0 with a step size of 0.5

Fig. 2

Effects of varying on the prevalence of drug abuse, starting from 0 up to 1.0 with a step size of 0.5

Conclusions

A mathematical model that incorporates imitation and limited rehabilitation has been formulated using nonlinear ordinary differential equations. It has been shown that the classical —threshold is not the key to control drug abuse within a population. In fact drug abuse problems may persist in the population even with subthreshold values of . This was shown to result, in particular when , and are high enough such that , and . Considerable effort should be directed towards reducing , and , this done by increasing the value of , and so as to avoid backward bifurcation. Also, results from the model application show that increasing and lead to an increase in the prevalence of drug abuse. Thus, communities should have suitable capacity for the treatment of drug abusers and specific health and/or education programs may be employed to reduce the imitation coefficient .

Limitations

Like in any model development, the model is not without limitations.

The model did not take into account contextual dynamics, such as drug supply chains or changes in interdiction.

Also, the study presented here ignored detailed social and economic characteristics.

Other initiation processes, not included in this work, for instance, initiation by self-conversion, drug supply chains etc. may form part of the author’s future research considerations.

Additional file 1. Appendices S1, S2.