An epidemiological modeling investigation of the long-term changing dynamics of the plague epidemics in Hong Kong

1.   Introduction

Plague is a zoonotic bacterial disease caused by the bacterium Yersinia pestis (or Y. pestis, formerly known as Pasteurella pestis) that affects humans and other mammals ^[1,2]. The disease is described as a disastrous cause of the black death that transpired around the middle of the 13th century and caused the death of millions of humans ^[1,2,3]. In recent years, major outbreaks occurred between 2010 and 2015 and affected over 3240 individuals, including more than 600 deaths ^[1]. The disease is endemic in Africa, Asia, South America, and Europe, with Peru, Madagascar, and the Democratic Republic of the Congo as the most affected countries, recently ^[1].

Plague is a zoonotic disease usually found in small mammals and their fleas and transmits between animals ^[4,5,6,7]. Humans can also be infected either by the bite of infected vector fleas, contact with infected items (or body fluids), or the intake of infected respiratory droplets ^[1,4]. The early 20th-century plague epidemics in Hong Kong were part of a global epidemic, with the Bombay outbreak in the 1890s identifying fleas as vectors for Yersinia pestis ^[4,8]. This identification was essential for developing effective strategies for plague prevention and management by targeting flea vectors and their rodent hosts. Plague can be a severe disease in humans, with a case-fatality ratio (CFR) ranging between 30–60% for the bubonic type and most fatal for the pneumonic type, especially when left untreated ^[1,2]. Until now, there has been no specific vaccine against plague infection. However, effective treatment against the disease's infection is currently available ^[1].

The third historically-recorded pandemic in China was introduced through Yunnan city in the 1850s and spread widely for the next few decades; however, it reached pandemic proportions in Hong Kong in 1894, killing over 2500 people and forcing the evacuation of roughly one-third of the city's population ^[3]. At various times from 1894 to 1904, the disease was reintroduced to the city via marine trade, killing thousands of people ^[3,9]. At the turn of the 20th century, most research on plague transmission involving rats, fleas, and humans was conducted by European and Japanese scientists in Hong Kong or using Hong Kong data ^[3,4,5,10] with emphasis on analyzing the correlation between rat infestations and plague and the many types of plague in rats ^[3,4].

In past centuries, plague dynamics have shifted dramatically between pandemics and years of relative inactivity. The 14th-century Black Death, a horrific pandemic that raced throughout Europe, serves as a stark illustration of the explosive nature of plague dynamics ^[11]. Thus, exploring long-term changing dynamics in plague epidemics could offer an intriguing insight into the intricate interaction of infectious agents, hosts, and the environment over prolonged durations.

Several epidemic modeling studies have been utilized to explore the dynamics of plague and related disease transmission (see, for instance, ^{[4,5,6,7,12,13,14,15,16,17]} and the references therein). In particular, Dean et al. ^[4] investigated the dynamics of plague transmission in Europe during the during the Black Death. They found that human ectoparasites, like body lice and human fleas, were more likely to have caused the rapidly developing epidemics than rats, especially in pre-industrial Europe. They also revealed that ectoparasite infections transpire with high infectivity in humans, which was in line with an experimental study by Boegler et al. ^[14]. Keeling et al. ^[5] employed a mechanistic model to study the transmission dynamics of bubonic plague in human, rat, and flea populations. They used sensitivity analysis to find some important parameters that are needed to control the plague. These include the equilibrium number of rat cases and the rate of infection during the first epizootic cycle. According to their study, they suggested that relatively small rodent populations can cause severe plague epidemics, which highlighted that timely intervention for vector control is a significant strategy. Nguyen et al. ^[6] studied the transmission dynamics of bubonic and pneumonic plague in Madagascar. They found that the estimated transmission rate from person to person was higher than the flea-to-human transmission rate. This further showed that person-to-person transmission should also be prioritized to effectively mitigate plague outbreaks.

