Chaos in a predator-prey-based mathematical model for illicit drug consumption

TL;DR

This paper introduces a predator-prey mathematical model for illicit drug consumption, revealing complex dynamics including periodic and chaotic behaviors, and demonstrates its potential to inform policy through data calibration and bifurcation analysis.

Contribution

It proposes a novel predator-prey inspired model for drug use dynamics, incorporating stabilizing and destabilizing effects, and applies bifurcation analysis to real-world data from Colorado and Washington.

Findings

Model exhibits limit cycles matching observed periodic data.

Increasing non-user growth rate leads to chaotic dynamics.

Model aligns well with empirical data from Colorado and Washington.

Abstract

Recently, a mathematical model describing the illicit drug consumption in a population consisting of drug users and non-users has been proposed. The model describes the dynamics of non-users, experimental users, recreational users, and addict users within a population. The aim of this work is to propose a modified version of this model by analogy with the classical predator-prey models, in particular considering non-users as prey and users as predator. Hence, our model includes a stabilizing effect of the growth rate of the prey, and a destabilizing effect of the predator saturation. Functional responses of Verhulst and of Holling type II have been used for modeling these effects. To forecast the marijuana consumption in the states of Colorado and Washington, we used data from Hanley (2013) and a genetic algorithm to calibrate the parameters in our model. Assuming that the population of non-users increases in proportion with the demography, and following the seminal works of Sir Robert May (1976), we use the growth rate of non-users as the main bifurcation parameter. For the state of Colorado, the model first exhibits a limit cycle, which agrees quite accurately with the reported periodic data in Hanley (2013). By further increasing the growth rate of non-users, the population then enters into two chaotic regions, within which the evolution of the variables becomes unpredictable. For the state of Washington, the model also exhibits a periodic solution, which is again in good agreement with observed data. A chaotic region for Washington is likewise observed in the bifurcation diagram. Our research confirms that mathematical models can be a useful tool for better understanding illicit drug consumption, and for guiding policy-makers towards more effective policies to contain this epidemic.

Tables

Table 1. Table 1: Interpretation of the parameters in NERA model

Parameter	Sociological Meaning
$r_1$	Influence rate of $E\left(t\right)$ on $N\left(t\right)$
$r_2$	Influence rate of $R\left(t\right)$ on $E\left(t\right)$
$r_3$	Rate at which recreational users change to addicts
${\alpha }_1$	Influence rate of $A\left(t\right)$ on $N\left(t\right)$
${\alpha }_2$	Influence rate of $R\left(t\right)$ on $N\left(t\right)$
${\gamma }_1$	Influence rate of $A\left(t\right)$ on $E\left(t\right)$
${\beta }_1$	Rate of moving in and out of the Nonuser category
${\beta }_2$	Rate at which experimental users quit drugs
${\beta }_3$	Rate at which recreational users quit drugs
${\beta }_4$	Rate at which addicts quit drugs

Table 2. Table 2: Parameter values for the consumption of Marijuana in Colorado

$r_1$	$r_2$	$r_3$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\beta }_1$	${\beta }_2$	${\beta }_3$	${\beta }_4$
0.44	0.193	0.029	0.103	0.043	0.031	0.042	0.016	0.052	0.047

Table 3. Table 3: Parameter values for the consumption of Marijuana in Washington

$r_1$	$r_2$	$r_3$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\beta }_1$	${\beta }_2$	${\beta }_3$	${\beta }_4$
0.38	0.142	0.034	0.099	0.112	0.032	0.015	0.03	0.066	0.039

Table 4. Table 4: LCEs of NERA model ( 1 ) for Colorado for various values of β 1 subscript 𝛽 1 \beta_{1} .

${\beta }_1$	LCE spectrum	Dynamics of the attractor	Hausdorff dimension
	($0,-,-,-$)	Periodic Motion (Limit Cycle)	$D=1$
	($0,0,-,-$)	3-Torus (Quasi-Periodic Motion)	$D=2$
	($0,-,-,-$)	Periodic Motion (Limit Cycle)	$D=1$
	($0,0,-,-$)	3-Torus (Quasi-Periodic Motion)	$D=2$

