Sequential data assimilation of the stochastic SEIR epidemic model for regional COVID-19 dynamics

Time-dependence of the contact parameter

An estimation of the time-dependence of the contact parameter is achieved via the model’s best fit. An approximative instantaneous negative log-likelihood L(t_k, β) of the contact parameter β at observation time t_k is obtained from the ensemble Kalman filter (see Model inference based on sequential data assimilation). Thus, by determining the minimum of L(t_k, β) with respect to β at time t_k we estimate the time-dependence of the best fit β*(t_k) (Fig. 2c). The black line represents the average time dependence for all 320 regions included in the analysis; standard deviations are indicated by the grey area. The results for the two example regions are given in their corresponding colors.

The non-pharmaceutical interventions to counter the spread of COVID-19 were implemented at slightly different points in time across Germany. In the majority of the regions, closings of schools and other educational institutions started on March 16, while large-scale contact bans were implemented on March 22. As a result of these social distancing measures will have an impact on the contact parameter, we expected to observe a related drop in the contact parameter over time. Before we present a corresponding analysis, it should be made clear that none of these measures can produce an immediate effect on the observed cases of infected individuals because of the latency period. Since sequential data assimilation will need several data points to adapt the model to the data, the related interval should be as long as possible to achieve a reliable estimate of the contact parameter. Therefore, we selected the average value of β*(t_k) over the three days from March 17th to March 19th as a pre-intervention value. The average over March 31st to April 2nd is taken as an estimate of the post-intervention value. To analyze the effect across regions, we computed average values β_pre (March 17-19) and β_post (March 31-April 2) of the relevant β*(t_k) for all regions. The resulting scatter plot indicates a clear reduction of the numerical value of the contact parameter from β_pre to β_post (Fig. 2d). The reduction is statistically significant (Wilcoxon test, p < 0.01).

Simulations with time-dependent contact parameter

The contact parameter β is the most critical parameter determining the dynamics of the stochastic SEIR model. After a time-resolved estimation of the best fit β_∗(t_k), we are able to generate simulations from an initial state to predict the future trajectory (Fig. 3). Simulations I begin with the first epidemic day in the corresponding region with greater than or equal to 30 cases. The initial number of infected I₀ is set to the observed number of cases y_obs(t₀), while the initial number of exposed is set to E₀ = g/a · I₀, which would hold at epidemic equilibrium for m > 0 (for m = 0 both E and I tend to zero with E/I ≈ g/a). The initial number of infected was disturbed by noise representing uncertainties in the initial model states. Forward iterations with the estimated time-varying contact parameter show that the slope of the epidemic curve is approximately reproduced by the model (Fig. 3a,c; grey lines indicate the ensemble of simulated trajectories; blue points are observed data).

• Download figure
• Open in new tab

Figure 3:

Simulations of the stochastic SEIR model for two example regions. Simulations I indicate an ensemble of 100 runs of the model with initial conditions from the first epidemic day with number of cases greater than or equal to 30 (grey: ensemble of trajectories; blue: observations). Simulations II start at March 26, using an ensemble size of 100 after data assimilation (grey: ensemble of trajectories; red: observations). (a) Cumulative cases of infected individuals over time for LK Köln. Daily reported new cases for Köln. (c) Cumulative cases for LK Münster. (d) Daily new cases for Münster. Date refers to report of case at RKI.

Simulation II starts at March 26 and exploits the full potential of sequential data assimilation. The sequential data assimilation approach via the ensemble Kalman filter (see Ensemble Kalman filter) is based on the forward modeling of an ensemble of trajectories. After each time step (1 day), the ensemble of trajectories is compared to the next observation and adjusted via a linear regression step. Thus, we obtained an adapted ensemble of internal model states for each epidemic day. Here we exploit this fact to run a forward simulation with initial conditions from the assimilated ensemble of internal model states. The corresponding forward simulations are close to the real time-evolution of the epidemics in the two example regions (Fig. 3a,c; grey lines indicate the ensemble of simulated trajectories; blue points are observed data). A related plot of the daily reported new cases indicates an approximately constant level of numbers of new cases for Köln (Fig. 3b) and slowly decreasing daily new cases for Münster (Fig. 3d); both predictions are in agreement with empirical observations.

Predictions for two different scenarios

The forward simulations discussed in the previous section demonstrated the predictive power of the SEIR model when using sequential data assimilation. In the next step, we generated simulations under two different scenarios. In scenario I, we started with the adapted ensemble of internal model states after data assimilation (April 16) and iterated the model forward with the mean contact parameter estimated from the period of April 14 to April 16, that is well after interventions were implemented (Fig. 4, green area). The simulations continue to match the time-course of infected cases for both example regions (Fig. 4a,b). Daily reported case numbers show a decline for both regions (Fig. 4c,d).

