Background

Constraint-based modeling uses mass balances, flux capacity, and reaction directionality constraints to predict fluxes through metabolism. Although transcriptional regulation and thermodynamic constraints have been integrated into constraint-based modeling, kinetic rate laws have not been extensively used.

Results

In this study, an in vivo kinetic parameter estimation problem was formulated and solved using multi-omic data sets for Escherichia coli. To narrow the confidence intervals for kinetic parameters, a series of kinetic model simplifications were made, resulting in fewer kinetic parameters than the full kinetic model. These new parameter values are able to account for flux and concentration data from 20 different experimental conditions used in our training dataset. Concentration estimates from the simplified kinetic model were within one standard deviation for 92.7% of the 790 experimental measurements in the training set. Gibbs free energy changes of reaction were calculated to identify reactions that were often operating close to or far from equilibrium. In addition, enzymes whose activities were positively or negatively influenced by metabolite concentrations were also identified. The kinetic model was then used to calculate the maximum and minimum possible flux values for individual reactions from independent metabolite and enzyme concentration data that were not used to estimate parameter values. Incorporating these kinetically-derived flux limits into the constraint-based metabolic model improved predictions for uptake and secretion rates and intracellular fluxes in constraint-based models of central metabolism.

Conclusions

This study has produced a method for in vivo kinetic parameter estimation and identified strategies and outcomes of kinetic model simplification. We also have illustrated how kinetic constraints can be used to improve constraint-based model predictions for intracellular fluxes and biomass yield and identify potential metabolic limitations through the integrated analysis of multi-omics datasets.

Copper (Cu) is an important enzyme co-factor that is also extremely toxic at high intracellular concentrations, making active efflux mechanisms essential for preventing Cu accumulation. Here, we have investigated the mechanistic role of metallochaperones in regulating Cu efflux. We have constructed a computational model of Cu trafficking and efflux based on systems analysis of the Cu stress response of Halobacterium salinarum. We have validated several model predictions via assays of transcriptional dynamics and intracellular Cu levels, discovering a completely novel function for metallochaperones. We demonstrate that in addition to trafficking Cu ions, metallochaperones also function as buffers to modulate the transcriptional responsiveness and efficacy of Cu efflux. This buffering function of metallochaperones ultimately sets the upper limit for intracellular Cu levels and provides a mechanistic explanation for previously observed Cu metallochaperone mutation phenotypes.

Author Summary

Copper (Cu) toxicity is a problem of medical, agricultural, and environmental significance. Cu toxicity severely inhibits growth of plant roots significantly affecting their morphology; Cu overload also accounts for some of the most common metal-metabolism abnormalities and neuropsychiatric problems including Wilson's and Menkes diseases. There is a large body of literature on how Cu enters and exits the cell; the kinetic and structural details of Cu translocation between trafficking, sensing, metabolic, and pumping proteins; and phenotypes associated with defects in metalloregulatory and efflux functions. Although the role of metallochaperones in Cu-cytotoxicity has been poorly studied, it has been observed that in animals deletion of metallochaperones results in elevated intracellular Cu levels along with overexpression of the P1-type ATPase efflux pump, ultimately causing malformation with high mortality. These observations are mechanistically explained by a predictive model of the Cu circuit in Halobacterium salinarum, which serves as an excellent model system for Cu trafficking and regulation in organisms with multiple chaperones. Constructed through iterative modeling and experimentation, this model accurately recapitulates known dynamical properties of the Cu circuit and predicts that intracellular Cu-buffering emerges as a consequence of the interplay of paralogous metallochaperones that traffic and allocate Cu to distinct targets.

Rationale

Activation state-dependent secretion of alpha-1 proteinase inhibitor (A1PI) by monocytes and macrophages was first reported in 1985. Since then, monocytes and tissue macrophages have emerged as key sentinels of infection and tissue damage via activation of self-assembling pattern recognition receptors (inflammasomes), which trigger inflammation and cell death in a caspase-1 dependent process. These studies examine the relationship between A1PI expression in primary monocytes and monocytic cell lines, and inflammatory cytokine expression in response to inflammasome directed stimuli.

Methods

IL-1 β expression was examined in lung macrophages expressing wild type A1PI (A1PI-M) or disease-associated Z isoform A1PI (A1PI-Z). Inflammatory cytokine release was evaluated in THP-1 monocytic cells or THP-1 cells lacking the inflammasome adaptor ASC, transfected with expression vectors encoding A1PI-M or A1PI-Z. A1PI-M was localized within monocytes by immunoprecipitation in hypotonic cell fractions. Cell-free titration of A1PI-M was performed against recombinant active caspase-1 in vitro.

Results

IL-1 β expression was elevated in lung macrophages expressing A1PI-Z. Overexpression of A1PI-M in THP-1 monocytes reduced secretion of IL-1β and TNF-α. In contrast, overexpression of A1PI-Z enhanced IL-1β and TNF- α secretion in an ASC dependent manner. A1PI-Z-enhanced cytokine release was inhibited by a small molecule caspase-1 inhibitor but not by high levels of exogenous wtA1PI. Cytosolic localization of A1PI-M in monocytes was not diminished with microtubule-inhibiting agents. A1PI-M co-localized with caspase-1 in gel-filtered cytoplasmic THP-1 preparations, and was co-immunoprecipitated with caspase 1 from nigericin-stimulated THP-1 cell lysate. Plasma-derived A1PI inhibited recombinant caspase-1 mediated conversion of a peptide substrate in a dose dependent manner.

