Analysis of a mathematical model for the spreading of the monkeypox virus with constant proportional-Caputo derivative operator

1.   Introduction

Presently, people in many countries face a series of epidemics, each caused by a different type of emerging or reemerging virus. Monkeypox ($\mathbb{MPOX}$) disease is one of these urgent situations that need to be monitored appropriately since the World Health Organization (WHO) has recently reported a resurgence of the $\mathbb{MPOX}$ outbreak, with the disease reemerging and spreading across multiple countries worldwide. Numerous cases and clusters have been identified simultaneously in diverse geographical areas ^[1]. The $\mathbb{MPOX}$ virus is the cause of a zoonotic illness. The spreading of this virus to humans mainly comes through bites or scratches from wild animals, such as rodents and primates. Human-to-human transmission arises through respiratory droplets or contact with body fluids when a person touches a lesion on an infected individual or their items. The symptoms of $\mathbb{MPOX}$ patients frequently resemble those of smallpox but are less severe. It begins with a slight fever, chills, weariness, muscle aches, exhaustion, and headache. The $\mathbb{MPOX}$ can cause swelling in the lymph nodes, while smallpox does not cause lymphadenopathy. The time between the onset of symptoms and the rash completely recovering typically lasts 14–28 days ^[2,3,4]. Currently, there is no vaccine that directly protects against $\mathbb{MPOX}$ disease. However, due to the similarity between the virus that causes $\mathbb{MPOX}$ disease and smallpox disease, vaccination against smallpox is the best choice. It effectively prevents $\mathbb{MPOX}$ disease by as much as 80–85 percent ^[5]. Mathematical models of epidemics have long been used to help us understand the effects of numerous disease transmissions in the real world. Understanding the virus and its transmission dynamics is critical for establishing efficient prevention and control strategies. Therefore, many researchers have focused on developing the $\mathbb{MPOX}$ epidemic model for various reasons, including potential human outbreaks. Some examples of interesting works, such as the analysis of the $\mathbb{MPOX}$ disease with the impact of vaccination, were studied in ^[6,7] A quarantine class and an enlightenment campaign were incorporated into the $\mathbb{MPOX}$ model, which can be found in ^[8], while the factor of isolated humans was added to the model, as seen in ^[9,10,11]. In ^[12], pox-like analysis was investigated under the factor of recovery with permanent immunity. Monkeypox analysis using data and statistical tools were examined in ^[13,14].

Fractional calculus is an area of pure mathematics that generalizes the concepts of differentiation and integration involving non-integer derivatives and integrals. It is a valuable and efficient tool for describing complex dynamical systems and simulating real-world problems, particularly in epidemiology, since it has hereditary properties and describes memory in the context of numerous diseases. Besides, fractional models provide a more realistic illness trajectory. Fractional differential operators enhance epidemiological modeling by combining memory, heterogeneity, and nonlocal dynamics. This leads to more accurate predictions, a better understanding of diseases, and improved control strategies for illnesses such as monkeypox, COVID-19, HIV, and so on. Fractional derivatives, which can take any real number as an order, are highly flexible and nonlocal. These characteristics make them more dependable than classical derivatives for approximating real data since they take into account global interactions and memory effects. Researchers have defined the many forms of fractional derivative operators in fractional calculus. Each form is created to capture distinct characteristics of the fractional derivative idea. For instance, Riemann-Liouville ($\mathbb{RL}$, one of the original definitions of fractional derivative), Caputo ^[15,16], Caputo-Katugampola ^[17], Caputo-Fabrizio ^[18,19], Atangana-Baleanu ^[18,19], Hilfer ^[20], constant proportional-Caputo $(\mathbb{CPC})$ ^[21], and $(k, \psi)$-Hilfer proportional ^[22]. The $\mathbb{CPC}$ derivative is a modern development in the fractional derivative, which is designed in ^[21]. It is a hybrid concept between Caputo's type and the proportional derivatives to gain a novel type of fractional calculus and potential applications for modeling actual data from various problems. In addition, the $\mathbb{CPC}$ operator is a more comprehensive framework than the Caputo fractional derivative operator, for which supporting research can be found in ^[23,24]. This operator offers a significant advantage by effectively addressing current challenges in different issues that traditional operators cannot adequately evaluate. Several studies ^[25,26,27] highlight that the $\mathbb{CPC}$ operator provides a more practical and accurate approach for examining mathematical models applied to real-world problems using real data compared to classical and fractional-order operators.

The application of fractional calculus to mathematical modeling is a novel approach that has captured the interest of many researchers across diverse scientific domains. It is particularly intriguing due to the memory effect, a unique and powerful characteristic of fractional-order models, which remains supereminent over the classical models due to the diverse properties of fractional operators that yield more accurate and reliable data. Many fractional models consist of differential equations carefully designed by researchers to address the pressing problems within the global environment. Extensive research has been conducted on studying and analyzing fractional-order models for disease transmission dynamics. El-Mesady et al. ^[28] designed a fractional-order vaccination mathematical model for tuberculosis incorporating a susceptible class with an underlying ailment. Peter et al. ^[29] developed a fractional mathematical model for studying measles infection with double-dose vaccination. The work ^[30] constructed a mathematical model under the Atangana-Baleanu-Caputo derivative to examine meningitis with treatment and vaccination dynamics. The authors ^[31] presented the pneumococcal pneumonia infection model using fractional-order derivatives in the sense of the Caputo-Fabrizio operator. Peter ^[32] studied the transmission dynamics of a Brucellosis model under the Caputo-Fabrizio fractional operator. Additionally, $\mathbb{MPOX}$ virus infection is one of the research areas that has gotten a lot of attention and produced some fascinating results since the disease resurfaced after a long hiatus. We refer the reader to several previously interesting works about the $\mathbb{MPOX}$ fractional model. For example, Peter et al. ^[33] presented and established the dynamical behavior of the $\mathbb{MPOX}$ virus model by using both classical and differential equations via the Caputo-Fabrizio fractional derivative. Ngungu et al. ^[34] studied the dynamics of the $\mathbb{MPOX}$ virus spreading with a non-pharmaceutical intervention using real-time data with the Caputo-Fabrizio operator. Wireko et al. ^[35] used fractal-fractional operators $(\mathbb{FFO}s)$ to explore the biological behavior of the $\mathbb{MPOX}$ disease. Sudsutad et al. ^[36] investigated the theoretical analysis for the transmission of the $\mathbb{MPOX}$ virus fractional model under the $\mathbb{FFO}s$ involving the Atangana-Baleanu sense. The study by El-Mesady et al. ^[37] looked into how the $\mathbb{MPOX}$ virus spreads in human hosts and rodent populations using a Caputo fractional-order nonlinear model. The $\mathbb{MPOX}$ virus is established by applying a deterministic mathematical model in the context of the Atangana-Baleanu fractional derivative that depends on the generalized Mittag-Leffler kernel ^[38]. Zhang et al. ^[39] studied a deterministic Caputo fractional-order mathematical model of Marburg-$\mathbb{MPOX}$ virus co-infection transmission. Liu et al. ^[40] analyzed the dynamics of a $\mathbb{MPOX}$ disease with the impact of vaccination utilizing a fractional mathematical model. For more works, we refer readers to see ^[41,42]. However, even though many researchers employ fractional differential systems to understand real-world phenomena, finding exact solutions to such systems through manual methods is still a formidable challenge. Hence, several efficient techniques are produced to find the approximated solution of fractional differential systems, such as the predictor-corrector ^[43,44], Adams-Bashforth ^[43,45], Newton polynomial ^[43,45], and the decomposition formula method ^[46,47]. The last introduced powerful knowledge is Ulam stability, an essential technique in mathematical analysis and other related sciences. Ulam stability is vital to maintaining the stability and reliability of solutions under minor variations. It ensures that solutions remain valid and usable even though conditions change marginally. There are various kinds of Ulam stability commonly used, like Ulam-Hyers ($\mathbb{UH}$) stability and Ulam-Hyers-Rassias ($\mathbb{UHR}$) stability; see the history and its application in ^{[48,49,50,51]}.

