Mathematical Modeling and backward bifurcation in monkeypox disease under real observed data

We propose a mathematical model to analyze the monkeypox disease in the context of the known cases of the USA epidemic. We formulate the model and obtain their essential properties. The equilibrium points are found and their stability is demonstrated. We prove that the model is locally asymptotical stable (LAS) at disease free equilibrium (DFE) under R0<1. The presence of an endemic equilibrium is demonstrated, and the phenomena of backward bifurcation is discovered in the monkeypox disease model. In the monkeypox infectious disease model, the parameters that lead to backward bifurcation are θr, τ1, and ξr. When R0>1, we determine the model’s global asymptotical stability (GAS). To parameterize the model using real data, we obtain the real value of the model parameters and compute R1=0.5905. Additionally, we do a sensitivity analysis on the parameters in R0. We conclude by presenting specific numerical findings.


Introduction
Animals, especially rodents and non-human primates, can pass the virus that causes monkeypox to people. Though often less severe, it resembles smallpox clinically. In Africa, particularly in the central and western areas, the illness has emerged as a serious public health concern since it mostly affects residents who live close to tropical rainforests. Monkeypox has emerged as the most significant Orthopoxvirus threat public health since the eradication of smallpox in 1980. Numerous animal species, including dormice, rope and tree squirrels, monkey species, Gambian pouched rats, and others have been shown to have the virus. To comprehend the virus's natural history and find the reservoirs where it lives, more study is needed. The first verified case of monkeypox was found in a nine-month-old baby in the Congo in 1970, and ever since then, the number of human cases has progressively climbed in a number of African nations. The severity of monkeypox epidemics varies, and the disease burden is poorly characterized. Nigeria reported around 200 verified cases and 500 suspected cases in 2017, with a 3 percent death rate. In order to stop the spread of the illness and lessen its negative effects on public health, it is essential to continue monitoring and researching monkeypox outbreaks [1][2][3].
reported the most instances in the US, and on August 4, 2022, the CDC declared a national public health emergency. The earliest state to report at least one case was Wyoming, and the first death in the US was on August 30, 2022 [7][8][9][10].
Mathematical modeling has been considered useful for understanding the physical problems, engineering and telecommunication problems [11][12][13]. Also, various kinds of applications of the mathematical models and their application can be seen in the recent work [14][15][16]. Application of the mathematical models in fluid related problems can be seen in [17,18]. The Various mathematical models for different disease dynamics have been presented in the literature, as seen in sources such as [19][20][21][22][23][24][25]. For instance, [19] presents a fractional order model for TB infection, while [20] studies the formulation of viral kinetics in HIV-1 infection. The authors proposed a mathematical study to analyze the HIV/AIDS disease using the actual data of Indonesia [26]. The malaria mathematical model under the seasonal factor of the mosquitoes has been analyzed in [27]. A mathematical study has been organized to understand the COVID-19 and the anxiety factors [28], modeling and simulation of the epidemic model using under mass action mechanism [29], and three species food chain dynamics [30]. COVID-19 epidemic [22], breast cancer [23], Lassa fever infection [24], and polio dynamics [25], the other infectious diseases [31][32][33] are also recent applications of mathematical modeling. Mathematical models for monkeypox disease have also been proposed in the literature. For example, in [34], the authors considered the monkeypox disease and study the result through a dynamical model. Monkeypox virus interaction between humans and animals is considered through mathematical modeling in [6], while [35] considered a mathematical study that explore the human transmission mechanism of the monkeypox virus. The treatment of monkeypox and the vaccine intervention is studied using compartmental mathematical models in [36], while [37] explored the monkeypox dynamics under an imperfect vaccine. In [38], the authors consider the two routes of disease transmission, such as human-to-human and rodents-to-human transmission of the monkeypox virus, and provide a detailed analysis of the disease without actual cases. The compartment mathematical model is presented to help understand the complex process of the disease. The writers of [39] employed mathematical modeling by segmenting the populations of people and animals into three subclasses each, and they produced comprehensive mathematical findings without using actual data. Our work aims to update the existing literature by considering a new mathematical model that studies the monkeypox disease with real data and considers both the possible routes of transmission such as humanto-human and rodents-to-human. To our knowledge, this work will be the first to study the mathematical modeling of monkeypox disease using real data. We hope that our work will be useful for disease elimination efforts in the USA.
We divide section-wise the work as follows: We give the formulation of the model and explain each route of transmission and their flow rate from each class in Section ''Model Construction''. The equilibrium points and analysis of the model are discussed in Section ''Equilibrium points and their analysis''. Estimation of the model parameters and the sensitivity analysis of the basic reproduction number is given in Section ''Parameters estimations''. We show the graphical results in Section ''Numerical method''. Summary regarding the work and recommendations is provided in Section ''Conclusion''.

