Towards kinetic modeling of genome-scale metabolic networks without sacrificing stoichiometric,thermodynamic and physiological constraints

Mathematical modeling is an essential tool for the comprehensive understanding of cell metabolism and its interactions with the environmental and process conditions. Recent developments in the construction and analysis of stoichiometric models made it possible to define limits on steady-state metabolic behavior using flux balance analysis. However, detailed information on enzyme kinetics and enzyme regulation is needed to formulate kinetic models that can accurately capture the dynamic metabolic responses. The use of mechanistic enzyme kinetics is a difficult task due to uncertainty in the kinetic properties of enzymes. Therefore, the majority of recent works considered only mass action kinetics for reactions in metabolic networks. Herein, we applied the optimization and risk analysis of complex living entities (ORACLE) framework and constructed a large-scale mechanistic kinetic model of optimally grown Escherichia coli. We investigated the complex interplay between stoichiometry, thermodynamics, and kinetics in determining the flexibility and capabilities of metabolism. Our results indicate that enzyme saturation is a necessary consideration in modeling metabolic networks and it extends the feasible ranges of metabolic fluxes and metabolite concentrations. Our results further suggest that enzymes in metabolic networks have evolved to function at different saturation states to ensure greater flexibility and robustness of cellular metabolism.     

Bringing metabolic networks to life: integration of kinetic, metabolic, and proteomic data

Background

Translating a known metabolic network into a dynamic model requires reasonable guesses of all enzyme parameters. In Bayesian parameter estimation, model parameters are described by a posterior probability distribution, which scores the potential parameter sets, showing how well each of them agrees with the data and with the prior assumptions made.

Results

We compute posterior distributions of kinetic parameters within a Bayesian framework, based on integration of kinetic, thermodynamic, metabolic, and proteomic data. The structure of the metabolic system (i.e., stoichiometries and enzyme regulation) needs to be known, and the reactions are modelled by convenience kinetics with thermodynamically independent parameters. The parameter posterior is computed in two separate steps: a first posterior summarises the available data on enzyme kinetic parameters; an improved second posterior is obtained by integrating metabolic fluxes, concentrations, and enzyme concentrations for one or more steady states. The data can be heterogenous, incomplete, and uncertain, and the posterior is approximated by a multivariate log-normal distribution. We apply the method to a model of the threonine synthesis pathway: the integration of metabolic data has little effect on the marginal posterior distributions of individual model parameters. Nevertheless, it leads to strong correlations between the parameters in the joint posterior distribution, which greatly improve the model predictions by the following Monte-Carlo simulations.

Conclusion

We present a standardised method to translate metabolic networks into dynamic models. To determine the model parameters, evidence from various experimental data is combined and weighted using Bayesian parameter estimation. The resulting posterior parameter distribution describes a statistical ensemble of parameter sets; the parameter variances and correlations can account for missing knowledge, measurement uncertainties, or biological variability. The posterior distribution can be used to sample model instances and to obtain probabilistic statements about the model's dynamic behaviour.     

Effect of rpoS gene knockout on the metabolism of Escherichia coli during exponential growth phase and early stationary phase based on gene expressions, enzyme activities and intracellular metabolite concentrations

The RNA polymerase sigma factor, encoded by rpoS gene, controls the expression of a large number of genes in Escherichia coli under stress conditions. The present study investigated the growth characteristics and metabolic pathways of rpoS gene knockout mutant of E. coli growing in LB media under aerobic condition. The analyses were made based on gene expressions obtained by DNA microarray and RT-PCR, enzyme activities and intracellular metabolite concentrations at the exponential and early stationary phases of growth. Although the glucose utilization pattern of the mutant was similar to the parent strain, the mutant failed to utilize acetate throughout the cultivation period. Microarray data indicated that the expression levels of several important genes of acetate metabolism such as acs, aceAB, cysDEK, fadR, etc. were significantly altered in the absence of rpoS gene. Interestingly, there was an increased activity of TCA cycle during the exponential growth phase, which was gradually diminished at the onset of stationary phase. Moreover, rpoS mutation had profound effect on the expression of several other genes of E. coli metabolic pathways that were not described earlier. The changes in the gene expressions, enzyme activities and intracellular metabolite concentrations of the rpoS mutant are discussed in details with reference to the major metabolic pathways of E. coli.     

