Approximate Bayesian computation

Approximate Bayesian Computation of radiocarbon and paleoenvironmental record shows population resilience on Rapa Nui (Easter Island)

Examining how past human populations responded to environmental and climatic changes is a central focus of the historical sciences. The use of summed probability distributions (SPD) of radiocarbon dates as a proxy for estimating relative population sizes provides a widely applicable method in this research area. Paleodemographic reconstructions and modeling with SPDs, however, are stymied by a lack of accepted methods for model fitting, tools for assessing the demographic impact of environmental or climatic variables, and a means for formal multi-model comparison. These deficiencies severely limit our ability to reliably resolve crucial questions of past human-environment interactions. We propose a solution using Approximate Bayesian Computation (ABC) to fit complex demographic models to observed SPDs. Using a case study from Rapa Nui (Easter Island), a location that has long been the focus of debate regarding the impact of environmental and climatic changes on its human population, we find that past populations were resilient to environmental and climatic challenges. Our findings support a growing body of evidence showing stable and sustainable communities on the island. The ABC framework offers a novel approach for exploring regions and time periods where questions of climate-induced demographic and cultural change remain unresolved.

Summed probability distributions of radiocarbon dates can be used to estimate past demography, but methods to test for associations with environmental change are lacking. Here, DiNapoli et al. propose an approach using Approximate Bayesian Computation and illustrate it in a case study of Rapa Nui.

Novel genomic approaches to study antagonistic coevolution between hosts and parasites

Host-parasite coevolution is ubiquitous, shaping genetic and phenotypic diversity and the evolutionary trajectory of interacting species. With the advances of high throughput sequencing technologies applicable to model and non-model organisms alike, it is now feasible to study in greater detail (a) the genetic underpinnings of coevolution, (b) the speed and type of dynamics at coevolving loci, and (c) the genomic consequences of coevolution. This review focuses on three recently developed approaches that leverage information from host and parasite full genome data simultaneously to pinpoint coevolving loci and draw inference on the coevolutionary history. First, co-genome-wide association study (co-GWAS) methods allow pinpointing the loci underlying host-parasite interactions. These methods focus on detecting associations between genetic variants and the outcome of experimental infection tests or on correlations between genomes of naturally infected hosts and their infecting parasites. Second, extensions to population genomics methods can detect genes under coevolution and infer the coevolutionary history, such as fitness costs. Third, correlations between host and parasite population size in time are indicative of coevolution, and polymorphism levels across independent spatially distributed populations of hosts and parasites can reveal coevolutionary loci and infer coevolutionary history. We describe the principles of these three approaches and discuss their advantages and limitations based on coevolutionary theory. We present recommendations for their application to various host (prokaryotes, fungi, plants, and animals) and parasite (viruses, bacteria, fungi, and macroparasites) species. We conclude by pointing out methodological and theoretical gaps to be filled to extract maximum information from full genome data and thereby to shed light on the molecular underpinnings of coevolution.

Imprints of Past Habitat Area Reduction on Extant Taxonomic, Functional, and Phylogenetic Composition

Past environmental changes have shaped the evolutionary and ecological diversity of extant organisms. Specifically, climatic fluctuations have made environmental conditions alternatively common or rare over time. Accordingly, most taxa have undergone restriction of their distribution to local refugia during habitat contraction, from which they could expand when suitable habitat became more common. Assessing how past restrictions in refugia have shaped species distributions and genetic diversity has motivated much research in evolutionary biology and biogeography. But there is still lack of clear synthesis on whether and how the taxonomic, functional and phylogenetic composition of extant multispecies assemblages retains the imprint of past restriction in refugia. We devised an original eco-evolutionary model to investigate the temporal dynamics of a regional species pool inhabiting a given habitat today, and which have experienced habitat reduction in the past. The model includes three components: (i) a demographic component driving stochastic changes in population sizes and extinctions due to habitat availability, (ii) a mutation and speciation component representing how divergent genotypes emerge and define new species over time, and (iii) a trait evolution component representing how trait values have changed across descendants over time. We used this model to simulate dynamics of multispecies assemblages that occupied a restricted refugia in the past and could expand their distribution subsequently. We characterized the past restriction in refugia in terms of two parameters representing the ending time of past refugia, and the extent of habitat restriction in the refugia. We characterized extant patterns of taxonomic, functional and phylogenetic diversity depending on these parameters. We found that extant relative abundances reflect the lasting influence of more recent refugia on demographic dynamics, while phylogenetic composition reflects the influence of more ancient habitat change. Extant functional diversity depends on the interplay between diversification dynamics and trait evolution, offering new options to jointly infer current trait adaptation and past trait evolution dynamics.

Ten-year projection of white-nose syndrome disease dynamics at the southern leading-edge of infection in North America

Predicting the emergence and spread of infectious diseases is critical for the effective conservation of biodiversity. White-nose syndrome (WNS), an emerging infectious disease of bats, has resulted in high mortality in eastern North America. Because the fungal causative agent Pseudogymnoascus destructans is constrained by temperature and humidity, spread dynamics may vary by geography. Environmental conditions in the southern part of the continent are different than the northeast, where disease dynamics are typically studied, making it difficult to predict how the disease will manifest. Herein, we modelled WNS pathogen spread in Texas based on cave densities and average dispersal distances of hosts, projecting these results out to 10 years. We parameterized a predictive model of WNS epidemiology and its effects on bat populations with observed cave environmental data. Our model suggests that bat populations in northern Texas will be more affected by WNS mortality than southern Texas. As such, we recommend prioritizing the preservation of large overwintering colonies of bats in north Texas through management actions. Our model illustrates that infectious disease spread and infectious disease severity can become uncoupled over a gradient of environmental variation and highlight the importance of understanding host, pathogen and environmental conditions across a breadth of environments.

Data-driven inference of high-accuracy isostable-based dynamical models in response to external inputs

Isostable reduction is a powerful technique that can be used to characterize behaviors of nonlinear dynamical systems using a basis of slowly decaying eigenfunctions of the Koopman operator. When the underlying dynamical equations are known, previously developed numerical techniques allow for high-order accuracy computation of isostable reduced models. However, in situations where the dynamical equations are unknown, few general techniques are available that provide reliable estimates of the isostable reduced equations, especially in applications where large magnitude inputs are considered. In this work, a purely data-driven inference strategy yielding high-accuracy isostable reduced models is developed for dynamical systems with a fixed point attractor. By analyzing steady-state outputs of nonlinear systems in response to sinusoidal forcing, both isostable response functions and isostable-to-output relationships can be estimated to arbitrary accuracy in an expansion performed in the isostable coordinates. Detailed examples are considered for a population of synaptically coupled neurons and for the one-dimensional Burgers' equation. While linear estimates of the isostable response functions are sufficient to characterize the dynamical behavior when small magnitude inputs are considered, the high-accuracy reduced order model inference strategy proposed here is essential when considering large magnitude inputs.

Toward Efficient Bayesian Approaches to Inference in Hierarchical Hidden Markov Models for Inferring Animal Behavior

The study of animal behavioral states inferred through hidden Markov models and similar state switching models has seen a significant increase in popularity in recent years. The ability to account for varying levels of behavioral scale has become possible through hierarchical hidden Markov models, but additional levels lead to higher complexity and increased correlation between model components. Maximum likelihood approaches to inference using the EM algorithm and direct optimization of likelihoods are more frequently used, with Bayesian approaches being less favored due to computational demands. Given these demands, it is vital that efficient estimation algorithms are developed when Bayesian methods are preferred. We study the use of various approaches to improve convergence times and mixing in Markov chain Monte Carlo methods applied to hierarchical hidden Markov models, including parallel tempering as an inference facilitation mechanism. The method shows promise for analysing complex stochastic models with high levels of correlation between components, but our results show that it requires careful tuning in order to maximize that potential.

