Delay calibration method and device in spatial metapopulation model of infectious disease

Abstract

The application provides a delay calibration method and device in an infectious disease space collection population model, and belongs to the technical field of infectious disease transmission modeling and risk assessment; the application adopts a delay time discretization method to convert continuous distribution into a discrete probability vector; an infected person allocation algorithm based on matrix operation and a time axis translation algorithm based on matrix operation are designed; the efficiency is significantly improved, and the time complexity is reduced from to; experimental verification shows that the processing speed of a 100-million-population city is at least improved by 7 times, and the processing speed of a single 100-million-population city can be improved by more than 10 times; the whole algorithm is based on matrix operation under a determined dimension variable, so it is easier to use GPU to accelerate operation, and then the algorithm can be extended to a real-time epidemic monitoring and early warning system in the future.

Description

Technical Field

[0001] This invention belongs to the field of infectious disease transmission modeling and risk assessment technology, specifically involving an efficient calibration algorithm for case reporting delay in an infectious disease spatial ensemble population model through probability distribution discretization and matrix operations. Background Technology

[0002] Current infectious disease emergency early warning systems typically include steps such as case reporting, laboratory confirmation, and epidemiological investigation. Because each step inevitably involves time consumption, the system generally suffers from delays. These delays can be caused by multiple factors, such as an inadequate system construction leading to complex and lengthy reporting procedures; a shortage of laboratory reagents limiting the efficiency of pathogen detection; and outdated models and methods resulting in distorted epidemiological investigation results.

[0003] Infectious disease dynamics models are an important method for reconstructing the actual spread of diseases. These models abstractly describe the spatiotemporal processes of pathogen transmission using mathematical tools. By reconstructing the interaction mechanisms among biological, social, behavioral, and environmental factors related to transmission, they provide a reasoning and computational method for studying how these factors influence each other. This allows for practical applications such as infectious disease transmission prediction, hazard effect assessment, intervention optimization, and counterfactual reasoning. Spatial ensemble population models of infectious diseases are a type of infectious disease dynamics model capable of studying the mechanisms of cross-regional disease transmission. These models offer a good trade-off between realistic representation and computational efficiency, playing a crucial role in the early stages of disease dynamics when it was primarily driven by large-scale population movement. However, the lag in current infectious disease emergency warning systems leads to significant delays in case reporting data, affecting the accuracy of parameter estimation in infectious disease dynamics models. Therefore, methods for calibrating infectious disease case reporting delays are crucial for accurately estimating epidemiological parameters and, consequently, accurately assessing the current spread of infectious diseases.

[0004] In existing research on infectious disease transmission modeling and parameter estimation, the delayed calibration of infectious disease case reporting data is one of the key factors affecting the accuracy of inversion of the true transmission trend of infectious diseases. Our research indicates that existing methods address this issue from both the infectious disease model design and epidemiological parameter estimation stages:

[0005] 1. Infectious Disease Model Design: When building an infectious disease model, the lag of observed data should be considered in the model design as variables or parameters. For example, using time-delay differential difference equations to model the dynamics of infectious disease transmission can more accurately describe various time delay phenomena in the transmission of infectious diseases (such as incubation period, delayed immunity, etc.). However, this method has drawbacks such as high solution difficulty and strong parameter sensitivity. Alternatively, the medical treatment process can be incorporated into the model, and model parameters can be estimated using more stable and low-lag observed data such as hospital bed occupancy rate, intensive care unit (ICU) occupancy rate, and number of deaths. However, this method requires medical treatment data, which increases the modeling cost, and the accuracy of medical treatment data also needs secondary verification. Or, the "reporting rate" can be directly introduced as a parameter to be estimated, which intuitively describes the delay between the model's predicted data and the observed data. However, the increase in the number of parameters to be estimated may affect the accuracy and robustness of parameter estimation.