Table 5. Table 5: LCEs of NERA model ( 1 ) for Washington for various values of β 1 subscript 𝛽 1 \beta_{1} .

${\beta }_1$	LCE spectrum	Dynamics of the attractor	Hausdorff dimension
	($0,-,-,-$)	Periodic Motion (Limit Cycle)	$D=1$
	($0,0,-,-$)	3-Torus	$D=2$
	($0,-,-,-$)	Periodic Motion (Limit Cycle)	$D=1$

Equations

$$\frac{dX}{dt}=\beta X\left(1-X\right)$$

$$\frac{X}{h+X}Y$$

$$\frac{dN}{dt}$$

$$\frac{dE}{dt}$$

$$\frac{dR}{dt}$$

$$\frac{dA}{dt}$$

$$O(0,0,0,0)\text{mbox}and{I}_1(1,0,0,0).$$

$$I_2(N_2^{\ast },E_2^{\ast },0,0)=\left(\frac{{\beta }_{2}h}{r_1-{\beta }_2},\frac{r_1-{\beta }_2\left(1+h\right)}{{\left(r_1-{\beta }_2\right)}^2}{\beta }_{1}h,0,0\right)$$

$$r_1-{\beta }_2>0\text{mbox}and{r}_1-{\beta }_2\left(1+h\right)>0.$$

$$J_2(N_2^{\ast },0,0,A_2^{\ast })=\left(\frac{{\beta }_{4}h}{{\alpha }_1-{\beta }_4},0,0,\frac{{\alpha }_1-{\beta }_4\left(1+h\right)}{{\left({\alpha }_1+{\beta }_4\right)}^2}{\beta }_{1}h\right)$$

$${\alpha }_1-{\beta }_4>0\text{mbox}and{\alpha }_1-{\beta }_4\left(1+h\right)>0.$$

$$\frac{dN}{dt}$$

$$\frac{dE}{dt}$$

$$\frac{dR}{dt}$$

$$\frac{dA}{dt}$$

$$J=\left(\begin{array}{cccc}{J}_{11}& -r_1{p}_1(N)& -{\alpha }_2{p}_1(N)& -{\alpha }_1{p}_1(N)\text{vspace}6pt\text{hfill}\\ r_{1}E{p}_1^{\prime }(N)& J_{22}& -r_2{p}_2(E)& -{\gamma }_1{p}_2(E)\text{vspace}6pt\text{hfill}\\ {\alpha }_{2}R{p}_1^{\prime }(N)& r_{2}R{p}_2^{\prime }(E)& J_{33}& -r_3{p}_3(R)\text{vspace}6pt\text{hfill}\\ {\alpha }_{1}A{p}_1^{\prime }(N)& {\gamma }_{1}A{p}_2^{\prime }(E)& -r_{3}A{p}_3^{\prime }(R)& J_{44}\text{vspace}6pt\text{hfill}\end{array}\right)$$

$$J_{11}\left(N,E,R,A\right)$$

$$J_{22}\left(N,E,R,A\right)$$

$$J_{33}\left(N,E,R,A\right)$$

$$J_{44}\left(N,E,R,A\right)$$

$$g\left(0\right)={\beta }_1\text{mbox},g^{\prime }\left(0\right)=-{\beta }_1\text{mbox},p_i\left(0\right)=0\text{mbox},p_i^{\prime }\left(0\right)=\frac{1}{h},$$

$$g\left(1\right)=0\text{mbox},g^{\prime }\left(1\right)=-{\beta }_1\text{mbox},p_i\left(1\right)=\frac{1}{h+1}\text{mbox},p_i^{\prime }\left(1\right)=\frac{h}{{\left(h+1\right)}^2}.$$

$${\lambda }_1={\beta }_1\text{mbox};{\lambda }_2=-{\beta }_2\text{mbox};{\lambda }_3=-{\beta }_3\text{mbox};{\lambda }_4=-{\beta }_4.$$

$${\lambda }_1=-{\beta }_1\text{mbox};{\lambda }_2=-{\beta }_2+\frac{r_1}{h+1}\text{mbox};{\lambda }_3=-{\beta }_3+\frac{{\alpha }_2}{h+1}\text{mbox};{\lambda }_4=-{\beta }_4+\frac{{\alpha }_1}{h+1}.$$