Large-scale evaluation of in silico gene deletions in Saccharomyces cerevisiae

A large-scale in silico evaluation of gene deletions in Saccharomyces cerevisiae was conducted using a genome-scale reconstructed metabolic model. The effect of 599 single gene deletions on cell viability was simulated in silico and compared to published experimental results. In 526 cases (87.8%), the in silico results were in agreement with experimental observations when growth on synthetic complete medium was simulated. Viable phenotypes were correctly predicted in 89.4% (496 out of 555) and lethal phenotypes were correctly predicted in 68.2% (30 out of 44) of the cases considered. The in silico evaluation was solely based on the topological properties of the metabolic network which is based on well-established reaction stoichiometry. No interaction or regulatory information was accounted for in the in silico model. False predictions were analyzed on a case-by-case basis for four possible inadequacies of the in silico model: (1) incomplete media composition, (2) substitutable biomass components, (3) incomplete biochemical information, and (4) missing regulation. This analysis eliminated a number of false predictions and suggested a number of experimentally testable hypotheses. A genome-scale in silico model can thus be used to systematically reconcile existing data and fill in our knowledge gaps about an organism.     

An integrated analysis of the genome of the hyperthermophilic archaeon Pyrococcus abyssi

The hyperthermophilic euryarchaeon Pyrococcus abyssi and the related species Pyrococcus furiosus and Pyrococcus horikoshii, whose genomes have been completely sequenced, are presently used as model organisms in different laboratories to study archaeal DNA replication and gene expression and to develop genetic tools for hyperthermophiles. We have performed an extensive re-annotation of the genome of P. abyssi to obtain an integrated view of its phylogeny, molecular biology and physiology. Many new functions are predicted for both informational and operational proteins. Moreover, several candidate genes have been identified that might encode missing links in key metabolic pathways, some of which have unique biochemical features. The great majority of Pyrococcus proteins are typical archaeal proteins and their phylogenetic pattern agrees with its position near the root of the archaeal tree. However, proteins probably from bacterial origin, including some from mesophilic bacteria, are also present in the P. abyssi genome.     

Further insights into the isoenzyme composition and activity of glutamate dehydrogenase in Arabidopsis thaliana

Following the discovery that in Arabidopsis, a third isoenzyme of NADH-dependent glutamate dehydrogenase (GDH) is expressed in the mitochondria of the root companion cells, we have re-examined the GDH isoenzyme composition. By analyzing the NADH-GDH isoenzyme composition of single, double and triple mutants deficient in the expression of the three genes encoding the enzyme, we have found that the α, β and γ polypeptides that comprise the enzyme can be assembled into a complex combination of heterohexamers in roots. Moreover, we observed that when one or two of the three root isoenzymes were missing from the mutants, the remaining isoenzymes compensated for this deficiency. The significance of such complexity is discussed in relation to the metabolic and signaling function of the NADH-GDH enzyme. Although it has been shown that a fourth gene encoding a NADPH-dependent enzyme is present in Arabidopsis, we were not able to detect corresponding enzyme activity, even in the triple mutant totally lacking NADH-GDH activity.     

Differential gene expression profiles and real-time measurements of growth parameters in Saccharomyces cerevisiae grown in microliter-scale bioreactors equipped with internal stirring

