Search

Isogenic cells cultured together show heterogeneity in their proliferation rate. To determine the differences between fast and slow-proliferating cells, we developed a method to sort cells by proliferation rate, and performed RNA-seq on slow and fast proliferating subpopulations of pluripotent mouse embryonic stem cells (mESCs) and mouse fibroblasts. We found that slowly proliferating mESCs have a more naïve pluripotent character. We identified an evolutionarily conserved proliferation-correlated transcriptomic signature that is common to all eukaryotes: fast cells have higher expression of genes for protein synthesis and protein degradation. This signature accurately predicted growth rate in yeast and cancer cells, and identified lineage-specific proliferation dynamics during development, using C. elegans scRNA-seq data. In contrast, sorting by mitochondria membrane potential revealed a highly cell-type specific mitochondria-state related transcriptome. mESCs with hyperpolarized mitochondria are fast proliferating, while the opposite is true for fibroblasts. The mitochondrial electron transport chain inhibitor antimycin affected slow and fast subpopulations differently. While a major transcriptional-signature associated with cell-to-cell heterogeneity in proliferation is conserved, the metabolic and energetic dependency of cell proliferation is cell-type specific.

Generalist microbes have adapted to a multitude of environmental stresses through their integrated stress response system. Individual stress responses have been quantified by E. coli metabolism and expression (ME) models under thermal, oxidative and acid stress, respectively. However, the systematic quantification of cross-stress & cross-talk among these stress responses remains lacking. Here, we present StressME: the unified stress response model of E. coli combining thermal (FoldME), oxidative (OxidizeME) and acid (AcidifyME) stress responses. StressME is the most up to date ME model for E. coli and it reproduces all published single-stress ME models. Additionally, it includes refined rate constants to improve prediction accuracy for wild-type and stress-evolved strains. StressME revealed certain optimal proteome allocation strategies associated with cross-stress and cross-talk responses. These stress-optimal proteomes were shaped by trade-offs between protective vs. metabolic enzymes; cytoplasmic vs. periplasmic chaperones; and expression of stress-specific proteins. As StressME is tuned to compute metabolic and gene expression responses under mild acid, oxidative, and thermal stresses, it is useful for engineering and health applications. The modular design of our open-source package also facilitates model expansion (e.g., to new stress mechanisms) by the computational biology community.

Identifying plausible mechanisms for compartmentalization and accumulation of the organic intermediates of early metabolic cycles in primitive cells has been a major challenge in theories of life's origins. Here, we propose a mechanism, where positive membrane potentials elevate the concentration of the organic intermediates. Positive membrane potentials are generated by positively charged surfaces of protocell membranes due to accumulation of transition metals. We find that (i) positive membrane potentials comparable in magnitude to those of modern cells can increase the concentration of the organic intermediates by several orders of magnitude; (ii) generation of large membrane potentials destabilize ion distributions; (iii) violation of electroneutrality is necessary to induce nonzero membrane potentials; and (iv) violation of electroneutrality enhances osmotic pressure and diminishes reaction efficiency, resulting in an evolutionary driving force for the formation of lipid membranes, specialized ion channels, and active transport systems.

Sustaining a robust metabolic network requires a balanced and fully functioning proteome. In addition to amino acids, many enzymes require cofactors (coenzymes and engrafted prosthetic groups) to function properly. Extensively validated resource allocation models, such as genome-scale models of metabolism and gene expression (ME-models), have the ability to compute an optimal proteome composition underlying a metabolic phenotype, including the provision of all required cofactors. Here we apply the ME-model for Escherichia coli K-12 MG1655 to computationally examine how environmental conditions change the proteome and its accompanying cofactor usage. We found that: (1) The cofactor requirements computed by the ME-model mostly agree with the standard biomass objective function used in models of metabolism alone (M-models); (2) ME-model computations reveal non-intuitive variability in cofactor use under different growth conditions; (3) An analysis of ME-model predicted protein use in aerobic and anaerobic conditions suggests an enrichment in the use of peroxyl scavenging acids in the proteins used to sustain aerobic growth; (4) The ME-model could describe how limitation in key protein components affect the metabolic state of E. coli. Genome-scale models have thus reached a level of sophistication where they reveal intricate properties of functional proteomes and how they support different E. coli lifestyles.

Mathematical models of metabolic networks utilize simulation to study system-level mechanisms and functions. Various approaches have been used to model the steady state behavior of metabolic networks using genome-scale reconstructions, but formulating dynamic models from such reconstructions continues to be a key challenge. Here, we present the Mass Action Stoichiometric Simulation Python (MASSpy) package, an open-source computational framework for dynamic modeling of metabolism. MASSpy utilizes mass action kinetics and detailed chemical mechanisms to build dynamic models of complex biological processes. MASSpy adds dynamic modeling tools to the COnstraint-Based Reconstruction and Analysis Python (COBRApy) package to provide an unified framework for constraint-based and kinetic modeling of metabolic networks. MASSpy supports high-performance dynamic simulation through its implementation of libRoadRunner: the Systems Biology Markup Language (SBML) simulation engine. Three examples are provided to demonstrate how to use MASSpy: (1) a validation of the MASSpy modeling tool through dynamic simulation of detailed mechanisms of enzyme regulation; (2) a feature demonstration using a workflow for generating ensemble of kinetic models using Monte Carlo sampling to approximate missing numerical values of parameters and to quantify biological uncertainty, and (3) a case study in which MASSpy is utilized to overcome issues that arise when integrating experimental data with the computation of functional states of detailed biological mechanisms. MASSpy represents a powerful tool to address challenges that arise in dynamic modeling of metabolic networks, both at small and large scales.

The fitness landscape is a concept commonly used to describe evolution towards optimal phenotypes. It can be reduced to mechanistic detail using genome-scale models (GEMs) from systems biology. We use recently developed GEMs of Metabolism and protein Expression (ME-models) to study the distribution of Escherichia coli phenotypes on the rate-yield plane. We found that the measured phenotypes distribute non-uniformly to form a highly stratified fitness landscape. Systems analysis of the ME-model simulations suggest that this stratification results from discrete ATP generation strategies. Accordingly, we define "aero-types", a phenotypic trait that characterizes how a balanced proteome can achieve a given growth rate by modulating 1) the relative utilization of oxidative phosphorylation, glycolysis, and fermentation pathways; and 2) the differential employment of electron-transport-chain enzymes. This global, quantitative, and mechanistic systems biology interpretation of fitness landscape formed upon proteome allocation offers a fundamental understanding of bacterial physiology and evolution dynamics.