$$0<r_1-{\beta }_2\left(1+h\right)-h{r}_1\text{mbox}and0<r_1-{\beta }_2,$$

$$0<{\alpha }_1-{\beta }_4\left(1+h\right)-h{\alpha }_1\text{mbox}and0<{\alpha }_1-{\beta }_4.$$

$$0<h{r}_1<r_1-{\beta }_2\left(1+h\right),$$

$$0<h{\alpha }_1<{\alpha }_1-{\beta }_4\left(1+h\right).$$

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_2}{dt}$$

$$\frac{d{x}_3}{dt}$$

$$\frac{d{x}_4}{dt}$$

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_2}{dt}$$

$$\frac{d{x}_3}{dt}$$

$$\frac{d{x}_4}{dt}$$

Full text

I Introduction

Mathematical modeling of illicit drug consumption is a very difficult and complex problem. To this aim predator-prey models have been used at the end of the nineties. Then, in 2013 a model called NERA has been built to describe the dynamics of the non (N), experimental (E), recreational (R) and addict (A) user categories, respectively, within a given population. However, the original NERA model didn’t involve limitation in drug consumption and was consequently unable to transcribe the periodic evolution of each category. So, we have modified this model by analogy with the classical predator-prey models and while considering non-users (N) as prey and users (E, R and A) as predator. Then, using data from the state of Colorado and Washington and a genetic algorithm, we calibrated our predator-prey NERA models by estimating their parameters sets. This allowed us to account for the periodic evolution of each category. Then, by considering that the population of Nonusers increases in proportion of the demography we highlighted chaotic regions within which the evolution of the variables becomes unpredictable. Thus, it appears that our validated NERA model can be a precious tool in forecasting of illicit drug consumption and can be of substantial interest to policy-makers in the problematic of illicit drug consumption.

Since the beginning of the nineties, many continuous time models have been proposed to describe the dynamics of drug consumption Baveja1 ; Baveja2 ; Caulkins1 ; Caulkins2 ; Caulkins3 ; Gragnani1 ; Gragnani3 . They mainly consisted in second order nonlinear models involving two variables which were used for predicting the long term dynamics of addicts and dealers of a large drug market, like that of an entire country. Few years later, Gragnani et al. Gragnani2 proposed an extension to such models by adding a third ordinary differential equation to the original two-dimensional dynamical systems. This third variable represented the enforcement exerted by the authorities. Let’s notice on the one hand that these models already used limitations in the growth and decay of each variable (at least for the two first) and on the other hand that, Gragnani et al. Gragnani2 proved the existence of slow-fast limit cycles according to the singular perturbation theory as solution of their three-dimensional dynamical system. Recently, Dauhoo et al. Dauhoo have considered that drug users are generally classified into three main categories, depending on their consumption frequency and the control they have over the drug, i.e., Experimental (E), Recreational (R) and Addicts (A) users. By adding a fourth variable, i.e., the Non (N) users to these first three ones, they proposed the NERA model. Although this four-dimensional dynamical system takes into account the mutual influence that drug users (E, R and A) can have on non-users (N) and on each other it does not contain any limitation in the growth and decay of each variable. Thus, no oscillatory or chaotic regime could be observed. That’s the reason why we have proposed to modify Dauhoo’s NERA model Dauhoo by introducing a limitation in each “functional response”. In their paper, Dauhoo et al. Dauhoo wrote the following sentence: “Anyone could be a ‘prey’ to illicit drugs”. This led us to make the analogy with predator-prey models.

II General Predator-Prey NERA Model

