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

Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Confirmatory Results Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation Francesco Pinotti , Julien Thézé , Xavier Bailly , Guillaume Fournié doi: https://doi.org/10.1101/2025.11.25.690436 Francesco Pinotti 1 Université de Lyon, INRAE, VetAgro Sup, UMR EPIA , Clermont-Auvergne-Rhône-Alpes, Lyon, France 2 Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA , Saint-Genès-Champanelle, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: francesco.pinotti{at}inrae.fr Julien Thézé 1 Université de Lyon, INRAE, VetAgro Sup, UMR EPIA , Clermont-Auvergne-Rhône-Alpes, Lyon, France 2 Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA , Saint-Genès-Champanelle, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Xavier Bailly 1 Université de Lyon, INRAE, VetAgro Sup, UMR EPIA , Clermont-Auvergne-Rhône-Alpes, Lyon, France 2 Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA , Saint-Genès-Champanelle, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Guillaume Fournié 1 Université de Lyon, INRAE, VetAgro Sup, UMR EPIA , Clermont-Auvergne-Rhône-Alpes, Lyon, France 2 Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA , Saint-Genès-Champanelle, France 3 Department of Pathobiology and Population Sciences , Royal Veterinary College, London, United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Data/Code Preview PDF 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 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, they enable making epidemic forecasts, evaluate the impact of public health interventions and identify populations at heightened risk of infection 2 – 4 . Data used to inform transmission models come in various forms, e.g. incident cases, prevalence of infection and immune response, or eventually deaths 5 – 7 . Phylodynamic models work instead with genetic sequences and phylogenetic trees detailing ancestral relationships among sampled pathogens 8 . Despite the elegance and relative ease of implementation of transmission models, they often pose significant challenges in parameter estimation 9 . 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 10 . The traditional approach to parameter inference relies on the likelihood function p ( x |θ), which expresses the probability of observing data x given parameters θ for some underlying model. Inference via maximum likelihood seeks the parameter configurations that maximise p ( x |θ), while Bayesian inference uses the posterior distribution p (θ| x ), which is proportional to the likelihood. As with maximum likelihood, one may optimise the posterior or draw samples θ ∼ p (θ| x ) from it. Both operations can be carried out by repeatedly evaluating the likelihood p ( x |θ) for various θ values and fixed observed data x . However, such direct approach is feasible only for the simplest infectious disease models since p ( x |θ) 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 θ and the model’s output x . Direct evaluation of p (θ| 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 11 . SBI revolves around the idea that a single simulation that yields pseudo-data x given parameters θ, is equivalent to sampling from the likelihood distribution, i.e. x ∼ p ( x |θ). 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 12 . 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 13 . 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 13 . 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 11 . The synthetic likelihood method, for example, approximates the likelihood over summary statistics with a multivariate normal distribution 14 . Regression-based ABC methods instead assume an explicit parametric relationship between parameters and data that accounts also for heteroskedasticity 15 . 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 16 . 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 17 . Neural posterior estimation (NPE) is a recent addition to the SBI toolkit 18 , 19 . NPE aims to approximate the posterior distribution p (θ| 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 20 – 28 , 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 29 . 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 (θ| x ) from a collection of simulated samples drawn from the joint distribution p ( x , θ). This is achieved by training a Normalising Flow (NF) function q φ (θ| x ) with hyperparameters φ 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) 30 . The trained posterior estimator should be regarded as a function of both parameters θ 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 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 31 . 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 |θ) is known analytically, meaning that SBI is not necessary. Nonetheless, it allowed us to compare the estimated posterior with the true posterior p (θ| x ), which can be probed via standard MCMC methods. The transmission model was based on the SEIR model considered in 31 , 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 , θ). 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 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. Download figure Open in new tab 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 , θ). 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 31 , 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 32 . Download figure Open in new tab 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 31 . 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 33 and the CmdStanPy Python interface 34 . 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) 35 – 37 , 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. Download figure Open in new tab 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 , θ). 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 , θ) and draw 1000 samples from the estimated posterior for each simulation. The TARP criterion provides a sufficient condition regarding the global correctness of our posterior estimator, on average, over samples drawn from p ( x , θ). 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 l -C2ST 38 have been recently developed to deal with this question. The rationale behind l -C2ST is to train a classifier to distinguish between samples from the true posterior p (θ| x obs ) and the estimated posterior : 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 τ describing the ancestral relationships between 72 published EBOV sequences 39 . Following previous studies 40 , 41 , 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 41 , which carried a similar analysis on the same data, we ignored ρ and only inferred θ = ( T 0 , R 0 , σ, μ). Birth-Death models are the cornerstone of phylodynamics inference 42 . 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 8 , 43 . Notable software that can fit a BDEI model to genetic data include Phylodeep 44 and pyBDEI 45 , 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 41 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 46 . In ABC-LR, data x consists of a vector of tree summary statistics introduced in 41 (see also Table S3 ). In contrast, NPE can extract informative summaries automatically from a phylogeny τ, 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 τ with N = 72 leaves as a (2, N ) matrix 44 . We stress that the CBLV representation is lossless and considers all the information contained in the tree τ. Both ABC-LR and NPE use a training dataset consisting of 200000 samples from p ( x , θ) (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/σ and d I = 1/μ 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 ). Download figure Open in new tab 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 , θ). 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 0 compared to NPE. Inferred marginals align with previous estimates from 40 although we estimate a shorter latent period d E . We cannot make a direct comparison with the original MCMC analysis in terms of T 0 since it modelled the start of sampling differently. However, we can make a comparison in terms of the time T 1 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. l -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 41 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. Download figure Open in new tab 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 . 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 48 . 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 15 , which essentially corresponds to a multivariate normal posterior density. In contrast, NPE exploits Normalising Flows that can flexibly represent complex posterior distributions. Phylodeep 44 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 49 . 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 44 . Other types of networks may also be envisioned, such as recurrent neural networks for temporal data, graph neural networks and transformers 50 – 52 . 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 (θ| 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 |θ) and the likelihood-to-evidence ratio p ( x |θ)/ p ( x ), respectively 53 , 54 . 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 19 , 53 – 55 . Finally, NPE methods that are based on diffusion models and continuous NFs have also been proposed 56 , 57 . 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 58 . Amortisation paves the way for a range of simulation-based calibration tools that may be too time-consuming with other inference methods 59 . Rank-based SBC and related methods are relatively simple checks that can detect patterns of bias, overconfidence or underconfidence 35 , 36 . TARP works similarly, but provides necessary and sufficient conditions for a globally accurate posterior estimation 37 . Local methods like l -C2ST can help detect issues with specific problem instances 38 . 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 60 . 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 61 . 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 27 . 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 31 . 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 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 ≥ C i and D j ≥ 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 39 . Time series SEIR epidemic model The numbers of susceptibles S , exposed E , infectious I and recovered R evolve according to the following dynamical equations: Where N = 10 6 , σ −1 = 5. 3 days and μ −1 = 5. 61 days. The transmission rate is defined as , 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 and deaths d reported between t i −1 and t i follow a negative binomial distribution with means and , 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: where and . We retain only simulations where total cases fall between 100 and 1000000. Phylogenetic analysis of EBOV sequences We used BEAST v2.7.8 47 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 62 , 63 . We implemented the Birth-Death model using the BDMM-Prime BEAST package 43 , 64 , setting the sampling ratio to 0 before the oldest sample. We assume a strict molecular clock with a fixed rate 0.001984 yr −140 . We use the same prior specifications as in the BD 1 model in 40 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 65 and checked that the effective sample size of each estimated parameter was above 200. We obtained the MCC tree with TreeAnnotator v2.7.7 66 . Birth-death model with exposed individuals (BDEI model) According to this model, an infectious individual generates another infection with rate β = R 0 · μ or is removed with rate μ due to recovery or sampling. The sampling ratio ρ 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 σ. 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 ρ since it is well known that this parameter cannot be identified jointly with μ and 67 . Following 41 , we assume that ρ is approximately around 0.7 but still allow for some variability around this value; in practice, we draw ρ 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 68 package in BEAST v2.7.8. Further details on Neural Posterior Estimation and Normalising Flows Mathematically, the NPE estimator corresponds to the solution of the following optimisation problem: As explained in the Results section, we choose q φ (θ| x ) to be a Normalising Flow . An NF is defined as an invertible differentiable transformation V = F φ ( U ) with parameters φ 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: where is the Jacobian of the inverse transformation 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 φ ( u ) with u ∼ 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 φ (θ| x ), i.e. the approximate posterior. Technically speaking, the target variable V is just the parameters θ, while the data x represents some side-information on which the flow is conditioned to. Setting x = x obs in q φ (θ| x ) for some observed data x obs yields an estimate of the posterior p (θ| 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 φ (θ| x ). NPE automates the construction of summary statistics from data x by training a second neural network parameterised by ϕ that calculates an embedding y = G ϕ ( x ) that is then fed to the normalising flow F φ . Crucially, NPE learns G ϕ and F φ in a joint fashion by solving the following problem: 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 |θ), whether the latter is available or not. All SBC methods discussed here are based on N C calibration samples drawn from the joint distribution p ( x , θ) = p ( x |θ) p (θ), 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 posterior samples , where x i is the i -th calibration observation ( i = 1, …, N C ) simulated given θ = θ i . It then proceeds to measure the rank r i of θ i among values . If θ is one-dimensional, and ranges between 0 and N S . The rationale behind rank-based SBC is that if truly matches the posterior p (θ| x i ), then θ i and 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 is not perfectly calibrated and may thus expose signs of bias, overconfidence and/or underconfidence. If θ 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 , whereas other tests such as rank-based SBC only provide necessary but not sufficient conditions. As in rank-based SBC, TARP generates samples from for every sample (θ i , x i ). It then generates a reference point 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 θ within a hyper-ball centered around and extending up to θ i . 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 (labelled l = 0) from samples (θ, x ) ∼ p ( x , θ) (labelled l = 1). In practice, we recycle N existing draws (θ i , x i , l i = 0) from p ( x , θ) to generate data points where . Note that this operation does not require much overhead if inference is amortized. We then train a binary classifier on this augmented sample with 2 N C data points and let h ( x , θ) denote the classifier probability of class l = 1 for observation (θ, x ). Let us now consider a specific observation x obs and samples from the estimated posterior . If 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 , θ k ) = 0. 5. On the contrary, any deviation from chance-level predictions indicates some difference between estimated and true posterior. To measure this discrepancy, l -C2ST calculates the test statistic . The measured discrepancy may be compatible with noise and not reflect a genuine difference between and p (θ| x obs ). To tease these situations apart, must be compared against a null distribution of values obtained by training additional classifiers on data with permuted labels. This procedures yields a p-value that quantifies the probability of observing a discrepancy equal to or larger under the null hypothesis that . 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 from p (τ, θ) and calculate a vector of summary statistics x (τ) from each tree τ. 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 τ 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 θ i = m ( x i ) + σ( x i ) ς between parameters and summary statistics, where m ( x ) is the regression function, σ( x ) is the square root of the conditional variance of θ given x accounting for stochasticity, and ς is the residual. Once estimates and are available, samples from the approximate posterior are obtained by correcting samples θ i according to . Here we assume that both m ( x ) and σ( 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: where weights w penalise simulations that are far from x obs . These weights are calculated , where d max is the largest euclidean distance among the simulations in the reference table. Next, log σ 2 ( x ) is fitted to the logarithm of the empirical residuals by minimising the quantity: We apply the regression-based correction after mapping model parameters to the unbounded domain via a transformation θ’ = f (θ). 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 . Supplementary material View this table: View inline View popup Download powerpoint Table S1: Prior distributions for the SEIR model. View this table: View inline View popup Table S2: NF architectural details View this table: View inline View popup Table S3: summary statistics used in ABC. View this table: View inline View popup Table S4: Prior distributions for the BDEI model. The priors are the same as in 40 and 41 . Please note that we only sample ρ from its prior to construct training data but we do not infer this parameter. Download figure Open in new tab 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 ϕ ( 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 72 . Download figure Open in new tab 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 , where is the median posterior estimate and θ true is the true parameter. (B) Distribution of Relative Bias. The Relative Bias is calculated as . Results are based on 1000 simulations from p ( x , θ) and 10000 posterior samples for each simulation. Black bars indicate median values. Download figure Open in new tab 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 , θ). 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 , θ) and draw 5000 samples from the estimated posterior for each simulation. Download figure Open in new tab 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. Download figure Open in new tab 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 ). Download figure Open in new tab 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 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. Acknowledgements We thank David Abrial for reviewing the code and helping us set up an online code repository. Funder Information Declared European Union, https://ror.org/019w4f821 , 101151817 French National Research Agency French Ministry of Higher Education and Research Footnotes https://forge.inrae.fr/francesco.pinotti/SBI-EPI-PHYLO-tutorial 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 ). OpenUrl 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 ). OpenUrl Abstract / FREE Full Text 3. Pinotti , F. et al. Modelling the transmission dynamics of H9N2 avian influenza viruses in a live bird market . Nat. Commun . 15 , 3494 ( 2024 ). OpenUrl PubMed 4. ↵ Davies , N. G. et al. Age-dependent effects in the transmission and control of COVID-19 epidemics . Nat. Med . 26 , 1205 – 1211 ( 2020 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef PubMed 6. Huang , A. T. et al. Reconciling heterogeneous dengue virus infection risk estimates from different study designs . Proc. Natl. Acad. Sci . 122 , e2411768121 ( 2025 ). OpenUrl CrossRef PubMed 7. ↵ Parag , K. V. , Donnelly , C. A. & Zarebski , A. E. Quantifying the information in noisy epidemic curves . Nat. Comput. Sci . 2 , 584 – 594 ( 2022 ). OpenUrl PubMed 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 ). OpenUrl CrossRef PubMed 9. ↵ Swallow , B. et al. Challenges in estimation, uncertainty quantification and elicitation for pandemic modelling . Epidemics 38 , 100547 ( 2022 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef PubMed 11. ↵ Cranmer , K. , Brehmer , J. & Louppe , G. The frontier of simulation-based inference . Proc. Natl. Acad. Sci . 117 , 30055 – 30062 ( 2020 ). OpenUrl Abstract / FREE Full Text 12. ↵ Tavare , S. , Balding , D. J. , Griffiths , R. C. & Donnelly , P. Inferring Coalescence Times from DNA Sequence Data . Genetics 145 , 505 – 518 ( 1997 ). OpenUrl Abstract / FREE Full Text 13. ↵ Marin , J.-M. , Pudlo , P. , Robert , C. P. & Ryder , R. J. Approximate Bayesian computational methods . Stat. Comput . 22 , 1167 – 1180 ( 2012 ). OpenUrl CrossRef PubMed 14. ↵ Price , L. F. , Drovandi , C. C. , Lee , A. & Nott , D. J. Bayesian Synthetic Likelihood . J. Comput. Graph. Stat . 27 , 1 – 11 ( 2018 ). OpenUrl 15. ↵ Blum , M. G. B. & François , O. Non-linear regression models for Approximate Bayesian Computation . Stat. Comput . 20 , 63 – 73 ( 2010 ). OpenUrl CrossRef Web of Science 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 ). OpenUrl CrossRef 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 ). OpenUrl PubMed 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 ). OpenUrl 21. Vasist , M. et al. Neural posterior estimation for exoplanetary atmospheric retrieval . Astron. Astrophys . 672 , A147 ( 2023 ). OpenUrl 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 ). OpenUrl 23. Caspi , I. et al. Mutation rate, selection, and epistasis inferred from RNA virus haplotypes via neural posterior estimation . Virus Evol . 9 , vead033 ( 2023 ). OpenUrl CrossRef PubMed 24. Gonçalves , P. J. et al. Training deep neural density estimators to identify mechanistic models of neural dynamics . eLife 9 , e56261 ( 2020 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl 28. ↵ Fen , C. Fast Simulation-Based Bayesian Estimation of Heterogeneous and Representative Agent Models using Normalizing Flow Neural Networks . Preprint at doi: 10.48550/arXiv.2203.06537 ( 2022 ). OpenUrl CrossRef 29. ↵ Chatha , P. , Bu , F. , Regier , J. , Snitkin , E. & Zelner , J. Neural Posterior Estimation for Stochastic Epidemic Modeling . Preprint at doi: 10.48550/arXiv.2412.12967 ( 2024 ). OpenUrl CrossRef 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 ). OpenUrl 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef 36. ↵ Talts , S. , Betancourt , M. , Simpson , D. , Vehtari , A. & Gelman , A. Validating Bayesian Inference Algorithms with Simulation-Based Calibration . Preprint at doi: 10.48550/arXiv.1804.06788 ( 2020 ). OpenUrl CrossRef 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 ). OpenUrl 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 ). OpenUrl 39. ↵ Gire , S. K. et al. Genomic surveillance elucidates Ebola virus origin and transmission during the 2014 outbreak . Science 345 , 1369 – 1372 ( 2014 ). OpenUrl Abstract / FREE Full Text 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef 43. ↵ Vaughan , T. G. & Stadler , T. Bayesian Phylodynamic Inference of Multitype Population Trajectories Using Genomic Data . Mol. Biol. Evol . 42 , msaf130 ( 2025 ). OpenUrl PubMed 44. ↵ Voznica , J. et al. Deep learning from phylogenies to uncover the epidemiological dynamics of outbreaks . Nat. Commun . 13 , 3896 ( 2022 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl PubMed 47. ↵ Bouckaert , R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis . PLOS Comput. Biol . 15 , e1006650 ( 2019 ). OpenUrl CrossRef PubMed 48. ↵ Handbook of Approximate Bayesian Computation . ( Chapman and Hall/CRC, New York , 2018 ). doi: 10.1201/9781315117195 . OpenUrl CrossRef 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 ). OpenUrl CrossRef PubMed Web of Science 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl 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 ). OpenUrl 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 doi: 10.48550/arXiv.2011.01808 ( 2020 ). OpenUrl CrossRef 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl 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 ). OpenUrl Abstract / FREE Full Text 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 ). OpenUrl CrossRef PubMed Web of Science 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 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef PubMed 66. ↵ Drummond , A. J. & Rambaut , A. BEAST: Bayesian evolutionary analysis by sampling trees . BMC Evol. Biol . 7 , 214 ( 2007 ). OpenUrl CrossRef PubMed 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 ). OpenUrl CrossRef PubMed 68. ↵ Vaughan , T. G. ReMASTER: improved phylodynamic simulation for BEAST 2.7 . Bioinformatics 40 , btae015 ( 2024 ). OpenUrl PubMed 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 ). OpenUrl 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 ). View the discussion thread. Back to top Previous Next Posted November 28, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation Francesco Pinotti , Julien Thézé , Xavier Bailly , Guillaume Fournié bioRxiv 2025.11.25.690436; doi: https://doi.org/10.1101/2025.11.25.690436 Share This Article: Copy Citation Tools Simulation based-inference of epidemiological and phylodynamic models via Neural Posterior Estimation Francesco Pinotti , Julien Thézé , Xavier Bailly , Guillaume Fournié bioRxiv 2025.11.25.690436; doi: https://doi.org/10.1101/2025.11.25.690436 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41968) Biophysics (21460) Cancer Biology (18603) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19910) Epidemiology (2067) Evolutionary Biology (24327) Genetics (15614) Genomics (22513) Immunology (17741) Microbiology (40423) Molecular Biology (17193) Neuroscience (88646) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7647) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)