Comparative Cost-Effectiveness Analysis of Awareness and Treatment Strategies for Controlling Mpox in Diverse Epidemiological Settings

1. Introduction

Mpox, previously known as monkeypox, is a zoonotic viral infection caused by the monkeypox virus (MPXV), classified under the genus Orthopoxvirus in the Poxviridae family [1,2]. First identified in 1958 in laboratory monkeys in Copenhagen, Denmark, the disease was recognized as a human health concern in 1970, following the first confirmed human case in the Democratic Republic of Congo (DRC) [1,2,3]. In 2022, the World Health Organization (WHO) officially renamed the disease Mpox to mitigate the stigmatization associated with its former name [4]. Despite its classification as a neglected tropical disease, the resurgence of Mpox cases beyond endemic regions underscores the urgent need for targeted public health interventions, particularly in resource-limited settings where access to healthcare remains a challenge.

Epidemiologically, mpox remains endemic in several African countries, particularly in central and western Africa, where sporadic outbreaks have been documented since its discovery. The virus is characterized by two genetically distinct clades: Clade I (Congo Basin) and Clade II (West African). Clade I (with subclades Ia and Ib) is associated with a higher fatality rate for cases of approximately 10% and increased transmissibility, making it a critical public health concern in regions with limited healthcare infrastructure [2,5,6]. Clade II (with subclades IIa and IIb), specifically Clade IIb, is responsible for the global outbreak of 2022, exhibits a lower virulence but remains a significant epidemiological threat due to its potential for human-to-human transmission, particularly in densely populated urban centers [2,5,6]. These clades have spread significantly to Europe and North America; however, the threat posed by clade I remains low, while clade II mpox is still circulating at low levels.

The re-emergence of mpox beyond its endemic regions serves as a stark reminder of the global interconnectedness of infectious diseases. It underscores the importance of integrated surveillance systems, rapid response frameworks, and international cooperation to prevent and manage emerging zoonotic threats effectively. The 2022 outbreak marked a pivotal moment in mpox epidemiology, as the disease spread rapidly beyond its endemic regions to more than 100 countries, including Europe, the Americas, and Asia. This unprecedented spread highlighted the potential for global dissemination of the virus, driven by increased international travel, urbanization, and possibly asymptomatic transmission. The outbreak also underscored the urgent need for comprehensive surveillance, rapid diagnostic testing, and targeted public health interventions to curb the spread, particularly in settings with limited resources [7,8,9].

Moreover, mathematical modeling studies have been instrumental in elucidating the transmission dynamics of Mpox and identifying optimal control strategies. Ngungu et al. [10] developed a mathematical epidemiological model using real data from the United Kingdom to assess the impact of non-pharmaceutical interventions on the transmission dynamics of Mpox, highlighting the effectiveness of targeted control measures in mitigating infection spread. Furthermore, the dynamics of mpox transmission have been extensively studied through some novel modeling approaches, emphasizing the role of contaminated surfaces, asymptomatic infection, and co-infections with other diseases such as COVID-19. For instance, Hassan et al. [11] developed a deterministic model incorporating environmental transmission pathways and contaminated surfaces to emphasize the role of fomite transmission in sustaining outbreaks and to assess the stability of the disease-free and endemic equilibria. Their findings highlight that while effective control measures can reduce disease prevalence, the persistence of contaminated surfaces remains a significant risk factor for Mpox transmission. Similarly, Li et al. [7] assessed the impact of asymptomatic infection on disease spread in Nigeria, demonstrating that personal protection coupled with effective vaccination is crucial for mitigating outbreaks. Co-infection scenarios involving Mpox and other viral diseases such as COVID-19 have also been examined. Acheneje et al. [12] analyzed the concurrent transmission dynamics of both diseases, incorporating optimal control strategies aimed at reducing co-infection rates through social distancing, rodenticides, and vaccination. Their cost–benefit analysis underscores the importance of integrated control measures in minimizing both infection rates and associated costs.

Environmental factors have also been incorporated into mpox transmission models. Alshehri and Ullah [9] emphasized the impact of environmental viral concentration, showing that environmental contamination can significantly influence disease dynamics, particularly in densely populated urban settings. Their model recommends implementing optimal control strategies targeting environmental viral loads to effectively contain outbreaks.

Recent studies have further explored the application of optimal control theory to Mpox transmission dynamics and have assessed the effectiveness of intervention strategies such as vaccination, personal protective equipment, and isolation measures, underscoring the importance of context-specific control measures tailored to different epidemiological settings. Adepoju and Ibrahim [13] focused on vaccination and immunity loss following recovery, employing Pontryagin’s maximum principle to determine cost-effective control strategies. Rashid et al. [8] extended these efforts by developing a hybrid fractional Mpox model that integrates external factors, illustrating the benefits of targeted interventions in reducing the basic reproduction number and associated healthcare costs. Building on their team’s previous work on mathematical modeling approaches to investigate the transmission dynamics of the Mpox virus [14], Peter et al. [15] conducted a comparative analysis of cost-effective strategies for Mpox control, identifying the optimal combination of treatment and public awareness campaigns. Their findings indicate that while combined strategies are effective in reducing disease prevalence, targeted treatment-only interventions exhibit lower incremental cost-effectiveness ratios (ICERs), particularly in resource-limited settings like Nigeria and the DRC.

Building on these foundational studies and our previous work in [16], the present study develops and analyzes an optimal control model for Mpox transmission dynamics across four countries: Nigeria, Spain, Italy, and the DRC. By implementing a compartmental model and Pontryagin’s Maximum Principle, we aim to assess the cost-effectiveness of awareness, treatment, and combined control strategies through numerical simulations. The efficiency indices (EI), total cost (TC), and ICER for each intervention are evaluated, providing insights into the optimal public health strategies for mitigating Mpox outbreaks in diverse epidemiological contexts.

The remainder of this study is organized as follows: In Section 2, we develop a mathematical model to analyze the transmission dynamics of Mpox, extending it to include optimal control strategies focused on awareness, treatment, and a combination of both interventions. This section also derives the control reproduction number and characterizes the optimal control strategies using Pontryagin’s Maximum Principle. Section 3 presents numerical simulations to assess the impact of the proposed control strategies across four countries — Nigeria, Spain, Italy, and the Democratic Republic of Congo — evaluating the effectiveness of awareness-only, treatment-only, and combined control measures in reducing Mpox transmission. In Section 4, a cost-effectiveness analysis is conducted, involving the computation of Efficiency Indices (EI), Total Cost (TC), and Incremental Cost-Effectiveness Ratios (ICER) for each intervention scenario, identifying the most economically viable strategies in both resource-rich and resource-limited settings. Finally, in Section 5, we discuss the key findings, highlighting their implications for public health policy and identifying potential avenues for future research, particularly in optimizing intervention strategies for emerging zoonotic diseases in settings with limited healthcare resources.

2. Mathematical Model Formulation

The model we developed in this article takes into consideration the two strains (Clade I and Clade II) of the Mpox virus, which are based on two different infection pathways in the human population of the model. The total human population at time t, denoted by $N\left(t\right)$, is subdivided into susceptible humans $S\left(t\right)$, exposed humans $E\left(t\right)$, infectious humans based on Clade I $I_1\left(t\right)$, infectious humans based on Clade II $I_2\left(t\right)$, hospitalized humans $H\left(t\right)$, and recovered humans $R\left(t\right)$. So, the total human population is given by

$$N\left(t\right)=S\left(t\right)+E\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+H\left(t\right)+R\left(t\right),$$

(1)

while the mammal population is sub-divided into two compartments, susceptible $S_m\left(t\right)$ and infectious $I_m\left(t\right)$, such that the total mammal population is given by

$$M\left(t\right)=S_m\left(t\right)+I_m\left(t\right).$$

(2)

The total population of both human and mammal populations considered in the model dynamics is given as

$$N_T\left(t\right)=N\left(t\right)+M\left(t\right).$$

(3)

Thus, the transmission model with two strain infection dynamics is represented by the following system of non-linear ordinary differential equations:

$$\begin{array}{cc}\hfill \dot{S}& ={\Pi }_h-({\lambda }_h+{\mu }_h)S+ϵR\hfill \\ \hfill \dot{E}& ={\lambda }_{h}S-\gamma E_h-{\mu }_h{E}_h\hfill \\ \hfill \dot{I_1}& =(1-\theta )\gamma E_h-({\tau }_1+{\delta }_1+{\mu }_h)I_1\hfill \\ \hfill \dot{I_2}& =\theta \gamma E_h-({\tau }_2+{\delta }_2+{\mu }_h+\psi )I_2\hfill \\ \hfill \dot{H}& ={\tau }_1{I}_1+{\tau }_2{I}_2-({\delta }_3+{\mu }_h+\eta )H\hfill \\ \hfill \dot{R}& =\eta H+\psi I_2-(ϵ+{\mu }_h)R\hfill \\ \hfill \dot{S_m}& ={\Pi }_m-{\lambda }_m{S}_m-{\mu }_m{S}_m\hfill \\ \hfill \dot{I_m}& ={\lambda }_m{S}_m-{\mu }_m{I}_m\hfill \end{array}$$

(4)

under the following initial conditions:

$$S\left(0\right)>0,E\left(0\right)\ge 0,I_1\left(0\right)\ge 0,I_2\left(0\right)\ge 0,H\left(0\right)\ge 0,R\ge 0,S_m>0,I_m\ge 0.$$

(5)

The effective contact rate with an infectious human or infectious mammal is defined as:

$${\lambda }_h={\displaystyle \frac{{\beta }_h(I_1+\sigma I_2)}{N_h}},{\lambda }_m={\displaystyle \frac{\alpha {\beta }_{hm}(I_1+\sigma I_2)}{N_h}}+{\displaystyle \frac{{\beta }_m{I}_m}{N_m}}.$$

(6)

The schematic diagram of the model is presented in Figure 1 and the state variables and parameters of the model are described in Table 1. In the appendix, we present some mathematical analysis of model (4).

2.1. Development of Optimal Control Model

In this section, we aim to develop an optimal control model designed to reduce the incidence of Mpox in the selected countries Italy, Spain, the Democratic Republic of the Congo (DRC), and Nigeria. This model incorporates three distinct control variables:

• Awareness Campaign ${\mathbf{u}}_{\mathbf{1}}$: This control represents initiatives aimed at preventing Mpox through educational measures. It leverages social media, mainstream media, and civil society organizations, including religious groups, to disseminate information and promote preventive behaviors.
• Effective Treatment and Management ${\mathbf{u}}_{\mathbf{2}}$: This control focuses on specific antiviral treatments for Mpox, including tecovirimat (TPOXX), brincidofovir, and cidofovir, to improve patient outcomes and reduce disease transmission.

These control measures are integrated into the model 4 to optimize the mitigation strategies of Mpox in selected countries.

$$\begin{array}{cc}\hfill \dot{S}& ={\Pi }_h-(1-u_1){\lambda }_{h}S+{\mu }_{h}S+ϵR\hfill \\ \hfill \dot{E}& =(1-u_1){\lambda }_{h}S-\gamma E_h-{\mu }_h{E}_h\hfill \\ \hfill \dot{I_1}& =(1-\theta )\gamma E_h-({\tau }_1+{\delta }_1+{\mu }_h)I_1\hfill \\ \hfill \dot{I_2}& =\theta \gamma E_h-({\tau }_2+{\delta }_2+{\mu }_h+\psi )I_2\hfill \\ \hfill \dot{H}& ={\tau }_1{I}_1+{\tau }_2{I}_2-({\delta }_3+{\mu }_h)H-(\eta +u_2)H\hfill \\ \hfill \dot{R}& =(u_2+\eta )H+\psi I_2-(ϵ+{\mu }_h)R\hfill \\ \hfill \dot{S_m}& ={\Pi }_m-(1-u_3){\lambda }_m{S}_m-{\mu }_m{S}_m\hfill \\ \hfill \dot{I_m}& =(1-u_3){\lambda }_m{S}_m-{\mu }_m{I}_m\hfill \end{array}$$

(7)

under the following initial conditions:

$$S\left(0\right)>0,E\left(0\right)\ge 0,I_1\left(0\right)\ge 0,I_2\left(0\right)\ge 0,H\left(0\right)\ge 0,R\ge 0,S_m>0,I_m\ge 0.$$

(8)

The goal is to reduce the number of infected individuals for both Clade I and Clade II while keeping intervention costs minimal. The control measures include $u_1\left(t\right)$ for awareness campaigns, $u_2\left(t\right)$ for effective treatment. Mathematically, this objective is achieved by minimizing the following cost function, formulated as a quadratic function of the control variables to align with established standards in the literature [17,18,19,20]. Thus, the model’s objective functional is given by:

$$\mathcal{J}(u_1,u_2)={\int }_{t_0}^{t_f}\left(w_{1}E+w_2{I}_1+w_3{I}_2+w_{4}H+{\displaystyle \frac{1}{2}}\sum _{i=1}^2{\varphi }_i{u}_i^2\right).dt,$$

(9)

and

$$\mathcal{J}(u_1^{*},u_2^{*})=min\mathcal{J}(u_1,u_2):u_1,u_2\in \mathbb{U},$$

(10)

Thus, the control set is defined as

$$\mathbb{U}=\{u_1\left(t\right),u_2\left(t\right)\in {\mathbb{R}}^2\}$$

where $u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)$ are Lebesgue measurable functions satisfying

$$0\le u_1\left(t\right),u_2\left(t\right)\le 1,\mathrm{for}0\le t.$$

From Equation (9), $t_0$ and $t_f$ denote the initial and final time, respectively. The parameters $w_1,w_2,w_3,$ and $w_4$ represent weight constants corresponding to the exposed individuals, individuals infected with strain I, individuals infected with strain II, and hospitalized individuals, respectively. These weights reflect the relative importance of reducing specific population groups to effectively mitigate the spread of Mpox in selected countries. Additionally, ${\varphi }_i$ (for $i=1, 2$) are weights designed to capture the relative costs or efforts associated with implementing each time-dependent control strategy.

For the purpose of this study, we assume hypothetical values for the upper limits of the control functions:

$$u_{1max}=1,u_{2max}=1$$

However, these values are selected based on the assumption that they are realistically achievable given the available resources.

2.2. Control Reproduction Number

We establish a threshold parameter to control Mpox spread in selected countries using the next-generation matrix method. This approach is widely recognized and reliable, as demonstrated in related works [3,19,21,22]. To achieve this, model Equation (11) is simplified to focus on the compartments relevant to disease transmission.

$$\begin{array}{cc}\hfill \dot{E}& =(1-u_1){\lambda }_{h}S-\gamma E_h-{\mu }_h{E}_h\hfill \\ \hfill \dot{I_1}& =(1-\theta )\gamma E_h-({\tau }_1+{\delta }_1+{\mu }_h)I_1,\hfill \\ \hfill \dot{I_2}& =\theta \gamma E_h-({\tau }_2+{\delta }_2+{\mu }_h+\psi )I_2,\hfill \\ \hfill \dot{I_m}& ={\lambda }_m{S}_m-{\mu }_m{I}_m,\hfill \end{array}$$

(11)

$$F=\left[\begin{array}{cccc}0& (1-u_1){\beta }_h& (1-u_1){\beta }_h\sigma & 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& {\displaystyle \frac{{\beta }_{hm}\alpha {\mu }_h{\pi }_m}{{\pi }_h{\mu }_m}}& {\displaystyle \frac{{\beta }_{hm}\alpha \sigma {\mu }_h{\pi }_m}{{\pi }_h{\mu }_m}}& {\beta }_m\end{array}\right];V^{-1}=\left[\begin{array}{cccc}\gamma +{\mu }_h& 0& 0& 0\\ (-1+\theta )\gamma & k_1& 0& 0\\ -\theta \gamma & 0& k_2+\psi & 0\\ 0& 0& 0& {\mu }_m\end{array}\right]$$

for simplification ${\kappa }_1=({\delta }_1+{\mu }_h+{\tau }_1),{\kappa }_2=({\delta }_2+{\mu }_h+{\tau }_2),$ while the control reproduction number is computed as:

$$\begin{array}{cc}\hfill \rho \left(F{V}^{-1}\right)={\mathcal{R}}_c& ={\displaystyle \frac{(1-u_1){\beta }_h\gamma \left[{\kappa }_1\sigma \theta +\left(1-\theta \right)\left({\kappa }_2+\psi \right)\right]}{\left(\gamma +{\mu }_h\right){\kappa }_1\left({\kappa }_2+\psi \right)}}\hfill \end{array}$$

(12)

in the absence of control measure to mitigate the spread of Mpox in each of the countries, the control reproduction number computed becomes the basic reproduction number computed as:

$$\begin{array}{cc}\hfill {\mathcal{R}}_0& ={\displaystyle \frac{{\beta }_h\gamma \left[{\kappa }_1\sigma \theta +\left(1-\theta \right)\left({\kappa }_2+\psi \right)\right]}{\left(\gamma +{\mu }_h\right){\kappa }_1\left({\kappa }_2+\psi \right)}}\hfill \end{array}$$

(13)

In Figure 2, the effect of the control measure on the control reproduction number in Equation (12) is illustrated by plotting the transmission rate (${\beta }_h$) against the control measure ($u_1$) for selected countries. The visualization shows that as the transmission rate decreases and disease awareness increases, the reproduction number gradually declines over time.

2.3. Characterization of the optimal control

The optimal control variables $(u_1^{*},u_2^{*})$, associated with the three intervention strategies, along with their corresponding state variables $E^{*},I_1^{*},I_2^{*},H^{*}$, are determined by applying Pontryagin’s Maximum Principle [17,23]. The necessary conditions for optimality are derived accordingly, while the terminal conditions for the adjoint variables are specified using the transversality conditions as described by Hartl [18]. Thus, the Hamiltonian $\mathbb{H}$ is given by:

$$\begin{array}{cc}\hfill \mathbb{H}=& w_{1}E+w_2{I}_1+w_3{I}_2+w_{4}H+{\displaystyle \frac{1}{2}}({\varphi }_1{u}_1^2+{\varphi }_2{u}_2^2+{\varphi }_3{u}_3^2)\hfill \\ \hfill & +{\chi }_1\left[{\Pi }_h-(1-u_1)({\lambda }_h+{\mu }_h)S+ϵR\right]\hfill \\ \hfill & +{\chi }_2\left[(1-u_1){\lambda }_{h}S-\gamma E_h-{\mu }_h{E}_h\right]\hfill \\ \hfill & +{\chi }_3\left[(1-\theta )\gamma E_h-({\tau }_1+{\delta }_1+{\mu }_h)I_1\right]\hfill \\ \hfill & +{\chi }_4\left[\theta \gamma E_h-({\tau }_2+{\delta }_2+{\mu }_h+\psi )I_2\right]\hfill \\ \hfill & +{\chi }_5\left[{\tau }_1{I}_1+{\tau }_2{I}_2-({\delta }_3+{\mu }_h+\eta +u_2)H\right]\hfill \\ \hfill & +{\chi }_6\left[(\eta +u_2)H+\psi I_2-(ϵ+{\mu }_h)R\right]\hfill \\ \hfill & +{\chi }_7\left[{\Pi }_m-{\lambda }_m{S}_m-{\mu }_m{S}_m,\right]\hfill \\ \hfill & +{\chi }_8\left[{\lambda }_m{S}_m-{\mu }_m{I}_m\right]\hfill \end{array}$$

(14)

where ${\chi }_i,i=1,\dots ,8$ are the adjoint variables. Thus, the Pontryagins maximum principle [18] is applied.

The optimal solution for $u^{*}=(u_1,u_2)$ can be determined by applying the necessary conditions established by Pontryagin’s principle. This process entails differentiating $\mathbb{H}$ with respect to each control variable and solving the resulting equations. Through this approach, we obtain the expressions for the optimal control variables $u^{*}=(u_1^{*},u_2^{*})$ as follows:

$${\displaystyle \frac{\partial H}{\partial u_1}}=0,{\displaystyle \frac{\partial H}{\partial u_2}}=0$$

Then the following optimal conditions for $(u_1^{*},u_2^{*})$ are obtained as:

$$\begin{array}{cc}\hfill u_1=& {\displaystyle \frac{S\left[{\beta }_h(\sigma I_2+I_1)({\chi }_2-{\chi }_1)-{\mu }_{h}N{\chi }_1\right]}{N{\varphi }_1}}\hfill \\ \hfill u_2=& {\displaystyle \frac{H({\chi }_5-{\chi }_6)}{{\varphi }_2}}\hfill \end{array}$$

(15)

ensuring the optimal values obtained in Equation (15) is always in the interval $[0,1]$ the optimal value for each control is then defined as follows:

$$\begin{array}{cc}\hfill u_1^{*}=& min\left\{max\left(0,{\displaystyle \frac{S\left[{\beta }_h(\sigma I_2+I_1)({\chi }_2-{\chi }_1)-{\mu }_{h}N{\chi }_1\right]}{N{\varphi }_1}}\right),{u_1}_{max}\right\},\hfill \\ \hfill u_2^{*}=& min\left\{max\left(0,{\displaystyle \frac{H({\chi }_5-{\chi }_6)}{{\varphi }_2}}\right),{u_2}_{max}\right\}\hfill \end{array}$$

(16)

By substituting $u_1^{*},u_2^{*}$ into Equation (14), the optimal condition for the Hamiltonian $H^{*}$ is obtained. Furthermore, the adjoint system is derived by computing the partial derivatives of $\mathbb{H}$ with respect to each state variable

$${\displaystyle \frac{d{\chi }^{*}\left(t\right)}{dt}}=-{\displaystyle \frac{\partial H^{*}}{\partial x}}$$

leading to the following expressions:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_1^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial S}}=& -{\chi }_1\left[(1-u_1){\beta }_h{\displaystyle \frac{(I_2\sigma +I_1)S}{N^2}}-(1-u_1)\left({\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)}{N}}+{\mu }_h\right)\right]\hfill \\ \hfill & -{\chi }_2\left[(1-u_1){\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)}{N}}-(1-u_1){\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)S}{N^2}}\right]\hfill \\ & -{\chi }_7{\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)S_m}{N^2}}+{\chi }_8{\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)S_m}{N^2}}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_2^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial E}}=& -w_1-{\displaystyle \frac{{\chi }_1(1-u_1){\beta }_h(I_2\sigma +I_1)S}{N^2}}-{\chi }_2\left[-(1-u_1){\beta }_h{\displaystyle \frac{(I_2\sigma +I_1)S}{N^2}}-\gamma -{\mu }_h\right]\hfill \\ & -{\chi }_3(1-\theta )\gamma -{\chi }_4\theta \gamma -{\chi }_7\alpha {\beta }_{hm}{\displaystyle \frac{(I_2\sigma +I_1)S_m}{N^2}}\hfill \\ & +{\chi }_8\alpha {\beta }_{hm}{\displaystyle \frac{(I_2\sigma +I_1)S_m}{N^2}}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_3^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial I_1}}=& -w_2+{\chi }_1(1-u_1)\left({\displaystyle \frac{{\beta }_h}{N}}-{\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)}{N^2}}\right)S\hfill \\ & -{\chi }_2\left[{\displaystyle \frac{(1-u_1)S{\beta }_h}{N}}-(1-u_1){\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)S}{N^2}}\right]-{\chi }_3(-r_1-d_1-{\mu }_h)\hfill \\ & -{\chi }_5{r}_1+{\chi }_7\left[{\displaystyle \frac{\alpha {\beta }_{hm}}{N}}-{\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)}{N^2}}\right]S_m\hfill \\ & -{\chi }_8\left[{\displaystyle \frac{\alpha {\beta }_{hm}}{N}}-{\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)}{N^2}}\right]S_m\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_4^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial I_2}}=& -w_3+{\chi }_1(1-u_1)\left({\displaystyle \frac{{\beta }_h\sigma }{N}}-{\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)}{N^2}}\right)S\hfill \\ & -{\chi }_2(1-u_1)\left[{\displaystyle \frac{S{\beta }_h\sigma }{N}}-{\displaystyle \frac{{\beta }_h(I_2\sigma +I_1)S}{N^2}}\right]-{\chi }_4(-r_2-d_2-{\mu }_h-\psi )\hfill \\ & -{\chi }_5{r}_2-{\chi }_6\psi +{\chi }_7\left[{\displaystyle \frac{\alpha {\beta }_{hm}\sigma }{N}}-{\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)}{N^2}}\right]S_m\hfill \\ & -{\chi }_8\left[{\displaystyle \frac{\alpha {\beta }_{hm}\sigma }{N}}-{\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)}{N^2}}\right]S_m\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_5^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial H}}=& -w_4-{\chi }_1(1-u_1){\beta }_h{\displaystyle \frac{(I_2\sigma +I_1)S}{N^2}}+{\chi }_2(1-u_1){\beta }_h{\displaystyle \frac{(I_2\sigma +I_1)S}{N^2}}\hfill \\ & -{\chi }_5(-d_3-{\mu }_h-\eta -u_2)-{\chi }_6(\eta +u_2)\hfill \\ & -{\chi }_7(1-u_3){\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)S_m}{N^2}}+{\chi }_8(1-u_3){\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)S_m}{N^2}}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_6^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial R}}=& -{\chi }_1\left[{\displaystyle \frac{(1-u_1){\beta }_h(I_2\sigma +I_1)S}{N^2}}+ϵ\right]+{\chi }_2{\displaystyle \frac{(1-u_1){\beta }_h(I_2\sigma +I_1)S}{N^2}}\hfill \\ & -{\chi }_6(-ϵ-{\mu }_h)-{\chi }_7(1-u_3){\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)S_m}{N^2}}\hfill \\ & +{\chi }_8(1-u_3){\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)S_m}{N^2}}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_7^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial S_m}}=& -{\chi }_7\left[{\displaystyle \frac{(1-u_3){\beta }_m{I}_m{S}_m}{{(S_m+I_m)}^2}}-(1-u_3)\left({\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)}{N}}+{\displaystyle \frac{{\beta }_m{I}_m}{S_m+I_m}}\right)-{\mu }_m\right]\hfill \\ & -{\chi }_8\left[-{\displaystyle \frac{(1-u_3){\beta }_m{I}_m{S}_m}{{(S_m+I_m)}^2}}+(1-u_3)\left({\displaystyle \frac{\alpha {\beta }_{hm}(I_2\sigma +I_1)}{N}}+{\displaystyle \frac{{\beta }_m{I}_m}{S_m+I_m}}\right)\right]\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\chi }_8^{*}}{dt}}=-{\displaystyle \frac{\partial \mathbb{H}}{\partial I_m}}=& {\chi }_7(1-u_3)\left[{\displaystyle \frac{{\beta }_m}{S_m+I_m}}-{\displaystyle \frac{{\beta }_m{I}_m}{{(S_m+I_m)}^2}}\right]S_m\hfill \\ & -{\chi }_8\left[(1-u_3)\left({\displaystyle \frac{{\beta }_m}{S_m+I_m}}-{\displaystyle \frac{{\beta }_m{I}_m}{{(S_m+I_m)}^2}}\right)S_m-{\mu }_m\right]\hfill \end{array}$$

