Modeling the monkeypox infection using the Mittag–Leffler kernel

Muhammad Altaf Khan; Mutum Zico Meetei; Kamal Shah; Thabet Abdeljawad; Mohammad Y. Alshahrani

doi:10.1515/phys-2023-0111

Article Open Access

Modeling the monkeypox infection using the Mittag–Leffler kernel

Muhammad Altaf Khan , Mutum Zico Meetei , Kamal Shah , Thabet Abdeljawad and Mohammad Y. Alshahrani

Published/Copyright: September 20, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 21 Issue 1

Abstract

This article presents the mathematical formulation for the monkeypox infection using the Mittag–Leffler kernel. A detailed mathematical formulation of the fractional-order Atangana-Baleanu derivative is given. The existence and uniqueness results of the fractional-order system is established. The local asymptotical stability for the disease-free case, when ℛ 0 < 1 , is given. The global asymptotical stability is given when ℛ 0 > 1 . The backward bifurcation analysis for fractional system is shown. The authors give a numerical scheme, solve the model, and present the results graphically. Some graphical results are shown for disease curtailing in the USA.

Keywords: Monkeypox outbreak USA; fractional-order model; local and global stability; novel numerical scheme

1 Introduction

Monkeypox is a disease that is transmissible to the human population through animals although it is clinically less severe than smallpox. It displays signs that resemble those of smallpox. Monkeypox has replaced smallpox as the most notable Orthopoxvirus to public health since smallpox was eliminated in 1980 and smallpox vaccination was subsequently discontinued. Monkeypox, which mostly affects central and west Africa, has begun to enter towns and is regularly observed near tropical rainforests. Animals live on non-human primates and several rodent species as hosts. The monkeypox virus is able to infect several animal species, such as rope, primate, and tree squirrels. To understand the virus history in more detail, some more research work is needed. In addition, it is necessary to find out the reservoirs of the monkeypox virus and determine how it disseminates in the wild [1,2].

A serious condition that affects public health worldwide is monkeypox. The illness spread over the entire world in addition to central and western Africa. Contact with infected pet prairie dogs in the USA resulted in the world’s first monkeypox outbreak outside of Africa in 2003. These species had been kept with pouched rats and dormice brought from the Gambia. This outbreak spread monkeypox throughout the nation, resulting in over 70 cases. There have also been reports of monkeypox in September 2018 among people who came from Nigeria to Israel. Similarly, reports show evidence of the virus in the UK in the years 2018, 2019, 2021, and 2022. The infected cases in Singapore in 2019 and in the USA in 2021 have been recorded. The cases of monkeypox are also recorded in 2022 in many non-endemic countries [1].

Fractional calculus (FC) has got too much interest from researchers’ point of view due to its various prosperities, such as memory, heredity, crossover behavior, and many more. FC has been implied in scientific problems and many disease models in the literature and found suitable for disease dynamics due to its characteristics. For example, Guo and Li, in their study [3], used the FC to obtain results for the people who involve in online game and become addicted to it. The role of vaccination in the disease control has been studied using fractional-order derivative in the study by Baba et al. [4]. The discrete fractional derivative to study the COVID-19 epidemic has been used in the study by Abbes et al. [5]. Asamoah et al. [6] the authors formulated the listeriosis infection in fractional derivative. The cholera infection has been studied by George et al. [7] using the fractional-order system. To understand the liver disease dynamics, the authors formulated them in fractional derivative [8] and obtained their dynamical results. A single-route transmission model in fractional derivative is proposed in the study by Okyere and Ackora-Prah [9]. The HIV/AIDS disease under the fractional-order derivative is considered in the study by Farman et al. [10]. The Hepatitis C infection model using a non-singular kernel is suggested in the study by Evirgen et al. [11]. The rubella disease model using a fractional model based on Mittag–Leffler kernel is considered in the study by Koca [12]. The authors considered the fractional model for COVID-19 disease Bhatter et al. [13]. A mathematical model for diabetes using fractional derivative is considered in [14]. Other related work that used fractional operators in the study by Karaagac et al. the recent past are given in the study by Karaagac et al. [15–17].

In the literature, numerous models in terms of mathematics are reported to study the monkeypox disease. For instance, the researchers created a mathematical model for monkeypox viral transmission in the study by Madubueze et al. [18]. Somma et al. in their study [19] examined the mathematical modeling of the disease and obtained the results. Lasisi et al. [20] created a mathematical model to comprehend how the monkeypox virus spreads to people. Usman and Adamu [21] proposed to cure monkeypox and administer the vaccination using a compartmental mathematical model. Emeka et al. [22] examined the the dynamics of the monkeypox virus under incomplete vaccination. In the study by Peter et al. [23], a compartmental mathematical model is taken into account to comprehend the intricate nature of the illness process. The scientists looked into the transmission of the monkeypox virus from rodents to humans and from humans to humans, as well as a comprehensive examination of the illness. Allehiany et al. [24] have lately examined the dynamics of monkeypox under real data and their backward bifurcation. Section-wise details of this work are as follows: The description of the model is given in Section 2. The analysis of the arbitrary order model, its positivity and boundedness, as well as the existence of solutions and uniqueness, is given in Section 3. The equilibria of the model and their analysis are shown in Section 4. Numerical scheme and its application to the fractional-order monkeypox disease model are given in Section 5, while the findings are briefly discussed in Section 6.

2 Model construction

Allehiany et al. [24] considered the monkeypox disease in an integer-order derivative using the recent cases in the USA and obtained the dynamical results. In this work, we will consider the work in [24] by extending it to the fractional-order derivative in Atangana–Baleanu derivative. The population of the humans are divided into five, while rodents are divided into three. The groups that are involved in the human population including susceptible S h ( t ) , exposed E h ( t ) , people infected with monkeypox virus I h ( t ) , quarantined people, Q h ( t ) , and the people who get recovery from disease R h ( t ) . The total population of human represented by N ( t ) is determined as follows:

(1) N h ( t ) = S h ( t ) + E h ( t ) + I h ( t ) + Q h ( t ) + R h ( t ) .

The parameter Ψ h represents the healthy population’s recruitment rate, whereas ν h represents its natural death rate. At a rate ϖ 1 a person becomes ill with the virus after coming into contact with an infected rodent. Animal-to-human (zoonotic) transmission can occur by coming into contact with the blood, bodily fluids, cutaneous, or gastric wounds of infected animals. At a rate of ϖ 2 , the healthy people get infections while contacting the infected person. The route of transmission for ϖ 1 and ϖ 2 is shown by the force of infection as follows:

Γ ( t ) = ϖ 1 I r N r + ϖ 2 I h N h .

After the disease symptoms last 2–4 weeks, the person is identified as being infected with the virus, and hence, at a rate of η 1 , joins the infected class I h , while certain members of the people are quarantined at a rate of η 2 . The people recovered at the rates ψ and ϕ , respectively at infected and quarantined classes. People in the infected and quarantined class die from infection at a rate of ε 1 and ε 2 , respectively.

The animal population are divided into three groups, such as susceptible S r ( t ) , exposed E r ( t ) , and infected rodents I r ( t ) . The total population of rodents is shown by,

(2) N r ( t ) = S r ( t ) + E r ( t ) + I r ( t ) .

The parameters Ψ r and η r define the birth and natural death rate of the rodents population, respectively. The contact rate ϖ r by which a healthy rodent gets infected when it interacts with an infected rodent. The contact rate with the chance of a rodent developing an illness for each interaction with an infected rodent is represented by the parameter ϖ r . The force of infection of rodent population is shown as follows:

χ ( t ) = ϖ r I r N r .

The exposed animals become infected with a rate given by κ r . The following nonlinear system developed in the study by Allehiany et al. [24] based on the discussion above in integer order derivative is given as follows:

(3) d S h d t = Ψ h − Γ ( t ) S h − ν h S h , d E h d t = Γ ( t ) S h − ( η 1 + η 2 + ν h ) E h , d I h d t = η 1 E h − ( ψ + ν h + ε 1 ) I h , d Q h d t = η 2 E h − ( ϕ + ε 2 + ν h ) Q h , d R h d t = ψ I h + ϕ Q h − ν h R h , d S r d t = Ψ r − χ ( t ) S r − ν r S r , d E r d t = χ ( t ) S r − ( ν r + κ r ) E r , d I r d t = κ r E r − ν r I r ,

where

Γ ( t ) = ϖ 1 I r N r + ϖ 2 I h N h and χ ( t ) = ϖ r I r N r .

The initial conditions subject to system (3) are as follows:

(4) S h ( 0 ) ≥ 0 , E h ( 0 ) ≥ 0 , I h ( 0 ) ≥ 0 , Q h ( 0 ) ≥ 0 , R h ( 0 ) ≥ 0 , S r ( 0 ) ≥ 0 , E r ( 0 ) ≥ 0 , I r ( 0 ) ≥ 0 .