After composing all of the stories and being motivated by the above discussions, the $\mathbb{CPC}$ operator is particularly interesting in this aspect because it is a recent operator, and there has been little literature-based research on its application. Also, to the best of our current understanding, no studies using the proposed derivative operator have been conducted or published in the existing literature on the dynamics of the monkeypox virus as well as the benefit of using the $\mathbb{CPC}$ operator, as mentioned above, to address memory and hereditary properties for resulting in more accurate prediction and translation. Therefore, this study takes a comprehensive approach to analyzing the behavior of the $\mathbb{MPOX}$ model and investigating the factors influencing population changes, which may help control the monkeypox outbreak. We have developed a deterministic $\mathbb{MPOX}$ model by extending the works ^[9,33] in terms of incorporating the vaccinated-individuals compartment into the model and yielding a unique classical model. Later, we transform the model into a fractional-order system using the $\mathbb{CPC}$ operator to gain a better understanding of the virus dynamics. Our focus is on studying the dynamic behavior of this model, a task we accomplish by leveraging the well-known fixed-point theory of Banach's type to prove the solution's existence and uniqueness. We investigate the stability of equilibrium points with the help of the basic reproduction number. Moreover, we ensure the solutions exist by analyzing their stability via various Ulam stability. A decomposition formula for the $\mathbb{CPC}$ derivative technique is derived to obtain the numerical scheme, and some graphics in the numerical simulation are shown to visualize the system's behavior analysis. This work will help to fill the gap in the study of monkeypox transmittance using fractional derivatives and expand the scope of this study to benefit disease control.

The remaining sections of the paper are as follows: Section 2 introduces some concepts of the $\mathbb{CPC}$ operators. The $\mathbb{MPOX}$ model construction is presented in Section 3. Section 4 is dedicated to investigating model analysis, including the solution's positiveness and boundedness, the basic reproduction number, and the local stability analysis of the equilibrium points. Section 5 explores the existence theory for the proposed model; that is, the uniqueness result is verified using Banach's fixed-point theorem. Section 6 verifies various Ulam's type stability and their generalization. The numerical scheme derived from a decomposition formula for the $\mathbb{CPC}$ derivative is determined in Section 7. Finally, the results are discussed via some examples, and the summation of this discussion is given in the final part.

2.   Preliminaries

This section provides some fundamental definitions and properties of fractional calculus, which will be used to analyze the system throughout this work.

Definition 2.1. (^[52]). The Caputo fractional order derivative of a function $f$ with order $\alpha \in (0, 1)$ is provided by

Definition 2.2. (^[52]). Assume $f(t)$ is an integrable function. The $\mathbb{RL}$-integral of $\alpha > 0$ is given by

Definition 2.3. (^[21]). A proportional-Caputo $(\mathbb{CP})$ is a hybrid operator that combines the proportional operator and the Caputo fractional derivative as:

where $\mathcal{K}_0(\alpha, t) = \alpha C^{2\alpha}t^{1-\alpha}$ and $\mathcal{K}_1(\alpha, t) = (1-\alpha)t^{\alpha}$, for $C$ is constant and $\alpha \in (0, 1)$.

Moreover, as defined in the particular case where $\mathcal{K}_0$ and $\mathcal{K}_1$ are depending only on $\alpha$, the $\mathbb{CPC}$ operator can be defined by

where $\mathcal{K}_0(\alpha)$ and $\mathcal{K}_1(\alpha)$ are constants with respect to $t$.

Here, this study uses the specific case when $\mathcal{K}_0(\alpha) = \alpha C^{2\alpha} Q^{1-\alpha}$ and $\mathcal{K}_1(\alpha) = (1-\alpha) Q^\alpha$ where $C$ and $Q$ are constants.

Definition 2.4. (^[21]). The inverse operator of the $\mathbb{CPC}$ fractional derivative operator is provided by:

which satisfies the relation below:

3.   Formulation of the mathematical models

3.1.   The integer order of the $\mathbb{MPOX}$ model

In this subsection, we formulate a deterministic model of the dynamics of $\mathbb{MPOX}$ transmission across two groups: humans denoted by $N_h$ and rodents denoted by $N_r$. The human population has six different compartments: the susceptible ($S_h$), exposed ($E_h$), infectious ($I_h$), clinically ill ($C_h$), recovered ($R_h$), and vaccinated individuals ($V_h$). The rodent population has three compartments: the susceptible ($S_r$), exposed ($E_r$), and infectious ($I_r$). Thus, the total populations of humans and rodents are given by ${N}_h(t) = {S}_h(t) + {E}_h(t) + {I}_h(t) + {C}_h(t) + {R}_h(t)+ {V}_h(t)$ and ${N}_r(t) = {S}_r(t) + {E}_r(t) + {I}_r(t)$, respectively. The transmission diagram of the population flow among these compartments is displayed in Figure 1.

Here, the descriptions of the population changes of each group are presented in the following:

● Human groups: The amount of susceptible humans changes due to recruitment through birth rate or immigration $\Lambda_h$. The vaccinated population further increased the number in this class after the induced immunity waned at the rate $\eta$. This population class is rejected by the natural death per capital rate $\mu_h$ and by the population gaining infection after contact with infected humans or infected rodents at the rate $b_1$ or $b_2$, respectively. It is also reduced by vaccinating at the rate $\lambda$. The amount of exposed humans grows with the force of infection in the term of $\phi_h = {(b_1 I_h(t)+b_2 I_r(t))}/{N_h(t)}$ and declines by the disease progression rate $a_1$ and the natural death per capital rate $\mu_h$. The amount of infected individuals increases after transitioning from the exposed class with the disease progression rate $a_1$. The natural recovery rate $\omega$ reduces the population in the infectious class due to immunity. This population is also reduced by the natural death rate $\mu_h$, the disease-induced death rate $\delta_1$, and the clinically ill rate $a_2$ after moving to a clinically ill class of humans. The number of clinically ill humans grows due to those who seek medical assistance after becoming sick with a clinically ill rate $a_2$. These individuals are decreased by disease-induced death rate $\delta_2$ and natural death rate $\mu_h$. Moreover, moving to the recovered class at the rate $\nu$ can reduce the population in this class. The group in the recovered compartment increases with the recovery rate of clinically ill humans $\nu$ and the natural recovery rate due to immunity $\omega$, while it declines with the natural death rate $\mu_h$. The susceptible individuals enter into the vaccinated class with a rate $\lambda$. The vaccinated individuals are declined by the natural death rate and the waning induced immunity with a rate $\eta$.