(24)

With the transversal condition that ${\chi }_i\left(t_f\right)=0$ for $i=1,2,\dots ,8$

3. Numerical analysis

This section explores the impact of the two optimal control strategies in mitigating the spread of monkeypox within the human population. The primary focus is on evaluating the effectiveness of combined control measures in halting disease transmission. To achieve this, we conduct a numerical simulation of the monkeypox model, comparing scenarios with and without optimized control interventions to assess the influence of the control variables introduced earlier. However, the effects of different control strategy are represented by $u_1,$ and $u_2$. These strategies are systematically categorized into two groups: single control, and combine controls allowing for a structured analysis of the three possible control strategies considered in this research. Specifically, the strategies are defined as follows:

• Strategy 1: Implements only awareness strategy $u_1$.
• Strategy 2: Implements only recommended effective treatment strategy $u_2$.
• Strategy 3: Implements both awareness and recommended treatment strategies combine $u_1$ and $u_2$.

The model simulations utilize parameter values from existing literature [16] for selected countries, The optimality system integrates the model equations, adjoint equations, and control characterizations, solved in MATLAB R2022B using the forward-backward sweep method. This approach involves a 16-dimensional system of ordinary differential equations solved iteratively with the fourth-order Runge-Kutta method. Paramater values use for the numerical simulation are presented in Table 2 and Table 3 with a final time horizon of 100 days. The weight coefficients and cost functional are assigned as ($w_1=40,w_2=30,w_3=10,w_4=30,{\varphi }_1=10,000,{\varphi }_2=80,000$).

Here, we clarify the choice to choose our cost. The higher expenditure associated with awareness efforts, as compared to treatment and zoonotic control, is largely due to the extensive scope and continuity required for effective implementation. Public education campaigns often involve repeated communication through various channels such as radio, television, social media, and community outreach programs. These activities demand ongoing investments in logistics, personnel, training, and materials to ensure the message reaches diverse and widespread populations.While treatment and zoonotic control typically involve more direct, targeted actions such as administering medication to confirmed cases or managing animal reservoirs in specific areas which can be executed with more contained costs. As a result, promoting awareness on a large scale tends to incur greater overall financial commitment.

Thus, we assume that ${\varphi }_1>{\varphi }_2$, indicating that the cost of awareness is higher than the cost of treatment and zoonotic control.

Model (4) parameter values are cited from the literature for Spain, Italy, Nigeria, and the Democratic Republic of the Congo. The parameters $({\Pi }_m,{\delta }_1,{\delta }_1,{\delta }_3)$ are hypothetically chosen, while others are taken from the literature: ${\Pi }_h$ [24], $\mu$ [1], and $\gamma$ [2,3,25].

3.1. Optima Control Simulation

In the following sections, a detailed simulation is presented, along with an in-depth analysis and discussion of intervention strategies.

As depicted in Figure 3a Nigeria, with its large and densely populated regions, faces increased risk of rapid Mpox transmission, especially in urban centers where close contact is common. The No control scenario (in solid red line), which shows the highest peak and slowest decline, reflects what could happen in the absence of intervention widespread transmission and prolonged outbreaks. The Treatment Only control $u_2$ and Awareness Control Only $u_1$ strategies offer moderate benefits, but alone they are insufficient to quickly suppress Mpox outbreaks, likely due to logistical challenges to access healthcare care and the reach of treatment. Interestingly, treatment alone performs better than awareness, suggesting that strengthening clinical response infrastructure may yield faster short term impact in Nigeria. The most effective result comes from the combination of Treatment and Awareness (green line), which achieves the fastest decline in Mpox disease prevalence. This aligns with Nigeria’s need for comprehensive strategies that include not just clinical treatment, but strong public awareness campaigns to reduce risky behavior, encourage early reporting, and counter misinformation.

In light of the current Mpox situation in Italy characterized by low but ongoing reported cases, the control plot depicted in Figure 4a highlights the comparative effectiveness of various intervention strategies. The simulation indicates that while treatment alone moderately reduces disease burden, public awareness campaigns (black solid line) have a stronger initial impact by limiting exposure and transmission. Notably, the combination of treatment and awareness achieves the most substantial and rapid decline in Mpox prevalence. This outcome supports the implementation of integrated public health measures in Italy, reinforcing the need for both medical readiness and sustained risk communication to prevent a potential resurgence.

The optimal control analysis for the Spain population as depicted in Figure 5a reveals that combining treatment and public awareness significantly reduces Mpox disease prevalence more effectively than either strategy alone. Without any control measures (dotted line in red), the disease peaks higher and declines slowly, indicating a severe outbreak. Implementing treatment alone moderately reduces prevalence of Mpox disease in spain, but awareness campaigns alone outperform treatment by promoting early prevention. However, the synergy of both controls results in the fastest and steepest decline in cases, achieving near-eradication within 60 days. Therefore, integrated intervention is the most effective strategy for controlling Mpox in Spain. Furthermore, In Figure 5b both awareness and treatment controls were used at full strength for most of the 100 days. This helped reduce Mpox cases quickly and effectively. The controls dropped near the end as the disease was brought under control.

In Figure 6a, the plot illustrates the effectiveness of various control strategies on Mpox prevalence in the Democratic Republic of the Congo (DRC). The red curve, representing no control, shows a sustained high burden of disease, reflecting the current challenges faced by the DRC, such as limited healthcare access and weak surveillance systems. The black curve (awareness only) offers a better improvement in the prevalence of the Mpox disease, indicating the importance of public education and publicity alone has positive impact in curtailing the disease. The blue dashed line (treatment only) demonstrates slight improvement, emphasizing the importance of timely medical intervention. However, the green dotted curve (treatment and awareness) shows the most significant reduction in Mpox cases, highlighting that a combined strategy is the most effective. This supports the need for integrated policies that improve both public health education and access to treatment in DRC.

4. Cost Effectiveness Analysis