This deterministic model aims at transcribing into mathematical functions variations of the number of individuals of each group due to their interactions with others. We make the assumption that such interactions are mainly characterized by the influence that individuals of one group may exert on the others. Thus, we consider that people influenced by a group leaves the group to which they belong to join the group which has influenced them. As a consequence, some individuals disappear of one group and appear in another. So, by analogy with the predator-prey models used for a long time in Theoretical Ecology Maybook ; Scudo , these influences, which cause such appearances and disappearances, i.e., increases and decreases (variations) of the number of individuals of each group can be regarded as predation of on group on another. In the predator-prey model that we propose, Non-users (N) represent prey for all other groups (E, R and A). Then, Experimental-users (E) are predator of both Recreational-users (R) and Non-users (N) while Addict-users (A) are predator of all R, E and N-users. According to Kuznetsov Kuznetsov , for the model to be realistic it is necessary to introduce a “stabilizing effect of the competition among prey and a destabilizing effect of the predator saturation”. To this aim, we have used two different types of “functional responses” for the growth of prey (N) and for the growth of the predators (E, R, A). We have first considered that, in the absence of any predator, prey growth (N) must be limited. Such limitation or stabilizing effect is generally introduced by using the logistic law introduced by Verhulst Verhulst . Concerning the predators (E, R, A), the saturation of the predator rate (destabilizing effect) can be modeled with the classical Holling type II “functional response” Holling1 ; Holling2 . Of course, various other “functional responses” could have been chosen to this aim GinouxVolterra . Let’s notice that such saturation in the predator rate represents a limitation in the influence of each predator group on the others. At last, we consider that in the absence of its predators, the number of individuals of one group (E, R, and A except N) can decrease by a kind of “natural mortality” which can correspond to individuals leaving this group. In the most dramatic case of drug addict, this could be due to overdose.

II.1 Functional responses

In 1837, the Belgian biologist Pierre-François Verhulst proposed a model that took into account the limitation imposed by the increasing population size of a prey called X in absence of any predator. This model, called logistic law, can be written as follows:

[TABLE]

In 1926, the Italian mathematician Vito Volterra developed the very first predator-prey models. The formulation of the equation representing the predation is based on the méthode des rencontres (method of encounters) and on the hypothèse des équivalents (hypothesis of equivalents) vol26 ; vol27 ; vol28 ; vol31 . The former assumes that for predation to occur between a predatory species (X) and a prey species (Y), it is first necessary to have encounters between these two species and that the number of encounters between them is proportional to the product of the number of individuals composing them, that is $\alpha X(t)Y(t)$, the coefficient of proportionality $\alpha$ being equal to the probability of an encounter. The second hypothesis consists in assuming that “there is a constant ratio between the disappearances and appearances of individuals caused by the encounters”, that is, that predation of the preys is equivalent to increase of the predators. At the beginning, Volterra considers this increase as immediate. This means that predation is immediately transcribed in terms of growth of the predator species, whereas its effect naturally occurs with some delay. In his works Volterra vol31 also took this delay into account. This won’t be the case in this paper. In the following, we will consider that the effects of influence that individuals of one group (X) may exert on another (Y) are analogous to the effects of predation of (Y) on (X). The mathematical function corresponding to the modeling of a behavior such as influence or predation is called a “functional response”. The functional response proposed by Volterra to describe predation does not take into account any limitation, i.e. “satiety” of the predator and so, the predation rate is a “linear function” of the prey. Thus, from the beginning of the thirties various types of “function responses” were proposed while using nonlinear mathematical functions presenting an asymptotic behavior and so a limitation gau ; gau2 .

In the late 1950s, entomologist Crawford Stanley Holling Holling1 ; Holling2 developed two new functional responses for predation, also intended to describe a certain satiety of the predator (Y) with respect to its prey (X): Holling function of type II and Holling function of type III. In this paper, we will only use the Holling type II functional response which can be represented by:

[TABLE]

where $h$ represents half-saturation, that is, the value of the prey density $X=h$ for which the predation level reaches a value equal to half its maximum.

So, in this work we propose to used both logistic law and Holling type II functional response for modeling the influence that exert the predators A, R and E on each others and on the prey N (See Fig. 1).

II.2 Model equations

So, we have the following system of ordinary differential equations:

[TABLE]

