Data sources

We used the numbers of confirmed mpox cases by age and sex in both historical and current outbreaks from public sources^7,26,58,59. We used reported cases from 1 January to 18 August 2024 for the current outbreak in the DRC for its stable growth trend and data availability at the sub-provincial level (Kamituga and other health zones of South Kivu) during that period. For Burundi, which only reported cases from late July 2024 onwards, we used case data from 16 September to 27 October, roughly corresponding to the second half of the stable growth phase until its peak, to exclude potential influences of reporting bias or initial transients. Further details of the datasets are available in Extended Data Table 1. Weekly incidence data for the DRC and South Kivu were collected from public reports^17,18,60. Empirical contact survey data from Manicaland, Zimbabwe²⁴ were retrieved via the socialmixr R package⁶¹. Population age distributions in the DRC, Burundi and Manicaland from relevant years were collected from censuses (where available)^62,63 and World Population Prospects 2024⁶⁴. We also used synthetic contact matrix data published by Prem et al.⁶⁵.

Model structure

We assumed that the observed number of mpox cases by age or by age and sex follows a multinomial distribution with probabilities proportional to the elements of the dominant eigenvector of the next-generation matrix.

where \(v(K)\) is the normalized dominant eigenvector of the next-generation matrix K. We constructed the next-generation matrix as a product of transmissibility, β, age-dependent relative susceptibility, σ_a, and a contact matrix, C. We used age-disaggregated case data and contact matrices for the analysis of clade Ia (that is, C = {c_ab}), where c_ab represents the mean number of community contacts (conversational and/or physical, depending on the type of contact matrix used) an individual from age group b has per unit time with those from age group a. We used multiple contact matrices potentially relevant to mpox transmission (Extended Data Table 2). For each contact matrix, β was adjusted so that the dominant eigenvector of the next-generation matrix estimated for the historical Tshuapa dataset was equal to the previously estimated reproduction number of 0.82 (ref. ¹⁵).

For the analysis of clade Ib, we used age- and sex-disaggregated case data and expanded the contact matrix to account for heterosexual transmission dynamics in addition to community contacts. We stratified male and female populations into groups with high and low sexual activity. Following the approach described by Endo et al.²⁰, we further stratified cases in the high-activity groups by their route of exposure (community- versus sexually associated exposures) and allowed sexually acquired cases to have higher risks of onward sexual transmission because of a substantial variance in sexual behaviour⁶⁶. Namely, in a highly heterogeneous sexual network, transmissions through sexual contact tend to concentrate among individuals with the highest rates of sexual activity; this renders sexual exposure an indicator of even higher contact rates than the average among the high-activity groups. We assumed that this variance effect is negligible among low-activity groups (such that stratification by exposure route can be disregarded) and that their sexual contacts are reflected in the existing contact matrices, whose definition of contacts includes sexual contacts. Despite its simplicity, our approach of stratifying sexual contacts by two activity classes and two exposure routes could well capture the stationary growth pattern of mpox transmission (see Supplementary Tables 2 and 3 and Supplementary Figs. 10–12 for a proof-of-concept network model).

Let subscripts H and L denote high- and low-activity groups and M/m and F/f represent male and female cases with sexual (upper case) and community exposures (lower case), respectively. Vectors i_X represent the age-stratified number of mpox infections in category X. The generation-wise reproduction process of i_X is then described in a block matrix format (with each block an 8-by-8 matrix accounting for age group stratification) as:

where O and I are zero and unit matrices, respectively. S_MF and S_FM are block contact matrices representing heterosexual transmission from sexually acquired high-activity individuals. Similarly, Σ_MF and Σ_FM represent community-acquired transmission from high-activity heterosexual individuals. Again, note that we expect S_X to be larger than Σ_X because sexually acquired cases are more likely to have more sexual partners than the population average. P = (p_a) and Q = (q_a) are diagonal matrices representing the age-dependent proportions of high-activity individuals for males and females, respectively. C_X is a sex-specific contact matrix and is assumed to be half of the sex-aggregated contact matrix hereafter for simplicity. We define B_σ as a diagonal block matrix whose diagonal blocks are an 8-by-8 diagonal matrix of the age-specific susceptibility σ_a multiplied by transmissibility β.