Evaluating the cost dynamics of monkeypox management (Mpox) in this study involves analyzing two main controls: awareness (prevention) costs and treatment (curative) costs. Both are integral to a country’s public health strategy, but they have different structures and financial implications in selected countries. Awareness campaigns are primarily focused on prevention reducing the spread of the disease through education and behavior change. On the other hand, treatment costs are incurred when individuals are diagnosed and require medical intervention. These costs can be significant, especially for moderate to severe cases. The cost of hospitalization includes expenses related to inpatient care, bed space, feeding, and utilities. Symptomatic treatments, such as antipyretics, antivirals, and antibiotics to manage secondary bacterial infections, contribute to the overall cost of treatment.

In countries like Nigeria, Spain, Italy, and the Democratic Republic of Congo (DRC), the per-case cost of treating Mpox is generally higher than that of conducting awareness campaigns, mainly due to the financial burden of hospitalization and specialized healthcare services. Although awareness interventions are typically more affordable on a per-unit basis, their cumulative cost can exceed that of treatment when deployed extensively over time, especially in regions with low incidence rates. These economic considerations are reflected in the construction of the objective functional in Equation (9), where weights are assigned to each compartment as follows: $w_1=20$, $w_2=3000$, $w_3=1000$, and $w_4=3000$. To reflect realistic cost differences, the cost coefficients for awareness and treatment are hypothetically set such that awareness (${\varphi }_1=1200$) is less costly than treatment (${\varphi }_2=5000$). Accordingly, Figure 7 presents the efficiency indices for different control strategies—awareness alone, treatment alone, and combined control—highlighting their respective effectiveness in minimizing infections across the countries studied, independent of cost considerations.

However, in Figure 7 control $u_1$ is the most effective in all countries due to its number of averted infections, especially in DRC. However, Control $u_2$ adds value mainly in Spain and Nigeria but is weak alone in DRC. Combining both controls ($u_1{u}_2$) gives the highest efficiency, showing strong synergy in most countries.

4.1. Incremental Cost Effectiveness Ratio

The Incremental Cost-Effectiveness Ratio (ICER) is a metric used to assess the relative cost and health benefits of two competing interventions, particularly when resources are limited. It is computed using the following expression:

$$ICER={\displaystyle \frac{\mathrm{change}\mathrm{in}\mathrm{total}\cos \mathrm{t}\mathrm{of}\mathrm{intervention}}{\mathrm{change}\mathrm{in}\mathrm{total}\mathrm{infections}\mathrm{averted}\mathrm{by}\mathrm{interventions}}}$$

(25)

When comparing two strategies (L and M), and M proves to be more effective than L (i.e., $TIA\left(M\right)>TIA\left(L\right)$), ICER becomes a key indicator for evaluating their relative value. Specifically, ICER quantifies the additional cost required to gain one extra unit of effectiveness when transitioning from strategy L to strategy M.

The ICER for the baseline strategy L is calculated as the ratio of its total cost ($TC\left(L\right)$) to its total effectiveness ($TIA\left(L\right)$):

$$ICER\left(L\right)={\displaystyle \frac{TC\left(L\right)}{TIA\left(L\right)}}$$

(26)

For the more effective strategy M, the ICER is determined by assessing the cost and effectiveness differences between it and the baseline strategy L. This is calculated as:

$$ICER\left(M\right)={\displaystyle \frac{TC\left(M\right)-TC\left(L\right)}{TIA\left(M\right)-TIA\left(L\right)}}$$

(27)

It is essential to clearly define the column headers for Table ( Table 4, Table 5 and Table 6). Column (A) represents the total infection in the absence of any control measures. Column (B) denotes the total infection associated with the implementation of a specific control strategy. The Total Infection Averted (TIA), given by the difference ($A-B$), quantifies thetotal number of infection averted through the application of the control strategy.

To determine the most cost-effective control strategy between the combined control ($u_1{u}_2$) and the awareness-only strategy ($u_1$), we apply the incremental cost-effectiveness ratio (ICER) formulas as defined in Equations (26) and (27), using the data presented in Table 4. Specifically, for Nigeria, the ICER values are computed as follows:

$$ICER\left(u_1{u}_2\right)={\displaystyle \frac{TC\left(u_1{u}_2\right)}{TA\left(u_1{u}_2\right)}}={\displaystyle \frac{3.04\times {10}^6}{2.06\times {10}^{10}}}\approx 0.00015$$

$$ICER\left(u_1\right)={\displaystyle \frac{TC\left(u_1\right)-TC\left(u_1{u}_2\right)}{TA\left(u_1\right)-TA\left(u_1{u}_2\right)}}={\displaystyle \frac{5.84\times {10}^5-3.04\times {10}^6}{1.73\times {10}^{10}-2.06\times {10}^{10}}}\approx 0.00075$$

Since $ICER\left(u_1{u}_2\right)<ICER\left(u_1\right)$, the awareness-only strategy ($u_1$) is considered less cost-effective and is thus eliminated from further consideration in Table 4. This indicates that the combined application of awareness and treatment controls yields a greater reduction in infection at a lower cost per case averted, making it the more economically efficient option.

Similarly, to assess the cost-effectiveness of the control strategies in Spain, we apply Equations (26) and (27) using the corresponding values from Table 4. The ICER values for the combined strategy ($u_1{u}_2$) and the awareness-only strategy ($u_1$) are calculated as follows:

$$ICER\left(u_1{u}_2\right)={\displaystyle \frac{TC\left(u_1{u}_2\right)}{TIA\left(u_1{u}_2\right)}}={\displaystyle \frac{3.07\times {10}^6}{5.18\times {10}^9}}\approx 0.00059$$

$$ICER\left(u_1\right)={\displaystyle \frac{TC\left(u_1\right)-TC\left(u_1{u}_2\right)}{TIA\left(u_1\right)-TIA\left(u_1{u}_2\right)}}={\displaystyle \frac{5.89\times {10}^5-3.07\times {10}^6}{4.02\times {10}^9-5.18\times {10}^9}}\approx 0.00213$$

Given that $ICER\left(u_1{u}_2\right)<ICER\left(u_1\right)$, the awareness-only control strategy ($u_1$) is considered less cost-effective in the Spanish context. Therefore, it is excluded from further evaluation in Table 4, reaffirming the economic advantage of the combined control strategy in achieving greater infection reduction per unit cost.

For Italy, the cost-effectiveness of the combined control strategy ($u_1{u}_2$) compared to the awareness-only strategy ($u_1$) is evaluated using the ICER formulas. Using the data in Table 4, the computations are as follows:

$$ICER\left(u_1{u}_2\right)={\displaystyle \frac{TC\left(u_1{u}_2\right)}{TIA\left(u_1{u}_2\right)}}={\displaystyle \frac{3.06\times {10}^6}{4.89\times {10}^9}}\approx 0.00063$$

$$ICER\left(u_1\right)={\displaystyle \frac{TC\left(u_1\right)-TC\left(u_1{u}_2\right)}{TIA\left(u_1\right)-TIA\left(u_1{u}_2\right)}}={\displaystyle \frac{5.88\times {10}^5-3.06\times {10}^6}{4.19\times {10}^9-4.89\times {10}^9}}\approx 0.00319$$

Since $ICER\left(u_1{u}_2\right)<ICER\left(u_1\right)$, the awareness-only control strategy ($u_1$) is deemed less cost-effective for Italy and is thus excluded from Table 4. The combined control ($u_1{u}_2$) offers a more efficient use of resources in minimizing infections.

For the Democratic Republic of Congo (DRC), we analyze the ICER values to determine the more cost-effective strategy between ($u_1{u}_2$) and ($u_1$), using the values from Table 4:

$$ICER\left(u_1{u}_2\right)={\displaystyle \frac{TC\left(u_1{u}_2\right)}{TIA\left(u_1{u}_2\right)}}={\displaystyle \frac{3.04\times {10}^6}{1.60\times {10}^{10}}}\approx 0.00019$$

$$ICER\left(u_1\right)={\displaystyle \frac{TC\left(u_1\right)-TC\left(u_1{u}_2\right)}{TIA\left(u_1\right)-TIA\left(u_1{u}_2\right)}}={\displaystyle \frac{5.87\times {10}^5-3.04\times {10}^6}{1.58\times {10}^{10}-1.60\times {10}^{10}}}\approx 0.00087$$

As $ICER\left(u_1{u}_2\right)<ICER\left(u_1\right)$, the awareness-only strategy ($u_1$) is again considered less cost-effective in the DRC and is excluded from further analysis. The combined strategy ($u_1{u}_2$) is identified as the more cost-efficient approach for controlling Mpox infections. Thus, based on earlier ICER analyses obtained for each country, the awareness-only strategy ($u_1$) was eliminated and the result presented in Table 5.