The evolution of antimicrobial resistance (AMR) poses a persistent threat to global public health. Sequencing efforts have already yielded genome sequences for thousands of resistant microbial isolates and require robust computational tools to systematically elucidate the genetic basis for AMR. Here, we present a generalizable machine learning workflow for identifying genetic features driving AMR based on constructing reference strain-agnostic pan-genomes and training random subspace ensembles (RSEs). This workflow was applied to the resistance profiles of 14 antimicrobials across three urgent threat pathogens encompassing 288 Staphylococcus aureus, 456 Pseudomonas aeruginosa, and 1588 Escherichia coli genomes. We find that feature selection by RSE detects known AMR associations more reliably than common statistical tests and previous ensemble approaches, identifying a total of 45 known AMR-conferring genes and alleles across the three organisms, as well as 25 candidate associations backed by domain-level annotations. Furthermore, we find that results from the RSE approach are consistent with existing understanding of fluoroquinolone (FQ) resistance due to mutations in the main drug targets, gyrA and parC, in all three organisms, and suggest the mutational landscape of those genes with respect to FQ resistance is simple. As larger datasets become available, we expect this approach to more reliably predict AMR determinants for a wider range of microbial pathogens.

Response to acid stress is critical for Escherichia coli to successfully complete its life-cycle by passing through the stomach to colonize the digestive tract. To develop a fundamental understanding of this response, we established a molecular mechanistic description of acid stress mitigation responses in E. coli and integrated them with a genome-scale model of its metabolism and macromolecular expression (ME-model). We considered three known mechanisms of acid stress mitigation: 1) change in membrane lipid fatty acid composition, 2) change in periplasmic protein stability over external pH and periplasmic chaperone protection mechanisms, and 3) change in the activities of membrane proteins. After integrating these mechanisms into an established ME-model, we could simulate their responses in the context of other cellular processes. We validated these simulations using RNA sequencing data obtained from five E. coli strains grown under external pH ranging from 5.5 to 7.0. We found: i) that for the differentially expressed genes accounted for in the ME-model, 80% of the upregulated genes were correctly predicted by the ME-model, and ii) that these genes are mainly involved in translation processes (45% of genes), membrane proteins and related processes (18% of genes), amino acid metabolism (12% of genes), and cofactor and prosthetic group biosynthesis (8% of genes). We also demonstrated several intervention strategies on acid tolerance that can be simulated by the ME-model. We thus established a quantitative framework that describes, on a genome-scale, the acid stress mitigation response of E. coli that has both scientific and practical uses.

Genome-scale metabolic models (GEMs) are mathematically structured knowledge bases of metabolism that provide phenotypic predictions from genomic information. GEM-guided predictions of growth phenotypes rely on the accurate definition of a biomass objective function (BOF) that is designed to include key cellular biomass components such as the major macromolecules (DNA, RNA, proteins), lipids, coenzymes, inorganic ions and species-specific components. Despite its importance, no standardized computational platform is currently available to generate species-specific biomass objective functions in a data-driven, unbiased fashion. To fill this gap in the metabolic modeling software ecosystem, we implemented BOFdat, a Python package for the definition of a Biomass Objective Function from experimental data. BOFdat has a modular implementation that divides the BOF definition process into three independent modules defined here as steps: 1) the coefficients for major macromolecules are calculated, 2) coenzymes and inorganic ions are identified and their stoichiometric coefficients estimated, 3) the remaining species-specific metabolic biomass precursors are algorithmically extracted in an unbiased way from experimental data. We used BOFdat to reconstruct the BOF of the Escherichia coli model iML1515, a gold standard in the field. The BOF generated by BOFdat resulted in the most concordant biomass composition, growth rate, and gene essentiality prediction accuracy when compared to other methods. Installation instructions for BOFdat are available in the documentation and the source code is available on GitHub (https://github.com/jclachance/BOFdat).

S. aureus is classified as a serious threat pathogen and is a priority that guides the discovery and development of new antibiotics. Despite growing knowledge of S. aureus metabolic capabilities, our understanding of its systems-level responses to different media types remains incomplete. Here, we develop a manually reconstructed genome-scale model (GEM-PRO) of metabolism with 3D protein structures for S. aureus USA300 str. JE2 containing 854 genes, 1,440 reactions, 1,327 metabolites and 673 3-dimensional protein structures. Computations were in 85% agreement with gene essentiality data from random barcode transposon site sequencing (RB-TnSeq) and 68% agreement with experimental physiological data. Comparisons of computational predictions with experimental observations highlight: 1) cases of non-essential biomass precursors; 2) metabolic genes subject to transcriptional regulation involved in Staphyloxanthin biosynthesis; 3) the essentiality of purine and amino acid biosynthesis in synthetic physiological media; and 4) a switch to aerobic fermentation upon exposure to extracellular glucose elucidated as a result of integrating time-course of quantitative exo-metabolomics data. An up-to-date GEM-PRO thus serves as a knowledge-based platform to elucidate S. aureus' metabolic response to its environment.

Allosteric regulation has traditionally been described by mathematically-complex allosteric rate laws in the form of ratios of polynomials derived from the application of simplifying kinetic assumptions. Alternatively, an approach that explicitly describes all known ligand-binding events requires no simplifying assumptions while allowing for the computation of enzymatic states. Here, we employ such a modeling approach to examine the "catalytic potential" of an enzyme-an enzyme's capacity to catalyze a biochemical reaction. The catalytic potential is the fundamental result of multiple ligand-binding events that represents a "tug of war" among the various regulators and substrates within the network. This formalism allows for the assessment of interacting allosteric enzymes and development of a network-level understanding of regulation. We first define the catalytic potential and use it to characterize the response of three key kinases (hexokinase, phosphofructokinase, and pyruvate kinase) in human red blood cell glycolysis to perturbations in ATP utilization. Next, we examine the sensitivity of the catalytic potential by using existing personalized models, finding that the catalytic potential allows for the identification of subtle but important differences in how individuals respond to such perturbations. Finally, we explore how the catalytic potential can help to elucidate how enzymes work in tandem to maintain a homeostatic state. Taken together, this work provides an interpretation and visualization of the dynamic interactions and network-level effects of interacting allosteric enzymes.

Genome-scale models of metabolism and macromolecular expression (ME-models) explicitly compute the optimal proteome composition of a growing cell. ME-models expand upon the well-established genome-scale models of metabolism (M-models), and they enable a new fundamental understanding of cellular growth. ME-models have increased predictive capabilities and accuracy due to their inclusion of the biosynthetic costs for the machinery of life, but they come with a significant increase in model size and complexity. This challenge results in models which are both difficult to compute and challenging to understand conceptually. As a result, ME-models exist for only two organisms (Escherichia coli and Thermotoga maritima) and are still used by relatively few researchers. To address these challenges, we have developed a new software framework called COBRAme for building and simulating ME-models. It is coded in Python and built on COBRApy, a popular platform for using M-models. COBRAme streamlines computation and analysis of ME-models. It provides tools to simplify constructing and editing ME-models to enable ME-model reconstructions for new organisms. We used COBRAme to reconstruct a condensed E. coli ME-model called iJL1678b-ME. This reformulated model gives functionally identical solutions to previous E. coli ME-models while using 1/6 the number of free variables and solving in less than 10 minutes, a marked improvement over the 6 hour solve time of previous ME-model formulations. Errors in previous ME-models were also corrected leading to 52 additional genes that must be expressed in iJL1678b-ME to grow aerobically in glucose minimal in silico media. This manuscript outlines the architecture of COBRAme and demonstrates how ME-models can be created, modified, and shared most efficiently using the new software framework.

