Machine learning spatio-temporal epidemiological model to evaluate Germany-county-level COVID-19 risk

Public health safety is the cornerstone of economic and social stability. The CoronaVirus Disease 2019 (COVID-19) is an infectious disease caused by the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) which became a global pandemic that spread across the world in 2020. As of early March 2020, the number of reported cases in Germany experienced a rapid growth. Only severe restrictions on contacts and travel can ease the situation in the summer. However, the number of infected cases have rapidly increased again since mid-September, named popularly as the second-wave pandemic worldwide. By January 25, 2021, the pandemic has caused 2141 665 confirmed cases and 52 087 fatalities in Germany ⁶ . Many states are facing the dilemma challenges of public heath and economics in curbing the spread of the COVID-19. The medical resources of Germany (ICU beds, health workers and doctors) have not arrived at their ceiling, but the current phase of the pandemic is not encouraging. As a highly-infectious disease, COVID-19 is expected to remain spreading and be even worse to cause a higher number of infections and fatalities in the future. No well adapted medical treatments nor a sufficient vaccine allocation is available at this moment, so the non-pharmaceutical public health policies, such as, social distancing and travel restrictions can be implemented to alleviate the spread of epidemic currently. This policy has proven an effective containment in East Asia, e.g. China [1], South Korea and Japan. In central Europe this strategy worked as well at early stage [2]. In the presently worsening health crisis, the need to establish a multi-scale (nation-state-county) COVID-19 risk prediction model emerges, which can guide the evaluation and improvement on the public policies of governments (states) at the different regional levels. The multi-scale risk information can help to improve the non-pharmaceutical intervention strategies such as cross-state travel restriction, home isolating regulations and allocation of vaccination resources. The public can also benefit from the information, to make effective access of medical resources and help implementing or following the public policies effectively.

Evaluation of the risk of the on-going COVID-19 pandemic across the diverse geographical administrative units is challenging, because of the substantial heterogeneity rooted in the federal state political system with diverse economic composition, traffic conditions and individual political views. Our research presents the needed tool which can optimize our understanding of the importance of the local and regional (state-county-level) pandemic dynamics. Studies devoted to the dynamics of the epidemics based on mathematical models [3, 4] (also see reviews [5–7]) have been mainly focusing on time-domain analysis and ignored the influence of geographic distribution and the population migration [2], which in fact presents cross-county transfers on a daily-scale. A model which can predict the details of the spatio-temporal evolution of the COVID-19 cases is urgently needed to appropriately assess the local risk in a county of politically diverse multi-scale regions. Works have been concerned with spatial distribution or local evolutions were mainly based on historical statistical data [8–10]. Several recent studies share similar goals. Different countries [11–14] merited this as an urgent need across the world. From a practical perspective, an intellectual and efficient risk prediction model needs to be performed on a county basis, like with all Landkreise ⁷ in Germany, since people are highly mobile and connective within Germany and to all of Europe as well as all other continents.

Machine learning (ML), a branch of artificial intelligence (AI), efficiently integrates statistical and inference algorithms and thus offers the opportunity to uncover hidden structures of evolution in complex data and to describe it with finite dynamical parameters. Therefore, a combination of ML algorithms and spatio-temporal epidemiological models for risk prediction in the context of COVID-19 pandemic is a promising option. Although the state-of-the-art ML techniques have been deployed into applications of COVID-19 pandemic predictions [15–18], the highest concern for the governments and residents is how fast the COVID-19 spreads. Available ML models lack of capability of predicting the needed multi-scale evolution, since they routinely focus on this exhibit promising implications for individuals and on a national scale, and thus are impeded by the paucity of validation information and limiting privacy rules in Germany.

