In the absence of point-of-care gonorrhea diagnostics that report antibiotic susceptibility, gonorrhea treatment is empiric and determined by standardized guidelines. These guidelines are informed by estimates of resistance prevalence from national surveillance systems. We examined whether guidelines informed by local, rather than national, surveillance data could reduce the incidence of gonorrhea and increase the effective lifespan of antibiotics used in treatment guidelines. We used a transmission dynamic model of gonorrhea among men who have sex with men (MSM) in 16 U.S. metropolitan areas to determine whether spatially adaptive treatment guidelines based on local estimates of resistance prevalence can extend the effective lifespan of hypothetical antibiotics. The rate of gonorrhea cases in these metropolitan areas was 5,548 cases per 100,000 MSM in 2017. Under the current strategy of updating the treatment guideline when the prevalence of resistance exceeds 5%, we showed that spatially adaptive guidelines could reduce the annual rate of gonorrhea cases by 200 cases (95% uncertainty interval: 169, 232) per 100,000 MSM population while extending the use of a first-line antibiotic by 0.75 (0.55, 0.95) years. One potential strategy to reduce the incidence of gonorrhea while extending the effective lifespan of antibiotics is to inform treatment guidelines based on local, rather than national, resistance prevalence.

Mathematical models have come to play a key role in global pandemic preparedness and outbreak response: helping to plan for disease burden, hospital capacity, and inform nonpharmaceutical interventions. Such models have played a pivotal role in the COVID-19 pandemic, with transmission models-and, by consequence, modelers-guiding global, national, and local responses to SARS-CoV-2. However, these models have largely not accounted for the social and structural factors, which lead to socioeconomic, racial, and geographic health disparities. In this piece, we raise and attempt to clarify several questions relating to this important gap in the research and practice of infectious disease modeling: Why do epidemiologic models of emerging infections typically ignore known structural drivers of disparate health outcomes? What have been the consequences of a framework focused primarily on aggregate outcomes on infection equity? What should be done to develop a more holistic approach to modeling-based decision-making during pandemics? In this review, we evaluate potential historical and political explanations for the exclusion of drivers of disparity in infectious disease models for emerging infections, which have often been characterized as "equal opportunity infectors" despite ample evidence to the contrary. We look to examples from other disease systems (HIV, STIs) and successes in including social inequity in models of acute infection transmission as a blueprint for how social connections, environmental, and structural factors can be integrated into a coherent, rigorous, and interpretable modeling framework. We conclude by outlining principles to guide modeling of emerging infections in ways that represent the causes of inequity in infection as central rather than peripheral mechanisms.

Allostery is a form of protein regulation, where ligands that bind sites located apart from the active site can modify the activity of the protein. The molecular mechanisms of allostery have been extensively studied, because allosteric sites are less conserved than active sites, and drugs targeting them are more specific than drugs binding the active sites. Here we quantify the importance of allostery in genetic disease. We show that 1) known allosteric proteins are central in disease networks, contribute to genetic disease and comorbidities much more than non-allosteric proteins, and there is an association between being allosteric and involvement in disease; 2) they are enriched in many major disease types like hematopoietic diseases, cardiovascular diseases, cancers, diabetes, or diseases of the central nervous system; 3) variants from cancer genome-wide association studies are enriched near allosteric proteins, indicating their importance to polygenic traits; and 4) the importance of allosteric proteins in disease is due, at least partly, to their central positions in protein-protein interaction networks, and less due to their dynamical properties.

The ability to aggregate experimental data analysis and results into a concise and interpretable format is a key step in evaluating the success of an experiment. This critical step determines baselines for reproducibility and is a key requirement for data dissemination. However, in practice it can be difficult to consolidate data analyses that encapsulates the broad range of datatypes available in the life sciences. We present STENCIL, a web templating engine designed to organize, visualize, and enable the sharing of interactive data visualizations. STENCIL leverages a flexible web framework for creating templates to render highly customizable visual front ends. This flexibility enables researchers to render small or large sets of experimental outcomes, producing high-quality downloadable and editable figures that retain their original relationship to the source data. REST API based back ends provide programmatic data access and supports easy data sharing. STENCIL is a lightweight tool that can stream data from Galaxy, a popular bioinformatic analysis web platform. STENCIL has been used to support the analysis and dissemination of two large scale genomic projects containing the complete data analysis for over 2,400 distinct datasets. Code and implementation details are available on GitHub: https://github.com/CEGRcode/stencil.

Cortical circuits generate excitatory currents that must be cancelled by strong inhibition to assure stability. The resulting excitatory-inhibitory (E-I) balance can generate spontaneous irregular activity but, in standard balanced E-I models, this requires that an extremely strong feedforward bias current be included along with the recurrent excitation and inhibition. The absence of experimental evidence for such large bias currents inspired us to examine an alternative regime that exhibits asynchronous activity without requiring unrealistically large feedforward input. In these networks, irregular spontaneous activity is supported by a continually changing sparse set of neurons. To support this activity, synaptic strengths must be drawn from high-variance distributions. Unlike standard balanced networks, these sparse balance networks exhibit robust nonlinear responses to uniform inputs and non-Gaussian input statistics. Interestingly, the speed, not the size, of synaptic fluctuations dictates the degree of sparsity in the model. In addition to simulations, we provide a mean-field analysis to illustrate the properties of these networks.

Human immunodeficiency virus type 1 (HIV-1) continues to be a major cause of disease and premature death. As with all viruses, HIV-1 exploits a host cell to replicate. Improving our understanding of the molecular interactions between virus and human host proteins is crucial for a mechanistic understanding of virus biology, infection and host antiviral activities. This knowledge will potentially permit the identification of host molecules for targeting by drugs with antiviral properties. Here, we propose a data-driven approach for the analysis and prediction of the HIV-1 interacting proteins (VIPs) with a focus on the directionality of the interaction: host-dependency versus antiviral factors. Using support vector machine learning models and features encompassing genetic, proteomic and network properties, our results reveal some significant differences between the VIPs and non-HIV-1 interacting human proteins (non-VIPs). As assessed by comparison with the HIV-1 infection pathway data in the Reactome database (sensitivity > 90%, threshold = 0.5), we demonstrate these models have good generalization properties. We find that the 'direction' of the HIV-1-host molecular interactions is also predictable due to different characteristics of 'forward'/pro-viral versus 'backward'/pro-host proteins. Additionally, we infer the previously unknown direction of the interactions between HIV-1 and 1351 human host proteins. A web server for performing predictions is available at http://hivpre.cvr.gla.ac.uk/.

The complexity of biological processes such as cell differentiation is reflected in dynamic transitions between cellular states. Trajectory inference arranges the states into a progression using methodologies propelled by single-cell biology. However, current methods, all returning a best trajectory, do not adequately assess statistical significance of noisy patterns, leading to uncertainty in inferred trajectories. We introduce a tree dimension test for trajectory presence in multivariate data by a dimension measure of Euclidean minimum spanning tree, a test statistic, and a null distribution. Computable in linear time to tree size, the tree dimension measure summarizes the extent of branching more effectively than globally insensitive number of leaves or tree diameter indifferent to secondary branches. The test statistic quantifies trajectory presence and its null distribution is estimated under the null hypothesis of no trajectory in data. On simulated and real single-cell datasets, the test outperformed the intuitive number of leaves and tree diameter statistics. Next, we developed a measure for the tissue specificity of the dynamics of a subset, based on the minimum subtree cover of the subset in a minimum spanning tree. We found that tissue specificity of pathway gene expression dynamics is conserved in human and mouse development: several signal transduction pathways including calcium and Wnt signaling are most tissue specific, while genetic information processing pathways such as ribosome and mismatch repair are least so. Neither the tree dimension test nor the subset specificity measure has any user parameter to tune. Our work opens a window to prioritize cellular dynamics and pathways in development and other multivariate dynamical systems.

Dendrodendritic interactions between excitatory mitral cells and inhibitory granule cells in the olfactory bulb create a dense interaction network, reorganizing sensory representations of odors and, consequently, perception. Large-scale computational models are needed for revealing how the collective behavior of this network emerges from its global architecture. We propose an approach where we summarize anatomical information through dendritic geometry and density distributions which we use to calculate the connection probability between mitral and granule cells, while capturing activity patterns of each cell type in the neural dynamical systems theory of Izhikevich. In this way, we generate an efficient, anatomically and physiologically realistic large-scale model of the olfactory bulb network. Our model reproduces known connectivity between sister vs. non-sister mitral cells; measured patterns of lateral inhibition; and theta, beta, and gamma oscillations. The model in turn predicts testable relationships between network structure and several functional properties, including lateral inhibition, odor pattern decorrelation, and LFP oscillation frequency. We use the model to explore the influence of cortex on the olfactory bulb, demonstrating possible mechanisms by which cortical feedback to mitral cells or granule cells can influence bulbar activity, as well as how neurogenesis can improve bulbar decorrelation without requiring cell death. Our methodology provides a tractable tool for other researchers.

Genome-scale metabolic network reconstructions (GENREs) are valuable tools for understanding microbial metabolism. The process of automatically generating GENREs includes identifying metabolic reactions supported by sufficient genomic evidence to generate a draft metabolic network. The draft GENRE is then gapfilled with additional reactions in order to recapitulate specific growth phenotypes as indicated with associated experimental data. Previous methods have implemented absolute mapping thresholds for the reactions automatically included in draft GENREs; however, there is growing evidence that integrating annotation evidence in a continuous form can improve model accuracy. There is a need for flexibility in the structure of GENREs to better account for uncertainty in biological data, unknown regulatory mechanisms, and context-specificity associated with data inputs. To address this issue, we present a novel method that provides a framework for quantifying combined genomic, biochemical, and phenotypic evidence for each biochemical reaction during automated GENRE construction. Our method, Constraint-based Analysis Yielding reaction Usage across metabolic Networks (CANYUNs), generates accurate GENREs with a quantitative metric for the cumulative evidence for each reaction included in the network. The structuring of CANYUNs allows for the simultaneous integration of three data inputs while maintaining all supporting evidence for biochemical reactions that may be active in an organism. CANYUNs is designed to maximize the utility of experimental and annotation datasets and to ultimately assist in the curation of the reference datasets used for the automatic construction of metabolic networks. We validated CANYUNs by generating an E. coli K-12 model and compared it to the manually curated reconstruction iML1515. Finally, we demonstrated the use of CANYUNs to build a model by generating an E. coli Nissle CANYUNs model using novel phenotypic data that we collected. This method may address key challenges for the procedural construction of metabolic networks by leveraging uncertainty and redundancy in biological data.