Estimating the effect of social inequalities on the mitigation of COVID-19 across communities in Santiago de Chile

We study the spatio-temporal spread of SARS-CoV-2 in Santiago de Chile using anonymized mobile phone data from 1.4 million users, 22% of the whole population in the area, characterizing the effects of non-pharmaceutical interventions (NPIs) on the epidemic dynamics. We integrate these data into a mechanistic epidemic model calibrated on surveillance data. As of August 1, 2020, we estimate a detection rate of 102 cases per 1000 infections (90% CI: [95-112 per 1000]). We show that the introduction of a full lockdown on May 15, 2020, while causing a modest additional decrease in mobility and contacts with respect to previous NPIs, was decisive in bringing the epidemic under control, highlighting the importance of a timely governmental response to COVID-19 outbreaks. We find that the impact of NPIs on individuals' mobility correlates with the Human Development Index of comunas in the city. Indeed, more developed and wealthier areas became more isolated after government interventions and experienced a significantly lower burden of the pandemic. The heterogeneity of COVID-19 impact raises important issues in the implementation of NPIs and highlights the challenges that communities affected by systemic health and social inequalities face adapting their behaviors during an epidemic.

Competitive binding of STATs to receptor phospho-Tyr motifs accounts for altered cytokine responses

Cytokines elicit pleiotropic and non-redundant activities despite strong overlap in their usage of receptors, JAKs and STATs molecules. We use IL-6 and IL-27 to ask how two cytokines activating the same signaling pathway have different biological roles. We found that IL-27 induces more sustained STAT1 phosphorylation than IL-6, with the two cytokines inducing comparable levels of STAT3 phosphorylation. Mathematical and statistical modeling of IL-6 and IL-27 signaling identified STAT3 binding to GP130, and STAT1 binding to IL-27Rα, as the main dynamical processes contributing to sustained pSTAT1 levels by IL-27. Mutation of Tyr613 on IL-27Rα decreased IL-27-induced STAT1 phosphorylation by 80% but had limited effect on STAT3 phosphorgylation. Strong receptor/STAT coupling by IL-27 initiated a unique gene expression program, which required sustained STAT1 phosphorylation and IRF1 expression and was enriched in classical Interferon Stimulated Genes. Interestingly, the STAT/receptor coupling exhibited by IL-6/IL-27 was altered in patients with systemic lupus erythematosus (SLE). IL-6/IL-27 induced a more potent STAT1 activation in SLE patients than in healthy controls, which correlated with higher STAT1 expression in these patients. Partial inhibition of JAK activation by sub-saturating doses of Tofacitinib specifically lowered the levels of STAT1 activation by IL-6. Our data show that receptor and STATs concentrations critically contribute to shape cytokine responses and generate functional pleiotropy in health and disease.

Cryptic transmission of SARS-CoV-2 and the first COVID-19 wave in Europe and the United States

Given the narrowness of the initial testing criteria, the SARS-CoV-2 virus spread through cryptic transmission in January and February, setting the stage for the epidemic wave experienced in March and April, 2020. We use a global metapopulation epidemic model to provide a mechanistic understanding of the global dynamic underlying the establishment of the COVID-19 pandemic in Europe and the United States (US). The model is calibrated on international case introductions at the early stage of the pandemic. We find that widespread community transmission of SARS-CoV-2 was likely in several areas of Europe and the US by January 2020, and estimate that by early March, only 1 - 3 in 100 SARS-CoV-2 infections were detected by surveillance systems. Modeling results indicate international travel as the key driver of the introduction of SARS-CoV-2 with possible importation and transmission events as early as December, 2019. We characterize the resulting heterogeneous spatio-temporal spread of SARS-CoV-2 and the burden of the first COVID-19 wave (February-July 2020). We estimate infection attack rates ranging from 0.78%-15.2% in the US and 0.19%-13.2% in Europe. The spatial modeling of SARS-CoV-2 introductions and spreading provides insights into the design of innovative, model-driven surveillance systems and preparedness plans that have a broader initial capacity and indication for testing.

Neuronal circuits overcome imbalance in excitation and inhibition by adjusting connection numbers

The interplay between excitation and inhibition is crucial for neuronal circuitry in the brain. Inhibitory cell fractions in the neocortex and hippocampus are typically maintained at 15 to 30%, which is assumed to be important for stable dynamics. We have studied systematically the role of precisely controlled excitatory/inhibitory (E/I) cellular ratios on network activity using mice hippocampal cultures. Surprisingly, networks with varying E/I ratios maintain stable bursting dynamics. Interburst intervals remain constant for most ratios, except in the extremes of 0 to 10% and 90 to 100% inhibitory cells. Single-cell recordings and modeling suggest that networks adapt to chronic alterations of E/I compositions by balancing E/I connectivity. Gradual blockade of inhibition substantiates the agreement between the model and experiment and defines its limits. Combining measurements of population and single-cell activity with theoretical modeling, we provide a clearer picture of how E/I balance is preserved and where it fails in living neuronal networks.

Bayesian Inference in Y-Linked Two-Sex Branching Processes with Mutations: ABC Approach

A Y-linked two-sex branching process with mutations and blind choice of males is a suitable model for analyzing the evolution of the number of carriers of a Y-linked allele and its mutations. Such a model considers a two-sex monogamous population in which each female chooses her partner from among the male population without caring about his type (i.e., the allele he carries). In this work, we deal with the problem of estimating the main parameters of these models by developing Bayesian inference in a parametric framework. First, we consider as a sample scheme the observation of the total number of females and males up to some generation as well as the number of males of each genotype in the last generation. Subsequently, we introduce the information on the mutated males in only the last generation, obtaining in this way a second sample scheme. For both samples, we apply the Approximate Bayesian Computation (ABC) method to approximate the posterior distributions of the main parameters of the model. The accuracy of the procedure based on these samples is illustrated and discussed by way of simulated examples.

Estimated Methicillin-Resistant Staphylococcus aureus Decolonization in Intensive Care Units Associated With Single-Application Chlorhexidine Gluconate or Mupirocin

IMPORTANCE

Chlorhexidine gluconate (CHG) and mupirocin are widely used to decolonize patients with methicillin-resistant Staphylococcus aureus (MRSA) and reduce risks associated with infection in hospitalized populations. Quantifying the association of an application of CHG alone or in combination with mupirocin with risk of MRSA infection is important for studies evaluating alternative decolonization strategies or schedules and for identifying whether there is room for improved decolonizing agents.

OBJECTIVE

To estimate the proportion of patients with MRSA decolonized per application of CHG and mupirocin from existing population-level studies.

DESIGN, SETTING, AND PARTICIPANTS

A stochastic mathematical model of an 18-bed intensive care unit (ICU) in an academic medical center operating over 1 year was used to estimate parameters for the proportion of simulated patients with MRSA decolonized per application of CHG and mupirocin. The model was conducted using approximate bayesian computation with data from an existing meta-analysis of studies conducted from February 2005 through January 2015. Data were analyzed from January 2018 through November 2019.

EXPOSURE

A universal decolonization protocol for colonized patients in the ICU using CHG or CHG and mupirocin in combination was simulated.

MAIN OUTCOMES AND MEASURES

The proportion of patients with MRSA decolonized per application of CHG and mupirocin was estimated.

RESULTS

The estimated proportion of patients with MRSA decolonized per application of CHG was 0.15 (95% credible interval, 0.01-0.42), and the estimated proportion per application of mupirocin in conjunction with CHG was 0.15 (95% credible interval, 0.01-0.54). A lag in colonization detection was associated with decreases in the CHG estimate (0.11; 95% credible interval, 0.01-0.30) and mupirocin estimate (0.10; 95% credible interval, 0.00-0.34), which were sensitive to the value of the modeled contact rate between nurses and patients. A 1% increase in the value of this parameter was associated with a 0.73% increase in the estimated combined outcomes associated with CHG and mupirocin (95% CI: 0.71, 0.75). Gaps longer than 24 hours in the administration of decolonizing agents were associated with a decrease of within-ICU MRSA transmission. Compared with a mean (SD) of 1.23 (0.27) acquisitions per 1000 patient-days in scenarios with no decolonizing bathing, a bathing protocol administering CHG and mupirocin every 120 hours was associated with a mean (SD) acquisition rate of 1.03 (0.24) acquisitions per 1000 patient days, a 16.3% decrease (95% CI, 14.7%-18.0%; P > .001).

