We present a novel algorithm, implemented in the software ARGinfer, for probabilistic inference of the Ancestral Recombination Graph under the Coalescent with Recombination. Our Markov Chain Monte Carlo algorithm takes advantage of the Succinct Tree Sequence data structure that has allowed great advances in simulation and point estimation, but not yet probabilistic inference. Unlike previous methods, which employ the Sequentially Markov Coalescent approximation, ARGinfer uses the Coalescent with Recombination, allowing more accurate inference of key evolutionary parameters. We show using simulations that ARGinfer can accurately estimate many properties of the evolutionary history of the sample, including the topology and branch lengths of the genealogical tree at each sequence site, and the times and locations of mutation and recombination events. ARGinfer approximates posterior probability distributions for these and other quantities, providing interpretable assessments of uncertainty that we show to be well calibrated. ARGinfer is currently limited to tens of DNA sequences of several hundreds of kilobases, but has scope for further computational improvements to increase its applicability.

Gene set enrichment tests (a.k.a. functional enrichment analysis) are among the most frequently used methods in computational biology. Despite this popularity, there are concerns that these methods are being applied incorrectly and the results of some peer-reviewed publications are unreliable. These problems include the use of inappropriate background gene lists, lack of false discovery rate correction and lack of methodological detail. To ascertain the frequency of these issues in the literature, we performed a screen of 186 open-access research articles describing functional enrichment results. We find that 95% of analyses using over-representation tests did not implement an appropriate background gene list or did not describe this in the methods. Failure to perform p-value correction for multiple tests was identified in 43% of analyses. Many studies lacked detail in the methods section about the tools and gene sets used. An extension of this survey showed that these problems are not associated with journal or article level bibliometrics. Using seven independent RNA-seq datasets, we show misuse of enrichment tools alters results substantially. In conclusion, most published functional enrichment studies suffered from one or more major flaws, highlighting the need for stronger standards for enrichment analysis.

A long-standing practice in the treatment of cancer is that of hitting hard with the maximum tolerated dose to eradicate tumors. This continuous therapy, however, selects for resistant cells, leading to the failure of the treatment. A different type of treatment strategy, adaptive therapy, has recently been shown to have a degree of success in both preclinical xenograft experiments and clinical trials. Adaptive therapy is used to maintain a tumor's volume by exploiting the competition between drug-sensitive and drug-resistant cells with minimum effective drug doses or timed drug holidays. To further understand the role of competition in the outcomes of adaptive therapy, we developed a 2D on-lattice agent-based model. Our simulations show that the superiority of the adaptive strategy over continuous therapy depends on the local competition shaped by the spatial distribution of resistant cells. Intratumor competition can also be affected by fibroblasts, which produce microenvironmental factors that promote cancer cell growth. To this end, we simulated the impact of different fibroblast distributions on treatment outcomes. As a proof of principle, we focused on five types of distribution of fibroblasts characterized by different locations, shapes, and orientations of the fibroblast region with respect to the resistant cells. Our simulation shows that the spatial architecture of fibroblasts modulates tumor progression in both continuous and adaptive therapy. Finally, as a proof of concept, we simulated the outcomes of adaptive therapy of a virtual patient with four metastatic sites composed of different spatial distributions of fibroblasts and drug-resistant cell populations. Our simulation highlights the importance of undetected metastatic lesions on adaptive therapy outcomes.

Performing a cognitive task requires going through a sequence of functionally diverse stages. Although it is typically assumed that these stages are characterized by distinct states of cortical synchrony that are triggered by sub-cortical events, little reported evidence supports this hypothesis. To test this hypothesis, we first identified cognitive stages in single-trial MEG data of an associative recognition task, showing with a novel method that each stage begins with local modulations of synchrony followed by a state of directed functional connectivity. Second, we developed the first whole-brain model that can simulate cortical synchrony throughout a task. The model suggests that the observed synchrony is caused by thalamocortical bursts at the onset of each stage, targeted at cortical synapses and interacting with the structural anatomical connectivity. These findings confirm that cognitive stages are defined by distinct states of cortical synchrony and explains the network-level mechanisms necessary for reaching stage-dependent synchrony states.

Transcription factors (TFs) play an important role in regulating gene expression, thus the identification of the sites bound by them has become a fundamental step for molecular and cellular biology. In this paper, we developed a deep learning framework leveraging existing fully convolutional neural networks (FCN) to predict TF-DNA binding signals at the base-resolution level (named as FCNsignal). The proposed FCNsignal can simultaneously achieve the following tasks: (i) modeling the base-resolution signals of binding regions; (ii) discriminating binding or non-binding regions; (iii) locating TF-DNA binding regions; (iv) predicting binding motifs. Besides, FCNsignal can also be used to predict opening regions across the whole genome. The experimental results on 53 TF ChIP-seq datasets and 6 chromatin accessibility ATAC-seq datasets show that our proposed framework outperforms some existing state-of-the-art methods. In addition, we explored to use the trained FCNsignal to locate all potential TF-DNA binding regions on a whole chromosome and predict DNA sequences of arbitrary length, and the results show that our framework can find most of the known binding regions and accept sequences of arbitrary length. Furthermore, we demonstrated the potential ability of our framework in discovering causal disease-associated single-nucleotide polymorphisms (SNPs) through a series of experiments.

Although estrogen-receptor-positive (ER+) breast cancer is generally associated with favorable prognosis, clinical outcome varies substantially among patients. Genomic assays have been developed and applied to predict patient prognosis for personalized treatment. We hypothesize that the recurrence risk of ER+ breast cancer patients is determined by both genomic mutations intrinsic to tumor cells and extrinsic immunological features in the tumor microenvironment. Based on the Cancer Genome Atlas (TCGA) breast cancer data, we identified the 72 most common genomic aberrations (including gene mutations and indels) in ER+ breast cancer and defined sample-specific scores that systematically characterized the deregulated pathways intrinsic to tumor cells. To further consider tumor cell extrinsic features, we calculated immune infiltration scores for six major immune cell types. Many individual intrinsic features are predictive of patient prognosis in ER+ breast cancer, and some of them achieved comparable accuracy with the Oncotype DX assay. In addition, statistical learning models that integrated these features predicts the recurrence risk of patients with significantly better performance than the Oncotype DX assay (our optimized random forest model AUC = 0.841, Oncotype DX model AUC = 0.792, p = 0.04). As a proof-of-concept, our study indicates the great potential of genomic and immunological features in prognostic prediction for improving breast cancer precision medicine. The framework introduced in this work can be readily applied to other cancers.

Perception is often characterized computationally as an inference process in which uncertain or ambiguous sensory inputs are combined with prior expectations. Although behavioral studies have shown that observers can change their prior expectations in the context of a task, robust neural signatures of task-specific priors have been elusive. Here, we analytically derive such signatures under the general assumption that the responses of sensory neurons encode posterior beliefs that combine sensory inputs with task-specific expectations. Specifically, we derive predictions for the task-dependence of correlated neural variability and decision-related signals in sensory neurons. The qualitative aspects of our results are parameter-free and specific to the statistics of each task. The predictions for correlated variability also differ from predictions of classic feedforward models of sensory processing and are therefore a strong test of theories of hierarchical Bayesian inference in the brain. Importantly, we find that Bayesian learning predicts an increase in so-called "differential correlations" as the observer's internal model learns the stimulus distribution, and the observer's behavioral performance improves. This stands in contrast to classic feedforward encoding/decoding models of sensory processing, since such correlations are fundamentally information-limiting. We find support for our predictions in data from existing neurophysiological studies across a variety of tasks and brain areas. Finally, we show in simulation how measurements of sensory neural responses can reveal information about a subject's internal beliefs about the task. Taken together, our results reinterpret task-dependent sources of neural covariability as signatures of Bayesian inference and provide new insights into their cause and their function.

A central goal in sensory neuroscience is to understand the neuronal signal processing involved in the encoding of natural stimuli. A critical step towards this goal is the development of successful computational encoding models. For ganglion cells in the vertebrate retina, the development of satisfactory models for responses to natural visual scenes is an ongoing challenge. Standard models typically apply linear integration of visual stimuli over space, yet many ganglion cells are known to show nonlinear spatial integration, in particular when stimulated with contrast-reversing gratings. We here study the influence of spatial nonlinearities in the encoding of natural images by ganglion cells, using multielectrode-array recordings from isolated salamander and mouse retinas. We assess how responses to natural images depend on first- and second-order statistics of spatial patterns inside the receptive field. This leads us to a simple extension of current standard ganglion cell models. We show that taking not only the weighted average of light intensity inside the receptive field into account but also its variance over space can partly account for nonlinear integration and substantially improve response predictions of responses to novel images. For salamander ganglion cells, we find that response predictions for cell classes with large receptive fields profit most from including spatial contrast information. Finally, we demonstrate how this model framework can be used to assess the spatial scale of nonlinear integration. Our results underscore that nonlinear spatial stimulus integration translates to stimulation with natural images. Furthermore, the introduced model framework provides a simple, yet powerful extension of standard models and may serve as a benchmark for the development of more detailed models of the nonlinear structure of receptive fields.

Mouse is the most used model for studying the impact of microbiota on its host, but the repertoire of species from the mouse gut microbiome remains largely unknown. Accordingly, the similarity between human and mouse microbiomes at a low taxonomic level is not clear. We construct a comprehensive mouse microbiota genome (CMMG) catalog by assembling all currently available mouse gut metagenomes and combining them with published reference and metagenome-assembled genomes. The 41'798 genomes cluster into 1'573 species, of which 78.1% are uncultured, and we discovered 226 new genera, seven new families, and one new order. CMMG enables an unprecedented coverage of the mouse gut microbiome exceeding 86%, increases the mapping rate over four-fold, and allows functional microbiota analyses of human and mouse linking them to the driver species. Comparing CMMG to microbiota from the unified human gastrointestinal genomes shows an overlap of 62% at the genus but only 10% at the species level, demonstrating that human and mouse gut microbiota are largely distinct. CMMG contains the most comprehensive collection of consistently functionally annotated species of the mouse and human microbiome to date, setting the ground for analysis of new and reanalysis of existing datasets at an unprecedented depth.

One of the most difficult sensorimotor behaviors exhibited by flying animals is the ability to track another flying animal based on its sound emissions. From insects to mammals, animals display this ability in order to localize and track conspecifics, mate or prey. The pursuing individual must overcome multiple non-trivial challenges including the detection of the sounds emitted by the target, matching the input received by its (mostly) two sensors, localizing the direction of the sound target in real time and then pursuing it. All this has to be done rapidly as the target is constantly moving. In this project, we set to mimic this ability using a physical bio-mimetic autonomous drone. We equipped a miniature commercial drone with our in-house 2D sound localization electronic circuit which uses two microphones (mimicking biological ears) to localize sound signals in real-time and steer the drone in the horizontal plane accordingly. We focus on bat signals because bats are known to eavesdrop on conspecifics and follow them, but our approach could be generalized to other biological signals and other man-made signals. Using two different experiments, we show that our fully autonomous aviator can track the position of a moving sound emitting target and pursue it in real-time. Building an actual robotic-agent, forced us to deal with real-life difficulties which also challenge animals. We thus discuss the similarities and differences between our and the biological approach.

High-throughput experimentation has revolutionized data-driven experimental sciences and opened the door to the application of machine learning techniques. Nevertheless, the quality of any data analysis strongly depends on the quality of the data and specifically the degree to which random effects in the experimental data-generating process are quantified and accounted for. Accordingly calibration, i.e. the quantitative association between observed quantities and measurement responses, is a core element of many workflows in experimental sciences. Particularly in life sciences, univariate calibration, often involving non-linear saturation effects, must be performed to extract quantitative information from measured data. At the same time, the estimation of uncertainty is inseparably connected to quantitative experimentation. Adequate calibration models that describe not only the input/output relationship in a measurement system but also its inherent measurement noise are required. Due to its mathematical nature, statistically robust calibration modeling remains a challenge for many practitioners, at the same time being extremely beneficial for machine learning applications. In this work, we present a bottom-up conceptual and computational approach that solves many problems of understanding and implementing non-linear, empirical calibration modeling for quantification of analytes and process modeling. The methodology is first applied to the optical measurement of biomass concentrations in a high-throughput cultivation system, then to the quantification of glucose by an automated enzymatic assay. We implemented the conceptual framework in two Python packages, calibr8 and murefi, with which we demonstrate how to make uncertainty quantification for various calibration tasks more accessible. Our software packages enable more reproducible and automatable data analysis routines compared to commonly observed workflows in life sciences. Subsequently, we combine the previously established calibration models with a hierarchical Monod-like ordinary differential equation model of microbial growth to describe multiple replicates of Corynebacterium glutamicum batch cultures. Key process model parameters are learned by both maximum likelihood estimation and Bayesian inference, highlighting the flexibility of the statistical and computational framework.