In this study, we aim to extend previous research, such as those considered in ^[4,5], by incorporating demographic information, including birth and death rates, to enhance our understanding of the dynamics of plague epidemics in Hong Kong during the 1902–1904 outbreaks. Accurate modeling of disease transmission, especially over extended outbreak periods, is crucial for predicting disease spread and burden. Incorrect estimates of key epidemiological quantities, such as reproduction numbers and serial intervals, can lead to inaccurate predictions and responses to epidemics and pandemics. Our model considers the interplay of human, flea, and rodent populations to infer the long-term changing dynamics of the plague epidemic in Hong Kong. By fitting the model to historical plague mortality data from 1902 to 1904, we aim to elucidate the underlying drivers of the dynamics and identify effective measures for future epidemic preparedness. Our findings contribute to a deeper understanding of plague dynamics and inform strategies for mitigating the impact of future outbreaks. The paper is arranged as follows: the epidemic model is presented in Sections 2 and analyzed in Section 3 and sensitivity analysis are conducted in Section 4 brief discussion of results is provided in Section 5.

2.   Materials and methods

2.1.   Surveillance data

In this study, we analyzed the data for plague-induced mortality during nine outbreaks in Hong Kong between 1902 and 1904 ^[18]. The data were made available to the public through secondary sources, including published articles, books, and government reports. We digitized the data from these sources to cover the entire outbreak period. Unlike other historical records, such as the Eyam report, which exhibited two fatality peaks, the Hong Kong data provided a continuous record of weekly mortality figures. This allowed us to model the entire epidemic period more comprehensively. Figure 1 shows the plague weekly time series cases from the surveillance data for Hong Kong between 1902 and 1904 using the data obtained from ^[18].

Figure 1.  Observed weekly number of plague death cases for human in Hong Kong between 1902–1904.

DownLoad: Full-Size Img PowerPoint

2.2.   Effective reproduction number

The effective reproduction number ($R_{eff}$) is an important epidemiological metric that measures the average number of secondary cases caused by an infectious individual in a partially susceptible population. In Figure 2, we estimate $R_{eff}$ based on the provided case data and utilize the serial interval (the time between successive cases in a transmission chain) for pneumonic plague, which was estimated from earlier studies to be around 2 to 4 days ^[19,20]. These estimations provide insights into the transmission potential of plague over time, incorporating uncertainty and smoothing to enhance clarity as well as aid in developing effective interventions.

Figure 2.  Estimated $R_{eff}$ for plague dynamics based on the time series data and serial interval information. The gray-dashed lines represent the 95% confidence interval (CI), the dotted line indicates the estimated $R_{eff}$, and the solid black line shows the estimated median $R_{eff}$.

DownLoad: Full-Size Img PowerPoint

2.3.   Spatial distribution of plague deaths cases in Kowloon and Island, Hong Kong

In Figure 3, the time series spatial analysis illustrates the transmission dynamics of plague in Kowloon and Hong Kong Island during the 1902–1904 epidemics. Both locations exhibited similar transmission patterns, with rodents playing a significant role in disease spread. The spatial distribution of the plague outbreak significantly impacted its transmission across different geographical regions in Hong Kong. These patterns likely arose from a combination of factors, including human ecological conditions, migration, and vector dynamics. Plague transmission often propagated through trade routes, human movement, and environmental factors, leading to localized and clustered outbreaks. Understanding these spatial patterns is crucial for developing effective surveillance, control strategies, and forecasting future epidemics. By identifying these patterns, targeted interventions can be devised to curb disease transmission and protect vulnerable populations. These findings highlight the importance of incorporating spatial analysis into epidemiological modeling efforts to derive actionable insights that inform public health strategies. To further explore the spatial heterogeneity of plague transmission, identify high-risk areas, and uncover the key drivers of disease dynamics, future work will incorporate advanced spatial modeling techniques. This approach aims to enhance our ability to design effective mitigation plans and targeted interventions tailored to specific geographic settings, thereby improving our capacity to prevent and manage plague outbreaks and strengthen preparedness against unexpected epidemics or pandemics.

Figure 3.  Spatial distribution of plague situations in Kowloon and Island in Hong Kong from 1902 to 1904.

DownLoad: Full-Size Img PowerPoint

2.4.   Plague transmission model

We employed a dynamic model to investigate the long-term changing dynamics of plague epidemics in Hong Kong from 1902 to 1904, drawing on insights from previous studies ^[4,5,15]. Our model, derived from ^[4], explores the interplay of human, flea, and rodent populations and incorporates demographic factors to infer shifts in plague epidemics while considering the population's homogeneous nature and the random distribution of vectors. Specifically, we extended the standard susceptible-exposed-infectious-removed (SEIR) epidemic model by integrating demographic information in humans and incorporating the logistic growth rate in the rodent population, a novel approach for analyzing Hong Kong's plague dynamics. This modification is particularly relevant for studying the relatively long period, during which demographic parameters play a crucial role in understanding plague transmission dynamics.