CONCLUSIONS AND RELEVANCE

These findings suggest that there may be room for significant improvement in anti-MRSA disinfectants, including the compounds themselves and their delivery mechanisms. Despite the decolonization estimates found in this study, these agents are associated with robust outcomes after delays in administration, which may help in alleviating concerns over patient comfort and toxic effects.

Fixation patterns in simple choice reflect optimal information sampling

Simple choices (e.g., eating an apple vs. an orange) are made by integrating noisy evidence that is sampled over time and influenced by visual attention; as a result, fluctuations in visual attention can affect choices. But what determines what is fixated and when? To address this question, we model the decision process for simple choice as an information sampling problem, and approximate the optimal sampling policy. We find that it is optimal to sample from options whose value estimates are both high and uncertain. Furthermore, the optimal policy provides a reasonable account of fixations and choices in binary and trinary simple choice, as well as the differences between the two cases. Overall, the results show that the fixation process during simple choice is influenced dynamically by the value estimates computed during the decision process, in a manner consistent with optimal information sampling.

The epidemiologic impact and cost-effectiveness of new tuberculosis vaccines on multidrug-resistant tuberculosis in India and China

BACKGROUND

Despite recent advances through the development pipeline, how novel tuberculosis (TB) vaccines might affect rifampicin-resistant and multidrug-resistant tuberculosis (RR/MDR-TB) is unknown. We investigated the epidemiologic impact, cost-effectiveness, and budget impact of hypothetical novel prophylactic prevention of disease TB vaccines on RR/MDR-TB in China and India.

METHODS

We constructed a deterministic, compartmental, age-, drug-resistance- and treatment history-stratified dynamic transmission model of tuberculosis. We introduced novel vaccines from 2027, with post- (PSI) or both pre- and post-infection (P&PI) efficacy, conferring 10 years of protection, with 50% efficacy. We measured vaccine cost-effectiveness over 2027-2050 as USD/DALY averted-against 1-times GDP/capita, and two healthcare opportunity cost-based (HCOC), thresholds. We carried out scenario analyses.

RESULTS

By 2050, the P&PI vaccine reduced RR/MDR-TB incidence rate by 71% (UI: 69-72) and 72% (UI: 70-74), and the PSI vaccine by 31% (UI: 30-32) and 44% (UI: 42-47) in China and India, respectively. In India, we found both USD 10 P&PI and PSI vaccines cost-effective at the 1-times GDP and upper HCOC thresholds and P&PI vaccines cost-effective at the lower HCOC threshold. In China, both vaccines were cost-effective at the 1-times GDP threshold. P&PI vaccine remained cost-effective at the lower HCOC threshold with 49% probability and PSI vaccines at the upper HCOC threshold with 21% probability. The P&PI vaccine was predicted to avert 0.9 million (UI: 0.8-1.1) and 1.1 million (UI: 0.9-1.4) second-line therapy regimens in China and India between 2027 and 2050, respectively.

CONCLUSIONS

Novel TB vaccination is likely to substantially reduce the future burden of RR/MDR-TB, while averting the need for second-line therapy. Vaccination may be cost-effective depending on vaccine characteristics and setting.

A mesoscopic simulator to uncover heterogeneity and evolutionary dynamics in tumors

Increasingly complex in silico modeling approaches offer a way to simultaneously access cancerous processes at different spatio-temporal scales. High-level models, such as those based on partial differential equations, are computationally affordable and allow large tumor sizes and long temporal windows to be studied, but miss the discrete nature of many key underlying cellular processes. Individual-based approaches provide a much more detailed description of tumors, but have difficulties when trying to handle full-sized real cancers. Thus, there exists a trade-off between the integration of macroscopic and microscopic information, now widely available, and the ability to attain clinical tumor sizes. In this paper we put forward a stochastic mesoscopic simulation framework that incorporates key cellular processes during tumor progression while keeping computational costs to a minimum. Our framework captures a physical scale that allows both the incorporation of microscopic information, tracking the spatio-temporal emergence of tumor heterogeneity and the underlying evolutionary dynamics, and the reconstruction of clinically sized tumors from high-resolution medical imaging data, with the additional benefit of low computational cost. We illustrate the functionality of our modeling approach for the case of glioblastoma, a paradigm of tumor heterogeneity that remains extremely challenging in the clinical setting.

Digging into Toxoplasma gondii infections via soil: A quantitative microbial risk assessment approach

Soil has been identified as an important source of exposure to a variety of chemical and biological contaminants. Toxoplasma gondii is one of those potential biological contaminants associated with serious health effects in pregnant women and immunocompromised patients. Gardening or consumption of homegrown vegetables may present an important route of T. gondii infection via accidental ingestion of soil. In the Netherlands, there is quantitative information on the risk of T. gondii infection via meat products, but not on the risk of infection through soil. The objective of this study was to develop a quantitative microbial risk assessment (QMRA) model for estimating the risk associated with T. gondii exposure via accidental soil ingestion in the Netherlands. In order to obtain the needed information, a magnetic capture method for detection of T. gondii oocysts in soil samples was developed, and T. gondii DNA was detected using qPCR targeting the 529 bp repeat element. The method was shown to provide 95% probability of detection (95% CI: 88-100%) when at least 34 oocysts are present in 25 g of soil. T. gondii DNA was detected in 5 of 148 soil samples with interpretable results (3%, 95% CI: 1.5-7.7%). Results for 18 samples were not interpretable due to PCR inhibition. The estimated amount of oocysts presented in qPCR positive samples was quantified by a linear model, and the amount varied from 8 to 478 in 25 g of soil. The estimated incidence rate of T. gondii infection from the QMRA model via soil varied from 0.3 to 1.8 per 1000 individuals per day. Several data gaps (e.g., soil contamination/ingestion and oocysts viability) have been identified in this study, the structure of the model can be applied to obtain more accurate estimates of the risk of T. gondii infection via soil when data become available.

Markov chain Monte Carlo with Gaussian processes for fast parameter estimation and uncertainty quantification in a 1D fluid-dynamics model of the pulmonary circulation

The past few decades have witnessed an explosive synergy between physics and the life sciences. In particular, physical modelling in medicine and physiology is a topical research area. The present work focuses on parameter inference and uncertainty quantification in a 1D fluid-dynamics model for quantitative physiology: the pulmonary blood circulation. The practical challenge is the estimation of the patient-specific biophysical model parameters, which cannot be measured directly. In principle this can be achieved based on a comparison between measured and predicted data. However, predicting data requires solving a system of partial differential equations (PDEs), which usually have no closed-form solution, and repeated numerical integrations as part of an adaptive estimation procedure are computationally expensive. In the present article, we demonstrate how fast parameter estimation combined with sound uncertainty quantification can be achieved by a combination of statistical emulation and Markov chain Monte Carlo (MCMC) sampling. We compare a range of state-of-the-art MCMC algorithms and emulation strategies, and assess their performance in terms of their accuracy and computational efficiency. The long-term goal is to develop a method for reliable disease prognostication in real time, and our work is an important step towards an automatic clinical decision support system.

Point process models for sweat gland activation observed with noise