Biological sequence families contain many sequences that are very similar to each other because they are related by evolution, so the strategy for splitting data into separate training and test sets is a nontrivial choice in benchmarking sequence analysis methods. A random split is insufficient because it will yield test sequences that are closely related or even identical to training sequences. Adapting ideas from independent set graph algorithms, we describe two new methods for splitting sequence data into dissimilar training and test sets. These algorithms input a sequence family and produce a split in which each test sequence is less than p% identical to any individual training sequence. These algorithms successfully split more families than a previous approach, enabling construction of more diverse benchmark datasets.

Emerging clinical evidence suggests that thrombosis in the microvasculature of patients with Coronavirus disease 2019 (COVID-19) plays an essential role in dictating the disease progression. Because of the infectious nature of SARS-CoV-2, patients' fresh blood samples are limited to access for in vitro experimental investigations. Herein, we employ a novel multiscale and multiphysics computational framework to perform predictive modeling of the pathological thrombus formation in the microvasculature using data from patients with COVID-19. This framework seamlessly integrates the key components in the process of blood clotting, including hemodynamics, transport of coagulation factors and coagulation kinetics, blood cell mechanics and adhesive dynamics, and thus allows us to quantify the contributions of many prothrombotic factors reported in the literature, such as stasis, the derangement in blood coagulation factor levels and activities, inflammatory responses of endothelial cells and leukocytes to the microthrombus formation in COVID-19. Our simulation results show that among the coagulation factors considered, antithrombin and factor V play more prominent roles in promoting thrombosis. Our simulations also suggest that recruitment of WBCs to the endothelial cells exacerbates thrombogenesis and contributes to the blockage of the blood flow. Additionally, we show that the recent identification of flowing blood cell clusters could be a result of detachment of WBCs from thrombogenic sites, which may serve as a nidus for new clot formation. These findings point to potential targets that should be further evaluated, and prioritized in the anti-thrombotic treatment of patients with COVID-19. Altogether, our computational framework provides a powerful tool for quantitative understanding of the mechanism of pathological thrombus formation and offers insights into new therapeutic approaches for treating COVID-19 associated thrombosis.

Circadian clocks tick a rhythm with a nearly 24-hour period in a variety of organisms. In the clock proteins of cyanobacteria, KaiA, KaiB, and KaiC, known as a minimum circadian clock, the slow KaiB-KaiC complex formation is essential in determining the clock period. This complex formation, occurring when the C1 domain of KaiC hexamer binds ADP molecules produced by the ATPase activity of C1, is considered to be promoted by accumulating ADP molecules in C1 through inhibiting the ADP/ATP exchange (ADP release) rather than activating the ATP hydrolysis (ADP production). Significantly, this ADP/ATP exchange inhibition accelerates the complex formation together with its promotion, implying a potential role in the period robustness under environmental perturbations. However, the molecular mechanism of this simultaneous promotion and acceleration remains elusive because inhibition of a backward process generally slows down the whole process. In this article, to investigate the mechanism, we build several reaction models of the complex formation with the pre-binding process concerning the ATPase activity. In these models, six KaiB monomers cooperatively and rapidly bind to C1 when C1 binds ADP molecules more than a given threshold while stabilizing the binding-competent conformation of C1. Through comparison among the models proposed here, we then extract three requirements for the simultaneous promotion and acceleration: the stabilization of the binding-competent C1 by KaiB binding, slow ADP/ATP exchange in the binding-competent C1, and relatively fast ADP/ATP exchange occurring in the binding-incompetent C1 in the presence of KaiB. The last two requirements oblige KaiC to form a multimer. Moreover, as a natural consequence, the present models can also explain why the binding of KaiB to C1 reduces the ATPase activity of C1.

The understanding of bacterial gene function has been greatly enhanced by recent advancements in the deep sequencing of microbial genomes. Transposon insertion sequencing methods combines next-generation sequencing techniques with transposon mutagenesis for the exploration of the essentiality of genes under different environmental conditions. We propose a model-based method that uses regularized negative binomial regression to estimate the change in transposon insertions attributable to gene-environment changes in this genetic interaction study without transformations or uniform normalization. An empirical Bayes model for estimating the local false discovery rate combines unique and total count information to test for genes that show a statistically significant change in transposon counts. When applied to RB-TnSeq (randomized barcode transposon sequencing) and Tn-seq (transposon sequencing) libraries made in strains of Caulobacter crescentus using both total and unique count data the model was able to identify a set of conditionally beneficial or conditionally detrimental genes for each target condition that shed light on their functions and roles during various stress conditions.

Synchronization of neural oscillations is thought to facilitate communication in the brain. Neurodegenerative pathologies such as Parkinson's disease (PD) can result in synaptic reorganization of the motor circuit, leading to altered neuronal dynamics and impaired neural communication. Treatments for PD aim to restore network function via pharmacological means such as dopamine replacement, or by suppressing pathological oscillations with deep brain stimulation. We tested the hypothesis that brain stimulation can operate beyond a simple "reversible lesion" effect to augment network communication. Specifically, we examined the modulation of beta band (14-30 Hz) activity, a known biomarker of motor deficits and potential control signal for stimulation in Parkinson's. To do this we setup a neural mass model of population activity within the cortico-basal ganglia-thalamic (CBGT) circuit with parameters that were constrained to yield spectral features comparable to those in experimental Parkinsonism. We modulated the connectivity of two major pathways known to be disrupted in PD and constructed statistical summaries of the spectra and functional connectivity of the resulting spontaneous activity. These were then used to assess the network-wide outcomes of closed-loop stimulation delivered to motor cortex and phase locked to subthalamic beta activity. Our results demonstrate that the spatial pattern of beta synchrony is dependent upon the strength of inputs to the STN. Precisely timed stimulation has the capacity to recover network states, with stimulation phase inducing activity with distinct spectral and spatial properties. These results provide a theoretical basis for the design of the next-generation brain stimulators that aim to restore neural communication in disease.

Most microbes live in spatially structured communities (e.g., biofilms) in which they interact with their neighbors through the local exchange of diffusible molecules. To understand the functioning of these communities, it is essential to uncover how these local interactions shape community-level properties, such as the community composition, spatial arrangement, and growth rate. Here, we present a mathematical framework to derive community-level properties from the molecular mechanisms underlying the cell-cell interactions for systems consisting of two cell types. Our framework consists of two parts: a biophysical model to derive the local interaction rules (i.e. interaction range and strength) from the molecular parameters underlying the cell-cell interactions and a graph based model to derive the equilibrium properties of the community (i.e. composition, spatial arrangement, and growth rate) from these local interaction rules. Our framework shows that key molecular parameters underlying the cell-cell interactions (e.g., the uptake and leakage rates of molecules) determine community-level properties. We apply our model to mutualistic cross-feeding communities and show that spatial structure can be detrimental for these communities. Moreover, our model can qualitatively recapitulate the properties of an experimental microbial community. Our framework can be extended to a variety of systems of two interacting cell types, within and beyond the microbial world, and contributes to our understanding of how community-level properties emerge from microscopic interactions between cells.