Several antimalarial drugs exist, but differences between life cycle stages among malaria species pose challenges for developing more effective therapies. To understand the diversity among stages and species, we reconstructed genome-scale metabolic models (GeMMs) of metabolism for five life cycle stages and five species of Plasmodium spanning the blood, transmission, and mosquito stages. The stage-specific models of Plasmodium falciparum uncovered stage-dependent changes in central carbon metabolism and predicted potential targets that could affect several life cycle stages. The species-specific models further highlight differences between experimental animal models and the human-infecting species. Comparisons between human- and rodent-infecting species revealed differences in thiamine (vitamin B1), choline, and pantothenate (vitamin B5) metabolism. Thus, we show that genome-scale analysis of multiple stages and species of Plasmodium can prioritize potential drug targets that could be both anti-malarials and transmission blocking agents, in addition to guiding translation from non-human experimental disease models.

Deep-coverage metabolomic profiling has revealed a well-defined development of metabolic decay in human red blood cells (RBCs) under cold storage conditions. A set of extracellular biomarkers has been recently identified that reliably defines the qualitative state of the metabolic network throughout this metabolic decay process. Here, we extend the utility of these biomarkers by using them to quantitatively predict the concentrations of other metabolites in the red blood cell. We are able to accurately predict the concentration profile of 84 of the 91 (92%) measured metabolites (p < 0.05) in RBC metabolism using only measurements of these five biomarkers. The median of prediction errors (symmetric mean absolute percent error) across all metabolites was 13%. The ability to predict numerous metabolite concentrations from a simple set of biomarkers offers the potential for the development of a powerful workflow that could be used to evaluate the metabolic state of a biological system using a minimal set of measurements.

Progress in systems medicine brings promise to addressing patient heterogeneity and individualized therapies. Recently, genome-scale models of metabolism have been shown to provide insight into the mechanistic link between drug therapies and systems-level off-target effects while being expanded to explicitly include the three-dimensional structure of proteins. The integration of these molecular-level details, such as the physical, structural, and dynamical properties of proteins, notably expands the computational description of biochemical network-level properties and the possibility of understanding and predicting whole cell phenotypes. In this study, we present a multi-scale modeling framework that describes biological processes which range in scale from atomistic details to an entire metabolic network. Using this approach, we can understand how genetic variation, which impacts the structure and reactivity of a protein, influences both native and drug-induced metabolic states. As a proof-of-concept, we study three enzymes (catechol-O-methyltransferase, glucose-6-phosphate dehydrogenase, and glyceraldehyde-3-phosphate dehydrogenase) and their respective genetic variants which have clinically relevant associations. Using all-atom molecular dynamic simulations enables the sampling of long timescale conformational dynamics of the proteins (and their mutant variants) in complex with their respective native metabolites or drug molecules. We find that changes in a protein's structure due to a mutation influences protein binding affinity to metabolites and/or drug molecules, and inflicts large-scale changes in metabolism.

The costs and benefits of protein expression are balanced through evolution. Expression of un-utilized protein (that have no benefits in the current environment) incurs a quantifiable fitness costs on cellular growth rates; however, the magnitude and variability of un-utilized protein expression in natural settings is unknown, largely due to the challenge in determining environment-specific proteome utilization. We address this challenge using absolute and global proteomics data combined with a recently developed genome-scale model of Escherichia coli that computes the environment-specific cost and utility of the proteome on a per gene basis. We show that nearly half of the proteome mass is unused in certain environments and accounting for the cost of this unused protein expression explains >95% of the variance in growth rates of Escherichia coli across 16 distinct environments. Furthermore, reduction in unused protein expression is shown to be a common mechanism to increase cellular growth rates in adaptive evolution experiments. Classification of the unused protein reveals that the unused protein encodes several nutrient- and stress- preparedness functions, which may convey fitness benefits in varying environments. Thus, unused protein expression is the source of large and pervasive fitness costs that may provide the benefit of hedging against environmental change.

The costs and benefits of protein expression are balanced through evolution. Expression of un-utilized protein (that have no benefits in the current environment) incurs a quantifiable fitness costs on cellular growth rates; however, the magnitude and variability of un-utilized protein expression in natural settings is unknown, largely due to the challenge in determining environment-specific proteome utilization. We address this challenge using absolute and global proteomics data combined with a recently developed genome-scale model of Escherichia coli that computes the environment-specific cost and utility of the proteome on a per gene basis. We show that nearly half of the proteome mass is unused in certain environments and accounting for the cost of this unused protein expression explains >95% of the variance in growth rates of Escherichia coli across 16 distinct environments. Furthermore, reduction in unused protein expression is shown to be a common mechanism to increase cellular growth rates in adaptive evolution experiments. Classification of the unused protein reveals that the unused protein encodes several nutrient- and stress- preparedness functions, which may convey fitness benefits in varying environments. Thus, unused protein expression is the source of large and pervasive fitness costs that may provide the benefit of hedging against environmental change.

Escher is a web application for visualizing data on biological pathways. Three key features make Escher a uniquely effective tool for pathway visualization. First, users can rapidly design new pathway maps. Escher provides pathway suggestions based on user data and genome-scale models, so users can draw pathways in a semi-automated way. Second, users can visualize data related to genes or proteins on the associated reactions and pathways, using rules that define which enzymes catalyze each reaction. Thus, users can identify trends in common genomic data types (e.g. RNA-Seq, proteomics, ChIP)--in conjunction with metabolite- and reaction-oriented data types (e.g. metabolomics, fluxomics). Third, Escher harnesses the strengths of web technologies (SVG, D3, developer tools) so that visualizations can be rapidly adapted, extended, shared, and embedded. This paper provides examples of each of these features and explains how the development approach used for Escher can be used to guide the development of future visualization tools.

Altered metabolism in cancer cells has been viewed as a passive response required for a malignant transformation. However, this view has changed through the recently described metabolic oncogenic factors: mutated isocitrate dehydrogenases (IDH), succinate dehydrogenase (SDH), and fumarate hydratase (FH) that produce oncometabolites that competitively inhibit epigenetic regulation. In this study, we demonstrate in silico predictions of oncometabolites that have the potential to dysregulate epigenetic controls in nine types of cancer by incorporating massive scale genetic mutation information (collected from more than 1,700 cancer genomes), expression profiling data, and deploying Recon 2 to reconstruct context-specific genome-scale metabolic models. Our analysis predicted 15 compounds and 24 substructures of potential oncometabolites that could result from the loss-of-function and gain-of-function mutations of metabolic enzymes, respectively. These results suggest a substantial potential for discovering unidentified oncometabolites in various forms of cancers.