Combining real-time growth kinetics measurements with global gene expression analysis of microbial cultures is of significant value for high-throughput biological research. We have performed differential gene expression analysis in the eukaryotic model Saccharomyces cerevisiae grown in galactose and glucose media in 150 muL bioreactors equipped with sensors for in situ and real-time measurements of optical density (OD), pH, and dissolved oxygen (DO). The microbioreactors were fabricated from poly(dimethylsiloxane) (PDMS) and poly(methyl methacrylate) (PMMA) and equipped with internal magnetic ministirrers and evaporation compensation by water replacement. In galactose-grown cells, the core genes of the GAL operon GAL2, GAL1, GAL7, and GAL10 were upregulated at least 100-fold relative to glucose-grown cells. These differential gene expression levels were similar to those observed in large-scale culture vessels. The increasing rate at which complete genomic sequences of microorganisms are becoming available offers an unprecedented opportunity for comparative investigations of these organisms. Our results from S. cerevisiae cultures grown in instrumented microbioreactors show that it is possible to integrate high-throughput studies of growth physiology with global gene expression analysis of microorganisms.     

Development of the osteoblast phenotype in primary human osteoblasts in culture: comparison with rat calvarial cells in osteoblast differentiation.

In rat osteoblast-like cells, a time-dependent sequence of growth and differentiation-dependent genes has been identified and a model of osteoblast differentiation in culture suggested. We investigated the expression of the bone matrix-associated proteins osteonectin and procollagen I and of the bone cell phenotype-related proteins alkaline phosphatase and osteocalcin during cell culture in primary human osteoblast like cells. Primary human explant cultures from nine young healthy donors were established under highly standardized conditions. Cells in the second passage were analyzed on different days from day 1 to 32, comparing cells growing under the influence of ascorbate with controls. Gene expression was determined by Northern blot analysis or polymerase chain reaction. Osteocalcin expression was also investigated after 1,25-(OH)(2)D(3) stimulation. On the protein level, newly synthesized collagen I, alkaline phosphatase activity, and secretion of osteocalcin were analyzed at all time points. On comparing our findings to the pattern of gene expression suggested for the rat calvarial osteoblast system, we found a similar developmental sequence for the so-called "proliferation" as well as a similar, but lengthened, sequence for the "matrix maturation stage." During "matrix maturation," we found an ongoing proliferation despite increased alkaline phosphatase and decreased procollagen I gene expression. Our study, therefore, shows that in pHOB the gene expression profile proceeded to the "matrix maturation stage," as defined by Owen and colleagues, independent of ongoing proliferation. We were unable to observe the mineralization period as demonstrated by the missing increase of osteocalcin expression and lack of nodule formation in our human osteoblast model. In contrast to the rat system, we found a proliferation stimulating influence of ascorbate, suggesting species-specific differences in response to differentiation factors. From these data, we conclude that general considerations on physiology and pathophysiology of bone cell differentiation have to be confirmed in the human osteoblastic cell system.     

Fluid shear-induced mechanical signaling in MC3T3-E1 osteoblasts requires cytoskeleton-integrin interactions

Mechanical stimulation of bone induces new bone formation invivo and increases the metabolic activity and gene expression ofosteoblasts in culture. We investigated the role of the actin cytoskeleton and actin-membrane interactions in the transmission ofmechanical signals leading to altered gene expression in cultured MC3T3-E1 osteoblasts. Application of fluid shear to osteoblasts causedreorganization of actin filaments into contractile stress fibers andinvolved recruitment of₁-integrins and -actinin tofocal adhesions. Fluid shear also increased expression of two proteinslinked to mechanotransduction in vivo, cyclooxygenase-2 (COX-2) and theearly response gene product c-fos. Inhibition of actin stress fiberdevelopment by treatment of cells with cytochalasin D, by expression ofa dominant negative form of the small GTPase Rho, or by microinjectioninto cells of a proteolytic fragment of -actinin that inhibits-actinin-mediated anchoring of actin filaments to integrins at theplasma membrane each blocked fluid-shear-induced gene expression inosteoblasts. We conclude that fluid shear-induced mechanical signalingin osteoblasts leads to increased expression of COX-2 and c-Fos througha mechanism that involves reorganization of the actin cytoskeleton.Thus Rho-mediated stress fiber formation and the -actinin-dependentanchorage of stress fibers to integrins in focal adhesions may promotefluid shear-induced metabolic changes in bone cells.

     

Variation in the large-scale organization of gene expression levels in the hippocampus relates to stable epigenetic variability in behavior

Background