This study aims to establish such a new evaluation paradigm by introducing a spatio-temporal ML assisted epidemiological model. This dynamical model consists of cellular automaton (CA) [19, 20] plus a modified SUIR model [13, 21], where the dynamical parameters are extracted from the reported data of Robert Koch-Institut (RKI) by the Conv2LSTM Neural Network, a deep learning technique regularly applied in video analysis. The present SUIR model is more inclusive than previously used susceptible-infectious-recovery (SIR) models, as it considers explicitly the Undetected or Asymptomatic populations [22, 23], hidden due to the lack of testing and also the systematic uncertainties for the testing. We reproduces the evolution of infection maps from 25 February to 28 March 2020, and extracts the infectious parameters from RKI data in September 2020 with a neural network which was training on the data-set from CA-SUIR model simulations. The well-learned machine predicts county-level infection maps from September to November 2020, which is matched with the reported situation. With regard to the high-profile holiday, the different public health policies are evaluated in the end of the section, in which the entire, partial and none lock-down strategies are compared.

The full flow chart is assembled in figure 1. to present the multi-scale spatio-temporal simulation of epidemics, where as the first step real maps are projected onto the topological equivalent maps as explained in figure 8. The cumulative infected cases in the real profiles ⁸ are reassembled into a topological equivalent map with 25 × 25 sites in Germany, which contains 412 counties (see details in figure 8). The micro-states of each county, its number of infected cases, are driven by the SUIR and CA simultaneously in the model, which construct the infection map. With simulated data from the CA-SUIR model, the Conv2LSTM neural networks [24] are designed to learn to extract dynamical parameters of CA-SUIR model. Subsequently, we transfer the well-trained neural networks to make timely prediction and risk evaluation on RKI infection maps.

Zoom In Zoom Out Reset image size

2.1. Extracting dynamical parameters via deep learning

To prepare training data for the Conv2LSTM, we generate 30-day simulation maps as short movies by associating the SUIR model parameters in a 10% uncertainty range. The total data size is 10 000, in which 90% is treated as the training data-set and the rest is the testing data-set. The 30 day movie has 30 frames which is fed to the machine with continuous 7 day segments as slip windows [25]. The structure of the Conv2LSTM neural networks are exhibited in figure 10. The Conv2LSTM neural networks we built are demonstrated in figure 2, in which the time distributed convolutional neural networks (CNNs) layers and long short-term memory (LSTM) layers are combined to decode the information hidden in the short video. The details of the structure are arranged in figure 10. The time distributed CNNs could extract features (with temporal structure preserved) from sequential image inputs and then embed them into the following LSTM layers, which is able to efficiently drill evolution rules from the short video.

Zoom In Zoom Out Reset image size

2.2. Spatio-temporal CA-SUIR model

We construct a cellular automaton model for tracking the epidemics spatial evolution in a closed system (nation with fixed boundary as figure 1 shown). The underlying structure is a L × L cell grid, where L is the system size, and the number is from the topological map. The color of cell characterizes the population in minimum geographical unit, which is county in our case. It is an instructive example to set the cell size to a typical area of city, such as 10²–$10^3~\textrm{km}^2$, which can simulate the epidemics in the corresponding resolution approximately. The resolution could be naturally increased into community level, but it needs more detailed data which is laboriously accessible, thus the following discussion focuses on county level and corresponding resolution. This model adopts the Moore neighbor, and cells update their states by transition tensor $P_{m,n}^i(t)$ with 4 ×L × L shape, which means the possibility of transfer from city (m, n) to direction i at t time-step. The Pⁱ (t) are 4 matrices $(P^\textrm{u}(t),P^\textrm{d}(t),P^\textrm{l}(t),P^\textrm{r}(t))$ with same lattice size, in which different colors labels the direction of transfer. In the following simulation, we adopt the mean field approximation [26], which leads to the uniform transition possibility $P^\textrm{u}(t) = P^\textrm{d}(t) = P^\textrm{l}(t) = P^\textrm{r}(t) = \epsilon$. The transfer rate is introduced to describe the average movement of the residents on each site. To mimic the border control induced closed boundary, a transfer rate  = 0 has been used in the primary epidemic stage [27]. The evolution of the CA-SUIR model is as follows: Initialization. Set the position of infection burst (x, y) and generate population distribution S(t) at the L × L lattice. At t = 0 time-step, epidemic turns on and population begins to migrate; SUIR dynamics. At t time-step, evolve SUIR model $t^{^{\prime}}$ times on each site; movement. The parts of Undiagnosed people on site (m, n) migrate into the neighbor site with transition matrices P^k (t, m, n) at t + 1 time-step; update all sites synchronously until T time-step is reached.

