EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 2, Article Number 5913
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

1

Impact of Singular and Non-Singular Kernels on2

Crossover Monkeypox Mathematical Model3

N. H. Sweilam1,∗, S. M. Al-Mekhlafi2,3, A. Ahmed4, D. G. Mohamed4,4

E. M. Abo-Eldahab4
5

1 Mathematics Department, Faculty of Science, Cairo University, Giza, Egypt6

2 Mathematics Department, Faculty of Education, Sana’a University, Yemen7

3 Department of Engineering Mathematics and Physics, Future University in Egypt, Egypt8

4 Mathematics Department, Faculty of Science, Helwan University, Cairo, Egypt9

10

Abstract. This study presents three crossover models describing monkeypox disease that includes
Caputo, Mittag-Leffler, and Caputo-Fabrizio definitions. To represent the monkeypox disease,
three models of variable-order fractional, fractal-fractional, and stochastic, as well as their piecewise
derivatives are provided at three different time periods. To approximate these models, we use the
nonstandard Grünwald-Letnikov finite difference method to approximate the deterministic model
with a singular kernel and a nonsingular Mittag-Leffler kernel to approximate the deterministic
model using the Toufik-Atangana method. Moreover, we use the approximation of the integral
Caputo-Fabrizio and Lagrange polynomial of two steps to approximate the deterministic model
with a nonsingular exponential decay kernel. We implemented the Milstein method to approximate
the stochastic differential equation. An analysis of the suggested model’s stability is conducted.
The effectiveness of the procedures was confirmed, and the theoretical results were supported
through numerical testing and comparisons with actual data.

2020 Mathematics Subject Classifications: 65L05, 26A33, 39A5011

Key Words and Phrases: Monkeypox disease, fractal-fractional derivative, Milstein method,12

Atangana-Baleanu Caputo operator, Caputo-Fabrizio operator13

14

1. Introduction15

The Monkeypox virus is the cause of this uncommon viral zoonotic disease. Monkey-16

pox, first identified in 1958, is a pox-like disease caused by the virus known as variola. It17

leads to fever, rash, and lymphadenopathy, with severe cases in children, pregnant women,18

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5913

Email addresses: nsweilam@sci.cu.edu.eg (N. H. Sweilam),
sih.almikhlafi@su.edu.ye (S. M. Al-Mekhlafi),
aya.a.mahmoud@science.helwan.edu.eg (A. Ahmed),
doaa.gamal@science.helwan.edu.eg (D. G. Mohamed),
emad.aboeldahab@yahoo.com (E. M. Abo-Eldahab)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 2 of 41

and immunocompromised individuals. Monkeypox cases have been reported in central and19

western African nations before the 2022 outbreak, with most cases linked to travel to en-20

demic countries [1]. Since May 2022, cases have been found in non-endemic and endemic21

countries, with most cases recorded in sexual health services. The biology, spread, threat,22

vaccinations, and available therapies of the virus are examined [2]. Researchers have sug-23

gested a number of theoretical and mathematical investigations to examine the dynamics of24

monkeypox [3], including analyzing its ecological niche, geographic distribution, projected25

burden, and transmission dynamics [4]. They have also proposed models for treatment26

and vaccination [5], a stochastic model [6], A new model that takes into account transfer27

from person to person [7], and a game-theoretic model for assessing vaccination plans [8].28

Fractional calculus is being increasingly applied in science and engineering, particu-29

larly in modeling deathly epidemics [6]. Recent works have focused on optimal control30

problems [9], two-strain epidemic models [10], and there are numerous models using vari-31

ous fractional derivatives [11]. This trend is gaining traction in various fields. Researchers32

have developed various models for various diseases, including Covid-19, vaccination effi-33

cacy, cholera outbreaks, human African trypanosomiasis epidemics, hepatitis B virus, and34

dengue fever and zika virus coinfection [12]-[13]. Also, there are important applications of35

Atangana-Baleanu, Caputo-Fabrizo, Caputo derivatives can be found in [14], [15], [16].36

These models use integer and fractional derivatives, fractional derivatives, and math-37

ematical techniques to analyze and control various diseases. Epidemiology modelling is38

intended to benefit from the use of generalised fractional mathematical models. Its main39

advantages include improved data fitting, multistage, capturing factuality and memory40

effects. A powerful tool for taking into consideration the system’s memory and inherited41

traits is provided by fractional epidemic models, as opposed to integer-order derivative for42

models which either disregard or prevent this from occurred. Furthermore, in terms of43

data fitting, the fractional variant provides one extra degree of flexibility over the integer44

model. A fractional system has various benefits, including the opportunity to select the45

fractional ordered value that best describes a model in practice, memory, and the most46

precise data fitting. The fractional derivatives models are further reinforced by inherited47

traits, making them more appropriate for representing real-world phenomena [17, 18].48

In order to better represent the complex behaviours of certain real-world situations,49

piecewise calculus has been developed in recent years. Various real-world phenomena,50

such as infectious illnesses, heat transport, fluid dynamics, and intricate advection issues,51

exhibit crossover behaviours [19]. Interestingly, recent studies indicate that piecewise52

formulations in differential equations give a more realistic depiction of these processes than53

standard fractional order or integer order methods. Remarkably, the notions of piecewise54

derivatives and short memory in fractional calculus are very similar. Considerable recent55

work has been done in this area, including citations to important publications like [20].56

Our motivation stems from the complex and heterogeneous nature of infectious disease57

dynamics, where traditional integer-order models may not fully capture the underlying58

memory effects, spatial heterogeneity, and randomness in disease transmission.59

- Variable-order fractional derivatives allow us to incorporate time-dependent changes in60

disease progression, reflecting evolving immunity, interventions, and behavioral adapta-61



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 3 of 41

tions.62

- Fractal-fractional derivatives help model spatial heterogeneity and self-similar transmis-63

sion patterns, which are particularly relevant in disease spread across different community64

structures.65

- Stochastic derivatives account for randomness in infection rates, environmental varia-66

tions, and uncertainties in disease reporting, making the model more realistic in capturing67

real-world epidemic fluctuations.68

These advanced modeling techniques offer a more flexible and accurate framework com-69

pared to traditional approaches, ultimately improving predictive power and decision-70

making in epidemiological studies.71

In order to create three new crossover systems for the Monkypox model, which is in-72

troduced for the first time in this paper, with the variable-order and fractal-fractional73

derivatives defined by three definitions (Atangana-Baleanu, Caputo-Fabrizo, and Caputo74

derivatives), we will apply stochastic derivatives to piecewise differential equations. The75

stability analysis of the fundamental system will be covered. Furthermore, we will in-76

vestigate the behaviour of stochastic derivatives using the Milstein technique and solve77

a crossover model based on a non-singular kernel described by Atangana-Baleanu deriva-78

tives using Toufik-Atangana method (TAM). Additionally, we employ the approxima-79

tion of integral Caputo-Fabrizio and Lagrange polynomial of two steps to approximate80

the deterministic model with a nonsingular exponential decay kernel. We will use the81

Grünwald−Letnikov nonstandard finite difference technique (GL-NSFDM) to approximate82

a determinate model utilising a singular kernel. The outcomes of infected persons derived83

from the suggested models are contrasted with actual data, whereby the US monkey-84

pox data (June 13-September 16, 2022). Furthermore, A comparison between the three85

methods with the Milstein method for solving the three crossover models is presented.86

A number of numerical simulations utilising variable order and fractional orders will be87

shown.88

The following is the study’s structure: pertinent definitions of fractional and variable-89

order Caputo derivatives are given in Section 2. Also, a brief overview of stochastic90

equations and a numerical method for solving them was presented in this section. Section91

3 presents three distinct generalised crossover mathematical models for the monkeypox92

virus based on principles of variable, fractal-fractional order, and the notions stochastic93

in three intervals; these crossover mathematical models defined by Atangana-Baleanu,94

Caputo-Fabrizo, Caputo derivatives. Section 4 theoretical analysis of basic model is ex-95

amined. Section 5 illustrates constructing GL-NSFDM to approximate the deterministic96

model using a singular kernel and the TAM, which uses a nonsingular Mittag-Leffler ker-97

nel to approximate the deterministic model. Moreover, Utilising the approximation of98

integral Caputo-Fabrizio and Lagrange polynomial of two steps to approximate the deter-99

ministic model with a nonsingular exponential decay kernel. Numerical simulations of the100

three suggested models are shown in Section 5. Lastly, a summary of the study’s main101

conclusions and contributions is given in Section 6.102



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 4 of 41

2. Fundamental Definitions103

We provide a few key definitions for fractional that will be utilised in the rest of this104

research in the section that follows.105

Definition 2.1. The derivatives of order α of Caputo’s for f(t) are determined by [9]:106

C
0 D

α
t (f(t)) =

1

Γ(n− α)

∫ t

0
(t− s)n−α−1f (n)(s)ds, t > 0, (1)

where n = N, n− 1 < α ≤ n, and Γ(·) is a Gamma function.107

108

Definition 2.2. The Caputo fractal-fractional derivative of a function f(t) is defined by109

fractal order β and fractional order α and a power law kernel as follows [9]:110

FFC
0 Dα,β

t (f(t)) =
1

Γ(n− α)

∫ t

0
(t− s)n−α−1(

d

dsβ
f(s))ds, (2)

where (n− 1) < α, β ≤ n ∈ N.111

Definition 2.3. (Atangana - Baleanu Caputo fractional derivative)[21]112

ABC
a Dα

t (f(t)) =
AB(α)

1− α

∫ t

a
f

′
(s)Eα(−α

(t− s)α

1− α
)ds, (3)

α ∈ [0, 1], AB(α) = 1 − α + α
Γ(α) and the Mittag Leffler function is an entire function113

defined by114

Eα(z) =
∑∞

k=0
zk

Γ(kα+1) .115

Definition 2.4. The fractal-fractional derivative of f(t) with the Mittag-Leffler type kernel116

and an order of α, may be found by [22].117

FFAB
0 Dα,β

t (f(t)) =
AB(α)

(1− α)

∫ t

0

d

dtβ
f(s)Eα(−α

(t− s)α

1− α
)ds, (4)

such that: α ∈ [0, 1], AB(α) = 1− α + α
Γ(α) and Eα(z) =

∑∞
k=0

zk

Γ(kα+1) characterises the118

Mittag Leffler function.119

Definition 2.5. The form of Atangana - Baleanu fractal-fractional integral is [21]120

FFAB
0 Iα,βt (f(t)) =

αβ

AB(α)Γ(α)

∫ t

0
f(s)(t− s)α−1sβ−1ds+

1− α

AB(α)
βtβ−1f(t), (5)

here α ∈ [0, 1], AB(α) = 1−α+ α
Γ(α) and Eα(z) =

∑∞
k=0

zk

Γ(kα+1) characterises the Mittag121

Leffler function.122



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 5 of 41

Definition 2.6. The definition of the Caputo-Fabrizio fractional derivative[15].123

CF
0 Dα

t (f(t)) =
M(α)

1− α

∫ t

0

d

