Search

The collection of immunoglobulin genes in an individual’s germline, which gives rise to B cell receptors via recombination, is known to vary significantly across individuals. In humans, for example, each individual has only a fraction of the several hundred known V alleles. Furthermore, the currently-accepted set of known V alleles is both incomplete (particularly for non-European samples), and contains a significant number of spurious alleles. The resulting uncertainty as to which immunoglobulin alleles are present in any given sample results in inaccurate B cell receptor sequence annotations, and in particular inaccurate inferred naive ancestors. In this paper we first show that the currently widespread practice of aligning each sequence to its closest match in the full set of IMGT alleles results in a very large number of spurious alleles that are not in the sample’s true set of germline V alleles. We then describe a new method for inferring each individual’s germline gene set from deep sequencing data, and show that it improves upon existing methods by making a detailed comparison on a variety of simulated and real data samples. This new method has been integrated into the partis annotation and clonal family inference package, available at https://github.com/psathyrella/partis, and is run by default without affecting overall run time., Author summary Antibodies are an important component of the adaptive immune system, which itself determines our response to both pathogens and vaccines. They are produced by B cells through somatic recombination of germline DNA, which results in a vast diversity of antigen binding affinities across the B cell repertoire. We typically learn about the development of this repertoire, and its history of interaction with antigens, by sequencing large numbers of the DNA sequences from which antibodies are derived. In order to understand such data, it is necessary to determine the combination of germline V, D, and J genes that was rearranged to form each such B cell receptor sequence. This is difficult, however, because the immunoglobulin locus exhibits an extraordinary level of diversity across individuals—encompassing both allelic variation and gene duplication, deletion, and conversion—and because the locus’s large size and repetitive structure make germline sequencing very difficult. In this paper we describe a new computational method that avoids this difficulty by inferring each individual’s set of immunoglobulin germline genes directly from expressed B cell receptor sequence data.

Phylogenetic networks are rooted, directed, acyclic graphs that model reticulate evolutionary histories. Recently, statistical methods were devised for inferring such networks from either gene tree estimates or the sequence alignments of multiple unlinked loci. Bi-allelic markers, most notably single nucleotide polymorphisms (SNPs) and amplified fragment length polymorphisms (AFLPs), provide a powerful source of genome-wide data. In a recent paper, a method called SNAPP was introduced for statistical inference of species trees from unlinked bi-allelic markers. The generative process assumed by the method combined both a model of evolution for the bi-allelic markers, as well as the multispecies coalescent. A novel component of the method was a polynomial-time algorithm for exact computation of the likelihood of a fixed species tree via integration over all possible gene trees for a given marker. Here we report on a method for Bayesian inference of phylogenetic networks from bi-allelic markers. Our method significantly extends the algorithm for exact computation of phylogenetic network likelihood via integration over all possible gene trees. Unlike the case of species trees, the algorithm is no longer polynomial-time on all instances of phylogenetic networks. Furthermore, the method utilizes a reversible-jump MCMC technique to sample the posterior of phylogenetic networks given bi-allelic marker data. Our method has a very good performance in terms of accuracy and robustness as we demonstrate on simulated data, as well as a data set of multiple New Zealand species of the plant genus Ourisia (Plantaginaceae). We implemented the method in the publicly available, open-source PhyloNet software package., Author summary The availability of genomic data has revolutionized the study of evolutionary histories and phylogeny inference. Inferring evolutionary histories from genomic data requires, in most cases, accounting for the fact that different genomic regions could have evolutionary histories that differ from each other as well as from that of the species from which the genomes were sampled. In this paper, we introduce a method for inferring evolutionary histories while accounting for two processes that could give rise to such differences across the genomes, namely incomplete lineage sorting and hybridization. We introduce a novel algorithm for computing the likelihood of phylogenetic networks from bi-allelic genetic markers and use it in a Bayesian inference method. Analyses of synthetic and empirical data sets show a very good performance of the method in terms of the estimates it obtains.

Cancer cells share several metabolic traits, including aerobic production of lactate from glucose (Warburg effect), extensive glutamine utilization and impaired mitochondrial electron flow. It is still unclear how these metabolic rearrangements, which may involve different molecular events in different cells, contribute to a selective advantage for cancer cell proliferation. To ascertain which metabolic pathways are used to convert glucose and glutamine to balanced energy and biomass production, we performed systematic constraint-based simulations of a model of human central metabolism. Sampling of the feasible flux space allowed us to obtain a large number of randomly mutated cells simulated at different glutamine and glucose uptake rates. We observed that, in the limited subset of proliferating cells, most displayed fermentation of glucose to lactate in the presence of oxygen. At high utilization rates of glutamine, oxidative utilization of glucose was decreased, while the production of lactate from glutamine was enhanced. This emergent phenotype was observed only when the available carbon exceeded the amount that could be fully oxidized by the available oxygen. Under the latter conditions, standard Flux Balance Analysis indicated that: this metabolic pattern is optimal to maximize biomass and ATP production; it requires the activity of a branched TCA cycle, in which glutamine-dependent reductive carboxylation cooperates to the production of lipids and proteins; it is sustained by a variety of redox-controlled metabolic reactions. In a K-ras transformed cell line we experimentally assessed glutamine-induced metabolic changes. We validated computational results through an extension of Flux Balance Analysis that allows prediction of metabolite variations. Taken together these findings offer new understanding of the logic of the metabolic reprogramming that underlies cancer cell growth., Author summary Hallmarks describing common key events in initiation, maintenance and progression of cancer have been identified. One hallmark deals with rewiring of metabolic reactions required to sustain enhanced cell proliferation. The availability of molecular, mechanistic models of cancer hallmarks will mightily improve optimized personal treatment and new drug discovery. Metabolism is the only hallmark for which it is currently possible to derive large scale mathematical models, which have predictive ability. In this paper, we exploit a constraint-based model of the core metabolism required for biomass conversion of the most relevant nutrients—glucose and glutamine—to clarify the logic of control of cancer metabolism. We newly report that, when available oxygen is not sufficient to fully oxidize available glucose and glutamine carbons–a situation compatible with that observed under normal oxygen conditions in human and in cancer cells growing in vitro—utilization of glutamine by reductive carboxylation and conversion of glucose and glutamine to lactate confer advantage for biomass production. Redox homeostasis can be maintained through the use of different alternative pathways. In conclusion, this paper offers a logic interpretation to the link between metabolic rewiring and enhanced proliferation, which may offer new approaches to targeted drug discovery and utilization.

Metaproteomic studies adopt the common bottom-up proteomics approach to investigate the protein composition and the dynamics of protein expression in microbial communities. When matched metagenomic and/or metatranscriptomic data of the microbial communities are available, metaproteomic data analyses often employ a metagenome-guided approach, in which complete or fragmental protein-coding genes are first directly predicted from metagenomic (and/or metatranscriptomic) sequences or from their assemblies, and the resulting protein sequences are then used as the reference database for peptide/protein identification from MS/MS spectra. This approach is often limited because protein coding genes predicted from metagenomes are incomplete and fragmental. In this paper, we present a graph-centric approach to improving metagenome-guided peptide and protein identification in metaproteomics. Our method exploits the de Bruijn graph structure reported by metagenome assembly algorithms to generate a comprehensive database of protein sequences encoded in the community. We tested our method using several public metaproteomic datasets with matched metagenomic and metatranscriptomic sequencing data acquired from complex microbial communities in a biological wastewater treatment plant. The results showed that many more peptides and proteins can be identified when assembly graphs were utilized, improving the characterization of the proteins expressed in the microbial communities. The additional proteins we identified contribute to the characterization of important pathways such as those involved in degradation of chemical hazards. Our tools are released as open-source software on github at https://github.com/COL-IU/Graph2Pro., Author Summary In recent years, meta-omic (including metatranscriptomic and metaproteomic) techniques have been adopted as complementary approaches to metagenomic sequencing to study functional characteristics and dynamics of microbial communities, aiming at a holistic understanding of a community to respond to the changes in the environment. Currently, metaproteomic data are largely analyzed using the bioinformatics tools originally designed in bottom-up proteomics. In particular, recent metaproteomic studies employed a metagenome-guided approach, in which complete or fragmental protein-coding genes were first predicted from metagenomic sequences (i.e., contigs or scaffolds), acquired from the matched community samples, and predicted protein sequences were then used in peptide identification. A key challenge of this approach is that the protein coding genes predicted from assembled metagenomic contigs can be incomplete and fragmented due to the complexity of metagenomic samples and the short reads length in metagenomic sequencing. To address this issue, in this paper, we present a graph-centric approach that exploits the de bruijn graph structure reported by metagenome assembly algorithms to improve metagenome-guided peptide and protein identification in metaproteomics. We show that our method can identify much more peptides and proteins, improving the characterization of the proteins expressed in the microbial communities.

This paper presents a framework for modelling biological tissues based on discrete particles. Cell components (e.g. cell membranes, cell cytoskeleton, cell nucleus) and extracellular matrix (e.g. collagen) are represented using collections of particles. Simple particle to particle interaction laws are used to simulate and control complex physical interaction types (e.g. cell-cell adhesion via cadherins, integrin basement membrane attachment, cytoskeletal mechanical properties). Particles may be given the capacity to change their properties and behaviours in response to changes in the cellular microenvironment (e.g., in response to cell-cell signalling or mechanical loadings). Each particle is in effect an ‘agent’, meaning that the agent can sense local environmental information and respond according to pre-determined or stochastic events. The behaviour of the proposed framework is exemplified through several biological problems of ongoing interest. These examples illustrate how the modelling framework allows enormous flexibility for representing the mechanical behaviour of different tissues, and we argue this is a more intuitive approach than perhaps offered by traditional continuum methods. Because of this flexibility, we believe the discrete modelling framework provides an avenue for biologists and bioengineers to explore the behaviour of tissue systems in a computational laboratory., Author Summary Modelling is an important tool in understanding the behaviour of biological tissues. In this paper we advocate a new modelling framework in which cells and tissues are represented by a collection of particles with associated properties. The particles interact with each other and can change their behaviour in response to changes in their environment. We demonstrate how the propose framework can be used to represent the mechanical behaviour of different tissues with much greater flexibility as compared to traditional continuum based methods.