The ability to predict human phenotypes and identify biomarkers of disease from metagenomic data is crucial for the development of therapeutics for microbiome-associated diseases. However, metagenomic data is commonly affected by technical variables unrelated to the phenotype of interest, such as sequencing protocol, which can make it difficult to predict phenotype and find biomarkers of disease. Supervised methods to correct for background noise, originally designed for gene expression and RNA-seq data, are commonly applied to microbiome data but may be limited because they cannot account for unmeasured sources of variation. Unsupervised approaches address this issue, but current methods are limited because they are ill-equipped to deal with the unique aspects of microbiome data, which is compositional, highly skewed, and sparse. We perform a comparative analysis of the ability of different denoising transformations in combination with supervised correction methods as well as an unsupervised principal component correction approach that is presently used in other domains but has not been applied to microbiome data to date. We find that the unsupervised principal component correction approach has comparable ability in reducing false discovery of biomarkers as the supervised approaches, with the added benefit of not needing to know the sources of variation apriori. However, in prediction tasks, it appears to only improve prediction when technical variables contribute to the majority of variance in the data. As new and larger metagenomic datasets become increasingly available, background noise correction will become essential for generating reproducible microbiome analyses.

Metabolism is directly and indirectly fine-tuned by a complex web of interacting regulatory mechanisms that fall into two major classes. On the one hand, the expression level of the catalyzing enzyme sets the maximal theoretical flux level (i.e., the net rate of the reaction) for each enzyme-controlled reaction. On the other hand, metabolic regulation controls the metabolic flux through the interactions of metabolites (substrates, cofactors, allosteric modulators) with the responsible enzyme. High-throughput data, such as metabolomics and transcriptomics data, if analyzed separately, do not accurately characterize the hierarchical regulation of metabolism outlined above. They must be integrated to disassemble the interdependence between different regulatory layers controlling metabolism. To this aim, we propose INTEGRATE, a computational pipeline that integrates metabolomics and transcriptomics data, using constraint-based stoichiometric metabolic models as a scaffold. We compute differential reaction expression from transcriptomics data and use constraint-based modeling to predict if the differential expression of metabolic enzymes directly originates differences in metabolic fluxes. In parallel, we use metabolomics to predict how differences in substrate availability translate into differences in metabolic fluxes. We discriminate fluxes regulated at the metabolic and/or gene expression level by intersecting these two output datasets. We demonstrate the pipeline using a set of immortalized normal and cancer breast cell lines. In a clinical setting, knowing the regulatory level at which a given metabolic reaction is controlled will be valuable to inform targeted, truly personalized therapies in cancer patients.

The assembly of proteins into complexes and their interactions with other biomolecules are often vital for their biological function. While it is known that mutations at protein interfaces have a high potential to be damaging and cause human genetic disease, there has been relatively little consideration for how this varies between different types of interfaces. Here we investigate the properties of human pathogenic and putatively benign missense variants at homomeric (isologous and heterologous), heteromeric, DNA, RNA and other ligand interfaces, and at different regions in proteins with respect to those interfaces. We find that different types of interfaces vary greatly in their propensity to be associated with pathogenic mutations, with homomeric heterologous and DNA interfaces being particularly enriched in disease. We also find that residues that do not directly participate in an interface, but are close in three-dimensional space, show a significant disease enrichment. Finally, we observe that mutations at different types of interfaces tend to have distinct property changes when undergoing amino acid substitutions associated with disease, and that this is linked to substantial variability in their identification by computational variant effect predictors.

conceptual representations are critical for human cognition. Despite their importance, key properties of these representations remain poorly understood. Here, we used computational models of distributional semantics to predict multivariate fMRI activity patterns during the activation and contextualization of abstract concepts. We devised a task in which participants had to embed abstract nouns into a story that they developed around a given background context. We found that representations in inferior parietal cortex were predicted by concept similarities emerging in models of distributional semantics. By constructing different model families, we reveal the models' learning trajectories and delineate how abstract and concrete training materials contribute to the formation of brain-like representations. These results inform theories about the format and emergence of abstract conceptual representations in the human brain.

Glutamate transporters preserve the spatial specificity of synaptic transmission by limiting glutamate diffusion away from the synaptic cleft, and prevent excitotoxicity by keeping the extracellular concentration of glutamate at low nanomolar levels. Glutamate transporters are abundantly expressed in astrocytes, and previous estimates have been obtained about their surface expression in astrocytes of the rat hippocampus and cerebellum. Analogous estimates for the mouse hippocampus are currently not available. In this work, we derive the surface density of astrocytic glutamate transporters in mice of different ages via quantitative dot blot. We find that the surface density of glial glutamate transporters is similar in 7-8 week old mice and rats. In mice, the levels of glutamate transporters increase until about 6 months of age and then begin to decline slowly. Our data, obtained from a combination of experimental and modeling approaches, point to the existence of stark differences in the density of expression of glutamate transporters across different sub-cellular compartments, indicating that the extent to which astrocytes limit extrasynaptic glutamate diffusion depends not only on their level of synaptic coverage, but also on the identity of the astrocyte compartment in contact with the synapse. Together, these findings provide information on how heterogeneity in the spatial distribution of glutamate transporters in the plasma membrane of hippocampal astrocytes my alter glutamate receptor activation out of the synaptic cleft.

Third-generation sequencing technologies can generate very long reads with relatively high error rates. The lengths of the reads, which sometimes exceed one million bases, make them invaluable for resolving complex repeats that cannot be assembled using shorter reads. Many high-quality genome assemblies have already been produced, curated, and annotated using the previous generation of sequencing data, and full re-assembly of these genomes with long reads is not always practical or cost-effective. One strategy to upgrade existing assemblies is to generate additional coverage using long-read data, and add that to the previously assembled contigs. SAMBA is a tool that is designed to scaffold and gap-fill existing genome assemblies with additional long-read data, resulting in substantially greater contiguity. SAMBA is the only tool of its kind that also computes and fills in the sequence for all spanned gaps in the scaffolds, yielding much longer contigs. Here we compare SAMBA to several similar tools capable of re-scaffolding assemblies using long-read data, and we show that SAMBA yields better contiguity and introduces fewer errors than competing methods. SAMBA is open-source software that is distributed at https://github.com/alekseyzimin/masurca.

Classical mathematical models of tumor growth have shaped our understanding of cancer and have broad practical implications for treatment scheduling and dosage. However, even the simplest textbook models have been barely validated in real world-data of human patients. In this study, we fitted a range of differential equation models to tumor volume measurements of patients undergoing chemotherapy or cancer immunotherapy for solid tumors. We used a large dataset of 1472 patients with three or more measurements per target lesion, of which 652 patients had six or more data points. We show that the early treatment response shows only moderate correlation with the final treatment response, demonstrating the need for nuanced models. We then perform a head-to-head comparison of six classical models which are widely used in the field: the Exponential, Logistic, Classic Bertalanffy, General Bertalanffy, Classic Gompertz and General Gompertz model. Several models provide a good fit to tumor volume measurements, with the Gompertz model providing the best balance between goodness of fit and number of parameters. Similarly, when fitting to early treatment data, the general Bertalanffy and Gompertz models yield the lowest mean absolute error to forecasted data, indicating that these models could potentially be effective at predicting treatment outcome. In summary, we provide a quantitative benchmark for classical textbook models and state-of-the art models of human tumor growth. We publicly release an anonymized version of our original data, providing the first benchmark set of human tumor growth data for evaluation of mathematical models.

Reactivation of fetal-specific genes and isoforms occurs during heart failure. However, the underlying molecular mechanisms and the extent to which the fetal program switch occurs remains unclear. Limitations hindering transcriptome-wide analyses of alternative splicing differences (i.e. isoform switching) in cardiovascular system (CVS) tissues between fetal, healthy adult and heart failure have included both cellular heterogeneity across bulk RNA-seq samples and limited availability of fetal tissue for research. To overcome these limitations, we have deconvoluted the cellular compositions of 996 RNA-seq samples representing heart failure, healthy adult (heart and arteria), and fetal-like (iPSC-derived cardiovascular progenitor cells) CVS tissues. Comparison of the expression profiles revealed that reactivation of fetal-specific RNA-binding proteins (RBPs), and the accompanied re-expression of 1,523 fetal-specific isoforms, contribute to the transcriptome differences between heart failure and healthy adult heart. Of note, isoforms for 20 different RBPs were among those that reverted in heart failure to the fetal-like expression pattern. We determined that, compared with adult-specific isoforms, fetal-specific isoforms encode proteins that tend to have more functions, are more likely to harbor RBP binding sites, have canonical sequences at their splice sites, and contain typical upstream polypyrimidine tracts. Our study suggests that compared with healthy adult, fetal cardiac tissue requires stricter transcriptional regulation, and that during heart failure reversion to this stricter transcriptional regulation occurs. Furthermore, we provide a resource of cardiac developmental stage-specific and heart failure-associated genes and isoforms, which are largely unexplored and can be exploited to investigate novel therapeutics for heart failure.

Assessment of differential gene expression by qPCR is heavily influenced by the choice of reference genes. Although numerous statistical approaches have been proposed to determine the best reference genes, they can give rise to conflicting results depending on experimental conditions. Hence, recent studies propose the use of RNA-Seq to identify stable genes followed by the application of different statistical approaches to determine the best set of reference genes for qPCR data normalization. In this study, however, we demonstrate that the statistical approach to determine the best reference genes from commonly used conventional candidates is more important than the preselection of 'stable' candidates from RNA-Seq data. Using a qPCR data normalization workflow that we have previously established; we show that qPCR data normalization using conventional reference genes render the same results as stable reference genes selected from RNA-Seq data. We validated these observations in two distinct cross-sectional experimental conditions involving human iPSC derived microglial cells and mouse sciatic nerves. These results taken together show that given a robust statistical approach for reference gene selection, stable genes selected from RNA-Seq data do not offer any significant advantage over commonly used reference genes for normalizing qPCR assays.

Social learning, copying other's behavior without actual experience, offers a cost-effective means of knowledge acquisition. However, it raises the fundamental question of which individuals have reliable information: successful individuals versus the majority. The former and the latter are known respectively as success-based and conformist social learning strategies. We show here that while the success-based strategy fully exploits the benign environment of low uncertainly, it fails in uncertain environments. On the other hand, the conformist strategy can effectively mitigate this adverse effect. Based on these findings, we hypothesized that meta-control of individual and social learning strategies provides effective and sample-efficient learning in volatile and uncertain environments. Simulations on a set of environments with various levels of volatility and uncertainty confirmed our hypothesis. The results imply that meta-control of social learning affords agents the leverage to resolve environmental uncertainty with minimal exploration cost, by exploiting others' learning as an external knowledge base.

Accurate quantification of nerves in cancer specimens is important to understand cancer behaviour. Typically, nerves are manually detected and counted in digitised images of thin tissue sections from excised tumours using immunohistochemistry. However the images are of a large size with nerves having substantial variation in morphology that renders accurate and objective quantification difficult using existing manual and automated counting techniques. Manual counting is precise, but time-consuming, susceptible to inconsistency and has a high rate of false negatives. Existing automated techniques using digitised tissue sections and colour filters are sensitive, however, have a high rate of false positives. In this paper we develop a new automated nerve detection approach, based on a deep learning model with an augmented classification structure. This approach involves pre-processing to extract the image patches for the deep learning model, followed by pixel-level nerve detection utilising the proposed deep learning model. Outcomes assessed were a) sensitivity of the model in detecting manually identified nerves (expert annotations), and b) the precision of additional model-detected nerves. The proposed deep learning model based approach results in a sensitivity of 89% and a precision of 75%. The code and pre-trained model are publicly available at https://github.com/IA92/Automated_Nerves_Quantification.

