7.2
CiteScore
3.7
Impact Factor
Generic selectors
Exact matches only
Search in title
Search in content
Post Type Selectors
Search in posts
Search in pages
Filter by Categories
ABUNDANCE ESTIMATION IN AN ARID ENVIRONMENT
Case Study
Correspondence
Corrigendum
Editorial
Full Length Article
Invited review
Letter to the Editor
Original Article
Retraction notice
REVIEW
Review Article
SHORT COMMUNICATION
Short review
7.2
CiteScore
3.7
Impact Factor
Generic selectors
Exact matches only
Search in title
Search in content
Post Type Selectors
Search in posts
Search in pages
Filter by Categories
ABUNDANCE ESTIMATION IN AN ARID ENVIRONMENT
Case Study
Correspondence
Corrigendum
Editorial
Full Length Article
Invited review
Letter to the Editor
Original Article
Retraction notice
REVIEW
Review Article
SHORT COMMUNICATION
Short review
View/Download PDF

Translate this page into:

Full Length Article
10 2024
:36;
103329
doi:
10.1016/j.jksus.2024.103329

Global stability and modeling with a non-singular kernel for fractional order heroin epidemic model: Insights from different population studies

School of Physical and Mathematical Sciences, Faculty of Exact and Natural Sciences, Pontificia Universidad Católica del Ecuador, Av. 12 de Octubre 1076, Apartado, Quito 17-01-2184, Ecuador
Research Centre on Public Health (CESP), Department of Medicine and Surgery, University of Milan-Bicocca, Monza, Italy
Department of Mathematics, College of Science, King Khalid University, P.O. Box 9004, Abha 61413, Saudi Arabia
Department of Mathematics, University of the Punjab, Quaid-e-Azam Campus, Lahore, Pakistan
Department of Mathematics, Near East University, TRNC, Mersin 10, Nicosia 99138, Turkey
Department of Mathematical Sciences, University of South Africa, Florida 0003, South Africa
Department of Computer Science and Mathematics, Lebanese American University, 1107-2020, Beirut, Lebanon

⁎Corresponding author. talatn@unisa.ac.za (Talat Nazir)

Disclaimer:
This article was originally published by Elsevier and was migrated to Scientific Scholar after the change of Publisher.

Abstract

There has been a worldwide epidemic of heroin that has affected people, families, societies, and cultures across the world. Now, the heroin epidemic has transitioned from heroin abuse to the overuse of synthetic narcotics, which are widely accessible and inexpensively produced. In this work, a novel mathematical approach is applied to investigate the dynamics of the heroin epidemic model and its harmful effect on society with different population data. A heroin model has been constructed with the importance of a non-singular kernel in the sense of a generalized Mittag-Leffler kernel. The well-posedness of the proposed model is proven via fixed-point theory. To examine the heroin model, two equilibrium states have been determined. These equilibrium states are proven to be locally and globally asymptotically stable. To analyze the behavior of heroin, a basic reproduction number and sensitivity analysis are used to determine the impact of different parameters mathematically as well as through simulations. To find the approximate solution, we implement the Toufik–Atangana numerical method at different fractional order values. The sensitivity of the heroin model is carried out, and 3-D graphs show the significance of the parameter involved in the model. Finally, the numerical outcomes are presented with different values of fractional parameters.

Keywords

Heroin epidemic model
Atangana–Baleanu fractional derivative
Sensitivity analysis
Stability
Numerical solution

Data availability

Our manuscript has no associated data.

1

1 Introduction

The escalating intake of drugs and other dangerous substances is a severe concern. Heroin addiction affects not just the level of life for the vast majority of people but also, to a greater extent, the overall picture of world peace and economic prosperity (Xu, 2023; Glicksberg et al., 2023). Based on data from the United Nations (U.N.) World Drug Report, 35 million people have severe problems related to drug misuse, and only around a seventh are receiving a cure (Fabien, 2018). In addition, it was disclosed that drug addiction occurred globally in 2017, with 3 million cases of HIV and 6.2 million cases of severe hepatitis B viral infections. Metric records indicate that the detrimental physiological effects of drug use are more widespread and substantial than previously thought and that controlling the prevalence of harmful drugs is essential (Saha and Samanta, 2019; Tolomeo et al., 2021).

An indispensable mechanism for comprehending and treating drug addiction challenges is mathematically modeled (Raza et al., 2023). It has come to light that models are beneficial tools for forecasting drug users’ behavior patterns and offering practical recommendations for therapeutic approaches (Chou and D’Orsogna, 2022; Raza et al., 2022b). Narcotic drug misuse, including heroin consumption, is already a major societal issue. Wangari and Stone (2017) have investigated the spread of heroin addiction in society with the saturated treatment function. Zhang and Xing (2020) have formulated a reaction–diffusion heroin disease model with stability analysis. The authors in Huang and Liu (2013) proposed a distributed delay model of heroin and examined the dynamics in the population. In Liu and Wang (2016), the authors have modeled the heroin model with the saturating rate of non-linear occurrence and provided conditions for the global disease dynamics. Khan et al. (2021) have discussed preferable control factors for heroin spread with sensitivity.

A noteworthy advancement has recently been made in fractional calculus, using new derivative and integral operators with kernels (Butt et al., 2023a,c; Zafar et al., 2021). The creatively suggested element takes advantage of a modified version of the Mittag-Leffler function (MLF), this structure encompasses the keystone and its properties, which play a vital role in developing innovative methods to attain various fascinating characteristics, like the expanding variations and mean square deformation interphase that are found in significant scenarios. The authors in Zafar et al. (2021) developed a heroin model using a fractional operator. Weera et al. (2023) have considered a fractional heroin model for the numerical solutions. The groundbreaking fractional derivative operator, developed by Atangana and Baleanu in 2016, has widespread use in several scientific and technological fields. It has been demonstrated that utilizing the AB-fractional operator in simulation produces a disordered framework. for a brief time (Syam and Al-Refai, 2019; Bas and Ozarslan, 2018; Gomez-Aguilar et al., 2016; Owolabi and Atangana, 2019). The AB-fractional derivative is an excellent mathematical method for modeling progressively complex critical difficulties in the setting of Caputo since it has now been found that the MLF is a more significant and effective filtering procedure than the exponential and power laws.

In this work, we develop an epidemic model of heroin comprising addiction and enduring immunity, then explore the worldwide patterns triggered by the previous arguments. We used the new ABC-fractional method and Toufik–Atangana’s approach (Mekkaoui and Atangana, 2017) to create a fractional representation and evaluate numerical results. To the researchers knowledge, a few studies have used the ABC-fractional operator to analyze the heroin epidemic scenario. Additionally, the researchers claim that there is an absence of appropriate research that examines robust management of computational structures from the standpoint of ABC-fractional forms.

The paper’s structure is as follows: In Section 2, we delve into the dynamics of heroin, providing detailed explanations and introducing fundamental concepts in fractional calculus. Moving on to Section 3, we demonstrate the model’s existence and uniqueness using mathematical approaches. Section 4 is dedicated to analyzing the model’s stability analytically. Section 5 covers the numerical method and its corresponding results. Finally, Section 6 offers insights into the research’s overall conclusions.

2