● Rodent groups: The amount of susceptible rodents changes due to recruitment through birth rate $\Lambda_r$. These individuals are decreased by interaction with infected rodents at a rate $b_3$ and natural death per capital rate of rodents $\mu_r$. The amount of exposed rodents grows with the force of infection in the term $\phi_r = {b_3 I_r(t)}/{N_r(t)}$. It is declined by the natural death rate $\mu_r$ and the progression rate $a_3$, which transit from exposed rodents to infected rodents. The amount of infected rodents increases from the exposed rodent's transit to infectious rodents at the rate $a_3$. The infected rodents are reduced due to the natural death rate $\mu_r$.

Therefore, based on the above description, a nonlinear system of the integer-order of the $\mathbb{MPOX}$ model corresponding with nine ordinary differential equations, the so-called $\mathbb{MPOX}$ model, is shown below:

where the positive initial conditions are ${S}_h(0) = {S}_{h_0}$, ${E}_h(0) = {E}_{h_0}$, ${I}_h(0) = {I}_{h_0}$, ${C}_h(0) = {C}_{h_0}$, ${R}_h(0) = {R}_{h_0}$, ${V}_h(0) = {V}_{h_0}$, ${S}_r(0) = {S}_{r_0}$, ${E}_r(0) = {E}_{r_0}$, and ${I}_r(0) = {I}_{r_0}$. The details of each all positive parameter are included and displayed in Table 1.

3.2.   The $\mathbb{MPOX}$ model under the $\mathbb{CPC}$ derivative operator

As we know, fractional infectious models have specific hereditary properties and describe memory regarding disease dynamics because they employ fractional derivatives with a higher degree of freedom to account for more complex diseases and provide knowledge into the behavior and control of epidemic diseases. Consequently, this subsection further develops the integer order $\mathbb{MPOX}$ model (3.1) utilizing the $\mathbb{CPC}$ derivative operator as follows:

where $\mathbb{F}_{j} = \mathbb{F}_{j}(t, {S}_h(t), {E}_h(t), {I}_h(t), {C}_h(t), {R}_h(t), {V}_h(t), {S}_r(t), {E}_r(t), {I}_r(t))$, for $j = 1, 2, \ldots, 9$ and $\mathbb{F}_{j}$ for the suggested model are given by

with $\phi_h = {(b_1 I_h(t)+b_2 I_r(t))}/{N_h(t)}$, $\phi_r = {b_3{I}_r(t)}/{{N}_r(t)}$, $c_1 = \mu_h+\lambda$, $c_2 = \mu_h+a_1$, $c_3 = \omega+a_2+\mu_h+\delta_1$, $c_4 = \nu+\mu_h+\delta_2$, $c_5 = \mu_h+\eta$, and $c_6 = a_3+\mu_r$ with ${S}_{h_{0}} \geq 0$, ${E}_{h_{0}} \geq 0$, ${I}_{h_{0}} \geq 0$, ${C}_{h_{0}} \geq 0$, ${R}_{h_{0}} \geq 0$, ${V}_{h_{0}} \geq 0$, ${S}_{r_{0}} \geq 0$, ${E}_{r_{0}} \geq 0$, and ${I}_{r_{0}} \geq 0$. The suggested model (3.2) is said to be the $\mathbb{CPC}$-$\mathbb{MPOX}$ model.

4.   Model analysis

4.1.   Positiveness and boundedness

Since the variables in real-world phenomena have positive values, particularly the epidemic model based on the human population, all variables and parameters are considered to be positive. Moreover, the invariant region is essential in mathematical modeling as it guarantees that solutions remain biologically feasible and mathematically meaningful. In this part, we investigate the factors that guarantee the solutions of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) are positive, and we locate the invariant area, ensuring that the solution is bounded. To achieve this purpose, we define

Theorem 4.1. All solutions of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) under the initial conditions are bounded in $\mathbb{R}^9_+$.

Proof. From the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2), we obtain

If a set of the conditions $(S_{h_{0}}, E_{h_{0}}, I_{h_{0}}, C_{h_{0}}, R_{h_{0}}, V_{h_{0}}, S_{r_{0}}, E_{r_{0}}, I_{r_{0}}) \in \mathbb{R}^9_+$, then according to the above system (4.1), the solution $\mathbb{F}$ cannot avoid from the hyper-planes: ($S_h$, $E_h$, $I_h$, $C_h$, $R_h$, $V_h$, $S_r$, $E_r$, $I_r$) = $(0, 0, 0, 0, 0, 0, 0, 0, 0)$. Additionally, the vector field points into $\mathbb{R}^9_+$ on each hyperplane that is surrounded by the nonnegative constants. This means that $\mathbb{R}^9_+$ is the positively invariant region. □

Theorem 4.2. Let a positive set of solutions $(S_h(t), E_h(t), I_h(t), C_h(t), R_h(t), V_h(t), S_r(t), E_r(t), I_r(t))$ under the initial conditions. Then, there exists a domain $\Omega = \Omega_h \times \Omega_r \subset \mathbb{R}^6_+ \times \mathbb{R}^3_+$, that is positively invariant so that

Proof. Since $N_h(t) = S_h(t) + E_h(t)+ I_h(t)+ C_h(t)+ R_h(t)+ V_h(t)$, then from the proposed model (3.2), we have

In the case of no disease, we get ${_{}^{\mathbb{CPC}}}\mathfrak{D}_{0, t}^{\alpha} {N}_h(t) = \Lambda_h - \mu_h N_h(t)$. It implies that

By a particular comparison principle, we obtain

Then, we get

Thus, the feasible region of human groups (4.2) is completed.

Similarly, since $N_r(t) = S_r(t)+ E_r(t)+ I_r(t)$, then

It implies that

By a particular comparison principle, we have

Then, we have

Thus, the feasible region of rodent groups (4.3) is obtained. This shows the boundedness of the solutions for the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2). Hence, the region $\Omega = \Omega_h \times \Omega_r$, defined as (4.2) and (4.3), is positively invariant. We therefore conclude that the proposed model is epidemiologically feasible and well-posed in $\Omega$. The proof is complete. □

4.2.   The stability of the equilibria

The equilibrium points of the the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) can be calculated by setting the right-hand side equal to zero. This gives two possible positive equilibria, which are

● The $\mathbb{MPOX}$ free equilibrium ($\mathfrak{E}_{MFE}^0$); this point is defined when there is no disease in the population:

Note that: $c_1 c_5 - \lambda \eta = \mu_h (\mu_h+\eta+\lambda) > 0$, $N_h^0 = S_h^0+E_h^0+I_h^0+C_h^0+R_h^0+V_h^0$, and $N_r^0 = S_r^0+E_r^0+I_r^0$.

● The $\mathbb{MPOX}$ endemic equilibrium ($\mathfrak{E}_{MEE}^*$); this point is defined when there is disease transmission in the population:

where

Next, the basic reproduction number is conducted as an epidemic indicator. It is an essential fundamental epidemiological parameter used to assess the long-term dynamics of the epidemic, which is denoted by $\mathfrak{R}_0$. It is also defined as the expected number of secondary infections generated by a single infected individual throughout their infectious period. The basic reproduction number can be used as a disease control measure because it indicates the spread of a disease in terms of transmission within a population. It helps determine whether an outbreak or epidemic will occur or the infection will eventually disappear. Here, we will calculate it by using the next-generation technique as in ^[53]. For the proposed model, the disease-free state variables are $S_h, R_h, S_r$ while the infected state variables are $E_h, I_h, C_h, V_h, E_r, I_r$. Therefore,

where $\mathcal{F}$ and $\mathcal{V}$ denote the transmission and transitions matrices, respectively. Then, the next-generation matrix is provided as below:

This implies the spectral radius of (4.2) ($\mathfrak{R}_0$)

Remark 4.3. We can notice that:

● If $\mathfrak{R}_h^0 > 1$ and $\mathfrak{R}_r^0 > 1$, then $\mathfrak{R}_0 > 1$.

● If $\mathfrak{R}_h^0 < 1$ and $\mathfrak{R}_r^0 > 1$, then $\mathfrak{R}_0 = \mathfrak{R}_r^0 > 1$.

● If $\mathfrak{R}_h^0 > 1$ and $\mathfrak{R}_r^0 < 1$, then $\mathfrak{R}_0 = \mathfrak{R}_h^0 > 1$.

● If $\mathfrak{R}_h^0 < 1$ and $\mathfrak{R}_r^0 < 1$, then $\mathfrak{R}_0 < 1$.

For the following two theorems, the stability of the equilibrium points will be analyzed. To do this, the Jacobian matrix ($\mathcal{J}$) of the system is given by

Theorem 4.4. If $\mathfrak{R}_0 < 1$, then the $\mathbb{MPOX}$ free equilibrium of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is locally asymptotically stable with the necessary and sufficient criteria:

Proof. From the asumption $\mathfrak{R}_0 < 1$, it implies that $\mathfrak{R}_0^h < 1$ and $\mathfrak{R}_0^r < 1$. By applying matrix (4.4) at the point $\mathfrak{E}^0_{MFE}$, we have

Then, the characteristic equation is computed by $|\mathcal{J}(\mathfrak{E}^0_{MFE})-\theta \hat{I}| = 0$, and $\hat{I}$ is a identity matrix. This yields the eigenvalues as the following:

Obviously, the eigenvalues $\theta_1, \theta_2, \theta_3$, and $\theta_{4, 5}$ have negative real parts. Since $\mathfrak{R}_0^h < 1$ and $\mathfrak{R}_0^r < 1$, then the eigenvalues $\theta_{6, 7}$ and $\theta_{8, 9}$ also have negative real parts, respectively. This assures the assumption (4.5) with $0 < \alpha \leq 1$. Therefore, the point $\mathfrak{E}^0_{MFE}$ is locally asymptotically stable. □

For the following theorem, we will utilize the Routh-Hurwitz criterion to investigate the local stability of endemic equilibria. Here, we will determine the conditions under which the endemic equilibrium point is locally asymptotically stable. The Jacobian matrix (4.4) yields the following matrix for $\mathfrak{E}^*_{MEE}$:

We can easily see two negative eigenvalues; one is $\omega_1 = -\mu_h$, and another one is $\omega_2 = -c_4$. For the rest of the equilibrium points, we consider the reduced matrix

This yields the associated characteristic equation as follows:

with $\epsilon_k$ denoting the coefficients of $\omega^{7-k}, k = 1, 2, 3, \ldots, 7$ after resetting the polynomial equation in the formula form. Since the proof of the local asymptotic stable needs the negative real parts of all roots of (4.6), we use the Routh-Hurwitz criteria to obtain the conditions for the stability of $\mathfrak{E}^*_{MEE}$. We define the following notation, $h_1 = {(\epsilon_1\epsilon_2-\epsilon_3)}/{\epsilon_1}$, $h_2 = {(\epsilon_1\epsilon_4-\epsilon_5)}/{\epsilon_1}$, $h_3 = {(\epsilon_1\epsilon_6-\epsilon_7)}/{\epsilon_1}$, $g_1 = {(\epsilon_3 h_1-\epsilon_1 h_2)}/{h_1}$, $g_2 = {(\epsilon_5 h_1-\epsilon_1 h_3)}/{h_1}$, $g_3 = \epsilon_7$, $d_1 = {(h_2g_1-h_1g_2)}/{g_1}$, $d_2 = {(h_3g_1-h_1g_3)}/{g_1}$, $e_1 = {(g_2d_1-g_1d_2)}/{d_1}$, $e_2 = g_3$, and $f_1 = {(d_2 e_1-d_1e_2)}/{e_1}$.

Thus, the Hurwitz assumptions concerning (4.6), which ensure that all roots have negative real parts, are as follows: $(\mathcal{H}_{1}). \; (i) \; \epsilon_1 > 0$; $(ii) \; \epsilon_7 > 0$; $(iii) \; \epsilon_1 \epsilon_2 > \epsilon_3$; $(iv) \; \epsilon_1 \epsilon_2 \epsilon_3 +\epsilon_1\epsilon_5 > \epsilon^2_1\epsilon_4 + \epsilon^2_3$; $(v) \; h_2g_1 > h_1g_2$; $(vi) \; g_2d_1 > g_1d_2$; and $(vii) \; d_2e_1 > d_1e_2$. Therefore, we can conclude this part by the following theorem:

Theorem 4.5. Suppose that the necessary and sufficient assumption $(\mathcal{H}_{1})$ of Hurwitz criteria is satisfied, then the $\mathbb{MPOX}$ endemic equilibrium of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is locally asymptotically stable.

5.   Qualitative analysis

