arXiv CSopen access

Simulation-based inference for rapid Bayesian parameter estimation in epidemiological models: a comparison with MCMC

Simulation-based inference for rapid Bayesian parameter estimation in epidemiological models: a comparison with MCMC Alina Bazarova1,2,∗,† , Johann Fredrik Jadebeck3,6,† , Henrik Zunker4 , Carolina J. Klett-Tammen5 , Torben Heinsohn5 , Wolfgang Wiechert3,6 , Katharina Nöh3 , and Stefan Kesselheim1,2

arXiv:2606.27286v1 [cs.AI] 25 Jun 2026

1

Forschungszentrum Jülich, Jülich Supercomputing Centre, Jülich, Germany 2 Helmholtz AI, Germany 3 Forschungszentrum Jülich, Institute of Bio- and Geosciences, Jülich, Germany 4 German Aerospace Center, Institute of Software Technology, Department High-Performance Computing, Cologne, Germany 5 Helmholtz Centre for Infection Research, Braunschweig, Germany 6 RWTH Aachen University, Computational Systems Biology, Aachen, Germany ∗ †

Corresponding author: [email protected] These authors contributed equally to this work.

Abstract Mechanistic epidemiological models are widely used to support infectious disease forecasting and public-health decision making. Bayesian calibration of such models is commonly performed using Markov chain Monte Carlo (MCMC), which can become computationally expensive for high-dimensional nonlinear systems and repeated near-real-time analyses. Here, we investigate simulation-based inference (SBI) using neural posterior estimation as a scalable alternative for Bayesian calibration of a mechanistic SECIR epidemiological model using COVID-19 intensive care unit (ICU) occupancy data from Germany during 2020. We compared SBI and MCMC across multiple epidemic phases using both 31-day inference windows and a substantially more challenging 201-day reconstruction problem involving multiple transmission change points. Posterior agreement was evaluated quantitatively using Wasserstein distances and Kullback–Leibler divergences together with posterior predictive checks. Across the 31-day windows, SBI recovered posterior distributions in strong agreement with MCMC while accurately reproducing observed ICU trajectories. In the 201-day setting, SBI preserved the dominant posterior structure despite increased uncertainty. SBI, by combining CPU and GPU resources, substantially reduced computational runtime compared with MCMC, which was restricted to running on CPUs. Whereas MCMC required approximately 1000 seconds for the 31-day inference problems, SBI achieved comparable posterior and predictive performance in approximately 60–70 seconds on a single NVIDIA A100 GPU. For the 201-day inference problem, SBI required an average of 157 seconds, while the MCMC runs took over 19,000 seconds. Our results demonstrate that SBI provides a rapid and computationally efficient framework for Bayesian calibration of mechanistic epidemiological models, supporting repeated near-real-time inference and rapid outbreak analysis.

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. 1

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

2

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 transmissionparameter 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 (IN S ), infectious with symptoms (ISy ), infected severe (ISev ), infected critical (ICr ), 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 ICr compartment is directly linked to the observed ICU occupancy. 3

The ODE system is parameterized by transition times (T ), transition probabilities (µ), and transmission-related scaling factors (c, ρ, ξ). 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 ρ, and the relative infectiousness of non-symptomatic (ξIN S ) and symptomatic (ξISy ) 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 Table 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 c0 denote the baseline contact rate before the change in transmission conditions. If an intervention is introduced at time t1 > t0 , where t0 is the start of the simulation, the initial contact rate c0 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 δ ∈ (0, 1) using differentiable interpolation as described in [22]. The time-varying contact rate c(t) is thus defined as:   t ≤ t1 c0 , c(t) := b c(t), t ∈ (t1 , t1 + δ)   (1 − r)c0 , t ≥ t1 + δ

(1)

where ĉ(t) denotes the smooth interpolation between c0 and (1 − r)c0 . 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]. Table 1. Parameters of the SECIR compartmental model Name Symbol Initial contact rate c0 Reduction factor r Transmission probability on contact ρ Relative infectiousness non-symptomatic individuals ξIN S Risk infectiousness symptomatic individuals ξISy Population size N Time exposed TE Time infected non-symptomatic TIN S Time infected symptomatic TISy Time infected severe TISev Time infected critical TICr I Transition prob. from IN S to ISy µISy NS Transition prob. from ISy to ISev µIISev Sy

Dimension R≥0 [0, 1] [0, 1] [0, 1] [0, 1] R>0 R>0 R>0 R>0 R>0 R>0 [0, 1] [0, 1]

Transition prob. from ISev to ICr Transition prob. from ICr to D

[0, 1] [0, 1]

4

µIICr Sev µD ICr

ISy

Exposed, Not Infectious E

ϕρ

ξI

1 TE

Infectious, No Symptoms INS ISy

I +ξI I NS NS Sy Sy

1−µI

N −D

TI

µI TI

Infectious, Symptoms ISy I

NS

1−µISev Sy

TI

NS

Susceptible S

NS

NS

Sy

I

µISev

Recovered R

Sy

TI

Sy

I

1−µICr

1−µD I TI

Dead D

Cr

Infected, Critical ICr

µD I TI

TI

Cr

Cr

I

µICr TI

Cr

Sev

Sev

Sev

Infected, Severe ISev

Sev

Figure 1. Overview of the SECIR model compartment and transition structure. The model partitions the population into susceptible (S), exposed but not yet infectious (E), infectious with no symptoms (IN S ), infectious with symptoms (ISy ), infected severe (ISev ), infected critical (ICr ), recovered (R), and dead (D) compartments. Arrows indicate transitions between compartments, and edge labels denote the corresponding transition rates. Disease progression is driven by mean transition times T and transition probabilities µ, while infection depends on the effective contact rate, the transmission probability on contact, and the relative infectiousness of non-symptomatic and symptomatic individuals.

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 obs obs making. Let xobs 1:231 = (x1 , . . . , x231 ) denote the observed ICU occupancy time series over this period. For our inference experiments below, we consider different time windows xobs t0 :tend , where t0 is the respective starting time and tend 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 TE , TIN S , TISy , TISev , and TICr as well as the relative infectiousness of symptomatic individuals ξISy . In addition, the Bayesian model infers changes in the time-varying contact rate on a grid. Let 5

∆ be the minimum number of days for which a contact-rate level is valid, then for a time series  0 +1 starting at day t0 and ending after day tend the number of intervals on the grid is K = tend −t . ∆ For τk = t0 + k∆, k = 0, . . . , K the unsmoothed effective time-varying contact rate βpc (t) is given by the piecewise-constant function βpc (t) := ρ · c0 · (1 − rk ) for t ∈ [τk−1 , τk ),

k = 1, . . . , K,

(2)

where the rk 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 βpc (t) over a short transition as described above (Eq. (1)). As before, larger values for rk 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 c0 = 1 and ρ = 1, meaning that for each interval on the grid the effective contact rate is controlled by a single parameter rk . The full parameter vector of the Bayesian model is given by θ. θ = (TE , TIN S , TISy , TISev , TICr , ξISy , r1 , . . . , rK )

(3)

