Analysis of illegal drug transmission model using fractional delay differential equations

1.   Introduction

Illegal drug addiction is a significant issue that threatens a large percentage of the global population, leading to diseases like heart problems, a rise in the rate of toxins in the body, persistent stomach infections, dyspepsia, psychological issues, and in some circumstances, death. There are also some intellectual consequences of illegal drug usage, such as poor looks, excessive anxiety, tension, and irritability. Illegal drug addiction is a brain disorder with behavioral and physiological features that contribute to compulsive and persistent drug use despite the adverse effects. It is one of the most pressing concerns the government wants to address worldwide. Despite efforts to curb illegal drug usage, persons who use illegal drugs increased by 20 million people between 2015 and 2016 ^[1]. According to the world health organization (WHO), 271 million individuals took illegal drugs at least once in 2017, equivalent to around 5.6 % of the global population. In terms of fatalities, the Global Burden of Disease Study estimated that 585,000 deaths and 42 million years of potential life were destroyed globally in 2017 due to illegal drugs ^[2,3].

Many researchers have explored the mathematical modeling of drug consumption ^[4,5,6,7] and a large number of people interested in continuous-time and discrete-time modelling ^[8]. To explore the propagation of heroin epidemics, Ma et al. ^[9] created many heroin epidemic models. In South Africa, Nyabadza et al. ^[10] investigated the methamphetamine transmission model. Liu et al. ^[11] studied global stability and backward bifurcation in a synthetic drug transmission model that included therapy. Saha and Samanta ^[12] investigated the drug transmission model with optimum control.

Since, a fractional-order system has a memory effect, there are several studies on fractional-order epidemiological systems ^{[13,14,15,16,17]}. Fractional calculus generalizes the order of derivatives, replacing integer order with fractional order. During a systematic investigation, it is clear that the integer-order model is a particular instance of the fractional-order model. The fractional-order system's solution must converge to the integer-order system when the order approaches one ^[18]. Fractional order systems are more appropriate in many disciplines than integer order systems. An integer order approach cannot describe the memory-related phenomena influenced by genetic factors ^[19]. The fractional-order method appears to reflect the data obtained from real-life situations better. Diethelm ^[20] compared the fractional-order systems numerical solutions to the integer-order systems, concluding that the fractional-order approach gives a better intelligible interpretation. Recently, many studies have been investigated in a fractional-order framework ^[14,21,22]. Several numerical methods are proposed to solve the fractional differential and integral equations with delay, based on the predictor-corrector of Adams-Bashforth-Moulton, finite difference, B-linear spline, cubic spline, and integro quadratic spline interpolations ^{[23,24,25,26,27,28]}. The predictor-corrector method of Adams-Bashforth-Moulton method is employed in order to reduce computations ^[23]. Also, several numerical methods are proposed to solve the fractional differential equations with the nonlinear delay with Markovian switching signals ^[29,30,31].

As recovery from illegal drug addiction takes time and the risk of relapse is high, successful rehabilitation programs should be available and long-lasting to promote stable behavior change and long-term abstinence. Scientists have also discovered that chronic drug misuse affects the structure and chemistry of the brain, and these changes can continue for months or even years after the person stops taking drugs. Volkow et al. ^[32] provided the primary example of the long-term effects of addiction on the brain and heart. Crocq et al. discussed the cultural history of man's relationship with addictive drugs, as well as the historical roots of addiction research, which deals with addictive substances and their patterns of usage throughout history ^[33].