The aim of this article is to construct spatial models for the activation of sweat glands for healthy subjects and subjects suffering from peripheral neuropathy by using videos of sweating recorded from the subjects. The sweat patterns are regarded as realizations of spatial point processes and two point process models for the sweat gland activation and two methods for inference are proposed. Several image analysis steps are needed to extract the point patterns from the videos and some incorrectly identified sweat gland locations may be present in the data. To take into account the errors, we either include an error term in the point process model or use an estimation procedure that is robust with respect to the errors.

Bayesian genome scale modelling identifies thermal determinants of yeast metabolism

The molecular basis of how temperature affects cell metabolism has been a long-standing question in biology, where the main obstacles are the lack of high-quality data and methods to associate temperature effects on the function of individual proteins as well as to combine them at a systems level. Here we develop and apply a Bayesian modeling approach to resolve the temperature effects in genome scale metabolic models (GEM). The approach minimizes uncertainties in enzymatic thermal parameters and greatly improves the predictive strength of the GEMs. The resulting temperature constrained yeast GEM uncovers enzymes that limit growth at superoptimal temperatures, and squalene epoxidase (ERG1) is predicted to be the most rate limiting. By replacing this single key enzyme with an ortholog from a thermotolerant yeast strain, we obtain a thermotolerant strain that outgrows the wild type, demonstrating the critical role of sterol metabolism in yeast thermosensitivity. Therefore, apart from identifying thermal determinants of cell metabolism and enabling the design of thermotolerant strains, our Bayesian GEM approach facilitates modelling of complex biological systems in the absence of high-quality data and therefore shows promise for becoming a standard tool for genome scale modeling.

While temperature impacts the function of all cellular components, it’s hard to rule out how the temperature dependence of cell phenotypes emerged from the dependence of individual components. Here, the authors develop a Bayesian genome scale modelling approach to identify thermal determinants of yeast metabolism.

Hindsight is 2020 vision: a characterisation of the global response to the COVID-19 pandemic

BACKGROUND

The global impact of COVID-19 and the country-specific responses to the pandemic provide an unparalleled opportunity to learn about different patterns of the outbreak and interventions. We model the global pattern of reported COVID-19 cases during the primary response period, with the aim of learning from the past to prepare for the future.

METHODS

Using Bayesian methods, we analyse the response to the COVID-19 outbreak for 158 countries for the period 22 January to 9 June 2020. This encompasses the period in which many countries imposed a variety of response measures and initial relaxation strategies. Instead of modelling specific intervention types and timings for each country explicitly, we adopt a stochastic epidemiological model including a feedback mechanism on virus transmission to capture complex nonlinear dynamics arising from continuous changes in community behaviour in response to rising case numbers. We analyse the overall effect of interventions and community responses across diverse regions. This approach mitigates explicit consideration of issues such as period of infectivity and public adherence to government restrictions.

RESULTS

Countries with the largest cumulative case tallies are characterised by a delayed response, whereas countries that avoid substantial community transmission during the period of study responded quickly. Countries that recovered rapidly also have a higher case identification rate and small numbers of undocumented community transmission at the early stages of the outbreak. We also demonstrate that uncertainty in numbers of undocumented infections dramatically impacts the risk of multiple waves. Our approach is also effective at pre-empting potential flare-ups.

CONCLUSIONS

We demonstrate the utility of modelling to interpret community behaviour in the early epidemic stages. Two lessons learnt that are important for the future are: i) countries that imposed strict containment measures early in the epidemic fared better with respect to numbers of reported cases; and ii) broader testing is required early in the epidemic to understand the magnitude of undocumented infections and recover rapidly. We conclude that clear patterns of containment are essential prior to relaxation of restrictions and show that modelling can provide insights to this end.

Identifiability analysis for stochastic differential equation models in systems biology

Mathematical models are routinely calibrated to experimental data, with goals ranging from building predictive models to quantifying parameters that cannot be measured. Whether or not reliable parameter estimates are obtainable from the available data can easily be overlooked. Such issues of parameter identifiability have important ramifications for both the predictive power of a model, and the mechanistic insight that can be obtained. Identifiability analysis is well-established for deterministic, ordinary differential equation (ODE) models, but there are no commonly adopted methods for analysing identifiability in stochastic models. We provide an accessible introduction to identifiability analysis and demonstrate how existing ideas for analysis of ODE models can be applied to stochastic differential equation (SDE) models through four practical case studies. To assess structural identifiability, we study ODEs that describe the statistical moments of the stochastic process using open-source software tools. Using practically motivated synthetic data and Markov chain Monte Carlo methods, we assess parameter identifiability in the context of available data. Our analysis shows that SDE models can often extract more information about parameters than deterministic descriptions. All code used to perform the analysis is available on Github.

Sometimes hidden but always there: the assumptions underlying genetic inference of demographic histories

Demographic processes directly affect patterns of genetic variation within contemporary populations as well as future generations, allowing for demographic inference from patterns of both present-day and past genetic variation. Advances in laboratory procedures, sequencing and genotyping technologies in the past decades have resulted in massive increases in high-quality genome-wide genetic data from present-day populations and allowed retrieval of genetic data from archaeological material, also known as ancient DNA. This has resulted in an explosion of work exploring past changes in population size, structure, continuity and movement. However, as genetic processes are highly stochastic, patterns of genetic variation only indirectly reflect demographic histories. As a result, past demographic processes need to be reconstructed using an inferential approach. This usually involves comparing observed patterns of variation with model expectations from theoretical population genetics. A large number of approaches have been developed based on different population genetic models that each come with assumptions about the data and underlying demography. In this article I review some of the key models and assumptions underlying the most commonly used approaches for past demographic inference and their consequences for our ability to link the inferred demographic processes to the archaeological and climate records. This article is part of the theme issue 'Cross-disciplinary approaches to prehistoric demography'.

Computational evidence for hierarchically structured reinforcement learning in humans

Humans have the fascinating ability to achieve goals in a complex and constantly changing world, still surpassing modern machine-learning algorithms in terms of flexibility and learning speed. It is generally accepted that a crucial factor for this ability is the use of abstract, hierarchical representations, which employ structure in the environment to guide learning and decision making. Nevertheless, how we create and use these hierarchical representations is poorly understood. This study presents evidence that human behavior can be characterized as hierarchical reinforcement learning (RL). We designed an experiment to test specific predictions of hierarchical RL using a series of subtasks in the realm of context-based learning and observed several behavioral markers of hierarchical RL, such as asymmetric switch costs between changes in higher-level versus lower-level features, faster learning in higher-valued compared to lower-valued contexts, and preference for higher-valued compared to lower-valued contexts. We replicated these results across three independent samples. We simulated three models-a classic RL, a hierarchical RL, and a hierarchical Bayesian model-and compared their behavior to human results. While the flat RL model captured some aspects of participants' sensitivity to outcome values, and the hierarchical Bayesian model captured some markers of transfer, only hierarchical RL accounted for all patterns observed in human behavior. This work shows that hierarchical RL, a biologically inspired and computationally simple algorithm, can capture human behavior in complex, hierarchical environments and opens the avenue for future research in this field.

The Origins of Coca: Museum Genomics Reveals Multiple Independent Domestications from Progenitor Erythroxylum gracilipes

Coca is the natural source of cocaine as well as a sacred and medicinal plant farmed by South American Amerindians and mestizos. The coca crop comprises four closely related varieties classified into two species (Amazonian and Huánuco varieties within Erythroxylum coca Lam., and Colombian and Trujillo varieties within Erythroxylum novogranatense (D. Morris) Hieron.) but our understanding of the domestication and evolutionary history of these taxa is nominal. In this study, we use genomic data from natural history collections to estimate the geographic origins and genetic diversity of this economically and culturally important crop in the context of its wild relatives. Our phylogeographic analyses clearly demonstrate the four varieties of coca comprise two or three exclusive groups nested within the diverse lineages of the widespread, wild species Erythroxylum gracilipes; establishing a new and robust hypothesis of domestication wherein coca originated two or three times from this wild progenitor. The Colombian and Trujillo coca varieties are descended from a single, ancient domestication event in northwestern South America. Huánuco coca was domesticated more recently, possibly in southeastern Peru. Amazonian coca either shares a common domesticated ancestor with Huánuco coca, or it was the product of a third and most recent independent domestication event in the western Amazon basin. This chronology of coca domestication reveals different Holocene peoples in South America were able to independently transform the same natural resource to serve their needs; in this case, a workaday stimulant. [Erythroxylum; Erythroxylaceae; Holocene; Museomics; Neotropics; phylogeography; plant domestication; target-sequence capture.]