The temporal evolution parts are driven by the SUIR model and are demonstrated in figure 3, where the infectious rate λ controls the rate of spread which represents the probability of transmitting disease between a susceptible and an infectious individual. The diagnosed rate σ dictate the probability of latent infectious individuals becoming confirmed by testing. Removed rate γ contains the rate of mortality and recovery from the infectious. Considering the special features of COVID-19, we introduce a considerable population that could be infected but without any symptoms, or, they just have no access to perform a detection(with possibility p), which is concluded in a new compartment, undiagnosed (U). In contrast to the infected who is usually subjected to be isolated in medical facilities or quarantine at home, the undiagnosed population has normal mobility which leads to the growing-up infection; the infectivity could decrease after d days [28], eventually removed from the infection population, which means immunity will retain at least several months [29]. In a closed population with no birth or mortality per day, the modified SUIR model becomes,

$$\begin{align} \frac{d S}{d t} & = -\frac{S}{N}\lambda U \end{align} \tag{ 1 }$$

$$\begin{align} \frac{d U}{d t} & = \frac{S}{N}\lambda U -\sigma p U- \frac{\sigma(1-p)}{d}U, \end{align} \tag{ 2 }$$

$$\begin{align} \frac{d I}{d t} & = \sigma p U-\gamma I \end{align} \tag{ 3 }$$

$$\begin{align} \frac{d R}{d t} & = \gamma I + \frac{\sigma(1-p)}{d}U, \end{align} \tag{ 4 }$$

where $N = S + U + I + R$ is the total population and we adopt a redefinition as $\sigma \equiv \sigma p + \frac{\sigma(1-p)}{d}$ in the results.

Zoom In Zoom Out Reset image size

It should be noted that the choice of the spatial or temporal models are not unique. The modified SIR (or SEIR) models are routinely adopt in forecasting COVID-19 temporal evolution. Although it is not able to properly forecast the long-term situation [30], it offers a reasonable base-line for short-term prediction which could be used to train the deep learning machine sufficiently. On the other side, the temporal model could be improved by introducing more realistic factors, such as the risk of reinfection [31], the isolation measures for patients with mild symptoms [32] or the close contacts [33]. With regard to the spatial part, the philosophy is similar. The accurate tracking of the individuals [34] or agent-based models [8] could improve the performance of the prediction on dynamical transmission of COVID-19 in Germany.

3.1. Pandemic dynamics modeling, data-set generation and network capacity

With the daily number of cumulative cases in Germany (Data is from the open database of RKI ⁹ ), we determine the temporal parameters of the CA-SUIR model to be infectious rate λ = 0.293, diagnosed rate σ = 0.129, and removed rate γ = 0.628 by matching well with the first-wave burst out in March especially the rapidly increasing stage (from 25 February 2020 to 28 March 2020), as presented in figure 4. The initial undiagnosed cases remain same spatial distribution as the reported infection cases with U₀≃ 700, meanwhile the transfer rate connected counties is set as  = 0.08. Figure 4 suggests that there is an overestimation of the simulation at early stage (first 20 days), which could be naturally understood to be related with the underestimation in reality at that time since the limited detection ability and the supply of test reagents were insufficient [35]. The estimation to the uncertainty of undiagnosed population shown in figure 4, which is actually from the stability issues of RT-PCR testing [36] and unreported cases in practical situation.

Zoom In Zoom Out Reset image size

We train the Conv2LSTM to predict the dynamical parameters of the involved epidemiological model (see the learning curve in figure 11). The four sub-figures in figure 5 demonstrate the comparison between the prediction from the Conv2LSTM and ground-truth of testing data, they are the transfer rate , the infectious rate λ, the diagnosed rate σ and undiagnosed population U respectively. The well-trained machine extracts parameters from the RKI infection map of March 2020 as $\lambda = 0.26, \sigma = 0.09 , U_0 = 1352, \epsilon = 0.068$. They are numerically matched with the fitting results, in which the diagnosed rate and initial undiagnosed population maintain the consistency with fitting parameters in the relationship σU₀.

