Unfolding the Transmission Dynamics of Monkeypox Virus: An Epidemiological Modelling Analysis

Abstract

Monkeypox (mpox) is a zoonotic viral disease that has caused recurring outbreaks in West Africa. The current global mpox virus (mpoxv) epidemic in endemic and non-endemic areas has seriously threatened public health. In this study, we design an SEIR-based deterministic model that considers prodromal stage, differential infectivity, and hospitalisation to investigate the transmission behaviour of mpoxv, which could help enhance control interventions. The model is theoretically analyzed by computing essential epidemiological quantities/dynamics, such as the basic reproduction number, which estimates the number of secondary infections caused by a typical primary case in an entirely susceptible community. Stability of the model’s equilibrium states is examined to evaluate the transmission potential of the mpoxv. Furthermore, partial rank correlation coefficient was adopted for sensitivity analysis to determine the top-rank model’s parameters for controlling the spread of mpoxv. Moreover, numerical simulations and model predictions are performed and are used to evaluate the influence of some crucial model parameters that help in strengthening the prevention and control of mpoxv infection.

1. Introduction

Public health and the economy continue to face a serious worldwide threat from emerging diseases [1]. A contagious disease known as monkeypox (mpox) is caused by the monkeypox virus (mpoxv) [2]. It is a viral infection that can spread from animals to humans (also known as viral zoonosis) [2]. Other zoonotic diseases causing serious concerns to public health include Lassa fever [3] and query fever [4]. In 1970, the first human cases were found in the Democratic Republic of the Congo, Liberia, and Sierra Leone [2]. The unexpected large outbreaks of mpoxv in both endemic and non-endemic nations, as well as the prevalence of sexual transmission routes, have recently raised severe concerns globally [2,5,6].

Even though mpox is a zoonotic viral infection found primarily in the tropical rainforests of Central and West Africa, it is occasionally exported to non-endemic areas. Contact with an infected individual, animal, or material can transmit mpoxv to humans [2]. In Nigeria, the primary mpoxv transmission routes are zoonotic and human-to-human. Potential symptoms comprise swollen lymph nodes, fever, lymphadenopathy, headache, and a rash that develops blisters and then encrusts the entire body [2]. The incubation period lasts 5 to 21 days, and 9.8 (95% CI: 5.9–21.4) days is the estimated mean serial interval (i.e., the interval between an infected person’s symptoms and those of a subsequently infected person) [2,7,8,9].

Over 83,000 people were affected by the virus between 10 May and 29 December 2022 in 110 locations with >98% of those locations with no previous infection [10,11]. Additionally, a significant proportion has been linked to males, especially those who identify as gay or bisexual [2,6,12]. There is no evidence that the virus might spread via sexual contact. However, transmission via close contact, particularly sexual contact, has been reported [13,14].

Smallpox vaccines are also effective against mpoxv infection [2]. However, the vaccine’s inadequacy is a serious problem and might make the situation more contagious. Thus, it is essential to put greater effort into non-pharmaceutical intervention control strategies for prompt mpoxv control. On 21 July 2022, the World Health Organization (WHO) declared the virus a Public Health Emergency of International Concern (PHEIC) due to the recent rise in mpoxv in non-endemic countries linked to an epidemic in the United Kingdom [2,15]. The basic reproduction number (${\mathcal{R}}_0$, or the typical number of secondary mpoxv cases produced by a single case throughout the infectious period in a community without immunity or any intervention) of the mpoxv was predicted by early models to be >1 in populations of men who have sex with men (MSM) and ≤1 in other settings [2,7,16,17]. As a result, investigating the effectiveness of non-pharmaceutical intervention options to aid in transmission halt is critical.

To comprehend mpoxv transmission, considerable epidemiological models have been developed and employed, for instance, [5,12,18,19,20,21,22] and the references therein. In order to study the mpoxv transmission patterns and control strategy, Yuan et al. [5] developed a model taking into account the zoonotic transmission route. They discovered that effectively quarantining contacts with mpoxv in high-risk populations by post-exposure vaccination greatly slows the transmission. Additionally, in their similar research, Yuan et al. [18] proposed a new model to scrutinise the efficiency of vaccination and other control means at large gatherings. Their research demonstrates that, in the absence of public health management measures, the risk of a mpoxv outbreak persists at high levels during mass gathering occasions. They recommended that stringent contact tracing, testing, case identification, and isolation, as well as awareness campaigns, are crucial tactics to stop the spread of the mpoxv outbreaks. Endo et al. [12] developed a transmission model and applied it to empirical sexual partnership data to show that the heavy-tailed sexual partnership distribution could explain the ongoing spread of mpoxv among MSM population. They estimated ${\mathcal{R}}_0$ to be >1 in the MSM context, complicating outbreak containment. They proposed that enlightenment programs could facilitate the prevention and early detection of mpoxv among MSM communities.

Although clinically less severe than smallpox, the symptoms of the viral zoonosis mpox are comparable to those experienced by individuals with smallpox in the past. Mpox has taken over as the most significant orthopoxvirus for public health since smallpox was eradicated in 1980 and smallpox vaccinations were subsequently discontinued. Primarily affecting Central and West Africa, the virus has been spreading into cities and is frequently seen close to tropical rainforests. Numerous rodent species and non-human primates serve as hosts for animals. As a result, it is critical to investigate more information and validate facts about the so-called ailment mpox, which has inspired us to conduct this study, thereby establishing more facts about the ailment.

This study’s objective is to develop a novel dynamic model to simulate the spread of mpoxv as an emerging zoonosis. The model incorporates human-to-human and animal-to-human transmission to evaluate the impact of high- and low-risk populations, various infection stages, and isolation on overall transmission. The model is utilised to determine the most effective tactics for containing the mpoxv outbreak in Nigeria. The results would help with a more thorough investigation of the mechanisms of disease transmission within and between the animal reservoir and humans, which would help calculate the risk of outbreaks via various transmission pathways.

Our model distinguishes different phases of the infection, including asymptomatically infected (prodromal) and symptomatically infected (mild and severe) stages [2,11,23]. The mild stage is characterised by common signs and symptoms that can be treated if caught early, including fever, exhaustion, body pains, and very light rashes in some areas of the body. The severity of the infection, however, is dependent on the infected individual’s health and the method of exposure, with complex symptoms including extensive rashes that affect all parts of the body, including the face, mouth, eyes, soles of the feet, throat, and genital and anal areas.

This paper is summarised as follows. Following this introductory Section 1, the mpoxv model is formulated and theoretically analyzed in Section 2 and Section 3, respectively. Numerical simulations and sensitivity analysis are conducted in Section 4. We end the paper with concluding remarks in Section 5.

2. Methods

2.1. Mpox Outbreak Data

More often, we retrieved the epidemiological cases data for the mpoxv infection on a weekly basis for Nigeria issued by the Nigeria Center for Disease Control (NCDC) for the situation from January 1 through August 27 of 2022 [23] after laboratory confirmation and case characterisation of mpoxv from the NCDC. From the data, we derived the weekly cumulative incidence and examined the Nigerian mpoxv cases scenarios.

2.2. Mpoxv Model Description

The model developed in this study will describe epidemiologically the dynamics of mpoxv transmission by examining the mpoxv outbreak’s propagation patterns and prevention measures using a traditional SEIR-typed model. The proposed model accounts for high- and low-risk susceptible populations to investigate the transmission behaviour of the mpoxv. We specifically distinguished between the various infectious stages to evaluate the overall dynamics of mpoxv epidemic propagation. Hospitalisation is also taken into account by the model as an intervention strategy, which is vital for the prevention and management of mpoxv. Isolation is usually enforced to help curtail widespread disease outbreaks, where infected people are admitted to hospitals or isolation facilities [5,24,25]. The parameter $\rho$ designates the fraction of newly recruited people moving to a high-risk susceptible population, $S_h$, while $1-\rho$ represents the fraction of people moving to a low-risk susceptible population, $S_l$. The parameters ${\pi }_h$ and ${\pi }_r$ represent recruitment rates for humans and rodents, respectively. All other parameters are epidemiologically explained in

Table 1. Enunciation of variables and parameters of governing system (1).

Variable	Description
$N_h$	Overall human population
$S_h$	Population of susceptible humans at high-risk
$S_l$	Population of susceptible humans at low-risk
$E_h$	Exposed humans
P	Populations of infectious humans in the prodromal stage
$I_1$	Symptomatically mild infectious humans
$I_2$	Symptomatically severe infectious humans
H	Hospitalised humans
$R_h$	Recovered humans
$N_r$	Total rodent population
$S_r$	Susceptible rodents
$E_r$	Exposed rodents
$I_r$	Infectious rodents
$R_r$	Recovered rodents
Parameter
${\beta }_i(i=hh,rh,rr)$	Transmission rates/transmission probabilities
$\alpha ,\eta$	Modification parameters
${\pi }_h\left({\pi }_r\right)$	Human recruitment rate (rodents)
$\rho$	Fraction of newly recruited humans moving to $S_h$
v	Rate of reduction in infectiousness from $S_l$
${\sigma }_h\left({\sigma }_r\right)$	Progression rates
$\omega$	Rate at which mpoxv progress from P to $I_1$
$\theta$	Progression rate from mild to severe disease
$k_1,k_2$	Hospitalisation rates
${\delta }_m(m=i,h,r)$	mpoxv-induced death rates
${\tau }_j(j=1,2,3,r)$	Recovery rates
${\mu }_h\left({\mu }_r\right)$	Human (rodents) natural death rate

. We separated the total population of humans at time t, denoted by $N_h\left(t\right)$, into high-risk susceptible, $S_h\left(t\right)$; low-risk susceptible, $S_l\left(t\right)$; exposed, $E_h\left(t\right)$; prodromal (asymptomatically infected), $P\left(t\right)$; mild infectious, $I_1\left(t\right)$; severe infectious, $I_2\left(t\right)$; hospitalisation, $H\left(t\right)$; and recovered, $R_h\left(t\right)$. Hence,

$$N_h\left(t\right)=S_h+S_l+E_h+P+I_1+I_2+H+R_h.$$

Similarly, we separated the overall population of rodents at time t, indicated by $N_r\left(t\right)$, into $S_r\left(t\right)$, $E_r\left(t\right)$, $I_r\left(t\right)$, and $R_r\left(t\right)$, hence