Conclusions

Our results suggest that monocyte/macrophage-expressed A1PI-M antagonizes IL-1β secretion possibly via caspase-1 inhibition, a function which disease-associated A1PI-Z may lack. Therapeutic approaches which limit inflammasome responses in patients with A1PI deficiency, in combination with A1PI augmentation, may provide additional respiratory tissue-sparing benefits.

Regulation of the actin cytoskeleton is essential for epithelial cell polarity and protein trafficking within human uterine epithelium. The actin-binding protein cofilin is involved in regulation of actin dynamics by promoting actin branching and cytoskeleton reorganization. Dual immunohistochemical staining of cofilin and G-actin (represented by DNAse I staining) revealed cofilin-G-actin colocalization in the apical side of luminal epithelial cells of human uterine endometrium during the proliferative phase of the menstrual cycle. Interestingly, during the secretory phase of the menstrual cycle, cofilin was only present on the basolateral side. To determine whether the disease endometriosis causes a different pattern of actin remodeling, we investigated an established baboon model of induced endometriosis. The cofilin pattern in the secretory phase of baboons with endometriosis was similar to the proliferative phase in normal animals; cofilin was observed in the apical parts of luminal and glandular epithelium. A phosphatase regulating the activity of cofilin, slingshot (SSH1), revealed a similar staining pattern within these tissues. These patterns were confirmed through quantitative image analysis. Quantification of messenger RNA (mRNA) detected upregulated SSH1 and suggested a progesterone resistance-related pattern of nuclear steroid hormone receptors, but no change in membrane progesterone receptors (mPR alpha or mPR beta) was observed in endometriosis. Our data indicate that the severe dyssynchrony during menstrual cycle phases in endometriosis is connected with improper cytoskeleton rearrangements. We suggest that cofilin-mediated actin reorganization in uterine epithelial cells might be important in preparation for blastocyst implantation; dysregulation of this reorganization may lead to decreased fertility in endometriosis.

Regulation of the actin cytoskeleton is essential for epithelial cell polarity and protein trafficking within human uterine epithelium. The actin-binding protein cofilin is involved in regulation of actin dynamics by promoting actin branching and cytoskeleton reorganization. Dual immunohistochemical staining of cofilin and G-actin (represented by DNAse I staining) revealed cofilin-G-actin colocalization in the apical side of luminal epithelial cells of human uterine endometrium during the proliferative phase of the menstrual cycle. Interestingly, during the secretory phase of the menstrual cycle cofilin was only present on the basolateral side. To determine whether the disease endometriosis causes a different pattern of actin remodeling, we investigated an established baboon model of induced endometriosis. The cofilin pattern in the secretory phase of baboons with endometriosis was similar to the proliferative phase in normal animals; cofilin was observed in the apical parts of luminal and glandular epithelium. A phosphatase regulating the activity of cofilin, slingshot (SSH1), revealed a similar staining pattern within these tissues. These patterns were confirmed through quantitative image analysis. Quantification of mRNA detected upregulated SSH1 and suggested a progesterone resistance related pattern of nuclear steroid hormone receptors, but no change in membrane progesterone receptors (mPR alpha or mPR beta) was observed in endometriosis. Our data indicate that the severe dyssynchrony during menstrual cycle phases in endometriosis is connected with improper cytoskeleton rearrangements. We suggest that cofilin-mediated actin reorganization in uterine epithelial cells might be important in preparation for blastocyst implantation; dysregulation of this reorganization may lead to decreased fertility in endometriosis.

Predicting cellular responses to perturbations is an important task in systems biology. We report a new approach, RELATCH, which uses flux and gene expression data from a reference state to predict metabolic responses in a genetically or environmentally perturbed state. Using the concept of relative optimality, which considers relative flux changes from a reference state, we hypothesize a relative metabolic flux pattern is maintained from one state to another, and that cells adapt to perturbations using metabolic and regulatory reprogramming to preserve this relative flux pattern. This constraint-based approach will have broad utility where predictions of metabolic responses are needed.

Background

A33 is a type II integral membrane protein expressed on the extracellular enveloped form of vaccinia virus (VACV). Passive transfer of A33-directed monoclonal antibodies or vaccination with an A33 subunit vaccine confers protection against lethal poxvirus challenge in animal models. Homologs of A33 are highly conserved among members of the Orthopoxvirus genus and are potential candidates for inclusion in vaccines or assays targeting extracellular enveloped virus activity. One monoclonal antibody directed against VACV A33, MAb-1G10, has been shown to target a conformation-dependent epitope. Interestingly, while it recognizes VACV A33 as well as the corresponding variola homolog, it does not bind to the monkeypox homolog. In this study, we utilized a random phage display library to investigate the epitope recognized by MAb-1G10 that is critical for facilitating cell-to-cell spread of the vaccinia virus.

Results

By screening with linear or conformational random phage libraries, we found that phages binding to MAb-1G10 display the consensus motif CEPLC, with a disulfide bond formed between two cysteine residues required for MAb-1G10 binding. Although the phage motif contained no linear sequences homologous to VACV A33, structure modeling and analysis suggested that residue D115 is important to form the minimal epitope core. A panel of point mutants expressing the ectodomain of A33 protein was generated and analyzed by either binding assays such as ELISA and immunoprecipitation or a functional assessment by blocking MAb-1G10 mediated comet inhibition in cell culture.