microRNAs (miRNAs) are (18-22nt long) noncoding short (s)RNAs that suppress gene expression by targeting the 3' untranslated region of target mRNAs. This occurs through the seed sequence located in position 2-7/8 of the miRNA guide strand, once it is loaded into the RNA induced silencing complex (RISC). G-rich 6mer seed sequences can kill cells by targeting C-rich 6mer seed matches located in genes that are critical for cell survival. This results in induction of Death Induced by Survival gene Elimination (DISE), through a mechanism we have called 6mer seed toxicity. miRNAs are often quantified in cells by aligning the reads from small (sm)RNA sequencing to the genome. However, the analysis of any smRNA Seq data set for predicted 6mer seed toxicity requires an alternative workflow, solely based on the exact position 2-7 of any short (s)RNA that can enter the RISC. Therefore, we developed SPOROS, a semi-automated pipeline that produces multiple useful outputs to predict and compare 6mer seed toxicity of cellular sRNAs, regardless of their nature, between different samples. We provide two examples to illustrate the capabilities of SPOROS: Example one involves the analysis of RISC-bound sRNAs in a cancer cell line (either wild-type or two mutant lines unable to produce most miRNAs). Example two is based on a publicly available smRNA Seq data set from postmortem brains (either from normal or Alzheimer's patients). Our methods (found at https://github.com/ebartom/SPOROS and at Code Ocean: https://doi.org/10.24433/CO.1732496.v1) are designed to be used to analyze a variety of smRNA Seq data in various normal and disease settings.

Update of medRxiv. 2021 Jan 26;:

While acquired chemoresistance is recognized as a key challenge to treating many types of cancer, the dynamics with which drug sensitivity changes after exposure are poorly characterized. Most chemotherapeutic regimens call for repeated dosing at regular intervals, and if drug sensitivity changes on a similar time scale then the treatment interval could be optimized to improve treatment performance. Theoretical work suggests that such optimal schedules exist, but experimental confirmation has been obstructed by the difficulty of deconvolving the simultaneous processes of death, adaptation, and regrowth taking place in cancer cell populations. Here we present a method of optimizing drug schedules in vitro through iterative application of experimentally calibrated models, and demonstrate its ability to characterize dynamic changes in sensitivity to the chemotherapeutic doxorubicin in three breast cancer cell lines subjected to treatment schedules varying in concentration, interval between pulse treatments, and number of sequential pulse treatments. Cell populations are monitored longitudinally through automated imaging for 600-800 hours, and this data is used to calibrate a family of cancer growth models, each consisting of a system of ordinary differential equations, derived from the bi-exponential model which characterizes resistant and sensitive subpopulations. We identify a model incorporating both a period of growth arrest in surviving cells and a delay in the death of chemosensitive cells which outperforms the original bi-exponential growth model in Akaike Information Criterion based model selection, and use the calibrated model to quantify the performance of each drug schedule. We find that the inter-treatment interval is a key variable in determining the performance of sequential dosing schedules and identify an optimal retreatment time for each cell line which extends regrowth time by 40%-239%, demonstrating that the time scale of changes in chemosensitivity following doxorubicin exposure allows optimization of drug scheduling by varying this inter-treatment interval.

Beyond natural stimuli such as growth factors and stresses, the ability to experimentally modulate at will the levels or activity of specific intracellular signaling molecule(s) in specified cells within a tissue can be a powerful tool for uncovering new regulation and tissue behaviors. Here we perturb the levels of cAMP within specific cells of an epithelial monolayer to probe the time-dynamic behavior of cell-cell communication protocols implemented by the cAMP/PKA pathway and its coupling to the ERK pathway. The time-dependent ERK responses we observe in the perturbed cells for spatially uniform cAMP perturbations (all cells) can be very different from those due to spatially localized perturbations (a few cells). Through a combination of pharmacological and genetic perturbations, signal analysis, and computational modeling, we infer how intracellular regulation and regulated cell-cell coupling each impact the intracellular ERK response in single cells. Our approach reveals how a dynamic gap junction state helps sculpt the intracellular ERK response over time in locally perturbed cells.

The actin cortex is an active adaptive material, embedded with complex regulatory networks that can sense, generate, and transmit mechanical forces. The cortex exhibits a wide range of dynamic behaviours, from generating pulsatory contractions and travelling waves to forming organised structures. Despite the progress in characterising the biochemical and mechanical components of the actin cortex, the emergent dynamics of this mechanochemical system is poorly understood. Here we develop a reaction-diffusion model for the RhoA signalling network, the upstream regulator for actomyosin assembly and contractility, coupled to an active actomyosin gel, to investigate how the interplay between chemical signalling and mechanical forces regulates stresses and patterns in the cortex. We demonstrate that mechanochemical feedback in the cortex acts to destabilise homogeneous states and robustly generate pulsatile contractions. By tuning active stress in the system, we show that the cortex can generate propagating contraction pulses, form network structures, or exhibit topological turbulence.

Several important aspects related to SARS-CoV-2 transmission are not well known due to a lack of appropriate data. However, mathematical and computational tools can be used to extract part of this information from the available data, like some hidden age-related characteristics. In this paper, we present a method to investigate age-specific differences in transmission parameters related to susceptibility to and infectiousness upon contracting SARS-CoV-2 infection. More specifically, we use panel-based social contact data from diary-based surveys conducted in Belgium combined with the next generation principle to infer the relative incidence and we compare this to real-life incidence data. Comparing these two allows for the estimation of age-specific transmission parameters. Our analysis implies the susceptibility in children to be around half of the susceptibility in adults, and even lower for very young children (preschooler). However, the probability of adults and the elderly to contract the infection is decreasing throughout the vaccination campaign, thereby modifying the picture over time.

Developing measures for rapid and early detection of disease re-emergence is important to perform science-based risk assessment of epidemic threats. In the past few years, several early warning signals (EWS) from complex systems theory have been introduced to detect impending critical transitions and extend the set of indicators. However, it is still debated whether they are generically applicable or potentially sensitive to some dynamical characteristics such as system noise and rates of approach to critical parameter values. Moreover, testing on empirical data has, so far, been limited. Hence, verifying EWS performance remains a challenge. In this study, we tackle this question by analyzing the performance of common EWS, such as increasing variance and autocorrelation, in detecting the emergence of COVID-19 outbreaks in various countries. Our work illustrates that these EWS might be successful in detecting disease emergence when some basic assumptions are satisfied: a slow forcing through the transitions and not-fat-tailed noise. In uncertain cases, we observe that noise properties or commensurable time scales may obscure the expected early warning signals. Overall, our results suggest that EWS can be useful for active monitoring of epidemic dynamics, but that their performance is sensitive to certain features of the underlying dynamics. Our findings thus pave a connection between theoretical and empirical studies, constituting a further step towards the application of EWS indicators for informing public health policies.

One of the unique traits of membrane proteins is that a significant fraction of their hydrophobic amino acids is exposed to the hydrophobic core of lipid bilayers rather than being embedded in the protein interior, which is often not explicitly considered in the protein structure and function predictions. Here, we propose a characteristic and predictive quantity, the membrane contact probability (MCP), to describe the likelihood of the amino acids of a given sequence being in direct contact with the acyl chains of lipid molecules. We show that MCP is complementary to solvent accessibility in characterizing the outer surface of membrane proteins, and it can be predicted for any given sequence with a machine learning-based method by utilizing a training dataset extracted from MemProtMD, a database generated from molecular dynamics simulations for the membrane proteins with a known structure. As the first of many potential applications, we demonstrate that MCP can be used to systematically improve the prediction precision of the protein contact maps and structures.

Estimates of correlation between pairs of genes in co-expression analysis are commonly used to construct networks among genes using gene expression data. As previously noted, the distribution of such correlations depends on the observed expression level of the involved genes, which we refer to this as a mean-correlation relationship in RNA-seq data, both bulk and single-cell. This dependence introduces an unwanted technical bias in co-expression analysis whereby highly expressed genes are more likely to be highly correlated. Such a relationship is not observed in protein-protein interaction data, suggesting that it is not reflecting biology. Ignoring this bias can lead to missing potentially biologically relevant pairs of genes that are lowly expressed, such as transcription factors. To address this problem, we introduce spatial quantile normalization (SpQN), a method for normalizing local distributions in a correlation matrix. We show that spatial quantile normalization removes the mean-correlation relationship and corrects the expression bias in network reconstruction.

In many human cancers, the rate of cell growth depends crucially on the size of the tumor cell population. Low, zero, or negative growth at low population densities is known as the Allee effect; this effect has been studied extensively in ecology, but so far lacks a good explanation in the cancer setting. Here, we formulate and analyze an individual-based model of cancer, in which cell division rates are increased by the local concentration of an autocrine growth factor produced by the cancer cells themselves. We show, analytically and by simulation, that autocrine signaling suffices to cause both strong and weak Allee effects. Whether low cell densities lead to negative (strong effect) or reduced (weak effect) growth rate depends directly on the ratio of cell death to proliferation, and indirectly on cellular dispersal. Our model is consistent with experimental observations from three patient-derived brain tumor cell lines grown at different densities. We propose that further studying and quantifying population-wide feedback, impacting cell growth, will be central for advancing our understanding of cancer dynamics and treatment, potentially exploiting Allee effects for therapy.

Accurate pitch perception of harmonic complex tones is widely believed to rely on temporal fine structure information conveyed by the precise phase-locked responses of auditory-nerve fibers. However, accurate pitch perception remains possible even when spectrally resolved harmonics are presented at frequencies beyond the putative limits of neural phase locking, and it is unclear whether residual temporal information, or a coarser rate-place code, underlies this ability. We addressed this question by measuring human pitch discrimination at low and high frequencies for harmonic complex tones, presented either in isolation or in the presence of concurrent complex-tone maskers. We found that concurrent complex-tone maskers impaired performance at both low and high frequencies, although the impairment introduced by adding maskers at high frequencies relative to low frequencies differed between the tested masker types. We then combined simulated auditory-nerve responses to our stimuli with ideal-observer analysis to quantify the extent to which performance was limited by peripheral factors. We found that the worsening of both frequency discrimination and F0 discrimination at high frequencies could be well accounted for (in relative terms) by optimal decoding of all available information at the level of the auditory nerve. A Python package is provided to reproduce these results, and to simulate responses to acoustic stimuli from the three previously published models of the human auditory nerve used in our analyses.

Protein-protein interactions (PPIs) are ubiquitous biomolecular processes that are central to virtually all aspects of cellular function. Identifying small molecules that modulate specific disease-related PPIs is a strategy with enormous promise for drug discovery. The design of drugs to disrupt PPIs is challenging, however, because many potential drug-binding sites at PPI interfaces are "cryptic": When unoccupied by a ligand, cryptic sites are often flat and featureless, and thus not readily recognizable in crystal structures, with the geometric and chemical characteristics of typical small-molecule binding sites only emerging upon ligand binding. The rational design of small molecules to inhibit specific PPIs would benefit from a better understanding of how such molecules bind at PPI interfaces. To this end, we have conducted unbiased, all-atom MD simulations of the binding of four small-molecule inhibitors (SP4206 and three SP4206 analogs) to interleukin 2 (IL2)-which performs its function by forming a PPI with its receptor-without incorporating any prior structural information about the ligands' binding. In multiple binding events, a small molecule settled into a stable binding pose at the PPI interface of IL2, resulting in a protein-small-molecule binding site and pose virtually identical to that observed in an existing crystal structure of the IL2-SP4206 complex. Binding of the small molecule stabilized the IL2 binding groove, which when the small molecule was not bound emerged only transiently and incompletely. Moreover, free energy perturbation (FEP) calculations successfully distinguished between the native and non-native IL2-small-molecule binding poses found in the simulations, suggesting that binding simulations in combination with FEP may provide an effective tool for identifying cryptic binding sites and determining the binding poses of small molecules designed to disrupt PPI interfaces by binding to such sites.

The basic reproduction number (R0) of an infection determines the impact of its control. For many endemic infections, R0 is often estimated from appropriate country-specific seroprevalence data. Studies sometimes pool estimates from the same region for settings lacking seroprevalence data, but the reliability of this approach is unclear. Plausibly, indicator-based approaches could predict R0 for such settings. We calculated R0 for rubella for 98 settings and correlated its value against 66 demographic, economic, education, housing and health-related indicators. We also trained a random forest regression algorithm using these indicators as the input and R0 as the output. We used the mean-square error to compare the performances of the random forest, simple linear regression and a regional averaging method in predicting R0 using 4-fold cross validation. R0 was <5, 5-10 and >10 for 81, 14 and 3 settings respectively, with no apparent regional differences and in the limited available data, it was usually lower for rural than urban areas. R0 was most correlated with educational attainment, and household indicators for the Pearson and Spearman correlation coefficients respectively and with poverty-related indicators followed by the crude death rate considering the Maximum Information Coefficient, although the correlation for each was relatively weak (Pearson correlation coefficient: 0.4, 95%CI: (0.24,0.48) for educational attainment). A random forest did not perform better in predicting R0 than simple linear regression, depending on the subsets of training indicators and studies, and neither out-performed a regional averaging approach. R0 for rubella is typically low and using indicators to estimate its value is not straightforward. A regional averaging approach may provide as reliable an estimate of R0 for settings lacking seroprevalence data as one based on indicators. The findings may be relevant for other infections and studies estimating the disease burden and the impact of interventions for settings lacking seroprevalence data.

Converging evidence suggests the brain encodes time in dynamic patterns of neural activity, including neural sequences, ramping activity, and complex dynamics. Most temporal tasks, however, require more than just encoding time, and can have distinct computational requirements including the need to exhibit temporal scaling, generalize to novel contexts, or robustness to noise. It is not known how neural circuits can encode time and satisfy distinct computational requirements, nor is it known whether similar patterns of neural activity at the population level can exhibit dramatically different computational or generalization properties. To begin to answer these questions, we trained RNNs on two timing tasks based on behavioral studies. The tasks had different input structures but required producing identically timed output patterns. Using a novel framework we quantified whether RNNs encoded two intervals using either of three different timing strategies: scaling, absolute, or stimulus-specific dynamics. We found that similar neural dynamic patterns at the level of single intervals, could exhibit fundamentally different properties, including, generalization, the connectivity structure of the trained networks, and the contribution of excitatory and inhibitory neurons. Critically, depending on the task structure RNNs were better suited for generalization or robustness to noise. Further analysis revealed different connection patterns underlying the different regimes. Our results predict that apparently similar neural dynamic patterns at the population level (e.g., neural sequences) can exhibit fundamentally different computational properties in regards to their ability to generalize to novel stimuli and their robustness to noise-and that these differences are associated with differences in network connectivity and distinct contributions of excitatory and inhibitory neurons. We also predict that the task structure used in different experimental studies accounts for some of the experimentally observed variability in how networks encode time.

Metastatic cancer accounts for over 90% of all cancer deaths, and evaluations of metastasis potential are vital for minimizing the metastasis-associated mortality and achieving optimal clinical decision-making. Computational assessment of metastasis potential based on large-scale transcriptomic cancer data is challenging because metastasis events are not always clinically detectable. The under-diagnosis of metastasis events results in biased classification labels, and classification tools using biased labels may lead to inaccurate estimations of metastasis potential. This issue is further complicated by the unknown metastasis prevalence at the population level, the small number of confirmed metastasis cases, and the high dimensionality of the candidate molecular features. Our proposed algorithm, called Positive and unlabeled Learning from Unbalanced cases and Sparse structures (PLUS), is the first to use a positive and unlabeled learning framework to account for the under-detection of metastasis events in building a classifier. PLUS is specifically tailored for studying metastasis that deals with the unbalanced instance allocation as well as unknown metastasis prevalence, which are not considered by other methods. PLUS achieves superior performance on synthetic datasets compared with other state-of-the-art methods. Application of PLUS to The Cancer Genome Atlas Pan-Cancer gene expression data generated metastasis potential predictions that show good agreement with the clinical follow-up data, in addition to predictive genes that have been validated by independent single-cell RNA-sequencing datasets.

High throughput biology enables the measurements of relative concentrations of thousands of biomolecules from e.g. tissue samples. The process leaves the investigator with the problem of how to best interpret the potentially large number of differences between samples. Many activities in a cell depend on ordered reactions involving multiple biomolecules, often referred to as pathways. It hence makes sense to study differences between samples in terms of altered pathway activity, using so-called pathway analysis. Traditional pathway analysis gives significance to differences in the pathway components' concentrations between sample groups, however, less frequently used methods for estimating individual samples' pathway activities have been suggested. Here we demonstrate that such a method can be used for pathway-based survival analysis. Specifically, we investigate the pathway activities' association with patients' survival time based on the transcription profiles of the METABRIC dataset. Our implementation shows that pathway activities are better prognostic markers for survival time in METABRIC than the individual transcripts. We also demonstrate that we can regress out the effect of individual pathways on other pathways, which allows us to estimate the other pathways' residual pathway activity on survival. Furthermore, we illustrate how one can visualize the often interdependent measures over hierarchical pathway databases using sunburst plots.

Testing is recommended for all close contacts of confirmed COVID-19 patients. However, existing pooled testing methods are oblivious to the circumstances of contagion provided by contact tracing. Here, we build upon a well-known semi-adaptive pooled testing method, Dorfman's method with imperfect tests, and derive a simple pooled testing method based on dynamic programming that is specifically designed to use information provided by contact tracing. Experiments using a variety of reproduction numbers and dispersion levels, including those estimated in the context of the COVID-19 pandemic, show that the pools found using our method result in a significantly lower number of tests than those found using Dorfman's method. Our method provides the greatest competitive advantage when the number of contacts of an infected individual is small, or the distribution of secondary infections is highly overdispersed. Moreover, it maintains this competitive advantage under imperfect contact tracing and significant levels of dilution.

Temporal environmental variations affect diversity in communities of competing populations. In particular, the covariance between competition and environment is known to facilitate invasions of rare species via the storage effect. Here we present a quantitative study of the effects of temporal variations in two-species and in diverse communities. Four scenarios are compared: environmental variations may be either periodic (seasonal) or stochastic, and the dynamics may support the storage effect (global competition) or not (local competition). In two-species communities, coexistence is quantified via the mean time to absorption, and we show that stochastic variations yield shorter persistence time because they allow for rare sequences of bad years. In diverse communities, where the steady-state reflects a colonization-extinction equilibrium, the actual number of temporal niches is shown to play a crucial role. When this number is large, the same trends hold: storage effect and periodic variations increase both species richness and the evenness of the community. Surprisingly, when the number of temporal niches is small global competition acts to decrease species richness and evenness, as it focuses the competition to specific periods, thus increasing the effective fitness differences.

Tessellations emerge in many natural systems, and the constituent domains often contain regular patterns, raising the intriguing possibility that pattern formation within adjacent domains might be correlated by the geometry, without the direct exchange of information between parts comprising either domain. We confirm this paradoxical effect, by simulating pattern formation via reaction-diffusion in domains whose boundary shapes tessellate, and showing that correlations between adjacent patterns are strong compared to controls that self-organize in domains with equivalent sizes but unrelated shapes. The effect holds in systems with linear and non-linear diffusive terms, and for boundary shapes derived from regular and irregular tessellations. Based on the prediction that correlations between adjacent patterns should be bimodally distributed, we develop methods for testing whether a given set of domain boundaries constrained pattern formation within those domains. We then confirm such a prediction by analysing the development of 'subbarrel' patterns, which are thought to emerge via reaction-diffusion, and whose enclosing borders form a Voronoi tessellation on the surface of the rodent somatosensory cortex. In more general terms, this result demonstrates how causal links can be established between the dynamical processes through which biological patterns emerge and the constraints that shape them.

There is no single way to represent a task. Indeed, despite experiencing the same task events and contingencies, different subjects may form distinct task representations. As experimenters, we often assume that subjects represent the task as we envision it. However, such a representation cannot be taken for granted, especially in animal experiments where we cannot deliver explicit instruction regarding the structure of the task. Here, we tested how rats represent an odor-guided choice task in which two odor cues indicated which of two responses would lead to reward, whereas a third odor indicated free choice among the two responses. A parsimonious task representation would allow animals to learn from the forced trials what is the better option to choose in the free-choice trials. However, animals may not necessarily generalize across odors in this way. We fit reinforcement-learning models that use different task representations to trial-by-trial choice behavior of individual rats performing this task, and quantified the degree to which each animal used the more parsimonious representation, generalizing across trial types. Model comparison revealed that most rats did not acquire this representation despite extensive experience. Our results demonstrate the importance of formally testing possible task representations that can afford the observed behavior, rather than assuming that animals' task representations abide by the generative task structure that governs the experimental design.

The SARS-CoV-2 pandemic is a major concern all over the world and, as vaccines became available at the end of 2020, optimal vaccination strategies were subjected to intense investigation. Considering their critical role in reducing disease burden, the increasing demand outpacing production, and that most currently approved vaccines follow a two-dose regimen, the cost-effectiveness of delaying the second dose to increment the coverage of the population receiving the first dose is often debated. Finding the best solution is complex due to the trade-off between vaccinating more people with lower level of protection and guaranteeing higher protection to a fewer number of individuals. Here we present a novel extended age-structured SEIR mathematical model that includes a two-dose vaccination schedule with a between-doses delay modelled through delay differential equations and linear optimization of vaccination rates. By maintaining the minimum stock of vaccines under a given production rate, we evaluate the dose interval that minimizes the number of deaths. We found that the best strategy depends on an interplay between the vaccine production rate and the relative efficacy of the first dose. In the scenario of low first-dose efficacy, it is always better to vaccinate the second dose as soon as possible, while for high first-dose efficacy, the best strategy of time window depends on the production rate and also on second-dose efficacy provided by each type of vaccine. We also found that the rate of spread of the infection does not affect significantly the thresholds of the best window, but is an important factor in the absolute number of total deaths. These conclusions point to the need to carefully take into account both vaccine characteristics and roll-out speed to optimize the outcome of vaccination strategies.

Anthropogenic environments such as those created by intensive farming of livestock, have been proposed to provide ideal selection pressure for the emergence of antimicrobial-resistant Escherichia coli bacteria and antimicrobial resistance genes (ARGs) and spread to humans. Here, we performed a longitudinal study in a large-scale commercial poultry farm in China, collecting E. coli isolates from both farm and slaughterhouse; targeting animals, carcasses, workers and their households and environment. By using whole-genome phylogenetic analysis and network analysis based on single nucleotide polymorphisms (SNPs), we found highly interrelated non-pathogenic and pathogenic E. coli strains with phylogenetic intermixing, and a high prevalence of shared multidrug resistance profiles amongst livestock, human and environment. Through an original data processing pipeline which combines omics, machine learning, gene sharing network and mobile genetic elements analysis, we investigated the resistance to 26 different antimicrobials and identified 361 genes associated to antimicrobial resistance (AMR) phenotypes; 58 of these were known AMR-associated genes and 35 were associated to multidrug resistance. We uncovered an extensive network of genes, correlated to AMR phenotypes, shared among livestock, humans, farm and slaughterhouse environments. We also found several human, livestock and environmental isolates sharing closely related mobile genetic elements carrying ARGs across host species and environments. In a scenario where no consensus exists on how antibiotic use in the livestock may affect antibiotic resistance in the human population, our findings provide novel insights into the broader epidemiology of antimicrobial resistance in livestock farming. Moreover, our original data analysis method has the potential to uncover AMR transmission pathways when applied to the study of other pathogens active in other anthropogenic environments characterised by complex interconnections between host species.

Protein structure fundamentally underpins the function and processes of numerous biological systems. Fold recognition algorithms offer a sensitive and robust tool to detect structural, and thereby functional, similarities between distantly related homologs. In the era of accurate structure prediction owing to advances in machine learning techniques and a wealth of experimentally determined structures, previously curated sequence databases have become a rich source of biological information. Here, we use bioinformatic fold recognition algorithms to scan the entire AlphaFold structure database to identify novel protein family members, infer function and group predicted protein structures. As an example of the utility of this approach, we identify novel, previously unknown members of various pore-forming protein families, including MACPFs, GSDMs and aerolysin-like proteins.

Being permanently confronted with an uncertain world, brains have faced evolutionary pressure to represent this uncertainty in order to respond appropriately. Often, this requires visiting multiple interpretations of the available information or multiple solutions to an encountered problem. This gives rise to the so-called mixing problem: since all of these "valid" states represent powerful attractors, but between themselves can be very dissimilar, switching between such states can be difficult. We propose that cortical oscillations can be effectively used to overcome this challenge. By acting as an effective temperature, background spiking activity modulates exploration. Rhythmic changes induced by cortical oscillations can then be interpreted as a form of simulated tempering. We provide a rigorous mathematical discussion of this link and study some of its phenomenological implications in computer simulations. This identifies a new computational role of cortical oscillations and connects them to various phenomena in the brain, such as sampling-based probabilistic inference, memory replay, multisensory cue combination, and place cell flickering.

R is an increasingly preferred software environment for data analytics and statistical computing among scientists and practitioners. Packages markedly extend R's utility and ameliorate inefficient solutions to data science problems. We outline 10 simple rules for finding relevant packages and determining which package is best for your desired use. We begin in Rule 1 with tips on how to consider your purpose, which will guide your search to follow, where, in Rule 2, you'll learn best practices for finding and collecting options. Rules 3 and 4 will help you navigate packages' profiles and explore the extent of their online resources, so that you can be confident in the quality of the package you choose and assured that you'll be able to access support. In Rules 5 and 6, you'll become familiar with how the R Community evaluates packages and learn how to assess the popularity and utility of packages for yourself. Rules 7 and 8 will teach you how to investigate and track package development processes, so you can further evaluate their merit. We end in Rules 9 and 10 with more hands-on approaches, which involve digging into package code.

Stable isotope-assisted metabolic flux analysis (MFA) is a powerful method to estimate carbon flow and partitioning in metabolic networks. At its core, MFA is a parameter estimation problem wherein the fluxes and metabolite pool sizes are model parameters that are estimated, via optimization, to account for measurements of steady-state or isotopically-nonstationary isotope labeling patterns. As MFA problems advance in scale, they require efficient computational methods for fast and robust convergence. The structure of the MFA problem enables it to be cast as an equality-constrained nonlinear program (NLP), where the equality constraints are constructed from the MFA model equations, and the objective function is defined as the sum of squared residuals (SSR) between the model predictions and a set of labeling measurements. This NLP can be solved by using an algebraic modeling language (AML) that offers state-of-the-art optimization solvers for robust parameter estimation and superior scalability to large networks. When implemented in this manner, the optimization is performed with no distinction between state variables and model parameters. During each iteration of such an optimization, the system state is updated instead of being calculated explicitly from scratch, and this occurs concurrently with improvement in the model parameter estimates. This optimization approach starkly contrasts with traditional "shooting" methods where the state variables and model parameters are kept distinct and the system state is computed afresh during each iteration of a stepwise optimization. Our NLP formulation uses the MFA modeling framework of Wiechert et al. [1], which is amenable to incorporation of the model equations into an NLP. The NLP constraints consist of balances on either elementary metabolite units (EMUs) or cumomers. In this formulation, both the steady-state and isotopically-nonstationary MFA (inst-MFA) problems may be solved as an NLP. For the inst-MFA case, the ordinary differential equation (ODE) system describing the labeling dynamics is transcribed into a system of algebraic constraints for the NLP using collocation. This large-scale NLP may be solved efficiently using an NLP solver implemented on an AML. In our implementation, we used the reduced gradient solver CONOPT, implemented in the General Algebraic Modeling System (GAMS). The NLP framework is particularly advantageous for inst-MFA, scaling well to large networks with many free parameters, and having more robust convergence properties compared to the shooting methods that compute the system state and sensitivities at each iteration. Additionally, this NLP approach supports the use of tandem-MS data for both steady-state and inst-MFA when the cumomer framework is used. We assembled a software, eiFlux, written in Python and GAMS that uses the NLP approach and supports both steady-state and inst-MFA. We demonstrate the effectiveness of the NLP formulation on several examples, including a genome-scale inst-MFA model, to highlight the scalability and robustness of this approach. In addition to typical inst-MFA applications, we expect that this framework and our associated software, eiFlux, will be particularly useful for applying inst-MFA to complex MFA models, such as those developed for eukaryotes (e.g. algae) and co-cultures with multiple cell types.

Protein structure alignment algorithms are often time-consuming, resulting in challenges for large-scale protein structure similarity-based retrieval. There is an urgent need for more efficient structure comparison approaches as the number of protein structures increases rapidly. In this paper, we propose an effective graph-based protein structure representation learning method, GraSR, for fast and accurate structure comparison. In GraSR, a graph is constructed based on the intra-residue distance derived from the tertiary structure. Then, deep graph neural networks (GNNs) with a short-cut connection learn graph representations of the tertiary structures under a contrastive learning framework. To further improve GraSR, a novel dynamic training data partition strategy and length-scaling cosine distance are introduced. We objectively evaluate our method GraSR on SCOPe v2.07 and a new released independent test set from PDB database with a designed comprehensive performance metric. Compared with other state-of-the-art methods, GraSR achieves about 7%-10% improvement on two benchmark datasets. GraSR is also much faster than alignment-based methods. We dig into the model and observe that the superiority of GraSR is mainly brought by the learned discriminative residue-level and global descriptors. The web-server and source code of GraSR are freely available at www.csbio.sjtu.edu.cn/bioinf/GraSR/ for academic use.

RNA sequencing has been widely used as an essential tool to probe gene expression. While standard practices have been established to analyze RNA-seq data, it is still challenging to interpret and remove artifactual signals. Several biological and technical factors such as sex, age, batches, and sequencing technology have been found to bias these estimates. Probabilistic estimation of expression residuals (PEER), which infers broad variance components in gene expression measurements, has been used to account for some systematic effects, but it has remained challenging to interpret these PEER factors. Here we show that transcriptome diversity-a simple metric based on Shannon entropy-explains a large portion of variability in gene expression and is the strongest known factor encoded in PEER factors. We then show that transcriptome diversity has significant associations with multiple technical and biological variables across diverse organisms and datasets. In sum, transcriptome diversity provides a simple explanation for a major source of variation in both gene expression estimates and PEER covariates.

The functional near-infrared spectroscopy (fNIRS) can detect hemodynamic responses in the brain and the data consist of bivariate time series of oxygenated hemoglobin (oxy-Hb) and deoxygenated hemoglobin (deoxy-Hb) on each channel. In this study, we investigate oscillatory changes in infant fNIRS signals by using the oscillator decompisition method (OSC-DECOMP), which is a statistical method for extracting oscillators from time series data based on Gaussian linear state space models. OSC-DECOMP provides a natural decomposition of fNIRS data into oscillation components in a data-driven manner and does not require the arbitrary selection of band-pass filters. We analyzed 18-ch fNIRS data (3 minutes) acquired from 21 sleeping 3-month-old infants. Five to seven oscillators were extracted on most channels, and their frequency distribution had three peaks in the vicinity of 0.01-0.1 Hz, 1.6-2.4 Hz and 3.6-4.4 Hz. The first peak was considered to reflect hemodynamic changes in response to the brain activity, and the phase difference between oxy-Hb and deoxy-Hb for the associated oscillators was at approximately 230 degrees. The second peak was attributed to cardiac pulse waves and mirroring noise. Although these oscillators have close frequencies, OSC-DECOMP can separate them through estimating their different projection patterns on oxy-Hb and deoxy-Hb. The third peak was regarded as the harmonic of the second peak. By comparing the Akaike Information Criterion (AIC) of two state space models, we determined that the time series of oxy-Hb and deoxy-Hb on each channel originate from common oscillatory activity. We also utilized the result of OSC-DECOMP to investigate the frequency-specific functional connectivity. Whereas the brain oscillator exhibited functional connectivity, the pulse waves and mirroring noise oscillators showed spatially homogeneous and independent changes. OSC-DECOMP is a promising tool for data-driven extraction of oscillation components from biological time series data.

Hypophosphatasia (HPP) is a rare inherited disorder characterized by defective bone mineralization and is highly variable in its clinical phenotype. The disease occurs due to various loss-of-function mutations in ALPL, the gene encoding tissue-nonspecific alkaline phosphatase (TNSALP). In this work, a data-driven and biophysics-based approach is proposed for the large-scale analysis of ALPL mutations-from nonpathogenic to severe HPPs. By using a pipeline of synergistic approaches including sequence-structure analysis, network modeling, elastic network models and atomistic simulations, we characterized allosteric signatures and effects of the ALPL mutations on protein dynamics and function. Statistical analysis of molecular features computed for the ALPL mutations showed a significant difference between the control, mild and severe HPP phenotypes. Molecular dynamics simulations coupled with protein structure network analysis were employed to analyze the effect of single-residue variation on conformational dynamics of TNSALP dimers, and the developed machine learning model suggested that the topological network parameters could serve as a robust indicator of severe mutations. The results indicated that the severity of disease-associated mutations is often linked with mutation-induced modulation of allosteric communications in the protein. This study suggested that ALPL mutations associated with mild and more severe HPPs can exert markedly distinct effects on the protein stability and long-range network communications. By linking the disease phenotypes with dynamic and allosteric molecular signatures, the proposed integrative computational approach enabled to characterize and quantify the allosteric effects of ALPL mutations and role of allostery in the pathogenesis of HPPs.

Composition and functions of microbial communities affect important traits in diverse hosts, from crops to humans. Yet, mechanistic understanding of how metabolism of individual microbes is affected by the community composition and metabolite leakage is lacking. Here, we first show that the consensus of automatically generated metabolic reconstructions improves the quality of the draft reconstructions, measured by comparison to reference models. We then devise an approach for gap filling, termed COMMIT, that considers metabolites for secretion based on their permeability and the composition of the community. By applying COMMIT with two soil communities from the Arabidopsis thaliana culture collection, we could significantly reduce the gap-filling solution in comparison to filling gaps in individual reconstructions without affecting the genomic support. Inspection of the metabolic interactions in the soil communities allows us to identify microbes with community roles of helpers and beneficiaries. Therefore, COMMIT offers a versatile fully automated solution for large-scale modelling of microbial communities for diverse biotechnological applications.

Implicit motor recalibration allows us to flexibly move in novel and changing environments. Conventionally, implicit recalibration is thought to be driven by errors in predicting the sensory outcome of movement (i.e., sensory prediction errors). However, recent studies have shown that implicit recalibration is also influenced by errors in achieving the movement goal (i.e., task errors). Exactly how sensory prediction errors and task errors interact to drive implicit recalibration and, in particular, whether task errors alone might be sufficient to drive implicit recalibration remain unknown. To test this, we induced task errors in the absence of sensory prediction errors by displacing the target mid-movement. We found that task errors alone failed to induce implicit recalibration. In additional experiments, we simultaneously varied the size of sensory prediction errors and task errors. We found that implicit recalibration driven by sensory prediction errors could be continuously modulated by task errors, revealing an unappreciated dependency between these two sources of error. Moreover, implicit recalibration was attenuated when the target was simply flickered in its original location, even though this manipulation did not affect task error - an effect likely attributed to attention being directed away from the feedback cursor. Taken as a whole, the results were accounted for by a computational model in which sensory prediction errors and task errors, modulated by attention, interact to determine the extent of implicit recalibration.

To date, genome-wide association studies have identified thousands of statistically-significant associations between genetic variants, and phenotypes related to a myriad of traits and diseases. A key goal for human-genetics research is to translate these associations into functional mechanisms. Popular gene-set analysis tools, like MAGMA, map variants to genes they might affect, and then integrate genome-wide association study data (that is, variant-level associations for a phenotype) to score genes for association with a phenotype. Gene scores are subsequently used in competitive gene-set analyses to identify biological processes that are enriched for phenotype association. By default, variants are mapped to genes in their proximity. However, many variants that affect phenotypes are thought to act at regulatory elements, which can be hundreds of kilobases away from their target genes. Thus, we explored the idea of augmenting a proximity-based mapping scheme with publicly-available datasets of regulatory interactions. We used MAGMA to analyze genome-wide association study data for ten different phenotypes, and evaluated the effects of augmentation by comparing numbers, and identities, of genes and gene sets detected as statistically significant between mappings. We detected several pitfalls and confounders of such "augmented analyses", and introduced ways to control for them. Using these controls, we demonstrated that augmentation with datasets of regulatory interactions only occasionally strengthened the enrichment for phenotype association amongst (biologically-relevant) gene sets for different phenotypes. Still, in such cases, genes and regulatory elements responsible for the improvement could be pinpointed. For instance, using brain regulatory-interactions for augmentation, we were able to implicate two acetylcholine receptor subunits involved in post-synaptic chemical transmission, namely CHRNB2 and CHRNE, in schizophrenia. Collectively, our study presents a critical approach for integrating regulatory interactions into gene-set analyses for genome-wide association study data, by introducing various controls to distinguish genuine results from spurious discoveries.

Update of bioRxiv. 2021 Apr 27;:

Although chemotherapy is a standard treatment for cancer, it comes with significant side effects. In particular, certain agents can induce severe muscle loss, known as cachexia, worsening patient quality of life and treatment outcomes. 5-fluorouracil, an anti-cancer agent used to treat several cancers, has been shown to cause muscle loss. Experimental data indicates a non-linear dose-dependence for muscle loss in mice treated with daily or week-day schedules. We present a mathematical model of chemotherapy-induced muscle wasting that captures this non-linear dose-dependence. Area-under-the-curve metrics are proposed to quantify the treatment's effects on lean mass and tumour control. Model simulations are used to explore alternate dosing schedules, aging effects, and morphine use in chemotherapy treatment with the aim of better protecting lean mass while actively targeting the tumour, ultimately leading to improved personalization of treatment planning and improved patient quality of life.

Clathrin-coated structures must assemble on cell membranes to internalize receptors, with the clathrin protein only linked to the membrane via adaptor proteins. These structures can grow surprisingly large, containing over 20 clathrin, yet they often fail to form productive vesicles, instead aborting and disassembling. We show that clathrin structures of this size can both form and disassemble spontaneously when adaptor protein availability is low, despite high abundance of clathrin. Here, we combine recent in vitro kinetic measurements with microscopic reaction-diffusion simulations and theory to differentiate mechanisms of stable vs unstable clathrin assembly on membranes. While in vitro conditions drive assembly of robust, stable lattices, we show that concentrations, geometry, and dimensional reduction in physiologic-like conditions do not support nucleation if only the key adaptor AP-2 is included, due to its insufficient abundance. Nucleation requires a stoichiometry of adaptor to clathrin that exceeds 1:1, meaning additional adaptor types are necessary to form lattices successfully and efficiently. We show that the critical nucleus contains ~25 clathrin, remarkably similar to sizes of the transient and abortive structures observed in vivo. Lastly, we quantify the cost of bending the membrane under our curved clathrin lattices using a continuum membrane model. We find that the cost of bending the membrane could be largely offset by the energetic benefit of forming curved rather than flat structures, with numbers comparable to experiments. Our model predicts how adaptor density can tune clathrin-coated structures from the transient to the stable, showing that active energy consumption is therefore not required for lattice disassembly or remodeling during growth, which is a critical advance towards predicting productive vesicle formation.

Focal sources (FS) are believed to be important triggers and a perpetuation mechanism for paroxysmal atrial fibrillation (AF). Detecting FS and determining AF sustainability in atrial tissue can help guide ablation targeting. We hypothesized that sustained rotors during FS-driven episodes indicate an arrhythmogenic substrate for sustained AF, and that non-invasive electrical recordings, like electrocardiograms (ECGs) or body surface potential maps (BSPMs), could be used to detect FS and AF sustainability. Computer simulations were performed on five bi-atrial geometries. FS were induced by pacing at cycle lengths of 120-270 ms from 32 atrial sites and four pulmonary veins. Self-sustained reentrant activities were also initiated around the same 32 atrial sites with inexcitable cores of radii of 0, 0.5 and 1 cm. FS fired for two seconds and then AF inducibility was tested by whether activation was sustained for another second. ECGs and BSPMs were simulated. Equivalent atrial sources were extracted using second-order blind source separation, and their cycle length, periodicity and contribution, were used as features for random forest classifiers. Longer rotor duration during FS-driven episodes indicates higher AF inducibility (area under ROC curve = 0.83). Our method had accuracy of 90.6±1.0% and 90.6±0.6% in detecting FS presence, and 93.1±0.6% and 94.2±1.2% in identifying AF sustainability, and 80.0±6.6% and 61.0±5.2% in determining the atrium of the focal site, from BSPMs and ECGs of five atria. The detection of FS presence and AF sustainability were insensitive to vest placement (±9.6%). On pre-operative BSPMs of 52 paroxysmal AF patients, patients classified with initiator-type FS on a single atrium resulted in improved two-to-three-year AF-free likelihoods (p-value < 0.01, logrank tests). Detection of FS and arrhythmogenic substrate can be performed from ECGs and BSPMs, enabling non-invasive mapping towards mechanism-targeted AF treatment, and malignant ectopic beat detection with likely AF progression.

New technologies, such as multiplex immunofluorescence microscopy (mIF), are being developed and used for the assessment and visualization of the tumor immune microenvironment (TIME). These assays produce not only an estimate of the abundance of immune cells in the TIME, but also their spatial locations. However, there are currently few approaches to analyze the spatial context of the TIME. Therefore, we have developed a framework for the spatial analysis of the TIME using Ripley's K, coupled with a permutation-based framework to estimate and measure the departure from complete spatial randomness (CSR) as a measure of the interactions between immune cells. This approach was then applied to epithelial ovarian cancer (EOC) using mIF collected on intra-tumoral regions of interest (ROIs) and tissue microarrays (TMAs) from 160 high-grade serous ovarian carcinoma patients in the African American Cancer Epidemiology Study (AACES) (94 subjects on TMAs resulting in 263 tissue cores; 93 subjects with 260 ROIs; 27 subjects with both TMA and ROI data). Cox proportional hazard models were constructed to determine the association of abundance and spatial clustering of tumor-infiltrating lymphocytes (CD3+), cytotoxic T-cells (CD8+CD3+), and regulatory T-cells (CD3+FoxP3+) with overall survival. Analysis was done on TMA and ROIs, treating the TMA data as validation of the findings from the ROIs. We found that EOC patients with high abundance and low spatial clustering of tumor-infiltrating lymphocytes and T-cell subsets in their tumors had the best overall survival. Additionally, patients with EOC tumors displaying high co-occurrence of cytotoxic T-cells and regulatory T-cells had the best overall survival. Grouping women with ovarian cancer based on both cell abundance and spatial contexture showed better discrimination for survival than grouping ovarian cancer cases only by cell abundance. These findings underscore the prognostic importance of evaluating not only immune cell abundance but also the spatial contexture of the immune cells in the TIME. In conclusion, the application of this spatial analysis framework to the study of the TIME could lead to the identification of immune content and spatial architecture that could aid in the determination of patients that are likely to respond to immunotherapies.

Glucose is central to many biological processes, serving as an energy source and a building block for biosynthesis. After glucose enters the cell, hexokinases convert it to glucose-6-phosphate (Glc-6P) for use in anaerobic fermentation, aerobic oxidative phosphorylation, and the pentose-phosphate pathway. We here describe a genetic screen in Saccharomyces cerevisiae that generated a novel spontaneous mutation in hexokinase-2, hxk2G238V, that confers resistance to the toxic glucose analog 2-deoxyglucose (2DG). Wild-type hexokinases convert 2DG to 2-deoxyglucose-6-phosphate (2DG-6P), but 2DG-6P cannot support downstream glycolysis, resulting in a cellular starvation-like response. Curiously, though the hxk2G238V mutation encodes a loss-of-function allele, the affected amino acid does not interact directly with bound glucose, 2DG, or ATP. Molecular dynamics simulations suggest that Hxk2G238V impedes sugar binding by altering the protein dynamics of the glucose-binding cleft, as well as the large-scale domain-closure motions required for catalysis. These findings shed new light on Hxk2 dynamics and highlight how allosteric changes can influence catalysis, providing new structural insights into this critical regulator of carbohydrate metabolism. Given that hexokinases are upregulated in some cancers and that 2DG and its derivatives have been studied in anti-cancer trials, the present work also provides insights that may apply to cancer biology and drug resistance.

In moving animal groups, social interactions play a key role in the ability of individuals to achieve coordinated motion. However, a large number of environmental and cognitive factors are able to modulate the expression of these interactions and the characteristics of the collective movements that result from these interactions. Here, we use a data-driven fish school model to quantitatively investigate the impact of perceptual and cognitive factors on coordination and collective swimming patterns. The model describes the interactions involved in the coordination of burst-and-coast swimming in groups of Hemigrammus rhodostomus. We perform a comprehensive investigation of the respective impacts of two interactions strategies between fish based on the selection of the most or the two most influential neighbors, of the range and intensity of social interactions, of the intensity of individual random behavioral fluctuations, and of the group size, on the ability of groups of fish to coordinate their movements. We find that fish are able to coordinate their movements when they interact with their most or two most influential neighbors, provided that a minimal level of attraction between fish exist to maintain group cohesion. A minimal level of alignment is also required to allow the formation of schooling and milling. However, increasing the strength of social interactions does not necessarily enhance group cohesion and coordination. When attraction and alignment strengths are too high, or when the heading random fluctuations are too large, schooling and milling can no longer be maintained and the school switches to a swarming phase. Increasing the interaction range between fish has a similar impact on collective dynamics as increasing the strengths of attraction and alignment. Finally, we find that coordination and schooling occurs for a wider range of attraction and alignment strength in small group sizes.

SARS-CoV-2 Spike (Spike) binds to human angiotensin-converting enzyme 2 (ACE2) and the strength of this interaction could influence parameters relating to virulence. To explore whether population variants in ACE2 influence Spike binding and hence infection, we selected 10 ACE2 variants based on affinity predictions and prevalence in gnomAD and measured their affinities and kinetics for Spike receptor binding domain through surface plasmon resonance (SPR) at 37°C. We discovered variants that reduce and enhance binding, including three ACE2 variants that strongly inhibited (p.Glu37Lys, ΔΔG = -1.33 ± 0.15 kcal mol-1 and p.Gly352Val, predicted ΔΔG = -1.17 kcal mol-1) or abolished (p.Asp355Asn) binding. We also identified two variants with distinct population distributions that enhanced affinity for Spike. ACE2 p.Ser19Pro (ΔΔG = 0.59 ± 0.08 kcal mol-1) is predominant in the gnomAD African cohort (AF = 0.003) whilst p.Lys26Arg (ΔΔG = 0.26 ± 0.09 kcal mol-1) is predominant in the Ashkenazi Jewish (AF = 0.01) and European non-Finnish (AF = 0.006) cohorts. We compared ACE2 variant affinities to published SARS-CoV-2 pseudotype infectivity data and confirmed that ACE2 variants with reduced affinity for Spike can protect cells from infection. The effect of variants with enhanced Spike affinity remains unclear, but we propose a mechanism whereby these alleles could cause greater viral spreading across tissues and cell types, as is consistent with emerging understanding regarding the interplay between receptor affinity and cell-surface abundance. Finally, we compared mCSM-PPI2 ΔΔG predictions against our SPR data to assess the utility of predictions in this system. We found that predictions of decreased binding were well-correlated with experiment and could be improved by calibration, but disappointingly, predictions of highly enhanced binding were unreliable. Recalibrated predictions for all possible ACE2 missense variants at the Spike interface were calculated and used to estimate the overall burden of ACE2 variants on Covid-19.

The human immune system consists of a highly intelligent network of billions of independent, self-organized cells that interact with each other. Machine learning (ML) is an artificial intelligence (AI) tool that automatically processes huge amounts of image data. Immunotherapies have revolutionized the treatment of blood cancer. Specifically, one such therapy involves engineering immune cells to express chimeric antigen receptors (CAR), which combine tumor antigen specificity with immune cell activation in a single receptor. To improve their efficacy and expand their applicability to solid tumors, scientists optimize different CARs with different modifications. However, predicting and ranking the efficacy of different "off-the-shelf" immune products (e.g., CAR or Bispecific T-cell Engager [BiTE]) and selection of clinical responders are challenging in clinical practice. Meanwhile, identifying the optimal CAR construct for a researcher to further develop a potential clinical application is limited by the current, time-consuming, costly, and labor-intensive conventional tools used to evaluate efficacy. Particularly, more than 30 years of immunological synapse (IS) research data demonstrate that T cell efficacy is not only controlled by the specificity and avidity of the tumor antigen and T cell interaction, but also it depends on a collective process, involving multiple adhesion and regulatory molecules, as well as tumor microenvironment, spatially and temporally organized at the IS formed by cytotoxic T lymphocytes (CTL) and natural killer (NK) cells. The optimal function of cytotoxic lymphocytes (including CTL and NK) depends on IS quality. Recognizing the inadequacy of conventional tools and the importance of IS in immune cell functions, we investigate a new strategy for assessing CAR-T efficacy by quantifying CAR IS quality using the glass-support planar lipid bilayer system combined with ML-based data analysis. Previous studies in our group show that CAR-T IS quality correlates with antitumor activities in vitro and in vivo. However, current manually quantified IS quality data analysis is time-consuming and labor-intensive with low accuracy, reproducibility, and repeatability. In this study, we develop a novel ML-based method to quantify thousands of CAR cell IS images with enhanced accuracy and speed. Specifically, we used artificial neural networks (ANN) to incorporate object detection into segmentation. The proposed ANN model extracts the most useful information to differentiate different IS datasets. The network output is flexible and produces bounding boxes, instance segmentation, contour outlines (borders), intensities of the borders, and segmentations without borders. Based on requirements, one or a combination of this information is used in statistical analysis. The ML-based automated algorithm quantified CAR-T IS data correlates with the clinical responder and non-responder treated with Kappa-CAR-T cells directly from patients. The results suggest that CAR cell IS quality can be used as a potential composite biomarker and correlates with antitumor activities in patients, which is sufficiently discriminative to further test the CAR IS quality as a clinical biomarker to predict response to CAR immunotherapy in cancer. For translational research, the method developed here can also provide guidelines for designing and optimizing numerous CAR constructs for potential clinical development. Trial Registration: ClinicalTrials.gov NCT00881920.

Host-microbiome interactions and the microbial community have broad impact in human health and diseases. Most microbiome based studies are performed at the genome level based on next-generation sequencing techniques, but metaproteomics is emerging as a powerful technique to study microbiome functional activity by characterizing the complex and dynamic composition of microbial proteins. We conducted a large-scale survey of human gut microbiome metaproteomic data to identify generalist species that are ubiquitously expressed across all samples and specialists that are highly expressed in a small subset of samples associated with a certain phenotype. We were able to utilize the metaproteomic mass spectrometry data to reveal the protein landscapes of these species, which enables the characterization of the expression levels of proteins of different functions and underlying regulatory mechanisms, such as operons. Finally, we were able to recover a large number of open reading frames (ORFs) with spectral support, which were missed by de novo protein-coding gene predictors. We showed that a majority of the rescued ORFs overlapped with de novo predicted protein-coding genes, but on opposite strands or in different frames. Together, these demonstrate applications of metaproteomics for the characterization of important gut bacterial species.

Epithelial tissues constitute an exotic type of active matter with non-linear properties reminiscent of amorphous materials. In the context of a proliferating epithelium, modeled by the quasistatic vertex model, we identify novel discrete tissue scale rearrangements, i.e. cellular rearrangement avalanches, which are a form of collective cell movement. During the avalanches, the vast majority of cells retain their neighbors, and the resulting cellular trajectories are radial in the periphery, a vortex in the core. After the onset of these avalanches, the epithelial area grows discontinuously. The avalanches are found to be stochastic, and their strength is correlated with the density of cells in the tissue. Overall, avalanches redistribute accumulated local spatial pressure along the tissue. Furthermore, the distribution of avalanche magnitudes is found to obey a power law, with an exponent consistent with sheer induced avalanches in amorphous materials. To understand the role of avalanches in organ development, we simulate epithelial growth of the Drosophila eye disc during the third instar using a computational model, which includes both chemical and mechanistic signaling. During the third instar, the morphogenetic furrow (MF), a ~10 cell wide wave of apical area constriction propagates through the epithelium. These simulations are used to understand the details of the growth process, the effect of the MF on the growth dynamics on the tissue scale, and to make predictions for experimental observations. The avalanches are found to depend on the strength of the apical constriction of cells in the MF, with a stronger apical constriction leading to less frequent and more pronounced avalanches. The results herein highlight the dependence of simulated tissue growth dynamics on relaxation timescales, and serve as a guide for in vitro experiments.

Circadian rhythms in a wide range of organisms are mediated by molecular mechanisms based on transcription-translation feedback. In this paper, we use bifurcation theory to explore mathematical models of genetic oscillators, based on Kim & Forger's interpretation of the circadian clock in mammals. At the core of their models is a negative feedback loop whereby PER proteins (PER1 and PER2) bind to and inhibit their transcriptional activator, BMAL1. For oscillations to occur, the dissociation constant of the PER:BMAL1 complex, [Formula: see text], must be ≤ 0.04 nM, which is orders of magnitude smaller than a reasonable expectation of 1-10 nM for this protein complex. We relax this constraint by two modifications to Kim & Forger's 'single negative feedback' (SNF) model: first, by introducing a multistep reaction chain for posttranscriptional modifications of Per mRNA and posttranslational phosphorylations of PER, and second, by replacing the first-order rate law for degradation of PER in the nucleus by a Michaelis-Menten rate law. These modifications increase the maximum allowable [Formula: see text] to ~2 nM. In a third modification, we consider an alternative rate law for gene transcription to resolve an unrealistically large rate of Per2 transcription at very low concentrations of BMAL1. Additionally, we studied extensions of the SNF model to include a second negative feedback loop (involving REV-ERB) and a supplementary positive feedback loop (involving ROR). Contrary to Kim & Forger's observations of these extended models, we find that, with our modifications, the supplementary positive feedback loop makes the oscillations more robust than observed in the models with one or two negative feedback loops. However, all three models are similarly robust when accounting for circadian rhythms (~24 h period) with [Formula: see text] ≥ 1 nM. Our results provide testable predictions for future experimental studies.

Understanding and characterising biochemical processes inside single cells requires experimental platforms that allow one to perturb and observe the dynamics of such processes as well as computational methods to build and parameterise models from the collected data. Recent progress with experimental platforms and optogenetics has made it possible to expose each cell in an experiment to an individualised input and automatically record cellular responses over days with fine time resolution. However, methods to infer parameters of stochastic kinetic models from single-cell longitudinal data have generally been developed under the assumption that experimental data is sparse and that responses of cells to at most a few different input perturbations can be observed. Here, we investigate and compare different approaches for calculating parameter likelihoods of single-cell longitudinal data based on approximations of the chemical master equation (CME) with a particular focus on coupling the linear noise approximation (LNA) or moment closure methods to a Kalman filter. We show that, as long as cells are measured sufficiently frequently, coupling the LNA to a Kalman filter allows one to accurately approximate likelihoods and to infer model parameters from data even in cases where the LNA provides poor approximations of the CME. Furthermore, the computational cost of filtering-based iterative likelihood evaluation scales advantageously in the number of measurement times and different input perturbations and is thus ideally suited for data obtained from modern experimental platforms. To demonstrate the practical usefulness of these results, we perform an experiment in which single cells, equipped with an optogenetic gene expression system, are exposed to various different light-input sequences and measured at several hundred time points and use parameter inference based on iterative likelihood evaluation to parameterise a stochastic model of the system.

The most common kind of cancer among women is breast cancer. Understanding the tumor microenvironment and the interactions between individual cells and cytokines assists us in arriving at more effective treatments. Here, we develop a data-driven mathematical model to investigate the dynamics of key cell types and cytokines involved in breast cancer development. We use time-course gene expression profiles of a mouse model to estimate the relative abundance of cells and cytokines. We then employ a least-squares optimization method to evaluate the model's parameters based on the mice data. The resulting dynamics of the cells and cytokines obtained from the optimal set of parameters exhibit a decent agreement between the data and predictions. We perform a sensitivity analysis to identify the crucial parameters of the model and then perform a local bifurcation on them. The results reveal a strong connection between adipocytes, IL6, and the cancer population, suggesting them as potential targets for therapies.

Proteins are typically represented by discrete atomic coordinates providing an accessible framework to describe different conformations. However, in some fields proteins are more accurately represented as near-continuous surfaces, as these are imprinted with geometric (shape) and chemical (electrostatics) features of the underlying protein structure. Protein surfaces are dependent on their chemical composition and, ultimately determine protein function, acting as the interface that engages in interactions with other molecules. In the past, such representations were utilized to compare protein structures on global and local scales and have shed light on functional properties of proteins. Here we describe RosettaSurf, a surface-centric computational design protocol, that focuses on the molecular surface shape and electrostatic properties as means for protein engineering, offering a unique approach for the design of proteins and their functions. The RosettaSurf protocol combines the explicit optimization of molecular surface features with a global scoring function during the sequence design process, diverging from the typical design approaches that rely solely on an energy scoring function. With this computational approach, we attempt to address a fundamental problem in protein design related to the design of functional sites in proteins, even when structurally similar templates are absent in the characterized structural repertoire. Surface-centric design exploits the premise that molecular surfaces are, to a certain extent, independent of the underlying sequence and backbone configuration, meaning that different sequences in different proteins may present similar surfaces. We benchmarked RosettaSurf on various sequence recovery datasets and showcased its design capabilities by generating epitope mimics that were biochemically validated. Overall, our results indicate that the explicit optimization of surface features may lead to new routes for the design of functional proteins.

Atomic force microscopy (AFM) can visualize the dynamics of single biomolecules under near-physiological conditions. However, the scanning tip probes only the molecular surface with limited resolution, missing details required to fully deduce functional mechanisms from imaging alone. To overcome such drawbacks, we developed a computational framework to reconstruct 3D atomistic structures from AFM surface scans, employing simulation AFM and automatized fitting to experimental images. We provide applications to AFM images ranging from single molecular machines, protein filaments, to large-scale assemblies of 2D protein lattices, and demonstrate how the obtained full atomistic information advances the molecular understanding beyond the original topographic AFM image. We show that simulation AFM further allows for quantitative molecular feature assignment within measured AFM topographies. Implementation of the developed methods into the versatile interactive interface of the BioAFMviewer software, freely available at www.bioafmviewer.com, presents the opportunity for the broad Bio-AFM community to employ the enormous amount of existing structural and modeling data to facilitate the interpretation of resolution-limited AFM images.

Identifying copy number variants (CNVs) can provide diagnoses to patients and provide important biological insights into human health and disease. Current exome and targeted sequencing approaches cannot detect clinically and biologically-relevant CNVs outside their target area. We present SavvyCNV, a tool which uses off-target read data from exome and targeted sequencing data to call germline CNVs genome-wide. Up to 70% of sequencing reads from exome and targeted sequencing fall outside the targeted regions. We have developed a new tool, SavvyCNV, to exploit this 'free data' to call CNVs across the genome. We benchmarked SavvyCNV against five state-of-the-art CNV callers using truth sets generated from genome sequencing data and Multiplex Ligation-dependent Probe Amplification assays. SavvyCNV called CNVs with high precision and recall, outperforming the five other tools at calling CNVs genome-wide, using off-target or on-target reads from targeted panel and exome sequencing. We then applied SavvyCNV to clinical samples sequenced using a targeted panel and were able to call previously undetected clinically-relevant CNVs, highlighting the utility of this tool within the diagnostic setting. SavvyCNV outperforms existing tools for calling CNVs from off-target reads. It can call CNVs genome-wide from targeted panel and exome data, increasing the utility and diagnostic yield of these tests. SavvyCNV is freely available at https://github.com/rdemolgen/SavvySuite.

Bayesian methods are widely used in the GWAS meta-analysis. But the considerable consumption in both computing time and memory space poses great challenges for large-scale meta-analyses. In this research, we propose an algorithm named SMetABF to rapidly obtain the optimal ABF in the GWAS meta-analysis, where shotgun stochastic search (SSS) is introduced to improve the Bayesian GWAS meta-analysis framework, MetABF. Simulation studies confirm that SMetABF performs well in both speed and accuracy, compared to exhaustive methods and MCMC. SMetABF is applied to real GWAS datasets to find several essential loci related to Parkinson's disease (PD) and the results support the underlying relationship between PD and other autoimmune disorders. Developed as an R package and a web tool, SMetABF will become a useful tool to integrate different studies and identify more variants associated with complex traits.

Many organism behaviors are innate or instinctual and have been "hard-coded" through evolution. Current approaches to understanding these behaviors model evolution as an optimization problem in which the traits of organisms are assumed to optimize an objective function representing evolutionary fitness. Here, we use a mechanistic birth-death dynamics approach to study the evolution of innate behavioral strategies in a simulated population of organisms. In particular, we performed agent-based stochastic simulations and mean-field analyses of organisms exploring random environments and competing with each other to find locations with plentiful resources. We find that when organism density is low, the mean-field model allows us to derive an effective objective function, predicting how the most competitive phenotypes depend on the exploration-exploitation trade-off between the scarcity of high-resource sites and the increase in birth rate those sites offer organisms. However, increasing organism density alters the most competitive behavioral strategies and precludes the derivation of a well-defined objective function. Moreover, there exists a range of densities for which the coexistence of many phenotypes persists for evolutionarily long times.

Our ability to perceive a stable visual world in the presence of continuous movements of the body, head, and eyes has puzzled researchers in the neuroscience field for a long time. We reformulated this problem in the context of hierarchical convolutional neural networks (CNNs)-whose architectures have been inspired by the hierarchical signal processing of the mammalian visual system-and examined perceptual stability as an optimization process that identifies image-defining features for accurate image classification in the presence of movements. Movement signals, multiplexed with visual inputs along overlapping convolutional layers, aided classification invariance of shifted images by making the classification faster to learn and more robust relative to input noise. Classification invariance was reflected in activity manifolds associated with image categories emerging in late CNN layers and with network units acquiring movement-associated activity modulations as observed experimentally during saccadic eye movements. Our findings provide a computational framework that unifies a multitude of biological observations on perceptual stability under optimality principles for image classification in artificial neural networks.

Infections caused by antibiotic-resistant bacteria have become more prevalent during past decades. Yet, it is unknown whether such infections occur in addition to infections with antibiotic-susceptible bacteria, thereby increasing the incidence of infections, or whether they replace such infections, leaving the total incidence unaffected. Observational longitudinal studies cannot separate both mechanisms. Using plasmid-based beta-lactam resistant E. coli as example we applied mathematical modelling to investigate whether seven biological mechanisms would lead to replacement or addition of infections. We use a mathematical neutral null model of individuals colonized with susceptible and/or resistant E. coli, with two mechanisms implying a fitness cost, i.e., increased clearance and decreased growth of resistant strains, and five mechanisms benefitting resistance, i.e., 1) increased virulence, 2) increased transmission, 3) decreased clearance of resistant strains, 4) increased rate of horizontal plasmid transfer, and 5) increased clearance of susceptible E. coli due to antibiotics. Each mechanism is modelled separately to estimate addition to or replacement of antibiotic-susceptible infections. Fitness costs cause resistant strains to die out if other strain characteristics are maintained equal. Under the assumptions tested, increased virulence is the only mechanism that increases the total number of infections. Other benefits of resistance lead to replacement of susceptible infections without changing the total number of infections. As there is no biological evidence that plasmid-based beta-lactam resistance increases virulence, these findings suggest that the burden of disease is determined by attributable effects of resistance rather than by an increase in the number of infections.