Phylogeographic analysis and species distribution modelling of the wood frog Batrachyla leptopus (Batrachylidae) reveal interglacial diversification in south western Patagonia

Background

The evolutionary history of southern South American organisms has been strongly influenced by Pleistocene climate oscillations. Amphibians are good models to evaluate hypotheses about the influence of these climate cycles on population structure and diversification of the biota, because they are sensitive to environmental changes and have restricted dispersal capabilities. We test hypotheses regarding putative forest refugia and expansion events associated with past climatic changes in the wood frog Batrachyla leptopus distributed along ∼1,000 km of length including glaciated and non-glaciated areas in southwestern Patagonia.

Methods

Using three mitochondrial regions (D-loop, cyt b, and coI) and two nuclear loci (pomc and crybA1), we conducted multilocus phylogeographic analyses and species distribution modelling to gain insights of the evolutionary history of this species. Intraspecific genealogy was explored with maximum likelihood, Bayesian, and phylogenetic network approaches. Diversification time was assessed using molecular clock models in a Bayesian framework, and demographic scenarios were evaluated using approximate Bayesian computation (ABC) and extended Bayesian skyline plot (EBSP). Species distribution models (SDM) were reconstructed using climatic and geographic data.

Results

Population structure and genealogical analyses support the existence of four lineages distributed north to south, with moderate to high phylogenetic support (Bootstrap > 70%; BPP > 0.92). The diversification time of B. leptopus’ populations began at ∼0.107 mya. The divergence between A and B lineages would have occurred by the late Pleistocene, approximately 0.068 mya, and divergence between C and D lineages was approximately 0.065 mya. The ABC simulations indicate that lineages coalesced at two different time periods, suggesting the presence of at least two glacial refugia and a postglacial colonization route that may have generated two southern lineages (p = 0.93, type I error: <0.094, type II error: 0.134). EBSP, mismatch distribution and neutrality indexes suggest sudden population expansion at ∼0.02 mya for all lineages. SDM infers fragmented distributions of B. leptopus associated with Pleistocene glaciations. Although the present populations of B. leptopus are found in zones affected by the last glacial maximum (∼0.023 mya), our analyses recover an older history of interglacial diversification (0.107–0.019 mya). In addition, we hypothesize two glacial refugia and three interglacial colonization routes, one of which gave rise to two expanding lineages in the south.

Inferring the Effect of Species Interactions on Trait Evolution

Models of trait evolution form an important part of macroevolutionary biology. The Brownian motion model and Ornstein–Uhlenbeck models have become classic (null) models of character evolution, in which species evolve independently. Recently, models incorporating species interactions have been developed, particularly involving competition where abiotic factors pull species toward an optimal trait value and competitive interactions drive the trait values apart. However, these models assume a fitness function rather than derive it from population dynamics and they do not consider dynamics of the trait variance. Here, we develop a general coherent trait evolution framework where the fitness function is based on a model of population dynamics, and therefore it can, in principle, accommodate any type of species interaction. We illustrate our framework with a model of abundance-dependent competitive interactions against a macroevolutionary background encoded in a phylogenetic tree. We develop an inference tool based on Approximate Bayesian Computation and test it on simulated data (of traits at the tips). We find that inference performs well when the diversity predicted by the parameters equals the number of species in the phylogeny. We then fit the model to empirical data of baleen whale body lengths, using three different summary statistics, and compare it to a model without population dynamics and a model where competition depends on the total metabolic rate of the competitors. We show that the unweighted model performs best for the least informative summary statistic, while the model with competition weighted by the total metabolic rate fits the data slightly better than the other two models for the two more informative summary statistics. Regardless of the summary statistic used, the three models substantially differ in their predictions of the abundance distribution. Therefore, data on abundance distributions will allow us to better distinguish the models from one another, and infer the nature of species interactions. Thus, our framework provides a conceptual approach to reveal species interactions underlying trait evolution and identifies the data needed to do so in practice. [Approximate Bayesian computation; competition; phylogeny; population dynamics; simulations; species interaction; trait evolution.]

CaliPro: A Calibration Protocol That Utilizes Parameter Density Estimation to Explore Parameter Space and Calibrate Complex Biological Models

Introduction

Mathematical and computational modeling have a long history of uncovering mechanisms and making predictions for biological systems. However, to create a model that can provide relevant quantitative predictions, models must first be calibrated by recapitulating existing biological datasets from that system. Current calibration approaches may not be appropriate for complex biological models because: 1) many attempt to recapitulate only a single aspect of the experimental data (such as a median trend) or 2) Bayesian techniques require specification of parameter priors and likelihoods to experimental data that cannot always be confidently assigned. A new calibration protocol is needed to calibrate complex models when current approaches fall short.

Methods

Herein, we develop CaliPro, an iterative, model-agnostic calibration protocol that utilizes parameter density estimation to refine parameter space and calibrate to temporal biological datasets. An important aspect of CaliPro is the user-defined pass set definition, which specifies how the model might successfully recapitulate experimental data. We define the appropriate settings to use CaliPro.

Results

We illustrate the usefulness of CaliPro through four examples including predator-prey, infectious disease transmission, and immune response models. We show that CaliPro works well for both deterministic, continuous model structures as well as stochastic, discrete models and illustrate that CaliPro can work across diverse calibration goals.

Conclusions

We present CaliPro, a new method for calibrating complex biological models to a range of experimental outcomes. In addition to expediting calibration, CaliPro may be useful in already calibrated parameter spaces to target and isolate specific model behavior for further analysis.

Estimating the establishment of local transmission and the cryptic phase of the COVID-19 pandemic in the USA

We use a global metapopulation transmission model to study the establishment of sustained and undetected community transmission of the COVID-19 pandemic in the United States. The model is calibrated on international case importations from mainland China and takes into account travel restrictions to and from international destinations. We estimate widespread community transmission of SARS-CoV-2 in February, 2020. Modeling results indicate international travel as the key driver of the introduction of SARS-CoV-2 in the West and East Coast metropolitan areas that could have been seeded as early as late-December, 2019. For most of the continental states the largest contribution of imported infections arrived through domestic travel flows.

Real-Time Forecasting of the COVID-19 Outbreak in Chinese Provinces: Machine Learning Approach Using Novel Digital Data and Estimates From Mechanistic Models

BACKGROUND

The inherent difficulty of identifying and monitoring emerging outbreaks caused by novel pathogens can lead to their rapid spread; and if left unchecked, they may become major public health threats to the planet. The ongoing coronavirus disease (COVID-19) outbreak, which has infected over 2,300,000 individuals and caused over 150,000 deaths, is an example of one of these catastrophic events.

OBJECTIVE

We present a timely and novel methodology that combines disease estimates from mechanistic models and digital traces, via interpretable machine learning methodologies, to reliably forecast COVID-19 activity in Chinese provinces in real time.

METHODS

Our method uses the following as inputs: (a) official health reports, (b) COVID-19-related internet search activity, (c) news media activity, and (d) daily forecasts of COVID-19 activity from a metapopulation mechanistic model. Our machine learning methodology uses a clustering technique that enables the exploitation of geospatial synchronicities of COVID-19 activity across Chinese provinces and a data augmentation technique to deal with the small number of historical disease observations characteristic of emerging outbreaks.