$$N_r\left(t\right)=S_r+E_r+I_r+R_r.$$

In light of the aforementioned details, we illustrated the mpoxv model in Figure 1. The following systems of non-linear ordinary differential equations are satisfied by the state variables and model parameters listed in Table 1:   

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_l}{\mathrm{d}t}& ={\pi }_h(1-\rho )-v{\lambda }_h{S}_l-{\mu }_h{S}_l,\hfill \\ \hfill \frac{\mathrm{d}{S}_h}{\mathrm{d}t}& ={\pi }_h\rho -{\lambda }_h{S}_h-{\mu }_h{S}_h,\hfill \\ \hfill \frac{\mathrm{d}{E}_h}{\mathrm{d}t}& =(v{S}_l+S_h){\lambda }_h-({\sigma }_h+{\mu }_h)E_h,\hfill \\ \hfill \frac{\mathrm{d}P}{\mathrm{d}t}& ={\sigma }_h{E}_h-(\omega +{\mu }_h)P,\hfill \\ \hfill \frac{\mathrm{d}{I}_1}{\mathrm{d}t}& =\omega P-(\theta +k_1+{\tau }_1+{\mu }_h)I_1,\hfill \\ \hfill \frac{\mathrm{d}{I}_2}{\mathrm{d}t}& =\theta I_1-(k_2+{\delta }_i+{\tau }_2+{\mu }_h)I_2,\hfill \\ \hfill \frac{\mathrm{d}H}{\mathrm{d}t}& =k_1{I}_1+k_2{I}_2-({\delta }_h+{\tau }_3+{\mu }_h)H,\hfill \\ \hfill \frac{\mathrm{d}{R}_h}{\mathrm{d}t}& ={\tau }_1{I}_1+{\tau }_2{I}_2+{\tau }_{3}H-{\mu }_h{R}_h,\hfill \\ \hfill \frac{\mathrm{d}{S}_r}{\mathrm{d}t}& ={\pi }_r-{\lambda }_r{S}_r-{\mu }_r{S}_r,\hfill \\ \hfill \frac{\mathrm{d}{E}_r}{\mathrm{d}t}& ={\lambda }_r{S}_r-({\sigma }_r+{\mu }_r)E_r,\hfill \\ \hfill \frac{\mathrm{d}{I}_r}{\mathrm{d}t}& ={\sigma }_r{E}_r-({\delta }_r+{\tau }_r+{\mu }_r)I_r,\hfill \\ \hfill \frac{\mathrm{d}{R}_r}{\mathrm{d}t}& ={\tau }_r{I}_r-{\mu }_r{R}_r,\hfill \end{array}$$

(1)

The forces of infection for humans and rats are provided here by the model (1) variables, respectively.

$$\begin{array}{c}\hfill {\lambda }_h=\frac{{\beta }_{hh}(\alpha P+\eta I_1+I_2)}{N_h}+{\beta }_{rh}\frac{I_r}{N_r}\mathrm{and}{\lambda }_r=\frac{{\beta }_{rr}{I}_r}{N_r}.\end{array}$$

(2)

where, $\frac{{\beta }_{hh}(\alpha P+\eta I_1+I_2)}{N_h}$, ${\beta }_{rh}\frac{I_r}{N_r}$, and $\frac{{\beta }_{rr}{I}_r}{N_r}$ indicate human-to-human, rodent-to-human, and rodent-to-rodent transmission routes, respectively.

2.3. Fundamental Properties of the Mpoxv Model

As model (1) scrutinises the transmission of mpoxv in humans and rodents, all the parameters and variables of the model are taken to be strictly positive. To explore the fundamental qualitative properties of model (1), the rate of change of the overall populations of humans $N_h^{\prime }\left(t\right)$ and rodents $N_r^{\prime }\left(t\right)$ are firstly determined. Hence

$$\begin{array}{c}\hfill \frac{d{N}_h}{dt}={\pi }_h-{\mu }_h{N}_h-{\delta }_i{I}_2-{\delta }_{h}H⩽{\pi }_h-{\mu }_h{N}_h,\end{array}$$

(3)

and

$$\begin{array}{c}\hfill \frac{d{N}_r}{dt}={\pi }_r-{\mu }_r{N}_r-{\delta }_i{I}_r⩽{\pi }_r-{\mu }_r{N}_r.\end{array}$$

(4)

2.3.1. Positivity of Model Solutions

To be epidemiologically significant, the mpoxv model (1) must show that all of the state variables are non-negative for all time $t>0$. Alternatively, solutions of model system (1) with positive initial data will remain positive for all times $t>0$.

Lemma 1.

Let the initial data $\Theta \left(0\right)\ge 0$, where $\Theta \left(t\right)=\left(S_l,S_h,E_h,P,I_1,I_2,H,R_h,S_r,E_r,I_r,R_r\right)$. Then the solutions $\Phi \left(t\right)$ of model (1) are positive for all $t>0$.

Proof. 

Let $t^1=sup\left\{t>0:\Theta \left(t\right)>0\in [0,t]\right\}.$ Thus, $t^1>0$. Since, ${\pi }_h>0,$ it follows from the first equation of (1) that

$$\frac{d{S}_l}{dt}\ge -\left({\lambda }_h\left(t\right)+{\mu }_h\right)S_l.$$

(5)

Applying comparison theorem in conjunction with separation of variables method, we have

$$S_l\left(t^1\right)\ge S_l\left(0\right)\exp \left[-\left({\int }_0^{t_1}{\lambda }_h\left(u\right)du+{\mu }_h{t}_1\right)\right]>0.$$

(6)

Hence, $S_l\left(t\right)>0$. Similarly, it can be shown that the remaining components of $\Theta \left(t\right)$, i.e., $S_h,E_h,P,I_1,I_2,H,R_h,S_r,E_r,I_r,R_r$ are all non-negative for all $t>0$. Hence, $\Theta \left(t\right)>0$ for all $t>0$.    □

2.3.2. Invariant Region

The following biological feasible region is considered,

$$\begin{array}{c}\hfill \Omega =\left\{(S_l,S_h,E_h,P,I_1,I_2,H,R_h,S_r,E_r,I_r,R_r)\in {\mathbb{R}}_{+}^{12}:N_h⩽\frac{{\pi }_h}{{\mu }_h},N_r⩽\frac{{\pi }_r}{{\mu }_r}\right\}.\end{array}$$

To ensure that any system solutions that begin in the area $\Omega$ remain in the area $\Omega$ for all non-negative times t (i.e., $t⩾0$), the variables $N_h$ and $N_r$ given in Equations (3) and (4) have been simplified. Until now, the region $\Omega$ has been positively invariant; thus, solutions only apply to this region are sufficient. Thus,  model (1) satisfies the results for a normal existence, uniqueness, and continuation as per previous works [26,27,28,29].

3. Analytical Results

3.1. Disease-Free Equilibrium and Basic Reproduction Number

Disease-free equilibrium (DFE) for governing system (1), denoted by ${\Gamma }^0$, is derived by solving the right side of model (1) to zero: $\frac{d{S}_l}{dt}=\frac{d{S}_h}{dt}=\frac{d{E}_h}{dt}=\frac{dP}{dt}==\frac{d{I}_1}{dt}=\frac{d{I}_2}{dt}=\frac{dH}{dt}=\frac{d{R}_h}{dt}=\frac{d{S}_r}{dt}=\frac{d{E}_r}{dt}=\frac{d{I}_r}{dt}=\frac{d{R}_r}{dt}=0$. This produces $S_h^0=\frac{\rho {\pi }_h}{{\mu }_h}$, $S_l^0=\frac{(1-\rho ){\pi }_h}{{\mu }_h}$, $E_h^0=P^0=I_1^0=I_2^0=H^{\ast }=R_h^0=0$, and $S_r^0=\frac{{\pi }_r}{{\mu }_r}$, $E_r^0=I_r^0=R_r^0=0$ for human and rodents, respectively. The DFE for model (1) is obtained below

$$\begin{array}{c}\hfill {\Gamma }^0=\{S_l^0,S_h^0,E^0,P^0,I_1^0,I_2^0,H^0,R_h^0,S_r^0,E_r^0,I_r^0,R_r^0\}=\left\{\frac{(1-\rho ){\pi }_h}{{\mu }_h},\frac{{\pi }_h\rho }{{\mu }_h},0,0,0,0,0,0,\frac{{\pi }_r}{{\mu }_r},0,0,0\right\}\end{array}$$

(7)

3.2. Basic Reproduction Number

The concept of the next-generation matrix (NGM) [30] is used to compute the basic reproduction number ($Ro$) for the governing model (1) [24,30,31,32]. The linear stability of $\Gamma 0$ is obtained by implementing NGM technique for model (1) [30,31], with F matrices, designating the infectious new terms, and V, designating the other transferring terms, with the general formula ${\mathcal{R}}_0={\rho }_{r}F{V}^{-1}$, where the ${\rho }_r$ in the relation characterises the spectral radius of the NGM.

Therefore, the ${\mathcal{R}}_0$ is computed as follows.

$$\mathcal{F}=\left[\begin{array}{c}(v{S}_l^0+S_h^0){\lambda }_h^0\\ 0\\ 0\\ 0\\ 0\\ S_r^0{\lambda }_r^0\\ 0\end{array}\right]\mathrm{and}\mathcal{V}=\left[\begin{array}{c}{a}_1{E}_h^0\\ -{\sigma }_h{E}_h^0+a_2{P}^0\\ -\omega P^0+a_3{I}_1^0\\ -\theta I_1^0+a_4{I}_2^0\\ -k_1{I}_1^0-k_2{I}_2^0+a_5{H}^0\\ a_6{E}_r^0\\ -{\sigma }_r{E}_r^0+a_7{I}_r^0\end{array}\right],$$

As a result, the definitions of the mpoxv infection and transition matrices are respectively computed as:

$$F=\left[\begin{array}{ccccccc}0& m_1& m_2& m_3& 0& 0& m_4\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& m_5\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right]\mathrm{and}V=\left[\begin{array}{ccccccc}{a}_1& 0& 0& 0& 0& 0& 0\\ -{\sigma }_h& a_2& 0& 0& 0& 0& 0\\ 0& -\omega & a_3& 0& 0& 0& 0\\ 0& 0& -{\theta }_1& a_4& 0& 0& 0\\ 0& 0& -k_1& -k_2& a_5& 0& 0\\ 0& 0& 0& 0& 0& a_6& 0\\ 0& 0& 0& 0& 0& -{\sigma }_r& a_7\end{array}\right],$$