COVID-19 vaccines have been approved for children of age five and older in many countries. However, there is an ongoing debate as to whether children should be vaccinated and at what priority. In this work, we use mathematical modeling and optimization to study how vaccine allocations to different age groups effect epidemic outcomes. In particular, we consider the effect of extending vaccination campaigns to include the vaccination of children. When vaccine availability is limited, we consider Pareto-optimal allocations with respect to competing measures of the number of infections and mortality and systematically study the trade-offs among them. In the scenarios considered, when some weight is given to the number of infections, we find that it is optimal to allocate vaccines to adolescents in the age group 10-19, even when they are assumed to be less susceptible than adults. We further find that age group 0-9 is included in the optimal allocation for sufficiently high values of the basic reproduction number.

The increasing availability of single-cell RNA-sequencing (scRNA-seq) data from various developmental systems provides the opportunity to infer gene regulatory networks (GRNs) directly from data. Herein we describe IQCELL, a platform to infer, simulate, and study executable logical GRNs directly from scRNA-seq data. Such executable GRNs allow simulation of fundamental hypotheses governing developmental programs and help accelerate the design of strategies to control stem cell fate. We first describe the architecture of IQCELL. Next, we apply IQCELL to scRNA-seq datasets from early mouse T-cell and red blood cell development, and show that the platform can infer overall over 74% of causal gene interactions previously reported from decades of research. We will also show that dynamic simulations of the generated GRN qualitatively recapitulate the effects of known gene perturbations. Finally, we implement an IQCELL gene selection pipeline that allows us to identify candidate genes, without prior knowledge. We demonstrate that GRN simulations based on the inferred set yield results similar to the original curated lists. In summary, the IQCELL platform offers a versatile tool to infer, simulate, and study executable GRNs in dynamic biological systems.

Resource patchiness and aqueous phase fragmentation in soil may induce large differences local growth conditions at submillimeter scales. These are translated to vast differences in bacterial age from cells dividing every thirty minutes in close proximity to plant roots to very old cells experiencing negligible growth in adjacent nutrient poor patches. In this study, we link bacterial population demographics with localized soil and hydration conditions to predict emerging generation time distributions and estimate mean bacterial cell ages using mechanistic and heuristic models of bacterial life in soil. Results show heavy-tailed distributions of generation times that resemble a power law for certain conditions, suggesting that we may find bacterial cells of vastly different ages living side by side within small soil volumes. Our results imply that individual bacteria may exist concurrently with all of their ancestors, resulting in an archive of bacterial cells with traits that have been gained (and lost) throughout time-a feature unique to microbial life. This reservoir of bacterial strains and the potential for the reemergence of rare strains with specific functions may be critical for ecosystem stability and function.

Network-based approaches are becoming increasingly popular for drug discovery as they provide a systems-level overview of the mechanisms underlying disease pathophysiology. They have demonstrated significant early promise over other methods of biological data representation, such as in target discovery, side effect prediction and drug repurposing. In parallel, an explosion of -omics data for the deep characterization of biological systems routinely uncovers molecular signatures of disease for similar applications. Here, we present RPath, a novel algorithm that prioritizes drugs for a given disease by reasoning over causal paths in a knowledge graph (KG), guided by both drug-perturbed as well as disease-specific transcriptomic signatures. First, our approach identifies the causal paths that connect a drug to a particular disease. Next, it reasons over these paths to identify those that correlate with the transcriptional signatures observed in a drug-perturbation experiment, and anti-correlate to signatures observed in the disease of interest. The paths which match this signature profile are then proposed to represent the mechanism of action of the drug. We demonstrate how RPath consistently prioritizes clinically investigated drug-disease pairs on multiple datasets and KGs, achieving better performance over other similar methodologies. Furthermore, we present two case studies showing how one can deconvolute the predictions made by RPath as well as predict novel targets.

Integration of multi-omics data with molecular interaction networks enables elucidation of the pathophysiology of Alzheimer's disease (AD). Using the latest genome-wide association studies (GWAS) including proxy cases and the STRING interactome, we identified an AD network of 142 risk genes and 646 network-proximal genes, many of which were linked to synaptic functions annotated by mouse knockout data. The proximal genes were confirmed to be enriched in a replication GWAS of autopsy-documented cases. By integrating the AD gene network with transcriptomic data of AD and healthy temporal cortices, we identified 17 gene clusters of pathways, such as up-regulated complement activation and lipid metabolism, down-regulated cholinergic activity, and dysregulated RNA metabolism and proteostasis. The relationships among these pathways were further organized by a hierarchy of the AD network pinpointing major parent nodes in graph structure including endocytosis and immune reaction. Control analyses were performed using transcriptomics from cerebellum and a brain-specific interactome. Further integration with cell-specific RNA sequencing data demonstrated genes in our clusters of immunoregulation and complement activation were highly expressed in microglia.

A variational autoencoder (VAE) is a machine learning algorithm, useful for generating a compressed and interpretable latent space. These representations have been generated from various biomedical data types and can be used to produce realistic-looking simulated data. However, standard vanilla VAEs suffer from entangled and uninformative latent spaces, which can be mitigated using other types of VAEs such as β-VAE and MMD-VAE. In this project, we evaluated the ability of VAEs to learn cell morphology characteristics derived from cell images. We trained and evaluated these three VAE variants-Vanilla VAE, β-VAE, and MMD-VAE-on cell morphology readouts and explored the generative capacity of each model to predict compound polypharmacology (the interactions of a drug with more than one target) using an approach called latent space arithmetic (LSA). To test the generalizability of the strategy, we also trained these VAEs using gene expression data of the same compound perturbations and found that gene expression provides complementary information. We found that the β-VAE and MMD-VAE disentangle morphology signals and reveal a more interpretable latent space. We reliably simulated morphology and gene expression readouts from certain compounds thereby predicting cell states perturbed with compounds of known polypharmacology. Inferring cell state for specific drug mechanisms could aid researchers in developing and identifying targeted therapeutics and categorizing off-target effects in the future.

Aging affects almost all aspects of an organism-its morphology, its physiology, its behavior. Isolating which biological mechanisms are regulating these changes, however, has proven difficult, potentially due to our inability to characterize the full repertoire of an animal's behavior across the lifespan. Using data from fruit flies (D. melanogaster) we measure the full repertoire of behaviors as a function of age. We observe a sexually dimorphic pattern of changes in the behavioral repertoire during aging. Although the stereotypy of the behaviors and the complexity of the repertoire overall remains relatively unchanged, we find evidence that the observed alterations in behavior can be explained by changing the fly's overall energy budget, suggesting potential connections between metabolism, aging, and behavior.

In vertebrates, the octopeptide angiotensin II (AngII) is an important in vivo regulator of the cardiovascular system. It acts mainly through two G protein-coupled receptors, AT1 and AT2. To better understand distinctive features of these receptors, we carried out a phylogenetic analysis that revealed a mirror evolution of AT1 and AT2, each one split into two clades, separating fish from terrestrial receptors. It also revealed that hallmark mutations occurred at, or near, the sodium binding site in both AT1 and AT2. Electrostatics computations and molecular dynamics simulations support maintained sodium binding to human AT1 with slow ingress from the extracellular side and an electrostatic component of the binding free energy around -3kT, to be compared to around -2kT for human AT2 and the δ opioid receptor. Comparison of the sodium binding modes in wild type and mutated AT1 and AT2 from humans and eels indicates that the allosteric control by sodium in both AT1 and AT2 evolved during the transition from fish to amniota. The unusual S7.46N mutation in AT1 is mirrored by a L3.36M mutation in AT2. In the presence of sodium, the N7.46 pattern in amniota AT1 stabilizes the inward orientation of N3.35 in the apo receptor, which should contribute to efficient N3.35 driven biased signaling. The M3.36 pattern in amniota AT2 favours the outward orientation of N3.35 and the receptor promiscuity. Both mutations have physiological consequences for the regulation of the renin-angiotensin system.

Studying similarities in protein molecules has become a fundamental activity in much of biology and biomedical research, for which methods such as multiple sequence alignments are widely used. Most methods available for such comparisons cater to studying proteins which have clearly recognizable evolutionary relationships but not to proteins that recognize the same or similar ligands but do not share similarities in their sequence or structural folds. In many cases, proteins in the latter class share structural similarities only in their binding sites. While several algorithms are available for comparing binding sites, there are none for deriving structural motifs of the binding sites, independent of the whole proteins. We report the development of SiteMotif, a new algorithm that compares binding sites from multiple proteins and derives sequence-order independent structural site motifs. We have tested the algorithm at multiple levels of complexity and demonstrate its performance in different scenarios. We have benchmarked against 3 current methods available for binding site comparison and demonstrate superior performance of our algorithm. We show that SiteMotif identifies new structural motifs of spatially conserved residues in proteins, even when there is no sequence or fold-level similarity. We expect SiteMotif to be useful for deriving key mechanistic insights into the mode of ligand interaction, predict the ligand type that a protein can bind and improve the sensitivity of functional annotation.

Bursting is one of the fundamental rhythms that excitable cells can generate either in response to incoming stimuli or intrinsically. It has been a topic of intense research in computational biology for several decades. The classification of bursting oscillations in excitable systems has been the subject of active research since the early 1980s and is still ongoing. As a by-product, it establishes analytical and numerical foundations for studying complex temporal behaviors in multiple timescale models of cellular activity. In this review, we first present the seminal works of Rinzel and Izhikevich in classifying bursting patterns of excitable systems. We recall a complementary mathematical classification approach by Bertram and colleagues, and then by Golubitsky and colleagues, which, together with the Rinzel-Izhikevich proposals, provide the state-of-the-art foundations to these classifications. Beyond classical approaches, we review a recent bursting example that falls outside the previous classification systems. Generalizing this example leads us to propose an extended classification, which requires the analysis of both fast and slow subsystems of an underlying slow-fast model and allows the dissection of a larger class of bursters. Namely, we provide a general framework for bursting systems with both subthreshold and superthreshold oscillations. A new class of bursters with at least 2 slow variables is then added, which we denote folded-node bursters, to convey the idea that the bursts are initiated or annihilated via a folded-node singularity. Key to this mechanism are so-called canard or duck orbits, organizing the underpinning excitability structure. We describe the 2 main families of folded-node bursters, depending upon the phase (active/spiking or silent/nonspiking) of the bursting cycle during which folded-node dynamics occurs. We classify both families and give examples of minimal systems displaying these novel bursting patterns. Finally, we provide a biophysical example by reinterpreting a generic conductance-based episodic burster as a folded-node burster, showing that the associated framework can explain its subthreshold oscillations over a larger parameter region than the fast subsystem approach.