2 Formulation and configuration of heroin model

This section outlines the creation of the heroin disease epidemic model (Raza et al., 2022a), and three subgroups comprise the whole population. U ( t ) represents susceptible persons, V ( t ) shows drug users, and W ( t ) displays non-drug users, which are the state variables in the model, respectively. This model considers the following assumptions:

  • The total population N is the sum of all the subclasses and its size is assumed to be fixed in the whole period.

  • All the biological parameters must be nonnegative or positive.

  • Homogeneous mixing is considered.

  • All of the individuals are supposed to be identically susceptible to drug habit.

Therefore, the heroin epidemic’s expansion is visualized in Fig. 1, and it is elucidated through the accompanying set of ODEs:

(1)
d U d t = η σ 1 U V N δ 1 U , d V d t = σ 1 U V N ρ V + σ 2 V W N δ 1 V , d W d t = ρ V σ 2 V W N δ 1 W , subject to initial conditions U 0 0 , V 0 0 , W 0 0 and N = U + V + W .
The heroin epidemic’s expansion within sub populations.
Fig. 1
The heroin epidemic’s expansion within sub populations.

Due to the complex structure of the population, we equalize the system (1) by subsidizing the state variables: S = U N , X = V N , Y = W N as follows

(2)
d S ( t ) d t = η σ 1 S X δ 1 S , t 0 , d X ( t ) d t = σ 1 S X ρ X + σ 2 X Y δ 1 X , t 0 , d Y ( t ) d t = ρ X σ 2 X Y δ 1 Y , t 0 , with initial conditions S 0 0 , X 0 0 , Y 0 0 and S ( t ) + X ( t ) + Y ( t ) = 1 . An elaborate explanation of the parameters employed in the study, along with a description of their scientific significance and relationship to certain physical occurrences to be provided in Table 1.

The traditional derivative neglects memory effects that exist in numerous biological structures. Therefore, We employ the Atangana–Baleanu fractional derivative in place of the conventional integer-order derivative in order to enlarge framework (2) for the present study. This transformation will enable us to track memory effects and get further insights into the dynamics of the heroin epidemic. Before continuing, we revisit key aspects of Atangana–Baleanu and Caputo derivatives (Butt et al., 2023b; Abdeljawad and Baleanu, 2017; Rashid et al., 2022; Butt et al., 2022).

Table 1 Interpretation of parameters of the heroin model.
Parameter Interpretation
η New drug users
σ 1 Probability rate of infection
σ 2 Non-drug consumers have a high probability of becoming drug consumers.
δ 1 Individuals who quit using drugs as a habit
δ 2 Rate of drug-related deaths
δ 3 Death rate for drug-related deaths following treatment
ρ Drug users’ transition to becoming non-users

Definition 1

Consider w 1 ( a , b ) , b > a , be a mapping and then, the Atangana–Baleanu fractional derivative in Caputo sense is expressed (Butt et al., 2023b) as:

(3)
a A B C D t α w ( t ) = M ( α ) 1 α a t w ̇ ( ξ ) E α [ α ( t ξ ) α 1 α ] d ξ , where, 0 α 1 and M ( α ) = 1 α + α Γ ( α ) . E α ( ν ) shows the Mittag-Leffler function and expressed as follows:
(4)
E α ( ν ) = β = 0 ν β 1 + α β , α , β .

Definition 2

The AB-fractional integral of the function w 1 ( a , b ) is written by Butt et al. (2023b):

(5)
a A B I t α w ( t ) = 1 α M ( α ) w ( t ) + α M ( α ) Γ ( α ) a t w ( ξ ) ( t ξ ) α 1 d ξ .

Lemma 1

For w 1 ( a , b ) , then both the AB-fractional derivative and integral operator satisfies the Newton–Leibniz identity (Abdeljawad and Baleanu, 2017):

(6)
a A B I t α a A B C D t α w ( t ) = w ( t ) w ( a ) .

Lemma 2

For two functions w , u Δ ( a , b ) , b > a , then the AB-fractional operator satisfies the subsequent variant (Abdeljawad and Baleanu, 2017; Rashid et al., 2022):

(7)
a A B C D t α w ( t ) a A B C D t α u ( t ) = Δ w ( t ) u ( t ) .

Therefore, the system (2) can be extended to a fractional framework of order α ( 0 , 1 ] , with constant transmission rates. Thus, the Atangana–Baleanu fractional derivative defines our proposed model:

(8)
0 A B C D t α S ( t ) = η σ 1 S X δ 1 S , 0 A B C D t α X ( t ) = σ 1 S X ρ X + σ 2 X Y δ 1 X , 0 A B C D t α Y ( t ) = ρ X σ 2 X Y δ 1 Y . The domain of the proposed model (8) is described by the set Ξ as:
(9)
Ξ : = { S ( t ) , X ( t ) , Y ( t ) R + 3 : 0 S ( t ) + X ( t ) + Y ( t ) η δ 1 } .
It is now sufficient to prove that the set of records mentioned in Tables 1 and (9) are positively constant. The preceding lemma will serve as the basis for the proof of (9).

Lemma 3

(Abdeljawad and Baleanu, 2017; Rashid et al., 2022) Suppose that w 1 ( a , b ) and assume 0 A B C D t α w ( t ) 1 ( a , b ) , α ( 0 , 1 ] . Then, we have w ( t ) = w ( a ) + 1 Γ ( α ) 0 A B C D t α w ( t ) ( z a ) α , when t [ 0 , z ] .

It is noticed that by Lemma 3, if w ( t ) [ 0 , δ ] , 0 A B C D t α w ( t ) [ 0 , δ ] and 0 A B C D t α w ( t ) 0 , for all t ( 0 , δ ] , α ( 0 , 1 ] , then the function w ( t ) is increasing, and if 0 A B C D α w ( t ) 0 , for all t ( 0 , δ ] , then the function w ( t ) is decreasing for all t ( 0 , δ ] .

Utilizing Lemma 3, we show the accumulation remains positively invariant, we have

(10)
0 A B C D t α S | S = 0 = η 0 , 0 A B C D t α X | X = 0 = 0 0 , 0 A B C D t α Y | Y = 0 = ρ X ( t ) 0 . All the solutions to (8) are positive and contained in R + 3 due to (10), and the set specified in (9) is positively invariant for framework. In addition, we proceed by adding all framework equations to illustrate the bounded nature of the solutions in the fractional system (8), that gives
(11)
0 A B C D t α N ( t ) = η δ 1 N .
Applying Laplace transform, we get L 0 A B C D t α N ( t ) + δ 1 N L ( η ) , L N [ ( 1 σ ) s α σ α 1 α ] s α 1 N ( 0 ) 1 α M ( α ) [ s α + α 1 α ] η s , which can be written as L N 1 σ α ( 1 α ) ( 1 σ ) s α 1 × { 1 α ( 1 σ ) M ( α ) [ 1 α + α s α 1 α ] η s + N ( 0 ) ( 1 σ ) s } , where σ = δ 1 ( 1 α ) M ( α ) .

On using the inverse Laplace transform, yields