ds
f(s)e(

−α
1−α

(t−s))ds, (6)

Consider α to be a fractional order constant with the value 0 < α ≤ 1 and M(α) =124

1− α+ α
Γ(α) .125

Definition 2.7. The fractional integral of the function f(t) with an exponential decay126

kernel is provided by [21]127

CF
0 Iαt (f(t)) =

α

M(α)

∫ t

0
f(s)ds+

1

M(α)
(1− α)f(t), (7)

In this case, M(α) represents a normalising function, M(0) = 1,M(1) = 1 and 0 < α ≤ 1.128

Definition 2.8. The Caputo-Fabrizio fractal-fractional derivative[22].129

FFCF
0 Dα,β

t (f(t)) =
M(α)

1− α

d

dtβ

∫ t

0
exp(− α

1− α
(t− s))f(s)ds, (8)

whereas M(0) = M(1) = 1, 0 < α ≤ 1, and β ≤ n ∈ N.130

Definition 2.9. With order α and an exponentially decaying type kernel, the fractal-131

fractional integral of f(t) is provided by [21].132

FFCF
0 Iα,βt (f(t)) =

αβ

M(α)

∫ t

0
sβ−1f(s)ds+

β(1− α)tβ−1f(t)

M(α)
, (9)

In this case, M(α) represents a normalising function, M(0) = 1,M(1) = 1 and 0 < α ≤ 1.133

2.1. Stochastic Differential Equation (SDE)134

Consider the real-valued process Y (t) with t ∈ [0, T ] satisfying the following a stochas-135

tic differential equation (SDE) [23].136

dY (t) = Φ(Y (t), t)dt+ σΨ(Y (t), t)∆B(t), (10)

where ∆B = B(tn+1)−B(tn) is the Brownian increment on [tn, tn+1].137

2.2. Milstein method138

By including a second-order ”correction” factor, which is obtained from the stochastic139

Taylor series expansion of Y (t) by using Ito’s lemma on the Φ(.) and Ψ(.) functions,140

Milstein technique improves the accuracy of the Euler-Maruyama method approximation.141

The following differential form is obtained using the Milstein technique [24]:142

Yn+1 − Yn = Φ(Yn, tn)h+Ψ(Yn, tn)∆Bn + 0.5Ψ′(Yn, tn)Ψ(Yn, tn)((∆Bn)
2 − h), (11)

where Ψ′(Y (t), t) denotes the derivative of Ψ(Y (t), t), ∆B = B(tn+1) − B(tn) is the143

Brownian increment on [tn, tn+1].144



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 6 of 41

3. The Crossover of Monkeypox Mathematical Model145

Here we modified the analysis of monkeypox virus using mathematics which given in146

[25] to a hybrid piecewise variable of fractional, fractal-fractional order, and stochastic147

Monkeypox disease mathematical model. Using the idea of a system of piecewise differen-148

tial equations, we will present three models; two models with nonsingular kernel and the149

third one with singular kernel. The impact of singular and non-singular kernels on the150

crossover monkeypox mathematical model primarily affects how the memory and hered-151

itary effects are incorporated into the model. Singular and non-singular kernels arise152

in fractional-order differential equations, which are increasingly used to model infectious153

diseases due to their ability to capture long-term memory effects.154

3.1. The Crossover of Monkeypox Mathematical Model with Nonsingular155

Kernel (CM1)156

A crossover monkeypox model considers both human-to-human and animal-to-human157

transmission pathways. The choice of kernel affects:158

• Epidemic Threshold (Basic Reproduction Number, R0): Singular kernels might lead159

to higher R0 due to stronger memory effects, whereas non-singular kernels may160

predict lower, more immediate transmission dynamics.161

• Peak Infection Time: Singular kernels delay the infection peak, while non-singular162

kernels lead to a more immediate peak.163

Applying the Atangana–Baleanu definition , it incorporates a nonsingular and nonlocal164

kernel, making it suitable for modeling memory effects in complex systems. Also, it gives165

rise to the idea of a piecewise system of differential equations. The Atangana–Baleanu166

variable order operator ’κ’ is used to extend the deterministic model in 0 < t ≤ t1, the167

fractal-fractional of the Atangana–Baleanu is employed in t1 < t ≤ t2, and the stochastic168

differential equation (SDE) in t2 < t ≤ Tf . The following is an expression of the new169

model:170

ABC
0 Dκ

t Sh = Λκ
h −

(βκ
1 Ir + βκ

2 Ih)Sh

Nh
− µκ

hSh + ϕκQh,

ABC
0 Dκ

t Eh =
(βκ

1 Ir + βκ
2 Ih)Sh

Nh
− (γκ1 + γκ2 + µκ

h)Eh,

ABC
0 Dκ

t Ih = γκ1Eh − (µκ
h + δκh + ρκ)Ih,

ABC
0 Dκ

t Qh = γκ2Eh − (ϕκ + τκ + µκ
h + δκh)Qh,

ABC
0 Dκ

t Rh = ρκIh + τκQh − µκ
hRh,

ABC
0 Dκ

t Sr = Λκ
r −

βκ
3SrIr
Nr

− µκ
rSr,

ABC
0 Dκ

t Er =
βκ
3SrIr
Nr

− (µκ
r + γκ3 )Er,



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 7 of 41

ABC
0 Dκ

t Ir = γκ3Er − (µκ
r + δκr )Ir, (12)

with initial conditions Sh0, Eh0, Ih0, Qh0, Rh0, Sr0, Er0, Ir0 at t = 0.171

Furthermore, if t1 < t ≤ t2, the model has the following expression:172

FFAB
0 Dα,β

t Sh = Λα
h − (βα

1 Ir + βα
2 Ih)Sh

Nh
− µα

hSh + ϕαQh,

FFAB
0 Dα,β

t Eh =
(βα

1 Ir + βα
2 Ih)Sh

Nh
− (γα1 + γα2 + µα

h)Eh,

FFAB
0 Dα,β

t Ih = γα1Eh − (µα
h + δαh + ρα)Ih,

FFAB
0 Dα,β

t Qh = γα2Eh − (ϕα + τα + µα
h + δαh )Qh,

FFAB
0 Dα,β

t Rh = ραIh + ταQh − µα
hRh,

FFAB
0 Dα,β

t Sr = Λα
r − βα

3 SrIr
Nr

− µα
r Sr,

FFAB
0 Dα,β

t Er =
βα
3 SrIr
Nr

− (µα
r + γα3 )Er,

FFAB
0 Dα,β

t Ir = γα3Er − (µα
r + δαr )Ir, (13)

with initial conditions173

Sh(t1) = Sh1, Eh(t1) = Eh1, Ih(t1) = Ih1, Sr(t1) = Sr1,

Rh(t1) = Rh1, Er(r1) = Er1, Qh(t1) = Qh1, Ir(t1) = Ir1.

Additionally, if t2 < t ≤ Tf , the model has the following expression:174

dSh = (Λh −
(β1Ir + β2Ih)Sh

Nh
− µhSh + ϕQh)dt+ σ1ShdB1,

dEh = (
(β1Ir + β2Ih)Sh

Nh
− (γ1 + γ2 + µh)Eh)dt+ σ2EhdB2,

dIh = (γ1Eh − (µh + δh + ρ)Ih)dt+ σ3IhdB3,

dQh = (γ2Eh − (ϕ+ τ + µh + δh)Qh)dt+ σ4QhdB4,

dRh = (ρIh + τQh − µhRh)dt+ σ5RhdB5,

dSr = (Λr −
β3SrIr
Nr

− µrSr)dt+ σ6SrdB6,

dEr = (
β3SrIr
Nr

− (µr + γ3)Er)dt+ σ7ErdB7,

dIr = (γ3Er − (µr + δr)Ir)dt+ σ8IrdB8. (14)
175

Sh(t2) = Sh2, Eh(t2) = Eh2, Ih(t2) = Ih2, Qh(t2) = Qh2,

Rh(t2) = Rh2, Sr(t2) = Sr2, Er(r2) = Er2, Ir(t2) = Ir2.

Table 1 provides explanations of the variables used in the model, whereas Table 2176

includes the parameters and their corresponding values.177



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 8 of 41

Table 1: Variables in the system.

The variable Description

Nh The population of humans (Sh + Eh + Ih +Qh +Rh).

Sh Susceptible class.

Eh Exposed class.

Ih Infected class.

Qh Isolation class.

Rh Recovery class.

Nr The population of rodent (Sr + Er + Ir).

Sr Susceptible class.

Er Exposed class.

Ir Infected class.

Table 2: The values of the model’s parameters.

Parameter Description Value

Λh The recruitment rate of humans 0.34857
Λr The recruitment rate of rodents 0.60822
β1 The contact rate between rodents and humans 0.00103563
β2 The contact rate between humans 0.99993476
β3 The contact rate between rodents 0.50838481
α1 The transmission rate of humans from exposed to infectious class 0.07296767
α2 The rate of identification of suspected case 0.00000213
α3 The transmission rate of rodents from exposed to infectious class 0.03353088
ϕ The proportion of humans not detected after diagnosis 0.33749119
τ The rate of transmission from isolation to recovered class 0.99128100
ρ The recovery rate 0.16786558
µh The natural death rates of humans 1/(79× 365)
µr The natural death rates of rodents 1/(5× 365)
δh The disease induced death rates of humans 0.18291202
δr The disease induced death rates of rodents 0.00000255

3.2. The Crossover of Monkeypox Mathematical Model with Exponential178

Decay Kernel (CM2)179

The Caputo-Fabrizio (CF) fractional derivative is an alternative fractional operator180

introduced to overcome certain limitations of classical fractional derivatives, particularly181

in handling memory effects with a nonsingular kernel. Applying the Caputo–Fabrizio def-182

inition (nonsingular kernel) gives rise to The concept of a system of piecewise differential183

equations. The Caputo–Fabrizio variable order operator ’κ’ is used to extend the deter-184



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 9 of 41

ministic model in 0 < t ≤ t1, the fractal-fractional of the Caputo–Fabrizio is employed185

in t1 < t ≤ t2, and the stochastic differential equation (SDE) in t2 < t ≤ Tf . Here is an186

expression for the new model:187

CF
0 Dκ

t Sh = Λκ
h −

(βκ
1 Ir + βκ

2 Ih)Sh

Nh
− µκ

hSh + ϕκQh,

CF
0 Dκ

t Eh =
(βκ

1 Ir + βκ
2 Ih)Sh

Nh
− (γκ1 + γκ2 + µκ

h)Eh,

CF
0 Dκ

t Ih = γκ1Eh − (µκ
h + δκh + ρκ)Ih,

CF
0 Dκ

t Qh = γκ2Eh − (ϕκ + τκ + µκ
h + δκh)Qh,

CF
0 Dκ

t Rh = ρκIh + τκQh − µκ
hRh,

CF
0 Dκ

t Sr = Λκ
r −

βκ
3SrIr
Nr

− µκ
rSr,

CF
0 Dκ

t Er =
βκ
3SrIr
Nr

− (µκ
r + γκ3 )Er,