Correlation patterns in multiple sequence alignments of homologous proteins can be exploited to infer information on the three-dimensional structure of their members. The typical pipeline to address this task, which we in this paper refer to as the three dimensions of contact prediction, is to (i) filter and align the raw sequence data representing the evolutionarily related proteins; (ii) choose a predictive model to describe a sequence alignment; (iii) infer the model parameters and interpret them in terms of structural properties, such as an accurate contact map. We show here that all three dimensions are important for overall prediction success. In particular, we show that it is possible to improve significantly along the second dimension by going beyond the pair-wise Potts models from statistical physics, which have hitherto been the focus of the field. These (simple) extensions are motivated by multiple sequence alignments often containing long stretches of gaps which, as a data feature, would be rather untypical for independent samples drawn from a Potts model. Using a large test set of proteins we show that the combined improvements along the three dimensions are as large as any reported to date., Author Summary Proteins are large molecules that living cells make by stringing together building blocks called amino acids or peptides, following their blue-prints in the DNA. Freshly made proteins are typically long, structure-less chains of peptides, but shortly afterwards most of them fold into characteristic structures. Proteins execute many functions in the cell, for which they need to have the right structure, which is therefore very important in determining what the proteins can do. The structure of a protein can be determined by X-ray diffraction and other experimental approaches which are all, to this day, somewhat labor-intensive and difficult. On the other hand, the order of the peptides in a protein can be read off from the DNA blue-print, and such protein sequences are today routinely produced in large numbers. In this paper we show that many similar protein sequences can be used to find information about the structure. The basic approach is to construct a probabilistic model for sequence variability, and then to use the parameters of that model to predict structure in three-dimensional space. The main technical novelty compared to previous contributions in the same general direction is that we use models more directly matched to the data.

Ethanol is one of the most widely used recreational substances in the world and due to its ubiquitous use, ethanol abuse has been the cause of over 3.3 million deaths each year. In addition to its effects, ethanol’s primary metabolite, acetaldehyde, is a carcinogen that can cause symptoms of facial flushing, headaches, and nausea. How strongly ethanol or acetaldehyde affects an individual depends highly on the genetic polymorphisms of certain genes. In particular, the genetic polymorphisms of mitochondrial aldehyde dehydrogenase, ALDH2, play a large role in the metabolism of acetaldehyde. Thus, it is important to characterize how genetic variations can lead to different exposures and responses to ethanol and acetaldehyde. While the pharmacokinetics of ethanol metabolism through alcohol dehydrogenase have been thoroughly explored in previous studies, in this paper, we combined a base physiologically-based pharmacokinetic (PBPK) model with a whole-body genome-scale model (WBM) to gain further insight into the effect of other less explored processes and genetic variations on ethanol metabolism. This combined model was fit to clinical data and used to show the effect of alcohol concentrations, organ damage, ALDH2 enzyme polymorphisms, and ALDH2-inhibiting drug disulfiram on ethanol and acetaldehyde exposure. Through estimating the reaction rates of auxiliary processes with dynamic Flux Balance Analysis, The PBPK-WBM was able to navigate around a lack of kinetic constants traditionally associated with PK modelling and demonstrate the compensatory effects of the body in response to decreased liver enzyme expression. Additionally, the model demonstrated that acetaldehyde exposure increased with higher dosages of disulfiram and decreased ALDH2 efficiency, and that moderate consumption rates of ethanol could lead to unexpected accumulations in acetaldehyde. This modelling framework combines the comprehensive steady-state analyses from genome-scale models with the dynamics of traditional PK models to create a highly personalized form of PBPK modelling that can push the boundaries of precision medicine., Author summary Alcohol is a widely used recreational drug in many parts of the world and it is often abused or misused, leading to the deaths of millions of people each year from driving under the influence and overdose. Additionally, the body breaks down alcohol into acetaldehyde, a carcinogen that has its own effects ranging from headaches and nausea to liver damage. The effects of ethanol and acetaldehyde vary due to genetic variations that create different forms of the enzymes responsible for breaking them down. Due to these differences, it is important to characterize how these changes affect the metabolism of alcohol and acetaldehyde. To capture these differences, we have created a new model that integrates the traditional pharmacokinetic model with a whole-body genome-scale model that can characterize different genetic variations. In addition, traditional models often require experimentally measured data, yet with this new framework we avoid this tedious process by mathematically solving the genome-scale model with the dynamic Flux Balance Analysis technique, allowing for gap filling. Through this model, we show that the whole-body genome-scale model demonstrates flexibility and robustness that has not been seen before in pharmacokinetic models. Our model combines advantages from both pharmacokinetic and genome-scale modelling and can be personalized to characterize individual reactions to other drugs and further precision medicine.

Individual cells show variability in their signaling dynamics that often correlates with phenotypic responses, indicating that cell-to-cell variability is not merely noise but can have functional consequences. Based on this observation, we reasoned that cell-to-cell variability under the same treatment condition could be explained in part by a single signaling motif that maps different upstream signals into a corresponding set of downstream responses. If this assumption holds, then repeated measurements of upstream and downstream signaling dynamics in a population of cells could provide information about the underlying signaling motif for a given pathway, even when no prior knowledge of that motif exists. To test these two hypotheses, we developed a computer algorithm called MISC (Motif Inference from Single Cells) that infers the underlying signaling motif from paired time-series measurements from individual cells. When applied to measurements of transcription factor and reporter gene expression in the yeast stress response, MISC predicted signaling motifs that were consistent with previous mechanistic models of transcription. The ability to detect the underlying mechanism became less certain when a cell’s upstream signal was randomly paired with another cell’s downstream response, demonstrating how averaging time-series measurements across a population obscures information about the underlying signaling mechanism. In some cases, motif predictions improved as more cells were added to the analysis. These results provide evidence that mechanistic information about cellular signaling networks can be systematically extracted from the dynamical patterns of single cells., Author summary Cells use molecular signaling networks to translate dynamically changing stimuli into appropriate downstream responses. Specialized network structures, or motifs, allow cells to properly decode a variety of temporal input signals. In this paper, we present an algorithm that infers signaling motifs from multiple examples of an upstream signal paired with its downstream response in a population of single cells. We compare the predictive power of single-cell versus averaged time-series traces and the incremental benefit of adding more single-cell traces to the algorithm. We use this approach to understand how yeast respond to environmental stresses.

The partial success of tumor immunotherapy induced by checkpoint blockade, which is not antigen-specific, suggests that the immune system of some patients contain antigen receptors able to specifically identify tumor cells. Here we focused on T-cell receptor (TCR) repertoires associated with spontaneous breast cancer. We studied the alpha and beta chain CDR3 domains of TCR repertoires of CD4 T cells using deep sequencing of cell populations in mice and applied the results to published TCR sequence data obtained from human patients. We screened peripheral blood T cells obtained monthly from individual mice spontaneously developing breast tumors by 5 months. We then looked at identical TCR sequences in published human studies; we used TCGA data from tumors and healthy tissues of 1,256 breast cancer resections and from 4 focused studies including sequences from tumors, lymph nodes, blood and healthy tissues, and from single cell dataset of 3 breast cancer subjects. We now report that mice spontaneously developing breast cancer manifest shared, Public CDR3 regions in both their alpha and beta and that a significant number of women with early breast cancer manifest identical CDR3 sequences. These findings suggest that the development of breast cancer is associated, across species, with biomarker, exclusive TCR repertoires., Author summary More than 1,100 cancer treatments are currently in clinical testing. Of these, 295 are immune-oncology medicines and vaccines. For Checkpoint Immunotherapy alone, 52,539 patients were enrolled in 2017. Of these treatments, 108 are for breast cancer, the leading cancer diagnosed in women. Although checkpoint immunotherapy has yet to be proven beneficial in most types of breast cancer, current observations of which specific T cell go into the tumor, and how these cells are connected to the response to treatment, have lead to a new understating of the role of the immune system. Yet, and in spite of this intense research into the role of T cells in tumors, very little is known about how clones of T cells behave during tumor formation. The paper we present here studies this clonal behavior over time. We study it during mammary tumor development in mice. We see that a set of T cell clones, which comes from a diversity of different nucleotide sequences, stands out in its behavior compared with other T cell clones. Importantly, we show that these mouse T cell clones have a strong connection to human breast cancer. The specific identification of a subset of T cell clones, together with the behavior of these clones over time, is critical to our understanding of the effect of modifications to the immune system, such as the modifications made by immunotherapy treatments, and to our understanding of how to better control the behavior of the immune system.

Fitting stochastic transmission models to electronic patient data can offer detailed insights into the transmission of healthcare-associated infections and improve infection control. Pathogen whole-genome sequencing may improve the precision of model inferences, but computational constraints have limited modelling applications predominantly to small datasets and specific outbreaks, whereas large-scale sequencing studies have mostly relied on simple rules for identifying/excluding plausible transmission. We present a novel approach for integrating detailed epidemiological data on patient contact networks in hospitals with large-scale pathogen sequencing data. We apply our approach to study Clostridioides difficile transmission using a dataset of 1223 infections in Oxfordshire, UK, 2007–2011. 262 (21% [95% credibility interval 20–22%]) infections were estimated to have been acquired from another known case. There was heterogeneity by sequence type (ST) in the proportion of cases acquired from another case with the highest rates in ST1 (ribotype-027), ST42 (ribotype-106) and ST3 (ribotype-001). These same STs also had higher rates of transmission mediated via environmental contamination/spores persisting after patient discharge/recovery; for ST1 these persisted longer than for most other STs except ST3 and ST42. We also identified variation in transmission between hospitals, medical specialties and over time; by 2011 nearly all transmission from known cases had ceased in our hospitals. Our findings support previous work suggesting only a minority of C. difficile infections are acquired from known cases but highlight a greater role for environmental contamination than previously thought. Our approach is applicable to other healthcare-associated infections. Our findings have important implications for effective control of C. difficile., Author summary Preventing infections spreading in hospitals is a major priority for healthcare systems globally. Mathematical models can be used with electronic hospital records to reconstruct how infections spread, which in turn can help guide infection control interventions. Sequencing the genetic code (DNA) of the bacteria that cause infections can help to follow when transmission is occurring, as bacterial DNA from two patients infected from the same source will likely be very similar. Our paper describes a new statistical method for combining hospital records and sequencing data to track infections. We use our approach to study over 1200 infections of Clostridioides difficile (C. diff.), which is a common cause of diarrhoea in hospitals. We show that only a minority of infections are acquired from other unwell patients, but the amount of spread varies by the subtype of C. diff involved. We also show that different C. diff subtypes survive in the hospital environment for longer than others and may need enhanced control strategies. We also can detect differences in spread at different hospitals and show that by the end of the study we had largely eliminated transmission of C. diff from unwell patients in our hospitals.