RESULTS

Our model is able to produce stable and accurate forecasts 2 days ahead of the current time and outperforms a collection of baseline models in 27 out of 32 Chinese provinces.

CONCLUSIONS

Our methodology could be easily extended to other geographies currently affected by COVID-19 to aid decision makers with monitoring and possibly prevention.

An automatic adaptive method to combine summary statistics in approximate Bayesian computation

To infer the parameters of mechanistic models with intractable likelihoods, techniques such as approximate Bayesian computation (ABC) are increasingly being adopted. One of the main disadvantages of ABC in practical situations, however, is that parameter inference must generally rely on summary statistics of the data. This is particularly the case for problems involving high-dimensional data, such as biological imaging experiments. However, some summary statistics contain more information about parameters of interest than others, and it is not always clear how to weight their contributions within the ABC framework. We address this problem by developing an automatic, adaptive algorithm that chooses weights for each summary statistic. Our algorithm aims to maximize the distance between the prior and the approximate posterior by automatically adapting the weights within the ABC distance function. Computationally, we use a nearest neighbour estimator of the distance between distributions. We justify the algorithm theoretically based on properties of the nearest neighbour distance estimator. To demonstrate the effectiveness of our algorithm, we apply it to a variety of test problems, including several stochastic models of biochemical reaction networks, and a spatial model of diffusion, and compare our results with existing algorithms.

A Synthetic Likelihood Solution to the Silent Synapse Estimation Problem

Functional features of synaptic populations are typically inferred from random electrophysiological sampling of small subsets of synapses. Are these samples unbiased? Here, we develop a biophysically constrained statistical framework to address this question and apply it to assess the performance of a widely used method based on a failure-rate analysis to quantify the occurrence of silent (AMPAR-lacking) synapses. We simulate this method in silico and find that it is characterized by strong and systematic biases, poor reliability, and weak statistical power. Key conclusions are validated by whole-cell recordings from hippocampal neurons. To address these shortcomings, we develop a simulator of the experimental protocol and use it to compute a synthetic likelihood. By maximizing the likelihood, we infer silent synapse fraction with no bias, low variance, and superior statistical power over alternatives. Together, this generalizable approach highlights how a simulator of experimental methodologies can substantially improve the estimation of physiological properties.

Phenotypic selection through cell death: stochastic modelling of O-6-methylguanine-DNA methyltransferase dynamics

Glioblastoma (GBM) is the most aggressive malignant primary brain tumour with a median overall survival of 15 months. To treat GBM, patients currently undergo a surgical resection followed by exposure to radiotherapy and concurrent and adjuvant temozolomide (TMZ) chemotherapy. However, this protocol often leads to treatment failure, with drug resistance being the main reason behind this. To date, many studies highlight the role of O-6-methylguanine-DNA methyltransferase (MGMT) in conferring drug resistance. The mechanism through which MGMT confers resistance is not well studied-particularly in terms of computational models. With only a few reasonable biological assumptions, we were able to show that even a minimal model of MGMT expression could robustly explain TMZ-mediated drug resistance. In particular, we showed that for a wide range of parameter values constrained by novel cell growth and viability assays, a model accounting for only stochastic gene expression of MGMT coupled with cell growth, division, partitioning and death was able to exhibit phenotypic selection of GBM cells expressing MGMT in response to TMZ. Furthermore, we found this selection allowed the cells to pass their acquired phenotypic resistance onto daughter cells in a stable manner (as long as TMZ is provided). This suggests that stochastic gene expression alone is enough to explain the development of chemotherapeutic resistance.

Reducing complexity and unidentifiability when modelling human atrial cells

Mathematical models of a cellular action potential (AP) in cardiac modelling have become increasingly complex, particularly in gating kinetics, which control the opening and closing of individual ion channel currents. As cardiac models advance towards use in personalized medicine to inform clinical decision-making, it is critical to understand the uncertainty hidden in parameter estimates from their calibration to experimental data. This study applies approximate Bayesian computation to re-calibrate the gating kinetics of four ion channels in two existing human atrial cell models to their original datasets, providing a measure of uncertainty and indication of potential issues with selecting a single unique value given the available experimental data. Two approaches are investigated to reduce the uncertainty present: re-calibrating the models to a more complete dataset and using a less complex formulation with fewer parameters to constrain. The re-calibrated models are inserted back into the full cell model to study the overall effect on the AP. The use of more complete datasets does not eliminate uncertainty present in parameter estimates. The less complex model, particularly for the fast sodium current, gave a better fit to experimental data alongside lower parameter uncertainty and improved computational speed. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

Reducing complexity and unidentifiability when modelling human atrial cells

Mathematical models of a cellular action potential (AP) in cardiac modelling have become increasingly complex, particularly in gating kinetics, which control the opening and closing of individual ion channel currents. As cardiac models advance towards use in personalized medicine to inform clinical decision-making, it is critical to understand the uncertainty hidden in parameter estimates from their calibration to experimental data. This study applies approximate Bayesian computation to re-calibrate the gating kinetics of four ion channels in two existing human atrial cell models to their original datasets, providing a measure of uncertainty and indication of potential issues with selecting a single unique value given the available experimental data. Two approaches are investigated to reduce the uncertainty present: re-calibrating the models to a more complete dataset and using a less complex formulation with fewer parameters to constrain. The re-calibrated models are inserted back into the full cell model to study the overall effect on the AP. The use of more complete datasets does not eliminate uncertainty present in parameter estimates. The less complex model, particularly for the fast sodium current, gave a better fit to experimental data alongside lower parameter uncertainty and improved computational speed. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

Probabilistic Updating of Structural Models for Damage Assessment Using Approximate Bayesian Computation

A novel probabilistic approach for model updating based on approximate Bayesian computation with subset simulation (ABC-SubSim) is proposed for damage assessment of structures using modal data. The ABC-SubSim is a likelihood-free Bayesian approach in which the explicit expression of likelihood function is avoided and the posterior samples of model parameters are obtained using the technique of subset simulation. The novel contributions of this paper are on three fronts: one is the introduction of some new stopping criteria to find an appropriate tolerance level for the metric used in the ABC-SubSim; the second one is the employment of a hybrid optimization scheme to find finer optimal values for the model parameters; and the last one is the adoption of an iterative approach to determine the optimal weighting factors related to the residuals of modal frequency and mode shape in the metric. The effectiveness of this approach is demonstrated using three illustrative examples.

Multilevel modelling as a tool to include variability and uncertainty in quantitative microbiology and risk assessment. Thermal inactivation of Listeria monocytogenes as proof of concept

Variability is inherent in biology and also substantial for microbial populations. In the context of food safety risk assessment, it refers to differences in the response of different bacterial strains (between-strain variability) and different cells (within-strain variability) to the same condition (e.g. inactivation treatment). However, its quantification based on empirical observations and its incorporation in predictive models is a challenge for both experimental design and (statistical) analysis. In this article we propose the use of multilevel models to quantify (different levels of) variability and uncertainty and include them in the predictions. As proof of concept, we analyse the microbial inactivation of Listeria monocytogenes to thermal treatments including different levels of variability (between-strain and within-strain) and uncertainty. The relationship between the microbial count and time was expressed using a (non-linear) Weibullian model. Moreover, we defined stochastic hypotheses to describe the different types of variation at the level of the kinetic parameters, as well as in the observations (microbial counts). The model parameters (kinetic parameters and variances) are estimated using Bayesian statistics. The multilevel approach was compared against an analogous, single-level model. The multilevel methodology shrinks extreme parameter estimates towards the mean according to uncertainty, thus mitigating overfitting. In addition, this approach enables to easily incorporate different levels of variation (between-strain and/or within-strain variability and/or uncertainty) in the predictions. On the other hand, multilevel (Bayesian) models are more complex to define, implement, analyse and communicate than single-level models. Nevertheless, their ability to incorporate different sources of variability in predictions make them very suitable for Quantitative Microbial Risk Assessment.