To further evaluate the cost-effectiveness of the implemented control strategies, we apply the incremental cost-effectiveness ratio (ICER) formulas defined in Equations (26) and (27). Using the data from Table 5, we compare the combined control strategy ($u_1{u}_2$) with the treatment only strategy ($u_2$). The resulting ICER values, alongside key efficiency and cost metrics, are presented in Table 6.

From Table 6, it is evident that the treatment only strategy ($u_2$) consistently exhibits a lower ICER value compared to the combined control ($u_1{u}_2$) across all countries. This indicates that, from a cost effectiveness perspective, implementing effective treatment and management control alone yields greater economic efficiency in preventing infections than the more resource intensive combined strategy ($u_1{u}_2$).

5. Conclusions

5.1. Discussion and Concluding Remarks

The analysis of cost-effectiveness for Mpox control strategies across Nigeria, Spain, Italy, and the DRC reveals critical insights into how resources can be optimally allocated for epidemic management in these regions. The results, drawn from Table 6, show that the effective treatment only strategy ($u_2$) is consistently more cost effective than the combined control strategy ($u_1{u}_2$) in all countries.

In Nigeria, where healthcare infrastructure may be stretched, the ICER for effective treatment only ($u_2$) is significantly lower at 0.00004 compared to the combined control strategy ($u_1{u}_2$), which stands at 0.00015. This indicates that prioritizing effective treatment rather than trying to combine it with awareness campaigns will yield more efficient use of limited resources. The cost effectiveness of effective treatment only becomes even more important when considering Nigeria’s resource constraints and the need for rapid, and impactful interventions that can be scaled effectively across the country.

Similarly, in Spain, with a relatively more robust healthcare system but still facing pressures from occasional outbreaks, the ICER for effective treatment only (0.00023) is lower than the combined control strategy (0.00059). This suggests that, even in a high income country, effective treatment interventions should take precedence as they offer the greatest cost benefit. Spain can benefit from maintaining a clear focus on treatment, which will allow for a more efficient allocation of resources without compromising the response to the epidemic.

In Italy, where the healthcare system is advanced but still vulnerable to the pressures of managing infectious disease outbreaks, the effective treatment only strategy ($u_2$) demonstrates a significant cost advantage, with an ICER of 0.00016, compared to the combined strategy ($u_1{u}_2$) with ICER of 0.00063. Given Italy’s extensive experience in the management of infectious diseases (as seen in its response to COVID-19), the findings suggest that even in well funded healthcare systems, effective treatment interventions are more cost effective and should remain the focal point of control efforts.

In DRC, a country with a fragile healthcare system and limited resources, the effective treatment only strategy is again more cost effective, with an ICER of 0.00004, much lower than the combined strategy at 0.00019. The findings for DRC highlight the importance of focusing on treatment as the primary response, particularly in low resource settings where cost efficiency is a key consideration. Given the challenges DRC faces with healthcare access and infrastructure, an effective treatment focused approach could potentially save more lives without overburdening the already limited healthcare capacity.

However, across these diverse settings ranging from low resource environments like Nigeria and DRC to more advanced systems in Spain and Italy the effective treatment only strategy proves to be the most cost effective option for controlling Mpox. This aligns with the realities each country faces in terms of healthcare infrastructure, financial constraints, and the need for effective measures to control epidemics.

5.2. Limitation of the Study and Future Work

One notable limitation of this study is the exclusion of zoonotic transmission dynamics in the modeling framework. Specifically, the model does not account for human interactions with animal reservoirs, such as those working in zoos, wildlife parks, veterinary clinics, pet ownership settings, and related environments. Incorporating a dedicated compartment to represent individuals who are regularly exposed to potential animal hosts would enhance the model’s realism and provide a more comprehensive understanding of Mpox transmission pathways. Future studies should consider integrating these dynamics to better inform control strategies, particularly in regions where human-animal interactions are frequent and significant.

Appendix A.1. Positivity and Boundedness

For model 4 to be epidemiologically meaningful, it is imperative to show that all its state variables are non-negative for all time (t), and that $\mathcal{D}$ is indeed bounded. Thus, we adopt the following theorem:

Theorem A1.

let the initial values for the model 4 be $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I_1\left(0\right)\ge 0,I_2\left(0\right)\ge 0,H\left(0\right)\ge 0,R\left(0\right)\ge 0,S_m\left(0\right)\ge 0,l_m\left(0\right)\ge 0$. Then the solutions $(S,E,I_1,I_{2}H,R,S_m{I}_m)$ are positive for all time $(t>0)$.

Proof. 

Let $t_1=supt>0:S>0,E>0,I_1>0,I_2>0,H>0,R>0,S_m>0,I_m>0\in [0,t].$ Thus, $t_1>0$ We have, from the first equation of model 4 that

$${\displaystyle \frac{dS}{dt}}=-({\lambda }_h+{\mu }_h)S,where{\lambda }_h={\displaystyle \frac{{\beta }_h(I_1+\sigma I_2)}{N_h}}$$

which can be written as:

$$\int {\displaystyle \frac{dS}{dt}}=-\int {\displaystyle \frac{dS}{dt}},$$

so that

$$S\left(t_1\right)=S\left(0\right)exp\left[-{\int }_0^{t_1}\lambda \left(u\right).du\right]>0,$$

Similarly, this can be shown that $E>0,I_1>0,I_2>0,H>0,R>0,S_m>0,I_m>0$. □

The SEIHR X SI formulated model is represented by differential equations in system 4, which is to be analyzed in a feasible region $\mathcal{D}$, and all state variables and parameters of the model are assumed to be positive $\forall t\ge 0$. The bounded region is obtained through the following lemma.

Lemma A1.

The compact $\mathcal{D}$ defined as (${\Delta }_1\times {\Delta }_2$) is a positively invariant set, which attracts all positive orbits in ${\mathcal{R}}_{+}^8$.

Proof. 

According to the system of Equation (4), we have: Proof. Since $N_h\left(t\right)=S\left(t\right)+E\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+R_h\left(t\right)$ and $N_m\left(t\right)=S_m\left(t\right)+I_m\left(t\right)$, then the derivative of $N_h\left(t\right)$ is given by

$${\displaystyle \frac{d{N}_h}{dt}}={\displaystyle \frac{dS}{dt}}+{\displaystyle \frac{dE}{dt}}+{\displaystyle \frac{d{I}_1}{dt}}+{\displaystyle \frac{d{1}_2}{dt}}+{\displaystyle \frac{dH}{dt}}+{\displaystyle \frac{dR}{dt}}$$

(A1)

and

$${\displaystyle \frac{d{N}_m}{dt}}={\displaystyle \frac{d{S}_m}{dt}}+{\displaystyle \frac{d{I}_m}{dt}}$$

(A2)

Simplifying eqn (A1) becomes:

$$\begin{array}{c}\hfill {\displaystyle \frac{d{N}_h}{dt}}={\Pi }_h-{\mu }_h{N}_h-({\delta }_1{I}_h+{\delta }_2{I}_2+{\delta }_{3}H\end{array}$$

(A3)

Since ${\delta }_1,{\delta }_2,{\delta }_3>0$ and $I_1\left(t\right)I_2\left(t\right),H\left(t\right),\ge 0$ for all $t\ge 0$, eqn (A3) becomes:

$$\begin{array}{c}\hfill {\displaystyle \frac{d{N}_h}{dt}}\le {\Pi }_h-{\mu }_h{N}_h\end{array}$$

(A4)

Integrating eqn (A4) and solving gives:

$$\begin{array}{cc}\hfill N_h\left(t\right)\le & {\displaystyle \frac{{\Pi }_h}{{\mu }_h}}+\left(N_h\left(0\right)-{\displaystyle \frac{{\Pi }_h}{{\mu }_h}}\right){\mathbf{e}}^{-{\mu }_{h}t}\hfill \\ \hfill \le & N_h\left(0\right){\mathbf{e}}^{-{\mu }_{h}t}+{\displaystyle \frac{{\Pi }_h}{{\mu }_h}}\left(1-{\mathbf{e}}^{-{\mu }_{h}t}\right)\hfill \end{array}$$

(A5)

Similarly, solving the second equation of (A2):

$$\begin{array}{cc}\hfill N_m\left(t\right)\le & {\displaystyle \frac{{\Pi }_m}{\varphi +\delta }}+\left(N_R\left(0\right)-{\displaystyle \frac{{\Pi }_m}{{\mu }_m}}\right){\mathbf{e}}^{-{\mu }_{m}t}\hfill \\ \hfill \le & N_m\left(0\right){\mathbf{e}}^{-{\mu }_{m}t}+{\displaystyle \frac{{\Pi }_m}{{\mu }_m}}\left(1-{\mathbf{e}}^{-{\mu }_{m}t}\right)\hfill \end{array}$$

(A6)

Therefore, as $t\ge 0,N_h\left(t\right)\le {\displaystyle \frac{{\Pi }_h}{{\mu }_h}}$ and $N_m\left(t\right)\le {\displaystyle \frac{{\Pi }_m}{{\mu }_m}}$. Hence:

$$\begin{array}{cc}\hfill {\Omega }_1=& \left\{(S,E,I_1{I}_{2}H,R_h)\in {\mathcal{R}}_{+}^6:N_h\le {\displaystyle \frac{{\Pi }_h}{{\mu }_h}}\right\}\hfill \\ \hfill {\Omega }_2=& \left\{(S_m,I_m)\in {\mathcal{R}}_{+}^2:N_m\le {\displaystyle \frac{{\Pi }_m}{{\mu }_m}}\right\}\hfill \\ \hfill \mathcal{D}=& \left\{{\Omega }_1\times {\Omega }_2|N_h\le {\displaystyle \frac{{\Pi }_h}{{\mu }_h}},N_m\le {\displaystyle \frac{{\Pi }_m}{{\mu }_m}}\right\}\hfill \end{array}$$

(A7)

□

This implies that ${\displaystyle \frac{{\Pi }_h}{{\mu }_h}}$ and ${\displaystyle \frac{{\Pi }_m}{{\mu }_m}}$ are the upper bond for human population $N_h\left(t\right)$ and mammal reservoirs, $N_m\left(t\right)$ respectively. Thus, from Lemma A1 the set $\Omega$ is positive invariant with respect to the system (4). Then we conclude that the set $\Omega$ is positively invariant and all solutions of system (4) are non-negative and epidemiologically well-posed.

Appendix A.2. Disease-Free Equilibrium

At the disease-free equilibrium (DFE), the compartment $E=I_1=I_2=H=R=l_m=0.$ The DFE of the model ${\mathcal{E}}_0=(S^0,E^0,I_1^0,I_2^0,H^0,R^0,S_m^0,l_m^0)$, is given as

$${\mathcal{E}}_0=({\displaystyle \frac{{\Pi }_h}{{\mu }_h}},0,0,0,0,0,{\displaystyle \frac{{\Pi }_m}{{\mu }_m}},0)$$

To compute the basic reproduction number ${\mathcal{R}}_0$, we used the next-generation matrix approach. The rate of transfer vectors into and out of the affected compartments f and v respectively, are given by

$$f=\left[\begin{array}{c}{\displaystyle \frac{S{\beta }_h(I_1+I_2\sigma )}{N_h}}\\ 0\\ 0\\ S_m({\displaystyle \frac{I_m{\beta }_m}{N_m}}+{\displaystyle \frac{\alpha {\beta }_{hm}(I_1+I_2\sigma )}{N_h}})\end{array}\right],v=\left[\begin{array}{c}E\left({\gamma }_h+{\mu }_h\right)\\ -E{\gamma }_h\left(1-\theta \right)+I_1\left({\delta }_1+{\mu }_h+{\tau }_1\right)\\ -E{\gamma }_h\theta +I_2\left({\delta }_2+{\mu }_h+\psi +{\tau }_2\right)\\ H\left({\delta }_3+\eta +{\mu }_h\right)-I_1{\tau }_1-I_2{\tau }_2\\ I_m{\mu }_m\end{array}\right]$$

The Jacobian matrix of f and v are given by;

$$F=\left[\begin{array}{ccccc}0& {\beta }_h& {\beta }_h\sigma & 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& {\displaystyle \frac{{\Pi }_m\alpha {\beta }_{hm}{\mu }_h}{{\Pi }_h{\mu }_m}}& {\displaystyle \frac{{\Pi }_m\alpha {\beta }_{hm}{\mu }_h\sigma }{{\Pi }_h{\mu }_m}}& 0& {\displaystyle \frac{{\Pi }_m{\beta }_m}{N_m{\mu }_m}}\end{array}\right],\mathrm{and}V=\left[\begin{array}{ccccc}{m}_0& 0& 0& 0& 0\\ -{\gamma }_h\left(1-\theta \right)& m_1& 0& 0& 0\\ -{\gamma }_h\theta & 0& m_3& 0& 0\\ 0& -{\tau }_1& -{\tau }_2& m_4& 0\\ 0& 0& 0& 0& {\mu }_m\end{array}\right]$$

So that,

$$\begin{array}{cc}\hfill \rho \left(F{V}^{-1}\right)={\mathcal{R}}_0& ={\displaystyle \frac{{\beta }_h{\gamma }_h\left[\sigma \theta \left({\delta }_1+{\mu }_h+{\tau }_1\right)+\left(1-\theta \right)\left({\delta }_2+{\mu }_h+\psi +{\tau }_2\right)\right]}{\left({\gamma }_h+{\mu }_h\right)\left({\delta }_1+{\mu }_h+{\tau }_1\right)\left({\delta }_2+{\mu }_h+\psi +{\tau }_2\right)}}\hfill \end{array}$$

(A8)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{{\beta }_h{\gamma }_h\left[\sigma \theta m_1+(1-\theta )m_2\right]}{m_0{m}_1{m}_2}}\hfill \end{array}$$

(A9)

where ${\mathcal{R}}_0$ is positive since all parameters are positive and $0\le \theta \le 1$.

Appendix A.3. Existence and uniqueness of solution

We show that eqn (4) has a unique solution $(S,E,I_1,I_2,H,R,S_m,l_m)\in {\mathbb{R}}_{+}^8$ under the condition eqn (5) using the Lipschitz condition.

Consider a general form of ODE $y^{\prime }=f(x,y)$ with $y\left(x_0\right)=y_0,x\in \mathbb{R}$. If f is continuous and has partial derivatives that are continuous in the region $\mathbb{R}$, then $f(x,y)$ is Lipschitz continuous.

We re-write our model 4 in the form $F\left(x\right)$ with solution $x\in {\mathbb{R}}_{+}^8$

$$\begin{array}{cc}\hfill F_1& ={\Pi }_h-({\lambda }_h+{\mu }_h)S+ϵR,\hfill \\ \hfill F_2& ={\lambda }_{h}S-\gamma E_h-{\mu }_h{E}_h\hfill \\ \hfill F_3& =(1-\theta )\gamma E_h-({\tau }_1+{\delta }_1+{\mu }_h)I_1,\hfill \\ \hfill F_4& =\theta \gamma E_h-({\tau }_2+{\delta }_2+{\mu }_h+\psi )I_2,\hfill \\ \hfill F_5& ={\tau }_1{I}_1+{\tau }_2{I}_2-({\delta }_3+{\mu }_h+\eta )H\hfill \\ \hfill F_6& =\eta H+\psi I_2-(ϵ+{\mu }_h)R,\hfill \\ \hfill F_7& ={\Pi }_m-{\lambda }_m{S}_m-{\mu }_m{S}_m,\hfill \\ \hfill F_8& ={\lambda }_m{S}_m-{\mu }_m{I}_m,\hfill \end{array}$$

(A10)

Lemma A2.

Let $F\left(x\right)$ be continuous differentiable at $x\in {\mathbb{R}}_{+}^8$.

Theorem A2.

The model (4) has a unique solution if $F\left(x\right)$ is continuous and have continuous partial derivative.

Proof. 

Some partial derivatives of $F\left(x\right)$ are given below

${\displaystyle \frac{\partial F_1}{\partial x_1}}=|-({\lambda }_h+{\mu }_h)|<\infty {\displaystyle \frac{\partial F_1}{\partial x_2}}=0<\infty ;{\displaystyle \frac{\partial F_1}{\partial x_3}}=|{\displaystyle \frac{{\beta }_h{x}_1}{N_h}}|<\infty {\displaystyle \frac{\partial F_1}{\partial x_4}}=\left|{\displaystyle \frac{{\beta }_h\sigma x_1}{N_h}}\right|<\infty$

${\displaystyle \frac{\partial F_1}{\partial x_5}}=0<\infty ;{\displaystyle \frac{\partial F_1}{\partial x_6}}=ϵ<\infty {\displaystyle \frac{\partial F_1}{\partial x_7}}=0<\infty ;{\displaystyle \frac{\partial F_1}{\partial x_8}}=0<\infty$