2.1 Model in fractional-order

We first present some related definitions here. For more details on the definition and its application, see [25].

Definition 1

For a function f , the Atangana–Baleanu derivative in Caputo sense is given as follows:

(5) D t σ ABC [ f ( t ) ] = G ( σ ) 1 − σ ∫ a 2 t f ′ ( x ) E σ [ ζ ( t − x ) σ ] d x ,

where ζ = − σ 1 − σ , σ ∈ [ 0 , 1 ] , a 2 > a 1 , f ∈ H 1 ( a 1 , a 2 ) , and G ( σ ) = 1 − σ + σ ∕ Γ ( σ ) .

The Laplace transform has been used to obtain the following integral for the Atangana–Baleanu derivative.

Definition 2

The integral related to Eq. (5) can be defined as follows:

(6) I t σ a 1 ABC [ f ( t ) ] = 1 − σ G ( σ ) f ( t ) + σ G ( σ ) Γ ( σ ) ∫ a 1 t f ( x ) ( t − x ) σ − 1 d x ,

where the function G ( σ ) is called the normalization function, and it holds for G ( 0 ) = 1 = G ( 1 ) .

The fractional-order system is superior to the integer-order systems, due to many properties such as the heredity memory effects and the crossover behavior that make it more significant compared to integer systems. One of the more interesting properties of the fractional-order system is that it provides reasonable fitting to the cases compared to the non-fractional system. With the above such advantages, we shall use the result given in Eq. (5) by applying it to our monkeypox infection model (3), and obtain the following fractional system:

(7) D t σ 0 ABC S h = Ψ h − Γ ( t ) S h − ν h S h , D t σ 0 ABC E h = Γ ( t ) S h − ( η 1 + η 2 + ν h ) E h , D t σ 0 ABC I h = η 1 E h − ( ψ + ν h + ε 1 ) I h , D t σ 0 ABC Q h = η 2 E h − ( ϕ + ε 2 + ν h ) Q h , D t σ 0 ABC R h = ψ I h + ϕ Q h − ν h R h , D t σ 0 ABC S r = Ψ r − χ ( t ) S r − ν r S r , D t σ 0 ABC E r = χ ( t ) S r − ( ν r + κ r ) E r , D t σ 0 ABC I r = κ r E r − ν r I r ,

where

Γ ( t ) = ϖ 1 I r N r + ϖ 2 I h N h , and χ ( t ) = ϖ r I r N r .

3 Model analysis

This section will explore the related mathematical properties associated with the system (3). We first demonstrate the model’s positivity and boundedness, then the existence and uniqueness solution of the system will be presented. The total population of human is

D t σ 0 ABC N h ( t ) = D t σ 0 ABC S h ( t ) + D t σ 0 ABC E h ( t ) + D t σ 0 ABC I h ( t ) + D t σ 0 ABC Q h ( t ) + D t σ 0 ABC R h ( t ) .

Furthermore, we obtain the following:

(8) D t σ 0 ABC N h ( t ) = Ψ h − ν h N h − ε 1 I h − ε 2 Q h , D t σ 0 ABC N h ( t ) ≤ Ψ h − ν h N h .

With the solution of Eq. (8) using the Laplace transform, we obtain

(9) N h ( t ) ≤ G ( σ ) G ( σ ) + ( 1 − σ ) ν h N h ( 0 ) + ( 1 − σ ) Ψ h G ( σ ) + ( 1 − σ ) ν h × E σ , 1 − σ ν h G ( σ ) + ( 1 − σ ) ν h t σ + σ Ψ h G ( σ ) + ( 1 − σ ) ν h E σ , σ + 1 × − σ ν h G ( σ ) + ( 1 − σ ) ν h t σ .

The asymptomatic nature of the Mittag–Leffler leads to the following, when t → ∞ :

(10) lim t → ∞ N h ( t ) ≤ Ψ h ν h .

In a similar way, we can obtain the total dynamics of the rodent population,

D t σ 0 ABC N r ( t ) = D t σ 0 ABC S r ( t ) + D t σ 0 ABC E r ( t ) + D t σ 0 ABC I r ( t ) ,

(11) D t σ 0 ABC N r = Ψ r − ν r N r .

Solution of Eq. (11) gives

(12) N r ( t ) ≤ G ( σ ) G ( σ ) + ( 1 − σ ) ν r N r ( 0 ) + ( 1 − σ ) Ψ r G ( σ ) + ( 1 − σ ) ν r × E σ , 1 − σ ν r G ( σ ) + ( 1 − σ ) ν r t σ + σ Ψ r G ( σ ) + ( 1 − σ ) ν r E σ , σ + 1 × − σ ν r G ( σ ) + ( 1 − σ ) ν r t σ .

It follows from Eq. (12) that

(13) lim t → ∞ N r ( t ) ≤ Ψ r ν r .

When t → ∞ , then Eqs (10) and (13) provide Ψ h ∕ ν h and Ψ r ∕ ν r , respectively. Thus, for any t ≥ 0 , the monkeypox infection model (7) possesses nonnegative solution. Therefore, any affiliated solutions to the model (7) shall remain positive for any t ≥ 0 . Thus, the system (7) is epidemiologically well-posed, and its dynamical characteristics may be explored in the following feasible region:

(14) Ω = Ω h × Ω r ⊆ R + 5 + R + 3 ,

where

Ω h = ϒ 1 ∈ R + 5 : ϒ 2 ≤ Ψ h ν h , Ω r = ϒ 3 ∈ R + 3 : ϒ 4 ≤ Ψ r ν r ,

where ϒ 1 = ( S h , E h , I h , Q h , R h ) , ϒ 2 = S h + E h + I h + Q h + R h , ϒ 3 = ( S r , E r , I r ) , and ϒ 4 = S r + E r + I r .

3.1 Boundedness and positive solution

We study the existence and uniqueness of the monkeypox infection model (7). To do this, we first made the following changes to the system (7):

(15) D t σ 0 ABC S h = L 1 ( t , W ) , D t σ 0 ABC E h = L 2 ( t , W ) , D t σ 0 ABC I h = L 3 ( t , W ) , D t σ 0 ABC Q h = L 4 ( t , W ) , D t σ 0 ABC R h = L 5 ( t , W ) , D t σ 0 ABC S r = L 6 ( t , W ) , D t σ 0 ABC E r = L 7 ( t , W ) , D t σ 0 ABC I r = L 8 ( t , W ) ,

where the kernels are represented as follows:

(16) L 1 ( t , W ) = Ψ h − Γ ( t ) S h − ν h S h , L 2 ( t , W ) = Γ ( t ) S h − ( η 1 + η 2 + ν h ) E h , L 3 ( t , W ) = η 1 E h − ( ψ + ν h + ε 1 ) I h , L 4 ( t , W ) = η 2 E h − ( ϕ + ε 2 + ν h ) Q h , L 5 ( t , W ) = ψ I h + ϕ Q h − ν h R h , L 6 ( t , W ) = Ψ r − χ ( t ) S r − ν r S r , L 7 ( t , W ) = χ ( t ) S r − ( ν r + κ r ) E r , L 8 ( t , W ) = κ r E r − ν r I r ,

and W = ( S h , E h , I h , Q h , R h , S r , E r , I r ) . Using the well-known fixed point theorem called the Banach fixed point theorem used as an application to the COVID-19 infection, see [26]. Applying the fractional integral to system (15), we have