Conclusions

These results confirm L118 as a component of the MAb-1G10 binding epitope, and further identify D115 as an essential residue. By defining the minimum conformational structure, as well as the conformational arrangement of a short peptide sequence recognized by MAb-1G10, these results introduce the possibility of designing small molecule mimetics that may interfere with the function of A33 in vivo. This information will also be useful for designing improved assays to evaluate the potency of monoclonal and polyclonal products that target A33 or A33-modulated EV dissemination.

Constraint-based models of metabolism have been used in a variety of studies on drug discovery, metabolic engineering, evolution, and multi-species interactions. These genome-scale models can be generated for any sequenced organism since their main parameters (i.e., reaction stoichiometry) are highly conserved. Their relatively low parameter requirement makes these models easy to develop; however, these models often result in a solution space with multiple possible flux distributions, making it difficult to determine the precise flux state in the cell. Recent research efforts in this modeling field have investigated how additional experimental data, including gene expression, protein expression, metabolite concentrations, and kinetic parameters, can be used to reduce the solution space. This mini-review provides a summary of the data-driven computational approaches that are available for reducing the solution space and thereby improve predictions of intracellular fluxes by constraint-based models.

The segmentation gene network in Drosophila embryo solves the fundamental problem of embryonic patterning: how to establish a periodic pattern of gene expression, which determines both the positions and the identities of body segments. The gap gene network constitutes the first zygotic regulatory tier in this process. Here we have applied the systems-level approach to investigate the regulatory effect of gap gene Kruppel (Kr) on segmentation gene expression. We acquired a large dataset on the expression of gap genes in Kr null mutants and demonstrated that the expression levels of these genes are significantly reduced in the second half of cycle 14A. To explain this novel biological result we applied the gene circuit method which extracts regulatory information from spatial gene expression data. Previous attempts to use this formalism to correctly and quantitatively reproduce gap gene expression in mutants for a trunk gap gene failed, therefore here we constructed a revised model and showed that it correctly reproduces the expression patterns of gap genes in Kr null mutants. We found that the remarkable alteration of gap gene expression patterns in Kr mutants can be explained by the dynamic decrease of activating effect of Cad on a target gene and exclusion of Kr gene from the complex network of gap gene interactions, that makes it possible for other interactions, in particular, between hb and gt, to come into effect. The successful modeling of the quantitative aspects of gap gene expression in mutant for the trunk gap gene Kr is a significant achievement of this work. This result also clearly indicates that the oversimplified representation of transcriptional regulation in the previous models is one of the reasons for unsuccessful attempts of mutant simulations.

Author Summary

Systems biology is aimed to develop an understanding of biological function or process as a system of interacting components. Here we apply the systems-level approach to understand how the blueprints for segments in the fruit fly Drosophila embryo arise. We obtain gene expression data and use the gene circuits method which allow us to reconstruct the segment determination process in the computer. To understand the system we need not only to describe it in detail, but also to comprehend what happens when certain stimuli or disruptions occur. Previous attempts to model segmentation gene expression patterns in a mutant for a trunk gap gene were unsuccessful. Here we describe the extension of the model that allows us to solve this problem in the context of Kruppel (Kr) gene. We show that remarkable alteration of gap gene expression patterns in Kr mutants can be explained by dynamic decrease of the activating effect of Cad on a target gene and exclusion of Kr from the complex network of gap gene interactions, that makes it possible for other interactions, in particular between hb and gt, to come into effect.

Phenotypic differences of genetically identical cells under the same environmental conditions have been attributed to the inherent stochasticity of biochemical processes. Various mechanisms have been suggested, including the existence of alternative steady states in regulatory networks that are reached by means of stochastic fluctuations, long transient excursions from a stable state to an unstable excited state, and the switching on and off of a reaction network according to the availability of a constituent chemical species. Here we analyse a detailed stochastic kinetic model of two-component system signalling in bacteria, and show that alternative phenotypes emerge in the absence of these features. We perform a bifurcation analysis of deterministic reaction rate equations derived from the model, and find that they cannot reproduce the whole range of qualitative responses to external signals demonstrated by direct stochastic simulations. In particular, the mixed mode, where stochastic switching and a graded response are seen simultaneously, is absent. However, probabilistic and equation-free analyses of the stochastic model that calculate stationary states for the mean of an ensemble of stochastic trajectories reveal that slow transcription of either response regulator or histidine kinase leads to the coexistence of an approximate basal solution and a graded response that combine to produce the mixed mode, thus establishing its essential stochastic nature. The same techniques also show that stochasticity results in the observation of an all-or-none bistable response over a much wider range of external signals than would be expected on deterministic grounds. Thus we demonstrate the application of numerical equation-free methods to a detailed biochemical reaction network model, and show that it can provide new insight into the role of stochasticity in the emergence of phenotypic diversity.

Author Summary