GC skew is a phenomenon observed in many bacterial genomes, wherein the two replication strands of the same chromosome contain different proportions of guanine and cytosine nucleotides. Here we demonstrate that this phenomenon, which was first discovered in the mid-1990s, can be used today as an analysis tool for the 15,000+ complete bacterial genomes in NCBI’s Refseq library. In order to analyze all 15,000+ genomes, we introduce a new method, SkewIT (Skew Index Test), that calculates a single metric representing the degree of GC skew for a genome. Using this metric, we demonstrate how GC skew patterns are conserved within certain bacterial phyla, e.g. Firmicutes, but show different patterns in other phylogenetic groups such as Actinobacteria. We also discovered that outlier values of SkewIT highlight potential bacterial mis-assemblies. Using our newly defined metric, we identify multiple mis-assembled chromosomal sequences in previously published complete bacterial genomes. We provide a SkewIT web app https://jenniferlu717.shinyapps.io/SkewIT/ that calculates SkewI for any user-provided bacterial sequence. The web app also provides an interactive interface for the data generated in this paper, allowing users to further investigate the SkewI values and thresholds of the Refseq-97 complete bacterial genomes. Individual scripts for analysis of bacterial genomes are provided in the following repository: https://github.com/jenniferlu717/SkewIT., Author summary Even though every guanine (G) is paired with a cytosine (C) in double-stranded DNA molecules, bacterial genomes have more G’s than C’s when we focus only on a single strand in the direction of replication, called the leading strand. This phenomenon, called GC skew, is so ubiquitous that it has been used reliably to identify the replication origin (the location from which DNA begins the process of replicating itself) in thousands of bacteria. Here we describe a new method that automatically captures the “skewness” of a genome by finding the origin and terminus of replication, and then reporting skewness as a single number. We calculated this value for over 15,000 genomes, and found that most phylogenetic groups have a characteristic amount of skewness. We also observed that an unusually low value for skewness sometimes indicated that the genome was incorrectly assembled. To assist others in this type of analysis, we developed a graphical tool to compute and display GC-skew for any genome of interest.

Millions of people worldwide develop foodborne illnesses caused by Salmonella enterica (S. enterica) every year. The pathogenesis of S. enterica depends on flagella, which are appendages that the bacteria use to move through the environment. Interestingly, populations of genetically identical bacteria exhibit heterogeneity in the number of flagella. To understand this heterogeneity and the regulation of flagella quantity, we propose a mathematical model that connects the flagellar gene regulatory network to flagellar construction. A regulatory network involving more than 60 genes controls flagellar assembly. The most important member of the network is the master operon, flhDC, which encodes the FlhD4C2 protein. FlhD4C2 controls the construction of flagella by initiating the production of hook basal bodies (HBBs), protein structures that anchor the flagella to the bacterium. By connecting a model of FlhD4C2 regulation to a model of HBB construction, we investigate the roles of various feedback mechanisms. Analysis of our model suggests that a combination of regulatory mechanisms at the protein and transcriptional levels induce bistable FlhD4C2 levels and heterogeneous numbers of flagella. Also, the balance of regulatory mechanisms that become active following HBB construction is sufficient to provide a counting mechanism for controlling the total number of flagella produced., Author summary Salmonella causes foodborne illnesses in millions of people worldwide each year. Flagella, which are appendages that the bacteria use to move through the environment, are a key factor in the infection process. Populations of genetically identical bacteria have been observed to contain both motile cells, generally with 6–10 flagella, and nonmotile cells, with no flagella. In this paper, we use mathematical models of the gene network that regulates flagellar construction to explore how the bacteria controls the number of flagella produced. We suggest that a bacterium must accumulate a threshold amount of a master regulator protein to initiate flagella production and failure to reach the threshold results in no flagella. Downstream mechanisms that impact the amount of master regulator protein are sufficient to determine how many flagella are produced.

Contributing to broader impacts is an important aspect of scientific research. Engaging practicing K-12 teachers as part of a research project can be an effective approach for addressing broader impacts requirements of grants, while also advancing researcher and teacher professional growth. Our focus is on leveraging teachers' professional expertise to develop science education materials grounded in emerging scientific research. In this paper, we describe ten simple rules for planning, implementing, and evaluating teacher engagement to support the broader impact goals of your research project. These collaborations can lead to the development of instructional materials or activities for students in the classroom or provide science research opportunities for teachers. We share our successes and lessons learned while collaborating with high school biology teachers to create technology-based, instructional materials developed from basic biological research. The rules we describe are applicable across teacher partnerships at any grade level in that they emphasize eliciting and respecting teachers' professionalism and expertise.

Understanding evolution requires detailed knowledge of genotype-phenotype maps; however, it can be a herculean task to measure every phenotype in a combinatorial map. We have developed a computational strategy to predict the missing phenotypes from an incomplete, combinatorial genotype-phenotype map. As a test case, we used an incomplete genotype-phenotype dataset previously generated for the malaria parasite’s ‘chloroquine resistance transporter’ (PfCRT). Wild-type PfCRT (PfCRT3D7) lacks significant chloroquine (CQ) transport activity, but the introduction of the eight mutations present in the ‘Dd2’ isoform of PfCRT (PfCRTDd2) enables the protein to transport CQ away from its site of antimalarial action. This gain of a transport function imparts CQ resistance to the parasite. A combinatorial map between PfCRT3D7 and PfCRTDd2 consists of 256 genotypes, of which only 52 have had their CQ transport activities measured through expression in the Xenopus laevis oocyte. We trained a statistical model with these 52 measurements to infer the CQ transport activity for the remaining 204 combinatorial genotypes between PfCRT3D7 and PfCRTDd2. Our best-performing model incorporated a binary classifier, a nonlinear scale, and additive effects for each mutation. The addition of specific pairwise- and high-order-epistatic coefficients decreased the predictive power of the model. We evaluated our predictions by experimentally measuring the CQ transport activities of 24 additional PfCRT genotypes. The R2 value between our predicted and newly-measured phenotypes was 0.90. We then used the model to probe the accessibility of evolutionary trajectories through the map. Approximately 1% of the possible trajectories between PfCRT3D7 and PfCRTDd2 are accessible; however, none of the trajectories entailed eight successive increases in CQ transport activity. These results demonstrate that phenotypes can be inferred with known uncertainty from a partial genotype-phenotype dataset. We also validated our approach against a collection of previously published genotype-phenotype maps. The model therefore appears general and should be applicable to a large number of genotype-phenotype maps., Author summary Biological macromolecules are built from chains of building blocks. The function of a macromolecule depends on the specific chemical properties of the building blocks that make it up. Macromolecules evolve through mutations that swap one building block for another. Understanding how biomolecules work and evolve therefore requires knowledge of the effects of mutations. The effects of mutations can be measured experimentally; however, because there are a vast number of possible combinations of mutations, it is often difficult to make enough measurements to understand biomolecular function and evolution. In this paper, we describe a simple method to predict the effects of mutations on biomolecules from a small number of measurements. This method works by appropriately averaging the effects of mutations seen in different contexts. We test the method by predicting the effects of mutations on a PfCRT—a macromolecule from the malarial parasite that confers drug resistance. We find that our method is fast and effective. Using a small number of measurements, we were able to gain insight into the evolutionary steps by which this macromolecule conferred drug resistance. To make this method accessible to other researchers, we have released it as an open-source software package: https://gpseer.readthedocs.io.

Over the past two decades, substantial amount of work has been conducted to characterize different odorant receptors, neuroanatomy and odorant response properties of the early olfactory system of Drosophila melanogaster. Yet many odorant receptors remain only partially characterized, and the odorant transduction process and the axon hillock spiking mechanism of the olfactory sensory neurons (OSNs) have yet to be fully determined. Identity and concentration, two key characteristics of the space of odorants, are encoded by the odorant transduction process. Detailed molecular models of the odorant transduction process are, however, scarce for fruit flies. To address these challenges we advance a comprehensive model of fruit fly OSNs as a cascade consisting of an odorant transduction process (OTP) and a biophysical spike generator (BSG). We model odorant identity and concentration using an odorant-receptor binding rate tensor, modulated by the odorant concentration profile, and an odorant-receptor dissociation rate tensor, and quantitatively describe the mechanics of the molecular ligand binding/dissociation of the OTP. We model the BSG as a Connor-Stevens point neuron. The resulting spatio-temporal encoding model of the Drosophila antenna provides a theoretical foundation for understanding the neural code of both odorant identity and odorant concentration and advances the state-of-the-art in a number of ways. First, it quantifies on the molecular level the spatio-temporal level of complexity of the transformation taking place in the antennae. The concentration-dependent spatio-temporal code at the output of the antenna circuits determines the level of complexity of olfactory processing in the downstream neuropils, such as odorant recognition and olfactory associative learning. Second, the model is biologically validated using multiple electrophysiological recordings. Third, the model demonstrates that the currently available data for odorant-receptor responses only enable the estimation of the affinity of the odorant-receptor pairs. The odorant-dissociation rate is only available for a few odorant-receptor pairs. Finally, our model calls for new experiments for massively identifying the odorant-receptor dissociation rates of relevance to flies., Author summary Identity and concentration, intrinsically embedded in the odorant space, are two key characteristics of olfactory coding that define the level of complexity of neural processing throughout the olfactory system in the fruit fly. In this paper we advance a theoretical foundation for understanding these two characteristics by quantifying mathematically the odorant space and devising a biophysical model of the olfactory sensory neurons (OSNs). To validate our modeling approach, we propose and apply an algorithm to estimate the affinity value and the dissociation rate, the two characteristics that define odorant identity, of multiple odorant-receptor pairs. We then evaluate our model with a multitude of odorant waveforms and demonstrate that the model output reproduces the temporal responses of OSNs obtained from in vivo electrophysiology recordings. Furthermore, we evaluate the model at the OSN population level and quantify on the molecular level the spatio-temporal level of complexity of the transformation taking place between the odorant space and the OSNs. The resulting concentration-dependent spatio-temporal code determines the level of complexity of the input space driving olfactory processing in the downstream neuropils. Lastly, our model demonstrates that the currently available data for OSN responses only enables estimation of affinity value. This calls for new experiments for massively identifying the odorant-receptor dissociation rates of relevance to flies.

