Author summary↩︎

Mechanistic epidemiological models played an important role during the COVID-19 pandemic by helping estimate disease transmission, forecast healthcare demand, and evaluate intervention strategies. These models require repeated calibration to surveillance data as epidemic conditions change. Bayesian inference provides a principled framework for such calibration while quantifying uncertainty, but standard approaches based on Markov chain Monte Carlo (MCMC) are often computationally expensive and difficult to apply in near-real-time settings.

In this study, we investigated simulation-based inference (SBI), a machine-learning-based approximate Bayesian approach, as a faster alternative for parameter estimation in a mechanistic SECIR epidemiological model using COVID-19 ICU occupancy data from Germany. We compared SBI directly with conventional MCMC across multiple epidemic phases and both short and long inference windows. SBI recovered posterior parameter distributions that closely matched MCMC while substantially reducing runtime from approximately 1000 seconds to about one minute on a single GPU for inference using 31-day time series. Furthermore, SBI scaled much better for inferring parameters from a 201-series, requiring approximately three minutes, compared to 19,000 seconds for MCMC. We additionally quantified posterior agreement using distributional similarity metrics and posterior predictive checks.

Our results show that SBI can provide rapid and scalable Bayesian calibration for epidemiological models, supporting repeated inference and rapid outbreak analysis in time-sensitive public-health settings.

Introduction↩︎

Mechanistic epidemiological models are widely used to describe infectious-disease dynamics and to support public-health decision making. Representing key biological and epidemiological processes explicitly, such models have been used to estimate transmission dynamics, project epidemic trajectories, assess intervention strategies, and anticipate healthcare-system burden [1]–[4]. During the COVID-19 pandemic, mechanistic models played a central role in forecasting hospital and intensive care unit (ICU) occupancy and in evaluating the potential effects of non-pharmaceutical interventions, [5], [6].

Repeated calibration is particularly important [2], [7], [8], as transmission patterns, contact behavior, testing practices, healthcare utilization, and intervention measures may change rapidly during an ongoing outbreak and model parameters inferred from one epidemic phase may not be accurate for longer than a few days or weeks. Public-health decisions based on projected hospital or ICU occupancy therefore require inference methods that are capable of updating parameter estimates and projections, including uncertainty, on timescales compatible with operational decision-making. Delays in model calibration risk reducing the value of model-based forecasts, especially when epidemic dynamics change rapidly [9].

Bayesian inference provides a principled framework for calibrating epidemiological models because it combines prior knowledge, observational data, and uncertainty into interpretable inferences, thereby enabling forecasts [10]. In practice, Bayesian calibration is commonly performed using Markov chain Monte Carlo (MCMC) methods, which generate samples from the posterior distribution through repeated evaluations of the likelihood. MCMC-based inference has been applied successfully to a wide range of infectious-disease modeling problems, including influenza, Ebola, and COVID-19 transmission [11]–[13].

Despite its strong theoretical guarantees, MCMC is computationally expensive for complex epidemiological models [12], [14], sometimes taking months for detailed analyses [15]. Modern epidemic simulators often involve high-dimensional parameter spaces, nonlinear dynamical systems, stochastic processes, or computationally intensive numerical solvers [16]. As a result, posterior sampling with MCMC, which may require thousands to millions of likelihood evaluations, becomes prohibitively slow for repeated analyses. This computational burden limits the applicability of MCMC for repeated near-real-time analyses during rapidly evolving outbreaks.

Simulation-based inference (SBI) has emerged as a promising alternative to MCMC within Bayesian inference for simulator-based models [16], [17]. Instead of relying on repeated likelihood evaluations, SBI learns a probabilistic mapping from simulated data to model parameters using neural density estimation. Neural posterior estimation (NPE), in particular, trains a conditional density estimator on simulated parameter-data pairs and subsequently evaluates the learned posterior approximation on observed data. Once trained, SBI can rapidly generate posterior samples, making SBI attractive for applications that require repeated inference and fast recalibration.

Recent studies have started exploring SBI for epidemiological applications, including transmission-parameter estimation, outbreak reconstruction, and forecasting in stochastic and compartmental epidemic settings [18]–[20]. However, the use of learned posterior approximations in these scenarios introduces new methodological challenges: the quality of the simulation design, prior specification, simulation budget, neural-network architecture, and training strategy. Moreover, accurate posterior predictive trajectories do not necessarily imply that the inferred posterior distribution over parameters is accurate, because different parameter combinations can produce similar epidemic dynamics. Careful validation against established Bayesian inference methods is therefore essential before SBI can be used as a substitute for conventional Bayesian approaches [21].

Here, we investigate NPE as a simulation-based alternative to MCMC for Bayesian calibration of a mechanistic (SECIR) model fitted to COVID-19 ICU occupancy data from Germany. We focus on ICU occupancy because of its direct relevance for healthcare-system burden as well as its lower sensitivity to reporting practice and testing intensity. The inference problem includes epidemiological disease progression parameters as well as time-varying effective contact-rate parameters that capture changes in transmission conditions over time.

We compare SBI and MCMC across multiple epidemic phases and inference horizons. First, we evaluate three operationally relevant 31-day time windows representing distinct phases of the German COVID-19 epidemic in 2020. These shorter time windows reflect an operational setting in which models may need to be recalibrated repeatedly, as new data arrive. Second, we consider a substantially more challenging 201-day reconstruction problem involving multiple transition change points and a higher-dimensional parameter space. We assess agreement between SBI and MCMC using posterior predictive checks as well as quantitative comparisons of marginal posterior distributions based on Wasserstein distance (WD) and Kullback–Leibler divergence (KLD). By combining predictive evaluation with posterior comparison, we aim to determine whether SBI can provide a computationally efficient approximation to MCMC-based Bayesian calibration while preserving the posterior structure relevant for epidemiological interpretation and uncertainty quantification.

Materials and methods↩︎

The epidemiological model↩︎