Synthetic Lethality (SL) is currently defined as a type of genetic interaction in which the loss of function of either of two genes individually has limited effect in cell viability but inactivation of both genes simultaneously leads to cell death. Given the profound genomic aberrations acquired by tumor cells, which can be systematically identified with -omics data, SL is a promising concept in cancer research. In particular, SL has received much attention in the area of cancer metabolism, due to the fact that relevant functional alterations concentrate on key metabolic pathways that promote cellular proliferation. With the extensive prior knowledge about human metabolic networks, a number of computational methods have been developed to predict SL in cancer metabolism, including the genetic Minimal Cut Sets (gMCSs) approach. A major challenge in the application of SL approaches to cancer metabolism is to systematically integrate tumor microenvironment, given that genetic interactions and nutritional availability are interconnected to support proliferation. Here, we propose a more general definition of SL for cancer metabolism that combines genetic and environmental interactions, namely loss of gene functions and absence of nutrients in the environment. We extend our gMCSs approach to determine this new family of metabolic synthetic lethal interactions. A computational and experimental proof-of-concept is presented for predicting the lethality of dihydrofolate reductase (DHFR) inhibition in different environments. Finally, our approach is applied to identify extracellular nutrient dependences of tumor cells, elucidating cholesterol and myo-inositol depletion as potential vulnerabilities in different malignancies.

