Mathematical assessment of monkeypox disease with the impact of vaccination using a fractional epidemiological modeling approach

This present paper aims to examine various epidemiological aspects of the monkeypox viral infection using a fractional-order mathematical model. Initially, the model is formulated using integer-order nonlinear differential equations. The imperfect vaccination is considered for human population in the model formulation. The proposed model is then reformulated using a fractional order derivative with power law to gain a deeper understanding of disease dynamics. The values of the model parameters are determined from the cumulative reported monkeypox cases in the United States during the period from May 10th to October 10th, 2022. Besides this, some of the demographic parameters are evaluated from the population of the literature. We establish sufficient conditions to ensure the existence and uniqueness of the model’s solution in the fractional case. Furthermore, the stability of the endemic equilibrium of the fractional monkeypox model is presented. The Lyapunov function approach is used to demonstrate the global stability of the model equilibria. Moreover, the fractional order model is numerically solved using an efficient numerical technique known as the fractional Adams-Bashforth-Moulton method. The numerical simulations are conducted using estimated parameters, considering various values of the fractional order of the Caputo derivative. The finding of this study reveals the impact of various model parameters and fractional order values on the dynamics and control of monkeypox.


Basic definitions
Fractional differential operators in mathematical modeling are well-known and it has been successful in the fields like science, engineering, and epidemiology.We give a few fundamental definitions along with some proprieties that be will used in the rest of the study 30,31 .

Definition 1
The left and right Caputo fractional derivative operator of the function f ∈ L ∞ (R) ∩ C(R) is defined as follows Definition 2 For x ∈ R , the generalized Mittag-Leffler function E α,β (x) is defined by and satisfies the following property given in (1) , α, β > 0,

Formulation of the integer order monkeypox model
The total human population N h is divided into six different compartments i.e S h susceptible humans, E h exposed humans, I h infected humans, C h clinically ill humans, V h vaccinated humans and R h recovered individuals.Simi- larly, N r represents the overall non-human population which is further divided into four different subgroups i.e. S r susceptible, E r exposed, I r infected and R r recovered animals respectively.The susceptible human class is initiated through recruited rate φ h .This population is declined by the natural mortality rate µ h and by the population gaining infection after interacting with infected humans (or animals) at the contact rate β 1 or β 2 respectively.The susceptible human group is also decreased by vaccinating at the rate α h .This class is further increased by the vaccinated population after waning the induced immunity at the rate η .Thus, the dynamics of the susceptible human class are modeled via the following equation The class of exposed individuals is increased by joining the newly infected individuals with the force of infection . Moreover, this population class is declined by the transmission rate β and the mortality rate µ h .Thus, we obtained the following differential equation describing the dynamics of exposed human class The number of infected people rises at a rate of β after the transition from being exposed class.This class is decreased by the natural death rate µ h , the disease-induced death rate δ 1 , the recovered population at the ω 1 and by the people who became clinically ill and join next class at the rate γ .Thus, the resulting dynamics can be modeled via the following equation The people in I h class become clinically (or critically) ill at the rate γ and join C h class.These individuals are reduced by natural mortality rate µ h and disease-induced death rate δ 2 .Further, it is decreased by joining the recovered class at the rate ρ .Hence, we obtained the following equation describing the dynamics of C h population The susceptible individuals are vaccinated at the rate α h .The vaccinated individuals are decreased due to natural mortality rate µ h and due to the waning of induced immunity at a rate η .Thus, we have Finally, the class of recovered humans is formed due to the transition of individuals from I h and C h at the recovery rates ρ and ω 1 respectively.Further, it is decreased due to the natural death rate µ h .Thus, we obtain The recruitment rate of susceptible animals is φ r and the natural death rate in all animal populations is denoted by µ r .The force of infection in this case is given by (3) E α,β (x) = x E α,α+β (x) + 1 Ŵ(β) .
where β 3 is the effective contact rate causing the disease transmission rate among animals.The flow from exposed to the infected animal compartment is denoted by ε .The death rate induced due to infection in infected animals is δ 3 .The parameter ω 2 shows the recovery rate of infected animals.Thus, the dynamics of monkeypox in animals are given by the following system Fractional monkeypox model in the Caputo sense.This section presents the extension of the integer order monkeypox model to the fractional case.The well-known Caputo derivative having order 0 < ϑ ≤ 1 is utilized to formulate the generalized fractional epidemic model.Fractional epidemic models used fractional derivatives to capture more complex and non-local effects in the disease transmission dynamics.Such epidemic models have a greater degree of freedom and offer insights into the behavior and control of infectious diseases for a particular set of data 22 .The proposed model for transmission of monkeypox model is described by the following system where,

