Mathematical Modeling of Leukemia within Stochastic Fractional Delay Differential Equations

1  Introduction

Leukemia is one of the very significant global public health problems that cut across all demographics. This disease results from abnormal or premature white blood cells, which grow unsymmetrical in the blood and interfere with the normal functioning of healthy cells. This therefore impairs immunity to infections and general well-being [1]. In [2], the authors analyzed new pharmacodynamic parameters associated with ibrutinib responses in chronic lymphocytic leukemia by prospective study in real-world patients. Mathematical modeling is integrated to predict treatment outcomes. In [3], the authors develop nonlinear Ordinary Differential Equation (ODE) models describing the population dynamics of leukemic cells as functions of differences in feedback configurations and kinetic properties like self-renewal differentiation and division probabilities, proliferation, and death rates. The proposal extends to how these factors impact the course of leukemia. In [4], the authors attempted to mathematically describe the progression dynamics of chronic myeloid leukemia by providing a mathematical model of cloned hematopoiesis through nonlinear systems of differential equations. In [5], the authors integrated pharmacokinetic-pharmacodynamic (PKPD) models used to assess the clonal reduction potential of promising candidate drugs compared to high-dose cytarabine in the consolidation therapy of Acute Myeloid Leukemia (AML). Since the goal is to discover better alternatives. In [6], the authors described a simple model that describes an interface between leukemic cells and the body’s autologous immune response in the chronic phase of chronic myelogenous leukemia (CML). The model attempts to capture the dynamic behavior of the growth of leukemic cells coupled with the immune system response over time. In [7], the authors studied a mathematical model describing Chimeric Antigen Receptor T (CAR-T) cells, leukemia tumors, and B cell competition. All interactions are studied in detail to understand the essence of their behavior. In [8], the authors presented and discussed an autologous tumor-immune response model for CML. The paper advances toward a mathematical modeling understanding of the dynamics of CML. In [9], the authors looked at a mathematical myeloid leukemia model, emphasizing the existence and stability of trivial and nontrivial equilibrium points. Stability analysis is given for the equilibrium states. In [10], the authors developed a mathematical model to analyze the interactions between leukemia stem cells with the bone marrow microenvironment. This model enables the simulation of some dynamics in the course of chronic myeloid leukemia progression dynamics. In [11], the authors presented a novel two-parameter discrete distribution by combining Poisson and Quasi-Shanker distributions. This new distribution offers enhanced modeling capabilities. In [12], the authors analyzed Leukemia as a blood cancer characterized by an excess number of white blood cells. If identified at any stage with accuracy, the treatment would be effective. In [13], the authors contributed to the improvement of further knowledge, and accuracies in diagnostic methods, and also play an important role in present times for the diagnosis and treatment of acute leukemia. In [14], the authors demonstrated how to infer from experimental measurements with live cell fluorescence labeling and flow cytometry on the growth regime and division strategy of leukemia cell populations an analytical model for which cell growth and division rates depend on powers of the size. In [15], the authors develop a mathematical model to describe the acute lymphoblastic leukemia (ALL) behavior, which includes the evolution of a leukemic clone during the treatment process. In [16], the authors implemented a new model, simulating the evolution of cancerous cells in patients with chronic lymphocytic leukemia receiving chemotherapy treatment, to capture the dynamics of treatment response and disease progression. In [17], the authors introduced a nonlinear leukemia dynamical system using a piecewise modified ABC fractional-order derivative. Conduct both theoretical and computational analyses of the model, specifically regarding crossover effects. In [18], the authors explored Chronic lymphocytic leukemia (CLL) protein-protein interaction networks using a new approach where both statistical thermodynamics and systems biology are brought together to identify proteins integral to the onset of the disease. In [19], the authors proposed an approach to the pattern recognition of ALL through the application of advanced computational deep learning techniques. Advanced computational deep learning, in this case, will focus on improving the outcome of leukemia diagnosis with complex algorithms aimed at detecting and classifying abnormal cell patterns. In [20], the authors discussed optimal control in a reaction-diffusion leukemia immune model that captures the interactions and dynamics between leukemia cells, normal cells, and CAR-T cells, including how to control leukemia growth, enhance CAR-T therapy efficacy, and preserve healthy cell populations. In [21], the authors provided a comprehensive introduction to the mechanistic mathematical and computational modeling of blood cell formation depicting aspects of both normal and pathological processes, involving basic principles and applications of modeling techniques concerning disorders related to hematopoiesis. Alsakaji et al. studied a stochastic tumor-immune interaction model with external treatments and time delays, which is an optimal control problem [22]. Rihan et al. studied a fractional order delay differential model of a tumor-immune system with vaccine efficacy: stability, bifurcation, and control [23].