Automatic characterization of fluorescent labeling in intact mammalian tissues remains a challenge due to the lack of quantifying techniques capable of segregating densely packed nuclei and intricate tissue patterns. Here, we describe a powerful deep learning-based approach that couples remarkably precise nuclear segmentation with quantitation of fluorescent labeling intensity within segmented nuclei, and then apply it to the analysis of cell cycle dependent protein concentration in mouse tissues using 2D fluorescent still images. First, several existing deep learning-based methods were evaluated to accurately segment nuclei using different imaging modalities with a small training dataset. Next, we developed a deep learning-based approach to identify and measure fluorescent labels within segmented nuclei, and created an ImageJ plugin to allow for efficient manual correction of nuclear segmentation and label identification. Lastly, using fluorescence intensity as a readout for protein concentration, a three-step global estimation method was applied to the characterization of the cell cycle dependent expression of E2F proteins in the developing mouse intestine.

Cells can measure shallow gradients of external signals to initiate and accomplish a migration or a morphogenetic process. Recently, starting from mathematical models like the local-excitation global-inhibition (LEGI) model and with the support of empirical evidence, it has been proposed that cellular communication improves the measurement of an external gradient. However, the mathematical models that have been used have over-simplified geometries (e.g., they are uni-dimensional) or assumptions about cellular communication, which limit the possibility to analyze the gradient sensing ability of more complex cellular systems. Here, we generalize the existing models to study the effects on gradient sensing of cell number, geometry and of long- versus short-range cellular communication in 2D systems representing epithelial tissues. We find that increasing the cell number can be detrimental for gradient sensing when the communication is weak and limited to nearest neighbour cells, while it is beneficial when there is long-range communication. We also find that, with long-range communication, the gradient sensing ability improves for tissues with more disordered geometries; on the other hand, an ordered structure with mostly hexagonal cells is advantageous with nearest neighbour communication. Our results considerably extend the current models of gradient sensing by epithelial tissues, making a step further toward predicting the mechanism of communication and its putative mediator in many biological processes.

