
Abstract
Introduction: what only mathematical models can reveal
An operational model of problematic drug use epidemics
The crucial parameters to be estimated externally and the scenario parameters
Some scenario analyses
Further developments
References
ABSTRACT
A modified version of a model for the human immunodeficiency virus/acquired immunodeficiency syndrome (HIV/AIDS) epidemic, proposed at the beginning of the 1990s and recently generalized, is presented as a representation of a problematic drug use epidemic. The model can be used both to estimate relevant epidemic macroparameters and to carry out scenario analyses. The model is of the “moverstayer” type and allows for heterogeneous risk behaviour among the susceptible population. Such models treat the “susceptibles” as subdivided into two groups: the “stayers”, that is, the individuals who are considered not to be at risk of infection, and the “movers”, who are at risk. Because of the interactions between infectious individuals (pushers) and the susceptibles, some of the latter may pass into the drug user “compartments” and begin a drug user “career”. The model described here comprises two different stages of hidden drug use. The first (light use) stage, which can be more strictly defined, is the initial (or nonproblematic) stage of drug use, following which light drug users can either stop using drugs or pass to hard drug use (or death). Other relevant stages taken into account are the therapy stage, the recidivist use stage (either light or hard use) and the temporary nonuse stage (the latter are the visible stages). A simpler version of the model is studied using a Markov hybrid approximation and some “what if” scenario analyses are obtained by simulation. A more complex and realistic model is also outlined for possible further development.
The drug problem and its consequences for society represent a complicated field for research. Policy makers and researchers are seeking answers to a number of questions concerning drug use, its consequences and related costs. They are examining the extent of drug use, the consequences of drug use, which policies are effective, the costs of such policies and the cost of drug use to society. The con sequences comprise, among many others, the adverse effects of infectious diseases and the costs to society of drugrelated criminality. It is thus important to understand and measure drug use and how it responds to drug control interventions [1]. The present article will contribute to that objective by introducing a simple compartmental epidemic model of drug use that incorporates various effects on initiation into new use or relapse into recidivist use of drugs.
Compartmental models represent a powerful mathematical tool that is well established in modelling the spread of diseases in a population. They provide a framework in which numbers of people in different “compartments” (each one homogeneous with respect to some specified characteristic) and the relationships between such compartments, modelling the dynamics of the population, can be described in mathematical terms.
Dynamic compartmental modelling of epidemic processes, producing either operational or transmission models, occurs through the usual representation of the dynamics involved by means of a system of stochastic or deterministic differential (or difference) equations; this is the case for both the operational and the transmission models. The main difference between the two kinds of model lies in the fact that transmission models take into account the dynamic processes at the microlevel, modelling the interactions between individuals belonging to the different subgroups involved in the epidemic, whereas operational models work on macrovariables or indicators suitable for use in estimating the size of the phenomena or monitoring the impact of various interventions, modelled by suitable scenario parameters. Many models of the two types have been developed to study the human immunodeficiency virus/acquired immunodeficiency syndrome (HIV/AIDS) epidemic and they can also be used, with some modifications, to represent problematic drug use epidemics.
Some examples of transmission models are given in Dietz [2], Hadeler [3], Kretzschmar and Dietz [4] in relation to the HIV/AIDS epidemic and in Billard and Dayananda [5] and in Behrens [1] for problematic drug use epidemics. Those models are complex because the researchers introduce into the modelling process a detailed formalization of the interactions existing, or assumed to exist, between a large number of subjects involved in the epidemic process. This is the case of the analysis of HIV transmission across risk groups, such as drug users and heterosexuals, when some specific hypothesis is formulated as regards the contact pattern [4, 6] or when the time at which the contact occurs is assumed somehow to affect the probability of transmission of the infection [4]. As already mentioned, such models can be helpful when specific contact or transmission patterns need to be thoroughly analysed, but they may be extremely cumbersome when the epidemic as a whole is being examined. Detailed model structures contain all sorts of parameters regulating every single interaction in the simulation process and seldom are enough data available for reliable estimation of all those many parameters.
A more efficient way of obtaining a picture of a simulated epidemic is by using a “simple” operational model. In contrast to transmission models, simple models do not attempt to include all the possible group or individual interactions into the modelling structure, but summarize the dynamics of the epidemic in terms of certain nonlinear interactions and group the infected individuals in chains of “compartments”. Most of the parameters controlling the dynamics in such systems are derived from epidemiological studies, external to the model, and their values are based simply on specialized studies or on the literature; only a limited number of “internal” parameters are usually left to be estimated by fitting the existing data or to be used for scenario analyses. Typically, the set of internal parameters includes some form of control of the transmission and of the size of the core group; other internal parameters may have different origins and interpretations, depending on the design of each individual model. The ability of simple models to describe an epidemic correctly is theoretically much more limited when compared with that of complex transmission models. In general, however, simple models turn out to be much more tractable, both because of the limited number of parameters required for their functioning and because of the quality of their output. In fact, while the large number of parameters in the transmission models often results in their unidentifiability, parameter estimates in simple models are usually reliable enough for the models to be used successfully in complex scenario analyses and under other difficult conditions.
Most of the simple models have a stochastic basis, although their output results in some deterministic forecast curves. Such models, defined as “hybrid” models in Bailey [79], are the result of a deterministic approximation of stochastic models, using mean times of stay in compartments with conditional Poisson arrivals and maximum likelihood (ML) or numerical estimates of internal para meters. Their use, however theoretically disputable, is more and more frequent among applied researchers because of their technical simplicity and the clarity and immediate possibility of use and interpretation of their results: their output consists of deterministic curves of epidemic indicators (incidence, prevalence, dynamics of the core group and so on), while the parameter ML estimation procedure is derived directly from the distributional hypotheses of the original stochastic model. One such model for the HIV/AIDS epidemic was proposed at the beginning of the 1990s and has recently been generalized. It makes it possible to obtain scenario analyses with ease [10, 11]. That model, when used in conjunction with suitable backcalculation methods, can reduce uncertainties in the estimates of incidence curves [12] and has also been used for indirectly estimating the prevalence of injecting drug use in Italy [13]. A modified version of the model is presented below for use in relation to epidemics of problematic drug use. (The model will be presented from a mathematical point of view elsewhere [14].)
The graphs provided in figure I describe the main features of the proposed model. The model is of the “moverstayer” type [10, 15, 16] and allows for heterogeneous risk behaviour among the susceptible population. Such models treat the susceptible population as subdivided into two groups: the “stayers”, that is the individuals who are considered not to be at risk of infection (such models are suit able for scenario analyses in order to assess the impact of various proportions of vaccinated persons on the probability of extinction of a given epidemic) and the “movers”, who are at risk. Because of the interactions between infectious individuals (for the present study these are taken to be problematic drug users who are also pushers)** and the susceptibles, or as a result of the pressure of the black market on the susceptibles, some of the latter may pass into the drug user “compartments” and begin a drug user “career”. As in the model proposed by Behrens and colleagues [1], the model described here comprises two different stages of hidden drug use. The first (light drug use) stage, which can be more strictly defined, is the initial (or nonproblematic) stage of drug use, following which light drug users can either stop using drugs or pass to hard drug use (or death). The other arrows in the diagram show all the other possible transitions in a drug user career. The curves connecting the drug use (infectious) compartments and the susceptible (or temporary nonuse) compartments indicate the possible interactions that may produce transitions from susceptibles (or temporary nonuse) to infectious (nonlinear terms in the equations); the other possible transitions from susceptibles (or temporary nonuse) to infectious are induced by the pressure of the black market and are assumed to be represented in linear terms in the equations.
In order to formulate the corresponding equations (either deterministic or stochastic), some further hypotheses must be explored and the known and unknown parameters described. A first approximation may be through a Markov model, possibly a marked Markov process. In such a case the length of stay in each compartment is assumed to be distributed exponentially and the results are useful in providing a first qualitative insight into the epidemic process. A more realistic approximation is by using semiMarkov processes. In such a case the length of stay in each compartment may be assumed to be distributed differently with respect to an exponential variable, possibly depending also on covariates describing the individual experience of drug use. Such a study is much more complex, but suitable mathematical and simulation techniques can be used. In the present article, only the Markov model is considered and used to make some scenario analyses.
On the basis of figure I, it is immediately possible to formulate the equations of the model either in the form of deterministic (continuous or discrete) equations or as stochastic (continuous or discrete) equations. The discrete stochastic equations will be reported elsewhere [14]; only the simulation programme and some scenarios aimed at evaluating the impact of different parameters on the epidemic are considered in the present article. The state variables used in the model (with the exception of S ( t ), which is the proportion of stayers at time t , represent the incidence per one million inhabitants; the unit time for the simulation runs is taken to equal one week. The equation for the proportion of stayers, S , is derived from the hypothesis [10] that the new entrances in the susceptible compartment are divided into stayers and movers according to constant proportions S _{o} and 1S _{o} (stationarity), with 0 < S _{o } < 1, even if other hypotheses can be incidentally incorporated into the model as well as possible transitions from the movers to the stayers compartments as a result of prevention campaigns or law enforcement activities. This further development is outlined briefly in the last part of this article.
A qualitative analysis of the epidemic (transient and asymptotic behaviours), in particular as regards the incidence of drug use, that is the transitions from susceptible to drug user, can be conducted by analysing jointly the equation for X and the equation for S . The method is analogous to that used in Rossi [10] and will be considered elsewhere [14].
As regards the parameters and the distributions of the length of stay, some are already available from study of the latency period [17]. Some can be derived using therapy data already available in various places. The demographic parameters regulating the dynamics of the susceptible population, µ_{01} , µ_{10} and π_{17} , are assumed to be known and are countryspecific. The other µ parameters can be estimated externally using the information from mortality studies among drug users, which is available for most countries of the European Union [18]. The parameters µ_{23} and µ_{34} (natural history parameters) can be estimated on the basis of data available from study of the latency period [19, 20]. The parameters µ_{45} , µ_{46} , µ_{54} , µ_{56} , µ_{61} and µ_{65} (therapy parameters) can be obtained, at least as regards their order of magnitude, from therapy data available in most countries; the values of all these parameters for Italy (order of magnitude) are given in table 1. All the other parameters, µ_{12} , µ_{26} , v_{12} , v_{13} , v_{15} , v_{26} , v_{36} and v_{56} , can be used as scenario parameters, as well as the parameter “initial proportion of stayers”, S _{o} .
In the present article it is assumed (following Billard and Dayananda [5]) that most pushers are hard drug users (basic vscenarios). It has been suggested elsewhere that most pushers are soft drug users [21], so some simulation runs will be devoted to obtaining scenarios based on this latter hypothesis (alternative vscenarios) by modifying the order of magnitude of the v parameters reported in table 1 (10^{5} instead of 10^{6} and vice versa).
All the µ and π parameters represent transition rates per person of the origin compartment per week, the n parameters are rates per week per pair. Using the simulation procedure described below, some impact analyses can be conducted to evaluate the influence of the scenario parameters on the course of the epidemic. It is also possible to make some further scenario analyses to evaluate the impact of increased efficacy of therapy services by taking the scenario parameters as fixed and using the therapy parameters as variable. Similarly, the natural history parameters can be used to make further “what if” scenario analyses (the results of such analyses will be presented elsewhere).
The simulation procedure, used to obtain scenario analyses, is written in Splus for personal computers. All the parameters can be modified at the beginning of each run. The standard output comprises graphs of the prevalence curves in each compartment and of the incidence curves of major interest. It is possible to choose the total simulation time, which is measured in weeks.
Some scenario analyses have been conducted to study the impact of S _{o} and some other parameters on the qualitative behaviour of the epidemic. The results are reported in graphic form below. As expected, the effect of S _{o} is higher with respect to the other parameters [11] considered for the simulation runs reported here.
Some scenarios have been obtained based on various hypotheses concerning the parameters. The results are reported in graphic form and a summary table (table 2) describes the macrocharacteristics of the different epidemics using the location and the size of the peaks for the incidence and prevalence curves, the value of the cumulative curve of deaths at the end of the period taken into account and the location of the maximum of the proportion of stayers, which represents the time of the beginning of the endemic phase (saturation effect).
Scenario 1 is obtained using the parameters listed in table 1, with S _{o} = 0.98 (this means that the size of the population at risk is taken to be about 20,000 per one million inhabitants), µ_{12} = 10^{5} (this parameter measures the strength of the black market: its value means that the black market can be assumed to induce about one person of the risk group per month to try drugs at the beginning of the epidemic) and µ_{26} = 0.004 (this is the parameter regulating spontaneous exits from the nonaddicted user population: its value means that the proportion of nonaddicts who stop using drugs is about half those who proceed to addictive behavour). Various graphs representing the curves of interest are shown in figures IIVII.

Figure III. Scenario 1: Prevalence curve of light users 
Figure IV. Scenario 1: Incidence curve from light users
to heavy users 
Figure V. Scenario 1: Prevalence curve of hard users 
Figure VI. Scenario 1: Prevalence curve of clients of
therapy services 
Figure VII. Scenario 1: Proportion of stayers 
This scenario is obtained using the parameters listed in table 1, with S _{o} = 0.99 (here the size of the population at risk is taken to be about 10,000 per one million inhabitants), µ_{12} = 10^{5} (as in scenario 1) and µ_{26} = 0.004 (as in scenario 1). Graphs representing the curves of interest are shown in figures VIIIXIII.
Figure VIII. Scenario 2: Incidence curve from susceptibles to light users 
Figure IX. Scenario 2: Prevalence curve of light users 
Figure X. Scenario 2: Incidence curve from light users to hard users 
Figure XI. Scenario 2: Prevalence curve of hard users 
Figure XII. Scenario 2: Prevalence curve of clients of therapy
services 
Figure XIII. Scenario 2: Proportion of stayers 
Scenario 3 is obtained using the parameters listed in table 1, with S _{o}_{ } = 0.98 (as in scenario 1), µ_{12} = 10^{6} (this parameter measures the strength of the black market: its value means that the black market can be assumed to induce about one person of the risk group per 10 months to try drugs at the beginning of the epidemic) and µ_{26} = 0.004 (as in scenario 2). Graphs representing the curves of interest are shown in figures XIVXIX.
Figure XIV. Scenario 3: Incidence curve from users 
Figure XV. Scenario 3: Prevalence curve of light susceptibles to
light users 
Figure XVI. Scenario 3: Incidence curve from light users 
Figure XVII. Scenario 3: Prevalence curve of clients of therapy
services 
Figure XVIII. Scenario 3: Prevalence curve of clients of therapy
services 
Figure XIX. Scenario 3: Proportion of stayers 
This scenario is obtained using the parameters listed in table 1, with S _{o} = 0.98 (as in scenario 3), µ_{12} = 10^{5} (as in scenarios 1 and 2) and µ_{26} = 0.0004 (this is the parameter regulating spontaneous exits from the nonaddicted user population: its value means that the proportion of nonaddicts who stop using drugs is about 1/20 of those who proceed to addictive behaviour). Graphs representing the curves of interest are shown in figures XXXXV.
Figure XX. Scenario 4: Incidence curve from susceptibles
to light users 
Figure XXI. Scenario 4: Prevalence curve of light users 
Figure XXII. Scenario 4: Incidence curve from light
users to hard users 
Figure XXIII. Scenario 4: Prevalence curve of hard users 
Figure XXIV. Scenario 4: Prevalence curve of clients of
therapy services 
Figure XXV. Scenario 4: Proportion of stayers 
The present scenario is obtained using the parameters listed in table 1, but modifying the order of magnitude of the n parameters (10^{5} instead of 10^{6} and vice versa). This means that it is assumed that there is higher possibility that a nonaddicted individual can induce a susceptible to try drugs than an addicted one can. The values S _{o} = 0.98, µ_{12} = 10^{5} and µ_{26} = 0.004 are the same as in scenario 1. Graphs representing the curves of interest are shown in figures XXVIXXXI.
Figure XXVI. Scenario 5: Incidence curve
from susceptibles to light users 
Figure XXVII. Scenario 5: Prevalence curve
of light users 
Figure XXVIII. Scenario 5: Incidence curve
from light users to hard users 
Figure XXIX. Scenario 5: Prevalence curve of
hard users 
Figure XXX. Scenario 5: Prevalence curve of
clients of therapy services 
Figure XXXI. Scenario 5: Proportion of
stayers 
The present scenario is obtained using the parameters listed in table 1, but modifying the order of magnitude of the v parameters (10^{5} instead of 10^{6} and vice versa), as in scenario 5, and with S _{o} = 0.99, µ_{12} = 10^{5} and µ_{26} = 0.004 (as in scenario 2). Graphs representing the curves of interest are shown in figures XXXIIXXXVII.
Figure XXXII. Scenario 6: Incidence curve from
susceptibles to light users 
Figure XXXIII. Scenario 6: Prevalence curve of light
users 
Figure XXXIV. Scenario 6: Incidence curve from from
light users to hard users 
Figure XXXV. Scenario 6: Prevalence curve of hard users 
Figure XXXVI. Scenario 6: Prevalence curve of clients of
therapy services 
Figure XXXVII. Scenario 6: Proportion of stayers 
An analysis of the graphs shows clearly that the parameter with the highest impact on the course of the epidemic is S _{o} , which is a measure of the size of the group of susceptibles who are at risk of infection (the core group). The bigger the core group, the faster the evolution of the epidemic and the higher the prevalence and incidence curves. The influence of the parameter that measures the pressure of the black market, µ_{12} , appears to be less important. As regards the comparison between the basic hypothesis for the n parameters, that is, that pushers are mainly hard users, and the alternative hypothesis, that pushers are mainly light users, on the basis of the results of the simulations, it appears that the alternative scenarios show faster evolutions of epidemics. Table 2 summarizes the results in terms of the macroindicators (location and size of the peaks for the incidence and prevalence curves, value of the cumulative curve of deaths at the end of the period taken into account and the location of the maximum of the proportion of stayers) described above. All the values related to size represent rates per one million inhabitants.
It can be observed that the initial proportion of stayers, S _{o} is the parameter with the highest impact: there are higher differences in the macroparameters between scenario 1 and scenario 2 than between scenarios 1 and 3 or between scenarios 1 and 4.
Higher values of S _{0} correspond to larger and fastergrowing epidemics, whereas the parameter related to the pressure of the market is not as important, at least in the range of values examined here. A minor effect also arises from parameter µ_{26} , which models transitions from light use to nonuse (spontaneous cessation). As mentioned above, alternative scenarios correspond essentially to fastergrowing epidemics, other things beings equal.
Further analyses will be reported elsewhere [14].
The use of suitable markers (marked processes) might make it possible to incorporate further descriptions of each individual involved in an epidemic; for in stance, markers may take into account the number of incarcerations, the number of nonfatal overdoses or the number of failed therapy interventions. Possible repeated therapy interventions should be incorporated into the model in order to obtain a more realistic picture of the career of the problematic drug user; there is general agreement that the time spent in the therapy compartment for the first therapy episode is distributed differently with respect to times related to the following episodes. The complexity of such analysis is evident and the interpretation of the results may be quite uncertain and unreliable owing to lack of information about the statistical distribution of the length of stay in the various compartments. It is therefore preferable to wait until the therapy data sets comprise reliable and complete information that will make possible external estimation of the parameters of interest.
In order to obtain more realistic models, the transition parameters should not be taken as constant, but should be represented as functions taking into account the history of drug use for a given individual, represented by statistical variables taken as known (covariates), and the history of policy interventions (availability of services, law enforcement activities and so on), represented by other covariates (timedependent) and, possibly, by latent variables. This would result in a realistic but at the same time unreliable and intractable transmission model of very limited use. Some years ago, Kaplan [22] published a paper entitled “Can bad models suggest good policies?” in which he clearly explained how simple models, taking into account only the main peculiarities of the phenomenon of interest and neglecting minor effects, which may simply mask the relevant behaviours, are much more useful in order to assess the main consequences of policy interventions than complex models. The title itself is a good indication of the usefulness of simple mathematical models in addressing complex policy issues. Finally, if it were to be necessary to evaluate the impact of some kind of prevention campaign or intervention relating to problematic drug use among susceptibles, using the moverstayer model proposed above, that issue could be dealt with easily by introducing the possibility that a mover become a stayer during the period of the campaign (or soon after). This can be simply modelled by assuming that during a short period a given proportion of movers become stayers and modifying the equations of the model to incorporate that fact. Analysis of the new equations would then make it possible to state that the efficacy of a primary prevention intervention of that length is higher at the beginning of the epidemic and decreases in the following period [14]. The effect of a secondary prevention intervention can also be taken into account in a similar manner by modifying the equations of the model. Analysis of the new equations would make it possible to state that the efficacy of a secondary prevention intervention of the given length is lower at the beginning of the epidemic and increases in the following period. The model may thus always be used to forecast the effect of such types of intervention and to determine the best combination of the two [14]. Similarly, the effect of law enforcement, which may influence movers, pushing a proportion of them to become stayers as a result of information about the adverse consequences of using illegal drugs [1], can be modelled assuming some parameters depending either on the number of addicts assisted by the healthcare services or the number of incarcerated addicts or both.
Further generalizations might concern a more realistic approach to modelling the length of stay in the various compartments, taking into account the heterogeneity of individual behaviours. Most of these issues will be addressed in future contributions.
*The present article was funded in part by the European Monitoring Centre on Drugs and Drug Addiction in the framework of the European Network to Develop Policy Relevant Models and SocioEconomic Analyses of Drug Use, Consequences and Interventions. The simulation procedure used to obtain scenario analyses was developed by Alessandro Ghizzoni.
**From surveys conducted among military conscripts, reported in the annual report on the state of the drug problem in Italy for the year 1999, published by the National Focal Point, it appears that, in terms of the reasons for drug use, the two most mentioned factors were curiosity (more than 40 per cent) and peer group pressure (more than 30 per cent).