Despite sharing the same genes, identical twins demonstrate substantial variability in behavioral traits and in their risk for disease. Epigenetic factors–DNA and chromatin modifications that affect levels of gene expression without affecting the DNA sequence–are thought to be important in establishing this variability. Epigenetically-mediated differences in the levels of gene expression that are associated with individual variability traditionally are thought to occur only in a gene-specific manner. We challenge this idea by exploring the large-scale organizational patterns of gene expression in an epigenetic model of behavioral variability.

Methodology/Findings

To study the effects of epigenetic influences on behavioral variability, we examine gene expression in genetically identical mice. Using a novel approach to microarray analysis, we show that variability in the large-scale organization of gene expression levels, rather than differences in the expression levels of specific genes, is associated with individual differences in behavior. Specifically, increased activity in the open field is associated with increased variance of log-transformed measures of gene expression in the hippocampus, a brain region involved in open field activity. Early life experience that increases adult activity in the open field also similarly modifies the variance of gene expression levels. The same association of the variance of gene expression levels with behavioral variability is found with levels of gene expression in the hippocampus of genetically heterogeneous outbred populations of mice, suggesting that variation in the large-scale organization of gene expression levels may also be relevant to phenotypic differences in outbred populations such as humans. We find that the increased variance in gene expression levels is attributable to an increasing separation of several large, log-normally distributed families of gene expression levels. We also show that the presence of these multiple log-normal distributions of gene expression levels is a universal characteristic of gene expression in eurkaryotes. We use data from the MicroArray Quality Control Project (MAQC) to demonstrate that our method is robust and that it reliably detects biological differences in the large-scale organization of gene expression levels.

Conclusions

Our results contrast with the traditional belief that epigenetic effects on gene expression occur only at the level of specific genes and suggest instead that the large-scale organization of gene expression levels provides important insights into the relationship of gene expression with behavioral variability. Understanding the epigenetic, genetic, and environmental factors that regulate the large-scale organization of gene expression levels, and how changes in this large-scale organization influences brain development and behavior will be a major future challenge in the field of behavioral genomics.     

The continuous hydrolysis of geniposide to genipin using immobilized β-glucosidase on calcium alginate gel

Summary The reaction velocity of immobilized -glucosidase was approximated by the first-order reaction kinetics. A plug flow reactor was used for continuous hydrolysis of geniposide with this immobilized enzyme. The activity of this immobilized enzyme was retained 100% for 600 h. The amount of genipin formed by using the immobilized enzyme was 17 fold that formed using the native enzyme without reuse. Using immobilized enzyme, purity and yield of genipin, which is a hydrolyzate of geniposide, was improved comparing with the native enzyme.     

Monte-Carlo Modeling of the Central Carbon Metabolism of Lactococcus lactis: Insights into Metabolic Regulation

Metabolic pathways are complex dynamic systems whose response to perturbations and environmental challenges are governed by multiple interdependencies between enzyme properties, reactions rates, and substrate levels. Understanding the dynamics arising from such a network can be greatly enhanced by the construction of a computational model that embodies the properties of the respective system. Such models aim to incorporate mechanistic details of cellular interactions to mimic the temporal behavior of the biochemical reaction system and usually require substantial knowledge of kinetic parameters to allow meaningful conclusions. Several approaches have been suggested to overcome the severe data requirements of kinetic modeling, including the use of approximative kinetics and Monte-Carlo sampling of reaction parameters. In this work, we employ a probabilistic approach to study the response of a complex metabolic system, the central metabolism of the lactic acid bacterium Lactococcus lactis, subject to perturbations and brief periods of starvation. Supplementing existing methodologies, we show that it is possible to acquire a detailed understanding of the control properties of a corresponding metabolic pathway model that is directly based on experimental observations. In particular, we delineate the role of enzymatic regulation to maintain metabolic stability and metabolic recovery after periods of starvation. It is shown that the feedforward activation of the pyruvate kinase by fructose-1,6-bisphosphate qualitatively alters the bifurcation structure of the corresponding pathway model, indicating a crucial role of enzymatic regulation to prevent metabolic collapse for low external concentrations of glucose. We argue that similar probabilistic methodologies will help our understanding of dynamic properties of small-, medium- and large-scale metabolic networks models.     