It is a surprising fact that genetically identical bacteria, living in identical conditions, can develop in completely different ways: for example, one subpopulation might grow very fast and another very slowly. These different phenotypes are thought to be one reason why bacteria that cause disease can survive antibiotic treatment or become persistent. This diversity of behaviour is usually attributed to the existence of multiple stable phenotypic states, or to the coexistence of one stable state with another unstable excited state, or finally to the possibility of the whole biochemical system that controls the phenotype being switched on and off. In this paper we describe a different scenario that leads to phenotypic diversity in two-component system signalling, a very common mechanism that bacteria use to sense external signals and control their response to changes in their environment. We use probability theory and equation-free computational analysis to calculate the average number of molecules of each chemical species present in the two-component system and hence show that sporadic production of either of two key chemical components required for signalling can delay the response to the external signal in some bacterial cells and so lead to the emergence of two distinct cell populations.

Shewanella oneidensis MR-1 is a facultative anaerobe that derives energy by coupling organic matter oxidation to the reduction of a wide range of electron acceptors. Here, we quantitatively assessed the lactate and pyruvate metabolism of MR-1 under three distinct conditions: electron acceptor-limited growth on lactate with O2, lactate with fumarate, and pyruvate fermentation. The latter does not support growth but provides energy for cell survival. Using physiological and genetic approaches combined with flux balance analysis, we showed that the proportion of ATP produced by substrate-level phosphorylation varied from 33% to 72.5% of that needed for growth depending on the electron acceptor nature and availability. While being indispensable for growth, the respiration of fumarate does not contribute significantly to ATP generation and likely serves to remove formate, a product of pyruvate formate-lyase-catalyzed pyruvate disproportionation. Under both tested respiratory conditions, S. oneidensis MR-1 carried out incomplete substrate oxidation, whereby the tricarboxylic acid (TCA) cycle did not contribute significantly. Pyruvate dehydrogenase was not involved in lactate metabolism under conditions of O2 limitation but was required for anaerobic growth, likely by supplying reducing equivalents for biosynthesis. The results suggest that pyruvate fermentation by S. oneidensis MR-1 cells represents a combination of substrate-level phosphorylation and respiration, where pyruvate serves as an electron donor and an electron acceptor. Pyruvate reduction to lactate at the expense of formate oxidation is catalyzed by a recently described new type of oxidative NAD(P)H-independent d-lactate dehydrogenase (Dld-II). The results further indicate that pyruvate reduction coupled to formate oxidation may be accompanied by the generation of proton motive force.

The physiology of ethanologenic Escherichia coli grown anaerobically in alkali-pretreated plant hydrolysates is complex and not well studied. To gain insight into how E. coli responds to such hydrolysates, we studied an E. coli K-12 ethanologen fermenting a hydrolysate prepared from corn stover pretreated by ammonia fiber expansion. Despite the high sugar content (∼6% glucose, 3% xylose) and relatively low toxicity of this hydrolysate, E. coli ceased growth long before glucose was depleted. Nevertheless, the cells remained metabolically active and continued conversion of glucose to ethanol until all glucose was consumed. Gene expression profiling revealed complex and changing patterns of metabolic physiology and cellular stress responses during an exponential growth phase, a transition phase, and the glycolytically active stationary phase. During the exponential and transition phases, high cell maintenance and stress response costs were mitigated, in part, by free amino acids available in the hydrolysate. However, after the majority of amino acids were depleted, the cells entered stationary phase, and ATP derived from glucose fermentation was consumed entirely by the demands of cell maintenance in the hydrolysate. Comparative gene expression profiling and metabolic modeling of the ethanologen suggested that the high energetic cost of mitigating osmotic, lignotoxin, and ethanol stress collectively limits growth, sugar utilization rates, and ethanol yields in alkali-pretreated lignocellulosic hydrolysates.

Genome-scale network reconstructions are useful tools for understanding cellular metabolism, and comparisons of such reconstructions can provide insight into metabolic differences between organisms. Recent efforts toward comparing genome-scale models have focused primarily on aligning metabolic networks at the reaction level and then looking at differences and similarities in reaction and gene content. However, these reaction comparison approaches are time-consuming and do not identify the effect network differences have on the functional states of the network. We have developed a bilevel mixed-integer programming approach, CONGA, to identify functional differences between metabolic networks by comparing network reconstructions aligned at the gene level. We first identify orthologous genes across two reconstructions and then use CONGA to identify conditions under which differences in gene content give rise to differences in metabolic capabilities. By seeking genes whose deletion in one or both models disproportionately changes flux through a selected reaction (e.g., growth or by-product secretion) in one model over another, we are able to identify structural metabolic network differences enabling unique metabolic capabilities. Using CONGA, we explore functional differences between two metabolic reconstructions of Escherichia coli and identify a set of reactions responsible for chemical production differences between the two models. We also use this approach to aid in the development of a genome-scale model of Synechococcus sp. PCC 7002. Finally, we propose potential antimicrobial targets in Mycobacterium tuberculosis and Staphylococcus aureus based on differences in their metabolic capabilities. Through these examples, we demonstrate that a gene-centric approach to comparing metabolic networks allows for a rapid comparison of metabolic models at a functional level. Using CONGA, we can identify differences in reaction and gene content which give rise to different functional predictions. Because CONGA provides a general framework, it can be applied to find functional differences across models and biological systems beyond those presented here.