The stochastic fractional delayed methodology has significant scientific benefits because it combines stochastic processes, fractional calculus, and temporal delays to mimic complex systems more precisely. This approach uses fractional derivatives to capture memory effects, which is important in fields like biology and economics since it accounts for past effects on current conditions. Temporal delays are incorporated to provide a realistic depiction of the lag between cause and effect, while the stochastic component handles the inherent randomness in systems. Together, these elements enhance the models’ predictive power and robustness across multiple domains, rendering them more realistic representations of real-world dynamics.

For the computational modeling of the study on leukemia, stochastic fractional delay differential equations (SFDDEs) were used and implemented in MATLAB software with fractions calculus and stochastic modules. To solve the computational challenges, adaptive step-sizing and parallel processing were added while maximizing speed but maintaining stable solutions. The model was validated by comparing it with analytical solutions for simplified cases, ensuring the results are robust and reliable.

The structure of the paper is as follows: Section 1 provides an overview and a detailed examination of leukemia-like diseases that have been reported in the literature. Section 2 examines the construction of the stochastic fractional delayed leukemia disease model, the resulting mathematical analysis, and the local and global levels of the model’s equilibria, reproduction number, and stability analysis. The reproduction number sensitivity that we obtain from the Section 3 SFDDEs for the system. In Section 4, the stochastic fractional delayed model was examined using the Grunwald Letnikov Nonstandard Finite Difference (GL-NSFD) approach. In Section 5, the positivity and boundedness of the GL-NSFD were analyzed. In Section 6, the explicit focus is on numerical simulations and the results displayed. Final opinions offer a comprehensive synopsis of the work completed under Section 7.

2  Model Formulation

In this paper, we discussed a model that shows how leukemia spreads through three subgroups of cells. These groups are susceptible cells (S), infected cells (I), and immune cells (W). The susceptible group (S) includes cells that can get leukemia but don’t have it yet. The infected group (I) has cells that have leukemia and can spread it. The immune group (W) includes cells that have recovered from leukemia and are now immune, meaning they can’t get it again. Further, let Λ be the rate at which the susceptible blood cells enter the circulatory blood from compartments like bone marrow, lymph nodes, and thymus. Parameters μ, γ, and b are the natural mortality rates of susceptible blood cells, infected cells, and immune cells, respectively. The parameter β is the infection rate of susceptible blood cells. The rate at which the infected cells are recovered due to encounters with immune cells is denoted by ε.

A flow diagram of the model has been provided in Fig. 1. The basic structure for the investigations is the deterministic model introduced in [1], whereby the first-order temporal derivatives are substituted by fractional Caputo derivatives of order α. We believe that such an adjustment provides a more appropriate picture of the diffusion of leukemia disease. In the model analysis, Caputo fractional derivatives are employed to include memory on the fact that the system depends on previous user behavior, and it seems to store this information in long-term memory about other cases of leukemia. From this modification, one can see that the model is well suited to representing the processes of leukemia disease, especially over some time due to the non-linearity. The random or stochastic disturbances are incorporated to provide for the random nature and variation that is present in the systems under study. These stochastic components can facilitate and consider the complex spread of disease, as well as the unpredictable dynamics visible in the actual utilization of leukemia disease.

Dtα0c[S]=Λα−βαS(t−τ)I(t−τ)e−(μα+γα)τ−μαS+σ1αSdB(t),(1)

Dtα0c[I]=βαS(t−τ)I(t−τ)e−(μα+γα)τ−(γα+εα)I+σ2αIdB(t),(2)

Dtα0c[W]=εαI−bαW+σ3αWdB(t).(3)

Figure 1: Flow diagram of leukemia disease

The nonnegative (initial) conditions for the system (1)–(3) are