Evaluating and tracking wound size is a fundamental metric for the wound assessment process. Good location and size estimates can enable proper diagnosis and effective treatment. Traditionally, laboratory wound healing studies include a collection of images at uniform time intervals exhibiting the wounded area and the healing process in the test animal, often a mouse. These images are then manually observed to determine key metrics -such as wound size progress- relevant to the study. However, this task is a time-consuming and laborious process. In addition, defining the wound edge could be subjective and can vary from one individual to another even among experts. Furthermore, as our understanding of the healing process grows, so does our need to efficiently and accurately track these key factors for high throughput (e.g., over large-scale and long-term experiments). Thus, in this study, we develop a deep learning-based image analysis pipeline that aims to intake non-uniform wound images and extract relevant information such as the location of interest, wound only image crops, and wound periphery size over-time metrics. In particular, our work focuses on images of wounded laboratory mice that are used widely for translationally relevant wound studies and leverages a commonly used ring-shaped splint present in most images to predict wound size. We apply the method to a dataset that was never meant to be quantified and, thus, presents many visual challenges. Additionally, the data set was not meant for training deep learning models and so is relatively small in size with only 256 images. We compare results to that of expert measurements and demonstrate preservation of information relevant to predicting wound closure despite variability from machine-to-expert and even expert-to-expert. The proposed system resulted in high fidelity results on unseen data with minimal human intervention. Furthermore, the pipeline estimates acceptable wound sizes when less than 50% of the images are missing reference objects.