We use a mechanistic compartmental model to describe the dynamics of COVID-19 infections and ICU occupancy in Germany in 2020. The model structure follows the SECIR-formulation introduced in [22] and is implemented in the MEmilio framework [23]. The total population is partitioned into eight epidemiological compartments: susceptible (\(S\)), exposed (\(E\)), infectious with no symptoms (\(I_{NS}\)), infectious with symptoms (\(I_{Sy}\)), infected severe (\(I_{Sev}\)), infected critical (\(I_{Cr}\)), recovered (\(R\)), and dead (\(D\)) (see Fig. 1). We use ordinary differential equations (ODEs) to describe the transitions of individuals between these compartments. The model is well suited for modeling the early phase of the COVID-19 pandemic, since represents a largely susceptible population without prior immunity, separates non-symptomatic and symptomatic population, and explicitly tracks disease states, thereby providing important indicators for decision-making. In particular, the \(I_{Cr}\) compartment is directly linked to the observed ICU occupancy.

The ODE system is parameterized by transition times (\(T\)), transition probabilities (\(\mu\)), and transmission-related scaling factors (\(c, \rho, \xi\)). Transition times determine the average duration spent in epidemiological states, while transition probabilities determine the fraction of individuals progressing to the next disease state. Transition is governed by the effective contact rate \(c\) and the transmission probability upon contact \(\rho\), and the relative infectiousness of non-symptomatic (\(\xi_{I_{NS}}\)) and symptomatic (\(\xi_{I_{Sy}}\)) individuals. Individuals without symptoms either recover or progress to symptomatic infection. Symptomatic individuals then either recover or progress to severe and critical disease states, where they either recover or die. The full list of the 15 model parameters is given in 1).

To represent changes in transmission conditions over time, the model allows time-dependent reduction to modify the baseline contact rate upon behavioral changes or non-pharmaceutical interventions (NPIs). Let \(c_0\) denote the baseline contact rate before the change in transmission conditions. If an intervention is introduced at time \(t_1 > t_0\), where \(t_0\) is the start of the simulation, the initial contact rate \(c_0\) is modified by a reduction factor \(r\), where \(r\) ranges from 0 (no reduction of contacts) to 1 (complete suppression of contacts). To avoid discontinuities in the ODE system, changes in the contact rate are not applied instantaneously. Instead, the contact rate is smoothed over a short transition interval \(\delta \in (0, 1)\) using differentiable interpolation as described in [22]. The time-varying contact rate \(c(t)\) is thus defined as:

\[\begin{align} \label{eq:damping} {c}(t) := \begin{cases} c_0, & t \leq t_1 \\ \widehat{c}(t), & t \in (t_1, t_1 + \delta) \\ (1-r)c_0, & t \geq t_1 + \delta \end{cases} \end{align}\tag{1}\] where \(\hat{c}(t)\) denotes the smooth interpolation between \(c_0\) and \((1-r)c_0\). The large number of simulations required for both SBI as well as for the likelihood evaluations for MCMC were performed using the C++ implementation of the SECIR model in the MEmilio framework [23].

Data↩︎

The reported ICU occupancy of COVID-19 patients [24] in Germany from April 24 to December 10 2020 serves as the empirical basis for comparing SBI and MCMC powered inference. We use the ICU occupancy rather than reported daily case count because case counts are more strongly affected by changes in testing intensity, reporting delays, and under-reporting [25], [26]. Furthermore, ICU occupancy is also directly relevant for healthcare-system burden and public-health decision making. Let \(x^{obs}_{1:231}=(x^{obs}_1,\dots,x_{231}^{obs})\) denote the observed ICU occupancy time series over this period. For our inference experiments below, we consider different time windows \(x^{obs}_{t_0:t_{end}}\), where \(t_0\) is the respective starting time and \(t_{end}\) is the end time.

Bayesian parametrization↩︎

Estimated parameters↩︎

The Bayesian model was designed to infer epidemiological parameters together with a time-varying effective contact rate. Specifically, the Bayesian model estimates the five transmission-time parameters \(T_E\), \(T_{I_{NS}}\), \(T_{I_{Sy}}\), \(T_{I_{Sev}}\), and \(T_{I_{Cr}}\) as well as the relative infectiousness of symptomatic individuals \(\xi_{I_{Sy}}\). In addition, the Bayesian model infers changes in the time-varying contact rate on a grid. Let \(\Delta\) be the minimum number of days for which a contact-rate level is valid, then for a time series starting at day \(t_0\) and ending after day \(t_{end}\) the number of intervals on the grid is \(K=\left\lfloor{\tfrac{t_{end}-t_0+1}{\Delta}}\right\rfloor\). For \(\tau_k=t_0+k\Delta, k=0,\dots,K\) the unsmoothed effective time-varying contact rate \(\beta_{\mathrm{pc}}(t)\) is given by the piecewise-constant function \[\label{eq:effective95dampening} {\beta_{\mathrm{pc}}}(t) := \rho \cdot c_0 \cdot (1-r_k) \, \text{for}\, t\in[\tau_{k-1}, \tau_{k}) ,\qquad k=1,\dots,K,\tag{2}\] where the \(r_k\) are reduction factors for each of the \(K\) intervals. The effective contact rate used in the ODE system is obtained by smoothing transitions between consecutive levels of \(\beta_{\mathrm{pc}}(t)\) over a short transition as described above (Eq. 1 ). As before, larger values for \(r_k\) correspond to stronger reduction in the effective contact rate (1 representing reduction to 0 contacts and 0 representing no reduction of contacts). To avoid non-identifiability of the baseline contact rate and the transition probability on contact, we fix \(c_0=1\) and \(\rho=1\), meaning that for each interval on the grid the effective contact rate is controlled by a single parameter \(r_k\).

The full parameter vector of the Bayesian model is given by \(\theta\). \[\theta=(T_E, T_{I_{NS}}, T_{I_{Sy}}, T_{I_{Sev}},T_{I_{Cr}}, \xi_{I_{Sy}}, r_1,\dots,r_K)\]

Model Initialization↩︎