Antigen recognition by T-cells is guided by the T-cell receptor (TCR) heterodimer formed by α and β chains. A huge diversity of TCR sequences should be maintained by the immune system in order to be able to mount an effective response towards foreign pathogens, so, due to cooperative binding of α and β chains to the pathogen, any constraints on chain pairing can have a profound effect on immune repertoire structure, diversity and antigen specificity. By integrating available structural data and paired chain sequencing results we were able to show that there are almost no constraints on pairing in TCRαβ complexes, allowing naive T-cell repertoire to reach the highest possible diversity. Additional analysis reveals that the specific choice of contacting amino acids can still have a profound effect on complex conformation. Moreover, antigen-driven selection can distort the uniform landscape of chain pairing, while small, yet significant, differences in the pairing can be attributed to various specialized T-cell subsets such as MAIT and iNKT T-cells, as well as other TCR sets specific to certain antigens., Author summary In the present paper we study chain pairing preferences in the T-cell receptor (TCR) heterodimer complex. The TCR molecule is formed by α and β chains and binding of both of these chains to an antigen presented by the major histocompatibility complex (MHC) molecule is required in order to trigger an immune response against foreign pathogens and neoantigens. We show that chain pairing in the TCR complex is nearly random ensuring a highly diverse set of TCRs required for recognition of a vast set of antigens. Our results also show that chain pairing preferences can nevertheless influence TCR complex geometry and biases in TCR chain pairing can be used to identify antigen-driven selection or selection towards specialized subsets of T-cells such as mucosal-associated and natural killer invariant T-cells.

Social insect colonies exhibit colony-level phenotypes such as social immunity and task coordination, which are produced by the individual phenotypes. Mapping the genetic basis of such phenotypes requires associating the colony-level phenotype with the genotypes in the colony. In this paper, we examine alternative approaches to DNA extraction, library construction, and sequencing for genome wide association studies (GWAS) of colony-level traits using a population sample of Cataglyphis niger ants. We evaluate the accuracy of allele frequency estimation from sequencing a pool of individuals (pool-seq) from each colony using either whole-genome sequencing or reduced representation genomic sequencing. Based on empirical measurement of the experimental noise in sequenced DNA pools, we show that reduced representation pool-seq is drastically less accurate than whole-genome pool-seq. Surprisingly, normalized pooling of samples did not result in greater accuracy than un-normalized pooling. Subsequently, we evaluate the power of the alternative approaches for detecting quantitative trait loci (QTL) of colony-level traits by using simulations that account for an environmental effect on the phenotype. Our results can inform experimental designs and enable optimizing the power of GWAS depending on budget, availability of samples and research goals. We conclude that for a given budget, sequencing un-normalized pools of individuals from each colony provides optimal QTL detection power., Author summary Genomic mapping techniques are used to map phenotypes to genotypes. Mapping is of general interest in any biological system, including fundamental studies of biological traits, clinical studies of genetic predisposition to disease, and agro- and bio-technological studies of domesticated plants and animals. Typically, such studies associate phenotypic measurements of individuals with their genotypes. Here we evaluate methodological approaches for genomic mapping of phenotypes that are expressed at the level of a group rather than that of individuals. We demonstrate that genomic sequencing of a DNA pool from multiple samples provides increased statistical power within a limited budget. Our results facilitate more efficient use of resources in genomic mapping studies that investigate group-level phenotypes.

Cell-to-cell variability generates subpopulations of drug-tolerant cells that diminish the efficacy of cancer drugs. Efficacious combination therapies are thus needed to block drug-tolerant cells via minimizing the impact of heterogeneity. Probabilistic models such as Bliss independence have been developed to evaluate drug interactions and their combination efficacy based on probabilities of specific actions mediated by drugs individually and in combination. In practice, however, these models are often applied to conventional dose-response curves in which a normalized parameter with a value between zero and one, generally referred to as fraction of cells affected (fa), is used to evaluate the efficacy of drugs and their combined interactions. We use basic probability theory, computer simulations, time-lapse live cell microscopy, and single-cell analysis to show that fa metrics may bias our assessment of drug efficacy and combination effectiveness. This bias may be corrected when dynamic probabilities of drug-induced phenotypic events, i.e. induction of cell death and inhibition of division, at a single-cell level are used as metrics to assess drug efficacy. Probabilistic phenotype metrics offer the following three benefits. First, in contrast to the commonly used fa metrics, they directly represent probabilities of drug action in a cell population. Therefore, they deconvolve differential degrees of drug effect on tumor cell killing versus inhibition of cell division, which may not be correlated for many drugs. Second, they increase the sensitivity of short-term drug response assays to cell-to-cell heterogeneities and the presence of drug-tolerant subpopulations. Third, their probabilistic nature allows them to be used directly in unbiased evaluation of synergistic efficacy in drug combinations using probabilistic models such as Bliss independence. Altogether, we envision that probabilistic analysis of single-cell phenotypes complements currently available assays via improving our understanding of heterogeneity in drug response, thereby facilitating the discovery of more efficacious combination therapies to block drug-tolerant cells., Author summary Resistance to therapy due to tumor cell heterogeneity poses a major challenge to the use of cancer drugs. Cell-to-cell variability generates subpopulations of drug-tolerant cells that diminish therapeutic efficacy, even in populations of cells that are scored as highly sensitive based on drug potency. Overcoming such heterogeneity and blocking subpopulations of drug-tolerant cells motivate efforts toward identifying efficacious combination therapies. The success of these efforts depends on our ability to distinguish how heterogeneous populations of cells respond to individual drugs, and how these responses are influenced by combined drug interactions. In this paper, we propose mathematical and experimental frameworks to evaluate time-dependent drug interactions based on probabilistic metrics that quantify drug-induced tumor cell killing or inhibition of division at a single-cell level. These metrics can reveal heterogeneous drug responses and their changes with time and drug combinations. Thus, they have important implications for designing efficacious combination therapies, especially those designed to block or overcome drug-tolerant subpopulations of cancer cells.

Identification of induced pluripotent stem (iPS) progenitor cells, the iPS forming cells in early stage of reprogramming, could provide valuable information for studying the origin and underlying mechanism of iPS cells. However, it is very difficult to identify experimentally since there are no biomarkers known for early progenitor cells, and only about 6 days after reprogramming initiation, iPS cells can be experimentally determined via fluorescent probes. What is more, the ratio of progenitor cells during early reprograming period is below 5%, which is too low to capture experimentally in the early stage. In this paper, we propose a novel computational approach for the identification of iPS progenitor cells based on machine learning and microscopic image analysis. Firstly, we record the reprogramming process using a live cell imaging system after 48 hours of infection with retroviruses expressing Oct4, Sox2 and Klf4, later iPS progenitor cells and normal murine embryonic fibroblasts (MEFs) within 3 to 5 days after infection are labeled by retrospectively tracing the time-lapse microscopic image. We then calculate 11 types of cell morphological and motion features such as area, speed, etc., and select best time windows for modeling and perform feature selection. Finally, a prediction model using XGBoost is built based on the selected six types of features and best time windows. Our model allows several missing values/frames in the sample datasets, thus it is applicable to a wide range of scenarios. Cross-validation, holdout validation and independent test experiments show that the minimum precision is above 52%, that is, the ratio of predicted progenitor cells within 3 to 5 days after viral infection is above 52%. The results also confirm that the morphology and motion pattern of iPS progenitor cells is different from that of normal MEFs, which helps with the machine learning methods for iPS progenitor cell identification., Author summary Identification of induced pluripotent stem (iPS) progenitor cells could provide valuable information for studying the origin and underlying mechanism of iPS cells. However, it is very difficult to identify experimentally since there are no biomarkers known for early progenitor cells, and only after about 6 days of induction, iPS cells can be experimentally determined via fluorescent probes. What is more, the percentage of the progenitor cells during the early induction period is below 5%, too low to capture experimentally in early stage. In this work, we proposed an approach for the identification of iPS progenitor cells, the iPS forming cells, based on machine learning and microscopic image analysis. The aim is to help biologists to enrich iPS progenitor cells during the early stage of induction, which allows experimentalists to select iPS progenitor cells with much higher probability, and furthermore to study the biomarkers which trigger the reprogramming process.