Genome-scale metabolic models have proven useful for answering fundamental questions about metabolic capabilities of a variety of microorganisms, as well as informing their metabolic engineering. However, only a few models are available for oxygenic photosynthetic microorganisms, particularly in cyanobacteria in which photosynthetic and respiratory electron transport chains (ETC) share components. We addressed the complexity of cyanobacterial ETC by developing a genome-scale model for the diazotrophic cyanobacterium, Cyanothece sp. ATCC 51142. The resulting metabolic reconstruction, iCce806, consists of 806 genes associated with 667 metabolic reactions and includes a detailed representation of the ETC and a biomass equation based on experimental measurements. Both computational and experimental approaches were used to investigate light-driven metabolism in Cyanothece sp. ATCC 51142, with a particular focus on reductant production and partitioning within the ETC. The simulation results suggest that growth and metabolic flux distributions are substantially impacted by the relative amounts of light going into the individual photosystems. When growth is limited by the flux through photosystem I, terminal respiratory oxidases are predicted to be an important mechanism for removing excess reductant. Similarly, under photosystem II flux limitation, excess electron carriers must be removed via cyclic electron transport. Furthermore, in silico calculations were in good quantitative agreement with the measured growth rates whereas predictions of reaction usage were qualitatively consistent with protein and mRNA expression data, which we used to further improve the resolution of intracellular flux values.

Author Summary

Cyanobacteria have been promoted as platforms for biofuel production due to their useful physiological properties such as photosynthesis, relatively rapid growth rates, ability to accumulate high amounts of intracellular compounds and tolerance to extreme environments. However, development of a computational model is an important step to synthesize biochemical, physiological and regulatory understanding of photoautotrophic metabolism (either qualitatively or quantitatively) at a systems level, to make metabolic engineering of these organisms tractable. When integrated with other genome-scale data (e.g., expression data), numerical simulations can provide experimentally testable predictions of carbon fluxes and reductant partitioning to different biosynthetic pathways and macromolecular synthesis. This work is the first to computationally explore the interactions between components of photosynthetic and respiratory systems in detail. In silico predictions obtained from model analysis provided insights into the effects of light quantity and quality upon fluxes through electron transport pathways, alternative pathways for reductant consumption and carbon metabolism. The model will not only serve as a platform to develop genome-scale metabolic models for other cyanobacteria, but also as an engineering tool for manipulation of photosynthetic microorganisms to improve biofuel production.

Shewanella oneidensis MR-1 sequentially utilizes lactate and its waste products (pyruvate and acetate) during batch culture. To decipher MR-1 metabolism, we integrated genome-scale flux balance analysis (FBA) into a multiple-substrate Monod model to perform the dynamic flux balance analysis (dFBA). The dFBA employed a static optimization approach (SOA) by dividing the batch time into small intervals (i.e., ∼400 mini-FBAs), then the Monod model provided time-dependent inflow/outflow fluxes to constrain the mini-FBAs to profile the pseudo-steady-state fluxes in each time interval. The mini-FBAs used a dual-objective function (a weighted combination of “maximizing growth rate” and “minimizing overall flux”) to capture trade-offs between optimal growth and minimal enzyme usage. By fitting the experimental data, a bi-level optimization of dFBA revealed that the optimal weight in the dual-objective function was time-dependent: the objective function was constant in the early growth stage, while the functional weight of minimal enzyme usage increased significantly when lactate became scarce. The dFBA profiled biologically meaningful dynamic MR-1 metabolisms: 1. the oxidative TCA cycle fluxes increased initially and then decreased in the late growth stage; 2. fluxes in the pentose phosphate pathway and gluconeogenesis were stable in the exponential growth period; and 3. the glyoxylate shunt was up-regulated when acetate became the main carbon source for MR-1 growth.

Author Summary

This study integrates two modeling approaches, a Monod kinetic model and genome-scale flux balance analysis, to analyze the dynamic metabolism of an environmentally important bacterium (S. oneidensis MR-1). The modeling results reveal that MR-1 metabolism is suboptimal for biomass growth, while MR-1 continuously reprograms the intracellular flux distributions in adaption to nutrient conditions. This innovative dFBA framework can be widely used to investigate transient cell metabolisms in response to environmental variations. Furthermore, the dFBA is able to simulate metabolite-labeling dynamics in 13C-tracer experiments, and thus can serve as a springboard to advanced 13C-assisted dynamic metabolic flux analysis by using labeled proteinogenic amino acids to improve flux results.

Background

Carbon-13 (13C) analysis is a commonly used method for estimating reaction rates in biochemical networks. The choice of carbon labeling pattern is an important consideration when designing these experiments. We present a novel Monte Carlo algorithm for finding the optimal substrate input label for a particular experimental objective (flux or flux ratio). Unlike previous work, this method does not require assumption of the flux distribution beforehand.

Results

Using a large E. coli isotopomer model, different commercially available substrate labeling patterns were tested computationally for their ability to determine reaction fluxes. The choice of optimal labeled substrate was found to be dependent upon the desired experimental objective. Many commercially available labels are predicted to be outperformed by complex labeling patterns. Based on Monte Carlo Sampling, the dimensionality of experimental data was found to be considerably less than anticipated, suggesting that effectiveness of 13C experiments for determining reaction fluxes across a large-scale metabolic network is less than previously believed.

Conclusions

While 13C analysis is a useful tool in systems biology, high redundancy in measurements limits the information that can be obtained from each experiment. It is however possible to compute potential limitations before an experiment is run and predict whether, and to what degree, the rate of each reaction can be resolved.

Background

Despite the availability of numerous complete genome sequences from E. coli strains, published genome-scale metabolic models exist only for two commensal E. coli strains. These models have proven useful for many applications, such as engineering strains for desired product formation, and we sought to explore how constructing and evaluating additional metabolic models for E. coli strains could enhance these efforts.