(12)
N ( t ) η δ 1 η δ 1 ( 1 σ ) d d t 0 t E α σ α ( 1 α ) ( 1 σ ) ( t x ) α d x + 1 ( 1 α ) E α σ α ( 1 α ) ( 1 σ ) t α N ( 0 ) . Here, the MLF is represented as E α 1 , α 2 . Taking into account the notion that the MLF has asymptotic properties, thus E α 1 , α 2 ( β ) q = 1 θ β q / Γ ( α 2 α 1 q ) + O ( | 1 β 1 + θ | ) , for | β | , α 1 π 2 < | a r g ( β ) | π . N ( ζ ) η / δ 1 is easy to grasp since ζ tends to . As a consequence, (9) demonstrates the biological viability in domain of model (8).

3

3 Existence and uniqueness of heroin model

Model (8) provides a framework for analysis that incorporates the dynamics of the heroin epidemic, and the ABC-fractional derivative system may be represented as

(13)
0 A B C D t α S ( t ) = Ω 1 ( t , S ) , 0 A B C D t α X ( t ) = Ω 2 ( t , X ) , 0 A B C D t α Y ( t ) = Ω 3 ( t , Y ) , where the kernels are defined as
(14)
Ω 1 ( t , S ) = η σ 1 S X δ 1 S , Ω 2 ( t , X ) = σ 1 S X ρ X + σ 2 X Y δ 1 X , Ω 3 ( t , Y ) = ρ X σ 2 X Y δ 1 Y ,
supplements with ICs S ( 0 ) = S 0 , X ( 0 ) = X 0 , Y ( 0 ) = Y 0 .

Below is a discussion of the alternatives for fractional-order systems, including their existence and uniqueness of solutions. We employ the renowned Banach fixed-point theorem (Banach, 1922) to establish the existence of a solution to problem (8). For an in-depth examination of some useful results of fixed point and compressions, you can turn to Ref. Panda (2020), Şahin et al. (2023), Şahin and Alagöz (2023), Aggarwal et al. (2023).

Theorem 1

Theorem 1 Banach, 1922

Let ( X , d ) be a complete metric space and mapping f : X X satisfies d ( f x , f y ) α d ( x , y ) , for all x , y X , where α [ 0 , 1 ) is a constant. Then there exists a unique x X such that x = f x . Moreover, for any x 0 X , the iterative sequence x n + 1 = f x n converges to x .

We shall now show the existence and uniqueness of the solution by following the procedures listed below. The formulation (8) is handled using the AB-fractional integral:

(15)
S ( t ) S ( 0 ) = 1 α M ( α ) Ω 1 ( t , S ) + α M ( α ) Γ ( α ) 0 t Ω 1 ( ζ , S ) ( t ζ ) α 1 d ζ , X ( t ) X ( 0 ) = 1 α M ( α ) Ω 2 ( t , X ) + α M ( α ) Γ ( α ) 0 t Ω 2 ( ζ , X ) ( t ζ ) α 1 d ζ , Y ( t ) Y ( 0 ) = 1 α M ( α ) Ω 3 ( t , Y ) + α M ( α ) Γ ( α ) 0 t Ω 3 ( ζ , Y ) ( t ζ ) α 1 d ζ . Let us assume that for any Banach space, H ( I ) = [ 0 , t ] of real-valued continuous functions on I , there exists an interval I = [ 0 , t ] such that W = H ( I ) × H ( I ) × H ( I ) . The corresponding norm of these mappings is as described below: S = sup t I | S ( t ) | , X = sup t I | X ( t ) | , Y = sup t I | Y ( t ) | .

The contraction and the Lipschitz assumption provide the cornerstones of our subsequent theorem.

Theorem 2

If the kernels are Ω i , i = 1 , 2 , 3 in (8), then there exists ξ i , i = 1 , 2 , 3 , such that

(16)
Ω 1 ( t , S ) Ω 1 ( t , S ˆ ) ξ 1 S ( t ) S ˆ ( t ) , Ω 2 ( t , X ) Ω 2 ( t , X ˆ ) ξ 2 X ( t ) X ˆ ( t ) , Ω 3 ( t , Y ) Ω 3 ( t , Y ˆ ) ξ 3 Y ( t ) Y ˆ ( t ) are contraction mappings for ξ i [ 0 , 1 ) , i = 1 , 2 , 3 .

Proof

We shall commence from the initial category denoted as S . In this context, let S and S ˆ represent two mappings. It is essential to consider the following factors:

(17)
Ω 1 ( t , S ) Ω 1 ( t , S ˆ ) = σ 1 S X δ 1 S ( σ 1 S ˆ X δ 1 S ˆ ) , = σ 1 X ( S S ˆ ) δ 1 ( S S ˆ ) , ( σ 1 m 1 + δ 1 ) S S ˆ , ξ 1 S S ˆ , where ξ 1 = ( σ 1 m 1 + δ 1 ) , S = sup t I | S ( t ) | = m 1 , X = sup t I | X ( t ) | = m 2 , Y = sup t I | Y ( t ) | = m 3 . Therefore, S satisfies the Lipschitz condition, and it is a contraction if 0 ( σ 1 m 1 + δ 1 ) < 1 . For the next scenarios, the Lipschitz criteria are given in the same manner.

Now for t = t n , n = 1 , 2 , , we present the iterative relationship of (15) that follows:

(18)
S n ( t ) ( 0 ) = 1 α M ( α ) Ω 1 ( t , S n 1 ) + α M ( α ) Γ ( α ) 0 t Ω 1 ( ζ , S n 1 ) ( t ζ ) α 1 d ζ , X n ( t ) ( 0 ) = 1 α M ( α ) Ω 2 ( t , X n 1 ) + α M ( α ) Γ ( α ) 0 t Ω 2 ( ζ , X n 1 ) ( t ζ ) α 1 d ζ , Y n ( t ) ( 0 ) = 1 α M ( α ) Ω 3 ( t , Y n 1 ) + α M ( α ) Γ ( α ) 0 t Ω 3 ( ζ , Y n 1 ) ( t ζ ) α 1 d ζ , along with appropriate starting condition S ( 0 ) = S 0 , X ( 0 ) = X 0 , Y ( 0 ) = Y 0 .

To determine the difference between the subsequent components, utilize the following equations:

(19)
1 n ( t ) = S n ( t ) S n 1 ( t ) = 1 α M ( α ) ( Ω 1 ( t , S n 1 ) Ω 1 ( t , S n 2 ) ) + α M ( α ) Γ ( α ) 0 t ( Ω 1 ( ζ , S n 1 ) Ω 1 ( ζ , S n 2 ) ) ( t ζ ) α 1 d ζ ,
(20)
2 n ( t ) = X n ( t ) X n 1 ( t ) = 1 α M ( α ) ( Ω 2 ( t , X n 1 ) Ω 2 ( t , X n 2 ) ) + α M ( α ) Γ ( α ) 0 t ( Ω 2 ( ζ , X n 1 ) Ω 2 ( ζ , X n 2 ) ) ( t ζ ) α 1 d ζ ,
(21)
3 n ( t ) = Y n ( t ) Y n 1 ( t ) = 1 α M ( α ) ( Ω 3 ( t , Y n 1 ) Ω 3 ( t , Y n 2 ) ) + α M ( α ) Γ ( α ) 0 t ( Ω 3 ( ζ , Y n 1 ) Ω 3 ( ζ , Y n 2 ) ) ( t ζ ) α 1 d ζ .
Applying the norm on (19)(21), we obtain
(22)
1 n ( t ) = S n ( t ) S n 1 ( t ) = 1 α M ( α ) Ω 1 ( t , S n 1 ) Ω 1 ( t , S n 2 ) + α M ( α ) Γ ( α ) 0 t Ω 1 ( ζ , S n 1 ) Ω 1 ( ζ , S n 2 ) ( t ζ ) α 1 d ζ ,
(23)
2 n ( t ) = X n ( t ) X n 1 ( t ) = 1 α M ( α ) Ω 2 ( t , X n 1 ) Ω 2 ( t , X n 2 ) + α M ( α ) Γ ( α ) 0 t Ω 2 ( ζ , X n 1 ) Ω 2 ( ζ , X n 2 ) ( t ζ ) α 1 d ζ ,
(24)
3 n ( t ) = Y n ( t ) Y n 1 ( t ) = 1 α M ( α ) Ω 3 ( t , Y n 1 ) Ω 3 ( t , Y n 2 ) + α M ( α ) Γ ( α ) 0 t Ω 3 ( ζ , Y n 1 ) Ω 3 ( ζ , Y n 2 ) ( t ζ ) α 1 d ζ .
Furthermore, Eq. (22) can be transformed into the following expressions: 1 n ( t ) = S n ( t ) S n 1 ( t ) , 1 α M ( α ) Ω 1 ( t , S n 1 ) Ω 1 ( t , S n 2 ) + α M ( α ) Γ ( α ) 0 t Ω 1 ( ζ , S n 1 ) Ω 1 ( ζ , S n 2 ) ( t ζ ) α 1 d ζ , 1 α M ( α ) L 1 S n 1 S n 2 + α M ( α ) Γ ( α ) 0 t ξ 1 S n 1 S n 2 ( t ζ ) α 1 d ζ , 1 α M ( α ) + t α M ( α ) Γ ( α ) ξ 1 1 ( n 1 ) ( t ) . As a result, we get
(25)
1 n ( t ) 1 α M ( α ) + t α M ( α ) Γ ( α ) ξ 1 1 ( n 1 ) ( t ) .
In this regard, the subsequent inequalities can be used to simplify all other representations of (23)(24):
(26)
2 n ( t ) 1 α M ( α ) + t α M ( α ) Γ ( α ) ξ 2 2 ( n 1 ) ( t ) ,
(27)
3 n ( t ) 1 α M ( α ) + t α M ( α ) Γ ( α ) ξ 3 3 ( n 1 ) ( t ) .

The theorem is now described in detail as follows.

Theorem 3

The suggested heroin fractional model (8) has exact solutions if the following suppositions exist: that is, we can workout P 0 α such that

(28)
1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ξ i < 1 , i = 1 , 2 , 3 .

Proof

Utilizing Eqs. (25)(27), we have

(29)
1 n ( t ) S ( 0 ) [ ( 1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ) ξ 1 ] n , 2 n ( t ) X ( 0 ) [ ( 1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ) ξ 2 ] n , 3 n ( t ) Y ( 0 ) [ ( 1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ) ξ 3 ] n . This establishes the continued existence of the previously indicated strategies. Moreover, we proceed as directed to confirm that the function stated above is a consequence of (8):
(30)
S ( t ) S ( 0 ) = S n 1 a 1 n ̃ ( t ) , X ( t ) X ( 0 ) = X n 1 a 2 n ̃ ( t ) , Y ( t ) Y ( 0 ) = Y n 1 a 3 n ̃ ( t ) .
Therefore, we have
(31)
a 1 n ̃ ( t ) 1 α M ( α ) Ω 1 ( t , S n ) Ω 1 ( t , S n 1 ) + α M ( α ) Γ ( α ) 0 t Ω 1 ( ζ , S n ) Ω 1 ( ζ , S n 1 ) ( t ζ ) α 1 d ζ 1 α M ( α ) ξ 1 S n S n 1 + t α M ( α ) Γ ( α ) ξ 1 S n S n 1 .
By continuing the process recursively, gives
(32)
a 1 n ̃ ( t ) 1 α M ( α ) + t α M ( α ) Γ ( α ) n + 1 ξ 1 n S n S n 1 n .
Setting t α = P 0 α yields a 1 n ̃ ( t ) 1 α M ( α ) + P 0 α M ( α ) Γ ( α ) n + 1 ξ 1 n S n S n 1 n . Considering a 1 n ̃ ( t ) 0 . Applying limit as n to (35), it is noted that a 1 n ̃ ( t ) 0 for
(33)
1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ξ 1 < 1 .
In a similar way, we show that a 2 n ̃ ( t ) 0 , and a 3 n ̃ ( t ) 0 , then
(34)
1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ξ i < 1 , i = 1 , 2 , 3 .

The existence of fractional model (8) is guaranteed by Theorems 2 and 3 with the implementation of Banach fixed point theorem. Our following result will demonstrate the singularity of the solutions.

Theorem 4

The suggested heroin fractional model (8) has a unique solution, if

(35)
1 α M ( α ) + P 0 α M ( α ) Γ ( α ) ξ i < 1 , i = 1 , 2 , 3 .

Proof

Assume that S ( t ) ˆ , X ( t ) ˆ , and Y ( t ) ˆ are other solutions to (8). Then

(36)
S ( t ) S ( t ) ˆ = 1 α M ( α ) ( Ω 1 ( t , S ) Ω 1 ( t , S ˆ ) ) + α M ( α ) Γ ( α ) 0 t ( Ω 1 ( ζ , S ) Ω 1 ( ζ , S ˆ ) ) ( t ζ ) α 1 d ζ . Taking norm on both sides, then we have
(37)
S ( t ) S ( t ) ˆ 1 α M ( α ) ξ 1 S S ˆ + t α M ( α ) Γ ( α ) ξ 1 S S ˆ .
Since 1 1 α M ( α ) P 0 α M ( α ) Γ ( α ) > 0 , we determine that S ( t ) S ( t ) ˆ = 0 . Hence, we find that S ( t ) = S ( t ) ˆ . Similarly, we can demonstrate that X ( t ) = X ( t ) ˆ and Y ( t ) = Y ( t ) ˆ . With this, we conclude the proof. □

4

4 Qualitative perspectives of the heroin model

4.1

4.1 Equilibrium points and basic reproduction number

The proposed fractional heroin epidemic model (8) has two equilibrium points in the feasible region of the model (9). To determine the equilibrium states of the suggested system, we put 0 A B C D t α S ( t ) = 0 , 0 A B C D t α X ( t ) = 0 , 0 A B C D t α Y ( t ) = 0 , which gives that

(38)
η σ 1 S X δ 1 S = 0 , σ 1 S X ρ X + σ 2 X Y δ 1 X = 0 , ρ X σ 2 X Y μ 1 Y = 0 . After simplifying the above system of equations for the state variables, the suggested fractional model has heroin free and heroin persistence equilibrium points which are stated as follows:
(39)
H 0 = S 0 X 0 Y 0 = η δ 1 0 0
and
(40)
H 1 = S 1 X 1 Y 1 = ρ ( ρ + δ 1 ) σ 2 ( 1 + σ 2 ) ρ σ 1 η δ 1 S 1 σ 1 S 1 σ 2 X 1 δ 1 ρ X 1 .