Precise identification of target sites of RNA-binding proteins (RBP) is important to understand their biochemical and cellular functions. A large amount of experimental data is generated by in vivo and in vitro approaches. The binding preferences determined from these platforms share similar patterns but there are discernable differences between these datasets. Computational methods trained on one dataset do not always work well on another dataset. To address this problem which resembles the classic "domain shift" in deep learning, we adopted the adversarial domain adaptation (ADDA) technique and developed a framework (RBP-ADDA) that can extract RBP binding preferences from an integration of in vivo and vitro datasets. Compared with conventional methods, ADDA has the advantage of working with two input datasets, as it trains the initial neural network for each dataset individually, projects the two datasets onto a feature space, and uses an adversarial framework to derive an optimal network that achieves an optimal discriminative predictive power. In the first step, for each RBP, we include only the in vitro data to pre-train a source network and a task predictor. Next, for the same RBP, we initiate the target network by using the source network and use adversarial domain adaptation to update the target network using both in vitro and in vivo data. These two steps help leverage the in vitro data to improve the prediction on in vivo data, which is typically challenging with a lower signal-to-noise ratio. Finally, to further take the advantage of the fused source and target data, we fine-tune the task predictor using both data. We showed that RBP-ADDA achieved better performance in modeling in vivo RBP binding data than other existing methods as judged by Pearson correlations. It also improved predictive performance on in vitro datasets. We further applied augmentation operations on RBPs with less in vivo data to expand the input data and showed that it can improve prediction performances. Lastly, we explored the predictive interpretability of RBP-ADDA, where we quantified the contribution of the input features by Integrated Gradients and identified nucleotide positions that are important for RBP recognition.

Protozoan parasites cause diverse diseases with large global impacts. Research on the pathogenesis and biology of these organisms is limited by economic and experimental constraints. Accordingly, studies of one parasite are frequently extrapolated to infer knowledge about another parasite, across and within genera. Model in vitro or in vivo systems are frequently used to enhance experimental manipulability, but these systems generally use species related to, yet distinct from, the clinically relevant causal pathogen. Characterization of functional differences among parasite species is confined to post hoc or single target studies, limiting the utility of this extrapolation approach. To address this challenge and to accelerate parasitology research broadly, we present a functional comparative analysis of 192 genomes, representing every high-quality, publicly-available protozoan parasite genome including Plasmodium, Toxoplasma, Cryptosporidium, Entamoeba, Trypanosoma, Leishmania, Giardia, and other species. We generated an automated metabolic network reconstruction pipeline optimized for eukaryotic organisms. These metabolic network reconstructions serve as biochemical knowledgebases for each parasite, enabling qualitative and quantitative comparisons of metabolic behavior across parasites. We identified putative differences in gene essentiality and pathway utilization to facilitate the comparison of experimental findings and discovered that phylogeny is not the sole predictor of metabolic similarity. This knowledgebase represents the largest collection of genome-scale metabolic models for both pathogens and eukaryotes; with this resource, we can predict species-specific functions, contextualize experimental results, and optimize selection of experimental systems for fastidious species.

Most normative models in computational neuroscience describe the task of learning as the optimisation of a cost function with respect to a set of parameters. However, learning as optimisation fails to account for a time-varying environment during the learning process and the resulting point estimate in parameter space does not account for uncertainty. Here, we frame learning as filtering, i.e., a principled method for including time and parameter uncertainty. We derive the filtering-based learning rule for a spiking neuronal network-the Synaptic Filter-and show its computational and biological relevance. For the computational relevance, we show that filtering improves the weight estimation performance compared to a gradient learning rule with optimal learning rate. The dynamics of the mean of the Synaptic Filter is consistent with spike-timing dependent plasticity (STDP) while the dynamics of the variance makes novel predictions regarding spike-timing dependent changes of EPSP variability. Moreover, the Synaptic Filter explains experimentally observed negative correlations between homo- and heterosynaptic plasticity.

Emerging evidence suggests that host-microbe interaction in the cervicovaginal microenvironment contributes to cervical carcinogenesis, yet dissecting these complex interactions is challenging. Herein, we performed an integrated analysis of multiple "omics" datasets to develop predictive models of the cervicovaginal microenvironment and identify characteristic features of vaginal microbiome, genital inflammation and disease status. Microbiomes, vaginal pH, immunoproteomes and metabolomes were measured in cervicovaginal specimens collected from a cohort (n = 72) of Arizonan women with or without cervical neoplasm. Multi-omics integration methods, including neural networks (mmvec) and Random Forest supervised learning, were utilized to explore potential interactions and develop predictive models. Our integrated analyses revealed that immune and cancer biomarker concentrations were reliably predicted by Random Forest regressors trained on microbial and metabolic features, suggesting close correspondence between the vaginal microbiome, metabolome, and genital inflammation involved in cervical carcinogenesis. Furthermore, we show that features of the microbiome and host microenvironment, including metabolites, microbial taxa, and immune biomarkers are predictive of genital inflammation status, but only weakly to moderately predictive of cervical neoplastic disease status. Different feature classes were important for prediction of different phenotypes. Lipids (e.g. sphingolipids and long-chain unsaturated fatty acids) were strong predictors of genital inflammation, whereas predictions of vaginal microbiota and vaginal pH relied mostly on alterations in amino acid metabolism. Finally, we identified key immune biomarkers associated with the vaginal microbiota composition and vaginal pH (MIF), as well as genital inflammation (IL-6, IL-10, MIP-1α).

Estimating the changes of epidemiological parameters, such as instantaneous reproduction number, Rt, is important for understanding the transmission dynamics of infectious diseases. Current estimates of time-varying epidemiological parameters often face problems such as lagging observations, averaging inference, and improper quantification of uncertainties. To address these problems, we propose a Bayesian data assimilation framework for time-varying parameter estimation. Specifically, this framework is applied to estimate the instantaneous reproduction number Rt during emerging epidemics, resulting in the state-of-the-art 'DARt' system. With DARt, time misalignment caused by lagging observations is tackled by incorporating observation delays into the joint inference of infections and Rt; the drawback of averaging is overcome by instantaneously updating upon new observations and developing a model selection mechanism that captures abrupt changes; the uncertainty is quantified and reduced by employing Bayesian smoothing. We validate the performance of DARt and demonstrate its power in describing the transmission dynamics of COVID-19. The proposed approach provides a promising solution for making accurate and timely estimation for transmission dynamics based on reported data.

We examine the structure of the visual motion projected on the retina during natural locomotion in real world environments. Bipedal gait generates a complex, rhythmic pattern of head translation and rotation in space, so without gaze stabilization mechanisms such as the vestibular-ocular-reflex (VOR) a walker's visually specified heading would vary dramatically throughout the gait cycle. The act of fixation on stable points in the environment nulls image motion at the fovea, resulting in stable patterns of outflow on the retinae centered on the point of fixation. These outflowing patterns retain a higher order structure that is informative about the stabilized trajectory of the eye through space. We measure this structure by applying the curl and divergence operations on the retinal flow velocity vector fields and found features that may be valuable for the control of locomotion. In particular, the sign and magnitude of foveal curl in retinal flow specifies the body's trajectory relative to the gaze point, while the point of maximum divergence in the retinal flow field specifies the walker's instantaneous overground velocity/momentum vector in retinotopic coordinates. Assuming that walkers can determine the body position relative to gaze direction, these time-varying retinotopic cues for the body's momentum could provide a visual control signal for locomotion over complex terrain. In contrast, the temporal variation of the eye-movement-free, head-centered flow fields is large enough to be problematic for use in steering towards a goal. Consideration of optic flow in the context of real-world locomotion therefore suggests a re-evaluation of the role of optic flow in the control of action during natural behavior.

A critical step in studying biological features (e.g., genetic variants, gene families, metabolic capabilities, or taxa) is assessing their diversity and distribution among a sample of individuals. Accurate assessments of these patterns are essential for linking features to traits or outcomes of interest and understanding their functional impact. Consequently, it is of crucial importance that the measures employed for quantifying feature diversity can perform robustly under any evolutionary scenario. However, the standard measures used for quantifying and comparing the distribution of features, such as prevalence, phylogenetic diversity, and related approaches, either do not take into consideration evolutionary history, or assume strictly vertical patterns of inheritance. Consequently, these approaches cannot accurately assess diversity for features that have undergone recombination or horizontal transfer. To address this issue, we have devised RecPD, a novel recombination-aware phylogenetic-diversity statistic for measuring the distribution and diversity of features under all evolutionary scenarios. RecPD utilizes ancestral-state reconstruction to map the presence / absence of features onto ancestral nodes in a species tree, and then identifies potential recombination events in the evolutionary history of the feature. We also derive several related measures from RecPD that can be used to assess and quantify evolutionary dynamics and correlation of feature evolutionary histories. We used simulation studies to show that RecPD reliably reconstructs feature evolutionary histories under diverse recombination and loss scenarios. We then applied RecPD in two diverse real-world scenarios including a preliminary study type III effector protein families secreted by the plant pathogenic bacterium Pseudomonas syringae and growth phenotypes of the Pseudomonas genus and demonstrate that prevalence is an inadequate measure that obscures the potential impact of recombination. We believe RecPD will have broad utility for revealing and quantifying complex evolutionary processes for features at any biological level.

Highly polymorphic interaction of KIR3DL1 and KIR3DS1 with HLA class I ligands modulates the effector functions of natural killer (NK) cells and some T cells. This genetically determined diversity affects severity of infections, immune-mediated diseases, and some cancers, and impacts the course of immunotherapies, including transplantation. KIR3DL1 is an inhibitory receptor, and KIR3DS1 is an activating receptor encoded by the KIR3DL1/S1 gene that has more than 200 diverse and divergent alleles. Determination of KIR3DL1/S1 genotypes for medical application is hampered by complex sequence and structural variation, requiring targeted approaches to generate and analyze high-resolution allele data. To overcome these obstacles, we developed and optimized a model for imputing KIR3DL1/S1 alleles at high-resolution from whole-genome SNP data. We designed the model to represent a substantial component of human genetic diversity. Our Global imputation model is effective at genotyping KIR3DL1/S1 alleles with an accuracy ranging from 88% in Africans to 97% in East Asians, with mean specificity of 99% and sensitivity of 95% for alleles >1% frequency. We used the established algorithm of the HIBAG program, in a modification named Pulling Out Natural killer cell Genomics (PONG). Because HIBAG was designed to impute HLA alleles also from whole-genome SNP data, PONG allows combinatorial diversity of KIR3DL1/S1 with HLA-A and -B to be analyzed using complementary techniques on a single data source. The use of PONG thus negates the need for targeted sequencing data in very large-scale association studies where such methods might not be tractable.

Investigation of the diversity of malaria parasite antigens can help prioritize and validate them as vaccine candidates and identify the most common variants for inclusion in vaccine formulations. Studies of vaccine candidates of the most virulent human malaria parasite, Plasmodium falciparum, have focused on a handful of well-known antigens, while several others have never been studied. Here we examine the global diversity and population structure of leading vaccine candidate antigens of P. falciparum using the MalariaGEN Pf3K (version 5.1) resource, comprising more than 2600 genomes from 15 malaria endemic countries. A stringent variant calling pipeline was used to extract high quality antigen gene 'haplotypes' from the global dataset and a new R-package named VaxPack was used to streamline population genetic analyses. In addition, a newly developed algorithm that enables spatial averaging of selection pressure on 3D protein structures was applied to the dataset. We analysed the genes encoding 23 leading and novel candidate malaria vaccine antigens including csp, trap, eba175, ama1, rh5, and CelTOS. Our analysis shows that current malaria vaccine formulations are based on rare haplotypes and thus may have limited efficacy against natural parasite populations. High levels of diversity with evidence of balancing selection was detected for most of the erythrocytic and pre-erythrocytic antigens. Measures of natural selection were then mapped to 3D protein structures to predict targets of functional antibodies. For some antigens, geographical variation in the intensity and distribution of these signals on the 3D structure suggests adaptation to different human host or mosquito vector populations. This study provides an essential framework for the diversity of P. falciparum antigens to be considered in the design of the next generation of malaria vaccines.