Results

We used the genomic information from 16 E. coli strains to generate an E. coli pangenome metabolic network by evaluating their collective 76,990 ORFs. Each of these ORFs was assigned to one of 17,647 ortholog groups including ORFs associated with reactions in the most recent metabolic model for E. coli K-12. For orthologous groups that contain an ORF already represented in the MG1655 model, the gene to protein to reaction associations represented in this model could then be easily propagated to other E. coli strain models. All remaining orthologous groups were evaluated to see if new metabolic reactions could be added to generate a pangenome-scale metabolic model (iEco1712_pan). The pangenome model included reactions from a metabolic model update for E. coli K-12 MG1655 (iEco1339_MG1655) and enabled development of five additional strain-specific genome-scale metabolic models. These additional models include a second K-12 strain (iEco1335_W3110) and four pathogenic strains (two enterohemorrhagic E. coli O157:H7 and two uropathogens). When compared to the E. coli K-12 models, the metabolic models for the enterohemorrhagic (iEco1344_EDL933 and iEco1345_Sakai) and uropathogenic strains (iEco1288_CFT073 and iEco1301_UTI89) contained numerous lineage-specific gene and reaction differences. All six E. coli models were evaluated by comparing model predictions to carbon source utilization measurements under aerobic and anaerobic conditions, and to batch growth profiles in minimal media with 0.2% (w/v) glucose. An ancestral genome-scale metabolic model based on conserved ortholog groups in all 16 E. coli genomes was also constructed, reflecting the conserved ancestral core of E. coli metabolism (iEco1053_core). Comparative analysis of all six strain-specific E. coli models revealed that some of the pathogenic E. coli strains possess reactions in their metabolic networks enabling higher biomass yields on glucose. Finally the lineage-specific metabolic traits were compared to the ancestral core model predictions to derive new insight into the evolution of metabolism within this species.

Conclusion

Our findings demonstrate that a pangenome-scale metabolic model can be used to rapidly construct additional E. coli strain-specific models, and that quantitative models of different strains of E. coli can accurately predict strain-specific phenotypes. Such pangenome and strain-specific models can be further used to engineer metabolic phenotypes of interest, such as designing new industrial E. coli strains.

Recent evidence from serum metabolomics indicates that specific metabolic disturbances precede β-cell autoimmunity in humans and can be used to identify those children who subsequently progress to type 1 diabetes. The mechanisms behind these disturbances are unknown. Here we show the specificity of the pre-autoimmune metabolic changes, as indicated by their conservation in a murine model of type 1 diabetes. We performed a study in non-obese prediabetic (NOD) mice which recapitulated the design of the human study and derived the metabolic states from longitudinal lipidomics data. We show that female NOD mice who later progress to autoimmune diabetes exhibit the same lipidomic pattern as prediabetic children. These metabolic changes are accompanied by enhanced glucose-stimulated insulin secretion, normoglycemia, upregulation of insulinotropic amino acids in islets, elevated plasma leptin and adiponectin, and diminished gut microbial diversity of the Clostridium leptum group. Together, the findings indicate that autoimmune diabetes is preceded by a state of increased metabolic demands on the islets resulting in elevated insulin secretion and suggest alternative metabolic related pathways as therapeutic targets to prevent diabetes.

Author Summary

We have recently found that distinct metabolic disturbances precede β-cell autoimmunity in children who later progress to type 1 diabetes (T1D). Here we performed a murine study using non-obese diabetic (NOD) mice that recapitulated the protocol used in human, followed up by independent studies where NOD mice were studied in relation to risk of diabetes progression. We found that young female NOD mice who later progress to autoimmune diabetes exhibit the same lipidomic pattern as prediabetic children. These metabolic changes are accompanied by enhanced glucose-stimulated insulin secretion, upregulation of insulinotropic amino acids in islets, elevated plasma leptin and adiponectin, and diminished gut microbial diversity of the Clostridium leptum subgroup. The metabolic phenotypes observed in our study could be relevant as end points for studies investigating T1D pathogenesis and/or responses to interventions. By proceeding from a clinical study via metabolomics and modeling to an experimental model using a similar study design, then evolving further to tissue-specific studies, we hereby also present a conceptually novel approach to reversed translation that may be useful in future therapeutic studies in the context of prevention and treatment of T1D as well as of other diseases characterized by long prodromal periods.

What governs the concentrations of metabolites within living cells? Beyond specific metabolic and enzymatic considerations, are there global trends that affect their values? We hypothesize that the physico-chemical properties of metabolites considerably affect their in-vivo concentrations. The recently achieved experimental capability to measure the concentrations of many metabolites simultaneously has made the testing of this hypothesis possible. Here, we analyze such recently available data sets of metabolite concentrations within E. coli, S. cerevisiae, B. subtilis and human. Overall, these data sets encompass more than twenty conditions, each containing dozens (28-108) of simultaneously measured metabolites. We test for correlations with various physico-chemical properties and find that the number of charged atoms, non-polar surface area, lipophilicity and solubility consistently correlate with concentration. In most data sets, a change in one of these properties elicits a ∼100 fold increase in metabolite concentrations. We find that the non-polar surface area and number of charged atoms account for almost half of the variation in concentrations in the most reliable and comprehensive data set. Analyzing specific groups of metabolites, such as amino-acids or phosphorylated nucleotides, reveals even a higher dependence of concentration on hydrophobicity. We suggest that these findings can be explained by evolutionary constraints imposed on metabolite concentrations and discuss possible selective pressures that can account for them. These include the reduction of solute leakage through the lipid membrane, avoidance of deleterious aggregates and reduction of non-specific hydrophobic binding. By highlighting the global constraints imposed on metabolic pathways, future research could shed light onto aspects of biochemical evolution and the chemical constraints that bound metabolic engineering efforts.