where $a_1={\sigma }_h+{\mu }_h$, $a_2=\omega +{\mu }_h$, $a_3=\theta +k_1+{\tau }_1+{\mu }_h$, $a_4=k_2+{\delta }_i+{\tau }_2+{\mu }_h$, $a_5={\delta }_h+{\tau }_3+{\mu }_h$, $a_6={\sigma }_r+{\mu }_r$, $a_7={\delta }_r+{\tau }_r+{\mu }_r$, $m_1=\left(v(1-\rho )+\rho \right){\beta }_{hh}\alpha$, $m_2=\left(v(1-\rho )+\rho \right){\beta }_{hh}\eta$, $m_3=\left(v(1-\rho )+\rho \right){\beta }_{hh}$, $m_4=\frac{{\beta }_{rh}{\pi }_r}{{\mu }_r}$, and $m_5=\frac{{\beta }_{rr}{\pi }_r}{{\mu }_r}$. Straightforward computation produces

$$V^{-1}=\left[\begin{array}{ccccccc}{a_1}^{-1}& 0& 0& 0& 0& 0& 0\\ \frac{{\sigma }_h}{a_2{a}_1}& {a_2}^{-1}& 0& 0& 0& 0& 0\\ \frac{\omega {\sigma }_h}{a_3{a}_2{a}_1}& \frac{\omega }{a_3{a}_2}& {a_3}^{-1}& 0& 0& 0& 0\\ \frac{\theta \omega {\sigma }_h}{a_4{a}_3{a}_2{a}_1}& \frac{\theta \omega }{a_4{a}_3{a}_2}& \frac{\theta }{a_4{a}_3}& {a_4}^{-1}& 0& 0& 0\\ \frac{\omega {\sigma }_h\left(k_1{a}_4+k_2\theta \right)}{a_4{a}_3{a}_2{a}_1{a}_5}& \frac{\omega \left(k_1{a}_4+k_2\theta \right)}{a_4{a}_3{a}_2{a}_5}& \frac{k_1{a}_4+k_2\theta }{a_4{a}_3{a}_5}& \frac{k_2}{a_4{a}_5}& {a_5}^{-1}& 0& 0\\ 0& 0& 0& 0& 0& {a_6}^{-1}& 0\\ 0& 0& 0& 0& 0& \frac{{\sigma }_r}{a_6{a}_7}& {a_7}^{-1}\end{array}\right].$$

and

$$F·V^{-1}=\left[\begin{array}{ccccccc}\frac{m_1{\sigma }_h}{a_1{a}_2}+\frac{m_2\omega {\sigma }_h}{a_1{a}_2{a}_3}+\frac{m_3\theta \omega {\sigma }_h}{a_1{a}_2{a}_3{a}_4}& \frac{m_1}{a_2}+\frac{m_2\omega }{a_3{a}_2}+\frac{m_3\theta \omega }{a_3{a}_2{a}_4}& \frac{m_2}{a_3}+\frac{m_3\theta }{a_3{a}_4}& \frac{m_3}{a_4}& 0& \frac{m_4{\sigma }_r}{a_6{a}_7}& \frac{m_4}{a_7}\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& \frac{m_5{\sigma }_r}{a_6{a}_7}& \frac{m_5}{a_7}\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right].$$

Hence, the reproduction numbers for human and rodent populations are respectively obtained as:

$$\begin{array}{c}\hfill {\mathcal{R}}_0^h=\frac{{\sigma }_h\left(\omega a_4{m}_2+\omega m_3\theta +a_3{a}_4{m}_1\right)}{a_4{a}_3{a}_2{a}_1},\mathrm{and}{\mathcal{R}}_0^r=\frac{m_5{\sigma }_r}{a_6{a}_7},\end{array}$$

So that the ${\mathcal{R}}_0$ is now given as

$$\begin{array}{c}\hfill {\mathcal{R}}_0=max\{{\mathcal{R}}_0^h,{\mathcal{R}}_0^r\}.\end{array}$$

(8)

where $m_1=\left(v(1-\rho )+\rho \right){\beta }_{hh}\alpha$, $m_2=\left(v(1-\rho )+\rho \right){\beta }_{hh}\eta$, $m_3=\left(v(1-\rho )+\rho \right){\beta }_{hh}$, $m_4=\frac{{\beta }_{rh}{\pi }_r}{{\mu }_r}$, and $m_5=\frac{{\beta }_{rr}{\pi }_r}{{\mu }_r}$.

By considering [30] and the local stability of the DFE of model (1), we state the following useful result. Similar results can be found in [33].

Theorem 1.

The DFE of model (1) is locally-asymptotically stable whenever ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$, with ${\mathcal{R}}_0=max\{{\mathcal{R}}_0^h,{\mathcal{R}}_0^r\}$.

The epidemiological essence of the results above implies that a measly intake of mpoxv infection would not result in a substantial outbreak when ${\mathcal{R}}_0<1$. For the mpoxv containment measure, the prerequisite for making ${\mathcal{R}}_0<1$ is suitable but unnecessary. As a result, when ${\mathcal{R}}_0$ is decreased to less than 1, the disease eventually disappears, and when ${\mathcal{R}}_0$ is increased to more than 1, the disease lingers. Hence, efficacious intervention strategies are mandated to mitigate the illness effectively [31].

3.3. Endemic Equilibrium and Its Stability

3.3.1. Endemic Equilibrium

Invasion of mpoxv indicates that at least one infected compartment is not empty. The vector field of (1) is set to zero, and we performed some algebraic calculations, which resulted in an endemic equilibrium (EE) state. Thus, the EE is calculated below.

Suppose the region $\Gamma$ at EE is given by

$$\begin{array}{c}\hfill {\Gamma }^{\ast }=\{S_l^{\ast },S_h^{\ast },E_h^{\ast },P^{\ast },I_1^{\ast },I_2^{\ast },H^{\ast },R_h^{\ast },S_r^{\ast },E_r^{\ast },I_r^{\ast },R_r^{\ast }\}.\end{array}$$

(9)

In forms of ${\lambda }_h^{\ast }$ and ${\lambda }_r^{\ast }$, the EE becomes computed as follows.

$$\begin{array}{cc}\hfill & S_l^{\ast }=\frac{(1-\rho ){\pi }_h}{v{\lambda }_h+{\mu }_h},S_h^{\ast }=\frac{\rho {\pi }_h}{{\lambda }_h+{\mu }_h},E_h^{\ast }=\frac{{\lambda }_h{\pi }_h\left((1-\rho )v{\lambda }_h+(1-\rho )v{\mu }_h+\rho v{\lambda }_h+\rho {\mu }_h\right)}{a_1\left(v{{\lambda }_h}^2+v{\lambda }_h{\mu }_h+{\lambda }_h{\mu }_h+{{\mu }_h}^2\right)},\hfill \\ \hfill & P^{\ast \ast }=\frac{{\lambda }_h{\sigma }_h{\pi }_h\left((1-\rho )v{\lambda }_h+(1-\rho )v{\mu }_h+\rho v{\lambda }_h+\rho {\mu }_h\right)}{a_2{a}_1\left(v{{\lambda }_h}^2+v{\lambda }_h{\mu }_h+{\lambda }_h{\mu }_h+{{\mu }_h}^2\right)},\hfill \\ \hfill & I_1^{\ast }=\frac{\omega {\lambda }_h{\sigma }_h{\pi }_h\left((1-\rho )v{\lambda }_h+(1-\rho )v{\mu }_h+\rho v{\lambda }_h+\rho {\mu }_h\right)}{a_3{a}_2{a}_1\left(v{{\lambda }_h}^2+v{\lambda }_h{\mu }_h+{\lambda }_h{\mu }_h+{{\mu }_h}^2\right)},\hfill \\ \hfill & I_2^{\ast }=\frac{\omega \theta {\lambda }_h{\sigma }_h{\pi }_h\left((1-\rho )v{\lambda }_h+(1-\rho )v{\mu }_h+\rho v{\lambda }_h+\rho {\mu }_h\right)}{a_4{a}_3{a}_2{a}_1\left(v{{\lambda }_h}^2+v{\lambda }_h{\mu }_h+{\lambda }_h{\mu }_h+{{\mu }_h}^2\right)},\hfill \\ \hfill & H^{\ast }=\frac{\left(\theta k_2+a_4{k}_1\right)\omega {\lambda }_h{\sigma }_h{\pi }_h\left((1-\rho )v{\lambda }_h+(1-\rho )v{\mu }_h+\rho v{\lambda }_h+\rho {\mu }_h\right)}{a_4{a}_3{a}_2{a}_1\left(v{{\lambda }_h}^2+v{\lambda }_h{\mu }_h+{\lambda }_h{\mu }_h+{{\mu }_h}^2\right)a_5},\hfill \\ \hfill & R_h^{\ast }=\frac{\left(\theta a_5{\tau }_2+\theta k_2{\tau }_3+a_4{a}_5{\tau }_1+a_4{k}_1{\tau }_3\right)\omega {\lambda }_h{\sigma }_h{\pi }_h\left((1-\rho )v{\lambda }_h+(1-\rho )v{\mu }_h+\rho v{\lambda }_h+\rho {\mu }_h\right)}{a_4{a}_3{a}_2{a}_1\left(v{{\lambda }_h}^2+v{\lambda }_h{\mu }_h+{\lambda }_h{\mu }_h+{{\mu }_h}^2\right){\mu }_h{a}_5},\hfill \\ \hfill & S_r^{\ast }=\frac{{\pi }_r}{{\lambda }_r+{\mu }_r},E_r^{\ast }=\frac{{\lambda }_r{\pi }_r}{a_6\left({\lambda }_r+{\mu }_r\right)},I_r^{\ast }=\frac{{\sigma }_r{\lambda }_r{\pi }_r}{a_7{a}_6\left({\lambda }_r+{\mu }_r\right)},\mathrm{and}{R}_r^{\ast }=\frac{{\tau }_r{\sigma }_r{\lambda }_r{\pi }_r}{{\mu }_r{a}_7{a}_6\left({\lambda }_r+{\mu }_r\right)}.\hfill \end{array}$$

