Modeling Mpox Transmission and Burden Across VariousVaccination Levels in the DRC Using an Age-Structured Model

1. Introduction

Mpox is a viral disease caused by the mpox virus (MPXV), which belongs to the Orthopoxvirus genus. This virus mainly affects animals such as monkeys, squirrels, and rodents [1]. Additionally, Mpox is classified as a zoonotic disease, meaning it can be transmitted from animals to humans. However, transmission is not limited to animal-to-human interactions; human-to-human transmission has also been documented. Researchers suggest that Mpox spreads among humans through close contact with infected individuals, particularly through respiratory droplets or direct contact with skin lesions and other bodily fluids. In addition, environmental transmission is possible through contact with contaminated objects or surfaces. Common symptoms of Mpox infection include fever, headache, muscle pain, and fatigue [2].

Although MPOX cases are found in numerous countries worldwide, they are most commonly found in Africa. The first monkeypox case was identified in 1958 following the occurrence of outbreaks in a group of monkeys that were being studied. Globally, there have been a total of 88,122 reported cases, with 30,555 of them occurring in the United States [3].

Since 2017, there has been a resurgence of Mpox in Nigeria, the Democratic Republic of Congo (DRC), and other parts of Africa, attributed to the waning immunity to smallpox vaccines and the accumulation of cohorts that have never been vaccinated against smallpox [4]. In 2022, a global outbreak of Mpox resulted in more than 91,000 confirmed cases in 115 countries, establishing chains of human-to-human transmission and leading to a renewed focus on vaccination as a preventive measure [5]. In 2023, there were 14,626 suspected cases of clade I Mpox in the Democratic Republic of Congo, more than double the 6,216 cases recorded in 2020 [6]. The sexual transmission of clade I Mpox was observed for the first time, with outbreaks among MSM (men who have sex with men) networks and commercial sex workers [6]. However, unlike the clade IIb outbreak, where more than 95% infections were recorded among adult men [7], cases of the ongoing DRC epidemic have been concentrated in children. From 1 January to 26 May 2024, 50% of the confirmed cases were under 15 years old, reflecting the age patterns seen in the MPOX dynamics [6]. In addition, preliminary evidence from the DRC indicates that children under 15 years old accounted for 86% (254 out of 296) of reported deaths in 2024, suggesting a potentially higher severity of Mpox infections in younger age groups [6].

Age-specific variations in disease severity are crucial for understanding the burden of Mpox and for designing targeted interventions to reduce overall morbidity and mortality. Recent studies have provided age-specific severity estimates for clade I Mpox, highlighting a significantly higher risk of severe outcomes among younger populations, particularly children under 15 years old [8]. Furthermore, it remains unclear to what extent age-related severity differences can be explained by the concentration of vaccine coverage in older age groups [8]. Public health policies addressing the Mpox outbreak require adequate disease modeling strategies. This need has been emphasized during the COVID-19 pandemic, where mathematical models of disease transmission have played a critical role in shaping public health policies [9]. However, Mpox has historically been a neglected zoonotic disease, resulting in limited development of mathematical models for its transmission, particularly those assessing the impact of vaccination.

Various mathematical models have been developed to understand Mpox dynamics before, during, and after the 2022-2023 outbreak. Some of these models specifically focus on the transmission of Mpox from animals to humans [10,11,12,13,14,15,16], while others examine sexual transmission or spread within high-risk populations [17,18,19,20]. Additional models have explored the impact of smallpox or Mpox vaccination on curbing the virus’s spread [21,22], while others have investigated the co-interaction between Mpox and other diseases [17,23]. Despite these contributions, the impact of vaccination on MPOX transmission and its overall burden remains a critical and pressing concern.

In this study, we developed a mathematical age-structured model with two subgroups: individuals under 15 years old and those over 15. The model emphasizes the impact of vaccination, as the upcoming campaign primarily targets children under 15 [24]. We analyzed trends in MPOX cases and deaths averted in the DRC under various vaccination coverage levels.

2. Method

2.1. Model Formulation

In this study, we developed a deterministic compartmental model of Mpox stratified by age group, infection status and vaccination status to describe the impact of a vaccine on Mpox transmission. The mathematical model comprises twelve epidemiological states depending on the individual’s health, vaccination, and age status. The total population at time t denoted by $N\left(t\right)$ is divided into two subpopulations, meaning under-15 and over-15 years old individuals, which are represented by subscripts 1 and 2, respectively. Each subpopulation, with size $N_1$ or $N_2$ is further subdivided into six compartments that are Susceptible ($S_1$), Vaccinated ($V_1$), Exposed ($E_1$), Infectious ($I_1$), Detected infectious ($I_{d_1}$), and Recovered ($R_1$). Thus, the total population size for the under-15 years old individuals is given by:

$$N_1=S_1+V_1+E_1+I_1+I_{d1}+R_1.$$

Similarly, the over-15 year subpopulation denoted by $N_2$ is also further subdivided into six compartments: Susceptible ($S_2$), Vaccinated ($V_2$), Exposed ($E_1$), Infectious ($I_2$), Detected infectious ($I_{d_2}$), and Recovered ($R_2$). The total population size for the over-15 years old individuals is given by:

$$N_2=S_2+V_2+E_2+I_2+I_{d2}+R_2.$$

Therefore, the total population at time t, denoted by $N\left(t\right)$, is

$$N\left(t\right)=N_1\left(t\right)+N_2\left(t\right).$$

Newborns enter the system at a constant rate ($\Lambda$), joining the susceptible compartment, $S_1$. From $S_1$, individuals may be vaccinated at a rate ${\nu }_t$, transitioning to the vaccinated compartment $V_1$. As they age, individuals in $V_1$ move to $V_2$ at a constant ageing rate $\zeta$. Vaccinated individuals in $V_1$ and $V_2$ lose immunity over time at a rate $\omega$, returning to the susceptible compartments $S_1$ and $S_2$, respectively.

From the susceptible compartments ($S_1$ and $S_2$), individuals may become exposed ($E_1$ and $E_2$) through effective contact with confirmed infectious individuals ($I_{d1}$ and $I_{d2}$) or undetected infectious individuals ($I_1$ and $I_2$) at rates $(1-{\Phi }_1){\lambda }_1$ and $(1-{\Phi }_2){\lambda }_2$, respectively. Here, ${\lambda }_1$ and ${\lambda }_2$ represent the force of infection for under-15 and over-15 individuals (${\lambda }_1={\lambda }_2$), respectively, while ${\Phi }_1$ and ${\varphi }_2$ account for the reduction in the force of infection due to age-specific biological, behavioral, and epidemiological factors. Vaccinated individuals ($V_1$ and $V_2$) may also become exposed ($E_1$ and $E_2$) through effective contact with infected, at rates $(1-\eta )(1-{\Phi }_1){\lambda }_1$ and $(1-\eta )(1-{\Phi }_2){\lambda }_2$, where $\eta$ represents vaccine efficacy.

Once exposed, individuals progress to the infectious compartments ($I_1$ and $I_2$) at a rate $\alpha$. A fraction of infectious individuals may be detected through contact tracing and testing and move to the confirmed infectious compartments ($I_{d1}$ and $I_{d2}$) at rates ${\rho }_d$. Individuals in these compartments may recover at rates ${\gamma }_{d1}$ and ${\gamma }_{d2}$, transitioning to the recovered compartments ($R_1$ and $R_2$), or die at rates ${\delta }_{d1}$ and ${\delta }_{d2}$. Similarly, those in the undetected infectious compartments ($I_1$ and $I_2$) may recover at rates ${\gamma }_{u1}$ and ${\gamma }_{u2}$ or die at rates ${\delta }_{u1}$ and ${\delta }_{u2}$.