Due to covid pandemic and lockdown, unemployment rates, assault, burglary, and drug addiction are increasing ^[34]. As a result, drug transmission is a significant danger, and it is vital to establish mechanisms to prevent drug transmission in society. To the best of our knowledge, a fractional delayed model is, therefore, provided for the first time to determine the delay between a drug user and therapy (to optimize the number of days with cost-effectiveness for a rehabilitation program). Fractional-order models have an additional parameter (order of the derivative), which is beneficial for numerical simulations. Hence, changing the order of the derivative can stabilize (destabilize) some unstable (stable) systems for particular parameter values around their equilibrium points. In a fractional-order framework, we investigated the drug transmission model with the impact of the awareness campaign's psychological treatment. A drug user should have time before relapsing into physiological addiction owing to the therapy's effects and self-control. According to Duffy's Napa Valley Rehab, thirty-five percent of individuals who are in treatment for fewer than 90 days reported relapse in the following year, compared to 17 percent of those in treatment for 90 days or more who reported relapse the next year. But, the long-term rehabilitation approach has drawbacks such as cost and conflict with job and household obligations. Motivated by the work of other dynamical systems with time delay, ^[35,36,37] a fractional delayed model is introduced to estimate the time lag between a drug-addicted person and their treatment.

In Section 2, we review some fundamental definitions and theories of fractional-order differential equations, followed by a description and formulation of the drug model in Section 3. We also discuss whether the solution of the proposed system exists, is unique, bounded, and non-negative in the next Section 4. The equilibrium points and drug abuser number are calculated in Section 5. In Section 6 local and global stability of equilibrium points (both drug-free and drug addiction) systematically is discussed. The delayed model is described in Section 7, and the stability of the delayed model is discussed in Section 8. Finally, numerical simulations are performed in Section 9 to validate the findings, followed by result and discussion; and conclusions of the whole work.

2.   Preliminaries

Fractional derivatives are robust tools for describing memory and heredity characteristics in many systems and processes. Different approaches to differentiating and integrating arbitrary orders have been proposed, each with unique characteristics. The Riemann-Liouville (R-L) and Caputo derivatives are the most commonly used definitions ^[38,39,40]. Several papers have been presented in subsequent years examining the concept and implementations of innovative fractional-order derivatives obtained by substituting the R-L and Caputo derivatives with a non-singular kernel ^[41,42]. Strong mathematical reasons demonstrate that these non-singular kernel variants have significant faults that should make them unsuitable for deployment ^[43]. The definition with non-singular kernels violates the fundamental theorem of fractional calculus ^[43]. As a result, the suggested model employs Caputo's definition, which is defined as:

Definition 1. Suppose that $\eta > 0,$ $t > 0$. Then for $n\in \mathbb{N}$

is called the Caputo fractional derivative of order $\eta$. When dealing with real-world problems, the Caputo derivative is particularly useful since it allows conventional initial and boundary conditions and the derivative of a constant is zero, which is not the case with the Riemann–Liouville fractional derivative ^[38].

Definition 2. The Mittag-Leffler function arises in the solution of fractional differential equations ^[44], introduced by G. M. Mittag-Leffler, is given by

The two-parameter Mittag-Leffler function is defined by the series expansion:

where $\eta > 0$ and $\gamma > 0$. G. M. Mittag Leffler developed the Laplace transform of the Mittag-Leffler function ^[45]:

The Laplace transform of the Caputo fractional derivative is defined as ^[39]

Lemma 1. Consider the fractional differential equation

with $0 < \eta\leq1$ and $y:[t_0, \infty]\times\Omega\rightarrow \mathbb{R}^n, \Omega \in \mathbb{R}^n$. Then a unique solution of above equation on $[t_0, \infty]\times\Omega$ exists provided $y(t, x)$ obeys the locally Lipchitz condition with respect to $x$ ^[38].

Lemma 2. (Generalized Mean Value Theorem ^[46]) Let $f(t) \in C[a, b]$ and $\dfrac{d^\eta f(t)}{dt^\eta} \in C[a, b]$ for $0 < \eta \leq 1$, then

Remark 1. If $f(t) \in C[0, b]$ and $\dfrac{d^\eta f(t)}{dt^\eta} \in C[0, b]$ for $0 < \eta \leq 1$. It is clear that if $\dfrac{d^\eta f(t)}{dt^\eta} \geq 0 \; \forall \; t \in (0, b]$, then $f(t)$ is non-decreasing and if $\dfrac{d^\eta f(t)}{dt^\eta} \leq 0\; \forall\; t \in (0, b]$, then $f(t)$ is non-increasing for all $t \in [0, b]$ ^[46].