The ammonia-oxidizing bacterium Nitrosomonas europaea has been widely recognized as an important player in the nitrogen cycle as well as one of the most abundant members in microbial communities for the treatment of industrial or sewage wastewater. Its natural metabolic versatility and extraordinary ability to degrade environmental pollutants (e.g., aromatic hydrocarbons such as benzene and toluene) enable it to thrive under various harsh environmental conditions. Constraint-based metabolic models constructed from genome sequences enable quantitative insight into the central and specialized metabolism within a target organism. These genome-scale models have been utilized to understand, optimize, and design new strategies for improved bioprocesses. Reduced modeling approaches have been used to elucidate Nitrosomonas europaea metabolism at a pathway level. However, genome-scale knowledge about the simultaneous oxidation of ammonia and pollutant metabolism of N. europaea remains limited. Here, we describe the reconstruction, manual curation, and validation of the genome-scale metabolic model for N. europaea, iGC535. This reconstruction is the most accurate metabolic model for a nitrifying organism to date, reaching an average prediction accuracy of over 90% under several growth conditions. The manually curated model can predict phenotypes under chemolithotrophic and chemolithoorganotrophic conditions while oxidating methane and wastewater pollutants. Calculated flux distributions under different trophic conditions show that several key pathways are affected by the type of carbon source available, including central carbon metabolism and energy production.

Adaptive sequential behavior is a hallmark of human cognition. In particular, humans can learn to produce precise spatiotemporal sequences given a certain context. For instance, musicians can not only reproduce learned action sequences in a context-dependent manner, they can also quickly and flexibly reapply them in any desired tempo or rhythm without overwriting previous learning. Existing neural network models fail to account for these properties. We argue that this limitation emerges from the fact that sequence information (i.e., the position of the action) and timing (i.e., the moment of response execution) are typically stored in the same neural network weights. Here, we augment a biologically plausible recurrent neural network of cortical dynamics to include a basal ganglia-thalamic module which uses reinforcement learning to dynamically modulate action. This "associative cluster-dependent chain" (ACDC) model modularly stores sequence and timing information in distinct loci of the network. This feature increases computational power and allows ACDC to display a wide range of temporal properties (e.g., multiple sequences, temporal shifting, rescaling, and compositionality), while still accounting for several behavioral and neurophysiological empirical observations. Finally, we apply this ACDC network to show how it can learn the famous "Thunderstruck" song intro and then flexibly play it in a "bossa nova" rhythm without further training.

Biomolecular condensates formed via liquid-liquid phase separation (LLPS) play a crucial role in the spatiotemporal organization of the cell material. Nucleic acids can act as critical modulators in the stability of these protein condensates. To unveil the role of RNA length in regulating the stability of RNA binding protein (RBP) condensates, we present a multiscale computational strategy that exploits the advantages of a sequence-dependent coarse-grained representation of proteins and a minimal coarse-grained model wherein proteins are described as patchy colloids. We find that for a constant nucleotide/protein ratio, the protein fused in sarcoma (FUS), which can phase separate on its own-i.e., via homotypic interactions-only exhibits a mild dependency on the RNA strand length. In contrast, the 25-repeat proline-arginine peptide (PR25), which does not undergo LLPS on its own at physiological conditions but instead exhibits complex coacervation with RNA-i.e., via heterotypic interactions-shows a strong dependence on the length of the RNA strands. Our minimal patchy particle simulations suggest that the strikingly different effect of RNA length on homotypic LLPS versus RBP-RNA complex coacervation is general. Phase separation is RNA-length dependent whenever the relative contribution of heterotypic interactions sustaining LLPS is comparable or higher than those stemming from protein homotypic interactions. Taken together, our results contribute to illuminate the intricate physicochemical mechanisms that influence the stability of RBP condensates through RNA inclusion.

Fitness conferred by the same allele may differ between genotypes and environments, and these differences shape variation and evolution. Changes in amino acid propensities at protein sites over the course of evolution have been inferred from sequence alignments statistically, but the existing methods are data-intensive and aggregate multiple sites. Here, we develop an approach to detect individual amino acids that confer different fitness in different groups of species from combined sequence and phylogenetic data. Using the fact that the probability of a substitution to an amino acid depends on its fitness, our method looks for amino acids such that substitutions to them occur more frequently in one group of lineages than in another. We validate our method using simulated evolution of a protein site under different scenarios and show that it has high specificity for a wide range of assumptions regarding the underlying changes in selection, while its sensitivity differs between scenarios. We apply our method to the env gene of two HIV-1 subtypes, A and B, and to the HA gene of two influenza A subtypes, H1 and H3, and show that the inferred fitness changes are consistent with the fitness differences observed in deep mutational scanning experiments. We find that changes in relative fitness of different amino acid variants within a site do not always trigger episodes of positive selection and therefore may not result in an overall increase in the frequency of substitutions, but can still be detected from changes in relative frequencies of different substitutions.

Biliary ducts collect bile from liver lobules, the smallest functional and anatomical units of liver, and carry it to the gallbladder. Disruptions in this process caused by defective embryonic development, or through ductal reaction in liver disease have a major impact on life quality and survival of patients. A deep understanding of the processes underlying bile duct lumen formation is crucial to identify intervention points to avoid or treat the appearance of defective bile ducts. Several hypotheses have been proposed to characterize the biophysical mechanisms driving initial bile duct lumen formation during embryogenesis. Here, guided by the quantification of morphological features and expression of genes in bile ducts from embryonic mouse liver, we sharpened these hypotheses and collected data to develop a high resolution individual cell-based computational model that enables to test alternative hypotheses in silico. This model permits realistic simulations of tissue and cell mechanics at sub-cellular scale. Our simulations suggest that successful bile duct lumen formation requires a simultaneous contribution of directed cell division of cholangiocytes, local osmotic effects generated by salt excretion in the lumen, and temporally-controlled differentiation of hepatoblasts to cholangiocytes, with apical constriction of cholangiocytes only moderately affecting luminal size.