Natural mortality affects all compartments at rates ${\mu }_1$ and ${\mu }_2$. Furthermore, individuals in compartments associated with the under-15 subpopulation transition to the corresponding over-15 subpopulation compartment through ageing at the rate $\zeta$.

When developing the mathematical model, we made some assumptions or comments, which are as follows.

(i)

Vaccination is administered to under-15 years old individuals that are susceptible. The model doesn’t consider vaccination of detected and confirmed infectious individuals.

(ii)

We assume that there is homogeneous mixing among the population, which means that every individual in the community is equally likely to mix and acquire infections from each member when they make contact.

(iii)

Since the Mpox outbreak has persisted for a long time, we include the vital dynamics (birth and natural death) in the model.

(iv)

The proportion of under-15s in the total Mpox cases remains constant at 0.5 throughout the modeling period.

Using the flow diagram in Figure 1 and the above description, we obtain the following system of equations:

$$\begin{array}{cc}\hfill \dot{S_1}& =\Lambda -(1-{\Phi }_1){\lambda }_1{S}_1+\omega V_1-(\nu \left(t\right)+{\mu }_1+\zeta )S_1,\hfill \\ \hfill \dot{V_1}& =\nu \left(t\right)S_1-(1-\eta )(1-{\Phi }_1){\lambda }_1{V}_1-(\omega +{\mu }_1+\zeta )V_1,\hfill \\ \hfill {\dot{E}}_1& =(1-{\Phi }_1){\lambda }_1{S}_1+(1-\eta )(1-{\Phi }_1){\lambda }_1{V}_1-(\alpha +{\mu }_1+\zeta )E_1,\hfill \\ \hfill {\dot{I}}_1& =\alpha E_1-({\rho }_d+{\gamma }_{u1}+{\delta }_{u1}+{\mu }_1+\zeta )I_1,\hfill \\ \hfill {\dot{I}}_{d1}& ={\rho }_d{I}_1-({\gamma }_{d1}+{\delta }_{d1}+{\mu }_1+\zeta )I_{d1},\hfill \\ \hfill {\dot{R}}_1& ={\gamma }_{u1}{I}_1+{\gamma }_{d1}I{d}_1-({\mu }_1+\zeta )R_1,\hfill \\ \hfill {\dot{S}}_2& =\zeta S_1-(1-{\Phi }_2){\lambda }_2{S}_2+\omega V_2-{\mu }_2{S}_2,\hfill \\ \hfill {\dot{V}}_2& =\zeta V_1-(1-\eta )(1-{\Phi }_2){\lambda }_2{V}_2-(\omega +{\mu }_2)V_2,\hfill \\ \hfill {\dot{E}}_2& =\zeta E_1+(1-{\Phi }_2){\lambda }_2{S}_2+(1-\eta )(1-{\Phi }_2){\lambda }_2{V}_2-(\alpha +{\mu }_2)E_2,\hfill \\ \hfill {\dot{I}}_2& =\zeta I_1+\alpha E_2-({\rho }_d+{\gamma }_{u2}+{\delta }_{u2}+{\mu }_2)I_2,\hfill \\ \hfill {\dot{I}}_{d2}& =\zeta I{d}_1+{\rho }_d{I}_2-({\gamma }_{d2}+{\delta }_{d2}+{\mu }_2)I_{d2},\hfill \\ \hfill {\dot{R}}_2& =\zeta R_1+{\gamma }_{u2}{I}_2+{\gamma }_{d2}I{d}_2-{\mu }_2{R}_2,\hfill \\ \hfill \dot{N_v}& =\nu \left(t\right)S_1.\hfill \end{array}$$

(1)

In (1), ${\lambda }_1$ and ${\lambda }_2$ are the forces of infection for under-15 and over-15 individuals respectively, which are defined by:

$${\lambda }_1={\displaystyle \frac{{\beta }_{u1}{I}_1+{\beta }_{d1}{I}_{d1}+{\beta }_{u2}{I}_2+{\beta }_{d2}{I}_{d2}}{N}},{\lambda }_2={\lambda }_1,$$

where $N=N_1+N_2,N_1=S_1+V_1+E_1+I_1+I_{d1}+R_1,N_2=S_2+V_2+E_2+I_2+I_{d2}+R_2$ and ${\Phi }_2$ is a modification parameter that accounts for the fact that over 15 years are less likely to be infected.

The proportion $p_1$ of MPOX cases among those under 15 can be approximated using the expression:

$$p_1\approx {\displaystyle \frac{(1-{\Phi }_1)P_1}{(1-{\Phi }_1)N_1+(1-{\Phi }_2)P_2}},$$

(2)

where $P_1=0.235$ and $P_2=0.765$ are the proportions of under-15 and over-15 years old individuals, respectively, in DRC [25]. Thus, from (2), we deduce the expression of ${\Phi }_2$, which is:

$${\Phi }_2\approx 1-{\displaystyle \frac{(1-p_1)(1-{\Phi }_1)P_1}{p_1{P}_2}}.$$

(3)

With $p_1\approx 0.5$ as reported by [6],

$${\Phi }_2\approx 1-{\displaystyle \frac{(1-{\Phi }_1)P_1}{P_2}}.$$

(4)

From (4), if $0<{\Phi }_1\le 1$, one can notice that $0<{\Phi }_2\le 1$ since $P_1<P_2$.

The last equation of (1) is added to determine the cumulative number of vaccinated individuals by time t ($N_v\left(t\right)$), which is used to dynamically adjust the vaccination rate. This helps to ensure that the vaccination coverage over a given period reaches a target value (see section 2.2.2).

From (1), it follows that the dynamics of the whole population is governed by the equation: $\dot{N}={\dot{N}}_1+{\dot{N}}_2=\Lambda -{\mu }_1{N}_1-{\mu }_2{N}_2-{\delta }_{d1}{I}_{d1}-{\delta }_{u1}{I}_1-{\delta }_{d2}{I}_{d2}-{\delta }_{u2}{I}_2\le \Lambda -{\mu }_1{N}_1.$ Applying the comparison theorem, we obtain $N\left(t\right)\le {\displaystyle \frac{\Lambda }{{\mu }_1}}+(N_1\left(0\right)-{\displaystyle \frac{\Lambda }{{\mu }_1}})e^{-{\mu }_{1}t}$. Thus, ${lim}_{t\to \infty }N\left(t\right)\le {\displaystyle \frac{\Lambda }{{\mu }_1}}$. Let $y=(S_1,V_1,E_1,I_1,I_{d1},R_1,S_2,V_2,E_2,I_2,I_{d2},R_2)$ and $\Omega =\{y\in {\mathbb{R}}_{+}^{12}:y_{1\le i\le 12}\le {\displaystyle \frac{\Lambda }{{\mu }_1}}\}$. Thus, $\Omega$ is a forward invariant set. The model is analyzed in this set.

Table 1. State variables and their description.