[0006] 2. Epidemiological Parameter Estimation: Based on the Bayesian parameter estimation framework, the estimated parameters in the infectious disease model are updated by calibrating the time when the predicted data is observed and matching it with the actual observed data. For example, the daily number of newly recorded infections is first predicted using the infectious disease dynamics model, and then a corresponding reporting delay time is assigned to each infected person who should be observed, thereby calibrating the actual number of infectious disease cases observed each day. This type of method has the advantages of simplicity, intuitiveness, efficiency, and easy extensibility, and can be applied to most infectious disease dynamics models without increasing model complexity. Its disadvantage is high time complexity.

[0007] While both of the aforementioned methods have their advantages, the accuracy and efficiency of these methods will face even greater challenges for infectious disease models, such as spatial ensemble population models, which involve multiple regions and have higher complexity. This invention can be categorized as a delayed calibration method for infectious disease case report data (Category II). It is specifically designed for the characteristics of spatial ensemble population models of infectious diseases, exhibiting lower time complexity and significantly improving parameter estimation efficiency while maintaining the accuracy of parameter estimation.

[0008] Traditional methods require iterating through each region, each day, and each newly recorded infected individual, sampling and calibrating the observation time of each patient based on a probability distribution function. This iterative process has significant time complexity, severely impacting computational efficiency. In particular, the algorithm's iteration process is population-scale dependent, making its shortcomings even more pronounced when dealing with large populations, thus hindering its ability to handle large-scale, multi-regional, and long-term infectious disease modeling tasks. Summary of the Invention

[0009] In view of this, the present invention designs a reporting delay calibration algorithm with lower time complexity based on the estimation idea of ​​the frequentist school, and provides a delay calibration method and device in the spatial set population model of infectious diseases. It can efficiently update the number of daily reported infected persons through matrix operations while ensuring the accuracy of parameter estimation, thereby improving the algorithm's running efficiency.

[0010] A delay calibration method in an infectious disease spatial ensemble population model includes: First, fitting the report delay time probability distribution to obtain the probability distribution function. ; Indicates a time delay, in days;

[0011] The second step is to fit the probability distribution function. Discretization yields the discrete distribution law. ;

[0012] The third step is to construct the probability distribution vector:

[0013] The probability distribution vector is obtained by arranging the probabilities in the discrete probability law in sequence. This refers to the time delay between infection and diagnosis for each infected person. For 0 days, 1 day, , The probability of the day, then the probability of the day. element for:

[0014] ;

[0015] in, This represents the time interval from infection to diagnosis for each recorded infected individual; T represents the observation time, in days.

[0016] Step 4: Calculate the distribution matrix of infected individuals:

[0017] According to the infectious disease dynamics model, let the number of newly recorded infected persons in region i on day t within the observation period be denoted as . The set of its reporting delay times is equivalent to a random variable. Extracted from One sample; according to the ideas of the frequentist school, Defined as:

[0018]

[0019] in, ; Indicates the total number of regions; This means that the median of the sample is equal to The number of recorded infections in region i The number of people at each location is determined by the probability distribution vector. The estimated result is:

[0020]

[0021] The first The number of newly recorded infections in all regions of the day is aggregated into a vector, denoted as . ; A suitable delayed reporting time is assigned to each newly added recorded infected person through matrix operations, and the reporting time for the first... Daily distribution matrix of infected individuals:

[0022]

[0023] Right now Each element This indicates that among the newly recorded infections in region i on day t, the reporting delay period is [number] days. The number of people;

[0024] Step 5: Update the infected person count matrix:

[0025] Construct a translation matrix with T rows and T columns. for:

[0026]

[0027] Using translation matrix Will Shifting the elements in the matrix to the right by t units yields the adjusted distribution matrix as follows:

[0028]

[0029] in, Representation matrix The power of t, that is, for Right multiply by t matrices To achieve The effect of shifting the element to the right by column t;

[0030] At this point, the number of newly recorded infected individuals on day t is updated using matrix operations:

[0031] ;

[0032] Among them, matrix This matrix represents the actual number of infected individuals that can be observed.

[0033] Traversing observation time Complete the task of calibrating the delayed reporting of infectious disease cases.