• Download figure
• Open in new tab

Figure 4:

Model predictions for COVID-19 after data assimilation in comparison to data (black lines). In scenario I (green area), an assimilated ensemble of internal model states starts the forecast with contact parameter β_post (continuation of social distancing interventions). In scenario II (red area), the equivalent forecast is generated with contact parameter β_pre (termination of interventions). Predictions for cumulative case numbers in Heinsberg. (b) Predictions of daily new cases in Heinsberg. (c) Predictions of cumulative cases for Warendorf. (d) Predictions of daily new cases for Warendorf. Date refers to report of case at RKI.

In scenario II, we assumed that all governmental intervention measures had been terminated. Therefore, we used the estimated contact parameters from the period of March 17 to March 19. Again, we started simulations with the adapted ensemble of internal states after sequential data assimilation (Fig. 4, red area). For both example regions, we observe a strong increase in infected cases under scenario II (Fig. 4c,d). This dramatic increase can be seen most clearly in the plot of daily numbers of new cases (for more examples see Supplementary Information Appendix).

RKI data on COVID-19 in Germany

The Robert Koch Institute (RKI), the central scientific institution in the field of biomedicine within the portfolio of the Federal Ministry of Health, provides daily access to the number of confirmed cases, deaths, and recovered patients, broken down into 412 counties, six age groups, and by gender. As they are official records, only cases certified by doctors or labs according to a strict medical protocol in accordance with the Infection Protection Act are entered into the data base. The exact time of an infection is usually not known. The associated time stamp refers to the date on which the local health authority became aware of the case and recorded it electronically. As records are passed from the physician or lab via local and state health authorities to the RKI, there is a delay of several days before cases are reported on the website. Thus, the statistics relating to the most recent three or four days are incomplete and cannot be interpreted; retrospective updates and corrections are made for all days of the pandemic spell as they become available. Also the date at which the case is reported at the RKI (i.e., the Erkrankungsdatum/date of illness) is not more recent than the date at which the infection occurred. We use data up to and including April 30, as reported on June 22, 2020; they are included as part of the supplement.

SEIR model and basic reproductive number

The SEIR epidemic model is a four-compartment model with susceptibles (individuals who are able to contract the disease), exposed (individuals who are infected but not yet infectious), infectious, and recovered (individuals who are immune). The model is typically formulated as a system of ordinary differential equations (ODE), i.e., where the total number of individuals N = S + E + I + R is constant under temporal evolution due to . The ODE system, Eq. (1-4), has a non-trivial equilibrium point, denoted as epidemic equilibrium (S₀, E₀, I₀, R₀), where the number of susceptibles S₀ at equilibrium is related to the basic reproductive number [8]. Since we aim at a short-term description of the system, we neglect birth and death processes here, equivalent to the limit m → 0, we obtain

We use a numerical values of g = 1/3 per day, equivalent to an average infectious period of D = 3 days [16]. As a consequence, the critical condition for disease containment R₀ < 1 is obtained for β < β_crit = 1/3 per day in our model. The median latency period has been estimated as 5.2 days [13]. Here we used a numerical value of a = 1/5.2 per day for the rate parameter of the exposed individuals.

The stochastic SEIR model

While the classical model is formulated as a system of ordinary differential equations, we are exploring its application to a relatively low number of cases in the early phase of the current epidemics on the regional level with population sizes from 10⁵ to 10⁶. Therefore, we use the stochastic SEIR model in the form of a master equation [9], which is particularly useful for modeling small numbers of infected individuals occuring in smaller regions or at the beginning of an epidemic.

The demographic SEIR model contains four variables denoted by X = (S, E, I, R)^T ∈ ℕ⁴, representing the number of individuals in each of the four classes with a constant population size N = S + E + I + R. The transition rate of the ODE compartmental model translates into transition probabilities in the master equation formulation for the evolution of the model’s conditional probability density, that is, with transition probabilities given in Table 1 and initial condition X₀. The individual trajectories for the model’s temporal evolution can be simulated exactly and numerically efficiently [9] using Gillespie’s algorithm [11].

Model inference based on sequential data assimilation

Publicly available data on the cumulative number of infected individuals is used to infer the model states X = (S, E, I, R)^T and the contact parameter β of the stochastic SEIR model. Note that the cumulative number of infected individuals corresponds to Y = I + R in the SEIR model.