Pluripotent stem cells retain the developmental timing of their species of origin in vitro, an observation that suggests the existence of a cell-intrinsic developmental clock, yet the nature and machinery of the clock remain a mystery. We hypothesize that one possible component may lie in species-specific differences in the kinetics of transcriptional responses to differentiation signals. Using a liquid-handling robot, mouse and human pluripotent stem cells were exposed to identical neural differentiation conditions and sampled for RNA-sequencing at high frequency, every 4 or 10 minutes, for the first 10 hours of differentiation to test for differences in transcriptomic response rates. The majority of initial transcriptional responses occurred within a rapid window in the first minutes of differentiation for both human and mouse stem cells. Despite similarly early onsets of gene expression changes, we observed shortened and condensed gene expression patterns in mouse pluripotent stem cells compared to protracted trends in human pluripotent stem cells. Moreover, the speed at which individual genes were upregulated, as measured by the slopes of gene expression changes over time, was significantly faster in mouse compared to human cells. These results suggest that downstream transcriptomic response kinetics to signaling cues are faster in mouse versus human cells, and may offer a partial account for the vast differences in developmental rates across species., Author summary For regenerative medicine, a slow developmental clock operating in human pluripotent stem cells imposes severe obstacles to many pluripotent stem cell therapies. Producing the desired cell types with the necessary level of physiological maturity in clinically-relevant time frames will be essential if the therapy is aimed at being patient-specific. Thus, understanding what controls species-specific rates of differentiation and identifying ways to accelerate it would have a broad impact on improving the therapeutic value of human pluripotent stem cells. In this paper, we measured the immediate-early gene expression responses to neural differentiation signals in mouse and human pluripotent stem cells. We quantitatively summarized the gene expression patterns using a flexible model that allowed us to identify shared trends and dynamic changes across species. We found that responses occurred equally early across species, but expression patterns were prolonged in human cells and the rate of upregulation in gene expression was significantly faster in mouse cells. Our results represent a significant contribution in furthering the knowledge of early developmental transcriptional programs.

Phenotypic plasticity is a vital strategy for plants to deal with changing conditions by inducing phenotypes favourable in different environments. Understanding how natural selection acts on variation in phenotypic plasticity in plants is therefore a central question in ecology, but is often ignored in modelling studies. Here we present a new modelling approach that allows for the analysis of selection for variation in phenotypic plasticity as a response strategy. We assess selection for shade avoidance strategies of Arabidopsis thaliana in response to future neighbour shading signalled through a decrease in red:far-red (R:FR) ratio. For this, we used a spatially explicit 3D virtual plant model that simulates individual Arabidopsis plants competing for light in different planting densities. Plant structure and growth were determined by the organ-specific interactions with the light environment created by the vegetation structure itself. Shade avoidance plastic responses were defined by a plastic response curve relating petiole elongation and lamina growth to R:FR perceived locally. Different plasticity strategies were represented by different shapes of the response curve that expressed different levels of R:FR sensitivity. Our analyses show that the shape of the selected shade avoidance strategy varies with planting density. At higher planting densities, more sensitive response curves are selected for than at lower densities. In addition, the balance between lamina and petiole responses influences the sensitivity of the response curves selected for. Combining computational virtual plant modelling with a game theoretical analysis represents a new step towards analysing how natural selection could have acted upon variation in shade avoidance as a response strategy, which can be linked to genetic variation and underlying physiological processes., Author summary Plants are able to respond to changes in the environment. Particularly, plants can show different structural traits (e.g. stem height and leaf size) in different planting densities. These trait changes are the result of so-called plastic responses that can be induced by changes in the light spectrum. Although a great part of the physiological processes underlying these plastic responses have been identified, it remains unclear how these plastic responses, and variations therein, can be the result of selection. In this paper we analyse selection on different plastic responses within different dynamic competitive environments. We use a 3D virtual plant model that simulates realistic plant growth based on light absorption, photosynthesis and specific light signals that induce changes in leaf growth. The 3D model simulated various competitive vegetation stands consisting of plants with different plastic responses at different planting densities. We conclude that selection in different densities results in different plastic responses. We advocate that our modelling approach allows for analyses related to selection on plastic responses itself, instead of on specific trait values in different environments, and is therefore an important new step forward in understanding the role of plastic responses for plant performance.

The origin of avian flight is one of the most controversial debates in Paleontology. This paper investigates the wing performance of Caudipteryx, the most basal non-volant dinosaur with pennaceous feathered forelimbs by using modal effective mass theory. From a mechanical standpoint, the forced vibrations excited by hindlimb locomotion stimulate the movement of wings, creating a flapping-like motion in response. This shows that the origin of the avian flight stroke should lie in a completely natural process of active locomotion on the ground. In this regard, flapping in the history of evolution of avian flight should have already occurred when the dinosaurs were equipped with pennaceous remiges and rectrices. The forced vibrations provided the initial training for flapping the feathered wings of theropods similar to Caudipteryx., Author summary The origin of avian flight in the perspective of mechanics has been investigated for the first time. We reported the first evidence for flapping hypothesis based on principle of physical modeling. This is significant because using modal effective mass method and reconstructed Caudipteryx, the most basal non-volant winged dinosaur, we captured significant and negligible modes and realized that resonance oscillation of Caudipteryx wings could occur as the running speed approached to the primary frequencies. Such forced vibrations induced by legs' motions during running trained the Caudipteryx and the other feathered dinosaurs to flap their wings.

Basal gene expression levels have been shown to be predictive of cellular response to cytotoxic treatments. However, such analyses do not fully reveal complex genotype- phenotype relationships, which are partly encoded in highly interconnected molecular networks. Biological pathways provide a complementary way of understanding drug response variation among individuals. In this study, we integrate chemosensitivity data from a large-scale pharmacogenomics study with basal gene expression data from the CCLE project and prior knowledge of molecular networks to identify specific pathways mediating chemical response. We first develop a computational method called PACER, which ranks pathways for enrichment in a given set of genes using a novel network embedding method. It examines a molecular network that encodes known gene-gene as well as gene-pathway relationships, and determines a vector representation of each gene and pathway in the same low-dimensional vector space. The relevance of a pathway to the given gene set is then captured by the similarity between the pathway vector and gene vectors. To apply this approach to chemosensitivity data, we identify genes whose basal expression levels in a panel of cell lines are correlated with cytotoxic response to a compound, and then rank pathways for relevance to these response-correlated genes using PACER. Extensive evaluation of this approach on benchmarks constructed from databases of compound target genes and large collections of drug response signatures demonstrates its advantages in identifying compound-pathway associations compared to existing statistical methods of pathway enrichment analysis. The associations identified by PACER can serve as testable hypotheses on chemosensitivity pathways and help further study the mechanisms of action of specific cytotoxic drugs. More broadly, PACER represents a novel technique of identifying enriched properties of any gene set of interest while also taking into account networks of known gene-gene relationships and interactions., Author summary Gene expression levels have been used to study the cellular response to drug treatments. However, analysis of gene expression without considering gene interactions cannot fully reveal complex genotype-phenotype relationships. Biological pathways reveal the interactions among genes, thus providing a complementary way of understanding the drug response variation among individuals. In this paper, we aim to identify pathways that mediate the chemical response of each drug. We used the recently generated CTRP pharmacogenomics data and CCLE basal expression data to identify these pathways. We showed that using the prior knowledge encoded in molecular networks substantially improves pathway identification. In particular, we integrate genes and pathways into a large heterogeneous network in which links are protein-protein interactions and gene-pathway affiliations. We then project this heterogeneous network onto a low-dimensional space, which enables more precise similarity measurements between pathways and drug-response-correlated genes. Extensive experiments on two benchmarks show that our method substantially improved the pathway identification performance by using the molecular networks. More importantly, our method represents a novel technique of identifying enriched properties of any gene set of interest while also taking into account networks of known gene-gene relationships and interactions.

Mutation signatures in cancer genomes reflect endogenous and exogenous mutational processes, offering insights into tumour etiology, features for prognostic and biologic stratification and vulnerabilities to be exploited therapeutically. We present a novel machine learning formalism for improved signature inference, based on multi-modal correlated topic models (MMCTM) which can at once infer signatures from both single nucleotide and structural variation counts derived from cancer genome sequencing data. We exemplify the utility of our approach on two hormone driven, DNA repair deficient cancers: breast and ovary (n = 755 samples total). We show how introducing correlated structure both within and between modes of mutation can increase accuracy of signature discovery, particularly in the context of sparse data. Our study emphasizes the importance of integrating multiple mutation modes for signature discovery and patient stratification, and provides a statistical modeling framework to incorporate additional features of interest for future studies., Author summary Over time DNA accumulates mutations from a variety of sources. Some mutations result from external mutagens, such as UV radiation, while others result from processes occurring within the cell itself. Each of these sources can impart characteristic patterns of mutations on the genome, known as mutation signatures, which can be detected using computational techniques. Loss of DNA repair mechanisms can leave specific mutation signatures in the genomes of cancer cells. To identify cancers with broken DNA-repair processes, accurate methods are needed for detecting mutation signatures and, in particular, their activities or probabilities within individual cancers. In this paper, we introduce a class of statistical modeling methods used for natural language processing, known as “topic models”, that outperform standard methods for signature analysis. We show that topic models that incorporate signature probability correlations across cancers perform best, while jointly analyzing multiple mutation types improves robustness to low mutation counts.

Filamentous actin (F-actin) and non-muscle myosin II motors drive cell motility and cell shape changes that guide large scale tissue movements during embryonic morphogenesis. To gain a better understanding of the role of actomyosin in vivo, we have developed a two-dimensional (2D) computational model to study emergent phenomena of dynamic unbranched actomyosin arrays in the cell cortex. These phenomena include actomyosin punctuated contractions, or "actin asters" that form within quiescent F-actin networks. Punctuated contractions involve both formation of high intensity aster-like structures and disassembly of those same structures. Our 2D model allows us to explore the kinematics of filament polarity sorting, segregation of motors, and morphology of F-actin arrays that emerge as the model structure and biophysical properties are varied. Our model demonstrates the complex, emergent feedback between filament reorganization and motor transport that generate as well as disassemble actin asters. Since intracellular actomyosin dynamics are thought to be controlled by localization of scaffold proteins that bind F-actin or their myosin motors we also apply our 2D model to recapitulate in vitro studies that have revealed complex patterns of actomyosin that assemble from patterning filaments and motor complexes with microcontact printing. Although we use a minimal representation of filament, motor, and cross-linker biophysics, our model establishes a framework for investigating the role of other actin binding proteins, how they might alter actomyosin dynamics, and makes predictions that can be tested experimentally within live cells as well as within in vitro models., Author summary Recent genetic and mechanical studies of embryonic development reveal a critical role for intracellular scaffolds in generating the shape of the embryo and constructing internal organs. In this paper we developed computer simulations of these scaffolds, composed of filamentous actin (F-actin), a rod-like protein polymer, and mini-thick filaments, composed of non-muscle myosin II, forming a two headed spring-like complex of motor proteins that can walk on, and remodel F-actin networks. Using simulations of these dynamic interactions, we can carry out virtual experiments where we change the physics and chemistry of F-actin polymers, their associated myosin motors, and cross-linkers and observe the changes in scaffolds that emerge. For example, by modulating the motor stiffness of the myosin motors in our model we can observe the formation or loss of large aster-like structures. Such fine-grained control over the physical properties of motors or filaments within simulations are not currently possible with biological experiments, even where mutant proteins or small molecule inhibitors can be targeted to specific sites on filaments or motors. Our approach reflects a growing adoption of simulation methods to investigate microscopic features that shape actomyosin arrays and the mesoscale effects of molecular scale processes. We expect predictions from these models will drive more refined experimental approaches to expose the many roles of actomyosin in development.