Zoom In Zoom Out Reset image size

3.2. Machine learning-assisted risk prediction

To assess the COVID-19 risk in realistic situation timely, we deploy the well-trained neural network into the current cases in transfer learning manner. Transfer learning is a ML method where a model developed for a task is generally reused as the starting point for a model on a second task. In our case, we transfer the pre-trained models into the on-going second-wave COVID-19 epidemics in Germany, extracting the dynamical parameters of the spatio-temporal pandemics from the 30 day RKI infection maps with beginning from September 8, 2020, which yields the dynamical parameters as $\lambda = 0.21, \sigma = 0.167, U_0 = 7426, \epsilon = 0.09$.

Figure 6 describes the homogeneity between simulation and RKI infection maps from September to November 2020 in Germany. The former cumulative infection maps shown in figure 6(a) are generated from the SUIR model with the ML parameters during 3 months. The latter maps shown in figure 6(b) are projected from the RKI data, in which we present the cumulative infection cases in real profiles until the end of November. The top 5 biggest cities, Berlin (13^∘23E, 52^∘31N), Hamburg (10^∘00E, 53^∘33N), München (11^∘35E, 48^∘03N), Köln (6^∘58E, 50^∘56N) and Frankfurt (8^∘41E, 50^∘07N) demonstrate a worse situation than the other sites. The slight mismatching on the map of November is due to the actual soft lock-down policy has been executed from 2 November. The functional regulations will contain the pandemic moderately, which will be examined under different policies in the following part.

Zoom In Zoom Out Reset image size

3.3. Public policies evaluation

In public policy, three different degrees of restriction rules, named as the unlocked, soft-locked and locked strategies are evaluated. The last one is the lock-down policy came into play first in Wuhan [37, 38] which shows strictest restriction to the mobility of residents compared with the other policies. The restrictive rules take effect in different places with differentiated implementation, which includes but are not limited to office, campus, school, restaurant and cross-county travel. How is the effect of implementation? It is of crucial concern for both governments and residents who are enduring the social-economical pressure in pandemic [11, 39–41]. The comprehensive influences at county-level are averaged as an equivalent effect to the transfer rate, which means the effect of policies are coarse graining to the county-level. Reasonably, it could be improved if we collect more precise data below county-level, e.g. individual or community-level [42].

The differences among them in simulation are the transfer rate and susceptible population: the most strict locked rules has transfer rate  = 0, the soft-locked with half transfer  = 0.045, and the rule with no additional restriction has  = 0.09 as default; meanwhile, the containment policies could result in a consequence of infectious population decreasing that effectively deplete the susceptible population [1]. It leads to the reduction of the susceptible population, reaching to 1/2, 1/8 for the locked and soft-locked rules respectively, in which the 1/8 is from the rough estimation to the necessary out time. These three strategies are implemented from the same time point 25 November 2020, and their 30 day simulation results are shown in figure 7. In our evaluation to the policies, 2.3 million people would be infected at Christmas of 2020 in the worst situation, where daily increase is approximately 80 thousand, which means the total mortality will reach 34.5 thousand approximately (estimate with mortality rate ~1.5% [43]). With regard to the soft-locked policy, the infection number reduces to 1.4 million with 21 thousand fatalities and the corresponding daily increasing cases are 20 thousand. In the most strict rules, although it restricts the mobility, the number of infection cases will be 1.07 million and the daily increase is only 4 thousand. Compared with the reported cases shown in figure 7(d), they achieve the consistency in geographic distribution. Considering the reported cases from RKI, cumulative cases 1612 648 with daily increasing 25 533 and 29 182 fatalities, the lock-down policies have clearly made a difference. It can effectively reduce the chance of exposure to the infectious population, which reflects in the reduction of the transfer rate and effective susceptible population in our model. It is also found that the efficiency of the vaccination or allocation of the vaccine could be embedded into our CA-SUIR model through the similar method that reduction to the susceptible population [1, 44], which is mathematically equivalent with introducing a suppression factor to the susceptible population as −νS. The factor ν could include a time dependent effect to show the progress of the vaccination in different counties. As a matter of fact, our framework has the alignment with reality which suggests a succinct scheme for quantifying the effects of the public policies. The quantification of policy effects is conceivably shown in the variation of the dynamical parameters or populations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we construct a multi-scale spatio-temporal epidemiological model, which enables the prediction on the state-county-level COVID-19 risk over 412 Landkreise in Germany. Reasonably, it is convenient to replace the module by other elaborate spatial or temporal models for reaching a more accurate prediction in our framework. Regarding the ML module, transfer learning helps us to extract the dynamical parameters from the RKI infection maps through the well-trained Conv2LSTM neural networks.