Genomic analysis of the natural history of attention-deficit/hyperactivity disorder using Neanderthal and ancient Homo sapiens samples

Attention-deficit/hyperactivity disorder (ADHD) is an impairing neurodevelopmental condition highly prevalent in current populations. Several hypotheses have been proposed to explain this paradox, mainly in the context of the Paleolithic versus Neolithic cultural shift but especially within the framework of the mismatch theory. This theory elaborates on how a particular trait once favoured in an ancient environment might become maladaptive upon environmental changes. However, given the lack of genomic data available for ADHD, these theories have not been empirically tested. We took advantage of the largest GWAS meta-analysis available for this disorder consisting of over 20,000 individuals diagnosed with ADHD and 35,000 controls, to assess the evolution of ADHD-associated alleles in European populations using archaic, ancient and modern human samples. We also included Approximate Bayesian computation coupled with deep learning analyses and singleton density scores to detect human adaptation. Our analyses indicate that ADHD-associated alleles are enriched in loss of function intolerant genes, supporting the role of selective pressures in this early-onset phenotype. Furthermore, we observed that the frequency of variants associated with ADHD has steadily decreased since Paleolithic times, particularly in Paleolithic European populations compared to samples from the Neolithic Fertile Crescent. We demonstrate this trend cannot be explained by African admixture nor Neanderthal introgression, since introgressed Neanderthal alleles are enriched in ADHD risk variants. All analyses performed support the presence of long-standing selective pressures acting against ADHD-associated alleles until recent times. Overall, our results are compatible with the mismatch theory for ADHD but suggest a much older time frame for the evolution of ADHD-associated alleles compared to previous hypotheses.

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

Motivated by the rapid spread of coronavirus disease 2019 (COVID-19) in mainland China, we use a global metapopulation disease transmission model to project the impact of travel limitations on the national and international spread of the epidemic. The model is calibrated on the basis of internationally reported cases and shows that, at the start of the travel ban from Wuhan on 23 January 2020, most Chinese cities had already received many infected travelers. The travel quarantine of Wuhan delayed the overall epidemic progression by only 3 to 5 days in mainland China but had a more marked effect on the international scale, where case importations were reduced by nearly 80% until mid-February. Modeling results also indicate that sustained 90% travel restrictions to and from mainland China only modestly affect the epidemic trajectory unless combined with a 50% or higher reduction of transmission in the community.

Flexible model selection for mechanistic network models

Network models are applied across many domains where data can be represented as a network. Two prominent paradigms for modelling networks are statistical models (probabilistic models for the observed network) and mechanistic models (models for network growth and/or evolution). Mechanistic models are better suited for incorporating domain knowledge, to study effects of interventions (such as changes to specific mechanisms) and to forward simulate, but they typically have intractable likelihoods. As such, and in a stark contrast to statistical models, there is a relative dearth of research on model selection for such models despite the otherwise large body of extant work. In this article, we propose a simulator-based procedure for mechanistic network model selection that borrows aspects from Approximate Bayesian Computation along with a means to quantify the uncertainty in the selected model. To select the most suitable network model, we consider and assess the performance of several learning algorithms, most notably the so-called Super Learner, which makes our framework less sensitive to the choice of a particular learning algorithm. Our approach takes advantage of the ease to forward simulate from mechanistic network models to circumvent their intractable likelihoods. The overall process is flexible and widely applicable. Our simulation results demonstrate the approach's ability to accurately discriminate between competing mechanistic models. Finally, we showcase our approach with a protein-protein interaction network model from the literature for yeast (Saccharomyces cerevisiae).

A practical guide to pseudo-marginal methods for computational inference in systems biology

For many stochastic models of interest in systems biology, such as those describing biochemical reaction networks, exact quantification of parameter uncertainty through statistical inference is intractable. Likelihood-free computational inference techniques enable parameter inference when the likelihood function for the model is intractable but the generation of many sample paths is feasible through stochastic simulation of the forward problem. The most common likelihood-free method in systems biology is approximate Bayesian computation that accepts parameters that result in low discrepancy between stochastic simulations and measured data. However, it can be difficult to assess how the accuracy of the resulting inferences are affected by the choice of acceptance threshold and discrepancy function. The pseudo-marginal approach is an alternative likelihood-free inference method that utilises a Monte Carlo estimate of the likelihood function. This approach has several advantages, particularly in the context of noisy, partially observed, time-course data typical in biochemical reaction network studies. Specifically, the pseudo-marginal approach facilitates exact inference and uncertainty quantification, and may be efficiently combined with particle filters for low variance, high-accuracy likelihood estimation. In this review, we provide a practical introduction to the pseudo-marginal approach using inference for biochemical reaction networks as a series of case studies. Implementations of key algorithms and examples are provided using the Julia programming language; a high performance, open source programming language for scientific computing (https://github.com/davidwarne/Warne2019_GuideToPseudoMarginal).

Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art

Stochasticity is a key characteristic of intracellular processes such as gene regulation and chemical signalling. Therefore, characterizing stochastic effects in biochemical systems is essential to understand the complex dynamics of living things. Mathematical idealizations of biochemically reacting systems must be able to capture stochastic phenomena. While robust theory exists to describe such stochastic models, the computational challenges in exploring these models can be a significant burden in practice since realistic models are analytically intractable. Determining the expected behaviour and variability of a stochastic biochemical reaction network requires many probabilistic simulations of its evolution. Using a biochemical reaction network model to assist in the interpretation of time-course data from a biological experiment is an even greater challenge due to the intractability of the likelihood function for determining observation probabilities. These computational challenges have been subjects of active research for over four decades. In this review, we present an accessible discussion of the major historical developments and state-of-the-art computational techniques relevant to simulation and inference problems for stochastic biochemical reaction network models. Detailed algorithms for particularly important methods are described and complemented with Matlab^® implementations. As a result, this review provides a practical and accessible introduction to computational methods for stochastic models within the life sciences community.

Inference of coevolutionary dynamics and parameters from host and parasite polymorphism data of repeated experiments

There is a long-standing interest in understanding host-parasite coevolutionary dynamics and associated fitness effects. Increasing amounts of genomic data for both interacting species offer a promising source to identify candidate loci and to infer the main parameters of the past coevolutionary history. However, so far no method exists to perform the latter. By coupling a gene-for-gene model with coalescent simulations, we first show that three types of biological costs, namely, resistance, infectivity and infection, define the allele frequencies at the internal equilibrium point of the coevolution model. These in return determine the strength of selective signatures at the coevolving host and parasite loci. We apply an Approximate Bayesian Computation (ABC) approach on simulated datasets to infer these costs by jointly integrating host and parasite polymorphism data at the coevolving loci. To control for the effect of genetic drift on coevolutionary dynamics, we assume that 10 or 30 repetitions are available from controlled experiments or several natural populations. We study two scenarios: 1) the cost of infection and population sizes (host and parasite) are unknown while costs of infectivity and resistance are known, and 2) all three costs are unknown while populations sizes are known. Using the ABC model choice procedure, we show that for both scenarios, we can distinguish with high accuracy pairs of coevolving host and parasite loci from pairs of neutrally evolving loci, though the statistical power decreases with higher cost of infection. The accuracy of parameter inference is high under both scenarios especially when using both host and parasite data because parasite polymorphism data do inform on costs applying to the host and vice-versa. As the false positive rate to detect pairs of genes under coevolution is small, we suggest that our method complements recently developed methods to identify host and parasite candidate loci for functional studies.

On the use of approximate Bayesian computation Markov chain Monte Carlo with inflated tolerance and post-correction