[0034] Ideally, in the first step, we assume that the time interval from infection to diagnosis for each recorded infected person is a continuous random variable. The probability distribution function is denoted as .

[0035] Ideally, the probability distribution function is fitted using the observed data. The parameters determine its expression.

[0036] Ideally, the discrete distribution law obtained in the second step is... for:

[0037] .

[0038] Ideally, in step four, for Round each value to the nearest integer.

[0039] Preferably, the matrix Each element in the equation is rounded down.

[0040] Ideally, a sparse matrix storage method should be used for storing the matrix.

[0041] A delay calibration device includes a first module, a second module, a third module, a fourth module, and a fifth module;

[0042] The first module is used to perform the first step, fitting the probability distribution of the report delay time to obtain the probability distribution function. ; Indicates a time delay;

[0043] The second module is used to perform the second step, fitting the probability distribution function. Discretization yields the discrete distribution law. ;

[0044] The third module is used to perform the third step, constructing a probability distribution vector;

[0045] The fourth module is used to perform the fourth step, calculating the distribution matrix of infected individuals;

[0046] The fifth module is used to perform the fifth step, updating the number of infected persons matrix, thereby completing the task of calibrating the delayed reporting of infectious disease cases.

[0047] The present invention has the following beneficial effects:

[0048] This invention provides a delay calibration method and apparatus for a spatial ensemble population model of infectious diseases. It employs a delay time discretization method to convert a continuous distribution into a discrete probability vector; it designs an infected person allocation algorithm based on matrix operations and a time axis translation algorithm based on matrix operations; significantly improving efficiency, with a time complexity of approximately [missing value]. Down to Where N is the city population, L is the number of cities, and T is the observation duration. Experiments show that the processing speed for 100 cities with a population of over one million is improved by at least 7 times, and the processing speed for a single city with a population of over one million can be improved by more than 10 times. The entire algorithm is based on matrix operations under variables of a fixed dimension, making it easier to utilize GPUs for accelerated computation, and thus it can be extended to real-time epidemic monitoring and early warning systems in the future. Attached Figure Description

[0049] Figure 1 This is a flowchart of the present invention;

[0050] Figure 2 This represents the time required for a single iteration in the parameter estimation algorithm compared to the traditional traversal calibration method.

[0051] Figure 3 The results are parameter estimations for real-world cases;

[0052] Figure 4 Basic reproduction number of real cases The estimation results. Detailed Implementation

[0053] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0054] Current infectious disease emergency early warning systems typically report cases in integer time units of "days." To improve computational efficiency, infectious disease dynamics models are usually discretized and iteratively solved using integer step sizes of 1. To match the infectious disease dynamics model, we assume that the delay time for case reporting also occurs in integer multiples of "days." In the infectious disease spatial ensemble population model, let E represent the population in the exposure state. This represents the recorded infected individuals that will be observed. Therefore, the number of newly recorded infected individuals in the i-th region on day t can be represented as those from the state. Transition to state The population, the number of this population is , , It is the total number of regions. , This refers to the number of observation days. Due to the lag in infectious disease emergency early warning systems, the actual number of infected individuals that can be observed in the i-th region on day t is not equal to... Therefore, the main objective of this invention is to calibrate the time when these infected individuals were observed, and to calculate the actual number of recorded infections that can be observed each day and in each region, so as to match the actual observation data for parameter estimation. Let the matrix... This represents a matrix indicating the actual number of infected individuals that can be observed, where each element... This represents the number of infected individuals that can actually be observed in the i-th region on day t. The traditional method is... Each infected person is assigned a reporting delay time. and accumulate them into the matrix. middle:

[0055]

[0056] Within this framework, all regions L, all observation days T, and all actually observable infected individuals are processed cyclically. This allows us to obtain the actual number of infections that can be observed daily in each region. However, this method involves at least three nested loops, resulting in high time complexity. This invention transforms the above-mentioned traversal calibration process into a matrix operation with a single loop, significantly improving efficiency while maintaining effectiveness. See the detailed process below. Figure 1 The specific steps are as follows.