To initialize the epidemiological compartments for inference starting at time \(t_0\) we use the initialization procedure described in [27]. This procedure combines reported case data, ICU occupancy data and population data. Reported confirmed cases [28] are mapped to the exposed, (a)symptomatic infectious, severe, critical, recovered and deceased compartments by shifting the case time series according to assumed transition times and weighting the resulting contributions with the corresponding disease progression probabilities. ICU occupancy was obtained from the DIVI intensive care registry [24]. Population data [29] is used to set the susceptible compartment as the remaining population after all other compartments have been initialized.

Prior distributions↩︎

Because ICU occupancy provides information on only one of the eight epidemiological compartments, prior information is required to constrain weakly identifiable disease-progression and transmission parameters. We therefore used literature-informed priors for the epidemiological parameters and bounded uniform priors for the time-varying reduction factors [30]–[32]. The prior distributions are summarized in 2, while the priors on the reduction factors for 201-day window are summarized in , Spec B.

Prior and posterior predictive checks↩︎

Prior and posterior predictive checks were used to assess the consistency of the Bayesian model before and after conditioning on the observed ICU occupancy data.

For prior predictive checks, 1000 parameter samples were drawn from the prior distribution and propagated through the SECIR model to generate prior predictive trajectories. These trajectories were compared with the observed ICU occupancy data to evaluate whether the prior distributions generated epidemiologically plausible epidemic dynamics.

For posterior predictive evaluation, 1000 posterior samples were drawn from the inferred posterior distribution. Each sampled parameter vector was propagated through the SECIR model to generate one simulated ICU occupancy trajectory. The resulting ensemble of simulated trajectories was used to approximate the posterior predictive distribution. For each time point, we computed the posterior predictive mean and the pointwise 95 % credible interval. As a quantitative measure of predictive performance, we calculated the root mean squared error (RMSE) between the trajectory generated from the posterior predictive mean trajectory and the observed ICU occupancy data.

Prior predictive checks were used to assess the suitability of the prior specification, whereas posterior predictive checks were used to evaluate whether the inferred posterior distributions reproduced the observed epidemic dynamics.

Markov chain Monte Carlo reference inference↩︎

To compute baseline inferences for comparison we use Markov chain Monte Carlo (MCMC). MCMC is especially suited to create such baseline inferences as it has strong theoretical guarantees that samples will be asymptotically distributed according the targeted posterior distribution.

Observation model and likelihood↩︎

For MCMC inference, we used an explicit likelihood function relating simulated ICU occupancy to observed ICU occupancy. Following [33], we assumed Student’s t-distributed observation errors with \(\nu=4\) degrees of freedom. Student’s t-distribution is widely used to construct likelihoods for noisy observations as it’s less sensitive to outliers. For an observation window from \(t_0\) to \(T\), the likelihood was \[p(x^{obs}_{t_0:T}|\theta,t_{0:T})= \prod_{i=t_0}^{T} p(x^{obs}_i|\theta, t_i)= \prod_{i=t_0}^{T} \text{Student-}t_{\nu}(x^{obs}_i|x^{sim}(\theta, t_i), \sigma_i)\] where \(x^{sim}(\theta, t_i)\) denotes the simulated ICU occupancy at day \(t_i\) and the observation noise modeled as \(\sigma_i=\sqrt{|x^{sim}(\theta, t_i)|}\).

MCMC reference inference↩︎

We used MCMC as the reference Bayesian inference method against which SBI was compared. MCMC is suitable for this purpose because, under standard regularity conditions and sufficient convergence, it provides samples from the targeted posterior distribution. Posterior sampling was performed using the differential-evolution Metropolis-Z algorithm (DE-MZ) [34]. DE-MZ is an adaptive, gradient-free MCMC method that uses the history of the MCMC chains to adapt proposals/samples to the scale and correlation structure of the posterior during sampling. Both DE-MZ and SBI do not require the evaluation of the likelihood gradient, which is computationally expensive for the SECIR ODE system.

We used the DE-MZ implementation provided by PyMC v5.0.0 [35]. Tuning is implemented by PyMC. As a rule, for \(n\) samples we want to draw, we sample an additional \(n\) samples beforehand for tuning and discard them once tuning is done. Convergence was assessed using rank-normalized \(\hat{R}\) diagnostic [36] and check that the value is close to 1.

Simulation Based Inference↩︎

Simulation-Based Inference (SBI) leverages Artificial Intelligence to construct an approximate Bayesian framework. Here, we use the Neural Posterior Estimation (NPE) method [37], which approximates the posterior distribution of the model parameters.

Amortized Neural Posterior Estimation.↩︎

We used simulation-based inference with Neural Posterior Estimation (NPE) to approximate the posterior distribution over model parameters. NPE learns a conditional density estimator from simulated parameter–data pairs and subsequently generates approximate posterior samples for observed data without requiring repeated likelihood evaluations.

The objective is to approximate the Bayesian posterior distribution

\[p(\theta \mid X) \propto L(X \mid \theta)\pi(\theta),\]

where \(L(\cdot)\) denotes the likelihood function and \(\pi(\cdot)\) the prior distribution. Rather than evaluating the likelihood explicitly, NPE trains a neural density estimator using simulations generated from the prior predictive distribution.

Training data were generated by sampling parameter vectors \(\theta_j \sim \pi(\theta)\) and propagating them through the SECIR simulator to obtain synthetic ICU occupancy trajectories \(\hat{X}_j = f(\theta_j)\), where \(f\) denotes the epidemiological simulator. This produced a set of simulated parameter–data pairs \((\theta_j,\hat{X}_j)\) used to train the neural density estimator.

The network learns a conditional posterior approximation \(\tilde{p}(\theta \mid \hat{X})\), where the simulated trajectory serves as conditioning information and the distribution over model parameters is represented directly. In this way, the network learns a mapping from simulated epidemic trajectories to posterior distributions over parameters.

The approach is amortized because the posterior estimator is trained once using a large collection of simulated parameter–data pairs and can subsequently be applied to any new observation generated by the same simulator and prior distribution. After training, the observed ICU occupancy data \(X\) are provided as conditioning information and parameter samples are drawn from the learned posterior approximation \(\tilde{p}(\theta \mid X)\). A schematic overview of the procedure is shown in Fig. 2.

Neural Density Estimator and Trajectory Embedding↩︎