Metabolic network reconstructions represent valuable scaffolds for '-omics' data integration and are used to computationally interrogate network properties. However, they do not explicitly account for the synthesis of macromolecules (i.e., proteins and RNA). Here, we present the first genome-scale, fine-grained reconstruction of Escherichia coli's transcriptional and translational machinery, which produces 423 functional gene products in a sequence-specific manner and accounts for all necessary chemical transformations. Legacy data from over 500 publications and three databases were reviewed, and many pathways were considered, including stable RNA maturation and modification, protein complex formation, and iron-sulfur cluster biogenesis. This reconstruction represents the most comprehensive knowledge base for these important cellular functions in E. coli and is unique in its scope. Furthermore, it was converted into a mathematical model and used to: (1) quantitatively integrate gene expression data as reaction constraints and (2) compute functional network states, which were compared to reported experimental data. For example, the model predicted accurately the ribosome production, without any parameterization. Also, in silico rRNA operon deletion suggested that a high RNA polymerase density on the remaining rRNA operons is needed to reproduce the reported experimental ribosome numbers. Moreover, functional protein modules were determined, and many were found to contain gene products from multiple subsystems, highlighting the functional interaction of these proteins. This genome-scale reconstruction of E. coli's transcriptional and translational machinery presents a milestone in systems biology because it will enable quantitative integration of '-omics' datasets and thus the study of the mechanistic principles underlying the genotype-phenotype relationship.

Recent advances in reconstruction and analytical methods for signaling networks have spurred the development of large-scale models that incorporate fully functional and biologically relevant features. An extended reconstruction of the human Toll-like receptor signaling network is presented herein. This reconstruction contains an extensive complement of kinases, phosphatases, and other associated proteins that mediate the signaling cascade along with a delineation of their associated chemical reactions. A computational framework based on the methods of large-scale convex analysis was developed and applied to this network to characterize input-output relationships. The input-output relationships enabled significant modularization of the network into ten pathways. The analysis identified potential candidates for inhibitory mediation of TLR signaling with respect to their specificity and potency. Subsequently, we were able to identify eight novel inhibition targets through constraint-based modeling methods. The results of this study are expected to yield meaningful avenues for further research in the task of mediating the Toll-like receptor signaling network and its effects.

The study of dynamic functions of large-scale biological networks has intensified in recent years. A critical component in developing an understanding of such dynamics involves the study of their hierarchical organization. We investigate the temporal hierarchy in biochemical reaction networks focusing on: (1) the elucidation of the existence of "pools" (i.e., aggregate variables) formed from component concentrations and (2) the determination of their composition and interactions over different time scales. To date the identification of such pools without prior knowledge of their composition has been a challenge. A new approach is developed for the algorithmic identification of pool formation using correlations between elements of the modal matrix that correspond to a pair of concentrations and how such correlations form over the hierarchy of time scales. The analysis elucidates a temporal hierarchy of events that range from chemical equilibration events to the formation of physiologically meaningful pools, culminating in a network-scale (dynamic) structure-(physiological) function relationship. This method is validated on a model of human red blood cell metabolism and further applied to kinetic models of yeast glycolysis and human folate metabolism, enabling the simplification of these models. The understanding of temporal hierarchy and the formation of dynamic aggregates on different time scales is foundational to the study of network dynamics and has relevance in multiple areas ranging from bacterial strain design and metabolic engineering to the understanding of disease processes in humans.

Public gene expression databases are a rapidly expanding resource of organism responses to diverse perturbations, presenting both an opportunity and a challenge for bioinformatics workflows to extract actionable knowledge of transcription regulatory network function. Here, we introduce a five-step computational pipeline, called iModulonMiner, to compile, process, curate, analyze, and characterize the totality of RNA-seq data for a given organism or cell type. This workflow is centered around the data-driven computation of co-regulated gene sets using Independent Component Analysis, called iModulons, which have been shown to have broad applications. As a demonstration, we applied this workflow to generate the iModulon structure of Bacillus subtilis using all high-quality, publicly-available RNA-seq data. Using this structure, we predicted regulatory interactions for multiple transcription factors, identified groups of co-expressed genes that are putatively regulated by undiscovered transcription factors, and predicted properties of a recently discovered single-subunit phage RNA polymerase. We also present a Python package, PyModulon, with functions to characterize, visualize, and explore computed iModulons. The pipeline, available at https://github.com/SBRG/iModulonMiner, can be readily applied to diverse organisms to gain a rapid understanding of their transcriptional regulatory network structure and condition-specific activity.

Immunization through repeated direct venous inoculation of Plasmodium falciparum (Pf) sporozoites (PfSPZ) under chloroquine chemoprophylaxis, using the PfSPZ Chemoprophylaxis Vaccine (PfSPZ-CVac), induces high-level protection against controlled human malaria infection (CHMI). Humoral and cellular immunity contribute to vaccine efficacy but only limited information about the implicated Pf-specific antigens is available. Here, we examined Pf-specific antibody profiles, measured by protein arrays representing the full Pf proteome, of 40 placebo- and PfSPZ-immunized malaria-naïve volunteers from an earlier published PfSPZ-CVac dose-escalation trial. For this purpose, we both utilized and adapted supervised machine learning methods to identify predictive antibody profiles at two different time points: after immunization and before CHMI. We developed an adapted multitask support vector machine (SVM) approach and compared it to standard methods, i.e. single-task SVM, regularized logistic regression and random forests. Our results show, that the multitask SVM approach improved the classification performance to discriminate the protection status based on the underlying antibody-profiles while combining time- and dose-dependent data in the prediction model. Additionally, we developed the new fEature diStance exPlainabilitY (ESPY) method to quantify the impact of single antigens on the non-linear multitask SVM model and make it more interpretable. In conclusion, our multitask SVM model outperforms the studied standard approaches in regard of classification performance. Moreover, with our new explanation method ESPY, we were able to interpret the impact of Pf-specific antigen antibody responses that predict sterile protective immunity against CHMI after immunization. The identified Pf-specific antigens may contribute to a better understanding of immunity against human malaria and may foster vaccine development.

Generalist microbes have adapted to a multitude of environmental stresses through their integrated stress response system. Individual stress responses have been quantified by E. coli metabolism and expression (ME) models under thermal, oxidative and acid stress, respectively. However, the systematic quantification of cross-stress & cross-talk among these stress responses remains lacking. Here, we present StressME: the unified stress response model of E. coli combining thermal (FoldME), oxidative (OxidizeME) and acid (AcidifyME) stress responses. StressME is the most up to date ME model for E. coli and it reproduces all published single-stress ME models. Additionally, it includes refined rate constants to improve prediction accuracy for wild-type and stress-evolved strains. StressME revealed certain optimal proteome allocation strategies associated with cross-stress and cross-talk responses. These stress-optimal proteomes were shaped by trade-offs between protective vs. metabolic enzymes; cytoplasmic vs. periplasmic chaperones; and expression of stress-specific proteins. As StressME is tuned to compute metabolic and gene expression responses under mild acid, oxidative, and thermal stresses, it is useful for engineering and health applications. The modular design of our open-source package also facilitates model expansion (e.g., to new stress mechanisms) by the computational biology community.