Basic reproduction number is a principle parameter of infectious disease models and usually represented by R 0 . This is referred to as the estimated count of secondary infections generated by one infected person within a population entirely susceptible to the infection. For the proposed heroin fractional epidemic model, R 0 is developed utilizing the upcoming generation matrix strategy (Ahmad et al., 2021; Riaz et al., 2024; Ain et al., 2024).

The suggested system can be written as

(41)
0 A B C D t α Q = A ( Q ) B ( Q ) , where Q = ( X , Y ) and
(42)
A = σ 1 S X + σ 2 X Y σ 2 X Y , B = ( ρ + δ 1 ) X ρ X δ 1 Y .
Following are the Jacobian matrices of (42), expressed as A ̄ = σ 1 S + σ 2 Y σ 2 X σ 2 Y σ 2 X , B ̄ = ( ρ + δ 1 ) 0 ρ δ 1 . The spectral radius of A ̄ B ̄ 1 ( H 0 ) is R 0 . Thus, we have
(43)
R 0 = Φ ( A ̄ B ̄ 1 ) = σ 1 η δ 1 ( ρ + δ 1 ) .

4.2

4.2 Stability analysis

This section discusses the local and global behavior of the considered model (8) at the heroin-free and heroin persistence equilibrium states that by imposing constraints on the R 0 .

4.2.1

4.2.1 Stability of heroin free equilibrium point H 0

Theorem 5

Theorem 5 Butt et al., 2023a

Heroin free equilibrium point H 0 of system (8) is locally asymptotically stable if R 0 < 1 and unstable when R 0 > 1 .

Proof

The Jacobian matrix J 0 ( S 0 , X 0 , Y 0 ) for proposed fractional model (8) at heroin free state H 0 is computed as:

(44)
J 0 ( H 0 ) = δ 1 σ 1 η δ 1 0 0 σ 1 η δ 1 ( ρ + δ 1 ) 0 0 ρ δ 1 . The Jacobian matrix J 0 ( H 0 ) has the following eigenvalues:
(45)
λ 1 = δ 1 ,
(46)
λ 2 = δ 1 ,
(47)
λ 3 = σ 1 η δ 1 ( ρ + δ 1 ) .
The first two eigenvalues λ 1 and λ 2 are clearly less than zero, as δ 1 is positive. Now, consider Eq. (47), such as:
(48)
λ 3 = ( ρ + δ 1 ) [ σ 1 η δ 1 ( ρ + δ 1 ) 1 ] = ( ρ + δ 1 ) [ R 0 1 ] .
Thus, λ 3 is negative if and only if R 0 is less than one. Therefore, all eigenvalues of J 0 ( H 0 ) are negative, and consequently, H 0 is locally asymptotically stable which completes the proof. □

Theorem 6

Theorem 6 Butt et al., 2023a

Heroin free equilibrium point H 0 of system (8) is globally asymptotically stable (GAS) if R 0 < 1 and unstable when R 0 > 1 .

Proof

To examine the global behavior of heroin free equilibrium point, we choose a Lyapunov function as:

(49)
F 0 = S S 0 S 0 ln S S 0 + X + Y . Taking derivative of F 0 with respect to time as:
(50)
0 A B C D t α F 0 = 1 S 0 S 0 A B C D t α S + 0 A B C D t α X + 0 A B C D t α Y .
Employing model Eqs. (8), we yields
(51)
0 A B C D t α F 0 = 1 S 0 S [ η σ 1 S X δ 1 S ] + σ 1 S X + σ 2 X Y δ 1 X σ 2 X Y δ 1 Y .
Above equality, after some calculations takes the form:
(52)
0 A B C D t α F 0 = δ 1 S ( S S 0 ) 2 + σ 1 S 0 X δ 1 ( X + Y ) ,
where σ 1 S 0 X 0 . Therefore
(53)
0 A B C D t α F 0 δ 1 S ( S S 0 ) 2 δ 1 ( X + Y ) .
Using the fact that the result of first term is positive when R 0 1 , then
(54)
0 A B C D t α F 0 0 ,
for all ( S , X , Y ) R + 3 . It is observed that 0 A B C D t α F 0 0 if and only if S = S 0 , X = X 0 = 0 , and Y = Y 0 = 0 . By setting S = S 0 , X = Y = 0 in (8), we have S η δ 1 when t . Applying LaSalles invariance principle (Lassalle, 1976), every solution of the proposed model having initial conditions approaches to H 0 , as t . Thus H 0 is globally asymptotically stable in the region Ξ . As a consequence, the growth of epidemic disappear from the population. □

4.2.2

4.2.2 Stability of heroin persistence equilibrium point H 1

Theorem 7

Theorem 7 Butt et al., 2023a

Heroin persistence equilibrium point H 1 of system (8) is locally asymptotically stable if R 0 > 1 and unstable when R 0 < 1 .

Proof

The Jacobian matrix J 1 ( S 1 , X 1 , Y 1 ) for the considered fractional model (8) at heroin persistence H 1 is computed as:

(55)
J 1 ( H 1 ) = σ 1 X 1 δ 1 σ 1 S 1 0 σ 1 X 1 σ 1 S 1 + σ 2 Y 1 ( ρ + δ 1 ) σ 2 X 1 0 ρ σ 2 Y 1 σ 2 X 1 δ 1 . The eigenvalues of the Jacobian matrix J 1 ( H 1 ) is obtained by solving the following equation:
(56)
| J 1 ( H 1 ) λ I | = 0 ,
where I is the identity matrix of 3-by-3. Therefore, we obtain the characteristic equation as:
(57)
( σ 1 X 1 δ 1 λ ) [ ( σ 1 S 1 + σ 2 Y 1 ( ρ + δ 1 ) λ ) ( σ 2 X 1 δ 1 λ ) σ 2 X 1 ρ + σ 2 2 X 1 Y 1 ] + σ 1 S 1 X 1 ( σ 2 X 1 δ 1 λ ) = 0 .
One eigenvalue of J 1 ( H 1 ) is given by
(58)
λ 1 = ( σ 1 X 1 + δ 1 ) < 0
and, the remaining two eigenvalues correspond to the solutions of the following equation:
(59)
λ 2 + ( e + a b ) + ( c d a b + a e ) = 0 ,
where a = ( σ 2 X 1 + δ 1 ) > 0 , b = ( σ 1 S 1 + σ 2 Y 1 + σ 1 2 S 1 Y 1 ) > 0 , c = σ 2 X 1 ρ > 0 , d = σ 2 2 X 1 Y 1 > 0 , e = ρ + δ 1 > 0 . Therefore, according to the Routh–Hurwitz criteria for second degree polynomials, the following conditions are validated, if R 0 > 1 . e + a b > 0 , c d a b + a e > 0 . Thus, heroin persistence equilibrium point is locally asymptotically stable. □

Theorem 8

Theorem 8 Butt et al., 2023a