The posterior distribution was approximated using a Masked Autoregressive Flow (MAF) [38], a normalizing-flow density estimator commonly used for neural posterior estimation. In the neural posterior estimation framework, MAF models the conditional posterior distribution autoregressively as \[p(\theta \mid \hat{X})= \prod_{i=1}^{d} p(\theta_i \mid \theta_{1}, \hat{X}),\] where \(\hat{X}\) denotes the simulated ICU occupancy trajectory used as conditioning information and \(d\) is the dimension of the parameter vector. Additional mathematical details on normalizing flows, MAF, and the underlying Masked Autoencoder for Density Estimation (MADE) architecture are provided in the Supplementary Information, .

Because the simulated ICU occupancy trajectories are time series, we employed a Convolutional Neural Network (CNN) embedding to extract a lower-dimensional representation before conditioning the density estimator. Rather than providing the full trajectory directly to the MAF, the embedding network maps the simulated trajectory \(\hat{X}\) to a learned feature representation \(h(\hat{X})\), which is subsequently used as conditioning information. CNNs are well suited for extracting informative local and multiscale features from one-dimensional temporal data [39]. The final architecture consisted of a four-layer CNN embedding network followed by a MAF density estimator with ten hidden features and two flow transformations.

Training was performed using the Sequential Neural Posterior Estimation (SNPE-C) objective [37], [40]. The loss combines a contrastive (atomic) term with a maximum-likelihood regularization term, \[\mathcal{L}_{\mathrm{atomic}} + \mathcal{L}_{\mathrm{MLE}},\] where the regularization term reduces density leakage outside the support of the prior distribution. Additional details on the SNPE-C objective are provided in the Supplementary Information ().

SBI implementation↩︎

Parameter inference was performed using the sbi Python library (version 0.22.0; [41]). Within the memilio framework, we developed a dedicated SBI wrapper compatible with the simulate_for_sbi interface and parameter priors specified through MultipleIndependent distributions.

The simulated ICU occupancy trajectories were embedded using a convolutional neural network (CNN) constructed via the CNNEmbedding function. To assess the effect of the embedding architecture, we evaluated CNNs with between two and four convolutional layers and output channel sizes ranging from 6 to 12 channels per layer. Candidate architectures included channel configurations such as \([6,12,12]\) and \([6,12,12,12]\). The resulting embedding was used as conditioning information for a MAF density estimator with ten hidden features and two flow transformations. Posterior inference was performed using the SNPE-C algorithm.

To investigate the sensitivity of SBI to training hyperparameters, we evaluated simulation budgets of 20000, 50000, and 100000 forward simulations in combination with batch sizes of 450, 900, 1800, 3600, 7200, 14400, 22500, 45000, and 90000.

To improve training performance for large simulation datasets, we implemented a custom PyTorch Dataset that uses vectorized tensor indexing and avoids repeated preprocessing while preserving the standard SNPE-C training and inference workflow. Training was performed on a GPU using the custom dataloader.

After training, the observed ICU occupancy trajectory was supplied as conditioning information to the learned posterior estimator, enabling sampling from the approximate posterior distribution over model parameters.

Comparison of posterior distributions↩︎

To quantify agreement between posterior distributions obtained via simulation-based inference (SBI) and Markov chain Monte Carlo (MCMC), we compared the marginal posterior distributions of each inferred parameter using the 1-Wasserstein distance and the Kullback–Leibler (KL) divergence.

For a given scalar parameter \(\theta_j\), let \(p(\theta_j)\) denote the marginal posterior density inferred by SBI and \(q(\theta_j)\) the corresponding marginal posterior density inferred by MCMC. Since both posteriors are represented by samples, continuous density estimates \(\hat{p}(\theta_j)\) and \(\hat{q}(\theta_j)\) were obtained using Gaussian kernel density estimation (KDE). The densities were evaluated on a uniform grid spanning the combined support of both sample sets and normalized numerically.

0.0.0.1 Wasserstein distance.

The first-order Wasserstein distance (Earth Mover’s Distance) between the two marginal posterior distributions was computed as \[W_1(p,q)=\inf_{\gamma \in \Pi(p,q)} \mathbb{E}_{(x,y)\sim\gamma} \left[|x-y|\right],\] where \(\Pi(p,q)\) denotes the set of all joint distributions with marginals \(p\) and \(q\). In practice, this quantity was computed directly from the posterior samples using a standard empirical estimator.

0.0.0.2 Kullback–Leibler divergence.

The KL divergence from the SBI posterior to the MCMC posterior was defined as \[D_{\mathrm{KL}}(p(\theta_j)||q(\theta_j))=\int \hat{p}(\theta_j) \log \frac{\hat{p}(\theta_j)}{\hat{q}(\theta_j)} d\theta_j .\]

The integral was approximated numerically on the discretized grid, \[D_{\mathrm{KL}}(p(\theta_j)||q(\theta_j)) \approx \sum_{k=1}^{K} \hat{p}(\theta_{j,k}) \log \frac{\hat{p}(\theta_{j,k})}{\hat{q}(\theta_{j,k})} \Delta\theta ,\] where \(\Delta\theta\) denotes the grid spacing. To improve numerical stability, density values were clipped to a small positive constant before evaluating the logarithm.

0.0.0.3 Symmetric KL divergence.

Because the KL divergence is asymmetric, we additionally computed the symmetric KL divergence, \[D_{\mathrm{sym}}\!\left(p(\theta_j),q(\theta_j)\right) = \frac{1}{2} \left[ D_{\mathrm{KL}}\!\left(p(\theta_j)\|q(\theta_j)\right) + D_{\mathrm{KL}}\!\left(q(\theta_j)\|p(\theta_j)\right) \right].\] All metrics were computed independently for each inferred parameter. For each inference configuration, we report the mean and standard deviation of the resulting distances across 16 independent SBI runs.

Results↩︎

Experimental setups↩︎