where $\beta_{1}$ is the growth rate of the population of the prey ($N$) in the absence of any predator ($E,R,A$), $\beta_{i}$ with $i=2,3,4$ are the “natural mortality” of each predator ($E,R,A$) in the absence of all others and $r_{i}$ with $i=1,2,3$ is the predation rate of $A$ on $R$, $R$ on $E$ and $E$ on $N$ respectively. $\alpha_{1}$ and $\gamma_{1}$ represent the predation rate of $A$ on $N$ and $E$ respectively while $\alpha_{2}$ is that of $R$ on $N$. Thus, in this four-dimensional dynamical system, a set of 11 positive parameters: $(\beta_{1},\beta_{2},\beta_{3},\beta_{4},r_{1},r_{2},r_{3},\alpha_{1},\alpha_{2},\gamma_{1},h)$ is used.

Remark. Let’s notice that for each Holling type II functional response a different half saturation $h$ could have been chosen. Nevertheless, the aim of this work is to propose the most simple and consistent model for illicit drug consumption.

The sociological meaning of each parameter used in this model (1) is given in Table 1.

II.3 Dynamic aspects

Due to the presence of the Holling type II functional responses, the determination of the fixed points of this four-dimensional dynamical system (1) is not trivial while using the classical nullclines method. Nevertheless, by posing $E=R=A=0$ two obvious fixed points can be easily found:

[TABLE]

By posing $R=A=0$, a third fixed point can be also obtained:

[TABLE]

It follows that this fixed point $I_{2}$ is positive provided that:

[TABLE]

Then, by posing $E=R=0$, a fourth fixed point can be also obtained:

[TABLE]

It follows that this fixed point $J_{2}$ is positive provided that:

[TABLE]

Finally, while posing $A=0$, a fifth fixed point $I_{3}(N_{3}^{*},E_{3}^{*},R_{3}^{*},0)$ can be determined but its expression is too long to be written here.

Remark. Let’s notice that all fixed points depend on parameters $(\beta_{1})$.

Following the works of Freedman Freed , the system (1) may be written as:

[TABLE]

where $g(N)=\beta_{1}(1-N)$, $p_{1}(N)=\dfrac{N}{h+N}$, $p_{2}(E)=\dfrac{E}{h+E}$ and $p_{3}(R)=\dfrac{R}{h+R}$. Such a formulation will simplify the computation of the eigenvalues of the functional Jacobian matrix (5) presented below.

II.3.1 Functional Jacobian matrix

The Jacobian matrix of system (4) reads:

[TABLE]

where the diagonal terms read:

[TABLE]

Let’s notice that for $i=1,2,3$ we have for $N=0$ and then for $N=1$:

[TABLE]

II.3.2 Eigenvalues at O(0,0,0,0)

Taking into account the above conditions (6), the functional Jacobian matrix (5) is diagonal. Thus, the four eigenvalues are:

[TABLE]

It follows that the origin $O$ is a saddle.

II.3.3 Eigenvalues at $I_{1}$(1,0,0,0)

Taking into account the above conditions, the functional Jacobian matrix (5) is diagonal. Thus, the four eigenvalues are:

[TABLE]

According to conditions (2-3), both $\lambda_{2}$ and $\lambda_{3}$ are positive. So, whatever the values of the parameters, it follows that $I_{1}$ is also a saddle.

II.3.4 Eigenvalues at $I_{2}$($\rm{\textit{N}_{2}^{*}}$, $\rm{\textit{E}_{2}^{*}}$, 0, 0) and at $J_{2}$($\rm{\textit{N}_{2}^{*}}$, 0, 0, $\rm{\textit{A}_{2}^{*}}$)

Although the four eigenvalues evaluated at $I_{2}$ and $J_{2}$ are too long to be expressed, two of them contain a square root. So, according to the choice of the parameters, these eigenvalues may be complex conjugate. Thus, if we assume that the expression in the square root is negative, we can look for the sign of the remaining part which can be considered as the real part of these eigenvalues. Such a real part is positive provided that:

[TABLE]

[TABLE]

Combining these conditions with the previous one (2-3), we find that

[TABLE]

It follows that $I_{2}$ and $J_{2}$ are a saddle-foci (two eigenvalues are complex conjugate with positive real parts and the two others eigenvalues are real).

II.3.5 Eigenvalues at $I_{3}$($\rm{\textit{N}_{3}^{*}}$, $\rm{\textit{E}_{3}^{*}}$, $\rm{\textit{R}_{3}^{*}}$, 0)