The CA-SUIR model reproduces the epidemics in March and September acknowledged as the first- and second-wave of COVID-19 pandemic. Even though the prediction focuses on the infectious dynamics in the current paper, it is conveniently feasible to derive recovery and mortality rate from the model under the guarantee of accurate data sources confronted. The simulations perform properly to describe and predict the data in the initial 30 day phase, while it shows a tendency towards a faster smoothing-out of the pronounced local fine structure and persistent hot spots, as compared to the data. It could be understood that the observed infection cases concentrated in the big cities, such as Berlin, Hamburg, München, Köln and Frankfurt than the other sites in the simulation, which reflects high population density enhances the diffusion of the COVID-19. The prediction could be evidently improved by introducing more abundant interactions across and within the counties. Some practical improvements can be implemented in future works, such as replacing the geographic lattice to traffic networks [5, 42, 45–47], introducing the higher resolution lattices beyond the county-level which also means higher-order interactions among the residents [12, 20, 42, 48]. The well-trained machine generates the infection maps which are consistent with the RKI data in 3 months. With regard to the slight uncertainty in extracting dynamical parameters, it could be avoided by training the neural networks on more heterogeneous data-set. We will present the further improvement in our future papers.

The ML-assisted CA-SUIR model indicates that knowledge of risk evaluation in the county-level can be exhibited to the public directly and timely, meanwhile, the 7 day windows could forecast the current situation robustly. However, we could not ignore the limitation of this study. The main shortage of our framework is highly dependent on the accuracy of the CA-SUIR model which is oversimplified for the complicated COVID-19 situation. Our study was also limited by the absence of accurate mobility data which unavoidably affects the feasibility. The assessment to the different public health policies suggests that the transmission modes of coronavirus could be shaped by the non-pharmaceutical interventions dramatically. The other valuable application of our evaluation system is to improve allocation of medical resources which are usually unequal in different counties, some of them are suffering from the epidemic that far exceeds their medical capacity.

The authors thank Esteban Vargas, Maria Babarossa and Jan Fuhrmann for useful discussions and comments. The work on this research is supported by the BMBF through the ErUM-Data funding and the Samson AG AI Grant (L W, K Z), by the Walter Greiner Gesellschaft zur Foerderung der physikalischen Grundlagenforschung e.V. through the Judah M. Eisenberg Laureatus Chair at Geothe Universitaet Frankfurt am Main (H S), by the National Natural Science Foundation of China, Grant No. 11875002 (T X, Y J), by the Zhuobai and Bajian Programs of Beihang University (Y J). The authors also thank the donation of NVIDIA GPUs from NVIDIA Corporation.

The data that support the findings of this study are openly available at the following URL/DOI: https://covid-19-geohub-deutschland-esridech.hub.arcgis.com/search?collection=Dataset and https://github.com/Anguswlx/Conv2LSTM-COVID19.