We compared SBI with MCMC of the SECIR model to COVID-19 ICU occupancy data from Germany. Inference was performed for three 31-day inference windows and one 201-day inference window representing distinct epidemic phases, including a declining post-first-wave phase of June 2020 with low incidence summer period and two increasing autumn 2020 phases, started at offsets of 0, 160, and and 200 days after April 24, 2020. The extended 201-day window was used to assess performance in a long and higher-dimensional reconstruction problem involving multiple changes in transmission conditions.

For the 31-day windows, the selected SBI configuration used 50000 simulations and a batch size of 14400. For the 201-day inference problem, which involved a higher-dimensional parameter space and multiple change points, the selected configuration used 100000 simulations and a batch size of 1800. Each SBI configuration was repeated 16 times. Posterior predictive performance, posterior agreement with MCMC, and runtime are summarized in Table 3.

Prior predictive performance↩︎

Prior predictive checks were particularly important for the 201-day inference problem because the larger number of change-point parameters substantially increased the flexibility of the model. While most epidemiological parameters were constrained by literature-informed priors, broad prior bounds on the change-point parameters produced highly dispersed trajectories, including implausible epidemic dynamics and partially non-identifiable parameter configurations (Fig. , left). As a result, a substantial proportion of simulations occupied regions that were inconsistent with the observed epidemic trajectory.

The sensitivity of the model to change-point parameters is amplified over the longer 201-day reconstruction window, where relatively small parameter changes can lead to markedly different epidemic trajectories. For SBI, priors that support implausible epidemic dynamics impose an avoidable burden on training. While principled methods are being developed to meaningfully constrain priors [42], here we adopted tighter prior bounds for the change-point parameters in both the MCMC and SBI analyses. The revised specification retained sufficient variability to capture a broad range of plausible epidemic scenarios while excluding unrealistic dynamics.

Under the final prior specification, the observed ICU occupancy trajectory was well covered by the prior predictive envelope (Fig. , right), indicating that the priors were neither overly restrictive nor unrealistically diffuse.

Posterior predictive performance↩︎

Figure 3 shows representative posterior predictive checks for the selected inference settings. As a quantitative summary, Table 3 reports the mean RMSE and standard deviation across 16 independently trained SBI models.

4pt

Across the three 31-day windows, SBI reproduced the observed ICU trajectories well. Across the selected short-window configurations, the offset 0 window achieved the lowest RMSE of \(66.41 \pm 26.16\). This window corresponds to the declining phase after the first epidemic wave. The offset 160 and offset 200 windows yielded RMSE values of \(85.88 \pm 54.78\) and \(98.79 \pm 66.22\), respectively. These later windows involved more rapidly changing epidemic dynamics and showed broader posterior predictive uncertainty, but the posterior predictive trajectories still followed the observed ICU occupancy patterns. Compared with the MCMC reference posterior predictive distributions, which achieved RMSE values of 31.8, 14.4, and 18.2 for the offset 0, 160, and 200 windows, respectively, SBI generally produced broader credible intervals while retaining similar overall trajectory shapes and temporal trends.

The 201-day inference problem was substantially more challenging because the model had to reconstruct both the decline after the first wave and the subsequent increase during the second wave within a single inference window. The selected 201-day SBI configuration yielded an RMSE of \(325.72 \pm 135.78\). This error was larger than for the 31-day windows, as expected from the longer reconstruction horizon, the higher-dimensional parameter space, and the need to infer multiple change-point parameters simultaneously. Nevertheless, the posterior predictive trajectories captured the dominant temporal structure of the observed ICU data, including the prolonged low-incidence period and the sharp autumn increase in ICU occupancy. The corresponding MCMC posterior predictive distribution achieved an RMSE of 23.1 and produced narrower credible intervals, although both approaches reproduced the major epidemic phases observed in the ICU occupancy data.

Agreement between SBI and MCMC posteriors↩︎

We next compared the posterior distributions inferred by SBI with the corresponding MCMC reference posteriors. Agreement was assessed visually using marginal posterior distributions, see Fig. 4, and quantitatively using the first Wasserstein distance, KL divergence, and symmetric KL divergence, Tab 3.

For the 31-day inference windows, SBI recovered marginal posterior distributions that were broadly consistent with the MCMC reference posteriors (Fig. 4). The strongest agreement was observed for the offset 0 window (top panel). In this setting, the selected SBI configuration achieved an average first-order Wasserstein distance of \(0.75 \pm 0.51\) and a symmetric KL divergence of \(9.58 \pm 13.2\). The posterior modes inferred by SBI and MCMC were closely aligned for most epidemiological and reduction factor parameters.

Posterior agreement remained good for the offset 160 and offset 200 windows, although discrepancies increased for some parameters. In these later windows, SBI generally produced broader marginal posterior distributions than MCMC. This broadening was most visible for selected disease-duration and reduction factors. For offset 160, the selected SBI configuration yielded a Wasserstein distance of \(2.37 \pm 1.43\) and a symmetric KL divergence of \(75.79 \pm 31.35\). For offset 200, the corresponding values were \(1.89 \pm 0.45\) and \(46.44 \pm 9.58\). Despite these larger divergences, the dominant posterior regions remained aligned between SBI and MCMC.

The 201-day inference problem involved a substantially larger parameter space and showed broader posterior uncertainty than the 31-day analyses, . Nevertheless, SBI preserved the main posterior structure observed in the MCMC reference. For the selected 201-day configuration, the Wasserstein distance was \(0.35 \pm 0.35\) and the symmetric KL divergence was \(72.14 \pm 24.67\). Several reduction factor parameters showed substantial overlap between SBI and MCMC, whereas larger discrepancies were observed for parameters near the end of the observation window. These late change points were less strongly constrained by the data and therefore showed broader posterior support under SBI.

Overall, SBI recovered posterior distributions that were consistent with the MCMC reference posteriors across all inference settings. While the SBI marginals were generally smoother and somewhat broader than the corresponding MCMC distributions, the dominant posterior regions and parameter trends were preserved.

Effect of simulation budget and batch size↩︎

We evaluated the effect of simulation budget and batch size on SBI performance using posterior predictive RMSE, posterior distance metrics, runtime, and variability across repeated runs (Table 3, ).