Model construction
This work extend the work proposed in [39] by which the authors formulated a mathematical model for monkeypox disease by considering the population of humans into three classes, such as susceptible, infected, and recovered individuals while the population of animals is also distributed into three classes such as susceptible, infected, and recovered rodents. The authors analyzed the model and obtained their mathematical results without real statistics. Updating the work in [39], we consider studying the monkeypox virus under an extended mathematical model with real observed data from the US. We added the exposed and the quarantined populations and used the recent outbreak data in this new mathematical model proposed in this study. The model consists two different populations, that is the humans and the animal (rodents). The human population consists of five classes, that is, susceptible, ℎ ( ), exposed ℎ ( ), individuals infected with monkeypox virus ℎ ( ), individuals that are quarantined ℎ ( ), and people who are recovered from the disease are ℎ ( ). We divide the population of rodents into three subclasses; that is the susceptible, exposed, and infected rodents, respectively denoted by ( ), ( ), and ( ). Individuals are infected through animal to human, and humanto-human contact, and these two route of transmission are considered in the present work. So, the total population of humans and rodents are respectively ℎ ( ) = ℎ ( ) + ℎ ( ) + ℎ ( ) + ℎ ( ) + ℎ ( ) and ( ) = ( ) + ( ) + ( ). The dynamics of these two populations can be shown through the following nonlinear system of evolutionary differential equation: where ( ) = 1 + 2 ℎ ℎ , and ( ) = , are the force of infection, whereas the initial conditions are In the monkeypox system developed in (1), the parameter ℎ denotes the recruitment rate of the healthy population, and its natural mortality rate is given by ℎ . Two transmission routes possibly exist to infect humans: direct contact of healthy people with infected rodents, and human-to-human transmission that causes infection in human populations. The rate 1 is defined as the effective contact rate at which a human is infected by the virus per contact with an infected rodent. Through contact with the blood, body fluids, cutaneous or gastrointestinal wounds of diseased animals, animal-to-human (zoonotic) transfer can take place. Dormice, Gambian pouched rats, rope and tree squirrels, numerous monkey species, and other animals have all demonstrated signs of the monkeypox virus infection throughout Africa. Despite rodents being the most plausible candidate, the natural reservoir for the monkeypox virus has still not been found. Raw meat and products obtained from infected animals are major risk factors for infection. Living near or in woods can expose people to infected animals inadvertently or inadequately.
The rate 2 is defined as the effective contact rate among humans for transmission of the monkeypox virus. Close interaction with infected skin conditions, respiratory droplets, or newly contaminated objects can all result in the spread of the disease from person to person. Since droplet respiratory particle transmission frequently needs longterm face-to-face interaction, health professionals, family, and friends,  as well as other close connections of active patients are more at risk. Lastly, when a healthy person comes into having intimate touch with someone who has the monkeypox virus, they become infected and, after a period of 2-4 weeks when disease symptoms appear, they are identified as a member of the infected class ℎ at a rate of 1 . Some individuals in the exposed class are quarantined at a rate of 2 , depending on the severity of their symptoms, which may or may not indicate a severe infection. Recovery from infection occurs in the infected and quarantined classes at rates of and , respectively. Individuals in the infected class die due to the infection at a rate of 1 , while those in the quarantined class have a disease death rate of 2 . The population of animals (rodents) increases through the birth rate and decreases due to the natural death rate 3 . The parameter 3 is the contact rate between infected and uninfected rodents, and it represents the likelihood of a rodent contracting an infection for each encounter with an infected rodent. Exposed animals become infected at a rate given by . The interaction of humans and rodents population and their parameters transfer from each compartment is shown in Fig. 1.