It is also critical to understand that the model considers different transmission pathways, rat to human and rat to rat via flea vector. Although the compartment of the flea vector was not explicitly present in the model, it is assumed to be included mechanistically as a reservoir that transfers the plague disease from dead rats to healthy rats or humans ^[4], we also noted that the higher the number of fleas in the population, the more plague infection. The ratio $\frac{\beta_{r}}{\beta_{h}}$ represents the rate of transmission from infected deceased rat (infected-dead rat or rat carcass) to susceptible rat and from infected rat carcass to susceptible humans, which is given by ^[4] based on the parameter values from Table 2. According to previous reports ^[15], the biting rate of an infectious flea is approximately 5 times more likely to cause an infection in susceptible rats than in susceptible humans. For more detail regarding the rat-flea model" formulation, see ^[4]. The total human population at time $t$, represented by $N_h(t)$, is split into sub-compartments of susceptible, $S_h(t)$, infected, $I_h(t)$, recovered, $R_h(t)$, and plagued-deceased, $D_h(t)$, humans, such that

The total flea (vector) population at time $t$, represented by $N_{f}(t)$, is split into two compartments: the average number of fleas living on a rat ($H_f(t)$) and the number of free infectious fleas ($F_f(t)$) that are searching for a host. Hence, we have

The total rat population at time $t$, represented by $N_{r}(t)$, is divided into sub-compartments of susceptible $S_r(t)$, infected $I_r(t)$, recovered, $R_r(t)$, and plague deceased, $D_r(t)$, rats. Thus,

These definitions and the incorporation of demographic parameters allow for a comprehensive analysis of plague dynamics over the studied period. For more detailed information regarding the rat-flea model formulation, please refer to ^[4].

The plague model is epidemiologically described in Eq (2.1). The state variables and parameters listed in Tables 1 and 2, respectively, all of which are assumed to be positive, satisfy the next set of non-linear ordinary differential equations (ODEs).

Table 1.  Interpretation of the state variables of the model (2.1).

Variable	Interpretation
$N_{h}$	Total population of humans
$S_h$	Susceptible humans with risk of plague infection
$I_{h}$	Infected humans with symptoms of plague infection
$R_h$	Recovered humans
$D_h$	Plague deceased-humans
$N_{f}$	Total population of fleas
$H_f$	Average number of fleas living on a rat
$F_f$	Number of free infectious fleas searching for a host
$N_{r}$	Total population of rats
$S_{r}$	Susceptible rats
$I_{r}$	Plague-infected rats
$R_r$	Recovered rats
$D_r$	Plague deceased-rats

---

 | Show Table

DownLoad: CSV

Table 2.  Table of summary for the model parameters given in model (2.1). The ranges for the transmission rates $\beta_h$ and $\beta_r$ were derived from ^[4]. All other parameter values were chosen based on our simulation process.

Parameter	Notation	Baseline value	Remark/Unit	Source(s)
Transmission rate for humans	$\beta_h$	$0.1001 (0.01-0.2)$	$\text{variable}$	[4]
Transmission rate for rodents	$\beta_r$	$0.5001 (0.02-0.75)$	$\text{variable}$	[4]
Average infectious period in humans	$1/\gamma_h$	$10 (2.38 - 66.67)$	$\text{day}$	[4]
Average infectious period in rats	$1/\gamma_r$	$5.2 (1.59 - 100)$	$\text{day}$	[4]
Humans recovery rate	$g_h$	$0.4 (0.01-0.8)$	$\text{dimensionless}$	[4]
Rats recovery rate	$g_r$	$0.1 (0.01-0.5)$	$\text{dimensionless}$	[4]
Flea growth rate	$r_f$	$0.0084 (0.0001-0.05)$	$\text{day}^{-1}$	[5]
Fleas carrying capacity	$K_f$	$6 (3.29-11.17)$	$\text{dimensionless}$	[4,5]
Flea searching efficiency	$a$	$0.004 (0-1)$	$\text{dimensionless}$	[5]
Natural death rate of fleas	$1/\mu_f$	$5.0 (3.23 - 100)$	$\text{day}$	[4]
Natural death rate of rats	$1/\mu_r$	$5.0 (4.0 - 50.0)$	$\text{day}$	[5]
Natural death rate of humans	$1/\mu_h$	$1/(35\times 365) (1/(32\times365) - 1/(45\times365))$	$\text{day}^{-1}$	[32,33]
Recruitment rate of humans	$\pi_h$	$42.27 (32.88 - 46.23)$	$\text{person \cdot day}^{-1}$	[32,33]

