Stability and bifurcation analysis of a heroin model with diffusion, delay and nonlinear incidence rate.

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.

Introduction

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 all-time 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 (Wang et al. 2019). 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), and 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; 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.

Fig. 1

The schematic diagram of the system (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 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”).

Basic results

Consider . The initial conditions for the model system (2) are:with and where be the space of functions which are continuous from to 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:

Theorem 1

All solutions of system (2) with initial conditions (3) are positive and bounded for

Proof

It can easily be shown that all the solutions of (2) with initial conditions (3) are defined on and remain positive By adding the equations of system (2), we get and So, Hence and i.e. N(t) is bounded, namely Hence are bounded.

Straightforward computation shows that for b sufficiently large, then system (2) has two positive equilibriums , whereand is the positive root of Eq. (5)where,which can be writtenSo the condition to get two positive critical points can be obtained for large values of b (from Eq. (5)):

Stability and bifurcation analysis of system (2)

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

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 iswithandWhen , Eq. (6) becomeswhereNow for , is locally asymptotically stable if the condition (H): , , is satisfied based on the the Routh–Hurwitz criteria, i.e., we haveThe stability condition for has been discussed in the next sub-section.

Hopf bifurcation

For , we assume that is a root of Eq. (6). Then separating real and imaginary part we have,Thus we can obtainwhereLet , then Eq. (8) becomesSuppose that , then Eq. (9) has at least one positive root say , i.e., Eq. (8) has a positive root and Eq. (7) has a pair of purely imaginary roots . For , we haveDifferentiating Eq. (7) with respect to , we getFurther,where .

Obviously, if and , then . Based on the discussion above and the Hopf bifurcation theorem in Upadhyay and Kumari (2019), we have the following results.

Theorem 2

For system (2), if the conditions (H) hold, then the positive equilibrium is locally asymptotically stable when ; system (2) undergoes a Hopf bifurcation at when and a family of periodic solutions bifurcate from .

Properties of Hopf bifurcation

Introduce a new perturbation parameter with , then is the Hopf bifurcation value of system (2). Let , , , and normalize by . Then, system (2) becomes the following equation in the phase space where , : , F: , withA and B are defined by the following form:andwhere , , , 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 thatChoosingand is the Dirac delta function. For , define the operator andThen, system (2) is equivalent toFor , define and a bilinear productwith .

Based on discussion above, we know that are eigenvalues of A(0) and . Suppose that and are the corresponding eigenfunctions. By direct calculation we obtainFrom , we haveUsing 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 , , and as follows:with and can be solved byThus, we havewhere determines the direction of the Hopf bifurcation; determines the stability of the bifurcating periodic solutions and 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).

Theorem 3

For system (2), if (), then the Hopf bifurcation is supercritical (subcritical); if (), then the bifurcating periodic solutions are stable (unstable); if (), then the period of the bifurcating periodic solutions increase (decrease).

The reaction diffusion system

We consider the following reaction diffusion system associated to the reaction system (1)with the Dirichlet boundary conditionswhere is an open bounded domain in with smooth boundary . The positive constants and are called the diffusion terms. The initial dataare assumed to be nonnegative and uniformly bounded on .

Local existence and positivity of solutions

Since the nonlinearitiesare continuously differentiable on , 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 , where the three dimensional identity matrix, is the Laplacian operator and denotes the transposition. Under these assumptions, the following local existence result is well known (see Smoller 2012; Henry 2006).

Proposition 1

The system (15)–(17) admits a unique, classical local solution on . If then

Since the reactions are quasi-positive, i.e.then the nonnegativity of the solutions is preserved by application of classical results on invariant regions.

Global existence of solutions

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) which we summarize in our case (i.e. ) as follows

Proposition 2

We suppose that the functions and h are of polynomial growth and satisfyfor all where ; and are positive reels. Then all solutions of system (15) with positive Dirichlet boundary conditions and positive initial data are positive and global in time (i.e.).

Since for our system all conditions of the above Proposition (2) are satisfied except (20) which can be verified, using the positivity of the solutions for all times, as followsThen we have the following proposition

Proposition 3

All solutions of system (15) with Dirichlet boundary conditions (16) and positive initial data (17) are positive and global in time (i.e. ).

Turing instability (DDI)

We say that a system of the formexhibits 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 “self-organization” 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 formwhere 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:

Proposition 4

A sufficient condition for diffusion-driven instability is that either

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

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 iswhere are given by (4). The largest diagonal element of A is which is positive under the following conditionthat is

Theorem 4

A sufficient condition for diffusion-driven instability for the reaction diffusion system (15) is given by (21) for sufficiently small diffusion .

Example 1

For , , , , , , 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 conditionandrespectively. 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 and , or Condition-2: (23) and and

Numerical simulation

In the previous sections, the conditions for stability and Hopf-bifurcation of the system (2) have been derived analytically. To illustrate our obtained analytical results in a better way, some numerical simulations have been performed. The parameters are chosen in such a way so that it satisfies the conditions for the existence of the Hopf-bifurcation. The set of parametric values are as follows: , , , , , . Then,from which we obtain the positive equilibrium . With this set of parameters, we can verify that the conditions (H) and are satisfied by some complex computations and obtain and the critical value of delay . Based on Theorem (2), we conclude that is seen to be stable for less value of the delay i.e., (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., , we get a stable periodic solution where the Hopf-bifurcation occurs and there we get a stable limit cycle around the equilibrium point (Fig. 3). For large value of delay () 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.

Fig. 2

Time series solution of the system (2) for when equilibrium point is stable. Parameters are written in the text

Fig. 3

Time series solution of the system (2) for when equilibrium point is stable periodic. Parameters are written in the text

Fig. 4

Time series solution of the system (2) for when equilibrium point has increasing oscillation. Parameters are written in the text

Fig. 5

One parameter bifurcation diagram for S with respect to

Fig. 6

One parameter bifurcation diagram for with respect to

Fig. 7

One parameter bifurcation diagram for with respect to

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 . Stable dynamics has been observed after solving the above system. In Fig. 8, the reaction-diffusion system shows stable dynamics at . On the other hand for the increased value of time delay , the system admits oscillatory dynamics as shown in Fig. 9.

Fig. 8

Hopf bifurcation in the reaction–diffusion system with delay, when , , , , , , with time lag

Fig. 9

Hopf bifurcation in the reaction–diffusion system with delay, when , , , , , , with time lag

Conclusions

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, 2003), 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, 1992), 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 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 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.

For sufficiently small diffusion for with respect to those of S and , or

For sufficiently large diffusion for S with respect to those of and , or

For sufficiently large diffusion for with respect to those of and S.

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.