Dynamics of muscle glycogenolysis modeled with pH time course computation and pH-dependent reaction equilibria and enzyme kinetics

Cellular metabolites are moieties defined by their specific binding constants to H+, Mg2+, and K+ or anions without ligands. As a consequence, every biochemical reaction in the cytoplasm has an associated proton stoichiometry that is generally noninteger- and pH-dependent. Therefore, with metabolic flux, pH is altered in a medium with finite buffer capacity. Apparent equilibrium constants and maximum enzyme velocities, which are functions of pH, are also altered. We augmented an earlier mathematical model of skeletal muscle glycogenolysis with pH-dependent enzyme kinetics and reaction equilibria to compute the time course of pH changes. Analysis shows that kinetics and final equilibrium states of the closed system are highly constrained by the pH-dependent parameters. This kinetic model of glycogenolysis, coupled to creatine kinase and adenylate kinase, simulated published experiments made with a cell-free enzyme mixture to reconstitute the network and to synthesize PCr and lactate in vitro. Using the enzyme kinetic and thermodynamic data in the literature, the simulations required minimal adjustments of parameters to describe the data. These results show that incorporation of appropriate physical chemistry of the reactions with accurate kinetic modeling gives a reasonable simulation of experimental data and is necessary for a physically correct representation of the metabolic network. The approach is general for modeling metabolic networks beyond the specific pathway and conditions presented here.     

Prediction of temporal gene expression. Metabolic opimization by re-distribution of enzyme activities.

A computational approach is used to analyse temporal gene expression in the context of metabolic regulation. It is based on the assumption that cells developed optimal adaptation strategies to changing environmental conditions. Time-dependent enzyme profiles are calculated which optimize the function of a metabolic pathway under the constraint of limited total enzyme amount. For linear model pathways it is shown that wave-like enzyme profiles are optimal for a rapid substrate turnover. For the central metabolism of yeast cells enzyme profiles are calculated which ensure long-term homeostasis of key metabolites under conditions of a diauxic shift. These enzyme profiles are in close correlation with observed gene expression data. Our results demonstrate that optimality principles help to rationalize observed gene expression profiles.     

Dosage, deletions and dominance: simple models of the evolution of gene expression

Dominance of the wild-type allele over spontaneous null mutations, such as deletions, can be explained in terms of the effects of changes in enzyme dose on the flux of metabolic pathways. If ever increasing levels of enzyme activity have ever decreasing effects on the flux of the biochemical pathway, then halving of dosage will always have a lesser effect on flux than half the effect of complete removal of gene activity. Furthermore, if gene expression rates are high, then halving of dose can have a negligible effect on flux and dominance will be strong. Given that strong dominance appears to be common, this leaves open the issue of why enzyme activity levels are so high that a halving of expression rates is of minimal effect. Why produce so much surplus enzyme? One explanation, suggested by Haldane, is that selection favoured high expression levels as a defence against mutation. We model this scenario formally and show that protection from mutation is an extremely weak force determining expression levels. The selective coefficients are only of the order of the mutation rate. However, if we suppose a linear mapping of flux with fitness and a monotonic cost to increased gene expression, it follows simply that here exists an optimal level of gene expression. By contrast to the mutational model, doubling of gene expression rates when the system is distant from the optimum is associated with extremely high selective coefficients (orders of magnitude higher than the mutation rate). When the cost of gene expression is slight the optimal rate of expression is such that strong dominance will follow.     

Prediction of missing enzyme genes in a bacterial metabolic network. Reconstruction of the lysine-degradation pathway of Pseudomonas aeruginosa

