Abstract
Despite reporting very few mpox cases in early 2023, Mainland China observed a surge of over 500 cases during the summer. Amid ambiguous prevention strategies and stigma surrounding mpox transmission, the epidemic silently escalated. This study aims to quantify the scale of the mpox epidemic and assess the transmission dynamics of the virus by estimating the effective reproduction number (Re) during its early phase. Publicly available data were aggregated to obtain daily mpox incidences in Mainland China, and the Re was estimated using an exponential growth model. The mean Re value was found to be 1.76 (95% credible interval: 1.51–2.06), suggesting a case doubling time of approximately 2 weeks. This estimate was compared with Re values from 16 other countries’ national outbreaks in 2022 that had cumulative incidences exceeding 700 symptomatic cases by the end of that year. The Re estimates for these outbreaks ranged from 1.17 for Portugal to 2.88 for Colombia. The pooled mean Re was 1.66 (95% credible interval: 1.40–1.92), which aligns closely with the Re for Mainland China. These findings underscore the need for immediate and effective control measures including targeted vaccination campaigns to mitigate the further spread and impact of the epidemic.
1. Introduction
On 9 August 2023, the China Centers for Disease Control and Prevention (CDC) announced 491 cases in Mainland China in July, representing a marked escalation from 106 cases in June and a five-fold increase (1). Although concerning, this surge was largely anticipated by health experts (2, 3). Amid ambiguous prevention strategies and the stigmatization of specific patient groups—especially men who have sex with men (MSM)—the mpox epidemic intensified unnoticed.
Mpox is a viral infection mainly transmitted through direct skin-to-skin contact (4). The global mpox outbreak of 2022–2023 was primarily driven by sexual transmission within MSM networks (5). With limited spillover to the broader general population, the 2022 outbreaks predominantly affected high-risk MSM communities (6, 7), making these communities primary targets for subsequent vaccination campaigns (8, 9).
The global transmission of mpox was evident in 2022, but it was not until 2023 that Western Pacific nations, including Japan, South Korea, and Taiwan, saw significant case counts (10-13). By the end of 2022, many of these countries had minimal or no local mpox transmissions, possibly due to ongoing coronavirus disease 2019 (COVID-19)-related travel and social restrictions (14). Given the considerable MSM population in Mainland China, the uncontrollable spread of mpox presents a significant public health concern. An effective mpox vaccine exists, however, challenges such as limited accessibility and disease-related stigma have hampered its broad dissemination (2), underscoring the need to understand mpox transmission dynamics in Mainland China and predict the evolution of future outbreaks.
Epidemics typically follow a consistent transmission cycle: introduction, escalation, peak, and eventual decline (15). The initial stages, often marked by exponential growth, offer important insights into understanding of disease spread dynamics. It is crucial to estimate the effective reproduction number, Re, which indicates the average number of secondary transmissions per primary case, during the period in which interventions, behavioral shifts, or acquired immunity have not yet made an impact.
Prior Re evaluations for the 2022 national mpox outbreaks have revealed variability (16-19). Using a sub-exponential growth model, Du and colleagues (16) estimated that the Re ranged from 1.02 in Portugal to 1.95 in the United States. Conversely, other reports have noted higher Re values. For example, for the United Kingdom, Re values of 2.32 and 3.14 were reported (17, 18), compared with a value of 1.18 estimated by Du et al. (16). Whereas these studies (17, 18) applied variations of a simple susceptible–exposed–infected–recovered model, they also relied on case counts that were classified by the date of reporting. Such an approach can result in an upward bias in estimates because it may reflect not only actual transmission dynamics but also the recognition of an outbreak and the discovery of infection clusters. Administrative batching in case reporting can further bias the estimation of Re. Nevertheless, a recent
UK mpox study (19) employed a detailed, more complex susceptible–exposed– infected–recovered model that factored in potential biases, generating Re values that ranged from 5.16 to 3.02 during the initial stage of the outbreak. Given diverse factors, including cultural differences, differences in outbreak response, and challenges with data collection, identifying how these estimates align with Re in Mainland China is crucial.
In this study, we analyzed data from the World Health Organization (WHO) and other governmental sources to enhance our understanding of mpox transmission in Mainland China (20-22). By focusing on case counts solely attributed to Mainland China and estimating the Re during the initial phase of the outbreak, we aimed to compare counts with estimates from 2022 and project potential case counts for the near future. Our ultimate goal was to inform ongoing public health interventions, and potentially optimize efficiency of these interventions, while also contributing to global epidemic preparedness.
2. Methods
Data collection
This study used three primary data sources. First, daily incidences of confirmed mpox cases in Mainland China, Taiwan, and the Hong Kong Special Administrative Region were extracted from the WHO global mpox dashboard (2022-2023) as of 15 August 2023 (20). These incidences were categorized by the date of symptom onset, followed by the date of diagnosis and the date of reporting. Second, a dataset of the daily incidences that was confirmed in Taiwan, categorized either by the date of symptom onset or the date of reporting, was retrieved from the National Infectious Disease Statistics System of the Taiwan Centers for Disease Control (CDC) on 23 August 2023 (21). Third, a line list of cases that were confirmed in Hong Kong was compiled using press releases from the Center for Health Protection of the Department of Health of the Government of Hong Kong Special Administrative Region on 23 August 2023 (22). The WHO data recorded 859 total cases in Mainland China, Taiwan, and Hong Kong. These cases were classified as follows: 615 by the date of symptom onset, 225 by the date of reporting, and 19 by the date of diagnosis. The datasets from Taiwan and Hong Kong recorded 292 and 33 cases, respectively; these cases were all classified by the date of symptom onset. Notably, one patient in Hong Kong lacked a definitive date for symptom onset. The report indicated that this patient began showing symptoms at the end of June. The patient was reported in August 2023 and was unlikely to have been accounted in the WHO dataset; therefore, the patient was omitted from our analysis.
To extract cases that were confirmed solely in Mainland China, the combined daily counts from the Taiwan and Hong Kong datasets were subtracted from the incidences recorded in the WHO dataset. In the process, two assumptions regarding data decomposition in the WHO dataset were investigated.
First, we assumed that all counts from Taiwan were reported by the date of symptom onset in the WHO dataset. However, upon subtracting these counts from the WHO data, the resulting daily incidence yielded predominantly negative values before 1 June 2023. Subsequently, only one date (21 June 2023) appeared with a negative incidence, when the WHO dataset recorded a daily count of two cases, whereas Taiwan and Hong Kong reported three and no cases, respectively. Therefore, we deduced that our initial assumption regarding the integration of Taiwan incidence data into the WHO dataset was likely different.
Thus, we hypothesized that all case counts from Taiwan were incorporated into the WHO dataset by their respective dates of reporting to the WHO. Upon examining the WHO reporting pattern, it became evident that most cases that were defined by their reporting dates—specifically May 2, 9, 16, 23, 30; June 6, 20, 27; and July 4—were recorded at equidistant intervals of 7 days. A single deviation from this pattern was noted: the inclusion of June 8 with the exclusion of June 13 (Supplementary Figure 1). Assuming that the reporting date in Taiwan had a one-day lag to its confirmation date, the case counts in Taiwan aligned almost perfectly with the counts of cases by reporting date in the WHO dataset (Supplementary Figure 2). Although incidences from the Taiwan dataset predominantly fell below the WHO dataset counts, matches were observed for the last five dates: June 8, 20, and 27, and July 4 and 11. This observation led us to the conclusion that the latter approach of imputing data from Taiwan to the WHO dataset was plausible.
Moreover, the upper bound of the exponential growth phase was set to tup = 5 july 2023, indicating the peak of observed daily incidence of symptomatic cases at 21. The exponential growth phase was investigated by considering a 45-day period and placing a lower bound at tlo= 21 May 2023. The robustness of the estimates was explored in the sensitivity analysis by varying the dates for both boundaries, which ranged from 16 to 26 May 2023 for the lower bound and 30 June to 10 July 2023 for the upper bound. The duration of the exponential growth phase that was considered in the sensitivity analysis was thus varied between 15 and 55 days.
Statistical framework
Reporting delay distribution
Given that all cases with unidentified symptom onset dates were back-projected to their estimated date of symptom onset in the subsequent analysis, estimating the reporting delay distribution was crucial. This distribution measured the period from the date of symptom onset to either the date of reporting or the date of diagnosis. For this estimation, we utilized the Hong Kong line list data, which is composed of 33 records, and a cutoff time of T = 23 August 2023.
Each of the reporting delays (time to diagnosis or time to reporting) was fit to the generalized gamma distribution (GGD), which is defined by three parameters: shape (Q), location (M), and scale (S) .The choice of GGD was guided by its flexibility and its capability to represent three commonly used distributions: gamma (when Q = S), Weibull (when Q = 1), and log-normal (when Q = 0) (23). The mathematical inference for the GGD was based on the following: if the random variable follows a GGD, then the transformation where Q exp (Q − 2z), Where z = (log(x) − M)/S, adheres to a gamma distribution with shape Q − 2 and scale 1.
In particular, Oi represented the date of symptom onset, and Di (where Oi≤ Di) was the date of reporting or diagnosis, extracted from the Hong Kong line list data, with I= 1,…,33. We assumed that the respective times of symptom onset (oi) and reporting or diagnosis (di) were uniformly distributed within their intervals as follows:
Where the symbol “∼” implies “is distributed as.” The difference di − oi, representing the observed delay, conforms to the following truncated gamma distribution:
where: yi = Q exp (Q −2 Zi), Yi = Q exp (Q−2Zi), and zi = (log(di − oi) − M)/S.Zi=(log(T−oi)−M)/S.
The likelihood is represented by the form:
where λ and Г are the probability density function and cumulative density function of the gamma distribution, respectively. The parameters are assumed to have weakly-informative priors: log Q, M, log S ∼ Normal(mean = 0,SD =1).
Reconstructed incidence
Given that not all of the cases were categorized by their date of symptom onset, it was essential to back-project the cases identified solely by their date of diagnosis or reporting to ascertain their dates of symptom onset.
Therefore, nt was the number of cases that were symptomatic on day t, counting from day 1 on tlo. Similarly,
and 
represented the cases with unknown symptom onset dates but with identified reporting and diagnosis dates, respectively. Any non-zero count 
(where ∗: = {r,d}) was then back-projected using sampling from a multinomial distribution as follows:
In the above, 
describes the discretized reporting delay distribution as follows:
The final equation in set (5) indicates the likelihood of a case being reported on day but manifesting symptoms before day 1 (i.e., before tlo). To compute the overall incidence, the following summation was employed:
Effective reproduction number (Re) during the initial phase of the epidemic Earlier, we defined the exponential growth phase using the date range tlo ≤ t ≤ tup. The Re for the initial phase of the epidemic was estimated by analyzing the case counts (ct) with the presumption of the exponential growth. The growth rate was estimated by incorporating the negative binomial likelihood as follows (24):
where tlo ≤ t ≤ tup and the parameters { i0,r, ϕ } were supported by the following weakly-informative priors:
Re was computed based on Wallinga and Lipsitch’s formula (25), which links Re,r, and gamma-distributed generation time with the mean μand SD σ:
For this study, we adopted the values μ=12.5 days and σ=5.7 days, as estimated in the study by Guzzetta et al. (26).
Extrapolation of case counts beyond 5 July 2023
Acknowledging the potential for variations in reported cases due to potential delays or modifications in reporting, we extrapolated the case counts from 5 July to 1 September 2023 under the assumption of sustained exponential growth. The case counts 
, which encompass both previously reported and extrapolated cases based on symptom onset date followed the following equation:
Here, “∼”denotes sampling from the negative binomial distribution. The cumulative counts of cases for July and August were derived by aggregating 
values across their respective date ranges.
For computing case counts based on the date of reporting (nt), the counts 
were adjusted by considering the reporting delay. For every time, the following sampling was conducted:
where 
refers to (5). Thereafter:
Finally, the aggregated count of the number of cases (nt) for July and August was determined using summation over the respective date ranges.
Global trends across the 2022 outbreaks
Mpox incidence data for 2022 was collated from the WHO dashboard (20). From this dataset, we selected 16 countries that confirmed over 700 symptomatic cases by the end of the year. Most of these nations provided daily incidence counts categorized by the date of symptom onset, the date of diagnosis, or the date of reporting. However, six counties in the WHO Region of Americas (AMRO)—Argentina, Brazil, Chile, Colombia, Mexico, and Peru—presented only aggregated case counts without any subtyping. For each country, we identified the exponential growth phase as a 45-day window that led up to the epidemic’s peak,which was represented by the highest daily case count for 2022. In situations in which multiple peaks were recorded, we selected the earliest occurrence. In our sensitivity analysis, we considered alternating windows of 30 and 60 days.
For simplicity, we applied the reporting delay distribution from the Hong Kong data—a reasonable approach given its close alignment to a globally reported median of 6 days, an inter-quartile range spanning 4–9 days (20), and the absence of more detailed global data. Using the methods described earlier, we then estimated the Re for each of the selected countries.
2.1. Technical details
We adopted the Bayesian framework to fit each realization of the daily incidence classified by the date of symptom onset (6) to models (7)–(9) using Markov Chain Monte Carlo (MCMC) sampling techniques. The Bayesian estimation process was performed using Stan software (27). The inherent structure of the reconstructed incidence was challenging given its reliance on sampling back-projected counts from a multinomial distribution (4). Because these counts, which were represented as integers, were not immediately compatible with the efficient use of Hamilton Monte Carlo sampling in Stan, we employed an approach similar to that previously described (28). For each realization of (4), we simulated the MCMC chain using 1,000 iterations for the tuning-in process and retained only a single posterior for further analysis. To construct the posterior distribution, we used 4,000 different chains that resulted in 4,000 posterior draws. The convergence of the MCMC simulation was inspected visually and was checked implicitly by fitting the model only for symptomatic cases. In doing so, the back-projection was not needed. One thousand iterations that were used for tuning-in and four chains that were composed of 1,000 posterior draws resulted in an R-hat statistic below 1.1 (29). The code is made available on http://github.com/aakhmetz/Mpox-in-MainlandChina-2023.
3. Results
Using the Hong Kong data, we estimated a mean delay of 3.4 days (95% credible interval [CrI]: 2.7–4.4) from symptom onset to case diagnosis and a mean delay of 6.0 days (95% CrI: 5.1–7.0) from symptom onset to case reporting. Upon integrating these delays into the case back-projection (Figure 1), the Re was estimated at 1.76 (95% CrI: 1.49–2.06). Correspondingly, the doubling time was 15 days (95% CrI: 11–21). Projecting the incidence of mpox past 5 July 2023 to cover the months of July and August, the mean count of cases that were categorized by the symptom onset date was 1,080 (95% CrI: 650–1,720) for July and 5,200 (95% CrI: 1,820–12,060) for August.
The exponential growth phase (dashed lines) represents the lower bound on 21 May 2023 and the upper bound on 5 July 2023.
Adjusting these numbers for the reporting delay, the mean numbers of reported cases were 800 (95% CrI: 540–1,180) and 3,850 (95% CrI: 1,490–8,290) for July and August, respectively.
The analysis of national outbreaks from selected countries in 2022 revealed a range in Re values. Portugal had the lowest mean Re value at 1.17 (95% CrI: 1.07–1.30), whereas Colombia had the highest value at 2.88 (95% CrI: 2.46–3.34), as illustrated in Figure 2. The pooled mean Re was 1.71 (95% CrI: 1.48–1.95), with a between-country variance (heterogeneity) of τ 2=0.19. However, upon evaluating the data quality provided to the WHO by individual countries, it became evident that the nations within the WHO AMRO—with the exceptions of the United States and Canada—provided only aggregated case counts without subtyping them by date of symptom onset, date of diagnosis, or date of reporting. Such reporting likely introduced higher variability and potential bias in the Re values. By excluding these countries from the analysis, the pooled mean Re dropped slightly to 1.66 (95% CrI: 1.40–1.92), and the between-country variance declined to τ2=0.15.Notably, the Re, value for Mainland China remained similar to the pooled average regardless of whether AMRO countries were included in the analysis.
The Re values are shown in the first column on the right, and a period of 45 days is indicated in the last column. Countries that did not report case counts subtyped by the date of symptom onset, diagnosis date, and reporting date were excluded from the calculation of the pooled mean and are shown in dashed black.
The sensitivity analysis confirmed the robustness of our findings but highlighted some variability. When considering only cases with a known symptom onset date, the Re for Mainland China remained consistent at 1.75 (95% CrI: 1.48–2.05), suggesting that most cases classified by date of reporting were predominantly linked to Taiwan. Varying the time window for the identified exponential growth phase—in which the lower bound (tup) varied between 16 May 2023 and 26 June 2023, and the upper bound (tlo) varied between 30 June 2023 and 10 July 2023—the median posterior Re ranged from 1.58 to 1.72 (95% CrI 1.38–2.05; Supplementary Figure 3). The posterior predictions for reported cases ranged from 650 to 840 (95% CrI: 480–1,400) for July and 2,150 to 3,540 (95% CrI: 1,110–9,690) for August (Supplementary Figure 4, Supplementary Figure 5). For the selected 2022 outbreaks, adjusting the length of the exponential growth phase to either 30 or 60 days did not critically alter the prior estimates (Supplementary Figure 6). Importantly, when we focused solely on cases with a known date of symptom onset, the estimated Re remained unchanged, indicating the minimal impact of case records with unknown symptom onset on the Re for the selected countries.
4. Discussion
In our analysis, we estimated that the Re of the mpox epidemic in Mainland China was within the range of 1.6–1.8. This range suggests a doubling time of approximately 2 weeks. Importantly, our estimate was based on the early phase of the epidemic, in which an evident exponential growth of cases was observed and control measures had likely not yet been implemented. Given the relatively low incidence of mpox and no deaths reported during the epidemic’s early phase, significant behavioral changes in at-risk populations were likely not yet in effect, implying that mpox transmission occurred largely unhindered within affected communities during this early phase.
The Re value estimated in our study closely aligns with values derived from selected national outbreaks in 2022, suggesting the potentially universal character of the outbreaks. However, an accurate estimation of infections remains challenging given the likely underreporting of mpox cases. Several factors might contribute to this under-ascertainment. One reason may be the societal stigmatization associated with the disease that leads to individuals’ reluctance to disclose their infection for fear of facing societal prejudice or discrimination. Another possibility is that infections can result from intimate encounters outside of a committed relationship, leading individuals to withhold their health status because of the fear of repercussions in their personal lives.
Consequently, most cases in the dataset may be either severe or detected through active surveillance measures such as those implemented at entry points into Mainland China.
Our study has several limitations. First, we did not differentiate between indigenous and imported infections in our dataset, a factor that may have influenced the estimation of Re. However, by focusing exclusively on the exponential phase in which case growth is most evident, the influence of imported cases on our estimates was likely minimal. Second, the actual method of imputing cases from Taiwan in the WHO dataset may have differed from ours or could have been modified by the WHO over time. For example, the observed pattern of data imputations by date of reporting that were aggregated every 7 days was no longer feasible at the end of July and August. Third, the reporting delay distribution was estimated using Hong-Kong line list data, which may differ from the data from Mainland China given the differences in healthcare systems, case ascertainment, and case investigation protocols. Fourth, the results of comparisons of Re values across different countries must be interpreted with caution.
The proactive control measures adopted by some nations versus the more passive strategies employed by others can introduce variances. For instance, the outbreak in Portugal evolved at a slower pace given the government’s proactive position. The case counts were almost steady over 2–3 months before a decline was observed. Although the Re values in our study and that of Du et al. (16) were close to 1, a granular look at the case counts combined with a time-varying Re revealed a brief Re peak of 2.7 on 10 May 2022 that then gradually declined to 1 and dropped below this threshold by the end of June 2022 (30). Our mathematical framework successfully captured this average trend across Portugal and other countries, and a pooled mean Re of 1.66 summarized our overall findings. Finally, the WHO dataset retrieved on 15 August 2023 contained incomplete case counts given that some cases were yet to be reported. At the time of writing the results of this study, a snapshot for 12 September 2023 became available.
The Re estimate only slightly changed to 1.63 (95% CrI: 1.43–1.86); however, the case counts for May and June were updated and revealed possible long reporting delays (Supplementary Figure 1A and Supplementary Figure 7). Although mpox may be associated with a lower mortality rate compared with diseases such as smallpox, its morbidity is notably high, and mpox is associated with health complications such as pneumonia, confusion, and eye infections (31). Mpox disproportionally affects people living with human immunodeficiency virus (HIV) (32). The transmission of mpox throughout Mainland China in the summer of 2023 was largely unconstrained given that regulatory guidelines were first published only at the end of July. This timeline put vulnerable populations at heightened risk of infections and subsequent complications. Nevertheless, the potential impact of efficient vaccination campaigns and information programs aimed at raising awareness of the disease cannot be overlooked; such interventions could significantly reduce case numbers and prevent a larger future outbreak. In contrast to Mainland China, Taiwan managed to stabilize its mpox situation in 2023, maintaining an Re value of approximately one. Drawing from its previous experience—specifically, its management of the hepatitis A virus (HAV) outbreak among MSM in 2015–2017 (33)—the Taiwan government promptly initiated a vaccination campaign, targeted at high-risk groups that was bolstered by awareness campaign. The strategy not only ensured low vaccine hesitancy but also prevented an escalating outbreak. Mainland China could adopt a similar approach with vaccination campaigns focused on high-risk social groups coupled with the dissemination of accurate health messages. Avoiding stigmatization of the disease, an oversight during the early stages of the oubreak (2, 3), could further contribute to higher vaccination rates.
5. Conclusions
In this study, the mean Re value estimated for the early phase of the mpox outbreak in Mainland China was between 1.6 and 1.8, closely aligning with the pooled mean Re from selected outbreaks in 2022. This alignment suggests that Mainland China could experience a similar trajectory of case growth as observed in other countries previously affected by mpox. These findings underscore the urgency for implementing immediate and effective control measures, including targeted vaccination coupled with information campaigns to reduce stigmatization of vulnerable groups associated with the disease.
Data Availability
All data produced are available online at the github
Financial support
This study was supported by the National Science and Technology Council, Taiwan (NSTC #111-2314-B-002-289).
Conflict of interest
The authors declare that there are no known competing financial or personal relationships that could have appeared to influence the work reported in this paper.
Ethical standards
This study was based on publicly available data and did not require ethical approval.
Data availability statement
All data used for this study can be found at: https://github.com/aakhmetz/Mpox-in-MainlandChina-2023/tree/main/data
Figures
Daily incidences reported by the World Health Organization (A), Taiwan Centers for Disease Control (B), and the Government of the Hong Kong Special Administrative Region (C). The snapshot dates are indicated in the legends. Each case count is subtyped by the date of symptom onset, date of diagnosis, and date of reporting.
The concordance in case counts between the World Health Organization dataset (blue; as of 15 August 2023) and the Taiwan Centers for Disease Control dataset (pink; as of 23 August 2023) that were aggregated by date of reporting. The dates of reporting for Taiwan were aggregated to the dates of reporting from the WHO dataset. Taiwan cases beyond the latest date (11 July 2023) were omitted.
Effective reproduction number (Re ) during the initial phase of the epidemic in Mainland China, depending on the time window that identified the exponential phase. (A) shows the posterior median Re value, and (B) and (C) show the lower (2.5th percentile) and upper (97.5th percentile) bounds of the 95% credible interval, respectively. The horizontal axis indicates the lower bound, and the vertical axis indicates the upper bound that were used for the exponential phase window. The date is indicated as a month/day in 2023.
(A) shows the posterior median of number of cases, and (B) and (C) show the lower (2.5th percentile) and upper (97.5th percentile) bounds of the 95% credible interval, respectively. The horizontal axis indicates the lower bound date, and the vertical axis indicates the upper bound that were used for the exponential phase window. The date is indicated as a month/day in 2023.
(A) shows the posterior median of number of cases, and (B) and (C) show the lower (2.5th percentile) and upper (97.5th percentile) bounds of the 95% credible interval, respectively. The horizontal axis indicates the lower bound date, and the vertical axis indicates the upper bound that were used for the exponential phase window. The date is indicated as a month/day in 2023.
The Re values are shown in the first column on the right, whereas a period of 30 days (A) or 60 days (B) is indicated in the last column. Countries that did not report case counts subtyped by the date of symptom onset, diagnosis date, or reporting date, were excluded from the calculation of the pooled mean and are shown in dashed black.
Each case count is subtyped by the date of symptom onset, date of diagnosis, and date of reporting.
Acknowledgements
We thank Anahid Pinchis from Edanz (www.edanz.com/ac) for editing a draft of this manuscript.
Footnotes
E-mails: s8111127{at}gmail.com (P.W.)
Correcting the typographical mistakes with figures' references (Word issue) + updating the abstract
References