Concerning this last point $I_{3}$, the analytical analysis of stability is no more possible and it becomes then necessary to choose a parameter set.

Remark. Let’s notice that the number of real positive fixed points depends on the choice of parameters.

II.3.6 Bifurcation analysis

Since $I_{2}$, $J_{2}$ and $I_{3}$ do not depend on parameters $(\beta_{4},r_{3},\alpha_{1},\gamma_{1})$, bifurcations can occur in these models (for both states of Colorado and Washington). Following the works of May May1976 , we propose to choose the parameter $\beta_{1}$, i.e., the growth rate of the population of the prey ($N$) or the rate of moving in and out of the Nonuser category, as bifurcation parameter. Let’s notice that this choice is based on the assumption according to which the population of Nonusers increases in proportion of the demography. Then, for the same reasons as previously, an analytical analysis would be difficult even impossible. So, in the next section we will set all the parameters except $\beta_{1}$ and then we will use a bifurcation diagram to determine the values of the bifurcation parameters.

II.3.7 Existence of bounded solutions

Following the works of Freedman Freed and while posing $N=x_{1}$, $E=x_{2}$, $R=x_{3}$ and $A=x_{4}$, the system (1) can be rewritten as follows:

[TABLE]

where $g(x_{1})=\beta_{1}(1-x_{1})$, $p_{i}(x_{i})=\dfrac{x_{i}}{h+x_{i}}$ with $i=1,2,3$.

Moreover, analysis of experimental data available on the prevalence of marijuana in the population of 21+ in the states of Colorado and Washington Hanley has shown that the influence of $A$ on $N$ and $E$ as well as that of $R$ on $N$ are in fact very weak. So, we will consider that $\alpha_{1}\ll 1$, $\alpha_{2}\ll 1$ and $\gamma_{1}\ll 1$. Thus, under these assumptions, model (8) reads:

[TABLE]

Then, according to Theorem 2.1 stated by Freedman (Freed, , p. 72), “all solutions initiating in the nonnegative cone are bounded and eventually enter a certain attracting set described below.”

III Applications of Predator-Prey NERA Model

In order to perform numerical experiments for forecasting the marijuana consumption in the states of Colorado and Washington beyond the year of the I – 502 implementation, we used data from Hanley Hanley . The Washington State Institute for Public Policy conducted a benefit-cost analysis of the implementation of I – 502, which legalizes recreational marijuana use for adults within the two states Hanley . The data collected in the latter report is the result of an analysis of population-level data in order to monitor four key indicators of marijuana use, namely, current marijuana use, lifetime marijuana use, marijuana abuse or dependency and age of initiation prior to implementation of I – 502. In the case of our NERA model, the first three categories correspond to the Recreational, Experimental and Addict category respectively Dauhoo . The report highlights the importance of examining trends in that manner will allow them to monitor whether the implementation of I – 502 appears to affect these key indicators of marijuana use over time. Our numerical experiments aim to forecast the marijuana consumption in the states of Colorado and Washington beyond the year of the I – 502 implementation using data from Hanley Hanley . The NERA predator-prey model is calibrated by estimating the parameters $r_{1}$, $r_{2}$, $r_{3}$, $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$, $\beta_{1}$, $\beta_{2}$, $\beta_{3}$ and $\beta_{4}$. Parameter $h$ has been arbitrarily chosen equal to $1/2$ as usually done in theoretical ecology Scudo . To this aim, we used these data and a genetic algorithm as explained in Dauhoo et al. Dauhoo . The fitness function used has been chosen to minimize the error that our model generates. MATLAB Optimtool is used and the fitness function is inserted in the genetic algorithm. We hence obtain the set of values in Tab. 2 and Tab. 3 corresponding to the NERA predator-prey model for Colorado and Washington. Obviously, since the parameters sets of both NERA predator-prey models have been calibrated starting from the data from Hanley Hanley , numerical integration of Eqs. (1) with parameters sets from Tab. 2 and 3 are in perfect agreement with the observed data. Thus, it confirms that the NERA predator-prey models can be a precious tool in forecasting of illicit drug consumption in the population of the state considered.