(10)

Epidemiologically, the presence of EE demonstrates that the mpoxv circulates and prevails in the community by indicating that at least one of the infected classes in the model is not empty.

3.3.2. Global Stability Analysis of the Endemic Equilibrium

Herein, we examined the interior-of-the-feasible-region solutions for the governing model that converge to the unique EE provided by ${\Gamma }^{\ast }$ only if ${\mathcal{R}}_0>1$. At ${\Gamma }^{\ast }$, the mpoxv will spread over the neighbourhood and endure. With the help of the Lyapunov function (LF) way, which can be possible by generating the LF from the governing system to confirm the general stability of the EE; see [26,27,34] and [34] for examples of other modelling assessments that have extensively employed this strategy.

Theorem 2.

The EE, ${\Gamma }^{\ast \ast }$, is globally-asymptotically stable (GAS) in the Ω region only if ${\mathcal{R}}_0>1$, and the condition (i) $\left(1-\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right)\left(1-\frac{P{\lambda }_h^{\ast }}{P^{\ast }{\lambda }_h}\right)\ge 0$, (ii) $\left(2+\frac{P}{P^{\ast }}\right)\le \frac{I_1^{\ast }P}{I_1{P}^{\ast }}+\frac{H}{H^{\ast }}+\frac{H^{\ast }{I}_1}{H{I}_1^{\ast }}$, and (iii) $\left(2+\frac{I_1}{I_1^{\ast }}\right)\le \frac{I_2^{\ast }{I}_1}{I_2{I}_1^{\ast }}+\frac{H}{H^{\ast }}+\frac{H^{\ast }{I}_2}{H{I}_2^{\ast }}$, hold.

For the proof of the above Theorem 2, see below.

Proof. 

Considering the technique as in [26,27,28,35,36] by constructing an LF as follows:

$$\begin{array}{cc}\hfill \Xi \left(t\right)=& {\zeta }_1\left(S_l-S_l^{\ast }-S_l^{\ast }ln\frac{S_l}{S_l^{\ast }}\right)+{\zeta }_2\left(S_h-S_h^{\ast }-S_h^{\ast }ln\frac{S_h}{S_h^{\ast }}\right)+{\zeta }_3\left(E_h-E_h^{\ast }-E_h^{\ast }ln\frac{E_h}{E_h^{\ast }}\right)+\hfill \\ \hfill & {\zeta }_4\left(P-P^{\ast }-P^{\ast }ln\frac{P}{P^{\ast }}\right)+{\zeta }_5\left(I_1-I_1^{\ast }-I_1^{\ast }ln\frac{I_1}{I_1^{\ast }}\right)+{\zeta }_6\left(I_2-I_2^{\ast }-I_2^{\ast }ln\frac{I_2}{I_2^{\ast }}\right)+\hfill \\ \hfill & {\zeta }_7\left(H-H^{\ast }-H^{\ast }ln\frac{H}{H^{\ast }}\right)+{\zeta }_8\left(S_r-S_r^{\ast }-S_r^{\ast }ln\frac{S_r}{S_r^{\ast }}\right)+{\zeta }_9\left(E_r-E_r^{\ast }-E_r^{\ast }ln\frac{E_r}{E_r^{\ast }}\right)+\hfill \\ \hfill & {\zeta }_{10}\left(I_r-I_r^{\ast }-I_r^{\ast }ln\frac{I_r}{I_r^{\ast }}\right).\hfill \end{array}$$

(11)

implies the derivative of the LF (11), along with the solutions of (1), are given by

$$\begin{array}{cc}\hfill \dot{\Xi }\left(t\right)=& {\zeta }_1\left(1-\frac{S_l^{\ast }}{S_l}\right)\dot{S_l}+{\zeta }_2\left(1-\frac{S_h^{\ast }}{S_h}\right)\dot{S_h}+{\zeta }_3\left(1-\frac{E^{\ast }}{E}\right)\dot{E}+{\zeta }_4\left(1-\frac{P^{\ast }}{P}\right)\dot{P}+{\zeta }_5\left(1-\frac{I_1^{\ast }}{I_1}\right)\dot{I_1}+\hfill \\ \hfill & {\zeta }_6\left(1-\frac{I_2^{\ast }}{I_2}\right)\dot{I_2}+{\zeta }_7\left(1-\frac{H^{\ast }}{H}\right)\dot{H}+{\zeta }_8\left(1-\frac{S_r^{\ast }}{S_r}\right)\dot{S_r}+{\zeta }_9\left(1-\frac{E_r^{\ast }}{E_r}\right)\dot{E_r}+{\zeta }_{10}\left(1-\frac{I_r^{\ast }}{I_r}\right)\dot{I_r}.\hfill \end{array}$$

(12)

By direct computation from Equation (12), we have

$$\begin{array}{cc}\hfill {\zeta }_1\left(1-\frac{S_l^{\ast }}{S_l}\right)\dot{S_l}& ={\zeta }_1\left(1-\frac{S_l^{\ast }}{S_l}\right)\left(\pi (1-\rho )-v{\lambda }_h{S}_l-{\mu }_h{S}_l\right)\hfill \\ \hfill & ={\zeta }_1\left(1-\frac{S_l^{\ast }}{S_l}\right)\left(v{\lambda }_h^{\ast }{S}_l^{\ast }+{\mu }_h{S}_l^{\ast }-v{\lambda }_h{S}_l-{\mu }_h{S}_l\right)\hfill \\ \hfill & ={\zeta }_{1}v{\lambda }_h^{\ast }{S}_l^{\ast }\left(1-\frac{S_l^{\ast }}{S_l}\right)\left(1-\frac{{\lambda }_h{S}_l}{{\lambda }_h^{\ast }{S}_l^{\ast }}\right)-{\zeta }_1{\mu }_h\frac{{(S_l-S_l^{\ast })}^2}{S_l}\hfill \\ \hfill & \le {\zeta }_{1}v{\lambda }_h^{\ast }{S}_l^{\ast }\left(1-\frac{{\lambda }_h{S}_l}{{\lambda }_h^{\ast }{S}_l^{\ast }}-\frac{S_l^{\ast }}{S_l}+\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right),\hfill \end{array}$$

(13)

and

$$\begin{array}{cc}\hfill {\zeta }_2\left(1-\frac{S_h^{\ast }}{S_h}\right)\dot{S_h}& ={\zeta }_2\left(1-\frac{S_h^{\ast }}{S_h}\right)\left(\pi \rho -{\lambda }_h{S}_h-{\mu }_h{S}_h\right)\hfill \\ \hfill & ={\zeta }_2\left(1-\frac{S_h^{\ast }}{S_h}\right)\left({\lambda }_h^{\ast }{S}_h^{\ast }+{\mu }_h{S}_h^{\ast }-{\lambda }_h{S}_h-{\mu }_h{S}_h\right)\hfill \\ \hfill & ={\zeta }_2{\lambda }_h^{\ast }{S}_h^{\ast }\left(1-\frac{S_h^{\ast }}{S_h}\right)\left(1-\frac{{\lambda }_h{S}_h}{{\lambda }_h^{\ast }{S}_h^{\ast }}\right)-{\zeta }_2{\mu }_h\frac{{(S_h-S_h^{\ast })}^2}{S_h}\hfill \\ \hfill & \le {\zeta }_2{\lambda }_h^{\ast }{S}_h^{\ast }\left(1-\frac{{\lambda }_h{S}_h}{{\lambda }_h^{\ast }{S}_h^{\ast }}-\frac{S_h^{\ast }}{S_h}+\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right),\hfill \end{array}$$

(14)

and

$$\begin{array}{cc}\hfill {\zeta }_3\left(1-\frac{E_h^{\ast }}{E_h}\right)\dot{E_h}& ={\zeta }_3\left(1-\frac{E_h^{\ast }}{E_h}\right)\left(v\lambda S_l+\lambda S_h-a_1{E}_h\right)\hfill \\ \hfill & ={\zeta }_3\left(1-\frac{E_h^{\ast }}{E_h}\right)\left(v{\lambda }_h{S}_l+{\lambda }_h{S}_h-(v{\lambda }_h^{\ast }{S}_l^{\ast }+{\lambda }_h^{\ast }{S}_h^{\ast })\frac{E_h}{E_h^{\ast }}\right)\hfill \\ \hfill & ={\zeta }_{3}v{\lambda }_h^{\ast }{S}_l^{\ast }\left(1-\frac{E_h^{\ast }}{E_h}\right)\left(\frac{{\lambda }_h{S}_l}{{\lambda }_h^{\ast }{S}_l^{\ast }}-\frac{E_h}{E_h^{\ast }}\right)+{\zeta }_3{\lambda }_h^{\ast }{S}_h^{\ast }\left(1-\frac{E_h^{\ast }}{E_h}\right)\left(\frac{{\lambda }_h{S}_h}{{\lambda }_h^{\ast }{S}_h^{\ast }}-\frac{E_h}{E_h^{\ast }}\right)\hfill \\ \hfill & ={\zeta }_{3}v{\lambda }_h^{\ast }{S}_l^{\ast }\left(\frac{\lambda S_l}{{\lambda }^{\ast }{S}_l^{\ast }}-\frac{E_h}{E_h^{\ast }}-\frac{{\lambda }_h{S}_l{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_l^{\ast }{E}_h}+1\right)+{\zeta }_3{\lambda }_h^{\ast }{S}_h^{\ast }\left(\frac{{\lambda }_h{S}_h}{{\lambda }_h^{\ast }{S}_h^{\ast }}-\frac{E_h}{E_h^{\ast }}-\frac{{\lambda }_h{S}_h{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_h^{\ast }{E}_h}+1\right),\hfill \end{array}$$

(15)

and

$$\begin{array}{cc}\hfill {\zeta }_4\left(1-\frac{P^{\ast }}{P}\right)\dot{P}& ={\zeta }_4\left(1-\frac{P^{\ast }}{P}\right)\left({\sigma }_h{E}_h-a_{2}P\right)\hfill \\ \hfill & ={\zeta }_4\left(1-\frac{P^{\ast }}{P}\right)\left({\sigma }_h{E}_h-{\sigma }_h{E}_h^{\ast }\frac{P}{P^{\ast }}\right)\hfill \\ \hfill & ={\zeta }_4{\sigma }_h{E}_h^{\ast }\left(1-\frac{P^{\ast }}{P}\right)\left(\frac{E_h}{E_h^{\ast }}-\frac{P}{P^{\ast }}\right)\hfill \\ \hfill & ={\zeta }_4{\sigma }_h{E}_h^{\ast }\left(\frac{E_h}{E_h^{\ast }}-\frac{P}{P^{\ast }}-\frac{P^{\ast }{E}_h}{P{E}_h^{\ast }}+1\right);\hfill \end{array}$$

(16)

similarly,

$$\begin{array}{cc}\hfill {\zeta }_5\left(1-\frac{I_1^{\ast }}{I_1}\right)\dot{I_1}& ={\zeta }_5\omega P^{\ast }\left(\frac{P}{P^{\ast }}-\frac{I_1}{I_1^{\ast }}-\frac{I_1^{\ast }P}{I_1{P}^{\ast }}+1\right),\hfill \end{array}$$

(17)

and

$$\begin{array}{cc}\hfill {\zeta }_6\left(1-\frac{I_2^{\ast }}{I_2}\right)\dot{I_2}& ={\zeta }_6\theta I_1^{\ast }\left(\frac{I_1}{I_1^{\ast }}-\frac{I_2}{I_2^{\ast }}-\frac{I_2^{\ast }{I}_1}{I_{1}2I_1^{\ast }}+1\right),\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\zeta }_7\left(1-\frac{H^{\ast }}{H}\right)\dot{H}& ={\zeta }_7{k}_1{I}_1^{\ast }\left(\frac{I_1}{I_1^{\ast }}-\frac{H}{H^{\ast }}-\frac{H^{\ast }{I}_1}{H{I}_1^{\ast }}+1\right)+{\zeta }_7{k}_2{I}_2^{\ast }\left(\frac{I_2}{I_2^{\ast }}-\frac{H}{H^{\ast }}-\frac{H^{\ast }{I}_2}{H{I}_2^{\ast }}+1\right),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\zeta }_8\left(1-\frac{S_r^{\ast }}{S_r}\right)\dot{S_r}& ={\zeta }_8{\lambda }_r^{\ast }{S}_r^{\ast }\left(1-\frac{{\lambda }_r{S}_r}{{\lambda }_r^{\ast }{S}_r^{\ast }}-\frac{S_r^{\ast }}{S_r}+\frac{{\lambda }_r}{{\lambda }_r^{\ast }}\right),\hfill \end{array}$$