In order to simulate the population evolution on lattice, the topological mapping is adopted in our CA computation. We firstly embed the map in a rectangle area whose length and width are chosen as the maximum value of the Germany map, as in figure 8(left). Then the area is segmented into uniform sites of L × L. And each site is labelled by its row and column numbers as B_mn , where $1\leqslant m, n\leqslant L$. If a county occupies N sites, all kinds of the population, such as susceptible, undiagnosed, infected and removed, are assigned equally into each site. And if there is more than one county or county part on a site, the populations of the site is set as the sum of all these counties or county parts. Because the geographic distance is not important in our simulation, which focus on the population distributions, we set each site as a square when plotting the distribution in figure 8(right). In the current simulation, we confine all the population in Germany and thus set the out-of-Germany transfer possibility at the border as zero.

According to the initial infection cases distribution and the locations of main cities of Germany, we set the initial undiagnosed population distribution U₀(m, n) to be proportional to the population distribution of the cities N(m, n). In the simulations for March we choose the transfer rate  = 0.08 as a constant, the population of transfer on site is proportional to the rate in each time-step. With the fitted parameters from above SUIR model, we generate the simulated evolution of the cumulative infection cases map in figure 9, which demonstrates the high consistency with the RKI data map, named as real profiles. In which figure 9(a) shows that, as of the end of March, the cumulative infection cases are presented in dynamical simulation with fitting parameters at 2020-03-05, 2020-03-15 and 2020-03-25 respectively. Correspondingly, figure 9(b) shows the cumulative infection cases in real profiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this part, we discuss the fitting details based on the mathematical simplification to the SUIR model. In the large susceptible population limit and ignoring the recovery and mortality cases for the uncertainty at early stage, there are only two equations are relevant,

$$\begin{align} \frac{d U}{d t} & = (\lambda -\sigma p- \frac{\sigma (1-p)}{d})U \end{align} \tag{ 5 }$$

$$\begin{align} \frac{d (I+R)}{d t} & = (\sigma p + \frac{\sigma (1-p)}{d} )U, \end{align} \tag{ 6 }$$

where we redefine the $\sigma p + \frac{\sigma (1-p)}{d}$ as σ. The solutions are clearly derived as,

$$\begin{align} U(t) & = U_0 Exp[(\lambda -\sigma)t] \end{align} \tag{ 7 }$$

$$\begin{align} I+R & = I_0+R_0-\frac{\sigma U_0}{\lambda -\sigma}+\frac{\sigma U(t)}{\lambda -\sigma}. \end{align} \tag{ 8 }$$

The I + R fitting is easy now with $I+R = a_* +b_* Exp[c_* t]$. This means we can not obtain all the independent parameters in the solution, but the following redefinition could be helpful,

$$\begin{align} a_* & = I_0+R_0-\frac{\sigma U_0}{\lambda -\sigma} \end{align} \tag{ 9 }$$

$$\begin{align} b_* & = \frac{\sigma U_0}{\lambda -\sigma} \end{align} \tag{ 10 }$$

$$\begin{align} c_* & = \lambda -\sigma, \end{align} \tag{ 11 }$$

or equivalently only the following combinations could be obtained:

$$\begin{align} &I_0+R_0 = a_*+b_* \end{align} \tag{ 12 }$$

$$\begin{align} &\sigma U_0 = b_* c_* \end{align} \tag{ 13 }$$

$$\begin{align} &\lambda -\sigma = c_*. \end{align} \tag{ 14 }$$

As for the product σU₀, which means it is difficult to determine the σ and U₀ individually. It will introduce the inevitable uncertainty to the two parameters, while the confidence to the multiplicative variable could be reserved.

The nuts-and-bolts of the neural network are: Convolution layers within TimeDistributed wrapper, from Convolution 2D (64, kernel size = 3 × 3, 'ReLU', 'same' padding) to Convolution 2D (64, kernel size = 3 × 3, strides = 2, 'ReLU', 'same' padding) to Convolution 2D (64, kernel size = 3 × 3, 'ReLU', 'same' padding) to Convolution 2D (64, kernel size = 3 × 3, 'ReLU', 'same' padding) and Flatten; for the LSTM part, 256 LSTM cells with 'sequence = False' and activation = 'tanh' are used; dense layers contains 256 neurons with 'ReLU' active function, before the output is the 'sigmoid' function to fit the normalized targets.