(17) S h ( t ) − S h ( 0 ) = 1 − σ G ( σ ) L 1 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 1 ( k , W ) ( t − x ) σ − 1 d x , E h ( t ) − E h ( 0 ) = 1 − σ G ( σ ) L 2 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 2 ( k , W ) ( t − x ) σ − 1 d x , I h ( t ) − I h ( 0 ) = 1 − σ G ( σ ) L 3 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 3 ( k , W ) ( t − x ) σ − 1 d x , Q h ( t ) − Q h ( 0 ) = 1 − σ G ( σ ) L 4 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 4 ( k , W ) ( t − x ) σ − 1 d x , R h ( t ) − R h ( 0 ) = 1 − σ G ( σ ) L 5 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 5 ( k , W ) ( t − x ) σ − 1 d x , S r ( t ) − S r ( 0 ) = 1 − σ G ( σ ) L 6 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 6 ( k , W ) ( t − x ) σ − 1 d x , E r ( t ) − E r ( 0 ) = 1 − σ G ( σ ) L 7 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 7 ( k , W ) ( t − x ) σ − 1 d x , I r ( t ) − I r ( 0 ) = 1 − σ G ( σ ) L 8 ( t , W ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 8 ( k , W ) ( t − x ) σ − 1 d x .

Let the set B = U ( J ) × ⋯ × U ( J ) , where U ( J ) denotes the Banach space of real-valued continuous functions defined on an interval J = [ 0 , T ] , with the corresponding norm defined by ‖ ( S h , E h , I h , Q h , R h ) ‖ = ‖ S h ‖ + ‖ E h ‖ + ‖ I h ‖ + ‖ Q h ‖ + ‖ R h ‖ , and ‖ ( S r , E r , I r ) ‖ = ‖ S r ‖ + ‖ E r ‖ + ‖ I r ‖ , where

‖ S h ‖ = sup t ∈ J ∣ S h ( t ) ∣ = b 1 , ‖ E h ‖ = sup t ∈ J ∣ E h ( t ) ∣ = b 2 , ‖ I h ‖ = sup t ∈ J ∣ I h ( t ) ∣ = b 3 , ‖ Q h ‖ = sup t ∈ J ∣ Q h ( t ) ∣ = b 4 , ‖ R h ‖ = sup t ∈ J ∣ R h ( t ) ∣ = b 5 , ‖ S r ‖ = sup t ∈ J ∣ S r ( t ) ∣ = b 6 , ‖ E r ‖ = sup t ∈ J ∣ E r ( t ) ∣ = b 7 , ‖ I r ‖ = sup t ∈ J ∣ I r ( t ) ∣ = b 8 .

Theorem 1

(Lipschitz condition and contraction) For each of the kernels L 1 , … , L 8 in Eq. (15), there exists M i for i = 1 , 2 , … , 8 , such that

(18) ‖ L 1 ( t , S h ) − L 1 ( t , S h 1 ) ‖ ≤ M 1 ‖ S h − S h 1 ‖ , ‖ L 2 ( t , E h ) − L 2 ( t , E h 1 ) ‖ ≤ M 2 ‖ E h − E h 1 ‖ , ‖ L 3 ( t , I h ) − L 3 ( t , I h 1 ) ‖ ≤ M 3 ‖ I h − I h 1 ‖ , ‖ L 4 ( t , Q h ) − L 4 ( t , Q h 1 ) ‖ ≤ M 4 ‖ Q h − Q h 1 ‖ , ‖ L 5 ( t , R h ) − L 5 ( t , R h 1 ) ‖ ≤ M 5 ‖ R h − R h 1 ‖ , ‖ L 6 ( t , S r ) − L 6 ( t , S r 1 ) ‖ ≤ M 6 ‖ S r − S r 1 ‖ , ‖ L 7 ( t , E r ) − L 7 ( t , E r 1 ) ‖ ≤ M 7 ‖ E r − E r 1 ‖ , ‖ L 8 ( t , I r ) − L 8 ( t , I r 1 ) ‖ ≤ M 8 ‖ I r − I r 1 ‖

and are contractions for 0 < M i < 1 , i = 1 , 2 , … , 8 .

Proof

We give the result first for S h equation of model (7),

(19) ‖ L 1 ( t , S h ) − L 1 ( t , S h 1 ) ‖ = ‖ Ψ h − ϖ 1 I r N r + ϖ 2 I h N h S h − ν h S h − Ψ h + ϖ 1 I r N r + ϖ 2 I h N h S h 1 + ν h S h 1 ‖ , ≤ ‖ Ψ h − ( ϖ 1 I r + ϖ 2 I h ) S h − ν h S h − Ψ h + ( ϖ 1 I r + ϖ 2 I h ) S h 1 + ν h S h 1 ‖ , ≤ ‖ ( ϖ 1 I r + ϖ 2 I h ) ( S h 1 − S h ) + ν h ( S h 1 − S h ) ‖ , ≤ ‖ ( ϖ 1 I r + ϖ 2 I h + ν h ) ( S h 1 − S h ) ‖ , ≤ ‖ ( ϖ 1 sup t ∈ J ‖ I r ( t ) ‖ + ϖ 2 sup t ∈ J ‖ I h ( t ) ‖ + ν h ) ( S h 1 − S h ) ‖ , ≤ ( ϖ 1 b 8 + ϖ 2 b 3 + ν h ) ‖ ( S h 1 − S h ) ‖ , ≤ M 1 ‖ ( S h 1 − S h ) ‖ ,

where M 1 = ( ϖ 1 b 7 + ϖ 2 b 3 + ν h ) . L 1 ( t , S h ) holds the Lipschitz property with condition M 1 . Furthermore, we obtain contraction if 0 < M 1 < 1 . Using the aforementioned approach, we can determine

‖ L 2 ( t , E h ) − L 2 ( t , E h 1 ) ‖ ≤ M 2 ‖ E h 1 − E h ‖ , ‖ L 3 ( t , I h ) − L 3 ( t , I h 1 ) ‖ ≤ M 3 ‖ I h 1 − I h ‖ , ‖ L 4 ( t , Q h ) − L 4 ( t , Q h 1 ) ‖ ≤ M 4 ‖ Q h 1 − Q h ‖ , ‖ L 5 ( t , R h ) − L 5 ( t , R h 1 ) ‖ ≤ M 5 ‖ R h 1 − R h ‖ , ‖ L 6 ( t , S r ) − L 6 ( t , S r 1 ) ‖ ≤ M 6 ‖ S r 1 − S r ‖ , ‖ L 7 ( t , E r ) − L 7 ( t , E r 1 ) ‖ ≤ M 7 ‖ E r 1 − R r ‖ , ‖ L 8 ( t , I r ) − L 8 ( t , I r 1 ) ‖ ≤ M 8 ‖ I r 1 − I r ‖ ,

where M 2 = P 1 , M 3 = P 2 , M 4 = P 3 , M 5 = ν h , M 6 = ( ϖ r b 8 + ν r ) , m 7 = P 4 , and M 8 = ν r .

When t = t n , n = 1 , 2 , … , we can have the recursive form for system (15):

(20) S h n ( t ) = 1 − σ G ( σ ) L 1 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 1 ( k , W n − 1 ) ( t − x ) σ − 1 d x , E h n ( t ) = 1 − σ G ( σ ) L 2 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 2 ( k , W n − 1 ) ( t − x ) σ − 1 d x , I h n ( t ) = 1 − σ G ( σ ) L 3 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 3 ( k , W n − 1 ) ( t − x ) σ − 1 d x , Q h n ( t ) = 1 − σ G ( σ ) L 4 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 4 ( k , W n − 1 ) ( t − x ) σ − 1 d x , R h n ( t ) = 1 − σ G ( σ ) L 5 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 5 ( k , W n − 1 ) ( t − x ) σ − 1 d x , S r n ( t ) = 1 − σ G ( σ ) L 6 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 6 ( k , W n − 1 ) ( t − x ) σ − 1 d x , E r n ( t ) = 1 − σ G ( σ ) L 7 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 7 ( k , W n − 1 ) ( t − x ) σ − 1 d x , I r n ( t ) = 1 − σ G ( σ ) L 8 ( t , W n − 1 ) + σ G ( σ ) Γ ( σ ) ∫ 0 t L 8 ( k , W n − 1 ) ( t − x ) σ − 1 d x ,

where the initial conditions are shown in Eq. (4). The difference among the successive terms in Eq. (20) is

(21) B 1 n ( t ) = S h n ( t ) − S h n − 1 ( t ) = 1 − σ G ( σ ) ( L 1 ( t , S h n − 1 ) − L 1 ( t , S h n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 1 ( t , S h n − 1 ) − L 1 ( t , S h n − 2 ) ) ( t − x ) σ − 1 d x , B 2 n ( t ) = E h n ( t ) − E h n − 1 ( t ) = 1 − σ G ( σ ) ( L 2 ( t , E h n − 1 ) − L 2 ( t , E h n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 2 ( t , E h n − 1 ) − L 2 ( t , E h n − 2 ) ) ( t − x ) σ − 1 d x , B 3 n ( t ) = I h n ( t ) − I h n − 1 ( t ) = 1 − σ G ( σ ) ( L 3 ( t , I h n − 1 ) − L 3 ( t , I h n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 3 ( t , I h n − 1 ) − L 3 ( t , I h n − 2 ) ) ( t − x ) σ − 1 d x , B 4 n ( t ) = Q h n ( t ) − Q h n − 1 ( t ) = 1 − σ G ( σ ) ( L 4 ( t , Q h n − 1 ) − L 4 ( t , Q h n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 4 ( t , Q h n − 1 ) − L 4 ( t , Q h n − 2 ) ) ( t − x ) σ − 1 d x , B 5 n ( t ) = R h n ( t ) − R h n − 1 ( t ) = 1 − σ G ( σ ) ( L 5 ( t , R h n − 1 ) − L 5 ( t , R h n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 5 ( t , R h n − 1 ) − L 5 ( t , R h n − 2 ) ) ( t − x ) σ − 1 d x , B 6 n ( t ) = S r n ( t ) − S r n − 1 ( t ) = 1 − σ G ( σ ) ( L 6 ( t , S r n − 1 ) − L 6 ( t , S r n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 6 ( t , S r n − 1 ) − L 6 ( t , S r n − 2 ) ) ( t − x ) σ − 1 d x , B 7 n ( t ) = E r n ( t ) − E r n − 1 ( t ) = 1 − σ G ( σ ) ( L 7 ( t , E r n − 1 ) − L 7 ( t , E r n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 7 ( t , E r n − 1 ) − L 7 ( t , E r n − 2 ) ) ( t − x ) σ − 1 d x , B 8 n ( t ) = I r n ( t ) − I r n − 1 ( t ) = 1 − σ G ( σ ) ( L 8 ( t , I r n − 1 ) − L 8 ( t , I r n − 2 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 8 ( t , I r n − 1 ) − L 8 ( t , I r n − 2 ) ) ( t − x ) σ − 1 d x ,

Taking norm of Eq. (21), and consider their first equation

(22) ∥ B 1 n ( t ) ∥ = ∥ S h n ( t ) − S h n − 1 ( t ) ∥ ≤ 1 − σ G ( σ ) ∥ ( L 1 ( t , S h n − 1 ) − L 1 ( t , S h n − 2 ) ) ∥ + σ G ( σ ) Γ ( σ ) ∫ 0 τ ∥ ( L 1 ( t , S h n − 1 ) − L 1 ( t , S h n − 2 ) ) ∥ ( t − x ) σ − 1 d x ≤ 1 − σ G ( σ ) M 1 ∥ S h n − 1 − S h n − 2 ∥ + σ G ( σ ) Γ ( σ ) M 1 ∫ 0 τ ∥ S h n − 1 − S h n − 2 ∥ ( t − x ) σ − 1 d x ≤ M 1 ∥ B 1 ( n − 1 ) ( t ) ∥ 1 − σ G ( σ ) + t σ G ( σ ) Γ ( σ ) .

We have the following:

(23) ∥ B 1 n ( t ) ∥ ≤ M 1 G ¯ ∥ B 1 ( n − 1 ) ( t ) ∥ ,

where G ¯ = 1 − σ G ( σ ) + t σ G ( σ ) Γ ( σ ) . With the same method, we obtain the following for the rest of equations:

(24)□ ∥ B 2 n ( t ) ∥ ≤ M 2 G ¯ ∥ B 2 ( n − 1 ) ( t ) ∥ , ∥ B 3 n ( t ) ∥ ≤ M 3 G ¯ ∥ B 3 ( n − 1 ) ( t ) ∥ , ∥ B 4 n ( t ) ∥ ≤ M 4 G ¯ ∥ B 4 ( n − 1 ) ( t ) ∥ , ∥ B 5 n ( t ) ∥ ≤ M 5 G ¯ ∥ B 5 ( n − 1 ) ( t ) ∥ , ∥ B 6 n ( t ) ∥ ≤ M 6 G ¯ ∥ B 6 ( n − 1 ) ( t ) ∥ , ∥ B 7 n ( t ) ∥ ≤ M 7 G ¯ ∥ B 7 ( n − 1 ) ( t ) ∥ , ∥ B 8 n ( t ) ∥ ≤ M 8 G ¯ ∥ B 8 ( n − 1 ) ( t ) ∥ .

Theorem 2

The fractional monkeypox infection model (7)can have a solution if we can find K 0 that can satisfy the inequality

1 − σ G ( σ ) + K 0 σ G ( σ ) Γ ( σ ) M i , i = 1 , 2 , … 8 .

Proof

Using Eqs. (23) and (24), we have

(25) ∥ B 1 n ( t ) ∥ ≤ M 1 G ¯ ∥ B 1 ( n − 1 ) ( t ) ∥ , ∥ B 2 n ( t ) ∥ ≤ M 2 G ¯ ∥ B 2 ( n − 1 ) ( t ) ∥ , ∥ B 3 n ( t ) ∥ ≤ M 3 G ¯ ∥ B 3 ( n − 1 ) ( t ) ∥ , ∥ B 4 n ( t ) ∥ ≤ M 4 G ¯ ∥ B 4 ( n − 1 ) ( t ) ∥ , ∥ B 5 n ( t ) ∥ ≤ M 5 G ¯ ∥ B 5 ( n − 1 ) ( t ) ∥ , ∥ B 6 n ( t ) ∥ ≤ M 6 G ¯ ∥ B 6 ( n − 1 ) ( t ) ∥ , ∥ B 7 n ( t ) ∥ ≤ M 7 G ¯ ∥ B 7 ( n − 1 ) ( t ) ∥ , ∥ B 8 n ( t ) ∥ ≤ M 8 G ¯ ∥ B 8 ( n − 1 ) ( t ) ∥ .

The existence of the model solution is confirmed by Theorem 1, and now we shall show that the functions S h ( t ) , … , I r ( t ) show the solutions to Eq. (7). Assume that the following hold:

(26) S h ( t ) − S h ( 0 ) = S h n ( t ) − b 1 n ( t ) , E h ( t ) − E h ( 0 ) = E h n ( t ) − b 2 n ( t ) , I h ( t ) − I h ( 0 ) = I h n ( t ) − b 3 n ( t ) , Q h ( t ) − Q h ( 0 ) = Q h n ( t ) − b 4 n ( t ) , R h ( t ) − R h ( 0 ) = R h n ( t ) − b 5 n ( t ) , S r ( t ) − S r ( 0 ) = S r n ( t ) − b 6 n ( t ) , E r ( t ) − E r ( 0 ) = E r n ( t ) − b 7 n ( t ) , I r ( t ) − I r ( 0 ) = I r n ( t ) − b 8 n ( t ) .

It follows from Eq. (26) that

(27) ∥ b 1 n ( t ) ∥ ≤ 1 − σ G ( σ ) ∥ ( L 1 ( σ , S h n ) − L 1 ( σ , S h n − 1 ) ) ∥ + σ G ( σ ) Γ ( σ ) ∫ 0 σ ∥ ( L 1 ( σ , S h n ) − L 1 ( τ , S h n − 1 ) ) ∥ ( σ − x ) σ − 1 d x , ≤ 1 − σ G ( σ ) M 1 ∥ S h n − S h n − 1 ∥ + σ n G ( σ ) Γ ( σ ) M 1 ∥ S h n − S h n − 1 ∥ .

We have the following after repeating the process:

(28) ∥ b 1 n ( t ) ∥ ≤ [ G ¯ ] n + 1 M 1 n ∥ S h n − S h n − 1 ∥ n ,

when t = K 0 σ , we obtain

(29) ∥ b 1 n ( t ) ∥ ≤ 1 − σ G ( σ ) + M 0 σ G ( σ ) Γ ( σ ) n + 1 M 1 n ∥ S h n − S h n − 1 ∥ n , ∥ b 1 n ( t ) ∥ → 0 .

Using limit on both sides of Eq. (29), we obtain

(30) G ¯ M 1 < 1 .

In similar way, it can be shown that ‖ b 2 n ‖ → 0 , ‖ b 3 n ‖ → 0 , … , ‖ b 8 n ‖ → 0 ,

(31) G ¯ M i < 1 , i = 1 , 2 , … , 8 .

Theorems 1 and 2 ensure the existence of the monkeypox infection fractional model (7) using the fixed point theorem. Now, in the following theorem, we shall show the uniqueness.□

Theorem 3

(Solution uniqueness) The fractional system (7)possesses a unique solution provided that

(32) G ¯ M i < 1 , i = 1 , 2 , … , 8 .

Proof

Consider that S h 1 , E h 1 , I h 1 , Q h 1 , R h 1 , S r 1 , E r 1 , and I r 1 , are another set of solutions of system (7), then

(33) S h ( t ) − S h 1 ( t ) = 1 − σ G ( σ ) ( L 1 ( t , S h ) − L 1 ( t , S h 1 ) ) + σ G ( σ ) Γ ( σ ) ∫ 0 t ( L 1 ( t , S h ) − L 1 ( t , S h 1 ) ) ( t − x ) σ − 1 d x .

Applying the norm on both sides of Eq. (33), we have

(34) ‖ S h ( t ) − S h 1 ( t ) ‖ ≤ 1 − σ G ( σ ) M 1 ‖ S h ( t ) − S h 1 ( t ) ‖ + t σ G ( σ ) Γ ( σ ) M 1 ∣ S h ( t ) − S h 1 ( t ) ∣ .

Since ( 1 − M 1 G ¯ ) > 0 , we obtain ‖ S h ( t ) − S h 1 ( t ) ‖ = 0 . So, S h ( t ) = S h 1 ( t ) . Repeating the aforementioned process, we can have E h ( t ) = E h 1 ( t ) ,…, I r ( t ) = I r 1 ( t ) .□

4 Equilibria and their analysis

We shall investigate the equilibrium points such as the disease-free (DFE) and the endemic equilibrium (EE) associated with the model (7) in the present portion. The DFE of the system (7) shall be denoted by Δ 0 , and is obtained as follows:

D t σ 0 ABC S h = 0 , D t σ 0 ABC E h = 0 , D t σ 0 ABC I h = 0 = 0 , D t σ 0 ABC R h = 0 , D t σ 0 ABC R h = 0 , D t σ 0 ABC S r = 0 , D t σ 0 ABC E r = 0 , D t σ 0 ABC I r = 0 .

So, we obtain

Δ 0 = ( S h 0 , 0 , 0 , 0 , 0 , S r 0 , 0 , 0 ) = Ψ h ν h , 0 , 0 , 0 , 0 , Ψ r ν r , 0 , 0 .

The DFE Δ 0 is useful in the computation of the threshold quantity, say ℛ 0 . The well-known method described in the study by Van den Driessche and Watmough [27] will be used to obtain ℛ 0 for the system (7). We have the matrices

F = 0 ϖ 2 0 S h 0 ϖ 1 S r 0 0 0 0 0 0 0 0 ϖ r 0 0 0 0 , and V = P 1 0 0 0 − η 1 P 2 0 0 0 0 P 4 0 0 0 − κ r ν r ,

where P 1 = ( ν h + η 1 + η 2 ) , P 2 = ( ν h + ε 1 + ψ ) , P 3 = ( ν h + ε 2 + ϕ ) , and P 4 = ( ν r + κ r ) . While using ρ ( F V − 1 ) , we can obtain the basic reproduction number ℛ 0 for the fractional system (7), which is provided by

(35) ℛ 0 = max ϖ 2 η 1 P 1 P 2 , ϖ r κ r P 4 ν r , ℛ 0 = max { ℛ 1 , ℛ 2 } .

We shall present the following theorem to show the locally asymptotically stable (LAS) of model (7).

Theorem 4

The monkeypox fractional system (7) at Δ 0 is LAS if ℛ 1 < 1 and σ ∈ [ 0 , 1 ] , and all the associated eigenvalues λ k for k = 1 , … , 8 hold

∣ arg ( λ ( k ) ) ∣ > σ π 2 .

Proof

The following Jacobian matrix is obtained at the monkeypox-free equilibrium Δ 0 :

J ( Δ 0 ) = − ν h 0 − ϖ 2 0 0 0 0 − ν r ϖ 1 Ψ h ν h Ψ r 0 − P 1 ϖ 2 0 0 0 0 ν r ϖ 1 Ψ h ν h Ψ r 0 η 1 − P 2 0 0 0 0 0 0 η 2 0 − P 3 0 0 0 0 0 0 ψ ϕ − ν h 0 0 0 0 0 0 0 0 − ν r 0 − ϖ r 0 0 0 0 0 0 − P 4 ϖ r 0 0 0 0 0 0 κ r − ν r .

In J ( Δ 0 ) , we have the eigenvalues that clearly have a negative real part: − ( ν h + ε 2 + ϕ ) , − ν h , − ν h , − ν r . The following characteristics equation will be used to determine the final four eigenvalues:

(36) λ 4 + f 1 λ 3 + f 2 λ 2 + f 3 λ + f 4 = 0 ,

where

f 1 = P 1 + P 2 + P 4 + ν r , f 2 = P 1 P 2 ( 1 − ℛ 1 ) + P 4 η r ( 1 − ℛ 2 ) + ( P 1 + P 2 ) ( ν r + P 4 ) , f 3 = ( P 1 + P 2 ) P 4 ν r ( 1 − ℛ 2 ) + P 1 P 2 ( 1 − ℛ 1 ) ( ν r + P 4 ) , f 4 = P 4 P 1 P 2 ν r ( 1 − ℛ 1 ) ( 1 − ℛ 2 ) .

The coefficients f j , where j = 2 , 3 , 4 , in Eq. (36) shall be shown to be positive. f 1 > 0 while f k , where j = 2 , 3 , 4 can be positive if ℛ 0 < 1 . Furthermore, to obtain eigenvalues with negative real parts, the coefficients f j > 0 must satisfy the Routh–Hurtwiz criteria. For Eq. (36), it is required to prove F = f 1 f 2 f 3 > f 3 2 + f 1 2 f 4 . This condition is satisfied as follows:

F = ( P 1 + P 2 ) ( P 4 + ν r ) [ P 1 ( F 3 P 2 + P 2 2 ( 1 − ℛ 1 ) ( P 4 + ν r ) + P 4 ( 1 − ℛ 2 ) ν r ( P 4 + ν r ) ) ] + ( P 1 + P 2 ) ( P 4 + ν r ) [ F 1 P 4 ( 1 − ℛ 2 ) ν r + F 2 P 1 2 ] > 0 ,

where

F 1 = P 2 ( P 4 + ν r ) + P 4 ( 1 − ℛ 2 ) ν r + P 2 2 , F 2 = P 2 ( 1 − ℛ 1 ) ( P 4 + η r ) + P 4 ( 1 − ℛ 2 ) ν r + P 2 2 ( 1 − ℛ 1 ) 2 , F 3 = ( 1 − ℛ 1 ) ( P 4 2 + ν r 2 ) + 2 P 4 ν r ( 1 − ℛ 1 ℛ 2 ) .

The condition F > 0 ensures that the polynomial of order four will give four eigenvalues with negative real parts. Hence, the monkeypox model (3) under equilibrium Δ 0 is LAS if ℛ 0 < 1 .□

4.1 Endemic equilibria and backward bifurcation

The EE of the monkeypox model (3) is denoted by Δ 1 , Δ 1 = ( S h * , E h * , I h * , Q h * , R h * , S r * , E r * , I r * ) , and is calculated as follows:

(37) S h * = Ψ h ν h + Γ * , E h * = Γ * S h * P 1 , I h * = η 1 E h * P 2 , Q h * = η 2 E h * P 3 , R h * = ψ I h * + ϕ Q h * ν h , S r * = Ψ r χ * + ν r , E r * = χ * S r * P 4 , I r * = E r * κ r ν r ,

Using Eq. (37) into Γ * ,

Γ * = ϖ 1 I r * N r * + ϖ 2 I h * N h * , χ * = ϖ r I r * N r * ,

we obtain

χ * = χ * ϖ r κ r ( P 4 + χ * ) ν r + χ * κ r .

Furthermore, we put χ * into expression Γ * and finally obtain the following result:

a 1 Γ * 2 + a 2 Γ * + a 3 = 0 ,

where

a 1 = ϖ r P 4 ( P 3 η 1 ( ν h + ψ ) + P 2 ( η 2 ( ν h + ϕ ) + P 3 ν h ) ) , a 2 = ϖ 1 P 4 ν r [ P 3 η 1 ( ν h + η 1 ψ ) + P 2 η 2 ( ν h + ϕ ) ] ( 1 − ℛ 2 ) + P 3 ν h [ P 2 ( ϖ r ( P 1 P 4 − ϖ 1 κ r ) + ϖ 1 P 4 ν r ) − ϖ 2 η 1 ϖ r P 4 ] , a 3 = ϖ 1 P 1 P 2 P 3 P 4 ν h ν r ( 1 − ℛ 2 ) .

The existence of the unique endemic equilibria and the possible existence of the backward bifurcation in system (7) are summarized in the following theorem.

Theorem 5

The fractional system (7) has:

When a 3 < 0 iff ℛ 2 > 1 , we obtain unique EE;
When a 2 < 0 and a 3 = 0 , or a 2 2 − 4 a 1 a 3 = 0 , we obtain unique EE;
When a 3 > 0 , a 2 < 0 , and a 2 2 − 4 a 1 a 3 > 0 , we obtain two endemic equilibria.

The case (i) of Theorem (5) gives a clear indication of the unique EE for the fractional system (7) if ℛ 2 > 1 . For backward bifurcation, case (iii) has enough conditions. In the procedure to obtain the result for the backward bifurcation, in Eq. (7), we choose a 2 2 − 4 a 1 a 3 = 0 and then solve for the critical values of ℛ 2 , denoted by ℛ 2 c , given by

ℛ 2 c = 1 − a 2 2 4 a 1 ϖ 1 P 1 P 2 P 3 P 4 ν h ν r .

The backward bifurcation would occur for the values of ℛ 2 , such that ℛ 2 c < ℛ 2 < 1 . The bifurcation diagram is shown in Figure 1 by considering the values of the parameters shown in the numerical section, except ϖ r = 0.0007028 , η 1 = 0.606 , and κ r = 0.000399 . With the list of parameter values, we obtain ℛ 2 = 0.5404 < 1 .

Figure 1

Backward bifurcation in monkeypox model.

4.2 Globally asymptotically stable (GAS) of EE

The GAS of the fractional system (7) will be carried out here. The following is given for (7) at EE:

(38) Ψ h = Γ * ( t ) S h * + ν h S h * , Γ ( t ) * S h * = P 1 E h * , η 1 E h * = P 2 I h * , η 2 E h * = P 3 Q h * , Ψ r = χ * ( t ) S r * + ν r S r * , χ ( t ) * S r * = P 4 E r * , κ r E r * = ν r I r * .

We will use the result in Eq. (38) later in the proof of the following theorem:

Theorem 6

The fractional system (7) is GAS if ℛ 0 > 1 .

Proof

We consider the Lyapunov function as follows:

(39) ℒ ( t ) = S h − S h * − S h * ln S h S h * + E h − E h * − E h * ln E h E h * + ϖ 2 I h * S h * η 1 E h * I h − I h * − I * ln I h I h * + ϖ 1 I r * S h * η 2 E h * Q h − Q h * − Q h * ln Q h Q h * + S r − S r * − S r * ln S r S r * + E r − E r * − E r * ln E r E r * + ϖ r I r * S r * κ r E r * I r − I r * − I r * ln I r I r * .

We obtain the following while taking the time derivative of Eq. (39):

(40) D t σ 0 ABC ℒ = 1 − S h * S h D t σ 0 ABC S h + 1 − E h * E h D t σ 0 ABC E h + ϖ 2 I h * S h * η 1 E h * 1 − I h * I h D t σ 0 ABC I h + ϖ 1 I r * S h * η 2 E h * 1 − Q h * Q h D t σ 0 ABC Q h + 1 − S r * S r D t σ 0 ABC S r + 1 − E r * E r D t σ 0 ABC E r + ϖ r I r * S r * κ r E r * 1 − I r * I r D t σ 0 ABC I r .

Direct calculation of the terms in Eq. (40) are obtained as follows:

(41) 1 − S h * S h D t σ 0 ABC S h = 1 − S h * S h Ψ h − ϖ 1 I r N r + ϖ 2 I h N h S h − ν h S h , ≤ 1 − S h * S h [ Ψ h − ( ϖ 1 I r + ϖ 2 I h ) S h − ν h S h ] , ≤ 1 − S h * S h [ ( ϖ 1 I r * + ϖ 2 I h * ) S h * + ν h S h * − ( ϖ 1 I r + ϖ 2 I h ) S h − ν h S h ] , ≤ ν h S h * 2 − S h * S h − S h S h * + ϖ 1 I r * S h * 1 − S h * S h − I r S h I r * S h * + I r I r * + ϖ 2 I h * S h * 1 − S h * S h − I h S h I h * S h * + I h I h * , ≤ ϖ 1 I r * S h * 1 − S h * S h − I r S h I r * S h * + I r I r * + ϖ 2 I h * S h * 1 − S h * S h − I h S h I h * S h * + I h I h * ,

(42) 1 − E h * E h D t σ 0 ABC E h = 1 − E h * E h [ ( ϖ 1 I r + ϖ 2 I h ) S h − P 1 E h ] , = 1 − E h * E h ( ϖ 1 I r + ϖ 2 I h ) S h − ( ϖ 1 I r * + ϖ 2 I h * ) S h * E h E h * , = ϖ 1 I r * S h * 1 − E h E h * + I r S h I r * S h * − I r S h E h * I r * S h * E h + ϖ 2 I h * S h * 1 − E h E h * + I h S h I h * S h * − I h S h E h * I h * S h * E h ,

(43) ϖ 2 I h * S h * η 1 E h * 1 − I h * I h D t σ 0 ABC I h = ϖ 2 I h * S h * η 1 E h * 1 − I h * I h [ η 1 E h − P 2 I h ] , = ϖ 2 I h * S h * E h * 1 − I h * I h E h − E h * I h * I h , = ϖ 2 I h * E h * 1 − I h I h * − E h I h * E h * I h + E h E h * ,

(44) ϖ 1 I r * S h * η 2 E h * 1 − Q h * Q h D t σ 0 ABC Q h = ϖ 1 I r * S h * η 2 E h * 1 − Q h * Q h [ η 2 E h − P 3 Q h ] , = ϖ 1 I r * S h * E h * 1 − Q h * Q h E h − E h * Q h * Q h , = ϖ 1 I r * S h * 1 + E h E h * − E h Q h * E h * Q h − Q h Q h * .

(45) 1 − S r * S r D t σ 0 ABC S r ≤ 1 − S r * S r [ Ψ r − ϖ r I r S r − ν r S r ] , = 1 − S r * S r [ ϖ r I r * S r * + ν r S r * − ϖ r I r S r − ν r S r ] , = ν r S r * 2 − S r S r * − S r * S r + ϖ r I r * S r * 1 − S r * S r − I r S r I r * S r * + I r I r * ,

(46) 1 − E r * E r D t σ 0 ABC E r ≤ 1 − E r * E r [ ϖ r I r S r − P 4 E r ] , ≤ 1 − E r * E r ϖ r I r S r − ϖ r I r * S r * E r E r * , ≤ ϖ r I r * S r * 1 − E r E r * + I r S r I r * S r * − I r S r E r * E r I r * S r * ,

(47) ϖ r I r * S r * κ r E r * 1 − I r * I r D t σ 0 ABC I r ≤ ϖ r I r * S r * κ r E r * 1 − I r * I r [ κ r E r − ν r I r ] , ≤ ϖ r I r * S r * E r * 1 − I r * I r E r − E r * I r * I r , ≤ ϖ r I r * S r * 1 + E r E r * + I r E r I r * E r * − I r I r * .

Using Eqs. (41)–(47) in Eq. (40), and after some simplification, we achieve the following:

(48) D t σ 0 ABC ℒ = ϖ 1 I r * S h * 3 − S h * S h − Q h Q h * + I r I r * − E h Q h * E h * Q h − I r S h E h * E h I r * S h * + ϖ 2 I h * S h * 3 − S h * S h − E h I h * I h E h * − I h S h E h * E h I h * S h * + ϖ r I r * S r * 3 − S r * S r − E r I r * I r E r * − I r S r E r * E r I r * S r * .

The following are nonnegative due to the property of arithmetic geometric mean:

3 − S h * S h − E h I h * I h E h * − I h S h E h * E h I h * S h * ≤ 0 , 3 − S r * S r − E r I r * I r E r * − I r S r E r * E r I r * S r * ≤ 0 ,

and D t σ 0 ABC ℒ ≤ 0 if

(49) 3 − S h * S h − Q h Q h * + I r I r * − E h Q h * E h * Q h − I r S h E h * E h I r * S h * ≤ 0 .

So, the subset for D t σ 0 ABC ℒ ≤ 0 is Δ 1 which is invariant and largest. Thus, the EE Δ 1 is GAS if ℛ 0 > 1 .□

5 Numerical scheme

In this section, we present the numerical algorithm for the numerical solution of the monkeypox fractional model (7). This method is based on the Atangana–Baleanu fractional derivative for the numerical solution of the fractional-order models. We follow the procedure given in the study by Toufik and Atangana [28] and derive the algorithm first in the general case and later for our model equations. Let us consider the following nonlinear fractional differential equations:

(50) D t σ 0 ABC χ ( t ) = f ( t , χ ( t ) ) , χ ( 0 ) = χ 0 .

We arrive at the following equation, after utilizing the fundamental theorem of FC:

(51) χ ( t ) = χ ( 0 ) + ( 1 − σ ) G ( σ ) f ( t , χ ( t ) ) + σ Γ ( σ ) G ( σ ) ∫ 0 t f ( x , χ ( x ) ) ( t − x ) σ − 1 d x .

At t = t j + 1 , and further approximating the function f ( x , χ ( x ) ) and simplifying, we finally obtain

(52) χ j + 1 = χ ( 0 ) + ( 1 − σ ) G ( σ ) f ( t j , χ ( t j ) ) + σ G ( σ ) ∑ m = 0 j h σ f ( t m , χ ( t m ) ) Γ ( σ + 2 ) B j , m 1 − h σ f ( t m − 1 , χ ( t m − 1 ) ) Γ ( α + 2 ) B j , m 2 ,

where

(53) B j , m 1 = [ ( j + 1 − m ) σ ( j − m + 2 + σ ) − ( j − m ) σ ( j − m + 2 + 2 σ ) ] , B j , m 2 = [ ( j − m + 1 ) σ + 1 − ( j − m ) σ ( j − m + 1 + σ ) ] .

We shall adapt the scheme shown in Eqs. (52)–(54), apply it to the fractional model (7), and obtain:

S h j + 1 = S h ( 0 ) + ( 1 − σ ) G ( σ ) f 1 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 1 ( D 1 ) B j , m 1 − f 1 ( D 2 ) B j , m 2 } , E h j + 1 = E h ( 0 ) + ( 1 − σ ) G ( σ ) f 2 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 2 ( D 1 ) B j , m 1 − f 2 ( D 2 ) B j , m 2 } , I h j + 1 = I h ( 0 ) + ( 1 − σ ) G ( σ ) f 3 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 3 ( D 1 ) B j , m 1 − f 3 ( D 2 ) B j , m 2 } , Q h j + 1 = Q h ( 0 ) + ( 1 − σ ) G ( σ ) f 4 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 4 ( D 1 ) B j , m 1 − f 4 ( D 2 ) B j , m 2 } ,

(54) R h j + 1 = R h ( 0 ) + ( 1 − σ ) G ( σ ) f 5 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 5 ( D 1 ) B j , m 1 − f 5 ( D 2 ) B j , m 2 } , S r j + 1 = S r ( 0 ) + ( 1 − σ ) G ( σ ) f 6 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 6 ( D 1 ) B j , m 1 − f 6 ( D 2 ) B j , m 2 } , E r j + 1 = E r ( 0 ) + ( 1 − σ ) G ( σ ) f 7 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 7 ( D 1 ) B j , m 1 − f 7 ( D 2 ) B j , m 2 } , I r j + 1 = I r ( 0 ) + ( 1 − σ ) G ( σ ) f 8 ( t j , S h j , E h j , I h j , Q h j , R h j , S r j , E r j , I r j ) + σ G ( σ ) ∑ m = 0 j h σ Γ ( σ + 2 ) { f 8 ( D 1 ) B j , m 1 − f 8 ( D 2 ) B j , m 2 } ,

where

D 1 = t m , S h ( t m ) , E h ( t m ) , I h ( t m ) , Q h ( t m ) , R h ( t m ) , S r ( t m ) , E r ( t m ) , I r ( t m ) , D 2 = t m − 1 , S h ( t m − 1 ) , E h ( t m − 1 ) , I h ( t m − 1 ) , Q h ( t m − 1 ) , R h ( t m − 1 ) , S r ( t m − 1 ) , I r ( t m − 1 ) ,

and

f 1 ( D 1 ) = Ψ h − Γ ( t ) S h − ν h S h , f 2 ( D 1 ) = Γ ( t ) S h − ( η 1 + η 2 + ν h ) E h , f 3 ( D 1 ) = η 1 E h − ( ψ + ν h + ε 1 ) I h , f 4 ( D 1 ) = η 2 E h − ( ϕ + ε 2 + ν h ) Q h , f 5 ( D 1 ) = ψ I h + ϕ Q h − ν h R h , f 6 ( D 1 ) = Ψ r − χ ( t ) S r − ν r S r , f 7 ( D 1 ) = χ ( t ) S r − ( ν r + κ r ) E r , f 8 ( D 1 ) = κ r E r − ν r I r .

5.1 Numerical simulation

Here, we numerically analyze the fractional-order system (7) using the parameter values given in the study by Allehiany et al. [24] as follows: Ψ h = ν h × N h ( 0 ) , ν h = 1 76.4 × 365 , ϖ 1 = 0.000041552 , ϖ 2 = 0.6307 , η 1 = 0.0306 , η 2 = 0.0571 , ε 1 = 0.3356 , ε 2 = 0.04149 , ψ = 0.0369 , ϕ = 0.02951 , Ψ r = ν r × N r ( 0 ) , ν r = 1 ∕ ( 5 × 365 ) , ϖ r = 0.1028 , and κ r = 0.0799 . The time unit is considered in days. We simulate the model using the fractional scheme shown in above section and present the graphical results. The stability and convergence of the fractional-order σ for its many values on the dynamics of human and rodent populations are shown in Figures 2 and 3. The solution behaviors of the human and rodent compartments show that the results are converge and stable for the proposed values of σ .

Figure 2

The graph shows the human population. (a–d) denote susceptible, exposed, infected, and quarantined population, respectively.

$Figure 3 Numerical result of the recovered human and the rodents compartments when σ = 1 , 0.99 , 0.98 , 0.97 \sigma =1,0.99,0.98,0.97 . Subgraph (a) shows the recovered human population, while (b–d), respectively, show the susceptible, exposed, and infected rodent population.$

Figure 3

Numerical result of the recovered human and the rodents compartments when σ = 1 , 0.99 , 0.98 , 0.97 . Subgraph (a) shows the recovered human population, while (b–d), respectively, show the susceptible, exposed, and infected rodent population.

Figure 4 shows the behavior of the human compartments when ϖ 1 is varied and σ = 0.96 is fixed. From graphical result, we can see that there is decrease in the cases when varying ϖ 1 . By follows the guidelines such as, eliminating the shelter, food sources, and water for the rodents, a better decrease in the future cases will be observed.

$Figure 4 Numerical results of the human compartments when σ = 0.96 \sigma =0.96 and ϖ 1 = 0.000041552 , 0.000031552 , 0.000021552 , 0.000011552 {\varpi }_{1}=0.000041552,0.000031552,0.000021552,0.000011552 , whereas (a–c) show exposed, infected, and quarantined people, respectively.$

Figure 4

Numerical results of the human compartments when σ = 0.96 and ϖ 1 = 0.000041552 , 0.000031552 , 0.000021552 , 0.000011552 , whereas (a–c) show exposed, infected, and quarantined people, respectively.

Figure 5 describes the dynamics of human compartments with the variation in the contact parameter ϖ 2 and σ = 0.96 fixed. Decreasing the contact between human to human, the number of monkeypox infected cases are decreased. Rapid case identification and surveillance are crucial for epidemic containment. Intimate contact with ill patients is the big risk that generate the infected cases in the disease outbreak. Healthcare workers and family members are more at risk for infection. While treating individuals with a monkeypox virus infection that has been suspected or confirmed, or when handling specimens from such patients, health workers should adhere to the prescribed infection control methods.

$Figure 5 Numerical result of the human compartments when σ = 0.96 \sigma =0.96 and ϖ 1 = 0.000041552 , 0.000031552 , 0.000021552 , 0.000011552 {\varpi }_{1}=0.000041552,0.000031552,0.000021552,0.000011552 , whereas (a–c) show, the exposed, infected, and quarantined individuals, respectively.$

Figure 5

Numerical result of the human compartments when σ = 0.96 and ϖ 1 = 0.000041552 , 0.000031552 , 0.000021552 , 0.000011552 , whereas (a–c) show, the exposed, infected, and quarantined individuals, respectively.

The parameter η 1 and its impact on the human compartment is shown in Figure 6. Decreasing the parameter η 1 , the number of infected humans decreased. An infected person with the monkeypox virus shall be isolated and also their close contact with other healthy people shall be minimized in order to control the infection spread further in the human population.

$Figure 6 Variation in η 1 {\eta }_{1} and σ = 0.96 \sigma =0.96 fixed, and its impact on the population. Subgraphs (a–c) show the exposed, infected, and quarantined people, respectively.$

Figure 6

Variation in η 1 and σ = 0.96 fixed, and its impact on the population. Subgraphs (a–c) show the exposed, infected, and quarantined people, respectively.

The parameter η 2 that causes humans to be quarantined when it is identified that they have a risk of being infected is shown in Figure 7. The result in Figure 7 indicates that the number of exposed, infected, and quarantined population decreased when the quarantine rate increased.

$Figure 7 Variation in η 2 {\eta }_{2} and σ = 0.96 \sigma =0.96 fixed, and their impact on population. Subgraphs (a–c) show, respectively, the exposed, infected, and the quarantined people.$

Figure 7

Variation in η 2 and σ = 0.96 fixed, and their impact on population. Subgraphs (a–c) show, respectively, the exposed, infected, and the quarantined people.

6 Conclusion

In this work, a fractional model in the Atangana–Baleanu derivative is proposed and obtained the dynamics of the monkeypox disease with the real data in the USA. We formulated the model for the monkeypox infection in Atangana–Baleanu derivative. The existence and uniqueness of the system are explored briefly. The local asymptotical result for the fractional system was obtained and discussed. We presented the LAS of the fractional system for ℛ 0 < 1 . The endemic equilibria and their existence for fractional system were presented. The backward bifurcation analysis for fractional system was explored. The GAS for endemic case has been shown when ℛ 0 > 1 .

The values of the parameters obtained from the real data in the study by Allehiany et al. [24] are used to perform the numerical simulation for the fractional system. We solved the fractional system and presented the results graphically. The findings suggest that minimizing interaction between rats and people by removing their access to food, water, and shelter, among other things, will reduce infection cases. Furthermore, by routinely disposing waste within or outside the house, the risk of rats can be reduced. It is also possible to reduce human-to-human transmission to reduce the number of cases in the future. Avoid handling any clothing, linens, blankets, or other items that have come into contact with an infected person or animal. Divvy up the healthy people from the monkeypox victims. Wash your hands with soap and water thoroughly after coming into touch with any infected individuals or animals. Stay away from any animals that could have the infection.

Acknowledgments

The authors extend their appreciation to the Deanship of Scientific Research at the King Khalid University for funding this work through large Groups project under grant number RGP.2/227/43. Kamal Shah and Thabet Abdeljawad are thankful to Prince Sultan University for funding the publication of this manuscript and support through the theoretical and applied sciences research lab.

Funding information: This work was funded through large Groups project under grant number RGP.2/227/43 by the Deanship of Scientific Research at the King Khalid University.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors declare that there is no potential conflict of interests regarding the publications of this article.
Data availability statement: Data are available on the reference shown inside the manuscript.

References

[1] World Health Organization. https://www.who.int/news-room/fact-sheets/detail/monkeypox, note= accessed on March 05, 2023. Search in Google Scholar

[2] Jezek Z, Szczeniowski M, Paluku K, Mutombo M, Grab B. Human monkeypox: confusion with chickenpox. Acta Tropica. 1988;45(4):297–307. Search in Google Scholar

[3] Guo Y, Li T. Fractional-order modeling and optimal control of a new online game addiction model based on real data. Commun Nonlinear Sci Numer Simul. 2023:121:107221. 10.1016/j.cnsns.2023.107221Search in Google Scholar

[4] Baba IA, Humphries UW, Rihan FA. Role of vaccines in controlling the spread of COVID-19: a fractional-order model. Vaccines. 2023;11(1):145. 10.3390/vaccines11010145Search in Google Scholar PubMed PubMed Central

[5] Abbes A, Ouannas A, Shawagfeh N, Jahanshahi H. The fractional-order discrete COVID-19 pandemic model: stability and chaos. Nonlinear Dynamics. 2023;111(1):965–83. 10.1007/s11071-022-07766-zSearch in Google Scholar PubMed PubMed Central

[6] Asamoah JK, Addai E, Arthur YD, Okyere E. A fractional mathematical model for listeriosis infection using two kernels. Decision Anal J. 2023;6:100191. 10.1016/j.dajour.2023.100191Search in Google Scholar

[7] George R, Mohammadi K, Mohammadi H, Ghorbanian R, Rezpour S, Duc A. The study of cholera transmission using an SIRZ fractional-order mathematical model. Fractals. 2023. 10.1142/S0218348X23400534Search in Google Scholar

[8] Azeem M, Farman M, Abukhaled M, Nisar KS, Akgul A. Epidemiological analysis of human liver model with fractional operator. Fractals. 2023;31(4):2340047. 10.1142/S0218348X23400479Search in Google Scholar

[9] Okyere S, Ackora-Prah J. Modeling and analysis of monkeypox disease using fractional derivatives. Results Eng. 2023;17:100786. 10.1016/j.rineng.2022.100786Search in Google Scholar PubMed PubMed Central

[10] Farman M, Akgü l A, Tekin MT, Akram MM, Ahmad A, Mahmoud EE, et al. Fractal fractional-order derivative for HIV/AIDS model with Mittag–Leffler kernel. Alexandr Eng J. 2022;61(12):10965–80. 10.1016/j.aej.2022.04.030Search in Google Scholar

[11] Evirgen F, Ucar E, OOzdemir N, Altun E, Abdeljawad T. The impact of nonsingular memory on the mathematical model of Hepatitis C virus. Fractals. 2023;31(4):2340065. 10.1142/S0218348X23400650Search in Google Scholar

[12] Koca I. Analysis of rubella disease model with non-local and non-singular fractional derivatives. Int J Optim Control Theories Appl (IJOCTA). 2018;8(1):17–25. 10.11121/ijocta.01.2018.00532Search in Google Scholar

[13] Bhatter S, Jangid K, Abidemi A, Owolabi K, Purohit S, et al. A new fractional mathematical model to study the impact of vaccination on COVID-19 outbreaks. Decision Analytics J. 2023;6:100156. 10.1016/j.dajour.2022.100156Search in Google Scholar

[14] Karaagac B, Owolabi KM, Pindza E. A computational technique for the Caputo fractal-fractional diabetes mellitus model without genetic factors. Int J Dyn Control. 2023;11:1–18. 10.1007/s40435-023-01131-7Search in Google Scholar PubMed PubMed Central

[15] Khan H, Alzabut J, Baleanu D, Alobaidi G, Rehman MU. Existence of solutions and a numerical scheme for a generalized hybrid class of n-coupled modified ABC-fractional differential equations with an application. AIMS Math. 2023;8(3):6609–25. 10.3934/math.2023334Search in Google Scholar

[16] Hussain S, Tunccc O, Rahman G, Khan H, Nadia E. Mathematical analysis of stochastic epidemic model of MERS-corona & application of ergodic theory. Math Comput Simulat. 2023;207:130–50. 10.1016/j.matcom.2022.12.023Search in Google Scholar PubMed PubMed Central

[17] Khan H, Alzabut J, Shah A, He ZY, Etemad S, Rezapour S, et al. On fractal-fractional waterborne disease model: a study on theoretical and numerical aspects of solutions via simulations. Fractals. 2023;31:2340055. 10.1142/S0218348X23400558Search in Google Scholar

[18] Madubueze C, Onwubuya IO, Nkem GN, Chazuka Z. On the transmission dynamics of the monkeypox virus in the presence of environmental transmission. Front Appl Math Stat. 2022;8:1061546. 10.3389/fams.2022.1061546Search in Google Scholar

[19] Somma SA, Akinwande NI, Chado UD. A mathematical model of monkey pox virus transmission dynamics. Ife J Sci. 2019;21(1):195–204. 10.4314/ijs.v21i1.17Search in Google Scholar

[20] Lasisi N, Akinwande N, Oguntolu F. Development and exploration of a mathematical model for transmission of monkey-pox disease in humans. Math Models Eng. 2020;6(1):23–33. 10.21595/mme.2019.21234Search in Google Scholar

[21] Usman S, Adamu II. Modeling the transmission dynamics of the monkeypox virus infection with treatment and vaccination interventions. J Appl Math Phys. 2017;5(12):2335. 10.4236/jamp.2017.512191Search in Google Scholar

[22] Emeka P, Ounorah M, Eguda F, Babangida B. Mathematical model for monkeypox virus transmission dynamics. Epidemiol Open Access. 2018;8(3):1000348. Search in Google Scholar

[23] Peter OJ, Kumar S, Kumari N, Oguntolu FA, Oshinubi K, Musa R. Transmission dynamics of Monkeypox virus: a mathematical modelling approach. Model Earth Syst Environ. 2022;8(3):3423–34. 10.1007/s40808-021-01313-2Search in Google Scholar PubMed PubMed Central

[24] Allehiany F, DarAssi MH, Ahmad I, Khan MA, Tag-eldin EM. Mathematical Modeling and backward bifurcation in monkeypox disease under real observed data. Results Phys. 2023;50:106557. 10.1016/j.rinp.2023.106557Search in Google Scholar PubMed PubMed Central

[25] Atangana A, Baleanu D. New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Thermal Sci. 2016;20(2):763–9. 10.2298/TSCI160111018ASearch in Google Scholar

[26] Panda SK. Applying fixed point methods and fractional operators in the modelling of novel coronavirus 2019-nCoV/SARS-CoV-2. Results Phys. 2020;19:103433. 10.1016/j.rinp.2020.103433Search in Google Scholar PubMed PubMed Central

[27] Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002;180(1–2):29–48. 10.1016/S0025-5564(02)00108-6Search in Google Scholar PubMed

[28] Toufik M, Atangana A. New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models. Europ Phys J Plus. 2017;132:1–16. 10.1140/epjp/i2017-11717-0Search in Google Scholar

Received: 2023-07-05

Revised: 2023-08-21

Accepted: 2023-08-28

Published Online: 2023-09-20

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2023-0111

Keywords for this article

Monkeypox outbreak USA; fractional-order model; local and global stability; novel numerical scheme

Creative Commons

BY 4.0