
Heroin addiction is a continuously progressing phenomenon that represents a major problem for world's the public health; this indicates that the development of new methodologies to address the issue at the international level is a crucial priority. To study its transmission dynamics, a new stochastic fractional delayed heroin model based on stochastic fractional delay differential equations (SFDDEs) was developed to focus on the positive aspects of randomness and memory effects. The positive, boundedness, existence, and uniqueness of the model were studied rigorously. The equilibria (i.e., heroin-free equilibrium and the present equilibrium, which gives a clue about both eradication and persistence cases), reproduction number, and sensitivity of parameters were analyzed. The local and global stability of the new model was studied around its steady states. Also, well-known theorems are presented to investigate the extinction and persistence of heroin. The Grunwald-Letnikove non-standard finite difference (GL-NSFD) method was used for the efficient computational analysis of the stochastic fractional delayed model. For the dynamical consistency of the model, the positivity and boundedness of an efficient method were studied rigorously. The given study focuses on delay strategies and fractional calculus that could be useful in formulating specific measures for regulating addiction. Moreover, the simulated results support the theoretical analysis of the model and validate it.
Citation: Feliz Minhós, Ali Raza, Umar Shafique. An efficient computational analysis for stochastic fractional heroin model with artificial decay term[J]. AIMS Mathematics, 2025, 10(3): 6102-6127. doi: 10.3934/math.2025278
[1] | Wedad Albalawi, Muhammad Imran Liaqat, Fahim Ud Din, Kottakkaran Sooppy Nisar, Abdel-Haleem Abdel-Aty . Well-posedness and Ulam-Hyers stability results of solutions to pantograph fractional stochastic differential equations in the sense of conformable derivatives. AIMS Mathematics, 2024, 9(5): 12375-12398. doi: 10.3934/math.2024605 |
[2] | Zhengwen Yin, Yuanshun Tan . Threshold dynamics of stochastic SIRSW infectious disease model with multiparameter perturbation. AIMS Mathematics, 2024, 9(12): 33467-33492. doi: 10.3934/math.20241597 |
[3] | Aisha F. Fareed, Emad A. Mohamed, Mokhtar Aly, Mourad S. Semary . A novel numerical method for stochastic conformable fractional differential systems. AIMS Mathematics, 2025, 10(3): 7509-7525. doi: 10.3934/math.2025345 |
[4] | Omar Kahouli, Saleh Albadran, Zied Elleuch, Yassine Bouteraa, Abdellatif Ben Makhlouf . Stability results for neutral fractional stochastic differential equations. AIMS Mathematics, 2024, 9(2): 3253-3263. doi: 10.3934/math.2024158 |
[5] | Wedad Albalawi, Muhammad Imran Liaqat, Kottakkaran Sooppy Nisar, Abdel-Haleem Abdel-Aty . Qualitative study of Caputo Erdélyi-Kober stochastic fractional delay differential equations. AIMS Mathematics, 2025, 10(4): 8277-8305. doi: 10.3934/math.2025381 |
[6] | Muhammad Imran Liaqat, Fahim Ud Din, Wedad Albalawi, Kottakkaran Sooppy Nisar, Abdel-Haleem Abdel-Aty . Analysis of stochastic delay differential equations in the framework of conformable fractional derivatives. AIMS Mathematics, 2024, 9(5): 11194-11211. doi: 10.3934/math.2024549 |
[7] | Wedad Albalawi, Muhammad Imran Liaqat, Fahim Ud Din, Kottakkaran Sooppy Nisar, Abdel-Haleem Abdel-Aty . Significant results in the pth moment for Hilfer fractional stochastic delay differential equations. AIMS Mathematics, 2025, 10(4): 9852-9881. doi: 10.3934/math.2025451 |
[8] | Zafar Iqbal, Muhammad Aziz-ur Rehman, Muhammad Imran, Nauman Ahmed, Umbreen Fatima, Ali Akgül, Muhammad Rafiq, Ali Raza, Ali Asrorovich Djuraev, Fahd Jarad . A finite difference scheme to solve a fractional order epidemic model of computer virus. AIMS Mathematics, 2023, 8(1): 2337-2359. doi: 10.3934/math.2023121 |
[9] | Yuanfu Shao . Dynamics and optimal harvesting of a stochastic predator-prey system with regime switching, S-type distributed time delays and Lévy jumps. AIMS Mathematics, 2022, 7(3): 4068-4093. doi: 10.3934/math.2022225 |
[10] | Tahir Khan, Fathalla A. Rihan, Muhammad Bilal Riaz, Mohamed Altanji, Abdullah A. Zaagan, Hijaz Ahmad . Stochastic epidemic model for the dynamics of novel coronavirus transmission. AIMS Mathematics, 2024, 9(5): 12433-12457. doi: 10.3934/math.2024608 |
Heroin addiction is a continuously progressing phenomenon that represents a major problem for world's the public health; this indicates that the development of new methodologies to address the issue at the international level is a crucial priority. To study its transmission dynamics, a new stochastic fractional delayed heroin model based on stochastic fractional delay differential equations (SFDDEs) was developed to focus on the positive aspects of randomness and memory effects. The positive, boundedness, existence, and uniqueness of the model were studied rigorously. The equilibria (i.e., heroin-free equilibrium and the present equilibrium, which gives a clue about both eradication and persistence cases), reproduction number, and sensitivity of parameters were analyzed. The local and global stability of the new model was studied around its steady states. Also, well-known theorems are presented to investigate the extinction and persistence of heroin. The Grunwald-Letnikove non-standard finite difference (GL-NSFD) method was used for the efficient computational analysis of the stochastic fractional delayed model. For the dynamical consistency of the model, the positivity and boundedness of an efficient method were studied rigorously. The given study focuses on delay strategies and fractional calculus that could be useful in formulating specific measures for regulating addiction. Moreover, the simulated results support the theoretical analysis of the model and validate it.
Humans have been aware of heroin's stimulating qualities for a very long time. The primary sources of heroin are smack and chakra. Even in current times, there is proof that these poppies were cultivated. Aeon's ago, poppies were harvested for their opium-rich sap, which was widely used. In particular, opium poppies were widespread throughout the eighteenth century. In 1897, a German pharmaceutical company announced that heroin was now a medicinal product and started marketing it to treat various illnesses like morphine addiction and tuberculosis [1]. In [2], authors provided numerical solutions for the traditional White and Comiskey model (NFD-WCM) based on regressive fractional order derivatives. The NDF-WCM has genuine and precise solutions of NFD-WCM. Levenberg-Marquard backpropagation (LMB), a controlled probabilistic determining technique using neural networks (NNs), provides the solutions to the fractional NFD-WCM. In [3], the authors focused on an instance of the heroin pandemic with limited balanced network technology. In [4], authors examined a model of the transitioning of heroin into a global epidemic, wherein the opium plants' blooming and fruiting seasons influence heroin availability. In [5], authors focused on the dynamic characteristics of a stochastic model for a worldwide heroin epidemic with Leonardy sounds. In [6], authors presented and examined a probabilistic model with a propagated delay to describe heroin use in an inconsistent setting. In [7], authors emphasized the nonlinear age-space structured heroin propagation model's computational boundary reliability. The primary nonlinear age-space structured model was discretized spatially to generate a semi-discrete structure. In [8], authors committed to researching the dynamics of a heroin-cocaine pandemic model constructed according to a certain age. In [9], authors created a unique R-D heroin epidemic model that included an enclosure for everlasting vaccination and a segment for readmission. In [10], authors explained the dynamics of a deterministic model of illegal opioid use disorder (IOUD). In [11], authors discussed White and Comiskey's standard mathematical model, which illustrates heroin epidemics and provides mathematical and numerical evidence for the existence of stable equilibrium. In [12], authors utilized an Atangana–Blaenau fractional-order derivative in the Caputo practical sense to create a heroin epidemic model that closely depicted real-world issues and was endowed with both permanent immunity and recovery. In [13], authors proved that both susceptibility and recovery are age-dependent, and the current age-structured heroin epidemic model was developed. In [14], authors examined internal competition among drug users for a small amount of drugs in a heroin epidemic model with treatment age. The impact of a drug user's treatment duration before stopping treatment was examined, and the model under consideration was capable of backward bifurcation, indicating the potential for two endemic equilibriums. In [15], authors investigated an epidemic model of heroin use that takes into account only the age distribution of current users. In [16], authors evaluated the heroin epidemic model's overall characteristics, including its nonlinear distribution function and time-distributed latency. In [17], authors created a model of the HIV/heroin co-infection pandemic to explain how the two diseases co-transmit within the same community. In [18], authors showed that the dissemination of knowledge regarding the consequences of heroin is a necessary condition for an individual's behavioral reaction. Also, a mathematical model of the heroin pandemic was analyzed, incorporating preventative information and treatment as control interventions. In [19], authors examined in detail a novel spatiotemporal model that incorporates the concept of the Laplacian operator into the cocaine-heroin epidemiology model with essential dynamics. In [20], authors employed various demographic data, and an innovative mathematical technique was used to examine the dynamics of the heroin epidemic model and its detrimental effects on society. The authors studied different fractional models like COVID-19, Hepatitis B, HIV/HCV co-infection, and hepatitis C using different techniques [21–25].
Using a stochastic delayed technique offers benefits, in science by combining stochastic processes, fractional calculus, and temporal delays to create more precise models of complex systems. This advanced strategy utilizes derivatives to account for memory effects, which are vital in fields like biology and economics where past events heavily influence present circumstances. By including delays this approach paints a realistic picture of the time lag between cause and effect, thus improving the accuracy of system portrayals. The stochastic aspect deals with uncertainty. Randomness is present in these systems, ensuring that models can better handle unpredictability. Together, these components enhance the accuracy and reliability of models, resulting in accurate representations of real-world phenomena.
The paper is structured into seven successive sections to explore significant areas of the research. Section 1 presents an adequate literature review of infections comparable to heroin. This section provides background information and related findings to set the context for the study. Section 2 examines the complexity of the delayed model, especially the stochastic fractional delayed heroin model. This section also provides mathematical discussions on local and global stabilities' assessment of the model's equilibria. Furthermore, the section includes an analysis of the reproduction number as well as the stability analysis of different levels of the model. Section 3 goes further in establishing the sensitivity of the reproduction number for the system. This sensitivity analysis is conducted using a stochastic fractional delay differential equation, which gives a more profound insight into the alterations of variables across the system. In section 4, the stochastic fractional delayed model is analyzed using the stochastic GL-NSFD method. This section revolves around the utilization of this method in analyzing the model as well as its feasibility. Section 5 focuses on the analysis of the non-negative and boundedness of the stochastic GL-NSFD method. Achieving these properties is important to support the generalizability of the model in capturing real-world scenarios. Section 6 is devoted to the computation and shows how the numerical simulations, derived in the previous sections, can be implemented. It presents a detailed analysis of the findings, which focuses on the implications of the research in the context of practical application. Finally, Section 7 contains the findings of the study, its originality, and future research recommendations.
This section analyzes the dynamics of the heroin epidemic model, as presented in [1]. The model comprises many significant state variables to represent the different groups of the population affected by the disease. In specific terms, U(t) defines the probability density function of the group of persons most vulnerable to heroin usage based on their propensity to get addicted. The variable V(t) depicts the quantitative future state of the affected community in connection with continual drug abuse, namely active heroin users. W(t), representing non-drug users, is used for evaluating the exposure to risk between drug takers and current heroin users. Figure 1 shows the flow in and out of these states over time, creating a series of flows that is useful to give insights into how the epidemic unfolds. The flow map informs a representation that characterizes susceptible persons, drug users, and non-drug users in the model community and facilitates the interpretation of our ability to capture historically reported heroin epidemic situations.
Several parameters that characterize rates of transmission determine the dynamics of the epidemic and are essential for the assessment of the epidemic's progress. The main parameters are (See Table 1):
Parameters | Description |
Λ | The rate at which new people begin using heroin |
μ1 | The rate at which current drug users discontinue their regular use |
β1 | The risk of infection for drug users |
β3 | The pace at which non-drug users enter the drug-user compartment |
ρ | The mortality rate among drug users |
e−μ1τ | The artificial delay parameter |
These parameters are incorporated into the model using a set of linked delay differential equations (DDEs) with first-order nonlinearity, together with the long-term effects of treatment and quitting, giving a complete framework for studying the heroin epidemic model.
dU(t)dt=Λ−β1U(t−τ)V(t−τ)e−μ1τ−μ1U(t),dV(t)dt=β1U(t−τ)V(t−τ)e−μ1τ−ρV(t)+β3V(t)W(t)−μ1V(t),dW(t)dt=ρV(t)−β3V(t)W(t)−μ1W(t).} | (1) |
The non-negative (initial) conditions for the system (1) are U(0)≥0,V(0)≥0,W(0)≥0, t≥0,τ<t and N=U+V+W.
The delayed deterministic model (1) is the basic epidemiological model, in which the first-order temporal derivatives are replaced by Caputo fractional derivatives of order α. In our opinion, this will provide a better representation of heroin viral infection. The model analysis is given in part by Caputo fractional derivatives to account for memory effects. These derivatives also inform how the system relies on previous states and saves it in long-term memory about past infections. With this modification, the model is now able to capture a complex nonlinear response of infection dynamics over time. Stochastic perturbations in the model help to manage this natural unpredictability and variability of biological systems. These stochastic elements are included in the model to simulate and evaluate the complexities of disease transmission more intuitively, accounting for the inherent randomness and unpredictable nature of biological processes within real outbreaks.
0cDαt[U(t)]=Λα−βα1U(t−τ)V(t−τ)e−μα1τ−μα1U(t)+σ1U(t)d(B(t))0cDαt[V(t)]=βα1U(t−τ)V(t−τ)e−μα1τ−ραV(t)+βα3V(t)W(t)−μα1V(t)+σ2V(t)d(B(t))0cDαt[W(t)]=ραV(t)−βα3V(t)W(t)−μα1W(t)+σ3W(t)d(B(t))}. | (2) |
This strategy indicates how the random disturbances in the system can be represented by the stochastic fluctuations σi;i=1,2,3. Let B(t) be Brownian motion, a continuous stochastic process for time t>0. The parameter τ acts as a time delay to the system provided that τ<t; also, the effect of the delayed feedback is visible only after a time interval.
Preliminaries. In the conceptualization of Caputo, the following foundational preliminary definitions are crucial for a thorough understanding of the fractional derivative concept:
Definition 1. For a function q∈Cn, the Caputo fractional derivative of order α∈(m−1,m),m∈N is
0cDαtq(t)=1Γ(m−α)∫t0qm(T)dT(t−T)α+1−m. | (3) |
Definition 2. For the function q(t), the expression describes the equivalent fractional integral with order α>0.
Iαtq(t)=1Γ(α)∫t0(t−T)α−1q(T)dT. | (4) |
Where "Γ" is the gamma function displayed.
The purpose of this section is to verify the existence and uniqueness of solutions for the stochastic fractional delayed model given in (2). We adapt the notion of fractional integral to this specific case when the system is subject to those exact initial conditions to establish these properties. In this benchmark, σi=0;i=1,2,3. Given that there is no interacting deterministic feature, this assumption helps reduce the stochastic noise components to zero and makes us able to focus exclusively on the deterministic features of any system. After that, CL can be combined with applicable mathematical tools to guarantee the existence and uniqueness of fractal delayed model outcomes.
U(t)=U0+1Γ(α)∫t0(t−u)α−1P1(u,U)duV(t)=V0+1Γ(α)∫t0(t−u)α−1P2(u,V)duW(t)=W0+1Γ(α)∫t0(t−u)α−1P3(u,W)du}. | (5) |
In system (5), the functions defined under the integral are as follows:
P1(t,U)=Λα−βα1U(t)V(t)e−μα1τ−μα1U(t)P2(t,V)=βα1U(t)V(t)e−μα1τ−ραV(t)+βα3V(t)W(t)−μα1V(t)P3(t,W)=ραV(t)−βα3V(t)W(t)−μα1W(t)}. | (6) |
There is also a general assumption that the limiting functions U(t),V(t),andW(t) are strictly non-negative. For these functions of asymptotic behavior, positive constants ε1,ε2,andε3 should be set. In such a manner, as the period (t) undergoes an increasing scale, the values of U(t),V(t),andW(t) tend to approach particular positive values below which no greater number than a predetermined one can be found (those will be denoted by ε1,ε2,andε3).
Such that
‖U\left(t\right)‖\le {\mathsf{ε}}_{1},‖V\left(t\right)‖\le {\mathsf{ε}}_{2}\;\mathrm{a}\mathrm{n}\mathrm{d}\;‖W\left(t\right)‖\le {\mathsf{ε}}_{3} . |
Theorem 1. Let {\mathcal{P}}_{i} be functions for i = \mathrm{1, 2}, 3 with 0\le \mathcal{W} = max\{\mathrm{1, 2}, 3\} < 1 , satisfying 0\le \mathcal{W} < 1 , and assuming {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 . Then each function {\mathcal{P}}_{i} satisfies the Lipschitz condition and is a contraction mapping.
Proof. First of all, we consider the function {\mathcal{P}}_{1} as our main concern. To further understand the associations and consequences of the following functions, let's define and analyze them in the context of U and {U}_{1} .
‖{\mathcal{P}}_{1}\left(t,U\right)-{\mathcal{P}}_{2}\left(t,{U}_{1}\right)‖ = ‖{\beta }_{1}^{\alpha }\left(U-{U}_{1}\right)V{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\left(U-{U}_{1}\right)‖ . |
‖{\mathcal{P}}_{1}\left(t,U\right)-{\mathcal{P}}_{2}\left(t,{U}_{1}\right)‖\le ‖{\beta }_{1}^{\alpha }\left(U-{U}_{1}\right)V{e}^{-{\mu }_{1}^{\alpha }\tau }‖+‖{\mu }_{1}^{\alpha }\left(U-{U}_{1}\right)‖ . |
‖{\mathcal{P}}_{1}\left(t,U\right)-{\mathcal{P}}_{2}\left(t,{U}_{1}\right)‖\le \left({\beta }_{1}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }‖V‖+{\mu }_{1}^{\alpha }\right)‖U-{U}_{1}‖ . |
‖{\mathcal{P}}_{1}\left(t,U\right)-{\mathcal{P}}_{2}\left(t,{U}_{1}\right)‖\le \left({\beta }_{1}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }{\mathsf{ε}}_{2}+{\mu }_{1}^{\alpha }\right)‖U-{U}_{1}‖ . |
‖{\mathcal{P}}_{1}\left(t,U\right)-{\mathcal{P}}_{2}\left(t,{U}_{1}\right)‖\le {\xi }_{1}‖U-{U}_{1}‖ . | (7) |
The equation {\xi }_{1} = \left({\beta }_{1}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }{\mathsf{ε}}_{2}+{\mu }_{1}^{\alpha }\right) confirms that function {\xi }_{1} satisfies the Lipschitz condition. To ensure that the Lipschitz criteria apply for the functions {\mathcal{P}}_{i} , where i = \mathrm{2, 3} , we may employ an analogous method that takes advantage of the nature of these functions. Additionally, the fact that these functions are contractions is confirmed under the specific conditions \mathcal{W} = max\{\mathrm{1, 2}, 3\} < 1 and \mathcal{W} = max\{\mathrm{1, 2}, 3\} < 1 . This implies the presence of contractile properties in the system, which are required for the solutions to stabilize and converge. It is also noteworthy that the system, as shown by Eq (6), has a consistent expression. For the overall validity and stability of the model to be maintained, the solutions must be continuous for them to behave smoothly over the designated domain.
\left.\begin{array}{c}{U}_{n}\left(t\right) = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}{\left(t-u\right)}^{\alpha -1}{\mathcal{P}}_{1}\left(u,{U}_{n-1}\right)du\\ {V}_{n}\left(t\right) = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}{\left(t-u\right)}^{\alpha -1}{\mathcal{P}}_{2}\left(u,{V}_{n-1}\right)du\\ {W}_{n}\left(t\right) = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}{\left(t-u\right)}^{\alpha -1}{\mathcal{P}}_{3}\left(u,{W}_{n-1}\right)du\end{array}\right\} . | (8) |
The representation of the variation between the two terms in (8) is
\left.\begin{array}{c}{\psi }_{n-1}\left(t\right) = \left({U}_{n}\left(t\right)-{U}_{n-1}\left(t\right)\right) = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}\left({\mathcal{P}}_{1}\left(u,{U}_{n-1}\right)-{\mathcal{P}}_{1}\left(u,{U}_{n-2}\right)\right)du\\ {\varphi }_{n-1}\left(t\right) = \left({V}_{n}\left(t\right)-{V}_{n-1}\left(t\right)\right) = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}\left({\mathcal{P}}_{2}\left(u,{V}_{n-1}\right)-{\mathcal{P}}_{2}\left(u,{V}_{n-2}\right)\right)du\\ {\vartheta }_{n-1}\left(t\right) = \left({W}_{n}\left(t\right)-{W}_{n-1}\left(t\right)\right) = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}\left({\mathcal{P}}_{3}\left(u,{W}_{n-1}\right)-{\mathcal{P}}_{3}\left(u,{W}_{n-2}\right)\right)du\end{array}\right\} . | (9) |
Therefore, we have
{U}_{n}\left(t\right) = \sum _{i = 0}^{n}{\psi }_{i}\left(t\right),{V}_{n}\left(t\right) = \sum _{i = 0}^{n}{\varphi }_{i}\left(t\right),{W}_{n}\left(t\right) = \sum _{i = 0}^{n}{\vartheta }_{i}\left(t\right). | (10) |
Let,
‖{\psi }_{n}\left(t\right)‖ = ‖{U}_{n}\left(t\right)-{U}_{n-1}\left(t\right)‖ . |
‖{\psi }_{n}\left(t\right)‖ = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}\left({\mathcal{P}}_{1}\left(u,{U}_{n-1}\right)-{\mathcal{P}}_{1}\left(u,{U}_{n-2}\right)\right)du . |
‖{\psi }_{n}\left(t\right)‖ = \frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}‖{U}_{n}\left(t\right)-{U}_{n-1}\left(t\right)‖du . |
‖{\psi }_{n}\left(t\right)‖ = \frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}{\psi }_{n-1}\left(t\right)du . | (11) |
When the remaining equations of system (9) are solved using the same approach, the following outcomes are obtained:
\left.\begin{array}{c}‖{\varphi }_{n}\left(t\right)‖ = \frac{{\xi }_{2}}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}{\varphi }_{n-1}\left(t\right)du\\ ‖{\vartheta }_{n}\left(t\right)‖ = \frac{{\xi }_{3}}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}{\vartheta }_{n-1}\left(t\right)du\end{array}\right\} . | (12) |
As required.
Theorem 2. \left(i\right) There exists a uniform function defined in system (9).
(ⅱ) If there is a {t}_{*} > 1 such that \frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)} < 1, then the model system (2) has at least one solution, for i = \mathrm{12, 3} and {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3.
Proof. The functions U\left(t\right), V\left(t\right), \;\mathrm{a}\mathrm{n}\mathrm{d}\;W\left(t\right) are bounded, and each kernel {\mathcal{P}}_{i} for i = \mathrm{12, 3}, fulfills the Lipschitz condition, leading to the derivation of the following relations.
\left.\begin{array}{c}‖{\psi }_{n}\left(t\right)‖\le ‖U\left(0\right)‖{‖\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)‖}^{n}\\ ‖{\varphi }_{n}\left(t\right)‖\le ‖V\left(0\right)‖{‖\frac{{\xi }_{2}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)‖}^{n}\\ ‖{\vartheta }_{n}\left(t\right)‖\le ‖W\left(0\right)‖{‖\frac{{\xi }_{3}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)‖}^{n}\end{array}\right\} . | (13) |
In the system (13) demonstration, the function given in (10) is demonstrated to exist and be uniform.
To prove (ⅱ), one must show that the system of solutions of (2) is reached by U\left(t\right), V\left(t\right)\;\mathrm{a}\mathrm{n}\mathrm{d}\;W\left(t\right) . Specifically, we define {A}_{n}\left(t\right), {B}_{n}\left(t\right), \;\mathrm{a}\mathrm{n}\mathrm{d}\;{C}_{n}\left(t\right) as terms that function after n variations are performed, so that
\left.\begin{array}{c}U\left(t\right)-U\left(0\right) = {U}_{n}\left(t\right)-{A}_{n}\left(t\right)\\ V\left(t\right)-V\left(0\right) = {V}_{n}\left(t\right)-{B}_{n}\left(t\right)\\ W\left(t\right)-W\left(0\right) = {W}_{n}\left(t\right)-{C}_{n}\left(t\right)\end{array}\right\} . | (14) |
Utilizing the triangle inequality in conjunction with {\xi }_{1} Lipschitz condition, we conclude that
‖{A}_{n}\left(t\right)‖ = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}\left({\mathcal{P}}_{1}\left(u,{U}_{n-1}\right)-{\mathcal{P}}_{1}\left(u,{U}_{n-2}\right)\right)du . |
‖{A}_{n}\left(t\right)‖\le \frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}‖{U}_{n}\left(t\right)-{U}_{n-1}\left(t\right)‖ . | (15) |
Repetitively carrying out the procedure in (15), we get
‖{A}_{n}\left(t\right)‖\le {‖\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)‖}^{n+1}{\mathsf{ε}}_{1} . | (16) |
Next, at {t}_{*} , one acquires
‖{A}_{n}\left(t\right)‖\le {‖\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left({t}_{*}\right)‖}^{n+1}{\mathsf{ε}}_{1} . | (17) |
Assuming n\to \infty as the limit.
\underset{n\to \infty }{\mathrm{lim}}‖{A}_{n}\left(t\right)‖\le \underset{n\to \infty }{\mathrm{lim}}{‖\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left({t}_{*}\right)‖}^{n+1}{\mathsf{ε}}_{1} . | (18) |
By applying the hypothesis \frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left({t}_{*}\right) < 1 , we have from (18) that
\underset{n\to \infty }{\mathrm{lim}}‖{A}_{n}\left(t\right)‖ = 0 . | (19) |
By using the same process as for n→∞, we get
\left.\begin{array}{c}‖{B}_{n}\left(t\right)‖\to 0\\ ‖{C}_{n}\left(t\right)‖\to 0\end{array}\right\} . | (20) |
Therefore, there is certainly a single solution.
Theorem 3. If \left(1-\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)\right) > 0 and {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3, then the system (2) has a unique solution.
Proof. Consider that another collection of solutions to (2) is represented by the sets {U}_{1}, {V}_{1}, \;\mathrm{a}\mathrm{n}\mathrm{d}\;{W}_{1} .
‖U\left(t\right)-{U}_{1}\left(t\right)‖ = \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_{0}^{t}\left({\mathcal{P}}_{1}\left(u,U\right)-{\mathcal{P}}_{1}\left(u,{U}_{1}\right)\right)du . |
‖U\left(t\right)-{U}_{1}\left(t\right)‖\le \frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}‖U\left(t\right)-{U}_{1}\left(t\right)‖ . | (21) |
If the terms in (21), are rearranged, one gets
\left(1-\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)\right)‖U\left(t\right)-{U}_{1}\left(t\right)‖\le 0 . | (22) |
By applying the hypothesis \left(1-\frac{{\xi }_{1}}{\mathrm{\Gamma }\left(\alpha \right)}\left(t\right)\right) > 0 , we have from (22) that
‖U\left(t\right)-{U}_{1}\left(t\right)‖ = 0 . | (23) |
It follows from this because U\left(t\right) = {U}_{1}\left(t\right) .
Applying the identical process to every solution for i = \mathrm{2, 3}, we arrive at:
V\left(t\right) = {V}_{1}\left(t\right) , W\left(t\right) = {W}_{1}\left(t\right) . | (24) |
Hence, the theorem is proved.
Theorem 4. For {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 , the shown stochastic fractional delayed heroin model (2) has a positive solution in {\mathbb{R}}^{+3} based on its initial conditions.
Proof. Whenever the states are possibly constrained by the initial conditions, the states must remain non-negative throughout the system. This means that the system's variables and parameters, unless specified, cannot be negative at any point in the evolutionary framework mentioned above. Non-negativity is important for convolution's sake as negative states will not make much real sense or even be well-defined mathematically in many cases in most physical models. We get
{\left.{}_{0}{}^{c}{D}_{t}^{\alpha }\left[U\left(t\right)\right]\right|}_{U = 0} = {\mathrm{\Lambda }}^{\alpha }\ge 0 , {\left.{}_{0}{}^{c}{D}_{t}^{\alpha }\left[V\left(t\right)\right]\right|}_{V = 0} = 0\ge 0 , {\left.{}_{0}{}^{c}{D}_{t}^{\alpha }\left[W\left(t\right)\right]\right|}_{W = 0} = {\rho }^{\alpha }V\ge 0 . | (25) |
The positive solution is achieved based on the given stochastic fractional delayed heroin model as described in the system (2) if the initial condition is located in the feasible region. This means that even in the presence of stochastic randomness and fractional time intervals, the system is possible and will not go beyond the realms of physical or biological reality to arrive at a positive outcome. The stability of the model and the accuracy with which it describes the dynamics of heroin addicts' behavior are thus retained; this means that the position of initial parameters in the allowable region is defined correctly.
Theorem 5. For any time t, the system (2) in feasible region \mathcal{H} = \left\{\left(U\left(t\right), V\left(t\right), W\left(t\right)\right)\in {\mathbb{R}}^{+3};0 < N\left(t\right)\le \frac{{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }}, \forall t\ge 0, \tau < t\right\} (where N\left(t\right) = U\left(t\right)+V\left(t\right)+W\left(t\right) is the total human population) with the initial condition is bounded assuming that {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 .
Proof. The total sum of human populations can be written as;
{}_{0}{}^{c}{D}_{t}^{\alpha }N\left(t\right)\le {\mathrm{\Lambda }}^{\alpha }-{\mu }_{1}^{\alpha }N\left(t\right) . | (26) |
After resolving the inequality above, we obtain
N\left(t\right)\le \frac{{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }}+\left(N\left(0\right)-\frac{{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }}\right){e}^{-\mu t} . | (27) |
Using Grown's inequality
\underset{t\to \infty }{\mathrm{lim}}\mathit{Sup}N\left(t\right)\le \frac{{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }} . | (28) |
Therefore, the epidemiologically feasible region for the propagation of heroin is provided by (28).
\mathcal{H} = \left\{\left(U\left(t\right),V\left(t\right),W\left(t\right)\right)\in {\mathbb{R}}^{+3};0 < N\le \frac{{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }}, \forall t\ge 0,\tau < t\right\} . | (29) |
From an epidemiological perspective, when analyzing and predicting the epidemic characteristics of heroin addiction, it is reasonable to use the stochastic fractional delayed heroin model (2). This model is suitable and stable to analogize the dissemination of heroin use among people, providing adequate points of view on the key parameters of the distribution course. According to the report made in (29), this positive invariance means that the structure of the model can preserve the epidemiological parameters' integrity in its usage.
This section examines the intricate stages of the stochastic fractional delayed heroin model (2) in connection with the dynamics of heroin infection. It provides a comprehensive examination of the system's behavior under different conditions, highlighting the heroin-present equilibrium a state in which the infection persists, and the heroin-free equilibrium a state in which the disease is eradicated. By assessing these equilibrium states, the section aims to shed information on the underlying mechanisms that limit the infection's tendency to spread throughout a population. It also considers the effects of stochastic components, fractional orders, and delays on the stability and transition between these equilibria, giving insight into the effectiveness of various intervention strategies.
Therefore,
\text{Heroin-free equilibrium} = HFE = {\mathcal{O}}^{0} = \left({U}_{0},{V}_{0},{W}_{0}\right) = \left(\frac{{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }},\mathrm{0,0},0\right) . | (30) |
\text{Heroin-present equilibrium} = HPE = {\mathcal{O}}^{\mathcal{*}} = \left({U}^{*},{V}^{*},{W}^{*}\right) . | (31) |
{U}^{*} = \frac{{\rho }^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)-{\beta }_{3}^{\alpha }\left(1+{\beta }_{3}^{\alpha }\right)}{{\beta }_{1}^{\alpha }{\rho }^{\alpha }} , {V}^{*} = \frac{{\mathrm{\Lambda }}^{\alpha }-{\mu }_{1}^{\alpha }{U}^{*}}{{\beta }_{1}^{\alpha }{U}^{*}{e}^{-{\mu }_{1}^{\alpha }\tau }} , {W}^{*} = \frac{{\beta }_{3}^{\alpha }{V}^{*}+{\mu }_{1}^{\alpha }}{{\rho }^{\alpha }{V}^{*}} . |
In epidemiology, the reproduction number holds great significance. This implies whether or not the illness is common among people in general. In a population where the reproduction number is less than one, disease may be averted; if it is greater than one, disease already exists. To determine the reproduction number, we apply next-generation methodology. \mathcal{L} therefore, becomes the matrix of transmission, and \mathcal{M} the matrix of transition.
\mathcal{L} = \left[\begin{array}{cc}\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }}& 0\\ 0& 0\end{array}\right] , \mathcal{M} = \left[\begin{array}{cc}\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)& 0\\ -{\rho }^{\alpha }& {\mu }_{1}^{\alpha }\end{array}\right] . |
The largest eigenvalue of the matrix in heroin-free equilibrium, also known as the spectral radius or reproduction number, is as follows:
{R}_{0} = \frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} . | (32) |
Theorem 6. Heroin-free equilibrium (30) is locally asymptotically stable for \alpha \in (0, 1) if {R}_{0} < 1 assuming that {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 .
Proof. The stochastic fractional delayed heroin model (2) is linearized around (30) to get a 3\times 3 dimensional Jacobian matrix with negative real components and eigenvalues:
{\lambda }_{1} = -{\mu }_{1}^{\alpha },{\lambda }_{2} = -{\mu }_{1}^{\alpha } . {\lambda }_{3} = \frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }}{{\mu }_{1}^{\alpha }}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right) . |
Therefore, if {R}_{0} < 1 , then the heroin-free equilibrium of the given stochastic fractional delayed heroin model (2) is stable locally. Otherwise, (30) is unstable in the local sense if {R}_{0} > 1 .
Theorem 7. The stochastic fractional delayed heroin model (2) is globally asymptotically stable (GAS) at heroin-free equilibrium, {\mathcal{O}}^{0} if {R}_{0} < 1 with the assumption of {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 .
Proof. Define the Lyapunov function \mathcal{L}:\mathcal{H}\to \mathbb{R} as, {\mathcal{O}}^{0} = \left({U}_{0}, {V}_{0}, {W}_{0}\right)
\mathcal{L}\left(t\right) = \left[U-{U}_{0}-{U}_{0}log\frac{U}{{U}_{0}}\right]+{V}_{0}+{W}_{0} . |
{}_{0}{}^{c}{D}_{t}^{\alpha }\mathcal{L}\left(t\right) = \left[\frac{U-{U}_{0}}{U}\right]{D}_{t}^{\alpha }U+{D}_{t}^{\alpha }V+{D}_{t}^{\alpha }W. |
\begin{align} {}_{0}{}^{c}{D}_{t}^{\alpha }\mathcal{L}\left(t\right) =& \left[\frac{U-{U}_{0}}{U}\right]\left({\mathrm{\Lambda }}^{\alpha }-{\beta }_{1}^{\alpha }U\left(t\right)V\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }-{\mu }_{1}^{\alpha }U\left(t\right)\right)+\\&\qquad \left({\beta }_{1}^{\alpha }U\left(t\right)V\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }-{\rho }^{\alpha }V\left(t\right)+{\beta }_{3}^{\alpha }V\left(t\right)W\left(t\right)-{\mu }_{1}^{\alpha }V\left(t\right)\right) \end{align} |
+\left({\rho }^{\alpha }V\left(t\right)-{\beta }_{3}^{\alpha }V\left(t\right)W\left(t\right)-{\mu }_{1}^{\alpha }W\left(t\right)\right) . |
{}_{0}{}^{c}{D}_{t}^{\alpha }\mathcal{L}\left(t\right)\le -{\mathrm{\Lambda }}^{\alpha }\frac{{\left(U-{U}_{0}\right)}^{2}}{U{U}_{0}}-{\mu }_{1}^{\alpha }V\left(t\right)\left[1-\frac{{\beta }_{1}^{\alpha }U\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }}\right]-{\mu }_{1}^{\alpha }W\left(t\right) . |
This implies that {}_{0}{}^{c}{D}_{t}^{\alpha }\mathcal{L}\le 0 if {R}_{0} < 1 and {}_{0}{}^{c}{D}_{t}^{\alpha }\mathcal{L} = 0 if U\left(t\right) = {U}_{0}, V\left(t\right) = W\left(t\right) = 0 . Therefore, {\mathcal{O}}^{0} is globally asymptotically stable.
Theorem 8. Heroin-present equilibrium (31) is locally asymptotically stable for \alpha \in (0, 1) if {R}_{0} > 1 with the assumption that {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 .
Proof. The stochastic fractional delayed heroin model (2) is linearized around (31) to get a 3\times 3 dimensional Jacobian matrix with negative real components and eigenvalues
{\lambda }_{1} = -\left({\beta }_{1}^{\alpha }{V}^{*}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right) . |
The second-order polynomial follows as
{\lambda }^{2}+\left(E+AB\right)\lambda +\left(C-DAB+AE\right) = 0 . |
Here,
\begin{gathered} A = \left({\beta }_{3}^{\alpha }{V}^{*}+{\mu }_{1}^{\alpha }\right) > 0 , B = {\beta }_{1}^{\alpha }{U}^{*}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{W}^{*}+{\beta }_{1}^{\alpha }{U}^{*}{V}^{*}{e}^{-{\mu }_{1}^{\alpha }\tau } > 0 , C = {\rho }^{\alpha }{\beta }_{1}^{\alpha }{V}^{*}{e}^{-{\mu }_{1}^{\alpha }\tau } > 0 , \\ D = {\beta }_{3}^{\alpha }{W}^{*}{V}^{*} > 0 , E = \left({\mu }_{1}^{\alpha }+{\rho }^{\alpha }\right) > 0 . \end{gathered} |
Which is the second-order polynomial where the coefficients of the polynomial are positive with E+AB > 0, CDAB+AE > 0 . So, by the Routh-Hurwitz Criterion for the second-degree polynomial, the coefficient of the characteristic equation is positive with the constraint {R}_{0} > 1. Hence, the heroin-present equilibrium of the given system (2) is locally stable. If {R}_{0} < 1, then Routh Hurwitz's condition for stability is violated. Thus, (31) is locally unstable.
Theorem 9. The stochastic fractional delayed heroin model (2) is globally asymptotically stable (GAS) at heroin-present equilibrium, ( {\mathcal{O}}^{\mathcal{*}} ) if {R}_{0} > 1 with assumption of {\sigma }_{i} = 0;i = \mathrm{1, 2}, 3 .
Proof. Define the Lyapunov function Z:\mathcal{H}\to \mathbb{R} as,
Z = {k}_{1}\left(U-{U}^{*}-{U}^{*}\mathit{ln}\left(\frac{U}{{U}^{*}}\right)\right)+{k}_{2}\left(V-{V}^{*}-{V}^{*}\mathit{ln}\left(\frac{V}{{V}^{*}}\right)\right)+{k}_{3}\left(W-{W}^{*}-{W}^{*}\mathit{ln}\left(\frac{W}{{W}^{*}}\right)\right) . |
Given positive constants {k}_{i}(i = \mathrm{1, 2}, 3) , we can express the following equation:
{}_{0}{}^{c}{D}_{t}^{\alpha }Z = {k}_{1}\left[\frac{U-{U}^{*}}{U}\right]{D}_{t}^{\alpha }U+{k}_{2}\left[\frac{V-{V}^{*}}{V}\right]{D}_{t}^{\alpha }V+{k}_{3}\left[\frac{W-{W}^{*}}{W}\right]{D}_{t}^{\alpha }W . |
{}_{0}{}^{c}{D}_{t}^{\alpha }Z = -{k}_{1}{\mathrm{\Lambda }}^{\alpha }\frac{{\left(U-{U}^{*}\right)}^{2}}{U{U}^{*}}-{k}_{2}\left({\beta }_{1}^{\alpha }U\left(t\right)V\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }V\left(t\right)W\left(t\right)\right)\frac{{\left(V-{V}^{*}\right)}^{2}}{V{V}^{*}}-{k}_{3}{\rho }^{\alpha }V\left(t\right)\frac{{\left(W-{W}^{*}\right)}^{2}}{W{W}^{*}} . |
If we choose {k}_{i}\;where\;\left(i = \mathrm{1, 2}, 3\right)
{}_{0}{}^{c}{D}_{t}^{\alpha }Z = -{\mathrm{\Lambda }}^{\alpha }\frac{{\left(U-{U}^{*}\right)}^{2}}{U{U}^{*}}-\left({\beta }_{1}^{\alpha }U\left(t\right)V\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }V\left(t\right)W\left(t\right)\right)\frac{{\left(V-{V}^{*}\right)}^{2}}{V{V}^{*}}-{\rho }^{\alpha }V\left(t\right)\frac{{\left(W-{W}^{*}\right)}^{2}}{W{W}^{*}} . |
{}_{0}{}^{c}{D}_{t}^{\alpha }Z\le 0 , for {R}_{0} > 1 and {}_{0}{}^{c}{D}_{t}^{\alpha }Z = 0 if and only if U = {U}^{*}, V = {V}^{*}, W = {W}^{*} . Hence Lasalle's invariance principle \left(31\right) is globally asymptotically stable.
Definition. Let B\left(t\right) be a Brownian motion and Y\left(t\right) be an Itô drift-diffusion process that satisfies the stochastic fractional delay differential equation:
{}_{0}{}^{c}{D}_{t}^{\alpha }\left[V\left(t\right)\right] = {\beta }_{1}^{\alpha }U\left(t\right)V\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }V\left(t\right)W\left(t\right)-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)V\left(t\right)+{\sigma }_{2}V\left(t\right)d\left(B\left(t\right)\right) . |
{}_{0}{}^{c}{D}_{t}^{\alpha }\left[V\left(t\right)\right]\le {\beta }_{1}^{\alpha }U\left(t\right)V\left(t\right){e}^{-{\mu }_{1}^{\alpha }\tau }-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)V\left(t\right)+{\sigma }_{2}V\left(t\right)d\left(B\left(t\right)\right) . |
If f\left(V, t\right)\in {C}^{2}\left({\mathbb{R}}^{2}, \mathbb{ }\mathbb{R}\right), then f(V\left(t\right), t) is also an Itô drift-diffusion process, which satisfies as follows:
{}_{0}{}^{c}{D}_{t}^{\alpha }\left(f\left(V\left(t\right),t\right)\right) = {}_{0}{}^{c}{D}_{t}^{\alpha }f\left(V\left(t\right),t\right)dt+{}_{0}{}^{c}{D}_{t}^{\alpha }\left(\left(V\left(t\right),t\right)\right)dB\left(t\right)+\frac{1}{2}{}_{0}{}^{c}{{D}^{2}}_{t}^{\alpha }\left(\left(V\left(t\right),t\right)\right){dB\left(t\right)}^{2}. |
Let us introduce {R}_{o}^{s} = {\mathrm{R}}_{\mathrm{o}}^{\mathrm{d}}-\frac{{\mathrm{\sigma }}_{2}^{2}}{2\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}.
Lemma. The unique solution of the system (2) lies in the region \mathcal{H} if it satisfies \left(U\left(0\right), V\left(0\right), W\left(0\right)\right)\in {\mathbb{R}}^{+3} .
Definition. The infected individuals will be extinct in the system (2), if \underset{t\to \infty }{\mathit{lim}}V\left(\mathrm{t}\right) = 0 , \forall t\ge 0 .
Theorem 10. If {R}_{o}^{d} < 1 and {\mathrm{\sigma }}_{2}^{2} < \frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} , then the infected individuals of the system (2) exponentially tend to zero.
Proof. Let us consider the initial data \left(U\left(0\right), V\left(0\right), W\left(0\right)\right)\in {\mathbb{R}}^{+3} and the system (2) admits the solution as \left(U\left(t\right), V\left(t\right), W\left(t\right)\right) , with σ and c being the randomness and drift, respectively, if it satisfies the stochastic fractional delay differential equation:
{}_{0}{}^{c}{D}_{t}^{\alpha }\left[V\left(t\right)\right]\le \left(\frac{{\beta }_{1}^{\alpha }U{e}^{-{\mu }_{1}^{\alpha }\tau }}{\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right)dt+{\sigma }_{2}cV\left(t\right)d\left(B\left(t\right)\right) . |
By using the Itô's lemma with f\left(V\right) = \mathit{ln}\left(V\right) , we have
{}_{0}{}^{c}{D}_{t}^{\alpha }ln\left(V\right) = {}_{0}{}^{c}{D}_{t}^{\alpha }\left(V\right)dV+\frac{1}{2}{}_{0}{}^{c}{{D}^{2}}_{t}^{\alpha }\left(V\right){V}^{2}{\mathrm{\sigma }}^{2}dt . |
{}_{0}{}^{c}{D}_{t}^{\alpha }\mathit{ln}\left(V\right) = \frac{1}{V}dV+\frac{1}{2}(-\frac{1}{{V}^{2}}){V}^{2}{\mathrm{\sigma }}^{2}dt . |
{}_{0}{}^{c}{D}_{t}^{\alpha }\mathit{ln}\left(V\right)\le (\frac{{\beta }_{1}^{\alpha }U{e}^{-{\mu }_{1}^{\alpha }\tau }}{\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)-\frac{1}{2}{\mathrm{\sigma }}_{2}^{2})dt+{\sigma }_{2}cd\left(B\left(t\right)\right) . |
\mathrm{ln}\left(V\right)\le \mathrm{ln}\left(V\left(0\right)\right)\le {\int }_{0}^{t}\left(\frac{{\beta }_{1}^{\alpha }U{e}^{-{\mu }_{1}^{\alpha }\tau }}{\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)-\frac{1}{2}{\mathrm{\sigma }}_{2}^{2}\right)dt+{\int }_{0}^{\mathrm{t}}{\sigma }_{2}cd\left(B\left(t\right)\right) . |
Notice that, W\left(t\right) = {\int }_{0}^{\mathrm{t}}{\sigma }_{2}cd\left(B\left(t\right)\right) with W\left(0\right) = 0 .
If {\mathrm{\sigma }}_{2}^{2} > \frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} ,
\mathrm{ln}\left(V\right)\ge \left(\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)-\frac{1}{2}\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}\right)t+W\left(t\right)+\mathit{ln}V\left(0\right) , |
\frac{\mathrm{ln}\left(V\right)}{\mathrm{t}}\ge \left(\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{2{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right)+\frac{\mathrm{W}\left(\mathrm{t}\right)}{\mathrm{t}}+\frac{\mathrm{ln}V\left(0\right)}{\mathrm{t}} , |
\underset{\mathrm{t}\to \mathrm{\infty }}{\mathrm{lim}}\frac{\mathrm{ln}\left(V\right)}{\mathrm{t}}\ge \left(\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{2{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right) > 0 , \text{with} \; \underset{\mathrm{t}\to \mathrm{\infty }}{\mathrm{lim}}\frac{\mathrm{W}\left(\mathrm{t}\right)}{\mathrm{t}} = 0. |
If {\mathrm{\sigma }}_{2}^{2} < \frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} , then
\mathrm{ln}\left(V\left(\mathrm{t}\right)\right)\le \left(\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)-\frac{1}{2}{\mathrm{\sigma }}_{2}^{2}\right)t+\mathrm{W}\left(\mathrm{t}\right)+\mathrm{ln}V\left(0\right), |
\frac{\mathrm{ln}\left(V\right)}{\mathrm{t}}\le \left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\left(\frac{{\beta }_{1}^{\alpha }{\mathrm{\Lambda }}^{\alpha }{e}^{-{\mu }_{1}^{\alpha }\tau }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}-1\right)+\frac{W\left(t\right)}{t}+\frac{\mathit{ln}V\left(0\right)}{t}, |
\underset{t\to \infty }{\mathit{lim}}sup\frac{\mathit{ln}\left(V\right)}{t}\le \left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\left({R}_{o}^{d}-1\right) , |
when {R}_{o}^{d} < 1 , we get \underset{t\to \infty }{\mathit{lim}}sup\frac{\mathit{ln}\left(V\right)}{t}\le 0 , \underset{t\to \infty }{lim}V\left(t\right) = 0 , as desired.
{R}_{o}^{d} = {R}_{o}^{s}-\frac{{\sigma }_{2}^{2}}{2\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} < 1 . |
In this section, we examine the behavior of model parameters concerning reproduction number {R}_{0} . Examine the transmission and spread of disease with the sensitive analysis of the model. Preliminary: The formalized sensitivity index of a variable ℮ depends on a parameter \varrho :
{\mathcal{E}}_{\varrho }^{℮} = \frac{\varrho }{℮}\times \frac{\partial ℮}{\partial \varrho } . |
In spatial terms, determine the sensitive indices of parameters concerning the reproduction number {R}_{0} .
{\mathbb{V}}_{{\mathrm{\Lambda }}^{\alpha }} = \frac{{\mathrm{\Lambda }}^{\alpha }}{{R}_{0}}\times \frac{\partial {R}_{0}}{\partial {\mathrm{\Lambda }}^{\alpha }} = 1 > 0 , {\mathbb{V}}_{{\beta }_{1}^{\alpha }} = \frac{{\beta }_{1}^{\alpha }}{{R}_{0}}\times \frac{\partial {R}_{0}}{\partial {\beta }_{1}^{\alpha }} = 1 > , {\mathbb{V}}_{{\rho }^{\alpha }} = \frac{{\rho }^{\alpha }}{{R}_{0}}\times \frac{\partial {R}_{0}}{\partial {\rho }^{\alpha }} = -\frac{1}{\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} > 0 , |
{\mathbb{V}}_{{\mu }_{1}^{\alpha }} = \frac{{\mu }_{1}^{\alpha }}{{R}_{0}}\times \frac{\partial {R}_{0}}{\partial {\mu }_{1}^{\alpha }} = -\frac{2{\mu }_{1}^{\alpha }+{\rho }^{\alpha }}{{\mu }_{1}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} < 0 . |
The uncertainty signs and values of the sensitivity indicators are presented in Tables 2 and 3. In Table 2, the positive signs of the uncertainty indicators of {R}_{0} correspond to the more sensitive parameters, whereas the non-positive signs indicate less sensitive parameters. Furthermore, if we control the transmission parameter {\beta }_{1}, which is very sensitive, then the disease will be controlled in the population (See Figure 2).
Parameters | Signs |
\mathbf{\Lambda } | Positive |
{\boldsymbol{\beta }}_{\mathbf{1}} | Positive |
\boldsymbol{\rho } | Negative |
{\boldsymbol{\mu }}_{\mathbf{1}} | Negative |
Parameters | Values |
\mathbf{\Lambda } | 1 |
{\boldsymbol{\beta }}_{\mathbf{1}} | 1 |
\boldsymbol{\rho } | -1.92 |
{\boldsymbol{\mu }}_{\mathbf{1}} | -3.92 |
This section presents a numerical method for the stochastic fractional delayed heroin model (2) that underlies it. Another instance of the stochastic fractional delayed system is provided:
\left.\begin{array}{c}{}_{0}{}^{c}{D}_{t}^{\alpha }\left[U\left(t\right)\right]{|}_{t = {t}_{n}} = {\mathrm{\Lambda }}^{\alpha }-{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }-{\mu }_{1}^{\alpha }{U}_{n+1}+{\sigma }_{1}{U}_{n}d\left(B\left(t\right)\right)\\ {}_{0}{}^{c}{D}_{t}^{\alpha }\left[V\left(t\right)\right]{|}_{t = {t}_{n}} = {\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }-{\rho }^{\alpha }{V}_{n+1}+{\beta }_{3}^{\alpha }{V}_{n}{W}_{n}-{\mu }_{1}^{\alpha }{V}_{n+1}+{\sigma }_{2}{V}_{n}d\left(B\left(t\right)\right)\\ {}_{0}{}^{c}{D}_{t}^{\alpha }\left[W\left(t\right)\right]{|}_{t = {t}_{n}} = {\rho }^{\alpha }{V}_{n}-{\beta }_{3}^{\alpha }{V}_{n}{W}_{n+1}-{\mu }_{1}^{\alpha }{W}_{n+1}+{\sigma }_{3}{W}_{n}d\left(B\left(t\right)\right)\end{array}\right\} . | (33) |
The first approach is the Grunwald-Letnikove approach, or GL, as follows:
{}_{0}{}^{c}{D}_{t}^{\alpha }\upsilon \left(t\right){|}_{t = {t}_{n}} = \frac{1}{{\left(\mathcal{K}\left(\Delta t\right)\right)}^{\alpha }}\left({\mathcal{K}}_{n+1}-\sum _{i = 1}^{n+1}{\nu }_{i}{\mathcal{K}}_{n+1-i}-{\rho }_{n+1}{\mathcal{K}}_{0}\right) . | (34) |
Here, {\nu }_{i} = {\left(-1\right)}^{i-1}\left(\begin{array}{c}\alpha \\ i\end{array}\right){\nu }_{1} = \alpha .
{\rho }_{i} = \frac{{i}^{-\alpha }}{\mathrm{\Gamma }\left(1-\alpha \right)};i = \mathrm{1,2},3,\dots n+1 . |
Now, the outcome that follows helps to verify some other hypotheses.
NSFD rules are added to the GL approach, making the discrete model for susceptible humans as follows:
\frac{1}{{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }}\left({U}_{n+1}-\sum _{i = 1}^{n+1}{\nu }_{i}{U}_{n+1-i}-{\rho }_{n+1}{U}_{0}\right) = {\mathrm{\Lambda }}^{\alpha }-{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }-{\mu }_{1}^{\alpha }{U}_{n+1}+{\sigma }_{1}{U}_{n}\Delta {B}_{n}. | (35) |
{U}_{n+1} = \frac{\sum _{i = 1}^{n+1}{\nu }_{i}{U}_{n+1-i}+{\rho }_{n+1}{U}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mathrm{\Lambda }}^{\alpha }+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\sigma }_{1}{U}_{n}\Delta {B}_{n}}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)} . | (36) |
Additionally, the latent and breaking out heroin model GL-NSFD scheme is as follows:
{V}_{n+1} = \frac{\sum _{i = 1}^{n+1}{\nu }_{i}{V}_{n+1-i}+{\rho }_{n+1}{V}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{3}^{\alpha }{V}_{n}{W}_{n}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\sigma }_{2}{V}_{n}\Delta {B}_{n}}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)} . | (37) |
{W}_{n+1} = \frac{\sum _{i = 1}^{n+1}{\nu }_{i}{W}_{n+1-i}+{\rho }_{n+1}{W}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }{V}_{n}\right)+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\sigma }_{3}{W}_{n}\Delta {B}_{n}}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)} . | (38) |
The positivity and boundedness of the solution for the systems (2) are confirmed by the following theorem.
Theorem 11. Suppose that {U}_{0}\ge 0, {V}_{0}\ge 0, {W}_{0}\ge 0, {\mathrm{\Lambda }}^{\alpha }\ge 0, {\beta }_{1}^{\alpha }\ge 0, {\mu }_{1}^{\alpha }\ge 0, {\rho }^{\alpha }\ge 0, {d}_{1}^{\alpha }\ge 0, {m}_{1}^{\alpha }\ge 0, then {U}_{n}\ge 0, {V}_{n}\ge 0, {W}_{n}\ge 0 for all n = \mathrm{1, 2}, 3, \dots .. and \Delta {B}_{n} = 0 .
Proof. For this, by using the induction method, we get,
For n = 0
{U}_{1} = \frac{{\nu }_{1}{U}_{0}+{\rho }_{1}{U}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mathrm{\Lambda }}^{\alpha }}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{V}_{1} = \frac{{\nu }_{1}{V}_{0}+{\rho }_{1}{V}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{U}_{1}{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{0}{W}_{0}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{W}_{1} = \frac{{\nu }_{1}{W}_{0}+{\rho }_{1}{W}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }{V}_{0}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)}\ge 0 . |
For n = 1
{U}_{2} = \frac{{\nu }_{1}{U}_{1}+{\nu }_{2}{U}_{0}+{\rho }_{2}{U}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mathrm{\Lambda }}^{\alpha }}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{V}_{2} = \frac{{\nu }_{1}{V}_{1}+{\nu }_{2}{V}_{0}+{\rho }_{2}{V}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{U}_{2}{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{1}{W}_{1}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{W}_{2} = \frac{{\nu }_{1}{W}_{1}+{\nu }_{2}{W}_{0}+{\rho }_{2}{W}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }{V}_{0}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)}\ge 0 . |
For n = 2
{U}_{3} = \frac{{\nu }_{1}{U}_{2}+{\nu }_{2}{U}_{1}+{\nu }_{3}{U}_{0}+{\rho }_{3}{U}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mathrm{\Lambda }}^{\alpha }}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{V}_{3} = \frac{{\nu }_{1}{V}_{2}+{\nu }_{2}{V}_{1}+{\nu }_{3}{V}_{0}+{\rho }_{3}{V}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{U}_{3}{V}_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{2}{W}_{2}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{W}_{3} = \frac{{\nu }_{1}{W}_{2}+{\nu }_{2}{W}_{1}+{\nu }_{3}{W}_{0}+{\rho }_{3}{W}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }{V}_{0}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)}\ge 0 . |
Suppose that for n = 1, 2, 3, ......, n-1 , {U}_{n}\ge 0 , {V}_{n}\ge 0 , and {W}_{n}\ge 0 .
Thus for n = n ,
{U}_{n+1} = \frac{\sum _{i = 1}^{n+1}{\nu }_{i}{U}_{n+1-i}+{\rho }_{n+1}{U}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mathrm{\Lambda }}^{\alpha }}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{V}_{n+1} = \frac{\sum _{i = 1}^{n+1}{\nu }_{i}{V}_{n+1-i}+{\rho }_{n+1}{V}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{3}^{\alpha }{V}_{n}{W}_{n}}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)}\ge 0 . |
{W}_{n+1} = \frac{\sum _{i = 1}^{n+1}{\nu }_{i}{W}_{n+1-i}+{\rho }_{n+1}{W}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }{V}_{n}\right)}{1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)}\ge 0 . |
As required.
Theorem 12. Suppose that {U}_{0}+{V}_{0}+{W}_{0} = 1, {\mathrm{\Lambda }}^{\alpha }\ge 0, {\beta }_{1}^{\alpha }\ge 0, {\mu }_{1}^{\alpha }\ge 0, {\rho }^{\alpha }\ge 0, {d}_{1}^{\alpha }\ge 0, {m}_{1}^{\alpha }\ge 0 and {\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\ge 0, then {U}_{n}, {V}_{n}, {W}_{n} are all bounded for all n = \mathrm{1, 2}, 3, \dots.n, and \Delta {B}_{n} = 0 .
Proof. For this,
{U}_{n+1}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mu }_{1}^{\alpha }{U}_{n+1}+{V}_{n+1}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right){V}_{n+1}+{W}_{n+1} |
+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right){W}_{n+1} |
= \sum\limits _{i = 1}^{n+1}{\nu }_{i}{U}_{n+1-i}+{\rho }_{n+1}{U}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\mathrm{\Lambda }}^{\alpha }+\sum\limits _{i = 1}^{n+1}{\nu }_{i}{V}_{n+1-i}+{\rho }_{n+1}{V}_{0} |
+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }{\beta }_{3}^{\alpha }{V}_{n}{W}_{n}+\sum _{i = 1}^{n+1}{\nu }_{i}{W}_{n+1-i}+{\rho }_{n+1}{W}_{0}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }{V}_{n}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{n+1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{n+1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{n+1} |
= \sum _{i = 1}^{n+1}{\nu }_{i}\left({U}_{n+1-i}+{V}_{n+1-i}+{W}_{n+1-i}\right)+{\rho }_{n+1}\left({U}_{0}+{V}_{0}+{W}_{0}\right)+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{n+1}{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{n}{W}_{n}\right) . |
Next, we use the induction method to evaluate the further iteration. Then,
For n = 0
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{1} |
= {\nu }_{1}\left({U}_{0}+{V}_{0}+{W}_{0}\right)+{\rho }_{1}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{1}{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{0}{W}_{0}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{1} |
= {\nu }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{1}{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{0}{W}_{0}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{1} = {\Upsilon }_{1} . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{1}\le {\nu }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{1}{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{0}{W}_{0}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{1}\le {\nu }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{1}{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{0}{W}_{0}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{1}\le {\nu }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{1}{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{0}{W}_{0}\right) . |
{U}_{1}\le \frac{{\Upsilon }_{1}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{0}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right)} . |
{V}_{1}\le \frac{{\Upsilon }_{1}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right)} . |
{W}_{1}\le \frac{{\Upsilon }_{1}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right)} . |
{U}_{1}\le {\Upsilon }_{1} , {V}_{1}\le {\Upsilon }_{1} , {W}_{1}\le {\Upsilon }_{1} . |
For n = 1
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{2} |
= {\nu }_{1}\left({U}_{1}+{V}_{1}+{W}_{1}\right)+{\nu }_{2}\left({U}_{0}+{V}_{0}+{W}_{0}\right)+{\rho }_{2}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{U}_{2}{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{V}_{1}{W}_{1}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{2} |
= \alpha \left(1\right)+\alpha \left({\Upsilon }_{1}+{\Upsilon }_{1}+{\Upsilon }_{1}\right)+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{1}{\Upsilon }_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{2}{\Upsilon }_{3}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{2} |
\le \alpha +3\alpha {\Upsilon }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{1}{\Upsilon }_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{2}{\Upsilon }_{3}\right) = {\Upsilon }_{2} . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{2}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{2} |
\le \alpha +3\alpha {\Upsilon }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{1}{\Upsilon }_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{2}{\Upsilon }_{3}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{2}\le \alpha +3\alpha {\Upsilon }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{1}{\Upsilon }_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{2}{\Upsilon }_{3}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{2}\le \alpha +3\alpha {\Upsilon }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{1}{\Upsilon }_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{2}{\Upsilon }_{3}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{2}\le \alpha +3\alpha {\Upsilon }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{1}{\Upsilon }_{2}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{2}{\Upsilon }_{3}\right) . |
{U}_{2}\le \frac{{\Upsilon }_{2}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right)} . |
{V}_{2}\le \frac{{\Upsilon }_{2}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right)} . |
{W}_{2}\left(t\right)\le \frac{{\Upsilon }_{2}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right)} . |
{U}_{2}\le {\Upsilon }_{2} , {V}_{2}\le {\Upsilon }_{2} , {W}_{2}\le {\Upsilon }_{2} . |
Now, consider that
{U}_{n}\le {\Upsilon }_{n},{V}_{n}\le {\Upsilon }_{n},{W}_{n}\le {\Upsilon }_{n} . |
Here,
{\Upsilon }_{n} = \alpha +3\alpha \left({\Upsilon }_{n-1},+{\Upsilon }_{n-2},\dots .{\Upsilon }_{2}+{\Upsilon }_{1}\right)+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\alpha +3\alpha {\Upsilon }_{1}+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)} |
+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
For n = n
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{n+1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{n+1}+\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{n+1} |
= \sum\limits _{i = 1}^{n+1}{\nu }_{i}\left({U}_{n+1-i}+{V}_{n+1-i}+{W}_{n+1-i}\right)+{\rho }_{n+1}\left({U}_{0}+{V}_{0}+{W}_{0}\right) |
+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
= {\nu }_{1}\left({U}_{n}+{V}_{n}+{W}_{n}\right)+{\nu }_{2}\left({U}_{n-1}+{V}_{n-1}+{W}_{n-1}\right)+{\nu }_{3}\left({U}_{n-2}+{V}_{n-2}+{W}_{n-2}\right)+\dots |
+{\nu }_{n}\left({U}_{1}+{V}_{1}+{W}_{1}\right)+{\nu }_{n+1}\left({U}_{0}+{V}_{0}+{W}_{0}\right)+{\rho }_{n+1}+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
\le \alpha \left(3{\Upsilon }_{n}\right)+\alpha \left(3{\Upsilon }_{n-1}\right)+\alpha \left(3{\Upsilon }_{n-2}\right)+\dots +\alpha \left(3{\Upsilon }_{2}\right)+\alpha \left(3{\Upsilon }_{1}\right)+\alpha \left(1\right)+\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)} |
+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
\le \alpha +\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+3\alpha \left({\Upsilon }_{n-1},+{\Upsilon }_{n-2},\dots .{\Upsilon }_{2}+{\Upsilon }_{1}\right)+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) = {\Upsilon }_{n+1} . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right){U}_{n+1} |
\le \alpha +\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+3\alpha \left({\Upsilon }_{n-1},+{\Upsilon }_{n-2},\dots .{\Upsilon }_{2}+{\Upsilon }_{1}\right)+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right){V}_{n+1} |
\le \alpha +\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+3\alpha \left({\Upsilon }_{n-1},+{\Upsilon }_{n-2},\dots .{\Upsilon }_{2}+{\Upsilon }_{1}\right)+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right){W}_{n+1} |
\le \alpha +\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}+3\alpha \left({\Upsilon }_{n-1},+{\Upsilon }_{n-2},\dots .{\Upsilon }_{2}+{\Upsilon }_{1}\right)+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\mathrm{\Lambda }}^{\alpha }+{\beta }_{1}^{\alpha }{\Upsilon }_{n}{\Upsilon }_{n+1}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\beta }_{3}^{\alpha }{\Upsilon }_{n+1}{\Upsilon }_{n+2}\right) . |
{U}_{n+1}\le \frac{{\Upsilon }_{n+1}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\beta }_{1}^{\alpha }{V}_{n}{e}^{-{\mu }_{1}^{\alpha }\tau }+{\mu }_{1}^{\alpha }\right)\right)} , {V}_{n+1}\le \frac{{\Upsilon }_{n+1}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({\rho }^{\alpha }+{\mu }_{1}^{\alpha }\right)\right)} , {W}_{n+1}\le \frac{{\Upsilon }_{n+1}}{\left(1+{\left(\mathcal{K}\left(h\right)\right)}^{\alpha }\left({d}_{1}^{\alpha }+{m}_{1}^{\alpha }\right)\right)} . |
{U}_{n+1}\le {\Upsilon }_{n+1} , {V}_{n+1}\le {\Upsilon }_{n+1} , {W}_{n+1}\le {\Upsilon }_{n+1} . |
The simulations' parameters are explained in this section. The primary features of the simulated graphs are examined using the set of parametric variables given in Table 4. Furthermore, these graphs are created at a time when heroin is broadly used in the human population and finally achieves a stable, endemic form. Four appropriate values of α are selected at the current equilibrium to examine the dynamics of the uninfected humans. The heroin epidemic has demonstrated that the rates at which infections spread vary throughout nations. Since every person or group of people has distinct physical surroundings, immune systems, health concerns, and other factors, these rates have biological significance.
Parameters | Values | Source [1] |
\mathbf{\Lambda } | 0.5 | Estimated |
{\boldsymbol{\beta }}_{\mathbf{1}} | 1.7 | Fitted |
{\boldsymbol{\mu }}_{\mathbf{1}} | 0.5 | Estimated |
{\boldsymbol{\beta }}_{\mathbf{3}} | 0.8 | Fitted |
\boldsymbol{\rho } | 0.02 | Fitted |
\boldsymbol{b} | 0.22 | Estimated |
{\boldsymbol{\sigma }}_{\boldsymbol{i}} | 0\le {\sigma }_{i}\le 1 | Fitted |
This section provides a detailed explanation of the comparison of the stochastic fractional delayed heroin epidemic model at different values of the fractional order \alpha and the effects of delay tactics \tau on susceptible, drug, and non-drug user populations. Figure 3 shows the behavior of the susceptible population over time for various fractional orders \alpha . The fractional order \alpha affects the memory and hereditary properties of the system. A lower \alpha indicates that the system has a stronger memory effect, which means that past states have a more significant influence on current dynamics. Therefore, Figure 3 shows that for lower values of \alpha , the susceptible population decreases more slowly, indicating that individuals remain at risk for longer periods. As \alpha increases, the decrease becomes more rapid, suggesting a quicker transition from susceptibility, potentially due to a faster adoption of protective behaviors or interventions. Figure 4 illustrates the dynamics of drug users over time for different fractional orders \alpha . A lower \alpha could indicate a slower decline or even a persistence in the number of drug users due to the long-term effects of addiction or the slow impact of rehabilitation efforts. Higher values of \alpha may show a quicker reduction in drug user numbers, indicating that the system is more responsive to current interventions or that the effects of past states (such as prior addiction) are fading faster. So, Figure 4 demonstrates that with increasing \alpha , the number of drug users decreases more rapidly, suggesting more effective rehabilitation or prevention efforts as the memory effect of addiction wanes. Figure 5 displays the behavior of non-drug users over time for various values of \alpha . The behavior of non-drug users is influenced by both the rate of recovery from drug use and the resistance to becoming drug users. For lower \alpha , non-drug users might increase more slowly, indicating a prolonged influence of past drug use or societal factors that delay recovery. As \alpha increases, the population of non-drug users could rise more quickly, reflecting a stronger impact of present conditions, such as effective awareness campaigns or policies that reduce drug use. Figure 5 shows a more significant increase in non-drug users at higher \alpha values, indicating more successful prevention and recovery efforts as the past influences diminish. Figure 6 shows the effect of a delay tactic \left(\tau = 10\right) on the susceptible population for different values of \alpha . The delay tactic represents an intervention or a time lag in the response to the epidemic (e.g., delayed introduction of a drug prevention program). With \left(\tau = 10\right) , the system responds with a lag to changes in the population or intervention efforts. Figure 6 shows that for lower \alpha , the effect of the delay is more pronounced, leading to prolonged periods of susceptibility. As \alpha increases, the system may recover more quickly from the delay, reflecting a reduced impact of the past. This indicates that higher fractional orders \alpha could mitigate the negative effects of delayed interventions, leading to a quicker reduction in the susceptible population. Figure 7 depicts the effect of a delay tactic \left(\tau = 10\right) on the drug user population for different values of \alpha . The delay tactic here could represent a delay in the implementation of drug rehabilitation programs or policies. For lower \alpha , the delay might result in a slower decrease or even a temporary increase in drug users, as the past states continue to influence the present. Higher \alpha values might show that the system adjusts more quickly, reducing the number of drug users even with the delay. Figure 7 shows that for higher \alpha , the negative impact of the delay is less severe, suggesting that systems with weaker memory effects (higher \alpha) are more resilient to delays in intervention. Figure 8 shows how the delay tactic \left(\tau = 10\right) affects the non-drug user population for different \alpha values. The delay could slow down the transition from drug user to non-drug user or delay the effects of prevention programs. For lower \alpha , the increase in non-drug users might be slower or less pronounced, reflecting the lingering effects of past drug use or the delayed response to interventions. As \alpha increases, the population of non-drug users might increase more rapidly despite the delay, indicating that the system is less affected by past states. The graph shows that higher \alpha values result in a more robust increase in non-drug users, even in the presence of delayed interventions, highlighting the importance of reducing the system's memory effect in drug prevention and recovery strategies. The graphs provide a comprehensive view of how the fractional order \alpha and delay tactics \tau influence the dynamics of the susceptible, drug, and non-drug user populations in a stochastic heroin epidemic model. Lower fractional orders indicate stronger memory effects, leading to slower responses to changes and interventions, while higher fractional orders result in quicker adjustments and more effective epidemic control. The delay tactics illustrate the critical timing of interventions, showing that systems with weaker memory effects are more resilient to delays.
The study of the stochastic fractional delayed heroin epidemic model is useful in understanding the phenomenon of epidemic growth of heroin addiction within the population. For this, the population of humans has been categorized into three groups: the probability density function of individuals most vulnerable to heroin usage based on their propensity for addiction, the quantitative future state of the affected community due to continual drug abuse, namely active heroin users and non-drug users. This categorization is used for evaluating risk exposure between drug takers and current heroin users. The dynamical properties of the stochastic fractional delay differential equations for the model were studied rigorously. A detailed stability analysis shows that the stability of the system depend on the basic reproduction number {R}_{0} . When {R}_{0} < 1 , the unique free equilibrium is stable, suggesting that addiction can be eliminated if conditions are right. On the other hand, if {R}_{0} > 1 , the present equilibrium is stable, implying that the spread of addiction will continue but can be controlled. It also examines how delay tactics affect the spread of addiction. The stochastic GL-NSFD scheme was implemented for numerical solutions of the model. Well-known theorems are presented for the positivity and boundedness of the proposed numerical method. Therefore, this stochastic fractional delayed model can be used as a useful tool in describing and modeling addiction impacts in communities. The delay strategy is of particular interest in providing methodical suggestions concerning the elaboration of successful health policies and preventive measures. In addition, this model can be used as the basis for further research studies of other diseases where a complex dynamic interaction occurs to evaluate the effectiveness of various control strategies.
All authors made substantial contributions to the study, participated in writing, and approved the manuscript. Also, read and approved a manuscript with the given study.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The datasets analyzed during the current study are available from the corresponding author upon reasonable request. No pharmaceutical study is based on the heroin epidemic model. Our focus is a non-pharmaceutical strategy to study the transmission dynamics of the heroin epidemic model. We take a model and data from reference number [1] and the model extends into stochastic with fractional delayed formulation.
This research was partially supported by the Fundação para a Ciência e Tecnologia, FCT, under the project https://doi.org/10.54499/UIDB/04674/2020.
The authors affirm they have no conflicts of interest to disclose concerning the current study.
[1] |
A. Raza, Y. M. Chu, M. Y. Bajuri, A. Ahmadian, N. Ahmed, M. Rafiq, et al., Dynamical and nonstandard computational analysis of heroin epidemic model, Results Phys., 34 (2022), 105245. https://doi.org/10.1016/j.rinp.2022.105245 doi: 10.1016/j.rinp.2022.105245
![]() |
[2] |
W. Weera, T. Botmart, S. Zuhra, Z. Sabir, M. A. Z. Raja, S. Ben Said, A neural study of the fractional Heroin epidemic model, Comput., Mater. Con., 74 (2022), 4453–4467. https://doi.org/10.32604/cmc.2023.033232 doi: 10.32604/cmc.2023.033232
![]() |
[3] |
Y. Zhou, C. Tian, Z. Ling, Dynamical behavior of the Heroin epidemic model on a finite weighted network, Bull. Malays. Math. Sci. Soc., 46 (2023), 175. https://doi.org/10.1007/s40840-023-01568-1 doi: 10.1007/s40840-023-01568-1
![]() |
[4] |
H. Jiang, L. Chen, F. Wei, Q. Zhu, Survival analysis and probability density function of switching heroin model, Math. Biosci. Eng., 20 (2023), 13222–13249. https://doi.org/10.3934/mbe.2023590 doi: 10.3934/mbe.2023590
![]() |
[5] |
Y. Wei, J. Zhan, J. Guo, Asymptotic behaviors of a heroin epidemic model with nonlinear incidence rate influenced by stochastic perturbations, J. Appl. Anal. Comput., 14 (2024), 1060–1077. https://doi.org/10.11948/20230323 doi: 10.11948/20230323
![]() |
[6] |
S. Djilali, Y. Chen, S. Bentout, Dynamics of a delayed nonlocal reaction–diffusion heroin epidemic model in a heterogenous environment, Math. Methods Appl. Sci., 48 (2024), 273–307. https://doi.org/10.1002/mma.10327 doi: 10.1002/mma.10327
![]() |
[7] |
X. Liu, M. Zhang, Z. W. Yang, Numerical threshold stability of a nonlinear age-structured reaction diffusion heroin transmission model, Appl. Numer. Math., 204 (2024), 291–311. https://doi.org/10.1016/j.apnum.2024.06.016 doi: 10.1016/j.apnum.2024.06.016
![]() |
[8] |
A. Chekroun, M. N. Frioui, T. Kuniya, T. M. Touaoula, Mathematical analysis of an age structured heroin-cocaine epidemic model, Discrete Cont. Dyn. Syst. -Ser. B, 25 (2020), 4449–4477. https://doi.org/10.3934/dcdsb.2020107 doi: 10.3934/dcdsb.2020107
![]() |
[9] |
L. Zhang, Y. Xing, Stability analysis of a reaction‐diffusion Heroin epidemic model, Complexity, 2020 (2020), 1–16. https://doi.org/10.1155/2020/3781425 doi: 10.1155/2020/3781425
![]() |
[10] |
S. Cole, S. Wirkus, Modeling the dynamics of heroin and illicit opioid use disorder, treatment, and recovery, Bul. Math. Biol., 84 (2022), 48. https://doi.org/10.1007/s11538-022-01002-w doi: 10.1007/s11538-022-01002-w
![]() |
[11] |
Z. U. A. Zafar, H. Rezazadeh, M. Inc, K. S. Nisar, T. A. Sulaiman, A. Yusuf, Fractional order heroin epidemic dynamics, Alex. Eng. J., 60 (2021), 5157–5165. https://doi.org/10.1016/j.aej.2021.04.039 doi: 10.1016/j.aej.2021.04.039
![]() |
[12] |
S. Rashid, F. Jarad, A. G. Ahmad, K. M. Abualnaja, New numerical dynamics of the heroin epidemic model using a fractional derivative with Mittag-Leffler kernel and consequences for control mechanisms, Results Phys., 35 (2022), 105304. https://doi.org/10.1016/j.rinp.2022.105304 doi: 10.1016/j.rinp.2022.105304
![]() |
[13] |
A. Khan, G. Zaman, R. Ullah, N. Naveed, Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age, AIMS Math., 6 (2021), 1377–1394. https://doi.org/10.3934/math.2021086 doi: 10.3934/math.2021086
![]() |
[14] |
S. Bentout, S. Djilali, B. Ghanbari, Backward, Hopf bifurcation in a heroin epidemic model with treat age, Int. J. Model., Simul., Sci. Comput., 12 (2021), 2150018. https://doi.org/10.1142/S1793962321500185 doi: 10.1142/S1793962321500185
![]() |
[15] |
A. Din, Y. Li, Controlling heroin addiction via age-structured modeling, Adv. Differ. Equ., 2020 (2020), 1–17. https://doi.org/10.1186/s13662-020-02983-5 doi: 10.1186/s13662-020-02983-5
![]() |
[16] |
S. Djilali, S. Bentout, T. M. Touaoula, A. Tridane, S. Kumar, Global behavior of Heroin epidemic model with time distributed delay and nonlinear incidence function, Results Phys., 31 (2021), 104953. https://doi.org/10.1016/j.rinp.2021.104953 doi: 10.1016/j.rinp.2021.104953
![]() |
[17] |
X. C. Duan, X. Z. Li, M. Martcheva, Coinfection dynamics of heroin transmission and HIV infection in a single population, J. Biol. Dyn., 14 (2020), 116–142. https://doi.org/10.1080/17513758.2020.1726516 doi: 10.1080/17513758.2020.1726516
![]() |
[18] |
P. T. Sowndarrajan, L. Shangerganesh, A. Debbouche, D. F. Torres, Optimal control of a heroin epidemic mathematical model, Optimization, 71 (2022), 3107–3131. https://doi.org/10.1080/02331934.2021.2009823 doi: 10.1080/02331934.2021.2009823
![]() |
[19] | A. Zinihi, M. R. S. Ammi, M. Ehrhardt, A. Bachir, Dynamical analysis of a cocaine-heroin epidemiological model with spatial distributions, arXiv Preprint, 2024. https://doi.org/10.48550/arXiv.2405.15532 |
[20] |
M. Vivas-Cortez, A. Bakar, M. S. Alqarni, N. Raza, T. Nazir, M. Farman, Global stability and modeling with a non-singular kernel for fractional order heroin epidemic model: insights from different population studies, J. King Saud Univ. -Sci., 36 (2024), 103329. https://doi.org/10.1016/j.jksus.2024.103329 doi: 10.1016/j.jksus.2024.103329
![]() |
[21] |
P. A. Naik, M. Farman, A. Zehra, K. S. Nisar, E. Hincal, Analysis and modeling with fractal-fractional operator for an epidemic model with reference to COVID-19 modeling, Partial Differ. Equ. Appl. Math., 10 (2024), 100663. https://doi.org/10.1016/j.padiff.2024.100663 doi: 10.1016/j.padiff.2024.100663
![]() |
[22] |
A. Zehra, P. A. Naik, A. Hasan, M. Farman, K. S. Nisar, F. Chaudhry, et al., Physiological and chaos effect on dynamics of neurological disorder with memory effect of fractional operator: a mathematical study, Comput. Meth. Prog. Biomed., 250 (2024), 108190. https://doi.org/10.1016/j.cmpb.2024.108190 doi: 10.1016/j.cmpb.2024.108190
![]() |
[23] |
P. A. Naik, B. M. Yeolekar, S. Qureshi, N. Manhas, M. Ghoreishi, M. Yeolekar, et al., Global analysis of a fractional‐order Hepatitis B virus model under immune response in the presence of cytokines, Adv. Theory Simul., 7 (2024), 2400726. https://doi.org/10.1002/adts.202400726 doi: 10.1002/adts.202400726
![]() |
[24] |
P. A. Naik, B. M. Yeolekar, S. Qureshi, M. Yeolekar, A. Madzvamuse, Modeling and analysis of the fractional-order epidemic model to investigate mutual influence in HIV/HCV co-infection, Nonlinear Dyn., 112 (2024), 11679–11710. https://doi.org/10.1007/s11071-024-09653-1 doi: 10.1007/s11071-024-09653-1
![]() |
[25] |
P. A. Naik, M. Yavuz, S. Qureshi, K. M. Owolabi, A. Soomro, A. H. Ganie, Memory impacts in hepatitis C: a global analysis of a fractional-order model with an effective treatment, Comput. Meth. Prog. Biomed., 254 (2024), 108306. https://doi.org/10.1016/j.cmpb.2024.108306 doi: 10.1016/j.cmpb.2024.108306
![]() |
Parameters | Description |
\mathbf{\Lambda } | The rate at which new people begin using heroin |
{\boldsymbol{\mu }}_{\mathbf{1}} | The rate at which current drug users discontinue their regular use |
{\boldsymbol{\beta }}_{\mathbf{1}} | The risk of infection for drug users |
{\boldsymbol{\beta }}_{\mathbf{3}} | The pace at which non-drug users enter the drug-user compartment |
\boldsymbol{\rho } | The mortality rate among drug users |
{\boldsymbol{e}}^{-{\boldsymbol{\mu }}_{\mathbf{1}}\boldsymbol{\tau }} | The artificial delay parameter |
Parameters | Signs |
\mathbf{\Lambda } | Positive |
{\boldsymbol{\beta }}_{\mathbf{1}} | Positive |
\boldsymbol{\rho } | Negative |
{\boldsymbol{\mu }}_{\mathbf{1}} | Negative |
Parameters | Values |
\mathbf{\Lambda } | 1 |
{\boldsymbol{\beta }}_{\mathbf{1}} | 1 |
\boldsymbol{\rho } | -1.92 |
{\boldsymbol{\mu }}_{\mathbf{1}} | -3.92 |
Parameters | Values | Source [1] |
\mathbf{\Lambda } | 0.5 | Estimated |
{\boldsymbol{\beta }}_{\mathbf{1}} | 1.7 | Fitted |
{\boldsymbol{\mu }}_{\mathbf{1}} | 0.5 | Estimated |
{\boldsymbol{\beta }}_{\mathbf{3}} | 0.8 | Fitted |
\boldsymbol{\rho } | 0.02 | Fitted |
\boldsymbol{b} | 0.22 | Estimated |
{\boldsymbol{\sigma }}_{\boldsymbol{i}} | 0\le {\sigma }_{i}\le 1 | Fitted |
Parameters | Description |
\mathbf{\Lambda } | The rate at which new people begin using heroin |
{\boldsymbol{\mu }}_{\mathbf{1}} | The rate at which current drug users discontinue their regular use |
{\boldsymbol{\beta }}_{\mathbf{1}} | The risk of infection for drug users |
{\boldsymbol{\beta }}_{\mathbf{3}} | The pace at which non-drug users enter the drug-user compartment |
\boldsymbol{\rho } | The mortality rate among drug users |
{\boldsymbol{e}}^{-{\boldsymbol{\mu }}_{\mathbf{1}}\boldsymbol{\tau }} | The artificial delay parameter |
Parameters | Signs |
\mathbf{\Lambda } | Positive |
{\boldsymbol{\beta }}_{\mathbf{1}} | Positive |
\boldsymbol{\rho } | Negative |
{\boldsymbol{\mu }}_{\mathbf{1}} | Negative |
Parameters | Values |
\mathbf{\Lambda } | 1 |
{\boldsymbol{\beta }}_{\mathbf{1}} | 1 |
\boldsymbol{\rho } | -1.92 |
{\boldsymbol{\mu }}_{\mathbf{1}} | -3.92 |
Parameters | Values | Source [1] |
\mathbf{\Lambda } | 0.5 | Estimated |
{\boldsymbol{\beta }}_{\mathbf{1}} | 1.7 | Fitted |
{\boldsymbol{\mu }}_{\mathbf{1}} | 0.5 | Estimated |
{\boldsymbol{\beta }}_{\mathbf{3}} | 0.8 | Fitted |
\boldsymbol{\rho } | 0.02 | Fitted |
\boldsymbol{b} | 0.22 | Estimated |
{\boldsymbol{\sigma }}_{\boldsymbol{i}} | 0\le {\sigma }_{i}\le 1 | Fitted |