Published papers

In the post-genomic era, the large-scale, systematic, and functional analysis of all cellular components using transcriptomics, proteomics, and metabolomics, together with bioinformatics for the analysis of the massive amount of data generated by these "omics" methods are the focus of intensive research activities. As a consequence of these developments, systems biology, whose goal is to comprehend the organism as a complex system arising from interactions between its multiple elements, becomes a more tangible objective. Mathematical modeling of microorganisms and subsequent computer simulations are effective tools for systems biology, which will lead to a better understanding of the microbial cell and will have immense ramifications for biological, medical, environmental sciences, and the pharmaceutical industry.In this review, we describe various types of mathematical models (structured, unstructured, static, dynamic, etc.), of microorganisms that have been in use for a while, and others that are emerging. Several biochemical/cellular simulation platforms to manipulate such models are summarized and the E-Cell system 1 developed in our laboratory is introduced. Finally, our strategy for building a "whole cell metabolism model", including the experimental approach, is presented. © 2004 Elsevier B.V. All rights reserved.

All the archaeal genomes sequenced to date contain a single Type 2 RNase H gene. We found that the genome of a halophilic archaeon, Halobacterium sp. NRC-1, contains an open reading frame with similarity to Type 1 RNase H. The protein encoded by the Vng0255c gene, possessed amino acid sequence identities of 33 % with Escherichia coli RNase HI and 34 % with a Bacillus subtilis RNase HI homologue. The B. subtilis RNase HI homologue, however, lacks amino acid sequences corresponding to a basic protrusion region of the E. coli RNase HI, and the Vng0255c has the similar deletion. As this deletion apparently conferred a complete loss of RNase H activity on the B. subtilis RNase HI homologue protein, the Vng0255c product was expected to exhibit no RNase H activity. However, the purified recombinant Vng0255c protein specifically cleaved an RNA strand of the RNA/DNA hybrid in vitro, and when the Vng0255c gene was expressed in an E. coli strain MIC2067 it could suppress the temperature-sensitive growth defect associated with the loss of RNase H enzymes of this strain. These results in vitro and in vivo strongly indicate that the Halobacterium Vng0255c is the first archaeal Type 1 RNase H. This enzyme, unlike other Type 1 RNases H, was able to cleave an Okazaki fragment-like substrate at the junction between the 3′-side of ribonucleotide and 5′-side of deoxyribonucleotide. It is likely that the archaeal Type 1 RNase H plays a role in the removal of the last ribonucleotide of the RNA primer from the Okazaki fragment during DNA replication.

Motivation: Many important problems in cell biology require the dense nonlinear interactions between functional modules to be considered. The importance of computer simulation in understanding cellular processes is now widely accepted, and a variety of simulation algorithms useful for studying certain subsystems have been designed. Many of these are already widely used, and a large number of models constructed on these existing formalisms are available. A significant computational challenge is how we can integrate such sub-cellular models running on different types of algorithms to construct higher order models. Results: A modular, object-oriented simulation meta-algorithm based on a discrete-event scheduler and Hermite polynomial interpolation has been developed and implemented. It is shown that this new method can efficiently handle many components driven by different algorithms and different timescales. The utility of this simulation framework is demonstrated further with a 'composite' heat-shock response model that combines the Gillespie - Gibson stochastic algorithm and deterministic-differential equations. Dramatic improvements in performance were obtained without significant accuracy drawbacks. A multi-timescale demonstration of coupled harmonic oscillators is also shown. © Oxford University Press 2004; all rights reserved.

Hox genes play a critical role in morphogenesis of the early embryo along the anteroposterior axis. In mammals, 39 Hox genes with extensive homology are organized into 13 paralogous groups, forming four clusters on four separate chromosomes. The genes within each cluster are arranged in a 3' to 5' direction and expressed in a temporally and spatially coordinated manner along the anteroposterior axis in the vertebrae, limbs and viscera, including the gastrointestinal tract, but little is known about their spatial expression in the adult gastrointestinal tract. We used the quantitative polymerase chain reaction (PCR) intercalater method with SYBR Green to quantify human Hox gene expression in the adult gastrointestinal tract tissue: esophagus, stomach, duodenum, jejunum, ileum, ileocecum, cecum, ascending colon, transverse colon, descending colon and rectum. Hox gene expression was normalized to glyceraldehyde-3-phosphate-dehydrogenase (GAPDH) gene expression. The spatial expression pattern was analyzed by the multivariate method. The expression level of all 39 Hox genes could be measured in a reproducible manner. Genes with higher expression in the foregut-derived segments tended to have lower expression in hindgut-derived segments, whereas those with low expression in the former tended to have higher in the latter. Principal components analysis and permax analysis revealed a position-specific expression pattern of Hox genes along the anteroposterior axis of the adult gastrointestinal tract. The pattern recapitulates the expression pattern in the embryonic gastrointestinal tract. We suggest that Hox genes may play a pivotal role in the position-specific regenerative process of intestinal epithelial cells.