The transcriptional regulatory network (TRN) of E. coli consists of thousands of interactions between regulators and DNA sequences. Regulons are typically determined either from resource-intensive experimental measurement of functional binding sites, or inferred from analysis of high-throughput gene expression datasets. Recently, independent component analysis (ICA) of RNA-seq compendia has shown to be a powerful method for inferring bacterial regulons. However, it remains unclear to what extent regulons predicted by ICA structure have a biochemical basis in promoter sequences. Here, we address this question by developing machine learning models that predict inferred regulon structures in E. coli based on promoter sequence features. Models were constructed successfully (cross-validation AUROC > = 0.8) for 85% (40/47) of ICA-inferred E. coli regulons. We found that: 1) The presence of a high scoring regulator motif in the promoter region was sufficient to specify regulatory activity in 40% (19/47) of the regulons, 2) Additional features, such as DNA shape and extended motifs that can account for regulator multimeric binding, helped to specify regulon structure for the remaining 60% of regulons (28/47); 3) investigating regulons where initial machine learning models failed revealed new regulator-specific sequence features that improved model accuracy. Finally, we found that strong regulatory binding sequences underlie both the genes shared between ICA-inferred and experimental regulons as well as genes in the E. coli core pan-regulon of Fur. This work demonstrates that the structure of ICA-inferred regulons largely can be understood through the strength of regulator binding sites in promoter regions, reinforcing the utility of top-down inference for regulon discovery.

Likelihood ratios are frequently utilized as basis for statistical tests, for model selection criteria and for assessing parameter and prediction uncertainties, e.g. using the profile likelihood. However, translating these likelihood ratios into p-values or confidence intervals requires the exact form of the test statistic's distribution. The lack of knowledge about this distribution for nonlinear ordinary differential equation (ODE) models requires an approximation which assumes the so-called asymptotic setting, i.e. a sufficiently large amount of data. Since the amount of data from quantitative molecular biology is typically limited in applications, this finite-sample case regularly occurs for mechanistic models of dynamical systems, e.g. biochemical reaction networks or infectious disease models. Thus, it is unclear whether the standard approach of using statistical thresholds derived for the asymptotic large-sample setting in realistic applications results in valid conclusions. In this study, empirical likelihood ratios for parameters from 19 published nonlinear ODE benchmark models are investigated using a resampling approach for the original data designs. Their distributions are compared to the asymptotic approximation and statistical thresholds are checked for conservativeness. It turns out, that corrections of the likelihood ratios in such finite-sample applications are required in order to avoid anti-conservative results.

The number of published metagenome assemblies is rapidly growing due to advances in sequencing technologies. However, sequencing errors, variable coverage, repetitive genomic regions, and other factors can produce misassemblies, which are challenging to detect for taxonomically novel genomic data. Assembly errors can affect all downstream analyses of the assemblies. Accuracy for the state of the art in reference-free misassembly prediction does not exceed an AUPRC of 0.57, and it is not clear how well these models generalize to real-world data. Here, we present the Residual neural network for Misassembled Contig identification (ResMiCo), a deep learning approach for reference-free identification of misassembled contigs. To develop ResMiCo, we first generated a training dataset of unprecedented size and complexity that can be used for further benchmarking and developments in the field. Through rigorous validation, we show that ResMiCo is substantially more accurate than the state of the art, and the model is robust to novel taxonomic diversity and varying assembly methods. ResMiCo estimated 7% misassembled contigs per metagenome across multiple real-world datasets. We demonstrate how ResMiCo can be used to optimize metagenome assembly hyperparameters to improve accuracy, instead of optimizing solely for contiguity. The accuracy, robustness, and ease-of-use of ResMiCo make the tool suitable for general quality control of metagenome assemblies and assembly methodology optimization.

Simulating the spread of infectious diseases in human communities is critical for predicting the trajectory of an epidemic and verifying various policies to control the devastating impacts of the outbreak. Many existing simulators are based on compartment models that divide people into a few subsets and simulate the dynamics among those subsets using hypothesized differential equations. However, these models lack the requisite granularity to study the effect of intelligent policies that influence every individual in a particular way. In this work, we introduce a simulator software capable of modeling a population structure and controlling the disease's propagation at an individualistic level. In order to estimate the confidence of the conclusions drawn from the simulator, we employ a comprehensive probabilistic approach where the entire population is constructed as a hierarchical random variable. This approach makes the inferred conclusions more robust against sampling artifacts and gives confidence bounds for decisions based on the simulation results. To showcase potential applications, the simulator parameters are set based on the formal statistics of the COVID-19 pandemic, and the outcome of a wide range of control measures is investigated. Furthermore, the simulator is used as the environment of a reinforcement learning problem to find the optimal policies to control the pandemic. The obtained experimental results indicate the simulator's adaptability and capacity in making sound predictions and a successful policy derivation example based on real-world data. As an exemplary application, our results show that the proposed policy discovery method can lead to control measures that produce significantly fewer infected individuals in the population and protect the health system against saturation.

When judging the average value of sample stimuli (e.g., numbers) people tend to either over- or underweight extreme sample values, depending on task context. In a context of overweighting, recent work has shown that extreme sample values were overly represented also in neural signals, in terms of an anti-compressed geometry of number samples in multivariate electroencephalography (EEG) patterns. Here, we asked whether neural representational geometries may also reflect a relative underweighting of extreme values (i.e., compression) which has been observed behaviorally in a great variety of tasks. We used a simple experimental manipulation (instructions to average a single-stream or to compare dual-streams of samples) to induce compression or anti-compression in behavior when participants judged rapid number sequences. Model-based representational similarity analysis (RSA) replicated the previous finding of neural anti-compression in the dual-stream task, but failed to provide evidence for neural compression in the single-stream task, despite the evidence for compression in behavior. Instead, the results indicated enhanced neural processing of extreme values in either task, regardless of whether extremes were over- or underweighted in subsequent behavioral choice. We further observed more general differences in the neural representation of the sample information between the two tasks. Together, our results indicate a mismatch between sample-level EEG geometries and behavior, which raises new questions about the origin of common psychometric distortions, such as diminishing sensitivity for larger values.

The clustering of platelet glycoprotein receptors with cytosolic YxxL and YxxM motifs, including GPVI, CLEC-2 and PEAR1, triggers activation via phosphorylation of the conserved tyrosine residues and recruitment of the tandem SH2 (Src homology 2) domain effector proteins, Syk and PI 3-kinase. We have modelled the clustering of these receptors with monovalent, divalent and tetravalent soluble ligands and with transmembrane ligands based on the law of mass action using ordinary differential equations and agent-based modelling. The models were experimentally evaluated in platelets and transfected cell lines using monovalent and multivalent ligands, including novel nanobody-based divalent and tetravalent ligands, by fluorescence correlation spectroscopy. Ligand valency, receptor number, receptor dimerisation, receptor phosphorylation and a cytosolic tandem SH2 domain protein act in synergy to drive receptor clustering. Threshold concentrations of a CLEC-2-blocking antibody and Syk inhibitor act in synergy to block platelet aggregation. This offers a strategy for countering the effect of avidity of multivalent ligands and in limiting off-target effects.

