# Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation

Francesco Pinotti, Julien Thézé, Xavier Bailly, Guillaume Fournié

<figure><img src="image_1.png" /></figure>

## To cite this version:

Francesco Pinotti, Julien Thézé, Xavier Bailly, Guillaume Fournié. Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation. 2026. hal-05477583

HAL Id: hal-05477583

https://hal.inrae.fr/hal-05477583v1

Preprint submitted on 26 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation

Francesco Pinotti<fnref n="1" />,<fnref n="2" />, Julien Thézé<fnref n="1" />,<fnref n="2" />, Xavier Bailly<fnref n="1" />,<fnref n="2" />, Guillaume Fournié<fnref n="1" />,<fnref n="2" />,<fnref n="3" />

## Abstract

Mathematical models are essential for understanding and forecasting infectious disease dynamics, yet parameter inference remains challenging when likelihoods are intractable or unknown. Simulation-based inference (SBI) offers a powerful alternative by leveraging model simulations in place of explicit likelihood evaluations. However, traditional SBI approaches such as Approximate Bayesian Computation often rely on ad-hoc summary statistics and suffer from reduced efficiency in high-dimensional settings. Recent advances in SBI address these limitations by employing flexible statistical or machine-learning surrogates to approximate likelihoods or posterior distributions, improving scalability and accuracy. Neural Posterior Estimation (NPE), in particular, employs neural density estimators to learn the posterior directly from simulated data, automatically extracting informative representations without requiring handcrafted statistics. Despite its promise, NPE has seen limited application in infectious disease epidemiology and phylodynamics. In this study, we assess the ability of NPE to infer parameters in mechanistic models of infectious disease dynamics. Using data from the 2014 Ebola virus outbreak in Sierra Leone, we present two case studies reflecting common outbreak-analysis scenarios. First, we fit a compartmental Susceptible-Exposed-Infectious-Recovered transmission model to time series of reported cases and deaths. Second, we fit a Birth & Death phylodynamic model with exposed and infectious hosts to an early Ebola virus phylogeny. In both applications, NPE produces accurate and reliable posterior estimates. We further show that the amortized nature of NPE facilitates model calibration and criticism, enabling rapid assessment of model fit across multiple datasets. Our results highlight the potential of NPE as a flexible, scalable, and statistically principled tool for epidemiological and phylodynamic inference. Finally, we provide detailed online tutorials illustrating the NPE workflow in both applications.

## Introduction

Mathematical models are central tools for infectious disease epidemiology. One of their most appealing features, which was apparent already in the seminal work of Kermack and McKendrick<fnref n="1" />, is the ability to explain the shape and tempo of real-world outbreaks from a handful of basic mechanistic processes, e.g. transmission, latency and recovery. Moreover,

<footnote n="1">¹Université de Lyon, INRAE, VetAgro Sup, UMR EPIA, Clermont-Auvergne-Rhône-Alpes, Lyon, France</footnote>

<footnote n="2">²Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA, Saint-Genès-Chamanelle, France</footnote>

<footnote n="3">³Department of Pathobiology and Population Sciences, Royal Veterinary College, London, United Kingdom</footnote>
they enable making epidemic forecasts, evaluate the impact of public health interventions and identify populations at heightened risk of infection²⁻⁴. Data used to inform transmission models come in various forms, e.g. incident cases, prevalence of infection and immune response, or eventually deaths⁵⁻⁷. Phylodynamic models work instead with genetic sequences and phylogenetic trees detailing ancestral relationships among sampled pathogens⁸. Despite the elegance and relative ease of implementation of transmission models, they often pose significant challenges in parameter estimation⁹. These challenges are even more pronounced for the most detailed models, e.g. microsimulations or agent-based models among others, which provide a granular representation of the target population with possibly complex interactions among its constituent elements¹⁰.

The traditional approach to parameter inference relies on the likelihood function $p(x|\theta)$, which expresses the probability of observing data $x$ given parameters $\theta$ for some underlying model. Inference via maximum likelihood seeks the parameter configurations that maximise $p(x|\theta)$, while Bayesian inference uses the posterior distribution $p(\theta|x)$, which is proportional to the likelihood. As with maximum likelihood, one may optimise the posterior or draw samples $\theta \sim p(\theta|x)$ from it. Both operations can be carried out by repeatedly evaluating the likelihood $p(x|\theta)$ for various $\theta$ values and fixed observed data $x$. However, such direct approach is feasible only for the simplest infectious disease models since $p(x|\theta)$ is often unknown or intractable. The likelihood being a function of both the model and the data, its accessibility depends on the relationship between the input $\theta$ and the model's output $x$. Direct evaluation of $p(\theta|x)$ is often difficult or even impossible, particularly for detailed, large-scale epidemiological or phylodynamic models. This limitation has motivated the development of simulation-based inference (SBI) methods, which bypass likelihood computation by learning the relationship between simulated data and parameters.

(SBI) has now become a popular and reliable alternative to traditional likelihood-based approaches¹¹. SBI revolves around the idea that a single simulation that yields pseudo-data $x$ given parameters $\theta$, is equivalent to sampling from the likelihood distribution, i.e. $x \sim p(x|\theta)$. Repeated simulations can thus be used to explore the likelihood function, even when it is intractable or unknown. This property is particularly appealing, because performing simulations is often much easier than explicitly evaluating the likelihood.

Approximate Bayesian Calculation (ABC) is perhaps the most widely used and one of the earliest methods in the SBI toolkit¹². The rejection-ABC algorithm couples simulation with a rejection step, which retains only simulated pseudo-data, and the associated parameters, that are sufficiently similar to observed data. Despite their conceptual simplicity and theoretical guarantees, all ABC-like algorithms share some common limitations. First, rejection is not straightforward to calibrate and invariably adds some approximation to posterior estimates¹³. Second, ABC suffers from the curse of dimensionality, which can be partly mitigated by reducing the data to a smaller set of summary statistics $y$. However, most ABC algorithms offer no guidance regarding the choice of the most relevant summary statistics, which is generally informed by expert knowledge. Moreover, replacing the data with summary statistics removes some useful information contained in the data, further exposing inference to a potentially significant and uncontrolled degree of approximation¹³.
Recent developments in SBI have focused on alleviating some of these limitations while making inference more sample-efficient and scalable. The general strategy is to replace the likelihood or the posterior distribution with a surrogate function that is easier to approximate and analyse¹¹. The synthetic likelihood method, for example, approximates the likelihood over summary statistics with a multivariate normal distribution¹⁴. Regression-based ABC methods instead assume an explicit parametric relationship between parameters and data that accounts also for heteroskedasticity¹⁵. Bayesian Optimization for Likelihood-Inference (BOLFI) uses a Gaussian Process model to approximate the relationship among parameters and the discrepancy between model simulations and data¹⁶. Emulation and history matching is yet another method that uses a statistical model to describe the relationship between simulations and parameters with a focus on efficiently removing ‘implausible’ parameter regions¹⁷.

Neural posterior estimation (NPE) is a recent addition to the SBI toolkit¹⁸,¹⁹. NPE aims to approximate the posterior distribution $p(\theta|x)$ with a flexible function based on one or more neural networks trained on simulated data. NPE thus learns the probabilistic mapping between data and parameters and automatically extracts informative representations of the data that serve as implicit summary statistics, effectively removing the need for handcrafted statistics or ad-hoc distance functions. NPE has been applied to complex inverse problems in physics, evolution, neuroscience and economics²⁰⁻²⁸, but its potential in epidemiology and phylodynamics has remained largely unexplored. To our knowledge, only a single study has explored the potential of NPE in the context of infectious disease dynamics²⁹. While promising, this work assumed that the posterior distribution is a multivariate normal, which may bias inference in some applications, and did not consider typical sources of data that are available during outbreaks.

The aim of this study is to assess the performance of NPE for parameter inference in mechanistic models of infectious disease dynamics. We illustrate this through two case studies that are highly relevant to epidemiological outbreak analyses. Both applications use data from the 2014 Ebola virus (EBOV) outbreak in Sierra Leone. In the first example, we use NPE to fit a Susceptible-Exposed-Infectious-Removed (SEIR) model to a time series of reported cases and deaths. In the second example we fit a phylodynamic model to an early phylogeny of EBOV sequences. Finally, we show how the amortized property of NPE enables effective calibration and criticism of epidemiological models.

# Results

## Neural Posterior Estimation

NPE aims to approximate the posterior distribution $p(\theta|x)$ from a collection of simulated samples $\{(\theta_i, x_i)\}_{i=1,...,N_{sim}}$ drawn from the joint distribution $p(x, \theta)$. This is achieved by training a Normalising Flow (NF) function $q_{\phi}(\theta|x)$ with hyperparameters $\phi$ on these simulations. NFs are highly expressive functions that not only can approximate arbitrary probability distributions with good accuracy, but that are also fast to evaluate and to sample from (see Methods for further details on NFs and training)³⁰. The trained posterior estimator $\hat{q}(\theta|x)$ should be regarded as a function of both parameters $\theta$ and data $x$. This means that if we are
interested in performing inference from a specific observation $x_{obs}$, we only need to set
$x = x_{obs}$ in $\hat{q}(\theta|x)$ to target the corresponding posterior distribution.