[0057] The first step is to fit the probability distribution of the reporting delay time. It is assumed that the time interval from infection to diagnosis for each recorded infected person is a random variable. According to existing research, the reporting delay time in the early stages of COVID-19 transmission can be well fitted by the gamma distribution. Since the gamma distribution is a continuous random variable distribution model, this invention also assumes... It is a continuous random variable, denoted as . Then, the probability distribution function can be fitted using the observed data. The parameters.

[0058] The second step is to design a discrete distribution. Since the smallest time unit in current infectious disease emergency early warning systems and common infectious disease models is an integer "day," the continuous random variable will be... Discretizing integers can speed up the efficiency of subsequent algorithms. Specifically, by calculating random variables... In The integral of the probability density function within a left-closed, right-open interval with radius 0.5, centered at [a_1], is used to approximate the random variable. exist The probability of proximity. To conform to the definition of the distribution law, reasonable constraints need to be imposed on the boundary, ultimately constructing a law that can approximate... The discrete distribution law of it is denoted as :

[0059]

[0060] The third step is to construct the probability distribution vector. This involves arranging the probabilities in the probability distribution law sequentially to obtain the probability distribution vector. This refers to the time delay between infection and diagnosis for each infected person. For 0 days, 1 day, , The probability of the day, then the probability of the day. element for:

[0061] ;

[0062] The fourth step is to calculate the distribution matrix of infected individuals. Based on the infectious disease dynamics model, the number of newly recorded infected individuals in region i on day t is... Their reported delay times constitute a set equivalent to random variables. Extracted from One sample. According to the frequentist school of thought, It can also be defined as:

[0063]

[0064] in, This means that the median of the sample is equal to The number of items, within the research framework of the invention, corresponds to a delayed reporting time of [number]. The number of people. Therefore, at this time, the region Recorded infections in The number of people at a location can be determined using probability distribution vectors. The estimated result is:

[0065]

[0066] The first The number of newly recorded infections in all regions of the day is aggregated into a vector, denoted as . Therefore, matrix operations can be used to assign an appropriate delayed reporting time to each newly recorded infected person, and to calculate the... Daily distribution matrix of infected individuals:

[0067]

[0068] Right now Each element This indicates that among the newly recorded infections in region i on day t, the reporting delay period is [number] days. The number of people. Note that since the number of people should be an integer, the approximate equal sign is used here, and the integer part will be rounded off in the final output.

[0069] Step 5: Update the infected person count matrix. Assume the current time is t, and in order to update the infected person distribution matrix at this time... Accumulated to In the middle, it is also necessary to The elements in the matrix are shifted to the right by t units. Construct a translation matrix with T rows and T columns. for:

[0070]

[0071] Clearly, any matrix right multiplication The power of t is equivalent to shifting all elements of the matrix to the right by t units, filling missing elements with 0s, and discarding redundant elements. The resulting adjusted distribution matrix is ​​as follows:

[0072]

[0073] in, Representation matrix The power of t, that is, for Right multiply by t matrices To achieve The effect of shifting the element to the right by column t;

[0074] At this point, the number of newly recorded infected individuals on day t is updated using matrix operations:

[0075] ;

[0076] Among them, matrix This matrix represents the actual number of infected individuals that can be observed.

[0077] Traversing observation time Complete the task of calibrating the delayed reporting of infectious disease cases.

[0078] Thus, this invention utilizes matrix calculations to replace element-wise iteration, significantly improving computational efficiency, especially when dealing with large urban populations. When the size is large, the efficiency improvement brought by the method of this invention is more significant. It is worth noting that in our experiments, we used the computing engine "odin" to accelerate the simulation of the disease transmission process (FitzJohn RG, Knock ES, Whittles LK, et al. Reproducible parallel inference and simulation of stochastic state space models using odin, dust, and mcstate[J]. Wellcome OpenResearch, 2021, 5: 288.), and this engine currently only supports the calculation of matrices of a certain dimension. Therefore, our algorithm can be used as part of the infectious disease dynamics model as input into the algorithm package "odin" for accelerated calculation. Figure 2The time required for a single iteration of the method of this invention and the traditional method are plotted for different total number of cities L and different population sizes N in each city. It can be seen that replacing iterative loops with matrix operations significantly improves the efficiency of parameter estimation.