In this part, an analysis of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) will be investigated applying Banach's contraction mapping principle ^[54]. First, the Banach space on $[0, T]$ for all continuous real-valued functions is given by $\mathbb{X} = \mathcal{C}(\mathcal{I}\times {{\mathbb{R}}^{9}}, \mathbb{R})$ under the norm $\left\| \mathcal{Y} \right\| = \left\| ({{S}_{h}}, {{E}_{h}}, {{I}_{h}}, {{C}_{h}}, {{R}_{h}}, {{V}_{h}}, {{S}_{r}}, {{E}_{r}}, {{I}_{r}})\right\| = \left\| {{S}_{h}} \right\|+\left\| {{E}_{h}} \right\|+\left\| {{I}_{h}} \right\|+\left\| {{C}_{h}} \right\|+\left\| {{R}_{h}} \right\|+\left\| {{V}_{h}} \right\|+\left\| {{S}_{r}} \right\|+\left\| {{E}_{r}} \right\|+\left\| {{I}_{r}} \right\|$, where ${{S}_{h}}$, ${{E}_{h}}$, ${{I}_{h}}$, ${{C}_{h}}$, ${{R}_{h}}$, ${{V}_{h}}$, ${{S}_{r}}$, ${{E}_{r}}$, ${{I}_{r}}\in \mathbb{X}$ and $\left\|{{S}_{h}} \right\| = \sup_{t\in[0, T]}\left| {{S}_{h}}(t) \right| = {{\mathfrak{B}}_{{{h}_{1}}}}$, $\left\| {{E}_{h}} \right\| = \sup_{t\in[0, T]}\left| {{E}_{h}}(t) \right| = {{\mathfrak{B}}_{{{h}_{2}}}}$, $\left\| {{I}_{h}} \right\| = \sup_{t\in[0, T]}\left| {{I}_{h}}(t) \right| = {{\mathfrak{B}}_{{{h}_{3}}}}$, $\left\| {{C}_{h}} \right\| = \sup_{t\in[0, T]}\left| {{C}_{h}}(t) \right| = {{\mathfrak{B}}_{{{h}_{4}}}}$, $\left\| {{R}_{h}} \right\| = \sup_{t\in[0, T]}\left| {{R}_{h}}(t) \right| = {{\mathfrak{B}}_{{{h}_{5}}}}$, $\left\| {{V}_{h}} \right\| = \sup_{t\in[0, T]}\left| {{V}_{h}}(t) \right| = {{\mathfrak{B}}_{{{h}_{6}}}}$, $\left\| {{S}_{r}} \right\| = \sup_{t\in[0, T]}\left| {{S}_{r}}(t) \right| = {{\mathfrak{B}}_{{{r}_{1}}}}$, $\left\| {{E}_{r}} \right\| = \sup_{t\in[0, T]}\left| {{E}_{r}}(t) \right| = {{\mathfrak{B}}_{{{r}_{2}}}}$, and $\left\| {{I}_{r}} \right\| = \sup_{t\in[0, T]}\left| {{I}_{r}}(t) \right| = {{\mathfrak{B}}_{{{r}_{3}}}}$. Next, assume $\mathbb{F}\in \mathbb{X}$ and $\mathcal{Y}\in \mathcal{C}([0, T], \mathbb{R})$, and the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) can be presented as

where

when ${{\mathbb{G}}_{i}}, \, \, i = 1, \, \, 2, \ldots, \, \, 9$ are given by (3.3). Then, the problem (5.1) can be written applying Definition 2.4 as below

From the problem (5.1), the Eq (5.2) can be presented in the integral form:

Define an operator $\mathcal{F}:\mathbb{X}\to \mathbb{X}$ where $\mathcal{F} = ({{\mathcal{F}}_{1}}, {{\mathcal{F}}_{2}}, {{\mathcal{F}}_{3}}, {{\mathcal{F}}_{4}}, {{\mathcal{F}}_{5}}, {{\mathcal{F}}_{6}}, {{\mathcal{F}}_{7}}, {{\mathcal{F}}_{8}}, {{\mathcal{F}}_{9}})$. In consideration of (5.3) to (5.11), we obtain

Next, the transformation of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) to $\mathcal{Y} = \mathcal{F}\mathcal{Y}$ that is a fixed point problem, will be later applied with a fixed point theory to prove that the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) has a solution.

Theorem 5.1. Suppose that $\mathbb{F}\in \mathbb{X}$ corresponding with the assumption $(A_1)$ is as follows:

$(A_{1})$ There is a constant $\mathcal{N}_{\max} = \max\{ \mathcal{N}_{1}, \mathcal{N}_{2}, \mathcal{N}_{3}, \mathcal{N}_{4}, \mathcal{N}_{5}, \mathcal{N}_{6}, \mathcal{N}_{7}, \mathcal{N}_{8}, \mathcal{N}_{9}\} > 0$, such that

where $i = 1, 2, 3, \ldots, 9$ and ${S}_h$, ${E}_h$, ${I}_h$, ${C}_h$, ${R}_h$, ${V}_h$, ${S}_r$, ${E}_r$, ${I}_r \in \mathbb{X}$, and $t \in [0, T]$.

If

then the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) has a unique solution.

Proof. Define a bounded, closed, and convex subset ${{D}_{{{r}_{a}}}} : = \{({{S}_{h}}, {{E}_{h}}, {{I}_{h}}, {{C}_{h}}, {{R}_{h}}, {{V}_{h}}, {{S}_{r}}, {{E}_{r}}, {{I}_{r}}) \in \mathbb{X}:\|({{S}_{h}}, {{E}_{h}}, {{I}_{h}}, {{C}_{h}}, {{R}_{h}}, {{V}_{h}}, {{S}_{r}}, {{E}_{r}}, {{I}_{r}})\|\le {{r}_{a}}\}$ with a radius $r_a$ defined as

where ${{\mathcal{K}}_{\max}} = \max \{{{S}_{{{h}_{0}}}}$, ${{E}_{{{h}_{0}}}}$, ${{I}_{{{h}_{0}}}}$, ${{C}_{{{h}_{0}}}}$, ${{R}_{{{h}_{0}}}}$, ${{V}_{{{h}_{0}}}}$, ${{S}_{{{r}_{0}}}}$, ${{E}_{{{r}_{0}}}}$, ${{I}_{{{r}_{0}}}}\}$ and $\mathbb{F}_{\max }^{*} = \max\{\mathbb{F}_{1}^{*}$, $\mathbb{F}_{2}^{*}$, $\mathbb{F}_{3}^{*}$, $\mathbb{F}_{4}^{*}$, $\mathbb{F}_{5}^{*}$, $\mathbb{F}_{6}^{*}$, $\mathbb{F}_{7}^{*}$, $\mathbb{F}_{8}^{*}$, $\mathbb{F}_{9}^{*}\}$. Let $\sup_{t\in [0, T]}|{{\mathbb{F}}_{i}}(s, 0)| = \mathbb{F}_{i}^{*} < +\infty$, $i = 1$, $2$, $3$, $\ldots$, $9$. The process is divided into two parts.

Step I. We prove that $\mathcal{F} D_{r_{a}} \subset D_{r_{a}}$.

For any $({{S}_{h}(t)}$, ${{E}_{h}(t)}$, ${{I}_{h}(t)}$, ${{C}_{h}(t)}$, ${{R}_{h}(t)}$, ${{V}_{h}(t)}$, ${{S}_{r}(t)}$, ${{E}_{r}(t)}$, ${{I}_{r}(t)})\in {{D}_{{{r}_{a\, }}}}$, $t\in[0, T]$, we have

In the same process, we also have

This yields that

Therefore, $\mathcal{F} D_{r_{a}} \subset D_{r_{a}}$.