S(0)≥0,I(0)≥0,W(0)≥0,t≥0,τ<t,

and

N=S+I+W.

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 so that 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 deep and thorough understanding of the fractional derivative concept:

Definition 1: For a function 𝓆∈Cn, the Caputo fractional derivative of order α∈(𝓂−1,𝓂),𝓂∈N is

Dtα0c(t)=1Γ(𝓂−α)∫0t𝓆𝓂(T)dT(t−T)α+1−𝓂.

Definition 2: For thefunction 𝓆(t), the expression describes the equivalent fractional integral with order α>0.

Itα𝓆(t)=1Γ(α)∫0t(t−T)α−1𝓆(T)dT.

where “Γ” is the gamma function displayed.

2.1 Existence and Uniqueness of the Stochastic Fractional Delayed Model

This section of the paper establishes the stochastic fractional delayed model’s existence and uniqueness. Apply the fractional integral in this case, assuming σiα=0;i=1,2,3, under the system (1)–(3) starting condition.

S(t)=S0+1Γ(α)∫0t(t−s)α−1ℏ1(s,S)dS,(4)

I(t)=I0+1Γ(α)∫0t(t−s)α−1ℏ2(s,I)ds,(5)

W(t)=W0+1Γ(α)∫0t(t−s)α−1ℏ3(s,W)ds.(6)

The functions defined under the integral in system (4)–(6) are

ℏ1(t,F)=Λα−βαSIe−(μα+γα)τ−μαS,(7)

ℏ2(t,I)=βαSIe−(μα+γα)τ−(γα+εα)I,(8)

ℏ3(t,P)=εαI−bαW.(9)

Furthermore, it is assumed that ℰ1,ℰ2, and ℰ3 exist as positive constants and that S(t),I(t), and W(t) are non-negative limiting functions, such that

‖S(t)‖≤ℰ1,‖I(t)‖≤ℰ2,‖W(t)‖≤ℰ3.

Theorem 1: The function ℏi for i=1,2,3 fulfills the Lipshitz requirement and are contraction mappings by assuming σiα=0;i=1,2,3, if the condition 0≤𝒲=max{1,2,3}< 1 is true.

Proof: First, we consider the function ℏ1. Examine the subsequent S and S1:

‖ℏ1(t,S)−ℏ2(t,S1)‖=‖βα(S−S1)Ie−(μα+γα)τ+μα(S−S1)‖,‖ℏ1(t,S)−ℏ2(t,S1)‖≤‖βα(S−S1)Ie−(μα+γα)τ‖+‖μα(S−S1)‖,‖ℏ1(t,S)−ℏ2(t,S1)‖≤(βα(S−S1)e−(μα+γα)τ‖I‖+μα)‖S−S1‖,‖ℏ1(t,S)−ℏ2(t,S1)‖≤(βαe−(μα+γα)τℰ2+μα)‖S−S1‖,‖ℏ1(t,S)−ℏ2(t,S1)‖≤ξ1‖S−S1‖.(10)

In this instance, ξ1=(βαe−(μα+γα)τℰ2+μα). Lipshitz’s condition is satisfied. This method may also be used to meet the Lipschitz criteria for ℏi,i=2,3. Additionally, if 0≤𝒲=max{1,2,3}<1, then the functions are contractions. Furthermore, there is constant writing in the system (7)–(9).

Sn(t)=1Γ(α)∫0t(t−s)α−1ℏ1(s,Sn−1)ds,(11)

In(t)=1Γ(α)∫0t(t−s)α−1ℏ2(s,In−1)ds,(12)

Wn(t)=1Γ(α)∫0t(t−s)α−1ℏ3(s,Wn−1)ds.(13)

The two terms’ variation is represented by the following in (11)–(13):

ψn−1(t)=(Sn(t)−Sn−1(t))=1Γ(α)∫0t(ℏ1(s,Sn−1)−ℏ1(s,Sn−2))ds,(14)

φn−1(t)=(In(t)−In−1(t))=1Γ(α)∫0t(ℏ2(s,In−1)−ℏ2(s,In−2))ds,(15)

ϑn−1(t)=(Wn(t)−Wn−1(t))=1Γ(α)∫0t(ℏ3(s,Wn−1)−ℏ3(s,Wn−2))ds.(16)