[0079] This invention is an improvement upon the literature "Li R, Pei S, Chen B, et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)[J]. Science, 2020, 368(6490): 489-493.", constructing a suitable infectious disease transmission modeling and parameter estimation framework to test the application effect of this invention. The improvements mainly include the following two parts:

[0080] 1. Use binomial and multinomial distributions to describe the cross-regional movement of people in different epidemic states and the transfer process between different cells. Compared with the traditional Poisson distribution, the binomial and multinomial distributions are more suitable for describing the transfer process between cells of a limited population and for modeling the early cross-regional spread of disease. Blindly using the Poisson distribution can easily lead to data overshoot, thus overestimating the number of people in the mobile population who are capable of transmitting the disease, especially in the early stages of the disease when the number of infected people is small, which may lead to large errors.

[0081] 2. Parameter estimation is performed using the particle Markov chain Monte Carlo (pMCMC) algorithm. Compared to the traditional ensemble Kalman filter (EnKF) algorithm, pMCMC is more suitable for parameter estimation problems in nonlinear systems such as infectious disease transmission.

[0082] Table 1 Initialization settings for simulation experiment parameters

[0083]

[0084] Based on the aforementioned infectious disease transmission modeling and parameter estimation framework, we constructed a city network using 375 cities in China as nodes and conducted experimental verification based on real observational data. The Li R literature provides population movement data and the number of newly reported COVID-19 cases per day for 14 days, from January 10th to January 23rd, 2020. As of January 23rd, 2020, the national COVID-19 monitoring report showed a cumulative total of 801 cases nationwide, with 454 cases reported in Wuhan. The correctness and reliability of this data have been verified in the original paper. Therefore, it is suitable as a control method for this invention. Taking Wuhan as the starting point of disease transmission, Table 1 summarizes the specific parameter ranges and priors. The final parameter estimation results are as follows: Figure 3 As shown.

[0085] The results show that the parameters obtained by both methods generally exhibit a unimodal distribution, and the peak positions of most parameters are very close. Since the effectiveness of the traditional method in this case has been verified in LiR literature, this indirectly validates the effectiveness of the method presented in this paper. It should be noted that the two methods have different effects on... The estimation results differ somewhat; both models exhibit two local extrema, one close to 0 and the other close to 0.5. The actual percentage of infected individuals is recorded. It cannot be 0; it is more likely to be 0.5. The result obtained by the corresponding method of this invention is more reliable for this parameter.

[0086] Due to the basic reproduction number It is one of the most important epidemiological parameters. Based on the parameter estimation results, we calculate the basic reproduction number under the simulation scenario. :

[0087]

[0088] Figure 4 Showing The distribution of. From The distribution also shows that the results obtained by the model corresponding to this invention exhibit a unimodal distribution, while the competing model has two local peaks. Since the values ​​of the two main peaks are very close, this indirectly verifies the accuracy of our method. However, based on real-world conditions... Typically, only one definite value can exist, therefore our method is more reliable.

[0089] Finally, we summarized all the parameters and The median and 95% confidence interval (Table 2) of the two methods, although there are some differences, are very close. This indirectly confirms that our method has made significant improvements in algorithm efficiency and robustness while ensuring accuracy and effectiveness.

[0090] Table 2 Comparison of parameter estimation results of the present invention and traditional methods in real-world cases.

[0091]

[0092] Based on the above-described delay calibration method, this invention also provides a delay calibration device, comprising a first module, a second module, a third module, a fourth module, and a fifth module.

[0093] The first module is used to fit the probability distribution of the report delay time, obtaining the probability distribution function. ; Indicates a time delay;

[0094] The second module is used to fit the probability distribution function. Discretization yields the discrete distribution law. ;

[0095] The third module is used to construct the probability distribution vector;