Author Summary

What governs the identity and concentrations of metabolites within living cells? The first part of this question has received much attention. Organisms were found to qualitatively prefer hydrophilic and charged metabolites, a phenomenon that was explained to be a result of constraints imposed by contemporary as well as archaic metabolism. However, among the metabolites that are used, a quantitative preference has never been analyzed systematically. Here we use the most comprehensive data sets of metabolite concentrations available to explore such trends. We find that in various organisms and growth conditions, living cells minimize the concentrations of non-polar, un-charged metabolites. More specifically, metabolites' hydrophobicity alters concentrations by two orders of magnitudes on average and explains up to half of the variation of metabolite concentrations within cells. We suggest that this can be attributed to an evolutionary pressure to avoid an unspecific hydrophobic effect: the preference of hydrophobic surfaces in an aqueous environment to adhere to other hydrophobic surfaces. Our findings shed light on the evolution of the internal makeup of living cells and can assist in establishing metabolic models that support synthetic biology and metabolic engineering efforts.

The use of computational models in metabolic engineering has been increasing as more genome-scale metabolic models and computational approaches become available. Various computational approaches have been developed to predict how genetic perturbations affect metabolic behavior at a systems level, and have been successfully used to engineer microbial strains with improved primary or secondary metabolite production. However, identification of metabolic engineering strategies involving a large number of perturbations is currently limited by computational resources due to the size of genome-scale models and the combinatorial nature of the problem. In this study, we present (i) two new bi-level strain design approaches using mixed-integer programming (MIP), and (ii) general solution techniques that improve the performance of MIP-based bi-level approaches. The first approach (SimOptStrain) simultaneously considers gene deletion and non-native reaction addition, while the second approach (BiMOMA) uses minimization of metabolic adjustment to predict knockout behavior in a MIP-based bi-level problem for the first time. Our general MIP solution techniques significantly reduced the CPU times needed to find optimal strategies when applied to an existing strain design approach (OptORF) (e.g., from ∼10 days to ∼5 minutes for metabolic engineering strategies with 4 gene deletions), and identified strategies for producing compounds where previous studies could not (e.g., malate and serine). Additionally, we found novel strategies using SimOptStrain with higher predicted production levels (for succinate and glycerol) than could have been found using an existing approach that considers network additions and deletions in sequential steps rather than simultaneously. Finally, using BiMOMA we found novel strategies involving large numbers of modifications (for pyruvate and glutamate), which sequential search and genetic algorithms were unable to find. The approaches and solution techniques developed here will facilitate the strain design process and extend the scope of its application to metabolic engineering.

Background

Rhodobacter sphaeroides is one of the best studied purple non-sulfur photosynthetic bacteria and serves as an excellent model for the study of photosynthesis and the metabolic capabilities of this and related facultative organisms. The ability of R. sphaeroides to produce hydrogen (H2), polyhydroxybutyrate (PHB) or other hydrocarbons, as well as its ability to utilize atmospheric carbon dioxide (CO2) as a carbon source under defined conditions, make it an excellent candidate for use in a wide variety of biotechnological applications. A genome-level understanding of its metabolic capabilities should help realize this biotechnological potential.

Results

Here we present a genome-scale metabolic network model for R. sphaeroides strain 2.4.1, designated iRsp1095, consisting of 1,095 genes, 796 metabolites and 1158 reactions, including R. sphaeroides-specific biomass reactions developed in this study. Constraint-based analysis showed that iRsp1095 agreed well with experimental observations when modeling growth under respiratory and phototrophic conditions. Genes essential for phototrophic growth were predicted by single gene deletion analysis. During pathway-level analyses of R. sphaeroides metabolism, an alternative route for CO2 assimilation was identified. Evaluation of photoheterotrophic H2 production using iRsp1095 indicated that maximal yield would be obtained from growing cells, with this predicted maximum ~50% higher than that observed experimentally from wild type cells. Competing pathways that might prevent the achievement of this theoretical maximum were identified to guide future genetic studies.

Conclusions

iRsp1095 provides a robust framework for future metabolic engineering efforts to optimize the solar- and nutrient-powered production of biofuels and other valuable products by R. sphaeroides and closely related organisms.