Therefore, we have

Sn(t)=∑i=0nψi(t),(17)

In(t)=∑i=0nφi(t),(18)

Wn(t)=∑i=0nϑi(t).(19)

Let,

‖ψn(t)‖=‖Sn(t)−Sn−1(t)‖,‖ψn(t)‖=1Γ(α)∫0t(ℏ1(s,Sn−1)−ℏ1(s,Sn−2))ds,‖ψn(t)‖=ξ1Γ(α)∫0t‖Sn(t)−Sn−1(t)‖ds,‖ψn(t)‖=ξ1Γ(α)∫0tψn−1(t)ds.(20)

The remaining equations in the system (15) and (16) might be solved using the same method to get

‖φn(t)‖=ξ2Γ(α)∫0tφn−1(t)ds,(21)

‖ϑn(t)‖=ξ3Γ(α)∫0tϑn−1(t)ds.(22)

As required. □

Theorem 2: Show that (i) The system (14)–(16) has a specified uniform function. (ii) If there is a t∗>1 such that ξ1Γ(α)<1. Assuming σiα=0;i=1,2,3, the model system has at least one solution if ξ1Γ(α)<1 for i=1,2,3.

Proof: Since each kernel ℏ1 for i=1,2,3, satisfies the Lipschitz condition and the functions S(t),I(t), and W(t) are bounded, the following relations may be determined.

‖ψn(t)‖≤‖S(0)‖‖ξ1Γ(α)(t)‖n,(23)

‖ϑn(t)‖≤‖I(0)‖‖ξ2Γ(α)(t)‖n,(24)

‖ψn(t)‖≤‖W(0)‖‖ξ3Γ(α)(t)‖n.(25)

In the system (23)–(25) figure, it is demonstrated that the function given in (17)–(19) exists and is consistent. It is necessary to show that S(t),I(t), and W(t) converge to the system of solutions of (1)–(3) to prove (ii). To do this, we define the remaining terms after n changes as An(t),Bn(t), and Cn(t). Thus,

S(t)−S(0)=Sn(t)−An(t),(26)

I(t)−I(0)=In(t)−Bn(t),(27)

W(t)−W(0)=Wn(t)−Cn(t).(28)

By using the ξ1 Lipschitz condition and the triangle inequality, we obtain to determine that

‖An(t)‖=1Γ(α)∫0t(ℏ1(s,Sn−1)−ℏ1(s,Sn−2))ds,‖An(t)‖≤ξ1Γ(α)‖Sn(t)−Sn−1(t)‖.(29)

When the process in (29) is repeated, we obtain

‖An(t)‖≤‖ξ1Γ(α)(t)‖n+1ℰ1.(30)

Next, at t∗, one acquires

‖An(t)‖≤‖ξ1Γ(α)(t∗)‖n+1ℰ1.(31)

Assuming n→∞ as the limit,

limn→∞‖An(t)‖≤limn→∞‖ξ1Γ(α)(t∗)‖n+1ℰ1.(32)

By applying the hypothesis ξ1Γ(α)(t∗)<1, we have from (32) yield.

limn→∞‖An(t)‖=0.(33)

By using the same process as for n → ∞, we get

‖Bn(t)‖→0,(34)

‖Cn(t)‖→0.(35)

As a result, there must be only one solution.

As desired. □

Theorem 3: If (1−ξ1Γ(α)(t))>0 for the assumption σiα=0;i=1,2,3, then the system (1)–(3) has a unique solution.

Proof: Consider that another collection of solutions to (1)–(3) is represented by the sets S1,I1, and W1.

‖S(t)−S1(t)‖=1Γ(α)∫0t(ℏ1(s,Sn−1)−ℏ1(s,Sn−2))ds,‖S(t)−S1(t)‖≤ξ1Γ(α)‖S(t)−S1(t)‖.(36)

If the terms in (36), are rearranged, one gets

(1−ξ1Γ(α)(t))‖S(t)−S1(t)‖≤0,(37)

By applying the hypothesis (1−ξ1Γ(α)(t))>0, we have from (37) yield.

‖S(t)−S1(t)‖=0.(38)

It follows from this because S(t)=S1(t). Applying the identical process to every solution for i=2,3, we arrive at