CF
0 Dκ

t Ir = γκ3Er − (µκ
r + δκr )Ir, (15)

with initial conditions Sh0, Eh0, Ih0, Qh0, Rh0, Sr0, Er0, Ir0 at t = 0.188

Furthermore, if t1 < t ≤ t2, the model has the following expression:189

FFCF
0 Dα,β

t Sh = Λα
h − (βα

1 Ir + βα
2 Ih)Sh

Nh
− µα

hSh + ϕαQh,

FFCF
0 Dα,β

t Eh =
(βα

1 Ir + βα
2 Ih)Sh

Nh
− (γα1 + γα2 + µα

h)Eh,

FFCF
0 Dα,β

t Ih = γα1Eh − (µα
h + δαh + ρα)Ih,

FFCF
0 Dα,β

t Qh = γα2Eh − (ϕα + τα + µα
h + δαh )Qh,

FFCF
0 Dα,β

t Rh = ραIh + ταQh − µα
hRh,

FFCF
0 Dα,β

t Sr = Λα
r − βα

3 SrIr
Nr

− µα
r Sr,

FFCF
0 Dα,β

t Er =
βα
3 SrIr
Nr

− (µα
r + γα3 )Er,

FFCF
0 Dα,β

t Ir = γα3Er − (µα
r + δαr )Ir, (16)

with initial conditions190

Sh(t1) = Sh1, Eh(t1) = Eh1, Ih(t1) = Ih1, Sr(t1) = Sr1,

Rh(t1) = Rh1, Er(r1) = Er1, Qh(t1) = Qh1, Ir(t1) = Ir1.

Additionally, if t2 < t ≤ Tf , the model has the following expression:191

dSh = (Λh −
(β1Ir + β2Ih)Sh

Nh
− µhSh + ϕQh)dt+ σ1ShdB1,



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 10 of 41

dEh = (
(β1Ir + β2Ih)Sh

Nh
− (γ1 + γ2 + µh)Eh)dt+ σ2EhdB2,

dIh = (γ1Eh − (µh + δh + ρ)Ih)dt+ σ3IhdB3,

dQh = (γ2Eh − (ϕ+ τ + µh + δh)Qh)dt+ σ4QhdB4,

dRh = (ρIh + τQh − µhRh)dt+ σ5RhdB5,

dSr = (Λr −
β3SrIr
Nr

− µrSr)dt+ σ6SrdB6,

dEr = (
β3SrIr
Nr

− (µr + γ3)Er)dt+ σ7ErdB7,

dIr = (γ3Er − (µr + δr)Ir)dt+ σ8IrdB8. (17)

192

Sh(t2) = Sh2, Eh(t2) = Eh2, Ih(t2) = Ih2, Qh(t2) = Qh2,

Rh(t2) = Rh2, Sr(t2) = Sr2, Er(r2) = Er2, Ir(t2) = Ir2.

3.3. The Crossover of Monkeypox Mathematical Model with Singular193

Kernel (CM3)194

The Caputo definition of a fractional derivative is one of the most commonly used195

formulations in fractional calculus, particularly in applications involving physical and en-196

gineering problems.197

By using the Caputo formulation (singular kernel), the concept of a piecewise differen-198

tial equation system is introduced. In 0 < t ≤ t1, the model of deterministic behaviour is199

extended using the Caputo variable order operator, and in t1 < t ≤ t2, fractal-fractional200

of caputo variable order is used where the range of t2 < t ≤ tf contains an expansion of a201

stochastic differential equation (SDE). The following is an expression for the new model:202

C
0 D

κ
t Sh = Λκ

h −
(βκ

1 Ir + βκ
2 Ih)Sh

Nh
− µκ

hSh + ϕκQh,

C
0 D

κ
t Eh =

(βκ
1 Ir + βκ

2 Ih)Sh

Nh
− (γκ1 + γκ2 + µκ

h)Eh,

C
0 D

κ
t Ih = γκ1Eh − (µκ

h + δκh + ρκ)Ih,
C
0 D

κ
t Qh = γκ2Eh − (ϕκ + τκ + µκ

h + δκh)Qh,
C
0 D

κ
t Rh = ρκIh + τκQh − µκ

hRh,

C
0 D

κ
t Sr = Λκ

r −
βκ
3SrIr
Nr

− µκ
rSr,

C
0 D

κ
t Er =

βκ
3SrIr
Nr

− (µκ
r + γκ3 )Er,

C
0 D

κ
t Ir = γκ3Er − (µκ

r + δκr )Ir, (18)

with initial conditions Sh0, Eh0, Ih0, Qh0, Rh0, Sr0, Er0, Ir0 at t = 0.203



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 11 of 41

Furthermore, if t1 < t ≤ t2, the model has the following expression:204

FFC
0 Dα,β

t Sh = Λα
h − (βα

1 Ir + βα
2 Ih)Sh

Nh
− µα

hSh + ϕαQh,

FFC
0 Dα,β

t Eh =
(βα

1 Ir + βα
2 Ih)Sh

Nh
− (γα1 + γα2 + µα

h)Eh,

FFC
0 Dα,β

t Ih = γα1Eh − (µα
h + δαh + ρα)Ih,

FFC
0 Dα,β

t Qh = γα2Eh − (ϕα + τα + µα
h + δαh )Qh,

FFC
0 Dα,β

t Rh = ραIh + ταQh − µα
hRh,

FFC
0 Dα,β

t Sr = Λα
r − βα

3 SrIr
Nr

− µα
r Sr,

FFC
0 Dα,β

t Er =
βα
3 SrIr
Nr

− (µα
r + γα3 )Er,

FFC
0 Dα,β

t Ir = γα3Er − (µα
r + δαr )Ir, (19)

with initial conditions205

Sh(t1) = Sh1, Eh(t1) = Eh1, Ih(t1) = Ih1, Sr(t1) = Sr1,

Rh(t1) = Rh1, Er(r1) = Er1, Qh(t1) = Qh1, Ir(t1) = Ir1.

Additionally, if t2 < t ≤ Tf , the model has the following expression:206

dSh = (Λh −
(β1Ir + β2Ih)Sh

Nh
− µhSh + ϕQh)dt+ σ1ShdB1,

dEh = (
(β1Ir + β2Ih)Sh

Nh
− (γ1 + γ2 + µh)Eh)dt+ σ2EhdB2,

dIh = (γ1Eh − (µh + δh + ρ)Ih)dt+ σ3IhdB3,

dQh = (γ2Eh − (ϕ+ τ + µh + δh)Qh)dt+ σ4QhdB4,

dRh = (ρIh + τQh − µhRh)dt+ σ5RhdB5,

dSr = (Λr −
β3SrIr
Nr

− µrSr)dt+ σ6SrdB6,

dEr = (
β3SrIr
Nr

− (µr + γ3)Er)dt+ σ7ErdB7,

dIr = (γ3Er − (µr + δr)Ir)dt+ σ8IrdB8. (20)
207

Sh(t2) = Sh2, Eh(t2) = Eh2, Ih(t2) = Ih2, Sr(t2) = Sr2,

Rh(t2) = Rh2, Er(r2) = Er2, Qh(t2) = Qh2, Ir(t2) = Ir2.

4. Theoretical Analysis of Model208

To obtain the equilibrium points of the crossover models, we put all derivatives equal209

to zero then we have two points as follows:210



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 12 of 41

Disease-free equilibrium [26], when there is no sickness, and given as follows:211

ϵ0 = (
Λh

µh
, 0, 0, 0, 0,

Λr

µr
, 0, 0). (21)

Endemic equilibrium [26] when the virus continues to spread across the population,212

and given as follows:213

ϵ1 = (S∗∗
h , E∗∗

h , I∗∗h , Q∗∗
h , R∗∗

h , S∗∗
r , E∗∗

r , I∗∗r )
214

S∗∗
h =

k1k3Λh

µhk1k3 − α2ΦΦh + k1k3Φh
,

E∗∗
h =

k3ΛhΦh

µhk1k3 − α2ΦΦh + k1k3Φh
,

I∗∗h =
k3α1ΛhΦh

k2(µhk1k3 − α2ΦΦh + k1k3Φh)
,

Q∗∗
h =

α2ΛhΦh

µhk1k3 − α2ΦΦh + k1k3Φh
,

R∗∗
h =

(k3ρα1 + k2τα2)ΛhΦh

µhk2(µhk1k3 − α2ΦΦh + k1k3Φh)
,

S∗∗
r =

λr

µr +Φr
,

E∗∗
r =

λr

k4(µr +Φr)
,

I∗∗r =
Φrα3λr

k4k5(µr +Φr)
.

(22)

The formulas k1 = µh +α2 +α1, k2 = ρ+ δh +µh, k3 = δh + τ +µh +Φ, k4 = α3 +µr,215

k5 = δr + µr, Φh =
β1I∗r+β2T ∗∗

h
Nh

, Φr =
β3I∗∗r
Nr

.216

217

Reproduction number [26]218

R0 =
α1β2

(α1 + α2 + µh)(µh + δh + ρ)
. (23)

Stability of endemic equilibrium point219

The Jacobian matrix about the endemic equilibria ϵ1 is given as:220

J(ϵ1) =

−(
β1Ir+β2Ih

Nh
) − µh 0 − β2Sh

Nh
ϕ 0 0 0 − β1Sh

Nh
β1Ir+β2Ih

Nh
−α1 − α2 − µh

β2Sh
Nh

0 0 0 0
β1Sh
Nh

0 α1 −µh − δh − γ 0 0 0 0 0
0 α2 0 −ϕ − τ − δh − µh 0 0 0 0
0 0 γ τ −µh 0 0 0

0 0 0 0 0 − β3Ir
Nr

− µr 0 − β3Sr
Nr

0 0 0 0 0
β3Ir
Nr

−µr − α3
β3Sr
Nr

0 0 0 0 0 0 α3 µr − δr


,



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 13 of 41

then the characteristic equation of Jacobian is given as:221

1

NhNr
[(−x− µh) (−ϕα2 (Irβ1 + Ihβ2) (x+ ρ+ δh + µh) + (−x− τ − ϕ− δh − µh) (Shα1β2 (x+ µh)

− (x+ α1 + α2 + µh) (x+ ρ+ δh + µh) (Irβ1 + Ihβ2 +Nh(x+ µh)))) (Srα3β3(x+ µr)

−(x+ α3 + µr) (Irβ3 +Nr(x+ µr)) (x+ (µr + δr)))] = 0. (24)

This, when the polynomial is converted to standard form, may be further expressed as222

follows:223

x8 + C1x
7 ++C2x

6 + C3x
5 + C4x

4 + C5x
3 + C6x

2 + C7x+ C8 = 0,

where C ′
is are the coefficients of x8−i; i = 1, 2, . . . 8. Therefore, the Monkeypox endemic224

equilibrium (MFE) is asymptotically stable if R0 > 1 and C1 > 0, C1C2 > C3, C1C2C3 +225

C0C1C5 > C0C
2
3 + C2

1C4, P
∗Q > PQ∗, MQ∗ > P ∗N, M∗N > MN∗, XN∗ > TM∗ such226