State Variable	Description
$S_1\left(S_2\right)$	Under 15 (over 15) Susceptible individuals
$V_1\left(V_2\right)$	Under 15 (over 15) Vaccinated individuals
$E_1\left(E_2\right)$	Under 15 (over 15) Exposed individuals
$I_1\left(I_2\right)$	Under 15 (over 15) Infectious individuals
$I_{d_1}\left(I_{d_2}\right)$	Under 15 (over 15) Detected infectious individuals
$R_1\left(R_2\right)$	Under 15 (over 15) Recovered individuals

Table 1. State variables and their description.

State Variable	Description
$S_1\left(S_2\right)$	Under 15 (over 15) Susceptible individuals
$V_1\left(V_2\right)$	Under 15 (over 15) Vaccinated individuals
$E_1\left(E_2\right)$	Under 15 (over 15) Exposed individuals
$I_1\left(I_2\right)$	Under 15 (over 15) Infectious individuals
$I_{d_1}\left(I_{d_2}\right)$	Under 15 (over 15) Detected infectious individuals
$R_1\left(R_2\right)$	Under 15 (over 15) Recovered individuals

Table 2. Description of the model parameters.

Parameter	Description	Unit
${\Lambda }_0$	Birth rate in DRC	${\mathrm{Individuals}}^{-1}$${\mathrm{week}}^{-1}$
$\zeta$	Age-based transition rate	${\mathrm{Individuals}}^{-1}$${\mathrm{week}}^{-1}$
${\mu }_1\left({\mu }_2\right)$	Natural death rate of under 15 (over 15) individuals	${\mathrm{week}}^{-1}$
$\nu$	Vaccination Rate	${\mathrm{week}}^{-1}$
$\omega$	Waning rate of the vaccine efficacy	${\mathrm{week}}^{-1}$
$\eta$	Vaccine efficacy	dimensionless
${\Phi }_1$(${\Phi }_2$)	Percent reduction in the force of infection due to intervention measures (excluding vaccination) and acquired immunity.	dimensionless
${\beta }_{u_1}\left({\beta }_{u_2}\right)$	Transmission rates of under 15 (over 15) infectious individuals	${\mathrm{week}}^{-1}$
${\beta }_{d_1}\left({\beta }_{d_2}\right)$	Transmission rates of under 15 (over 15) detected infectious individuals	${\mathrm{week}}^{-1}$
$1/\alpha$	Latent period	${\mathrm{week}}^{-1}$
${\rho }_d$	Detection rate of infectious individuals	${\mathrm{week}}^{-1}$
${\delta }_{u_1}\left({\delta }_{u_2}\right)$	Mpox-induced death rate for under 15 (over 15) infectious individuals	${\mathrm{week}}^{-1}$
${\delta }_{d_1}\left({\delta }_{d_2}\right)$	Mpox-induced death rate for under 15 (over 15) detected infectious individuals	${\mathrm{week}}^{-1}$
${\gamma }_{u_1}\left({\gamma }_{u_2}\right)$	Recovery rate of under 15 (over 15) infectious individuals	${\mathrm{week}}^{-1}$
${\gamma }_{d_1}\left({\gamma }_{d_2}\right)$	Recovery rate of under 15 (over 15) detected infectious individuals	${\mathrm{week}}^{-1}$

Table 2. Description of the model parameters.

Parameter	Description	Unit
${\Lambda }_0$	Birth rate in DRC	${\mathrm{Individuals}}^{-1}$${\mathrm{week}}^{-1}$
$\zeta$	Age-based transition rate	${\mathrm{Individuals}}^{-1}$${\mathrm{week}}^{-1}$
${\mu }_1\left({\mu }_2\right)$	Natural death rate of under 15 (over 15) individuals	${\mathrm{week}}^{-1}$
$\nu$	Vaccination Rate	${\mathrm{week}}^{-1}$
$\omega$	Waning rate of the vaccine efficacy	${\mathrm{week}}^{-1}$
$\eta$	Vaccine efficacy	dimensionless
${\Phi }_1$(${\Phi }_2$)	Percent reduction in the force of infection due to intervention measures (excluding vaccination) and acquired immunity.	dimensionless
${\beta }_{u_1}\left({\beta }_{u_2}\right)$	Transmission rates of under 15 (over 15) infectious individuals	${\mathrm{week}}^{-1}$
${\beta }_{d_1}\left({\beta }_{d_2}\right)$	Transmission rates of under 15 (over 15) detected infectious individuals	${\mathrm{week}}^{-1}$
$1/\alpha$	Latent period	${\mathrm{week}}^{-1}$
${\rho }_d$	Detection rate of infectious individuals	${\mathrm{week}}^{-1}$
${\delta }_{u_1}\left({\delta }_{u_2}\right)$	Mpox-induced death rate for under 15 (over 15) infectious individuals	${\mathrm{week}}^{-1}$
${\delta }_{d_1}\left({\delta }_{d_2}\right)$	Mpox-induced death rate for under 15 (over 15) detected infectious individuals	${\mathrm{week}}^{-1}$
${\gamma }_{u_1}\left({\gamma }_{u_2}\right)$	Recovery rate of under 15 (over 15) infectious individuals	${\mathrm{week}}^{-1}$
${\gamma }_{d_1}\left({\gamma }_{d_2}\right)$	Recovery rate of under 15 (over 15) detected infectious individuals	${\mathrm{week}}^{-1}$

Table 3. Fixed parameter values of the model.

Parameter	Value	Reference
${\Lambda }_0$	$7.94\times {10}^{-4}$${\mathrm{Individuals}}^{-1}$. ${\mathrm{week}}^{-1}$	[26]
N	$109\times {10}^6$Individuals	[26]
$\omega$	$(1/32)$${\mathrm{week}}^{-1}$	[27]
$\eta$	0.8	[28]
$\zeta$	$(1/(15\times 52\left)\right)$${\mathrm{week}}^{-1}$	Calculated
${\mu }_1\left({\mu }_2\right)$	$3.5\times {10}^{-4}$ ($2.80\times {10}^{-4}$) ${\mathrm{week}}^{-1}$	Calculated
$\alpha$	$0.7$${\mathrm{week}}^{-1}$	[29]
${\delta }_{u_1}\left({\delta }_{u_2}\right)$	$0.0055$ (0.004) ${\mathrm{week}}^{-1}$	Assumed
${\delta }_{d_1}\left({\delta }_{d_2}\right)$	0.004 (0.003) ${\mathrm{week}}^{-1}$	${\mathrm{Calibrated}}^{*}$
${\gamma }_{u_1}\left({\gamma }_{u_2}\right)$	1/3.5 (1/3) ${\mathrm{week}}^{-1}$	Assumed
${\gamma }_{d_1}\left({\gamma }_{d_2}\right)$	1/3 (1/2.5) ${\mathrm{week}}^{-1}$	[29]

$*$${\delta }_{d_1}$ and ${\delta }_{d_2}$ are calibrated to ensure the model accurately aligns with the observed MPOX-related death incidence (our world in data).

Table 3. Fixed parameter values of the model.

