Access the full text.
Sign up today, get DeepDyve free for 14 days.
C. Holling (1959)
The Components of Predation as Revealed by a Study of Small-Mammal Predation of the European Pine SawflyThe Canadian Entomologist, 91
(1996)
Numerical calculation of Lyapunov exponents
J. Ginoux (2017)
The paradox of Vito Volterra’s predator-prey modelLettera Matematica, 5
J. Caulkins (1993)
Preface: Mathematical models of drug markets and drug policyMathematical and Computer Modelling, 17
J. Caulkins (1993)
Local Drug Markets' Response to Focused Police EnforcementOper. Res., 41
Noure-Roukayya Adam, M. Dauhoo, O. Kavian (2015)
An Analysis of the Dynamical Evolution of Experimental, Recreative and Abusive Marijuana Consumption in the States of Colorado and Washington Beyond the Implementation of I–502The Journal of Mathematical Sociology, 39
(1993)
Collapsing street markets for illicit drugs : the benefits of being decisive
(1991)
Hierachies of dynamical systems, in A Chaotic Hierarchy, edited by G
V. Volterra (1928)
Variations and Fluctuations of the Number of Individuals in Animal Species living togetherIces Journal of Marine Science, 3
A. Gragnani, S. Rinaldi, G. Feichtinger (1997)
Dynamics of Drug Consumption: a Theoretical ModelSocio-economic Planning Sciences, 31
C. Holling (1959)
Some Characteristics of Simple Types of Predation and ParasitismThe Canadian Entomologist, 91
A. Wolf, J. Swift, H. Swinney, J. Vastano (1985)
Determining Lyapunov exponents from a time seriesPhysica D: Nonlinear Phenomena, 16
M. Rosenzweig (1971)
Paradox of Enrichment: Destabilization of Exploitation Ecosystems in Ecological TimeScience, 171
(1965)
A case of the existence of a denumerable set of periodic motions
J. Marsden, M. McCracken, Springer New, York Berlin (1976)
The Hopf Bifurcation and Its Applications
G. Gauze, G. Teissier (1935)
Vérifications expérimentales de la théorie mathématique de la lutte pour la vie
(2013)
Legalization of recreational marijuana in Washington: Monitoring trends in use prior to the implementation of I-502
M. Dauhoo, B. Korimboccus, S. Issack (2013)
On the dynamics of illicit drug consumption in a given populationIma Journal of Applied Mathematics, 78
F. Scudo, J. Ziegler (1978)
The Golden Age of Theoretical Ecology: 1923–1940
V. Volterra, M. Brelot (1931)
Leçons sur la théorie mathématique de la lutte pour la vie
R. May, A. McLean (1977)
Theoretical Ecology: Principles and Applications
(1942)
Abzweigung einer periodischen Lösung von einer stationären Lösung eines Differentialsystems
H. Poincaré
méthodes nouvelles de la mécanique céleste
V. Volterra
Fluctuations in the Abundance of a Species considered MathematicallyNature, 118
F. Nakajima (2008)
The paradox of enrichment
J. Eckmann, D. Ruelle (1985)
Ergodic theory of chaos and strange attractorsReviews of Modern Physics, 57
A. Gragnani, S. Rinaldi, G. Feichtinger (1995)
Slow-Fast limit cycles in controlled drug markets
Alok Baveja, R. Batta, J. Caulkins, M. Karwan (1993)
Modeling the response of illicit drug markets to local enforcementSocio-economic Planning Sciences, 27
R. May (1976)
Simple mathematical models with very complicated dynamicsNature, 261
H. Freedman, J. So (1985)
Global stability and persistence of simple food chainsBellman Prize in Mathematical Biosciences, 76
(1928)
Variation and ﬂuctuations of the number of individuals in animal species living together . Translated by Miss Mary Evelyn Wells (“Journal du Conseil international l’exploration de la mer”
J. Caulkins, R. Larson, Thomas Rich (1993)
Geography's impact on the success of focused local drug enforcement operationsSocio-economic Planning Sciences, 27
鈴木 良明 (2005)
The Golden Age of Theoretical Ecology:1923-1940 : P18～33 (Biomathematics Kyoto Summer School)
(1976)
An English translation, with comments, is included as Section 5 in Marsden et al
(1967)
The existence of a denumerable set of periodic motions in four-dimensional space in an extended neighborhood of a saddle-focus
V. Volterra
[Book Reviews]Nature, 119
(1838)
Notice sur la loi que suit la population dans son accroissement
G. Gause (1971)
The struggle for existence
1,2 3 4 4 5,6,7 Jean-Marc Ginoux , Roomila Naeck , Yusra Bibi Ruhomally , Muhammad Zaid Dauhoo , Matjaˇz Perc Laboratoire LIS, CNRS, UMR 7020, Universit´e de Toulon, BP 20132, F-83957 La Garde cedex, France Departament de Matema`tiques, Universitat Auto`noma de Barcelona, 08193 Bellaterra, Barcelona, Catalonia, Spain PSASS, Ecoparc de Sologne, Domaine de Villemorant, Neug-sur Beuvron 41210, France Department of Mathematics, University of Mauritius, Reduit, Mauritius Faculty of Natural Sciences and Mathematics, University of Maribor, Koroˇska cesta 160, SI-2000 Maribor, Slovenia Center for Applied Mathematics and Theoretical Physics, University of Maribor, Mladinska 3, SI-2000 Maribor, Slovenia and Complexity Science Hub Vienna, Josefst¨adterstraße 39, A-1080 Vienna, Austria 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 modiﬁed 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 eﬀect of the growth rate of the prey, and a destabilizing eﬀect of the predator saturation. Functional responses of Verhulst and of Holling type II have been used for modeling these eﬀects. 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 ﬁrst 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 conﬁrms that mathematical models can be a useful tool for better understanding illicit drug consumption, and for guiding policy-makers towards more eﬀective policies to contain this epidemic. I. INTRODUCTION Mathematical modeling of illicit drug consumption is a very diﬃcult 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 modiﬁed 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 [1–5, 12, 14]. 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. [13] proposed an extension to such models by adding a third ordinary diﬀerential 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 ﬁrst) and on the other hand that, Gragnani et al. [13] 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. [6] have considered that drug users are generally classiﬁed 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 ﬁrst three ones, they proposed the NERA model. Although this four-dimensional dynamical system takes into account the mutual inﬂuence that drug users (E, R and A) can have on non-users (N) and on each other it does arXiv:1907.12973v1 [math.DS] 29 Jul 2019 2 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 [6] by introducing a limitation in each “functional response”. In their paper, Dauhoo et al. [6] 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 inﬂuence that individuals of one group may exert on the others. Thus, we consider that people inﬂuenced by a group leaves the group to which they belong to join the group which has inﬂuenced 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 [22, 27], these inﬂuences, 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 [20], for the model to be realistic it is necessary to introduce a “stabilizing eﬀect of the competition among prey and a destabilizing eﬀect of the predator saturation”. To this aim, we have used two diﬀerent types of “functional responses” for the growth of prey (N) and for the growth of the predators (E, R, A). We have ﬁrst considered that, in the absence of any predator, prey growth (N) must be limited. Such limitation or stabilizing eﬀect is generally introduced by using the logistic law introduced by Verhulst [30]. Concerning the predators (E, R, A), the saturation of the predator rate (destabilizing eﬀect) can be modeled with the classical Holling type II “functional response” [16, 17]. Of course, various other “functional responses” could have been chosen to this aim [11]. Let’s notice that such saturation in the predator rate represents a limitation in the inﬂuence 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. A. Functional responses In 1837, the Belgian biologist Pierre-Franc¸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: dX = βX (1 − X) dt In 1926, the Italian mathematician Vito Volterra developed the very ﬁrst predator-prey models. The formulation of the equation representing the predation is based on the m´ ethode des rencontres (method of encounters) and on the hypoth`ese des ´equivalents (hypothesis of equivalents) [31–34]. The former assumes that for predation to occur between a predatory species (X) and a prey species (Y), it is ﬁrst 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 αX(t)Y (t), the coeﬃcient of proportionality α 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 eﬀect naturally occurs with some delay. In his works Volterra [34] also took this delay into account. This won’t be the case in this paper. In the following, we will consider that the eﬀects of inﬂuence that individuals of one group (X) may exert on another (Y) are analogous to the eﬀects of predation of (Y) on (X). The mathematical function corresponding to the modeling of a behavior such as inﬂuence 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 [9, 10]. 3 In the late 1950s, entomologist Crawford Stanley Holling [16, 17] 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: h + X 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 inﬂuence that exert the predators A, R and E on each others and on the prey N (See Fig. 1). B. Model equations So, we have the following system of ordinary diﬀerential equations: dN N N N = β N (1 − N) − r E − α A − α R, 1 1 1 2 dt h + N h + N h + N dE N E E = r E − r R − β E − γ A, 1 2 2 1 dt h + N h + E h + E (1) dR E R N = r R − β R − r A + α R, 2 3 3 2 dt h + E h + R h + N dA R N E = r A − β A + α A + γ A, 3 4 1 1 dt h + R h + N h + E where β is the growth rate of the population of the prey (N) in the absence of any predator (E, R, A), β with 1 i i = 2, 3, 4 are the “natural mortality” of each predator (E, R, A) in the absence of all others and r with i = 1, 2, 3 is the predation rate of A on R, R on E and E on N respectively. α and γ represent the predation rate of A on N 1 1 and E respectively while α is that of R on N. Thus, in this four-dimensional dynamical system, a set of 11 positive parameters: (β , β , β , β , r , r , r , α , α , γ , h) is used. 1 2 3 4 1 2 3 1 2 1 Remark. Let’s notice that for each Holling type II functional response a diﬀerent 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 I. TABLE I: Interpretation of the parameters in NERA model Parameter Sociological Meaning r Inﬂuence rate of E(t) on N(t) r Inﬂuence rate of R(t) on E(t) r Rate at which recreational users change to addicts α Inﬂuence rate of A(t) on N(t) α Inﬂuence rate of R(t) on N(t) γ Inﬂuence rate of A(t) on E(t) β Rate of moving in and out of the Nonuser category β Rate at which experimental users quit drugs β Rate at which recreational users quit drugs β Rate at which addicts quit drugs 4 4 FIG. 1: Schematic representation of the predator-prey NERA model. 5 C. Dynamic aspects Due to the presence of the Holling type II functional responses, the determination of the ﬁxed 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 ﬁxed points can be easily found: O(0, 0, 0, 0) and I (1, 0, 0, 0). By posing R = A = 0, a third ﬁxed point can be also obtained: β h r − β (1 + h) 2 1 2 ∗ ∗ I (N , E , 0, 0) = , β h, 0, 0 2 1 2 2 2 r − β 1 2 (r − β ) 1 2 It follows that this ﬁxed point I is positive provided that: r − β > 0 and r − β (1 + h) > 0. (2) 1 2 1 2 Then, by posing E = R = 0, a fourth ﬁxed point can be also obtained: β h α − β (1 + h) 4 1 4 ∗ ∗ J (N , 0, 0, A ) = , 0, 0, β h 2 1 2 2 α − β 1 4 (α + β ) 1 4 It follows that this ﬁxed point J is positive provided that: α − β > 0 and α − β (1 + h) > 0. (3) 1 4 1 4 ∗ ∗ ∗ Finally, while posing A = 0, a ﬁfth ﬁxed point I (N , E , R , 0) can be determined but its expression is too long 3 3 3 to be written here. Remark. Let’s notice that all ﬁxed points depend on parameters (β ). Following the works of Freedman [8], the system (1) may be written as: dN = Ng (N) − (r E + α A + α R) p (N) , 1 1 2 1 dt dE = E (−β + r p (N)) − (r R + γ A) p (E) , 2 1 1 2 1 2 dt (4) dR = R (−β + r p (E) + α p (N)) − r Ap (E) , 3 2 2 2 1 3 3 dt dA = A [−β + r p (R) + α p (N) + γ p (E)] , 4 3 3 1 1 1 2 dt N E R where g(N) = β (1 − N), p (N) = , p (E) = and p (R) = . Such a formulation will simplify 1 1 2 3 h + N h + E h + R the computation of the eigenvalues of the functional Jacobian matrix (5) presented below. 1. Functional Jacobian matrix The Jacobian matrix of system (4) reads: J −r p (N) −α p (N) −α p (N) 11 1 1 2 1 1 1 r Ep (N) J −r p (E) −γ p (E) 1 22 2 2 1 2 J = (5) ′ ′ α Rp (N) r Rp (E) J −r p (R) 2 2 33 3 3 1 2 ′ ′ ′ α Ap (N) γ Ap (E) −r Ap (R) J 1 1 3 44 1 2 3 6 where the diagonal terms read: ′ ′ J (N, E, R, A) = g (N) + Ng (N) − (r E + α A + α R) p (N) , 11 1 1 2 J (N, E, R, A) = −β + r p (N) − γ up (E) , 22 2 1 1 1 J (N, E, R, A) = −β + r p (E) + α p (N) − r Ap (R) , 33 3 2 2 2 1 3 J (N, E, R, A) = −β + r p (R) + α p (N) + γ p (E) . 44 4 3 3 1 1 1 2 Let’s notice that for i = 1, 2, 3 we have for N = 0 and then for N = 1: ′ ′ g (0) = β , g (0) = −β , p (0) = 0, p (0) = , 1 1 i (6) 1 h ′ ′ g (1) = 0, g (1) = −β , p (1) = , p (1) = . 1 i h + 1 (h + 1) 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: λ = β ; λ = −β ; λ = −β ; λ = −β . 1 1 2 2 3 3 4 4 It follows that the origin O is a saddle. 3. Eigenvalues at I1(1,0,0,0) Taking into account the above conditions, the functional Jacobian matrix (5) is diagonal. Thus, the four eigenvalues are: r α α 1 2 1 λ = −β ; λ = −β + ; λ = −β + ; λ = −β + . 1 1 2 2 3 3 4 4 h + 1 h + 1 h + 1 According to conditions (2-3), both λ and λ are positive. So, whatever the values of the parameters, it follows 2 3 that I is also a saddle. ∗ ∗ ∗ ∗ 4. Eigenvalues at I (N , E , 0, 0) and at J (N , 0, 0, A ) 2 2 2 2 2 2 Although the four eigenvalues evaluated at I and J are too long to be expressed, two of them contain a square root. 2 2 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: 0 < r − β (1 + h) − hr and 0 < r − β , 1 2 1 1 2 0 < α − β (1 + h) − hα and 0 < α − β . 1 4 1 1 4 Combining these conditions with the previous one (2-3), we ﬁnd that 0 < hr < r − β (1 + h) , 1 1 2 (7) 0 < hα < α − β (1 + h) . 1 1 4 It follows that I and J are a saddle-foci (two eigenvalues are complex conjugate with positive real parts and the 2 2 two others eigenvalues are real). 7 ∗ ∗ ∗ 5. Eigenvalues at I (N , E , R , 0) 3 3 3 3 Concerning this last point I , 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 ﬁxed points depends on the choice of parameters. 6. Bifurcation analysis Since I , J and I do not depend on parameters (β , r , α , γ ), bifurcations can occur in these models (for both 2 2 3 4 3 1 1 states of Colorado and Washington). Following the works of May [23], we propose to choose the parameter β , 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 diﬃcult even impossible. So, in the next section we will set all the parameters except β and then we will use a bifurcation diagram to determine the values of the bifurcation parameters. 7. Existence of bounded solutions Following the works of Freedman [8] and while posing N = x , E = x , R = x and A = x , the system (1) can be 1 2 3 4 rewritten as follows: dx = x g (x ) − (r x + α x + α x ) p (x ) , 1 1 1 2 1 4 2 3 1 1 dt dx = x (−β + r p (x )) − (r x + γ x ) p (x ) , 2 2 1 1 1 2 3 1 4 2 2 dt (8) dx = x (−β + r p (x ) + α p (x )) − r x p (x ) , 3 3 2 2 2 2 1 1 3 4 3 3 dt dx = x [−β + r p (x ) + α p (x ) + γ p (x )] , 4 4 3 3 3 1 1 1 1 2 2 dt where g(x ) = β (1 − x ), p (x ) = with i = 1, 2, 3. 1 1 1 i i h + x Moreover, analysis of experimental data available on the prevalence of marijuana in the population of 21+ in the states of Colorado and Washington [15] has shown that the inﬂuence of A on N and E as well as that of R on N are in fact very weak. So, we will consider that α ≪ 1, α ≪ 1 and γ ≪ 1. Thus, under these assumptions, model (8) 1 2 1 reads: dx = x g (x ) − r x p (x ) , 1 1 1 2 1 1 dt dx = x [−β + r p (x )] − r x p (x ) , 2 2 1 1 1 2 3 2 2 dt (9) dx = x [−β + r p (x )] − r x p (x ) , 3 3 2 2 2 3 4 3 3 dt dx = x [−β + r p (x )] . 4 4 3 3 3 dt Then, according to Theorem 2.1 stated by Freedman [8, 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 [15]. The Washington 8 State Institute for Public Policy conducted a beneﬁt-cost analysis of the implementation of I – 502, which legalizes recreational marijuana use for adults within the two states [15]. 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 ﬁrst three categories correspond to the Recreational, Experimental and Addict category respectively [6]. The report highlights the importance of examining trends in that manner will allow them to monitor whether the implementation of I – 502 appears to aﬀect 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 [15]. The NERA predator-prey model is calibrated by estimating the parameters r , r , r , α , α , α , β , β , β and β . Parameter h has been 1 2 3 1 2 3 1 2 3 4 arbitrarily chosen equal to 1/2 as usually done in theoretical ecology [27]. To this aim, we used these data and a genetic algorithm as explained in Dauhoo et al. [6]. The ﬁtness function used has been chosen to minimize the error that our model generates. MATLAB Optimtool is used and the ﬁtness 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 [15], numerical integration of Eqs. (1) with parameters sets from Tab. 2 and 3 are in perfect agreement with the observed data. Thus, it conﬁrms that the NERA predator-prey models can be a precious tool in forecasting of illicit drug consumption in the population of the state considered. TABLE II: Parameter values for the consumption of Marijuana in Colorado r r r α α α β β β β 1 2 3 1 2 3 1 2 3 4 0.44 0.193 0.029 0.103 0.043 0.031 0.042 0.016 0.052 0.047 TABLE III: Parameter values for the consumption of Marijuana in Washington r r r α α α β β β β 1 2 3 1 2 3 1 2 3 4 0.38 0.142 0.034 0.099 0.112 0.032 0.015 0.03 0.066 0.039 FIG. 2: Bifurcation diagram of model (1) u as function of β . max 1 A. NERA model for Colorado Still using experimental data [15], we set the parameters from Tab. II where β is the bifurcation parameter. By varying continuously β from 0.02 to 0.8 (all other parameters are those given in Tab. 2), we determine the values for 1 9 which bifurcations occur by plotting the bifurcation diagram presented in Fig. 2. We observe in this diagram (Fig. 2), that for the value of β ∈ [0.02, 0.31] the solution of the NERA model (1) for Colorado is a limit cycle (LC) as exempliﬁed in Fig. 3a. Then, still increasing parameter β up to the ﬁrst bifurcation value β = 0.3175, the solution becomes chaotic (C) as highlighted in Fig. 3b. Such chaotic feature persists up to the second bifurcation value β = 0.37 starting from which the solution becomes again a limit cycle (LC) (see Fig 3c). Then, starting from the third bifurcation value β = 0.55 a chaotic attractor (C) appears again (see Fig. 3d). (a) β = 0.30 (b) β = 0.35 1 1 (c) β = 0.40 (d) β = 0.60 1 1 FIG. 3: Phase portraits of model (1) in the (N, E, R)-space for various values β . The algorithm developed by Marco Sandri [26] for Mathematica has been used to perform the numerical cal- culation of the Lyapunov characteristics exponents (LCE) of the NERA predator-prey model (1) for Colorado with the parameters from Tab. 2 and with β ∈ [0.02, 0.8]. As an example for β = 0.3, 0.35, 0.4 and 0.6, San- 1 1 dri’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 [19], a classiﬁcation of (autonomous) continuous-time attractors of dynamical system (1) on the basis of their Lyapunov spec- trum, together with their Hausdorﬀ 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 and involving the two algorithms pro- posed by Wolf et al. [35] and Eckmann and Ruelle [7] (see https://fr.mathworks.com/matlabcentral/ﬁleexchange/233- let). Results obtained by both algorithms are consistent. 10 TABLE IV: LCEs of NERA model (1) for Colorado for various values of β . β LCE spectrum Dynamics of the attractor Hausdorﬀ dimension 0.02 < β < 0.31 (0,−,−,−) Periodic Motion (Limit Cycle) D = 1 0.31 < β < 0.37 (0, 0,−,−) 3-Torus (Quasi-Periodic Motion) D = 2 0.37 < β < 0.55 (0,−,−,−) Periodic Motion (Limit Cycle) D = 1 0.55 < β < 0.8 (0, 0,−,−) 3-Torus (Quasi-Periodic Motion) D = 2 B. NERA model for Washington Now, we set the parameters from Tab. III where β is still the bifurcation parameter. By varying continuously β 1 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. FIG. 4: Bifurcation diagram of model (1) u as function of β . max 1 We observe from this diagram (Fig. 4) that for the value of β ∈ [0.02, 0.334] the solution of the NERA model (1) for Washington is a limit cycle (LC) in the (N, E)-plane as exempliﬁed in Fig. 5a. Then, still increasing parameter β up to the ﬁrst bifurcation value β = 0.335, the solution becomes chaotic (C) as highlighted in Fig. 5b & Fig. 5c. Such chaotic feature persists up to the second bifurcation value β = 0.357 starting from which the solution becomes again a limit cycle (see Fig 5d). Still using the algorithm developed by Marco Sandri [26] we numerically compute the Lyapunov characteristics exponents (LCE) of the NERA model (1) for Washington with the parameters from Tab. 3 and with β ∈ [0.02, 0.36]. In this case, for β = 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 classiﬁcation of attractors of dynamical system (1) on the basis of its Lyapunov spectrum, together with its Hausdorﬀ 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. 11 (a) β = 0.27 (b) β = 0.34 1 1 (c) β = 0.345 (d) β = 0.3485 1 1 FIG. 5: Phase portraits of model (1) in the (N, E, R)-space for various values β . TABLE V: LCEs of NERA model (1) for Washington for various values of β . β LCE spectrum Dynamics of the attractor Hausdorﬀ dimension 0.2 < β < 0.33 (0,−,−,−) Periodic Motion (Limit Cycle) D = 1 0.3375 < β < 0.44 (0, 0,−,−) 3-Torus D = 2 0.44 < β < 0.45 (0,−,−,−) Periodic Motion (Limit Cycle) D = 1 IV. DISCUSSION By considering that drug users are classiﬁed into four main categories: non (N), experimental (E), recreational (R) and addicts (A) users, Dauhoo et al. [6] have proposed the NERA model. Nevertheless, although this four- dimensional dynamical system took into account the mutual inﬂuence 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 modiﬁed 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 eﬀect” of the growth rate of the preys (N) and a “destabilizing eﬀect” of the predators (E, R and A) saturation. Functional responses of Verhulst and Holling type II have been used for modeling these eﬀects. Then, in order to perform numerical experiments for 12 forecasting the marijuana consumption in the states of Colorado and Washington beyond the year of the I – 502 implementation, we used data from Hanley [15] and a genetic algorithm as explained in Dauhoo et al. [6]. 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 [27]. 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 [23], we chose the parameter β , 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 β ∈ [0.02, 0.31], the solution is a limit cycle conﬁrming thus the behavior of the observed data from Hanley [15]. 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 β up to the ﬁrst bifurcation value β = 0.3175, we have shown that the solution becomes quasi periodic (3-Torus). Such chaotic attractor persists up to the second bifurcation value β = 0.37 starting from which the solution becomes again a limit cycle. When parameter β = 0.55 reaches the third bifurcation value a chaotic attractor appears again. These results have been conﬁrmed 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 β ∈ [0.02, 0.26] the solution is a limit cycle conﬁrming thus the behavior of the observed data from Hanley [15]. Then, still increasing parameter β up to the ﬁrst bifurcation value β = 0.34, the solution becomes quasi periodic (3-Torus). Such chaotic attractor persists up to the second bifurcation value β = 0.357 starting from which the solution becomes again a limit cycle. These results have been also conﬁrmed by the computation of the Lyapunov characteristics exponents. Thus, we have shown that an increase of the population of non-users (N) leads ﬁrst 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 [25]. We have thus conﬁrmed 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. Acknowledgments Matjaˇz Perc was supported by the Slovenian Research Agency (Grants J1-7009, J4-9302, J1-9112 and P5-0027) [1] A. Baveja, R. Batta, J. P. Caulkins & M. H. Karwan, Modeling the response of illicit drug markets to local enforcement, Socio-Economic Planning Sciences, 27(2), (1993) p. 73-89. [2] A. Baveja, R. Batta, J. P. Caulkins & M. H. Karwan, Collapsing street markets for illicit drugs: the beneﬁts of being decisive. Working Paper 93-39, Heinz School of Public Policy and Management, Carnegie Mellon University, Pittsburgh, PA (1993). [3] J. P. Caulkins (Ed.), Mathematical Models of Drug Markets and Drug Policy, Special issue of Mathematical and Computer Modelling 17, (1993) p. 1-115. [4] J. P. Caulkins, Local drug markets’ response to focused police enforcement, Operations Res. 41(5), (1993) p. 848-863. [5] J. P. Caulkins, R. C. Larson & T. F. Rich, Geography’s impact on the success of focused local drug enforcement operations, Socio-Economic Planning Sci. 27(1), (1993) p. 119-130. [6] M. Dauhoo, B. Korimboccus & S. Issack, On the dynamics of illicit drug consumption in a given population, IMA Journal of Applied Mathematics, 78(3), (2013) p. 432-448. [7] J.P. Eckmann & D. Ruelle, Ergodic Theory of Chaos and Strange Attractors, Rev. Mod. Phys., 57, (1985) p. 617-656. [8] H.I. Freedman & J.W.H.So, Global Stability and Persistence of Simple Food Chains, Mathematical Biosciences, 76(1) (1985) 69–86. [9] G.F. Gause, The struggle for existence, Williams and Wilkins, Baltimore, (1935). [10] G.F. Gause, V´eriﬁcations exp´erimentales de la th´eorie math´ematique de la lutte pour la vie, Hermann, Paris, (1935). [11] J.M. Ginoux, The paradox of Vito Volterra’s predator-prey model, Lettera Matematica International, Springer-Verlag, Milan (2017) p. 1-7. 13 [12] A. Gragnani, G. Feichtinger & S. Rinaldi, Dynamics of drug consumption: a theoretical model, International Insitute for Applied Systems Analysis, Laxenburg, Austria, Working Paper 94-77 (1994). [13] A. Gragnani, S. Rinaldi & G. Feichtinger, Slow-fast limit cycles in controlled drug markets, Proceedings of 3rd European Control Conference, Rome, Italy, Vol 4 (1995) p. 3031-3034. [14] A. Gragnani, S. Rinaldi & G. Feichtinger, Dynamics of drug consumption: a theoretical model, Socio-Economic Planning Sciences, 31(2), (1997) p. 127-137. [15] S. Hanley, Legalization of Recreational Marijuana in Washington: Monitoring Trends in Use Prior to the Implementa- tion of I-502, http://www.wsipp.wa.gov/ReportFile/1540/Wsipp Legalization-of-Recreational-Marijuana- in-Washington- Monitoring-Trends-in-Use-Prior-to-the-Implementation-of-I-502 Full-Report.pdf. [16] C.S. Holling, The components of predation as revealed by a study of small-mammal predation of the European pine sawﬂy, Can. Entomol. 91, (1959) p. 293-320 [17] C.S. Holling, Some characteristics of simple types of predation and parasitism, Can. Entomol. 91, (1959) p. 385-398. [18] E. Hopf, Abzweigung einer periodischen L¨osung von einer stationa¨ren L¨osung eines Diﬀerentialsystems, Berichten der Mathematisch-Physischen Klasse der S¨achsischen Akademie der Wissenschaften zu Leipzig XCIV, (1942) p. 1-22 (1942). An English translation, with comments, is included as Section 5 in Marsden et al. [1976]. [19] M. Klein and G. Baier Hierachies of dynamical systems, in A Chaotic Hierarchy, edited by G. Baier and M. Klein. Singapore: World Scientiﬁc (1991). [20] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, Berlin, 1998. [21] J.E. Marsden & M. McCracken, The Hopf Bifurcation and Its Applications, Springer-Verlag, New York, (1976). [22] R.M. May, Theoretical Ecology: Principles and Applications, Saunders, New York, (1976). [23] R.M. May, Simple mathematical models with very complicated dynamics, Nature 261, (1976) p. 459-467. [24] H. Poincar´e, H. M´ethodes Nouvelles de la M´ecanique C´eleste, Paris : Gauthier-Villars, volume I, II & III (1892-93-99). [25] M. Rosenzweig, The Paradox of Enrichment, Science 171, (1971), p. 385-387. [26] M. Sandri Numerical Calculation of Lyapunov Exponents, The Mathematica Journal, 6(3), (1996) p. 78-84. [27] F.M. Scudo & J.R. Ziegler, The Golden Age of Theoretical Ecology 1923-1940, Lecture Notes in Biomathematic, vol. 22, Springer-Verlag, Berlin-Heidelberg-New York 1978. [28] L.P. Shilnikov L.P. A case of the existence of a denumerable set of periodic motions,Sov. Math. Dokl., 6 (1965), p. 163-166. [29] L.P. Shilnikov The existence of a denumerable set of periodic motions in four-dimensional space in an extended neighborhood of a saddle-focus, Sov. Math. Dokl., 8(1) (1967), p.54-58. [30] P.F. Verhulst, Notice sur la loi que suit la population dans son accroissement, Corresp. Math. Phys., X (1838) 113–121. [31] V. Volterra, “Variazioni e ﬂuttuazioni del numero d’individui in specie animali conviventi,” Mem. Acad. Lincei III, 6 , 31–113 (1926). [32] V. Volterra, Fluctuations in the abundance of a species considered mathematically (“Nature”, vol. CXVIII, 1926 , pp. 558—560). Sotto lo stesso titolo furono poi pubblicate due lettere, una del Lotka e una del Volterra (Ibidem, vol. CXIX, 1927, pp. 12–13). [33] V. Volterra, Variation and ﬂuctuations of the number of individuals in animal species living together. Translated by Miss Mary Evelyn Wells (“Journal du Conseil international l’exploration de la mer”, Copenhague, vol. III, n. I, 1928, pp. 3–51). [34] V. Volterra, Le¸cons sur la Th´eorie Math´ematique de la Lutte pour la Vie, Gauthier-Villars, Paris, (1931). [35] A. Wolf, J.B. Swift, H.L. Swinney & J.A. Vastano, Determining Lyapunov Exponents from a Time Series, Physica D, 16, (1985) p. 285-317.
Nonlinear Sciences – arXiv (Cornell University)
Published: Jul 29, 2019
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.