In the present study, we combine sequential data assimilation for the model states with an approximate log-likelihood function for the contact parameter [20]. The basic algorithmic idea is to propagate an ensemble of M model forecasts using Gillespie’s algorithm up to the next available observation point t_k. The forecast ensemble is denoted by with n ∈ {1, …, M}. We used an ensemble size of M = 100 in this study. The reported cumulative number of infected individuals y (t) is then used via a linear regression approach to obtain the adjusted model states . This step is implemented via the ensemble Kalman filter in its formulation of [21, 20]. While the forecast ensemble is used to compute the temporary negative log-likelihood L(t_k, β) of the model’s contact parameter β at time t_k, the adjusted model states serve as starting values for the next Gillespie prediction cycle.

The above algorithm is run with a fixed range of contact parameters β ∈ [β_min, β_max] and for a fixed time window [t_initial, t_final] of available data points y_obs(t_k). The best fit contact parameter β_∗(t_k) at any time any t_k is the one that minimises the temporary negative log-likelihood function, that is, with L(t_k, β) defined by (13) below.

Ensemble Kalman filter

The reported cumulative number of infected individuals y_obs(t_k) is linked to the model states X = (S, E, I, R)^T via i.e., H = (0, 0, 1, 1). As an initial condition, we set I(t₀) equal to the reported number of infected cases and R(t₀) = 0 so that y_obs(t₀) = I(t₀)+R(t₀), and E(t₀) = g/a·I(t₀) with additive noise. We assume that the errors in the observed y_obs(t_k) is additive Gaussian with mean zero and variance ρ. We set ρ = 10 in our experiments. The analysis step of the ensemble Kalman filter is now based on the empirical mean and the empirical covariance matrix of the forecast ensemble. These two quantities are used to quantify the forecast uncertainty. Combining the forecast uncertainty with the assumed data error model leads to the linear regression formula [21, 20] with the Kalman gain defined by

The resulting analysis ensemble , can be mapped back onto the integers ℕ⁴ if needed.

Model evidence

The model’s negative log-likelihood at an observation time t_k is approximated by

Note that the first contribution penalises the data misfit while the second penalises model uncertainty. The smaller the negative log-likelihood, the better the chosen model parameter β fits the data y_obs(t_k) at time t_k. The best parameter fit over a time window [t_min, t_max] is defined as the value of β which minimises the cumulative negative log-likelihood that is,

Parameter recovery from simulated data

To test the inference scheme, we simulated data for 20 days. In Figure 5a, the black line indicates the evolution of the SEIR model’s predicted cumulative numbers of infected individuals, Y (t_k) = HX(t_k) = I(t_k) + R(t_k). As in the real data, red dots represent the daily number of reported cases. In the simulation, the contact rate was chosen as β_true = 0.6. In the following, we analyzed whether this true value could be recovered using the inference procedures described above.

• Download figure
• Open in new tab

Figure 5:

Parameter recovery analysis. (a) Simulated data with b = 0.6. (b) Negative log-likelihood L_cum(β) indicates a minimum at about the true parameter value.

We varied the contact rate β and determined the cumulative negative log likelihood values L_cum(β), Eq. (14). The position of the minimum of L_cum(β) indicates the best estimate for the numerical value of the underlying contact rate β_∗, Eq. (15). The position of the minimum turns out to be close to the true value, β_∗ ≈ β_true = 0.6 (Fig. 5b). Thus, parameter recovery can be demonstrated for a relatively short time series of 10 observations, which represents a typical data-set in the early phase of newly emerging epidemics. Next, we apply our inference scheme to real data.

Application to empirical data

Since the parameter inference was successful with simulated data, the next step was an application to empirical observations. We applied the inference framework to two regional data sets from the RKI data base. As an example, we selected the COVID-19 time series for Köln (RKI data, population size N = 1, 085, 664), which includes 27 days of observations with more than 30 cases and is plotted in Figure 6a. The parameter estimation yields an estimate for the contact rate of β_∗ ≈ 0.7 (Fig. 6b). Thus, an analysis of the negative log-likelihood function produced qualitatively similar results for the simulated SEIR time series and the empirical data for a representative region. In the main text, we carry out an estimation of the time-resolved instantaneous optimal parameter values β_∗(t_k), Eq. (7), using the instantaneous negative log-likelihood function L(t_k, β), Eq. (13).

• Download figure
• Open in new tab

Figure 6:

Contact parameter estimates for real data. (a) Data for Köln; date refers to report of case at RKI. (b) Negative log-likelihood L_total(β) for Köln give a minimum at β_∗ ≈ 0.7.

We found that our results were relatively insensitive to the choice of the measurement error variance ρ appearing in (12) and (13). At the same time, we emphasise that the errors in the reported cumulative numbers of infected individuals are complex, may vary over time, and will certainly impact on the inferred parameters. The same applies to the unknown initial model states X(t₀) = (S(t₀), E(t₀), I(t₀), R(t₀))^T and their uncertainties.