Parameter	Value	Reference
${\Lambda }_0$	$7.94\times {10}^{-4}$${\mathrm{Individuals}}^{-1}$. ${\mathrm{week}}^{-1}$	[26]
N	$109\times {10}^6$Individuals	[26]
$\omega$	$(1/32)$${\mathrm{week}}^{-1}$	[27]
$\eta$	0.8	[28]
$\zeta$	$(1/(15\times 52\left)\right)$${\mathrm{week}}^{-1}$	Calculated
${\mu }_1\left({\mu }_2\right)$	$3.5\times {10}^{-4}$ ($2.80\times {10}^{-4}$) ${\mathrm{week}}^{-1}$	Calculated
$\alpha$	$0.7$${\mathrm{week}}^{-1}$	[29]
${\delta }_{u_1}\left({\delta }_{u_2}\right)$	$0.0055$ (0.004) ${\mathrm{week}}^{-1}$	Assumed
${\delta }_{d_1}\left({\delta }_{d_2}\right)$	0.004 (0.003) ${\mathrm{week}}^{-1}$	${\mathrm{Calibrated}}^{*}$
${\gamma }_{u_1}\left({\gamma }_{u_2}\right)$	1/3.5 (1/3) ${\mathrm{week}}^{-1}$	Assumed
${\gamma }_{d_1}\left({\gamma }_{d_2}\right)$	1/3 (1/2.5) ${\mathrm{week}}^{-1}$	[29]

$*$${\delta }_{d_1}$ and ${\delta }_{d_2}$ are calibrated to ensure the model accurately aligns with the observed MPOX-related death incidence (our world in data).

2.2. DRC MPOX Data and Parameter Estimation Procedure

2.2.1. DRC Mpox Data

Data on the daily MPOX cases in DRC were obtained from MPOX data respiratory by Our World in Data (https://ourworldindata.org/mpox, accessed on 15 November 2024). These data were converted in weekly cases and span from April 30, 2023, to September 22, 2024. The start date considered for the modeling aligns with the initial establishment of a clearer upward trend, signifying the onset of sustained transmission dynamics suitable for modeling. It also ensures sufficient data points for early transmission patterns while avoiding noise from sporadic cases in the preceding months. In the DRC, where MPOX is endemic or elsewhere, the reported confirmed cases are infectious individuals who have been tested (especially with PCR). Hence, we assume that the weekly number of confirmed MPOX cases corresponds to newly detected infectious individuals from the $ODEs$ (${\rho }_d{I}_1+{\rho }_d{I}_2+\zeta I_{d_1}$). Demographic data (population size, birth rate and life expectancy) for the country were obtained from [26].

2.2.2. Model Fiting and Parameter Estimation Procedure

We consider a mixed population where under-15 and over-15-year-old individuals interact, with the transition from the under-15 subpopulation to the over-15 subpopulation governed by the parameter $\zeta$. We set $V_1\left(0\right)=E_1\left(0\right)=R_1\left(0\right)=V_2\left(0\right)=E_2\left(0\right)=R_2\left(0\right)=0$; $I_1\left(0\right)=I_2\left(0\right)=(I_{d1}\left(0\right)+I_{d2}\left(0\right))/\left(2{\rho }_d\right)$; $S_1\left(0\right)=N_1-I_1\left(0\right)-I_{d1}\left(0\right)$ and $S_2\left(0\right)=N_2-I_2\left(0\right)-I_{d2}\left(0\right)$. The expression $I_{d1}\left(0\right)+I_{d2}\left(0\right)=I_d\left(0\right)$ represents the number of reported MPOX cases in DRC on April 30, 2024. For simplicity, we assume $I_{d2}\left(0\right)=I_{d2}\left(0\right)=I_d\left(0\right)/2$.

The solutions of (1) were obtained using the MatlabR2024b built-in function $ODE45$. The minimization was achieved through the built-in constraint minimization algorithm fminsearchbnd. The fitting was carried out using a nonlinear least squares technique, finding the best set of parameter values that minimizes the root mean square differences between the observed weekly confirmed MPOX cases and the predicted weekly cases from the model (1). Thus, we estimated the best values of the unknown parameters (${\beta }_{u_1}$, ${\beta }_{u_2}$, ${\beta }_{d_1}$, ${\beta }_{d_2}$, ${\rho }_d$, ${\Phi }_1$, and ${\Phi }_2$) in the model (1) using a cross-validation procedure of parameter estimation [30], where 90% of the observations were considered to estimate values of the seven unknown model parameters with an $RMS{E}_1$ value and the remaining observations were used to validate the model with an $RMS{E}_2$ value. We simulated 2,000 different starting values of the seven parameters using the resample function in Matlab. The final values considered for the seven unknown parameters were those related to both the smallest values of $RMS{E}_2$ and $RMS{E}_1$ to guarantee a good fit of the model and a good predictive power [30]. The fixed parameters used in the model-fitting process were obtained from the literature as presented in Table 3, while other fixed parameters were calculated, calibrated or assumed).

Thus, we estimated the best values of the unknown parameters (${\beta }_{u_1}$, ${\beta }_{u_2}$, ${\beta }_{d_1}$, ${\beta }_{d_2}$, ${\rho }_d$, ${\Phi }_1$, and ${\Phi }_2$) in the model (1) using a cross-validation procedure of parameter estimation [30], where 90% of the observations were considered to estimate values of the seven unknown model parameters with an $RMS{E}_1$ value and the remaining observations were used to validate the model with an $RMS{E}_2$ value. We simulated 2,000 different starting values of the seven parameters using the resample function in Matlab. The final values considered for the seven unknown parameters were those related to both the smallest values of $RMS{E}_2$ and $RMS{E}_1$ to guarantee a good fit of the model and a good predictive power [30]. The fixed parameters used in the model-fitting process were obtained from the literature and presented in Table 3, while other fixed parameters were calculated, calibrated or assumed.

In the fitting procedure, we introduce vaccination in the model from January to December, 2025, which targets under-15 individuals with LC16m8 vaccine following [31]. Before January 2025, all vaccination parameters ($\omega$, $\nu$, and $\eta$) were set to zero. Our model did not account for the vaccination campaign, which started in October 2024, since it prioritized health workers and other high-risk groups (adults) due to limited vaccine availability (265,000 doses) [24]. We define a vaccinated individual as one who has received at least one dose of the LC16m8 vaccine. We consider different vaccination coverage levels, $C_{\mathrm{target}}$ ($0\%$, $10\%$, $20\%$, $30\%$, and $40\%$). We implemented a dynamic vaccination rate, $\nu \left(t\right)$, designed to ensure that the vaccination coverage reaches the predefined target, $C_{\mathrm{target}}$, over the specified vaccination period, $T_{vacc}$ ($T_{vacc}=50$ weeks: January 12 - December 31, 2025). This approach explicitly accounts for the interplay between vaccination and other dynamic factors affecting the susceptible population.

The total number, $N_v$, of individuals vaccinated by the end of the modeling period, T (April 30, 2023 - December 31, 2025; $T=140$ weeks), is proportional to the initial susceptible population, $S_1\left(0\right)$, and the target coverage, $C_{\mathrm{target}}$:

$${\int }_0^T\nu \left(t\right)S_1\left(t\right)dt=C_{\mathrm{target}}·S_1\left(0\right),$$

(5)

where $\nu \left(t\right)$ is the time-dependent vaccination flow rate, $S_1\left(t\right)$ is the number of susceptibles at time t and $S_1\left(0\right)$ is the initial number of susceptibles.

To achieve $C_{\mathrm{target}}$ over the modeling period, the vaccination rate $\nu \left(t\right)$ is dynamically adjusted. At any time t during the modeling period ($0\le t\le T$), the vaccination rate is computed as:

$$\nu \left(t\right)=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}t<t_{\mathrm{start}},\hfill \\ {\displaystyle \frac{N_{vr}\left(t\right)}{(T-t+1)·S_1\left(t\right)}},\hfill & \mathrm{if}{t}_{\mathrm{start}}\le t\le T\hfill \end{array}\right.$$

(6)

where

• $t_{\mathrm{start}}$ is start time of the vaccination campaign. Since our modeling period started on April 30, 2023, and we consider the vaccination starting in the second week of January 2025, $t_{\mathrm{start}}$ = 89 weeks).
• $N_{vr}\left(t\right)$ is the remaining number of individuals to be vaccinated at time t, calculated as:

  $$N_{vr}\left(t\right)=N_v\left(T\right)-N_v\left(t\right)=C_{\mathrm{target}}·S_1\left(0\right)-N_v\left(t\right);$$

  (7)

  $N_v\left(t\right)$ is the total number of under-15 individuals vaccinated at time t, determined by the ODEs (see equations (1)).
• $(T-t+1)$ spreads the vaccination over the remaining weeks in the modeling period.

To quantify the vaccine impact, we determined the vaccine effectiveness in terms of effective reproduction number and weekly and total cases and deaths averted for one year (January to December 2025) considering the above vaccination levels ($10\%$, $20\%$, $30\%$, and $40\%$) combined with non-pharmaceutical interventions targeting under-15s (Figure 3, Figure 4, Figure 5 and Figure 6).

3. Results

3.1. Analytical Results

3.1.1. Basic and Control Reproduction Numbers

Before determining the control reproduction number, the disease-free equilibrium ($DFE$) of the model is given by:

$N_0=(S_1^0,V_1^0,E_1^0,I_1^0,I_{d_1}^0,R_1^0,S_2^0,V_2^0,E_2^0,I_2^0,I_{d_2}^0,R_2^0)=(S_1^0,V_1^0,0,0,0,0,S_2^0,V_2^0,0,0,0)$, where

$S_1^0={\displaystyle \frac{\Lambda ({\mu }_1+\omega +\zeta )}{({\mu }_1+\zeta )({\mu }_1+\omega +\nu +\zeta )}}$, $V_1^0={\displaystyle \frac{\Lambda \nu }{({\mu }_1+\zeta )({\mu }_1+\omega +\nu +\zeta )}}$,

$S_2^0={\displaystyle \frac{\Lambda \zeta \left[({\mu }_2+\omega )({\mu }_1+\omega +\zeta )+\nu \omega \right]}{{\mu }_2({\mu }_2+\omega )({\mu }_1+\zeta )({\mu }_1+\omega +\nu +\zeta )}}$, $V_2^0={\displaystyle \frac{\Lambda \nu \zeta }{({\mu }_2+\omega )({\mu }_1+\zeta )({\mu }_1+\omega +\nu +\zeta )}}$.

In the absence of vaccination, $\omega =\nu =0$ and $S_1^0={\displaystyle \frac{\Lambda }{({\mu }_1+\zeta )}}$, $V_1^0=0$, $S_2^0={\displaystyle \frac{\Lambda \zeta }{{\mu }_2({\mu }_1+\zeta )}}$, $V_2^0=0$.

To determine the control reproduction number, we use the next-generation approach. The new infection matrix F and the transition matrix V, all evaluated at DFE, are given by:

$F=\left(\begin{array}{cccccc}0& {\beta }_{u_1}{M}_1& {\beta }_{d_1}{M}_1& 0& {\beta }_{u_2}{M}_1& {\beta }_{d_2}{M}_1\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& {\beta }_{u_1}{M}_2& {\beta }_{d_1}{M}_2& 0& {\beta }_{u_2}{M}_2& {\beta }_{d_2}{M}_2\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$,

$V=\left(\begin{array}{cccccc}\alpha +{\mu }_1+\zeta & 0& 0& 0& 0& 0\\ -\alpha & A_1& 0& 0& 0& 0\\ 0& -{\rho }_d& A_2& 0& 0& 0\\ -\zeta & 0& 0& \alpha +{\mu }_2& 0& 0\\ 0& -\zeta & 0& -\alpha & A_3& 0\\ 0& 0& -\zeta & 0& -{\rho }_d& A_4\end{array}\right)$,

where,

$A_1={\mu }_1+{\gamma }_{u_1}+{\delta }_{u_1}+{\rho }_d+\zeta$, $A_2={\mu }_1+{\gamma }_{d_1}+{\delta }_{d_1}+\zeta$, $A_3={\mu }_2+{\gamma }_{u_2}+{\delta }_{u_2}+{\rho }_d$, $A_4={\mu }_2+{\gamma }_{d_2}+{\delta }_{d_2}$,

$M_{01}={\displaystyle \frac{(1-{\Phi }_1)(S_1^0+(1-\eta )V_1^0)}{N^0}}$, $M_{02}={\displaystyle \frac{(1-{\Phi }_2)(S_2^0+(1-\eta )V_2^0)}{N^0}}$.

$$V^{-1}=\left(\begin{array}{cccccc}{\displaystyle \frac{1}{\alpha +{\mu }_1+\zeta }}& 0& 0& 0& 0& 0\\ {\displaystyle \frac{\alpha }{A_1(\alpha +{\mu }_1+\zeta )}}& {\displaystyle \frac{1}{A_1}}& 0& 0& 0& 0\\ {\displaystyle \frac{{\rho }_d\alpha }{A_1{A}_2(\alpha +{\mu }_1+\zeta )}}& {\displaystyle \frac{{\rho }_d}{A_1{A}_2}}& {\displaystyle \frac{1}{A_2}}& 0& 0& 0\\ {\displaystyle \frac{\zeta }{(\alpha +{\mu }_1+\zeta )(\alpha +{\mu }_2)}}& 0& 0& {\displaystyle \frac{1}{\alpha +{\mu }_2}}& 0& 0\\ X_2& {\displaystyle \frac{\zeta }{A_1{A}_3}}& 0& {\displaystyle \frac{\alpha }{A_3(\alpha +{\mu }_2)}}& {\displaystyle \frac{1}{A_3}}& 0\\ X_1& {\displaystyle \frac{\zeta {\rho }_d\left(A_2+A_3\right)}{A_2{A}_1{A}_3{A}_4}}& {\displaystyle \frac{\zeta }{A_2{A}_4}}& {\displaystyle \frac{{\rho }_d\alpha }{A_3{A}_4(\alpha +{\mu }_2)}}& {\displaystyle \frac{{\rho }_d}{A_3{A}_4}}& {\displaystyle \frac{1}{A_4}}\end{array}\right)$$

$X_1={\displaystyle \frac{\zeta {\rho }_d\alpha \left(A_2{A}_1+A_2(\alpha +{\mu }_2)+A_3(\alpha +{\mu }_2)\right)}{A_2{A}_1{A}_3{A}_4(\alpha +{\mu }_2)(\alpha +{\mu }_1+\zeta )}}$, $X_2={\displaystyle \frac{\zeta \alpha \left(A_1+\alpha +{\mu }_2\right)}{A_1{A}_3(\alpha +{\mu }_1+\zeta )(\alpha +{\mu }_2)}}$.

The largest eigenvalue of $F\times V^{-1}$ is the control reproduction number of the model given by:

$$R_{cv}=R_{c{v}_1}+R_{c{v}_2},$$

(8)

where,

$R_{c{v}_1}={\displaystyle \frac{\alpha \left(A_5{\rho }_d\zeta {\beta }_{d_2}+A_6{\rho }_d{\beta }_{d_1}+A_6{A}_2{\beta }_{u_1}+A_7{A}_2\zeta {\beta }_{u_2}\right)M_{01}}{A_1{A}_2{A}_3{A}_4(\alpha +{\mu }_2)(\alpha +{\mu }_1+\zeta )}}$ and $R_{c{v}_2}={\displaystyle \frac{\alpha \left({\rho }_d{\beta }_{d_2}+A_4{\beta }_{u_2}\right)M_{02}}{A_3{A}_4(\alpha +{\mu }_2)}}$,

with $A_5=(\alpha +{\mu }_2)(A_2+A_3)+A_1{A}_2$, $A_6=A_4(\alpha +{\mu }_2)A_3$, $A_7=A_4(A_1+{\mu }_2+\alpha )$.

To evaluate the numerical value of the control reproduction in the absence of vaccination ($R_c^0$), we set $\omega =\nu =\eta =0$ in the expression of $R_c$. The control reproduction number in the absence of vaccination is given by:

$$R_c=R_{c_1}+R_{c_2},$$

(9)

with

$R_{c_1}={\displaystyle \frac{(1-{\Phi }_1)\alpha \left(A_5{\rho }_d\zeta {\beta }_{d_2}+A_6{\rho }_d{\beta }_{d_1}+A_6{A}_2{\beta }_{u_1}+A_7{A}_2\zeta {\beta }_{u_2}\right){\mu }_2}{A_1{A}_2{A}_3{A}_4(\alpha +{\mu }_2)(\alpha +{\mu }_1+\zeta )({\mu }_2+\zeta )}}$,

and

$R_{c_2}={\displaystyle \frac{(1-{\Phi }_2)\alpha \left({\rho }_d{\beta }_{d_2}+A_4{\beta }_{u_2}\right)\zeta }{A_3{A}_4(\alpha +{\mu }_2)({\mu }_2+\zeta )}}$.

It is worth mentioning that in the absence of vaccination and other control measures, $\nu =\omega =\eta ={\Phi }_1={\rho }_d={\Phi }_2=0$. Hence, $V_1^0=V_2^0=0$. Replacing these particular values into the expression of the control reproduction number yields the basic reproduction number, $R_0$, of the model given by:

$$R_0=R_{01}+R_{02},$$

(10)

where

$R_{01}={\displaystyle \frac{\alpha {\mu }_2}{{\mu }_2+\zeta }}\left({\displaystyle \frac{{\beta }_{u_1}}{y_1(\alpha +{\mu }_1+\zeta )}}\right)$ and $R_{02}={\displaystyle \frac{\alpha \zeta {\beta }_{u_2}}{(\alpha +{\mu }_2)({\mu }_2+\zeta )y_2}}$,

with $y_1={\mu }_1+{\gamma }_{u_1}+{\delta }_{u_1}+\zeta$, $y_2={\mu }_2+{\gamma }_{u_2}+{\delta }_{u_2}$ and $y_3=\alpha +{\mu }_1$.

It is important to highlight that $R_{02}$ is the contribution of the over 15-year-old individuals involving the age-transition rate, $\zeta$.

3.1.2. Effective Reproduction

The effective reproduction number, $R_{ev}$ is given by:

$$R_{ev}\left(t\right)=R_{e{v}_1}\left(t\right)+R_{e{v}_2}\left(t\right),$$

(11)

where

$R_{e{v}_1}\left(t\right)={\displaystyle \frac{\alpha \left(A_5{\rho }_d\zeta {\beta }_{d_2}+A_6{\rho }_d{\beta }_{d_1}+A_6{A}_2{\beta }_{u_1}+A_7{A}_2\zeta {\beta }_{u_2}\right)M_1\left(t\right)}{A_1{A}_2{A}_3{A}_4(\alpha +{\mu }_2)(\alpha +{\mu }_1+\zeta )}}$,

$R_{e{v}_2}\left(t\right)={\displaystyle \frac{\alpha \left({\rho }_d{\beta }_{d_2}+A_4{\beta }_{u_2}\right)M_2\left(t\right)}{A_3{A}_4(\alpha +{\mu }_2)}}$, with $M_1\left(t\right)={\displaystyle \frac{(1-{\Phi }_1)(S_1\left(t\right)+(1-\eta )V_1\left(t\right))}{N\left(t\right)}}$,

$M_2\left(t\right)={\displaystyle \frac{(1-{\Phi }_2)(S_2\left(t\right)+(1-\eta )V_2\left(t\right)))}{N\left(t\right)}}$.

The effective reproduction number in the absence of vaccination, $\widehat{R_e}$ is given by:

$$R_e\left(t\right)=R_{e_1}\left(t\right)+R_{e_2}\left(t\right),$$

(12)

where

$R_{e_1}\left(t\right)={\displaystyle \frac{(1-{\Phi }_1)\alpha \left(A_5{\rho }_d\zeta {\beta }_{d_2}+A_6{\rho }_d{\beta }_{d_1}+A_6{A}_2{\beta }_{u_1}+A_7{A}_2\zeta {\beta }_{u_2}\right)S_1\left(t\right)}{A_1{A}_2{A}_3{A}_4(\alpha +{\mu }_2)(\alpha +{\mu }_1+\zeta )N\left(t\right)}}$, and

$R_{e_2}\left(t\right)={\displaystyle \frac{(1-{\Phi }_2)\alpha \left({\rho }_d{\beta }_{d_2}+A_4{\beta }_{u_2}\right)S_2\left(t\right)}{A_3{A}_4(\alpha +{\mu }_2)N\left(t\right)}}$; the $A_i$ are as above.

3.2. Stability Analysis of the DFE and Possible Extinction of Mpox

Theorem 1.

The Disease Free Equilibrium DFE is Globally asymptotically stable whenever ${\mathcal{R}}_0<1$.

Proof. 

We prove this theorem by following the same approach in [32]. We first set $N={\displaystyle \frac{\Lambda }{{\mu }_1}}$ into (1) to get the limiting system (13):

$$\begin{array}{cc}\hfill \dot{S_1}& =\Lambda -(1-{\Phi }_1){\widehat{\lambda }}_1{S}_1+\omega V_1-(\nu +{\mu }_1+\zeta )S_1,\hfill \\ \hfill \dot{V_1}& =\nu S_1-(1-\eta )(1-{\Phi }_1){\widehat{\lambda }}_1{V}_1-(\omega +{\mu }_1+\zeta )V_1,\hfill \\ \hfill {\dot{E}}_1& =(1-{\Phi }_1){\widehat{\lambda }}_1{S}_1+(1-\eta )(1-\varphi ){\widehat{\lambda }}_1{V}_1-(\alpha +{\mu }_1+\zeta )E_1,\hfill \\ \hfill {\dot{I}}_1& =\alpha E_1-A_1{I}_1,\hfill \\ \hfill {\dot{I}}_{d1}& ={\rho }_d{I}_1-A_2{I}_{d1},\hfill \\ \hfill {\dot{R}}_1& ={\gamma }_{u1}{I}_1+{\gamma }_{d1}I{d}_1-({\mu }_1+\zeta )R_1,\hfill \\ \hfill {\dot{S}}_2& =\zeta S_1-(1-{\Phi }_2){\widehat{\lambda }}_2{S}_2+\omega V_2-{\mu }_2{S}_2,\hfill \\ \hfill {\dot{V}}_2& =\zeta V_1-(1-\eta )(1-{\Phi }_2){\widehat{\lambda }}_2{V}_2-(\omega +{\mu }_2)V_2,\hfill \\ \hfill {\dot{E}}_2& =\zeta E_1+(1-{\Phi }_2){\widehat{\lambda }}_2{S}_2+(1-\eta )(1-\varphi ){\lambda }_2{V}_2-(\alpha +{\mu }_2)E_2,\hfill \\ \hfill {\dot{I}}_2& =\alpha E_2+\zeta I_1-A_3{I}_2,\hfill \\ \hfill {\dot{I}}_{d2}& =\zeta I{d}_1+{\rho }_d{I}_2-A_4{I}_{d2},\hfill \\ \hfill {\dot{R}}_2& =\zeta R_1+{\gamma }_{u2}{I}_2+{\gamma }_{d2}I{d}_2-{\mu }_2{R}_2,\hfill \end{array}$$