For the 31-day inference windows, configurations trained with 20000 simulations generally showed higher RMSE and greater variability across repeated runs. Increasing the simulation budget to 100000 did not consistently improve posterior predictive performance or agreement with the MCMC reference posteriors. Across the three 31-day windows, 50000 simulations provided stable posterior recovery while maintaining substantially lower computational cost. Among the tested batch sizes, a batch size of 14400 yielded consistently low RMSE, good agreement with MCMC, and stable runtimes, and was therefore selected for the subsequent analyses.

For the 201-day inference problem, the optimal configuration differed from that observed for the 31-day windows. Because this setting involved 21 inferred parameters and multiple reduction factors, posterior recovery was more sensitive to the training setup. The configuration with 100000 simulations and batch size 1800 achieved the most stable posterior agreement with MCMC across repeated runs while maintaining competitive posterior predictive performance. Larger batch sizes did not improve posterior agreement or RMSE.

Across experiments, very small batch sizes increased both runtime and variability across repeated runs. For the 31-day inference windows, the smallest tested batch size (900) required approximately 160–175 seconds depending on the inference window, whereas larger batch sizes yielded substantially lower and more stable runtimes (Table 3). In several runs, the increased runtime was partly attributable to inefficient posterior sampling from the trained neural posterior, where low acceptance rates increased the time required to generate posterior samples. Very large batch sizes did not consistently improve posterior agreement with MCMC or posterior predictive accuracy. Overall, the optimal hyperparameter configuration differed between the 31-day and 201-day inference settings, indicating that performance depended jointly on simulation budget, batch size, and inference dimensionality.

Consistency across inference windows↩︎

We additionally compared inferred parameters across the short- and long-window analyses to assess whether consistent posterior structures were recovered across different epidemic phases. For the epidemiological parameters that are not explicitly time dependent in the model, the exposed period and the durations of severe and critical infection showed substantial overlap across the three 31-day windows, with similar posterior regions inferred by both SBI and MCMC. In contrast, the symptomatic infection duration shifted toward larger values in the offset 160 and offset 200 windows compared with offset 0, and this pattern was observed for both inference methods. The relative infectiousness of symptomatic individuals remained broadly consistent between offset 0 and offset 160, but shifted toward larger values in the offset 200 window.

The 201-day inference provided an additional consistency check because it spans the same calendar period covered by the shorter 31-day windows within a single joint inference problem. Early reduction factor parameters in the 201-day posterior exhibited similar contact-reduction behavior to the offset 0 inference window, while later reduction factors reflected the stronger transmission dynamics inferred in the offset 160 window. Although the posterior distributions were generally broader in the 201-day setting, reflecting the larger number of inferred reduction factors and stronger parameter correlations, the dominant transmission patterns remained consistent between the short- and long-window analyses.

Runtime comparison↩︎

Details for all MCMC runs are reported in and all \(\hat{R}\) values are summarized in Table . SBI substantially reduced runtime compared with the MCMC reference inference. For the 31-day windows, MCMC required approximately 1077 s per inference run, achieving a maximum rank-normalized \(\hat{R}\) value of 1.0119 across all inferred parameters. All MCMC runs were performed on an AMD EPYC 9334 32-Core Processor, whereas SBI experiments were performed on a single NVIDIA A100 GPU. In contrast, the selected SBI configuration required \(71.47 \pm 19.22\) seconds for the offset 0 window, \(60.83 \pm 3.09\) seconds for the offset 160 window, and \(60.07 \pm 3.97\) seconds for the offset 200 window.

The runtime advantage was also observed for the 201-day inference problem. The MCMC reference inference required 19102 seconds, whereas the selected SBI configuration required \(157.08 \pm 42.28\) seconds. Thus, even in the higher-dimensional long-window setting, SBI remained computationally tractable and substantially faster than MCMC. Despite the high run-time, MCMC failed to achieve low rank-normalized \(\hat{R}\) (maximum \(\hat{R}=1.605\)) for the parameters.

Overall, SBI reduced inference time by approximately one order of magnitude for the 31-day windows and more than two orders of magnitude for the 201-day inference problem while maintaining comparable posterior predictive performance and broad agreement with the MCMC posterior distributions.

Discussion↩︎

In this study, we evaluated SBI with NPE as a computationally efficient alternative to MCMC for Bayesian calibration of a mechanistic SECIR epidemiological model fitted to COVID-19 ICU occupancy data from Germany. Across multiple epidemic phases and inference horizons, SBI recovered posterior distributions that were broadly consistent with the MCMC reference. By combining posterior predictive checks with quantitative comparisons of posterior distributions, we found that SBI preserved the dominant posterior structure while substantially reducing computational cost. These results show that NPE is capable to provide a practical approximation to MCMC-based Bayesian calibration for the SECIR ICU occupancy model considered here.

A central finding of this work is that SBI maintained good posterior predictive performance across qualitatively different epidemic regimes. The three 31-day windows included both declining and rapidly increasing ICU occupancy trajectories, yet the inferred SBI posteriors generated posterior predictive simulations that closely followed the observed data in all cases. This is particularly relevant for operational epidemiological modelling, where models must be recalibrated repeatedly as transmission dynamics change. The 201-day inference problem was substantially more demanding because the model had to jointly capture the decline following the first epidemic wave, the low incidence summer period, and the subsequent autumn increase within a single inference task. Although posterior uncertainty increased and posterior predictive error was larger in this setting, the dominant temporal features of the epidemic trajectory were still reproduced, indicating that SBI remained informative even in a substantially higher-dimensional inference problem.

Due to the unidentifiability of the parameters that together make up the contact rate, the inferred reduction factors should be interpreted as changes in the effective contact rate rather than direct estimates of specific interventions. The reduction factors therefore summarize the combined effects of non-pharmaceutical interventions, behavioral adaptation, seasonal contact patterns, and other unobserved drivers of transmission [5], [6], [43]. As such, they are epidemiologically informative beyond trajectory reconstruction alone, since they can indicate shifts in transmission conditions and support timely model recalibration and assessment of epidemic dynamics. Consistent with this interpretation, the offset 0 window suggested sustained transmission suppression following the first epidemic wave. The offset 160 window indicated that transmission was already being reduced despite continued increases in ICU occupancy, reflecting delays between transmission changes and their impact on critical care demand. The offset 200 window was consistent with stronger transmission control and slower epidemic growth, although uncertainty increased for the final change point. Similar qualitative transmission patterns were recovered in the 201-day reconstruction.