Model Initialization To initialize the epidemiological compartments for inference starting at time t0 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 Table 2, while the priors on the reduction factors for 201-day window are summarized in Prior bounds for contact-reduction (damping) parameters in the 201-day inference problem. Lower and upper bounds of the uniform prior distributions assigned to the 15 time-varying contact-reduction parameters (r0 , . . . , r14 ). Specification A corresponds to the initial prior configuration, while Specification B denotes the revised prior used after prior predictive checking to improve the plausibility of simulated epidemic trajectories and reduce prior support for unrealistic transmission dynamics, as show in Fig. S1, 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

6

Table 2. Priors for the Bayesian model. The median and 95% intervals cover reasonable values for the SECIR model. Prior Description TE ∼ Uniform(2. 67, 4. 0) median 3.36; 95% interval [2.7, 3.97] TIN S ∼ LogNormal(ln(10), 0.2) median 10; 95% interval [6.76, 14.80] TISy ∼ LogNormal(ln(10), 0.2) median 10; 95% interval [6.76, 14.80] TISev ∼ LogNormal(ln(10), 0.2) median 10; 95% interval [6.76, 14.80] TICr ∼ LogNormal(ln(10), 0.2) median 10; 95% interval [6.76, 14.80] ξISy ∼ Uniform(0.01, 0.9) median 0.46; 95% interval [0.03, 0.88] ind

rk ∼ Uniform(0, 1)

median 0.5; 95% interval [0.025, 0.975]

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 ν = 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 t0 to T , the likelihood was p(xobs t0 :T |θ, t0:T ) =

T Y

p(xobs i |θ, ti ) =

i=t0

T Y

sim Student-tν (xobs (θ, ti ), σi ) i |x

(4)

i=t0

where xsim p (θ, ti ) denotes the simulated ICU occupancy at day ti and the observation noise modeled as σi = |xsim (θ, ti )|. 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 7

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 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(θ | X) ∝ L(X | θ)π(θ),

(5)

