Numerical Treatments for a Crossover Cholera Mathematical Model Combining Different Fractional Derivatives Based on Nonsingular and Singular Kernels

1  Introduction

Cholera continues to pose a major public health challenge, particularly in endemic regions where sanitation and access to clean water are limited. According to Ali et al. [1], cholera remains a significant burden in many parts of the world, with periodic outbreaks contributing to considerable morbidity and mortality. Mathematical modeling has become an essential tool in understanding and predicting the spread of infectious diseases like cholera, offering insights that inform public health interventions and optimize resource allocation. Huppert and Katriel [2] emphasize the importance of mathematical frameworks in epidemiology, highlighting their capacity to simulate complex transmission dynamics and forecast outbreak patterns.

Recent advancements in mathematical modeling have incorporated fractional calculus to better represent the memory effects and hereditary properties inherent in disease transmission. Fractional-order models provide a more flexible and accurate description of epidemic processes compared to classical integer-order models, capturing the influence of past states on current disease progression. Lemos-Paião et al. [3] developed a cholera model with optimal control strategies, demonstrating how fractional dynamics can significantly enhance our understanding of disease spread and control measures. Similarly, the work by Lemos-Paião et al. [4] introduced an epidemic model for cholera that utilizes optimal control to manage treatment strategies, further illustrating the utility of fractional models in epidemiological applications.

Baleanu et al. [5] expanded on these approaches by proposing a generalized fractional model for real cholera outbreaks, incorporating stochastic components to account for the unpredictable nature of disease dynamics. Their model successfully captured the complexities observed in actual outbreaks, validating the effectiveness of fractional and stochastic approaches in improving predictive accuracy.

Fractional calculus has gained significant attention in recent years due to its ability to describe memory effects and hereditary properties in various physical and biological systems. Unlike classical integer-order derivatives, fractional derivatives provide a more generalized framework for modeling complex dynamical processes. In particular, special functions play a crucial role in solving fractional differential equations, offering analytical and numerical tools for applications in engineering, physics, and mathematical biology. The work presented in [6] explores the fundamental role of special functions in fractional models, highlighting their applicability in diverse scientific domains. Incorporating these special functions allows for more accurate representations of anomalous diffusion, viscoelastic materials, and population dynamics. This perspective aligns with our study, where fractional-order derivatives are employed to enhance the mathematical modeling of disease transmission and epidemiology, providing deeper insights into system behavior.

Building upon these foundational studies, this paper introduces a novel crossover cholera model that integrates multiple fractional derivatives with both singular and nonsingular kernels, bridging deterministic and stochastic processes to better reflect real-world conditions. The proposed model leverages three distinct types of fractional derivatives:

•   Caputo-Fabrizio fractional derivative (nonsingular kernel) [7], known for its smooth memory effects and exponential decay properties.

•   Caputo proportional constant fractional derivative (singular kernel) [8], capturing abrupt changes and reflecting more localized influences in the disease dynamics.

•   Atangana-Baleanu fractional derivative (nonsingular kernel) [9], offering a more generalized representation of memory effects through Mittag-Leffler functions.

These derivatives are applied across four distinct time intervals, reflecting the various stages of cholera progression, from initial infection to outbreak resolution. By incorporating stochastic differential equations (SDEs) alongside fractional derivatives, the model accounts for the inherent randomness and uncertainty present in disease transmission, environmental factors, and population mobility.

This crossover model not only refines the predictive accuracy of cholera outbreaks but also extends the applicability of fractional calculus to broader epidemiological contexts. By addressing the limitations of traditional models and incorporating stochastic elements, the proposed framework contributes to the ongoing effort to mitigate the global impact of cholera through improved forecasting and intervention strategies.

These derivatives are applied across four distinct time intervals, reflecting the various stages of cholera progression, from initial infection to outbreak resolution. By incorporating stochastic differential equations (SDEs) alongside fractional derivatives, the model accounts for the inherent randomness and uncertainty present in disease transmission, environmental factors, and population mobility.

To ensure accurate numerical approximations, the model employs a combination of advanced numerical techniques:

•   The approximation of Caputo proportional constant fractional derivative with Grünwald-Letnikov nonstandard finite difference method (CPC-GLNFDM) is applied to systems with singular kernels to achieve stable and accurate approximations [10].