Distant metastasis-free survival (DMFS) curves are widely used in oncology. They are classically analyzed using the Kaplan-Meier estimator or agnostic statistical models from survival analysis. Here we report on a method to extract more information from DMFS curves using a mathematical model of primary tumor growth and metastatic dissemination. The model depends on two parameters, α and μ, respectively quantifying tumor growth and dissemination. We assumed these to be lognormally distributed in a patient population. We propose a method for identification of the parameters of these distributions based on least-squares minimization between the data and the simulated survival curve. We studied the practical identifiability of these parameters and found that including the percentage of patients with metastasis at diagnosis was critical to ensure robust estimation. We also studied the impact and identifiability of covariates and their coefficients in α and μ, either categorical or continuous, including various functional forms for the latter (threshold, linear or a combination of both). We found that both the functional form and the coefficients could be determined from DMFS curves. We then applied our model to a clinical dataset of metastatic relapse from kidney cancer with individual data of 105 patients. We show that the model was able to describe the data and illustrate our method to disentangle the impact of three covariates on DMFS: a categorical one (Führman grade) and two continuous ones (gene expressions of the macrophage mannose receptor 1 (MMR) and the G Protein-Coupled Receptor Class C Group 5 Member A (GPRC5a) gene). We found that all had an influence in metastasis dissemination (μ), but not on growth (α).

Identifying plausible mechanisms for compartmentalization and accumulation of the organic intermediates of early metabolic cycles in primitive cells has been a major challenge in theories of life's origins. Here, we propose a mechanism, where positive membrane potentials elevate the concentration of the organic intermediates. Positive membrane potentials are generated by positively charged surfaces of protocell membranes due to accumulation of transition metals. We find that (i) positive membrane potentials comparable in magnitude to those of modern cells can increase the concentration of the organic intermediates by several orders of magnitude; (ii) generation of large membrane potentials destabilize ion distributions; (iii) violation of electroneutrality is necessary to induce nonzero membrane potentials; and (iv) violation of electroneutrality enhances osmotic pressure and diminishes reaction efficiency, resulting in an evolutionary driving force for the formation of lipid membranes, specialized ion channels, and active transport systems.

Sustaining a robust metabolic network requires a balanced and fully functioning proteome. In addition to amino acids, many enzymes require cofactors (coenzymes and engrafted prosthetic groups) to function properly. Extensively validated resource allocation models, such as genome-scale models of metabolism and gene expression (ME-models), have the ability to compute an optimal proteome composition underlying a metabolic phenotype, including the provision of all required cofactors. Here we apply the ME-model for Escherichia coli K-12 MG1655 to computationally examine how environmental conditions change the proteome and its accompanying cofactor usage. We found that: (1) The cofactor requirements computed by the ME-model mostly agree with the standard biomass objective function used in models of metabolism alone (M-models); (2) ME-model computations reveal non-intuitive variability in cofactor use under different growth conditions; (3) An analysis of ME-model predicted protein use in aerobic and anaerobic conditions suggests an enrichment in the use of peroxyl scavenging acids in the proteins used to sustain aerobic growth; (4) The ME-model could describe how limitation in key protein components affect the metabolic state of E. coli. Genome-scale models have thus reached a level of sophistication where they reveal intricate properties of functional proteomes and how they support different E. coli lifestyles.

The availability of bacterial transcriptomes has dramatically increased in recent years. This data deluge could result in detailed inference of underlying regulatory networks, but the diversity of experimental platforms and protocols introduces critical biases that could hinder scalable analysis of existing data. Here, we show that the underlying structure of the E. coli transcriptome, as determined by Independent Component Analysis (ICA), is conserved across multiple independent datasets, including both RNA-seq and microarray datasets. We subsequently combined five transcriptomics datasets into a large compendium containing over 800 expression profiles and discovered that its underlying ICA-based structure was still comparable to that of the individual datasets. With this understanding, we expanded our analysis to over 3,000 E. coli expression profiles and predicted three high-impact regulons that respond to oxidative stress, anaerobiosis, and antibiotic treatment. ICA thus enables deep analysis of disparate data to uncover new insights that were not visible in the individual datasets.

The fitness landscape is a concept commonly used to describe evolution towards optimal phenotypes. It can be reduced to mechanistic detail using genome-scale models (GEMs) from systems biology. We use recently developed GEMs of Metabolism and protein Expression (ME-models) to study the distribution of Escherichia coli phenotypes on the rate-yield plane. We found that the measured phenotypes distribute non-uniformly to form a highly stratified fitness landscape. Systems analysis of the ME-model simulations suggest that this stratification results from discrete ATP generation strategies. Accordingly, we define "aero-types", a phenotypic trait that characterizes how a balanced proteome can achieve a given growth rate by modulating 1) the relative utilization of oxidative phosphorylation, glycolysis, and fermentation pathways; and 2) the differential employment of electron-transport-chain enzymes. This global, quantitative, and mechanistic systems biology interpretation of fitness landscape formed upon proteome allocation offers a fundamental understanding of bacterial physiology and evolution dynamics.

Mathematical models of metabolic networks utilize simulation to study system-level mechanisms and functions. Various approaches have been used to model the steady state behavior of metabolic networks using genome-scale reconstructions, but formulating dynamic models from such reconstructions continues to be a key challenge. Here, we present the Mass Action Stoichiometric Simulation Python (MASSpy) package, an open-source computational framework for dynamic modeling of metabolism. MASSpy utilizes mass action kinetics and detailed chemical mechanisms to build dynamic models of complex biological processes. MASSpy adds dynamic modeling tools to the COnstraint-Based Reconstruction and Analysis Python (COBRApy) package to provide an unified framework for constraint-based and kinetic modeling of metabolic networks. MASSpy supports high-performance dynamic simulation through its implementation of libRoadRunner: the Systems Biology Markup Language (SBML) simulation engine. Three examples are provided to demonstrate how to use MASSpy: (1) a validation of the MASSpy modeling tool through dynamic simulation of detailed mechanisms of enzyme regulation; (2) a feature demonstration using a workflow for generating ensemble of kinetic models using Monte Carlo sampling to approximate missing numerical values of parameters and to quantify biological uncertainty, and (3) a case study in which MASSpy is utilized to overcome issues that arise when integrating experimental data with the computation of functional states of detailed biological mechanisms. MASSpy represents a powerful tool to address challenges that arise in dynamic modeling of metabolic networks, both at small and large scales.

The evolution of antimicrobial resistance (AMR) poses a persistent threat to global public health. Sequencing efforts have already yielded genome sequences for thousands of resistant microbial isolates and require robust computational tools to systematically elucidate the genetic basis for AMR. Here, we present a generalizable machine learning workflow for identifying genetic features driving AMR based on constructing reference strain-agnostic pan-genomes and training random subspace ensembles (RSEs). This workflow was applied to the resistance profiles of 14 antimicrobials across three urgent threat pathogens encompassing 288 Staphylococcus aureus, 456 Pseudomonas aeruginosa, and 1588 Escherichia coli genomes. We find that feature selection by RSE detects known AMR associations more reliably than common statistical tests and previous ensemble approaches, identifying a total of 45 known AMR-conferring genes and alleles across the three organisms, as well as 25 candidate associations backed by domain-level annotations. Furthermore, we find that results from the RSE approach are consistent with existing understanding of fluoroquinolone (FQ) resistance due to mutations in the main drug targets, gyrA and parC, in all three organisms, and suggest the mutational landscape of those genes with respect to FQ resistance is simple. As larger datasets become available, we expect this approach to more reliably predict AMR determinants for a wider range of microbial pathogens.

