key: cord-0062110-c57mwc7s
authors: Kundu, Soumen; Kumari, Nitu; Kouachi, Said; Kundu, Piu
title: Stability and bifurcation analysis of a heroin model with diffusion, delay and nonlinear incidence rate
date: 2021-04-20
journal: Model Earth Syst Environ
DOI: 10.1007/s40808-021-01164-x
sha: c100c6a501280353413fdcb0b07c5fd817fa7fb1
doc_id: 62110
cord_uid: c57mwc7s

As we all know, the use of heroin and other drugs in Europe and more specifically in Ireland and the resulting prevalence are well documented. A huge population is still dying using heroin every day. This may happen due to, several reasons like, excessive use of painkiller, lack of awareness etc. It has also inspired mathematical modelers to develop dynamical systems predicting the use of heroin in long run. In this work, the effect of heroin in Europe has been discussed by constructing a suitable mathematical model. Our model describes the process of treatment for heroin users by consolidating a sensible utilitarian structure that speaking to the restricted accessibility of treatment. In the treatment time frame, because of the discretion of the medication clients, some kind of time delay called immunity delay might be found. The effect of immunity delay on the system’s stability has been examined. The existence of positive solution and its boundedness has been established. Also, the local stability of the interior equilibrium point has been studied. Taking the immunity delay as the key parameter, the condition for Hopf-bifurcation has been studied. Using normal form theory and center manifold theorem, we have likewise talked about the direction and stability of delay induced Hopf-bifurcation. The corresponding reaction diffusion system with Dirichlet boundary condition has been considered and the Turing instability has been studied. Obtained solutions have also been plotted by choosing a suitable value of the parameters as the support of our obtained analytical results.

Addiction to drugs is an extremely common phenomenon among drug abusers. According to the world drug report 2018, global heroin and opium production reached an alltime high and drug-related deaths increased by 60% between 2000 and 2015 (Djilali et al. 2017) . Especially, in recent years, expanded utilization of heroin is an issue of concern in many parts of the world because that it can cross the blood-brain barrier within 15-20 s, which make the heroin users more prone to catch the human immunodeficiency virus (HIV) and the other blood-borne viruses diseases (Liu et al. 2019a, b; Yang et al. 2016b) . The extended use of heroin habituation and addiction in population which is named as the heroin epidemic presents some phenomena of epidemics such as rapid diffusion and clear geographic boundaries . Therefore, mathematical modeling has been universally used to describe heroin addiction in the epidemic way since 1981-1983 in Ireland (Kelly et al. 2003) .

The incidence rate is an important factor for epidemic mathematical modeling. There are some literature for heroin epidemic models with the bilinear incidence rate (Liu et al. 2019b; Fang et al. 2015; Huang and Liu 2013; Liu and Zhang 2011; Wang et al. 2011; Wangari and Stone 2017; Zhang and Wang 2019; Maji 2021) . In 1979, an exponential model has been simplified from Kermack and McKendrick's disease model by MacKintosh and Stewart (1979) . Their model MacKintosh and Stewart (1979) describes the extent of heroin use over a large area in an epidemic fashion. An ODE has been formulated by White and Comiskey (2007) to study the dynamics of disease modeling to the drug-using career. For this, the population has been divided into three categories, namely: susceptible, heroin users, and heroin users undergoing treatment. Wang et al. (2011) proposed a heroin epidemic model with bilinear incidence rate and investigated the stability of drug-free equilibrium and the endemic equilibrium. Liu and Zhang (2011) and Huang and Liu (2013) investigated a heroin epidemic model with bilinear incidence rate and distributed delays, respectively. Liu and Zhang (2011) showed that the basic reproduction number characterizes the disease transmission dynamics. Huang and Liu (2013) investigated the global stability of the heroin epidemic model by using the constructing appropriate Lyapunov function which improved the obtained results by Liu et al. in the literature Liu and Zhang (2011) . However, as stated in Fang et al. (2015) , all the heroin epidemic models did not consider the influence of the treat-age for heroin users during the treatment. Based on this fact, Fang et al. (2015) formulated a heroin epidemic model bilinear incidence rate and treat-age and studied the global asymptotic properties of the proposed model. Wangari and Stone (2017) proposed a heroin epidemic model with bilinear incidence rate and saturated treatment function and derived a condition ensuring global stability of the model by making use of the geometric approach. Considering that it usually takes some time for a heroin user to be susceptible. Zhang and Wang (2019) developed a heroin model with time delay and saturated treatment function based on the model proposed by Wangari and Stone (2017) . They have derived the sufficient conditions for local stability and existence of Hopf bifurcation of the model by taking the time delay due to the period used to cure heroin users as the bifurcation parameter and investigated the direction and stability of the Hopf bifurcation by using the normal form theory and the manifold center theorem.