Subject to non-negative initial conditions
Parameter estimation and data fitting.The objective of the present section is to accomplish the parameter estimation of the model.The estimation procedure is performed in two ways.Most of the model parameters are estimated from the real statistics of infected cases reported from 10 May to the end of October 2022 in the recent outbreak in the USA 1,3 .While some of the demographic parameters (the natural death and birth rate) are estimated from the USA population 32 .The estimation of parameters from the real disease data is significant for carrying out the numerical simulations more realistically.The estimation process of the proposed model (6) is conducted using the well-known nonlinear least square technique.This procedure involves the following main steps: 1. Define an objective function for quantifying the difference between the predicted and the actual data.This function is mathematically expressed as the sum of squared residuals (SSR) between the model predictions and the corresponding actual data points.2. Define or set the model parameters' initial values that need to be estimated.These values can be based on prior knowledge or initial guesses.3. Utilizing the initial parameter values set in step 2, the epidemic model is simulated in order to generate model predictions.4. Using a 'lsqcurvefit' optimization technique, for the minimization of objective function.
( www.nature.com/scientificreports/ 5. Set a convergence criterion to stop the iterative estimation scheme.This criterion is based on a maximum number of iterations, a minimal change in parameter estimates, or attaining a predefined threshold for the corresponding objective function.6.If the estimated parameter values do not adequately fit the real data curve or do not meet the convergence criterion, reset the initial parameter estimates and re-execute steps 3-5 until a reasonable agreement between the model simulation and the actual data is obtained.
By iteratively minimizing the objective function via a standard nonlinear least squares technique, the parameter estimates are refined, which leads to improved agreement between the model-predicted and observed data.The fitted and estimated values of the monkeypox epidemic model's parameters are listed in Table 1, while Fig. 1 illustrates the model's good fit to the observed data.

Basic analysis of the fractional monkeypox model
In this section, we address some of the basic and necessary mathematical analysis of the fractional monkeypox epidemic model (6).We proceed as follows  µ h and N r 0 (t) ≤ φ r µ r then the feasible domain for the given system is given by ⌣ R r 0 be positive, then all solution are will be positive for t > 0 .Considering the first equation in (6) we have then Since φ h + η V h ≥ 0 therefore, it implies By using the Laplace transform we have Applying inverse Laplace transform we have Since, quantities on right hand side are positive so we have S h ≥ 0 for t ≥ 0 .In similar approach, we have E h , I h , C h , V h , R h , S r , E r , I r and R r ≥ 0 ∀t > 0 corresponding to any non-negative initial values.Thus, the solu- tion remain in R 6 + × R 4 + for all t > 0 with non-negative initial values.Next to we prove the boundedness of the fractional model solution.The total human individuals is given by such that Taking Laplace transform on both sides we have By applying inverse Laplace transform we have In similar way, we can prove that N r (t) is bounded.Therefore, this established the notion of the set as required.So, we conclude that the region is epidemiological feasible and well posed in .Let T be a positive real number and consider J = [0, T] .We denote the set of all continuous function defined on J by C 0 e (J) with norm as �X� = sup {|X(t)|, t ∈ J}.Consider the system (6), (7) as following initial value problem (I.V.P.) where X(t) = (S h (t), E h (t), I h (t), C h (t), V h (t), R h (t), S r (t), E r (t), I r (t), R r (t)) represents the compartments and F is continuous function defined as follows Using property (2.1) in proposition we have The Picard iterations lead as follows Vol:.( 1234567890 9) can be written as follows

Lemma 1
The vector F(t, X(t)) described in (10) fulfills the well-known Lipschitz condition in the variable X on a set Vol.:(0123456789) Combining all these we have Lemma 2 Suppose we have (12) then the I.V.P. (6), ( 7) have unique solution X(t) ∈ C 0 e (J).
Proof For the required result, the fixed point theory and Picard-Lindelöf is used.Solution of the (6), ( 7) is considered as X(t) = W(X(t)), where W is the Picard operator defined as W : C 0 e J, R 10 Further, it leads to If ψ Ŵ(ϑ+1) W < 1 , then W gives a contraction and therefore, the problem (6), ( 7) has a unique solution.
Model equilibria and the basic reproduction number R 0 .The model disease free equilibrium (DFE) is given by Vol:.( 1234567890)  33 we compute the basic reproduction number R 0 .Let x = (E h , I h , C h , V h , E r , I r ) then we have where Then we considered the Jacobian of the Linearized system at disease free state which is given by Where, F and V are 6 × 6 matrices calculated as . Also, F and V contains the linear and non-linear terms of infected compartments.So, the reproduction number is calculated as Existence of endemic equilibrium point.The endemic equilibrium (EE) of the fractional monkeypox model is given by and is defined as follows ( 13) Moreover, by using ( 16) in (17) we have Global stability at EE.In the system (6 , therefore, we have h 1 = α 3 I h + α 4 I r and r 1 = α 5 I r .Therefore, the system (6) becomes Results for system (19) at steady states are calculated as follows Proof Consider the following nonlinear Lyapunov function given by (20) The Caputo fractional derivative of (20) implies ( 16) * r = Vol.:(0123456789) r , E * r , I * r .Thus, by LaSalle's invariance principle the EE is GAS in whenever, R 0 > 1.

Interpretations of R 0 versus model parameters
In this section, we analyze the impact of some model parameters on R 0 .The basic reproduction number R 0 is a key biological parameter that plays an important role in describing the disease dynamics and its possible control.In these graphical results, we consider the disease transmission rates β 1 , β 2 , vaccination rate α h , recovery rates ρ, ω 1 and vaccine waning rate η .The respective graphical results are shown in Figs.2a, 3, 4, 5 and 6a, while the corresponding contour plots in each case are shown in Figs.2b, 3, 4, 5 and 6b. Figure 2 shows the combine impact of β 1 and α h on R 0 .It shows that the R 0 decreases with a decrease in transmission coefficient β 1 and an increase in vaccination rate α h .Further, impact of β 2 along with α h upon R 0 is demonstrated in the Fig. 3a with the corresponding counter plot 3b.It can be observed that with the increase in vaccination rate and reduction in contact rate β the basic reproduction reduces to a value less than 1.Similarly, the impact of ρ and β 2 , ω 1 and β 2 , η and α h are presented in Figs.4a, 5a, and 6a respectively.We can observe that with the decrease in disease transmission rate β 2 and enhancement in vaccination rate α h and recovery rate ω the basic reproduction decreases significantly to a value less than unity.

Numerical results of the fractional model
In this section, we first investigate the numerical solution of the Caputo fractional monkeypox model.The generalized fractional Adams-Bash forth-Moulton approach 34 , an effective iterative method is used to solve the model numerically.The simulation for different model parameters and different values of ϑ are carried out using the generated numerical technique and the values for model parameters are listed in Table 1.The subsequent subsection contains the iterative solution.www.nature.com/scientificreports/Numerical technique.A concise numerical method for the approximate solution of the monkeypox transmission model in the Caputo sense is shown in this subsection.For this purpose, the fractional Adams-Bashforth-Moulton is used.The system (6) can be recomposed in the subsequent problem to get the desired scheme: where + , and M t, v(t) shows a continuous real valued function.Using the concept of integral in Caputo case, the aforementioned Eq. ( 22) is considered as follows: In order to use the procedure in 34 , a uniform grid on [0, T ] with step size h = T N , N ∈ N , where t u = uh, u = 0, 1 . . ., N is considered.Thus, the proposed model in (6) can be treated as:    Simulation of the fractional monkeypox model.In this section, we applied the iterative technique developed in the preceding section in order to simulate the monkeypox Caputo epidemic model (6).The simulation results are illustrated by taking the baseline values of estimated parameters listed in Table 1.In the simulation results, we mainly focus on the dynamics of infected human and animal populations.Initially, in the simulation, we analyze the dynamics of model state variables for various values of ϑ ∈ (0, 1] to study the impact of memory index on disease incidence.Further, the influence of variation in the recovery rate ω 1 on the behavior of the infected human population is predicted for many values of fractional order ϑ .Moreover, we considered different scenarios by decreasing the effective contact rates among susceptible humans and infected animals (i.e., β 1 , β 2 ) and the parameter β 3 showing the effective contact rate between susceptible and infected animals.In addition, these simulations are performed for various non-integer values of the parameter ϑ .The detailed discussion of the aforementioned simulation is described in the following subparts.Impact of recovery rate ω 1 and memory index ( ϑ).In this subpart, we analyzed the impact of variation in recov- ery rate ω 1 on the dynamics of cumulative cases of exposed E h , infected I h , and clinically infected C h human under different values of memory index ( ϑ ).The parameter showing recovery rate is perturbed with different rates (with baseline and with 10% , 30% , and 50% increase) depending on the treatment strength provided to the infected human population.The resulting simulation is accomplished in Figs. 9, 10 and 11 with subplots (a-c) for three values of fractional order ϑ .It is observed that with an increase in parameter ω 1 , the infected population in all aforementioned classes reduces showing the impact of treatment on the disease incidence.Additionally, it can observed that in the case of smaller values of ϑ , the decline in the infected population is slightly faster.www.nature.com/scientificreports/Impact of disease transmission rates and memory index ( ϑ).This section describes the influence of the infection transmission rates β 1 , β 2 , and β 3 on the dynamics of infected human and animal population classes.The graphi- cal results are initially performed for the baseline values of the aforementioned parameters given in Table 1.Subsequently, the parameter values are decreased with 10%, 20%, 40% , and 50% to the estimated values.Further- more, the simulation results are presented for three values of fractional order ϑ .These resulting simulations for the human population are shown in Figs. 12, 13 and 14 with respective subplots from (a-c) while the dynamics of the infected animal population are demonstrated in Figs. 12, 13 and 14d, e.In all plots, it is evident that the cumulative infected population significantly decreases as the infection transmission rates are reduced.Moreover, it can be seen from Figs. 13

Conclusion
In this manuscript, we analyzed the dynamics of a novel monkeypox infection with the case study of the recently reported outbreak.The model is first formulated using an integer-order nonlinear system of ten differential equations.The human population was divided into six distinct subgroups, while the animal population was divided into four classes.To estimate the model parameters, we fit the model to the actual cases of the 2022 monkeypox outbreak in the USA.Moreover, keeping the importance of the fractional modeling approach, the integer case model is extended to fractional order via the well-known Caputo operator.In the initial stage, we presented a comprehensive theoretical analysis of the fractional monkeypox model, including the existence and uniqueness of the solutions.The existence and stability analysis of the model equilibria are provided.Furthermore, the fractional model is solved numerically using the fractional Adams-Bashforth-Moulton approach.The simulation results illustrate the impact of disease incidence graphically, considering both the baseline values of the parameters and different values of the fractional order of the Caputo operator.Moreover, the model is simulated by increasing the parameter ω 1 (the recovery rate) and decreasing the disease transmission coefficients β 1 , β 2 , β 3 with different rates to their baseline values.The implementation of a real-data set of monkeypox infections makes this study more visible and important within the existing literature.with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates (e) Baseline values of β 1 , β 2 , β 3 with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates Baseline values of β 1 , β 2 , β 3 with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates Baseline values of β 1 , β 2 , β 3 with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates (e) Baseline values of β 1 , β 2 , β 3 with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates (a)  with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates (e)

Figure 1 .
Figure 1.Model predicted curve (solid line) to the observed data for the case when ϑ = 1.

Figure 2 .
Figure2.The subplot (a) illustrate the impact of β 1 (disease transmission coefficient relative to I h ) and α h (vaccination rate ) on R 0 , where (b) shows the respective contour plot.

2 Figure 3 .
Figure 3.The subplot (a) analyze the impact of β 2 (disease transmission coefficient relative to I r ) and α h (vaccination rate ) on R 0 , where (b) shows the respective contour plot.

Figure 4 .Figure 5 .Figure 6 .
Figure 4.The subplot (a) demonstrate the impact of β 2 (disease transmission coefficient relative to I r ) and ρ (recovery rate of C h ) on R 0 , where (b) shows the respective contour plot.

Figure 7 .
Figure 7. Simulations of only human classes in fractional monkeypox epidemic model (6) for various values of fractional order ϑ.
https://doi.org/10.1038/s41598-023-40745-xwww.nature.com/scientificreports/Impact of memory index ( ϑ).In the first case of simulation, we illustrate the impact of only fractional order ϑ of the Caputo operator on the model dynamics.The graphical results are shown in Figs.7a-e and 8a-c.The dynamics of only human population classes are shown in Fig.7while Fig.8illustrates the dynamical aspects of only animal subpopulation.The simulation in both cases is performed for four values of ϑ = 1, 0.95, 0.90, 0.85 .The time length considered in the simulation is taken to 500 days.It is observed that under the specific conditions, the solution trajectories in both population classes are converged to the steady state for all values of fractional order ϑ.

Figure 8 .
Figure 8. Simulations of only non-human classes in fractional monkeypox epidemic model (6) for various values of fractional order ϑ.

Baseline values of β 1 ,
β 2 , β 3 with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates

Table 1 .
Fitted and estimated values of the model parameters.
Baseline values of β 1 , β 2 , β 3 with 10% increase in contact rates with 20% increase in contact rates with 40% increase in contact rates with 50% increase in contact rates(b)