[0096] The fourth module is used to calculate the distribution matrix of infected individuals;

[0097] The fifth module is used to update the number of infected persons matrix, thereby completing the task of calibrating the delayed reporting of infectious disease cases.

[0098] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A delayed calibration method in a spatial ensemble population model of infectious diseases, characterized in that, include: The first step is to fit the probability distribution of the report delay time to obtain the probability distribution function. ; Indicates a time delay, in days; The second step is to fit the probability distribution function. Discretization yields the discrete distribution law. ; The third step is to construct the probability distribution vector: The probability distribution vector is obtained by arranging the probabilities in the discrete probability law in sequence. This refers to the time delay between infection and diagnosis for each infected person. For 0 days, 1 day, , The probability of the day, then the probability of the day. element for: ； in, This represents the time interval from infection to diagnosis for each recorded infected person; T represents the observation time, in days. Step 4: Calculate the distribution matrix of infected individuals: According to the infectious disease dynamics model, let the number of newly recorded infected persons in region i on day t within the observation period be denoted as . The set of its reporting delay times is equivalent to a random variable. Extracted from One sample; according to the ideas of the frequentist school, Defined as: ； in, ; Indicates the total number of regions; This means that the median of the sample is equal to The number of recorded infections in region i The number of people at each location is determined by the probability distribution vector. The estimated result is: ； The first The number of newly recorded infections in all regions of the day is aggregated into a vector, denoted as . ; A delayed reporting time is assigned to each newly added recorded infected person through matrix operations, and the reporting time for the [number]th [case] is calculated. Daily distribution matrix of infected individuals: ； Right now Each element This indicates that among the newly recorded infections in region i on day t, the reporting delay period is [number] days. The number of people; Step 5: Update the infected person count matrix: Construct a translation matrix with T rows and T columns. for: ； Using translation matrix Will Shifting the elements in the matrix to the right by t units yields the adjusted distribution matrix as follows: ； in, Representation matrix The power of t, that is, for Right multiply by t matrices To achieve The effect of shifting the element to the right by column t; At this point, the number of newly recorded infected individuals on day t is updated using matrix operations: ; Among them, matrix This matrix represents the actual number of infected individuals that can be observed. Traversing observation time Complete the task of calibrating the delayed reporting of infectious disease cases.

2. The delayed calibration method in a spatial ensemble population model of an infectious disease as described in claim 1, characterized in that, In the first step, it is assumed that the time interval from infection to diagnosis for each recorded infected person is a continuous random variable. The probability distribution function is denoted as .

3. The delay calibration method in a spatial ensemble population model of infectious diseases as described in claim 2, characterized in that, Fitting a probability distribution function using observation data The parameters determine its expression.

4. The delayed calibration method in a spatial ensemble population model of an infectious disease as described in claim 3, characterized in that, In the second step, the discrete distribution law is obtained. for: 。 5. The delayed calibration method in a spatial ensemble population model of an infectious disease as described in claim 4, characterized in that, In the fourth step, for Round each value to the nearest integer.

6. The delay calibration method in a spatial ensemble population model of an infectious disease as described in claim 5, characterized in that, matrix Each element in the equation is rounded down.

7. The delayed calibration method in a spatial ensemble population model of an infectious disease as described in claim 6, characterized in that, For matrix storage, a sparse matrix storage method is adopted.

8. An apparatus for implementing the delay calibration method according to any one of claims 1-7, characterized in that, It includes Module 1, Module 2, Module 3, Module 4, and Module 5; The first module is used to perform the first step, fitting the probability distribution of the report delay time to obtain the probability distribution function. ; Indicates a time delay; The second module is used to perform the second step, fitting the probability distribution function. Discretization yields the discrete distribution law. ; The third module is used to perform the third step, constructing a probability distribution vector; The fourth module is used to perform the fourth step, calculating the distribution matrix of infected individuals; The fifth module is used to perform the fifth step, updating the number of infected persons matrix, thereby completing the task of calibrating the delayed reporting of infectious disease cases.

Images