•   The Toufik-Atangana method was utilized for models involving nonsingular Mittag-Leffler kernels, ensuring convergence and minimizing computational error [11].

•   Caputo-Fabrizio integral approximation and Lagrange polynomial interpolation applied to nonsingular exponential decay kernels, capturing the long-term behavior of the disease [11].

•   The Milstein method is used for stochastic differential equations, providing precise approximations for random fluctuations in transmission dynamics [12].

The authors confirm that this study is original and has not been discussed before, and it is presented in this paper for the first time.

The structure of this study is as follows: Section 2 introduces key definitions of fractional derivatives, provides a brief overview of stochastic differential equations (SDEs), and outlines the numerical methods used to solve them. Section 3 presents the foundational model along with its mathematical analysis. Section 4 introduces a novel crossover model to describe the dynamics of cholera spread. This model integrates three fractional derivatives with both singular and nonsingular kernels, as well as SDEs across four distinct time intervals. The fractional derivatives employed include the Caputo-Fabrizio derivative (nonsingular kernel), the Caputo proportional constant derivative (singular kernel), and the Atangana-Baleanu derivative (nonsingular kernel). Section 5 details the application of four numerical methods to solve the proposed system and analyze its stability. Section 6 presents numerical simulations based on the proposed model, illustrating its performance and accuracy. Section 7 concludes the study by summarizing the key findings and highlighting the primary contributions of the research.

2  Preliminaries and Notations

In this section, we go over some of the most important fractional calculus definitions that are used in the remainder of this research.

Definition 1. Let Ω=[a,b], the value of α∈C, ℜ(α)>0,  −∞<a<b<∞, reference [13] defines the Riemann-Liouville’s derivatives of order α on the left and right sides, respectively, for a continuous function f(t).

Dtαaf(t)=1Γ(−α+n)(ddt)n∫atf(s)(−s+t)α−n+1ds,  t>a,(1)

Dbαtf(t)=1Γ(−α+n)(−ddt)n∫tbf(s)(−t+s)α−n+1ds,  t<b,(2)

where n=1+[ℜ(α)].

Definition 2. Let Ω=[a,b], the value of α∈C, ℜ(α)>0, −∞<a<b<∞. For a function that is continuous f(t), reference [13] defines the Riemann-Liouville’s integrals of order α on the left and right sides, respectively.

Itαaf(t)=1Γ(α)[∫atf(s)(−s+t)−1+αds],    t>a,  (3)

Ibαtf(t)=1Γ(α)[∫tbf(s)(−s+t)−1+αds],     t<b, (4)

where 0<α<1.

Definition 3. Let Ω=[a,b], −∞<a<b<∞, α∈C then reference [13] defines the left-side and right-side Caputo’s derivatives of order α for a function f(t), respectively, f∈ACn[a,b], i.e., f(k) is absolutely continuous for all k≤n, and f(n) is integrable and ACn[a,b] refers to the space of absolutely continuous functions with n derivatives on the interval [a,b].

(CDa+αf)(t)=(aCDtαf)(t)=1Γ(−α+n)∫atf(n)(s)(−ξ+t)α−n+1ds,  t>a,   (5)

(CDb−αf)(x)=(tCDbαf)(t)=(−1)nΓ(α+n)∫tbf(n)(s)(−t+s)α−n+1ds,   t<b.   (6)

Definition 4. [8] For α∈(0,1) and 𝒵(t)∈𝒲21(0,1), the Caputo-Febrizio fractional derivative of 𝒵(t) of order α is given by:

0CFDtα𝒵(t)=𝒩(α)1−α∫0tddτ𝒵(κ)exp[−α1−α(−κ+t)]dκ,(7)

the normalizing function 𝒩(α) ensures that 𝒩(0) =1 and 𝒩(1) = 1. The function 𝒵(t)’s fractional integral with exponential decay kernel is provided by:

0CFDtα𝒵(t)=1−α𝒩(α)𝒵(t)+α𝒩(α)∫0t𝒵(κ)dκ,(8)

where 𝒲21(0,1) is the Sobolev space of functions with square-integrable first derivatives, n=1+[ℜ(α)], ℜ(α)∉N0 and N0 is the set of non-negative integers.

Definition 5. The Caputo proportional fractional hybrid operator (CP) is generally defined as [8]:

Dtα0 CP𝒵(t)=(∫0t(𝒢1(s,α)𝒵(s)+𝒢0(s,α)𝒵′(s))(−s+t)−αds)1Γ(−α+1),=(𝒢1(t,α)𝒵(t)+𝒢0(t,α)𝒵′(t))(t−αΓ(−α+1)),(9)

where 𝒢0(α,t)=αt(−α+1), 𝒢1(α,t)=(−α+1)tα, 0<α<1.

Definition 6.The Caputo proportional constant fractional hybrid operator (CPC) can be defined as [8]:

Dtα0   CPC𝒵(t)=(∫0t(−s+t)−α(𝒢1(α)𝒵(s)+𝒢0(α)𝒵′(s))ds)1Γ(−α+1)=𝒢1(α)0RLIt1−α𝒵(t)+𝒢0(α)0CDtα𝒵(t),(10)

where 𝒢0(α)=α𝒯  (−α+1), 𝒢1(α)=(−α+1)𝒯α, 𝒯 is a constant.

Definition 7. According to Caputo, the Atangana Baleanu fractional derivative of the function 𝒵(t) is determined by the formula that follows [8]:

0ABCDtα𝒵(t)=𝒜ℬ(α)−α+1∫at𝒵˙(κ)Eα[−α−α+1(−κ+t)α]dκ,(11)

and the Riemann-Liouville Atangana Baleanu fractional derivative of the function 𝒵(t) is determined by the formula that follows:

0ABRDtα𝒵(t)=𝒜ℬ(α)−α+1ddτ∫at𝒵(κ)Eα[−α−α+1(−κ+t)α]dκ.(12)

Here, 𝒜ℬ(α)=−α+1+αΓ(α), 𝒵(t)∈𝒲21(0,1), and α∈(0,1). The fractional integral is represented by the following equation, which utilizes the Mittag-Leffler kernel:

0ABCItα𝒵(t)=−α+1𝒜ℬ(α)𝒵(κ)+α𝒜ℬ(α)Γ(α)+∫0t𝒵(κ)[(−κ+t)−1+α]dκ,(13)

where n=1+[ℜ(α)], ℜ(α)∉𝒩0.

Milstein method

SDE can be approximated numerically using the Milstein method in mathematics. It bears the name of Milstein, who published it in 1974 [12]. To solve the following (SDE):

d𝒳(t)=dW(t)g(𝒳(t))+f(𝒳(t))dt,𝒳(0)=𝒳0.(14)

Using the Milstein method like this:

𝒳1+n=𝒳n+f(𝒳n)Δt+g(𝒳n)ΔW+0.5g(𝒳n)g′(𝒳n)((ΔWn)2−Δt),(15)

where dΔWn describes the standard Brownian motion.

3  Model Setup

In this section, we introduce the basic cholera disease model as outlined in [3]. This model incorporates the following fundamental populations: quarantined individuals (𝒬), infectious individuals (ℐ), susceptible individuals (𝒮), and recovered individuals (ℛ). Additionally, the model considers the population of free bacteria in the environment (ℬ). This inclusion is critical, as cholera transmission occurs when a healthy individual consumes contaminated water containing these bacteria, which are subsequently removed from the aquatic environment.

The parameters of the model are defined in Table 1, and the mathematical representation of the model is presented as follows:

𝒮˙(t)=−βℬ(t)𝒮(t)ν+ℬ(t)+ωℛ(t)−μ𝒮(t)+Λ,ℐ˙(t)=βℬ(t)𝒮(t)ν+ℬ(t)−(δ+μ+λ1)ℐ(t),𝒬˙(t)=−(ε+μ+λ2)𝒬(t)+δℐ(t),ℛ˙(t)=−(μ+ω)ℛ(t)+ε𝒬(t),ℬ˙(t)=−ρℬ(t)𝒮(t)ν+ℬ(t)−dℬ(t)+ηℐ(t).(16)

3.1 Analysis of the Model

In this study, we assume non-negative initial conditions for system (16):

𝒬(0)=𝒬0≥0,ℐ(0)=ℐ0≥0,𝒮(0)=𝒮0≥0,ℛ(0)=ℛ0≥0,ℬ(0)=ℬ0≥0.

3.2 Solutions’ Boundedness and Positivity

In the following, the domain Ω described in the equations provides the biological constraints and the feasible region for the model variables. Here’s the interpretation of each component:

1. Human population domain Ωℋ