(20)

and

$$\begin{array}{cc}\hfill {\zeta }_9\left(1-\frac{E_r^{\ast }}{E_r}\right)\dot{E_r}& ={\zeta }_9{\lambda }_r^{\ast }{S}_r^{\ast }\left(\frac{{\lambda }_r{S}_r}{{\lambda }_r^{\ast }{S}_r^{\ast }}-\frac{E_r}{E_r^{\ast }}-\frac{{\lambda }_r{S}_r{E}_r^{\ast }}{{\lambda }_r^{\ast }{S}_r^{\ast }{E}_r}+1\right),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\zeta }_{10}\left(1-\frac{I_r^{\ast }}{I_r}\right)\dot{I_r}& ={\zeta }_{10}{\sigma }_r{E}_r^{\ast }\left(\frac{E_r}{E_r^{\ast }}-\frac{I_r}{I_r^{\ast }}-\frac{I_r^{\ast }{E}_r}{I_r{E}_r^{\ast }}+1\right).\hfill \end{array}$$

(22)

Substituting ${\zeta }_1={\zeta }_2={\zeta }_3={\zeta }_7={\zeta }_8={\zeta }_9=1$, ${\zeta }_4=\frac{{\lambda }^{\ast }}{{\sigma }_h{E}_h^{\ast }}(v{S}_l^{\ast }+S_h^{\ast })$, and ${\zeta }_5=\frac{k_1{I}_1^{\ast }}{\omega P^{\ast }}$, ${\zeta }_6=\frac{k_2{I}_2^{\ast }}{\theta I_1^{\ast }}$, ${\zeta }_{10}=\frac{{\lambda }_r^{\ast }{S}_r^{\ast }{I}_2^{\ast }}{{\sigma }_r{E}_r^{\ast }}$ and Equations (13)–(22) into Equation (12), we have