The bilinear incidence rate is more appropriate for communicable diseases such as H3N2 and SARS. Based on the limitation of the heroin epidemic models above, there have some heroin epidemic models with nonlinear incidence rates been proposed in recent years. Yang et al. (2016a) proposed a heroin epidemic model with a nonlinear incidence rate and investigated the stability of the equilibria for the model. Considering the time needed to return to an untreated drug user, Liu and Wang (2016b) studied a delayed multi-group heroin epidemic model with relapse phenomenon and nonlinear incidence rate, and derived the sufficient conditions for the global stability of the equilibria with the help of suitable Lyapunov functionals. Djilali et al. (2017) presented the global dynamics of a heroin epidemic model with a very general nonlinear incidence rate. Recently, Ma et al. (2017) developed the following heroin model with nonlinear incidence rate as following:

where S(t), U 1 (t) and U 2 (t) denote numbers of the susceptible population, heroin users not in treatment and heroin users in treatment at time t, respectively. b is the recruitment rate of the susceptible population; SU 2 1 is the infection rate of the heroin users not in treatment, which means that a susceptible person will become a heroin-addict once he contacts with heroin users twice; d is the natural death rate of the population; treatment has been taken by the heroin users at a rate ; sometimes heroin users neglect the treatment and reuse heroin ( ); after the treatment the heroin users becomes the susceptible ones at the rate of unit. Ma et al. (2017) showed that system (1) exhibits numerous kinds of bifurcation and the results obtained have certain effect to control the heroin spreading. The schematic diagram of our model (1) has been shown in Fig. 1 .

As stated in Ma et al. (2017) , 'once someone addicts to the drug, it is difficult to abstain from it just depending on self-control because of the natural facts'. Therefore, there should be a period for heroin users in treatment before they become susceptible ones since they have certain self-control ability. In general, delay differential equations show more complex dynamics rather than ordinary differential equations, which have been shown in some work about traditional epidemic models (Liu and Wang 2016a; Sirijampa et al. 2018; Tipsri and Chinviriyasit 2015; Haldar et al. 2020; Chhetri et al. 2020) and computer virus propagation models (Keshri and Mishra 2014; Ren et al. 2012 ; Upadhyay and Kumari 2019; Zhao and Bi 2017). Motivated

(1)

The schematic diagram of the system (1) by the above, we consider the following heroin epidemic model with time delay:

where is the immunity delay due to the self-control of the drug users in treatment. Our object of this work is to study the effect time delay on the dynamics of our assumed model. The model has extremely complicated kinetic features, which are mainly reflected in this article. The paper is organized into several sections and subsections. In "Basic results", it has been shown that solutions of system (2) are positive and bounded for a positive initial condition. The conditions for the existence of Hopf-bifurcation and local stability ("Local stability") of the interior equilibrium point have been derived in "Hopf bifurcation". Also, in this section, taking the delays as key parameters, different cases have been discussed for the occurrence of Hopf-bifurcation. Further, in "Properties of Hopf bifurcation", the direction of Hopf-bifurcation and the stability of the periodic solution are examined. In "The reaction diffusion system", we considered the reaction diffusion system associated to system (2), where in "Local existence and positivity of solutions", we proved the local existence and positivity of solutions with Dirichlet boundary conditions and in "Global existence of solutions", we proved the global existence, in time of the solutions. In "Turing instability (DDI)", we studied the Turing instability of the reaction diffusion system but with homogeneous Neumann conditions. Some numerical simulations have been done for our expository results in "Numerical simulation". The paper ends with a discussion of our work ("Conclusions").

The initial conditions for the model system (2) are:

be the space of functions which are continuous from [− , 0] to ℝ 3 + with supremum norm, called Banach Space. By the elementary theory of functional differential equation, we can say that system (2) with initial conditions (3) admits a unique solution.Now, to show the system (2) has positive and bounded solutions, we state the following result:

(2)

Theorem 1 All solutions of system (2) with initial conditions (3) are positive and bounded for t ⩾ 0.