Once the NF is trained, switching to a different observation comes with limited computational
cost and is nearly immediate, a property called *amortisation*. As we will show later, NPE
enables large simulation studies and powerful calibration strategies that would be
computationally prohibitive for non-amortised SBI methods.

Another advantage is that NPE can work directly with raw data, bypassing the need to select and calculate summary statistics, as required for standard ABC. Instead, NPE can automatically construct informative summary statistics using a dedicated neural network (see Methods).

## 2014 Ebola virus outbreak in Sierra Leone: multivariate time series data

In this section we aimed to infer epidemiological parameters of the 2014 Ebola virus
outbreak in Sierra Leone from reported cases and deaths reported in³¹. We fitted a
deterministic transmission model to cases {$c_i$}$_{i=1,...,M}$ and deaths {$d_i$}$_{i=1,...,M}$ reported at times
{$t_i$}$_{i=1,...,M}$ (see Methods). This is an example where the likelihood function $p(x|\theta)$ is known
analytically, meaning that SBI is not necessary. Nonetheless, it allowed us to compare the
estimated posterior $\hat{q}(\theta|x)$ with the true posterior $p(\theta|x)$, which can be probed via standard
MCMC methods.

The transmission model was based on the SEIR model considered in³¹, though our analysis
deviates slightly: first, the model was fitted to newly reported rather than cumulative cases.
Second, we assumed that reported cases followed a Negative Binomial rather than a
Poisson distribution to account for overdispersion in reporting. We inferred 5 parameters: the
duration $T_0$ between the occurrence of the primary case and the first reported cases, the
basic reproduction number $R_0$, the infection fatality rate $f_{death}$, the impact of public health
control interventions $k_{control}$ and the shape parameter $s$ of the reporting model. The prior
distribution is shown in Table S1.

We trained a Masked Autoregressive Flow-based (MAF, a specific type of NF) NPE
estimator on $N_{sim} = 50000$ draws from $p(x, \theta)$. Each sample $x$ was a matrix (2, $M$) where
each row represents a different time series of length $M$. The full network architecture is
shown in Fig. S1. Once training was done (see learning details in Table S2), the NPE
estimator $\hat{q}(\theta|x)$ was extremely fast at drawing posterior samples for any observed data $x_{obs}$,
and was readily specialized to any observation $x_{obs}$ (amortization). We first exploited these
properties to test the ability of NPE to recover model parameters from 1000 simulated
observations. Fig. 1 shows that NPE estimates are accurate almost everywhere in
parameter space, except perhaps when the shape parameter $s$ is large.
<figure><img src="image_3.png" /><figcaption>**Fig. 1: Simulation study for the SEIR model.** Estimated vs true parameter values. These results are based on 1000 draws from the joint distribution $p(x, \theta)$.</figcaption></figure>

Next, we applied NPE to real data from the 2014 Ebola outbreak in Sierra Leone and compared the NPE estimates with the posterior distribution estimated via MCMC. Fig. 2 shows that the NPE posterior (red) closely matched the true posterior (turquoise), further showcasing the potential of NPE. Our analysis recovered central estimates of $R_0$, $T_0$, $f_{death}$ and $k_{control}$ reported in³¹, while their C.I.s were much narrower. This discrepancy was not specific to the NPE posterior, which we have shown to match the true posterior distribution. Rather, the use of cumulative cases in the original study is likely to have yielded overconfident estimates³².
<figure><img src="image_4.png" /><figcaption>**Fig. 2: Comparison of NPE and MCMC inference results with time series data.** Diagonal plots show posterior marginal distributions obtained using NPE (red) and MCMC (turquoise). Lines represent central estimates (solid) and 95% C.I. (dotted) from³¹. Please note that there is no reference estimate for *s* as the original model assumed Poisson reporting. Off-diagonal plots show 90% and 50% high posterior density regions (outer and inner circles) obtained from either method. The black rectangle shows the reference estimates. Both distributions are calculated from 10000 posterior samples. MCMC samples were obtained using *Stan* software³³ and the *CmdStanPy* Python interface³⁴.</figcaption></figure>

Simulations from the posterior predictive distribution were close to the observed time series of reported cases and deaths, demonstrating the quality of our fit (Fig. 3A,B). Amortized inference enabled additional statistically rigorous model checks that would be computationally expensive if not prohibitive with MCMC. Most simulation-based calibration (SBC) methods, for example, require repeated parameter inference on simulated data. Here we considered two popular methods, rank-based SBC and Tests of Accuracy with Random Points (TARP)³⁵⁻³⁷, which test whether the estimated posterior displayed the correct coverage properties. Rank-based SBC checks the order statistics (ranks) of true parameters
with respect to posterior samples: if the estimated posterior is well calibrated, the ranks are uniformly distributed (Methods). Fig. 3C shows that the empirical rank distributions were almost indistinguishable from a uniform distribution for all parameters. By definition, rank-based SBC assesses the coverage property of marginal distributions only. TARP assesses the expected coverage probability of ball-shaped regions in parameter space. Fig. 3D shows that the expected coverage probability estimated via TARP closely matched the credibility level.

<figure><img src="image_5.png" /><figcaption>**Fig. 3: Model checking for SEIR model.** (A,B) Red lines represent simulated time series from the posterior predictive distribution $p(x|x_{obs})$ which is conditioned on EBOV data $x_{obs}$ (black). Panels A and B correspond to reported cases and deaths, respectively. (C) Cumulative distribution of ranks for individual parameters as part of rank-based SBC analysis. The latter is based on 1000 draws from the joint $p(x, \theta)$. For each simulation, the rank is calculated using 1000 samples from the estimated posterior. All cumulative distributions fall within the 95% C.I. of the corresponding uniform null distribution. (D) Expected coverage probability vs credibility level as calculated from the TARP analysis. We use 1000 simulations from the joint $p(x, \theta)$ and draw 1000 samples from the estimated posterior for each simulation.</figcaption></figure>

The TARP criterion provides a sufficient condition regarding the global correctness of our posterior estimator, on average, over samples drawn from $p(x, \theta)$. While encouraging, however, it does not inform about the local correctness for any specific observation $x_{obs}$.
Local diagnostics based on classifier 2-sample (C2ST) tests like I-C2ST³⁸ have been recently developed to deal with this question. The rationale behind I-C2ST is to train a classifier to distinguish between samples from the true posterior $p(\theta|x_{obs})$ and the estimated posterior $\hat{q}(\theta|x_{obs})$: if the two distributions are indeed similar, the classifier cannot tell them apart. The test proceeded by estimating the deviation of the classifier from chance-level for some observation $x_{obs}$ and then constructed a classical p-value. In the case of EBOV case data, we found a small discrepancy with $p = 0.21$, further supporting the correctness of our estimates.

## 2014 Ebola virus outbreak in Sierra Leone: Phylodynamics inference from sequence data

To further illustrate the relevance of NPE, we analyzed a phylogenetic dataset from the same EBOV outbreak in Sierra Leone. With this purpose, we fitted a mechanistic model to time-scaled phylogenetic trees $\tau$ describing the ancestral relationships between 72 published EBOV sequences³⁹. Following previous studies⁴⁰,⁴¹, we used a Birth-Death model with Exposed and Infectious classes (BDEI model). Briefly, this model is a stochastic extension of the SEIR epidemiological model, taking into account the removal of infectious hosts due to the isolation of cases following their active finding and sampling. Contrary to the SEIR model from the previous section, this BDEI model assumed an infinite supply of susceptible hosts, which is only accurate during the early phases of an outbreak. Unknown parameters include $T_0$, $R_0$ (as above), the latent period (1/σ), the infectious period (1/μ) and the proportion ρ of sampled infectious hosts. However, as explained in the Methods and following⁴¹, which carried a similar analysis on the same data, we ignored ρ and only inferred $\theta = (T_0, R_0, \sigma, \mu)$.

Birth-Death models are the cornerstone of phylodynamics inference⁴². However, Multi-type Birth-Death models such as the BDEI model are not easily tractable, especially when parameters vary in time as in our case (see Methods). Some BEAST packages offer functionalities to fit this kind of model, but they all make numerical approximations to evaluate the likelihood⁸,⁴³. Notable software that can fit a BDEI model to genetic data include Phylodeep⁴⁴ and pyBDEI⁴⁵, but these do not implement the exact same model–e.g they assume sampling to be always active ($T_0 = 0$). Therefore, we compared inference results from NPE with a linear regression correction ABC method (ABC-LR, see Methods). This is similar to the ABC-LASSO approach used in⁴¹ on the same dataset, but uses a Variance Inflation Factor (VIF) criterion to remove highly collinear data features rather than penalizing model complexity, as done in⁴⁶.