(13)

where

$${\widehat{\lambda }}_1={\displaystyle \frac{{\mu }_1({\beta }_{u1}{I}_1+{\beta }_{d1}{I}_{d1}+{\beta }_{u2}{I}_2+{\beta }_{d2}{I}_{d2})}{\Lambda }},{\widehat{\lambda }}_2={\widehat{\lambda }}_1.$$

Now, let us consider

$$J=F-V=\left(\begin{array}{cccccc}-(\alpha +{\mu }_1+\zeta )& {\beta }_{u_1}{M}_1& {\beta }_{d_1}{M}_1& 0& {\beta }_{u_2}{M}_1& {\beta }_{d_2}{M}_1\\ \alpha & -A_1& 0& 0& 0& 0\\ 0& {\rho }_d& -A_2& 0& 0& 0\\ \zeta & {\beta }_{u_1}{M}_2& {\beta }_{d_1}{M}_2^0& -({\mu }_2+\alpha )& {\beta }_{u_2}{M}_2& {\beta }_{d_2}{M}_2\\ 0& \zeta & 0& \alpha & -A_2& 0\\ 0& 0& \zeta & 0& {\rho }_d& -A_4\end{array}\right).$$

□

Assume ${\mathcal{R}}_0<1$, then the DFE is locally asymptotically stable, and the principal eigenvalue of J is negative. Notice that J is a quasi-positive matrix. Thus, by Perron-Frobenus theorem, J has a negative eigenvalue, say $-\lambda$$(\lambda >0)$ with a corresponding eigenvector strictly positive, U that satisfy $U^{T}J=-\lambda U^T$. Next, let Y denote the vector whose coordinates are the infected compartments, i.e., $Y={(E_1,I_1,I_{d1},E_2,I_2,I_{d2})}^T$ and define $H\left(Y\left(t\right)\right)=U^{T}Y\left(t\right)$. Based on (13),

$$Y\ge 0\Rightarrow \dot{H}\left(Y\left(t\right)\right)=U^T\dot{Y}\left(t\right)\le U^{T}JY\left(t\right)=-\lambda U^{T}Y\left(t\right)=-\lambda H\left(Y\left(t\right)\right).$$

Then $F\left(X\right(t\left)\right)\to 0$ and consequently $Y\left(t\right)\to 0$ as $t\to +\infty$. By invoking the limiting argument, one can notice that the DFE attracts all solutions of (13) in $\Omega$, completing the proof. The importance of this theorem in terms of the disease extension is that it guarantees the elimination of the disease whenever the basic reproduction threshold ${\mathcal{R}}_0$ is smaller than unity.

3.3. Numerical Results

The results from the model fitting demonstrate a strong agreement between the observed weekly numbers of confirmed cases (represented by red dots) and the predicted numbers of confirmed cases generated by the model (depicted as the blue curve; Figure 2(a)). Similarly, the predicted number of deaths aligns closely with the observed reported deaths (Figure 2(b)). The root mean square error ($RMS{E}_1$) for the model fit is 208.16. Additionally, the cross-validation process, performed on the remaining 10% of observations, yields an $RMS{E}_2$ value of 29.06, indicating the model’s robust predictive capability.

3.3.1. MPOX Dynamics in the Absence of Vaccination

The basic reproduction number ($R_0$) for MPOX was 1.73, while the control reproduction number ($R_c$) was 1.48. The peak in the number of confirmed cases occurred on September 29, 2024, with approximately 411 cases (Figure 2(a) and Table 4) and 4 reported deaths (Figure 2(b)). The detection rate of confirmed cases was low (0.0001), likely attributable to limited testing efforts.

Furthermore, we examine the evolution of the number of confirmed MPOX cases and deaths and the effective reproduction number over time in the absence of vaccination. The estimated model parameters, presented in Table 4, were used to predict these dynamics by December 2025. Figure 2 illustrates a declining trend in the number of confirmed MPOX cases and deaths, as well as a similar downward trend in the effective reproduction number during this period (Figure 2(c)). This highlights the gradual decline of Mpox transmission in the DRC. However, the projected increase in the effective reproduction number ($R_e$) starting around July 2025 raises concerns about a potential re-emergence of the disease.

3.3.2. Impact of Vaccination on MPOX Transmission Dynamics

We simulated the model to evaluate the impact of vaccination on MPOX dynamics in the DRC by predicting the effective reproduction number, the weekly numbers of confirmed MPOX cases and deaths under two scenarios: without vaccination and with increasing vaccination coverage levels (10%, 20%, 30%, and 40%). This simulation provides insights into the weekly trends of confirmed MPOX cases and deaths across these vaccination levels. Figure 3 illustrates the weekly incidence of MPOX-confirmed cases and deaths averted under different vaccination levels with the baseline scenario corresponding to no vaccination. The first two graphs (Figure 3(a) and (b)) of this figure depict the weekly number of confirmed cases and deaths in the absence of vaccination, showing a consistent decline throughout 2025. The weekly number of confirmed cases begins at approximately 300 in January 2025 and steadily declines to near zero by the end of the year. Similarly, the weekly number of reported deaths starts at around 3 in January 2025 and also diminishes to negligible levels by December 2025. The number of confirmed cases averted increases with higher vaccination coverage, peaking mid-year for all scenarios. Specifically, at 40% vaccination coverage, the weekly number of cases averted reaches a maximum of approximately 15 (3(c)). The weekly number of deaths averted shows a similar trend (3(d)).

These results emphasize that increasing vaccination coverage significantly reduces both the number of confirmed cases and deaths, particularly during the first half of 2025. The higher the vaccination level, the greater the averted burden of disease, highlighting the critical role of vaccination in controlling MPOX dynamics.

Figure 3. Weekly incidence of Mpox cases and deaths averted across combined vaccination and control measure (${\Phi }_1$) levels.

Figure 3. Weekly incidence of Mpox cases and deaths averted across combined vaccination and control measure (${\Phi }_1$) levels.

Combined Effect of Vaccination and Increasing Control Measures Levels (${\Phi }_1$)

Since under-15 individuals are the most affected by MPOX, non-pharmaceutical interventions targeting behavioral factors could be crucial in reducing MPOX transmission and its overall burden. In this study, we also examine the impact of varying levels of the reduction factor ${\Phi }_1$ of the force of infection ${\lambda }_1$ among under-15 individuals, in combination with different vaccination coverage levels, on the effective reproduction number, and the total number of cases and deaths averted by December 2025 (Figure 4, Figure 5 and Figure 6).