Proof It can easily be shown that all the solutions of (2) with initial conditions (3) are defined on ℝ 3 + and remain positive ∀ t ≥ 0. By adding the equations of system (2), we get

Straightforward computation shows that for b sufficiently large, then system (2) has two positive equilibriums E * (S * , U 1 * , U 2 * ) , where and U 2 * is the positive root of Eq. (5) where, which can be written So the condition to get two positive critical points can be obtained for large values of b (from Eq. (5)):

The subject of bifurcation is a vast area to do research. In a system of differential equations Hopf-bifurcation happens when the complex conjugate set of eigenvalues of a linearised system become purely imaginary at a fixed point.

,

So with these conditions, a system that contains at least two equations may have Hopf bifurcation. Hopf-bifurcation occurs at a point where a system changes state from stable to unstable, i.e. it is a local bifurcation in which a fixed point of a dynamical system loses stability. In this section by taking the delay as bifurcation parameters the necessary conditions for the Hopf-bifurcation are derived.

Local stability of an equilibrium point means that the equilibrium is stable to small perturbations, i.e., if we push the situation a bit out from the equilibrium point then the situation will return back on its own. Here we shall study the local stability with the help of Routh-Hurwitz criteria (Kundu et al. 2017 ). The characteristic equation of the linear section of system (2) at E * (S * , U 1 * , U 2 * ) is with and

0 is satisfied based on the the Routh-Hurwitz criteria, i.e., we have

A 0 =a 33 (a 12 a 21 − a 11 a 22 ) + a 11 a 23 a 32 ,

A 1 =a 11 a 22 + a 11 a 33 + a 22 a 33 − a 12 a 21 − a 23 a 32 ,

The stability condition for ≠ 0 has been discussed in the next sub-section.

For ≠ 0 , we assume that = i ( > 0) is a root of Eq. (6). Then separating real and imaginary part we have,

Suppose that a 0 < 0 , then Eq. (9) has at least one positive root say v 0 , i.e., Eq. (8) has a positive root 0 = √ v 0 and Eq. (7) has a pair of purely imaginary roots ±i 0 . For 0 , we have Differentiating Eq. (7) with respect to , we get

Based on the discussion above and the Hopf bifurcation theorem in Upadhyay and Kumari (2019), we have the following results. (2), if the conditions (H) hold, then the positive equilibrium E * (S * , U 1 * , U 2 * ) is locally asymptotically stable when ∈ [0, 0 ) ; system (2) undergoes a Hopf bifurcation at E * (S * , U 1 * , U 2 * ) when = 0 and a family of periodic solutions bifurcate from E * (S * , U 1 * , U 2 * ).

Introduce a new perturbation parameter = − 0 with ∈ R , then = 0 is the Hopf bifurcation value of

and normalize by t → (t∕ ) . Then, system (2) becomes the following equation in the phase space

A and B are defined by the following form:

and where a 13 = −2 U 1 * , a 14 = − S * , a 24 = 2 U 1 * , a 25 = S * and the form of other elements in the matrix A and B are given before the Eq. (7). From the Riesz representation theorem, there exists ( , ) of bounded variance such that Choosing (10) Based on discussion above, we know that i 0 are eigenvalues of A(0) and A * (0) . Suppose that ( ) = (1, 2 , 3 ) T e i 0 and * (s) = V(1, * 2 , * 3 )e i 0 s are the corresponding eigenfunctions. By direct calculation we obtain From ⟨ * , ⟩ = 1 , we have

Using the algorithm given in Hassard et al. (1981) and the similar computation process to that in Guo et al. (2019) , Sirijampa et al. (2018) and Li (2017) , we can obtain the expressions of g 20 , g 11 , g 02 and g 21 as follows:

with E 1 and E 2 can be solved by Thus, we have where 2 determines the direction of the Hopf bifurcation;

2 determines the stability of the bifurcating periodic solutions and T 2 determines the period of the bifurcating periodic solutions. In conclusion, we can obtain the following results based on the fundamental results about Hopf bifurcation in the literature Hassard et al. (1981) . (2), if 2 > 0 ( 2 < 0 ), then the Hopf bifurcation is supercritical (subcritical); if 2 < 0 ( 2 > 0 ), then the bifurcating periodic solutions are stable (unstable); g 20 =2V 0 [a 13 2 + a 14 2 2 +̄ * 2 (a 24 2 + a 25 2 2 )], g 11 =V 0 [a 13 ( 2 +̄2) + 2a 14 2̄2 +̄ * 2 (a 24 ( 2 +̄2) + 2a 25 2̄2 )],