Approximate Bayesian computation enables inference for complicated probabilistic models with intractable likelihoods using model simulations. The Markov chain Monte Carlo implementation of approximate Bayesian computation is often sensitive to the tolerance parameter: low tolerance leads to poor mixing and large tolerance entails excess bias. We propose an approach that involves using a relatively large tolerance for the Markov chain Monte Carlo sampler to ensure sufficient mixing and post-processing the output, leading to estimators for a range of finer tolerances. We introduce an approximate confidence interval for the related post-corrected estimators and propose an adaptive approximate Bayesian computation Markov chain Monte Carlo algorithm, which finds a balanced tolerance level automatically based on acceptance rate optimization. Our experiments show that post-processing-based estimators can perform better than direct Markov chain Monte Carlo targeting a fine tolerance, that our confidence intervals are reliable, and that our adaptive algorithm leads to reliable inference with little user specification.

An ABC Method for Whole-Genome Sequence Data: Inferring Paleolithic and Neolithic Human Expansions

Species generally undergo a complex demographic history consisting, in particular, of multiple changes in population size. Genome-wide sequencing data are potentially highly informative for reconstructing this demographic history. A crucial point is to extract the relevant information from these very large data sets. Here, we design an approach for inferring past demographic events from a moderate number of fully sequenced genomes. Our new approach uses Approximate Bayesian Computation, a simulation-based statistical framework that allows 1) identifying the best demographic scenario among several competing scenarios and 2) estimating the best-fitting parameters under the chosen scenario. Approximate Bayesian Computation relies on the computation of summary statistics. Using a cross-validation approach, we show that statistics such as the lengths of haplotypes shared between individuals, or the decay of linkage disequilibrium with distance, can be combined with classical statistics (e.g., heterozygosity and Tajima's D) to accurately infer complex demographic scenarios including bottlenecks and expansion periods. We also demonstrate the importance of simultaneously estimating the genotyping error rate. Applying our method on genome-wide human-sequence databases, we finally show that a model consisting in a bottleneck followed by a Paleolithic and a Neolithic expansion is the most relevant for Eurasian populations.

TopoFilter: a MATLAB package for mechanistic model identification in systems biology

Background

To develop mechanistic dynamic models in systems biology, one often needs to identify all (or minimal) representations of the biological processes that are consistent with experimental data, out of a potentially large set of hypothetical mechanisms. However, a simple enumeration of all alternatives becomes quickly intractable when the number of model parameters grows. Selecting appropriate dynamic models out of a large ensemble of models, taking the uncertainty in our biological knowledge and in the experimental data into account, is therefore a key current problem in systems biology.

Results

The TopoFilter package addresses this problem in a heuristic and automated fashion by implementing the previously described topological filtering method for Bayesian model selection. It includes a core heuristic for searching the space of submodels of a parametrized model, coupled with a sampling-based exploration of the parameter space. Recent developments of the method allow to balance exhaustiveness and speed of the model space search, to efficiently re-sample parameters, to parallelize the search, and to use custom scoring functions. We use a theoretical example to motivate these features and then demonstrate TopoFilter’s applicability for a yeast signaling network with more than 250’000 possible model structures.

Conclusions

TopoFilter is a flexible software framework that makes Bayesian model selection and reduction efficient and scalable to network models of a complexity that represents contemporary problems in, for example, cell signaling. TopoFilter is open-source, available under the GPL-3.0 license at https://gitlab.com/csb.ethz/TopoFilter. It includes installation instructions, a quickstart guide, a description of all package options, and multiple examples.

Reverse engineering directed gene regulatory networks from transcriptomics and proteomics data of biomining bacterial communities with approximate Bayesian computation and steady-state signalling simulations

BACKGROUND

Network inference is an important aim of systems biology. It enables the transformation of OMICs datasets into biological knowledge. It consists of reverse engineering gene regulatory networks from OMICs data, such as RNAseq or mass spectrometry-based proteomics data, through computational methods. This approach allows to identify signalling pathways involved in specific biological functions. The ability to infer causality in gene regulatory networks, in addition to correlation, is crucial for several modelling approaches and allows targeted control in biotechnology applications.

METHODS

We performed simulations according to the approximate Bayesian computation method, where the core model consisted of a steady-state simulation algorithm used to study gene regulatory networks in systems for which a limited level of details is available. The simulations outcome was compared to experimentally measured transcriptomics and proteomics data through approximate Bayesian computation.

RESULTS

The structure of small gene regulatory networks responsible for the regulation of biological functions involved in biomining were inferred from multi OMICs data of mixed bacterial cultures. Several causal inter- and intraspecies interactions were inferred between genes coding for proteins involved in the biomining process, such as heavy metal transport, DNA damage, replication and repair, and membrane biogenesis. The method also provided indications for the role of several uncharacterized proteins by the inferred connection in their network context.

CONCLUSIONS

The combination of fast algorithms with high-performance computing allowed the simulation of a multitude of gene regulatory networks and their comparison to experimentally measured OMICs data through approximate Bayesian computation, enabling the probabilistic inference of causality in gene regulatory networks of a multispecies bacterial system involved in biomining without need of single-cell or multiple perturbation experiments. This information can be used to influence biological functions and control specific processes in biotechnology applications.

Supermodeling: The Next Level of Abstraction in the Use of Data Assimilation

Data assimilation (DA) is a key procedure that synchronizes a computer model with real observations. However, in the case of overparametrized complex systems modeling, the task of parameter-estimation through data assimilation can expand exponentially. It leads to unacceptable computational overhead, substantial inaccuracies in parameter matching, and wrong predictions. Here we define a Supermodel as a kind of ensembling scheme, which consists of a few sub-models representing various instances of the baseline model. The sub-models differ in parameter sets and are synchronized through couplings between the most sensitive dynamical variables. We demonstrate that after a short pretraining of the fully parametrized small sub-model ensemble, and then training a few latent parameters of the low-parameterized Supermodel, we can outperform in efficiency and accuracy the baseline model matched to data by a classical DA procedure.

Unknotting the interactive effects of learning processes on cultural evolutionary dynamics

Forms of non-random copying error provide sources of inherited variation yet their effects on cultural evolutionary dynamics are poorly understood. Focusing on variation in granny and reef knot forms, we present a mathematical model that specifies how these variant frequencies are affected by non-linear interactions between copying fidelity, mirroring, handedness and repetition biases. Experiments on adult humans allowed these effects to be estimated using approximate Bayesian computation and the model is iterated to explain the prevalence of granny over reef knots in the wild. Our study system also serves to show conditions under which copying fidelity drives heterogeneity in cultural variants at equilibrium, and that interaction between unbiased forms of copying error can skew cultural variation.

Estimation of the within-herd transmission rates of bovine viral diarrhoea virus in extensively grazed beef cattle herds

Many research groups have developed mathematical models to simulate the dynamics of BVDV infections in cattle herds. However, most models use estimates for within-herd BVDV transmission rates that are either based on expert opinion or adapted from other dairy herd simulation models presented in the literature. There is currently little information on the transmission rates for BVDV in extensively grazed beef herds partly due to the logistical challenges in obtaining longitudinal data of individual animal’s seroconversion, and it may not be appropriate to apply the same transmission rates from intensive dairy herds given the significant differences in herd demographics and management. To address this knowledge gap, we measured BVDV antibody levels in 15 replacement heifers in each of 75 New Zealand beef breeding farms after their first calving and again at pregnancy scanning or weaning to check for seroconversion. Among these, data from 9 farms were used to infer the within-herd BVDV transmission rate with an approximate Bayesian computation method. The most probable within-herd BVDV transmission rate was estimated as 0.11 per persistently infected (PI) animal per day with a 95% highest posterior density interval between 0.03 and 0.34. This suggests that BVDV transmission in extensively grazed beef herds is generally slower than in dairy herds where the transmission rate has been estimated at 0.50 per PI animal per day and therefore may not be sufficient to ensure that all susceptible breeding females gain adequate immunity to the virus before the risk period of early pregnancy for generating new PI calves.