In-silico methods for the prediction of epitopes can support and improve workflows for vaccine design, antibody production, and disease therapy. So far, the scope of B cell and T cell epitope prediction has been directed exclusively towards peptidic antigens. Nevertheless, various non-peptidic molecular classes can be recognized by immune cells. These compounds have not been systematically studied yet, and prediction approaches are lacking. The ability to predict the epitope activity of non-peptidic compounds could have vast implications; for example, for immunogenic risk assessment of the vast number of drugs and other xenobiotics. Here we present the first general attempt to predict the epitope activity of non-peptidic compounds using the Immune Epitope Database (IEDB) as a source for positive samples. The molecules stored in the Chemical Entities of Biological Interest (ChEBI) database were chosen as background samples. The molecules were clustered into eight homogeneous molecular groups, and classifiers were built for each cluster with the aim of separating the epitopes from the background. Different molecular feature encoding schemes and machine learning models were compared against each other. For those models where a high performance could be achieved based on simple decision rules, the molecular features were then further investigated. Additionally, the findings were used to build a web server that allows for the immunogenic investigation of non-peptidic molecules (http://tools-staging.iedb.org/np_epitope_predictor). The prediction quality was tested with samples from independent evaluation datasets, and the implemented method received noteworthy Receiver Operating Characteristic-Area Under Curve (ROC-AUC) values, ranging from 0.69-0.96 depending on the molecule cluster.

Spider venom GDPD-like phospholipases D (SicTox) have been identified to be one of the major toxins in recluse spider venom. They are divided into two major clades: the α clade and the β clade. Most α clade toxins present high activity against lipids with choline head groups such as sphingomyelin, while activities in β clade toxins vary and include preference for substrates containing ethanolamine headgroups (Sicarius terrosus, St_βIB1). A structural comparison of available structures of phospholipases D (PLDs) reveals a conserved aromatic cage in the α clade. To test the potential influence of the aromatic cage on membrane-lipid specificity we performed molecular dynamics (MD) simulations of the binding of several PLDs onto lipid bilayers containing choline headgroups; two SicTox from the α clade, Loxosceles intermedia αIA1 (Li_αIA) and Loxosceles laeta αIII1 (Ll_αIII1), and one from the β clade, St_βIB1. The simulation results reveal that the aromatic cage captures a choline-headgroup and suggest that the cage plays a major role in lipid specificity. We also simulated an engineered St_βIB1, where we introduced the aromatic cage, and this led to binding with choline-containing lipids. Moreover, a multiple sequence alignment revealed the conservation of the aromatic cage among the α clade PLDs. Here, we confirmed that the i-face of α and β clade PLDs is involved in their binding to choline and ethanolamine-containing bilayers, respectively. Furthermore, our results suggest a major role in choline lipid recognition of the aromatic cage of the α clade PLDs. The MD simulation results are supported by in vitro liposome binding assay experiments.

In eukaryotes, the cell volume is observed to be strongly correlated with the nuclear volume. The slope of this correlation depends on the cell type, growth condition, and the physical environment of the cell. We develop a computational model of cell growth and proteome increase, incorporating the kinetics of amino acid import, protein/ribosome synthesis and degradation, and active transport of proteins between the cytoplasm and the nucleoplasm. We also include a simple model of ribosome biogenesis and assembly. Results show that the cell volume is tightly correlated with the nuclear volume, and the cytoplasm-nucleoplasm transport rates strongly influence the cell growth rate as well as the cell/nucleus volume ratio (C/N ratio). Ribosome assembly and the ratio of ribosomal proteins to mature ribosomes also influence the cell volume and the cell growth rate. We find that in order to regulate the cell growth rate and the cell/nucleus volume ratio, the cell must optimally control groups of kinetic and transport parameters together, which could explain the quantitative roles of canonical growth pathways. Finally, although not explicitly demonstrated in this work, we point out that it is possible to construct a detailed proteome distribution using our model and RNAseq data, provided that a quantitative cell division mechanism is known.

Collective living systems regularly achieve cooperative emergent functions that individual organisms could not accomplish alone. The rafts of fire ants (Solenopsis invicta) are often studied in this context for their ability to create aggregated structures comprised entirely of their own bodies, including tether-like protrusions that facilitate exploration of and escape from flooded environments. While similar protrusions are observed in cytoskeletons and cellular aggregates, they are generally dependent on morphogens or external gradients leaving the isolated role of local interactions poorly understood. Here we demonstrate through an ant-inspired, agent-based numerical model how protrusions in ant rafts may emerge spontaneously due to local interactions. The model is comprised of a condensed structural network of agents that represents the monolayer of interconnected worker ants, which floats on the water and gives ant rafts their form. Experimentally, this layer perpetually contracts, which we capture through the pairwise contraction of all neighboring structural agents at a strain rate of [Formula: see text]. On top of the structural layer, we model a dispersed, on-lattice layer of motile agents that represents free ants, which walk on top of the floating network. Experimentally, these self-propelled free ants walk with some mean persistence length and speed that we capture through an ant-inspired phenomenological model. Local interactions occur between neighboring free ants within some radius of detection, R, and the persistence length of freely active agents is tuned through a noise parameter, η as introduced by the Vicsek model. Both R and η where fixed to match the experimental trajectories of free ants. Treadmilling of the raft occurs as agents transition between the structural and free layers in accordance with experimental observations. Ultimately, we demonstrate how phases of exploratory protrusion growth may be induced by increased ant activity as characterized by a dimensionless parameter, [Formula: see text]. These results provide an example in which functional morphogenesis of a living system may emerge purely from local interactions at the constituent length scale, thereby providing a source of inspiration for the development of decentralized, autonomous active matter and swarm robotics.

During cell devision, maintaining the epigenetic information encoded in histone modification patterns is crucial for survival and identity of cells. The faithful inheritance of the histone marks from the parental to the daughter strands is a puzzle, given that each strand gets only half of the parental nucleosomes. Mapping DNA replication and reconstruction of modifications to equivalent problems in communication of information, we ask how well enzymes can recover the parental modifications, if they were ideal computing machines. Studying a parameter regime where realistic enzymes can function, our analysis predicts that enzymes may implement a critical threshold filling algorithm which fills unmodified regions of length at most k. This algorithm, motivated from communication theory, is derived from the maximum à posteriori probability (MAP) decoding which identifies the most probable modification sequence based on available observations. Simulations using our method produce modification patterns similar to what has been observed in recent experiments. We also show that our results can be naturally extended to explain inheritance of spatially distinct antagonistic modifications.

Single-cell RNA sequencing (scRNA-seq) methods are typically unable to quantify the expression levels of all genes in a cell, creating a need for the computational prediction of missing values ('dropout imputation'). Most existing dropout imputation methods are limited in the sense that they exclusively use the scRNA-seq dataset at hand and do not exploit external gene-gene relationship information. Further, it is unknown if all genes equally benefit from imputation or which imputation method works best for a given gene. Here, we show that a transcriptional regulatory network learned from external, independent gene expression data improves dropout imputation. Using a variety of human scRNA-seq datasets we demonstrate that our network-based approach outperforms published state-of-the-art methods. The network-based approach performs particularly well for lowly expressed genes, including cell-type-specific transcriptional regulators. Further, the cell-to-cell variation of 11.3% to 48.8% of the genes could not be adequately imputed by any of the methods that we tested. In those cases gene expression levels were best predicted by the mean expression across all cells, i.e. assuming no measurable expression variation between cells. These findings suggest that different imputation methods are optimal for different genes. We thus implemented an R-package called ADImpute (available via Bioconductor https://bioconductor.org/packages/release/bioc/html/ADImpute.html) that automatically determines the best imputation method for each gene in a dataset. Our work represents a paradigm shift by demonstrating that there is no single best imputation method. Instead, we propose that imputation should maximally exploit external information and be adapted to gene-specific features, such as expression level and expression variation across cells.

Alcohol Use Disorder (AUD) is a chronic, relapsing syndrome diagnosed by a heterogeneous set of behavioral signs and symptoms. There are no laboratory tests that provide direct objective evidence for diagnosis. Microarray and RNA-Seq technologies enable genome-wide transcriptome profiling at low costs and provide an opportunity to identify biomarkers to facilitate diagnosis, prognosis, and treatment of patients. However, access to brain tissue in living patients is not possible. Blood contains cellular and extracellular RNAs that provide disease-relevant information for some brain diseases. We hypothesized that blood gene expression profiles can be used to diagnose AUD. We profiled brain (prefrontal cortex, amygdala, and hypothalamus) and blood gene expression levels in C57BL/6J mice using RNA-seq one week after chronic intermittent ethanol (CIE) exposure, a mouse model of alcohol dependence. We found a high degree of preservation (rho range: [0.50, 0.67]) between blood and brain transcript levels. There was small overlap between blood and brain DEGs, and considerable overlap of gene networks perturbed after CIE related to cell-cell signaling (e.g., GABA and glutamate receptor signaling), immune responses (e.g., antigen presentation), and protein processing / mitochondrial functioning (e.g., ubiquitination, oxidative phosphorylation). Blood gene expression data were used to train classifiers (logistic regression, random forest, and partial least squares discriminant analysis), which were highly accurate at predicting alcohol dependence status (maximum AUC: 90.1%). These results suggest that gene expression profiles from peripheral blood samples contain a biological signature of alcohol dependence that can discriminate between CIE and Air subjects.

Insertions and deletions (indels) in human genomes are associated with a wide range of phenotypes, including various clinical disorders. High-throughput, next generation sequencing (NGS) technologies enable the detection of short genetic variants, such as single nucleotide variants (SNVs) and indels. However, the variant calling accuracy for indels remains considerably lower than for SNVs. Here we present a comparative study of the performance of variant calling tools for indel calling, evaluated with a wide repertoire of NGS datasets. While there is no single optimal tool to suit all circumstances, our results demonstrate that the choice of variant calling tool greatly impacts the precision and recall of indel calling. Furthermore, to reliably detect indels, it is essential to choose NGS technologies that offer a long read length and high coverage coupled with specific variant calling tools.

Hippocampal sharp wave/ripple oscillations are a prominent pattern of collective activity, which consists of a strong overall increase of activity with superimposed (140 - 200 Hz) ripple oscillations. Despite its prominence and its experimentally demonstrated importance for memory consolidation, the mechanisms underlying its generation are to date not understood. Several models assume that recurrent networks of inhibitory cells alone can explain the generation and main characteristics of the ripple oscillations. Recent experiments, however, indicate that in addition to inhibitory basket cells, the pattern requires in vivo the activity of the local population of excitatory pyramidal cells. Here, we study a model for networks in the hippocampal region CA1 incorporating such a local excitatory population of pyramidal neurons. We start by investigating its ability to generate ripple oscillations using extensive simulations. Using biologically plausible parameters, we find that short pulses of external excitation triggering excitatory cell spiking are required for sharp/wave ripple generation with oscillation patterns similar to in vivo observations. Our model has plausible values for single neuron, synapse and connectivity parameters, random connectivity and no strong feedforward drive to the inhibitory population. Specifically, whereas temporally broad excitation can lead to high-frequency oscillations in the ripple range, sparse pyramidal cell activity is only obtained with pulse-like external CA3 excitation. Further simulations indicate that such short pulses could originate from dendritic spikes in the apical or basal dendrites of CA1 pyramidal cells, which are triggered by coincident spike arrivals from hippocampal region CA3. Finally we show that replay of sequences by pyramidal neurons and ripple oscillations can arise intrinsically in CA1 due to structured connectivity that gives rise to alternating excitatory pulse and inhibitory gap coding; the latter denotes phases of silence in specific basket cell groups, which induce selective disinhibition of groups of pyramidal neurons. This general mechanism for sequence generation leads to sparse pyramidal cell and dense basket cell spiking, does not rely on synfire chain-like feedforward excitation and may be relevant for other brain regions as well.

The disease burden attributable to opportunistic pathogens depends on their prevalence in asymptomatic colonisation and the rate at which they progress to cause symptomatic disease. Increases in infections caused by commensals can result from the emergence of "hyperinvasive" strains. Such pathogens can be identified through quantifying progression rates using matched samples of typed microbes from disease cases and healthy carriers. This study describes Bayesian models for analysing such datasets, implemented in an RStan package (https://github.com/nickjcroucher/progressionEstimation). The models converged on stable fits that accurately reproduced observations from meta-analyses of Streptococcus pneumoniae datasets. The estimates of invasiveness, the progression rate from carriage to invasive disease, in cases per carrier per year correlated strongly with the dimensionless values from meta-analysis of odds ratios when sample sizes were large. At smaller sample sizes, the Bayesian models produced more informative estimates. This identified historically rare but high-risk S. pneumoniae serotypes that could be problematic following vaccine-associated disruption of the bacterial population. The package allows for hypothesis testing through model comparisons with Bayes factors. Application to datasets in which strain and serotype information were available for S. pneumoniae found significant evidence for within-strain and within-serotype variation in invasiveness. The heterogeneous geographical distribution of these genotypes is therefore likely to contribute to differences in the impact of vaccination in between locations. Hence genomic surveillance of opportunistic pathogens is crucial for quantifying the effectiveness of public health interventions, and enabling ongoing meta-analyses that can identify new, highly invasive variants.

Tick paralysis resulting from bites from Ixodes holocyclus and I. cornuatus is one of the leading causes of emergency veterinary admissions for companion animals in Australia, often resulting in death if left untreated. Availability of timely information on periods of increased risk can help modulate behaviors that reduce exposures to ticks and improve awareness of owners for the need of lifesaving preventative ectoparasite treatment. Improved awareness of clinicians and pet owners about temporal changes in tick paralysis risk can be assisted by ecological forecasting frameworks that integrate environmental information into statistical time series models. Using an 11-year time series of tick paralysis cases from veterinary clinics in one of Australia's hotspots for the paralysis tick Ixodes holocyclus, we asked whether an ensemble model could accurately forecast clinical caseloads over near-term horizons. We fit a series of statistical time series (ARIMA, GARCH) and generative models (Prophet, Generalised Additive Model) using environmental variables as predictors, and then combined forecasts into a weighted ensemble to minimise prediction interval error. Our results indicate that variables related to temperature anomalies, levels of vegetation moisture and the Southern Oscillation Index can be useful for predicting tick paralysis admissions. Our model forecasted tick paralysis cases with exceptional accuracy while preserving epidemiological interpretability, outperforming a field-leading benchmark Exponential Smoothing model by reducing both point and prediction interval errors. Using online particle filtering to assimilate new observations and adjust forecast distributions when new data became available, our model adapted to changing temporal conditions and provided further reduced forecast errors. We expect our model pipeline to act as a platform for developing early warning systems that can notify clinicians and pet owners about heightened risks of environmentally driven veterinary conditions.

The dynamics of a spreading disease and individual behavioral changes are entangled processes that have to be addressed together in order to effectively manage an outbreak. Here, we relate individual risk perception to the adoption of a specific set of control measures, as obtained from an extensive large-scale survey performed via Facebook-involving more than 500,000 respondents from 64 countries-showing that there is a "one-to-one" relationship between perceived epidemic risk and compliance with a set of mitigation rules. We then develop a mathematical model for the spreading of a disease-sharing epidemiological features with COVID-19-that explicitly takes into account non-compliant individual behaviors and evaluates the impact of a population fraction of infectious risk-deniers on the epidemic dynamics. Our modeling study grounds on a wide set of structures, including both synthetic and more than 180 real-world contact patterns, to evaluate, in realistic scenarios, how network features typical of human interaction patterns impact the spread of a disease. In both synthetic and real contact patterns we find that epidemic spreading is hindered for decreasing population fractions of risk-denier individuals. From empirical contact patterns we demonstrate that connectivity heterogeneity and group structure significantly affect the peak of hospitalized population: higher modularity and heterogeneity of social contacts are linked to lower peaks at a fixed fraction of risk-denier individuals while, at the same time, such features increase the relative impact on hospitalizations with respect to the case where everyone correctly perceive the risks.

A central problem in developmental and synthetic biology is understanding the mechanisms by which cells in a tissue or a Petri dish process external cues and transform such information into a coherent response, e.g., a terminal differentiation state. It was long believed that this type of positional information could be entirely attributed to a gradient of concentration of a specific signaling molecule (i.e., a morphogen). However, advances in experimental methodologies and computer modeling have demonstrated the crucial role of the dynamics of a cell's gene regulatory network (GRN) in decoding the information carried by the morphogen, which is eventually translated into a spatial pattern. This morphogen interpretation mechanism has gained much attention in systems biology as a tractable system to investigate the emergent properties of complex genotype-phenotype maps. In this study, we apply a Markov chain Monte Carlo (MCMC)-like algorithm to probe the design space of three-node GRNs with the ability to generate a band-like expression pattern (target phenotype) in the middle of an arrangement of 30 cells, which resemble a simple (1-D) morphogenetic field in a developing embryo. Unlike most modeling studies published so far, here we explore the space of GRN topologies with nodes having the potential to perceive the same input signal differently. This allows for a lot more flexibility during the search space process, and thus enables us to identify a larger set of potentially interesting and realizable morphogen interpretation mechanisms. Out of 2061 GRNs selected using the search space algorithm, we found 714 classes of network topologies that could correctly interpret the morphogen. Notably, the main network motif that generated the target phenotype in response to the input signal was the type 3 Incoherent Feed-Forward Loop (I3-FFL), which agrees with previous theoretical expectations and experimental observations. Particularly, compared to a previously reported pattern forming GRN topologies, we have uncovered a great variety of novel network designs, some of which might be worth inquiring through synthetic biology methodologies to test for the ability of network design with minimal regulatory complexity to interpret a developmental cue robustly.

Lymphocytes have been described to perform different motility patterns such as Brownian random walks, persistent random walks, and Lévy walks. Depending on the conditions, such as confinement or the distribution of target cells, either Brownian or Lévy walks lead to more efficient interaction with the targets. The diversity of these motility patterns may be explained by an adaptive response to the surrounding extracellular matrix (ECM). Indeed, depending on the ECM composition, lymphocytes either display a floating motility without attaching to the ECM, or sliding and stepping motility with respectively continuous or discontinuous attachment to the ECM, or pivoting behaviour with sustained attachment to the ECM. Moreover, on the long term, lymphocytes either perform a persistent random walk or a Brownian-like movement depending on the ECM composition. How the ECM affects cell motility is still incompletely understood. Here, we integrate essential mechanistic details of the lymphocyte-matrix adhesions and lymphocyte intrinsic cytoskeletal induced cell propulsion into a Cellular Potts model (CPM). We show that the combination of de novo cell-matrix adhesion formation, adhesion growth and shrinkage, adhesion rupture, and feedback of adhesions onto cell propulsion recapitulates multiple lymphocyte behaviours, for different lymphocyte subsets and various substrates. With an increasing attachment area and increased adhesion strength, the cells' speed and persistence decreases. Additionally, the model predicts random walks with short-term persistent but long-term subdiffusive properties resulting in a pivoting type of motility. For small adhesion areas, the spatial distribution of adhesions emerges as a key factor influencing cell motility. Small adhesions at the front allow for more persistent motility than larger clusters at the back, despite a similar total adhesion area. In conclusion, we present an integrated framework to simulate the effects of ECM proteins on cell-matrix adhesion dynamics. The model reveals a sufficient set of principles explaining the plasticity of lymphocyte motility.

Supervised machine learning applications in health care are often limited due to a scarcity of labeled training data. To mitigate the effect of small sample size, we introduce a pre-training approach, Patient Contrastive Learning of Representations (PCLR), which creates latent representations of electrocardiograms (ECGs) from a large number of unlabeled examples using contrastive learning. The resulting representations are expressive, performant, and practical across a wide spectrum of clinical tasks. We develop PCLR using a large health care system with over 3.2 million 12-lead ECGs and demonstrate that training linear models on PCLR representations achieves a 51% performance increase, on average, over six training set sizes and four tasks (sex classification, age regression, and the detection of left ventricular hypertrophy and atrial fibrillation), relative to training neural network models from scratch. We also compared PCLR to three other ECG pre-training approaches (supervised pre-training, unsupervised pre-training with an autoencoder, and pre-training using a contrastive multi ECG-segment approach), and show significant performance benefits in three out of four tasks. We found an average performance benefit of 47% over the other models and an average of a 9% performance benefit compared to best model for each task. We release PCLR to enable others to extract ECG representations at https://github.com/broadinstitute/ml4h/tree/master/model_zoo/PCLR.

Visual inspection of stimulus-induced gamma oscillations (30-70 Hz) often reveals a non-sinusoidal shape. Such distortions are a hallmark of non-linear systems and are also observed in mean-field models of gamma oscillations. A thorough characterization of the shape of the gamma cycle can therefore provide additional constraints on the operating regime of such models. However, the gamma waveform has not been quantitatively characterized, partially because the first harmonic of gamma, which arises because of the non-sinusoidal nature of the signal, is typically weak and gets masked due to a broadband increase in power related to spiking. To address this, we recorded local field potential (LFP) from the primary visual cortex (V1) of two awake female macaques while presenting full-field gratings or iso-luminant chromatic hues that produced huge gamma oscillations with prominent peaks at harmonic frequencies in the power spectra. We found that gamma and its first harmonic always maintained a specific phase relationship, resulting in a distinctive shape with a sharp trough and a shallow peak. Interestingly, a Wilson-Cowan (WC) model operating in an inhibition stabilized mode could replicate this shape, but only when the inhibitory population operated in the super-linear regime, as predicted recently. However, another recently developed model of gamma that operates in a linear regime driven by stochastic noise failed to produce salient harmonics or the observed shape. Our results impose additional constraints on models that generate gamma oscillations and their operating regimes.

As sequence and structure comparison algorithms gain sensitivity, the intrinsic interconnectedness of the protein universe has become increasingly apparent. Despite this general trend, β-trefoils have emerged as an uncommon counterexample: They are an isolated protein lineage for which few, if any, sequence or structure associations to other lineages have been identified. If β-trefoils are, in fact, remote islands in sequence-structure space, it implies that the oligomerizing peptide that founded the β-trefoil lineage itself arose de novo. To better understand β-trefoil evolution, and to probe the limits of fragment sharing across the protein universe, we identified both 'β-trefoil bridging themes' (evolutionarily-related sequence segments) and 'β-trefoil-like motifs' (structure motifs with a hallmark feature of the β-trefoil architecture) in multiple, ostensibly unrelated, protein lineages. The success of the present approach stems, in part, from considering β-trefoil sequence segments or structure motifs rather than the β-trefoil architecture as a whole, as has been done previously. The newly uncovered inter-lineage connections presented here suggest a novel hypothesis about the origins of the β-trefoil fold itself-namely, that it is a derived fold formed by 'budding' from an Immunoglobulin-like β-sandwich protein. These results demonstrate how the evolution of a folded domain from a peptide need not be a signature of antiquity and underpin an emerging truth: few protein lineages escape nature's sewing table.

Gmxapi provides an integrated, native Python API for both standard and advanced molecular dynamics simulations in GROMACS. The Python interface permits multiple levels of integration with the core GROMACS libraries, and legacy support is provided via an interface that mimics the command-line syntax, so that all GROMACS commands are fully available. Gmxapi has been officially supported since the GROMACS 2019 release and is enabled by default in current versions of the software. Here we describe gmxapi 0.3 and later. Beyond simply wrapping GROMACS library operations, the API permits several advanced operations that are not feasible using the prior command-line interface. First, the API allows custom user plugin code within the molecular dynamics force calculations, so users can execute custom algorithms without modifying the GROMACS source. Second, the Python interface allows tasks to be dynamically defined, so high-level algorithms for molecular dynamics simulation and analysis can be coordinated with loop and conditional operations. Gmxapi makes GROMACS more accessible to custom Python scripting while also providing support for high-level data-flow simulation algorithms that were previously feasible only in external packages.

Inferring the dynamics of pathogen transmission during an outbreak is an important problem in infectious disease epidemiology. In mathematical epidemiology, estimates are often informed by time series of confirmed cases, while in phylodynamics genetic sequences of the pathogen, sampled through time, are the primary data source. Each type of data provides different, and potentially complementary, insight. Recent studies have recognised that combining data sources can improve estimates of the transmission rate and the number of infected individuals. However, inference methods are typically highly specialised and field-specific and are either computationally prohibitive or require intensive simulation, limiting their real-time utility. We present a novel birth-death phylogenetic model and derive a tractable analytic approximation of its likelihood, the computational complexity of which is linear in the size of the dataset. This approach combines epidemiological and phylodynamic data to produce estimates of key parameters of transmission dynamics and the unobserved prevalence. Using simulated data, we show (a) that the approximation agrees well with existing methods, (b) validate the claim of linear complexity and (c) explore robustness to model misspecification. This approximation facilitates inference on large datasets, which is increasingly important as large genomic sequence datasets become commonplace.

Applying information theoretic measures to neuronal activity data enables the quantification of neuronal encoding quality. However, when the sample size is limited, a naïve estimation of the information content typically contains a systematic overestimation (upward bias), which may lead to misinterpretation of coding characteristics. This bias is exacerbated in Ca2+ imaging because of the temporal sparsity of elevated Ca2+ signals. Here, we introduce methods to correct for the bias in the naïve estimation of information content from limited sample sizes and temporally sparse neuronal activity. We demonstrate the higher accuracy of our methods over previous ones, when applied to Ca2+ imaging data recorded from the mouse hippocampus and primary visual cortex, as well as to simulated data with matching tuning properties and firing statistics. Our bias-correction methods allowed an accurate estimation of the information place cells carry about the animal's position (spatial information) and uncovered the spatial resolution of hippocampal coding. Furthermore, using our methods, we found that cells with higher peak firing rates carry higher spatial information per spike and exposed differences between distinct hippocampal subfields in the long-term evolution of the spatial code. These results could be masked by the bias when applying the commonly used naïve calculation of information content. Thus, a bias-free estimation of information content can uncover otherwise overlooked properties of the neural code.

While aerobic glycolysis, or the Warburg effect, has for a long time been considered a hallmark of tumor metabolism, recent studies have revealed a far more complex picture. Tumor cells exhibit widespread metabolic heterogeneity, not only in their presentation of the Warburg effect but also in the nutrients and the metabolic pathways they are dependent on. Moreover, tumor cells can switch between different metabolic phenotypes in response to environmental cues and therapeutic interventions. A framework to analyze the observed metabolic heterogeneity and plasticity is, however, lacking. Using a mechanistic model that includes the key metabolic pathways active in tumor cells, we show that the inhibition of phosphofructokinase by excess ATP in the cytoplasm can drive a preference for aerobic glycolysis in fast-proliferating tumor cells. The differing rates of ATP utilization by tumor cells can therefore drive heterogeneity with respect to the presentation of the Warburg effect. Building upon this idea, we couple the metabolic phenotype of tumor cells to their migratory phenotype, and show that our model predictions are in agreement with previous experiments. Next, we report that the reliance of proliferating cells on different anaplerotic pathways depends on the relative availability of glucose and glutamine, and can further drive metabolic heterogeneity. Finally, using treatment of melanoma cells with a BRAF inhibitor as an example, we show that our model can be used to predict the metabolic and gene expression changes in cancer cells in response to drug treatment. By making predictions that are far more generalizable and interpretable as compared to previous tumor metabolism modeling approaches, our framework identifies key principles that govern tumor cell metabolism, and the reported heterogeneity and plasticity. These principles could be key to targeting the metabolic vulnerabilities of cancer.

Nucleosomes are recognized as key regulators of transcription. However, the relationship between slow nucleosome unwrapping dynamics and bulk transcriptional properties has not been thoroughly explored. Here, an agent-based model that we call the dynamic defect Totally Asymmetric Simple Exclusion Process (ddTASEP) was constructed to investigate the effects of nucleosome-induced pausing on transcriptional dynamics. Pausing due to slow nucleosome dynamics induced RNAPII convoy formation, which would cooperatively prevent nucleosome rebinding leading to bursts of transcription. The mean first passage time (MFPT) and the variance of first passage time (VFPT) were analytically expressed in terms of the nucleosome rate constants, allowing for the direct quantification of the effects of nucleosome-induced pausing on pioneering polymerase dynamics. The mean first passage elongation rate γ(hc, ho) is inversely proportional to the MFPT and can be considered to be a new axis of the ddTASEP phase diagram, orthogonal to the classical αβ-plane (where α and β are the initiation and termination rates). Subsequently, we showed that, for β = 1, there is a novel jamming transition in the αγ-plane that separates the ddTASEP dynamics into initiation-limited and nucleosome pausing-limited regions. We propose analytical estimates for the RNAPII density ρ, average elongation rate v, and transcription flux J and verified them numerically. We demonstrate that the intra-burst RNAPII waiting times tin follow the time-headway distribution of a max flux TASEP and that the average inter-burst interval [Formula: see text] correlates with the index of dispersion De. In the limit γ→0, the average burst size reaches a maximum set by the closing rate hc. When α≪1, the burst sizes are geometrically distributed, allowing large bursts even while the average burst size [Formula: see text] is small. Last, preliminary results on the relative effects of static and dynamic defects are presented to show that dynamic defects can induce equal or greater pausing than static bottle necks.

Is it possible to learn and create a first Hidden Markov Model (HMM) without programming skills or understanding the algorithms in detail? In this concise tutorial, we present the HMM through the 2 general questions it was initially developed to answer and describe its elements. The HMM elements include variables, hidden and observed parameters, the vector of initial probabilities, and the transition and emission probability matrices. Then, we suggest a set of ordered steps, for modeling the variables and illustrate them with a simple exercise of modeling and predicting transmembrane segments in a protein sequence. Finally, we show how to interpret the results of the algorithms for this particular problem. To guide the process of information input and explicit solution of the basic HMM algorithms that answer the HMM questions posed, we developed an educational webserver called HMMTeacher. Additional solved HMM modeling exercises can be found in the user's manual and answers to frequently asked questions. HMMTeacher is available at https://hmmteacher.mobilomics.org, mirrored at https://hmmteacher1.mobilomics.org. A repository with the code of the tool and the webpage is available at https://gitlab.com/kmilo.f/hmmteacher.

Brain networks store new memories using functional and structural synaptic plasticity. Memory formation is generally attributed to Hebbian plasticity, while homeostatic plasticity is thought to have an ancillary role in stabilizing network dynamics. Here we report that homeostatic plasticity alone can also lead to the formation of stable memories. We analyze this phenomenon using a new theory of network remodeling, combined with numerical simulations of recurrent spiking neural networks that exhibit structural plasticity based on firing rate homeostasis. These networks are able to store repeatedly presented patterns and recall them upon the presentation of incomplete cues. Storage is fast, governed by the homeostatic drift. In contrast, forgetting is slow, driven by a diffusion process. Joint stimulation of neurons induces the growth of associative connections between them, leading to the formation of memory engrams. These memories are stored in a distributed fashion throughout connectivity matrix, and individual synaptic connections have only a small influence. Although memory-specific connections are increased in number, the total number of inputs and outputs of neurons undergo only small changes during stimulation. We find that homeostatic structural plasticity induces a specific type of "silent memories", different from conventional attractor states.

Biocatalysis is a promising approach to sustainably synthesize pharmaceuticals, complex natural products, and commodity chemicals at scale. However, the adoption of biocatalysis is limited by our ability to select enzymes that will catalyze their natural chemical transformation on non-natural substrates. While machine learning and in silico directed evolution are well-posed for this predictive modeling challenge, efforts to date have primarily aimed to increase activity against a single known substrate, rather than to identify enzymes capable of acting on new substrates of interest. To address this need, we curate 6 different high-quality enzyme family screens from the literature that each measure multiple enzymes against multiple substrates. We compare machine learning-based compound-protein interaction (CPI) modeling approaches from the literature used for predicting drug-target interactions. Surprisingly, comparing these interaction-based models against collections of independent (single task) enzyme-only or substrate-only models reveals that current CPI approaches are incapable of learning interactions between compounds and proteins in the current family level data regime. We further validate this observation by demonstrating that our no-interaction baseline can outperform CPI-based models from the literature used to guide the discovery of kinase inhibitors. Given the high performance of non-interaction based models, we introduce a new structure-based strategy for pooling residue representations across a protein sequence. Altogether, this work motivates a principled path forward in order to build and evaluate meaningful predictive models for biocatalysis and other drug discovery applications.

The massive assessment of immune evasion due to viral mutations that increase COVID-19 susceptibility can be computationally facilitated. The adaptive cytotoxic T response is critical during primary infection and the generation of long-term protection. Here, potential HLA class I epitopes in the SARS-CoV-2 proteome were predicted for 2,915 human alleles of 71 families using the netMHCIpan EL algorithm. Allele families showed extreme epitopic differences, underscoring genetic variability of protective capacity between humans. Up to 1,222 epitopes were associated with any of the twelve supertypes, that is, allele clusters covering 90% population. Next, from all mutations identified in ~118,000 viral NCBI isolates, those causing significant epitope score reduction were considered epitope escape mutations. These mutations mainly involved non-conservative substitutions at the second and C-terminal position of the ligand core, or total ligand removal by large recurrent deletions. Escape mutations affected 47% of supertype epitopes, which in 21% of cases concerned isolates from two or more sub-continental areas. Some of these changes were coupled, but never surpassed 15% of evaded epitopes for the same supertype in the same isolate, except for B27. In contrast to most supertypes, eight allele families mostly contained alleles with few SARS-CoV-2 ligands. Isolates harboring cytotoxic escape mutations for these families co-existed geographically within sub-Saharan and Asian populations enriched in these alleles according to the Allele Frequency Net Database. Collectively, our findings indicate that escape mutation events have already occurred for half of HLA class I supertype epitopes. However, it is presently unlikely that, overall, it poses a threat to the global population. In contrast, single and double mutations for susceptible alleles may be associated with viral selective pressure and alarming local outbreaks. The integration of genomic, geographical and immunoinformatic information eases the surveillance of variants potentially affecting the global population, as well as minority subpopulations.

Antimicrobial resistance presents a significant health care crisis. The mutation F98Y in Staphylococcus aureus dihydrofolate reductase (SaDHFR) confers resistance to the clinically important antifolate trimethoprim (TMP). Propargyl-linked antifolates (PLAs), next generation DHFR inhibitors, are much more resilient than TMP against this F98Y variant, yet this F98Y substitution still reduces efficacy of these agents. Surprisingly, differences in the enantiomeric configuration at the stereogenic center of PLAs influence the isomeric state of the NADPH cofactor. To understand the molecular basis of F98Y-mediated resistance and how PLAs' inhibition drives NADPH isomeric states, we used protein design algorithms in the osprey protein design software suite to analyze a comprehensive suite of structural, biophysical, biochemical, and computational data. Here, we present a model showing how F98Y SaDHFR exploits a different anomeric configuration of NADPH to evade certain PLAs' inhibition, while other PLAs remain unaffected by this resistance mechanism.

Successful navigation requires the ability to compute one's location and heading from incoming multisensory information. Previous work has shown that this multisensory input comes in two forms: body-based idiothetic cues, from one's own rotations and translations, and visual allothetic cues, from the environment (usually visual landmarks). However, exactly how these two streams of information are integrated is unclear, with some models suggesting the body-based idiothetic and visual allothetic cues are combined, while others suggest they compete. In this paper we investigated the integration of body-based idiothetic and visual allothetic cues in the computation of heading using virtual reality. In our experiment, participants performed a series of body turns of up to 360 degrees in the dark with only a brief flash (300ms) of visual feedback en route. Because the environment was virtual, we had full control over the visual feedback and were able to vary the offset between this feedback and the true heading angle. By measuring the effect of the feedback offset on the angle participants turned, we were able to determine the extent to which they incorporated visual feedback as a function of the offset error. By further modeling this behavior we were able to quantify the computations people used. While there were considerable individual differences in performance on our task, with some participants mostly ignoring the visual feedback and others relying on it almost entirely, our modeling results suggest that almost all participants used the same strategy in which idiothetic and allothetic cues are combined when the mismatch between them is small, but compete when the mismatch is large. These findings suggest that participants update their estimate of heading using a hybrid strategy that mixes the combination and competition of cues.

Traditional (genome-scale) metabolic models of cellular growth involve an approximate biomass "reaction", which specifies biomass composition in terms of precursor metabolites (such as amino acids and nucleotides). On the one hand, biomass composition is often not known exactly and may vary drastically between conditions and strains. On the other hand, the predictions of computational models crucially depend on biomass. Also elementary flux modes (EFMs), which generate the flux cone, depend on the biomass reaction. To better understand cellular phenotypes across growth conditions, we introduce and analyze new classes of elementary vectors for comprehensive (next-generation) metabolic models, involving explicit synthesis reactions for all macromolecules. Elementary growth modes (EGMs) are given by stoichiometry and generate the growth cone. Unlike EFMs, they are not support-minimal, in general, but cannot be decomposed "without cancellations". In models with additional (capacity) constraints, elementary growth vectors (EGVs) generate a growth polyhedron and depend also on growth rate. However, EGMs/EGVs do not depend on the biomass composition. In fact, they cover all possible biomass compositions and can be seen as unbiased versions of elementary flux modes/vectors (EFMs/EFVs) used in traditional models. To relate the new concepts to other branches of theory, we consider autocatalytic sets of reactions. Further, we illustrate our results in a small model of a self-fabricating cell, involving glucose and ammonium uptake, amino acid and lipid synthesis, and the expression of all enzymes and the ribosome itself. In particular, we study the variation of biomass composition as a function of growth rate. In agreement with experimental data, low nitrogen uptake correlates with high carbon (lipid) storage.