Genome-wide association studies (GWAS) identify genetic variants associated with traits or diseases. GWAS never directly link variants to regulatory mechanisms. Instead, the functional annotation of variants is typically inferred by post hoc analyses. A specific class of deep learning-based methods allows for the prediction of regulatory effects per variant on several cell type-specific chromatin features. We here describe "DeepWAS", a new approach that integrates these regulatory effect predictions of single variants into a multivariate GWAS setting. Thereby, single variants associated with a trait or disease are directly coupled to their impact on a chromatin feature in a cell type. Up to 61 regulatory SNPs, called dSNPs, were associated with multiple sclerosis (MS, 4,888 cases and 10,395 controls), major depressive disorder (MDD, 1,475 cases and 2,144 controls), and height (5,974 individuals). These variants were mainly non-coding and reached at least nominal significance in classical GWAS. The prediction accuracy was higher for DeepWAS than for classical GWAS models for 91% of the genome-wide significant, MS-specific dSNPs. DSNPs were enriched in public or cohort-matched expression and methylation quantitative trait loci and we demonstrated the potential of DeepWAS to generate testable functional hypotheses based on genotype data alone. DeepWAS is available at https://github.com/cellmapslab/DeepWAS.

Next-generation sequencing based base-by-base distance measures have become an integral complement to epidemiological investigation of infectious disease outbreaks. This study introduces PANPASCO, a computational pan-genome mapping based, pairwise distance method that is highly sensitive to differences between cases, even when located in regions of lineage specific reference genomes. We show that our approach is superior to previously published methods in several datasets and across different Mycobacterium tuberculosis lineages, as its characteristics allow the comparison of a high number of diverse samples in one analysis-a scenario that becomes more and more likely with the increased usage of whole-genome sequencing in transmission surveillance.

Response to acid stress is critical for Escherichia coli to successfully complete its life-cycle by passing through the stomach to colonize the digestive tract. To develop a fundamental understanding of this response, we established a molecular mechanistic description of acid stress mitigation responses in E. coli and integrated them with a genome-scale model of its metabolism and macromolecular expression (ME-model). We considered three known mechanisms of acid stress mitigation: 1) change in membrane lipid fatty acid composition, 2) change in periplasmic protein stability over external pH and periplasmic chaperone protection mechanisms, and 3) change in the activities of membrane proteins. After integrating these mechanisms into an established ME-model, we could simulate their responses in the context of other cellular processes. We validated these simulations using RNA sequencing data obtained from five E. coli strains grown under external pH ranging from 5.5 to 7.0. We found: i) that for the differentially expressed genes accounted for in the ME-model, 80% of the upregulated genes were correctly predicted by the ME-model, and ii) that these genes are mainly involved in translation processes (45% of genes), membrane proteins and related processes (18% of genes), amino acid metabolism (12% of genes), and cofactor and prosthetic group biosynthesis (8% of genes). We also demonstrated several intervention strategies on acid tolerance that can be simulated by the ME-model. We thus established a quantitative framework that describes, on a genome-scale, the acid stress mitigation response of E. coli that has both scientific and practical uses.

Pseudomonas aeruginosa, a main cause of human infection, can gain resistance to the antibiotic aztreonam through a mutation in NalD, a transcriptional repressor of cellular efflux. Here we combine computational analysis of clinical isolates, transcriptomics, metabolic modeling and experimental validation to find a strong association between NalD mutations and resistance to aztreonam-as well as resistance to other antibiotics-across P. aeruginosa isolated from different patients. A detailed analysis of one patient's timeline shows how this mutation can emerge in vivo and drive rapid evolution of resistance while the patient received cancer treatment, a bone marrow transplantation, and antibiotics up to the point of causing the patient's death. Transcriptomics analysis confirmed the primary mechanism of NalD action-a loss-of-function mutation that caused constitutive overexpression of the MexAB-OprM efflux system-which lead to aztreonam resistance but, surprisingly, had no fitness cost in the absence of the antibiotic. We constrained a genome-scale metabolic model using the transcriptomics data to investigate changes beyond the primary mechanism of resistance, including adaptations in major metabolic pathways and membrane transport concurrent with aztreonam resistance, which may explain the lack of a fitness cost. We propose that metabolic adaptations may allow resistance mutations to endure in the absence of antibiotics and could be targeted by future therapies against antibiotic resistant pathogens.

T cells use their T-cell receptors (TCRs) to scan other cells for antigenic peptides presented by MHC molecules (pMHC). If a TCR encounters a pMHC, it can trigger a signalling pathway that could lead to the activation of the T cell and the initiation of an immune response. It is currently not clear how the binding of pMHC to the TCR initiates signalling within the T cell. One hypothesis is that conformational changes in the TCR lead to further downstream signalling. Here we investigate four different TCRs in their free state as well as in their pMHC bound state using large scale molecular simulations totalling 26 000 ns. We find that the dynamical features within TCRs differ significantly between unbound TCR and TCR/pMHC simulations. However, apart from expected results such as reduced solvent accessibility and flexibility of the interface residues, these features are not conserved among different TCR types. The presence of a pMHC alone is not sufficient to cause cross-TCR-conserved dynamical features within a TCR. Our results argue against models of TCR triggering involving conserved allosteric conformational changes.

Horizontal gene transfer (HGT) has changed the way we regard evolution. Instead of waiting for the next generation to establish new traits, especially bacteria are able to take a shortcut via HGT that enables them to pass on genes from one individual to another, even across species boundaries. The tool Daisy offers the first HGT detection approach based on read mapping that provides complementary evidence compared to existing methods. However, Daisy relies on the acceptor and donor organism involved in the HGT being known. We introduce DaisyGPS, a mapping-based pipeline that is able to identify acceptor and donor reference candidates of an HGT event based on sequencing reads. Acceptor and donor identification is akin to species identification in metagenomic samples based on sequencing reads, a problem addressed by metagenomic profiling tools. However, acceptor and donor references have certain properties such that these methods cannot be directly applied. DaisyGPS uses MicrobeGPS, a metagenomic profiling tool tailored towards estimating the genomic distance between organisms in the sample and the reference database. We enhance the underlying scoring system of MicrobeGPS to account for the sequence patterns in terms of mapping coverage of an acceptor and donor involved in an HGT event, and report a ranked list of reference candidates. These candidates can then be further evaluated by tools like Daisy to establish HGT regions. We successfully validated our approach on both simulated and real data, and show its benefits in an investigation of an outbreak involving Methicillin-resistant Staphylococcus aureus data.