Transcriptional coactivator with PDZ-binding motif (TAZ) is known to bind to a variety of transcription factors to control cell differentiation and organ development. However, its role in uterine physiology has not yet been described. To study its regulation during the unique process of differentiation of fibroblasts into decidual cells (decidualization), we utilized the human uterine fibroblast (HuF) in vitro cell model. Immunocytochemistry data demonstrated that the majority of the TAZ protein is localized in the nucleus. Treatment of HuF cells with the embryonic stimulus cytokine interleukin 1 beta in the presence of steroid hormones (estradiol-17 beta and medroxyprogesterone acetate) for 13 days did not cause any apparent TAZ mRNA changes but resulted in a significant TAZ protein decline (approximately 62%) in total cell lysates. Analysis of cytosolic and nuclear extracts revealed that the decline of total TAZ was caused primarily by a drop of TAZ protein levels in the nucleus. TAZ was localized on the peroxisome proliferator-activated receptor response element site (located at position −1200 bp relative to the transcription start site) of the genomic region of decidualization marker insulin-like growth factor-binding protein 1 (IGFBP1) in HuF cells as detected by chromatin immunoprecipitation. TAZ is also present in human endometrium tissue as confirmed by immunohistochemistry. During the secretory phase of the menstrual cycle, specific TAZ staining particularly diminishes in the stroma, suggesting its participation during the decidualization process, as well as implantation. During early baboon pregnancy, TAZ protein expression remains minimal in the endometrium close to the implantation site. In summary, the presented evidence shows for the first time to date TAZ protein in the human uterine tract, its downregulation during in vitro decidualization, and its localization on the IGFBP1 promoter region, all of which indicate its presence in the uterine differentiation program during pregnancy.

Transcriptional regulator TAZ is present in the uterine differentiation program during pregnancy as evidenced by its decline during in vitro human decidualization induced by embryonic stimulus and its minimal expression at baboon implantation sites.

In biological systems that undergo processes such as differentiation, a clear concept of progression exists. We present a novel computational approach, called Sample Progression Discovery (SPD), to discover patterns of biological progression underlying microarray gene expression data. SPD assumes that individual samples of a microarray dataset are related by an unknown biological process (i.e., differentiation, development, cell cycle, disease progression), and that each sample represents one unknown point along the progression of that process. SPD aims to organize the samples in a manner that reveals the underlying progression and to simultaneously identify subsets of genes that are responsible for that progression. We demonstrate the performance of SPD on a variety of microarray datasets that were generated by sampling a biological process at different points along its progression, without providing SPD any information of the underlying process. When applied to a cell cycle time series microarray dataset, SPD was not provided any prior knowledge of samples' time order or of which genes are cell-cycle regulated, yet SPD recovered the correct time order and identified many genes that have been associated with the cell cycle. When applied to B-cell differentiation data, SPD recovered the correct order of stages of normal B-cell differentiation and the linkage between preB-ALL tumor cells with their cell origin preB. When applied to mouse embryonic stem cell differentiation data, SPD uncovered a landscape of ESC differentiation into various lineages and genes that represent both generic and lineage specific processes. When applied to a prostate cancer microarray dataset, SPD identified gene modules that reflect a progression consistent with disease stages. SPD may be best viewed as a novel tool for synthesizing biological hypotheses because it provides a likely biological progression underlying a microarray dataset and, perhaps more importantly, the candidate genes that regulate that progression.

Author Summary

We present a novel computational approach, Sample Progression Discovery (SPD), to discover biological progression underlying a microarray dataset. In contrast to the majority of microarray data analysis methods which identify differences between sample groups (normal vs. cancer, treated vs. control), SPD aims to identify an underlying progression among individual samples, both within and across sample groups. We validated SPD's ability to discover biological progression using datasets of cell cycle, B-cell differentiation, and mouse embryonic stem cell differentiation. We view SPD as a hypothesis generation tool when applied to datasets where the progression is unclear. For example, when applied to a microarray dataset of cancer samples, SPD assumes that the cancer samples collected from individual patients represent different stages during an intrinsic progression underlying cancer development. The inferred relationship among the samples may therefore indicate a trajectory or hierarchy of cancer progression, which serves as a hypothesis to be tested. SPD is not limited to microarray data analysis, and can be applied to a variety of high-dimensional datasets. We implemented SPD using MATLAB graphical user interface, which is available at http://icbp.stanford.edu/software/SPD/.

Purpose

To evaluate the impact of point-of-care (POC) testing for sexually transmitted infections (STIs) on reported awareness of test results and STI risk-reduction behaviors in adolescents.

Methods

Women ages 14–21 were recruited from the Emergency Department or Teen Health Clinic for this longitudinal study and were tested for STIs. Baseline demographics, risk behaviors, treatment, POC tests (wet mount and rapid antigen tests for Trichomonas vaginalis), and other STI test results (available 24–48 hours post-visit) were measured. These were compared to subject’s report of test results, abstinence, partner discussion, and partner testing during a post-visit telephone contact.

Results

Of 294 women, 155 (53%) were contacted: 65 (42%) had a positive STI test; 28 (43%) were POC positive; and 52 (33.5%) believed their STI results were positive. A positive POC test increased the proportion of subjects aware of being positive for an STI (89 vs. 21%, p<.01). Post-visit, 62% reported abstinence, 82% discussed testing with her partner, and 48% reported partner testing. Predictors of abstinence included a positive POC test (adjusted odds ratio (AOR) 4.6 confidence interval (CI) 1.5–13.6, prior abstinence of > 14 days (AOR 3.9, CI 1.7–9.0) and black race (AOR 3.5, CI 1.2–9.7). Women who believed their STI results were positive were more likely to report partner discussion (odds ratio (OR) 3.0, CI 1.0–8.8) and partner testing (OR 5.1, CI 2.4–11.2).

Conclusions

Awareness of STI results increases with POC testing. Effective communication of results can increase patient understanding and compliance with risk-reduction strategies which may impact the STI epidemic.