Figure 4. Effective reproduction number of MPOX across combined vaccination and control measure (${\Phi }_1$) levels.

Figure 4. Effective reproduction number of MPOX across combined vaccination and control measure (${\Phi }_1$) levels.

Figure 4 shows a decreasing trend in the effective reproduction number ($R_e$) across all scenarios, except for the baseline scenario, where an increase in $R_e$ is observed around July 2025, even with vaccination (Figure 4(a)). In scenarios with simulated improvements in control measure efficacy targeting the under-15s (Figure 4(b), (c), and (d)), the downward trend in $R_e$ becomes consistent with, however, a delay in the peak in MPOX cases and deaths ($R_e=1$) from January 2025 (Figure 4(b)) to June 2025 (for $40\%$ vaccination coverage) or later in September 2025 (for $10\%$ vaccination coverage) (Figure 4(d)). These findings emphasize the critical importance of adopting an integrated approach to Mpox control.

Figure 5 shows the trend in the total number of MPOX cases averted across varying vaccination coverage levels (10%, 20%, 30%, and 40%) and different control measure levels (${\Phi }_1$). The total number of MPOX cases without vaccination or control measures is 4602, serving as the baseline for comparison. With the baseline value of ${\Phi }_1$ (${\Phi }_1=3\%$), the number of cases averted increases with vaccination coverage, ranging from approximately 100 cases averted at 10% coverage to over 500 cases averted at 40% coverage. With an improved control measure efficacy, the number of cases averted increases substantially. For example, for ${\Phi }_1=10\%$, nearly 1000 cases are averted compared to the baseline scenario at 40% vaccination coverage and 1500 cases averted for an increased control measure efficacy (${\Phi }_1=15\%$) at the same vaccination level (40%). With the highest control measure efficacy considered (${\Phi }_1=20\%$), the number of averted cases reaches its maximum, surpassing 1800 cases at 40% vaccination coverage.

Figure 5. Total Mpox cases averted across combined vaccination and control measure (${\Phi }_1$) levels.

Figure 5. Total Mpox cases averted across combined vaccination and control measure (${\Phi }_1$) levels.

Figure 6 shows the total number of MPOX deaths averted across varying vaccination coverage levels (10%, 20%, 30%, and 40%) and different control measure efficacy levels (${\Phi }_1$=3%, 10%, 15%, and 20%). The total number of deaths in the absence of vaccination with baseline control measures is 51.

With the baseline control efficacy (${\Phi }_1$=3%), the number of deaths averted increases proportionally with vaccination coverage, ranging from approximately 1 death averted at 10% coverage to over 4 deaths averted at 40% coverage. With improved control efficacy, the number of deaths averted increases significantly. At 40% vaccination coverage, approximately 10, 15 and 17 deaths are averted for ${\Phi }_1$=10%, 15%, and 20%, respectively. These results highlight the synergistic effect of increasing vaccination coverage and improving control measure efficacy on reducing MPOX-related cases and deaths. The findings underscore the critical importance of combining vaccination campaigns with effective public health interventions to minimize MPOX transmission and burden.

Figure 6. Total Mpox Deaths averted across combined vaccination and control measure (${\Phi }_1$) levels.

Figure 6. Total Mpox Deaths averted across combined vaccination and control measure (${\Phi }_1$) levels.

The weekly numbers of under-15s that need to be vaccinated to achieve coverage levels of 10%, 20%, 30%, and 40% are provided as a supplementary file. A graph of these numbers for the first two months (January and February 2025) (Figure A1) shows an initial increase in the number of under-15s to be vaccinated in the first week of the campaign, followed by a steady allocation throughout the remainder of the vaccination campaign.

4. Discussion

This study uses a mathematical age-structured model to evaluate the impact of vaccination strategies on Mpox dynamics in the Democratic Republic of Congo (DRC), focusing on those under 15 years old as the primary target of upcoming vaccination campaigns. The model stratifies the total population into two subgroups according to age status (under-15 and over-15 year old subpopulation). It is fitted to weekly numbers of confirmed MPOX cases and deaths using the nonlinear least square technique. Our findings contribute to existing efforts to control the disease, with important implications for public health interventions.

Our study shows that vaccination can significantly reduce Mpox cases and deaths, particularly at higher coverage levels. For instance, at 40% coverage, the model predicts up to 1,800 cases and 17 deaths averted, emphasizing the critical role of vaccination campaigns in mitigating disease transmission and burden. These findings align with previous studies that highlighted the effectiveness of vaccination in reducing Mpox transmission during outbreaks ([21] and [22]). To be effective, the vaccination campaign will need to target under-15s in highly affected provinces in the DRC such as South-Kivu and Equateur [7].

Our analysis also highlights the synergistic benefits of combining vaccination with non-pharmaceutical interventions (NPIs). The study reveals a possible re-emergence of MPOX even with vaccination combined with no improved control measures targeting under-15s. Enhancing NPIs, such as public awareness campaigns, behavioral modifications, and improved hygiene practices, can address the potential rebound in $R_e$ and avert a re-emergence of Mpox in the population. Active surveillance systems should be expanded to rapidly detect and respond to new cases. Enhanced behavioral measures, such as reducing contact rates and improving awareness, mainly targeting under-15 individuals, further amplified the effectiveness of vaccination, supporting the conclusions from [23], who modeled the interplay between vaccination and behavioral interventions in zoonotic disease control. Similarly, [17] demonstrated that co-interventions, including vaccination, significantly reduce the overall burden of zoonotic diseases like Mpox.

Despite these encouraging results, our study underscores some challenges and gaps in current Mpox control efforts. First, the disproportionately low reduction of the force of infection linked to control measures (${\Phi }_1$=3.4%) among under 15-year-old individuals, as observed in this study, calls for targeted interventions focusing on this vulnerable age group. As noted by [23] and corroborated by recent findings [8], the need for age-specific severity estimates remains critical for optimizing vaccination strategies.

Another limitation is the lack of data on cross-immunity and LC16m8 vaccine waning immunity, which could influence long-term vaccination outcomes. Studies by [21] have highlighted the importance of understanding vaccine efficacy over time to ensure sustainable Mpox control. Future research should focus on addressing these gaps to improve model predictions and inform policy decisions.

Overall, this study emphasizes the need for integrated approaches combining vaccination with other interventions to effectively control Mpox in endemic regions. Public health efforts can maximize the impact of interventions and mitigate the disease’s burden by targeting the most affected populations and addressing existing data gaps. Our findings provide a valuable framework for policymakers and researchers to design and implement evidence-based Mpox control strategies.

5. Conclusions

This study developed and applied an age-structured mathematical model to evaluate the impact of vaccination on Mpox transmission and burden in the DRC. Our findings highlight the critical role of vaccination in reducing Mpox cases and deaths. Simulations revealed that increasing vaccination coverage from 10% to 40% substantially reduces disease transmission and burden, with the number of averted cases and deaths increasing significantly under higher coverage levels and enhanced non-pharmaceutical interventions. These results emphasize the need for targeted vaccination campaigns focusing on under-15-year-old individuals and the importance of integrating vaccination with other control measures to maximize public health outcomes. This work provides valuable insights into vaccination strategy design and its implications for controlling Mpox in endemic regions.