Noting that for t ≥ 1, \({{\bf{i}}}_{{{\rm{m}}}_{{\rm{H}}}}=P({{\bf{i}}}_{{{\rm{m}}}_{{\rm{H}}}}+{{\bf{i}}}_{{{\rm{m}}}_{{\rm{L}}}})\) and \({{\bf{i}}}_{{{\rm{f}}}_{{\rm{H}}}}=Q({{\bf{i}}}_{{{\rm{f}}}_{{\rm{H}}}}+{{\bf{i}}}_{{{\rm{f}}}_{{\rm{L}}}})\), the above relationship can be simplified into

where \({{\bf{i}}}_{{\rm{m}}}={{\bf{i}}}_{{{\rm{m}}}_{{\rm{H}}}}+{{\bf{i}}}_{{{\rm{m}}}_{{\rm{L}}}}\) and \({{\bf{i}}}_{{\rm{f}}}={{\bf{i}}}_{{{\rm{f}}}_{{\rm{H}}}}+{{\bf{i}}}_{{{\rm{f}}}_{{\rm{L}}}}\).

We assumed proportionate mixing to model S_MF, S_FM, Σ_MF and Σ_FM; further details on parameterization are described in the next section.

We constructed a multinomial likelihood for the observed number of male and female cases by age group:

where \([{{\bf{i}}}_{{{\rm{M}}}_{{\rm{H}}}}+{{\bf{i}}}_{{\rm{m}}};{{\bf{i}}}_{{{\rm{F}}}_{{\rm{H}}}}+{{\bf{i}}}_{{\rm{f}}}]\) is a vector obtained by sex-wise aggregation of the normalized dominant eigenvector of the next-generation matrix in equation (3).

Model fitting

Our model fitting process was twofold. We first fitted and validated the model of community contact transmission using the age distributions of cases from the historical and current clade Ia outbreaks in the DRC. We used this model to estimate age-dependent susceptibility parameters σ_a for the community contact matrix. We then constructed an expanded contact matrix including sexual contacts and calibrated it to the age–sex distributions of cases from clade Ib outbreaks in the DRC and Burundi to estimate additional parameters that describe transmission dynamics of clade Ib through community and sexual contact routes. All of the parameters used in the model are summarized in Extended Data Table 3.

Analysis of clade Ia outbreaks

As in the ‘Model structure’, we defined the next-generation matrix K = {k_ab} = βσ_ac_ab. The contact matrix c_ab was constructed from existing datasets (Extended Data Table 2), with adjustment for age structures between the source and modelled populations where different (using the density correction method described by Arregui et al.⁶⁷). Based on data availability, we used six age groups (0–4, 5–9, 10–19, 20–29, 30–39 and 40+ years) for the 2011–2015 Tshuapa data and eight age groups (0–4, 5–9, 10–14, 15–19, 20–29, 30–39, 40–49 and 50+ years) for the 2024 DRC endemic provinces data. When the age groups in the case data were coarser than the contact matrix, the elements were merged accordingly⁶⁸. Transmissibility β is a nuisance scaling parameter and was not estimated. We assumed that children aged 0–4 years may have potentially higher susceptibility to infection and that adults born before cessation of the smallpox vaccine in 1980¹¹ may have lower susceptibility. To account for this, we parameterized σ_a as

Here, ε represents the net vaccine protection among the vaccinated cohorts (that is, the product of the vaccine effectiveness and coverage in the historical campaign, also referred to as the effective vaccine coverage). The susceptibility for the age group that contains both vaccinated and unvaccinated cohorts (ages 30–39 years in the historical Tshuapa dataset and ages 40–49 years in the 2024 outbreak datasets) was specified as the weighted average between 1 − ε and 1, based on the proportion of individuals born in or before 1980 in those age groups⁶⁴ as of 2013 (the midpoint of the Tshuapa dataset period) or 2024. Given the high historical vaccination coverage in the DRC (89–97%)²⁹, the estimated vaccine protection ε would provide an approximation (or at least a reasonable lower bound) for the effectiveness against clade I mpox. We estimated σ₀ and ε using the Bayesian importance sampling method⁶⁹, based on the multinomial likelihood in equation (1), and obtained median estimates and 95% credible intervals.