Single-cell RNA sequencing (scRNA-seq) has been widely applied to discover new cell types by detecting sub-populations in a heterogeneous group of cells. Since scRNA-seq experiments have lower read coverage/tag counts and introduce more technical biases compared to bulk RNA-seq experiments, the limited number of sampled cells combined with the experimental biases and other dataset specific variations presents a challenge to cross-dataset analysis and discovery of relevant biological variations across multiple cell populations. In this paper, we introduce a method of variance-driven multitask clustering of single-cell RNA-seq data (scVDMC) that utilizes multiple single-cell populations from biological replicates or different samples. scVDMC clusters single cells in multiple scRNA-seq experiments of similar cell types and markers but varying expression patterns such that the scRNA-seq data are better integrated than typical pooled analyses which only increase the sample size. By controlling the variance among the cell clusters within each dataset and across all the datasets, scVDMC detects cell sub-populations in each individual experiment with shared cell-type markers but varying cluster centers among all the experiments. Applied to two real scRNA-seq datasets with several replicates and one large-scale droplet-based dataset on three patient samples, scVDMC more accurately detected cell populations and known cell markers than pooled clustering and other recently proposed scRNA-seq clustering methods. In the case study applied to in-house Recessive Dystrophic Epidermolysis Bullosa (RDEB) scRNA-seq data, scVDMC revealed several new cell types and unknown markers validated by flow cytometry. MATLAB/Octave code available at https://github.com/kuanglab/scVDMC., Author summary scRNA-seq enables detailed profiling of heterogeneous cell populations and can be used to reveal lineage relationships or discover new cell types. In the literature, there has been little effort directed towards developing computational methods for cross-population transcriptome analysis of multiple single-cell populations. The cross-cell-population clustering problem is different from the traditional clustering problem because single-cell populations can be collected from different patients, different samples of a tissue, or different experimental replicates. The accompanying biological and technical variation tends to dominate the signals for clustering the pooled single cells from the multiple populations. In this work, we have developed a multitask clustering method to address the cross-population clustering problem. The method simultaneously clusters each individual cell population and controls variance among the cell-type cluster centers within each cell population and across the cell populations. We demonstrate that our multitask clustering method significantly improves clustering accuracy and marker discovery in three public scRNA-seq datasets and also apply the method to an in-house Recessive Dystrophic Epidermolysis Bullosa (RDEB) dataset. Our results make it evident that multitask clustering is a promising new approach for cross-population analysis of scRNA-seq data.

Simultaneous changes in ion concentrations, glutamate, and cell volume together with exchange of matter between cell network and vasculature are ubiquitous in numerous brain pathologies. A complete understanding of pathological conditions as well as normal brain function, therefore, hinges on elucidating the molecular and cellular pathways involved in these mostly interdependent variations. In this paper, we develop the first computational framework that combines the Hodgkin–Huxley type spiking dynamics, dynamic ion concentrations and glutamate homeostasis, neuronal and astroglial volume changes, and ion exchange with vasculature into a comprehensive model to elucidate the role of glutamate uptake in the dynamics of spreading depolarization (SD)—the electrophysiological event underlying numerous pathologies including migraine, ischemic stroke, aneurysmal subarachnoid hemorrhage, intracerebral hematoma, and trauma. We are particularly interested in investigating the role of glutamate in the duration and termination of SD caused by K+ perfusion and oxygen-glucose deprivation. Our results demonstrate that glutamate signaling plays a key role in the dynamics of SD, and that impaired glutamate uptake leads to recovery failure of neurons from SD. We confirm predictions from our model experimentally by showing that inhibiting astrocytic glutamate uptake using TFB-TBOA nearly quadruples the duration of SD in layers 2-3 of visual cortical slices from juvenile rats. The model equations are either derived purely from first physical principles of electroneutrality, osmosis, and conservation of particles or a combination of these principles and known physiological facts. Accordingly, we claim that our approach can be used as a future guide to investigate the role of glutamate, ion concentrations, and dynamics cell volume in other brain pathologies and normal brain function., Author summary Pathological conditions such as seizure, migraine, traumatic brain injury, and stroke are associated with extreme changes in ion concentrations and glutamate, cell swelling, and heavy exchange of matter between neurons, glia, and vasculature. However, current experimental tools are capable of measuring only a few of these variables, which necessitates the development of biophysically relevant models. This study provides a comprehensive computational framework derived from first physical principles and physiological facts that enables us to investigate a wide range of key variables in SD caused by K+ perfusion and oxygen-glucose deprivation. While we use SD as an example, our approach can be extended to other pathological conditions and normal brain function.

Neural circuit development requires that synapses be formed between appropriate neurons. In addition, for a hierarchical network, successful development involves a sequencing of developmental events. It has been suggested that one mechanism that helps speed up development of proper connections is an early overproduction of synapses. Using a computational model of synapse development, such as adaptive synaptogenesis, it is possible to study such overproduction and its role in speeding up development; it is also possible to study other outcomes of synapse overproduction that are seemingly new to the literature. With a fixed number of neurons, adaptive synaptogenesis can control the speed of synaptic development in two ways: by altering the rate constants of the adaptive processes or by altering the initial number of rapidly but non-selectively accrued synapses. Using either mechanism, the simulations reveal that synapse overproduction appears as an unavoidable concomitant of rapid adaptive synaptogenesis. However, the shortest development times, which always produces the greatest amount of synapse overproduction, reduce adult performance by three measures: energy use, discrimination error rates, and proportional neuron allocation. Thus, the results here lead to the hypothesis that the observed speed of neural network development represents a particular inter-generational compromise: quick development benefits parental fecundity while slow development benefits offspring fecundity., Author summary The brain processes information by building and adjusting neural connections, and as an organism learns about the environment, it incorporates associations between features of the environment into these adjustable connections. Organisms, however, are not born with these connections; synapses must develop from a connectionless set of neurons. An important aspect of normal brain development is the timely development of the brain and of hierarchically arranged brain regions. In this paper, we show that an unsupervised algorithm, called adaptive synaptogenesis, builds neurons with stable connections and allows these neurons to appropriately discriminate patterns. Here, we show that larger rates of synaptogenesis and synaptic modification speed up development of stable connections. However, such a speed-up produces larger, steady-state energy-costs. We conjecture that these synaptic modification rates have evolved as a compromise via natural selection. In this regard, the evolution-theoretic idea of parent-offspring conflict arises as a developmental consideration; that is, time to stable connectivity must balance with subsequent neural performance. Importantly, we illustrate that the widely-observed phenomenon of synapse overproduction during development can be understood as a compromise between speedy and efficient development.

A complex disease has, by definition, multiple genetic causes. In theory, these causes could be identified individually, but their identification will likely benefit from informed use of anticipated interactions between causes. In addition, characterizing and understanding interactions must be considered key to revealing the etiology of any complex disease. Large-scale collaborative efforts are now paving the way for comprehensive studies of interaction. As a consequence, there is a need for methods with a computational efficiency sufficient for modern data sets as well as for improvements of statistical accuracy and power. Another issue is that, currently, the relation between different methods for interaction inference is in many cases not transparent, complicating the comparison and interpretation of results between different interaction studies. In this paper we present computationally efficient tests of interaction for the complete family of generalized linear models (GLMs). The tests can be applied for inference of single or multiple interaction parameters, but we show, by simulation, that jointly testing the full set of interaction parameters yields superior power and control of false positive rate. Based on these tests we also describe how to combine results from multiple independent studies of interaction in a meta-analysis. We investigate the impact of several assumptions commonly made when modeling interactions. We also show that, across the important class of models with a full set of interaction parameters, jointly testing the interaction parameters yields identical results. Further, we apply our method to genetic data for cardiovascular disease. This allowed us to identify a putative interaction involved in Lp(a) plasma levels between two ‘tag’ variants in the LPA locus (p = 2.42 ⋅ 10−09) as well as replicate the interaction (p = 6.97 ⋅ 10−07). Finally, our meta-analysis method is used in a small (N = 16,181) study of interactions in myocardial infarction., Author summary Interaction between organic molecules forms the basis of all biological systems. The availability of high-throughput genotyping and sequencing platforms enables us to cost-effectively genotype a large number of individuals. For sufficiently large datasets it is possible to reconstruct the genetic dependencies that underlie complex traits and diseases. However, there is a need for efficient statistical methodologies that can tackle the large sample size and computational resources required to study interaction. In this work we provide theory that reduces the required computational resources, and enable multiple research groups to effectively combine their results.