Importantly, our evaluation was not limited to posterior predictive checks. Similar trajectory predictions can arise from different combinations of epidemiological and transmission parameters, particularly in partially non-identifiable compartmental models fitted to a single data set. We therefore compared SBI and MCMC posteriors directly using marginal posterior distributions, Wasserstein distances, and KL divergences. This provides a stricter assessment of SBI than predictive performance alone. Across the 31-day windows, SBI recovered the major high-density regions of the MCMC reference posterior, although the SBI marginals were often smoother and broader. The broadening was especially visible in the later epidemic windows and the 201-day inference problem, consistent with the behavior expected from amortized approximate inference, where the neural density estimator may trade posterior sharpness for smoother uncertainty representations, [16], [17]. Nevertheless, the dominant posterior structure and the posterior predictive behavior were preserved across the analyzed scenarios.

The broader SBI posteriors were also reflected in the posterior predictive checks, where SBI generally produced wider credible intervals than the MCMC reference despite similar posterior predictive means. One possible explanation is the difference in how uncertainty is incorporated during inference. The MCMC reference conditions parameter estimates through an explicit likelihood and observation error model, whereas the SBI approach was trained on deterministic model simulations and did not explicitly model observation noise. Consequently, SBI may retain posterior mass over a wider range of parameter combinations that generate plausible epidemic trajectories, leading to broader predictive uncertainty after propagation through the SECIR model. This effect was most apparent in the later epidemic windows, where uncertainty in transmission-related parameters is amplified by the nonlinear epidemic dynamics.

The comparison between time windows provides an additional consistency check of the inferred epidemiological parameters. Several parameters associated with disease progression showed substantial overlap across the different observation windows, despite representing distinct phases of the epidemic, suggesting that key disease-progression time scales were recovered similarly across changing epidemic conditions. Where differences between windows were observed, they appeared consistently in both SBI and MCMC, indicating that they likely reflected features of the data and model structure rather than artifacts of the SBI approximation.

At the same time, these differences should be interpreted with caution. Because the model is fitted to ICU occupancy alone and time-varying transmission was represented through effective contact-rate reduction factors, individual epidemiological parameters may be only partially identifiable and can be challenging to disentangle from changes in behavior, seasonality, reporting, healthcare practice, or residual model mismatch. Thus, changes in inferred symptomatic duration or symptomatic infectiousness should not be interpreted as direct evidence for biological changes in the disease progression. Rather, they should be understood as effective parameter changes within the calibrated model that help explain the observed ICU trajectory together with the inferred contact-rate parameters.

The 201-day reconstruction provides further support for this interpretation. Because it covers the periods represented by the shorter 31-day windows within a single inference problem, it allows direct comparison between short- and long-window estimates. The long-window inference recovered transmission patterns that were qualitatively consistent with those identified in the shorter windows, although posterior uncertainty increased because of the larger number of inferred change points and stronger parameter correlations. Taken together, these findings indicate that the long-window inference provides a coherent, albeit less sharply identified, reconstruction of the epidemic phases captured by the shorter-window analyses.

SBI provided a substantial computational advantage over MCMC across both the short- and long-window inference problems. SBI reduced wall-clock inference time by approximately one order of magnitude for the 31-day windows and by more than two orders of magnitude for the 201-day reconstruction. The computational advantage became particularly pronounced in the higher-dimensional long-window setting, where MCMC required several hours and failed to achieve satisfactory convergence for all parameters. In contrast, SBI remained computationally tractable while maintaining broad agreement with the MCMC posterior and posterior predictive distributions. This reduction in runtime is particularly relevant for public-health applications in which epidemiological models must be recalibrated repeatedly as new data become available. Faster approximate posterior inference can facilitate more frequent model updates, rapid sensitivity analyses, and timely evaluation of alternative epidemic scenarios.

However, our runtime comparison should be interpreted as a practical wall-clock comparison rather than a pure algorithmic benchmark, since the SBI workflow leveraged both CPU-based simulations and GPU-accelerated neural network training, whereas the MCMC reference was executed entirely on CPU hardware.

Our results also highlight several methodological considerations for applying SBI to epidemiological models. First, prior specification proved particularly important for the high-dimensional 201-day inference problem. Broad priors over multiple reduction parameters generated highly variable and partly implausible epidemic trajectories, reducing the relevance of the simulated training data and impairing posterior learning. This behavior is consistent with previous observations that SBI performance depends strongly on the quality of the prior predictive distribution and on sufficient coverage of posterior-relevant regions of parameter space during training [40], [44]. Careful prior checking is therefore remains essential, particularly for long-window reconstruction problems involving many time-varying parameters

Second, we observed that increasing the number of simulations beyond 50,000 did not consistently improve posterior agreement or predictive accuracy in the 31-day inference problems. Instead, performance depended jointly on simulation budget, batch size, and inference complexity. Interestingly, the optimal configuration differed between the short- and long-window settings: while the 31-day analyses benefited from relatively large batch sizes, the substantially more complex 201-day inference required smaller batch sizes together with larger simulation budgets to achieve stable posterior recovery. One possible explanation is that smaller batches introduce additional stochasticity during optimization, which may improve generalization of the neural density estimator in higher-dimensional inference problems, consistent with observations from the deep-learning literature [45], [46]. These findings suggest that hyperparameter selection for SBI in epidemiological applications may need to be adapted to the temporal scale and dimensionality of the inference problem rather than transferred directly between settings.

The study leaves several points for future research. First, the analysis was performed using a single mechanistic SECIR model and one national ICU occupancy time series. Further work is needed to assess the extent to which the conclusions generalize to other deterministic model structures, simulators, countries, pathogens, or surveillance settings. Second, future studies could investigate whether combining ICU occupancy with additional epidemiological data sources improves parameter identifiability and posterior stability.