that:227

P =
C1C2 − C0C3

C1
, Q =

C1C4 − C0C5

C1

R =
C1C6 − C0C7

C1
, S = C8

P ∗ =
PC3 − C1Q

P
, Q∗ =

PC5 − C1R

P

R∗ =
PC7 − C1S

P
, M =

P ∗Q− PQ∗

P ∗

N =
P ∗R− PR∗

P ∗ , T =
PS

P ∗

M∗ =
MQ∗ − P ∗N

M
, N∗ =

MR∗ − P ∗T

M

X =
M∗N −MN∗

M∗ .

5. Numerical Schemes for Crossover Models228

5.1. Numerical Scheme with Nonsingular Kernel229

In this section, we will present the numerical methods to solve the crossover models230

with nonsingular kernel (12)-(14) and (15)-(17) and the crossover model with singular231

kernel (18)-(20) as follows:232

5.1.1. Toufik-Atangana method233

Consider the general formula of the variable-order Atangana-Baleanu differential equation234

given in (0, t1] as follows:235

ABC
0 Dκ

t Y (t) = f(t, Y (t)), Y (0) = Y0, κ = α(t) ∈ (0, 1]. (25)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 14 of 41

By integrate (25) using (5). Then by using the approximation of the variable-order236

Atangana-Baleanu integral equation and Lagrange polynomial of two step we have the237

following formula to approximate the variable-order Atangana-Baleanu differential equa-238

tion [21]:239

Yp+1 = Y0 +
1− κ

AB(κ)
f(tp, Yp) +

κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

f(tp, Yp)[(−s+ 1 + p)κ

(−s+ p+ κ+ 2)− (−s+ p)κ(−s+ p+ 2κ+ 2)]− κhκ

AB(k)Γ(κ+ 2)
p∑

s=0

[(−s+ p+ 1)κ+1 − (−s+ p)κ(−s+ p+ κ+ 1)]f(ts−1, Ys−1). (26)

When the suggested numerical method (26) is applied to system (12), the outcomes are240

Shp+1 = Sh0 +
1− κ

AB(κ)
(Λh −

(β1Irp + β2Ihp)Shp

Nhp
− µhShp + ϕQhp) +

κhκ

AB(κ)Γ(κ+ 2)

×
p∑

s=0

(Λh −
(β1Irs + β2Ihs)Shs

Nhs
− µhShs + ϕQhs)[(−s+ p+ 1)κ(−s+ κ+ 2 + p)

− (−s+ p)κ(−s+ 2κ+ 2 + p)]− κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

[
Λh −

(β1Ir(s−1) + β2Ih(s−1))

Nh(s−1)

Sh(s−1) − µhSh(s−1) + ϕQh(s−1)

]
[(−s+ p+ 1)κ+1 − (−s+ p)κ(−s+ κ+ 1 + p)],

Ehp+1 = Eh0 +
1− κ

AB(κ)
(
(β1Irp + β2Ihp)Shp

Nhp
− (γ1 + γ2 + µh)Ehp) +

κhκ

AB(κ)Γ(κ+ 2)

×
p∑

s=0

(
(β1Irs + β2Ihs)Shs

Nhs
− (γ1 + γ2 + µh)Ehs)[(−s+ p+ 1)κ(−s+ p+ 2 + κ)

− (−s+ p)κ(−s+ 2 + 2κ+ p)]− κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

[
(β1Ir(s−1) + β2Ih(s−1))

Nh(s−1)

∗Sh(s−1) − (γ1 + γ2 + µh)Eh(s−1)

]
[(−s+ p+ 1)κ+1 − (−s+ p)κ(−s+ κ+ 1 + p)],

Ihp+1 = Ih0 +
1− κ

AB(κ)
(γ1Ehp − (µh + δh + ρ)Ihp) +

κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

[γ1Ehs

−(µh + δh + ρ)Ihs] [(−s+ p+ 1)κ(−s+ κ+ 2 + p)− (−s+ p)κ(−s+ 2κ+ 2 + p)]

− κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

(γ1Eh(s−1) − (µh + δh + ρ)Ih(s−1))
[
(−s+ p+ 1)κ+1

−(−s+ p)κ(−s+ 1 + κ+ p)] ,

Qhp+1 = Qh0 +
1− κ

AB(κ)
(γ2Ehp − (ϕ+ τ + µh + δh)Qhp) +

κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

[γ2Ehs



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 15 of 41

−(ϕ+ τ + µh + δh)Qhs] [(−s+ p+ 1)κ(−s+ p+ κ+ 2)− (−s+ p)κ(−s+ p+ 2κ+ 2)]

− κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

(γ2Eh(s−1) − (ϕ+ τ + µh + δh)Qh(s−1))[(−s+ p+ 1)κ+1

− (−s+ p)κ(−s+ p+ 1 + κ)],

Rhp+1 = Rh0 +
1− κ

AB(κ)
(ρIhp + τQhp − µhRhp) +

κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

(ρIhs + τQhs − µhRhs)

[(−s+ p+ 1)κ(−s+ 2 + κ+ p)− (−s+ p)κ(−s+ 2 + 2κ+ p)]− κhκ

AB(κ)Γ(κ+ 2)

×
p∑

s=0

(ρIh(s−1) + τQh(s−1) − µhRh(s−1))[(−s+ 1 + p)κ+1 − (−s+ p)κ(−s+ κ+ 1 + p)],

Srp+1 = Sr0 +
1− κ

AB(κ)
(Λr −

β3SrpIrp
Nrp

− µrSrp) +
hκκ

Γ(2 + κ)AB(κ)

p∑
s=0

[
Λh −

β3SrsIrs
Nrs

−µrSrs] [(−s+ p+ 1)κ(−s+ κ+ 2 + p)− (−s+ p)κ(−s+ 2κ+ 2 + p)]

− κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

(Λh −
β3Sr(s−1)Ir(s−1)

Nr(s−1)
− µrSr(s−1))[(−s+ p+ 1)κ+1

− (−s+ p)κ(−s+ p+ κ+ 1)],

Erp+1 = Er0 +
1− κ

AB(κ)
(
β3SrpIr
Nrp

− (µr + γ3)Erp) +
κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

[
β3SrsIrs

Nrs

−(µr + γ3)Ers] [(−s+ p+ 1)κ(−s+ κ+ 2 + p)− (−s+ 2)κ(−s+ 2κ+ 2 + p)]

− κhκ

AB(κ)Γ(κ+ 2)

p∑
s=0

(
β3Sr(s−1)Ir(s−1)

Nr(s−1)
− (µr + γ3)Er(s−1))[(−s+ p+ 1)κ+1

− (−s+ p)κ(−s+ p+ κ+ 1)],

Irp+1 = Ir0 +
1− κ

AB(κ)
(γ3Erp − (µr + δr)Irp) +

α(t)hκ

AB(κ)Γ(κ+ 2)

p∑
s=0

[γ3Ers − (µr + δr)Irs]

[(−s+ p+ 1)κ(−s+ p+ 2 + κ)− (−s+ p)κ(−s+ p+ 2κ+ 2)]− κhκ

AB(κ)Γ(κ+ 2)

×
p∑

s=0

(γ3Er(s−1) − (µr + δr)Ir(s−1))[(−s+ p+ 1)κ+1 − (−s+ p)κ(−s+ p+ κ+ 1)].

(27)

Now consider the general formula for fractal-fractional Atangana-Baleanu deferential241

equation in the (t1, t2] given as follows:242

FFABC
0 Dα,β

t Y (t) = f(t, Y (t)), Y (0) = Y0, α, β ∈ (0, 1]. (28)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 16 of 41

Using [21] the approximation of (28) is given as follows:243

Yp+1 = Y0 +
1− α

AB(α)
tβ−1
p f(tp, Yp) +

αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
p f(tp, Yp) [(−s+ p+ 1)α(−s+ p+ α+ 2)

−(−s+ p)α(−s+ p+ 2α+ 2)]− αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
s−1 f(ts−1, Ys−1)

[
(−s+ p+ 1)α+1

−(−s+ α+ p+ 1)(−s+ p)α] . (29)

Utilising the proposed numerical technique TAM (5.1.1) in the system (13), we obtain:244

Shp+1 = Sh0 +
1− α

AB(α)
tβ−1
p (Λh −

(β1Irp + β2Ihp)Shp

Nhp
− µhShp + ϕQhp) +

αβhα

AB(α)Γ(α+ 2)

×
p∑

s=0

tβ−1
p (Λh −

(β1Irs + β2Ihs)Shs

Nhs
− µhShs + ϕQhs)[(−s+ α+ 2 + p)(−s+ p+ 1)α

− (−s+ p)α(−s+ 2α+ 2 + p)]− αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
s−1

[
Λh −

(β1Ir(s−1) + β2Ih(s−1))

Nh(s−1)

Sh(s−1) − µhSh(s−1) + ϕQh(s−1)

]
[(−s+ 1 + p)α+1 − (−s+ p)α(−s+ α+ 1 + p)],

Ehp+1 = Eh0 +
1− α

AB(α)
tβ−1
p (

(β1Irp + β2Ihp)Shp

Nhp
− (γ1 + γ2 + µh)Ehp) +

αβhα

AB(α)Γ(α+ 2)

×
p∑

s=0

tβ−1
p (

(β1Irs + β2Ihs)Shs

Nhs
− (γ1 + γ2 + µh)Ehs)[(−s+ p+ 1)α(−s+ α+ 2 + p)

− (−s+ p)α(−s+ 2α+ 2 + p)]− αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
s−1 (

(β1Ir(s−1) + β2Ih(s−1))Sh(s−1)

Nh(s−1)

− (γ1 + γ2 + µh)Eh(s−1))[(−s+ p+ 1)α+1 − (−s+ α+ 1 + p)(−s+ p)α],

Ihp+1 = Ih0 +
1− α

AB(α)
tβ−1
p (γ1Ehp − (µh + δh + ρ)Ihp) +

αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
p [γ1Ehs

−(µh + δh + ρ)Ihs] [(−s+ p+ 1)α(−s+ p+ α+ 2)− (−s+ p)α(−s+ p+ 2α+ 2)]

− αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
s−1 (γ1Eh(s−1) − (µh + δh + ρ)Ih(s−1))[(−s+ p+ 1)α+1

− (−s+ p)α(−s+ p+ 1 + α)],

Qhp+1 = Qh0 +
1− α

AB(α)
tβ−1
p (γ2Ehp − (ϕ+ τ + µh + δh)Qhp) +

αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
p [γ2Ehs

−(ϕ+ τ + µh + δh)Qhs] [(−s+ p+ 1)α(−s+ 2 + α+ p)− (−s+ p+ 2 + 2α)(−s+ p)α]

− αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
s−1 (γ2Eh(s−1) − (ϕ+ τ + µh + δh)Qh(s−1))[(−s+ 1 + p)α+1



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 17 of 41

− (−s+ p)α(−s+ α+ 1 + p)],

Rhp+1 = Rh0 +
1− α