Mitochondrial fatty-acid beta-oxidation (mFAO) plays a central role in mammalian energy metabolism. Multiple severe diseases are associated with defects in this pathway. Its kinetic structure is characterized by a complex wiring of which the functional implications have hardly been explored. Repetitive cycles of reversible reactions, each cycle shortening the fatty acid by two carbon atoms, evoke competition between intermediates of different chain lengths for a common set of ‘promiscuous’ enzymes (enzymes with activity towards multiple substrates). In our validated kinetic model of the pathway, substrate overload causes a steep and detrimental flux decline. Here, we unravel the underlying mechanism and the role of enzyme promiscuity in it. Comparison of alternative model versions elucidated the role of promiscuity of individual enzymes. Promiscuity of the last enzyme of the pathway, medium-chain ketoacyl-CoA thiolase (MCKAT), was both necessary and sufficient to elicit the flux decline. Subsequently, Metabolic Control Analysis revealed that MCKAT had insufficient capacity to cope with high substrate influx. Next, we quantified the internal metabolic regulation, revealing a vicious cycle around MCKAT. Upon substrate overload, MCKAT’s ketoacyl-CoA substrates started to accumulate. The unfavourable equilibrium constant of the preceding enzyme, medium/short-chain hydroxyacyl-CoA dehydrogenase, worked as an amplifier, leading to accumulation of upstream CoA esters, including acyl-CoA esters. These acyl-CoA esters are at the same time products of MCKAT and inhibited its already low activity further. Finally, the accumulation of CoA esters led to a sequestration of free CoA. CoA being a cofactor for MCKAT, its sequestration limited the MCKAT activity even further, thus completing the vicious cycle. Since CoA is also a substrate for distant enzymes, it efficiently communicated the ‘traffic jam’ at MCKAT to the entire pathway. This novel mechanism provides a basis to explore the role of mFAO in disease and elucidate similar principles in other pathways of lipid metabolism., Author summary Burning of fat is essential to provide energy to our bodies. However, fat metabolism is also involved in well-known diseases, such as metabolic syndrome and type 2 diabetes. In severe, inherited deficiencies of enzymes in the fatty-acid oxidation pathway, children run the risk of life-threatening low glucose levels or ‘hypoglycemia’. This typically occurs during overnight fasting, when the body depends on fat stores for its energy supply. To understand the sequence of events leading to hypoglycemia and eventually provide novel therapies, we simulated the process with a validated computer model. We found that when the supply of fat to the liver was increased beyond a critical value, the rate of fatty-acid oxidation declined. Under these circumstances, intermediate metabolites accumulated, very much like the situation in a traffic jam. The risk of a metabolic traffic jam was enhanced in fatty-acid oxidation disorders. In the body this would lead to insufficient energy to maintain blood glucose levels. In this paper we unravel the vicious cycle that leads to the detrimental flux decline and identify a critical role for the last enzyme of the pathway, MCKAT. Our results suggests that MCKAT or its regulators may be promising drug targets.

Cells use surface receptors to estimate concentrations of external ligands. Limits on the accuracy of such estimations have been well studied for pairs of ligand and receptor species. However, the environment typically contains many ligands, which can bind to the same receptors with different affinities, resulting in cross-talk. In traditional rate models, such cross-talk prevents accurate inference of concentrations of individual ligands. In contrast, here we show that knowing the precise timing sequence of stochastic binding and unbinding events allows one receptor to provide information about multiple ligands simultaneously and with a high accuracy. We show that such high-accuracy estimation of multiple concentrations can be realized with simple structural modifications of the familiar kinetic proofreading biochemical network diagram. We give two specific examples of such modifications. We argue that structural and functional features of real cellular biochemical sensory networks in immune cells, such as feedforward and feedback loops or ligand antagonism, sometimes can be understood as solutions to the accurate multi-ligand estimation problem., Author summary Cells live in chemically complex environments with many different chemical ligands around them. Can cells estimate concentrations of more ligands than they have receptor types? In this paper, we show that, surprisingly, the answer is “yes”, and the estimation can be implemented with simple biochemical components already present in many cells. Therefore, cells may “know” a lot more about their environment and thus may be able to implement more complex and accurate response strategies than was previously thought.

The genomic information of microbes is a major determinant of their phenotypic properties, yet it is largely unknown to what extent ecological associations between different species can be explained by their genome composition. To bridge this gap, this study introduces two new genome-wide pairwise measures of microbe-microbe interaction. The first (genome content similarity index) quantifies similarity in genome composition between two microbes, while the second (microbe-microbe functional association index) summarizes the topology of a protein functional association network built for a given pair of microbes and quantifies the fraction of network edges crossing organismal boundaries. These new indices are then used to predict co-occurrence between reference genomes from two 16S-based ecological datasets, accounting for phylogenetic relatedness of the taxa. Phylogenetic relatedness was found to be a strong predictor of ecological associations between microbes which explains about 10% of variance in co-occurrence data, but genome composition was found to be a strong predictor as well, it explains up to 4% the variance in co-occurrence when all genomic-based indices are used in combination, even after accounting for evolutionary relationships between the species. On their own, the metrics proposed here explain a larger proportion of variance than previously reported more complex methods that rely on metabolic network comparisons. In summary, results of this study indicate that microbial genomes do indeed contain detectable signal of organismal ecology, and the methods described in the paper can be used to improve mechanistic understanding of microbe-microbe interactions., Author summary It is still unknown to what extent ecological associations between microbes, as measured by co-occurrence of different taxa in 16S rRNA surveys, can be explained, or predicted, using composition and structure of microbial genomes alone. Here I introduce two new genome-wide, pairwise indices for quantifying the propensity of microbial species to interact with each other. The first measure quantifies similarity in genome composition between two microbes. The second measure summarizes the topology of a protein functional association network built for a given pair of microbes and quantifies the fraction of network edges crossing organismal boundaries. I then study the ability of two newly proposed and two previously reported indices to explain variation in microbial co-occurrence. All four measures are significantly correlated with co-occurrence of microbes even when accounting for evolutionary relationships between the species. One of the newly developed indices outperforms previously proposed ones and explains up to 3.5% of the variance in co-occurrence. In summary, the indices described here are able to detect ecological associations between species using only their genomic information; however, additional methods are needed to provide more reliable genomic tools for microbial ecology.

Accurate prediction of active sites is an important tool in bioinformatics. Here we present an improved structure based technique to expose active sites that is based on large changes of solvent accessibility accompanying normal mode dynamics. The technique which detects EXPOsure of active SITes through normal modEs is named EXPOSITE. The technique is trained using a small 133 enzyme dataset and tested using a large 845 enzyme dataset, both with known active site residues. EXPOSITE is also tested in a benchmark protein ligand dataset (PLD) comprising 48 proteins with and without bound ligands. EXPOSITE is shown to successfully locate the active site in most instances, and is found to be more accurate than other structure-based techniques. Interestingly, in several instances, the active site does not correspond to the largest pocket. EXPOSITE is advantageous due to its high precision and paves the way for structure based prediction of active site in enzymes., Author Summary In this paper, we present an improved technique to predict active sites in enzymes. Our technique is based on changes of solvent accessibility that accompany normal mode dynamics. We assert the technique strength using several enzyme datasets with known catalytic residues. We show the technique successfully locates the active site in most cases, and consistently surpasses the accuracy of other techniques. We show how the technique is advantageous and paves the way for high precision prediction of active sites.

Cross-feeding, a relationship wherein one organism consumes metabolites excreted by another, is a ubiquitous feature of natural and clinically-relevant microbial communities and could be a key factor promoting diversity in extreme and/or nutrient-poor environments. However, it remains unclear how readily cross-feeding interactions form, and therefore our ability to predict their emergence is limited. In this paper we developed a mathematical model parameterized using data from the biochemistry and ecology of an E. coli cross-feeding laboratory system. The model accurately captures short-term dynamics of the two competitors that have been observed empirically and we use it to systematically explore the stability of cross-feeding interactions for a range of environmental conditions. We find that our simple system can display complex dynamics including multi-stable behavior separated by a critical point. Therefore whether cross-feeding interactions form depends on the complex interplay between density and frequency of the competitors as well as on the concentration of resources in the environment. Moreover, we find that subtly different environmental conditions can lead to dramatically different results regarding the establishment of cross-feeding, which could explain the apparently unpredictable between-population differences in experimental outcomes. We argue that mathematical models are essential tools for disentangling the complexities of cross-feeding interactions., Author Summary Simple environments, even those used in laboratory experimental evolution, have proven vastly richer than originally thought, capable of generating and supporting genetic and phenotypic diversity. This was not foreseen by Gause’s seminal competitive exclusion theory, which predicted that simple single niche environments cannot support diversity. We now know that cross-feeding interactions can be a major driver of diversity maintenance in simple environments. Cross-feeding, a relationship wherein one organism consumes metabolites excreted by another, is a ubiquitous feature of natural and clinically-relevant microbial communities and even tumour cell populations. However, it remains unclear how readily such relationships form, and therefore our ability to predict their emergence is limited. Here we developed a mathematical model of cross-feeding and find that this system can display complex dynamics including multi-stable behaviour separated by a critical point. Therefore, the emergence of cross-feeding depends on complex interplay between density and frequency of competitors. Moreover we predict that small changes in environmental conditions can cause abrupt and irreversible shifts from cross-feeding permissive to cross-feeding prohibitive states. We argue that mathematical models are essential tools for disentangling the complexities of cross-feeding interactions.

Across the nervous system, certain population spiking patterns are observed far more frequently than others. A hypothesis about this structure is that these collective activity patterns function as population codewords–collective modes–carrying information distinct from that of any single cell. We investigate this phenomenon in recordings of ∼150 retinal ganglion cells, the retina’s output. We develop a novel statistical model that decomposes the population response into modes; it predicts the distribution of spiking activity in the ganglion cell population with high accuracy. We found that the modes represent localized features of the visual stimulus that are distinct from the features represented by single neurons. Modes form clusters of activity states that are readily discriminated from one another. When we repeated the same visual stimulus, we found that the same mode was robustly elicited. These results suggest that retinal ganglion cells’ collective signaling is endowed with a form of error-correcting code–a principle that may hold in brain areas beyond retina., Author Summary Neurons in most parts of the nervous system represent and process information in a collective fashion, yet the nature of this collective code is poorly understood. An important constraint placed on any such collective processing comes from the fact that individual neurons’ signaling is prone to corruption by noise. The information theory and engineering literatures have studied error-correcting codes that allow individual noise-prone coding units to “check” each other, forming an overall representation that is robust to errors. In this paper, we have analyzed the population code of one of the best-studied neural systems, the retina, and found that it is structured in a manner analogous to error-correcting schemes. Indeed, we found that the complex activity patterns over ~150 retinal ganglion cells, the output neurons of the retina, could be mapped onto collective code words, and that these code words represented precise visual information while suppressing noise. In order to analyze this coding scheme, we introduced a novel quantitative model of the retinal output that predicted neural activity patterns more accurately than existing state-of-the-art approaches.