The metabolic network is an important biological network which consists of enzymes and chemical compounds. However, a large number of metabolic pathways remains unknown, and most organism-specific metabolic pathways contain many missing enzymes. We present a novel method to identify the genes coding for missing enzymes using available genomic and chemical information from bacterial genomes. The proposed method consists of two steps: (a) estimation of the functional association between the genes with respect to chromosomal proximity and evolutionary association, using supervised network inference; and (b) selection of gene candidates for missing enzymes based on the original candidate score and the chemical reaction information encoded in the EC number. We applied the proposed methods to infer the metabolic network for the bacteria Pseudomonas aeruginosa from two genomic datasets: gene position and phylogenetic profiles. Next, we predicted several missing enzyme genes to reconstruct the lysine-degradation pathway in P. aeruginosa using EC number information. As a result, we identified PA0266 as a putative 5-aminovalerate aminotransferase (EC 2.6.1.48) and PA0265 as a putative glutarate semialdehyde dehydrogenase (EC 1.2.1.20). To verify our prediction, we conducted biochemical assays and examined the activity of the products of the predicted genes, PA0265 and PA0266, in a coupled reaction. We observed that the predicted gene products catalyzed the expected reactions; no activity was seen when both gene products were omitted from the reaction.     

Regulation of the human apolipoprotein AIV gene expression in transgenic mice

The apolipoprotein (Apo) AI-CIII-AIV gene cluster has a complex pattern of gene expression that is modulated by both gene- and cluster-specific cis-acting elements. In particular the regulation of Apo AIV expression has been previously studied in vivo and in vitro including several transgenic mouse lines but a complete, consistent picture of the tissue-specific controls is still missing. We have analysed the role of the Apo AIV 3' flanking sequences in the regulation of gene expression using both in vitro and in vivo systems including three lines of transgenic mice. The transgene consisted of a human fragment containing 7 kb of the 5' flanking region, the Apo AIV gene itself and 6 kb of the 3' flanking region (-7+6 Apo AIV). Accurate analysis of the Apo AIV mRNA levels using quantitative PCR and Northern blots showed that the 7+6 kb Apo AIV fragment confers liver-specific regulation in that the human Apo AIV transgene is expressed at approximately the same level as the endogenous mouse Apo AIV gene. In contrast, the intestinal regulation of the transgene did not follow, the pattern observed with the endogenous gene although it produced a much higher intestinal expression following the accepted human pattern. Therefore, this animal model provides an excellent substrate to design therapeutic protocols for those metabolic derangements that may benefit from variations in Apo AIV levels and its anti-atherogenic effect.     

Aberrant Behaviours of Reaction Diffusion Self-organisation Models on Growing Domains in the Presence of Gene Expression Time Delays

Turing’s pattern formation mechanism exhibits sensitivity to the details of the initial conditions suggesting that, in isolation, it cannot robustly generate pattern within noisy biological environments. Nonetheless, secondary aspects of developmental self-organisation, such as a growing domain, have been shown to ameliorate this aberrant model behaviour. Furthermore, while in-situ hybridisation reveals the presence of gene expression in developmental processes, the influence of such dynamics on Turing’s model has received limited attention. Here, we novelly focus on the Gierer–Meinhardt reaction diffusion system considering delays due the time taken for gene expression, while incorporating a number of different domain growth profiles to further explore the influence and interplay of domain growth and gene expression on Turing’s mechanism. We find extensive pathological model behaviour, exhibiting one or more of the following: temporal oscillations with no spatial structure, a failure of the Turing instability and an extreme sensitivity to the initial conditions, the growth profile and the duration of gene expression. This deviant behaviour is even more severe than observed in previous studies of Schnakenberg kinetics on exponentially growing domains in the presence of gene expression (Gaffney and Monk in Bull. Math. Biol. 68:99–130, 2006). Our results emphasise that gene expression dynamics induce unrealistic behaviour in Turing’s model for multiple choices of kinetics and thus such aberrant modelling predictions are likely to be generic. They also highlight that domain growth can no longer ameliorate the excessive sensitivity of Turing’s mechanism in the presence of gene expression time delays. The above, extensive, pathologies suggest that, in the presence of gene expression, Turing’s mechanism would generally require a novel and extensive secondary mechanism to control reaction diffusion patterning.