In this research endeavor, we undertake a comprehensive analysis of a compartmental model for the monkeypox disease, leveraging the Atangana–Baleanu fractional derivative framework. Our primary objective is to investigate the effectiveness of a range of control strategies in containing the transmission of this infectious ailment. The parameterization of the model is executed meticulously via the application of the maximum likelihood estimation technique. Our study involves a rigorous mathematical analysis of the considered model, which encompasses an exploration of the existence and uniqueness of solutions, as well as the establishment of conditions ensuring the compactness and continuity of these solutions. Subsequently, we embark on an extensive stability analysis of the model, complemented by the computation of both the effective and basic reproduction numbers. These calculations are instrumental in illuminating the long-term behavior of the epidemic. Additionally, we perform a sensitivity analysis of the basic reproduction number to discern the influence of various factors on disease transmission dynamics. To derive our numerical results, we implement the Adams–Bashforth predictor–corrector algorithm tailored for the Atangana–Baleanu fractional derivatives. We employ this numerical technique to facilitate the simulation of the model under a spectrum of fractional-order values, offering a visual representation of our findings. Our study underscores the pivotal roles of infection awareness, vaccination campaigns, and effective treatment in significantly curtailing disease transmission, thus contributing valuable insight to the field of epidemiology.

The monkeypox disease caused by the virus MPVX (the monkeypox virus) was detected for the very first time in the year 1958, but it was confused with smallpox until the beginning of the 1970s when it was clearly identified as monkeypox. Since the year 2002, specifically after May, at least 15 000 reports of monkeypox-infected cases have been recorded in every continent except for Antarctica. The infection is transmitted through direct contact. It can spread from infected rodents to humans and from infected humans to symptomatic ones. Since the outbreak of monkeypox, in several parts of Africa, researchers have started working on its transmission dynamics (TeWinkel, 2019; Bankuru , 2020; Leandry and Mureithi, 2023). According to the observation, the population of synanthropic rodents has rapidly grown in Africa, which has thereby led to the spread of the infection in various regions. Monkeypox, according to the WHO, is a variant of Arthur pox virus belonging to the same family as that of smallpox; however, the noted symptoms of monkeypox in comparison with smallpox are not that contagious and extreme. The most evidently observed symptoms of monkeypox include headaches, muscular pain, lymph nodes becoming swollen, fever, tiredness, rashes, and lesions (Usman and Adamu, 2017). Normally, the rate of mortality due to the infection is 10% of the infected population. A major number of deaths are recorded among children who are young primarily within the age group of 10  years. Five to twenty days is generally considered as the incubation time period. A total number of 86 516 confirmed cases of monkeypox have been recorded all over 113 countries as of March 18, 2023. In a span of 1 year, from May 2022 to 2023, the infection has widely spread across various regions of Asia, America, Africa, Europe, and Oceania. This was the first instance of the infection out-breaking rapidly in countries outside Central and West Africa. As of now, no absolutely clear and effective treatment methods have been found to control the transmission of the infection; however, vaccination does play a significant role in disease prevention as observed in a considerable number of cases. Antiviruses like vaccinia immune globulin, Tecovirimat, and Brincidofovir are proven effective in controlling the disease spread to a significant extent. Earlier, the vaccination given for smallpox had served 85% efficiently in containing the monkeypox infection; however, due to the successful eradication of smallpox worldwide, this control strategy can no longer be implemented. Due to lack of attention, not much study has been done on the transmission dynamics of this disease in the past.

Fractional calculus has garnered notable attention in recent years, emerging as a pivotal domain of mathematics distinguished by its rigorous applications across diverse facets of mathematical epidemiology (Atangana and Owolabi, 2018; Atangana and Qureshi, 2019; Biswas , 2022; Li , 2021; Baleanu , 2018). This prominence is attributed to the prevalence of fractional-order differential equations, which leverage non-local operators (Toufik and Atangana, 2017; Solís-Pérez , 2018). Such operators are quintessential and widely adopted tools in the realm of mathematical modeling. The salient distinction of fractional-order systems is their capacity to encompass a broader spectrum of problems, offering enhanced clarity and precision compared to traditional integer-order systems when investigating an array of real-world challenges within the purview of mathematical epidemiology (Liu , 2020; Majee , 2023). Furthermore, an inherent advantage of fractional operators, in contrast to their ordinary counterparts (Pandey 2013; Prasad and Bali, 2023; Saifuddin , 2016, 2017; Samanta, 2017; Samanta 2013; Stauch 2011), is their intrinsic hereditary properties, rendering them particularly appealing in scientific applications (Abro, 2020; Vijayalakshmi and Roselyn Besi, 2022).

Generally, when we consider fractional calculus, it is just an extended version of ordinary calculus (Agaba 2017; Ahmed 2021; Agusto, 2017; Ayinla 2021; Biswas 2017a; Cai 2017) where the order of the derivatives or integrals can be taken arbitrary as real or complex values (Yu , 2023; Huang and Yu, 2023; Parmar , 2023; Tasman, 2015; Zuo 2022). The motivation behind replacing ordinary calculus with fractional operators is that the latter can prove to be much more efficient in modeling various real-world problems specifically when the dynamics of the system keep changing with respect to the inherent constraints defined for the given system (Jan , 2019; Chinnathambi , 2021; Kumar , 2019) Also, fractional derivatives can be local and non-local in nature. Local fractional derivatives have recently found great relevance in mathematical modeling primarily due to their ability to retain certain significant characteristics of ordinary derivatives; however, these local operators do tend to lose the memory property, which is a usual hereditary characteristic of non-local fractional-order derivatives making them more significant (Panhwer , 2022; Abro , 2023). That is why, we have used the Atangana–Baleanu fractional derivative for our analysis as it is non-local in nature and has Mittag–Leffler function as its kernel; it is nonsingular as well (Aslam , 2021; Atangana and Koca, 2016; Qureshi and Jan, 2021; Srivastava , 2021; Roy and Pop, 2021).

Mathematical modeling involving fractional calculus has been a very significant area of study for scientists and epidemiologists working in this field (Din 2022; Du 2013; Lia 2017; Moore 2019; Muhammad Altaf and Atangana, 2019). Some of the prominent research in this domain can be found in the works of (Alzubaidi , 2023; Somma , 2019; Adel , 2023; Adom-Konadu, 2023). Mathematical models described by a system of fractional differential equations tend to possess a greater degree of freedom than ordinary differential equations (Dye and Wolpert, 1988; Geweke, 1991; Haario 2006; Haario 2001; Khan and Atangana, 2022; Laine, 2008; Lenhart and Workman, 2007; Okosun and Makinde, 2014); that is why they are considered more efficient in studying the dynamics of any particular disease for a given set of data as well as for capturing memory effects (Baleanu , 2020; Ali , 2021). In the present times, a wide range of fractional derivatives have been defined and employed for the mathematical modeling of diseases primarily due to their much more advanced characteristics, which helps in a better analysis of various practical real-world scenarios (Podlubny, 1999; Samko , 1993). Previously, various mathematical and statistical techniques had been used to discuss the impact of non-local operators and the fading memory effect by employing fractional differential equations. Various operators with fractional order have been implemented using local and non-local kernels (Jajarmi and Baleanu, 2018; Odibat, 2006, 2010). The Caputo and Riemann–Liouville fractional operators have one drawback; the kernel function used a singular despite being non-local. This singularity can be a hindrance when solving real-world problems. To overcome this drawback, Atangana and Baleanu defined a fractional order with a new operator, which involves the Mittag–Leffler function as its kernel, which is nonsingular in nature. So, in addition to all the properties of Caputo and Riemann–Liouville fractional derivatives, the kernel being nonsingular helps in a better understanding of the system's dynamics while dealing with practical mathematical models (Rahman , 2021; Rezapour , 2020).

A lot of work has been previously done by researchers to study the dynamics of monkeypox infection using both deterministic and fractional systems of differential equations. Okyere and Ackora-Prah (2023) in their work showed the impacts of using fractional-order derivatives on various compartments of their model under study. The influence of fractional derivatives in the dynamics of the monkeypox disease was highlighted, and the numerical simulation showed the conditions under which the infection could be eradicated in a span of five days. Adel (2023) in their work had investigated the dynamics of a novel fractional-order monkeypox model with optimal control. Their results provided a better insight into various previously used prevention and control strategies and suggested some new strategies for controlling disease transmission. Peter (2022, 2023) emphasized the mathematical modeling of the monkeypox disease using fractional derivatives; they determined the equilibrium points, conducted the stability analysis, and formulated various optimal control strategies for the prevention and control of the infection. They offered numerous parameters for controlling the infection, which will be of tremendous help to the community's policymakers. In their research, Alharbi (2022) described a novel fractional model of monkeypox and did the stability analysis to get a better understanding of the crucial features of the monkeypox virus, thereby helping the decision-makers in containing the disease transmission.

The distinctive contribution of our research is grounded in the novel approach we have employed to formulate a compartmental model for the monkeypox disease. In this innovative model, we subdivide the human population into distinct classes, predicated on their awareness or lack thereof regarding the infection. This stratification based on the awareness factor represents a noteworthy departure from conventional modeling practices, where the human population has traditionally been treated as a uniform entity in prior studies. Consequently, our research underscores the paramount significance of incorporating awareness dynamics into the modeling of the monkeypox disease. Furthermore, we have introduced vaccination as a pivotal control parameter within our model to assess its efficacy in curbing the transmission of the monkeypox disease. Within this novel framework, we have meticulously crafted a compartmental model for monkeypox, employing fractional derivatives as a mathematical foundation. Our research extends beyond model formulation, delving deeply into the intricate dynamics of this infectious disease. Our exploration advances the understanding of monkeypox transmission dynamics by introducing these distinctive elements, which have hitherto been underrepresented in the literature.

The structural organization of the remaining part of this paper is as follows. In Sec. II, we expound upon the foundational tenets and essential concepts of fractional calculus, providing the necessary theoretical underpinning for our subsequent modeling endeavors. Section III encompasses the pivotal formulation of our model. Initially, we construct a deterministic monkeypox model utilizing a system of ordinary differential equations. Subsequently, we refine this model by introducing fractional-order differential equations, employing the Atangana–Baleanu derivative framework [as outlined in the works of Baleanu (2018), Atangana and Owolabi (2018), and Danane (2020)]. Section IV delves into a comprehensive exploration of the fundamental properties intrinsic to our model. These properties encompass the determination of the invariant region, the establishment of conditions guaranteeing the existence and uniqueness of solutions, and the assurance of compactness and continuity of the model. Section V is dedicated to an exhaustive analysis of the model's equilibrium points, computation of both basic and effective reproduction numbers, and a detailed stability analysis. In Sec. VI, we elucidate the intricacies of model calibration, followed by an in-depth sensitivity analysis presented in Sec. VII. Section VIII offers a comprehensive presentation of our numerical simulations, accompanied by graphical representations of our results, thus providing a visual dimension to our findings. Finally, in Sec. IX, we draw a definitive conclusion, summarizing the key outcomes and insights derived from our research endeavors.

Here, we shall discuss some fundamental definitions and concepts of fractional calculus (Abdeljawad , 2020; Kumar , 2022).

Definition 1. Let us consider a function h(t) and a scalar q such that 0 < q < 1; then, we can define the Atangana–Baleanu derivative in the Caputo sense (ABC derivative) in the following manner:
ABC 0 D t q = G ( q ) 1 q 0 t E q [ q ( t z ) q 1 q ] h ( z ) d z .

The normalization function is denoted by G ( q ), and it satisfies G ( 0 ) = G ( 1 ) = 1. Here, E q denotes the Mittag–Leffler function of a single parameter defined in the following manner: E q ( u ) = k = 0 u k Γ ( q k + 1 ), where 0 < q < 1.

Definition 2. Let us consider a function h(t) and a scalar q such that 0 < q < 1, we can define the Atangana–Baleanu integral that corresponds to the ABC derivative in the following manner:
A B 0 I t q = 1 q G ( q ) h ( t ) + q G ( q ) Γ ( q ) 0 t f ( z ) ( t z ) q 1 d z .
Definition 3. Here, we shall give the definition of the Laplace transform of the Atangana–Baleanu fractional derivative taken in the Caputo sense.
L { 0 ABC D t q h ( t ) } ( s ) = G ( q ) 1 q s q L { h ( t ) } ( s ) s q 1 h ( 0 ) s q + q 1 q .

Here, we shall develop a compartmental model to describe the dynamics of monkeypox disease. For this, we have divided the human population into nine compartments, namely, susceptible unaware humans SU, susceptible aware humans SA, exposed humans EH, asymptomatically infected humans IA, symptomatically infected humans IH, hospitalized humans HH, recovered humans RH, unaware humans that took the vaccination VU, and vaccinated aware humans VH. Now coming to the rodent population, we have subcategorized it into two compartments, which are susceptible vector (rodents) SV and infected vector (rodents) IV. So, we have constructed our model with a total of 11 state variables.

Humans and rodents are recruited at a constant rate of πH and πV, respectively, into the population. λ denotes the fraction of susceptible humans that are aware. β , β 1 , β 2, and β3 denote the effective contact rate between rodents and unaware humans, unaware humans and other humans, aware humans and other humans, and infected rodents and other symptomatic rodents, respectively. Here, ρ1 denotes the movement rate of asymptomatically infected humans to either the hospitalized class or recovered class. θ2 denotes fraction of asymptomatically infected humans that are hospitalized and the remaining gets recovered. η1 denotes the rate of infection of exposed humans. The fraction of these infected humans that become symptomatic is given by θ1, and the remaining ones move into the asymptomatically infected class. Next, α1 gives the rate of treatment given to symptomatic infected human population of which a fraction η2 gets recovered and the remaining are further hospitalized. The rate of vaccination given to the human population is denoted by u1 and the efficiency rate of the vaccine by σ1. θ3 is the fraction of aware humans that are vaccinated. r denotes the rate of recovery of hospitalized humans. μh and μv are the natural mortality rates of the human and vector population, respectively, whereas δ and delta1 denote the disease induced death rate of the human and the rodent population, respectively.