In ABC-LR, data $x$ consists of a vector of tree summary statistics introduced in⁴¹ (see also Table S3). In contrast, NPE can extract informative summaries automatically from a phylogeny $\tau$, provided it is cast into an appropriate vector/matrix format $x$. Here we choose to use the Compact Bijective Ladderized Vectorization (CBLV) method to represent a binary rooted tree $\tau$ with $N = 72$ leaves as a $(2, N)$ matrix⁴⁴. We stress that the CBLV representation is lossless and considers all the information contained in the tree $\tau$. Both ABC-LR and NPE use a training dataset consisting of 200000 samples from $p(x, \theta)$ (see Table S2 for details on NPE training).
We first applied NPE and ABC-LR to 1000 simulated trees. Fig. 4 shows that both methods perform well at recovering the true parameters. $T_0$ and $R_0$ were the easiest to recover, while larger values of $d_E = 1/\sigma$ and $d_I = 1/\mu$ are underestimated. Importantly, both NPE and ABC-LR show similar patterns of bias, suggesting that NPE did not introduce any obvious artefact during inference. We observed similar performances in terms of mean relative error and mean relative bias (Fig. S2). NPE slightly outperformed ABC-LR in terms of parameters $R_0$, $d_E$ and $d_I$, while ABC-LR performed better on $T_0$, especially at low values. Finally, we verified that the NPE posterior estimator is well calibrated by running rank-based SBC and TARP (Fig. S3).
<figure><img src="image_6.png" /><figcaption>**Fig. 4: Simulation study for the BDEI model.** Estimated vs true parameter values. The left and right columns correspond to NPE and ABC-LR estimates, respectively. These results are based on 1000 draws from the joint distribution $p(x, \theta)$.</figcaption></figure>

Having validated methods on simulated data, we next applied NPE and ABC-LR to the empirical sequences from the 2014 Ebola virus outbreak in Sierra Leone. In practice, the
phylogenetic tree is not directly observed but should instead be inferred from sequence data,
ideally jointly with epidemiological mechanistic parameters θ. Here, we adopted a common
two-step strategy in which we first estimated the posterior distribution of the phylogenetic
tree and then fitted the BDEI model. We first used BEAST to draw 100 trees from the
posterior distribution using a simpler model with no latency (Methods). We then applied NPE
and ABC-LR to each posterior tree and aggregated posterior samples in order to account for
uncertainty in tree reconstruction. The results are shown in Fig. 5. Both methods yield similar
marginal and bivariate posterior distributions. ABC-LR estimates a slightly shorter infectious
period and larger T₀ compared to NPE. Inferred marginals align with previous estimates
from⁴⁰ although we estimate a shorter latent period d_E. We cannot make a direct comparison
with the original MCMC analysis in terms of T₀ since it modelled the start of sampling
differently. However, we can make a comparison in terms of the time T₁ to the first reported
sample (Fig. S4). NPE and ABC-LR indicate a more recent EBOV introduction than MCMC.
This might also explain why the latent period estimated via NPE and ABC-LR is also shorter
compared to the MCMC estimate. I-C2ST did not flag any issue with calibration when
applied to individual trees (p > 0.05 for all 100 posterior trees). The ABC posterior reported
in⁴¹ differed from our estimates, most likely because the initial analysis relied on a single,
maximum clade credibility (MCC) tree. We found that NPE and ABC-LR posterior
distributions were quite different from each other when considering the MCC tree (Fig. S5).
In particular, NPE estimates were significantly different from those reported in Fig. 5. Upon
closer investigation, we found that most clades in the MCC tree were poorly supported and
that the MCC tree deviated significantly from both posterior trees and trees drawn from the
prior-predictive distribution (Fig. S6). This suggests that such a summary tree is not a good
illustration of the tree posterior distribution, which impacts the efficiency of inference of the
BDEI model.
<figure><img src="image_7.png" /><figcaption>**Fig. 5: Comparison of NPE and ABC-LR inference results with EBOV sequence data.** Diagonal plots show posterior marginal distributions obtained using NPE (red) and ABC-LR (turquoise). Lines represent central estimates (solid) and 95% C.I. (dotted) from$^{40}$. Off-diagonal plots show 90% and 50% high posterior density regions (outer and inner circles) obtained from either method. The black rectangle shows the reference estimates. We do not show any reference estimate for $T_0$ since the original model used a different parameterisation (please check Fig. S4 for further details). Results are based on 10000 posterior samples obtained from 100 posterior trees estimated via BEAST$^{47}$.</figcaption></figure>

Discussion

We showed that Neural Posterior Estimation yields reliable Bayesian estimates in the context of real outbreak analyses characterised by noisy, hierarchical and multivariate data. We considered two models related to the 2014 EBOV outbreak in Sierra Leone where data consisted either of reported cases or genetic sequences. In the former case, we showed that NPE nearly recovered the true posterior distribution. In the second example the likelihood function was not available and we tested NPE against regression ABC, another popular simulation-based inference method. Both methods yielded similar estimates on both
simulated and posterior trees, which hints at NPE as a valid and competitive simulation-based inference method.

From a computational perspective, simulation and training steps are the main bottlenecks. We showed that NPE and ABC-LR are both efficient in terms of simulation requirements, using only 50000 and 200000 simulations in the SEIR and BDEI applications, respectively. Other SBI algorithms, e.g. rejection ABC or sequential variants of ABC, are typically much less efficient and may require millions of simulations⁴⁸. Training time depends mostly on neural network size, the number of training simulations and the training settings. Here we were able to complete training on an NVIDIA A40 gpu in about 1 hour in both SEIR and BDEI applications. After training, the NPE estimator is extremely fast to sample from (it took around 3s to draw 1000000 samples on our server) and can be readily applied to new observations thanks to its amortised nature.

NPE and regression ABC share the same general idea in that they both make a parametric assumption about the shape of the posterior distribution. Our implementation of regression ABC makes use of weighted heteroskedastic linear regression¹⁵, which essentially corresponds to a multivariate normal posterior density. In contrast, NPE exploits Normalising Flows that can flexibly represent complex posterior distributions. Phylodeep⁴⁴ also uses neural networks to recover mechanistic parameters from raw trees via regression. However, the method focuses on retrieving point estimates rather than a joint posterior distribution, although confidence intervals can still be obtained. Moreover, Phylodeep is restricted to a handful of popular Birth-Death models.

NPE can use a dedicated neural network to automatically learn summary statistics. This is a considerable advantage over ABC, since selecting appropriate summary statistics requires substantial care. Nonetheless, one must still choose an appropriate neural network for feature extraction, which is not a straightforward task either. One strategy would be to train multiple candidate architectures and choose the one with the smallest error on validation data. Smaller architectures require less data but risk underfitting, while more complex and larger architectures risk overfitting. Overfitting means that the inferred posterior estimator yields good results on training data but fails to generalise to new data points, e.g. those used for validation. Here we did not perform a systematic investigation as it was beyond the scope of our study, and relied on heuristic considerations and validation error as guiding principles. A general recommendation is to adopt networks with the appropriate inductive biases, i.e. networks that are wired to recognise and recover the correlational structure of the data. For example, convolutional neural networks can learn short-scale correlations that characterise, for example, image data⁴⁹. Both our examples make use of stacks of convolutional and fully-connected layers (Fig. S1), which can tackle data that are inherently hierarchical, as expected from time series and phylogenetic trees. The feature extractor in the phylodynamics example is indeed based on published networks such as Phylodeep⁴⁴. Other types of networks may also be envisioned, such as recurrent neural networks for temporal data, graph neural networks and transformers⁵⁰⁻⁵². It should be stressed that similar considerations apply also to the choice of NF type (here we used MAF), the number and size of layers and NF type-specific hyperparameters. Ultimately, users should select architectures whose complexity is commensurate with the number of training simulations.

NPE directly targets the posterior distribution $p(\theta|x)$, which has the advantage of making sampling as fast as a single forward pass of the corresponding neural network. Other neural-based methods like Neural Likelihood Estimation (NLE) and Neural Ratio Estimation (NRE) target instead the likelihood $p(x|\theta)$ and the likelihood-to-evidence ratio $p(x|\theta)/p(x)$, respectively⁵³,⁵⁴. These targets might be easier to learn than the posterior, but sampling still
requires a round of MCMC. Moreover, NLE cannot use a feature extractor, limiting its scalability to large data. Sequential versions of NPE, NLE and NRE can increase the accuracy of these algorithms on specific observations with less simulations, but the resulting posterior estimators are not amortized¹⁹,⁵³⁻⁵⁵. Finally, NPE methods that are based on diffusion models and continuous NFs have also been proposed⁵⁶,⁵⁷. These methods put less constraints on the network architecture than discrete NFs (like the MAF considered here), which can only accommodate invertible functions with a simple Jacobian structure as building blocks.

We showed that NPE can be easily incorporated within a principled Bayesian framework⁵⁸. Amortisation paves the way for a range of simulation-based calibration tools that may be too time-consuming with other inference methods⁵⁹. Rank-based SBC and related methods are relatively simple checks that can detect patterns of bias, overconfidence or underconfidence³⁵,³⁶. TARP works similarly, but provides necessary and sufficient conditions for a globally accurate posterior estimation³⁷. Local methods like I-C2ST can help detect issues with specific problem instances³⁸.