3.   Formulation of mathematical model

According to the substance use status, the illegal drug addiction model divides the entire population into four divisions. $S(t)$ stands for susceptible and refers to someone who is in danger of abusing any substance (illegal drug). Persons who use the substance (illegal drug) in any form occasionally are categorized as drug users, $I(t)$, individuals who use illegal drugs regularly in any form are defined as a drug abuser, $A(t)$, and those who quit using illegal drugs either via abstinence or rehabilitation are labeled as recovered, $R(t)$. Individuals enter the susceptible population at a rate of $\Pi$ through births and migration, according to the model. Susceptible people become drug users as a result of their interactions with drug users. The contact rate per capita is $\beta$ in this scenario. The natural recovery rate of drug users is $\delta$, while the rate at which people relapse to drug usage is $\rho$. The natural death rate is represented by $\mu$, while the rate at which all drug users abuse drugs is $\sigma$. The death rates of drug users and drug abusers are $r$ and $\theta$, respectively. All variables and parameters considered in the proposed model are defined in Table 1. We formulate the following model, which is a system of nonlinear fractional order differential equations with non-negative initial conditions, based on some assumptions and the flow diagram (Figure 1):

where $\eta \in (0, 1]$ is order of derivative for the proposed model, due to which memory property is included in proposed model. In the L.H.S. of the system, the human population is in fractional time with the dimension $t^{-\eta}$. Every constant has power $\eta$ in the R.H.S. of the system to maintain the system dimensionally balance since death rate, birth rate, etc. always have dimension $t^{-1}$. Hence, the use of parameter, $\eta$, instead of integer parameter can apparently lead to better results since one has an extra freedom degree.

4.   Model analysis

First, we have verified the existence and uniqueness of the solution and then demonstrate that it is non-negative and bounded.

4.1.   Existence and uniqueness

Theorem 1. There is a unique solution of model (3.1) along with non-negative initial conditions for $t \geq 0$. Furthermore, all the solutions are bounded.

Proof. Let $N(t)$ be the total population at any time, then

By taking Laplace transformation of Eq (4.1):

Now, by taking inverse Laplace transformation:

It is known that \(E_\eta(-x)\) is completely monotonic for \(x > 0\) if \(0 < \eta \leq 1\) and \(0\leq E_\eta(-\mu^\eta t^\eta)\leq 1\) ^[44,47].

According to Lemma 1, the solution of the proposed model (3.1) along with initial conditions for $t > 0$ exists, unique and bounded.

4.2.   Positivity of solution

As population cannot be negative, so, this section demonstrates the positivity of solutions. Consider that the condition (positivity of solution) is failed at some point to demonstrate $(S, I, A, R) \geq 0$. Let

Since

Then,

The solution of the proposed model is non-negative from Remark 1 of Lemma 2.

4.3.   Invariant set

The set $\Omega$ is invariant if for all the initial conditions are in $\Omega$, the solution of model (3.1) remains in $\Omega$. As a consequence, a positively invariant set will have positive solutions.

Hence, from Eqs 4.2 and 4.3

is positively invariant region.

5.   Equilibrium points

● Illegal drug-free equilibrium: When there are no illegal drug users (illegal drug users and drug abusers).

● Illegal drug addiction equilibrium: When illegal drug users are not zero.

5.1.   Illegal drug-free equilibrium

The illegal drug-free equilibrium $(E_0)$ is determined when $I = 0; A = 0; R = 0$. Therefore, the illegal drug free steady state is given by $\left(\frac{\Pi^\eta}{\mu^\eta}, 0, 0, 0\right)$.

5.2.   Drug abuser number