---

 | Show Table

DownLoad: CSV

It is worth mentioning that the susceptible flea population follows a logistic growth and is determined by the intrinsic growth rates, $r_f H_f(1-\frac{H_f}{K_f})$. To further gain more epidemiological insight into the plague dynamics, we computed a basic reproduction number ($R_0$), a key epidemiological metric that indicates the average number of secondary infections produced by one infected individual in a completely susceptible population It is used in evaluating disease transmission dynamics. For more details regarding $R_0$ computation, see Appendix A.

2.5.   Fitting framework

Following previous approaches ^{[21,22,23,24,25,26,27,28,29]}, we modeled the reported plague mortality cases in Hong Kong from 1902 to 1904 as a partially observed Markov process. The weekly reported cases, $\Psi_i$, were assumed to follow an over-dispersed Poisson process, modeled using a negative binomial (NB) distribution:

where $\tau$ is an over-dispersion parameter to be estimated. The likelihood function for the $i$-th week of the $j$-th year, $L_{i, j}(\cdot)$, represents the probability of the observed number of cases given the simulated cases $\Gamma_i$.

The overall log-likelihood, $l$, for the time series is:

where $\Phi$ is the total number of weeks in an epidemic wave, and $\Upsilon$ represents the parameter vector to be estimated. We used iterated filtering and plug-and-play likelihood-based inference to find the maximum likelihood estimates (MLE) of $\Upsilon$, implemented using the $\texttt{POMP}$ package in $\texttt{R}$ ^[30]. A fixed-time-step Euler-multinomial algorithm was used to simulate the ODE system.

To compare model performances, we used the Bayesian information criterion (BIC):

where $U$ is the number of data points, $V$ is the number of parameters, and $\hat{\Upsilon}$ denotes the estimated parameters ^[31].

For the fitting process, we calibrated the model using both human and rodent mortality data from 1902 to 1904. Although the human and rodent epidemics exhibited similar temporal trends, they differed in magnitude. To account for seasonal variations, we employed a cubic spline approach. The number of nodes in the cubic spline, $n_m$, was varied to find the optimal value that minimized the BIC. This approach allowed us to reconstruct the $\Gamma_i$ series, revealing that the 1902–1904 plague dynamics in human and rodent populations shared similar patterns with seasonal variations. It is important to note that while seasonal variations were accounted for, explicit fluctuations in rodent population size were not modeled directly ^[27].

3.   Results

3.1.   Long-term shifting dynamics of the plague outbreaks

In Figure 4, we present the model-fitting results for plague-induced mortalities in Hong Kong between 1902 and 1904, encompassing both human and rodent populations. Using biologically plausible parameters, our epidemic model effectively captures the temporal patterns of the plague epidemics during this period. The fitting results closely align with the surveillance data, indicating that our model can accurately assess the transmission dynamics of the plague in Hong Kong and potentially in other regions. The parameters used in the model are detailed in Table 2.

Figure 4.  Fitting results of the plague cases in Hong Kong from 1902–1904 for humans and rat populations using cubic spline functions. The black circles represent the reported number of plague death cases, and the red solid lines indicate the medians of the simulation. The grey-shaded regions denote the 95% confidence interval from the simulation. The blue dashed line represents the human transmission factor, and the green dashed line denotes rat susceptibility for the plague epidemics in Hong Kong from 1902–1904. The parameters used for the fitting results are summarized in Table 2. The inset panel shows the BIC as a function of $n_m$, the number of nodes in the cubic spline. The minimum BIC is attained at $n_m = 7$ in each case.

DownLoad: Full-Size Img PowerPoint

To account for seasonal variations in the transmission factor and rodent susceptibility (which refers to the likelihood of rodents becoming infected when exposed to Yersinia pestis, influencing disease spread among the rodent population), we employed a cubic spline function with uniformly distributed nodes over the study period. The inset panels in Figure 4(a) display the profiles as functions of the number of nodes ($n_m$). The optimal model, with the smallest BIC, was achieved with $n_m = 7$. Consequently, the main panels show the model fitting results with $n_m = 7$.