III.1 NERA model for Colorado

Still using experimental data Hanley , we set the parameters from Tab. 2 where $\beta_{1}$ is the bifurcation parameter. By varying continuously $\beta_{1}$ from 0.02 to 0.8 (all other parameters are those given in Tab. 2), we determine the values for which bifurcations occur by plotting the bifurcation diagram presented in Fig. 2.

We observe in this diagram (Fig. 2), that for the value of $\beta_{1}\in[0.02,0.31]$ the solution of the NERA model (1) for Colorado is a limit cycle (LC) as exemplified in Fig. 3a. Then, still increasing parameter $\beta_{1}$ up to the first bifurcation value $\beta_{1}^{b_{1}}=0.3175$, the solution becomes chaotic (C) as highlighted in Fig. 3b. Such chaotic feature persists up to the second bifurcation value $\beta_{1}^{b_{2}}=0.37$ starting from which the solution becomes again a limit cycle (LC) (see Fig 3c). Then, starting from the third bifurcation value $\beta_{1}^{b_{3}}=0.55$ a chaotic attractor (C) appears again (see Fig. 3d).

The algorithm developed by Marco Sandri Sandri for Mathematica${}^{\mbox{\scriptsize{\textregistered}}}$ has been used to perform the numerical calculation of the Lyapunov characteristics exponents (LCE) of the NERA predator-prey model (1) for Colorado with the parameters from Tab. 2 and with $\beta_{1}\in[0.02,0.8]$. As an example for $\beta_{1}=0.3,0.35,0.4$ and $0.6$, Sandri’s algorithm has provided respectively the following LCEs $(0,-0.0022,-0.013,-0.029)$, $(0,0,-0.0010,-0.065)$, $(0,-0.0011,-0.0011,-0.066)$ and $(0,0,-0.0017,-0.085)$. Then, according to the works of Klein and Baier KleinBaier , a classification of (autonomous) continuous-time attractors of dynamical system (1) on the basis of their Lyapunov spectrum, together with their Hausdorff dimension is presented in Tab. 4. LCEs values have been also computed with the Lyapunov Exponents Toolbox (LET) developed by Pr. Steve Siu for MatLab${}^{\mbox{\scriptsize{\textregistered}}}$ and involving the two algorithms proposed by Wolf et al. Wolf and Eckmann and Ruelle Eckmann (see https://fr.mathworks.com/matlabcentral/fileexchange/233-let). Results obtained by both algorithms are consistent.

III.2 NERA model for Washington

Now, we set the parameters from Tab. 3 where $\beta_{1}$ is still the bifurcation parameter. By varying continuously $\beta_{1}$ from 0.02 to 0.36 (all other parameters are those given in Tab. 3), we determine the values for which bifurcations occur by plotting the bifurcation diagram presented in Fig. 4.

We observe from this diagram (Fig. 4) that for the value of $\beta_{1}\in[0.02,0.334]$ the solution of the NERA model (1) for Washington is a limit cycle (LC) in the ($N,E$)-plane as exemplified in Fig. 5a. Then, still increasing parameter $\beta_{1}$ up to the first bifurcation value $\beta_{1}^{b_{1}}=0.335$, the solution becomes chaotic (C) as highlighted in Fig. 5b & Fig. 5c. Such chaotic feature persists up to the second bifurcation value $\beta_{1}^{b_{2}}=0.357$ starting from which the solution becomes again a limit cycle (see Fig 5d).

Still using the algorithm developed by Marco Sandri Sandri we numerically compute the Lyapunov characteristics exponents (LCE) of the NERA model (1) for Washington with the parameters from Tab. 3 and with $\beta_{1}\in[0.02,0.36]$. In this case, for $\beta_{1}=0.27,0.34,0.345$ and $0.3485$, Sandri’s algorithm provides respectively the following LCEs $(0,-0.0078,-0.0089,-0.019)$, $(0,0,-0.010,-0.017)$, $(0,0,-0.0092,-0.016)$ and $(0,0,-0.011,-0.016)$.

Then, as previously, a classification of attractors of dynamical system (1) on the basis of its Lyapunov spectrum, together with its Hausdorff dimension is presented in Tab. 5. Here again, LCEs values have been also computed with the Lyapunov Exponents Toolbox (LET) and results obtained by both algorithms are consistent.

IV Discussion

By considering that drug users are classified into four main categories: non (N), experimental (E), recreational (R) and addicts (A) users, Dauhoo et al. Dauhoo have proposed the NERA model. Nevertheless, although this four-dimensional dynamical system took into account the mutual influence that drug users (E, R and A) can have on non-users (N) and on each other, it did not contain any limitation in the growth and decay of each variable. Thus, no oscillatory or chaotic regime could be observed. So, the aim of this work was to propose a modified version of this NERA model by analogy with the classical predator-prey models and while considering non-users (N) as prey and users (E, R and A) as predator. Thus, this new model included a “stabilizing effect” of the growth rate of the preys (N) and a “destabilizing effect” of the predators (E, R and A) saturation. Functional responses of Verhulst and Holling type II have been used for modeling these effects. Then, in order to perform numerical experiments for forecasting the marijuana consumption in the states of Colorado and Washington beyond the year of the I – 502 implementation, we used data from Hanley Hanley and a genetic algorithm as explained in Dauhoo et al. Dauhoo . Thus, the NERA predator-prey model has been calibrated by estimating all the parameters except $h$ which has been arbitrarily chosen equal to $1/2$ as usually done in theoretical ecology Scudo . We hence obtained two parameters sets corresponding to the NERA predator-prey model for the state of Colorado and Washington. Following the works of May May1976 , we chose the parameter $\beta_{1}$, i.e., the growth rate of the population of the prey ($N$) or the rate of moving in and out of the Nonuser category, as bifurcation parameter. This choice was based on the assumption that the population of Nonusers increases in proportion of the demography. A stability and bifurcation analysis of the NERA model for Colorado and for Washington has therefore been performed.

Concerning the NERA model for Colorado, the bifurcation diagram has shown that when the value of parameter $\beta_{1}\in[0.02,0.31]$, the solution is a limit cycle confirming thus the behavior of the observed data from Hanley Hanley . So, the number of individuals of each group N, E, R and A oscillates in a deterministic way with a period and amplitude that can be numerically computed. Then, still increasing parameter $\beta_{1}$ up to the first bifurcation value $\beta_{1}^{b_{1}}=0.3175$, we have shown that the solution becomes quasi periodic (3-Torus). Such chaotic attractor persists up to the second bifurcation value $\beta_{1}^{b_{2}}=0.37$ starting from which the solution becomes again a limit cycle. When parameter $\beta_{1}^{b_{3}}=0.55$ reaches the third bifurcation value a chaotic attractor appears again. These results have been confirmed by the computation of the Lyapunov characteristics exponents.

Concerning the NERA model for Washington, the bifurcation diagram has shown that when the value of parameter $\beta_{1}\in[0.02,0.26]$ the solution is a limit cycle confirming thus the behavior of the observed data from Hanley Hanley . Then, still increasing parameter $\beta_{1}$ up to the first bifurcation value $\beta_{1}^{b_{1}}=0.34$, the solution becomes quasi periodic (3-Torus). Such chaotic attractor persists up to the second bifurcation value $\beta_{1}^{b_{2}}=0.357$ starting from which the solution becomes again a limit cycle. These results have been also confirmed by the computation of the Lyapunov characteristics exponents.

Thus, we have shown that an increase of the population of non-users (N) leads first to a periodic evolution of drug-users (E, R, A). But when this increase reaches a certain threshold for both states of Colorado and Washington, the evolution of all variables becomes unpredictable. Such result can be also interpreted by analogy with the so-called “paradox of enrichment” which states that increasing the food available to the prey caused the predator’s population to destabilize Rosenzweig . We have thus confirmed on one hand that the NERA predator-prey models can be a precious tool in forecasting of illicit drug consumption in the population of the state considered and, on the other hand that they can be of substantial interest to policy-makers in the problematic of illicit drug consumption.

Acknowledgements.

Matjaž Perc was supported by the Slovenian Research Agency (Grants J1-7009, J4-9302, J1-9112 and P5-0027)