The drug abuser number is defined as the number of newly addicted persons generated by a single addicted individual throughout the contagious period ($\Lambda = 10$ indicates a person with drug addiction would transmit it to an average of 10 additional people). The largest eigenvalue of $\zeta^*\xi^*{^{-1}}$ at $E_0$ is used to calculate the drug abuser number $\Lambda$ of given system ^[48].

where $\xi(t) = \xi_i^-(t)-\xi_i^+(t)$ and matrices $\zeta_i(t), \xi_i^-(t), \xi_i^+(t)$ are given by

At $E_0$, the Jacobian matrix of $\zeta(t)$ is given by

The Jacobian matrix of $\xi(t)$ is

The drug abuser number $\Lambda$ is obtained as

5.3.   Illegal drug addiction equilibrium

The illegal drug addiction equilibrium $(E_1)$ is determined when $I \neq 0; A \neq 0$. Therefore, the drug addiction steady state is given by $(S_1, I_1, A_1, R_1)$, where

For the existence of positive steady state $\Lambda$ must be greater than or equal to one. If $\Lambda = 1$ drug addiction steady state become drug free steady state. So, the drug addiction steady state is exists iff $\Lambda > 1$.

6.   Stability analysis

Theorem 2. The illegal drug free equilibrium point $E_0 = \left(\frac{\Pi^\eta}{\mu^\eta}, 0, 0, 0\right)$ of the system (3.1) is globally asymptotically stable under some restriction on parameters when $\Lambda < 1$; otherwise it is unstable.

Proof. Jacobian matrix at $E_0$ is

Clearly, one of the eigen-values is $-\mu^\eta$ lies in 2nd quadrant.

Now, for checking stability system we have to check rest of eigen-values. So for this we consider,

The equation of characteristic of above matrix is given by

By using Routh-Hurwitz criteria ^[49], if $P_1 > 0, P_3 > 0, P_1P_2-P_3 > 0$; then all the eigenvalues of the matrix will have negative real parts. Hence, the system (3.1) is stable whenever $\Lambda < 1$ ^[50]. Here $(I(t), A(t), R(t)) \to \left(0, 0, 0\right)$ as $t \to \infty$, so, $(S(t), I(t), A(t), R(t)) \to \left(\frac{\Pi^\eta}{\mu^\eta}, 0, 0, 0\right)$ as $t \to \infty$. So $E_0$ is globally asymptotically stable for $\Lambda < 1$ whenever $P_1 > 0, P_3 > 0, P_1P_2-P_3 > 0$.

Theorem 3. The illegal drug addiction steady state $E_1 = (S_1, I_1, A_1, R_1)$ exists and is locally asymptotically stable iff $\Lambda > 1$.

Proof. For illegal drug addiction equilibrium point $E_1 = (S_1, I_1, A_1, R_1)$ the Jacobian matrix of the given model (3.1) at is obtained as follows:

The characteristic equation corresponding to $J_1$ is

The proposed system is stable if the characteristics equation has roots with negative real parts. The following conditions are required for the stability of the system ^[49]:

As a result, $E_1 = (S_1, I_1, A_1, R_1)$ is locally asymptotically stable if the above conditions are met.

Theorem 4. The illegal drug addiction steady state $E_1 = (S_1, I_1, A_1, R_1)$ is globally asymptotically stable iff $\Lambda > 1$.

Proof. Let the positive definite Lyapunov function be defined by:

Now, the $\eta^{th}$ order derivative of Lyapunov function is given by ^[51]:

According to Gallegos and Duarte-Mermoud ^[52], the drug addiction steady state is globally asymptotically stable in $\Omega$ if $L < M$.

7.   Delayed model

An illegal drug user in treatment should have some time before relapsing into physiological addiction owing to the therapy's effects and self-control. Hence, a fractional delayed model is introduced to estimate the time lag between a drug-addicted person and their treatment. The fractional-order mathematical model subject to non-negative initial condition in the Banach space ${C}([-\tau, 0], {R}_{+}^{4})$ of continuous functions mapping the interval $[-\tau, 0]$ to ${R}_{+}^{4}$, is as follows:

All variables and parameters considered in the proposed model are defined in Table 1.

8.   Stability analysis of delayed model

The endemic equilibrium point of Eq (7.1) is $(S_1, I_1, A_1, R_1)$, and equations for $(S_1, I_1, A_1, R_1)$ are provided in Section 5. The effect of gradual increases in time delay on stability was investigated in this paper. The characteristic equation for system (7.1) at the equilibrium point $(S_1, I_1, A_1, R_1)$ is,

After solving, we get

The transcendental equation is known to have an unlimited number of complex roots. It is quite difficult to discern the indications of the roots when there is a time delay ($\tau$). As a result, we begin our research by determining that the temporal delay ($\tau$) is zero, and then develop stability criteria for $\tau > 0$.

At $\tau = 0$, the Eq (8.1) becomes Eq (6.1), and the condition of stability for Eq (6.1) is described in Section 6. Hence,

Theorem 5. The equilibrium point $(S_1, I_1, A_1, R_1)$ of system (7.1) is globally asymptotically stable for $\tau = 0$ iff condition (6.2) satisfied.

If $\tau > 0$, then the equation of characteristic (8.1) will have an indefinite number of roots. The existence of entirely imaginary solutions to the Eq (8.1) is a necessary for the variation in stability of $E_1$. Let $i\nu_0$ be the root of the Eq (8.1), then

Squaring and adding above two equations,

If $N_4 < 0;$ there must be two completely imaginary roots $i\nu_0$ of the Eq (8.3) exists. The Eq (8.2) is used to obtained $\tau$,

By taking derivative of the Eq (8.1) with respect to $\tau$,

Let us consider

$= \{G_1+G_2\}$, where

${\rm{If }}\ G_1, G_2 \ {\rm{are\ positive\ then}} \; \; \left[Re\left (\dfrac{d\lambda}{d\tau}\right)^{-1}\right]_{\lambda = i\nu_0} > 0$ which gives near $\tau = \tau_{\eta}$, there is a bifurcation of periodic solutions from $E_1 = (S_1, I_1, A_1, R_1)$. Hence,

Theorem 6. Suppose that $G_1, G_2 > 0$ for system (7.1), then $E_1 = (S_1, I_1, A_1, R_1)$ is asymptotically stable for $\tau \in [0, \tau_{\eta})$ and unstable when $\tau > \tau_{\eta}$ and exhibiting Hopf bifurcation at $E_1 = (S_1, I_1, A_1, R_1)$ when $\tau = \tau_{\eta}$.

9.   Numerical simulation

We exhibit numerical results in this part to show the dynamic nature of illegal drug propagation model, and to evaluate the analytical conclusions for various values of derivatives with time delays. The proposed model is solved using predictor-corrector method of Adams-Bashforth-Moulton which is implemented in MATLAB ^[23]. The variables and parameters used for evaluation are shown in Table 2 which were retrieved from ^[53]. It reveal that the equilibrium points for the proposed model are asymptotically stable for $\eta = 0.7, 0.8, 0.9, 1$, as shown in Figure 2. As $\eta$ is raised, the population converges quicker to its equilibria, as seen in Figure 2. The suggested model's stability is demonstrated by the trajectory of all populations, irrespective of order.

The behaviour of delayed model (7.1) for $\eta = 0.7, 0.8, 0.9, 1$ is shown in Figures 3–6. The plots show that the time delay for different order of derivative $\eta = 0.7, 0.8, 0.9, 1$, and demonstrate that the system is stable for $\tau < \tau_{\eta}$. The value time delay coefficient for the different order of derivative $\eta$ are shown in Table 3. Figures clearly indicate that when the order of derivative grows, the delay reduces, showing that the model's stability zone widens. According to Figure 7 drug abuser number increases as order ($\eta$). As the rate at which people relapse to drug use $(\rho)$ increases the drug abuser number increases which means the drug user population increases.

