Construction of effective reproduction number of infectious disease individuals based on spatiotemporal discriminant search model: take hand-foot-mouth disease as an example

Abstract

Objective

In order to facilitate the tracing of infectious diseases in a small area and to effectively carry out disease control and epidemiological investigations, this research proposes a novel spatiotemporal model to estimate effective reproduction number(Re)for infectious diseases, based on the fundamental concept of contact tracing.

Methods

This study utilizes the incidence of hand, foot, and mouth disease (HFMD) among children in Bishan District, Chongqing, China from 2015 to 2019. The study incorporates the epidemiological characteristics of HFMD and aims to construct a Spatiotemporal Correlation Discrimination of HFMD. Utilizing ARC ENGINE and C# programming for the creation of a spatio-temporal database dedicated to HFMD to facilitate data collection and analysis. The scientific validity of the proposed method was verified by comparing the effective reproduction number obtained by the traditional SEIR model.

Results

We have ascertained the optimal search radius for the spatiotemporal search model to be 1.5 km. Upon analyzing the resulting Re values, which range from 1.14 to 4.75, we observe a skewed distribution pattern from 2015 to 2019. The median and quartile Re value recorded is 2.42 (1.98, 2.72). Except for 2018, the similarity coefficient r of the years 2015, 2016, 2017, and 2019 were all close to 1, and p <0.05 in the comparison of the two models, indicating that the Re values obtained by using the search model and the traditional SEIR model are correlated and closely related. The results exhibited similarity between the Re curves of both models and the epidemiological characteristics of HFMD. Finally, we illustrated the regional distribution of Re values obtained by the search model at various time intervals on Geographic Information System (GIS) maps which highlighted variations in the incidence of diseases across different communities, neighborhoods, and even smaller areas.

Conclusion

The model comprehensively considers both temporal variation and spatial heterogeneity in disease transmission and accounts for each individual's distinct time of onset and spatial location. This proposed method differs significantly from existing mathematical models used for estimating Re in that it is founded on reasonable scientific assumptions and computer algorithms programming that take into account real-world spatiotemporal factors. It is particularly well-suited for estimating the Re of infectious diseases in relatively stable mobile populations within small geographical areas.

Introduction

The effective reproduction number (R_e) is an essential indicator for evaluating the scale and predicting the trend of infectious disease outbreaks [1]. Researchers use it to characterize the transmission behavior of a disease, while policymakers use it to develop effective mitigation strategies [2,3,4]. Unlike the basic reproduction number (R₀), R_e is better suited to practical situations as it represents the average number of secondary cases that may be produced by the primary case during its infectious period, and can therefore more accurately quantify a disease’s transmission capabilities [5, 6]. Common methods for obtaining the effective reproduction number include contact tracing [7], infectious disease transmission models and algorithms [8, 9], and estimation through exponential growth rates [10].

In recent years, effective mathematical modeling has become an influential method and scientific basis for the study of infectious disease management [11,12,13]. Such models are instrumental in simulating epidemic spread and predicting development trends. By constructing mathematical models, we can gain a better understanding of the characteristics and rules governing epidemic transmission, and provide a practical basis for scientific prevention and control. Such models also allow for timely formulation of response plans, swift action, and ultimately effective reduction of the harm caused by epidemics to protect public health and safety [14]. When it comes to obtaining R_e in the mathematical modeling of infectious disease dynamics, the most commonly used method is the classical effective reproduction number prediction based on SEIR models [15,16,17]. However, in recent years, there have been numerous research findings based on Bayesian inference methods using infectious disease outbreak dynamics [18, 19], the construction of infectious disease transmission matrices [20], and computer simulations to predict R_e [21]. For example, in recent years, there have been many research achievements in obtaining Re through the construction of differential equations and network mathematical models, such as approaches incorporating degree correlation [22, 23], weighted network models [24], and delayed differential equation models [25, 26]. However, due to the variability of infectious parameters and the diversity of network topologies, it is challenging to obtain accurate values that match the actual situation. Overall, there are several prevalent issues in the current methodologies for constructing R_e. For instance, one common challenge is the uncertainty associated with parameter estimations. The parameters utilized in the effective reproduction number model are derived through estimation procedures, which can be susceptible to uncertainties and errors. Another issue lies in the simplification of model assumptions. Many models rely on simplified assumptions, such as steady contact rates or fixed rates of epidemic spread. These assumptions may deviate from the actual scenario, potentially compromising the accuracy of R_e estimations. Moreover, a significant limitation is the inadequate consideration of spatial and temporal dynamics in disease transmission models. The transmission dynamics of actual epidemics are often influenced by spatial factors like geography and population movements, which are not comprehensively integrated into the existing models.