Figure 4 demonstrates that our model fits well with the reported plague cases from March 1902 to June 1904. The transmission factor ($R_0(t)$) and rodent susceptibility are time-dependent, as depicted by the right vertical axis for the human and rodent populations, respectively. The outbreak began around April 1902, peaked in July, and showed a declining trend, reaching control by June 1904.

These findings have significant public health implications. The accurate fit of the model underscores the importance of considering both human and rodent populations in understanding plague dynamics. The seasonal fluctuations indicate that increased rodent activity and elevated transmission rates during certain periods contribute to the epidemic's peaks. From a public health perspective, this implies that interventions such as rodent control, flea eradication, and public health campaigns should be intensified during late spring and early summer when the risk is highest. By incorporating seasonal dynamics into the model, public health authorities can better predict periods of increased risk and allocate resources more efficiently, leading to more timely and effective interventions. Implementing this preventive strategy can reduce the incidence of plague and mitigate its impact on public health ^[27]. Overall, our investigation highlights the importance of integrated strategies for pest control, continuous surveillance, and adaptable public health measures to effectively control and prevent plague outbreaks.

4.   Numerical simulations and sensitivity analysis

4.1.   Numerical results

As illustrated in Figure 5(a)–(c), our numerical simulations of the contour plots have uncovered a significant correlation between the $R_0$ and certain crucial epidemic parameters generated from the model. The observed trend indicates that an increase in both the transmission rate from rodents to humans ($\beta_r$) and the flea death rate ($\mu_f$) leads to a substantial rise in $R_0$ (see Figure 5. Similarly, in Figure 4(b), we found that an increase in $\beta_r$ and the rodent death rate ($\mu_r$) causes an increase in $R_0$, consequently increasing the epidemic's impact. However, Figure 5(c) revealed that an increase in $\mu_r$ and $\mu_f$ causes a decrease in $R_0$. This further demonstrates the critical epidemiological metric of these parameters in relation to $R_0$, which serves as an indicator of the potential for the plague to cause a serious epidemic. The results signify the influence of the transmission rate or contact rate of the disease, especially from rodents, fueling the likelihood of disease transmission per contact, logically contributing to an elevated $R_0$. Therefore, the observed increase in $R_0$ with rising $\beta_r$ and $\mu_f$, or $\beta_r$ and $\mu_r$, as well as the decrease in $R_0$ with rising $\mu_r$ and $\mu_f$, suggests a complex interplay between these epidemic parameters and the transmissibility of the plague in Hong Kong.

Figure 5.  Contour plots of the basic reproduction number ($R_0$, as a response function) in terms of key epidemic parameters of the model. (a): Contour plot showing the variation of $R_0$ with respect to the transmission rate from rodents to humans ($\beta_r$) on the x-axis and the death rate of fleas ($\mu_f$) on the y-axis. The plot indicates that $R_0$ will increase with an increase in $\beta_r$ and decrease in $\mu_f$. (b): Contour plot showing the variation of $R_0$ with respect to the transmission rate from rodents to humans ($\beta_r$) on the x-axis and the death rate of rodents ($\mu_r$) on the y-axis. The plot indicates that $R_0$ will increase with an increase in $\beta_r$ and decrease in $\mu_r$. (c): Contour plot showing the variation of $R_0$ with respect to the death rate of rodents ($\mu_r$) on the x-axis and the death rate of fleas ($\mu_f$) on the y-axis. The plot reveals that $R_0$ will decrease with an increase in both $\mu_r$ and $\mu_f$. All other parameter values used are given in Table 2.

DownLoad: Full-Size Img PowerPoint

4.2.   Sensitivity analysis

In this subsection, we performed a sensitivity analysis to determine the most sensitive parameter that should be prioritized in curtailing the spread of plague epidemics. We obtained the partial rank correlation coefficient (PRCC) for the sensitivity analysis by employing the same technique from previous studies ^{[21,22,34,35]}. In Figure 6, we show the results of the PRCC of the sensitivity analysis, which reveals key controllable epidemic parameters of the model (2.1). According to our results, the death rate of the flea vector ($\mu_f$) and infectious period in rats ($\gamma_r$) are the most sensitive parameters of the model that need emphasis for effective control of plague outbreaks. The PRCC of the basic reproduction number, $R_0$, was depicted in Figure 6 with the estimated parameters from the Table 2. The results further reveal that eliminating the plague vector by all means (e.g., adopting vector-borne disease control measures, such as the use of insecticides to kill the flea vector) is significant in mitigating plague epidemics.

Figure 6.  The partial rank correlation coefficient of the basic reproduction number, $R_0$, with respect to model (2.1) parameters. The dots represent the estimated correlation, and the bars are 95% CI. X-axis are the PRCC and Y-axis are the model's parameters. All the parameter values used are given in Table 2.

DownLoad: Full-Size Img PowerPoint

5.   Discussion

Infectious disease modeling has played a crucial role in understanding epidemics, providing an effective way of examining disease transmission and assessing disease control and prevention strategies ^[6]. Plague is described as a notorious and calamitous disease that affects humans and mammals ^[4]. We employed a deterministic model that includes human, flea, and rodent populations, incorporating demographic information to explore the dynamics of plague transmission in Hong Kong during the 1902–1904 outbreaks. The model was analyzed and fitted using plague mortality data from Hong Kong from 1902 to 1904. Our investigation into the long-term impact of the changing dynamics during epidemics revealed intriguing insights into the complex interactions that occur over extended periods between infectious organisms, hosts, and the environment. We highlighted that similar future outbreaks could be accurately analyzed using the same model framework.

To accurately capture the transmission dynamics of the plague during the 1902–1904 epidemic in Hong Kong, we utilized both reported human mortality data and raw rodent mortality data as part of our model fitting process. The human mortality data, systematically collected from the streets by public health authorities, provided reliable weekly counts of fatalities. In contrast, the rodent mortality data, though raw and less precise by contemporary standards, were gathered through systematic searches and inspections conducted during the outbreak. By employing a cubic spline function, we were able to calibrate our model using these historical records, which provided robust estimates of rodent mortality trends in relation to human cases. It is important to note that the spline function modeled seasonal trends rather than explicit fluctuations in rodent population size. This approach aligns with established methodologies in historical epidemiological modeling, where the integration of reconstructed data plays a crucial role in understanding epidemic dynamics ^[4,5]. The inclusion of both human and rodent data allowed us to derive meaningful insights into the interaction between rodent populations and the plague pathogen, thereby reinforcing the validity of our findings and their implications for public health interventions.

In Figure 4, we demonstrated the correlation between the model and the reported cases of the plague epidemics in Hong Kong, spanning from April 1902 to June 1904. From the model fitting results, we assessed the dynamics of human and rodent populations using a time-varying transmission component and the susceptibility of the rodent population. Our results suggest that the reported epidemics were driven by a combination of human and rodent outbreaks, with a transmission rate that varied over time between fleas and humans, fleas and rodents, and vice versa. The results also indicate that the peak of plague epidemics in Hong Kong occurred during the midyear. Additionally, we found that the human transmission factor reached its maximum level between April and June annually, followed by a rapid decline. The susceptibility of rodents also peaked during the same period each year, then significantly decreased. Our findings provide a comprehensive overview of the annual occurrences of plague infection in both human and rodent populations. The escalating severity of the plague was likely due to a rise in the human transmission rate and heightened susceptibility of rodents each April during the 1902–1904 plague epidemic. We anticipate that transmission rates may be similar across different situations or regions.

However, additional data, especially concerning the rodent population (such as population size and the occurrence of illnesses), could be highly beneficial for further understanding the dynamics of the plague. According to our investigation, recurring trends in the data suggest that the fundamental causes of the plague epidemic may reflect similar tendencies throughout the outbreak period, potentially due to seasonal factors causing fluctuations as a result of changing conditions. We aim to explore more dynamics related to seasonal effects in future work, as understanding these effects is also crucial in understanding plague dynamics ^[36].

The proposed approach exhibits comparable goodness-of-fit in terms of log-likelihood values. The observed transmission among humans can be attributed to either an immediate outbreak among rodents, resulting in time-varying spillover, a consistent predominance of plague transmission, or a combination of both mechanisms. Our epidemiological data imply that plague transmission does not occur at a consistent rate; instead, it exhibits a punctuated equilibrium, with outbreaks happening in clusters followed by periods of silence. Several factors influence these fluctuating dynamics, including environmental circumstances, host immunity, genetic immunity evolution, and the presence of reservoirs and vectors ^[37].

Further numerical simulations were conducted, as depicted in Figures 5 and 6, to demonstrate the influence of key model parameters on overall plague transmission dynamics. These parameters include the mortality rate of flea vectors and the rate at which rodents transmit the disease. Our contour simulations (Figure 5) highlight the critical interplay between biological and ecological factors in driving plague epidemics. Our findings show that an increase in $\beta_r$ and a decrease in $\mu_f$ cause a substantial increase in $R_0$ (see Figure 5(a)), indicating the need to target rodent-to-human transmission and manage flea populations. Interventions to reduce flea populations remain vital despite their high turnover. In Figure 5(b), the decrease in $R_0$ with higher rodent mortality ($\mu_r$) and a decrease in $\beta_r$ reveals that killing rodents might be a good strategy but may inadvertently increase human exposure to infected fleas, thereby emphasizing the need for integrated pest management strategies for effective control. Also, in Figure 5(c), the decrease in $R_0$ with increasing $\mu_r$ and $\mu_f$ underscores the importance of rodent and flea longevity in sustaining the epidemic. These results highlight the necessity for comprehensive control measures to reduce both flea and rodent populations, informing proactive public health interventions to manage and prevent future plague outbreaks.

Furthermore, our sensitivity analysis results, presented in Figure 6, identify the death rate of the flea vector ($\mu_f$) and the infectious period in rodents ($\gamma_r$) as the most critical parameters for controlling plague outbreaks. Targeting these parameters through interventions such as insecticide use is essential for effective epidemic control, highlighting the importance of vector management in public health strategies. Person-to-person transmission is crucial in the spread of a plague epidemic, particularly in cases of pneumonic plague, as illustrated in previous modeling studies. However, our present study primarily focuses on flea vector and rodent transmission dynamics relevant to the 1902–1904 plague epidemics in Hong Kong. In future studies, we intend to include this critical transmission mechanism (including person-to-person) to enhance the robustness and relevance of our model in both historical and contemporary contexts.

In summary, the main epidemiological findings from our study are as follows: (i) Estimation of the human transmission factor and rodent susceptibility during the plague epidemics in Hong Kong from 1902 to 1904. This finding highlights that the escalating severity of the plague was likely due to a rise in the human transmission rate and heightened susceptibility of rodents each April during the 1902–1904 plague epidemic. We anticipate that transmission rates may be similar across different regions. (ii) To our knowledge, no previous modeling study has analyzed plague transmission in Hong Kong using mortality data from 1902 to 1904 to identify the main transmission drivers during this period. Previous studies did not appropriately incorporate vital demographic factors, which are key to modeling the long-term dynamics of disease transmission. (iii) Additionally, through sensitivity analysis, we identified the key parameters that should be prioritized for effectively controlling plague transmission; this is crucial for future plague outbreaks and could provide valuable insights to public health practitioners and policymakers in designing suitable intervention measures. Furthermore, we also identified some crucial intervention measures that might be helpful for future epidemics.

Our investigation revealed that the long-term changing dynamics of plague epidemics reflect a complex interaction of variables that transcend temporal and spatial boundaries. Unraveling this complexity not only improves our understanding of historical pandemics but also enhances existing strategies for disease prevention and containment in a rapidly changing environment. Therefore, we emphasize, in particular, that controlling and eliminating flea and rodent vectors is crucial for preventing potential future plague outbreaks.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article. During the preparation of the final manuscript, the authors used grammar checker "QuillBot" in order to check and correct language errors and typos in the manuscript. After using this tool, the authors reviewed and edited the content as needed.

Acknowledgments

The authors are grateful to the handling editor and the anonymous reviewers for their invaluable comments and suggestions, which greatly contributed to improving the quality and clarity of our manuscript.

Conflict of Interest

The authors declare that there is no conflict of interest. Daihai He is an Associate Editor for the journal but had no role in the handling, reviewing, or decision-making process for this manuscript in any capacity.

Authors' Contributions

SSM and DH conceived the study and drafted the first manuscript. SSM, SZ, WM, AC, and DH carried out the investigation. All authors critically revised the manuscript and gave final approval for publication.