The deterministic model employed to study the entire dynamics of the transmission of the monkeypox disease for our current research is represented by the system of ordinary differential equations (ODEs) in Eq. (1) as follows:
S U = λ π H β I V S U N H β 1 I H S U N H μ h S U u 1 σ 1 ( 1 θ 3 ) S U , S A = ( 1 λ ) π H β 2 I H S A N H μ h S A u 1 σ 1 θ 3 S A , E H = β I V S U N H + β 1 I H S U N H + β 2 I H S A N H ( η 1 + μ h ) E H , I A = ( 1 θ 1 ) η 1 E H ( ρ 1 + μ h ) I A , I H = θ 1 η 1 E H ( α 1 + δ + μ h ) I H , H H = ( 1 η 2 ) α 1 I H + θ 2 ρ 1 I A ( μ h + r ) H H , R H = η 2 α 1 I H + ( 1 θ 2 ) ρ 1 I A + r H H μ h R H , V U = u 1 σ 1 ( 1 θ 3 ) S U μ h V U , V A = u 1 σ 1 θ 3 S A μ H V A , S V = π V β 3 S V I V N V μ v S V , I V = β 3 S V I V N V ( μ v + δ 1 ) I V .
(1)
Here, NH is the total human population and NV is the total vector (rodent) population. They satisfy the following equations:
N H = Λ H μ h N H δ I H , N V = Λ V μ v N V δ 1 I V .
(2)
Next, we will modify the given deterministic system (1) in frame of ABC derivative in the following manner:
ABC 0 D t q S U = λ π H β I V S U N H β 1 I H S U N H μ h S U u 1 σ 1 ( 1 θ 3 ) S U , ABC 0 D t q S A = ( 1 λ ) π H β 2 I H S A N H μ h S A u 1 σ 1 θ 3 S A , ABC 0 D t q E H = β I V S U N H + β 1 I H S U N H + β 2 I H S A N H ( η 1 + μ h ) E H , ABC 0 D t q I A = ( 1 θ 1 ) η 1 E H ( ρ 1 + μ h ) I A , ABC 0 D t q I H = θ 1 η 1 E H ( α 1 + δ + μ h ) I H , ABC 0 D t q H H = ( 1 η 2 ) α 1 I H + θ 2 ρ 1 I A ( μ h + r ) H H , ABC 0 D t q R H = η 2 α 1 I H + ( 1 θ 2 ) ρ 1 I A + r H H μ h R H , ABC 0 D t q V U = u 1 σ 1 ( 1 θ 3 ) S U μ h V U , ABC 0 D t q V A = u 1 σ 1 θ 3 S A μ H V A , ABC 0 D t q S V = π V β 3 S V I V N V μ v S V , ABC 0 D t q I V = β 3 S V I V N V ( μ v + δ 1 ) I V .
(3)

The parameters of the above model (2) are assumed to be strictly non-negative. Also, since this model takes into consideration the living population, the state variables are also assumed to be non-negative at the initial time t = 0.

Lemma 1. The initial data are assumed as follows: P H ( 0 ) 0 , P V ( 0 ) 0, where P H ( t ) and P V ( t ) are given by
P H ( t ) = ( S U ( t ) , S A ( t ) , E H ( t ) , I A ( t ) , I H ( t ) , H H ( t ) , R H ( t ) , V U ( t ) , V A ( t ) ) , P V ( t ) = ( S V ( t ) , I V ( t ) ) .

So, under the ABC derivative, the two closed sets P H ( 0 ) 0 , P V ( 0 ) 0 are positively invariant for the model (2).

Proof. To obtain the dynamics of the model (2) for the human and vector population, we shall sum up the compartmental equations for these two populations separately in the following manner (Atangana and Owolabi, 2018):
ABC 0 D t q N H ( t ) = ABC 0 D t q S U ( t ) + ABC 0 D t q S A ( t ) + ABC 0 D t q E H ( t ) + ABC 0 D t q I A ( t ) + ABC 0 D t q I H ( t ) + ABC 0 D t q H H ( t ) + ABC 0 D t q R H ( t ) + ABC 0 D t q V U ( t ) + ABC 0 D t q V A ( t ) .
So,
ABC 0 D t q N V ( t ) = ABC 0 D t q S V ( t ) + ABC 0 D t q I V ( t ) .
Now, for the human population, we get
ABC 0 D t q N H ( t ) = π H μ h N H δ I H , ABC 0 D t q N H ( t ) π H μ h N H .
(4)
On Eq. (4), we shall apply Laplace transformation to obtain
L { 0 ABC D t q N H ( t ) + μ h N H } L { π H } , N H ( t ) ( N ( q ) N ( q ) + ( 1 q ) μ h N ( 0 ) + ( 1 q ) π H N ( q ) + ( 1 q ) μ h ) E q , 1 ( η 1 t q ) + q π H N ( q ) + ( 1 q ) μ h E q , q + 1 ( η 1 t q ) ) .
(5)
Here, the Mittag–Leffler (ML) function with two parameters q and p is denoted by E q , p and η 1 = q μ h N ( q ) + ( 1 q ) μ h. Considering the asymptomatic nature of the ML function, we deduce the following result:
lim t N H ( t ) π H μ h ; for all t > 0 .
In a similar manner, we can deduce
lim t N V ( t ) π V μ v ; for all t > 0 .

Now, for the given initial conditions and t > 0, the solutions of the model (3) lying in PH and PV continue to stay in PH and PV, respectively. This shows that the regions PH and PV will serve as an attractor for all the solutions belonging to R 9 + and R 2 + and hence is positively invariant.