Exploratory model development was done using the historical Tshuapa dataset. The resulting model, described above, was then validated using the 2024 DRC endemic provinces dataset. We estimated two parameters, σ₀ and ε, from the historical Tshuapa dataset alone and compared the model outputs (after accounting for a decade’s shift in the age of the vaccinated cohort) with the age distribution from the (hold-out) 2024 DRC endemic provinces dataset (Supplementary Fig. 1). The model outputs showed good overall concordance with the validation data. We obtained final parameter estimates from the joint estimation using both datasets for the subsequent analyses.

Analysis of clade Ib outbreaks

We parameterized the block components representing sexual contacts among high-activity individuals Σ_MF, Σ_FM, S_MF and S_FM in equation (3) as follows. We assumed that the age groups 15–19, 20–29, 30–39 and 40–49 years potentially contain high-activity populations (that is, p_a = q_a = 0 outside these age groups). We assumed a simple proportionate mixing assumption (that is, no age assortativity) where each high-activity individual makes sexual contacts at a specified rate over the infectious period of MPXV, which are randomly assigned to other high-activity individuals of the opposite sex, irrespective of age. We denote the mean sexual contact rates among high-activity females and males infected through community contact as v_F and v_M; similarly, we denote the mean sexual contact rates among those infected through sexual contact as w_F and w_M. Note that v_F, v_M, w_F and w_M are defined to have the same scale as the community contact matrix in relation to onward transmission (that is, one unit of v_F, v_M, w_F or w_M represents the amount of sexual contact that contributes to the transmission equivalent of one daily community contact). The entries of Σ_MF, Σ_FM, S_MF and S_FM are then modelled as:

where n_a and m_a are the relative population sizes for males and females by age group (we assumed even sex distribution^70,71; that is, \({\sum }_{{{a}}}{n}_{{{a}}}={\sum }_{{{a}}}{m}_{{{a}}}=1/2\)). Note that S_X and Σ_Χ are proportional, with a ratio of w_X:v_X. We assumed that infection through community contact has no association with one’s sexual contact rate (that is, v_F and v_M represent the mean sexual contact rates among the high-activity groups). Reciprocity of heterosexual contacts then requires:

We estimated the ratio between w_F and v_F assuming that the sexual network among high-activity groups is represented by a configuration network (that is, no degree assortativity; an alternative scenario is discussed in Supplementary Fig. 7), in which v_X corresponds to the mean degree and w_X to the mean neighbour degree (or degree-weighted average of degree)³⁰. We assumed that FSWs represent the behaviour of the majority of high-activity female individuals and used FSW survey data from the DRC⁷². Since both the mean degree and mean neighbour degree can be characterized by the mean μ and standard deviation, σ, of the sexual contact degree distribution, π(x), their ratio is estimated as:

We used the reported mean and s.d. for the weekly number of clients (mean = 17.61; s.d. = 12.0) and non-paying partners (mean = 5.51; s.d. = 22.0) per FSW⁷² to derive w_F/v_F = 2.17. Here we assumed that the numbers of clients and non-paying partners are uncorrelated (an alternative scenario is discussed in Supplementary Fig. 7); the sexual contact degree as a sum of these two numbers would then have a mean of 17.61 + 5.51 = 23.12 and a variance of 12.0² + 22.0² = 628, which were supplied to equation (8). This and the reciprocity requirement in equation (7) assure that v_M and v_F are derived from w_F, {p_a} and {q_a}. We also set boundary conditions for the total number of high-activity individuals (that is, \({\sum }_{{{a}}}{n}_{{{a}}}{p}_{{{a}}}\) and \({\sum }_{{{a}}}{m}_{{{a}}}{q}_{{{a}}}\)). Namely, we assumed that high-activity males account for 10% of the male population among the sexually active age groups and high-activity females account for 1% of the total female population (translated into ~2.2% among the sexually active age groups), based on engagement data on commercial sex in Africa (the proportion reporting paid sex in the past 12 months among males aged 15–59 years in Mozambique⁷³ and the estimated number of FSWs per capita in the DRC³⁶).