As any other inference method, NPE has some limitations. One issue is that NPE may not perform robustly in some cases, e.g. when the model is unable to generate simulations that resemble the data. In this case, NPE extrapolates beyond the training data. We have seen this phenomenon in the context of the MCC EBOV tree, which was quite different from training trees. On one hand, this example warns against using a single summary tree when phylogenetic uncertainty is high⁶⁰. On the other hand, it underlies the importance of verifying model adequacy via, e.g., prior predictive checks. Analogously, NPE may also struggle to estimate the posterior distribution when some observations lie far from training points⁶¹. Increasing the size of the training set can alleviate this issue, at the cost of increased simulation and training times.

This study considers relatively simple but relevant epidemic models that are commonly found in the literature. We believe that NPE will be even more useful in the context of more complex models, e.g. agent-based models that are often regarded as stochastic black-boxes. Indeed, NPE was shown to enable inference in the case of a model of social dynamics²⁷. Finally, we note that all our analyses focus on a past outbreak and are hence retrospective in nature. Nonetheless, we believe that NPE may also facilitate model calibration in real time. This means dealing with rapidly unfolding outbreaks where data (cases, sequences) accumulate on a daily basis and model predictions must be updated with the same frequency. One potential strategy would be to train an NF once on simulations of varying sizes and exploit the amortized property to rapidly generate posterior estimates without any need for further training.

## Methods