Equilibrium points and their analysis
In this section, we will look at the basic features of the suggested model (1). We will first show the positivity and boundedness of the model, the equilibria, and its locally asymptotic results (LAS).

Model positivity and boundedness
Lemma 1. Assuming that the initial values of the variables in system (1) and (0) > 0, then for any > 0, any solution associated with the model will be non-negative for any non-negative > 0.

Proof.
To prove this result, we start with the first equation of system (1), Integrating Eq. (8), we obtain where ℎ0 represents the initial population size before the epidemic and is obviously positive. Therefore, ℎ ( ) > 0. By using the same method as above to prove the positivity of the remaining equations in system (1), Furthermore, the positive solution provided in Lemma 1 is obviously bounded, and it can be obtained by evaluating Eqs. (4) and (6) as approaches infinity. □

Equilibria and threshold
We will find the equilibria of the system (1) by equating the rate of change to zero at both the DFE and EE. The DFE associated to the model (1) will be denoted by 0 , which is given by The DFE 0 is useful in computing the threshold quantity, denoted by  0 . This threshold quantity is defined as the average number of monkeypox viruses inserted into a population with no monkeypox virus, producing further secondary monkeypox infections. The population with no monkeypox virus is known as the disease-free state, while the presence of monkeypox in a population is known as the endemic equilibrium. The monkeypox virus can be eliminated from the population if the threshold  0 < 1. However, if  0 > 1, the virus will become endemic. The threshold quantity  0 of the model (1) can be calculated using the well-known approach given in [40] by obtaining the required matrices: where 1 = ( ℎ + 1 + 2 ), 2 = ( ℎ + 1 + ), 3 = ( ℎ + 2 + ), and The spectral radius ( −1 ) provides the desired threshold quantity of the system (1), given by We provide the following theorem to show that the monkeypox model (1) is LAS.
Proof. The following Jacobian matrix is obtained at the monkeypoxfree equilibrium 0 , , − ℎ , − ℎ , − are the eigenvalues of the matrix ( 0 ) that contains negative real part. While the remaining four eigenvalues can be obtained using the characteristics equation: where 1 = 1 + 2 + 4 + , To ensure the stability of the system described by Eq. (11), it is important to demonstrate that the coefficients are positive, where = 2, 3, 4. While 1 is greater than zero, , where = 2, 3, 4, can only be positive if  0 < 1. Moreover, to confirm that the coefficients are positive, we need to satisfy the Routh-Hurwitz criteria, which will provide the eigenvalues with a negative real part. Therefore, for Eq. (11), we need to establish the following condition to ensure stability: ϝ = 1 2 3 > 2 3 + 2 1 4 . This condition is satisfied below, The condition ϝ > 0 guarantees that the fourth-order polynomial yields eigenvalues that have negative real parts. Therefore, if  0 < 1, the monkeypox system (1) is LAS at the DFE 0 .