Growth rate and yield are fundamental features of microbial growth. However, we lack a mechanistic and quantitative understanding of the rate-yield relationship. Studies pairing computational predictions with experiments have shown the importance of maintenance energy and proteome allocation in explaining rate-yield tradeoffs and overflow metabolism. Recently, adaptive evolution experiments of Escherichia coli reveal a phenotypic diversity beyond what has been explained using simple models of growth rate versus yield. Here, we identify a two-dimensional rate-yield tradeoff in adapted E. coli strains where the dimensions are (A) a tradeoff between growth rate and yield and (B) a tradeoff between substrate (glucose) uptake rate and growth yield. We employ a multi-scale modeling approach, combining a previously reported coarse-grained small-scale proteome allocation model with a fine-grained genome-scale model of metabolism and gene expression (ME-model), to develop a quantitative description of the full rate-yield relationship for E. coli K-12 MG1655. The multi-scale analysis resolves the complexity of ME-model which hindered its practical use in proteome complexity analysis, and provides a mechanistic explanation of the two-dimensional tradeoff. Further, the analysis identifies modifications to the P/O ratio and the flux allocation between glycolysis and pentose phosphate pathway (PPP) as potential mechanisms that enable the tradeoff between glucose uptake rate and growth yield. Thus, the rate-yield tradeoffs that govern microbial adaptation to new environments are more complex than previously reported, and they can be understood in mechanistic detail using a multi-scale modeling approach.

Genome-scale metabolic models (GEMs) are mathematically structured knowledge bases of metabolism that provide phenotypic predictions from genomic information. GEM-guided predictions of growth phenotypes rely on the accurate definition of a biomass objective function (BOF) that is designed to include key cellular biomass components such as the major macromolecules (DNA, RNA, proteins), lipids, coenzymes, inorganic ions and species-specific components. Despite its importance, no standardized computational platform is currently available to generate species-specific biomass objective functions in a data-driven, unbiased fashion. To fill this gap in the metabolic modeling software ecosystem, we implemented BOFdat, a Python package for the definition of a Biomass Objective Function from experimental data. BOFdat has a modular implementation that divides the BOF definition process into three independent modules defined here as steps: 1) the coefficients for major macromolecules are calculated, 2) coenzymes and inorganic ions are identified and their stoichiometric coefficients estimated, 3) the remaining species-specific metabolic biomass precursors are algorithmically extracted in an unbiased way from experimental data. We used BOFdat to reconstruct the BOF of the Escherichia coli model iML1515, a gold standard in the field. The BOF generated by BOFdat resulted in the most concordant biomass composition, growth rate, and gene essentiality prediction accuracy when compared to other methods. Installation instructions for BOFdat are available in the documentation and the source code is available on GitHub (https://github.com/jclachance/BOFdat).

S. aureus is classified as a serious threat pathogen and is a priority that guides the discovery and development of new antibiotics. Despite growing knowledge of S. aureus metabolic capabilities, our understanding of its systems-level responses to different media types remains incomplete. Here, we develop a manually reconstructed genome-scale model (GEM-PRO) of metabolism with 3D protein structures for S. aureus USA300 str. JE2 containing 854 genes, 1,440 reactions, 1,327 metabolites and 673 3-dimensional protein structures. Computations were in 85% agreement with gene essentiality data from random barcode transposon site sequencing (RB-TnSeq) and 68% agreement with experimental physiological data. Comparisons of computational predictions with experimental observations highlight: 1) cases of non-essential biomass precursors; 2) metabolic genes subject to transcriptional regulation involved in Staphyloxanthin biosynthesis; 3) the essentiality of purine and amino acid biosynthesis in synthetic physiological media; and 4) a switch to aerobic fermentation upon exposure to extracellular glucose elucidated as a result of integrating time-course of quantitative exo-metabolomics data. An up-to-date GEM-PRO thus serves as a knowledge-based platform to elucidate S. aureus' metabolic response to its environment.

Allosteric regulation has traditionally been described by mathematically-complex allosteric rate laws in the form of ratios of polynomials derived from the application of simplifying kinetic assumptions. Alternatively, an approach that explicitly describes all known ligand-binding events requires no simplifying assumptions while allowing for the computation of enzymatic states. Here, we employ such a modeling approach to examine the "catalytic potential" of an enzyme-an enzyme's capacity to catalyze a biochemical reaction. The catalytic potential is the fundamental result of multiple ligand-binding events that represents a "tug of war" among the various regulators and substrates within the network. This formalism allows for the assessment of interacting allosteric enzymes and development of a network-level understanding of regulation. We first define the catalytic potential and use it to characterize the response of three key kinases (hexokinase, phosphofructokinase, and pyruvate kinase) in human red blood cell glycolysis to perturbations in ATP utilization. Next, we examine the sensitivity of the catalytic potential by using existing personalized models, finding that the catalytic potential allows for the identification of subtle but important differences in how individuals respond to such perturbations. Finally, we explore how the catalytic potential can help to elucidate how enzymes work in tandem to maintain a homeostatic state. Taken together, this work provides an interpretation and visualization of the dynamic interactions and network-level effects of interacting allosteric enzymes.

Genome-scale models of metabolism and macromolecular expression (ME-models) explicitly compute the optimal proteome composition of a growing cell. ME-models expand upon the well-established genome-scale models of metabolism (M-models), and they enable a new fundamental understanding of cellular growth. ME-models have increased predictive capabilities and accuracy due to their inclusion of the biosynthetic costs for the machinery of life, but they come with a significant increase in model size and complexity. This challenge results in models which are both difficult to compute and challenging to understand conceptually. As a result, ME-models exist for only two organisms (Escherichia coli and Thermotoga maritima) and are still used by relatively few researchers. To address these challenges, we have developed a new software framework called COBRAme for building and simulating ME-models. It is coded in Python and built on COBRApy, a popular platform for using M-models. COBRAme streamlines computation and analysis of ME-models. It provides tools to simplify constructing and editing ME-models to enable ME-model reconstructions for new organisms. We used COBRAme to reconstruct a condensed E. coli ME-model called iJL1678b-ME. This reformulated model gives functionally identical solutions to previous E. coli ME-models while using 1/6 the number of free variables and solving in less than 10 minutes, a marked improvement over the 6 hour solve time of previous ME-model formulations. Errors in previous ME-models were also corrected leading to 52 additional genes that must be expressed in iJL1678b-ME to grow aerobically in glucose minimal in silico media. This manuscript outlines the architecture of COBRAme and demonstrates how ME-models can be created, modified, and shared most efficiently using the new software framework.

Several antimalarial drugs exist, but differences between life cycle stages among malaria species pose challenges for developing more effective therapies. To understand the diversity among stages and species, we reconstructed genome-scale metabolic models (GeMMs) of metabolism for five life cycle stages and five species of Plasmodium spanning the blood, transmission, and mosquito stages. The stage-specific models of Plasmodium falciparum uncovered stage-dependent changes in central carbon metabolism and predicted potential targets that could affect several life cycle stages. The species-specific models further highlight differences between experimental animal models and the human-infecting species. Comparisons between human- and rodent-infecting species revealed differences in thiamine (vitamin B1), choline, and pantothenate (vitamin B5) metabolism. Thus, we show that genome-scale analysis of multiple stages and species of Plasmodium can prioritize potential drug targets that could be both anti-malarials and transmission blocking agents, in addition to guiding translation from non-human experimental disease models.