AB(α)
tβ−1
p (ρIhp + τQhp − µhRhp) +

hααβ

Γ(α+ 2)AB(α)

p∑
s=0

tβ−1
p [ρIhs + τQhs

−µhRhs] [(−s+ 2 + α+ p)(−s+ 1 + p)α − (−s+ 2 + 2α+ p)(−s+ p)α]− αβhα

AB(α)Γ(α+ 2)

×
p∑

s=0

tβ−1
s−1 (ρIh(s−1) + τQh(s−1) − µhRh(s−1))[−(−s+ p+ 1 + α)(−s+ p)α + (−s+ p+ 1)α+1],

Srp+1 = Sr0 +
1− α

AB(α)
tβ−1
p (Λr −

β3SrpIrp
Nrp

− µrSrp) +
αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
p (Λh −

β3SrsIrs
Nrs

− µrSrs)[(−s+ p+ 2 + α)(−s+ p+ 1)α − (−s+ 2 + 2α+ p)(−s+ p)α]− αβhα

AB(α)Γ(α+ 2)

×
p∑

s=0

tβ−1
s−1 (Λh −

β3Sr(s−1)Ir(s−1)

Nr(s−1)
− µrSr(s−1))[(−s+ 1 + p)α+1 − (−s+ p)α(−s+ 1 + α+ p)],

Erp+1 = Er0 +
1− α

AB(α)
tβ−1
p (

β3SrpIr
Nrp

− (µr + γ3)Erp) +
αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
p (

β3SrsIrs
Nrs

− (µr + γ3)Ers)[(−s+ 1 + p)α(−s+ 2 + α+ p)− (−s+ p)α(−s+ 2 + 2α+ p)]− αβhα

AB(α)Γ(α+ 2)

×
p∑

s=0

tβ−1
s−1 (

β3Sr(s−1)Ir(s−1)

Nr(s−1)
− (µr + γ3)Er(s−1))[(−s+ 1 + p)α+1 − (−s+ p)α(−s+ 1 + α+ p)],

Irp+1 = Ir0 +
1− α

AB(α)
tβ−1
p (γ3Erp − (µr + δr)Irp) +

αβhα

AB(α)Γ(α+ 2)

p∑
s=0

tβ−1
p (γ3Ers − (µr + δr)Irs)

[(−s+ p+ 1)α(−s+ 2 + α+ p)− (−s+ 2 + 2α+ p)(−s+ p)α]− αβhα

AB(α)Γ(α+ 2)

×
p∑

s=0

tβ−1
s−1 (γ3Er(s−1) − (µr + δr)Ir(s−1))[(−s+ 1 + p)α+1 − (−s+ p)α(−s+ 1 + α+ p)].

(30)

245

Remark 5.1. The stability and error analysis for this method are found in [21].246

During the third interval (t2, Tf ]), to approximate the system (14), we will use Milstein247

method (11). Then the explicit solution given as follows:248

Shp+1 = Shp + (Λh −
(β1Ir + β2Ih)Sh

Nh
− µhSh + ϕQh)h+ σ1Shp∆Bp + 0.5σ2

1Shp

((∆Bp)
2 − h),



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 18 of 41

Ehp+1 = Ehp + (
(β1Ir + β2Ih)Sh

Nh
− (γ1 + γ2 + µh)Eh)h+ σ2Ehp∆Bp + 0.5σ2

2Ehp

((∆Bp)
2 − h),

Ihp+1 = Ihp + (γ1Eh − (µh + δh + ρ)Ih)h+ σ3Ihp∆Bp + 0.5σ2
3Ihp((∆Bp)

2 − h),

Qhp+1 = Qhp + (γ2Eh − (ϕ+ τ + µh + δh)Qh)h+ σ4Qhp∆Bp + 0.5σ2
4Qhp((∆Bp)

2 − ϕh),

Rhp+1 = Rhp + (ρIh + τQh − µhRh)ϕh+ σ5 ∗Rhp∆Bp + 0.5σ2
5Rhp((∆Bp)

2 − h),

Srp+1 = Srp + (Λr −
β3SrIr
Nr

− µrSr)h+ σ6Srp∆Bp + 0.5σ2
6Srp((∆Bp)

2 − h),

Erp+1 = Erp + (
β3SrIr
Nr

− (µr + γ3)Er)h+ σ7Erp∆Bp + 0.5σ2
7Erp((∆Bp)

2 − h),

Irp+1 = Irp + (γ3Er − (µr + δr)Ir)h+ σ8Irp∆Bp + 0.5σ2
8Irp((∆Bp)

2 − h).

(31)

5.2. Numerical Scheme with Exponential Decay Kernel249

Consider the general formula of the variable-order Caputo Fabrizo variable-order dif-250

ferential equation in (0, t1] as follows:251

CF
0 Dκ

t Y (t) = f(t, Y (t)), Y (0) = Y0, κ = α(t) ∈ (0, 1]. (32)

First integrate (32) using (9). Then by using the approximation of the variable-order Ca-252

puto Fabrizo integral equation and Lagrange polynomial of two step we have the following253

formula to approximate the variable-order Caputo Fabrizo differential equation [21]:254

Yp+1 = Yp +
1− κ

M(κ)
[f(tp, Yp)− f(tp−1, Yp−1)] +

κ

M(κ)
[f(tp, Yp)

3

2
h− f(tp−1, Yp−1)

1

2
h].

(33)

Applying the recommended numerical approach (33) to approximate the system (17) yields255

the following results:256