The rest of the partial derivatives exist, are continuous, and are bounded in the same way as eqn (A10). Hence, by Theorem A2, the model (4) has a unique solution. □

Appendix A.3.1. Endemic equilibrium

The model (4) has an endemic equilibrium ${\mathcal{E}}_1=(S^{**},E^{**},I_1^{**},I_2^{**},H^{**},R^{**},S_m^{**},l_m^{**})$ with ${\mathcal{E}}_1=0$ is given as: $S_m^{**}={\displaystyle \frac{{\Pi }_m}{{\lambda }_m^{*}+{\mu }_m}},I_m^{**}={\displaystyle \frac{{\Pi }_m}{{\mu }_m({\lambda }_m^{*}+{\mu }_m)}},E^{**}={\displaystyle \frac{{\lambda }_h^{**}{S}_h^{**}}{m_0}},I_1^{**}={\displaystyle \frac{(1-\theta )\gamma {\lambda }_h^{**}{S}_h^{**}}{m_0{m}_1}}{I}_2^{**}={\displaystyle \frac{\theta \gamma {\lambda }_h^{**}{S}_h^{**}}{m_0{m}_2}}$

$$H^{**}={\displaystyle \frac{{\tau }_1(1-\theta )\gamma m_2{\lambda }_h^{**}{S}_h^{**}+{\tau }_2\theta \gamma m_1{\lambda }_h^{**}{S}_h^{**}}{m_0{m}_1{m}_2{m}_3}},R^{**}={\displaystyle \frac{\eta {\tau }_1(1-\theta )\gamma m_2{\lambda }_h^{**}{S}_h^{**}+\eta {\tau }_2\theta \gamma m_1{\lambda }_h^{**}{S}_h^{**}+\psi \theta m_1{m}_2\gamma {\lambda }_h^{**}{S}^{**}}{m_0{m}_1{m}_2{m}_3{m}_4}}$$

Substituting Equation (A.3.1) into the forces of infection ${\lambda }_h^{**}$ and ${\lambda }_m^{**}$ and simplifying, such that

$${\lambda }_h^{**}={\displaystyle \frac{{\beta }_h(I_1^{**}+\sigma I_2^{**})}{N_h}},{\lambda }_m={\displaystyle \frac{\alpha {\beta }_{hm}(I_1^{**}+\sigma I_2^{**})}{N_h}}+{\displaystyle \frac{{\beta }_m{I}_m^{**}}{N_m}}$$

we have,

$$p_1{\left({\lambda }_h^{**}\right)}^2+p_2\left({\lambda }_h^{**}\right)=0$$

where $m_0=(\gamma +{\mu }_h)$, $m_1=({\tau }_1+{\delta }_1+{\mu }_h)$, $m_2=({\tau }_2+{\delta }_2++{\mu }_h+\psi )$, $m_3={\delta }_3+{\mu }_h+\eta$, $m_4=(ϵ+{\mu }_h)$, $p_1=m_1{m}_2{m}_3{m}_4{\Pi }_h+m_2{m}_3{m}_4(1-\theta )\gamma {\Pi }_h+m_1{m}_3{m}_4\theta \gamma {\Pi }_h+m_0{m}_4{\tau }_1(1-\theta )\gamma {\Pi }_h+m_1{m}_4{\tau }_2\theta \gamma {\Pi }_h+\eta m_0{\tau }_1(1-\theta )\gamma {\Pi }_h+\eta m_1{\tau }_2\theta \gamma {\Pi }_h$,$p_2=M_0{m}_1{m}_2{m}_3{m}_4{\Pi }_h-{\beta }_h{m}_2{m}_3{m}_4(1-\theta )\gamma {\Pi }_h-{\beta }_h{m}_1{m}_3{m}_4\sigma \theta \gamma {\Pi }_h$.

Appendix A.4. Local stability of Mpox-free equilibrium

We analyze the local stability of Mpox free equilibrium of the model system (4) by using the control reproduction number ${\mathcal{R}}_0$ in the following theorem as described in [26].

Theorem A3.

The Mpox free equilibrium ${\mathcal{E}}_0$, of the model (4) is locally asymptotically stable in the biological feasible region $\mathcal{D}$ if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof. 

In order to prove Theorem A3, the Jacobian matrix of the model (4) was obtained at Mpox disease-free equilibrium ${\mathcal{E}}_0$, where $E=I_1=I_2=H=R=I_m=0$, and $N_h=S$ and $N_m=S_m$. as

$${\mathrm{J}}_{\left({\mathcal{E}}_0\right)}=\left[\begin{array}{cccccccc}-{\mu }_h& 0& -{\beta }_h& -\sigma {\beta }_h& 0& ϵ& 0& 0\\ \\ 0& -(\gamma +{\mu }_h)& {\beta }_h& \sigma {\beta }_h& 0& 0& 0& 0\\ \\ 0& (1-\theta )\gamma & -({\tau }_1+{\delta }_1+{\mu }_h)& 0& 0& 0& 0& 0\\ \\ 0& \theta \gamma & 0& -({\tau }_2+{\delta }_2+{\mu }_h+\psi )& 0& 0& 0& 0\\ \\ 0& 0& {\tau }_1& {\tau }_2& -({\delta }_3+{\mu }_h+\eta )& 0& 0& 0\\ \\ 0& 0& 0& \psi & \eta & -(ϵ+{\mu }_h)& 0& 0\\ \\ 0& 0& -{\displaystyle \frac{\alpha {\beta }_{hm}{\Pi }_m{\mu }_h}{{\mu }_m{\Pi }_h}}& -{\displaystyle \frac{\sigma \alpha {\beta }_{hm}{\Pi }_m{\mu }_h}{{\mu }_m{\Pi }_h}}& 0& 0& -{\mu }_m& -{\beta }_m\\ \\ 0& 0& {\displaystyle \frac{\alpha {\beta }_{hm}{\Pi }_m{\mu }_h}{{\mu }_m{\Pi }_h}}& {\displaystyle \frac{\sigma \alpha {\beta }_{hm}{\Pi }_m{\mu }_h}{{\mu }_m{\Pi }_h}}& 0& 0& 0& -({\mu }_m-{\beta }_m)\end{array}\right]$$

(A11)

From eqn (A11), it is sufficient to show that all the eigenvalues of ${\mathcal{E}}_0$ are negative. We obtain the first four eigenvalues as: $-{\mu }_h,-{\mu }_m-(ϵ+{\mu }_h),-({\delta }_3+\eta +{\mu }_h),-({\mu }_m-{\beta }_m)$, while the remaining eigenvalues can be obtained from the sub-matrix ${\mathrm{J}}_1\left({\mathcal{E}}_0\right)$, which is given by:

$${{\mathrm{J}}_1}_{\left({\mathcal{E}}_0\right)}=\left[\begin{array}{ccc}-(\gamma +{\mu }_h)& {\beta }_h& \sigma {\beta }_h\\ (1-\theta )\gamma & -({\tau }_1+{\delta }_1+{\mu }_h)& 0\\ \theta \gamma & 0& -({\tau }_2+{\delta }_2+{\mu }_h+\psi )\end{array}\right]=\left[\begin{array}{ccc}-m_0& {\beta }_h& \sigma {\beta }_h\\ (1-\theta )\gamma & -m_1& 0\\ \theta \gamma & 0& -m_2\end{array}\right]$$

(A12)

The solutions of the characteristic polynomial for eqn (A12) are given as

$$x^3+{\varphi }_1{x}^2+{\varphi }_{2}x+{\varphi }_3=0$$

(A13)

where

${\varphi }_1=m_2+m_1+m_0$

${\varphi }_2=m_0(m_1+m_2)+m_1{m}_2-{\beta }_h\gamma [(1-\theta )+\theta \sigma ]$

${\varphi }_3=m_1{m}_0{m}_2(1-{\mathcal{R}}_0)$

Applying Routh-Hurwitz criterion, the cubic eqn (A13) will have a roots with negative real parts if and only if ${\varphi }_1>0,{\psi }_3>0and{\varphi }_1{\varphi }_2>{\varphi }_3.$ Clearly, ${\varphi }_1>0and{\varphi }_3>0,(if{\mathcal{R}}_0<\infty )$. As a result, the disease free equilibrium $\left({\mathcal{E}}_0\right)$ is locally asymptotically stable if ${\mathcal{R}}_0<\infty )$

□