Heroin persistence equilibrium point H 1 of system (8) is globally asymptotically stable if R 0 > 1 and unstable when R 0 < 1 .

Proof

Let us consider a Lyapunov function candidate as:

(60)
F 1 = S S 1 S 1 ln S S 1 + X X 1 X 1 ln X X 1 + Y Y 1 Y 1 ln Y Y 1 . Differentiating with respect to time yields:
(61)
0 A B C D t α F 1 = 1 S 1 S 0 A B C D t α S + 1 X 1 X 0 A B C D t α X + 1 Y 1 Y 0 A B C D t α Y .
Utilizing the equations of model, then
(62)
0 A B C D t α F 1 = 1 S 1 S [ η ( σ 1 X δ 1 ) S 1 ( σ 1 X + δ 1 ) ( S S 1 ) ] + 1 X 1 X [ σ 1 S X + σ 2 X Y ( ρ + δ 1 ) X 1 ( ρ + δ 1 ) ( X X 1 ) ] + 1 Y 1 Y [ ρ X ( σ 2 X + δ 1 ) Y 1 ( σ 2 X + δ 1 ) ( Y Y 1 ) ] .
After some simplification, we reached at:
(63)
0 A B C D t α F 1 = Γ 1 Γ 2 ,
where Γ 1 = η + ( σ 1 X + δ 1 ) ( S 1 ) 2 S + ( σ 1 S + σ 2 Y ) X + ( ρ + δ 1 ) ( X 1 ) 2 X + ρ X + ( σ 2 X + δ 1 ) ( Y 1 ) 2 Y and Γ 2 = ( σ 1 X + δ 1 ) ( S S 1 ) 2 S + ( σ 1 X + δ 1 ) S 1 + η S 1 S + ( ρ + δ 1 ) ( X X 1 ) 2 X + ( ρ + δ 1 ) X 1 + ( σ 1 S + σ 2 Y ) X 1 X + ρ X + ( σ 2 X + δ 1 ) ( Y Y 1 ) 2 Y + ( σ 2 X + δ 1 ) Y 1 + ρ X Y 1 Y . Since all the models parameters are non-negative, thus we have 0 A B C D t α F 1 0 , for Γ 1 Γ 2 , and the equality 0 A B C D t α F 1 = 0 , satisfies if and only if S = S 1 , X = X 1 , Y = Y 1 . The heroin persistence equilibrium point H 1 exist, when the basic reproduction number R 0 > 1 , the resulting scenario establishes that the sole invariant set H 1 attains the status of being the most significant invariant set encompassed within
(64)
Ξ : = { ( S , X , Y ) R + 3 : 0 A B C D t α F 1 = 0 } .
Thus, in accordance with LaSalle’s invariance principle (Lassalle, 1976), it is affirmed that H 1 achieves global asymptotic stability within the set Ξ . This indicates that individuals afflicted by heroin addiction will disseminate in the population, ultimately leading to a pandemic. □

5

5 Numerical analysis

This section is devoted for the numerical results and simulations of heroin epidemic model using a well known approach constructed by Toufik and Atangana (Mekkaoui and Atangana, 2017). We now implement this scheme to each equation of the system (8) to get the numerical solutions of the proposed system, we obtain the following results:

(65)
S ( t ) S ( 0 ) = 1 α M ( α ) Ω 1 t , S ( t ) + α M ( α ) Γ ( α ) 0 t Ω 1 ( ζ , S ( ζ ) ) ( t ζ ) α 1 d ζ . At a given point t = t j + 1 , j = 0 , 1 , 2 , then the Eq. (65) becomes S ( t j + 1 ) S ( 0 ) = 1 α M ( α ) Ω 1 t j , S ( t j ) + α M ( α ) Γ ( α ) t n t n + 1 Ω 1 ( ζ , S ( ζ ) ) ( t j + 1 ζ ) α 1 d ζ . The above equation can be written as:
(66)
S ( t j + 1 ) S ( 0 ) = 1 α M ( α ) Ω 1 t j , S ( t j ) + α M ( α ) Γ ( α ) n = 0 j t n t n + 1 Ω 1 ( ζ , S ( ζ ) ) ( t j + 1 ζ ) α 1 d ζ .
Applying two-step Lagrange polynomial interpolation on the function Ω 1 ( ζ , S ( ζ ) ) in the interval [ t j , t j + 1 ] . Therefore, we obtain
(67)
S j + 1 = S 0 + 1 α M ( α ) Ω 1 t j , S j + α M ( α ) Γ ( α ) n = 0 j Ω 1 ( t n , S ( t n ) ) h t n t n + 1 ( ζ t n 1 ) ( t j + 1 ζ ) α 1 d ζ Ω 1 ( t n 1 , S ( t n 1 ) ) h t n t n + 1 ( ζ t n ) ( t j + 1 ζ ) α 1 d ζ ,
where Π n 1 = t n t n + 1 ( ζ t n 1 ) ( t j + 1 ζ ) α 1 d ζ , = h α + 1 α ( α + 1 ) [ j + 1 n α j n + 2 + α j n α j n + 2 + 2 α ] and Π n = t n t n + 1 ( ζ t n ) ( t j + 1 ζ ) α 1 d ζ , = h α + 1 α ( α + 1 ) [ j + 1 n α + 1 j n α j n + 1 + α ] . The above integrals substituted into the Eq. (67), then we have the solution S ( t ) as follows:
(68)
S j + 1 = S 0 + 1 α M ( α ) Ω 1 t j , S j + α M ( α ) Γ ( α ) n = 0 j [ h α Ω 1 ( t n , S n ) Γ ( α + 2 ) ( j + 1 n α j n + 2 + α j n α j n + 2 + 2 α ) h α Ω 1 ( t n 1 , S n 1 ) Γ ( α + 2 ) j + 1 n α + 1 j n α j n + 1 + α ] .
Similarly, for the other state variables, we evaluate the following schemes:
(69)
X j + 1 = X 0 + 1 α M ( α ) Ω 2 t j , X j + α M ( α ) Γ ( α ) n = 0 j [ h α Ω 2 ( t n , X n ) Γ ( α + 2 ) ( j + 1 n α j n + 2 + α j n α j n + 2 + 2 α ) h α Ω 2 ( t n 1 , X n 1 ) Γ ( α + 2 ) j + 1 n α + 1 j n α j n + 1 + α ] ,
(70)
Y j + 1 = Y 0 + 1 α M ( α ) Ω 3 t j , Y j + α M ( α ) Γ ( α ) n = 0 j [ h α Ω 3 ( t n , Y n ) Γ ( α + 2 ) ( j + 1 n α j n + 2 + α j n α j n + 2 + 2 α ) h α Ω 3 ( t n 1 , Y n 1 ) Γ ( α + 2 ) j + 1 n α + 1 j n α j n + 1 + α ] .
We conduct mathematical simulations across various fractional values of the parameter α illustrating the dynamic spread of the heroin model within the human population over time. Our observations reveal a seamless match between the results from our suggested method and the inherent dynamics of the heroin model, confirming the reliability of the approach. Additionally, we explore the impact of the infection probability rate and the chance of non-drug addicts transitioning to drug addiction among humans. Finally, we present an analysis of the reproduction number’s sensitivity to the involved factors.