I(t)=I1(t),(39)

W(t)=W1(t).(40)

Hence proved. □

Theorem 4: For the initial conditions with assumption σiα=0;i=1,2,3, prove that the stochastic fractional delayed model (1)–(3) has a positive solution in R+3.

Proof: The feasible situation must be non-negative across the system to be considered under the initial condition. We obtain

Dtα0c[S(t)]|S=0=Λα>0, Dtα0c[I(t)]|I=0=βαSIe−(μα+γα)τ>0, Dtα0c[W(t)]|W=0=εαI>0,

hence, the stochastic fractional delayed model (1)–(3) has a positive solution, when the initial condition falls inside the feasible region. □

Theorem 5: The system (1)–(4) in the feasible region 𝒢={(S(t),I(t),W(t))∈R+3;0<N≤Λαμα,∀t≥0,τ<t}; (where N(t)=S(t)+I(t)+W(t)) at any given time t. The initial condition is bounded with the assumption that σiα=0;i=1,2,3.

Proof: The total sum of plant population can be written as

Dtα0cN(t)≤Λα−μαN(t).

After resolving the inequality above, we obtain

N(t)≤Λαμα+(N(0)−Λαμα)e−μαt.

Using Grown’s inequality

limt→∞SupN(t)≤Λαμα.(41)

Therefore, the epidemiologically feasible region for the propagation of cassava mosaic disease is provided by (41).

𝒢={(S(t),I(t),W(t))∈R+3;0<N≤Λαμα,∀t≥0,τ<t}.(42)

The stochastic fractional delayed model (1)–(3) is both positively invariant and realistic from an epidemiological point of view concerning the transmission of leukemia disease (42). Hence, the system (1)–(3) is bounded under the initial conditions. □

2.2 Model Equilibria and Reproduction Number

In this part, we examine the different states of the stochastic fractional delayed model (1)–(3) with the propagation of leukemia disease dynamics. It allows the analysis of the system behavior and how the system behaves at the free and present state of the leukemia equilibrium. It also provided information on how the fractional order, delays, and stochastic parameters behave and control the transmission of leukemia disease using different numerical techniques. Therefore,

Leukemia Free Equilibrium=LFE=C0=(S0,I0,W0)=(Λαμα,0,0),(43)

Leukemia Present Equilibrium=LPE=C∗=(S∗,I∗,W∗),(44)

S∗=(γα+εα)βαe−(μα+γα)τ, I∗=Λαβαe−(μα+γα)τ−μα(γα+εα)βαe−(μα+γα)τ(γα+εα), W∗=εαbαI∗

In epidemiology, the basic reproduction number is an important parameter. This suggests whether or not the disease is prevalent in the general population. If the value of it is less, then the disease is not spreading in the population, otherwise, the disease is present in the population. For the evaluation of the reproduction number using the Next-generation method. The largest eigenvalue or spectral radius of, at leukemia-free equilibrium (43) called reproduction number as follows:

R0=Λαβαe−(μα+γα)τμα(γα+εα).(45)

Theorem 6: Assuming σiα=0;i=1,2,3, the leukemia-free equilibrium (43) is locally asymptotically stable for α∈(0,1) if R0<1.

Proof: By linearizing the stochastic fractional delayed model (1)–(3) about (43) a 3×3 dimensional Jacobian matrix with negative real components and eigenvalues is obtained.

λ1=−μα, λ2=−bα, λ3=−(γα+εα)(1−R0).

As a result, the leukemia-free equilibrium of the provided stochastic fractional delayed model (1)–(3) is locally stable if R0<1. If R0>1, then (43) is unstable in the local sense. □

Theorem 7: Assuming σiα=0;i=1,2,3, the stochastic fractional delayed model (1)–(3) is globally asymptotically stable (GAS) at leukemia-free equilibrium, C0 if R0<1.

Proof: Define the Lyapunov function ℒ:𝒢→R defined as

ℒ(t)=[S−S0−S0logSS0]+I+W,Dtα0cℒ(t)=[S−S0S]DtαS+DtαI+DtαW,Dtα0cℒ(t)=[S−S0S](Λα−βαSIe−(μα+γα)τ−μαS)+(βαSIe−(μα+γα)τ−(γα+εα)I)+(εαI−bαW),Dtα0cℒ(t)≤−Λα(S−S0)2SS0−γαI(1−βαSe−(μα+γα)τγα)−bαW.