We present a computational framework, called DISCERN (DIfferential SparsE Regulatory Network), to identify informative topological changes in gene-regulator dependence networks inferred on the basis of mRNA expression datasets within distinct biological states. DISCERN takes two expression datasets as input: an expression dataset of diseased tissues from patients with a disease of interest and another expression dataset from matching normal tissues. DISCERN estimates the extent to which each gene is perturbed—having distinct regulator connectivity in the inferred gene-regulator dependencies between the disease and normal conditions. This approach has distinct advantages over existing methods. First, DISCERN infers conditional dependencies between candidate regulators and genes, where conditional dependence relationships discriminate the evidence for direct interactions from indirect interactions more precisely than pairwise correlation. Second, DISCERN uses a new likelihood-based scoring function to alleviate concerns about accuracy of the specific edges inferred in a particular network. DISCERN identifies perturbed genes more accurately in synthetic data than existing methods to identify perturbed genes between distinct states. In expression datasets from patients with acute myeloid leukemia (AML), breast cancer and lung cancer, genes with high DISCERN scores in each cancer are enriched for known tumor drivers, genes associated with the biological processes known to be important in the disease, and genes associated with patient prognosis, in the respective cancer. Finally, we show that DISCERN can uncover potential mechanisms underlying network perturbation by explaining observed epigenomic activity patterns in cancer and normal tissue types more accurately than alternative methods, based on the available epigenomic data from the ENCODE project., Author Summary Certain genes can regulate other genes’ expression and activity levels to perform key biological processes in a cell. Understanding how genes affect expression levels among one another is a fundamental goal in molecular biology. New statistical techniques that analyze genome-wide mRNA expression data obtained from different individuals can improve our understanding of gene interactions. In this paper, we focus on identifying the genes involved in changes in the topology of gene networks between two distinct biological status, for example, between diseased tissue and normal tissue. Our method, called DISCERN, takes two expression datasets, each from a different condition, for example, cancer and normal tissue, and identifies the genes that are likely to change their neighbors in the estimated gene networks between disease and normal conditions. Most analysis methods that compare gene expression datasets from two conditions address the question of which genes are significantly differentially expressed between conditions. The DISCERN method addresses a distinct question concerning which genes are significantly rewired in the inferred gene-regulator network in disease tissues. We show that DISCERN successfully identifies known disease specific genes, for example, genes whose expression levels are significantly associated with survival time in cancer, genes that have known disease specific driver mutations, as well as genes that have epigenomic evidence of regulatory rewiring between cancer and normal tissue types.

Allostery is one of the pervasive mechanisms through which proteins in living systems carry out enzymatic activity, cell signaling, and metabolism control. Effective modeling of the protein function regulation requires a synthesis of the thermodynamic and structural views of allostery. We present here a structure-based statistical mechanical model of allostery, allowing one to observe causality of communication between regulatory and functional sites, and to estimate per residue free energy changes. Based on the consideration of ligand free and ligand bound systems in the context of a harmonic model, corresponding sets of characteristic normal modes are obtained and used as inputs for an allosteric potential. This potential quantifies the mean work exerted on a residue due to the local motion of its neighbors. Subsequently, in a statistical mechanical framework the entropic contribution to allosteric free energy of a residue is directly calculated from the comparison of conformational ensembles in the ligand free and ligand bound systems. As a result, this method provides a systematic approach for analyzing the energetics of allosteric communication based on a single structure. The feasibility of the approach was tested on a variety of allosteric proteins, heterogeneous in terms of size, topology and degree of oligomerization. The allosteric free energy calculations show the diversity of ways and complexity of scenarios existing in the phenomenology of allosteric causality and communication. The presented model is a step forward in developing the computational techniques aimed at detecting allosteric sites and obtaining the discriminative power between agonistic and antagonistic effectors, which are among the major goals in allosteric drug design., Author Summary The 50th anniversary of Monod-Changeux-Jacob seminal paper “Allosteric proteins and cellular control systems” became the hallmark of a new wave in the allostery studies and the turning point in our vision of allostery and its implications in protein engineering and drug design. Recent experimental and theoretical works clearly show relevance of allosteric phenomenon to drug design, unraveling advantages of allosteric drugs in comparison to traditional orthosteric compounds. Remarkable simplicity of allosteric effectors and, at the same time, their potentially high specificity is one of the most important traits. The non conserved nature of allosteric ligands is a basis for avoiding drug resistance, and existence of latent regulatory sites make them attractive drug targets. The model presented in this work provides a theoretical framework for the quantification of the causality and energetics of allosteric regulation, which is a prerequisite for design of effector molecules with required characteristics. The synthesis between the thermodynamics of allostery and the intrinsic atomic nature of proteins and their interactions with the allosteric effectors accomplished in this work is a small initial step in the long endeavor towards future allosteric drugs.

Development of high-throughput monitoring technologies enables interrogation of cancer samples at various levels of cellular activity. Capitalizing on these developments, various public efforts such as The Cancer Genome Atlas (TCGA) generate disparate omic data for large patient cohorts. As demonstrated by recent studies, these heterogeneous data sources provide the opportunity to gain insights into the molecular changes that drive cancer pathogenesis and progression. However, these insights are limited by the vast search space and as a result low statistical power to make new discoveries. In this paper, we propose methods for integrating disparate omic data using molecular interaction networks, with a view to gaining mechanistic insights into the relationship between molecular changes at different levels of cellular activity. Namely, we hypothesize that genes that play a role in cancer development and progression may be implicated by neither frequent mutation nor differential expression, and that network-based integration of mutation and differential expression data can reveal these “silent players”. For this purpose, we utilize network-propagation algorithms to simulate the information flow in the cell at a sample-specific resolution. We then use the propagated mutation and expression signals to identify genes that are not necessarily mutated or differentially expressed genes, but have an essential role in tumor development and patient outcome. We test the proposed method on breast cancer and glioblastoma multiforme data obtained from TCGA. Our results show that the proposed method can identify important proteins that are not readily revealed by molecular data, providing insights beyond what can be gleaned by analyzing different types of molecular data in isolation., Author Summary Identification of cancer-related genes is an important task, made more difficult by heterogeneity between samples and even within individual patients. Methods for identifying disease-related genes typically focus on individual data sets such as mutational and differential expression data, and therefore are limited to genes that are implicated by each data set in isolation. In this work we propose a method that uses protein interaction network information to integrate mutational and differential expression data on a sample-specific level, and combine this information across samples in ways that respect the commonalities and differences between distinct mutation and differential expression profiles. We use this information to identify genes that are associated with cancer but not readily identifiable by mutations or differential expression alone. Our method highlights the features that significantly predict a gene’s association with cancer, shows improved predictive power in recovering cancer-related genes in known pathways, and identifies genes that are neither frequently mutated nor differentially expressed but show significant association with survival.

The Gene Ontology (GO) provides biologists with a controlled terminology that describes how genes are associated with functions and how functional terms are related to one another. These term-term relationships encode how scientists conceive the organization of biological functions, and they take the form of a directed acyclic graph (DAG). Here, we propose that the network structure of gene-term annotations made using GO can be employed to establish an alternative approach for grouping functional terms that captures intrinsic functional relationships that are not evident in the hierarchical structure established in the GO DAG. Instead of relying on an externally defined organization for biological functions, our approach connects biological functions together if they are performed by the same genes, as indicated in a compendium of gene annotation data from numerous different sources. We show that grouping terms by this alternate scheme provides a new framework with which to describe and predict the functions of experimentally identified sets of genes., Author Summary Investigating how a set of genes might collectively work together to perform various cellular processes has become a routine part of many biological analyses. In such analyses, genes of interest are compared to sets of genes annotated to various biological functions (or pathways) defined within carefully curated databases. One of the most comprehensive and widely used resources of this type is the Gene Ontology (GO) database. The Gene Ontology database is comprised of two important elements: (1) the ontology itself, which provides a controlled vocabulary of terms describing genetic function and also specifies how these functional terms are related to one another via a hierarchical structure; and (2) the set of annotations made using GO that connect individual genes to different functional terms. In our paper we investigate a method for organizing functional terms that results from connecting terms based on shared gene annotations. We find that this alternate classification has an organization that is highly distinct from the Gene Ontology hierarchy, challenging the way we think about the relationships between different biological functions. Finally, we show that these alternate collections of terms are highly associated with published cancer gene signatures, demonstrating that this alternative organization of biological functions can highlight important relationships between cellular processes and has the potential to lead to new insights and discoveries.

In simple organisms like E.coli, the metabolic response to an external perturbation passes through a transient phase in which the activation of a number of latent pathways can guarantee survival at the expenses of growth. Growth is gradually recovered as the organism adapts to the new condition. This adaptation can be modeled as a process of repeated metabolic adjustments obtained through the resilencings of the non-essential metabolic reactions, using growth rate as selection probability for the phenotypes obtained. The resulting metabolic adaptation process tends naturally to steer the metabolic fluxes towards high growth phenotypes. Quite remarkably, when applied to the central carbon metabolism of E.coli, it follows that nearly all flux distributions converge to the flux vector representing optimal growth, i.e., the solution of the biomass optimization problem turns out to be the dominant attractor of the metabolic adaptation process., Author Summary In modeling metabolic networks, concepts like biomass optimization are often used to determine flux distributions of simple organisms such as E.coli. Although they often give good results in practice, they normally rely on heuristic considerations like “evolution has tuned metabolic fluxes to optimize growth, hence optimizing growth gives reasonable fluxes”. The main result of this paper is to show that metabolic adaptation naturally leads to optimal growth, in the sense that the flux distribution associated to optimal growth is the dominant attractor of the fitness landscape of the metabolic adaptation process.