Ωℋ={(𝒬,ℐ,𝒮,ℛ)∈R04:𝒬(t)+ℐ(t)+𝒮(t)+ℛ(t)≤Λμ}.(17)

Biological Meaning: This constraint ensures that the total population remains within a biologically feasible limit, accounting for recruitment and death rates.

2. Bacterial population domain Ωℬ

Ωℬ={ℬ∈R0+:ℬ(t)≤Λημd}.(18)

Biological Meaning: This constraint reflects the balance between the contribution of bacteria to the environment (from infected individuals) and its natural decay. It ensures that the bacterial population remains within realistic ecological limits.

3. Combined domain Ω

Ω=Ωℋ×Ωℬ.

Biological Meaning: This domain ensures that the model operates within biologically meaningful constraints for both human and bacterial populations.

Lemma 1. The solutions of system (16) maintain non-negative values for all t≥0, provided that the initial conditions (3.1) are non-negative and belong to the non-negative orthant of five-dimensional Euclidean space, ℛ0+5.

Proof: We have

{𝒮(t)˙|ξ(𝒮)=ωℛ(t)+Λ>0,ℐ(t)˙|ξ(ℐ)=𝒮(t)ℬ(t)βν+ℬ(t)>0,𝒬(t)˙|ξ(𝒬)=δℐ(t)>0,ℛ(t)˙|ξ(ℛ)=ϵ𝒬(t)>0,ℬ(t)˙|ξ(ℬ)=ηℐ(t)>0,

where ξ(ν)={ν(t)=0,ν={𝒬,ℐ,𝒮,ℛ,ℬ}}. Therefore, any solution of system (16) is such that (𝒬,ℐ,𝒮,ℛ,ℬ)∈R0+5 for all t≥0. ◻

3.3 Analysis of Stability and Equilibrium Points

3.3.1 Equilibrium Points

Equilibrium points represent steady-state conditions where the system remains unchanged over time. In the context of the cholera model, equilibrium points might represent scenarios where the disease has either been eradicated (disease-free equilibrium) or has reached a stable endemic state.

In the context of cholera models, the disease-free equilibrium (DFE) represents a state where the disease has been eradicated from the population.

It is produced when we set the L.H.S. system (16) to zero and put ℐ=0. We’ve

E0=(𝒬0,ℐ0,𝒮0,ℛ0,ℬ0)=(0,0,Λμ,0,0).(19)

The endemic equilibrium corresponds to a non-trivial solution where ℐ≠0 (i.e., there is a persistent infection in the population). We denote for the endemic equilibrium by E∗ and given as:

E∗=(𝒬∗,ℐ∗,𝒮∗,ℛ∗,ℬ∗).

3.3.2 Stability Analysis

First, we compute the fundamental reproduction factor R0 [14].

Proposition 1: Model (16)’s fundamental reproduction factor is:

R0=βΛη(δ+λ1+μ)(μνd+ρΛ).

Proof: Let ℱi(t) represent the rate at which new infections appear in the compartment corresponding to index i [14]. Define 𝒱i+(t) as the rate at which individuals enter the compartment associated with index i through all other means, while 𝒱i−(t) denotes the rate at which individuals leave the compartment with index i. Based on this, the matrices ℱ(t), 𝒱+(t), and 𝒱−(t), which are associated with model (16), are given by:

ℱ(t)=[0βℬ(t)𝒮(t)ν+ℬ(t)000],𝒱+(t)=[Λ+ωℛ(t)0δℐ(t)ϵ𝒬(t)ηℐ(t)],𝒱−(t)=[βℬ(t)𝒮(t)ν+ℬ(t)+μ𝒮(t)b1ℐ(t)b2𝒬(t)b3ℛ(t)dℬ(t)]

where:

b1=δ+λ1+μ,b2=ϵ+λ2+μ,b3=μ+ω.

Thus, since:

𝒱(t)=𝒱−(t)−𝒱+(t),

we have

[𝒮′(t)ℐ′(t)𝒬′(t)ℛ′(t)ℬ′(t)]=ℱ(t)−𝒱(t).

The Jacobian matrix of ℱ(t) is:

F=[00000ℬ(t)βν+ℬ(t)000𝒮(t)κβ(ν+ℬ(t))2000000000000000].

The Jacobian matrix of 𝒱(t) is:

V=[ℬ(t)βν+ℬ(t)+μ00−ω𝒮(t)κβ(ν+ℬ(t))20b10000−δb20000−ϵb300−η00c].

Using the disease-free equilibrium E0 (19), we get the matrices F0 and V0 as:

F0=[000000000Λβνμ000000000000000],V0=[μ00−ωΛβνμ0b10000−δb20000−ϵb300−η00c].

The model (16)’s fundamental reproduction factor is:

R0=ρ(F0V0−1)=βΛη(μνd+ρΛ)(δ+λ1+μ).

So, the proof is finished. ◻

Now we prove the disease-free equilibrium E0’s local stability.

Theorem 2. The model (16)’s disease-free equilibrium (DFE) E0 is:

1.   Locally and asymptotically stable, provided that

Λβη<(μνd+ρΛ)(λ1+δ+μ).

2.   If Λβη>(μνd+ρΛ)(λ1+δ+μ), the DFE is unstable.

Moreover, a critical case arises if:

βΛη=(μνd+ρΛ)(δ+λ1+μ).

Proof: The linearized system of model (16)’s characteristic polynomial is provided by:

p~(χ)=det(ℱ0−𝒱0−ℐ5χ)

To calculate the polynomial p’s roots, we possess that

|−χ−μ00ω−Λβκμ0−χ−b100Λβνμ0δ−χ−b20000ϵ−χ−b300η00−χ−d|=0.

That is:

χ=−μ,χ=−b2,χ=−b3,p~(χ)=χ2+(d+b1)χ+db1−Λβηνμ=0. ◻

According to Routh’s standard [15], the roots of the polynomial p~(χ) have a negative real part if all of its coefficients have the same signal. As a result, the local asymptotic stability of the (DFE) E0 is established, p~(χ)’s coefficients are:

p~1=1>0,p~2=d+b1>0,p~3=db1−Λβηνμ.

Thus, if

db1−Λβηνμ>0⇔Λβη<(μνd+ρΛ)(μ+δ+λ1).

E0 be locally asymptotically stable.

else if db1−Λβηνμ<0⇔Λβη>(μνd+ρΛ)(μ+δ+λ1), the DFE is unstable.

else db1=Ληνμ⇔βΛη=(μνd+ρΛ)(μ+δ+λ1), yields a crucial case.

Proposition 2: The equations b1=μ+δ+λ1, b2=ϵ+μ+λ2, and b3=μ+ω should be considered. δ,λ∗,ϵ, ω>0 is assumed. Model (16) possesses the equilibrium of endemic if R0>1.

E∗=(𝒬∗,ℐ∗,𝒮∗,ℛ∗,ℬ∗)=(δΛb3λ∗𝒟,b2Λb3λ∗𝒟,b1Λb2b3𝒟,δΛϵλ∗𝒟,(ηβ−b1ρ)b2Λb3λ∗𝒟βc),(20)

where

𝒟=b2b1b3(μ+λ∗)−ϵδωλ∗,λ∗=ℬ∗βℬ∗+ν,

which is possible if

ηβ>b1ρ.(21)

Proof: The rate of transmission must be entirely positive in order for this equilibrium to be possible:

βℬ∗(t)ℬ∗(t)+ν>0.

The fourth, third and second equilibrium equations of (16) are solved successively, and we discover

𝒮∗=b1λ∗ℐ∗,ℐ∗=b2δ𝒬∗,𝒬∗=b3ϵℛ∗.

Next, we get

𝒮∗=ℛ∗(b2b1b3λ∗ϵδ).

The final value of 𝒮∗ that was assessed is then substituted into the initial equilibrium equation to get

ϵδλ∗Λ−ℛ∗𝒟=0.

It yields the fourth element of (20) as well as the first three through back substitution. Lastly, the practical value of ℬ∗, which needs to be nonnegative, is given by the fifth equation, which yields (21). Currently, from

λ∗=βℬ∗(ν+ℬ∗).

After rearranging and replacing the value of ℬ∗, we get

λ∗{[(βη−ρb1)Λ+νβcb1]b3b2−βνcδωϵ}=[Λβη−(νμc)b1+Λρ]b2βb3.

Consequently,

λ∗=b2b1b3β(−1+R0)(Λρ+cνμ)b2b3(ηβ−b1ρ)Λ+cβν[b2b1b3−ϵδω].

Considering (21), and the reality that

b2b1b3−ϵδω=(λ1+δ+μ)(λ2+ϵ+μ)(μ+ω)−ϵδω>0.

Since all of the other coefficients and λ1,λ2≥0 are positive, the value of λ∗ above is positive if and only if R0>1. An endemic equilibrium (20) exists in the model (16) in this situation. This brings the proof to an end.

4  The Piecewise Mathematical Model

The piecewise cholera mathematical model presented in this study enhances conventional epidemic modeling by employing sophisticated, multi-stage fractional differential equations to accurately capture the dynamics of cholera transmission. The model divides the progression of the disease into four distinct periods, each representing a critical stage in the spread of the infection. To reflect the unique characteristics of disease propagation at each stage, a different fractional derivative operator is applied within each interval.

In the initial phase 0<t≤T1, the model (16) utilizes the Caputo-Fabrizio fractional derivative, which provides a smooth representation of early-stage dynamics. This derivative effectively characterizes the rapid initial spread of the disease, with minimal memory effects influencing transmission.

During the second phase T1<t≤T2, the model (16) incorporates a proportional constant hybrid Caputo operator. This stage reflects the cumulative impact of infection, accounting for evolving interactions and immune responses as the disease progresses.

In the advanced stages of cholera, when T2<t≤T3, the model described in system (16) incorporates the Atangana-Baleanu fractional derivative. This operator is particularly effective in capturing the complex and irregular transmission patterns that emerge in the later stages of an outbreak. As the disease progresses, its spread becomes more unpredictable, making the Atangana-Baleanu derivative a suitable choice for modeling these dynamics.

The stochastic cholera model for T3<t≤Tf offers a more realistic representation of disease transmission by incorporating environmental and epidemiological randomness. Unlike deterministic models, it accounts for unpredictable fluctuations in disease dynamics, leading to improved epidemic forecasting and more effective strategies for adaptive public health interventions.

By adopting this piecewise differential approach, the model offers a more realistic and comprehensive depiction of epidemic stages, significantly improving predictive insights into the progression and potential containment of cholera outbreaks.

The resulting system is expressed as follows:

{1ϑ1−αDtα0 CF𝒮(t)=−βℬ(t)𝒮(t)ν+ℬ(t)+ωℛ(t)−μ𝒮(t)+∇,1ϑ1−αDtα0 CFℐ(t)=βB(t)𝒮(t)ν+ℬ(t)−(δ+λ1+μ)ℐ(t),1ϑ1−αDtα0 CF𝒬(t)=−(λ2+ϵ+μ)𝒬(t)+δℐ(t),    0<t≤T1,1ϑ1−αDtα0 CFℛ(t)=−(μ+ω)ℛ(t)+ϵ𝒬(t),1ϑ1−αDtα0 CFℬ(t)=−dℬ(t)−ρℬ(t)𝒮(t)ν+ℬ(t)+ηℐ(t),(22)

with initial conditions

𝒬0≥0, ℐ0≥0,𝒮0≥0,ℛ0≥0,ℬ0≥0.(23)

{1ϑ1−αDtα0 CPC𝒮(t)=−βℬ(t)𝒮(t)ν+ℬ(t)+ωℛ(t)−μ𝒮(t)+∇,1ϑ1−αDtα0 CPCℐ(t)=βℬ(t)𝒮(t)ν+ℬ(t)−(δ+λ1+μ)ℐ(t),1ϑ1−αDtα0 CPC𝒬(t)=−(λ2+ϵ+μ)𝒬(t)+δℐ(t),T1<t≤T2,1ϑ1−αDtα0 CPCℛ(t)=−(μ+ω)ℛ(t)+ϵ𝒬(t),1ϑ1−αDtα0 CPCℬ(t)=−dℬ(t)−ρℬ(t)𝒮(t)ν+ℬ(t)+ηℐ(t),(24)

with initial conditions

𝒬(t1)≥0, ℐ(t1)≥0,𝒮(t1)≥0,ℛ(t1)≥0,ℬ(t1)≥0.(25)

{1ϑ1−αDtα0 ABC𝒮(t)=−βℬ(t)𝒮(t)ν+ℬ(t)+ωℛ(t)−μ𝒮(t)+∇,1ϑ1−αDtα0 ABCℐ(t)=βℬ(t)𝒮(t)ν+ℬ(t)−(δ+λ1+μ)ℐ(t),1ϑ1−αDtα0 ABC𝒬(t)=−(λ2+ϵ+μ)𝒬(t)+δℐ(t),    T2<t≤T3,1ϑ1−αDtα0 ABCℛ(t)=−(μ+ω)ℛ(t)+ϵ𝒬(t),1ϑ1−αDtα0 ABCℬ(t)=−dℬ(t)−ρℬ(t)𝒮(t)ν+ℬ(t)+ηℐ(t),(26)

with

𝒬(t2)=𝒮2≥0, ℐ(t2)=ℐ2≥0,𝒮(t2)=𝒬2≥0,ℛ(t2)=ℛ2≥0,ℬ(t2)=ℬ2≥0.(27)

{d𝒮(t)=(−βℬ(t)𝒮(t)ν+B(t)+ωR(t)−μ𝒮(t)+∇)dt+σ1𝒮(t)dW1,dtℐ(t)=(βℬ(t)𝒮(t)ν+ℬ(t)−(δ+λ1+μ)ℐ(t))dt+σ2ℐ(t)dW2,d𝒬(t)=(−(λ2+ϵ+μ)𝒬(t)+δℐ(t))dt+σ3𝒬(t)dW3,    T3<t≤Tf,dℛ(t)=(−(ω+μ)ℛ(t)+ϵ𝒬(t))dt+σ4ℛ(t)dW4,dℬ(t)=(−dℬ(t)−ρℬ(t)𝒮(t)ν+ℬ(t)+ηℐ(t))dt+σ5ℬ(t)dW5,(28)

where σi,  i=1,2,…,5, represents the real constants of the stochastic environment intensity, and dWi describes the standard Brownian motion. with

𝒬(t3)=𝒮3≥0, ℐ(t3)=ℐ3≥0,𝒮(t3)=𝒬3≥0,ℛ(t3)=ℛ3≥0,  ℬ(t3)=ℬ3≥0.(29)

5  Numerical Approach to Solving the Proposed Model

5.1 Numerical Scheme with Exponential Decay Kernel

Consider the hybrid fractional order derivatives equation that follows:

Dtα0 CF𝒵(t)=𝒫(t,𝒵(t)),    0<α≤1,𝒵(0)=𝒵0.

Rewriting the above equation as [11] is possible based on Caputo-Fabrizio’s definition integral (8):

𝒵(t)−𝒵(0)=−α+1𝒩(α)𝒫(t,𝒵(t)) +α𝒩(α)∫0t𝒫(τ,𝒵(κ))dκ,(30)

We write (30) at t1+k=(1+k)Δt,

𝒵(t1+k)−𝒵(0)=−α+1𝒩(α)𝒫(tk,𝒵(tk)) +α𝒩(α)∫0t1+k𝒫(tκ,𝒵(tκ))dκ,

and at tk=kΔt:

𝒵(tk)−𝒵(0)=−α+1𝒩(α)𝒫(t−1+k,𝒵(t−1+k)) +α𝒩(α)∫0tk𝒫(tκ,𝒵(tκ))dκ,

where the normalization function 𝒩(α) ensures that 𝒩(0)=1 and 𝒩(1)=1.

Using these equations’ various approaches, the following can be written [11]:

𝒵(t1+k)−𝒵(tk)=−α+1𝒩(α)[𝒫(tk,𝒵(tk))−𝒫(t−1+k,𝒵(t−1+k))]+α𝒩(α)∫tkt1+k𝒫(tκ,𝒵(tκ))dκ.

The following can be obtained by entering its Lagrange polynomial into the equation above:

𝒵1+k−𝒵k=−α+1𝒩(α)[𝒫(tk,𝒵(tk))−𝒫(t−1+k,𝒵(t−1+k))]+α𝒩(α)∫tkt1+k[𝒫(tk,𝒵(tk))Δt(κ−t−1+k)−𝒫(t−1+k,𝒵(t−1+k))Δt(κ−tk)]dκ,

and we get the following:

𝒵1+k−𝒵k=−α+1𝒩(α)[𝒫(tk,𝒵(tk))−𝒫(t−1+k,Z(t−1+k))]+α𝒩(α)[𝒫(tk,𝒵(tk))Δt∫tkt1+k(κ−t−1+k)dκ−𝒫(t−1+k,𝒵(t−1+k))Δt∫tkt1+k(κ−tk)dκ].