This implies that Dtα0cℒ≤0 if R0<1 and Dtα0cℒ=0 if S(t)=S0,I(t)=W(t)=0. Therefore, C0 is globally asymptotically stable. □

Theorem 8: Assuming σiα=0;i=1,2,3, the leukemia-present equilibrium (44) is locally asymptotically stable for α∈(0,1) if R0>1.

Proof: By linearizing the stochastic fractional delayed model (1)–(3) about (44) a 3×3 dimensional Jacobian matrix with negative real components and eigenvalues is obtained.

λ1=−bα,λ2+a1λ+a0=0,a1=μα+βαS∗e−(μα+γα)τ,a0=μα(γα+εα)(R0−1)+(γα+εα)(2μα+βαI∗e−(μα+γα)τ).

As a result, the leukemia-present equilibrium of the provided stochastic fractional delayed model (1)–(3) is locally stable if R0>1. If R0<1, then (1)–(3) is unstable in the local sense. □

Theorem 9: Assuming σiα=0;i=1,2,3,4, the stochastic fractional delayed model (1)–(4) is globally asymptotically stable (GAS) at leukemia present equilibrium, C∗ if R0>1.

Proof: Define the Lyapunov function Z:𝒢→R defined as

Z=k1(S−S∗−S∗ln⁡(SS∗))+k2(I−I∗−I∗ln⁡(II∗))+k3(W−W∗−W∗ln⁡(WW∗)).

Given positive constants ki(i=1,2,3), we can express the following equation:

Dtα0cZ=k1[S−S∗S]DtαS+k2[I−I∗I]DtαI+k3[W−W∗W]DtαW,

Dtα0cZ=−k1Λα(S−S∗)2SS∗−k2βαSIe−(μα+γα)τ(I−I∗)2II∗−k3εαI(W−W∗)2WW∗.

If we choose ki where (i=1,2,3)

Dtα0cZ=−Λα(F−F∗)2FF∗−βαSIe−(μα+γα)τ(I−I∗)2II∗−εαI(W−W∗)2WW∗.

Dtα0cZ≤0, for R0>1 and Dtα0cZ=0 if and only if S=S∗,I=I∗,W=W∗. Hence by Lasalle’s invariance principle (44) is globally asymptotical stable. □

3  Sensitivity Analysis

In this section, we examine the behavior of model parameters concerning reproduction number R0. Examine the transmission and spread of disease with the sensitive analysis of the model. Preliminary: The formalized sensitivity index of a variable ℮, that depends differentiable on a parameter ϱ:

ℰϱ℮=ϱ℮×∂℮∂ϱ.

In spatial terms, determine the sensitive indices of parameters concerning the reproduction number R0.

VΛα=ΛαR0×∂R0∂Λα=1>0, Vβα=βαR0×∂R0∂βα=1>, Vμα=μαR0×∂R0∂μα=−1μα<0,

Vγα=γαR0×∂R0∂γα=−1(γα+εα)<0, Vεα=εαR0×∂R0∂εα=−1(γα+εα)<0.

Tables 1 and 2 provide the values of the sensitivity indicators together with the uncertainty indicators. All the parameters (Λα,βα) have positive sensitivity values. This means that the parameters have raised the value of R0 it becomes possible to make it high, hence making the system out of control. These parameters could represent transmission rates or factors that enhance the spread of the condition. The parameters (μα,γα,εα) associated with negative sensitivity values suggest that a higher value of these parameters results in a lower R0. These parameters probably measure natural mortality rates of susceptible blood cells, infected cells, and infected cells recovered due to encounters with immune cells which decrease the persistence or spread of the condition. For instance, if μα=−2,γα=−1.996,εα=−1.9996, the model also illustrates that such parameters’ increase would imply a remarkably reduced reproduction number, which underlines the importance of recovery or intervention mechanisms. Fig. 2 shows how different parameters affect R0. When increased, the parameters with negative sensitivity like (μα,γα,εα) decrease R0, and therefore the spread of the disease can be minimized if either recovery improves, or transmission is less. Conversely, those factors that lead to higher values of sensitivity parameters (Λα,andβα) are indicative of factors that raise the value of R0.