Dynamics of the populations of heroin model at heroin persistence state at different values of α .
Fig. 2
Dynamics of the populations of heroin model at heroin persistence state at different values of α .
PRCC statistics regarding the significance of factors associated with R 0 .
Fig. 3
PRCC statistics regarding the significance of factors associated with R 0 .

5.1

5.1 Sensitivity analysis

Making decisions about effectively managing a disease necessitates careful consideration of the sensitivity analysis concept. The sensitivity analysis enables us to examine how variables fluctuate when the parameters in R 0 are altered. It highlights the model’s most sensitive and impactful parameters and their effects on R 0 (Butt et al., 2023a,b).

Definition 3

(Butt et al., 2023b) The normalized forward sensitivity index ( Υ ς ) of the basic reproduction number R 0 that depends on a parameter ς is provided below as:

(71)
Υ ς = ς R 0 R 0 ς .

To examine the sensitivity of R 0 , we compute its derivatives as follows: R 0 δ 1 = σ 1 η ( ρ + 2 δ 1 ) [ δ 1 ( ρ + δ 1 ) ] 2 , R 0 η = σ 1 δ 1 ( ρ + δ 1 ) , R 0 ρ = σ 1 η δ 1 ( ρ + δ 1 ) 2 , R 0 σ 1 = η δ 1 ( ρ + δ 1 ) . The normalized sensitivity indices of the involved parameters are obtained as: Υ δ 1 = δ 1 R 0 R 0 δ 1 = 2 δ 1 ρ + δ 1 , Υ η = η R 0 R 0 η = 1 , Υ ρ = ρ R 0 R 0 ρ = ρ ρ + δ 1 , Υ σ 1 = σ 1 R 0 R 0 σ 1 = 1 .

Table 2 and Fig. 3 depict the favorable influence of η and σ 1 on the threshold parameter R 0 . This indicates that higher values of these parameters result in an increase in R 0 . The sensitivity indices obtained clearly show 10% rise in the new drug users, and infection rate occurs to raise the value of R 0 by 10%, 10%, correspondingly, and can ultimately leading to a disease outbreak. Across the other perspective, R 0 is negatively impacted by δ 1 and ρ as shown by Fig. 4 (see Table 3).

Table 2 Sensitivity indices of R 0 against the parameters.
Sensitivity index Value
Υ δ 1 −1.3333
Υ η + 1 . 0000
Υ ρ −0.3333
Υ σ 1 + 1 . 0000
The behaviors of R 0 under different biological parameters.
Fig. 4
The behaviors of R 0 under different biological parameters.
Table 3 Initial conditions and parameter values of the heroin model.
Initial condition Value Source
S ( 0 ) 0.50 Raza et al. (2022a)
X ( 0 ) 0.30 Raza et al. (2022a)
Y ( 0 ) 0.20 Raza et al. (2022a)
Parameter Value Source
η 0.04 Raza et al. (2022a)
σ 1 0.02 for H 0 (0.20 for H 1 ) Raza et al. (2022a)
σ 2 0.03 Raza et al. (2022a)
δ 1 0.04 Raza et al. (2022a)
ρ 0.02 Raza et al. (2022a)

5.2

5.2 Numerical results and discussion

This subsection examines the effect of the fractional parameter α on the heroin epidemic disease in the model. A graphical outcome of the model application (8) that depicts the impacts of fractional order on the population in every section. We have illustrated the heroin dynamics of a model for α = 0 . 80 , 0 . 84 , 0 . 88 , 0 . 92 , 0 . 96 and 1 , respectively. For the heroin-free equilibrium, an elevation in the value of α from 0.80 corresponds to boost in the count of susceptible cases. However, the populations of drug users and non-drug users decline undeviatingly by increasing the value of α from 0.8. In contrast, the susceptible decreases with the rise of fractional parameters in heroin persistence point and the behavior of drug users and non-drug users is observed oppositely, as illustrated by Fig. 2. At the end, when contemplating the case where α = 1 , we arrived at the results of Raza et al. (2022a), where a different numerical approach was used. The behavior of the individuals is illustrated in Fig. 5 at different values of σ 1 by fixing the value of the fractional parameter. These results show a distinct trend, the proportion of susceptible people declines as the infection probability rate rises. The reason for this is that a greater number of people who were previously susceptible to infection are now either recovering or current drug users due to the increased likelihood of infection. On the other hand, as the probability of infection increases, so does the population of drug users. As there is a direct association between the probability of infection and the number of new drug users, this increase happens because a greater infection rate makes more susceptible people more likely to start using drugs. On the other side, the findings indicate a rise in non-drug users, which may be people who have abstained from drug use in the past but have recovered or those who have rejected using drugs despite the risk of infection. The fact that this group is expanding suggests that, despite greater dangers, some members of the population at risk can recover or avoid infection.

The behavior of the individuals is illustrated in Fig. 6 at different values of σ 2 by fixing the value of the fractional parameter. According to the analysis, fewer people are susceptible to heroin use when the chance rate, that is, the rate at which people are exposed to or attracted by the drug—increases. This occurs because there is a greater probability that those at risk will become drug users when the chance rate is higher. With a rise in the chance rate, there is a consistent rise in the population of drug users. This pattern shows that increasing exposure or temptation rates leads to more people shifting from being vulnerable to actually using heroin. As the chance rate rises, the proportion of non-drug users falls. This group comprises people who have either never used heroin or who have overcome their addiction. This group may decline as a result of a higher chance rate, making it more difficult for people to prevent or rehabilitate from habit.

Effect of σ 1 on the behavior of heroin model populations for α = 0 . 92 , 0 . 96 and 1 .
Fig. 5
Effect of σ 1 on the behavior of heroin model populations for α = 0 . 92 , 0 . 96 and 1 .
Effect of σ 2 on the behavior of heroin model populations for α = 0 . 96 , 0 . 98 and 1.00.
Fig. 6
Effect of σ 2 on the behavior of heroin model populations for α = 0 . 96 , 0 . 98 and 1.00.

After examining, it can be stated that the mechanisms of heroin addiction are shaped by both the infection probability rate and the chance rate. A multimodal strategy that incorporates lowering risk factors, improving recovery guides, and putting in place powerful educational activities is necessary for successful management and safeguarding against heroin addiction. Public health measures can more successfully stop the rapid growth of heroin addiction and assist those who are impacted in a successful recovery by considering these factors.

6

6 Conclusion

In this manuscript, we have developed a novel fractional epidemic model of heroin population consequences by applying the Atangana–Baleanu fractional derivative. The fundamental reproduction number has been determined to investigate the dynamical behavior of the fractional model. The effects of heroin can be controlled by using a reproduction number. The primary properties such as positiveness, boundedness, existence, and uniqueness of the solution have been validated. It has been concluded that the equilibrium states are locally and globally asymptotically stable. In the end, we have applied a computational approach that assesses the dynamics of the proposed system. The results of the numerical technique show that the susceptible people increase with the increase in fractional order, and the other two subpopulations show opposite behavior.

CRediT authorship contribution statement