10.   Result and disscussion

According to Duffy's Napa Valley Rehab, thirty-five percent of individuals in treatment for fewer than 90 days reported relapse in the following year, compared to 17 percent of those in treatment for 90 days or more who reported deterioration the next year. But, the long-term rehabilitation approach has drawbacks such as cost and conflict with job and household obligations. To the best of our knowledge, a fractional delayed model is, therefore, provided for the first time to determine the delay between a drug user and therapy (to optimize the number of days with cost-effectiveness for a rehabilitation program).

There are two equilibrium points for the proposed model: illegal drug-free and illegal drug addiction. The drug abuser number ($\Lambda$) is also obtained for the proposed model. The illegal drug-free equilibrium point is stable whenever $\Lambda < 1$, whereas the drug addiction equilibrium point exists and is globally asymptotically stable whenever $\Lambda > 1$. The proposed delayed model exhibits Hopf bifurcation at the drug addiction equilibrium point. The predictor-corrector method of Adams-Bashforth-Moulton to solve fractional differential equations with or without delay is employed in order to reduce computations. The different plots for delayed model show that the time delay for different order of derivative $\eta = 0.7, 0.8, 0.9, 1$, and demonstrate that the system is stable for $\tau < \tau_{\eta}$.

According to our research, an illicit drug user can decrease drug intake in 9.7 weeks. On the other hand, the integer-order model predicts that an illegal drug user may reduce drug intake in 7.9 weeks. Since the success rate of a long-term rehab program is high in reality, a fractional-order model has a greater success rate than an integer-order model. Most rehab programs are 30, 60, and 90 days. Since the success rate of 90 days rehab program is high, in contrast, the long-term rehabilitation approach also has some disadvantages such as cost and conflict with job and household obligations. Therefore, the current study concludes that treatment programs lasting 70 to 75 days are more successful in reducing drug addiction.

11.   Conclusions

This research presents an illegal drug transmission model of fractional differential equations to explain the dynamics of illegal drug transmission better while considering drug users' histories. It is clear from this study that if the drug abuser number $(\Lambda)$ is less than unity, drug propagation can be regulated, and the illegal drug-free equilibrium can be maintained worldwide asymptotically. The drug abusers number increases as the value of the order of derivative $(\eta)$ grows to a specific value of $\eta$. The number of people who continue to use illegal drugs occasionally and regularly rises over time so does the susceptible population decrease with time. The people who have decided to either quit using drugs or avoid them due to education or self-awareness of the consequences of using these substances, as well as the anti-drug campaign is one of the causes for the reduction or avoidance of illegal drug usage might be this. So, the recovered population increases. The fractional-order illegal drug transmission mathematical model includes the delay coefficient $(\tau)$. The delayed model's results indicate that after the delay reaches a certain point, the model will oscillate repeatedly. The delay reduces as the order of derivatives grows from 0 to 1, showing that the model's stability zone widens. In this case, we find that an illegal drug user can reduce drug consumption in 9.7 weeks, whereas the integer-order model predicts that an illegal drug user can reduce drug consumption in 7.9 weeks. Because the success rate of a long-term rehab program is high in reality, the success rate of a fractional-order model is higher than that of an integer order model.

The proposed model can be generalized in the future by taking into account a variety of different elements. Different approaches may be used to get the best value for the fractional order of derivatives. Real-life data from various places may be used to confirm the simulation results.

Acknowledgments

Komal Bansal is grateful to the CSIR for providing financial assistance with 1091/(CSIR-UGC NET DEC. 2018). Along with that author want to thanks INTI International University, Malaysia for guiding and supporting in providing resources for this research.

Conflict of interest

We declare that we have no financial or personal relationships with other people or organisations that could improperly influence our work, and that we have no professional or other personal interest in any product, service, or company that could be interpreted as influencing the position presented in the manuscript.