**EBOV Time series data.** Cumulative EBOV cases and deaths reported during the 2014 outbreak in Sierra Leone were obtained from³¹. $C_i(D_i)$ denotes cumulative cases (deaths) reported up to time $t_i$, where $i = 1, ..., M'$, $M'$ is the length of the original time series and times $t_1, t_{M'}$ correspond to dates 27/05/2014 and 20/08/2014, respectively. We then calculated newly reported cases and deaths as $c_1 = C_1$, $d_1 = D_1$ and $c_i = C_i - C_{i-1}$, $d_i = D_i - D_{i-1}$ for $i > 1$, skipping intermediate time points for which any of these differences
were negative. In this case, we search for the first time point $j > i$ for which both $C_j \ge C_i$ and $D_j \ge D_i$. The resulting time series $c_i, d_i$ are therefore shorter than the original data ($i = 1, ..., M < M'$). Cases and deaths are log-transformed according to $\log(1 + x)$ before being passed to neural networks. Moreover, the neural network z-scores each input separately.

**EBOV sequence data.** We used 72 full-length EBOV sequences collected from 26/05/2014 to 18/06/2014 in Sierra Leone and originally reported in³⁹.

**Time series SEIR epidemic model.** The numbers of susceptibles $S$, exposed $E$, infectious $I$ and recovered $R$ evolve according to the following dynamical equations:

$$ \frac{dS}{dt} = - \beta(t)S/N \qquad (1) $$

$$ \frac{dE}{dt} = \beta(t)S/N - \sigma E \qquad (2) $$

$$ \frac{dI}{dt} = \sigma E - \mu I \qquad (3) $$

$$ \frac{dR}{dt} = \mu I \qquad (4) $$

Where $N = 10^6$, $\sigma^{-1} = 5.3$ days and $\mu^{-1} = 5.61$ days. The transmission rate is defined as $\beta(t) = \mu \cdot R_0 \cdot e^{-k_{control}t}$, where parameters $R_0, k_{control}$ are to be estimated. The parameter $k_{control} > 0$ quantifies the effectiveness of control measures: the larger $k_{control}$, the faster the decay of the transmission rate. Infected individuals have a probability $f_{death}$ to die. The outbreak begins with a single infectious individual at time $t = 0$, but cases are only reported after a delay $T_0$ to be estimated.

We assume a negative binomial observation model for both cases and deaths. More in detail, the expected numbers of cases $c_i$ and deaths $d_i$ reported between $t_{i-1}$ and $t_i$ follow a negative binomial distribution with means $\sigma \int_{t_{i-1}}^{t_i} E(t)dt$ and $\mu f_{death} \int_{t_{i-1}}^{t_i} I(t)dt$, respectively. The shape parameter $s$ is the same for both data sources and has to be estimated. This observation model assumes that all cases are reported on average, but allows some uncertainty around the expected number of cases. The full likelihood function is then given by:

$$ \prod_{i=1}^{M} \binom{c_i+s-1}{c_i-1} p_i^s (1-p_i)^{c_i} \cdot \binom{d_i+s-1}{d_i-1} q_i^s (1-q_i)^{d_i}, \qquad (5) $$
$$ \text{where } p_i = s/(s + \sigma \int_{t_{i-1}}^{t_i} E(t)dt) \text{ and } q_i = s/(s + \mu f_{\text{death}} \int_{t_{i-1}}^{t_i} I(t)dt). $$

We retain only simulations where total cases fall between 100 and 1000000.

**Phylogenetic analysis of EBOV sequences.** We used BEAST v2.7.8⁴⁷ to infer and date a phylogeny of available genetic sequences. We used a Birth-Death model as the tree prior, a HKY substitution model with gamma-distributed heterogeneity and 4 rate categories⁶²,⁶³. We implemented the Birth-Death model using the BDMM-Prime BEAST package⁴³,⁶⁴, setting the sampling ratio to 0 before the oldest sample. We assume a strict molecular clock with a fixed rate 0.001984 yr⁻¹⁴⁰. We use the same prior specifications as in the BD 1 model in⁴⁰ for all parameters but the *origin* parameter (the time between the earliest sample and the start of the outbreak), whose prior we assume to be uniform between 0 and 1000 days for computational convenience. We run a single MCMC chain with 20000000 steps. We assessed convergence using Tracer v1.7.2⁶⁵ and checked that the effective sample size of each estimated parameter was above 200. We obtained the MCC tree with TreeAnnotator v2.7.7⁶⁶.

**Birth-death model with exposed individuals (BDEI model).** According to this model, an infectious individual generates another infection with rate $\beta = R_0 \cdot \mu$ or is removed with rate $\mu$ due to recovery or sampling. The sampling ratio $\rho$ denotes the probability that a removed infection is also sampled. Newly infected individuals enter the exposed compartment, they cannot be sampled and become infectious with rate $\sigma$. Moreover, we assume that sampling starts only after $t > T_0$ and set the sampling ratio to 0 for $t < T_0$. A simulation stops when exactly 72 samples have been collected, yielding a time-scaled phylogeny determined by the full transmission chain. We reject simulations that fizzle out before reaching the desired number of samples or produce more than 50000 cases. We do not try to infer the sampling ratio $\rho$ since it is well known that this parameter cannot be identified jointly with $\mu$ and⁶⁷. Following⁴¹, we assume that $\rho$ is approximately around 0.7 but still allow for some variability around this value; in practice, we draw $\rho$ from a narrow prior distribution centered around 0.7 (see Table S4) and use these values to generate training data. We run BDEI simulations using the ReMASTER⁶⁸ package in BEAST v2.7.8.

**Further details on Neural Posterior Estimation and Normalising Flows.** Mathematically, the NPE estimator $\hat{q}(\theta|x) = q_{\hat{\phi}}(\theta|x)$ corresponds to the solution $\hat{\phi}$ of the following optimisation problem:

$$ \hat{\phi} = \underset{\phi}{\operatorname{argmin}}_{\hat{\phi}} - \sum_{i=1}^{N_{\text{sim}}} \log q_{\phi}(\theta_i|x_i). \qquad (6) $$

As explained in the Results section, we choose $q_{\phi}(\theta|x)$ to be a Normalising Flow. An NF is defined as an invertible differentiable transformation $V = F_{\phi}(U)$ with parameters $\phi$ mapping a 'simple' random variable $U$ with density $p_U(u)$, e.g. a multivariate normal or uniform, into a
new variable $V$ with density $p_V(v)$. These densities are related by the change of variable formula:

$$p_V(v) = p_U(F_\phi^{-1}(v)) \left| \det \partial_u F_\phi^{-1}(v) \right|, \qquad (7)$$

where $\partial_u F_\phi^{-1}(v)$ is the Jacobian of the inverse transformation $F_\phi^{-1}$ that maps $V$ back into $U$. Under appropriate conditions, NFs enable both fast sampling and evaluation of the target distribution $q_V$. Sampling amounts to calculating $v = F_\phi(u)$ with $u \sim p_U$, while density evaluation exploits the change of variable formula (Eq. 7). Highly expressive NFs that are also computationally convenient can be constructed by stacking multiple simpler transformations with a triangular Jacobian.

Here we consider Masked Autoregressive Flows (MAF) as the backbone of NPE inference$^{69}$. MAF stacks multiple Masked Autoencoder for Density Estimation (MADE) layers. The MADE is a feed-forward neural network with a masked weight matrix. Masking refers to specific weights being set to 0 in order to enforce autoregressive dependencies among network outputs.

The target distribution in the NPE context is $q_V = q_\phi(\theta|x)$, i.e. the approximate posterior. Technically speaking, the target variable $V$ is just the parameters $\theta$, while the data $x$ represents some side-information on which the flow is conditioned to. Setting $x = x_{obs}$ in $q_\phi(\theta|x)$ for some observed data $x_{obs}$ yields an estimate of the posterior $p(\theta|x_{obs})$. This means that inference via NPE is amortized, i.e. it can be readily performed on any dataset $x_{obs}$ by simply plugging $x_{obs}$ in $q_\phi(\theta|x)$.

NPE automates the construction of summary statistics from data $x$ by training a second neural network parameterised by $\phi$ that calculates an embedding $y = G_\phi(x)$ that is then fed to the normalising flow $F_\phi$. Crucially, NPE learns $G_\phi$ and $F_\phi$ in a joint fashion by solving the following problem:

$$\hat{(\phi, \hat{\phi})} = \operatorname*{argmin}_{\phi, \phi} - \sum_{i=1}^{N_{sim}} \log q_{\phi}(\theta_i | G_{\phi}(x_i)), \qquad (8)$$

Minimising the loss function expressed in Eq. 6 and Eq. 8 is a non-convex optimisation problem with no analytical solution. Nonetheless, good solutions can be found by means of heuristic algorithms such as stochastic gradient descent.

**NPE implementation and training.** We used the NPE algorithm implemented in the sbi Python module$^{70}$. sbi provides a friendly interface and implements many neural-based conditional density estimation algorithms and SBC metrics discussed in this manuscript. It relies on the popular pytorch module$^{71}$ to build neural networks and can run on GPUs to speedup training.
**Simulation-based calibration.** SBC methods aim to assess the accuracy of posterior estimates (including MCMC) without directly accessing the likelihood function $p(x|\theta)$, whether the latter is available or not. All SBC methods discussed here are based on $N_c$ calibration samples $\{(\theta_i', x_i')\}_{i=1,...,N_c}$ drawn from the joint distribution $p(x, \theta) = p(x|\theta)p(\theta)$, which is straightforward to simulate from under our assumptions.

SBC methods may assess local or global accuracy of the estimated posterior depending on whether they refer to a specific observation or not. Rank-based SBC and TARP belong to the latter group. Rank-based SBC generates $N_s$ posterior samples $\theta_i^j \sim \hat{q}(\theta|x_i)$, $j = 1, ..., N_s$, where $x_i$ is the $i$-th calibration observation ($i = 1, ..., N_c$) simulated given $\theta = \theta_i$. It then proceeds to measure the rank $r_i$ of $\theta_i$ among values $\{\theta_i^j\}_{j=1,...,N_s}$. If $\theta$ is one-dimensional, $r_i = \sum_{j=1}^{N_s} I(\theta_i < \theta_i^j)$ and ranges between 0 and $N_s$. The rationale behind rank-based SBC is that if $\hat{q}(\theta|x_i)$ truly matches the posterior $p(\theta|x_i)$, then $\theta_i$ and $\theta_i^j$ are i.i.d. variables and the ranks $r_i$ should be uniformly distributed between 0 and $N_s$. Any deviation from the uniform distribution suggests that the posterior estimator $\hat{q}(\theta|x)$ is not perfectly calibrated and may thus expose signs of bias, overconfidence and/or underconfidence. If $\theta$ is multivariate, rank-based SBC can be applied to individual components separately.

TARP is another global calibration method that looks at the expected coverage of the posterior estimator. The TARP criterion provides a sufficient condition to assess the global accuracy of $\hat{q}(\theta|x)$, whereas other tests such as rank-based SBC only provide necessary but not sufficient conditions. As in rank-based SBC, TARP generates samples $\theta_i^j$, $j = 1, ..., N_s$ from $\hat{q}(\theta|x_i)$ for every sample $(\theta_i', x_i')$. It then generates a reference point $\theta_i^r$ from an arbitrary reference point sampling distribution (here we use the prior distribution to make these proposals) and calculates the fraction $f_i$ of samples $\theta_i^j$ within a hyper-ball centered around $\theta_i^r$ and extending up to $\theta_i^r$. Theory indicates that for well-calibrated posteriors, the distribution of fractions $f_i$ is uniform, indicating that the expected coverage probability matches the credibility level. The corresponding cumulative density should then look like a straight line with unit slope. As in rank-based SBC, any deviations from the uniform distribution may be indicative of patterns of bias, overconfidence, underconfidence, etc.

l-C2ST is a local calibration method. As explained in the Results section, l-C2ST trains a classifier to distinguish samples $(\theta, x) \sim \hat{q}(\theta|x) \cdot p(x)$ (labelled $l=0$) from samples $(\theta, x) \sim p(x, \theta)$ (labelled $l=1$). In practice, we recycle $N_c$ existing draws $(\theta_i', x_i', l_i=0)$ from $p(x, \theta)$ to generate data points $(\hat{\theta}_i', x_i', l_i=1)$ where $\hat{\theta}_i \sim \hat{q}(\theta|x_i)$. Note that this operation does not require much overhead if inference is amortized. We then train a binary classifier
on this augmented sample with $2N_c$ data points and let $h(x, \theta)$ denote the classifier probability of class $l=1$ for observation $(\theta, x)$. Let us now consider a specific observation $x_{obs}$ and samples $\{\theta_k\}_{k=1,...,N_s}$ from the estimated posterior $\hat{q}(\theta|x_{obs})$. If $\hat{q}(\theta|x_{obs}) = p(\theta|x_{obs})$, an ideal classifier (in the sense of Bayes-optimal) would not be able to tell whether these samples come from one distribution or the other and $h(x_{obs'}, \theta_k) = 0.5$. On the contrary, any deviation from chance-level predictions indicates some difference between estimated and true posterior. To measure this discrepancy, I-C2ST calculates the test statistic $\hat{t}(x_{obs}) = N_s^{-1} \sum_{k=1}^{N_s} (h(x_{obs'}, \theta_k) - 0.5)^2$. The measured discrepancy $\hat{t}(x_{obs})$ may be compatible with noise and not reflect a genuine difference between $\hat{q}(\theta|x_{obs})$ and $p(\theta|x_{obs})$. To tease these situations apart, $\hat{t}(x_{obs})$ must be compared against a null distribution of $\hat{t}$ values obtained by training additional classifiers on data with permuted labels. This procedure yields a p-value that quantifies the probability of observing a discrepancy equal to $\hat{t}(x_{obs})$ or larger under the null hypothesis that $\hat{q}(\theta|x_{obs}) = p(\theta|x_{obs})$. In this work we use $N_c = 5000$ and $N_c = 25000$ calibration samples in the SEIR and BDEI examples, respectively.

**ABC-regression with VIF selection of summary statistics.** Following$^{41}$, we use an ABC regression approach to fit the BDEI model to a given phylogeny. First, we generate $N_{sim} = 200000$ samples $\{(\theta_i', \tau_i)\}_{i=1,...,N_{sim}}$ from $p(\tau, \theta)$ and calculate a vector of summary statistics $x(\tau)$ from each tree $\tau$. We use the same summary statistics proposed in$^{41}$, with a few minor modifications (a full list of summary statistics is reported in Table S3). We reject simulations where any summary statistics is undefined and divide each summary statistics by its mean absolute deviance so that all summaries share the same scale. Second, we remove all summary statistics with zero variance and those that are highly collinear. Specifically, we remove all summary statistics with a VIF above the recommended threshold of 10. We adopt a greedy strategy, removing the summary statistics with the largest VIF, recalculating VIFs for the remaining statistics and repeating this procedure until no VIF is larger than 10. This procedure results in a reference table $X$ containing $N_{sim}$ rows and $N_{ss} = 31$ columns denoting the remaining summary statistics. Third, given an observed tree $\tau_{obs}$ with summaries $x_{obs'}$, we select the $N_{abc} = 10000$ simulations with the smallest euclidean distance $d_i = \|x_i - x_{obs}\|$ from $x_{obs}$ and use these in the subsequent regression step. Previous work showed that ABC-regression is robust with respect to the choice of $N_{abc}$, although enough samples should be included to estimate regression parameters reliably$^{41}$. ABC-regression assumes a parametric relationship $\theta_i = m(x_i) + \sigma(x_i)\zeta$ between parameters and summary statistics, where $m(x)$ is the regression function, $\sigma(x)$ is the square root of the conditional variance of $\theta$ given $x$ accounting for stochasticity, and $\zeta$ is the residual. Once estimates $\hat{m}(x)$ and $\hat{\sigma}(x)$ are available, samples $\theta_i^c$ from the approximate posterior are obtained by correcting samples $\theta_i$ according to
$\theta_i^c = \hat{m}(x_{obs}) + \frac{\hat{\sigma}(x_{obs})}{\hat{\sigma}(x_i)} (\theta_i - \hat{m}(x_i))$. Here we assume that both $m(x)$ and $\sigma(x)$ are linear functions. In practice, these functions are estimated in two steps. First, $m(x)$ is estimated via weighted least squares, i.e. by minimizing the weighted sum of residuals:

$$ \sum_{i=1}^{N_{abc}} w_i \cdot (\theta_i - m(x_i))^2, \qquad (9) $$

where weights $w_i$ penalise simulations that are far from $x_{obs}$. These weights are calculated $w_i = 1 - d_i^2/(d_{max} + 0.01)^2$, where $d_{max}$ is the largest euclidean distance among the simulations in the reference table. Next, $\log \sigma^2(x)$ is fitted to the logarithm of the empirical residuals $\hat{\epsilon}_i = \theta_i - \hat{m}(x_i)$ by minimising the quantity:

$$ \sum_{i=1}^{N_{abc}} w_i \cdot (\log \hat{\epsilon}_i^2 - \log \sigma^2(x_i))^2. \qquad (10) $$

We apply the regression-based correction after mapping model parameters to the unbounded domain via a transformation $\theta' = f(\theta)$. The transformation is either a logit or a logarithm function depending on parameter support and prior range bounds. In the phylodynamics example we apply a logarithm transformation to $R_0$ and a logit transformation to the remaining parameters. This transformation guarantees that the corrected parameters do not lie outside the desired range.

## Code availability

Code, data and relevant tutorials are available at https://forge.inrae.fr/francesco.pinotti/SBI-EPI-PHYLO-tutorial.

## Acknowledgements

We thank David Abrial for reviewing the code and helping us set up an online code repository.

## References

1. Kermack, W. O. & McKendrick, A. G. A Contribution to the Mathematical Theory of Epidemics. Proc. R. Soc. Lond. Ser. Contain. Pap. Math. Phys. Character 115, 700–721 (1927).

2. Bertozzi, A. L., Franco, E., Mohler, G., Short, M. B. & Sledge, D. The challenges of modeling and forecasting the spread of COVID-19. Proc. Natl. Acad. Sci. 117,
16732–16738 (2020).

3. Pinotti, F. et al. Modelling the transmission dynamics of H9N2 avian influenza viruses in a live bird market. *Nat. Commun.* **15**, 3494 (2024).

4. Davies, N. G. et al. Age-dependent effects in the transmission and control of COVID-19 epidemics. *Nat. Med.* **26**, 1205–1211 (2020).

5. Pritchard, E. et al. Detecting changes in population trends in infection surveillance using community SARS-CoV-2 prevalence as an exemplar. *Am. J. Epidemiol.* **193**, 1848–1860 (2024).

6. Huang, A. T. et al. Reconciling heterogeneous dengue virus infection risk estimates from different study designs. *Proc. Natl. Acad. Sci.* **122**, e2411768121 (2025).

7. Parag, K. V., Donnelly, C. A. & Zarebski, A. E. Quantifying the information in noisy epidemic curves. *Nat. Comput. Sci.* **2**, 584–594 (2022).

8. Featherstone, L. A., Zhang, J. M., Vaughan, T. G. & Duchene, S. Epidemiological inference from pathogen genomes: A review of phylodynamic models and applications. *Virus Evol.* **8**, veac045 (2022).

9. Swallow, B. et al. Challenges in estimation, uncertainty quantification and elicitation for pandemic modelling. *Epidemics* **38**, 100547 (2022).

10. Willem, L., Verelst, F., Bilcke, J., Hens, N. & Beutels, P. Lessons from a decade of individual-based models for infectious disease transmission: a systematic review (2006-2015). *BMC Infect. Dis.* **17**, 612 (2017).

11. Cranmer, K., Brehmer, J. & Louppe, G. The frontier of simulation-based inference. *Proc. Natl. Acad. Sci.* **117**, 30055–30062 (2020).

12. Tavare, S., Balding, D. J., Griffiths, R. C. & Donnelly, P. Inferring Coalescence Times from DNA Sequence Data. *Genetics* **145**, 505–518 (1997).

13. Marin, J.-M., Pudlo, P., Robert, C. P. & Ryder, R. J. Approximate Bayesian computational methods. *Stat. Comput.* **22**, 1167–1180 (2012).

14. Price, L. F., Drovandi, C. C., Lee, A. & Nott, D. J. Bayesian Synthetic Likelihood. *J. Comput. Graph. Stat.* **27**, 1–11 (2018).
15. Blum, M. G. B. & François, O. Non-linear regression models for Approximate Bayesian Computation. *Stat. Comput.* **20**, 63–73 (2010).

16. Gutmann, M. U. & Corander, J. Bayesian Optimization for Likelihood-Free Inference of Simulator-Based Statistical Models. *J. Mach. Learn. Res.* **17**, 1–47 (2016).

17. Andrianakis, I. et al. Bayesian History Matching of Complex Infectious Disease Models Using Emulation: A Tutorial and a Case Study on HIV in Uganda. *PLOS Comput. Biol.* **11**, e1003968 (2015).

18. Papamakarios, G. & Murray, I. Fast ε-free Inference of Simulation Models with Bayesian Conditional Density Estimation. in *Advances in Neural Information Processing Systems* vol. 29 (Curran Associates, Inc., 2016).

19. Lueckmann, J.-M. et al. Flexible statistical inference for mechanistic models of neural dynamics. in *Proceedings of the 31st International Conference on Neural Information Processing Systems* 1289–1299 (Curran Associates Inc., Red Hook, NY, USA, 2017).

20. Khullar, G., Nord, B., Ćiprijanović, A., Poh, J. & Xu, F. DIGS: deep inference of galaxy spectra with neural posterior estimation. *Mach. Learn. Sci. Technol.* **3**, 04LT04 (2022).

21. Vasist, M. et al. Neural posterior estimation for exoplanetary atmospheric retrieval. *Astron. Astrophys.* **672**, A147 (2023).

22. Crisostomi, M., Dey, K., Barausse, E. & Trotta, R. Neural posterior estimation with guaranteed exact coverage: The ringdown of GW150914. *Phys. Rev. D* **108**, 044029 (2023).

23. Caspi, I. et al. Mutation rate, selection, and epistasis inferred from RNA virus haplotypes via neural posterior estimation. *Virus Evol.* **9**, vead033 (2023).

24. Gonçalves, P. J. et al. Training deep neural density estimators to identify mechanistic models of neural dynamics. *eLife* **9**, e56261 (2020).

25. Rodrigues, P. L. C., Moreau, T., Louppe, G. & Gramfort, A. HNPE: leveraging global parameters for neural posterior estimation. in *Proceedings of the 35th International Conference on Neural Information Processing Systems* 13432–13443 (Curran Associates Inc., Red Hook, NY, USA, 2021).
26. Hashemi, M. et al. Amortized Bayesian inference on generative dynamical network models of epilepsy using deep neural density estimators. Neural Netw. 163, 178–194 (2023).

27. Dyer, J., Cannon, P., Farmer, J. D. & Schmon, S. M. Black-box Bayesian inference for agent-based models. J. Econ. Dyn. Control 161, 104827 (2024).

28. Fen, C. Fast Simulation-Based Bayesian Estimation of Heterogeneous and Representative Agent Models using Normalizing Flow Neural Networks. Preprint at https://doi.org/10.48550/arXiv.2203.06537 (2022).

29. Chatha, P., Bu, F., Regier, J., Snitkin, E. & Zelner, J. Neural Posterior Estimation for Stochastic Epidemic Modeling. Preprint at https://doi.org/10.48550/arXiv.2412.12967 (2024).

30. Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S. & Lakshminarayanan, B. Normalizing flows for probabilistic modeling and inference. J Mach Learn Res 22, 57:2617-57:2680 (2021).

31. Althaus, C. L. Estimating the Reproduction Number of Ebola Virus (EBOV) During the 2014 Outbreak in West Africa. PLoS Curr. 6, ecurrents.outbreaks.91afb5e0f279e7f29e7056095255b288 (2014).

32. King, A. A., Domenech de Cellès, M., Magpantay, F. M. G. & Rohani, P. Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to Ebola. Proc. Biol. Sci. 282, 20150347 (2015).

33. Stan Development Team. Stan Reference Manual.

34. cmdstanpy: Python interface to CmdStan.

35. Cook, S. R., Gelman, A. & Rubin, D. B. Validation of Software for Bayesian Models Using Posterior Quantiles. J. Comput. Graph. Stat. 15, 675–692 (2006).

36. Talts, S., Betancourt, M., Simpson, D., Vehtari, A. & Gelman, A. Validating Bayesian Inference Algorithms with Simulation-Based Calibration. Preprint at https://doi.org/10.48550/arXiv.1804.06788 (2020).

37. Lemos, P., Coogan, A., Hezaveh, Y. & Perreault-Levasseur, L. Sampling-based accuracy
testing of posterior estimators for general inference. in Proceedings of the 40th
International Conference on Machine Learning vol. 202 19256–19273 (JMLR.org,
Honolulu, Hawaii, USA, 2023).

38. Linhart, J., Gramfort, A. & Rodrigues, P. L-C2ST: Local Diagnostics for Posterior Approximations in Simulation-Based Inference. Adv. Neural Inf. Process. Syst. 36, 56384–56410 (2023).

39. Gire, S. K. et al. Genomic surveillance elucidates Ebola virus origin and transmission during the 2014 outbreak. Science 345, 1369–1372 (2014).

40. Stadler, T., Kühnert, D., Rasmussen, D. A. & du Plessis, L. Insights into the early epidemic spread of ebola in sierra leone provided by viral sequence data. PLoS Curr. 6, ecurrents.outbreaks.02bc6d927ecee7bbd33532ec8ba6a25f (2014).

41. Saulnier, E., Gascuel, O. & Alizon, S. Inferring epidemiological parameters from phylogenies using regression-ABC: A comparative study. PLOS Comput. Biol. 13, e1005416 (2017).

42. MacPherson, A., Louca, S., McLaughlin, A., Joy, J. B. & Pennell, M. W. Unifying Phylogenetic Birth–Death Models in Epidemiology and Macroevolution. Syst. Biol. 71, 172–189 (2022).

43. Vaughan, T. G. & Stadler, T. Bayesian Phylodynamic Inference of Multitype Population Trajectories Using Genomic Data. Mol. Biol. Evol. 42, msaf130 (2025).

44. Voznica, J. et al. Deep learning from phylogenies to uncover the epidemiological dynamics of outbreaks. Nat. Commun. 13, 3896 (2022).

45. Zhukova, A., Hecht, F., Maday, Y. & Gascuel, O. Fast and Accurate Maximum-Likelihood Estimation of Multi-Type Birth–Death Epidemiological Models from Phylogenetic Trees. Syst. Biol. 72, 1387–1402 (2023).

46. Danesh, G. et al. Quantifying transmission dynamics of acute hepatitis C virus infections in a heterogeneous population using sequence data. PLOS Pathog. 17, e1009916 (2021).

47. Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian
evolutionary analysis. *PLOS Comput. Biol.* **15**, e1006650 (2019).

48. *Handbook of Approximate Bayesian Computation*. (Chapman and Hall/CRC, New York, 2018). doi:10.1201/9781315117195.

49. Krizhevsky, A., Sutskever, I. & Hinton, G. E. ImageNet Classification with Deep Convolutional Neural Networks. in *Advances in Neural Information Processing Systems* vol. 25 (Curran Associates, Inc., 2012).

50. Hochreiter, S. & Schmidhuber, J. Long Short-Term Memory. *Neural Comput.* **9**, 1735-1780 (1997).

51. Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M. & Monfardini, G. The Graph Neural Network Model. *IEEE Trans. Neural Netw.* **20**, 61-80 (2009).

52. Vaswani, A. et al. Attention is All you Need. in *Advances in Neural Information Processing Systems* vol. 30 (Curran Associates, Inc., 2017).

53. Hermans, J., Begy, V. & Louppe, G. Likelihood-free MCMC with amortized approximate ratio estimators. in *Proceedings of the 37th International Conference on Machine Learning* vol. 119 4239-4248 (JMLR.org, 2020).

54. Papamakarios, G., Sterratt, D. & Murray, I. Sequential Neural Likelihood: Fast Likelihood-free Inference with Autoregressive Flows. in *Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics* 837-848 (PMLR, 2019).

55. Greenberg, D. S., Nonnenmacher, M. & Macke, J. Automatic Posterior Transformation for Likelihood-Free Inference. in (2019).

56. Wildberger, J. et al. Flow Matching for Scalable Simulation-Based Inference. *Adv. Neural Inf. Process. Syst.* **36**, 16837-16864 (2023).

57. Geffner, T., Papamakarios, G. & Mnih, A. Compositional Score Modeling for Simulation-Based Inference. in *Proceedings of the 40th International Conference on Machine Learning* 11098-11116 (PMLR, 2023).

58. Gelman, A. et al. Bayesian Workflow. Preprint at https://doi.org/10.48550/arXiv.2011.01808 (2020).
59. Lueckmann, J.-M., Boelts, J., Greenberg, D., Goncalves, P. & Macke, J. Benchmarking Simulation-Based Inference. in *Proceedings of The 24th International Conference on Artificial Intelligence and Statistics* 343–351 (PMLR, 2021).

60. O'Reilly, J. E. & Donoghue, P. C. J. The Efficacy of Consensus Tree Methods for Summarizing Phylogenetic Relationships from a Posterior Sample of Trees Estimated from Morphological Data. *Syst. Biol.* **67**, 354–362 (2018).

61. Ward, D., Cannon, P., Beaumont, M., Fasiolo, M. & Schmon, S. Robust Neural Posterior Estimation and Statistical Model Criticism. *Adv. Neural Inf. Process. Syst.* **35**, 33845–33859 (2022).

62. Stadler, T., Kühnert, D., Bonhoeffer, S. & Drummond, A. J. Birth-death skyline plot reveals temporal changes of epidemic spread in HIV and hepatitis C virus (HCV). *Proc. Natl. Acad. Sci. U. S. A.* **110**, 228–233 (2013).

63. Hasegawa, M., Kishino, H. & Yano, T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. *J. Mol. Evol.* **22**, 160–174 (1985).

64. Kühnert, D., Stadler, T., Vaughan, T. G. & Drummond, A. J. Phylodynamics with Migration: A Computational Framework to Quantify Population Structure from Genomic Data. *Mol. Biol. Evol.* **33**, 2102–2116 (2016).

65. Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. *Syst. Biol.* **67**, 901–904 (2018).

66. Drummond, A. J. & Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. *BMC Evol. Biol.* **7**, 214 (2007).

67. Louca, S., McLaughlin, A., MacPherson, A., Joy, J. B. & Pennell, M. W. Fundamental Identifiability Limits in Molecular Epidemiology. *Mol. Biol. Evol.* **38**, 4010–4024 (2021).

68. Vaughan, T. G. ReMASTER: improved phylodynamic simulation for BEAST 2.7. *Bioinformatics* **40**, btae015 (2024).

69. Papamakarios, G., Pavlakou, T. & Murray, I. Masked Autoregressive Flow for Density Estimation. in *Advances in Neural Information Processing Systems* vol. 30 (Curran
Associates, Inc., 2017).

70. Boelts, J. et al. sbi reloaded: a toolkit for simulation-based inference workflows. *J. Open Source Softw.* **10**, 7754 (2025).

71. Paszke, A. et al. PyTorch: an imperative style, high-performance deep learning library. in *Proceedings of the 33rd International Conference on Neural Information Processing Systems* 8026–8037 (Curran Associates Inc., Red Hook, NY, USA, 2019).

72. Glorot, X., Bordes, A. & Bengio, Y. Deep Sparse Rectifier Neural Networks. in *Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics* 315–323 (JMLR Workshop and Conference Proceedings, 2011).

## Supplementary material

<table><caption>Table S1: Prior distributions for the SEIR model.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Prior</td></tr></thead><tbody><tr><td>R<sub>0</sub></td><td>Basic reproduction number.</td><td>Uniform(1,4)</td></tr><tr><td>T<sub>0</sub></td><td>Delay between the occurrence of the primary case and the first cases being reported.</td><td>Uniform(0 d,120 d)</td></tr><tr><td>s</td><td>Shape of Negative Binomial reporting distribution.</td><td>Uniform(0.1,20)</td></tr><tr><td>f<sub>death</sub></td><td>Infection fatality rate.</td><td>Uniform(0.1,0.9)</td></tr><tr><td>k<sub>control</sub></td><td>Impact of control measures</td><td>Uniform(0,0.02)</td></tr></tbody></table>

<table><caption>Table S2: NF architectural details</caption><thead><tr><td>Parameter</td><td>SEIR example</td><td>BDEI example</td></tr></thead><tbody><tr><td>Normalise input features?</td><td>Yes</td><td>No</td></tr><tr><td>Number of MAF layers</td><td>5</td><td>5</td></tr><tr><td>Number of hidden units in each MAF layer</td><td>50</td><td>256</td></tr><tr><td>Dropout</td><td>0.05</td><td>0</td></tr><tr><td>Optimizer</td><td>Adam</td><td>Adam</td></tr></tbody></table>
<table><caption>Table S3: summary statistics used in ABC.</caption><tr><td>Learning rate</td><td>0.0005</td><td>0.0005</td></tr><tr><td>Proportion of training samples dedicated to estimate validation error</td><td>0.1</td><td>0.02</td></tr><tr><td>Number of iterations with no improvement in terms of validation loss before stopping</td><td>50</td><td>30</td></tr></table>

<table><thead><tr><td>Name</td><td>Description</td></tr></thead><tbody><tr><td>Maximum temporal leaf-root distance</td><td>Largest temporal distance between the root and the leaves.</td></tr><tr><td>Minimum temporal leaf-root distance</td><td>Smallest temporal distance between the root and the leaves.</td></tr><tr><td>Mean branch length</td><td>Mean branch length calculated from all branches.</td></tr><tr><td>Median branch length</td><td>Median branch length calculated from all branches.</td></tr><tr><td>Branch length standard deviation</td><td>Standard deviation of branch lengths calculated from all branch lengths</td></tr><tr><td>Mean branch length (leaves)</td><td>Mean branch length calculated from external branches only.</td></tr><tr><td>Median branch length (leaves)</td><td>Median branch length calculated from external branches only.</td></tr><tr><td>Branch length standard deviation (leaves)</td><td>Standard deviation of branch lengths calculated from external branches only.</td></tr><tr><td>Mean branch length (internal)</td><td>Mean branch length calculated from internal branches only. This is calculated over 3 consecutive time periods*.</td></tr><tr><td>Median branch length (internal)</td><td>Median branch length calculated from internal branches only. This is calculated over 3 consecutive time periods*.</td></tr><tr><td>Branch length standard deviation (internal)</td><td>Standard deviation of branch lengths calculated from internal branches only. This is calculated over 3 consecutive time periods*.</td></tr><tr><td>Mean branch length ratio</td><td>Ratio of mean length of internal and external branches. This is calculated over 3 consecutive time periods*.</td></tr></tbody></table>
<table><tr><td>Median branch length ratio</td><td>Ratio of median length of internal and external branches. This is calculated over 3 consecutive time periods*.</td></tr><tr><td>Branch length standard deviation ratio</td><td>Ratio of the standard deviation of internal and external branch lengths. This is calculated over 3 consecutive time periods*.</td></tr><tr><td>Maximum number of extant lineages</td><td>Maximum number of extant lineages at any time as seen in the LTT** plot.</td></tr><tr><td>Timing of maximum number of extant lineages</td><td>Time corresponding to the maximum number of extant lineages as seen in the LTT plot.</td></tr><tr><td>LTT Slope 1</td><td>Slope of the line connecting the origin of the LTT plot and the point in the LTT plot corresponding to the maximum number of extant lineages</td></tr><tr><td>LTT Slope 2</td><td>Slope of the line connecting the point in the LTT plot corresponding to the maximum number of extant lineages and its endpoint.</td></tr><tr><td>LTT Slope ratio</td><td>Ratio of LTT Slope 1 and LTT Slope 2.</td></tr><tr><td>Mean sampling interval</td><td>Mean time between consecutive sampling events, i.e. downward steps in the LTT plot.</td></tr><tr><td>Mean branching interval</td><td>Mean time between consecutive branching events, i.e. upward steps in the LTT plot. This is calculated over three different time periods*.</td></tr><tr><td>Colless index</td><td>Measure of tree imbalance. It is defined as the difference between the numbers of leaves in the left and right subtrees starting from an internal node, summed over all internal nodes.</td></tr><tr><td>Sackin index</td><td>Measure of tree imbalance. It is defined as the sum of the depths** of all leaves.</td></tr><tr><td>Maximum ladder size</td><td>Size of the largest ladder motif, divided by the number of leaves. A ladder motif is a chain of internal nodes each having exactly one leaf child.</td></tr><tr><td>Proportion of ladders</td><td>Sum of sizes of ladders, divided by the number of internal nodes.</td></tr><tr><td>Staircaseness 1</td><td>Proportion of internal nodes with more leaves on one subtree.</td></tr></table>
<table><tr><td>Staircaseness 2</td><td>Mean ratio of the minimum and maximum numbers of leaves on the left and right subtrees across all internal nodes.</td></tr><tr><td>Width ratio</td><td>Ratio between the maximum width and the maximum depth*** of a tree. The width count $W_k$ is the number of nodes in the tree whose depth is $k$. The width of a tree is the distance $k$ corresponding to the largest width count.</td></tr><tr><td>Width absolute difference</td><td>It is calculated as the sum of absolute differences between consecutive width counts, i.e. $\sum_k |W_{k+1} - W_k|$.</td></tr><tr><td>LTT coordinates</td><td>It is an interpolated version of the LTT plot. First, 20 equidistant points $u_l$, $l = 1, ..., 20$ are laid out over the range $[0, T]$, where $T$ is the maximum temporal depth of the tree. Second, the number $n_l$ of lineages at time $t = u_l$ is interpolated directly from the LTT plot using the interp function from the numpy module in Python. This statistic thus comprises 40 values.</td></tr></table>

* The time interval spanned by the tree is split into 3 intervals of the same length. An internal branch that overlaps with some interval, even if only partially, will be associated with that interval. Please note that a given branch may be associated with multiple intervals. A node belongs to some interval if it falls inside it.

** The Lineage-through-time (LTT) plot is a piecewise constant function denoting the number of extant lineages at any point in time along the tree. The LTT plot increases by 1 after every split in the tree and decreases by 1 after each leaf (a sampling event).

*** The depth of a node is the number of steps in the path connecting the node to the root of the tree.

<table><caption>**Table S4: Prior distributions for the BDEI model.** The priors are the same as in$^{40}$ and$^{41}$. Please note that we only sample $\rho$ from its prior to construct training data but we do not infer this parameter.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Prior</td></tr></thead><tbody><tr><td>$R_0$</td><td>Basic reproduction number.</td><td>Lognormal(0,1.25)</td></tr><tr><td>$T_0$</td><td>Delay between the occurrence of the primary case and the start of reporting.</td><td>Uniform(0 d,92 d)</td></tr><tr><td>$\sigma$</td><td>Rate of onset of infectiousness.</td><td>Gamma.shape=0.5, rate=0.167 d<sup>-1</sup>),</td></tr></tbody></table>
<table><tr><td></td><td></td><td>Restricted within the range [0.0385 d<sup>-1</sup>, 1 d<sup>-1</sup>]</td></tr><tr><td>μ</td><td>Recovery rate.</td><td>Gamma.shape=0.5, rate=0.167 d<sup>-1</sup>), Restricted within the range [0.0385 d<sup>-1</sup>, 1 d<sup>-1</sup>]</td></tr><tr><td>ρ</td><td>Sampling ratio.</td><td>Beta(70,30)</td></tr></table>