a 11 a 12 b 13 a 21 a 22 a 23 0 a 32 a 33 + b 33

, then the period of the bifurcating periodic solutions increase (decrease).

We consider the following reaction diffusion system associated to the reaction system (1) with the Dirichlet boundary conditions

where Ω is an open bounded domain in ℝ N with smooth boundary Ω . The positive constants D 1 , D 2 and D 3 are called the diffusion terms. The initial data are assumed to be nonnegative and uniformly bounded on Ω.

Since the nonlinearities are continuously differentiable on ℝ 3 + , then for any initial data in ∞ Ω , it is easy to check directly their Lipschitz continuity on bounded subsets of the domain of a fractional power of the operator I 3 D 1 , D 2 , D 3 t Δ , where I 3 the three dimensional identity matrix, Δ is the Laplacian operator and (.) t denotes the transposition. Under these assumptions, the following local existence result is well known (see Smoller 2012; Henry 2006) .

The system (15)-(17) admits a unique, classical local solution (S, U 1 ,

Since the reactions are quasi-positive, i.e.

then the nonnegativity of the solutions is preserved by application of classical results on invariant regions.

To prove global existence of system (15) with Dirichlet boundary conditions (16) and positive initial data (17), we apply an old general result of Morgan (1989) for m-components systems on the form (15) 

We say that a system of the form

exhibits Turing instability, or DDI, if the system without diffusion, i.e. has locally stable equilibrium state which becomes unstable in the presence of diffusion. Here the function F (we assume it is regular) describes the reaction dynamics and D is a diagonal matrix of diffusion coefficients (see Elragig 2013) . Diffusion-driven instability is considered one of the mechanisms implicated in the spatial organization in morphogenesis Turing (1992) .

In studying pattern formation in RD systems the key first step is to determine the Turing space for a given model, i.e. the parameter set for the mode on which pattern formation can be triggered Murray (1982) . This can then be followed by bifurcation analysis of specific pattern formations Hoyle and Hoyle (2006) . Turing instability, or diffusion driven instability (DDI), is a concept first proposed by Turing (1952) .

We suppose that the boundary conditions are of Neumann and homogenous. Imposing such boundary conditions is due to their neutral nature as they do not pump the space with any additional material and this makes "selforganization" plausible. Taking other boundary conditions can influence the predictions where this can drive forming different patterns, (see Murray 2003) .

To analyse DDI mathematically, we use linearized stability analysis and Fourier transform together to get an ordinary differential (reaction) system on the form where A describes the Jacobian matrix of F (i.e. the linearized reaction matrix) and supposed to be stable (i.e. all its eigenvalues have negative real parts) and û is the Fourier transform of u in space.

Qian and Murray (2001) , the authors presented a simple, practical, necessary and sufficient condition for diffusion-driven linear instability and parameter space determination in nonlinear reaction systems with three species is presented:

(i) the largest diagonal elements of A is greater than zero with corresponding diffusion coefficient very small; or

(ii) the smallest diagonal cofactors of A is less than zero with corresponding diffusion coefficient very large.

Let us apply Proposition (4), to system (15) with Neumann homogenous boundary conditions.The Jacobian of F at the critical point E * (S * , U 1 * , U 2 * ) is where S * , U 1 * , U 2 * are given by (4). The largest diagonal element of A is a 22 = 2 SU 1 * − ( + d) which is positive under the following condition that is (21) for sufficiently small diffusion D 2 .

Example 1 For b = 1 , = 0.003 , d = 0.01 , = 0.2 , = 0.005 , = 0.001 , the above inequality is satisfied for For the Large diffusions, Straightforward computation shows that if the cofactors of the first and third diagonal elements of A are respectively less then zero under the following condition and respectively. With some tedious but elementary algebra, we can show that the two above inequalities are analogous to (21) and are satisfied under analogous conditions. Applying Proposition 4, we have

Theorem 5 Each of the following conditions will lead to diffusion driven instability Condition-1: (22) and D 1 ≫ D 2 and D 1 ≫ D 3 , or Condition-2: (23) and D 3 ≫ D 1 and D 3 ≫ D 2 .

0.125 22 < U 1 * < 72. 321.

In the previous sections, the conditions for stability and Hopf-bifurcation of the system (2) (2), we conclude that E * is seen to be stable for less value of the delay i.e., ∈ [0, 0 ) (Fig. 2) . Now with the increased values of delays the oscillation of the individual also increases. Therefore, at the critical value of delay i.e., = 0 , we get a stable periodic solution where the Hopf-bifurcation occurs and there we get a stable limit cycle around the equilibrium point E * (Fig. 3) . For large value of delay ( > 0 ) the system becomes unstable (Fig. 4) . The one parameter bifurcation diagram for all the individuals with respect to the delay has been plotted in Figs. 5, 6 and 7 respectively.

We have performed numerical solution for the spatially extended system of Eq. (24). We here have shown the existence of Hopf bifurcation in the reaction-diffusion system numerically with respect to delay parameter . The hopf bifurcation occurs at 0 = 7.375 . Stable dynamics has been observed after solving the above system. In Fig. 8 , (2) for < 0 when E 1 * equilibrium point is stable. Parameters are written in the text the reaction-diffusion system shows stable dynamics at = 5 < 0 = 7.375 . On the other hand for the increased value of time delay = 9 > 0 = 7.375 , the system admits oscillatory dynamics as shown in Fig. 9 .

As we all know that in Europe and more specifically in Ireland heroin and other drugs are mostly used, the report has been mentioned in Kelly et al. (2003) , Fang et al. (2015) , Huang and Liu (2013) , Liu and Zhang (2011 ), Wang et al. (2011 ), Wangari and Stone (2017 , Zhang and Wang (2019) , Maji (2021) , White and Comiskey (2007) , Yang et al. (2016a) , Elragig (2013 ), Friedman (2008 , Guo et al. (2019) , Hassard et al. (1981 ), Henry (2006 , Hoyle and Hoyle (2006) , Keshri and Mishra (2014) , Li (2017) , Liu and Wang (2016a) , Sirijampa et al. (2018) , Tipsri and Chinviriyasit (2015) , Haldar et al. (2020) , Chhetri et al. (2020) , Liu and Wang (2016b) , Ma et al. (2017 ), Morgan (1989 , Mulone and Straughan (2009 ), Murray (1982 , Qian and Murray (2001) , Ren et al. (2012) , Upadhyay and Kumari (2019) , Zhao and Bi (2017) , Sekiguchi and Ishiwata (2010) , Smoller (2012 ), Sporer (1999 , Turing (1952, Wiessing and Hartnoll (1999) . At that place, people are very much addicted to drugs and it causes many preventable deaths. In fat cells, heroin is well liquefiable, for that reason, it crosses the blood-brain barrier rapidly ('within 15-20 s') and "affects the brain and central nervous system, which causes both the 'rush' experience by users and the harmfulness" (Kelly et al. 2003) . The death related to heroin is also related to the use of alcohol or other drugs (Sporer 1999). Nowadays heroin becomes as popular among the population that its treatment becomes typically a heavy one. The dynamics of infectious disease has been studied using the techniques used in mathematics and statistics; see, for example, Mulone and Straughan (2009) , Sekiguchi and Ishiwata (2010) , White and Comiskey (2007) and Zhang and Teng (2008) ; however, little has been done to apply this method to the heroin epidemics.

In this work, a heroin model with the incident rate in the form of SU 2 1 is discussed. In our model, we consider a susceptible person to turn into a heroin someone who is addicted when he/she is enticed twice. Also, a heroin-client without treatment can't enter the susceptible population unexpectedly. This hypothesis may be more in line with the truth. Comparing with the works in Wang et al. (2019) and Ma et al. (2017) , completely different dynamical properties are discovered in this paper according to discussions.

In the present study, we not only investigate the effect of the time delay parameter on the model system dynamics but also present the complete analysis of the model system including the positivity of solutions, boundedness, stability and the properties related to the Hopf bifurcation. In particular, we have ensured the existence and uniqueness of solutions of the proposed model using elementary theory of functional differential equations. Following the existence and uniqueness of solutions, the boundedness of solution has also been established using comparison theorem. Sufficient condition for the local stability of interior equilibrium point have been discussed by considering for different values of the delay parameter. Furthermore, the sufficient conditions for the local stability of the interior equilibrium point have been investigated for various range of delay parameter. Along with this, in each case, the occurrence of Hopf bifurcation has been obtained by considering the time delay as a bifurcation parameter and the critical value of the delay parameter has also been obtained. From the analysis, we conclude that if the values of time delay for system (2) is less that it's corresponding critical value then the equilibrium for the system remains in a stable state while if the value of the delay is greater than its corresponding critical value, then the equilibrium for system (2) losses its stability. The properties of Hopf bifurcation , such as direction and stability have also been studied by means of the center manifold and normal form theory. Our result, however, indicates that when the recruitment rate of the susceptible population (i.e. b) is large enough then the system has two steady-states which are stable and can also lose, each one, its stability (bifurcation in the case of the reaction system and Turing instability when the diffusion are added) due to diffusion in the following three cases:

-For sufficiently small diffusion for U 1 with respect to those of S and U 2 , or -For sufficiently large diffusion for S with respect to those of U 1 and U 2 , or -For sufficiently large diffusion for U 2 with respect to those of U 1 and S.

For graphical verification of obtained analyticaly results and to observe the effects of the immunity delay on the stability of the system, some numerical computations have been presented by taking particular set of parametric values. We have systematically justifies all the analytical conditions. As we know that the educational awareness programs related to treatment may play a major role in protecting the individuals against use of drug. Therefore, the present study provides an opportunity for further extension of the work via incorporating of public awareness factor for treatment and to observe the effect of public awareness on the eradication of using heroin and other drugs. 

Influence of gestation delay and the role of additional food in holling type III predator-prey systems: a qualitative and quantitative investigation

A heroin epidemic model: very general non linear incidence, treat-age, and global stability

On transients, lyapunov functions and turing instabilities

Global asymptotic properties of a heroin epidemic model with treat-age

Partial differential equations of parabolic type

Hopf bifurcation analysis in a predatorprey model with time delay and food subsidies

Modeling and analysis of a predator-prey type eco-epidemic system with time delay

CUP Archive Henry D (2006) Geometric theory of semilinear parabolic equations

Pattern formation: an introduction to methods

A note on global stability for a heroin epidemic model with distributed delay

Prevalence of opiate use in Ireland 2000-2001: a 3-source capture-recapture study. Stationery Office

Two time-delay dynamic model on the transmission of malicious signals in wireless sensor network

Qualitative analysis for a delayed three species predator-prey model in presence of cooperation among preys

A study on time-delay rumor propagation model with saturated control function

Hopf bifurcation of a delayed siqr epidemic model with constant input and nonlinear incidence rate

Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate

Global behaviour of a heroin epidemic model with distributed delays

Dynamics of a stochastic heroin epidemic model

Dynamics of a stochastic heroin epidemic model with bilinear incidence and varying population size

Bifurcation of a heroin model with nonlinear incidence rate

A mathematical model of a heroin epidemic: implications for control policies

Dynamical analysis of a fractional-order predator-prey model incorporating a constant prey refuge and nonlinear incident rate

Global existence for semilinear parabolic systems

A note on heroin epidemics

Parameter space for turing instability in reaction diffusion mechanisms: a comparison of models

A simple method of parameter space determination for diffusion-driven instability with three species

A delayed computer virus propagation model and its dynamics

Global dynamics of a discretized sirs epidemic model with time delay

Hopf bifurcation analysis of a delayed seir epidemic model with infectious force in latent and infected period

Shock waves and reaction-diffusion equations

The effect of time delay on the dynamics of an seir model with nonlinear incidence

Collected works of am turing: morphogenesis. Saunders

Philosophical the royal biological transactions society sciences

Discrete and data packet delays as determinants of switching stability in wireless sensor networks

Dynamics of a heroin epidemic model with very population

Analysis of an age-structured multigroup heroin epidemic model

Analysis of a heroin epidemic model with saturated treatment function

Heroin epidemics, treatment and ode modelling

Study to obtain comparable national estimates for problem drug use prevalence for all EU member states

Global dynamics of a heroin epidemic model with age structure and nonlinear incidence

Global dynamical analysis of a heroin epidemic model on complex networks

Global behavior and permanence of sirs epidemic model with time delay

Hopf bifurcation of a heroin model with time delay and saturated treatment function

Hopf bifurcation of a computer virus spreading model in the network with limited anti-virus ability

Acknowledgements We thank the anonymous referee for valuable suggestions. One of the authors Mrs. Piu Kundu wishes to thank DST-INSPIRE India for their financial support.

The authors declare that they do not have any conflicts of interest.