Figure 2: Sensitivity indices of reproduction number (R0)

4  Stochastic Fractional Delayed GL-NSFD Method

This section presents a numerical method for the stochastic fractional delayed model (1)–(3) that underlies it. Another instance of the stochastic fractional delayed system is provided.

Dtα0c[S]|t=tn=Λα−βαSn+1Ine−(μα+γα)τ−μαSn+1+σ1αSndB(t),(46)

Dtα0c[I]|t=tn=βαSn+1Ine−(μα+γα)τ−(γα+εα)In+1+σ2αIndB(t),(47)

Dtα0c[W]|t=tn=εαIn+1−bαWn+1+σ3αWndB(t).(48)

First, the Grunwald-Letnikove approach, or GL is explained:

Dtα0cυ(t)|t=tn=1(𝒦(Δt))α(𝒦n+1−∑i=1n+1νi𝒦n+1−i−ρn+1𝒦0),(49)

here,

νi=(−1)i−1(αi)ν1=α,

ρi=i−αΓ(1−α);i=1,2,3,…,n+1.

Now, the outcome that follows helps verify some other hypotheses.

NSFD rules are added to the GL approach, making the discrete model for susceptible cells as

1(𝒦(h))α(Sn+1−∑i=1n+1νiSn+1−i−ρn+1S0)=Λα−βαSn+1Ine−(μα+γα)τ−μαSn+1+σ1αSnΔBn,(50)

Sn+1=∑i=1n+1νiSn+1−i−ρn+1S0+(𝒦(h))α(Λα+σ1αSnΔBn)1+(𝒦(h))α(βαIne−(μα+γα)τ+μα).(51)

Additionally, the latent and breaking out susceptible GL-NSFD scheme as

In+1=∑i=1n+1νiIn+1−i+ρn+1I0+(𝒦(h))α(βαSn+1Ine−(μα+γα)τ+σ2αInΔBn)1+(𝒦(h))α(γα+εα),(52)

Wn+1=∑i=1n+1νiWn+1−i+ρn+1W0+(𝒦(h))α(εαIn+1+σ3αWnΔBn)1+(𝒦(h))α(bα).(53)

5  Positivity and Boundedness of Stochastic Fractional Delayed GL-NSFD

The positivity and boundedness of the solution for the systems (1)–(3) are confirmed by the following theorem.

Theorem 10: Suppose that S0≥0,I0≥0,W0≥0,Λα≥0,βα≥0,μα≥0,γα≥0,εα≥0,bα≥0 then Sn≥0,In≥0,Wn≥0, for all n=1,2,3,… with assumption of ΔBn=0.

Proof: For this, by using the Induction method we get

For n=0,

S1=ν1S0+ρ1S0+(𝒦(h))α(Λα)1+(𝒦(h))α(βαI0e−(μα+γα)τ+μα)≥0,

I1=ν1I0+ρ1I0+(𝒦(h))α(βαS1I0e−(μα+γα)τ)1+(𝒦(h))α(γα+εα)≥0,

W1=ν1W0+ρ1W0+(𝒦(h))α(εαI1)1+(𝒦(h))α(bα)≥0.

For n=1,

S2=ν1S1+ν2S0+ρ2S0+(𝒦(h))α(Λα)1+(𝒦(h))α(βαI1e−(μα+γα)τ+μα)≥0,

I2=ν1I1+ν2I0+ρ2I0+(𝒦(h))α(βαS2I1e−(μα+γα)τ)1+(𝒦(h))α(γα+εα)≥0,

W2=ν1W1+ν2W0+ρ2W0+(𝒦(h))α(εαI2)1+(𝒦(h))α(bα)≥0.

For n=2,

S3=ν1S2+ν2S1+ν3S0+ρ3S0+(𝒦(h))α(Λα)1+(𝒦(h))α(βαI2e−(μα+γα)τ+μα)≥0,

I3=ν1I2+ν2I1+ν3I0+ρ3I0+(𝒦(h))α(βαS3I2e−(μα+γα)τ)1+(𝒦(h))α(γα+εα)≥0,