The complete SfiI and I-CeuI physical maps of four Bacillus subtilis (natto) strains, which were previously isolated as natto (fermented soybean) starters, were constructed to elucidate the genome structure. Not only the similarity in genome size and organization but also the microheterogeneity of the gene context was revealed. No large-scale genome rearrangements among the four strains were indicated by mapping of the genes, including 10 rRNA operons (rrn) and relevant genes required for natto production, to the loci corresponding to those of the B. subtilis strain Marburg 168. However, restriction fragment length polymorphism and the presence or absence of strain-specific DNA sequences, such as the prophages SPβ, skin element, and PBSX, as well as the insertion element IS4Bsu1, could be used to identify one of these strains as a Marburg type and the other three strains as natto types. The genome structure and gene heterogeneity were also consistent with the type of indigenous plasmids harbored by the strains.

We previously reported a computational approach to infer alternative splicing patterns from Mus musculus full-length cDNA clones and microarray data. Although we predicted a large number of unreported splice variants, the general mechanisms regulating alternative splicing were yet unknown. In the present study, we compared alternative exons and constitutive exons in terms of splice-site strength and frequency of potential regulatory sequences. These regulatory features were further compared among five different species: Homo sapiens, M. musculus, Arabidopsis thaliana, Oryza sativa, and Drosophila melanogaster. Solid statistical validations of our comparative analyses indicated that alternative exons have (1) weaker splice sites and (2) more potential regulatory sequences than constitutive exons. Based on our observations, we propose a combinatorial model of alternative splicing mechanisms, which suggests that alternative exons contain weak splice sites regulated alternatively by potential regulatory sequences on the exons.

It still remains intriguing how signal specificity is achieved when different signals are relayed by the common intracellular signal transduction pathways. A well documented example for signal specificity determination is found in rat phaeochromocytoma PC12 cells where epidermal growth factor (EGF) stimulation produces a transient mitogen-activated protein kinase (MAPK) activation and leads to cell proliferation while nerve growth factor (NGF) initiates a sustained MAPK activation and induces cell differentiation. In this simulation, we demonstrated that NGF-induced sustained MARK activation may mainly depend on continual regeneration of NGF receptors and that the presence of a small pool of surface receptors is enough to maintain a sustained MAPK activation. On the other hand, MAPK activation is not significantly sensitive to the half-life of internalized receptors and the levels of NGF-specific MAPK phosphatase MAP kinase phosphatase-3 (MKP-3), though cytoplasmic persistence of internalized NGF-bound receptors and the MKP-3 dependent feedback control also contribute to the sustaining of MAPK activation. These results are consistent with the recent experimental evidence that persistent tyrosine receptor kinase A (TrkA) activity is necessary to maintain transcription in the differentiating PC12 cells (Chang et al. 2003) and a sustained Src kinase activity is detected in response to NGF stimulation (Gatti 2003). It is suggested that sustained or transient MAPK activation induced by different growth factor and neurotrophins, which is crucial to their signaling specificity, could be satisfactorily accounted for by their specific receptor turnover kinetics rather than by the activation of specific downstream signaling cascades.