<figure><img src="image_8.png" /><figcaption>**Fig. S1: Neural network architectures to automatically extract features from raw data.** Panels A and B show the architectures used in the SEIR and BDEI examples, respectively. Each schematics shows to the neural network that maps simulation output *x* to a low-dimensional vector *y* according to $y = G_{\phi}(x)$, where *φ* denotes the ensemble of neural network parameters. The output *y* is then fed to the NF component of the NPE algorithm. Dense blocks consist of artificial neurons that connect to all input features at once. Convolutional blocks' units instead use a small sliding window to perform a 1D or 2D convolution operation on input features. Max- and Average-Pooling layers also use sliding windows to summarise and downscale input signals. Dense and convolutional layers use a RELU activation function<sup>72</sup>.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. S2: Additional calibration metrics for the BDEI model.** (A) Distribution of Relative Error for NPE (red) and ABC-LR (turquoise) for individual parameters. The Relative Error is calculated as $|\hat{\theta} - \theta_{true}|/\theta_{true}$, where $\hat{\theta}$ is the median posterior estimate and $\theta_{true}$ is the true parameter. (B) Distribution of Relative Bias. The Relative Bias is calculated as $(\hat{\theta} - \theta_{true})/\theta_{true}$. Results are based on 1000 simulations from $p(x, \theta)$ and 10000 posterior samples for each simulation. Black bars indicate median values.</figcaption></figure>