Endemic equilibria and backward bifurcation
We obtain analytic expression for the endemic equilibria and to determine the possible existence of an unique EE. We denote the EE of system (1) by 1 , which is given by the values ( * ℎ , * ℎ , * ℎ , * ℎ , * ℎ , * , * , * ), and obtain the desired result in the following, Using (12) into the following, we get after some calculations the following result, * = * Further, substituting the result * into * , and after simplifying, we obtain, The coefficient 1 is positive and 3 is positive whenever  2 < 1. The coefficient 2 can determine the existence of the backward bifurcation in the monkeypox model. The following summary are provided: Backward bifurcation would occur for values of  2 such that  2 <  2 < 1. The bifurcation diagram is shown in Fig. 2, considering the parameter values shown in Table 1, except for = 0.0007028, 1 = 0.606, and = 0.000399. With these parameters and the rest from Table 1, we have  2 = 0.5404 < 1, indicating that the monkeypox model (1) exhibits backward bifurcation.
The occurrence of backward bifurcation in the monkeypox transmission model (1) has important epidemiological implications. It implies that the conventional criterion of  0 < 1 is no longer sufficient for disease eradication, although it is still necessary. In this case, disease eradication would be determined by the initial sizes of the subpopulation in the model (i.e., state variables). Therefore, the practicality of controlling monkeypox when  0 < 1 may depend on the starting sizes of the subpopulation. □

