Viral adaptation rate is negatively correlated with viral population size in 24 pediatric HIV infections (Spearman’s rank correlation: p < 0.01).

<p>The trend-line was estimated using a weighted regression analysis (weighted regression: b = -0.0054, p < 0.01). The inset illustrates the bootstrap distribution of the slope, estimated from weighted regression, which indicates that the slope is less than zero. Data points are labeled by color according to the disease progression category of each patient as follows: slow non-progressors (SNP, dark blue), moderate non-progressors (MNP, light blue), moderate progressors (MP, orange), and rapid progressors (RP, pink). Error bars representing the uncertainty in our estimate were obtained using the bootstrap procedure described in [<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004694#pcbi.1004694.ref025" target="_blank">25</a>]. Specifically, the error bars depict the lower and upper quartile estimates from 250 bootstrap samples.</p

This figure is adapted from Fig 2a in Grenfell et al. [41].

<p>A simple population genetics model predicts that absolute within-host viral adaptation rate varies non-linearly with host immune response, which has a opposing effects on viral population size and the strength of immune selection. The left-hand side of the curve (A) can explain the negative relationship observed in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004694#pcbi.1004694.g002" target="_blank">Fig 2</a>: a weak immune response corresponds to large viral population but lower selective pressure. The shaded parts of the curve indicated by B and C predicts an absence or a positive relationship, respectively, between viral adaptation rate and viral population size.</p

The number of high-frequency replacement polymorphisms (scaled by the number of codons in each alignment and the number of years of observation) is not correlated with viral population size (Spearman’s rank correlation; p > 0.05).

<p>If natural selection were weak compared to genetic drift then a negative correlation would be expected, due to an increased fixation of slightly deleterious mutations in populations of small size. The data points are labeled using the color scheme employed in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004694#pcbi.1004694.g002" target="_blank">Fig 2</a>.</p

A schematic diagram that outlines the method used to estimate the rate of molecular adaptation in serially-sampled populations.

<p>(A) Viral sequences sampled from an earlier time point (the outgroup alignment) are compared with sequences sampled at a later time point (the main alignment). Mutations on the internal branch leading to the later sample (dark grey) represent nucleotide fixations, while all remaining mutations (light grey) correspond to polymorphisms in the later sample. Replacement (non-synonymous; diamonds) and silent (synonymous; circles) mutations are distinguished. (B) A consensus of the sequences from the earlier time point is used to identify whether fixations and polymorphisms are ancestral or derived. In this example, mutation has occurred in 7 out of 9 sites in the main alignment. (C) Nucleotide site-frequencies (i.e. the frequency of each mutation in the main alignment) are calculated and probabilistically assigned to three site-frequency ranges for both silent and replacement changes. Under neutral evolution, the ratio of replacement to silent changes in the mid site-frequency range, ρ_m/ σ_m, is expected to equal to the corresponding ratio in the high site-frequency range (ρ_h/ σ_h). Excess replacement changes in the high site-frequency range thus represent adaptive substitutions driven by positive selection (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004694#pcbi.1004694.e002" target="_blank">eq 2</a>). Note that invariant sites in the alignment (i.e. sites 6 and 7 in panel B) are assigned as silent or replacement using the degeneracy of the genetic code (see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004694#pcbi.1004694.s006" target="_blank">S2 Table</a> for details). Further, the site-frequency of invariant sites is probabilistically assigned using a Beta-binomial model (see <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1004694#sec002" target="_blank">Materials and Methods</a>).</p

Time-dependent reproduction number.

Time-dependent reproduction number generated by models with the highest evidence calculated using the Laplace approximation (orange lines) and referenced TI (blue lines). Note, the fitting data in this example contains superspreading events (which leads to very high values of Rt on certain days) so is not representative of SARS-CoV-2 transmission generally.</p

Bias and variance.

Evaluating normalising constants is important across a range of topics in statistical learning, notably Bayesian model selection. However, in many realistic problems this involves the integration of analytically intractable, high-dimensional distributions, and therefore requires the use of stochastic methods such as thermodynamic integration (TI). In this paper we apply a simple but under-appreciated variation of the TI method, here referred to as referenced TI, which computes a single model’s normalising constant in an efficient way by using a judiciously chosen reference density. The advantages of the approach and theoretical considerations are set out, along with pedagogical 1 and 2D examples. The approach is shown to be useful in practice when applied to a real problem —to perform model selection for a semi-mechanistic hierarchical Bayesian model of COVID-19 transmission in South Korea involving the integration of a 200D density.</div

Log-evidence estimated by Laplace and referenced TI approximations.

In each section, model with the highest log-evidence estimated by Laplace or referenced TI method is indicated in bold. The credible intervals for log-evidence comes from calculating the quantiles of the integral from Eq 2, where the integral values were obtained from the spline interpolated using running means of the expecations per λ over all iterations.</p

Dengue type-specific global database: 1943-2013

<p>Each observation represents a unique occurrence of dengue (1 or more cases in a given year in a unique location). The fields are described as follows:</p> <p>SOURCE: internal file referencing system for source of information</p> <p>YEAR: year of dengue occurrence</p> <p>COUNTRY: country of occurrence</p> <p>REGION: region of occurrence</p> <p>ADMIN_LEVEL: whether the occurrence is recorded at the admin-0 level (country, assigned value of 0), or the admin-1 level (province, assigned value of 1)</p> <p>GAUL_AD0: Global administrative unit layer (GAUL) code for the country level</p> <p>GAUL_AD1: GAUL code for the province-level; value is -999 when occurrence is recorded at the country level</p> <p>X: longitudinal coordinate of the centroid of the administrative unit polygon (WGS 1984 datum)</p> <p>Y: latitudinal coordinate of the centroid of the administrative unit polygon (WGS 1984 datum)</p> <p>DEN1-DEN4: represent presence of reporting for each of the four DENV types; assigned value of 1=reported, 0=not reported</p> <p>N_TYPES: the cumulative number of DENV types ever reported in the given unique location by the year of the occurrence</p> <p> </p

A description of the age-stratified version of the SEIR transmission model and a definition of the priors used within the Bayesian analysis.

A description of the age-stratified version of the SEIR transmission model and a definition of the priors used within the Bayesian analysis.</p

Knot sequences.

Compartmental models that describe infectious disease transmission across subpopulations are central for assessing the impact of non-pharmaceutical interventions, behavioral changes and seasonal effects on the spread of respiratory infections. We present a Bayesian workflow for such models, including four features: (1) an adjustment for incomplete case ascertainment, (2) an adequate sampling distribution of laboratory-confirmed cases, (3) a flexible, time-varying transmission rate, and (4) a stratification by age group. Within the workflow, we benchmarked the performance of various implementations of two of these features (2 and 3). For the second feature, we used SARS-CoV-2 data from the canton of Geneva (Switzerland) and found that a quasi-Poisson distribution is the most suitable sampling distribution for describing the overdispersion in the observed laboratory-confirmed cases. For the third feature, we implemented three methods: Brownian motion, B-splines, and approximate Gaussian processes (aGP). We compared their performance in terms of the number of effective samples per second, and the error and sharpness in estimating the time-varying transmission rate over a selection of ordinary differential equation solvers and tuning parameters, using simulated seroprevalence and laboratory-confirmed case data. Even though all methods could recover the time-varying dynamics in the transmission rate accurately, we found that B-splines perform up to four and ten times faster than Brownian motion and aGPs, respectively. We validated the B-spline model with simulated age-stratified data. We applied this model to 2020 laboratory-confirmed SARS-CoV-2 cases and two seroprevalence studies from the canton of Geneva. This resulted in detailed estimates of the transmission rate over time and the case ascertainment. Our results illustrate the potential of the presented workflow including stratified transmission to estimate age-specific epidemiological parameters. The workflow is freely available in the R package HETTMO, and can be easily adapted and applied to other infectious diseases.</div