Finally, the present analysis demonstrates retroperspective calibration rather than prospective forecasting. While the computational efficiency of SBI is promising for near-real-time epidemiological modelling, operational deployment would require evaluation in sequential forecasting settings with repeated updates as new observations become available.

In summary, SBI with NPE provided a fast and accurate approximation to MCMC-based Bayesian calibration for a mechanistic SECIR model fitted to German COVID-19 ICU occupancy data. Across multiple epidemic phases, SBI reproduced the observed ICU dynamics and recovered the dominant posterior structure of the MCMC reference while substantially reducing computational cost. Although posterior uncertainty increased in the more challenging long-window setting and performance remained sensitive to prior specification, SBI retained useful posterior information even in this higher-dimensional reconstruction problem. These findings support SBI as a promising approach for rapid calibration of mechanistic epidemiological models, provided that prior predictive checks, posterior validation, and careful consideration of parameter identifiability are incorporated into the inference workflow. Because SBI is an amortized inference method, its practical advantages are expected to be greatest in settings that require repeated inference under the same model and prior specification, where the initial simulation and training costs can be distributed across multiple analyses.

Acknowledgments↩︎

This work was supported by the Initiative and Networking Fund of the Helmholtz Association (grant agreement number KA1-Co-08, Project LOKI-Pandemics). The authors gratefully acknowledge computing time on the supercomputer JURECA [47] at Forschungszentrum Jülich under grant no. loki.

Supplementary information↩︎

1 Mean Autoregressive Flows↩︎

Generally, normalizing flows learn the distribution of the data \(X\) \(p_X(x;\phi)\) through a mapping \(f_{\phi}: \mathbb{R}^n\rightarrow \mathbb{R}^n\) between the latent variables \(Z\) and the data \(X\), where \(\phi\) are the parameters of the mapping. Then, the change of variables formula is used \[p_{X}(x;\phi)=p_{Z}(f^{-1}_{\phi}(x))\det|\frac{\partial f^{-1}_{\phi}}{\partial x}|,\] where \(\det|\frac{\partial f^{-1}_{\phi}}{\partial x}|\) is a Jacobian of the transformation between \(Z\) and \(X\). Typically, \(f_{\phi}\) is a chain of the transformations applied to a simple probability density function \(p_Z\), such as the one of the normal distribution. The network is then train via minimizing the negative log-likelihood \[\mathcal{L}=-\log{p_X(x)}=-\log{p_{Z}(f^{-1}_{\phi}(x))}-\log{\det|\frac{\partial f^{-1}_{\phi}}{\partial x}|}\]

MAFs model the probability density function of the data in an autoregressive way \(p_X(x)=\prod_{i=1}^np(x_{i}|x_{1:i-1})\), where \(x=(x_1,\dots,x_n)\) and \(x_{1:i}=(x_1,\dots,x_i)\). \[p(x_i|x_{1:i})\sim \mathcal{N}(x_i|m_{i,\phi},\exp{s_{i,\phi}}), m_{i,\phi}=m_{i,\phi}(x_{1:i-1}), s_{i,\phi}=s_{i,\phi}(x_{1:i-1})\]

where \(m_{i,\phi}\) and \(s_{i,\phi}\) are neural networks with shared weights \(\phi\) which is performed by Masked Autoencoder for Density Estimation (MADE) [48]. The latter applies binary masks on weight matrices in such a way that each output \(i\) depends only on the inputs \(1:i\). Therefore the transformation from \(Z\) to \(X\) is as follows \[x_i=z_i\exp{s_{i,\phi}}+m_{i,\phi},\quad z_i\sim \mathcal{N}(0,1)\]

In the neural posterior estimation setting considered in this work, the variables transformed by the flow correspond to the epidemiological parameters \(\theta\), while the simulated observations \(\hat{X}\) are provided as conditioning information. These conditioning variables are supplied to the MADE network and influence the functions \(m_{i,\phi}\) and \(s_{i,\phi}\) that parameterize the autoregressive transformations, while not themselves being transformed by the flow. Consequently, the learned distribution is

\[p(\theta|\hat{X}) = \prod_{i=1}^{n} p(\theta_i|\theta_{1:i-1},\hat{X}).\]

2 SNPE-C Training Objective↩︎

In SBI framework the training loss consists of two components: an atomic loss, which performs contrastive training on proposal-sampled parameter atoms, and an additional maximum-likelihood (MLE) term. The MLE term here acts as a regularizer by encouraging the density estimator to assign high probability to prior-drawn samples, thereby reducing density leakage outside the bounded support of the prior. Note, that the function to learn here is \(p_X(\theta;\phi|x)\), corresponding to the posterior distribution of the parameters \(\theta\) given the simulated data \(x\). Atomic loss reads as follows \[\mathcal{L}_{atomic}=-\log\frac{p(x;\phi|\theta^k)}{\sum_{j=1}^{K}p_X(x;\phi|\theta^j)},\] where \(\{\theta^1,\dots,\theta^K\}\) are all sets of parameters \(\theta\) generated during training, while \(\theta^k\) is the one underlying the simulated trajectory \(x\). \[\mathcal{L}_{MLE}=-\log{p_{X}(\theta;\phi|x)}=-\log{p_{Z}(f^{-1}_{\phi}(\theta;x))}-\log{\det|\frac{\partial f^{-1}_{\phi}}{\partial \theta}|}\]

The final loss being \[\mathcal{L}=\mathcal{L}_{atomic}+\mathcal{L}_{MLE}\]

3pt

S5 Table.

Detailed posterior agreement, predictive performance, and runtime metrics for all SBI configurations. For each inference setting, simulation budget, batch size, and epidemic window, the table reports parameter-wise agreement between SBI and MCMC posterior distributions measured using the first Wasserstein distance, Kullback–Leibler (KL) divergence, and symmetric KL divergence. Metrics are provided separately for all inferred epidemiological parameters and contact-rate change-point parameters. In addition, posterior predictive performance (RMSE), runtime statistics, number of failed runs, and simulation budgets are reported. Values are summarized as mean and standard deviation across 16 repeated SBI runs.

References↩︎