Shp+1 = Shp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(Λh −

(β1Irp + β2Ihp)Shp

Nhp
− µhShp + ϕQhp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(Λh −

(β1Ir(p−1) + β2Ih(p−1))Sh(p−1)

Nh(p−1)
− µhSh(p−1) + ϕQh(p−1)),

Ehp+1 = Ehp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(
(β1Irp + β2Ihp)Shp

Nhp
− (γ1 + γ2 + µh)Ehp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(

(β1Ir(p−1) + β2Ih(p−1))Sh(p−1)

Nh(p−1)
− (γ1 + γ2 + µh)Eh(p−1)),

Ihp+1 = Ihp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(γ1Ehp − (µh + δh + ρ)Ihp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(γ1Ehp−1 − (µh + δh + ρ)Ihp−1),



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 19 of 41

Qhp+1 = Qhp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(γ2Ehp − (ϕ+ τ + µh + δh)Qhp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(γ2Eh(p−1) − (ϕ+ τ + µh + δh)Qh(p−1)),

Rhp+1 = Rhp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(ρIhp + τQhp − µhRhp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(ρIh(p−1) + τQh(p−1) − µhRh(p−1)),

Srp+1 = Srp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(Λr −

β3SrpIrp
Nrp

− µrSrp)

− (
−κ+ 1

M(κ)
+ h

κ

2M(κ)
)
β3Sr(p−1)Ir(p−1)

Nr(p−1)
− µrSr(p−1)),

Erp+1 = Erp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(
β3SrpIr
Nrp

− (µr + γ3)Erp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(

β3Sr(p−1)Ir(p−1)

Nr(p−1)
− (µr + γ3)Er(p−1)),

Irp+1 = Irp + (
1− κ

M(κ)
+ h

3κ

2M(κ)
)(γ3Erp − (µr + δr)Irp)

− (
1− κ

M(κ)
+

κ

2M(κ)
h)(γ3Er(p−1) − (µr + δr)Ir(p−1)).

(34)

Consider the general formula for fractal-fractional Caputo Fabrizo defferential equation257

in the (t1, t2] given as follows:258

FFCF
0 Dα,β

t Y (t) = f(t, Y (t)), Y (0) = Y0, α, β ∈ (0, 1]. (35)

Using [21] the approximation of (35) is given as follows:259

Yp+1 = Yp +
1− α

M(α)
[βtβ−1

p f(tp, Yp)− βtβ−1
p−1f(tp−1, Yp−1)] +

α

M(α)

[
βtβ−1

p f(tp, Yp)
3

2
h

−βtβ−1
p f(tp−1, Yp−1)

1

2
h

]
. (36)

Applying the recommended numerical approach (36) to the system (16) yields the following260

results:261

Shp+1 = Shp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(Λh −
(β1Irp + β2Ihp)Shp

Nhp
− µhShp + ϕQhp)− (

1− α

M(α)

+
α

2M(α)
h)(βtβ−1

p )(Λh −
(β1Ir(p−1) + β2Ih(p−1))Sh(p−1)

Nh(p−1)
− µhSh(p−1) + ϕQh(p−1)),

Ehp+1 = Ehp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(
(β1Irp + β2Ihp)Shp

Nhp
− (γ1 + γ2 + µh)Ehp)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 20 of 41

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(
(β1Ir(p−1) + β2Ih(p−1))Sh(p−1)

Nh(p−1)
− (γ1 + γ2 + µh)Eh(p−1)),

Ihp+1 = Ihp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(γ1Ehp − (µh + δh + ρ)Ihp)

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(γ1Ehp−1 − (µh + δh + ρ)Ihp−1),

Qhp+1 = Qhp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(γ2Ehp − (ϕ+ τ + µh + δh)Qhp)

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(γ2Eh(p−1) − (ϕ+ τ + µh + δh)Qh(p−1)),

Rhp+1 = Rhp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(ρIhp + τQhp − µhRhp)

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(ρIh(p−1) + τQh(p−1) − µhRh(p−1)),

Srp+1 = Srp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(Λr −
β3SrpIrp

Nrp
− µrSrp)

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(
β3Sr(p−1)Ir(p−1)

Nr(p−1)
− µrSr(p−1)),

Erp+1 = Erp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(
β3SrpIr
Nrp

− (µr + γ3)Erp)

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(
β3Sr(p−1)Ir(p−1)

Nr(p−1)
− (µr + γ3)Er(p−1)),

Irp+1 = Irp + (
1− α

M(α)
+

3α

2M(α)
h)(βtβ−1

p )(γ3Erp − (µr + δr)Irp)

− (
1− α

M(α)
+

α

2M(α)
h)(βtβ−1

p )(γ3Er(p−1) − (µr + δr)Ir(p−1)).

(37)

262

Remark 5.2. The stability and error analysis for this method are found in [21].263

The third interval (t2, Tf ]), to approximate the system (17), we will use Milstein264

method (11). Then the explicit solution given as follows:265

Shp+1 = Shp + (Λh −
(β1Ir + β2Ih)Sh

Nh
− µhSh + ϕQh)h+ σ1Shp∆Bp + 0.5σ2

1Shp

((∆Bp)
2 − h),

Ehp+1 = Ehp + (
(β1Ir + β2Ih)Sh

Nh
− (γ1 + γ2 + µh)Eh)h+ σ2Ehp∆Bp + 0.5σ2

2Ehp

((∆Bp)
2 − h),



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 21 of 41

Ihp+1 = Ihp + (γ1Eh − (µh + δh + ρ)Ih)h+ σ3Ihp∆Bp + 0.5σ2
3Ihp((∆Bp)

2 − h),

Qhp+1 = Qhp + (γ2Eh − (ϕ+ τ + µh + δh)Qh)h+ σ4Qhp∆Bp + 0.5σ2
4Qhp((∆Bp)

2 − h),

Rhp+1 = Rhp + (ρIh + τQh − µhRh)h+ σ5Rhp∆Bp + 0.5σ2
5Rhp((∆Bp)

2 − h)

Srp+1 = Srp + (Λr −
β3SrIr
Nr

− µrSr)h+ σ6Srp∆Bp + 0.5σ2
6Srp((∆Bp)

2 − h),

Erp+1 = Erp + (
β3SrIr
Nr

− (µr + γ3)Er)h+ σ7Erp∆Bp + 0.5σ2
7Erp((∆Bp)

2 − h),

Irp+1 = Irp + (γ3Er − (µr + δr)Ir)h+ σ8Irp∆Bp + 0.5σ2
8Irp((∆Bp)

2 − h).

(38)

5.3. Numerical Scheme with Singular Kernel266

To approximate the models (18)-(20). We will use GL-NSFDM [27] to solve (18) and267

(19). GL-NSFDM was selected for solving these system due to the following considerations:268

• Consistency with fractional operators: The GL method provides a direct discretiza-269

tion of fractional derivatives and aligns well with the singular memory effects inherent270

in Caputo and Riemann-Liouville derivatives.271

• Accuracy in long-term dynamics: Since monkeypox exhibits extended incubation272

and immunity effects, the GL method effectively captures the long-memory charac-273

teristics essential for modeling real-world epidemic data.274

• Stability considerations: While the GL method is conditionally stable, we ensured275

numerical stability by selecting an appropriate step size h and confirming that the276

fractional-order system satisfied the necessary stability conditions.277

Also, we use Milstein method (11) to approximate (20). For the stochastic version278

of the model, we employed the Milstein method to simulate the influence of random279

fluctuations in disease transmission. The key reasons for this choice include:280

• Higher-order accuracy: Compared to the Euler-Maruyama method, the Milstein281

scheme incorporates additional correction terms that improve accuracy, which is282

crucial when modeling stochastic perturbations in disease dynamics.283

• Stability in stochastic systems: Monkeypox outbreaks are subject to random environ-284

mental and demographic variations. The Milstein method provides better stability285

properties in capturing these fluctuations compared to lower-order stochastic solvers.286

• Numerical Robustness: Given that stochastic differential equations (SDEs) can ex-287

hibit high sensitivity to parameter changes, the Milstein method’s ability to handle288

noise-driven processes with higher precision ensures more reliable simulation results.289



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 22 of 41

Consider the general form of variable order fractional Caputo differential equation as290

follows:291

C
0 D

κ
t Y (t) = f(t, Y (t)), Y (0) = Y0, κ = α(t) ∈ (0, 1]. (39)

To approximate (39) using GL-NSFDM [27] as follows:292

f(tp, Yp) =
1

hκ
(Yp+1 −

p+1∑
i=1

wiYp+1−i − qp+1Y0). (40)

By applying (40) for system (12) in (0, t1], the following is the explicit solution:293

Shp+1 −
p+1∑
i=1

wiShp+1−i − Sh0qp+1 = hκ(Λκ
h −

(βκ
1 Irp + βκ

2 Ihp)Shp+1

Nhp

− µκ
hShp+1 + ϕκQhp),

Ehp+1 −
p+1∑
i=1

wiEhp+1−i − qp+1Eh0 = hκ(
(βκ

1 Irp + βκ
2 Ihp)Shp+1

Nhp

− (γκ1 + γκ2 + µκ
h)Ehp+1),

Ihp+1 −
p+1∑
i=1

wiIhp+1−i − qp+1Ih0 = hκ(γκ1Ehp+1 − (µκ
h + δκh + ρκ)Ihp+1),

Qhp+1 −
p+1∑
i=1

wiQhp+1−i − qp+1Qh0 = hκ(γκ2Ehp+1 − (ϕκ + τκ + µκ
h + δκh)Qhp+1),

Rhp+1 −
p+1∑
i=1

wiRhp+1−i − qp+1Rh0 = hκ(ρκIhp+1 + τκQhp+1 − µκ
hRhp+1),

Srp+1 −
p+1∑
i=1

wiSrp+1−i − qp+1Sr0 = hκ(Λκ
r −

βκ
3Srp+1Irp

Nrp
− µκ

rSrp+1),

Erp+1 −
p+1∑
i=1

wiErp+1−i − qp+1Er0 = hκ(
βκ
3Srp+1Irp

Nrp
− (µκ

r + γκ3 )Erp+1),

Irp+1 −
p+1∑
i=1

wiIrp+1−i − qp+1Ir0 = hκ(γκ3Erp+1 − (µκ
r + δκr )Irp+1). (41)

The explicit expressions may be derived using a comprehensive fundamental calculation.294

Shp+1 =
1

hκ(
(βκ

1 Ir+βκ
2 Ih)

Nh
+ µκ

h) + 1
[hκ(Λκ

h + ϕκQhp+1)

+

p+1∑
i=1

wiShp+1−i + qp+1Sh0

]
,



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 23 of 41

Ehp+1 =
1

1 + hκ(γκ1 + γκ2 + µκ
h)

[
βκ
1 Irp + βκ

2 Ihp)Shp+1

Nhp

+

p+1∑
i=1

wiEhp+1−i + qp+1Eh0

]
,

Ihp+1 =
1

1 + hκ(µκ
h + δκh + ρκ)

[
hκγκ1Ehp+1 +

p+1∑
i=1

wiIhp+1−i + qp+1Ih0

]
,

Qhp+1 =
1

1 + hκ(ϕκ + τκ + µκ
h + δκh)

[hκγκ2Ehp+1

+

p+1∑
i=1

wiQhp+1−i + qp+1Qh0

]
,

Rhp+1 =
1

1 + hκµκ
h

[
hκ(ρκIhp+1 + τκQhp+1) +

p+1∑
i=1

wiRhp+1−i + qp+1Rh0

]
,

Srp+1 =
1

1 + hκ(µκ
r +

βκ
3 Srp+1Irp

Nrp
)

[
hκΛκ

r +

p+1∑
i=1

wiSrp+1−i + qp+1Sr0

]
,

Erp+1 =
1

1 + hκ(µκ
r + γκ3 )

[
hκ

βκ
3Srp+1Irp

Nrp
+

p+1∑
i=1

wiErp+1−i + qp+1Er0

]
,

Irp+1 =
1

1 + hκ(µκ
r + δκr )

[
hκγκ3Erp+1 +

p+1∑
i=1

wiIrp+1−i + qp+1Ir0

]
. (42)

295

Consider the general formula for fractal-fractional Caputo defferential equation in the296

(t1, t2] given as follows:297

FFC
0 Dα,β

t Y (t) = f(t, Y (t)), Y (0) = Y0, α, β ∈ (0, 1]. (43)

Using the relation between fractal-fractional derivative and fractional order derivative [21]298

we have:299

tβ−1βf(tp, Yp) =
1

hα
(Yp+1 −

p+1∑
i=1

wiYp+1−i − qp+1Y0). (44)

Applying the recommended numerical approach (44) to approximate the system (19) yields300

the following results:301

1

hα
(Shp+1 −

p+1∑
i=1

wiShp+1−i − qp+1Sh0) = (βtβ−1)

[
Λα
h −

(βα
1 Irp + βα

2 Ihp)Shp+1

Nhp

−µα
hShp+1 + ϕαQhp] ,



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 24 of 41

1

hα
(Ehp+1 −

p+1∑
i=1

wiEhp+1−i − qp+1Eh0) = (βtβ−1)

[
(βα

1 Irp + βα
2 Ihp)Shp+1

Nhp

−(γα1 + γα2 + µα
h)Ehp+1] ,

1

hα
(Ihp+1 −

p+1∑
i=1

wiIhp−i+1 − Ih0qp+1) = (βtβ−1)(γα1Ehp+1 − (µα
h + δαh + ρα)Ihp+1),

1

hα
(Qhp+1 −

p+1∑
i=1

wiQhp−i+1 −Qh0qp+1) = (βtβ−1)(γα2Ehp+1 − (ϕα + τα + µα
h + δαh )Qhp+1),

1

hα
(Rhp+1 −

p+1∑
i=1

wiRhp−i+1 −Rh0qp+1) = (βtβ−1)(ραIhp+1 + ταQhp+1 − µα
hRhp+1),

1

hα
(Srp+1 −

p+1∑
i=1

wiSrp−i+1 − Sr0qp+1) = (βtβ−1)(Λα
r − βα

3 Srp+1Irp
Nrp

− µα
r Srp+1),

1

hα
(Erp+1 −

p+1∑
i=1

wiErp−i+1 − Er0qp+1) = (βtβ−1)(
βα
3 Srp+1Irp

Nrp
− (µα

r + γα3 )Erp+1),

1

hα
(Irp+1 −

p+1∑
i=1

wiIrp−i+1 − Ir0qp+1) = (βtβ−1)(γα3Erp+1 − (µα
r + δαr )Irp+1). (45)

To obtain the explicit expressions, a comprehensive calculation can be carried out.302

Shp+1 =
1

1 + βtβ−1hα(
(βα

1 Ir+βα
2 Ih)

Nh
+ µα

h)

[
βtβ−1hα(Λα

h + ϕαQhp+1)

+

p+1∑
i=1

wiShp+1−i + qp+1Sh0

]
,

Ehp+1 =
1

1 + βtβ−1hα(γα1 + γα2 + µα
h)

[
βtβ−1hα

(βα
1 Irp + βα

2 Ihp)Shp+1

Nhp

+

p+1∑
i=1

wiEhp+1−i + qp+1Eh0

]
,

Ihp+1 =
1

1 + βtβ−1hα(µα
h + δαh + ρα)

[
βtβ−1hαγα1Ehp+1 +

p+1∑
i=1

wiIhp+1−i + qp+1Ih0

]
,

Qhp+1 =
1

1 + βtβ−1hα(ϕα + τα + µα
h + δαh )

[
βtβ−1hαγα2Ehp+1

+

p+1∑
i=1

wiQhp+1−i + qp+1Qh0

]
,



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 25 of 41

Rhp+1 =
1

1 + βtβ−1hαµα
h

[
βtβ−1hα(ραIhp+1 + ταQhp+1) +

p+1∑
i=1

wiRhp+1−i + qp+1Rh0

]
,

Srp+1 =
1

1 + βtβ−1hα(µα
r +

βα
3 Srp+1Irp

Nrp
)

[
βtβ−1hαΛα

r +

p+1∑
i=1

wiSrp+1−i + qp+1Sr0

]
,

Erp+1 =
1

1 + βtβ−1hα(µα
r + γα3 )

[
βtβ−1hα

βα
3 Srp+1Irp

Nrp
+

p+1∑
i=1

wiErp+1−i + qp+1Er0

]
,

Irp+1 =
1

1 + βtβ−1hα(µα
r + δαr )

[
βtβ−1hαγα3Erp+1 +

p+1∑
i=1

wiIrp+1−i + qp+1Ir0

]
. (46)

In the third interval (t2, Tf ]), to approximate the system (20), we will use Milstein303

method (11). Then the explicit solution given as follows:304

Shp+1 = Shp + (Λh −
(β1Ir + β2Ih)Sh

Nh
− µhSh + ϕQh)h+ σ1Shp∆Bp + 0.5σ2

1Shp

((∆Bp)
2 − h),

Ehp+1 = Ehp + (
(β1Ir + β2Ih)Sh

Nh
− (γ1 + γ2 + µh)Eh)h+ σ2Ehp∆Bp + 0.5σ2

2Ehp

((∆Bp)
2 − h),

Ihp+1 = Ihp + (γ1Eh − (µh + δh + ρ)Ih)h+ σ3Ihp∆Bp + 0.5σ2
3Ihp((∆Bp)

2 − h),

Qhp+1 = Qhp + (γ2Eh − (ϕ+ τ + µh + δh)Qh)h+ σ4Qhp∆Bp + 0.5σ2
4Qhp((∆Bp)

2 − h),

Rhp+1 = Rhp + (ρIh + τQh − µhRh)h+ σ5Rhp∆Bp + 0.5σ2
5Rhp((∆Bp)

2 − h),

Srp+1 = Srp + (Λr −
β3SrIr
Nr

− µrSr)h+ σ6Srp∆Bp + 0.5σ2
6Srp((∆Bp)

2 − h),

Erp+1 = Erp + (
β3SrIr
Nr

− (µr + γ3)Er)h+ σ7Erp∆Bp + 0.5σ2
7Erp((∆Bp)

2 − h),

Irp+1 = Irp + (γ3Er − (µr + δr)Ir)h+ σ8Irp∆Bp + 0.5σ2
8Irp((∆Bp)

2 − h).

(47)

305

Remark 5.3. Regarding the stability for this method, you can find it in [27].306

6. Numerical Simulations307

This part aims to verify the provided experimental results as well as the analytical308

expressions that were produced in the preceding sections. From t=0 on June 13 to t=95309

on September 16, using the infected cases, 2022 (96 data points) [25]. Sh0 = 10000, Eh0 =310

50, Ih0 = 5.86, Qh0 = 1.14, Rh0 = 0, Sr0 = 1000, Er0 = 100, and Ir0 = 10 are the model’s311



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 26 of 41

starting values. Table 2 lists the values of this model’s parameters.312

Multiple behaviours are displayed by the dynamics.This enables us to analyse and predict313

how the disease will develop from the beginning to the end, giving us the chance to see a314

variety of behaviours from crossover to stochastic processes. We got the results of infection315

cases from utilizing the techniques (27)-(31), (34)-(38), and (42)-(47) for the provided316

crossover models (12)-(14), (15)-(17), and (18)-(20), which are compared with real-life317

information from the United States (June 13, 2022, to September 16, 2022). Figures318

1-4 show the numerical results of the first crossover (12)-(14) model with non-singular319

kernel with different values of α, β, and κ. In Figures 1-2, we contrasted real data with the320

outcomes of infected humans that were produced using the suggested model (12)-(14). We321

noted the good results we have at α = 0.97, 0.99, β = 0.99, and κ = 0.98−0.001sin(t/10)2,322

and α = 0.93, 0.97, β = 0.99, and κ = 0.98− 0.001cos(t/50). Figure 3 describes the effect323

of changing α in behavior of the solution. Also, Figure 4 shows the solution behavior for324

the considered model at different values of α, β and κ = 0.94− 0.01t.325

Figures 5-8 show the numerical results of the first crossover (15)-(17) model with non-326

singular kernel and distinct α, β, and κ values. Also, we observed from Figures 7-8, the327

good results appear when α = 0.93, 0.97, β = 0.99, and κ = 0.98 − 0.001sin(t/10)2, and328

α = 0.93, β = 0.99, with κ = 0.95− 0.001cos(t/50).329

Figures 9-12 show the numerical results of the first crossover (18)-(20) model with330

singular kernel and different values of α, β, and κ. Also, we noted that from Figures 11-12,331

the good results appear when α = 0.93, 0.95, β = 0.99, and κ = 0.98− 0.001sin(t/10)2.332

Table 3 shows the CPU time for three schemes (27)-(31), (34)-(38), and (42)-(47), we333

noted that the scheme (34)-(38) with singular kernel is the fastest one.334

We concluded that from the comparesion with real data the best result we have from335

the first crossover model (14)-(12) with a non-singular kernel. Moreover if we compare336

our result with [28], we have excellent results.337

Table 3: The CPU time at Tf = 100, α = 0.95, β = 0.99, κ = 0.94− 0.01t.

h CM1 CM2 CM3

1 0.0348147
0.2 0.8495 0.155779 0.1417802
0.1 3.3487 0.549053 0.2576076
0.06 7.5111 1.085631 0.408594



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 27 of 41

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.99

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.97

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.95

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.93

Real data

Figure 1: Numerical simulation of the piecewise system (12)-(14) at various α, β = 0.99 and κ = 0.98 −
0.001sin2(t/10)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 28 of 41

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.99

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.97

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.95

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.93

Real data

Figure 2: Numerical simulation of the piecewise system (12)-(14) at various α, β = 0.99 and κ = 0.95 −
0.001cos(t/50)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 29 of 41

0 20 40 60 80 100
0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

t

S
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

500

1000

1500

2000

2500

3000

3500

t

E
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

100

200

300

400

500

600

700

t

I h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
−2

−1.5

−1

−0.5

0

0.5

1

1.5

t

Q
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

500

1000

1500

2000

2500

3000

3500

4000

4500

t

R
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

100

200

300

400

500

600

700

800

900

1000

t

S
r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

100

200

300

400

500

600

700

t

E
r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

200

400

600

800

1000

1200

t

I r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

Figure 3: Numerical simulation of the piecewise system (12)-(14) at various values α , β=0.99, κ=0.94-0.01t,
σ1=0.01, σ2 = 0.05, = 0.02, σ3= 0.01, σ4 = 0.05, σ5 = 0.02, σ6 = 0.01, σ7 = 0.05 and σ8 = 0.02



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 30 of 41

0 10 20 30 40 50 60 70 80 90 100
0

2000

4000

6000

8000

10000

12000

t

S
h

 

 

α=1,β=0.99

α=0.90,β=0.97

α=0.80,β=0.95

α=0.75,β=0.93

0 20 40 60 80 100
0

500

1000

1500

2000

2500

3000

3500

t

E
h

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

0 20 40 60 80 100
−2

−1.5

−1

−0.5

0

0.5

1

1.5

t

Q
h

0 20 40 60 80 100
0

1000

2000

3000

4000

5000

6000

t

R
h

0 20 40 60 80 100
0

100

200

300

400

500

600

700

800

900

1000

t

S
r

0 20 40 60 80 100
100

200

300

400

500

600

700

t

E
r

0 20 40 60 80 100
0

200

400

600

800

1000

1200

t

I r

Figure 4: Numerical simulation of the piecewise system (12)-(14) at various values α , β, κ=0.94-0.01t,
σ1 = 0.01, σ2 = 0.05, = 0.02, σ3= 0.01, σ4 = 0.05, σ5 = 0.02, σ6 = 0.01, σ7 = 0.05 and σ8 = 0.02



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 31 of 41

0 20 40 60 80 100
0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

t

S
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

500

1000

1500

2000

2500

3000

3500

t

E
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

100

200

300

400

500

600

700

t

I h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

0.2

0.4

0.6

0.8

1

1.2

1.4

t

Q
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

1000

2000

3000

4000

5000

6000

7000

t

R
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

100

200

300

400

500

600

700

800

900

1000

t

S
r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

100

200

300

400

500

600

700

t

E
r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

200

400

600

800

1000

1200

1400

1600

t

I r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

Figure 5: Numerical simulation of the piecewise system (15)-(17) at various values α , β=0.99, κ=0.94-0.01t,
σ1=0.01, σ2 = 0.05, = 0.02, σ3= 0.01, σ4 = 0.05, σ5 = 0.02, σ6 = 0.01, σ7 = 0.05 and σ8 = 0.02



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 32 of 41

0 10 20 30 40 50 60 70 80 90 100
0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

t

S
h

 

 

α=1,β=0.98

α=0.90,β=0.97

α=0.80,β=0.95

α=0.75,β=0.93

0 20 40 60 80 100
0

500

1000

1500

2000

2500

3000

3500

t

E
h

0 20 40 60 80 100
0

100

200

300

400

500

600

700

t

I h

0 20 40 60 80 100
0

0.2

0.4

0.6

0.8

1

1.2

1.4

t

Q
h

0 20 40 60 80 100
0

1000

2000

3000

4000

5000

6000

t

R
h

0 20 40 60 80 100
0

100

200

300

400

500

600

700

800

900

1000

t

S
r

0 20 40 60 80 100
0

100

200

300

400

500

600

700

t

E
r

0 20 40 60 80 100
0

200

400

600

800

1000

1200

1400

t

I r

Figure 6: Numerical simulation of the piecewise system (15)-(17 ) at various values α , β, κ=0.94-0.01t,
σ1 = 0.01, σ2 = 0.05, = 0.02, σ3= 0.01, σ4 = 0.05, σ5 = 0.02, σ6 = 0.01, σ7 = 0.05 and σ8 = 0.02



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 33 of 41

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.99

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.97

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.95

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.93

Real data

Figure 7: Numerical simulation of the piecewise system (15)-(17) numerical simulation at various α, β = 0.99
and κ = 0.98− 0.001sin2(t/10)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 34 of 41

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.99

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.97

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.95

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.93

Real data

Figure 8: Numerical simulation of the piecewise system (15)-(17) at various values α, β = 0.99 and κ =
0.95− 0.001cos(t/50)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 35 of 41

0 20 40 60 80 100
0

0.5

1

1.5

2

2.5

3

3.5
x 10

4

t

S
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

2000

4000

6000

8000

10000

12000

t

E
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

500

1000

1500

2000

2500

t

I h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

0.2

0.4

0.6

0.8

1

1.2

1.4

t

Q
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

2000

4000

6000

8000

10000

12000

14000

16000

18000

t

R
h

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

500

1000

1500

2000

2500

t

S
r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

200

400

600

800

1000

1200

1400

1600

t

E
r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

0 20 40 60 80 100
0

500

1000

1500

2000

2500

3000

3500

4000

4500

t

I r

 

 

α=1,α(t)=0.94−0.01t

α=0.90,α(t)=0.94−0.01t

α=0.80,α(t)=0.94−0.01t

α=0.75,α(t)=0.94−0.01t

Figure 9: Numerical simulation of the piecewise system (18)-(20) at various values α , β=0.99, κ=0.94-0.01t,
σ1=0.01, σ2 = 0.05, = 0.02, σ3= 0.01, σ4 = 0.05, σ5 = 0.02, σ6 = 0.01, σ7 = 0.05 and σ8 = 0.02



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 36 of 41

0 10 20 30 40 50 60 70 80 90 100
0

0.5

1

1.5

2

2.5

3

3.5

4
x 10

4

t

S
h

 

 

α=1,β=0.98

α=0.90,β=0.97

α=0.80,β=0.95

α=0.75,β=0.93

0 20 40 60 80 100
0

2000

4000

6000

8000

10000

12000

14000

t

E
h

0 10 20 30 40 50 60 70 80 90 100
0

200

400

600

800

1000

1200

1400

1600

1800

2000

t

I h

0 20 40 60 80 100
0

0.2

0.4

0.6

0.8

1

1.2

1.4

t

Q
h

0 10 20 30 40 50 60 70 80 90 100
0

0.5

1

1.5

2

2.5
x 10

4

t

R
h

0 20 40 60 80 100 120
0

500

1000

1500

2000

2500

t

S
r

0 20 40 60 80 100
0

500

1000

1500

2000

t

E
r

0 20 40 60 80 100
0

1000

2000

3000

4000

5000

6000

t

I r

Figure 10: Numerical simulation of the piecewise system (18)-(20) at various values α , β, κ=0.94-0.01t,
σ1 = 0.01, σ2 = 0.05, = 0.02, σ3= 0.01, σ4 = 0.05, σ5 = 0.02, σ6 = 0.01, σ7 = 0.05 and σ8 = 0.02



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 37 of 41

0 10 20 30 40 50 60 70 80 90 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.99

Real data

0 10 20 30 40 50 60 70 80 90 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.97

Real data

0 10 20 30 40 50 60 70 80 90 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.95

Real data

0 10 20 30 40 50 60 70 80 90 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.93

Real data

Figure 11: Numerical simulation of the piecewise system (18)-(20) at various values α, β = 0.99 and κ =
0.98− 0.001sin2(t/10)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 38 of 41

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.99

Real data

0 10 20 30 40 50 60 70 80 90 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.97

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.95

Real data

0 20 40 60 80 100
0

100

200

300

400

500

600

t

I h

 

 

α=0.93

Real data

Figure 12: Numerical simulation of the piecewise system (18)-(20) at various values α, β = 0.99 and κ =
0.95− 0.001cos(t/50)



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 39 of 41

7. Conclusions338

This study presents three crossover models on Monkeypox disease which combine Ca-339

puto, Mittag-Leffler, and the Caputo-Fabrizio definitions. A combination of three models340

variable order fractional, fractal-fractional, and stochastic with their piecewise derivatives341

to describe a Monkypox disease are presented in three intervals of time. The nonstandard342

Grünwald−Letnikov finite difference method was used to approximate the deterministic343

models with a singular kernel and the Toufik-Atangana technique to approximate the344

deterministic model involves a nonsingular Mittag-Leffler kernel. Moreover, the approxi-345

mation of the integral Caputo-Fabrizio and Lagrange polynomials of two steps was used346

to approximate the deterministic model with a nonsingular exponential decay kernel. The347

Milstein method was implemented to approximate the stochastic differential equation.348

Analysis was done on the suggested model’s stability.349

Regarding the impact of these kernels on model outcomes:350

• Singular kernels (e.g., Caputo fractional derivative) introduce power-law memory351

effects, meaning past states significantly influence present disease dynamics. This352

results in a slower decay of past influences, prolonged epidemic waves, and stronger353

hysteresis effects in the system. As a consequence, the disease dynamics exhibit a354

more persistent infection tail and slower stabilization.355

• Non-singular kernels (e.g., Atangana-Baleanu and Caputo-Fabrizio fractional deriva-356

tives) employ exponentially decaying memory effects, where the impact of past states357

diminishes more rapidly. This leads to faster stabilization of the disease, shorter358

memory retention, and more immediate responses to changes in transmission dy-359

namics.360

From the numerical comparison between these three operators, we found that CM1 is361

the fastest one. Also, the GL-NSFDM is more efficient than others because it can work362

with big step sizes. Moreover, from comparison the results of three crossover model with363

real data, we found that the first crossover model with non singular kernel and ABC364

operator is the best one to describe the monkeybox model.365

The piecewise crossover differential equations, constructed with fractional and variable366

order operators alongside stochastic derivatives, have opened new avenues for researchers367

across various fields, enabling them to capture diverse behaviors over time intervals. Ap-368

plying this approach to real-world problems allow researchers to more accurately reflect369

reality.370

References371

[1] Centers for Disease Control and Prevention. U.S. Monkeypox case trends reported to372

CDC. https://www.cdc.gov/poxvirus/monkeypox/response/2022/mpx-trends.html,373

2022.374



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 40 of 41

[2] A. A. Aljabali, M. A. Obeid, M. B. Nusair, A. Hmedat, and M. M. Tambuwala.375

Monkeypox virus: an emerging epidemic. Microbial Pathogenesis, 173(Pt A):105794,376

2022.377

[3] J. G. Breman. Monkeypox: an emerging infection for humans. In Emerging Viruses378

in Human Populations, pages 45–67. John Wiley & Sons, Ltd, Chichester, 2000.379

[4] R. S. Levine, A. T. Peterson, K. L. Yorita, D. Carroll, I. K. Damon, and M. G.380

Reynolds. Ecological niche and geographic distribution of human monkeypox in381

Africa. PLoS ONE, 2(1):e176, 2007.382

[5] S. Usman and I. Isa Adamu. Modeling the transmission dynamics of the monkeypox383

virus infection with treatment and vaccination interventions. Journal of Applied384

Mathematics and Physics, 5(12):2335–2353, 2017.385

[6] A. Khan, Y. Sabbar, and A. Din. Stochastic modeling of the Monkeypox 2022 epi-386

demic with cross-infection hypothesis in a highly disturbed environment. Mathemat-387

ical Biosciences and Engineering, 19(12):13560–13581, 2022.388

[7] R. Grant, L.-L. Nguyen, and R. Breban. Modelling human-to-human transmission of389

monkeypox. Bulletin of the World Health Organization, 98(9):638–640, 2020.390

[8] S. V. Bankuru, S. Kossol, W. Hou, P. Mahmoudi, J. Rychtář, and D. Taylor. A391

game-theoretic model of Monkeypox to assess vaccination strategies. PeerJ, 8:e9272,392

2020.393

[9] N. H. Sweilam, S. M. AL-Mekhlafi, and A. Almutairi. Fractal fractional optimal394

control for a novel malaria mathematical model; a numerical approach. Results in395

Physics, 19:103446, 2020.396

[10] M. Y. Sahnoune, A. Ez-zetouni, K. Akdim, et al. Qualitative analysis of a fractional-397

order two-strain epidemic model with vaccination and general non-monotonic inci-398

dence rate. International Journal of Dynamics and Control, 11(4):1532–1543, 2023.399

[11] C. Xu, C. Aouiti, M. Liao, P. Li, and Z. Liu. Chaos control strategy for a fractional-400

order financial model. Advances in Difference Equations, 2020:573, 2020.401

[12] K. N. Nabi, H. Abboubakar, and P. Kumar. Forecasting of COVID-19 pandemic: from402

integer derivatives to fractional derivatives. Chaos, Solitons & Fractals, 141:110283,403

2020.404

[13] S. K. Biswas, U. Ghosh, and S. Sarkar. Mathematical model of Zika virus dynamics405

with vector control and sensitivity analysis. Infectious Disease Modelling, 5:23–41,406

2020.407

[14] S. Hasan, A. El-Ajou, S. Hadid, M. Al-Smadi, and S. Momani. Atangana-Baleanu408

fractional framework of reproducing kernel technique in solving fractional population409

dynamics system. Chaos, Solitons & Fractals, 133:109624, 2020.410

[15] R. P. Chauhan, S. Kumar, B. S. T. Alkahtani, and S. S. Alzaid. A study on frac-411

tional order financial model by using Caputo-Fabrizio derivative. Results in Physics,412

57:107335, 2024.413

[16] N. H. Sweilam, F. A. Rihan, and S. M. AL-Mekhlafi. A fractional-order delay dif-414

ferential model with optimal control for cancer treatment based on synergy between415

anti-angiogenic and immune cell therapies. Discrete and Continuous Dynamical Sys-416

tems - Series S, 13(9):2403–2424, 2020.417



N. H. Sweilam et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 5913 41 of 41

[17] J. Danane, Z. Hammouch, K. Allali, S. Rashid, and J. Singh. A fractional-order model418

of coronavirus disease 2019 (COVID-19) with governmental action and individual419

reaction. Mathematical Methods in the Applied Sciences, 46(7):7808–7824, 2023.420

[18] M. Zamir, F. Nadeem, T. Abdeljawad, and Z. Hammouch. A fractional multi-order421

model to predict the COVID-19 outbreak in Morocco. Applied and Computational422

Mathematics, 20(1):177–203, 2021.423

[19] K. Shah, T. Abdeljawad, and A. Ali. Mathematical analysis of the Cauchy type dy-424

namical system under piecewise equations with Caputo fractional derivative. Chaos,425

Solitons & Fractals, 161:112356, 2022.426

[20] X. P. Li, H. F. Alrihieli, E. A. Algehyne, M. A. Khan, M. Y. Alshahrani, Y. Alraey,427

and M. B. Riaz. Application of piecewise fractional differential equation to COVID-19428

infection dynamics. Results in Physics, 39:105685, 2022.429

[21] A. Atangana and S. İğret Araz. New Numerical Scheme with Newton Polynomial:430

Theory, Methods, and Applications. Academic Press, London, 2021.431

[22] S. Etemad, I. Avci, P. Kumar, D. Baleanu, and S. Rezapour. Some novel mathematical432

analysis on the fractal–fractional model of the AH1N1/09 virus and its generalized433

Caputo-type version. Chaos, Solitons & Fractals, 162:112511, 2022.434

[23] P. E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations,435

volume 23 of Stochastic Modelling and Applied Probability. Springer, Berlin, 1992.436

[24] G. N. Mil’shtejn. Approximate integration of stochastic differential equations. Theory437

of Probability and Its Applications, 19(3):557–562, 1975.438

[25] O. J. Peter, S. Kumar, N. Kumari, et al. Transmission dynamics of Monkeypox439

virus: a mathematical modelling approach. Modeling Earth Systems and Environ-440

ment, 8(3):3423–3434, 2022.441

[26] O. J. Peter, S. Kumar, N. Kumari, F. A. Oguntolu, K. Oshinubi, and R. Musa.442

Transmission dynamics of Monkeypox virus: a mathematical modelling approach.443

Modeling Earth Systems and Environment, 8(3):3423–3434, 2022.444

[27] A. J. Arenas, G. González-Parra, and B. M. Chen-Charpentier. Construction of445

nonstandard finite difference schemes for the SI and SIR epidemic models of fractional446

order. Mathematics and Computers in Simulation, 121:48–63, 2016.447

[28] P. Kumar, M. Vellappandi, Z. Khan, S. M. Sivalingam, A. Kaziboni, and V. Govin-448

daraj. A case study of monkeypox disease in the United States using mathematical449

modeling with real data. Mathematics and Computers in Simulation, 213:444–465,450

2023.451