Here, we will try to analyze the existence and uniqueness property of the solution for the given fractional-order system (3). In order to demonstrate this, we shall consider on the interval Y = [ 0 , T ], a Banach space M ( Z ) comprising of all real-valued functions that are continuous having supremum, which is defined as follows:
M K = K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) X K ( Z ) ,
with the norm
| | S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V | | = | | S U | | + | | S A | | + | | E H | | + | | I A | | + | | I H | | + | | H H | | + | | R H | | + | | V U | | + | | V A | | + | | S V | | + | | I V | |, where
| | S U | | = sup t Z | S U ( t ) | , | | S A | | = sup t Z | S A ( t ) | , | | E H | | = sup t Z | E A ( t ) | , | | I A | | = sup t Z | I A ( t ) | , | | I H | | = sup t Z | I H ( t ) | , | | H H | | = sup t Z | H H ( t ) | , | | R H | | = sup t Z | R H ( t ) | , | | V U | | = sup t Z | V U ( t ) | , | | V A | | = sup t Z | V A ( t ) | , | | S V | | = sup t Z | S V ( t ) | , | | I V | | = sup t Z | I V ( t ) | .
According to the definition of AB fractional integral, we will now integrate all the equations of the model (4) to get the following results:
S U ( t ) S U ( 0 ) = A B 0 I t q { λ π H β I V S U N H β 1 I H S U N H μ h S U u 1 σ 1 ( 1 θ 3 ) S U } , S A ( t ) S A ( 0 ) = A B 0 I t q { ( 1 λ ) π H β 2 I H S A N H μ h S A u 1 σ 1 θ 3 S A } , E H ( t ) E H ( 0 ) = A B 0 I t q { β I V S U N H + β 1 I H S U N H + β 2 I H S A N H ( η 1 + μ h ) E H } , I A ( t ) I A ( 0 ) = A B 0 I t q { ( 1 θ 1 ) η 1 E H ( ρ 1 I A + μ h ) I A } , I H ( t ) I H ( 0 ) = A B 0 I t q { θ 1 η 1 E H ( α 1 + δ + μ h ) I H } , H H ( t ) H H ( 0 ) = A B 0 I t q { ( 1 η 2 ) α 1 I H + θ 2 ρ 1 I A ( μ h + r ) H H } , R H ( t ) R H ( 0 ) = A B 0 I t q { η 2 α 1 I H + ( 1 θ 2 ) ρ 1 I A + r H H μ h R H } , V U ( t ) V U ( 0 ) = A B 0 I t q { u 1 σ 1 ( 1 θ 3 ) S U μ h V U } , V A ( t ) V A ( 0 ) = A B 0 I t q { u 1 σ 1 θ 3 S A μ H V A } , S V ( t ) S V ( 0 ) = A B 0 I t q { π V β 3 S V I V N V μ v S V } , I V ( t ) I V ( 0 ) = A B 0 I t q { β 3 S V I V N V μ v I V } .
(6)
According to the axiomatic definition of the Atangana–Baleanu fractional integral (4), we obtain
S U ( t ) S U ( 0 ) = 1 q G ( q ) Q 1 ( t , S U ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 1 ( z , S U ( z ) ) d z , S A ( t ) S A ( 0 ) = 1 q G ( q ) Q 2 ( t , S A ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 2 ( z , S A ( z ) ) d z , E H ( t ) E H ( 0 ) = 1 q G ( q ) Q 3 ( t , E H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 3 ( z , E H ( z ) ) d z , I A ( t ) I A ( 0 ) = 1 q G ( q ) Q 4 ( t , I A ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 4 ( z , I A ( z ) ) d z , I H ( t ) I H ( 0 ) = 1 q G ( q ) Q 5 ( t , I H ( t ) ) + q Q ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 5 ( z , I H ( z ) ) d z , H H ( t ) H H ( 0 ) = 1 q G ( q ) Q 6 ( t , H H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 6 ( z , H H ( z ) ) d z , R H ( t ) R H ( 0 ) = 1 q G ( q ) Q 7 ( t , R H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 7 ( z , R H ( z ) ) d z , V U ( t ) V U ( 0 ) = 1 q G ( q ) Q 8 ( t , V U ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 8 ( z , V U ( z ) ) d z , V A ( t ) V A ( 0 ) = 1 q G ( q ) Q 9 ( t , V A ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 9 ( z , V A ( z ) ) d z , S V ( t ) S V ( 0 ) = 1 q G ( q ) Q 10 ( t , S V ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 10 ( z , S V ( z ) ) d z , I V ( t ) I V ( 0 ) = 1 q G ( q ) Q 11 ( t , I V ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 11 ( z , I V ( z ) ) d z ,
(7)
where
Q 1 ( t , S U ( t ) ) = λ π H β I V S U N H β 1 I H S U N H μ h S U u 1 σ 1 ( 1 θ 3 ) S U , Q 2 ( t , S A ( t ) ) = ( 1 λ ) π H β 2 I H S A N H μ h S A u 1 σ 1 θ 3 S A , Q 3 ( t , E H ( t ) ) = β I V S U N H + β 1 I H S U N H + β 2 I H S A N H ( η 1 + μ h ) E H , Q 4 ( t , I A ( t ) ) = ( 1 θ 1 ) η 1 E H ( ρ 1 I A + μ h ) I A , Q 5 ( t , I H ( t ) ) = θ 1 η 1 E H ( α 1 + δ + μ h ) I H , Q 6 ( t , H H ( t ) ) = ( 1 η 2 ) α 1 I H + θ 2 ρ 1 I A ( μ h + r ) H H , Q 7 ( t , R H ( t ) ) = η 2 α 1 I H + ( 1 θ 2 ) ρ 1 I A + r H H μ h R H , Q 8 ( t , V U ( t ) ) = u 1 σ 1 ( 1 θ 3 ) S U μ h V U , Q 9 ( t , V A ( t ) ) = u 1 σ 1 θ 3 S A μ H V A , Q 10 ( t , S V ( t ) ) = π V β 3 S V I V N V μ v S V , Q 11 ( t , I V ( t ) ) = β 3 S V I V N V μ v I V .
A necessary as well as sufficient condition for the following functions S U ( t ) , S A ( t ) , E H ( t ) , I A ( t ) , I H ( t ) , H H ( t ) , R H ( t ) , V U ( t ) , V A ( t ) , S V ( t ) , I V ( t ) to have an upper bound is that Lipschitz condition has to be satisfied for the expressions:
Q 1 ( t , S U ( t ) ) , Q 2 ( t , S A ( t ) ) , Q 3 ( t , E H ( t ) ) , Q 4 ( t , I A ( t ) ) , Q 5 ( t , I H ( t ) ) , Q 6 ( t , H H ( t ) ) , Q 7 ( t , R H ( t ) ) , Q 8 ( t , V U ( t ) ) , Q 9 ( t , V A ( t ) ) , Q 10 ( t , S V ( t ) ) , Q 11 ( t , I V ( t ) ) .
Now, taking two functions, S U ( t ) and S ̂ U then
| | Q 1 ( t , S U ( t ) ) Q 1 ( t , S ̂ U ( t ) ) | | = | | ( β I V N H + β 1 I H N H + u 1 σ 1 ( 1 θ 3 ) + μ h ) ( S U S ̂ U ) | | | | ( β I V N H + β 1 I H N H + u 1 σ 1 ( 1 θ 3 ) + μ h ) | | | | ( S U S ̂ U ) | | c 1 | | ( S U S ̂ U ) | | ,
where c 1 = | | ( β I V N H + β 1 I H N H + u 1 σ 1 ( 1 θ 3 ) + μ h ) | |.
Therefore, we obtain the following result:
| | Q 1 ( t , S U ( t ) ) Q 1 ( t , S ̂ U ( t ) ) | | c 1 | | ( S U S ̂ U ) | | .
(8)
For the other state variables, doing a similar kind of simplification as done for SU, we will get the following results:
| | Q 2 ( t , S A ( t ) ) Q 2 ( t , S ̂ A ( t ) ) | | c 2 | | ( S A S ̂ A ) | | , | | Q 3 ( t , E H ( t ) ) Q 3 ( t , E ̂ H ( t ) ) | | c 3 | | ( E H E ̂ A ) | | , | | Q 4 ( t , I A ( t ) ) Q 4 ( t , I ̂ A ( t ) ) | | c 4 | | ( I A I ̂ A ) | | , | | Q 5 ( t , I H ( t ) ) Q 5 ( t , I ̂ H ( t ) ) | | c 5 | | ( I H I ̂ H ) | | , | | Q 6 ( t , H H ( t ) ) Q 6 ( t , H ̂ H ( t ) ) | | c 6 | | ( H H H ̂ H ) | | , | | Q 7 ( t , R H ( t ) ) Q 7 ( t , R ̂ H ( t ) ) | | c 7 | | ( R H R ̂ H ) | | , | | Q 8 ( t , V U ( t ) ) Q 8 ( t , V ̂ U ( t ) ) | | c 8 | | ( V U V ̂ U ) | | , | | Q 9 ( t , V A ( t ) ) Q 9 ( t , V ̂ A ( t ) ) | | c 9 | | ( V A V ̂ A ) | | , | | Q 10 ( t , S V ( t ) ) Q 10 ( t , S ̂ V ( t ) ) | | c 10 | | ( S V S ̂ V ) | | , | | Q 11 ( t , I V ( t ) ) Q 11 ( t , I ̂ V ( t ) ) | | c 11 | | ( I V I ̂ V ) | | .
(9)

Thus, observing the system of equations we got above, we can conclude that Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , Q 6 , Q 7 , Q 8 , Q 9, Q10, and Q11 do satisfy the Lipschitz condition. Here, c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , c 7 , c 8 , c 9 , c 10 , c 11 represent the corresponding Lipschitz constants for the equations given above.

On the system of Eqs. (6), we apply the principle of recursion to obtain
S U k ( t ) = 1 q G ( q ) Q 1 ( t , S U k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 1 ( z , S U k 1 ( z ) ) d z + S U ( 0 ) , S A k ( t ) = 1 q G ( q ) Q 2 ( t , S A k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 2 ( z , S A k 1 ( z ) ) d z + S A ( 0 ) , E H k ( t ) = 1 q G ( q ) Q 3 ( t , E H k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 3 ( z , E H k 1 ( z ) ) d z + E H ( 0 ) , I A k ( t ) = 1 q G ( q ) Q 4 ( t , I A k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 4 ( z , I A k 1 ( z ) ) d z + I A ( 0 ) , I H k ( t ) = 1 q G ( q ) Q 5 ( t , I H k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 5 ( z , I H k 1 ( z ) ) d z + I H ( 0 ) , H H k ( t ) = 1 q G ( q ) Q 6 ( t , H H k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 6 ( z , H H k 1 ( z ) ) d z + H H ( 0 ) , R H k ( t ) = 1 q G ( q ) Q 7 ( t , R H k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 7 ( z , R H k 1 ( z ) ) d z + R H ( 0 ) , V U k ( t ) = 1 q G ( q ) Q 8 ( t , V U k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 8 ( z , V U k 1 ( z ) ) d z + V U ( 0 ) , V A k ( t ) = 1 q G ( q ) Q 9 ( t , V A k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 9 ( z , V A k 1 ( z ) ) d z + V A ( 0 ) , S V k ( t ) = 1 q G ( q ) Q 10 ( t , S V k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 10 ( z , S V k 1 ( z ) ) d z + S V ( 0 ) , I V k ( t ) = 1 q G ( q ) Q 11 ( t , I V k 1 ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 11 ( z , I V k 1 ( z ) ) d z + I V ( 0 ) .
The initial conditions are considered as follows:
S U 0 ( t ) = S U ( 0 ) , S A 0 ( t ) = S A ( 0 ) , E H 0 ( t ) = E H ( 0 ) , I A 0 ( t ) = I A ( 0 ) , I H 0 ( t ) = I H ( 0 ) , H H 0 ( t ) = H H ( 0 ) , R H 0 ( t ) = R H ( 0 ) , V U 0 ( t ) = V U ( 0 ) , V A 0 ( t ) = V A ( 0 ) , S V 0 ( t ) = S V ( 0 ) , I V 0 ( t ) = I V ( 0 ) .
Taking the difference of successive terms, we obtain
Ω S U , k ( t ) = S U k ( t ) S U k 1 ( t ) , Ω S U , k ( t ) = 1 q G ( q ) [ Q 1 ( t , S U k 1 ( t ) ) Q 1 ( t , S U k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 1 ( z , S U k 1 ( z ) ) Q 1 ( t , S U k 2 ( z ) ) d z , Ω S A , k ( t ) = S A k ( t ) S A k 1 ( t ) , Ω S A , k ( t ) = 1 q G ( q ) [ Q 2 ( t , S A k 1 ( t ) ) Q 2 ( t , S A k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 2 ( z , S A k 1 ( z ) ) Q 2 ( t , S A k 2 ( z ) ) d z , Ω E H , k ( t ) = E H k ( t ) E H k 1 ( t ) , Ω E H , k ( t ) = 1 q G ( q ) [ Q 3 ( t , E H k 1 ( t ) ) Q 3 ( t , E H k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 3 ( z , E H k 1 ( z ) ) Q 3 ( t , E H k 2 ( z ) ) d z , Ω I A , k ( t ) = I A k ( t ) I A k 1 ( t ) , Ω I A , k ( t ) = 1 q G ( q ) [ Q 4 ( t , I A k 1 ( t ) ) Q 4 ( t , I A k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 4 ( z , I A k 1 ( z ) ) Q 4 ( t , I A k 2 ( z ) ) d z , Ω I H , k ( t ) = I H k ( t ) I H k 1 ( t ) , Ω I H , k ( t ) = 1 q G ( q ) [ Q 5 ( t , I H k 1 ( t ) ) Q 5 ( t , I H k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 5 ( z , I H k 1 ( z ) ) Q 5 ( t , I H k 2 ( z ) ) d z , Ω H H , k ( t ) = H H k ( t ) H H k 1 ( t ) , Ω H H , k ( t ) = 1 q G ( q ) [ Q 6 ( t , H H k 1 ( t ) ) Q 6 ( t , H H k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 6 ( z , H H k 1 ( z ) ) Q 6 ( t , H H k 2 ( z ) ) d z , Ω R H , k ( t ) = R H k ( t ) R H k 1 ( t ) , Ω R H , k ( t ) = 1 q G ( q ) [ Q 7 ( t , R H k 1 ( t ) ) Q 7 ( t , R H k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 7 ( z , R H k 1 ( z ) ) Q 7 ( t , R H k 2 ( z ) ) d z , Ω V U , k ( t ) = V U k ( t ) V U k 1 ( t ) , Ω V U , k ( t ) = 1 q G ( q ) [ Q 8 ( t , V U k 1 ( t ) ) Q 8 ( t , V U k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 8 ( z , V U k 1 ( z ) ) Q 8 ( t , V U k 2 ( z ) ) d z , Ω V A , k ( t ) = V A k ( t ) V A k 1 ( t ) , Ω V A , k ( t ) = 1 q G ( q ) [ Q 9 ( t , V A k 1 ( t ) ) Q 9 ( t , V A k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 9 ( z , V A k 1 ( z ) ) Q 9 ( t , V A k 2 ( z ) ) d z , Ω S V , k ( t ) = S V k ( t ) S V k 1 ( t ) , Ω S V , k ( t ) = 1 q G ( q ) [ Q 10 ( t , S V k 1 ( t ) ) Q 10 ( t , S V k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 10 ( z , S V k 1 ( z ) ) Q 10 ( t , S V k 2 ( z ) ) d z , Ω I V , k ( t ) = I V k ( t ) I V k 1 ( t ) , Ω I V , k ( t ) = 1 q G ( q ) [ Q 11 ( t , I V k 1 ( t ) ) Q 11 ( t , I V k 2 ( t ) ) ] + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 11 ( z , I V k 1 ( z ) ) Q 11 ( t , I V k 2 ( z ) ) d z .
After getting the above results, we conclude that
S U k ( t ) = r = 0 k Ω S U , r ( t ) , S A k ( t ) = r = 0 k Ω S A , r ( t ) , E H k ( t ) = r = 0 k Ω E H , r ( t ) , I A k ( t ) = r = 0 k Ω I A , r ( t ) , I H k ( t ) = r = 0 k Ω I H , r ( t ) , H H k ( t ) = r = 0 k Ω H H , r ( t ) , R H k ( t ) = r = 0 k Ω R H , r ( t ) , V U k ( t ) = r = 0 k Ω V U , r ( t ) , V A k ( t ) = r = 0 k Ω V A , r ( t ) , S V k ( t ) = r = 0 k Ω S V , r ( t ) , I V k ( t ) = r = 0 k Ω I V , r ( t ) .
By using the property of recursion, we have
Ω S U , k 1 ( t ) = S U k 1 ( t ) S U k 2 ( t ) , Ω S A , k 1 ( t ) = S A k 1 ( t ) S A k 2 ( t ) , Ω E H , k 1 ( t ) = E H k 1 ( t ) E H k 2 ( t ) , Ω I A , k 1 ( t ) = I A k 1 ( t ) I A k 2 ( t ) , Ω I H , k 1 ( t ) = I H k 1 ( t ) I H k 2 ( t ) , Ω H H , k 1 ( t ) = H H k 1 ( t ) H H k 2 ( t ) , Ω R H , k 1 ( t ) = R H k 1 ( t ) R H k 2 ( t ) , Ω V U , k 1 ( t ) = V U k 1 ( t ) V U k 2 ( t ) , Ω V A , k 1 ( t ) = V A k 1 ( t ) V A k 2 ( t ) , Ω S V , k 1 ( t ) = S V k 1 ( t ) S V k 2 ( t ) , Ω I V , k 1 ( t ) = I V k 1 ( t ) I V k 2 ( t ) .
Using the above set of equations along with Eqs. (8) and (9), we get to the following results:
| | Ω S U , k ( t ) | | 1 q G ( q ) c 1 | | Ω S U , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 1 0 t ( t z ) q 1 | | Ω S U , k 1 ( z ) | | d z , | | Ω S A , k ( t ) | | 1 q G ( q ) c 2 | | Ω S A , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 2 0 t ( t z ) q 1 | | Ω S A , k 1 ( z ) | | d z , | | Ω E H , k ( t ) | | 1 q G ( q ) c 3 | | Ω E H , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 3 0 t ( t z ) q 1 | | Ω E H , k 1 ( z ) | | d z , | | Ω I A , k ( t ) | | 1 q G ( q ) c 4 | | Ω I A , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 4 0 t ( t z ) q 1 | | Ω I A , k 1 ( z ) | | d z , | | Ω I H , k ( t ) | | 1 q G ( q ) c 5 | | Ω I H , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 5 0 t ( t z ) q 1 | | Ω I H , k 1 ( z ) | | d z , | | Ω H H , k ( t ) | | 1 q G ( q ) c 6 | | Ω H H , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 6 0 t ( t z ) q 1 | | Ω H H , k 1 ( z ) | | d z , | | Ω R H , k ( t ) | | 1 q G ( q ) c 7 | | Ω R H , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 7 0 t ( t z ) q 1 | | Ω R H , k 1 ( z ) | | d z , | | Ω V U , k ( t ) | | 1 q G ( q ) c 8 | | Ω V U , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 8 0 t ( t z ) q 1 | | Ω V U , k 1 ( z ) | | d z , | | Ω V A , k ( t ) | | 1 q G ( q ) c 9 | | Ω V A , k 1 ( t ) | | + q G ( q ) Γ ( q ) c 9 0 t ( t z ) q 1 | | Ω V A , k 1 ( z ) | | d z , | | Ω S V , k ( t ) | | 1 q G ( q ) c 10 | | Ω S V , k 1 ( t ) | | + q G ( q ) Γ ( p ) c 10 0 t ( t z ) q 1 | | Ω S V , k 1 ( z ) | | d z , | | Ω I V , k ( t ) | | 1 q G ( q ) c 11 | | Ω I V , k 1 ( t ) | | + q G ( q ) Γ ( p ) c 11 0 t ( t z ) q 1 | | Ω I V , k 1 ( z ) | | d z .
(10)
Theorem. The solution of the monkeypox fractional compartmental model (2) will be unique in nature if for 0 t T, it will satisfy the condition as follows:
1 q G ( q ) r i + x q G ( q ) Γ ( q ) r i < 1 ,
for j = 1 , 2 , 3 , , 11
Proof. Previously, we had established the boundedness of S U ( t ) , S A ( t ) , E H ( t ) , I A ( t ) , I H ( t ) , H H ( t ) , R H ( t ) , V U ( t ) , V A ( t ) , S V ( t ), and I V ( t ). We have also observed from Eqs. (4) and (6) that Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8, Q9, Q10, and Q11 evidently satisfy the Lipschitz criterion. By employing the set of equations given in Eq. (6) along with a recursive formula, we get
| | Ω S U , k ( t ) | | | | S U ( 0 ) | | [ 1 q G ( q ) c 1 + x q G ( q ) Γ ( q ) c 1 ] k , | | Ω S A , k ( t ) | | | | S A ( 0 ) | | [ 1 q G ( q ) c 2 + x q G ( q ) Γ ( q ) c 2 ] k , | | Ω E H , k ( t ) | | | | E H ( 0 ) | | [ 1 q G ( q ) c 3 + x q G ( q ) Γ ( q ) c 3 ] k , | | Ω I A , k ( t ) | | | | I A ( 0 ) | | [ 1 q G ( q ) c 4 + x q G ( q ) Γ ( q ) c 4 ] k , | | Ω I H , k ( t ) | | | | I H ( 0 ) | | [ 1 q G ( q ) c 5 + x q G ( q ) Γ ( q ) c 5 ] k , | | Ω H H , k ( t ) | | | | H H ( 0 ) | | [ 1 q G ( q ) c 6 + x q G ( q ) Γ ( q ) c 6 ] k , | | Ω R H , k ( t ) | | | | R H ( 0 ) | | [ 1 q G ( q ) c 7 + x q G ( q ) Γ ( q ) c 7 ] k , | | Ω V U , k ( t ) | | | | V U ( 0 ) | | [ 1 q G ( q ) c 8 + x q G ( q ) Γ ( q ) c 8 ] k , | | Ω V A , k ( t ) | | | | V A ( 0 ) | | [ 1 q G ( q ) c 9 + x q G ( q ) Γ ( q ) c 9 ] k , | | Ω S V , k ( t ) | | | | S V ( 0 ) | | [ 1 q G ( q ) c 10 + x q G ( q ) Γ ( q ) c 10 ] k , | | Ω I V , k ( t ) | | | | I V ( 0 ) | | [ 1 q G ( q ) c 11 + x q G ( q ) Γ ( q ) c 11 ] k .
(11)
Therefore, we can say that the above-specified sequences do exist and satisfy the following:
| | Ω S U , k ( t ) | | 0 , | | Ω S A , k ( t ) | | 0 , | | Ω E H , k ( t ) | | 0 , | | Ω I A , k ( t ) | | 0 , | | Ω I H , k ( t ) | | 0 , | | Ω H H , k ( t ) | | 0 , | | Ω R H , k ( t ) | | 0 , | | Ω V U , k ( t ) | | 0 , | | Ω V A , k ( t ) | | 0 , | | Ω S V , k ( t ) | | 0 , | | Ω I V , k ( t ) | | 0.
For any value of j, in view of triangle inequality and using the set of Eqs. (6), we deduce that
| | S U k + j ( t ) S H k ( t ) | | s = k + 1 k + j U 1 s = U 1 k + 1 U 1 k + j + 1 1 U 1 , | | S A k + j ( t ) E H k ( t ) | | s = k + 1 k + j U 2 s = U 2 k + 1 U 2 k + j + 1 1 U 2 , | | E H k + j ( t ) I A k ( t ) | | s = k + 1 k + j U 3 s = U 3 k + 1 U 3 k + j + 1 1 U 3 , | | I A k + j ( t ) I H k ( t ) | | s = k + 1 k + j U 4 s = U 4 k + 1 U 4 k + j + 1 1 U 4 , | | I H k + j ( t ) T H k ( t ) | | s = k + 1 k + j U 5 s = U 5 k + 1 U 5 k + j + 1 1 U 5 , | | H H k + j ( t ) P H k ( t ) | | s = k + 1 k + j U 6 s = U 6 k + 1 U 6 k + j + 1 1 U 6 , | | R H k + j ( t ) R H k ( t ) | | s = k + 1 k + j U 7 s = U 7 k + 1 U 7 k + j + 1 1 U 7 , | | V U k + j ( t ) S R k ( t ) | | s = k + 1 k + j U 8 s = U 8 k + 1 U 8 k + j + 1 1 U 8 , | | V A k + j ( t ) I R k ( t ) | | s = k + 1 k + j U 9 s = U 9 k + 1 U 9 k + j + 1 1 U 9 , | | S V k + j ( t ) S V k ( t ) | | s = k + 1 k + j U 10 s = U 10 k + 1 U 10 k + j + 1 1 U 10 , | | I V k + j ( t ) I V k ( t ) | | s = k + 1 k + j U 11 s = U 11 k + 1 U 11 k + j + 1 1 U 11 .

Here, the symbol U i is given by the expression U i = 1 q G ( q ) r i + x q G ( q ) Γ ( q ) r i < 1 according to the hypothesis.

So, from the above results, we can observe that the sequences S U k , S A k , E H k , I A k , I H k , H H k , R H k , V U k , V A k , S V k , I V k are Cauchy sequences in the Banach space K(Z) that we had defined earlier. So, by nature, they are uniformly convergent, and hence, as we apply on these sequences the limit as n we will obtain exactly one solution, which establishes the uniqueness criterion.

Therefore, through this we can finally conclude that the solution for the given system of fractional differential equations for the monkeypox model (3) exists and is unique.

Let us first define an operator F as follows: F = ( F 1 , F 2 , F 3 , F 4 , F 5 , F 6 , F 7 , F 8 , F 9 , F 10 , F 11 ) :

Here, F 1 , F 2 , F 3 , F 4 , F 5 , F 6 , F 7 , F 8 , F 9 , F 10 , F 11 are defined in the following manner:
F 1 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = S U ( 0 ) + 1 q G ( q ) Q 1 ( t , S U ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 1 ( z , S U ( z ) ) d z = A 1 + B 1 , F 2 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = S A ( 0 ) + 1 q G ( q ) Q 2 ( t , S A ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 2 ( z , S A ( z ) ) d z = A 2 + B 2 , F 3 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = E H ( 0 ) + 1 q G ( q ) Q 3 ( t , E H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 3 ( z , E H ( z ) ) d z = A 3 + B 3 , F 4 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = I A ( 0 ) + 1 q G ( q ) Q 4 ( t , I A ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 4 ( z , I A ( z ) ) d z = A 4 + B 4 , F 5 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = I H ( 0 ) + 1 q G ( q ) Q 5 ( t , I H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 5 ( z , I H ( z ) ) d z = A 5 + B 5 , F 6 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = H H ( 0 ) + 1 q G ( q ) Q 6 ( t , H H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 6 ( z , H H ( z ) ) d z = A 6 + B 6 , F 7 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = R H ( 0 ) + 1 q G ( q ) Q 7 ( t , R H ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 7 ( z , R H ( z ) ) d z = A 7 + B 7 , F 8 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = V U ( 0 ) + 1 q G ( q ) Q 8 ( t , V U ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 8 ( z , V U ( z ) ) d z = A 8 + B 8 , F 9 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = V A ( 0 ) + 1 q G ( q ) Q 9 ( t , V A ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 9 ( z , V A ( z ) ) d z = A 9 + B 9 , F 10 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = S V ( 0 ) + 1 q G ( q ) Q 10 ( t , S V ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 10 ( z , S V ( z ) ) d z = A 10 + B 10 , F 11 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) = I V ( 0 ) + 1 q G ( q ) Q 11 ( t , I V ( t ) ) + q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 11 ( z , I V ( z ) ) d z = A 11 + B 11 ,
where
A 1 = S U ( 0 ) + 1 q G ( q ) Q 1 ( t , S U ( t ) ) and B 1 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 1 ( z , S U ( z ) ) d z ,
A 2 = S A ( 0 ) + 1 q G ( q ) Q 2 ( t , S A ( t ) ) and B 2 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 2 ( z , S A ( z ) ) d z ,
A 3 = E H ( 0 ) + 1 q G ( q ) Q 3 ( t , E H ( t ) ) and B 3 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 3 ( z , E H ( z ) ) d z ,
A 4 = I A ( 0 ) + 1 q G ( q ) Q 4 ( t , I A ( t ) ) and B 4 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 4 ( z , I A ( z ) ) d z ,
A 5 = I H ( 0 ) + 1 q G ( q ) Q 5 ( t , I H ( t ) ) and B 5 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 5 ( z , I H ( z ) ) d z ,
A 6 = H H ( 0 ) + 1 q G ( q ) Q 6 ( t , H H ( t ) ) and B 6 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 6 ( z , H H ( z ) ) d z ,
A 7 = R H ( 0 ) + 1 q G ( q ) Q 7 ( t , R H ( t ) ) and B 7 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 7 ( z , R H ( z ) ) d z ,
A 8 = V U ( 0 ) + 1 q G ( q ) Q 8 ( t , V U ( t ) ) and B 8 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 8 ( z , V U ( z ) ) d z ,
A 9 = V A ( 0 ) + 1 q G ( q ) Q 9 ( t , V A ( t ) ) and B 9 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 9 ( z , V A ( z ) ) d z ,
A 10 = S V ( 0 ) + 1 q G ( q ) Q 10 ( t , S V ( t ) ) and B 10 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 10 ( z , S V ( z ) ) d z ,
A 11 = I V ( 0 ) + 1 q G ( q ) Q 11 ( t , I V ( t ) ) and B 11 = q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 11 ( z , I V ( z ) ) d z .
Now, we define two sets A and B as follows:
A = ( A 1 , A 2 , A 3 , A 4 , A 5 , A 6 , A 7 , A 8 , A 9 , A 10 , A 11 ) , B = ( B 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , B 9 , B 10 , B 11 ) .

We will now go on to prove the compactness and continuity of the set B = ( B 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , B 9 , B 10 , B 11 ) defined above.

The operators Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , Q 6 , Q 7 , Q 8 , Q 9 , Q 10 , Q 11 that we had defined previously are continuous in nature and possess a positively bounded modulus represented below. a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 , a 8 , a 9 , a 10 , a 11 and b 1 , b 2 , b 3 , b 4 , b 5 , b 6 , b 7 , b 8 , b 9 , b 10 , b 11 used in the expressions below are some real scalars,
| Q 1 ( t , S U ( t ) ) | a 1 | | S U | | + b 1 , | Q 2 ( t , S A ( t ) ) | a 2 | | S A | | + b 2 , | Q 3 ( t , E H ( t ) ) | a 3 | | E H | | + b 3 , | Q 4 ( t , I A ( t ) ) | a 4 | | I A | | + b 4 , | Q 5 ( t , T H ( t ) ) | a 5 | | I H | | + b 5 , | Q 6 ( t , H H ( t ) ) | a 6 | | H H | | + b 6 , | Q 7 ( t , R H ( t ) ) | a 7 | | R H | | + b 7 , | Q 8 ( t , V U ( t ) ) | a 8 | | V U | | + b 8 , | Q 9 ( t , V A ( t ) ) | a 9 | | V A | | + b 9 , | Q 10 ( t , S V ( t ) ) | a 10 | | S V | | + b 10 , | Q 11 ( t , I V ( t ) ) | a 11 | | I V | | + b 11 .
A closed subset S of is now taken into consideration in the form as follows:
S = { ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) : | | S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V | | R , R > 0 } .
Now, we will get for ( S H , E H , I A , I H , T H , P H , R H , S R , I R , S V , I V ) S, | | B 1 ( t , S U ( t ) ) | | = max t [ 0 , T ] | q G ( q ) Γ ( q ) 0 t ( t z ) q 1 Q 1 ( z , S U ( z ) ) d z | x q G ( q ) Γ ( q ) 0 t ( x z ) q 1 | Q 1 ( Z , S U ( z ) ) d z |, which gives
| | B 1 ( t , S U ( t ) ) | | x q G ( q ) Γ ( q ) ( a 1 R + b 1 ) .
In a similar manner, we obtain
| | B 2 ( t , S A ( t ) ) | | x q G ( q ) Γ ( q ) ( a 2 R + b 2 ) , | | B 3 ( t , E H ( t ) ) | | x q G ( q ) Γ ( q ) ( a 3 R + b 3 ) , | | B 4 ( t , I A ( t ) ) | | x q G ( q ) Γ ( q ) ( a 4 R + b 4 ) , | | B 5 ( t , I H ( t ) ) | | x q G ( q ) Γ ( q ) ( a 5 R + b 5 ) , | | B 6 ( t , H H ( t ) ) | | x q G ( q ) Γ ( q ) ( a 6 R + b 6 ) , | | B 7 ( t , R H ( t ) ) | | x q G ( q ) Γ ( q ) ( a 7 R + b 7 ) , | | B 8 ( t , V U ( t ) ) | | x q G ( q ) Γ ( q ) ( a 8 R + b 8 ) , | | B 9 ( t , V A ( t ) ) | | x q G ( q ) Γ ( q ) ( a 9 R + b 9 ) , | | B 10 ( t , S V ( t ) ) | | x q G ( q ) Γ ( q ) ( a 10 R + b 10 ) , | | B 11 ( t , I V ( t ) ) | | x q G ( q ) Γ ( q ) ( a 11 R + b 11 ) .
The maximum norm of | | B ( B 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , B 9 , B 10 , B 11 ) | | is given by
| | B ( B 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , B 9 , B 10 , B 11 ) | | x q G ( q ) Γ ( q ) [ ( a 1 + a 2 + a 3 + a 4 + a 5 + a 6 + a 7 + a 8 + a 9 + a 10 + a 11 ) R + ( b 1 + b 2 + b 3 + b 4 + b 5 + b 6 + b 7 + b 8 + b 9 + b 10 + b 11 ] = C .

C being a positive constant. Therefore, | | B ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) | | C.

This establishes the uniform bounded property of the operator B.

Now, we will try to show the equicontinuity of B. Let us take two points t 1 , t 2 [ 0 , T ] with t 1 < t 2. Then,
| B 1 ( t 2 , S U ) B 1 ( t 1 , S U ) | = q G ( q ) Γ ( q ) | 0 t 2 ( t z ) q 1 Q 1 ( z , S U ( z ) ) d z 0 t 1 ( t z ) q 1 Q 1 ( z , S U ( z ) ) d z | , | B 1 ( t 2 , S U ) B 1 ( t 1 , S U ) q G ( q ) Γ ( q ) ( 0 t 2 ( t z ) q 1 d z 0 t 1 ( t z ) q 1 d z ) [ a 1 R + b 1 ] ,
which will eventually give
| B 1 ( t 2 , S U ) B 1 ( t 1 , S U ) a 1 R + b 1 G ( q ) Γ ( q ) ( t 2 q t 1 q ) .
Similarly, we get
| B 2 ( t 2 , S A ) B 2 ( t 1 , S A ) a 2 R + b 2 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 3 ( t 2 , E H ) B 3 ( t 1 , E H ) a 3 R + b 3 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 4 ( t 2 , I A ) B 4 ( t 1 , I A ) a 4 R + b 4 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 5 ( t 2 , I H ) B 5 ( t 1 , I H ) a 5 R + b 5 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 6 ( t 2 , H H ) B 6 ( t 1 , H H ) a 6 R + b 6 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 7 ( t 2 , R H ) B 7 ( t 1 , R H ) a 7 R + b 7 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 8 ( t 2 , V U ) B 8 ( t 1 , V U ) a 8 R + b 8 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 9 ( t 2 , V A ) B 9 ( t 1 , V A ) a 9 R + b 9 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 10 ( t 2 , S V ) B 10 ( t 1 , S V ) a 10 R + b 10 G ( q ) Γ ( q ) ( t 2 q t 1 q ) , | B 11 ( t 2 , I V ) B 11 ( t 1 , I V ) a 11 R + b 11 G ( q ) Γ ( q ) ( t 2 q t 1 q ) .

As t 1 t 2, we get | | B ( B 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , B 9 , B 10 , B 11 ) ( t 2 ) B ( B 1 , B 2 , B 3 , B 4 , B 5 , B 6 , B 7 , B 8 , B 9 , B 10 , B 11 ) | | ( t 1 ) | | 0. Therefore, B does not depend on ( S H , E H , I A , I H , T H , P H , R H , S R , I R , S V , I V ) . This establishes the equicontinuity of the operator B.

In order to establish the compactness and continuity parts, we employ the following theorem:

Arzela Ascoli Theorem : Let us consider a region Ω be in , and let F be a family of complex-valued functions on Ω that is pointwise bounded. Then every sequence { f n } in F has a subsequence that converges to a continuous function on Ω, the convergence being uniform on compact subsets. Since B is completely continuous, according to the Arzela–Ascoli theorem, it is relatively compact. This completes the result.

Here, we shall perform a rigorous analysis of the model (2) to find the disease-free and endemic equilibrium points and further perform their stability analysis.

To find the disease-free equilibrium (DFE) point, we shall equate the right-hand sides of the equations of the model (3) to 0. The infected variables are also equated to zero to compute this equilibrium point. On doing this, the DFE point that we get for the system of fraction differential equations (3) is given by the following manner.

By equating the right-hand sides of all the eight differential equations in the model (1) as well as all the infected variables to zero, we shall obtain the DFE point as follows:
E 0 = ( S U 0 , S A 0 , E H 0 , I A 0 , I H 0 , H H 0 , R H 0 , V U 0 , V A 0 , S V 0 , I V 0 ) , E 0 = ( λ π H μ h + u 1 σ 1 ( 1 θ 3 ) , ( 1 λ ) π H μ h + u 1 σ 1 θ 3 , 0 , 0 , 0 , 0 , 0 , u 1 σ 1 ( 1 θ 3 ) λ π H μ h ( μ h + u 1 σ 1 ( 1 θ 3 ) ) , u 1 σ 1 θ 3 ( 1 λ ) π H μ h ( μ h + u 1 σ 1 θ 3 ) , π V μ v , 0 ) .

1. Basic reproduction number

A total number of secondary infections are caused by just one infected individual among the susceptible class, while the course of the infection is defined as the basic reproduction number (Carr, 1981; Castillo-Chavez and Song, 2004; Castillo-Chavez 2002; Diekman 1990). It is considered a very significant tool in determining whether the disease with persist or fade out eventually in the course of time depending on the value of the basic reproduction number (Diekman 1990; Van den Driessche and Watmough, 2002). For our current fractional-order model, the infected compartments are E H , I A , I H , H H, and IV. We shall denote the right-hand side of the infected classes of the model (3) using two matrices F and V as defined below such that
ABC 0 D t q ( y ) = F ( y ) V ( y ) .
For our fractional-order model (4),
F = ( β I V S U N V + β 1 I H S U N H + β 2 I H S A N H 0 0 0 β 3 S V I V N V ) , V = ( ( η 1 + μ h ) E H ( 1 θ 1 ) η 1 E H + ( ρ 1 I A + μ h ) I A thet a 1 η 1 E H + ( α 1 + δ + μ h ) I H ( 1 η 2 ) α 1 I H θ 2 ρ 1 I A + ( μ h + r ) H H ( μ v + δ 1 ) I V ) .
Now, we shall evaluate the partial derivatives of the matrices F and V obtained above with respect to the infected classes at the DFE E 0 to obtain the following matrices F and V:
F = ( 0 0 β 1 + β 2 0 β 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 β 3 )
and
V = ( η 1 + μ h 0 0 0 0 ( 1 θ 1 ) η 1 ρ 1 + μ h 0 0 0 θ 1 η 1 0 α 1 + δ + μ h 0 0 0 θ 2 ρ 1 ( 1 η 2 ) α 1 μ h + r 0 0 0 0 0 μ v + δ 1 ) .
After some rigorous calculations, the basic reproduction number is obtained by the largest positive eigenvalue of F V 1, which is
R 0 = max { R 0 V , R 0 H } = max { β 3 μ v + δ 1 , η 1 θ 1 ( β 1 + β 2 ) ( η 1 + μ h ) ( α 1 + δ + μ h ) } .

2. Effective reproduction number

To assess the transmission dynamics of the disease under investigation in a timely and scientifically rigorous manner, we employ the concept of the effective reproduction number, denoted as Rt. Similar to the basic reproduction number, R 0, we estimate Rt by multiplying the fraction of the susceptible host population by R 0. Estimating Rt effectively necessitates consideration of several key parameters, including the number of infected cases, the chosen serial time interval, and the time of symptom onset in infected individuals. Rt is a crucial metric for understanding the current state of the epidemic. When Rt equals 1, it signifies that the disease is in an endemic phase, where the number of infected cases remains stable over time. Conversely, when Rt is less than or equal to 1, it indicates that the number of infected cases is declining, eventually leading to the eradication of the disease from the population under specific conditions. Rt essentially represents the average number of secondary cases generated by a single primary case within a population at a given point in time, which allows us to assess the potential for disease transmission. It is important to note that Rt is dynamic and may change as the immunity of the population changes over time due to various factors, including vaccination and natural infection. In our analysis, we make the assumption that the distribution of the serial interval follows a gamma distribution with a mean of 8.5, which is a common approach in epidemiological modeling. In Fig. 1, we have calculated the daily values of the effective reproduction number, Rt, using the moving average method. This plot illustrates the changing dynamics of disease transmission over time. Our model, which incorporates specific control measures, suggests that the infection rate declines to below 1 within a few months. This decline in Rt indicates a reduction in the spread of the disease, indicating that, under the applied measures, it is possible to control and manage the epidemic over a defined time period.

FIG. 1.

Plot of the effective reproduction number Rt with respect to time using the moving average method.

FIG. 1.

Plot of the effective reproduction number Rt with respect to time using the moving average method.

Close modal

Theorem. The DFE point E0 is locally asymptotically stable if the reproduction number R 0 < 1 and is unstable if R 0 > 1 (Ghosh 2016; Pal 2015).

Proof. For the system (2), the Jacobian matrix J computed at the DFE point E0 is given as follows:
( μ h u 1 σ 1 ( 1 θ 3 ) 0 0 0 0 0 0 0 0 0 0 0 μ h u 1 σ 1 θ 3 0 0 0 0 0 0 0 0 0 0 0 η 1 μ h 0 0 0 0 0 0 0 0 0 0 ( 1 θ 1 ) η 1 ρ 1 μ h 0 0 0 0 0 0 0 0 0 θ 1 η 1 0 α 1 δ μ h 0 0 0 0 0 0 0 0 0 θ 2 ρ 1 ( 1 η 2 ) α 1 μ h r 0 0 0 0 0 0 0 0 ( 1 θ 2 ) ρ 1 η 2 α 1 r 0 0 μ h 0 0 u 1 σ 1 ( 1 θ 3 ) 0 0 0 0 0 μ h 0 0 0 0 0 u 1 σ 1 θ 3 0 0 0 0 0 μ h 0 0 0 0 0 0 0 0 0 0 0 0 μ v 0 0 0 0 0 0 0 0 0 0 0 μ v δ 1 ) .
The eigenvalues of the Jacobian matrix are
Λ 1 = μ h u 1 σ 1 ( 1 θ 3 ) , Λ 2 = μ h u 1 σ 1 θ 3 , Λ 3 = α 1 δ μ h , Λ 4 = η 1 μ h , Λ 5 = ρ 1 μ h , Λ 6 = μ h r , Λ 7 = μ v δ 1 , Λ 8 = μ h 2 ( 3 μ h ( 1 e t a 2 ) α 1 ) 2 , Λ 9 = μ h 2 + ( 3 μ h ( 1 e t a 2 ) α 1 ) 2 , Λ 10 = μ h , Λ 11 = μ v .

We can observe that the eigenvalues of the Jacobian matrix at the disease-free equilibrium point are all negative. The eigenvalues being all negative indicate that the disease-free equilibrium point E 0 is locally asymptotically stable whenever R 0 < 1. From this, we can conclude that the transmission of the monkeypox disease can be contained within the given population under the condition R0 being less than unity and the initial population size taken for the fractional-order system lies in the basin of attraction of the equilibrium point E 0 The proof is thus established.

Theorem. An attractive solution will exist for the fractional-order model (2) provided that the zero solution Φ ( t ) = 0 will satisfy the equation lim t y 0 = 0 for | | y 0 | | < ϵ.

Proof. We have previously shown that pure model (2) possesses a solution, which is unique. We can now define the zero solution as follows:
( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) ( t ) = 0 .

Thus, | | ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) | | ϵ lim t ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) ( t ) = 0.

This shows that the solution will be attractive in nature and therefore asymptotically stable.

The investigation into the stability of functional equations can be traced back to a seminal inquiry posed by Stanislaw Ulam ca. 1940, with the primary focus on the stability of homomorphisms within groups. Approximately 1 year later, D. H. Hyer embarked on the exploration of this inquiry by leveraging the concept of additive mappings between Banach spaces. He harnessed the “contraction mapping theorem” to establish foundational results, marking the inception of significant advancements in the examination of the problem initially proposed by Ulam. Subsequently, in 1978, another distinguished researcher named Rassias undertook this question and extended and generalized the proofs and outcomes originally put forth by Hyer. The concept of the Hyer–Ulam stability has seen extensive research application in the realms of both ordinary and fractional differential equations. In the context of mathematical modeling for diseases, the analysis of stability assumes paramount importance. Moreover, the Hyer–Ulam stability framework is highly regarded as an indispensable tool for scrutinizing the dynamics of nonlinear fractional models. This form of stability analysis proves exceptionally valuable when addressing a wide spectrum of practical problems for which obtaining an exact solution is a formidable challenge. In such scenarios, the Hyer–Ulam stability framework facilitates the approximation of solutions, enabling the dynamic analysis of fractional-order systems.

In order to determine the stability of the solution of the model (2), we consider a particular kind of differential inequality described below.

Theorem. Let us consider the matrix H as follows:
( h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 ) .
Now, for the fractional-order system, the unique solution that exists will be Hyer–Ulam's stable if the spectral radius of H lies strictly within the unit circle given by | h 1 + h 2 + h 3 + h 4 + h 5 + h 6 + h 7 + h 8 + h 9 + h 10 + h 11 | < 1, where h 1 , h 2 , h 3 , h 4 , h 5 , h 6 , h 7 , h 8 , h 9 , h 10 , h 11 are defined as follows:
h 1 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 1 , h 2 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 2 , h 3 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 3 , h 4 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 4 , h 5 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 5 , h 6 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 6 , h 7 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 7 , h 8 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 8 , h 9 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 9 , h 10 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 10 , h 11 = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 11 .

Proof. We have already shown that the fractional-order system (3) has a unique solution. Now, let us consider two different solutions of the fractional-order model (2) of the form:

( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) and ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ). By employing (Sec. IV C), we get
| | F 1 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 1 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | = 1 q G ( q ) | | Q 1 ( t , S U ( t ) ) Q 1 ( t , S ̃ U ( t ) ) | | + q G ( q ) Γ ( q ) 0 t | | Q 1 ( z , S U ( z ) ) Q 1 ( z , S ̃ U ( z ) ) | | ( t z ) q 1 d z .
On further simplifying, we get
| | F 1 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 1 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | 1 q G ( q ) c 1 | | S U S ̃ U | | + x q G ( q ) Γ ( q ) c 1 | | S U S ̃ U | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 1 | | S U S ̃ U | | .
We can rewrite the above as
| | F 1 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 1 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 1 | | S U S ̃ U | | .
Similarly,
| | F 2 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 2 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 2 | | S A S ̃ A | | , | | F 3 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 3 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 3 | | E H E ̃ H | | , | | F 4 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 4 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 4 | | I A I ̃ A | | , | | F 5 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 5 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 5 | | I H I ̃ H | | , | | F 6 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 6 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 6 | | H H H ̃ H | | , | | F 7 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 7 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 7 | | R H R ̃ H | | , | | F 8 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 8 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 8 | | V U V ̃ U | | , | | F 9 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 9 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 9 | | V A V ̃ A | | , | | F 10 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 10 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 10 | | S V S ̃ V | | , | | F 11 ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F 11 ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) c 11 | | I V I ̃ V | | .
Thus, we can conclude that | | F ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) C [ | | ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ] , where C denotes the maximum value of the Lipschitz constants, i.e., C = max { c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , c 7 , c 8 , c 9 , c 10 , c 11 }. So,
| | F ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) F ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | = ( 1 q G ( q ) + x q G ( q ) Γ ( q ) ) [ c 1 | | S U S ̃ U | | + c 2 | | S A S ̃ A | | + c 3 | | E H E ̃ H | | + c 4 | | I A I ̃ A | | + c 5 | | I H I ̃ H | | + c 6 | | H H H ̃ H | | + c 7 | | R H R ̃ H | | + c 8 | | V U V ̃ U | | + c 9 | | V A V ̃ A | | + c 10 | | S V S ̃ V | | + c 11 | | I V I ̃ V | | ] .
Linearizing the above norms in matrix form for ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) and ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ), we shall obtain
| | ( S U , S A , E H , I A , I H , H H , R H , V U , V A , S V , I V ) ( S ̃ U , S ̃ A , E ̃ H , I ̃ A , I ̃ H , H ̃ H , H ̃ H , V ̃ U , V ̃ A , S ̃ V , I ̃ V ) | | ( h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 1 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 2 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 3 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 4 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 5 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 6 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 7 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 8 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 9 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 10 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 h 11 ) , ( | | S U S ̃ U | | | | S A S ̃ A | | | | E H E ̃ H | | | | I A I ̃ A | | | | I H I ̃ H | | | | H H H ̃ H | | | | R H R ̃ H | | | | V U V ̃ U | | | | V A V ̃ A | | | | S V S ̃ V | | | | I V I ̃ V | | ) .

Keeping in mind, the definition of h1, h2, h3, h4, h5, h6, h7, h8, h9, h10, and h11, we can observe that the above expression will converge to zero. For the above matrix, the eigenvalues are as follows: e 1 = 0 , e 2 = 0 , e 3 = 0 , e 4 = 0 , e 5 = 0 , e 6 = 0 , e 7 = 0 , e 8 = 0 , e 9 = 0 , e 10 = 0 and e 11 = h 1 + h 2 + h 3 + h 4 + h 5 + h 6 + h 7 + h 8 + h 9 + h 10 + h 11. We get the spectral radius S = max { | e s | , s = 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 } = h 1 + h 2 + h 3 + h 4 + h 5 + h 6 + h 7 + h 8 + h 9 + h 10 + h 11 < 1. Therefore, we can say that the fractional-order monkeypox model is Hyer–Ulam's stable as a result of the discrete dichotomy shown above.

Maximum likelihood estimation (MLE) stands as a pivotal and widely endorsed statistical tool employed for parameter estimation, as documented by Myung (2003), Portet (2020), and Zeng and Lin (2007). Its robustness is particularly evident when dealing with nonlinear models and data sets that do not conform to the normal distribution. To employ MLE, one must have a likelihood function, which serves as the foundation for making statistical inferences. In practice, when armed with a set of observed data and a mathematical model of interest, the challenge lies in deriving a suitable probability density function that accurately represents the data generation process. This endeavor can be viewed as an inverse problem and can be effectively addressed through the formulation of a likelihood function. Detailed information on the calibration of the model using the maximum likelihood estimation can be found in the  Appendix, accompanied by a visual representation of the data set. In our specific application, we employ MLE to calibrate model (1) to the annual incidence data of monkeypox cases in the United States. The parameters under estimation include β, β1, β2, β3, η1, and η2. Figure 2 illustrates the curve representing the model fitting achieved through the maximum likelihood estimation technique. This fitting process is employed to estimate the aforementioned parameters based on the recently recorded cases of monkeypox infection in the United States. Furthermore, Table I provides a comprehensive listing of the parameter values employed within our monkeypox model. This meticulous calibration using MLE enables us to refine the model's parameters and align it more closely with the empirical data, thus enhancing our understanding of the dynamics of monkeypox transmission in the United States (Table II).

FIG. 2.

Fitting of the monkeypox model (1) with the cumulative data of the United States.

FIG. 2.

Fitting of the monkeypox model (1) with the cumulative data of the United States.

Close modal
FIG. 3.

Sensitivity of R 0 with respect to the parameters β, β1, α1, and η1 using 10 000 samples (significant parameters have p-value less than 0.05).

FIG. 3.

Sensitivity of R 0 with respect to the parameters β, β1, α1, and η1 using 10 000 samples (significant parameters have p-value less than 0.05).

Close modal
TABLE I.

Description of parameters in model (1) and their values used for numerical simulations. The parameters estimated are with a confidence level of 95%. The total size of the human and rodent population taken initially are given by N H ( 0 ) and N V ( 0 ), respectively.

Parameters Descriptions Values Sources
λ  Fraction of susceptible humans that are aware  0.25  Assumed 
πH  Constant rate of recruitment for humans  πH* N H ( 0 ) day 1   
μh  Human mortality rate  0.0105 day 1  The World Bank (2012)  
πV  Constant rate of recruitment for vector  0.2 day 1  Peter (2022)  
μv  Vector mortality rate  0.002 day 1  Peter (2022)  
δ  Death rate of humans due to monkeypox  0.003 25 day−1  Assumed 
δ1  Death rate of rodents due to monkeypox  0.0064 day−1  Assumed 
α1  Movement rate from symptomatic infected human population  0.08 day−1  Assumed 
β  Effective contact rate of rodents to unaware humans  0.000 24 day−1  Estimated 
β1  Effective contact rate of unaware humans to other humans  0.3380 day−1  Estimated 
β2  Effective contact rate of aware humans to other humans  0.000 002 3 day−1  Estimated 
β3  Effective contact rate of rodents to rodents  0.1073 day−1  Estimated 
ρ1  Movement rate from asymptomatically infected humans  0.001  Assumed 
η1  Movement rate from exposed humans  0.92 day−1  Estimated 
η2  Fraction of symptomatic humans that are recovered  0.43 day−1  Estimated 
θ1  Fraction of exposed humans that become symptomatic  0.99  Assumed 
θ2  Fraction of asymptomatically infected humans that are hospitalized  0.5  Assumed 
θ3  Fraction of aware humans that are vaccinated  0.05  Assumed 
u1  Vaccination rate of humans  0–1  Assumed 
σ1  Efficiency rate of vaccine  0–1  Assumed 
r  Recovery rate of hospitalized individuals  0.036  Assumed 
Parameters Descriptions Values Sources
λ  Fraction of susceptible humans that are aware  0.25  Assumed 
πH  Constant rate of recruitment for humans  πH* N H ( 0 ) day 1   
μh  Human mortality rate  0.0105 day 1  The World Bank (2012)  
πV  Constant rate of recruitment for vector  0.2 day 1  Peter (2022)  
μv  Vector mortality rate  0.002 day 1  Peter (2022)  
δ  Death rate of humans due to monkeypox  0.003 25 day−1  Assumed 
δ1  Death rate of rodents due to monkeypox  0.0064 day−1  Assumed 
α1  Movement rate from symptomatic infected human population  0.08 day−1  Assumed 
β  Effective contact rate of rodents to unaware humans  0.000 24 day−1  Estimated 
β1  Effective contact rate of unaware humans to other humans  0.3380 day−1  Estimated 
β2  Effective contact rate of aware humans to other humans  0.000 002 3 day−1  Estimated 
β3  Effective contact rate of rodents to rodents  0.1073 day−1  Estimated 
ρ1  Movement rate from asymptomatically infected humans  0.001  Assumed 
η1  Movement rate from exposed humans  0.92 day−1  Estimated 
η2  Fraction of symptomatic humans that are recovered  0.43 day−1  Estimated 
θ1  Fraction of exposed humans that become symptomatic  0.99  Assumed 
θ2  Fraction of asymptomatically infected humans that are hospitalized  0.5  Assumed 
θ3  Fraction of aware humans that are vaccinated  0.05  Assumed 
u1  Vaccination rate of humans  0–1  Assumed 
σ1  Efficiency rate of vaccine  0–1  Assumed 
r  Recovery rate of hospitalized individuals  0.036  Assumed 
TABLE II.

Initial population sizes in model (1) estimated with a confidence interval of 95%.

Variables Descriptions Values Sources
N H ( 0 )  Initial total human population  332 403 650  U.S. Department of Commerce (2023)  
S U ( 0 )  Initial susceptible unaware human population  2.29 × 10 7  Estimated 
S A ( 0 )  Initial susceptible aware human population  N H ( 0 ) S U ( 0 ) E H ( 0 ) I A ( 0 ) I H ( 0 ) V U ( 0 ) V A ( 0 )   
E H ( 0 )  Initial population of exposed humans  206  Estimated 
I A ( 0 )  Initial population of humans infected with asymptomatic monkeypox  Assumed 
I H ( 0 )  Initial population of humans infected with symptomatic monkeypox  Assumed 
H H ( 0 )  Initial population of hospitalized humans  Assumed 
R H ( 0 )  Initial population of recovered humans  Assumed 
N V ( 0 )  Initial total rodent population  S V ( 0 ) + I V ( 0 )   
S V ( 0 )  Initial susceptible rodent population  1.82 × 10 7  Estimated 
I V ( 0 )  Initial population of rodent that are infected with monkeypox  2.3041 × 10 4  Estimated 
Variables Descriptions Values Sources
N H ( 0 )  Initial total human population  332 403 650  U.S. Department of Commerce (2023)  
S U ( 0 )  Initial susceptible unaware human population  2.29 × 10 7  Estimated 
S A ( 0 )  Initial susceptible aware human population  N H ( 0 ) S U ( 0 ) E H ( 0 ) I A ( 0 ) I H ( 0 ) V U ( 0 ) V A ( 0 )   
E H ( 0 )  Initial population of exposed humans  206  Estimated 
I A ( 0 )  Initial population of humans infected with asymptomatic monkeypox  Assumed 
I H ( 0 )  Initial population of humans infected with symptomatic monkeypox  Assumed 
H H ( 0 )  Initial population of hospitalized humans  Assumed 
R H ( 0 )  Initial population of recovered humans  Assumed 
N V ( 0 )  Initial total rodent population  S V ( 0 ) + I V ( 0 )   
S V ( 0 )  Initial susceptible rodent population  1.82 × 10 7  Estimated 
I V ( 0 )  Initial population of rodent that are infected with monkeypox  2.3041 × 10 4  Estimated 

In this section, we conduct a sensitivity analysis for the given fractional-order model (2). The primary objective of this analysis is to ascertain the key parameters that exert the most significant influence on the control of disease transmission (Marino 2008). This approach is based on the works of Peter (2023), Baithalu (2023), and Zarin (2022) and combines two well-established techniques: Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) multivariate analysis. To perform this sensitivity analysis, we randomly sample 1000 parameter sets. Our focus is on identifying the parameters that have a substantial impact on the basic reproduction number, R 0. From the analysis of Fig. 3, it becomes evident that the parameters β, β1, α1, and η1 exhibit sensitivity concerning the basic reproduction number, R 0. Notably, parameters β, β1, and η1 demonstrate positive sensitivity indices, implying that increases in these parameters result in higher values of R 0. This outcome is counterproductive to our goal of disease transmission control, as it exacerbates the spread of the infection. Thus, elevating these parameters would not be an effective strategy for containment. Conversely, a negative sensitivity index is observed for parameter α1. This implies that as α1 increases, signifying improved treatment for symptomatic infected individuals, the value of R 0 decreases. This inverse correlation with α1 is a promising finding. In practical terms, enhancing α1 by providing more effective treatment to symptomatic individuals, leading to their transition into the recovered or hospitalized classes, results in a decline in the value of R 0. This reduction in the basic reproduction number is indicative of the potential for disease eradication within the environment. In summary, increasing α1 emerges as a scientifically sound and effective control strategy for curbing the spread of the monkeypox infection. Identifying these significant parameters provides valuable insights into the factors that have the most substantial influence on the dynamics of disease transmission within the population. This information is instrumental in designing effective control strategies and interventions to mitigate the spread of the disease.

FIG. 4.

Time series solutions of the monkeypox model (1) exhibiting the impacts of effective contact rate β on the spread of infection. (a) Simulation results for IA with respect to β. (b) Simulation results for IH with respect to β. (c) Simulation results for HH with respect to β.

FIG. 4.

Time series solutions of the monkeypox model (1) exhibiting the impacts of effective contact rate β on the spread of infection. (a) Simulation results for IA with respect to β. (b) Simulation results for IH with respect to β. (c) Simulation results for HH with respect to β.

Close modal

In this section, we delve into the analysis of various control strategies designed to mitigate the propagation of the monkeypox infection. Drawing upon the findings from our numerical simulations, we aim to identify the most efficacious single or combined control strategy that could contribute to the eradication of the infection while maintaining environmental safety. For our numerical experiments, we employed the Adams–Bashforth predictor–corrector method. The computational framework was adapted from the approach outlined in (Biswas , 2022; Kumar, 2022; Khan , 2019). This method allowed us to approximate solutions for the fractional-order system, specifically within the context of the Atangana–Baleanu derivative, providing us with a precise and detailed analysis of the dynamics of the monkeypox infection.

In Fig. 4, we investigate the ramifications of varying the contact rate parameter, β, between the rodent population and unaware humans, on the dynamics of three distinct human populations in the context of the monkeypox infection: asymptomatic infected individuals (IA), symptomatic infected individuals (IH), and hospitalized individuals (HH). At β = 0.024, we observe a substantial escalation in the number of infected and hospitalized humans approximately 30 days into the simulation. The maximum values for IA, IH, and HH reach 3500, 125 000, and 75 000, respectively. Upon decreasing the value of β, significant changes become evident in the number of cases across all three compartments. When β = 0.000 24, the number of infected and hospitalized cases approaches zero, indicating a substantial improvement compared to the previous scenario. Therefore, diminishing the contact rate between the rodent population and the unaware human population emerges as an efficacious strategy for disease transmission control.

FIG. 5.

Time series solutions of the monkeypox model (1) exhibiting the impacts of effective contact rate β1 on the spread of infection. (a) Simulation results for IA with respect to β1. (b) Simulation results for IH with respect to β1. (c) Simulation results for HH with respect to β1.

FIG. 5.

Time series solutions of the monkeypox model (1) exhibiting the impacts of effective contact rate β1 on the spread of infection. (a) Simulation results for IA with respect to β1. (b) Simulation results for IH with respect to β1. (c) Simulation results for HH with respect to β1.

Close modal

In Fig. 5, we analyze the impact of varying the contact rate parameter β1 on the dynamics of three distinct human populations concerning the monkeypox infection: asymptomatic infected individuals (IA), symptomatic infected individuals (IH), and hospitalized individuals (HH). When β 1 = 0.338, there is a noticeable upsurge in the number of infected and hospitalized individuals, with the peaks occurring roughly 20 days into the simulation. The maximum values for IA, IH, and HH reach approximately 35, 1270, and 760, respectively, after around 100 days. As β1 is progressively reduced, a marginal decline in the number of infected and hospitalized individuals is observed. For instance, at β 1 = 0.000 338, the reduction is modest. After 100 days, the maximum number of asymptomatic infected (IA), symptomatic infected (IH), and hospitalized (HH) individuals becomes approximately 31, 1100, and 680. Thus, reducing the contact rate between unaware individuals and other humans does result in a decrease in the number of infected cases, but the reduction is not substantial. Consequently, this strategy may not be deemed one of the most effective approaches for controlling the spread of monkeypox infection.

FIG. 6.

Time series solutions of the monkeypox model (1) exhibiting the impacts of vaccination (u1) on the spread of infection. (a) Simulation results for IA with respect to u1. (b) Simulation results for IH with respect to u1. (c) Simulation results for HH with respect to u1.

FIG. 6.

Time series solutions of the monkeypox model (1) exhibiting the impacts of vaccination (u1) on the spread of infection. (a) Simulation results for IA with respect to u1. (b) Simulation results for IH with respect to u1. (c) Simulation results for HH with respect to u1.

Close modal

Figure 6 illustrates the influence of vaccination parameter u1 on the dynamics of three distinct human populations within the context of monkeypox infection: asymptomatic infected individuals (IA), symptomatic infected individuals (IH), and hospitalized individuals (HH). When u 1 = 0, signifying the absence of vaccination, the model reveals a substantial upsurge in the number of infected and hospitalized humans approximately 30 days into the simulation. After 100 days, these populations reach 184, 7500, and 4000 for IA, IH, and HH, respectively. The introduction of vaccination into the population results in a notable alteration in disease transmission dynamics. As the vaccination rate increases, a discernible decline in the prevalence of cases within each of the three compartments is observed. At u 1 = 0.75, there is a subtle but significant reduction in the number of infected and hospitalized cases when compared to the scenario with no vaccination. After 100 days, the maximum number of cases in the IA, IH, and HH compartments reaches 2, 200, and 100, respectively. Hence, providing vaccination to both infected individuals and those requiring hospitalization can be considered an efficacious strategy for mitigating the spread of the monkeypox infection.

FIG. 7.

Time series solutions of the monkeypox model (1) exhibiting the impacts of variations in the fractional order q on the spread of infection. (a) Simulation results for IA with respect to q. (b) Simulation results for IH with respect to q. (c) Simulation results for HH with respect to q.

FIG. 7.

Time series solutions of the monkeypox model (1) exhibiting the impacts of variations in the fractional order q on the spread of infection. (a) Simulation results for IA with respect to q. (b) Simulation results for IH with respect to q. (c) Simulation results for HH with respect to q.

Close modal

In Fig. 7, we have conducted an investigation into the impact of varying the fractional-order parameter, q on the dynamics of monkeypox infection. Specifically, we have analyzed the changes in the numbers of asymptomatically infected individuals (IA), symptomatically infected individuals (IH), and hospitalized cases (HH). Our findings suggest that the choice of the fractional-order parameter q plays a crucial role in influencing the progression of the monkeypox infection. When q is set to 0.99, we observe a noteworthy upsurge in the number of infected and hospitalized cases approximately 30 days into the simulation. After 100 days, the recorded values show a substantial increase, reaching around 35 for IA, 1270 for IH, and 760 for HH. This rapid escalation in cases indicates the highly contagious nature of monkeypox in this fractional order. Conversely, as we gradually decrease the fractional order q, we notice a concomitant decrease in the number of infected and hospitalized cases. Notably, at q = 0.8, a slight surge in cases is observed after approximately 50 days, though this surge is relatively modest. Under this scenario, the peak numbers in the compartments for asymptomatically infected, symptomatically infected, and hospitalized cases are limited to 3, 180, and 50, respectively. These results signify a notable improvement in controlling the transmission of the infection when compared to the higher q value. This underscores the significance of reducing the fractional-order parameter (q) as a highly recommended strategy for mitigating the spread of the monkeypox disease. This strategy appears effective in reducing the overall number of infected and hospitalized cases, indicating its potential as a valuable intervention in the management and containment of this infectious disease.

FIG. 8.

Graphical representation of the cumulative monkeypox cases in the United States.

FIG. 8.

Graphical representation of the cumulative monkeypox cases in the United States.

Close modal

In the current investigation, we have employed an innovative compartmental model to examine the dynamics of monkeypox disease. Within this model, we have subdivided the human population into two distinct compartments, distinguishing between those who are informed about the infection (the “aware” population) and those who remain uninformed (the “unaware” population). This approach has allowed us to conduct a comprehensive analysis of the transmission dynamics. The primary focus of our study encompasses the human and rodent populations, as they are the key categories under scrutiny. Our model was rigorously parameterized using empirical data on monkeypox-infected cases in the United States. The graphical representations of our model's outputs reveal a notable alignment with the collected empirical data, affirming the reliability of our model for the estimation of various parameters. These parameters, in turn, facilitate enhanced predictions concerning the dynamics of the monkeypox infection. In the course of our research, we have computed both the basic and effective reproduction numbers and established conditions for the stability of the disease-free equilibrium point. Furthermore, our sensitivity analysis of the basic reproduction number has yielded insightful conclusions regarding the impact of various factors on the transmission of the infection within the population. Additionally, we have demonstrated the Hyer–Ulam stability for our model. After a thorough examination of the asymptotic behavior of the disease-free equilibrium point, we have reached the critical conclusion that the complete eradication of monkeypox from the system is only achievable when the basic reproduction number, R 0, remains strictly below one. This condition holds, assuming that the initial population size of infected humans and rodents falls within the basin of attraction of the disease-free equilibrium point.

Our research outcomes have elucidated that the implementation of diverse control strategies, such as the dissemination of infection awareness, widespread vaccination campaigns, and the provision of efficacious treatments for infected individuals, holds paramount significance in mitigating the monkeypox epidemic and curtailing its propagation within the population. Employing numerical simulations, we have discerned the critical role played by minimizing the effective contact rate between individuals who are informed about the infection and those who remain unaware, as well as between the human and rodent populations, in effectively modulating disease transmission dynamics. These findings have the potential to offer valuable insights to governmental and community policymakers. They may use this knowledge to refine existing policies and develop novel, evidence-based approaches for enhancing the control of monkeypox spread, ultimately safeguarding public health and well-being.

The authors express their gratitude to the associate editor and the learned reviewers whose comments and suggestions have helped to improve this paper. Research of Santanu Biswas is supported by Dr. D. S. Kothari Postdoctoral Fellowship under the University Grants Commission scheme [Ref. No. F.4-2/2006 (BSR)/MA/19-20/0057].

The authors have no conflicts to disclose.

A. Santanu Biswas and B. Humaira Aslam have contributed equally to this work.

Santanu Biswas: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Resources (lead); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Humaira Aslam: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (lead); Resources (equal); Writing – original draft (lead); Writing – review & editing (equal). Pankaj Kumar Tiwari: Conceptualization (supporting); Supervision (supporting); Writing – review & editing (supporting).

The data are taken from https://www.commerce.gov and https://data.worldbank.org. Some parameter values have been taken from the work of Peter (2022), and the remaining ones are either assumed or estimated.

The data sets generated during and/or analyzed during the current research work are available from the corresponding author upon reasonable request.

Suppose for any data vector u = ( u 1 , u 2 , u 3 , . , u m ) and a parameter vector v = ( v 1 , v 2 , v 3 . , v m ), let f(u|v) be a probability density function such that with respect to this the likelihood function can be defined as L ( v | u ) = f ( u | v ), which represents the likelihood function of the parameter variable v, providing that the recorded data u are already given. After collecting the necessary data and determining the likelihood function, the next step would be to reach certain statistical conclusions about the population and the study. Since different values of the parameter variable correspond to different probability distributions, we are just interested in finding or estimating those parameter values that are indexed by the probability distribution underlying the given data set. The home likelihood estimation technique is based on the principle of estimating the parameter vector value that will maximize the obtain likelihood function. On exploring the multidimensional parameter space, the parameter vectors obtained as a result are known as the maximum likelihood estimate and are expressed as v MLE = ( v 1 , MLE , v 2 , MLE , , v k , MLE ). So, basically, MLE is just a method of finding a suitable probability distribution that will make the observed set of data most probable See Fig. 8.

1.
Abdeljawad
,
T.
,
Hajji
,
M. A.
,
Al-Mdallal
,
Q. M.
, and
Jarad
,
F.
, “
Analysis of some generalized ABC-fractional logistic models
,”
Alexandria Eng. J.
59
(
4
),
2141
2148
(
2020
).
7.
Abro
,
K. A.
, “
Fractional characterization of fluid and synergistic effects of free convective flow in circular pipe through Hankel transform
,”
Phys. Fluids
32
(
12
),
123102
(
2020
).
2.
Abro
,
K. A.
,
Atangana
,
A.
, and
Gómez-Aguilar
,
J. F.
, “
A comparative analysis of plasma dilution based on fractional integro-differential equation: An application to biological science
,”
Int. J. Model. Simul.
43
(1),
1
10
(
2023
).
2.
Adel
,
W.
,
Elsonbaty
,
A.
,
Aldurayhim
,
A.
, and
El-Mesady
,
A.
, “
Investigating the dynamics of a novel fractional-order monkeypox epidemic model with optimal control
,”
Alexandria Eng. J.
73
,
519
542
(
2023
).
3.
Adom-Konadu
,
A.
,
Bonyah
,
E.
,
Sackitey
,
A. L.
,
Anokye
,
M.
, and
Asamoah
,
J. K. K.
, “
A fractional order monkeypox model with protected travelers using the fixed point theorem and Newton polynomial interpolation
,”
Healthcare Anal.
3
,
100191
(
2023
).
4.
Adom-Konadu
,
A.
,
Bonyah
,
E.
,
Sackitey
,
A. L.
,
Anokye
,
M.
, and
Asamoah
,
J. K.
, “
A fractional order monkeypox model with protected travelers using the fixed point theorem and Newton polynomial interpolation
,”
Healthcare Analytics
3
,
100191
(
2023
).
5.
Agaba
,
G.
,
Kyrychko
,
Y.
, and
Blyuss
,
K.
, “
Mathematical model for the impact of awareness on the dynamics of infectious diseases
,”
Math. Biosci.
286
,
22
30
(
2017
).
6.
Ahmed
,
I.
,
Modu
,
G. U.
,
Yusuf
,
A.
,
Kumam
,
P.
, and
Yusuf
,
I.
, “
A mathematical model of Coronavirus Disease (COVID-19) containing asymptomatic and symptomatic classes
,”
Results Phys.
21
,
103776
(
2021
).
8.
Agusto
,
F. B.
, “
Mathematical model of Ebola transmission dynamics with relapse and reinfection
,”
Math. Biosci.
283
,
48
59
(
2017
).
9.
Alharbi
,
R.
,
Jan
,
R.
,
Alyobi
,
S.
,
Altayeb
,
Y.
, and
Khan
,
Z.
, “
Mathematical modeling and stability analysis of the dynamics of monkeypox via fractional-calculus
,”
Fractals
30
(
10
),
2240266
(
2022
).
10.
Ali
,
Z.
,
Rabiei
,
F.
,
Shah
,
K.
, and
Khodadadi
,
T.
, “
Fractal-fractional order dynamical behavior of an HIV/AIDS epidemic mathematical model
,”
Eur. Phys. J. Plus
136
(
1
),
36
(
2021
).
11.
Alzubaidi
,
A. M.
,
Othman
,
H. A.
,
Ullah
,
S.
,
Ahmad
,
N.
, and
Alam
,
M. M.
, “
Analysis of Monkeypox viral infection with human to animal transmission via a fractional and Fractal-fractional operators with power law kernel
,”
Math. Biosci. Eng.
20
,
6666
6690
(
2023
).
12.
Aslam
,
M.
,
Murtaza
,
R.
,
Abdeljawad
,
T.
,
Rahman
,
G. U.
,
Khan
,
A.
,
Khan
,
H.
, and
Gulzar
,
H.
, “
A fractional order HIV/AIDS epidemic model with Mittag-Leffler kernel
,”
Adv. Differ. Equations
2021
(
1
),
15
.
13.
Atangana
,
A.
and
Koca
,
I.
, “
Chaos in a simple nonlinear system with Atangana-Baleanu derivatives with fractional order
,”
Chaos, Solitons Fractals
89
,
447
454
(
2016
).
14.
Atangana
,
A.
and
Owolabi
,
K. M.
, “
New numerical approach for fractional differential equations
,”
Math. Model. Nat. Phenom.
13
(
1
),
3
(
2018
).
15.
Atangana
,
A.
and
Qureshi
,
S.
, “
Modeling attractors of chaotic dynamical systems with fractal-fractional operators
,”
Chaos, Solitons Fractals
123
,
320
337
(
2019
).
16.
Ayinla
,
A. Y.
,
Othman
,
W. A. M.
, and
Rabiu
,
M.
, “
A mathematical model of the tuberculosis epidemic
,”
Acta Biotheor.
69
,
225
255
(
2021
).
17.
Baithalu
,
R.
,
Mishra
,
S. R.
, and
Ali Shah
,
N.
, “
Sensitivity analysis of various factors on the micropolar hybrid nanofluid flow with optimized heat transfer rate using response surface methodology: A statistical approach
,”
Phys. Fluids
35
(
10
),
102016
(
2023
).
18.
Baleanu
,
D.
,
Jajarmi
,
A.
, and
Hajipour
,
M.
, “
On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag-Leffler kernel
,”
Nonlinear Dyn.
94
,
397
414
(
2018
).
19.
Baleanu
,
D.
,
Jajarmi
,
A.
,
Mohammadi
,
H.
, and
Rezapour
,
S.
, “
A new study on the mathematical modelling of human liver with Caputo-Fabrizio fractional derivative
,”
Chaos, Solitons Fractals
134
,
109705
(
2020
).
20.
Bankuru
,
S. V.
,
Kossol
,
S.
,
Hou
,
W.
,
Mahmoudi
,
P.
,
Rychtář
,
J.
, and
Taylor
,
D.
, “
A game-theoretic model of Monkeypox to assess vaccination strategies
,”
PeerJ
8
,
e9272
(
2020
).
21.
Biswas
,
S.
, “
Mathematical modeling of visceral leishmaniasis and control strategies
,”
Chaos, Solitons Fractals
104
,
546
556
(
2017a
).
22.
Biswas
,
S.
,
Aslam
,
H.
,
Das
,
S.
, and
Ghosh
,
A.
, “
Chaotic behavior in a novel fractional order system with no equilibria
,” in
Nonlinear Dynamics and Applications: Proceedings of the ICNDA
(
Springer
,
2022
), pp.
1081
1092
.
23.
Cai
,
L.
,
Li
,
X.
,
Tuncer
,
N.
,
Martcheva
,
M.
, and
Lashari
,
A.
, “
Optimal control of a malaria model with asymptomatic class and superinfection
,”
Math. Biosci.
288
,
94
(
2017
).
24.
Carr
,
J.
,
Applications of Centre Manifold Theory
(
Springer-Verlag
,
New York
,
1981
).
25.
Castillo-Chavez
,
C.
and
Song
,
B.
, “
Dynamical models of tuberculosis and their applications
,”
Math. Biosci. Eng.
1
(
2
),
361
404
(
2004
).
26.
Castillo-Chavez
,
C.
,
Zhilan
,
F.
, and
Huang
,
W.
, “
On the computation of R0 and its role on global stability
,” in
Mathematical Approaches for Emerging and Re Emerging Infectious Diseases: An Introduction
(
Springer-Verlag
,
New York
,
2002
).
27.
Chinnathambi
,
R.
,
Rihan
,
F. A.
, and
Alsakaji
,
H. J.
, “
A fractional-order model with time delay for tuberculosis with endogenous reactivation and exogenous reinfections
,”
Math. Methods Appl. Sci.
44
(
10
),
8011
8025
(
2021
).
28.
Danane
,
J.
,
Allali
,
K.
, and
Hammouch
,
Z.
, “
Mathematical analysis of a fractional differential model of HBV infection with antibody immune response
,”
Chaos, Solitons Fractals
136
,
109787
(
2020
).
29.
Diekman
,
O.
,
Heesterbeek
,
J. A. P.
, and
Metz
,
J. A. J.
, “
On the definition and computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations
,”
J. Math. Biol.
28
,
365
382
(
1990
).
30.
Din
,
A.
,
Li
,
Y.
,
Khan
,
F. M.
,
Khan
,
Z. U.
, and
Liu
,
P.
, “
On analysis of fractional order mathematical model of Hepatitis B using Atangana-Baleanu Caputo (ABC) derivative
,”
Fractals
30
(
1
),
2240017
(
2022
).
31.
Dye
,
C.
and
Wolpert
,
D.
, “
Earthquakes, influenza and cycles of Indian kala-azar
,”
Trans. R. Soc. Trop. Med. Hyg.
82
(
6
),
843
850
(
1988
).
32.
Du
,
M.
,
Wang
,
Z.
, and
Hu
,
H.
, “
Measuring memory with the order of fractional derivative
,”
Sci. Rep.
3
(
1
),
3431
(
2013
).
33.
Geweke
,
J.
, “
Evaluating the Accuracy of Sampling-bases approach to calculation of posterior moment
,”
Report No. 148
(
Federal Reserve Bank of Minneapolis
,
1991
).
34.
Ghosh
,
K.
,
Samanta
,
S.
,
Biswas
,
S.
,
Rana
,
S.
,
Elmojtaba
,
I.
,
Kesh
,
D.
, and
Chattopadhyay
,
J.
, “
Stability and bifurcation analysis of an eco-epidemiological model with multiple delays
,”
Nonlinear Stud.
23
(
2
),
167
208
(
2016
).
35.
Haario
,
H.
,
Laine
,
M.
,
Mira
,
A.
, and
Saksman
,
E.
, “
DRAM: Efficient adaptive MCMC
,”
Stat. Comput.
16
(
4
),
339
354
(
2006
).
36.
Haario
,
H.
,
Saksman
,
E.
, and
Tamminen
,
J.
, “
An adaptive Metropolis algorithm
,”
Bernoulli
7
,
223
242
(
2001
).
37.
Huang
,
Y.
and
Yu
,
L.
, “
Mathematical model of a hydraulic retarder based on Rankine-vortex dynamics
,”
Phys. Fluids
35
(
10
),
107127
(
2023
).
38.
Jajarmi
,
A.
and
Baleanu
,
D.
, “
A new fractional analysis on the interaction of HIV with CD4+ T-cells
,”
Chaos, Solitons Fractals
113
,
221
229
(
2018
).
39.
Jan
,
R.
,
Khan
,
M. A.
,
Kumam
,
P.
, and
Thounthong
,
P.
, “
Modeling the transmission of dengue infection through fractional derivatives
,”
Chaos, Solitons Fractals
127
,
189
216
(
2019
).
40.
Khan
,
A.
,
Khan
,
H.
,
Gómez-Aguilar
,
J. F.
, and
Abdeljawad
,
T.
, “
Existence and Hyers-Ulam stability for a nonlinear singular fractional differential equations with Mittag-Leffler kernel
,”
Chaos, Solitons Fractals
127
,
422
427
(
2019
).
41.
Khan
,
M. A.
and
Atangana
,
A.
, “
Mathematical modeling and analysis of COVID-19: A study of new variant Omicron
,”
Physica A
599
,
127452
(
2022
).
42.
Kumar
,
D.
,
Singh
,
J.
,
Al Qurashi
,
M.
, and
Baleanu
,
D.
, “
A new fractional SIRS-SI malaria disease model with application of vaccines, antimalarial drugs, and spraying
,”
Adv. Differ. Equations
2019
(
1
),
19
.
43.
Kumar
,
S.
,
Chauhan
,
R. P.
,
Aly
,
A. A.
,
Momani
,
S.
, and
Hadid
,
S.
, “
A study on fractional HBV model through singular and non-singular derivatives
,”
Eur. Phys. J. Spec. Top.
231
(
10
),
1885
1904
(
2022
).
44.
Laine
,
M.
, “
Adaptive MCMC methods with applications in environmental and geophysical models
,” Ph.D. thesis (
Finnish Meteorological Institute
,
2008
).
45.
Leandry
,
L.
, “
An investigation on the monkeypox virus dynamics in human and rodent populations for a deterministic mathematical model
,” Informatics in Medicine Unlocked 41, 101325 (
2023
).
46.
Lenhart
,
S.
and
Workman
,
J.
,
Optimal Control Applied to Biological Models
(
Chapman and Hall/CRC
,
2007
).
47.
Li
,
X. P.
,
Al Bayatti
,
H.
,
Din
,
A.
, and
Zeb
,
A.
, “
A vigorous study of fractional order COVID-19 model via ABC derivatives
,”
Results Phys.
29
,
104737
(
2021
).
48.
Lia
,
Y.
,
Haq
,
F.
,
Shah
,
K.
,
Shahzad
,
M.
, and
Rahman
,
G.
, “
Numerical analysis of fractional order pine wilt disease model with bilinear incident rate
,”
J. Math. Comput. Sci.
17
,
420
428
(
2017
).
49.
Liu
,
K.
,
Wang
,
J.
,
Zhou
,
Y.
, and
O'Regan
,
D.
, “
Hyers-Ulam stability and existence of solutions for fractional differential equations with Mittag-Leffler kernel
,”
Chaos, Solitons Fractals
132
,
109534
(
2020
).
50.
Majee
,
S.
,
Jana
,
S.
, and
Kar
,
T. K.
, “
Dynamical analysis of monkeypox transmission incorporating optimal vaccination and treatment with cost-effectiveness
,”
Chaos
33
(
4
),
043103
(
2023
).
51.
Marino
,
S.
,
Hogue
,
I. B.
,
Ray
,
C. J.
, and
Kirschner
,
D.
, “
A methodology for performing global uncertainty and sensitivity analysis in systems biology
,”
J. Theor. Biol.
254
,
178
196
(
2008
).
52.
Moore
,
E. J.
,
Sirisubtawee
,
S.
, and
Koonprasert
,
S.
, “
A Caputo-Fabrizio fractional differential equation model for HIV/AIDS with treatment compartment
,”
Adv. Differ. Equations
2019
(
1
),
1
20
.
53.
Muhammad Altaf
,
K.
and
Atangana
,
A.
, “
Dynamics of Ebola disease in the framework of different fractional derivatives
,”
Entropy
21
(
3
),
303
(
2019
).
54.
Myung
,
I. J.
, “
Tutorial on maximum likelihood estimation
,”
J. Math. Psychol.
47
(
1
),
90
100
(
2003
).
55.
Odibat
,
Z.
, “
Approximations of fractional integrals and Caputo fractional derivatives
,”
Appl. Math. Comput.
178
(
2
),
527
533
(
2006
).
56.
Odibat
,
Z. M.
, “
Analytic study on linear systems of fractional differential equations
,”
Comput. Math. Appl.
59
(
3
),
1171
1183
(
2010
).
57.
Okosun
,
K.
and
Makinde
,
O.
, “
A co-infection model of malaria and cholera diseases with optimal control
,”
Math. Biosci.
258
,
19
32
(
2014
).
58.
Okyere
,
S.
and
Ackora-Prah
,
J.
, “
Modeling and analysis of monkeypox disease using fractional derivatives
,”
Results Eng.
17
,
100786
(
2023
).
59.
Pal
,
N.
,
Samanta
,
S.
,
Biswas
,
S.
,
Alquran
,
M.
,
Al-Khaled
,
K.
, and
Chattopadhyay
,
J.
, “
Stability and bifurcation analysis of a three-species food chain model with delay
,”
Int. J. Bifurcation Chaos
25
(
9
),
1550123
(
2015
).
60.
Pandey
,
A.
,
Mubayi
,
A.
, and
Medlock
,
J.
, “
Comparing vector-host and sir models for dengue transmission
,”
Math. Biosci.
246
,
252
259
(
2013
).
61.
Panhwer
,
L. A.
,
Abro
,
K. A.
, and
Memon
,
I. Q.
, “
Thermal deformity and thermolysis of magnetized and fractional Newtonian fluid with rheological investigation
,”
Phys. Fluids
34
(
5
),
053115
(
2022
).
62.
Parmar
,
D.
,
Rathish Kumar
,
B. V.
,
Krishna Murthy
,
S. V. S. S. N. V. G.
, and
Kumar
,
S.
, “
Numerical study of entropy generation in magneto-convective flow of nanofluid in porous enclosure using fractional order non-Darcian model
,”
Phys. Fluids
35
(
9
),
097142
(
2023
).
63.
Peter
,
O. J.
,
Abidemi
,
A.
,
Ojo
,
M. M.
, and
Ayoola
,
T. A.
, “
Mathematical model and analysis of monkeypox with control strategies
,”
Eur. Phys. J. Plus
138
(
3
),
242
(
2023
).
64.
Peter
,
O. J.
,
Kumar
,
S.
,
Kumari
,
N.
,
Oguntolu
,
F. A.
,
Oshinubi
,
K.
, and
Musa
,
R.
, “
Transmission dynamics of Monkeypox virus: A mathematical modelling approach
,”
Model. Earth Syst. Environ.
8
,
3423
3434
(
2022
).
65.
Peter
,
O. J.
,
Madubueze
,
C. E.
,
Ojo
,
M. M.
,
Oguntolu
,
F. A.
, and
Ayoola
,
T. A.
, “
Modeling and optimal control of monkeypox with cost-effective strategies
,”
Model. Earth Syst. Environ.
9
(
2
),
1989
2007
(
2023
).
66.
Peter
,
O. J.
,
Oguntolu
,
F. A.
,
Ojo
,
M. M.
,
Olayinka Oyeniyi
,
A.
,
Jan
,
R.
, and
Khan
,
I.
, “
Fractional order mathematical model of monkeypox transmission dynamics
,”
Phys. Scr.
97
(
8
),
084005
(
2022
).
67.
Peter
,
O. J.
,
Yusuf
,
A.
,
Ojo
,
M. M.
,
Kumar
,
S.
,
Kumari
,
N.
, and
Oguntolu
,
F. A.
, “
A mathematical model analysis of meningitis with treatment and vaccination in fractional derivatives
,”
Int. J. Appl. Comput. Math.
8
(
3
),
117
(
2022
).
68.
Podlubny
,
I.
, “
An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications
,”
Math. Sci. Eng.
198
,
340
(
1999
).
69.
Portet
,
S.
, “
A primer on model selection using the Akaike Information Criterion
,”
Infect. Dis. Modell.
5
,
111
128
(
2020
).
70.
Prasad
,
B.
and
Bali
,
R.
, “
Mathematical study of nanoparticle loaded in red blood cells for drug delivery in an artery with stenosis
,”
Phys. Fluids
35
(
9
),
091902
(
2023
).
71.
Qureshi
,
S.
and
Jan
,
R.
, “
Modeling of measles epidemic with optimized fractional order under Caputo differential operator
,”
Chaos, Solitons Fractals
145
,
110766
(
2021
).
72.
Rahman
,
M. U.
,
Arfan
,
M.
,
Shah
,
Z.
,
Kumam
,
P.
, and
Shutaywi
,
M.
, “
Nonlinear fractional mathematical model of tuberculosis (TB) disease with incomplete treatment under Atangana-Baleanu derivative
,”
Alexandria Eng. J.
60
(
3
),
2845
2856
(
2021
).
73.
Rezapour
,
S.
,
Etemad
,
S.
, and
Mohammadi
,
H.
, “
A mathematical analysis of a system of Caputo-Fabrizio fractional differential equations for the anthrax disease model in animals
,”
Adv. Differ. Equations
2020
(
1
),
481
.
74.
Roy
,
N. C.
and
Pop
,
I.
, “
Exact solutions of Stokes' second problem for hybrid nanofluid flow with a heat source
,”
Phys. Fluids
33
(
6
),
063603
(
2021
).
75.
Saifuddin
,
M.
,
Biswas
,
S.
,
Samanta
,
S.
,
Sarkar
,
S.
, and
Chattopadhyay
,
J.
, “
Complex dynamics of an eco-epidemiological model with different competition coefficients and weak Allee in the predator
,”
Chaos, Solitons Fractals
91
,
270
285
(
2016
).
76.
Saifuddin
,
M.
,
Samanta
,
S.
,
Biswas
,
S.
, and
Chattopadhyay
,
J.
, “
An eco-epidemiological model with different competition coefficients and strong-Allee in the prey
,”
Int. J. Bifurcation Chaos
27
(
8
),
1730027
(
2017
).
77.
Samanta
,
S.
, “
Effects of awareness program and delay in the epidemic outbreak
,”
Math. Methods Appl. Sci.
40
(
5
),
1679
1695
(
2017
).
78.
Samanta
,
S.
,
Rana
,
S.
,
Sharma
,
A.
,
Misra
,
A.
, and
Chattopadhyay
,
J.
, “
Effect of awareness programs by media on the epidemic outbreaks: A mathematical model
,”
Appl. Math. Comput.
219
(
12
),
6965
6977
(
2013
).
79.
Samko
,
S. G.
,
Kilbas
,
A. A.
, and
Marichev
,
O. I.
,
Fractional Integrals and Derivatives: Theory and Applications
(
Gordon and Breach
,
1993
).
80.
Somma
,
S. A.
,
Akinwande
,
N. I.
, and
Chado
,
U. D.
, “
A mathematical model of monkey pox virus transmission dynamics
,”
Ife J. Sci.
21
(
1
),
195
204
(
2019
).
81.
Solís-Pérez
,
J. E.
,
Gómez-Aguilar
,
J. F.
, and
Atangana
,
A.
, “
Novel numerical method for solving variable-order fractional differential equations with power, exponential and Mittag-Leffler laws
,”
Chaos, Solitons Fractals
114
,
175
185
(
2018
).
82.
Srivastava
,
H. M.
,
Jan
,
R.
,
Jan
,
A.
,
Deebani
,
W.
, and
Shutaywi
,
M.
, “
Fractional-calculus analysis of the transmission dynamics of the dengue infection
,”
Chaos
31
(
5
),
053130
(
2021
).
83.
Stauch
,
A.
,
Sarkar
,
R.
,
Picado
,
A.
,
Ostyn
,
B.
,
Sundar
,
S.
,
Rijal
,
S.
,
Boelaert
,
M.
,
Dujardin
,
J.
, and
Duerr
,
H.
, “
Visceral leishmaniasis in the Indian subcontinent: Modelling epidemiology and control
,”
PLoS Negl. Trop. Dis.
5
(
11
),
e1405
(
2011
).
84.
Tasman
,
H.
, “
An optimal control strategy to reduce the spread of malaria resistance
,”
Math. Biosci.
262
,
73
79
(
2015
).
85.
TeWinkel
,
R. E.
, “
Stability analysis for the equilibria of a monkeypox model,” Ph.D. dissertation
(
The University of Wisconsin-Milwaukee
,
2019
).
86.
Toufik
,
M.
and
Atangana
,
A.
, “
New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models
,”
Eur. Phys. J. Plus
132
,
1
16
(
2017
).
87.
The World Bank
, see https://data.worldbank.org/ for the human mortality rate per day in the United States (2023).
87.
U.S. Department of Commerce, see https://www.commerce.gov for the initial total human population in the United States (2023).
88.
Usman
,
S.
and
Adamu
,
I. I.
, “
Modeling the transmission dynamics of the monkeypox virus infection with treatment and vaccination interventions
,”
J. Appl. Math. Phys.
5
(
12
),
2335
(
2017
).
89.
Van den Driessche
,
P.
and
Watmough
,
J.
, “
Reproduction numbers and the sub-threshold endemic equilibria for compartmental models of disease transmission
,”
Math. Biosci.
180
,
29
48
(
2002
).
90.
Vijayalakshmi
,
G. M.
and
Roselyn Besi
,
P.
, “
ABC fractional order vaccination model for Covid-19 with self-protective measures
,”
Int. J. Appl. Comput. Math.
8
(
3
),
130
(
2022
).
91.
Yu
,
W.
,
Song
,
S.
, and
Choi
,
J. I.
, “
Numerical simulations of underwater explosions using a compressible multi-fluid model
,”
Phys. Fluids
35
(
10
),
106102
(
2023
).
92.
Zarin
,
R.
,
Khan
,
A.
,
Akgül
,
A.
, and
Akgül
,
E. K.
, “
Fractional modeling of COVID-19 pandemic model with real data from Pakistan under the ABC operator
,”
AIMS Math.
7
(
9
),
15939
15964
(
2022
).
93.
Zeng
,
D.
and
Lin
,
D. Y.
, “
Maximum likelihood estimation in semiparametric regression models with censored data
,”
J. R. Stat. Soc., Ser. B
69
(
4
),
507
564
(
2007
).
94.
Zuo
,
C.
,
Zhu
,
F.
, and
Ling
,
Y.
, “
Analyzing COVID-19 vaccination behavior using an SEIRM/V epidemic model with awareness decay
,”
Front. Public Health
10
,
817749
(
2022
).