$$\begin{array}{cc}\hfill \dot{\Xi }\left(t\right)\le & v{\lambda }_h^{\ast }{S}_l^{\ast }\left(2-\frac{S_l^{\ast }}{S_l}-\frac{E_h}{E_h^{\ast }}-\frac{{\lambda }_h{S}_l{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_l^{\ast }{E}_h}+\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right)+\hfill \\ \hfill & {\lambda }_h^{\ast }{S}_h^{\ast }\left(2-\frac{S_h^{\ast }}{S_h}-\frac{E_h}{E_h^{\ast }}-\frac{{\lambda }_h{S}_h{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_h^{\ast }{E}_h}+\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right)+\hfill \\ \hfill & v{\lambda }_h^{\ast }{S}_l^{\ast }\left(\frac{E_h}{E_h^{\ast }}-\frac{P}{P^{\ast }}-\frac{P^{\ast }{E}_h}{P{E}_h^{\ast }}+1\right)+\hfill \\ \hfill & {\lambda }_h^{\ast }{S}_h^{\ast }\left(\frac{E_h}{E_h^{\ast }}-\frac{P}{P^{\ast }}-\frac{P^{\ast }{E}_h}{P{E}_h^{\ast }}+1\right)+\hfill \\ \hfill & k_1{I}_1^{\ast }\left(\frac{P}{P^{\ast }}-\frac{I_1^{\ast }P}{I_1{P}^{\ast }}-\frac{H}{H^{\ast }}-\frac{H^{\ast }{I}_1}{H{I}_1^{\ast }}+2\right)+\hfill \\ \hfill & k_2{I}_2^{\ast }\left(\frac{I_1}{I_1^{\ast }}-\frac{I_2^{\ast }{I}_1}{I_2{I}_1^{\ast }}-\frac{H}{H^{\ast }}-\frac{H^{\ast }{I}_2}{H{I}_2^{\ast }}+2\right)+\hfill \\ \hfill & {\lambda }_r^{\ast }{S}_r^{\ast }\left(2-\frac{S_r^{\ast }}{S_r}-\frac{E_r}{E_r^{\ast }}-\frac{{\lambda }_r{S}_r{E}_r^{\ast }}{{\lambda }_r^{\ast }{S}_r^{\ast }{E}_r}+\frac{{\lambda }_r}{{\lambda }_r^{\ast }}\right)+\hfill \\ \hfill & {\lambda }_r^{\ast }{S}_r^{\ast }\left(\frac{E_r}{E_r^{\ast }}-\frac{I_r}{I_r^{\ast }}-\frac{I_r^{\ast }{E}_r}{I_r{E}_r^{\ast }}+1\right).\hfill \end{array}$$

(23)

Surmise that $u\left(x\right)=1-x+lnx$, implies, $x>0$ yields $u\left(x\right)\le 0$, and, if $x=1$, implies $u\left(x\right)=0$. Thus $x-1\ge ln\left(x\right)$ for any $x>0$ [3,26,27,28,37].

Employing these definitions, with some algebra from (23) and the conditions above, yields

$$\begin{array}{cc}\hfill & \left(2-\frac{S_l^{\ast }}{S_l}-\frac{E_h}{E_h^{\ast }}-\frac{{\lambda }_h{S}_l{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_l^{\ast }{E}_h}+\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right)\hfill \\ \hfill & =\left(-(1-\frac{{\lambda }_h}{{\lambda }_h^{\ast }})(1-\frac{P{\lambda }_h^{\ast }}{P^{\ast }{\lambda }_h})+3-\frac{S_l^{\ast }}{S_l}-\frac{{\lambda }_h{S}_l{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_l^{\ast }{E}_h}-\frac{P{\lambda }_h^{\ast }}{P^{\ast }{\lambda }_h}-\frac{E_h}{E_h^{\ast }}+\frac{P}{P^{\ast }}\right)\hfill \\ \hfill & \le \left(-(\frac{S_l^{\ast }}{S_l}-1)-(\frac{{\lambda }_h{S}_l{E}^{\ast }}{{\lambda }_h^{\ast }{S}_l^{\ast }{E}_h}-1)-(\frac{P{\lambda }_h^{\ast }}{P^{\ast }{\lambda }_h}-1)-\frac{E_h}{E_h^{\ast }}+\frac{P}{P^{\ast }}\right)\hfill \\ \hfill & \le \left(-ln\left(\frac{S_l^{\ast }}{S_l}\frac{\lambda S_l{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_l^{\ast }{E}_h}\frac{P{\lambda }^{\ast }}{P^{\ast }{\lambda }_h}\right)-\frac{E_h}{E_h^{\ast }}+\frac{P}{P^{\ast }}\right)\hfill \\ \hfill & =\left(\frac{P}{P^{\ast }}-ln\left(\frac{P}{P^{\ast }}\right)+ln\left(\frac{E_h}{E_h^{\ast }}\right)-\frac{E_h}{E_h^{\ast }}\right).\hfill \end{array}$$

(24)

Similarly,

$$\begin{array}{cc}\hfill & \left(2-\frac{S_h^{\ast }}{S_h}-\frac{E_h}{E_h^{\ast }}-\frac{\lambda S_h{E}_h^{\ast }}{{\lambda }_h^{\ast }{S}_h^{\ast }{E}_h}+\frac{{\lambda }_h}{{\lambda }_h^{\ast }}\right)\le \left(\frac{P}{P^{\ast }}-ln\left(\frac{P}{P^{\ast }}\right)+ln\left(\frac{E_h}{E_h^{\ast }}\right)-\frac{E_h}{E_h^{\ast }}\right).\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & \left(2-\frac{S_r^{\ast }}{S_r}-\frac{E_r}{E_r^{\ast }}-\frac{{\lambda }_r{S}_r{E}_r^{\ast }}{{\lambda }_r^{\ast }{S}_r^{\ast }{E}_r}+\frac{{\lambda }_r}{{\lambda }_r^{\ast }}\right)\le \left(\frac{I_r}{I_r^{\ast }}-ln\left(\frac{I_r}{I_r^{\ast }}\right)+ln\left(\frac{E_r}{E_r^{\ast }}\right)-\frac{E_r}{E_r^{\ast }}\right).\hfill \end{array}$$

(26)

From Equation (23), we also have

$$\begin{array}{cc}\hfill \frac{E}{E^{\ast }}-\frac{P}{P^{\ast }}-\frac{P^{\ast }E}{P{E}^{\ast }}+1& =\left(u\left(\frac{P^{\ast }E}{P{E}^{\ast }}\right)+\frac{E}{E^{\ast }}-ln\left(\frac{E}{E^{\ast }}\right)-\frac{P}{P^{\ast }}+ln\left(\frac{P}{P^{\ast }}\right)\right)\hfill \\ \hfill & \le \frac{E}{E^{\ast }}-ln\left(\frac{E}{E^{\ast }}\right)+ln\left(\frac{P}{P^{\ast }}\right)-\frac{P}{P^{\ast }}.\hfill \end{array}$$

(27)

Similarly,

$$\begin{array}{c}\hfill \frac{E_r}{E_r^{\ast }}+\frac{I_r}{I_r^{\ast }}-\frac{I_r^{\ast }{E}_r}{I_r{E}_r^{\ast }}+1\le \frac{-I_r}{I_r^{\ast }}+ln\left(\frac{I_r}{I_r^{\ast }}\right)+ln\left(\frac{E_r}{E_r^{\ast }}\right)-\frac{E_r}{E_r^{\ast }}.\end{array}$$

(28)

Hence,

$$\begin{array}{cc}\hfill \dot{\Xi }\left(t\right)\le & v{\lambda }_h^{\ast }{S}_l^{\ast }\left(\frac{P}{P^{\ast }}-ln\left(\frac{P}{P^{\ast }}\right)+ln\left(\frac{E_h}{E_h^{\ast }}\right)-\frac{E_h}{E_h^{\ast }}\right)+\hfill \\ \hfill & {\lambda }_h^{\ast }{S}_h^{\ast }\left(\frac{P}{P^{\ast }}-ln\left(\frac{P}{P^{\ast }}\right)+ln\left(\frac{E_h}{E_h^{\ast }}\right)-\frac{E_h}{E_h^{\ast }}\right)+\hfill \\ \hfill & v{\lambda }_h^{\ast }{S}_l^{\ast }\left(\frac{E_h}{E_h^{\ast }}-ln\left(\frac{E_h}{E_h^{\ast }}\right)+ln\left(\frac{P}{P^{\ast }}\right)-\frac{P}{P^{\ast }}\right)+\hfill \\ \hfill & {\lambda }^{\ast }{S}_h^{\ast }\left(\frac{E_h}{E_h^{\ast }}-ln\left(\frac{E_h}{E_h^{\ast }}\right)+ln\left(\frac{P}{P^{\ast }}\right)-\frac{P}{P^{\ast }}\right)+\hfill \\ \hfill & k_1{I}_1^{\ast }\left(\frac{P}{P^{\ast }}-\frac{I_1^{\ast }P}{I_1{P}^{\ast }}-\frac{H}{H^{\ast }}-\frac{H^{\ast }{I}_1}{H{I}_1}+2\right)+\hfill \\ \hfill & k_2{I}_2^{\ast }\left(\frac{I_1}{I_1^{\ast }}-\frac{I_2^{\ast }{I}_1}{I_2{I}_1^{\ast }}-\frac{H}{H^{\ast }}-\frac{H^{\ast }{I}_2}{H{I}_2}+2\right)+\hfill \\ \hfill & {\lambda }_r^{\ast }{S}_r^{\ast }\left(\frac{I_r}{I_r^{\ast }}-ln\left(\frac{I_r}{I_r^{\ast }}\right)+ln\left(\frac{E_r}{E_r^{\ast }}\right)-\frac{E_r}{E_r^{\ast }}\right)+\hfill \\ \hfill & {\lambda }_r^{\ast }{S}_r^{\ast }\left(\frac{-I_r}{I_r^{\ast }}+ln\left(\frac{I_r}{I_r^{\ast }}\right)-ln\left(\frac{E_r}{E_r^{\ast }}\right)+\frac{E_r}{E_r^{\ast }}\right)\hfill \end{array}$$

(29)

Equations (13)–(29) and conditions(i)–(iii) ensure that $\dot{\Xi \left(t\right)}\le 0$. Further, the equality $\frac{d\Xi }{dt}=0$ holds only if $S_l=S_l^{\ast }$, $S_h=S_h^{\ast }$$E_h=E_h^{\ast }$, $P=P^{\ast }$, $I_1=I_1^{\ast }$, $I_2=I_2^{\ast }$, $H=H^{\ast }$, $S_r=S_r^{\ast }$, $E_r=E_r^{\ast }$, and $I_r=I_r^{\ast }$. Thus, the EE state (10) serves as the only invariant positive set to the governing model (1) contained entirely in $\left\{(S_l,S_h,E_h,P,I_1,I_2,H,S_r,E_r,I_r)\in \Omega :S_l=S_l^{\ast },S_h=S_h^{\ast },E_h=E_h^{\ast },P=P^{\ast },I_1=I_1^{\ast },I_2=I_2^{\ast },H=H^{\ast },S_r=S_r^{\ast },E_r=E_r^{\ast },I_r=I_r^{\ast }\right\}$. Note that $R_h$ and $R_r$ were not included in the GAS analysis since individuals in these classes have already recovered from the mpoxv and are not of much epidemiological interest in this regard. Hence, it follows from LaSalle’s invariance principle [34] that every solution to the Equation (1) with initial conditions in $\Omega$ converge to EE points, ${\Gamma }^{\ast }$, as $t\to \infty$. Thus, the positive EE is globally asymptotically stable.    □

3.4. Bifurcation Analysis

The center manifold theory (CMT) is applied to examine the bifurcation analysis [38]. We look at the parameters in the model (1) that produce forward, or backward bifurcation [3,30,38]. The moment ${\mathcal{R}}_0$ crosses unity from below, a small positive asymptotically stable equilibrium arises, and the DFE loses stability [38]. In the subsequent paragraph, we examine whether forward (exchange of DFE and EE stability at ${\mathcal{R}}_0=1$) or backward bifurcation occurs in the current model (1). Thus, the results below follow.

Theorem 3.

Forward bifurcation occurs in the current mpoxv model (1) at ${\mathcal{R}}_0=1$ only if the bifurcation coefficients, A and B, are negative and positive, respectively.

Proof. 

The proof of the above theorem is based on the CMT [38,39]. Consider the expression ${\displaystyle \frac{dx}{dt}}=f(x,\psi )$, with $\psi$ representing the bifurcation parameter and f being continuously differentiable at least twice in both x and $\psi$. The DFE (${\Gamma }^0$) is the point $(x_0=0,\psi =0)$ and the local stability of the DFE changes at the point $(x_0,0)$ [30]. We now demonstrate that there exists a nontrivial equilibrium around the bifurcation point $(x_0,\psi )$.

Surmise that ${\beta }_{hh}={\beta }_{hh}^{\ast }$ is taken as a bifurcation parameter and taking ${\mathcal{R}}_0=1$ from Equation (3) implies by Theorem A1.1, the DFE, $E_0$, is locally stable only if ${\beta }_{hh}<{\beta }_{hh}^{\ast }$ and unstable when ${\beta }_{hh}>{\beta }_{hh}^{\ast }$. Here, ${\beta }_{hh}={\beta }_{hh}^{\ast }$ is a bifurcation value.

In a simpler way, let $S_l=x_1$, $S_h=x_2$, $E_h=x_3$, $P=x_4$, $I_1=x_5$, $I_2=x_6$, $H=x_7$, $R=x_8$, $S_r=x_9$, $E_r=x_{10}$, $I_r=x_{11}$, and $R_r=x_{12}$, so that $N_h=x_1+x_2+x_3+x_4+x_5+x_6+x_7+x=x_8$ and $N_r=x_9+x_{10}+x_{11}+x_{12}$. Additionally, by using the identical vector notation with $x={(x_1,x_2,\cdots ,x_{12})}^T$, the model (1) as $\frac{dx}{dt}=f\left(x\right)$ where $f={(f_1,f_2,\cdots ,f_{12})}^T$ is as follows:

$$\begin{array}{cc}\hfill f_1& ={\pi }_h(1-\rho )-v{\lambda }_h{x}_1-{\mu }_h{x}_1,\hfill \\ \hfill f_2& ={\pi }_h\rho -{\lambda }_h{x}_2-{\mu }_h{x}_2,\hfill \\ \hfill f_3& =(v{x}_1+x_2){\lambda }_h-a_1{x}_3,\hfill \\ \hfill f_4& ={\sigma }_h{x}_3-a_2{x}_4,\hfill \\ \hfill f_5& =\omega x_4-a_3{x}_5,\hfill \\ \hfill f_6& =\theta x_5-a_4{x}_6,\hfill \\ \hfill f_7& =k_1{x}_5+k_2{x}_6-a_5{x}_7,\hfill \\ \hfill f_8& ={\tau }_1{x}_5+{\tau }_2{x}_6+{\tau }_3{x}_7-{\mu }_h{x}_8,\hfill \\ \hfill f_9& ={\pi }_r-{\lambda }_r{x}_9-{\mu }_r{x}_9,\hfill \\ \hfill f_{10}& ={\lambda }_r{x}_9-a_6{x}_{10},\hfill \\ \hfill f_{11}& ={\sigma }_r{x}_{10}-a_7{x}_{11},\hfill \\ \hfill f_{12}& ={\tau }_r{x}_{11}-{\mu }_r{x}_{12},\hfill \end{array}$$

(30)

and each of the linked forces of infection is represented by

$$\begin{array}{cc}\hfill {\lambda }_h=& \frac{{\beta }_{\mathit{hh}}\left(\alpha x_4+\eta x_5+x_6\right)}{x_1+x_2+x_3+x_4+x_5+x_6+x_7+x_8}+\frac{{\beta }_{\mathit{rh}}{x}_{11}}{x_9+x_{10}+x_{11}+x_{12}},\mathrm{and}\hfill \\ \hfill {\lambda }_r=& \frac{{\beta }_{\mathit{rr}}{x}_{11}}{x_9+x_{10}+x_{11}+x_{12}}.\hfill \end{array}$$

(31)

The Jacobian matrix of system (30), computed from the relation $J\left({\Gamma }^0\right)=F-V$, where F and V are infections and transition matrices, respectively, evaluated at the DFE $\left({\Gamma }^0\right)$ with ${\beta }_{hh}={\beta }_{hh}^{\ast }$, yields:

$$J\left({\Gamma }^0\right)=\left[\begin{array}{cccccccccccc}-{\mu }_h& 0& 0& -n_1& -n_2& -n_3& 0& 0& 0& 0& -n_4& 0\\ 0& -{\mu }_h& 0& -n_5& -n_6& -n_7& 0& 0& 0& 0& -n_8& 0\\ 0& 0& -a_1& n_9& n_{10}& n_{11}& 0& 0& 0& 0& n_{12}& 0\\ 0& 0& {\sigma }_h& -a_2& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \omega & -a_3& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& \theta & -a_4& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& k_1& k_2& -a_5& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& {\tau }_1& {\tau }_2& {\tau }_3& -{\mu }_h& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -{\mu }_r& 0& -{\beta }_{\mathit{rr}}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& -a_6& {\beta }_{\mathit{rr}}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {\sigma }_r& -a_7& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& {\tau }_r& -{\mu }_r\end{array}\right],$$

(32)

where $n_1=-v\left(-1+\rho \right){\beta }_{\mathit{hh}}\alpha$, $n_2=-v\left(-1+\rho \right){\beta }_{\mathit{hh}}\eta$, $n_3=-v\left(-1+\rho \right){\beta }_{\mathit{hh}}$, $n_4=\frac{v{\pi }_h\left(1-\rho \right){\beta }_{\mathit{rh}}{\mu }_r}{{\mu }_h{\pi }_r}$, $n_5=\rho {\beta }_{\mathit{hh}}\alpha$, $n_6=\rho {\beta }_{\mathit{hh}}\eta$, $n_7=\rho {\beta }_{\mathit{hh}}$, $n_8=\frac{{\pi }_h\rho {\beta }_{\mathit{rh}}{\mu }_r}{{\mu }_h{\pi }_r}$, $n_9=-\left(\left(v-1\right)\rho -v\right)\alpha {\beta }_{\mathit{hh}}$, $n_{10}=-\left(\left(v-1\right)\rho -v\right){\beta }_{\mathit{hh}}\eta$, $n_{11}=-\left(\left(v-1\right)\rho -v\right){\beta }_{\mathit{hh}}$, and $n_{12}=-\frac{\left(\left(v-1\right)\rho -v\right){\beta }_{\mathit{rh}}{\pi }_h{\mu }_r}{{\mu }_h{\pi }_r}$.

The Jacobian matrix $J\left({\Gamma }^0\right)$ of the system (30) has a simple zero eigenvalue (other eigenvalues have negative real parts). Thus, the CMT [30,38] may be employed to scrutinise the dynamics of the model (30) located at ${\beta }_{hh}={\beta }_{hh}^{\ast }$ [38] . The following computations are carried out using the notations in [38].

Eigenvectors of $J{\left({\Gamma }^0\right)}_{{\beta }_{hh}={\beta }_{hh}^{\ast }}$: When ${\mathcal{R}}_0=1$, one can observe that the $J\left({\Gamma }^0\right)$ has a right eigenvector (corresponding to the zero eigenvalues), given by $\mathbf{w}={(w_1,w_2,\cdots ,w_{12})}^T$, where

$$\begin{array}{cc}\hfill & w_1=\frac{-{\sigma }_h(a_3{a}_4{n}_1+a_4{n}_2\omega +n_3\omega \theta )}{{\mu }_h{a}_2{a}_3{a}_4}{w}_3,w_2=\frac{-{\sigma }_h(a_3{a}_4{n}_5+a_4{n}_6\omega +n_7\omega \theta )}{{\mu }_h{a}_2{a}_3{a}_4}{w}_3,w_3>0,w_4=\frac{{\sigma }_h}{a_2}{w}_3,\hfill \\ \hfill & w_5=\frac{{\sigma }_h\omega }{a_2{a}_3}{w}_3,w_6=\frac{{\sigma }_h\omega \theta }{a_2{a}_3{a}_4}{w}_3,w_7=\frac{{\sigma }_h\omega (a_4{k}_1+\theta k_2)}{a_2{a}_3{a}_4{a}_5}{w}_3,w_8=\frac{{\sigma }_h\omega ({\tau }_1{a}_4{a}_5+{\tau }_2{a}_5\theta +{\tau }_3{a}_4{k}_1+{\tau }_3\theta k_2)}{a_2{a}_3{a}_4{a}_5{\mu }_h}{w}_3,\hfill \\ \hfill & w_9=0,w_{10}=0,w_{11}=0,\mathrm{and}{w}_{12}=0.\hfill \end{array}$$

In similar way, the components of the left eigenvector of $J\left({\Gamma }^0\right)$ (corresponding to the zero eigenvalues), indicated by $\mathbf{v}=(v_1,v_2,\cdots ,v_9)$, are given by

$$\begin{array}{cc}\hfill & v_1=0,v_2=0,v_3>0,v_4=\frac{a_1}{{\sigma }_h}{v}_3,v_5=\frac{(a_1{a}_2-n_9{\sigma }_h)}{\omega {\sigma }_h}{v}_3,v_6=\frac{(a_1{a}_2{a}_3-a_3{n}_9{\sigma }_h-\omega {\sigma }_h{n}_{10})}{\theta \omega {\sigma }_h}{v}_3,\hfill \\ \hfill & v_7=0,v_8=0,v_9=0,v_{10}=\frac{{\sigma }_r{n}_{12}}{a_6{a}_7-{\sigma }_r{\beta }_{rr}}{v}_3,v_{11}=\frac{a_6{n}_{12}}{a_6{a}_7-{\sigma }_r{\beta }_{rr}}{v}_3,\mathrm{and}{v}_{12}=0.\hfill \end{array}$$

Taking to account the free components (entry) are taken to be $v_3=1$ and $w_3=\frac{1}{A_1+A_2}$ respectively, where

$$\begin{array}{cc}\hfill A_1=& 1+\frac{a_1}{a_2}+\left(\frac{a_1{a}_2-n_9{\sigma }_h}{a_2{a}_3}\right)\hfill \end{array}$$

and

$$\begin{array}{c}\hfill A_2=\frac{a_1{a}_2{a}_3-a_3{n}_9{\sigma }_h-\omega {\sigma }_h{n}_{10}}{a_2{a}_3{a}_4},\end{array}$$

so that $\mathbf{v}·\mathbf{w}=1$ (in line with [38]).

It can be demonstrated that the related bifurcation coefficients, A and B, are given, respectively, by computing the non-zero partial derivatives of $f_i(i=1,\cdots ,9)$:

$$\begin{array}{c}\hfill A=\sum _{k,i,j=1}^9{v}_k{w}_i{w}_j\frac{{\partial }^2{f}_k(0,0)}{\partial x_i\partial x_j}=\frac{2{\mu }_h{\beta }_{hh}{g}_1}{{\pi }_h}(-g_2\rho (1-v)-g_{3}v+w_2)v_3,\end{array}$$

(33)

$$\begin{array}{c}\hfill B=\sum _{k,i=1}^9{v}_k{w}_i\frac{{\partial }^2{f}_k(0,0)}{\partial x_i\partial {\beta }_1}=(v(1-\rho )+\rho )(\alpha w_4+\eta w_5+w_6)v_3,\end{array}$$

(34)

where $g_1=\alpha w_4+\eta w_5+w_6$, $g_2=w_1+w_2+w_3+w_4+w_5+w_6+w_7+w_8$, and $g_3=w_2+w_3+w_4+w_5+w_6+w_7+w_8$.

Since $v_3>0$$w_3>0$, $w_4>0$, $w_5>0$, and $w_6>0$, it implies that $B\ge 0$, and hence we have forward bifurcation if and only if $A<0$ (i.e., when $g_2>0$ and $g_3>0$ since $g_1>0$), else backward bifurcation (i.e., when $A>0$). To prove the above numerically, we used the parameter values given in

Table 2. Table summarizing model parameter values (1).

Parameter	Baseline (Range)	Units	Sources
$N_h$	$1.8\times {10}^8(1.2\times {10}^8$–$2.2\times {10}^8)$	Persons	[3,40]
$N_r$	$0.01\times N_h$	Rodents	Assumed
${\mu }_h$	$0.000045(0.00003$–$0.00006)$	$$\mathrm{Per}\mathrm{day}$$	[3,40]
${\mu }_r$	$0.002(0.001$–$0.006)$	$$\mathrm{Per}\mathrm{day}$$	[41]
${\pi }_h$	$2500(1000$–$5000)$	$$\mathrm{Per}\mathrm{day}$$	[3,24,42]
${\pi }_r$	$0.5(0$–$1)$	$$\mathrm{Per}\mathrm{day}$$	[3,43,44]
$\rho$	$0.8(0$–$1)$	$$\mathrm{Per}\mathrm{day}$$	[45]
v	$0.045(0$–$1)$	Dimensionless	Estimated
${\beta }_{hh}$	$0.03(0.01$–$0.5)$	$$\mathrm{Per}\mathrm{day}$$	Estimated
${\beta }_{rh}$	$0.3045(0.1$–$0.8)$	$$\mathrm{Per}\mathrm{day}$$	Estimated
${\beta }_{rr}$	$0.025(0.01$–$0.5)$	$$\mathrm{Per}\mathrm{day}$$	Estimated
$\alpha$	$0.75(0$–$1)$	Dimensionless	[5]
$\eta$	$0.8(0$–$1)$	Dimensionless	[5]
${\sigma }_h$	$0.033(0.01$–$0.2)$	$$\mathrm{Per}\mathrm{day}$$	[44]
${\sigma }_r$	$0.0083(0.0023$–$0.031)$	$$\mathrm{Per}\mathrm{day}$$	[44]
$\omega$	$0.0042(0.002$–$1)$	$$\mathrm{Per}\mathrm{day}$$	[5]
$\theta$	$0.021(0.01$–$0.1)$	$$\mathrm{Per}\mathrm{day}$$	Estimated
$k_1$	$0.2(0.1$–$0.5)$	$$\mathrm{Per}\mathrm{day}$$	Estimated
$k_2$	$0.4(0.3$–$0.5)$	$$\mathrm{Per}\mathrm{day}$$	Estimated
${\tau }_1$	$0.048(0.001$–$0.075)$	$$\mathrm{Per}\mathrm{day}$$	[5]
${\tau }_2$	$0.05(0.001$–$0.075)$	$$\mathrm{Per}\mathrm{day}$$	[5]
${\tau }_3$	$0.056(0.001$–$0.075)$	$$\mathrm{Per}\mathrm{day}$$	[5]
${\tau }_r$	$0.083(0.01$–$0.25)$	$$\mathrm{Per}\mathrm{day}$$	[44]
${\delta }_i$	$0.0011(0.001$–$0.025)$	$$\mathrm{Per}\mathrm{day}$$	[5]
${\delta }_h$	$0.001(0.001$–$0.025)$	$$\mathrm{Per}\mathrm{day}$$	[5]
${\delta }_r$	$0.057(0.012$–$0.085)$	$$\mathrm{Per}\mathrm{day}$$	[44]

to verify that A ($-7.451088966\times {10}^{-6}$) and B ($4.799895897$) are negative and positive, respectively.

Hence, the mpoxv model (1) exhibits the phenomenon of forward bifurcation at ${\mathcal{R}}_0=1$. No EE appears if ${\mathcal{R}}_0<1$, and the DFE is the only local attractor, but if ${\mathcal{R}}_0>1$, then the EE exists. For this reason, there is a forward bifurcation because, in the neighbourhood of the bifurcation point, the disease prevalence is an increasing function of ${\mathcal{R}}_0$, which implies that the mpoxv model (1) has forward bifurcation property.    □

4. Numerical Results

4.1. Model Prediction

We used the prior approach outlined in [35] to validate model (1) utilizing mpoxv surveillance data produced and released by the NCDC in [23]. Using the R statistical software, Pearson’s chi-square and the least square method are used for fitting the data using the R (version 4.2.1). The model is validated using actual mpoxv scenarios for cases from 1 January to 27 August 2022 (specifically for 35 epidemiological weeks). The model suited well to the mpoxv situation report in Nigeria using suitable biological parameters according to the results. This shows that the model can be used to understand the spread of mpoxv in the community. The total number of cases for 35 epidemiological weeks in the fitting findings of mpoxv confirmed cases are depicted in Figure 2.

We obtained demographic biological parameters for mpoxv in Nigeria from [40]. We also computed other related demographic parameters, including $\pi$ and $\mu$. $\pi$ is evaluated as 8100 per day. Given that, as of last year, the predicted life expectancy in Nigeria was $60.87$ each year, indicating that $\mu =\frac{1}{60.87\times 365}$$\mathrm{per}\mathrm{day}$= $4.5\times {10}^{-5}$$\mathrm{per}\mathrm{day}$. Other than that, everything is set as in Table 2. It is crucial to comprehend that poor resource settings produced most of the mpoxv cases in Nigeria. Figure 2 shows the model-fitting result for 35 epidemiological weeks (i.e., from January through August 2022). The initial conditions employed are given: $S_l\left(0\right)=12\times {10}^7$, $S_h\left(0\right)=6\times {10}^7$, $E_h\left(0\right)=2000$, $P\left(0\right)=70$, $I_1\left(0\right)=2$, $I_2\left(0\right)=1$, $H\left(0\right)=1$, $R_h\left(0\right)=0$, $S_r\left(0\right)=S_h\left(0\right)\times {10}^{-2}$, $E_r\left(0\right)=800$, $I_r\left(0\right)=16$, $R_r\left(0\right)=2$. Moreover, from the prediction result, we observed that the mpoxv cases had reached an all-time peak since 2022, suggesting timely intervention to curtail the spread of the mpoxv outbreak in Nigeria.

4.2. Numerical Simulations

In this part, we analysed the model numerically to gain a profound understanding of the dynamic behaviour of the model. The classical model (1) is numerically simulated in this part using the biological parameters specified in Table 2. In order to solve the nonlinear model depicted in the model (1). The numerical findings are presented using first-order convergent numerical techniques. The numerical method used to solve differential equations of integer order is precise, conditionally stable, and convergent. Based on factors that are treated as variables, we were able to provide the graphical results.

We conducted numerical simulations to examine some key model parameters and evaluate the dynamics of each compartment. We examined the dynamics of the model compartments and found that several vital parameters substantially impact whether the mpoxv increases or decreases. These variables can also provide insight into how to lessen transmissions, such as by reducing the dangers to the most susceptible groups and the number of sick people. Additionally, understanding these quantities helps us better understand the contagiousness of a particular infection and how changes in the environment and certain interventions can affect it. In order to achieve this, the individual compartments’ dynamic features are represented in Figure 3. Simulation results showing the effect of the rate of infectiousness in humans $\theta$ moving from $I_1$ to $I_2$ are depicted in Figure 4a,b with some varying values, which shows increasing and decreasing patterns at some points. Moreover, the model was simulated to show the impact of hospitalisation rate parameters $k_1$ on $I_1$ and H, and $k_2$ on $I_1$ and H in Figure 4c–f, respectively. The transmission rates and probability are simulated with some varying values on the compartments $S_r,E_r,I_r$, in Figure 4g–i, respectively. The behaviour of the rate at which mpoxv progresses, $\omega$, from P to $I_1$, is shown n the compartments P to $I_1$ in Figure 4i,j.

4.3. Sensitivity Analysis

In this subsection, sensitivity analysis is scrutinised to determine the powerful impact and influence of various parameters on mpoxv transmission in Nigeria. The sensitivity analysis of (1) with ${\mathcal{R}}_0$ and infectious attack rate as response functions were revealed using the partial rank correlation coefficient (PRCC); see Figure 5. This technique has been widely used previously [46,47]. Our study outcomes imply ${\beta }_{hh}$ and $\omega$ parameters are the most sensitive parameters of the model, necessitating close monitoring to reduce mpoxv transmission in Nigeria and other countries with similar settings. The diagrammatic representation of the PRCC in terms of ${\mathcal{R}}_0$ and infectious attack rate is depicted in Figure 5. The parameter values in Table 1 are employed to perform the job.

5. Concluding Remarks

The current spread of mpoxv in endemic and non-endemic countries appears to be strikingly different from previous epidemics, with a disproportionate number of cases in MSM [2,12]. Mpoxv is primarily transmitted from animals to humans in most endemic countries, such as Nigeria [2,11,23]. In this research, we proposed an epidemiological modelling approach to examine the transmission of mpoxv considering high-and low-risk susceptible populations. The model also considered the prodromal phase, mild and severe infection stages for mpoxv, as well as the role of hospitalisation to assess the prevailing mpoxv infection in Nigeria. The model fitted well to the Nigerian case scenarios to assess the current mpoxv situation, which would enhance mpoxv control, especially in West Africa, where the disease has been endemic for a long time.

We rigorously analysed the model mathematically and obtained some crucial threshold dynamics, such as the ${\mathcal{R}}_0$, which can be used to assess the effect of mpoxv control strategies. Subsequent theoretical analysis also revealed that the DFE of the model is GAS when ${\mathcal{R}}_0$ is less than or equal to unity and unstable when it is greater than 1. We also found that the EE is GAS when ${\mathcal{R}}_0$ is greater than 1, indicating the potential for mpoxv to spread in a community.

Moreover, the plausibility of the occurrence of bifurcation in the current model was investigated theoretically, and we found that forward bifurcation likely exists (see Section 3.4 for more detail). Consequently, forward bifurcation in the present model indicates the DFE and EE exchange stability at ${\mathcal{R}}_0=1$. That is, forward bifurcation exists in the current model at ${\mathcal{R}}_0=1$, and an endemic equilibrium exists (which is GAS under certain conditions as stated in Theorem 2) when ${\mathcal{R}}_0>1$. By epidemiological implication, the infection will persist in a community if the average number of newly infected individuals is greater than 1. Hence, we observed that the stability behavior of the current model changes (from stable to unstable) around the endemic equilibrium when the bifurcation parameter ${\beta }_{hh}$ varies [48,49].

Results from our simulations suggest that the risk of a mpoxv outbreak remains high, especially in high-risk areas without intervention. This includes rural communities where the interaction of humans and rodents is more prevalent [7,17]. In particular, the broader populace could be at risk if transmission rates rise among high-risk groups. In the event of such a situation brought on by virus evolution, isolation of symptomatic cases may not be sufficient to stop the spread. Case detection and contact tracking in high-risk groups would become crucial techniques to avoid bridging from high-risk to low-risk groups. Additionally, areas with a higher percentage of MSM may face greater dangers and need stricter public health measures, such as contact tracing and isolation. Monitoring the prevalence of mpoxv in animals is essential to indicate the likelihood of mpoxv outbreaks since the virus’ transmission is tightly linked to seasonal fluctuations in animal host abundance and activity level [5].

We also investigated sensitivity analysis using the partial rank correlation coefficient to determine the powerful impact and influence of various model parameters on mpoxv transmission in Nigeria. Our sensitivity analysis results (with ${\mathcal{R}}_0^h$ and infection attack rate as response functions) show that the parameters ${\beta }_{hh}$ and $\omega$ are the most sensitive parameters. Indicating that reducing mpoxv transmission in Nigeria and other nations with similar settings requires constant monitoring as well as awareness campaigns, especially in poor resource areas.

Our findings also provided a deeper understanding of the trends and mechanisms that influence mpoxv spread in Nigeria. Environmental sanitisation, particularly in high-risk locations, is required since mpoxv is a zoonotic illness with humans living in low-resource settings being the most vulnerable. This will effectively minimise the morbidity and mortality of the disease. As the mpoxv transmission is presumably driven by ecological factors with rodents as the principal host; thus, we plan to investigate the effect of seasonality on the overall mpoxv transmission. We suggest that to earn sufficient mpoxv prevention and mitigation, extensive research on mpoxv, the supply of enough medical resources, and awareness programs are essential. Finally, our model can be extended to incorporate seasonality in order to investigate further the influence of seasonality and rodent-driven factors for the transmission of mpoxv in Nigeria.