Step II. We prove that $\mathcal{F}$ is a contraction.

Let $({{S}_{h}}$, ${{E}_{h}}$, ${{I}_{h}}$, ${{C}_{h}}$, ${{R}_{h}}$, ${{V}_{h}}$, ${{S}_{r}}$, ${{E}_{r}}$, ${{I}_{r}})\in {{D}_{{{r}_{a}}}}$ and $(S_{h}^{*}$, $E_{h}^{*}$, $I_{h}^{*}$, $C_{h}^{*}$, $R_{h}^{*}$, $V_{h}^{*}$, $S_{r}^{*}$, $E_{r}^{*}$, $I_{r}^{*})\in {{D}_{{{r}_{a}}}}$, for $t \in [0, T]$. We obtain that

This yields that,

Likewise, we obtain

Since $\mathcal{F} = ({{\mathcal{F}}_{1}}, {{\mathcal{F}}_{2}}, {{\mathcal{F}}_{3}}, {{\mathcal{F}}_{4}}, {{\mathcal{F}}_{5}}, {{\mathcal{F}}_{6}}, {{\mathcal{F}}_{7}}, {{\mathcal{F}}_{8}}, {{\mathcal{F}}_{9}})$ and ${{\mathcal{N}}_{\max }} > 0$, then

Under the assumption (5.12), we can conclude that $\mathcal{F}$ is a contraction. Then, we get $\mathcal{F}$ has a unique fixed point. Therefore, the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) has a unique solution. □

Next, we investigate some sufficient criteria of the Ulam's stability for the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2). The definitions of these types and some necessary remarks are provided below.

Assume that ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ is a constant and ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in \mathcal{C}([0, T], {{\mathbb{R}}^{+}})$. The inequalities are given:

where $t\in [0, T]$ and ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{\kappa }_{{{\mathbb{F}}_{{{\mathcal{Y}}_{j}}}}}})}^{T}}$, for $j = 1, 2, \ldots, 9$.

Definition 5.2. The $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is called $\mathbb{UH}$ stable if there exists a constant ${{\Phi }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ such that for every ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ and each solution ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ of (5.13), there exist a solution $\mathcal{Y} \in \mathbb{X}$ of (3.2) with

where $t\in [0, T]$ and ${{\Phi }_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{\Phi }_{{{\mathbb{F}}_{{{\mathcal{Y}}_{i}}}}}})}^{T}}$, $i = 1, 2, \ldots, 9$.

Definition 5.3. The $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is called generalized $\mathbb{UH}$ stable if there exists a function ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in \mathcal{C}([0, T], {{\mathbb{R}}^{+}})$, with ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}(0) = 0$ so that for ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ and for each solution ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ of (5.14), there exist a solution $\mathcal{Y} \in \mathbb{X}$ of (3.2) with

where $t\in [0, T],$ and ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}})}^{T}}, \; i = 1, 2, \ldots, 9$.

Definition 5.4. The $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is called $\mathbb{UHR}$ stable with respect to ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in \mathcal{C}([0, T], {{\mathbb{R}}^{+}})$ if there exists a number ${{\Omega }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ so that for every ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ and for each solution ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ of (5.15), there exists a solution $\mathcal{Y} \in \mathbb{X}$ of (3.2) with

where $t\in [0, T], {{\Omega }_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{\Omega }_{{{\mathbb{F}}_{{{\mathcal{Y}}_{i}}}}}})}^{T}},$ and ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{P}_{{{\mathbb{F}}_{\mathcal{Y}}}_{i}}})}^{T}}, \; i = 1, 2, \ldots, 9$.

Definition 5.5. The $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is called generalized $\mathbb{UHR}$ stable with respect to ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in \mathcal{C}([0, T], {{\mathbb{R}}^{+}})$ if there exists a number ${{\Omega }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$ so that for each solution ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ of (5.15), there exists a solution $\mathcal{Y} \in \mathbb{X}$ of (3.2) with

where $t\in [0, T], {{\Omega }_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{\Omega }_{{{\mathbb{F}}_{{{\mathcal{Y}}_{i}}}}}})}^{T}},$ and ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{P}_{{{\mathbb{F}}_{\mathcal{Y}}}_{i}}})}^{T}}$, $i = 1, 2, \ldots, 9$.

Remark 5.6. ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ is the solution of (5.13), if and only if, there is ${{\chi}_{\mathcal{Y}}}(0) = 0$ satisfied $\left| {{\chi }_{\mathcal{Y}}}(t) \right|\le {{\kappa }_{{{\mathbb{G}}_{\mathcal{Y}}}}}$, ${{\chi }_{\mathcal{Y}}} = \max {{({{\chi }_{{{\mathcal{Y}}_{i}}}})}^{T}}, i = 1, 2, \ldots, 9, \; {{\kappa }_{{{F}_{\mathcal{Y}}}}} > 0$, and $^{\mathbb{C}\mathbb{P}\mathbb{C}}\mathfrak{D}_{0, t}^{\alpha }\mathcal{Y}(t) = \mathbb{F}(t, \mathcal{Y}(t))+{{\chi }_{\mathcal{Y}}}(t)$.