We have earlier generated a mutant mouse in a course of making a transgenic line that exhibited interesting heterozygote phenotypes, which exhibited failure to thrive, severe bone deformities, and polycystic kidneys. This mutant mouse provided a clue to uncover a unique role of expressed pseudogenes. In this mutant the transgene was integrated into the vicinity of the expressing pseudogene of Makorin1 called Makorin1-p1. This insertion reduced transcription of the Makorin1-p1, resulting in destabilization of the Makorin1 mRNA in trans via a cis-acting RNA decay element within the 5′ region of Makorin1 that is homologous between Makorin1 and Makorin1-p1. These findings demonstrate a novel and specific regulatory role of an expressed pseudogene as well as functional significance for noncoding RNAs. Next, we developed an original algorithm to determine how many pseudogenes are expressed. Based on our examination 2-3% of human processed pseudogenes are expressed using the most strict criteria. Interestingly, the mouse has a much smaller proportion of expressed pseudogenes (0.5-1%). Pseudogenes are functionally less constrained, and have accumulated more mutations than translated genes. If they have some functions in gene regulation, this property would allow more rapid functional diversification than protein-coding genes. In addition, some genetic phenomena that exhibit incomplete penetrance might be attributed to "mutation" or "variation" of pseudogenes. © Springer-Verlag 2004.

For the purpose of analyzing the relation between the splice sites and the order of introns, we conducted the following analysis for the GT-AG and GC-AG splice site groups. First, the pre-mRNAs of H. sapiens, M. musculus, D. melanogaster, A. thaliana and O. sativa were sampled by mapping the full-length cDNA to the genomes. Next, the consensus sequences at different regions of pre-mRNAs were analyzed in the five species. We also investigated the mononucleotide and dinucleotide frequencies in the extensive regions around the 5′ splice sites (5′ss) and 3′ splice sites (3′ss). As a result, differential frequencies of nucleotides at the first 5′ss in both the GT-AG and GC-AG splice site groups were observed in A. thaliana and O. sativa pre-mRNAs. The trend, which indicates that GC 5′ss possess strong consensus sequences, was observed not only in mammalian pre-mRNAs but also in the pre-mRNAs of D. melanogaster, A. thaliana and O. sativa. Furthermore, we examined the consensus sequences of the constitutive and alternative splice sites. It was suggested that in the case of the alternative GC-AG introns, the tendency to have a weak consensus sequence at 5′ss is different between H. sapiens and M. musculus pre-mRNAs. © Imperial College Press.

ST0753, the orthologous gene of Type 1 RNase H found in a thermoacidophilic archaeon, Sulfolobus tokodaii, was analyzed. The recombinant ST0753 protein exhibited RNase H activity in both in vivo and in vitro assays. The protein expressed in an RNase H-deficient mutant Escherichia coli strain functioned to suppress the temperature-sensitive phenotype associated with the lack of RNase H. The in vitro characteristics of the gene's RNase H activity were similar to those of Halobacterium RNase HI, the first archaeal Type 1 RNase H to be characterized. Surprisingly, the S.tokodaii RNase HI cleaved not only the RNA strand of an RNA/DNA hybrid but also an RNA strand of an RNA/RNA duplex in the presence of Mn2+ or Co2+. The result of gel filtration column chromatography showed this double-stranded RNA-dependent RNase (dsRNase) activity was coincident with S.tokodaii RNase HI. A site-directed mutagenesis study of essential amino acids for RNase H activity indicated that this activity also affected dsRNase activity. A single amino acid replacement of Asp-125 by Asn resulted in loss of dsRNase activity but not RNase H activity, suggesting that amino acid residues required for dsRNase activity seemed slightly different from those of RNase H activity. Some reverse transcriptases from retroelements can cleave double-stranded RNA, and this activity requires the RNase H domain. Similarities in primary structure and biochemical characteristics between S.tokodaii RNase HI and reverse transcriptases imply that the S.tokodaii enzyme might be derived from the RNase H domain of reverse transcriptase. © Oxford University Press 2004; all rights reserved.

Background: Elucidating the dynamic behaviour of genetic regulatory networks is one of the most significant challenges in systems biology. However, conventional quantitative predictions have been limited to small networks because publicly available transcriptome data has not been extensively applied to dynamic simulation. Results: We present a microarray data-based semi-kinetic (MASK) method which facilitates the prediction of regulatory dynamics of genetic networks composed of recurrently appearing network motifs with reasonable accuracy. The MASK method allows the determination of model parameters representing the contribution of regulators to transcription rate from time-series microarray data. Using a virtual regulatory network and a Saccharomyces cerevisiae ribosomal protein gene module, we confirmed that a MASK model can predict expression profiles for various conditions as accurately as a conventional kinetic model. Conclusions: We have demonstrated the MASK method for the construction of dynamic simulation models of genetic networks from time-series microarray data, initial mRNA copy number and first-order degradation constants of mRNA. The quantitative accuracy of the MASK models has been confirmed, and the results indicated that this method enables the prediction of quantitative dynamics in genetic networks composed of commonly used network motifs, which cover considerable fraction of the whole network. © 2005 Yugi et al., licensee BioMed Central Ltd.