Global stability of the endemic equilibrium
We show the global asymptotical stability (GAS) of the model (1) in the absence of equation ℎ , as it is independent on the rest of the equations of the model (1). We first give endemic equilibria in the following at a steady-state at 1 , for system (1): The expression shown in (13) will be used later in the proof of the GAS of EE 1 .
Proof. Let us look at the Lyapunov function provided by, The time differentiation of Eq. (14) gives, The expressions on the right hand side of Eq. (15) are obtained using the equations from the model (1) as follows: ( ( ( ( Substituting Eqs. (16)- (22) into Eq. (15), and after simplifications, we get Here, using the property of arithmetic geometric mean, we have Thus, the largest invariant subset for which  ′ ≤ 0 is 1 . It follows from the results given in [43], that EE 1 is GAS if  0 > 1. □

Parameters estimations
This section describes the parameter estimation process for the monkeypox disease model (1) using real cases reported in the United States of America (USA) from May 2022 to March 15, 2023. The data was obtained from the CDC, as documented in [44]. The data consists of daily cases. For this data fitting process, we employ the nonlinear least square technique, which is a well-known curve fitting method for differential equation systems. To set the initial values of the variables given in system (1), we consider the population of the United States in 2022 to be ℎ (0) = 333, 287, 557, that representing the total human population. Assuming no disease is present, the initial healthy population is ℎ (0) = 333, 279, 356. The number of people who have been exposed to the disease is set according to the data fit, is ℎ (0) = 8000, and the first reported monkeypox case is ℎ (0) = 1. The number of individuals who are quarantined is ℎ ( ) = 200, while we assume that there is no recovered people from the disease is ℎ (0) = 0. Since there is no available data for the initial conditions of the animal (rodent) population, we assume that (0) = 9450. For the rodents, we set (0) = 9000, (0) = 400, and (0) = 50, based on the fitting of the model.
Among the model parameters, some can be calculated directly from the model equations, while the others are obtained during the model fitting experiment. For example, we can calculate the average lifespan of humans in the USA as 1∕(76.4 × 365) per day [41], and the natural birth rate using the formula ℎ = ℎ × ℎ (0), which yields an approximate value of ℎ ≈ 11951. Similarly, we assume the average lifespan of rodents to be = 1∕(5 × 365) per day, but this may vary depending on the species. The birth rate of rodents can be obtained using the relation = × (0). The remaining parameters are determined through fitting the monkeypox disease model (1) to experimental data, and their values are listed in Table 1, with time measured in days. Using these parameter values, we compute the threshold  1 ≈ 0.5905, which is consistent with the monkeypox outbreak in the USA [44], where the threshold remained below 1 throughout the outbreak. The graph in Fig. 3 shows the model fitting results for the number of cases.

Sensitivity analysis
This subsection focuses on the sensitivity analysis of the basic reproductive number. The factors that impact the dynamics of the Through this research, we can identify the most sensitive characteristics that significantly affect the basic reproduction and demonstrate how specific parameters can impact the growth or decline of the basic reproductive number. To obtain the sensitivity analysis of the basic reproductive number  0 , we use the procedure outlined in the literature [45].
Definition 1. The normalized forward sensitivity index of a variable under the parameter, can be defined through the following formula, Consider the above formula (25), we can have the analytical expression of  0 , It should be noted that  0 = max {  1 ,  2 } . For each parameters  0 and their numerical values are given in Table 1 are used to get the sensitivity indices. We shall have the sensitivity index of  1 with respect to the parameter 1 , given by Similarly, we follows the above formula and obtained the sensitivity indices,   Table 2 lists the sensitivity indices of the parameters in  1 . It reveals that 2 is the most sensitive parameter that can increase or decrease  1 . Additionally, 1 , 1 , 2 , , and ℎ are also sensitive parameters that can enhance  1 . Even small variations in the values of these parameters can significantly impact the basic reproduction number. Therefore, it is essential to focus on these transmission rates to reduce the future spread of the disease. The impact of these parameters is demonstrated graphically in the numerical simulation section.
Similarly, Table 3 shows the sensitivity indices of  2 , where and are the sensitive parameters.

Numerical method
This section presents the numerical approach used to solve the monkeypox model (1) numerically. In the past some practical problems under some real data and their analysis is given [46][47][48]. We will also provide a brief overview of the parameterized method's history and explain how to calculate the classical derivative when it is applicable. Notably, the considered approach come first from the work in [49] for classical derivative and later modified by the authors in [50]. The recent application of this proposed method can be seen in the recent work [51]. To obtain the necessary scheme using the method shown in [50], we consider the following general nonlinear equation: where (0) = 0 is the initial condition. Both sides integration of Eq. (28) gives the following, Take = +1 and = , and then after subtraction, we get, Within the interval , the approximation of function ( , ( )) is given by, After the replacement, we have The second term that have +1 and + 1 appeared in Eq. (32) were substituted, and their formulae were defined in [49] led to an implicit scheme that occasionally cause complications during simulation. So, we have where +1 = + ℎ. For the function +1 , we have F.M. Allehiany et al.

Fig. 5.
Numerical results of the human compartments with the variation in 1 . Subgraphs (a-c) explains the exposed, infected and quarantined population.

Numerical simulation
We conducted a numerical study of the monkeypox model using the parameter values listed in Table 1. Using the numerical method, we graphically present the results of the sensitive parameters that have a significant impact on the elimination of the disease. We first verify that the model is stable at its equilibrium points. As shown in Fig. 4(a-c), the model approaches its equilibrium points, and the solution behavior demonstrates a trend towards the disease-free equilibrium case. Fig. 4 depicts the rodent population and their dynamic behavior for the high endemic state. Fig. 5 illustrates the behavior of the human compartment for varying contact rates 1 . It shows that decreasing contact between humans and rodents can decrease the human infection rate. To prevent rodent infestation, it is crucial to follow guidelines and preventive measures, such as eliminating sources of food, water, and shelter for rats. This decreases the risk of human infection, and it is important to regularly and frequently dispose of trash both inside and outside the house. Additionally, to reduce the possibility of contracting germs or illnesses, it is recommended to thoroughly clean any areas where there are signs of rodent activity.
We present Fig. 6 which depicts the dynamics of the human compartments with varying human-to-human transmission coefficient 2 . Decreasing the transmission among humans, we can see a significant reduction in the number of infected cases. Early detection of cases and effective surveillance measures are critical for epidemic containment.  The most important risk factor for monkeypox virus infection during human outbreaks is intimate contact with infected patients. Healthcare people and those family members are considered to be at a higher risk of this virus. Therefore, health workers must strictly adhere to prescribed infection control measures when treating individuals with suspected or confirmed monkeypox virus infections or handling their specimens. Ideally, caregivers should be individuals who have received a smallpox vaccination.
Most human monkeypox cases have historically been caused by primary animal-to-human transmission. Thus, it is essential to keep wild animals, especially those that are sick or dead, and their byproducts such as meat, blood, and other bodily fluids, away from unprotected human contact. Any foods containing animal meat or other components must be thoroughly cooked before consumption. Infected individuals with monkeypox virus should be quarantined to reduce further spread of the infection among humans. Kissing, massaging, and prolonged face-to-face interaction are also potential transmission routes. Personto-person transmission of monkeypox is now occurring in Canada. The majority of monkeypox cases in Canada have been reported in men with multiple sexual partners, which is consistent with global trends. However, it is important to note that the risk of exposure to the virus is not limited to any particular group and having multiple sexual partners can increase the risk of infection [52].
The impact of the parameter 1 on the human compartment can be seen in Fig. 7. As the parameter 1 decreases, the number of infected humans also decreases. To control the further spread of the monkeypox virus in the human population, it is important to isolate infected individuals and minimize their close contact with healthy people. Fig. 8 shows the parameter coefficient , which is responsible for humans being infected by animals (rodents). A decrease in the variation of leads to a reduction in human infections. To prevent infection spread to humans, it is recommended to limit close contact with animals carrying the virus. Additionally, it is advised to avoid direct contact with animals that carry infections transmittable to humans, such as contagious rash or contaminated body fluids (respiratory secretions, saliva, scabs, crusts, or sores). Our results in Figure align with these recommendations. Fig. 9 shows the population compartments of rodents, with variations in the contact rate parameter . By reducing the contact rate of , the number of exposed and infected rodents decreases. By decreasing the population of rodents, we can decrease the risk of infection in the human population. Necessary steps can be taken to achieve this, such as cleaning areas where rodents are active, removing water and food sources, and eliminating places where rodents can shelter. Additionally, regular garbage disposal from homes can also help to control rodent populations.

Conclusion
In this study, we investigated the dynamics of the monkeypox outbreak in the USA using real statistical data. We formulated a model and discussed its basic results in detail, including the calculation of the basic Fig. 8. Numerical results of the human compartments with the variation in . Subgraphs (a-c) explains respectively, the exposed, infected and the quarantined population. Fig. 9. Graphical result of the rodents compartments with the variation in . Subgraphs (a) and (b) describe the exposed, and the infected rodent population.
reproduction number, which consists of the human and rodent basic reproduction numbers. We also obtained local asymptotic results for the monkeypox model, showing that the model is locally asymptotically stable if  0 < 1. Additionally, we analyzed the endemic equilibria and found the existence of a backward bifurcation phenomenon. We established the result for the backward bifurcation under certain parameters and found that the model is not globally asymptotically stable for the disease-free case due to this bifurcation. However, we were able to establish global asymptotic stability at the endemic equilibria, and found that it is globally asymptotically stable at 1 .
We utilized real experimental data obtained from CDC [44] to fit our model using a nonlinear least square curve fit technique. Some of the parameters were obtained from literature and existing equations of the model, while the rest were fitted to the data. The basic reproduction number associated with humans was accurately estimated to be  1 = 0.5905 using our method. We performed a sensitivity analysis of the model and identified the most sensitive parameters that can help in disease elimination. The numerical scheme was also solved, and the results showed that reducing contact between rodents and humans by eliminating their food source, water, and shelter, regularly disposing of trash, and decreasing human-to-human transmission can all lead to a decrease in future cases. We recommended avoiding contact with infected animals and humans, washing hands with soap and water after any contact, and separating infected patients from healthy individuals. We conclude that our study can help in managing and controlling the spread of monkeypox.

Declaration of competing interest
The authors declare that they have no potential conflict of interest regarding the publishing of this paper.

Data availability
Data will be made available on request.