All proteomes contain both proteins and polypeptide segments that don't form a defined three-dimensional structure yet are biologically active-called intrinsically disordered proteins and regions (IDPs and IDRs). Most of these IDPs/IDRs lack useful functional annotation limiting our understanding of their importance for organism fitness. Here we characterized IDRs using protein sequence annotations of functional sites and regions available in the UniProt knowledgebase ("UniProt features": active site, ligand-binding pocket, regions mediating protein-protein interactions, etc.). By measuring the statistical enrichment of twenty-five UniProt features in 981 IDRs of 561 human proteins, we identified eight features that are commonly located in IDRs. We then collected the genetic variant data from the general population and patient-based databases and evaluated the prevalence of population and pathogenic variations in IDPs/IDRs. We observed that some IDRs tolerate 2 to 12-times more single amino acid-substituting missense mutations than synonymous changes in the general population. However, we also found that 37% of all germline pathogenic mutations are located in disordered regions of 96 proteins. Based on the observed-to-expected frequency of mutations, we categorized 34 IDRs in 20 proteins (DDX3X, KIT, RB1, etc.) as intolerant to mutation. Finally, using statistical analysis and a machine learning approach, we demonstrate that mutation-intolerant IDRs carry a distinct signature of functional features. Our study presents a novel approach to assign functional importance to IDRs by leveraging the wealth of available genetic data, which will aid in a deeper understating of the role of IDRs in biological processes and disease mechanisms.

