Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

Chaos in a predator-prey-based mathematical model for illicit drug consumption 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 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. 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 [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 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. [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 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 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 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 [22, 27], 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 [20], 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 [30]. 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” [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 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. 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 first 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 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 αX(t)Y (t), the coefficient 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 effect 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 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 [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 influence 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 differential 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 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 I. TABLE I: Interpretation of the parameters in NERA model Parameter Sociological Meaning r Influence rate of E(t) on N(t) r Influence rate of R(t) on E(t) r Rate at which recreational users change to addicts α Influence rate of A(t) on N(t) α Influence rate of R(t) on N(t) γ Influence 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 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: O(0, 0, 0, 0) and I (1, 0, 0, 0). By posing R = A = 0, a third fixed 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 fixed 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 fixed 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 fixed point J is positive provided that: α − β > 0 and α − β (1 + h) > 0. (3) 1 4 1 4 ∗ ∗ ∗ Finally, while posing A = 0, a fifth fixed 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 fixed 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 find 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 fixed 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 difficult 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 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 α ≪ 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 benefit-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 first 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 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 [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 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 [15], 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. 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 exemplified in Fig. 3a. Then, still increasing parameter β up to the first 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 classification of (autonomous) continuous-time attractors of dynamical system (1) on the basis of their Lyapunov spec- trum, 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 and involving the two algorithms pro- posed by Wolf et al. [35] and Eckmann and Ruelle [7] (see https://fr.mathworks.com/matlabcentral/fileexchange/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 Hausdorff 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 exemplified in Fig. 5a. Then, still increasing parameter β up to the first 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 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. 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 Hausdorff 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 classified 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 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 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 confirming 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 first 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 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 β ∈ [0.02, 0.26] the solution is a limit cycle confirming thus the behavior of the observed data from Hanley [15]. Then, still increasing parameter β up to the first 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 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 [25]. 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. 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 benefits 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´erifications 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 sawfly, 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 Differentialsystems, 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 Scientific (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 fluttuazioni 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 fluctuations 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. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nonlinear Sciences arXiv (Cornell University)

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

Download PDF
13 pages

Loading next page...
 
/lp/arxiv-cornell-university/chaos-in-a-predator-prey-based-mathematical-model-for-illicit-drug-iYlDYp400o

References (38)

ISSN
0096-3003
eISSN
ARCH-3338
DOI
10.1016/j.amc.2018.10.089
Publisher site
See Article on Publisher Site

Abstract

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 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. 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 [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 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. [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 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 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 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 [22, 27], 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 [20], 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 [30]. 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” [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 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. 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 first 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 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 αX(t)Y (t), the coefficient 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 effect 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 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 [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 influence 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 differential 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 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 I. TABLE I: Interpretation of the parameters in NERA model Parameter Sociological Meaning r Influence rate of E(t) on N(t) r Influence rate of R(t) on E(t) r Rate at which recreational users change to addicts α Influence rate of A(t) on N(t) α Influence rate of R(t) on N(t) γ Influence 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 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: O(0, 0, 0, 0) and I (1, 0, 0, 0). By posing R = A = 0, a third fixed 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 fixed 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 fixed 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 fixed point J is positive provided that: α − β > 0 and α − β (1 + h) > 0. (3) 1 4 1 4 ∗ ∗ ∗ Finally, while posing A = 0, a fifth fixed 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 fixed 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 find 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 fixed 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 difficult 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 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 α ≪ 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 benefit-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 first 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 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 [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 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 [15], 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. 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 exemplified in Fig. 3a. Then, still increasing parameter β up to the first 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 classification of (autonomous) continuous-time attractors of dynamical system (1) on the basis of their Lyapunov spec- trum, 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 and involving the two algorithms pro- posed by Wolf et al. [35] and Eckmann and Ruelle [7] (see https://fr.mathworks.com/matlabcentral/fileexchange/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 Hausdorff 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 exemplified in Fig. 5a. Then, still increasing parameter β up to the first 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 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. 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 Hausdorff 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 classified 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 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 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 confirming 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 first 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 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 β ∈ [0.02, 0.26] the solution is a limit cycle confirming thus the behavior of the observed data from Hanley [15]. Then, still increasing parameter β up to the first 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 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 [25]. 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. 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 benefits 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´erifications 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 sawfly, 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 Differentialsystems, 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 Scientific (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 fluttuazioni 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 fluctuations 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.

Journal

Nonlinear Sciences – arXiv (Cornell University)

Published: Jul 29, 2019

There are no references for this article.