Miguel Vivas-Cortez: Data curation, Funding acquisition, Methodology. Abu Bakar: Formal analysis, Software, Validation. M.S. Alqarni: Conceptualization, Resource Management, Methodology. Nauman Raza: Investigation,Supervision, Validation. Talat Nazir: Supervision, Project administration, Data curation. Muhammad Farman: Methodology, Validation.

Acknowledgments

The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University, Saudi Arabia for funding this work through Large Research Project under grant number RGP2/113/45. The authors thank the anonymous referees for their valuable constructive comments and suggestions, which improved the quality of this paper in the present form.

Declaration of competing interest

The authors declare that there is no conflict of interest regarding the submission of this paper.

References

  1. , , . Integration by parts and its applications of a new nonlocal fractional derivative with Mittag-Leffler nonsingular kernel. J. Nonlinear Sci.. 2017;10(3):1098-1107.
    [CrossRef] [Google Scholar]
  2. , , , . On common fixed points of non-Lipschitzian semigroups in a hyperbolic metric space endowed with a graph. J. Anal. 2023:2023.
    [CrossRef] [Google Scholar]
  3. , . Mathematical analysis for the effect of voluntary vaccination on the propagation of Corona virus pandemic. Results Phys.. 2021;31:104917
    [CrossRef] [Google Scholar]
  4. , . A stochastic approach for co-evolution process of virus and human immune system. Sci. Rep.. 2024;14(1):10337.
    [Google Scholar]
  5. , . Sur les opérations dans les ensembles abstraits et leur application aux équations intégrales. Fund. Math.. 1922;3:133-181.
    [Google Scholar]
  6. , , . Real world applications of fractional models by Atangana - Aleanu fractional derivative. Chaos Solit. Fractals. 2018;116(2018):121-125.
    [CrossRef] [Google Scholar]
  7. , . Numerical analysis of Atangana-Baleanu fractional model to understand the propagation of a novel Corona virus pandemic. Alex. Eng. J.. 2022;61(9):7007-7027.
    [CrossRef] [Google Scholar]
  8. , . Dynamical analysis of a nonlinear fractional cervical cancer epidemic model with the nonstandard finite difference method. Ain Shams Eng. J. 2023102479
    [CrossRef] [Google Scholar]
  9. , . Investigating the fractional dynamics and sensitivity of an epidemic model with nonlinear convex rate. Results Phys.. 2023;54:107089
    [CrossRef] [Google Scholar]
  10. , . Optimally analyzed fractional coronavirus model with Atangana-Baleanu derivative. Results Phys.. 2023;53:106929
    [CrossRef] [Google Scholar]
  11. , , . A mathematical model of reward-mediated learning in drug addiction. Chaos. 2022;32(2):021102
    [CrossRef] [Google Scholar]
  12. , . United nations office on drugs and crime: World drug report 2017. SIRIUS-Zeitschrift für Strategische Analysen. 2018;2(1):85-86.
    [Google Scholar]
  13. , , , . Heroin and fentanyl in dallas county: A 5-year retrospective review of toxicological, seized drug, and demographical data. J. Forensic Sci.. 2023;68(1):222-232.
    [CrossRef] [Google Scholar]
  14. , . Atangana-Baleanu fractional derivative applied to electromagnetic waves in dielectric media. J. Electromagn. Waves Appl.. 2016;30(15):1937-1952.
    [CrossRef] [Google Scholar]
  15. , , . A note on global stability for a heroin epidemic model with distributed delay. Appl. Math. Lett.. 2013;26(7):687-691.
    [CrossRef] [Google Scholar]
  16. , . Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age. AIMS Math.. 2021;6(2):1377-1394.
    [CrossRef] [Google Scholar]
  17. , . The Stability of Dynamical Systems, Regional Conference Series in Applied Mathematics. Philadelphia, Pa, USA: SIAM; .
    [CrossRef]
  18. , , . Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate. J. Nonlinear Sci. Appl.. 2016;9(5):2149-2160.
    [CrossRef] [Google Scholar]
  19. , , . New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models. Eur. Phys. J. Plus. 2017;132(444):1-16.
    [CrossRef] [Google Scholar]
  20. , , . On the formulation of Adams–Bashforth scheme with Atangana-Baleanu-Caputo fractional derivative to model chaotic problems. Chaos. 2019;29:023111
    [CrossRef] [Google Scholar]
  21. , . Applying fixed point methods and fractional operators in the modelling of novel Coronavirus 2019-nCoV/SARS-CoV-2. Results Phys.. 2020;19:103433
    [CrossRef] [Google Scholar]
  22. , . New numerical dynamics of the heroin epidemic model using a fractional derivative with Mittag-Leffler kernel and consequences for control mechanisms. Results Phys.. 2022;35:105304
    [Google Scholar]
  23. , . Dynamical and nonstandard computational analysis of heroin epidemic model. Results Phys.. 2022;34:105245
    [CrossRef] [Google Scholar]
  24. , . Numerical simulations of the fractional-order SIQ mathematical model of Corona virus disease using the nonstandard finite difference scheme. Malays. J. Math. Sci.. 2022;16(3)
    [CrossRef] [Google Scholar]
  25. , . A numerical efficient splitting method for the solution of HIV time periodic reaction–diffusion model having spatial heterogeneity. Phys. A. 2023;609:128385
    [CrossRef] [Google Scholar]
  26. , . Fractional dynamics and sensitivity analysis of measles epidemic model through vaccination. Arab J. Basic Appl. Sci.. 2024;31(1):265-281.
    [Google Scholar]
  27. , , . Dynamics of an epidemic model with impact of toxins. Phys. A. 2019;527:121152
    [CrossRef] [Google Scholar]
  28. , , . On the approximation of fixed points for the class of mappings satisfying (CSC)-condition in Hadamard spaces. Carpathian Math. Publ.. 2023;15(2):495-506.
    [CrossRef] [Google Scholar]
  29. , , , . Some fixed-point results for the K F -iteration process in hyperbolic metric spaces. Symmetry. 2023;15(7):1360.
    [CrossRef] [Google Scholar]
  30. , , . Fractional differential equations with Atangana-Baleanu fractional derivative: Analysis and applications. Chaos Solit. Fractals X. 2019;2:100013
    [CrossRef] [Google Scholar]
  31. , , , , . Chronic heroin use disorder and the brain: Current evidence and future implications. Progr. Neuro-Psychopharmacol. Biol. Psychiat.. 2021;111:110148
    [CrossRef] [Google Scholar]
  32. , , . Analysis of a heroin epidemic model with saturated treatment function. J. Appl. Math. (2017)
    [CrossRef] [Google Scholar]
  33. , . A neural study of the fractional heroin epidemic model. Comput. Mater. Contin.. 2023;74(2):2.
    [CrossRef] [Google Scholar]
  34. , . Dynamical analysis of a heroin - cocaine epidemic model with nonlinear incidence and spatial heterogeneity. J. Biol. Dyn.. 2023;17(1):2189026
    [CrossRef] [Google Scholar]
  35. , . Fractional order heroin epidemic dynamics. Alex. Eng. J.. 2021;60(6):5157-5165.
    [CrossRef] [Google Scholar]
  36. , , . Stability analysis of a reaction–diffusion heroin epidemic model. Complexity 2020:1-16.
    [CrossRef] [Google Scholar]
Show Sections