Assuming the community transmission patterns are similar between clades Ia and Ib, we used the same σ_a estimated for clade Ia throughout the clade Ib analysis. This implicitly assumes uniform smallpox vaccine coverage across geographical settings (see Supplementary Fig. 7 for a sensitivity analysis). The rest of the parameters, {p_a}, {q_a}, w_F and w_M, were estimated using the Markov chain Monte Carlo method (No-U-Turn Sampler⁷⁴). We obtained 2,000 Markov chain Monte Carlo samples while ensuring an effective sample size of at least 500 for every parameter; most parameters achieved an effective sample size of ~1,000 or larger. Using these posterior samples, we estimated the proportion of R_eff attributable to sexual transmission as a measure of the relative contribution of sexual contacts to mpox transmission. We defined it as the relative reduction in the effective reproduction number (R_eff) by excluding sexual contacts from the next-generation matrix. Namely, this proportion was given as 1 − λ_C/λ_SC, where λ_SC and λ_C are the dominant eigenvalues of the next-generation matrix (equation (3)) with and without the sexual contact blocks (S_X and Σ_X), respectively.

Model averaging and outputs

We used model averaging over candidate models using different contact matrix data (Extended Data Table 2) with Watanabe–Akaike weights⁷⁵ based on the combined likelihood for datasets 2, 3 and 4 in Extended Data Table 1. We did not include the likelihood for sub-provincial datasets (datasets 5 and 6) or Burundi (dataset 7) in the weights for the main analysis due to overlapping cases and potentially different data collection and contexts, respectively. Alternatively, weights also using these likelihoods were included in our sensitivity analysis (see Supplementary Fig. 7). Models 5 and 6 using synthetic contact matrices had substantially worse likelihood than other models, resulting in effectively zero weights, and were therefore not included in the main analysis.

Frequency of transmission between age groups

We estimated the relative frequency of infector–infectee age pairs at equilibrium using the next-generation matrix and normalized dominant eigenvector:

We then aggregated the block elements of M to obtain a sex- and transmission-route-aggregated frequency matrix.

Projection of time evolution of the reproduction numbers

We projected the effective reproduction number (R_y) of mpox in the DRC and Burundi by shifting the age of the smallpox-immunized cohorts and modifying the contact matrices⁶⁷ according to the population projection⁶⁴. Referring to the previous estimate¹⁵, we assumed that the effective reproduction number of clade Ia in the Tshuapa province was 0.82 (95% confidence interval = 0.79–0.85) in 2015 (the midpoint of the estimation period 2013–2017). We then projected the effective reproduction number for clade Ia between 2010 and 2030 as the dominant eigenvalue of the next-generation matrix for each year, accounting for the following two factors. First, the relative susceptibility of age groups that include the immunized cohorts was adjusted according to the estimated net protection, ε (model averaging weights from Table 1), and the proportion of immunized cohorts (those born in or before 1980, per the population projection⁶⁴) within the age groups. Second, the community contact matrices were also adjusted to the projected population age distribution using the density correction method⁶⁷. For clade Ib, we projected the effective reproduction number in 2024 onward assuming that patterns of community contact transmission are shared between clades Ia and Ib. Namely, the R_y for clade Ib was projected by incorporating the additional sexual contact contributions described in equation (3) (estimated for the Kivus and Burundi) into the projected next-generation matrices for clade Ia, assuming time-invariant sexual contact patterns. Here we assumed a hypothetical pre-outbreak state and neglected potential infection- or mpox-vaccination-derived immunity, as well as any future interventions. To yield an uncertainty range for these projections, we considered uncertainty in the reference estimate¹⁵ (both the confidence interval and the estimation period being for 2013–2017 instead of the midpoint, assuming normally and uniformly distributed errors, respectively) in addition to the uncertainty in our parameter estimation.

Effective reproduction numbers under focused and mass vaccination scenarios

We modelled the change in the effective reproduction number of clade Ib under a combination of FSW-focused and mass vaccination scenarios. We used an mpox vaccine effectiveness estimate of 86%⁷⁶, which we assumed to protect against infection but not against onward transmission in the case of breakthrough infections. FSW-focused vaccination was assumed to be offered to female high-activity individuals of 20–39 years of age. For mass vaccination, all individuals aged 20–39 years, regardless of sexual activity level, were included. The baseline effective reproduction numbers without vaccines were matched to the projected values reflecting the ageing of the smallpox-immunized cohorts for each population setting.

Ethics and inclusion

This study used publicly available data only and did not require ethics approval. Per the Global Code of Conduct for Equitable Research Partnerships, local researchers were included in the study design, study implementation, data ownership and authorship of publications. The local relevance of the research has been discussed with and guided by local researchers.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.