Remark 5.7. ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ is the solution of (5.14), if and only if, there is ${{\psi }_{\mathcal{Y}}}\in \mathbb{X}$ satisfied $\left| {{\psi }_{\mathcal{Y}}}(t) \right|\le {{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}{{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}(t), \, \, {{\psi }_{\mathcal{Y}}} = \max {{({{\psi }_{{{\mathcal{Y}}_{i}}}})}^{T}}, \, \, {{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}} = \max {{({{P}_{{{\mathbb{F}}_{{{\mathcal{Y}}_{i}}}}}})}^{T}}$, $i = 1, 2, \ldots, 9$, and $^{\mathbb{C}\mathbb{P}\mathbb{C}}\mathfrak{D}_{0, t}^{\alpha }\mathcal{Y}(t) = \mathbb{F}(t, \mathcal{Y}(t))+{{\psi }_{\mathcal{Y}}}(t)$.

Lemma 5.8. Assume that ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ is a solution of (5.13), then

Proof. Let ${{Z}_{\mathcal{Y}}}$ be the solution of (5.13). From Remark (5.6), we have

Thus, the solution to (5.21) is provided below:

Therefore,

The proof is completed. □

Theorem 5.9. Assume the conditions in Theorem 5.1 and Lemma 5.8 are satisfied, then the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is $\mathbb{UH}$ stable.

Proof. Assume that ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in {{\mathbb{R}}^{+}}$ and ${{Z}_{\mathcal{Y}}}$ is a solution of (5.13). Let $\mathcal{Y}\in \mathbb{X}$ be a unique solution of (3.2), then

Setting

this yields that $\left| {{Z}_{\mathcal{Y}}}(t)-\mathcal{Y}(t) \right| \le {{\Phi }_{{{\mathbb{F}}_{\mathcal{Y}}}}}{{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}$. Hence, by Definition 5.2, the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is $\mathbb{UH}$ stable. □

Corollary 5.10. Setting ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}({{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}) = {{\Phi }_{{{\mathbb{F}}_{\mathcal{Y}}}}}{{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}$ with ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}(0) = 0$ in Theorem 5.9, by Definition 5.3, the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is $\mathbb{UH}$ stable.

The following conditions are necessary to prove the $\mathbb{UHR}$ and generalized $\mathbb{UHR}$ stability.

$(D_{1})$ There is an increasing function ${{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in \mathcal{C}([0, T], {{\mathbb{R}}^{+}})$ and a number ${{\lambda }_{{{\mathbb{F}}_{\mathcal{Y}}}}} > 0$, so that

Lemma 5.11. Let ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ be the solution of (5.15), then

Proof. Let ${{Z}_{\mathcal{Y}}}\in \mathbb{X}$ be a solution of (5.15). Applying Remark 5.7, we obtain

Hence, a solution of (5.23) is provided below:

Since,

then, the inequality (5.22) is achieved. □

Theorem 5.12. Assume the conditions in Theorem 5.1 and Lemma 5.11 are satisfied, then the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is $\mathbb{UHR}$ stable.

Proof. Let ${{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}\in {{\mathbb{R}}^{+}}$ and ${{Z}_{\mathcal{Y}}}$ be the solution of (5.15). Suppose that $\mathcal{Y}\in \mathbb{X}$ is a unique solution of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2), then

Setting

this yields that $\left| {{Z}_{\mathcal{Y}}}(t)-\mathcal{Y}(t) \right|\, \, \le {{\Phi }_{{{\mathbb{F}}_{\mathcal{Y}}}}}{{\kappa }_{{{\mathbb{F}}_{\mathcal{Y}}}}}{{P}_{{{\mathbb{F}}_{\mathcal{Y}}}}}(t)$. Hence, by Definition 5.4, the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is $\mathbb{UHR}$ stable. □

Corollary 5.13. Setting ${{\kappa }_{{{\mathbb{G}}_{\mathcal{Y}}}}} = 1$ in Theorem 5.12, by Definition 5.5, the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) is $\mathbb{UHR}$ stable.

6.   Numerical algorithm

This section derives the numerical schemes for solving the approximated solution of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) by utilizing the decomposition technique for the $\mathbb{CPC}$ derivative operator. Next, we create an approximation design for the $\mathbb{CPC}$ derivative operator via $\alpha\in (0, 1]$ of a function $f(t)$. We will make a sequence of $N+1$ equations under $N+1$ conditions for the fractional Cauchy problem in the context of $\mathbb{CPC}$ derivative operator ^[47]. A sequence $(f_{N})$ of the solutions to such systems eventually leads to the solution of the obtained problem.

Theorem 6.1. Assume that $N$ is a positive number and $f\in A{{C}^{2}}([0, T], \mathbb{R})$. Let

${{\mathcal{V}}_{i}}:[0, T]\to \mathbb{R}$ be functions defined by

Then,

where $\underset{N\to \infty }{\mathop{\lim }}\, \mathcal{E}_{tr}(t) = 0$ for $t\in [0, T]$.

Proof. Using Definition 2.3 and ${\mathcal{K}_{1}}(\alpha) = {{Q}^{\alpha }}\left(1-\alpha \right)$, ${\mathcal{K}_{0}}(\alpha) = \alpha {{Q}^{1-\alpha }}{{C}^{2\alpha }}$, $\alpha \in (0, 1]$, we have

Let $u = {{Q}^{\alpha }}\left(1-\alpha \right)f(s)+\alpha {{Q}^{1-\alpha }}{{C}^{2\alpha }}{f}^{\prime}(s)$ and $dv = {{\left(t-s \right)}^{-\alpha }}ds$. Using the integrating by part technique, yields that

Applying the generalized binomial theorem, it follows that

Plugging (6.5) into (6.4), it follows that

where

Hence, by direct calculation, we get

Let $u = {{s}^{i}}$ and $dv = \frac{d}{ds}\left[ {{Q}^{\alpha }}\left(1-\alpha \right)f(s)+\alpha {{Q}^{1-\alpha }}{{C}^{2\alpha }}{f}^{\prime}(s) \right]ds$. Using the integrating by part technique, it follows that

Now, we show that $\mathcal{E}_{tr}(t) \to 0$ as $N \to \infty$ for $t \in [0, T]$, and to show this, we give an upper bound for the error term. Since $\Gamma (x+\alpha)\sim \Gamma (x){{x}^{\alpha}}$ and $s/t < 1$,

Then,

Substituting (6.7) into (6.6) with $M^{*}(t) = \max_{s \in[0, T]}\left\vert \frac{d}{ds}\left[ {{Q}^{\alpha }}\left(1-\alpha \right)f(s)+\alpha {{Q}^{1-\alpha }}{{C}^{2\alpha }}{f}^{\prime}(s) \right] \right\vert$ implies the following upper bound:

The righthand side of (6.8) tends to zero for all $t\in (0, T)$ as $N\to \infty$. □

To obtain the numerical approximation of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2), we apply Theorem 6.1. Then,

and

where $A_{N}$ and $B_{N, i}$ are given by (6.1) under the conditions

7.   Results and discussion

This section uses the numerical algorithm from the previous section to obtain the numerical solutions of the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2). The values of the basic parameters are listed as in Table 2 with the initial condition: $S_h(0) = 3.4 \times 10^8$, $E_h(0) = 1000$, $I_h(0) = 100$, $C_h(0) = 100$, $R_h(0) = 1 \times 10^6$, $V_h(0) = 1 \times 10^6$, $S_r(0) = 9 \times 10^4$, $E_r(0) = 100$, and $I_r(0) = 100$. From the given data, we obtain $\mathfrak{R}_0 = 1.5379 \times 10^4 > 1, \epsilon_1 = 1.42043, \epsilon_2 = 0.69421, \epsilon_3 = 0.13750, \epsilon_4 = 0.01003, \epsilon_5 = 0.00018, \epsilon_6 = 9.32369 \times 10^{-9}, \epsilon_7 = 1.03004 \times 10^{-13}, h_1 = 0.59741, h_2 = 0.00990, h_3 = 9.32362 \times 10^{-9}, g_1 = 0.11395, g_2 = 0.00018, g_3 = 1.03004 \times 10^{-13}, d_1 = 0.00896, d_2 = 9.32308 \times 10^{-9}, e_1 = 0.00018, e_2 = 1.03004 \times 10^{-13}, f_1 = 9.31797 \times 10^{-9}$. Consequently, these calculated values satisfy the conditions (i)–(vii) in Theorem $4.5$. Thus, the endemic equilibrium point is locally asymptotically stable.

Here, Figure 2 shows the numerical solutions for the classical $\mathbb{MPOX}$ model (3.1) of each group of human and rodent populations, while Figure 3 expresses the numerical solutions utilizing the $\mathbb{CPC}$ fractional operator for the $\mathbb{CPC}$-$\mathbb{MPOX}$ model (3.2) with varies $\alpha = 0.995$, $0.985$, $0.975$, $0.965$, $0.955$. As shown in both of the mentioned figures, the findings demonstrate that the fractional order model follows the same trend as the traditional model but is more flexible and has a higher degree of freedom. At larger fractional orders, increasing and decreasing behavior converges to the classical model more quickly than small fractional orders. Furthermore, from Figure 3 observation, a slight adjustment in the fractional order $\alpha$ was found to cause only a minor change in the behavior. This indicates that variations in the fractional order have a negligible impact on the stability of the disease dynamics across different groups within the proposed model. In addition, the simulation takes place over a $500$-day period. Under certain conditions, solution trajectories in all population groups reach the steady-state for all various values of $\alpha$ as time passes. Figure 2a and Figure 3a indicate the behavior of the susceptible humans. It is observed that the amount of them decreases rapidly in a short time in the beginning and tends to a steady-state as time tends to infinity. Figure 2b, Figure 3b, Figure 2c, Figure 3c, Figure 2d, Figure 3d, Figure 2i, and Figure 3i indicate the behavior of the exposed, infectious, clinically ill humans, and the infected rodents, respectively. We can see that the populations of each group are quite similar in order to gradually increase before reaching a steady state as time passes. Figure 2e and Figure 3e indicate the behavior of the recovered humans, and we extended the period to 5000 days for trend clearly observation. We found that the behavior of this group is increasing and tends to the steady state in the end. The behavior of the vaccinated humans is shown in Figure 2f and Figure 3f, demonstrating rapid population growth in a very short time and then reaching the equilibrium point. The behavior of the susceptible rodents drops initially before tending to a steady state, as seen in Figure 2g and Figure 3g. While Figure 2h and Figure 3h indicate the behavior of the exposed rodents, the population climbs to a peak and then falls until reaching a stable state.

Additionally, we place greater emphasis on examining the influence of the immunity-induced recovery rate $\omega$ on the dynamics of three human compartments regarding the $\mathbb{MPOX}$ infection, which are exposed, infectious, and clinically ill individuals. The parameter representing the recovery rate is varied with the following values: $0.4670$ (baseline from Table 2), $0.5137$ ($10\%$ increase), 0.5604 ($20\%$ increase), 0.6538 ($40\%$ increase), and $0.8406$ ($80\%$ increase) under vary $\alpha = 1$, $0.9$ and $0.8$ as seen in Figures 4–6, respectively. These graphs present two significant observations. First, an increase in $\omega$ levels corresponds to a clear reduction in the amount of infected individuals across all groups. Second, lower $\alpha$ values lead to a slightly faster decline in the amount of infected humans.

Similarly, investigating the impact of the contact rates between infected individuals or infected rodents on exposed, infected, and clinically ill humans is also highly significant. In this context, the parameter value for the contact rate of infected individuals with the susceptible population is utilized as follows: 0.5701 (baseline from Table 2), $0.62711$ ($10\%$ increase), $0.68412$ ($20\%$ increase), $0.79814$ ($40\%$ increase), and $0.91216$ ($60\%$ increase). The contact rate of infected rodents with the susceptible population is utilized as follows: $0.2508$ (baseline from Table 2), $0.27588$ ($10\%$ increase), $0.30096$ ($20\%$ increase), $0.35112$ ($40\%$ increase), and $0.40128$ ($60\%$ increase). Figures 7–9 demonstrate the impact of varying contact rate parameters $b_1$ and $b_2$ under different values of $\alpha$: $1$, $0.9$, and $0.8$, respectively. The graphs indicate that the contact rate substantially influences the population sizes of the respective groups. An increase in the contact rate leads to a corresponding increase in the amount of individuals in the groups mentioned. Besides, a reduction in the value of $\alpha$ leads to a slightly faster decline in the population sizes of all groups.

Therefore, based on all of the simulations above in various scenarios, we found a key finding that the fractional order or the memory index, the immunity-induced recovery rate, and the contact rates between infected individuals or infected rodents are factors that play an important role in determining whether monkeypox infection levels increase or decrease in the human population. These parameters can serve as control measures for monkeypox transmission.

8.   Conclusions

This study presented a $\mathbb{CPC}$ fractional derivative model for the spread of the monkeypox virus, accounting for interactions between humans and rodents. The model was developed to analyze system behavior and identify parameters that can help control the disease. For the theoretical part, we proved the positiveness and boundedness of solutions and examined the equilibrium points and the basic reproduction number. We investigated the local asymptotically stable steady state. We verified the qualitative results of the suggested model, including the existence and uniqueness results, using the Banach contraction mapping principle. Various Ulam's stability was demonstrated to ensure the existing solutions. The exactness of the theoretical guarantee is confirmed via the numerical simulations in all diagrams utilizing a decomposition formula for a constant proportional Caputo derivative. For the numerical and graphical parts, the graphical results highlight the key advantage of our proposed $\mathbb{CPC}$-$\mathbb{MPOX}$ model. It proves to be more realistic and practical than the traditional model, and its flexibility enhances precision, enabling us to achieve superior outcomes compared to the classical approach. Using the value of the parameters in Table 2, it was shown by numeric calculation that $\mathfrak{R}_0 > 1$ and all conditions of Routh-Hurwitz criteria are satisfied. Hence, we concluded that the endemic equilibrium point is locally asymptotically stable, which supports Theorem 4.5. Furthermore, the suggested model is shown by rising the value parameters of the recovery rate due to immunity and the contact rates of susceptible humans with infected humans and infected rodents at various levels and different fractional orders to analyze the effect in population dynamics influenced by the spread of $\mathbb{MPOX}$ via graphical simulation. The study demonstrated that the factors in our case can be incorporated into the control mechanisms of monkeypox transmission. These findings may offer valuable insights into the prevention and control of $\mathbb{MPOX}$ outbreaks in the future.

For future work, applying the $\mathbb{CPC}$ derivatives operator to study and analyze other epidemic models in real-world situations can be a proper alternative technique. On the other hand, other fractional operators, such as piece-wise and stochastic operators, can be considered for $\mathbb{MPOX}$ virus transmission to investigate real-world situations more realistically.

Author contributions

Jutarat Kongson: Conceptualization, Methodology, Writing-original draft, Writing-review & editing, Formal analysis, Funding acquisition; Chatthai Thaiprayoon: Conceptualization, Methodology, Writing–original draft, Writing–review & editing, Software, Supervision; Weerawat Sudsutad: Conceptualization, Methodology, Writing-original draft, Writing-review & editing, Formal analysis. All authors have read and agreed to the published version of the manuscript.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in creating this article.

Acknowledgments

J. Kongson (jutarat_@go.buu.ac.th) and C. Thaiprayoon (chatthai@go.buu.ac.th) would like to thank Burapha University for providing bench space and support.

Funding

This work was financially supported by (i) Burapha University (BUU), (ii) Thailand Science Research and Innovation (TSRI), and (iii) National Science Research and Innovation Fund (NSRF) (Fundamental Fund: Grant no. 75/2567).

Conflict of interest

The authors declare no conflict of interest.