Detecting quantitative trait loci (QTL) and estimating QTL variances (represented by the squared QTL effects) are two main goals of QTL mapping and genome-wide association studies (GWAS). However, there are issues associated with estimated QTL variances and such issues have not attracted much attention from the QTL mapping community. Estimated QTL variances are usually biased upwards due to estimation being associated with significance tests. The phenomenon is called the Beavis effect. However, estimated variances of QTL without significance tests can also be biased upwards, which cannot be explained by the Beavis effect; rather, this bias is due to the fact that QTL variances are often estimated as the squares of the estimated QTL effects. The parameters are the QTL effects and the estimated QTL variances are obtained by squaring the estimated QTL effects. This square transformation failed to incorporate the errors of estimated QTL effects into the transformation. The consequence is biases in estimated QTL variances. To correct the biases, we can either reformulate the QTL model by treating the QTL effect as random and directly estimate the QTL variance (as a variance component) or adjust the bias by taking into account the error of the estimated QTL effect. A moment method of estimation has been proposed to correct the bias. The method has been validated via Monte Carlo simulation studies. The method has been applied to QTL mapping for the 10-week-body-weight trait from an F2 mouse population.

At-sea behaviour of seabirds have received significant attention in ecology over the last decades as it is a key process in the ecology and fate of these populations. It is also, through the position of top predator that these species often occupy, a relevant and integrative indicator of the dynamics of the marine ecosystems they rely on. Seabird trajectories are recorded through the deployment of GPS, and a variety of statistical approaches have been tested to infer probable behaviours from these location data. Recently, deep learning tools have shown promising results for the segmentation and classification of animal behaviour from trajectory data. Yet, these approaches have not been widely used and investigation is still needed to identify optimal network architecture and to demonstrate their generalization properties. From a database of about 300 foraging trajectories derived from GPS data deployed simultaneously with pressure sensors for the identification of dives, this work has benchmarked deep neural network architectures trained in a supervised manner for the prediction of dives from trajectory data. It first confirms that deep learning allows better dive prediction than usual methods such as Hidden Markov Models. It also demonstrates the generalization properties of the trained networks for inferring dives distribution for seabirds from other colonies and ecosystems. In particular, convolutional networks trained on Peruvian boobies from a specific colony show great ability to predict dives of boobies from other colonies and from distinct ecosystems. We further investigate accross-species generalization using a transfer learning strategy known as 'fine-tuning'. Starting from a convolutional network pre-trained on Guanay cormorant data reduced by two the size of the dataset needed to accurately predict dives in a tropical booby from Brazil. We believe that the networks trained in this study will provide relevant starting point for future fine-tuning works for seabird trajectory segmentation.

By classifying patients into subgroups, clinicians can provide more effective care than using a uniform approach for all patients. Such subgroups might include patients with a particular disease subtype, patients with a good (or poor) prognosis, or patients most (or least) likely to respond to a particular therapy. Transcriptomic measurements reflect the downstream effects of genomic and epigenomic variations. However, high-throughput technologies generate thousands of measurements per patient, and complex dependencies exist among genes, so it may be infeasible to classify patients using traditional statistical models. Machine-learning classification algorithms can help with this problem. However, hundreds of classification algorithms exist-and most support diverse hyperparameters-so it is difficult for researchers to know which are optimal for gene-expression biomarkers. We performed a benchmark comparison, applying 52 classification algorithms to 50 gene-expression datasets (143 class variables). We evaluated algorithms that represent diverse machine-learning methodologies and have been implemented in general-purpose, open-source, machine-learning libraries. When available, we combined clinical predictors with gene-expression data. Additionally, we evaluated the effects of performing hyperparameter optimization and feature selection using nested cross validation. Kernel- and ensemble-based algorithms consistently outperformed other types of classification algorithms; however, even the top-performing algorithms performed poorly in some cases. Hyperparameter optimization and feature selection typically improved predictive performance, and univariate feature-selection algorithms typically outperformed more sophisticated methods. Together, our findings illustrate that algorithm performance varies considerably when other factors are held constant and thus that algorithm selection is a critical step in biomarker studies.

Annotation of cells in single-cell clustering requires a homogeneous grouping of cell populations. There are various issues in single cell sequencing that effect homogeneous grouping (clustering) of cells, such as small amount of starting RNA, limited per-cell sequenced reads, cell-to-cell variability due to cell-cycle, cellular morphology, and variable reagent concentrations. Moreover, single cell data is susceptible to technical noise, which affects the quality of genes (or features) selected/extracted prior to clustering. Here we introduce sc-CGconv (copula based graph convolution network for single clustering), a stepwise robust unsupervised feature extraction and clustering approach that formulates and aggregates cell-cell relationships using copula correlation (Ccor), followed by a graph convolution network based clustering approach. sc-CGconv formulates a cell-cell graph using Ccor that is learned by a graph-based artificial intelligence model, graph convolution network. The learned representation (low dimensional embedding) is utilized for cell clustering. sc-CGconv features the following advantages. a. sc-CGconv works with substantially smaller sample sizes to identify homogeneous clusters. b. sc-CGconv can model the expression co-variability of a large number of genes, thereby outperforming state-of-the-art gene selection/extraction methods for clustering. c. sc-CGconv preserves the cell-to-cell variability within the selected gene set by constructing a cell-cell graph through copula correlation measure. d. sc-CGconv provides a topology-preserving embedding of cells in low dimensional space.

Cardio/cerebrovascular diseases (CVD) have become one of the major health issue in our societies. But recent studies show that the present pathology tests to detect CVD are ineffectual as they do not consider different stages of platelet activation or the molecular dynamics involved in platelet interactions and are incapable to consider inter-individual variability. Here we propose a stochastic platelet deposition model and an inferential scheme to estimate the biologically meaningful model parameters using approximate Bayesian computation with a summary statistic that maximally discriminates between different types of patients. Inferred parameters from data collected on healthy volunteers and different patient types help us to identify specific biological parameters and hence biological reasoning behind the dysfunction for each type of patients. This work opens up an unprecedented opportunity of personalized pathology test for CVD detection and medical treatment.