where L(·) denotes the likelihood function and π(·) 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 θj ∼ π(θ) and propagating them through the SECIR simulator to obtain synthetic ICU occupancy trajectories X̂j = f (θj ), where f denotes the epidemiological simulator. This produced a set of simulated parameter–data pairs (θj , X̂j ) used to train the neural density estimator. The network learns a conditional posterior approximation p̃(θ | 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 p̃(θ | 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

8

Figure 2. Workflow of amortized neural posterior estimation (NPE). Parameters θ are sampled from the prior distribution and propagated through the SECIR epidemiological model to generate simulated ICU occupancy trajectories X̂. The resulting parameter–trajectory pairs (θ, X̂) are used to train a neural density estimator that approximates the conditional posterior distribution p(θ|X). After training, the observed ICU occupancy data X are provided as conditioning information, and samples are drawn from the learned posterior approximation p̂(θ|X). Because the posterior estimator is trained once and subsequently reused for inference, the approach is amortized. as p(θ | X̂) =

d Y

p(θi | θ1 , X̂),

(6)

i=1

where 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, Mean Autoregressive Flows. 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 X̂ to a learned feature representation h(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, Latomic + LMLE , (7) 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 (SNPE-C Training Objective). 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, 9

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 20,000, 50,000, and 100,000 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 θj , let p(θj ) denote the marginal posterior density inferred by SBI and q(θj ) the corresponding marginal posterior density inferred by MCMC. Since both posteriors are represented by samples, continuous density estimates p̂(θj ) and q̂(θ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. Wasserstein distance. The first-order Wasserstein distance (Earth Mover’s Distance) between the two marginal posterior distributions was computed as W1 (p, q) =

inf

γ∈Π(p,q)

E(x,y)∼γ [|x − y|] ,

(8)

where Π(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. Kullback–Leibler divergence. posterior was defined as

The KL divergence from the SBI posterior to the MCMC Z

DKL (p(θj )||q(θj )) =

p̂(θj ) log

p̂(θj ) dθj . q̂(θj )

(9)

The integral was approximated numerically on the discretized grid, DKL (p(θj )||q(θj )) ≈

K X k=1

p̂(θj,k ) log

p̂(θj,k ) ∆θ, q̂(θj,k )

(10)

where ∆θ denotes the grid spacing. To improve numerical stability, density values were clipped to a small positive constant before evaluating the logarithm.

10

Symmetric KL divergence. Because the KL divergence is asymmetric, we additionally computed the symmetric KL divergence, 1 [DKL (p(θj )∥q(θj )) + DKL (q(θj )∥p(θj ))] . (11) 2 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. Dsym (p(θj ), q(θj )) =

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 50,000 simulations and a batch size of 14,400. For the 201-day inference problem, which involved a higher-dimensional parameter space and multiple change points, the selected configuration used 100,000 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. Prior predictive simulations for the 200-day inference setting under two alternative change-point prior specifications. Each panel shows 10,000 ICU occupancy trajectories (vertical axis) over 200 days (horizontal axis) generated from parameters sampled from the prior distribution and propagated through the forward simulator. The observed trajectory (black curve) is overlaid for reference. Broader prior bounds (left, Spec A in Table S1) result in highly dispersed and partially implausible epidemic dynamics, whereas tighter bounds (right, Spec B in Table S1) concentrate prior mass around epidemiologically plausible trajectories while retaining sufficient variability, 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. 11

Under the final prior specification, the observed ICU occupancy trajectory was well covered by the prior predictive envelope (Fig. Prior predictive simulations for the 200-day inference setting under two alternative change-point prior specifications. Each panel shows 10,000 ICU occupancy trajectories (vertical axis) over 200 days (horizontal axis) generated from parameters sampled from the prior distribution and propagated through the forward simulator. The observed trajectory (black curve) is overlaid for reference. Broader prior bounds (left, Spec A in Table S1) result in highly dispersed and partially implausible epidemic dynamics, whereas tighter bounds (right, Spec B in Table S1) concentrate prior mass around epidemiologically plausible trajectories while retaining sufficient variability, 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. offset 0, 31 days 900 14400 45000 offset 160, 31 days 900 14400 45000 offset 200, 31 days 900 14400 45000 offset 0, 201 days 450 1800 22500

Wasserstein

KL

KL sym

RMSE

Runtime

0.98 (1.14) 0.75 (0.51) 0.63 (0.37)

17.18 (42.13) 17.92 (24.98) 17.3 (24.38)

9.35 (21.37) 9.58 (13.2) 8.92 (12.36)

153.52 (384.73) 66.41 (26.16) 70.43 (62.9)

175.19 (102.72) 71.47 (19.22) 63.13 (13.32)

3.38 (2.5) 2.37 (1.43) 2.47 (1.34)

189.49 (69.15) 134.67 (73.38) 115.07 (23.34) 138.17 (53.13) 75.79 (31.35) 85.88 (54.78) 146.1 (54.44) 78.8 (31.63) 137.86 (103.24)

166.1 (27.65) 60.83 (3.09) 61.01 (6.67)

1.90 (0.49) 1.89 (0.45) 1.98 (0.34)

87.6 (5.59) 85.44 (12.33) 84.6 (13.6)

58.46 (18.64) 46.44 (9.58) 47.06 (10.82)

68.72 (37.52) 98.79 (66.22) 144.51 (140.1)

161.11 (27.65) 60.07 (3.97) 59.68 (5.02)

0.50 (0.28) 0.35 (0.35) 0.46 (0.18)

61.8 (56.6) 141.69 (47.11) 132.25 (44.12)

39.01 (42.59) 72.14 (24.67) 66.98 (22.65)

323.43 (126.72) 874.51 (1000.34) 325.72 (135.78) 157.08 (42.28) 397.93 (129.21) 86.55 (9.23)

Table 3. Comparison of posterior agreement, posterior predictive performance, and runtime across inference configurations. For each inference window and batch size, we report the mean and standard deviation across 16 repeated runs of the 1-Wasserstein distance, Kullback-Leibler divergence, symmetric Kullback-Leibler divergence between SBI and MCMC posterior distributions, posterior predictive RMSE, and total SBI runtime in seconds. The selected configurations used for the main analyses were 50,000 simulations with batch size 14,400 for the 31-day inference windows and 100,000 simulations with batch size 1,800 for the 201-day inference window. All SBI experiments were performed on a single NVIDIA A100 GPU. 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 ± 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 ± 54.78 and 98.79 ± 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

12

A) 31 days, offset 0

B) 31 days, offset 160

ICU occupancy

2500

2000

2000

1500

RMSE SBI = 33.4 1500 RMSE MCMC = 31.8

1000

1000

500 0

5

10

15

20

25

30

0

C) 31 days, offset 200

ICU occupancy

5

10

15

20

25

30

D) 201 days, offset 0

5000

4000

4500

3000

4000 RMSE SBI = 34.7 RMSE MCMC = 18.2

2000

3500

1000

3000

RMSE SBI = 48.7 RMSE MCMC = 14.4

0

5

10

MCMC 95% credible interval

15 Days

20

25

30

0

SBI 95% credible interval

RMSE SBI = 150.2 RMSE MCMC = 23.1

0

25

50

75 100 125 150 175 200 Days

MCMC mean

SBI mean

Data

Figure 3. Posterior predictive checks for the inferred SBI posterior distributions across the analyzed inference windows. For each setting, ICU occupancy trajectories were simulated using parameter samples drawn from the inferred posterior distribution. The solid line shows the posterior predictive mean trajectory, the shaded region indicates the 95% pointwise credible interval, and the black line corresponds to the observed ICU occupancy data. The panels include the three 31-day inference windows starting at offsets 0 (upper left), 160 (upper right), and 200 (lower left) days after April 24, 2020, as well as the extended 201-day (lower right) inference window. Across all settings, the posterior predictive distributions capture the dominant temporal dynamics of the observed epidemic trajectory.

13

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 ± 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, Comparison of marginal posterior distributions inferred using simulation-based inference (SBI, blue) and Markov chain Monte Carlo (MCMC, orange) for the extended 200-day inference window starting on April 24, 2020. The panels show posterior distributions for epidemiological parameters and all inferred change-point parameters. Despite the substantially increased dimensionality of the inference problem, SBI recovered posterior modes and overall parameter trends broadly consistent with the corresponding MCMC posteriors, while generally producing smoother and broader marginal distributions 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±0.51 and a symmetric KL divergence of 9.58±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 ± 1.43 and a symmetric KL divergence of 75.79 ± 31.35. For offset 200, the corresponding values were 1.89 ± 0.45 and 46.44 ± 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, Comparison of marginal posterior distributions inferred using simulation-based inference (SBI, blue) and Markov chain Monte Carlo (MCMC, orange) for the extended 200-day inference window starting on April 24, 2020. The panels show posterior distributions for epidemiological parameters and all inferred change-point parameters. Despite the substantially increased dimensionality of the inference problem, SBI recovered posterior modes and overall parameter trends broadly consistent with the corresponding MCMC posteriors, while generally producing smoother and broader marginal distributions. Nevertheless, SBI preserved the main posterior structure observed in the MCMC 14

reference. For the selected 201-day configuration, the Wasserstein distance was 0.35 ± 0.35 and the symmetric KL divergence was 72.14 ± 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, S5 Table). For the 31-day inference windows, configurations trained with 20,000 simulations generally showed higher RMSE and greater variability across repeated runs. Increasing the simulation budget to 100,000 did not consistently improve posterior predictive performance or agreement with the MCMC reference posteriors. Across the three 31-day windows, 50,000 simulations provided stable posterior recovery while maintaining substantially lower computational cost. Among the tested batch sizes, a batch size of 14,400 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 100,000 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 31day 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

15

Figure 4. Comparison of marginal posterior distributions inferred using simulation-based inference (SBI, blue) and Markov chain Monte Carlo (MCMC, orange). The panels show posterior distributions for epidemiological parameters and change-point parameters inferred from ICU occupancy data for the 31-day windows starting at offsets 0 (upper panel), 160 (middle panel), and 200 (lower panel) days after April 24, 2020. Across all settings, SBI recovered posterior modes and overall parameter trends broadly consistent with the MCMC reference posteriors.

16

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 shortand long-window analyses.

Runtime comparison Details for all MCMC runs are reported in MCMC runtime details and diagnostics for all runs. For each inference setting, we report the number of chains and samples, acceptance rate, runtime, and maximum R̂ value. The reported maximum R̂ and threshold counts refer to inferred model parameters. The 30 day runs converged nicely, with just the 31 day, offset 160 run having a single parameter with R̂ = 1.012 > 1.01, which can effectively be interpreted as convergence. Conversely, for the 200 day run no parameters converged despite the long runtime and large sample budget. As shown in Table S3, the R̂ values for the posterior predictive values are lower, between 1.009 and 1.250, indicating that the 16 chains mostly agree on the posterior predictive even though they do not agree on the parameters. To keep computations light, we thinned the 201 days MCMC run by a factor of 100, significantly reducing memory footprint and all R̂ values are summarized in Table Rank-normalized R̂ values for all parameters and posterior predictive ICU occupancy. The predictive quantities are reported as derived diagnostics of chain agreement for the model outputs. For the first day of simulation (day 0, day 160, day 200), each of the runs used the initialization process described in materials and methods. The 31 day inference windows show good convergence. For the 200 day run, mixing was harder for the parameters than for the posterior predictive values and despite the long runtime, MCMC did not achieve convergence. 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 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 ± 19.22 seconds for the offset 0 window, 60.83 ± 3.09 seconds for the offset 160 window, and 60.07 ± 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 19,102 seconds, whereas the selected SBI configuration required 157.08 ± 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 R̂ (maximum R̂ = 1.605) for the parameters. Overall, SBI reduced inference time by approximately one order of magnitude for the 31day 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.

17

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 18

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 19

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.

20

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.

References [1] Lingcai Kong, Mengwei Duan, Jin Shi, Jie Hong, Zhaorui Chang, and Zhijie Zhang. Compartmental structures used in modeling covid-19: a scoping review. Infectious Diseases of Poverty, 11(1):72, Jun 2022. ISSN 2049-9957. doi: 10.1186/s40249-022-01001-y. [2] Eduard Campillo-Funollet, James Van Yperen, Phil Allman, Michael Bell, Warren Beresford, Jacqueline Clay, Matthew Dorey, Graham Evans, Kate Gilchrist, Anjum Memon, Gurprit Pannu, Ryan Walkley, Mark Watson, and Anotida Madzvamuse. Predicting and forecasting the impact of local outbreaks of covid-19: use of seir-d quantitative epidemiological modelling for healthcare demand and capacity. International Journal of Epidemiology, 50(4):1103–1113, 07 2021. ISSN 0300-5771. doi: 10.1093/ije/dyab106. [3] Sen Pei, Sasikiran Kandula, and Jeffrey Shaman. Differential effects of intervention timing on covid-19 spread in the united states. Science Advances, 6(49):eabd6370, 2020. doi: 10.1126/ sciadv.abd6370. [4] Jonas Gilg, Johann F. Jadebeck, Mariama Jaiteh, David Kerkmann, Niklas Medinger, Shabaz Memon, Anna Clara Wendler, Moritz Zeumer, Henrik Zunker, Maximilian Franz Betz, Ralf Hannemann-Tamas, Jonas Heinicke, Julian Litz, Achim Basermann, Cas Cremers, Manuel Dahmen, Andreas Gerndt, Jens Henrik Göbbert, Björn Hagemeier, Carolina J. Klett-Tammen, Berit Lange, Katharina Nöh, Sarah Strassburger, Michael Meyer-Herrmann, and Martin Joachim Kühn. A full software stack for epidemic disease management: Unlocking the joint potential of software technology and supercomputing. Technical report, arxiv, 2026. URL https: //elib.dlr.de/224465/. [5] Jonas Dehning, Johannes Zierenberg, F. Paul Spitzner, Michael Wibral, Joao Pinheiro Neto, Michael Wilczek, and Viola Priesemann. Inferring change points in the spread of covid-19 reveals the effectiveness of interventions. Science, 369(6500):eabb9789, 2020. doi: 10.1126/ science.abb9789. [6] Seth Flaxman, Swapnil Mishra, Axel Gandy, H. Juliette T. Unwin, Thomas A. Mellan, Helen Coupland, Charles Whittaker, Harrison Zhu, Tresnia Berah, Jeffrey W. Eaton, Mélodie Monod, Pablo N. Perez-Guzman, Nora Schmit, Lucia Cilloni, Kylie E. C. Ainslie, Marc Baguelin, Adhiratha Boonyasiri, Olivia Boyd, Lorenzo Cattarino, Laura V. Cooper, Zulma Cucunubá, Gina Cuomo-Dannenburg, Amy Dighe, Bimandra Djaafara, Ilaria Dorigatti, Sabine L. van Elsland, Richard G. FitzJohn, Katy A. M. Gaythorpe, Lily Geidelberg, Nicholas C. Grassly, William D. Green, Timothy Hallett, Arran Hamlet, Wes Hinsley, Ben Jeffrey, Edward Knock, Daniel J. Laydon, Gemma Nedjati-Gilani, Pierre Nouvellet, Kris V. Parag, Igor Siveroni, Hayley A. Thompson, Robert Verity, Erik Volz, Caroline E. Walters, Haowei Wang, Yuanrong Wang, Oliver J. Watson, Peter Winskill, Xiaoyue Xi, Patrick G. T. Walker, Azra C. Ghani, Christl A. Donnelly, Steven Riley, Michaela A. C. Vollmer, Neil M. Ferguson, Lucy C. Okell,

21

Samir Bhatt, and Imperial College COVID-19 Response Team. Estimating the effects of non-pharmaceutical interventions on covid-19 in europe. Nature, 584(7820):257–261, Aug 2020. ISSN 1476-4687. doi: 10.1038/s41586-020-2405-7. [7] Sebastian Funk, Anton Camacho, Adam J. Kucharski, Rosalind M. Eggo, and W. John Edmunds. Real-time forecasting of infectious disease dynamics with a stochastic semi-mechanistic model. Epidemics, 22:56–61, 2018. ISSN 1755-4365. doi: 10.1016/j.epidem.2016.11.003. The RAPIDD Ebola Forecasting Challenge. [8] Cécile Viboud, Kaiyuan Sun, Robert Gaffey, Marco Ajelli, Laura Fumanelli, Stefano Merler, Qian Zhang, Gerardo Chowell, Lone Simonsen, and Alessandro Vespignani. The rapidd ebola forecasting challenge: Synthesis and lessons learnt. Epidemics, 22:13–21, 2018. ISSN 1755-4365. doi: 10.1016/j.epidem.2017.08.002. The RAPIDD Ebola Forecasting Challenge. [9] Nicholas G. Reich, Logan C. Brooks, Spencer J. Fox, Sasikiran Kandula, Craig J. McGowan, Evan Moore, Dave Osthus, Evan L. Ray, Abhinav Tushar, Teresa K. Yamana, Matthew Biggerstaff, Michael A. Johansson, Roni Rosenfeld, and Jeffrey Shaman. A collaborative multiyear, multimodel assessment of seasonal influenza forecasting in the united states. Proceedings of the National Academy of Sciences, 116(8):3146–3154, 2019. doi: 10.1073/pnas.1812594116. [10] Radu V. Craiu and Jeffrey S. Rosenthal. Bayesian computation via markov chain monte carlo. Annual Review of Statistics and Its Application, 1(Volume 1, 2014):179–201, 2014. ISSN 2326-831X. doi: 10.1146/annurev-statistics-022513-115540. [11] Nicholas P. Jewell, Theodoros Kypraios, Paul Neal, and Gareth O. Roberts. Bayesian analysis for emerging infectious diseases. Journal of the Royal Statistical Society: Series C (Applied Statistics), 58(2):317–336, 2009. [12] Michael Li, Jonathan Dushoff, and Benjamin M Bolker. Fitting mechanistic epidemic models to data: A comparison of simple markov chain monte carlo approaches. Stat Methods Med Res, 27(7):1956–1967, July 2018. [13] Thomas House, Ashley Ford, Shiwei Lan, Samuel Bilson, Elizabeth Buckingham-Jeffery, and Mark Girolami. Bayesian uncertainty quantification for transmissibility of influenza, norovirus and ebola using information geometry. J R Soc Interface, 13(121), August 2016. [14] Radu V. Craiu and Jeffrey S. Rosenthal. Bayesian computation via markov chain monte carlo. Annual Review of Statistics and Its Application, 1:179–201, 2014. doi: 10.1146/ annurev-statistics-022513-115540. [15] Lorenzo Contento, Noemi Castelletti, Elba Raimúndez, Ronan Le Gleut, Yannik Schälte, Paul Stapor, Ludwig Christian Hinske, Michael Hoelscher, Andreas Wieser, Katja Radon, Christiane Fuchs, Jan Hasenauer, and the KoCo19 study group. Integrative modelling of reported case numbers and seroprevalence reveals time-dependent test efficiency and infectious contacts. Epidemics, 43:100681, 2023. doi: 10.1016/j.epidem.2023.100681. [16] Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceedings of the National Academy of Sciences, 117(48):30055–30062, 2020. doi: 10.1073/ pnas.1912789117. [17] Jan-Matthis Lueckmann, Jan Boelts, David S. Greenberg, Pedro J. Gonçalves, and Jakob H. Macke. Benchmarking simulation-based inference, 2021. 22

[18] Maxwell H Wang and Jukka-Pekka Onnela. Flexible bayesian inference on partially observed epidemics. J Complex Netw, 12(2):cnae017, March 2024. [19] Francesco Pinotti, Julien Thézé, Xavier Bailly, and Guillaume Fournié. Simulation basedinference of epidemiological and phylodynamic models via neural posterior estimation. bioRxiv, 2025. doi: 10.1101/2025.11.25.690436. [20] Vincent Wieland, Nils Wassmuth, Lorenzo Contento, Martin Kühn, and Jan Hasenauer. Assessment of simulation-based inference methods for stochastic compartmental models in epidemiological research, 2026. [21] J.P. Manzano-Patrón, Michael Deistler, Cornelius Schröder, Theodore Kypraios, Pedro J. Gonçalves, Jakob H. Macke, and Stamatios N. Sotiropoulos. Uncertainty mapping and probabilistic tractography using simulation-based inference in diffusion mri: A comparison with classical bayes. Medical Image Analysis, 103:103580, 2025. ISSN 1361-8415. doi: 10.1016/j.media.2025.103580. [22] Martin J. Kühn, Daniel Abele, Tanmay Mitra, Wadim Koslow, Majid Abedi, Kathrin Rack, Martin Siggel, Sahamoddin Khailaie, Margrit Klitz, Sebastian Binder, Luca Spataro, Jonas Gilg, Jan Kleinert, Matthias Häberle, Lena Plötzke, Christoph D. Spinner, Melanie Stecher, Xiao Xiang Zhu, Achim Basermann, and Michael Meyer-Hermann. Assessment of effective mitigation and prediction of the spread of SARS-CoV-2 in Germany using demographic information and spatial resolution. Mathematical Biosciences, page 108648, 2021. ISSN 0025-5564. doi: https://doi.org/10.1016/j.mbs.2021.108648. [23] Julia Bicker, Carlotta Gerstein, David Kerkmann, Sascha Korf, René Schmieding, Anna Wendler, Henrik Zunker, Daniel Abele, Maximilian Betz, Khoa Nguyen, Lena Plötzke, Kilian Volmer, Agatha Schmidt, Nils Waßmuth, Patrick Lenz, Daniel Richter, Hannah Tritzschak, Ralf Hannemann-Tamas, Julian Litz, Paul Johannssen, Marielena Borges, Annika Jungklaus, Manuel Heger, Annalena Lange, Elisabeth Kluth, Kathrin Rack, Vincent Wieland, Jonas Arruda, Sebastian Binder, Margrit Klitz, Martin Siggel, Manuel Dahmen, Achim Basermann, Michael Meyer-Hermann, Jan Hasenauer, and Martin J. Kühn. Memilio – a high performance modular epidemics simulation software for multi-scale and comparative simulations of infectious disease dynamics, 2026. [24] Robert Koch-Institut. Intensivkapazitäten und covid-19-intensivbettenbelegung in deutschland, November 2025. [25] Hendrik Streeck, Bianca Schulte, Beate M. Kümmerer, Enrico Richter, Tobias Höller, Christine Fuhrmann, Eva Bartok, Ramona Dolscheid-Pommerich, Moritz Berger, Lukas Wessendorf, Monika Eschbach-Bludau, Angelika Kellings, Astrid Schwaiger, Martin Coenen, Per Hoffmann, Birgit Stoffel-Wagner, Markus M. Nöthen, Anna M. Eis-Hübinger, Martin Exner, Ricarda Maria Schmithausen, Matthias Schmid, and Gunther Hartmann. Infection fatality rate of sars-cov2 in a super-spreading event in germany. Nature Communications, 11(1), November 2020. ISSN 2041-1723. doi: 10.1038/s41467-020-19509-y. [26] Daniela Gornyk, Manuela Harries, Stephan Glöckner, Monika Strengert, Tobias Kerrinnes, JanaKristin Heise, Henrike Maaß, Julia Ortmann, Barbora Kessel, Yvonne Kemmling, Berit Lange, and Gérard Krause. Sars-cov-2 seroprevalence in germany. Deutsches Ärzteblatt international, December 2021. ISSN 1866-0452. doi: 10.3238/arztebl.m2021.0364.

23

[27] Henrik Zunker, René Schmieding, David Kerkmann, Alain Schengen, Sophie Diexer, Rafael Mikolajczyk, Michael Meyer-Hermann, and Martin J. Kühn. Novel travel time aware metapopulation models and multi-layer waning immunity for late-phase epidemic and endemic scenarios. PLOS Computational Biology, 20(12), Dez 2024. doi: 10.1371/journal.pcbi.1012630. [28] Robert Koch-Institut. Sars-cov-2 infektionen in deutschland, 2025. [29] Regionaldatenbank Deutschland. Fortschreibung des Bevölkerungsstandes: 12411-02-03-4 Bevölkerung nach Geschlecht und Altersgruppen (17) - Stichtag 31.12. - regionale Tiefe: Kreise und krfr. Städte, 2022. URL https://www.regionalstatistik.de/genesis/online? operation=statistic&levelindex=0&levelid=1646144362683&code=12411#abreadcrumb. [30] Julia Schilling, Ann-Sophie Lehfeld, Dirk Schumacher, Alexander Ullrich, Michaela Diercke, Silke Buda, Walter Haas, and RKI COVID-19 Study Group. Krankheitsschwere der ersten COVID-19-welle in deutschland basierend auf den meldungen gemäß infektionsschutzgesetz. Journal of Health Monitoring, 5(S11):2–20, 2020. doi: 10.25646/7169. [31] Oyungerel Byambasuren, Magnolia Cardona, Katy Bell, Justin Clark, Mary-Louise McLaws, and Paul Glasziou. Estimating the extent of asymptomatic COVID-19 and its potential for community transmission: Systematic review and meta-analysis. Journal of the Association of Medical Microbiology and Infectious Disease Canada, 5(4):223–234, 2020. doi: 10.3138/ jammi-2020-0030. [32] Wafa Dhouib, Jihen Maatoug, Imen Ayouni, Nawel Zammit, Rim Ghammem, Sihem Ben Fredj, and Hassen Ghannem. The incubation period during the pandemic of COVID19: A systematic review and meta-analysis. Systematic Reviews, 10(1):101, 2021. doi: 10.1186/s13643-021-01648-y. [33] Jonas Dehning, Sebastian B. Mohr, Sebastian Contreras, Philipp Dönges, Emil N. Iftekhar, Oliver Schulz, Philip Bechtle, and Viola Priesemann. Impact of the Euro 2020 championship on the spread of COVID-19. Nature Communications, 14(1):122, 2023. doi: 10.1038/s41467-022-35512-x. [34] C. J. F. ter Braak and J. A. Vrugt. Differential evolution markov chain with snooker updater and fewer chains. Statistics and Computing, 18(4):435–446, 2008. doi: 10.1007/s11222-008-9104-9. [35] Oriol Abril-Pla, Virgile Andreani, Colin Carroll, Larry Dong, Christopher J. Fonnesbeck, Maxim Kochurov, Ravin Kumar, Junpeng Lao, Christian C. Luhmann, Osvaldo A. Martin, Michael Osthege, Ricardo Vieira, Thomas Wiecki, and Robert Zinkov. PyMC: A modern and comprehensive probabilistic programming framework in Python. PeerJ Computer Science, 9 (e1516), 2023. doi: 10.7717/peerj-cs.1516. [36] Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. b for assessing convergence of Rank-normalization, folding, and localization: An improved R MCMC (with discussion). Bayesian Analysis, 16(2):667–718, 2021. doi: 10.1214/20-BA1221. [37] Jan-Matthis Lueckmann, Giacomo Bassetto, Theofanis Karaletsos, and Jakob H. Macke. Likelihood-free inference with emulator networks, 2019. [38] George Papamakarios, Theo Pavlakou, and Iain Murray. Masked autoregressive flow for density estimation, 2018.

24

[39] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998. doi: 10.1109/5.726791. [40] David S. Greenberg, Marcel Nonnenmacher, and Jakob H. Macke. Automatic posterior transformation for likelihood-free inference, 2019. [41] Alvaro Tejero-Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pedro J. Gonçalves, David S. Greenberg, and Jakob H. Macke. sbi: A toolkit for simulationbased inference. Journal of Open Source Software, 5(52):2505, 2020. doi: 10.21105/joss.02505. [42] Sarah A. Vollert, Christopher Drovandi, Cailan Jeynes-Smith, Luz V. Pascal, and Matthew P. Adams. Beyond data: Leveraging non-empirical information and expert knowledge in bayesian model calibration, 2025. [43] Henrik Zunker, Philipp Dönges, Patrick Lenz, Seba Contreras, and Martin J. Kühn. Riskmediated dynamic regulation of effective contacts de-synchronizes outbreaks in metapopulation epidemic models. Chaos, Solitons & Fractals, 199:116782, 2025. ISSN 0960-0779. doi: 10.1016/ j.chaos.2025.116782. [44] George Papamakarios and Iain Murray. Fast ϵ-free inference of simulation models with bayesian conditional density estimation, 2018. [45] Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang. On large-batch training for deep learning: Generalization gap and sharp minima. In International Conference on Learning Representations, 2017. URL https:// openreview.net/forum?id=H1oyRlYgg. [46] Dominic Masters and Carlo Luschi. Revisiting small batch training for deep neural networks, 2018. URL https://arxiv.org/abs/1804.07612. [47] Jülich Supercomputing Centre. JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre. Journal of large-scale research facilities, 7(A182), 2021. doi: 10.17815/jlsrf-7-182. URL http://dx.doi.org/10. 17815/jlsrf-7-182. [48] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle. Made: Masked autoencoder for distribution estimation, 2015. arXiv:1502.03509.

25

Supplementary information S1

Mean Autoregressive Flows

Generally, normalizing flows learn the distribution of the data X pX (x; ϕ) through a mapping fϕ : Rn → Rn between the latent variables Z and the data X, where ϕ are the parameters of the mapping. Then, the change of variables formula is used pX (x; ϕ) = pZ (fϕ−1 (x)) det |

∂fϕ−1 ∂x

|,

(12)

∂f −1

ϕ where det | ∂x | is a Jacobian of the transformation between Z and X. Typically, fϕ is a chain of the transformations applied to a simple probability density function pZ , such as the one of the normal distribution. The network is then train via minimizing the negative log-likelihood

L = − log pX (x) = − log pZ (fϕ−1 (x)) − log det |

∂fϕ−1 ∂x

|

(13)

QnMAFs model the probability density function of the data in an autoregressive way pX (x) = i=1 p(xi |x1:i−1 ), where x = (x1 , . . . , xn ) and x1:i = (x1 , . . . , xi ). p(xi |x1:i ) ∼ N (xi |mi,ϕ , exp si,ϕ ), mi,ϕ = mi,ϕ (x1:i−1 ), si,ϕ = si,ϕ (x1:i−1 )

(14)

where mi,ϕ and si,ϕ are neural networks with shared weights ϕ 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 xi = zi exp si,ϕ + mi,ϕ ,

zi ∼ N (0, 1)

(15)

In the neural posterior estimation setting considered in this work, the variables transformed by the flow correspond to the epidemiological parameters θ, while the simulated observations X̂ are provided as conditioning information. These conditioning variables are supplied to the MADE network and influence the functions mi,ϕ and si,ϕ that parameterize the autoregressive transformations, while not themselves being transformed by the flow. Consequently, the learned distribution is p(θ|X̂) =

n Y

p(θi |θ1:i−1 , X̂).

(16)

i=1

S2

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 pX (θ; ϕ|x), corresponding to the posterior distribution of the parameters θ given the simulated data x. Atomic loss reads as follows p(x; ϕ|θk ) Latomic = − log PK , j j=1 pX (x; ϕ|θ ) 26

(17)

where {θ1 , . . . , θK } are all sets of parameters θ generated during training, while θk is the one underlying the simulated trajectory x. LM LE = − log pX (θ; ϕ|x) = − log pZ (fϕ−1 (θ; x)) − log det |

∂fϕ−1 ∂θ

|

(18)

The final loss being L = Latomic + LM LE

27

(19)

Figure S1. Prior predictive simulations for the 200-day inference setting under two alternative change-point prior specifications. Each panel shows 10,000 ICU occupancy trajectories (vertical axis) over 200 days (horizontal axis) generated from parameters sampled from the prior distribution and propagated through the forward simulator. The observed trajectory (black curve) is overlaid for reference. Broader prior bounds (left, Spec A in Table S1) result in highly dispersed and partially implausible epidemic dynamics, whereas tighter bounds (right, Spec B in Table S1) concentrate prior mass around epidemiologically plausible trajectories while retaining sufficient variability

28

Figure S2. Comparison of marginal posterior distributions inferred using simulation-based inference (SBI, blue) and Markov chain Monte Carlo (MCMC, orange) for the extended 200-day inference window starting on April 24, 2020. The panels show posterior distributions for epidemiological parameters and all inferred change-point parameters. Despite the substantially increased dimensionality of the inference problem, SBI recovered posterior modes and overall parameter trends broadly consistent with the corresponding MCMC posteriors, while generally producing smoother and broader marginal distributions.

29

Table S1. Prior bounds for contact-reduction (damping) parameters in the 201-day inference problem. Lower and upper bounds of the uniform prior distributions assigned to the 15 time-varying contact-reduction parameters (r0 , . . . , r14 ). Specification A corresponds to the initial prior configuration, while Specification B denotes the revised prior used after prior predictive checking to improve the plausibility of simulated epidemic trajectories and reduce prior support for unrealistic transmission dynamics, as show in Fig. S1. reduction factor 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14

Lower (Spec A) 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.8 0.8 0.8 0.8 0.7 0.7 0.7

Upper (Spec A) 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

Lower (Spec B) 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.7 0.7 0.7

Upper (Spec B) 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

Table S2. MCMC runtime details and diagnostics for all runs. For each inference setting, we report the number of chains and samples, acceptance rate, runtime, and maximum R̂ value. The reported maximum R̂ and threshold counts refer to inferred model parameters. The 30 day runs converged nicely, with just the 31 day, offset 160 run having a single parameter with R̂ = 1.012 > 1.01, which can effectively be interpreted as convergence. Conversely, for the 200 day run no parameters converged despite the long runtime and large sample budget. As shown in Table S3, the R̂ values for the posterior predictive values are lower, between 1.009 and 1.250, indicating that the 16 chains mostly agree on the posterior predictive even though they do not agree on the parameters. To keep computations light, we thinned the 201 days MCMC run by a factor of 100, significantly reducing memory footprint. Run 31 days, offset 0 31 days, offset 160 31 days, offset 200 201 days, offset 0

Chains Draws / chain Total samples Runtime Acceptance rate Max parameter R̂ Parameters with R̂ > 1.01 16 16 16 16

100 000 100 000 100 000 1 000 000

1 600 000 1 600 000 1 600 000 16 000 000

1.08e3s 1.09e3s 1.07e3s 1.91e4s

30

0.391 0.394 0.418 0.526

1.001 1.012 1.008 1.606

0 1 0 21

Table S3. Rank-normalized R̂ values for all parameters and posterior predictive ICU occupancy. The predictive quantities are reported as derived diagnostics of chain agreement for the model outputs. For the first day of simulation (day 0, day 160, day 200), each of the runs used the initialization process described in materials and methods. The 31 day inference windows show good convergence. For the 200 day run, mixing was harder for the parameters than for the posterior predictive values and despite the long runtime, MCMC did not achieve convergence. Quantity

31 days, offset 0 31 days, offset 160 31 days, offset 200 201 days

Inferred model parameters reduction factor 0 reduction factor 1 reduction factor 2 reduction factor 3 reduction factor 4 reduction factor 5 reduction factor 6 reduction factor 7 reduction factor 8 reduction factor 9 reduction factor 10 reduction factor 11 reduction factor 12 reduction factor 13 reduction factor 14 risk of infection from symptomatic time exposed time infected critical time infected no sym time infected sev time infected sym

1.001 1.000 – – – – – – – – – – – – – 1.000 1.000 1.001 1.001 1.000 1.000

1.010 1.003 – – – – – – – – – – – – – 1.012 1.002 1.002 1.003 1.004 1.003

1.006 1.004 – – – – – – – – – – – – – 1.008 1.002 1.004 1.003 1.004 1.003

1.554 1.284 1.449 1.217 1.380 1.543 1.606 1.544 1.309 1.160 1.066 1.126 1.214 1.080 1.155 1.317 1.228 1.369 1.401 1.067 1.470

Posterior predictive ICU occupancy predicted ICU occupancy day 0 predicted ICU occupancy day 1 predicted ICU occupancy day 2 predicted ICU occupancy day 3 predicted ICU occupancy day 4 predicted ICU occupancy day 5 predicted ICU occupancy day 6 predicted ICU occupancy day 7 predicted ICU occupancy day 8 predicted ICU occupancy day 9 predicted ICU occupancy day 10 predicted ICU occupancy day 11 predicted ICU occupancy day 12 predicted ICU occupancy day 13 predicted ICU occupancy day 14 predicted ICU occupancy day 15 predicted ICU occupancy day 16 predicted ICU occupancy day 17 predicted ICU occupancy day 18 predicted ICU occupancy day 19 predicted ICU occupancy day 20 predicted ICU occupancy day 21 predicted ICU occupancy day 22 predicted ICU occupancy day 23 predicted ICU occupancy day 24 predicted ICU occupancy day 25 predicted ICU occupancy day 26 predicted ICU occupancy day 27 predicted ICU occupancy day 28 predicted ICU occupancy day 29 predicted ICU occupancy day 30 predicted ICU occupancy day 31

– 2.000 1.000 1.000 1.000 1.000 1.001 1.001 1.001 1.001 1.001 1.001 1.001 1.001 1.001 1.001 1.001 1.001 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 –

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

– 1.106 1.091 1.069 1.049 1.033 1.022 1.017 1.021 1.031 1.046 1.062 1.077 1.088 1.094 1.095 1.093 1.088 1.081 1.073 1.064 1.056 1.050 1.044 1.041 1.041 1.041 1.043 1.044 1.047 1.051 1.055

Continued on next page

31

Quantity predicted ICU occupancy day 32 predicted ICU occupancy day 33 predicted ICU occupancy day 34 predicted ICU occupancy day 35 predicted ICU occupancy day 36 predicted ICU occupancy day 37 predicted ICU occupancy day 38 predicted ICU occupancy day 39 predicted ICU occupancy day 40 predicted ICU occupancy day 41 predicted ICU occupancy day 42 predicted ICU occupancy day 43 predicted ICU occupancy day 44 predicted ICU occupancy day 45 predicted ICU occupancy day 46 predicted ICU occupancy day 47 predicted ICU occupancy day 48 predicted ICU occupancy day 49 predicted ICU occupancy day 50 predicted ICU occupancy day 51 predicted ICU occupancy day 52 predicted ICU occupancy day 53 predicted ICU occupancy day 54 predicted ICU occupancy day 55 predicted ICU occupancy day 56 predicted ICU occupancy day 57 predicted ICU occupancy day 58 predicted ICU occupancy day 59 predicted ICU occupancy day 60 predicted ICU occupancy day 61 predicted ICU occupancy day 62 predicted ICU occupancy day 63 predicted ICU occupancy day 64 predicted ICU occupancy day 65 predicted ICU occupancy day 66 predicted ICU occupancy day 67 predicted ICU occupancy day 68 predicted ICU occupancy day 69 predicted ICU occupancy day 70 predicted ICU occupancy day 71 predicted ICU occupancy day 72 predicted ICU occupancy day 73 predicted ICU occupancy day 74 predicted ICU occupancy day 75 predicted ICU occupancy day 76 predicted ICU occupancy day 77 predicted ICU occupancy day 78 predicted ICU occupancy day 79 predicted ICU occupancy day 80 predicted ICU occupancy day 81 predicted ICU occupancy day 82 predicted ICU occupancy day 83 predicted ICU occupancy day 84 predicted ICU occupancy day 85 predicted ICU occupancy day 86 predicted ICU occupancy day 87 predicted ICU occupancy day 88 predicted ICU occupancy day 89 predicted ICU occupancy day 90 predicted ICU occupancy day 91 predicted ICU occupancy day 92 predicted ICU occupancy day 93 predicted ICU occupancy day 94 predicted ICU occupancy day 95 predicted ICU occupancy day 96

31 days, offset 0 31 days, offset 160 31 days, offset 200 201 days – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

1.059 1.062 1.068 1.073 1.079 1.086 1.092 1.099 1.104 1.109 1.111 1.111 1.109 1.103 1.094 1.083 1.070 1.058 1.047 1.040 1.037 1.035 1.038 1.040 1.044 1.048 1.052 1.055 1.059 1.061 1.063 1.065 1.065 1.065 1.064 1.062 1.061 1.060 1.061 1.064 1.070 1.078 1.088 1.098 1.107 1.114 1.117 1.114 1.108 1.097 1.083 1.069 1.057 1.051 1.054 1.069 1.095 1.128 1.163 1.195 1.220 1.238 1.247 1.250 1.247

Continued on next page

32

Quantity predicted ICU occupancy day 97 predicted ICU occupancy day 98 predicted ICU occupancy day 99 predicted ICU occupancy day 100 predicted ICU occupancy day 101 predicted ICU occupancy day 102 predicted ICU occupancy day 103 predicted ICU occupancy day 104 predicted ICU occupancy day 105 predicted ICU occupancy day 106 predicted ICU occupancy day 107 predicted ICU occupancy day 108 predicted ICU occupancy day 109 predicted ICU occupancy day 110 predicted ICU occupancy day 111 predicted ICU occupancy day 112 predicted ICU occupancy day 113 predicted ICU occupancy day 114 predicted ICU occupancy day 115 predicted ICU occupancy day 116 predicted ICU occupancy day 117 predicted ICU occupancy day 118 predicted ICU occupancy day 119 predicted ICU occupancy day 120 predicted ICU occupancy day 121 predicted ICU occupancy day 122 predicted ICU occupancy day 123 predicted ICU occupancy day 124 predicted ICU occupancy day 125 predicted ICU occupancy day 126 predicted ICU occupancy day 127 predicted ICU occupancy day 128 predicted ICU occupancy day 129 predicted ICU occupancy day 130 predicted ICU occupancy day 131 predicted ICU occupancy day 132 predicted ICU occupancy day 133 predicted ICU occupancy day 134 predicted ICU occupancy day 135 predicted ICU occupancy day 136 predicted ICU occupancy day 137 predicted ICU occupancy day 138 predicted ICU occupancy day 139 predicted ICU occupancy day 140 predicted ICU occupancy day 141 predicted ICU occupancy day 142 predicted ICU occupancy day 143 predicted ICU occupancy day 144 predicted ICU occupancy day 145 predicted ICU occupancy day 146 predicted ICU occupancy day 147 predicted ICU occupancy day 148 predicted ICU occupancy day 149 predicted ICU occupancy day 150 predicted ICU occupancy day 151 predicted ICU occupancy day 152 predicted ICU occupancy day 153 predicted ICU occupancy day 154 predicted ICU occupancy day 155 predicted ICU occupancy day 156 predicted ICU occupancy day 157 predicted ICU occupancy day 158 predicted ICU occupancy day 159 predicted ICU occupancy day 160 predicted ICU occupancy day 161

31 days, offset 0 31 days, offset 160 31 days, offset 200 201 days – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – 1.004

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

1.241 1.231 1.219 1.207 1.194 1.182 1.171 1.161 1.151 1.142 1.133 1.123 1.112 1.099 1.084 1.069 1.053 1.040 1.030 1.025 1.024 1.026 1.032 1.039 1.046 1.053 1.060 1.065 1.069 1.072 1.074 1.076 1.077 1.077 1.078 1.079 1.082 1.084 1.087 1.090 1.095 1.101 1.104 1.108 1.112 1.115 1.117 1.117 1.115 1.115 1.113 1.109 1.104 1.097 1.089 1.080 1.071 1.062 1.052 1.042 1.033 1.026 1.020 1.018 1.019

Continued on next page

33

Quantity predicted ICU occupancy day 162 predicted ICU occupancy day 163 predicted ICU occupancy day 164 predicted ICU occupancy day 165 predicted ICU occupancy day 166 predicted ICU occupancy day 167 predicted ICU occupancy day 168 predicted ICU occupancy day 169 predicted ICU occupancy day 170 predicted ICU occupancy day 171 predicted ICU occupancy day 172 predicted ICU occupancy day 173 predicted ICU occupancy day 174 predicted ICU occupancy day 175 predicted ICU occupancy day 176 predicted ICU occupancy day 177 predicted ICU occupancy day 178 predicted ICU occupancy day 179 predicted ICU occupancy day 180 predicted ICU occupancy day 181 predicted ICU occupancy day 182 predicted ICU occupancy day 183 predicted ICU occupancy day 184 predicted ICU occupancy day 185 predicted ICU occupancy day 186 predicted ICU occupancy day 187 predicted ICU occupancy day 188 predicted ICU occupancy day 189 predicted ICU occupancy day 190 predicted ICU occupancy day 191 predicted ICU occupancy day 192 predicted ICU occupancy day 193 predicted ICU occupancy day 194 predicted ICU occupancy day 195 predicted ICU occupancy day 196 predicted ICU occupancy day 197 predicted ICU occupancy day 198 predicted ICU occupancy day 199 predicted ICU occupancy day 200 predicted ICU occupancy day 201 predicted ICU occupancy day 202 predicted ICU occupancy day 203 predicted ICU occupancy day 204 predicted ICU occupancy day 205 predicted ICU occupancy day 206 predicted ICU occupancy day 207 predicted ICU occupancy day 208 predicted ICU occupancy day 209 predicted ICU occupancy day 210 predicted ICU occupancy day 211 predicted ICU occupancy day 212 predicted ICU occupancy day 213 predicted ICU occupancy day 214 predicted ICU occupancy day 215 predicted ICU occupancy day 216 predicted ICU occupancy day 217 predicted ICU occupancy day 218 predicted ICU occupancy day 219 predicted ICU occupancy day 220 predicted ICU occupancy day 221 predicted ICU occupancy day 222 predicted ICU occupancy day 223 predicted ICU occupancy day 224 predicted ICU occupancy day 225 predicted ICU occupancy day 226

31 days, offset 0 31 days, offset 160 31 days, offset 200 201 days – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

1.004 1.004 1.003 1.003 1.002 1.002 1.002 1.002 1.002 1.002 1.002 1.002 1.002 1.003 1.003 1.002 1.002 1.002 1.002 1.001 1.001 1.001 1.001 1.001 1.001 1.000 1.000 1.001 1.001 – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – –

– – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – 1.003 1.003 1.003 1.003 1.003 1.003 1.003 1.002 1.002 1.002 1.002 1.002 1.002 1.003 1.003 1.003 1.003 1.002 1.002 1.001 1.001 1.001 1.001 1.001 1.001 1.001

1.023 1.029 1.037 1.044 1.057 1.064 1.070 1.076 1.080 1.080 1.079 1.077 1.072 1.064 1.054 1.043 1.032 1.025 1.023 1.020 1.018 1.016 1.015 1.016 1.015 1.016 1.019 1.021 1.023 1.022 1.019 1.014 1.009 1.006 1.005 1.009 1.016 1.025 1.034 – – – – – – – – – – – – – – – – – – – – – – – – – –

Continued on next page

34

Quantity predicted ICU occupancy day 227 predicted ICU occupancy day 228 predicted ICU occupancy day 229 predicted ICU occupancy day 230

31 days, offset 0 31 days, offset 160 31 days, offset 200 201 days – – – –

– – – –

1.001 1.000 1.000 1.001

– – – –

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 changepoint 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.

35