Background: Cell simulation, which aims to predict the complex and dynamic behavior of living cells, is becoming a valuable tool. In silico models of human red blood cell (RBC) metabolism have been developed by several laboratories. An RBC model using the E-Cell simulation system has been developed. This prototype model consists of three major metabolic pathways, namely, the glycolytic pathway, the pentose phosphate pathway and the nucleotide metabolic pathway. Like the previous model by Joshi and Palsson, it also models physical effects such as osmotic balance. This model was used here to reconstruct the pathology arising from hereditary glucose-6-phosphate dehydrogenase (G6PD) deficiency, which is the most common deficiency in human RBC. Results: Since the prototype model could not reproduce the state of G6PD deficiency, the model was modified to include a pathway for de novo glutathione synthesis and a glutathione disulfide (GSSG) export system. The de novo glutathione (GSH) synthesis pathway was found to compensate partially for the lowered GSH concentrations resulting from G6PD deficiency, with the result that GSSG could be maintained at a very low concentration due to the active export system. Conclusion: The results of the simulation were consistent with the estimated situation of real G6PD-deficient cells. These results suggest that the de novo glutathione synthesis pathway and the GSSG export system play an important role in alleviating the consequences of G6PD deficiency. © 2005 Nakayama et al; licensee BioMed Central Ltd.

For many decades microorganisms have been used for industrial purposes; traditional fermentations such as brewing and production of food additives, aroma molecules, organic acids and pharmaceutical-like antibiotics or recombinant proteins are instances of the industrial microorganism utilization. Therefore, microorganism modeling and simulation have been required for engineering purposes, because of demands for design, optimization and quality control of large-scale fermentation plants. Modeling has recently become more highly developed, aided by the deciphering of microorganism genomes, the completion of metabolic databases, the development of analytical methodologies and improvements in the performance of computers. This paper reviews past and recent metabolic simulation of microorganisms, and also discusses the metabolome analytical techniques and the construction of large-scale microorganism models which are now being developed in our group. © 2005 Springer Science+Business Media, Inc.

Background: Many computer studies have employed either dynamic simulation or metabolic flux analysis (MFA) to predict the behaviour of biochemical pathways. Dynamic simulation determines the time evolution of pathway properties in response to environmental changes, whereas MFA provides only a snapshot of pathway properties within a particular set of environmental conditions. However, owing to the large amount of kinetic data required for dynamic simulation, MFA, which requires less information, has been used to manipulate large-scale pathways to determine metabolic outcomes. Results: Here we describe a simulation method based on cooperation between kinetics-based dynamic models and MFA-based static models. This hybrid method enables quasi-dynamic simulations of large-scale metabolic pathways, while drastically reducing the number of kinetics assays needed for dynamic simulations. The dynamic behaviour of metabolic pathways predicted by our method is almost identical to that determined by dynamic kinetic simulation. Conclusion: The discrepancies between the dynamic and the hybrid models were sufficiently small to prove that an MFA-based static module is capable of performing dynamic simulations as accurately as kinetic models. Our hybrid method reduces the number of biochemical experiments required for dynamic models of large-scale metabolic pathways by replacing suitable enzyme reactions with a static module. © 2005 Yugi et al; licensee BioMed Central Ltd.

Motivation: Transcription start site selection and alternative splicing greatly contribute to diversifying gene expression. Recent studies have revealed the existence of alternative first exons, but most have involved mammalian genes, and as yet the regulation of usage of alternative first exons has not been clarified, especially in plants. Results: We systematically identified putative alternative first exon transcripts in rice, verified the candidates using RT-PCR, and searched for the promoter elements that might regulate the alternative first exons. As a result, we detected a number of unreported alternative first exons, some of which are regulated in a tissue-specific manner. © The Author 2005. Published by Oxford University Press. All rights reserved.