Moreover, the contact tracing method [27], which is the most basic method used by primary disease control personnel to obtain Re, is suitable for community or regional infectious disease outbreaks. This method involves identifying possible contacts, tracing the possible sources and transmission chains of infection, and interrupting disease transmission as quickly as possible to minimize Re. In this paper, the fundamental principle of contact tracing was incorporated into spatiotemporal mathematical modeling of an infectious disease. Using hand, foot and mouth disease (HFMD) in children as a case study, we explored a spatiotemporal infectious disease model based on computer programming and simulated infectious disease dynamics to calculate the R_e. This model was built upon the specific temporal and spatial positioning of each individual case and utilized a search calculation method, thus accommodating temporal variability and spatial imbalance.

Data and methods

Data sources

The data of HFMD cases in Bishan District of Chongqing from May 2015 to December 2019 were selected that including gender, age, onset time, diagnosis time, residential address, preschool institution address and contact information. The data source was Bishan District Center for Disease Control and Prevention (China Information System for Disease Control and Prevention [28]).

Study design

Development of a spatial-temporal database

We utilized ARC ENGINE and C# for the programming and development of HFMD spatio-temporal database and analysis software-PIADS. This system is equipped with a spatio-temporal database feature that can overlay various factors’ properties, enabling advanced analysis and search. The software’s analysis functions provide a solid foundation for detailed examination of the spatio-temporal factors surrounding HFMD.

For now, we mainly included attribute data for basic information, temporal information, and spatial location information in the spatio-temporal database, and will add information for infectious disease R_e impact factors in the future.

Algorithm and program implementation for spatiotemporal correlation discrimination of HFMD

In the case of HFMD, we proposed a time-related judgment function for HFMD. To compare the time correlation between two cases, the onset time of the i-th case is denoted as $\:{t}_{i1}$, the time of diagnosis as $\:{t}_{i2}$ the incubation period as$\:{\:t}_{q}$ and the period of post-onset infection as $\:{t}_{g}$. Similarly, for the j-th case, its onset time is denoted as $\:{t}_{j1}$. The time correlation discriminant function of the i-th case can be expressed as follows:

$$\:\left\{\begin{array}{c}{t}_{j1}-\left({t}_{i1}-{t}_{q}\right)\le\:{t}_{q}+{t}_{g}\exists\:{er}_{t}\left\{\right\}\\\:{t}_{j1}-\left({t}_{i1}-{t}_{q}\right)>{t}_{q}+{t}_{g}\:\nexists\:{er}_{t}\left\{\right\}\end{array}\right.\:$$

(1)

further simplification results in

$$\:\left\{\begin{array}{c}{t}_{j1}-{t}_{i1}\le\:{t}_{g}\exists\:{er}_{t}\left\{\right\}\\\:{t}_{j1}-{t}_{i1}>{t}_{g}\:\nexists\:{er}_{t}\left\{\right\}\end{array}\right.\:$$

(2)

$\:{\:\:\:er}_{t}\left\{\right\}$ refers to the time-dependent set of cases, with $\:t_g=1\;\left(\;\mathrm{week}\;\right)$under natural conditions where isolation measures are not implemented immediately after diagnosis of HFMD. In cases where isolation measures are taken immediately after diagnosis of HFMD, $\:{t}_{g}={t}_{2}-{t}_{1}$, where $\:{t}_{2}$ is the time of diagnosis and $\:{t}_{1}$ is the time of onset. Taking into consideration the infectious incubation period of the child, the time window of infectivity for the first generation of cases is $\:{t}_{g}+{t}_{q}$. The average incubation period of HFMD is 3–5 days, and patients typically seek medical attention 1–2 days after onset of symptoms. Therefore, the window period for infectivity in the first generation of cases is approximately 1 week.

In addition, a spatial correlation discriminant function is proposed. To compare the spatial correlation between two cases, the one with the earlier onset time is designated as the first-generation case, and the one with the later onset time is considered susceptible. Let$\:({x}_{i},{y}_{i})$ represent the spatial coordinate of the i first-generation case and ($\:{x}_{j},{y}_{j}$) represent the spatial coordinate of the susceptible person. $\:\xi\:$ denotes the distance of influence of infectious dynamics. The discriminant function for assessing spatial correlation between these two cases is as follows:

$$\:\left\{\begin{array}{c}\sqrt{{\left({x}_{j}-{x}_{i}\right)}^{2}+{\left({y}_{j}-{y}_{i}\right)}^{2}}\le\:\xi\:\exists\:{er}_{r}\left\{\right\}\\\:\sqrt{{\left({x}_{j}-{x}_{i}\right)}^{2}+{\left({y}_{j}-{y}_{i}\right)}^{2}}>\xi\:\nexists\:{er}_{r}\left\{\right\}\end{array}\right.$$

(3)

As shown above, $\:{er}_{r}\left\{\right\}$ refers to the distance-dependent collection of cases.

When both temporal correlation and spatial correlation criteria are simultaneously met, it is considered that there is a high likelihood of an infection chain between the two. The relationship between the first and second generation cases is thus determined as follows:

$$\:\left\{\begin{array}{c}{t}_{j1}-{t}_{i1}\le\:{t}_{g}\:\exists\:\:{er}_{t}\left\{\right\}\:and\sqrt{{\left({x}_{j}-{x}_{i}\right)}^{2}+{\left({y}_{i}-{y}_{i}\right)}^{2}}\le\:\xi\:\exists\:{\:er}_{r}\left\{\right\}\Rightarrow\:{\rm\:A}\{{B}_{1},\,{B}_{1}\cdots\:{B}_{\:n}\}\\\:{t}_{j1}-{t}_{i1}>{t}_{g}\:\nexists\:\:{er}_{t}\left\{\right\}\:or\sqrt{{\left({x}_{j}-{x}_{i}\right)}^{2}+{\left({y}_{j}-{y}_{i}\right)}^{2}}>\xi\:\:\nexists\:\:{er}_{r}\left\{\right\}\nRightarrow\:{\rm\:A}\{{B}_{1},\,{B}_{1}\cdots\:{B}_{\:n}\}\end{array}\right.$$

(4)

$\mathrm A\left\{B_1,\,B_1\;\cdots\;B_\eta\right\}$ denotes the collection of second-generation cases derived from case A of the first generation. $\:{B}_{\:n}$ represents the NTH second-generation case arising from case A of the first generation, and n corresponds to the number of second-generation cases originating from case A of the first generation.

By formulating a spatio-temporal correlation discriminant function and integrating it with the spatio-temporal database and analysis system for HFMD, it is possible to develop a search algorithm and program for analyzing the spatio-temporal correlation of infectious diseases. Essentially, this method entails contact tracing through time-space correlation searches, enabling the identification of both primary and secondary cases and subsequently calculating the average infection period (Re) of the primary cases. These represent pivotal stages in the establishment of a spatio-temporal correlation search model.

Implementation of the SEIR Model R_e algorithm and procedure

In order to verify the scientific validity of this method, we introduce the classic SEIR model R_e calculation formula for comparison. The formula is shown as follows:

$$\:{R}_{e}=\left(1+\frac{r}{{b}_{1}}\right)\left(1+\frac{r}{{b}_{2}}\right)$$

(5)

Where r is the growth rate, $\:{b}_{1},\,{b}_{2}$ are the removal rates of the latent population and the infected population, respectively. The incubation period and the infection period can be respectively referred to as: $\:{T}_{E}=\frac{1}{{b}_{1}}$, $\:{T}_{1}=\frac{1}{{b}_{2}}$ Thus, we can derive the following formula:

$$\:{R}_{e}=\left(1+r{T}_{E}\right)\left(1+r{T}_{1}\right)$$

(6)

Based on the formula above, we programmed and developed on the HFMD spatiotemporal database. Assuming that cases are removed (i.e. recovered) one week after infection, we calculated the SEIR model’s effective reproduction number (Re) by dividing the daily increase in cases by the number of people currently infected (i.e. cumulative infections minus removed cases).

Result

Determining the optimal search radius for a search model

Based on the aforementioned spatiotemporal search model, we selected cases of Hand, Foot, and Mouth Disease during its peak season (April to June) and used spatial discriminant values of distances 0.5 km, 1.0 km, 1.5 km, 2.0 km, 2.5 km, 3.0 km, and 3.5 km. Through computer calculations, we obtained specific R_e values for each case under the prescribed conditions. Taking the mean of these R_e values yields the following:

From Fig. 1, it can be observed that the data for each year shows the fastest increase around a radius of 1.5 km. As the search radius continues to increase, the growth rate of R_e begins to decline. The optimal critical value for the search model is the point at which the search radius reaches the maximum R_e while ensuring that the search range for first-generation cases does not overlap. When R_e reaches a certain threshold, the search ranges begin to overlap, resulting in a slower increase in R_e as the radius goes beyond 1.5 km.

Fig. 1

R_e values obtained with different search radius during the peak periods from 2015 to 2019

In addition, we can analyze the search model principle and determine that when the search radius is half the average distance between first-generation cases, the search range for first-generation cases can achieve maximum coverage with minimal overlap, in accordance with the optimal search radius for the model. For example, during the peak season from April to June, we can obtain the number of first-generation cases for the i-th infectious day, denoted as ni, and the distances between each pair of them, using computer programming. The average distance for the i-th day can be calculated as:

$$\:\stackrel{-}{{D}_{i}}=\frac{\sum\:_{j}^{{c}_{ni}^{2}}{d}_{j}}{{c}_{ni}^{2}}$$

(7)

In the equation, $\:\stackrel{-}{{D}_{i}}$represents the mean distance between first-generation cases for day i, $\:{c}_{ni}^{2}$ represents the sample size of the distances between first-generation cases for day i, $\:{d}_{j}$ represents the jth distance between first-generation cases for day i, and $\:ni$ represents the number of first-generation cases for day i.

The average value of the mean distance on peak days is:

$$\:\stackrel{-}{\stackrel{-}{{D}_{i}}}=\frac{\sum\:_{i}^{m}\sum\:_{j}^{{c}_{ni}^{2}}{d}_{j}}{m}$$

(8)

Where$\:\stackrel{-}{\stackrel{-}{{D}_{i}}}$is the average value of the mean distance between first-generation cases during peak days. m represents the number of days during the peak period.

The average distance between first-generation and secondary cases from 2015 to 2019 during the months of April to June is illustrated in the figure below. (Fig. 2). The cumulative distribution function of distance between first-generation and secondary cases from 2015 to 2019 during the same period is represented in the figure below (Fig. 3).

Fig. 2

Average distance between first-generation and secondary cases during peak periods from 2015 to 2019

Fig. 3

Cumulative distribution function of the distance samples between first-generation and secondary cases during peak periods from 2015 to 2019

The optimal search radius is half the mean distance between first-generation cases on a daily basis during the peak period.

$$\:{Radius}_{opt}=\frac{\sum\:_{i}^{m}\stackrel{-}{{D}_{i}}}{2 m}$$

(9)

Where $\:{Radius}_{opt}$ is the optimal search radius.

By calculation, We obtained the average distance and optimal search radius for first-generation cases during the peak period from 2015 to 2019 (Fig. 4). The optimal search radius is approximately 1.5 km, aligning with the prior conclusion. Therefore, it is advisable to suggest 1.5 km as the optimal search radius for spatiotemporal search models.

Fig. 4

Average distance between first-generation cases and optimal search radius during peak periods from 2015 to 2019

The R_e value obtained from a spatiotemporal data search model

The figure below (Fig. 5A) reveals that the fluctuation of the R_e value of HFMD in children, ranging between 1.14 and 4.75 from 2015 to 2019. Additionally, Fig. 5B illustrates a skewed distribution of the R_e reconstructed by the search model, with a median of 2.42 and a quartile range of (1.98, 2.72). The figure presents a scatterplot with statistical indicators, where the central line represents the median and the upper and lower quartiles are indicated respectively.

Fig. 5

(A) Monthly bar chart depicting R_e values from 2015 to 2019, generated using spatiotemporal data search model (B) Scatter plot of monthly R_e values based on search model

The R_e value obtained through programming calculation based on the SEIR model

Based on the aforementioned SEIR model and programming algorithm, daily R_e values can be obtained. By aggregating the daily values within a month, we obtain the monthly Re value. The specific results are as follows (Fig. 6A):

Fig. 6

(A) Bar chart of monthly R_e values from 2015 to 2019 based on SEIR model and programming algorithm (B) Scatter plot of monthly Re values based on SEIR model

The scatterplot with statistics was utilized to demonstrate (Fig. 6B). It can be observed that the SEIR model for HFMD in children had Re values fluctuating between 1.48 and 4.34 from 2015 to 2019, indicating a normal distribution with a mean and standard deviation of 2.88 ± 0.56.

Comparison of R_e values between search model and SEIR model

The two sets of curves in the following figure respectively represent the R_e values of spatiotemporal Search models and SEIR model (Fig. 7).

Fig. 7

Trend graph of R_e values of spatial-temporal search model and Re calculation based on SEIR model through computer programming from 2015 to 2019

As shown in the figure, during the period from 2015 to 2019, the two sets of curves exhibited similar trends and roughly matched the “bimodal” trend of HFMD incidence, with the first peak occurring from April to June and the second peak appearing later, from October to December, with the second peak lower than the first peak. Except for 2018, the curve changes in other years showed a trend of simultaneous increase and decrease.

We set the R_e values obtained by the two methods at the same point in time on a monthly basis as the data for correlation analysis (Fig. 8). The results showed that the 2015 search model did not follow a normal distribution and Spearman correlation coefficient was used, while Pearson correlation coefficient was used for other data that followed a normal distribution. Except for 2018, the similarity coefficient r of the years 2015, 2016, 2017, and 2019 were all close to 1, and p < 0.05, indicating that the R_e values obtained by using the search model and the traditional SEIR model are correlated and closely related (Table 1).

Fig. 8

Linear regression analysis graph of the correlation between R_e values of search model and SEIR model from 2015 to 2017

Table 1 Correlation analysis of monthly re-value of hand, foot and mouth disease from 2015 to 2019

Differential expression of R_e of search models on GIS map

To illustrate the differential expression of R_e obtained by search models on the time and spatial axes, we randomly selected spatial-temporal data of cases from the community of Bicheng and Biquan during the period of May 14, 2016 to July 14, 2016, which includes the spatial location and R_e value of each case. Using a period of two weeks as a unit, we plotted the kernel density map of R_e. The specific results are as follows (Fig. 9):

Fig. 9

The spatial distribution and regional variances of R_e values as acquired through search model

As time passes, the distribution range and color depth of the regional R_e values shift, indicating an increasing, decreasing and then increasing trend in disease transmission during the period. The color depth represents the magnitude of R_e values, and the search model is capable of demonstrating the variability in the spread of the epidemic during the same time period within the region. The black cross symbol denotes the distribution of cases, and it can be seen that the number of cases in a given area does not necessarily correspond to the magnitude of R_e values.

As such, it is not appropriate to solely use the mean of R_e values to represent the development of regional epidemics. The reality is that there are variances in the prevalence of diseases in different communities, neighborhoods, and even smaller areas. Similarly, the number of local cases cannot represent the trend of the development and spread of the disease.

Discussion

We base our evaluation of the severity of infectious diseases mainly on R₀ and case fatality rate [29, 30] However, in practice, we have tools such as vaccination and quarantine to intervene in the progression of the disease. As a result, R_e is a more accurate representation of the actual situation than R0. Additionally, the search model’s setting of isolating first-generation cases at home after diagnosis and rendering them non-infectious is also in alignment with this.

In epidemiological investigations of infectious diseases, R₀ and R_e are commonly acquired using various methods such as contact tracing, infectious disease dynamic model calculations, and estimation of exponential growth rates. Contact tracing is based on real case tracking data, but is only applicable to the early stages of infectious disease transmission. As the transmission chain expands, it becomes difficult to establish a comprehensive data collection system to complete this work. For example, during the COVID-19 pandemic, we can observe China collecting individuals’ health information and travel trajectories through mobile applications to screen and monitor individuals who may potentially be infected with the coronavirus. This approach requires the mobilization of significant social resources and is typically suitable for nationwide or large-scale infectious disease outbreaks. Additionally, the different socioeconomic status, demographics, healthcare resources, and lifestyles across regions lead to variations in the intensity and characteristics of disease transmission, resulting in regional differences in Re [31]. The Re obtained through traditional mathematical modeling mostly pertains to the overall situation of larger regions (provinces, cities, counties) [32,33,34]. However, it is equally essential to understand the transmission dynamics of smaller geographic levels, especially at the local level (communities, streets), to formulate targeted prevention and control strategies.

From the perspective of research methodology, the fundamental premise of this study is rooted in contact tracing, using HFMD as a case study. Combining the fact that HFMD outbreaks mostly occur in childcare facilities or communities, susceptible children are mainly active within the community or within the community-to-childcare facility range during the latent and infectious periods. It is assumed that confirmed cases are searched within a reasonable activity radius around the residence of the primary case, and individuals who meet the conditions for transmission chain formation in terms of time and space are considered as secondary cases. This setting is a scientific assumption based on actual situations and real data. The advantage of this study lies in the stable model, straightforward computation, and easy implementation. It is suitable for assessing and predicting legally notifiable infectious disease outbreaks among scattered children and preschool children, whose activity ranges are relatively fixed. Local public health departments can all achieve this goal. Furthermore, as it is based on scientific assumptions and calculations conducted on the basis of real data for each actual case, the resulting Re value can not only reflect the overall situation, but also can be the Re value of a community or a housing estate, which can reflect the regional differences in the epidemic at the same time period and make differential judgments on the epidemic trend of different communities.

From the perspective of model construction, this study extracts individual time and space factors from a large amount of actual case data from the disease control center. By utilizing GPS coordinates and computer programming calculations, the model is constructed based on two key factors: time nodes and spatial radius. The selection of time nodes takes into account the average incubation period and infectious period of HFMD, with the transmission interval between primary and secondary cases set at 2 weeks. The spatial radius is determined by setting the model’s optimal critical value as the search radius that simultaneously meets two conditions: reaching the maximum Re as the search radius continues to increase and ensuring that the search ranges of the first-generation cases do not overlap. Through comparison of the computational results, a spatial radius of 1.5 km was ultimately chosen as the optimal value, aligning with the characteristics of the daily activity range of local scattered children and preschool children. The Re values obtained from the 2015 to 2019 search models are all greater than 1, indicating a continuous transmission and epidemic trend of HFMD, which is also consistent with the actual situation. By plotting the Re kernel density maps, the regional differences in Re obtained by the search model can be reflected. Disease control and prevention departments can obtain the specific situation of the epidemic in different neighborhoods and different blocks, which facilitates fine-tuning epidemic control and prediction. In addition, we can see that the magnitude and trend of R_e obtained with the spatio-temporal search mode and R_e obtained with the classical SEIR model are in good agreement except for 2018. Similarly, in the similarity analysis, apart from 2018, the similarity coefficient r for the remaining four years is close to 1 (p < 0.05). Given that there is a strong correlation and a relatively close relationship between them, this provides further verification of the scientific validity of the model. The discrepancy in 2018 may be related to the following reasons: the number of HFMD cases throughout 2018 showed a significant increase compared to 2015 (1557), 2016 (2701), and 2017 (1395), reaching a total of 3507 cases for the year. The distribution throughout the year was uneven, with only around 20 cases per month in January to March, while peak periods in June-July and November reached 500–650 cases per month. This uneven distribution may be the fundamental reason for the large difference in obtaining Re values, especially during January to March in 2018. Therefore, we can speculate that the mechanisms for obtaining Re in the two models lead to differences when the case distribution is extremely uneven. However, this highlights the necessity of exploring the search model, as it can reflect the differences in temporal and spatial distribution as well as local and overall differences. Its limitations are as follows: like most models of infectious disease, it cannot distinguish between invisible infection or transmission caused by a carrier; Second, the data comes from data reported by the Centers for Disease Control and Prevention. Underreporting or unnormalized treatment cases can lead to increased errors in the model calculation results.