<figure><img src="image_10.png" /></figure>
**Fig. S3: Model checking for BDEI model.** (A) Cumulative distribution of ranks for individual parameters as part of rank-based SBC analysis. The latter is based on 1000 draws from the joint $p(x, \theta)$. For each simulation, the rank is calculated using 1000 samples from the estimated posterior. (B) Expected coverage probability vs credibility level as calculated from the TARP analysis. We use 5000 simulations from the joint $p(x, \theta)$ and draw 5000 samples from the estimated posterior for each simulation.

<figure><img src="image_11.png" /></figure>

**Fig. S4: Time to first sample.** Here we show the distribution of the delay $T_1$ between the occurrence of the primary case and the first reported sample in the phylodynamics example (filled). $T_1$ is different from $T_0$, which denotes the delay before reporting begins (solid lines, same as in Fig. 5). By construction, $T_1 > T_0$. $T_1$ enables a comparison with the reference analysis (black lines), which assumed the sampling ratio to be 0 before the first sample$^{40}$. Our analysis yields a shorter estimate of $T_1$ compared with the original analysis but the corresponding posterior distributions and means are compatible with each other. Results are based on 10000 simulations.
<figure><img src="image_12.png" /><figcaption>**Fig. S5: Comparison of NPE and ABC-LR inference results on EBOV MCC tree.** NPE and ABC-LR yield different predictions on the MCC tree. Importantly, both methods yield parameter estimates that are not compatible with original MCMC estimates (black dot: mean, black line: 95% C.I.). The MCC tree is not an appropriate summary of the posterior distribution over trees (see Fig. S6).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. S6: The MCC tree is not representative of posterior or simulated trees.** (A) Distributions of counts of extreme summary statistics per tree for different tree sets. We calculated the distributions of 31 summary</figcaption></figure>
statistics (Those retained after VIF selection) from 200000 BDEI trees drawn from the joint distribution p(x, θ). A
summary derived from some tree is said to be extreme if it falls outside the 95% C.I. of this distribution. The blue
line shows counts from 100 trees drawn from the joint p(x, θ) and acts as a reference expectation. The orange
line corresponds to the 100 BEAST posterior trees used in Fig. 5 in the main text. The red line denotes the
number of extreme summary statistics for the MCC tree. This figure shows that the MCC tree is an outlier that is
poorly represented by BDEI simulations. Posterior trees display some extreme features as well, though not nearly
as many as for the MCC tree. The MCC and a few posterior trees are shown in (B) and (C), respectively. Nodes in
the MCC tree are painted according to their posterior support. Most internal nodes are poorly supported.