The translation start site, immediately downstream from the start codon, is a dominant factor for gene expression in Escherichia coli. At present, no method exists to improve the expression level of cloned genes, since it remains difficult to find the best codon combination within the region. We determined the expression parameters that correspond to all sense codons within the first four codons using GFPuv which encodes a derivative of green fluorescent protein. Using a genetic algorithm (GA)-based computer program, these parameters were incorporated in a simple, static model for the prediction of translation efficiency, and optimized to the expression level for 137 randomly isolated GFPuv genes. The calculated initial translation index (ITI), also proven for the DsRed2 gene that encodes a red fluorescent protein, should provide a solution to overcome the gene expression problem in cloned genes whose expression is often inherently blocked at the translation process. The proposed method facilitates heterologous protein production in E. coli, the most commonly used host in biological and industrial fields. © 2004 Published by Elsevier B.V.

Increased research aimed at simulating biological systems requires sophisticated parameter estimation methods. All current approaches, including genetic algorithms, need pre-existing equations to be functional. A generalized approach to predict not only parameters but also biochemical equations from only observable time-course information must be developed and a computational method to generate arbitrary equations without knowledge of biochemical reaction mechanisms must be developed. We present a technique to predict an equation using genetic programming. Our technique can search topology and numerical parameters of mathematical expression simultaneously. To improve the search ability of numeric constants, we added numeric mutation to the conventional procedure. As case studies, we predicted two equations of enzyme-catalyzed reactions regarding adenylate kinase and phosphofructokinase. Our numerical experimental results showed that our approach could obtain correct topology and parameters that were close to the originals. The mean errors between given and simulation-predicted time-courses were 1.6 × 10-5% and 2.0 × 10-3%, respectively. Our equation prediction approach can be applied to identify metabolic reactions from observable time-courses. © 2004 Elsevier Ireland Ltd. All rights reserved.

Biological simulation serves to unify the basic elements of systems biology, namely, model selection, experimentation and model refinement. To select biochemical models for simulation, metabolome analysis can be performed using capillary electrophoresis or liquid chromatography coupled with mass spectrometry. In this manner, selected models can be elaborated with temporal/spatial gene and protein expression data obtained from model organisms such as Escherichia coli. The E. coli single gene deletion mutant library (KO collection) and His-tag/GFP-fusion single open reading frame clone expression library (ASKA) are powerful resources for this task. The integration of parallel experimental datasets into dynamic simulation tools forms the remaining challenge for the systematic analysis and elucidation of biological networks and holds promise for biotechnological applications. © 2005 Elsevier Ltd. All rights reserved.

The SCO2299 gene from Streptomyces coelicolor encodes a single peptide consisting of 497 amino acid residues. Its N-terminal region shows high amino acid sequence similarity to RNase HI, whereas its C-terminal region bears similarity to the CobC protein, which is involved in the synthesis of cobalamin. The SCO2299 gene suppressed a temperature-sensitive growth defect of an Escherichia coli RNase H-deficient strain, and the recombinant SCO2299 protein cleaved an RNA strand of RNA-DNA hybrid in vitro. The N-terminal domain of the SCO2299 protein, when overproduced independently, exhibited RNase H activity at a similar level to the full length protein. On the other hand, the C-terminal domain showed no CobC-like activity but an acid phosphatase activity. The full length protein also exhibited acid phosphatase activity at almost the same level as the C-terminal domain alone. These results indicate that RNase H and acid phosphatase activities of the full length SCO2299 protein depend on its N-terminal and C-terminal domains, respectively. The physiological functions of the SCO2299 gene and the relation between RNase H and acid phosphatase remain to be determined. However, the bifunctional enzyme examined here is a novel style in the Type 1 RNase H family. Additionally, S. coelicolor is the first example of an organism whose genome contains three active RNase H genes. © 2005 FEBS.

We have developed a comprehensive software suite for bioinformatics research of cDNAs; it is aimed at rapid characterization of the features of genes and the proteins they code. Methods implemented include the detection of translation initiation and termination signals, statistical analysis of codon usage, comparative study of amino acid composition, comparative modeling of the structures of product proteins, prediction of alternative splice forms, and metabolic pathway reconstruction.