EP1752899A1

EP1752899A1 - Method for simulating a production process of a substance

Info

Publication number: EP1752899A1
Application number: EP06016619A
Authority: EP
Inventors: Yoshihiro Usuda; Shintaro Iwatani; Kazuhiko Matsui
Original assignee: Ajinomoto Co Inc
Current assignee: Ajinomoto Co Inc
Priority date: 2005-08-09
Filing date: 2006-08-09
Publication date: 2007-02-14
Also published as: JP5011675B2; JP2007047994A; US20070038419A1

Abstract

A method for simulating a substance production process using cells based on differential equations concerning intracellular metabolites and gene expression comprising preparing the differential equations concerning intracellular metabolites and gene expression, assigning values for parameters in the differential equations required for the simulation and solving the differential equations with the assigned values, which is characterized by the following (a) to (d): (a) the differential equations include a differential equation including a specific growth rate of the cells, and the specific growth rate is represented as a function of time; (b) all or a part of parameters of the differential equations are represented as a function of the specific growth rate or a function of time; (c) the differential equations include a formation rate for formation of a cell component from an intracellular metabolite, and the formation rate is represented as a function of the specific growth rate or a function of time; and (d) the differential equations include an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, and the inflow rate and/or the outflow rate is represented as a function of the specific growth rate or a function of time.

Description

BACKGROUND OF THE INVENTION

The present invention relates to a method for simulating a production process of a substance, typically an amino acid or a nucleic acid, using cells of a microorganism or the like, and a program therefor.
The technique of constructing a mathematical equation model of biochemical reactions caused by intracellular enzymes to estimate intracellular dynamic behaviors of metabolites is called metabolic simulation, and there are many examples thereof (Ishii, N. et al., J. Biotechnol., 113:281-294, 2004) and proposed many methods therefor ( U.S. Patent Application No. 2002/0022947 , International Publication Nos. WO2004/081862 , WO03/07217 , WO02/55995 , WO02/05205 , Japanese Patent Application Laid-Open (KOKAI) No. 2003-180400 ). As an example of comparatively large scale metabolic simulation, Teusink et al. performed simulation of anaerobic ethanol fermentation of Saccharomyces cerevisiae considering metabolic routes branching from the glycolytic system for producing glycogen, trehalose, glycerol and succinic acid (Teusink, B. et al., Eur. J. Biochem., 267:5313-5329, 2000). For Escherichia coli, Chassagnole et al. constructs a central metabolic model under a condition of a constant growth rate to perform simulation (Chassagnole, C. et al., Biotechnol. Bioeng., 26:203-216, 2002). According to another report, a part of the gene expression of E. coli metabolic enzymes was modeled and combined with an enzymatic reaction model to perform simulation (Wang, J. et al., J. Biotechnol., 92:133-158, 2001; Schmid J.W. et al., Metab. Eng., 6:364-377, 2004). Further, in the report of Varner, construction of a large scale model in which gene expression is incorporated into a kinetic model of enzymes is conceptually disclosed, and a specific growth rate is expressed by an equation using saturation coefficients of precursors for the maximum specific growth rate (Varner, J.D., Biotechnol. Bioeng., 69:664-678, 2000). For other organisms, further detailed models have been reported. Jeong et al. constructed a model of sporulation process of Bacillus subtilis in batch culture using mathematical equations (Jeong et al., Biotechnol. Bioeng., 35:160-184, 1990). Tomita et al., Bioinformatics 15:72-84, 1999, also has reported on the simulation of 127 genes involved in transcription and translation of the Mycoplasma genitalium genome, as well as energy production and phospholipid synthesis of the microorganism.

SUMMARY OF THE INVENTION

In simulation of a production process of a substance, typically a nucleic acid or an amino acid, using cells of microorganisms or the like, enzymatic reactions from a substrate to an objective product are represented by mathematical equations using kinetic parameters in many cases. However, batch culture or semi-batch culture (fed batch culture) is often used for production of a useful substance, and therefore the growth rate of cells changes. In connection with it, various kinds of parameters in the cells also change. In addition, besides the production rates of a substance serving as a substrate, components present in the medium and an objective product, production rates of by-products such as amino acids, organic acids and carbon dioxide (CO₂) also change during the process of substance production. Therefore, a technique for performing metabolic simulation with sufficient precision in such a manner that the simulation should well fit to experimental data of such growth rates or by-products is desired. If an accurate metabolic simulation reflecting experimental data is enabled, it becomes possible to conduct experiments of amplification or deletion of a gene by a computer in a short time (in silico experiments). It is expected that it should greatly shorten the development period for improving substance production ability of cells.
The inventors of the present invention conducted various experiments in view of the aforementioned problems. Consequently, they found that they could achieve more accurate metabolic simulation of the production of a substance, typically an amino acid or a nucleic acid, by cells of microorganisms or the like. The enhanced accuracy pertained when (A) a specific growth rate, serving as an index of cell growth, was represented by a mathematical equation that used a time function based on measured values, and (B) time functions or functions using the specific growth rate as a variable were employed with various parameters, including outflow rates of intracellular metabolites into cell components and, further, uptake rates of intracellular metabolites from the outside of the cells or excretion rates of the same to the outside of the cells.
The present invention was accomplished based on the aforementioned findings and provides the following.

[1] A method for effecting a simulation of a substance-production process that uses cells, wherein said simulation is based on a set of differential equations that represent intracellular metabolites and gene expression, said method comprising the steps of:
1. (a) including a specific growth rate of cells, expressed as a differential equation, in the set of differential equations;
2. (b) assigning values for parameters in the set of differential equations, wherein at least one of said parameters is represented as a growth rate factor;
3. (c) incorporating in the set of differential equations a formation rate for formation of a cell component from an intracellular metabolite, wherein said formation rate is represented as a growth rate factor;
4. (d) incorporating in the set of differential equations an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, wherein said inflow rate and said outflow rate are represented, respectively, as a growth rate factor;
5. (e) solving the set of differential equations; and
6. (f) generating data representative of the substance-production process.
[2] The method according to [1], wherein the differential equations including the specific growth rate of the cells include the differential equations represented as the following equations (1) to (3): $d [Metabolite] / dt = V_{input} - V_{output} - µ [Metabolite]$
$d [mRNA] / dt = k_{tanscription} [P] - (k_{dRNA} + µ) [mRNA]$
$d [Protein] / dt = k_{translation} [mRNA] - (k_{dProtein} + µ) [Protein]$
wherein, in the equation 1, [Metabolite] represents an intracellular concentration of a metabolite, V_input represents the sum of rates of reactions producing the metabolite, V_output represents the sum of rates of reactions consuming the metabolite, and µ represents the specific growth rate;
in the equation 2, [mRNA] represents a concentration of mRNA, k_{transcription} represents a rate constant of transcription, [P] represents a promoter concentration, k_dRNA represents a rate constant of decomposition of mRNA, and µ represents the specific growth rate, and
in the equation 3, [Protein] represents a concentration of a protein, k_translation represents a rate constant of translation, k_dProtein represents a rate constant of decomposition of the protein, and µ represents the specific growth rate.
[3] The method according to [1], wherein the growth rate factor is a function of the specific growth rate or a function of time.
[4] The method according to [1], wherein the specific growth rate is represented as a function of time, and the function is obtained by preparing a mathematic equation from measurement data of the specific growth rate in the production process.
[5] The method according to [1], wherein the growth rate factor representing the formation rate is obtained by preparing a mathematic equation expressing measurement data of the formation rate in the production process.
[6] The method according to [1], wherein the growth rate factor representing the inflow rate and/or the outflow rate is obtained by preparing a mathematic equation expressing measurement data of the inflow rate and/or the outflow rate in the production process.
[7] The method according to [1], wherein the metabolite taken up into the cells is a substrate and/or an organic substance in a medium.
[8] The method according to [1], wherein the metabolite excreted out of the cells is an objective substance and/or a by-product.
[9] The method according to [8], wherein the metabolite excreted out of the cells is an amino acid, an organic acid and/or carbon dioxide.
[10] The method according to [9], wherein the metabolite excreted out of the cells is an amino acid or an organic acid.
[11] The method according to [1], wherein the parameters required for the simulation are a rate constant of transcription and/or a rate constant of translation.
[12] The method according to [1], wherein the cells are those of a microorganism having an amino acid producing ability and/or an organic acid producing ability.
[13] The method according to [12], wherein the microorganism is Escherichia coli.
[14] The method according to [1], wherein a composition of cell components itself is represented by a mathematical equation using the specific growth rate of the cells or the cells' equivalent index concerning the growth.
[15] A computer program product for effecting a simulation of a substance-production process that uses cells, wherein said simulation is based on a set of differential equations that represent intracellular metabolites and gene expression, comprising:
1. (a) computer code for including a specific growth rate of cells, expressed as a differential equation, in the set of differential equations;
2. (b) computer code for assigning values for parameters in the set of differential equations, wherein at least one of said parameters is represented as a growth rate factor;
3. (c) computer code for incorporating in the set of differential equations a formation rate for formation of a cell component from an intracellular metabolite, wherein said formation rate is represented as a growth rate factor;
4. (d) computer code for incorporating in the set of differential equations an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, wherein said inflow rate and said outflow rate are represented, respectively, as a growth rate factor;
5. (e) computer code for solving the set of differential equations; and
6. (f) computer code for generating data representative of the substance-production process.
[16] A system for effecting a simulation of a substance-production process that uses cells, wherein said simulation is based on a set of differential equations that represent intracellular metabolites and gene expression, comprising: a processor for processing information; and
a storing means, including:
1. (a) computer code for including a specific growth rate of cells, expressed as a differential equation, in the set of differential equations;
2. (b) computer code for assigning values for parameters in the set of differential equations, wherein at least one of said parameters is represented as a growth rate factor;
3. (c) computer code for incorporating in the set of differential equations a formation rate for formation of a cell component from an intracellular metabolite, wherein said formation rate is represented as a growth rate factor;
4. (d) computer code for incorporating in the set of differential equations an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, wherein said inflow rate and said outflow rate are represented, respectively, as a growth rate factor;
5. (e) computer code for solving the set of differential equations; and
6. (f) computer code for generating data representative of the substance-production process.

According to the method of the present invention, it becomes possible to perform simulation of a production process of a substance, typically an amino acid or a nucleic acid, for a production process using cells of a microorganism or the like, in which a growth rate markedly changes.

BRIEF EXPLANATION OF DRAWINGS

Fig. 1 shows a diagram for explaining modeling of gene expression.
Fig. 2 shows examples of a growth curve obtained by measurement of turbidity (OD) (A) and a time equation of a specific growth rate µ obtained from the curve (B).
Fig. 3 shows a diagram for explaining the material balance of metabolites.
Fig. 4 shows examples of a concentration change curve obtained by measurement of an extracellular acetic acid (AcOH) concentration (A) and a time equation of an excretion rate obtained from the curve (B).
Fig. 5 shows a flowchart of a process to be executed by the program of the present invention.
Fig. 6 shows a functional block diagram of a computer on which the program of the present invention is installed.
Fig. 7 shows a modeled central metabolic map of E. coli. Metabolites, enzymes and genes are indicated with abbreviations. The arrows represent conversion between metabolites by enzymatic reactions, for each of which name of enzyme (capital letters) and gene encoding the enzyme (small letters in italics) are indicated. The proteins in the boxes are transcription factors, and the broken lines represent interactions with an effector. The dotted lines represent outflow from intracellular metabolites to cell components. The thick black frame represents the cell membrane, the substances outside the cells are extracellular substances, and the boxes on the cell membrane indicate uptake or excretion.
Figs. 8A and 8B show results of metabolic simulation of the E. coli central metabolic model: cell volume (growth) (A), extracellular glucose (B), G6P: (C), FDP (D), GA3P (E), PEP (F), PYR (G), 6PGC (H), R5P (I), ACCoA (J), CIT (K), AKG (L), SUCCoA (M), FUM (N), OAA (0), and extracellular CO₂ (P). For the extracellular glucose (B) and extracellular CO₂ (P), actual measured values are plotted as broken lines. The horizontal axes represent the culture time (minute) (the same shall apply to Figs. 9A, 9B, 10A, 10B, 11A, 11B, 12A and 12B).
Figs. 9A and 9B show results of gene expression simulation of the E. coli central metabolic model: ptsHI mRNA (A), EI (B), HPr (C), ptsG mRNA (D), IICBGlc (E), IICBGlc activity (F), pfkA mRNA (G), PFKA (H), PFKA activity (I), zwf mRNA (J) G6PD (K), G6PD activity (L), gltA mRNA (M), CS (N), CS activity (0), ppc mRNA (P), PPC (Q), and PPC activity (R). The values of metabolic fluxes at 315 minutes and 495 minutes converted into enzymatic reaction rates are indicated with closed circles.
Figs. 10A and 10B show results of gene expression simulation under a condition that RNA polymerase and ribosome concentrations are independent from µ: ptsHI mRNA (A), EI (B), HPr (C), ptsG mRNA (D), IICBGlc (E), IICBGlc activity (F), pfkA mRNA (G), PFKA (H), PFKA activity (I), zwf mRNA (J), G6PD(K), G6PD activity (L), gltA mRNA (M), CS (N), CS activity (0), ppc mRNA (P), PPC (Q), and PPC activity (R). The values of metabolic fluxes at 315 minutes and 495 minutes converted into enzymatic reaction rates are indicated with closed circles.
Figs. 11A and 11B show results of metabolic simulation when synthesis rates for cell components are set at 0: cell volume (growth) (A), extracellular glucose (B), G6P: (C), FDP (D), GA3P (E), PEP (F), PYR (G), 6PGC (H), R5P (I), ACCoA (J), CIT (K), AKG (L), SUCCoA (M), FUM (N), OAA (O), and extracellular CO₂ (P). For the extracellular glucose (B) and extracellular CO₂ (P), actual measured values are plotted as broken lines.
Figs. 12A and 12B show results of metabolic simulation when uptake and excretion of metabolites are set at 0: cell volume (growth) (A), extracellular glucose (B), G6P: (C), FDP (D), GA3P (E), PEP (F), PYR (G), 6PGC (H), R5P (I), ACCoA (J), CIT (K), AKG (L), SUCCoA (M), FUM (N), OAA (O), and extracellular CO₂ (P). For the extracellular glucose (B) and extracellular CO₂ (P), actual measured values are plotted as broken lines.

DETAILED DESCRIPTION OF THE INVENTION

Hereafter, the present invention will be explained in detail.
In this description, the phrase "substance production process using cells" and similar terminology refers to a process of biochemical reactions from a substrate to a product that is caused as sequential enzymatic reactions, using cells to produce an objective product.
The simulation method of the present invention is based on differential equations, which may be the same as a conventional simulation methodology except that specific conditions according to the present invention are employed. Conventional simulation methods comprise preparing differential equations for intracellular metabolites and gene expression, assigning values for parameters in the differential equations required for the simulation and solving the differential equations with the assigned values.
Differential equations can usually be prepared by incorporating mathematical equations for expression control into metabolic simulation.
Metabolic simulation is a technique for describing time-dependent dynamic changes by expressing intracellular biochemical reactions with mathematical equations, describing changes of substances caused by the reactions with differential equations, and solving them by numerical computation. The mathematical equations used for this purpose are often nonlinear ordinary differential equations (ODE), and this process of preparing mathematical equations is generally called "modeling." Typically, many nonlinear differential equations are solved by using a computer.
In this regard, the phrase "biochemical reactions" refers to processes for converting an intracellular metabolite by an enzymatic reaction. Data on such reactions are stored in databases for many organisms. For example, KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.ad.jp/kegg/; Kanehisa, M. et al., Nucleic Acids Res., 32:277-280, 2004) can be referred to. As for E. coli, EcoCyc (Encyclopedia of Escherichia coli K12 Genes and Metabolism, http://ecocyc.org/, Keseler et al. (2005) Nucleic Acids Res., 33, D334-D337, 2005) is known. For describing these biochemical enzymatic reactions with mathematical equations, dynamic equations based on Michaelis-Menten type reaction formulas are often used (Segel I.H., Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady-State Enzyme Systems, John Wiley & Sons, 1975). Parameters of respective enzymes can be collected from literatures. For example, for E. coli, Chassagnole et al. described enzymatic reactions of the glycolytic pathway and pentose phosphate pathway from glucose to acetyl-CoA with mathematical equations using kinetic parameters (Chassagnole, C. et al., Biotechnol. Bioeng., 26, 203-216, 2002).
It is gene expression that determines amounts of intracellular enzymes, and this process is realized through the processes of transcription of a gene into mRNA and translation of mRNA into a protein. By incorporating this gene expression into metabolic simulation, it becomes possible to describe more detailed intracellular behaviors. Examples include description of expression control with mathematical equations and incorporation of them into metabolic simulation. For example, Wang et al. described expression control of the sucrose and glycerol uptake system of E. coli with mathematical equations and combined them with an enzymatic reaction model of the glycolytic pathway to perform simulation (Wang, J. et al., J. Biotechnol., 92:133-158, 2001), and Schmid et al. modeled expression control of a tryptophan biosynthesis pathway gene (trp operon) and combined it with a central metabolic model of E. coli to perform analysis (Schmid J.W. et al., Metab. Eng., 6:364-377, 2004). Differential equations concerning intracellular metabolites and gene expression are prepared as described above.
Concentration change and gene expression of a metabolite can be generally represented by one or more differential equations. Thus, one commonly represents each intracellular metabolite with an equation 1, considering the enzymatic reaction rate for each intracellular metabolite and the dilution effect due to growth: $d [Metabolite] / dt = V_{input} - V_{output} - μ [Metabolite]$
In the equation 1, [Metabolite] represents the concentration of an intracellular metabolite, V_input represents the sum of rates of enzymatic reactions for producing the metabolite, V_output represents the sum of rates of enzymatic reactions consuming the metabolite, and µ represents the specific growth rate.
Gene expression is generally described in terms of the two stages of transcription and translation, i.e., as changes in concentrations of mRNA produced by transcription of a gene and protein produced by translation of mRNA (Fig. 1).
For mRNA produced by transcription of various genes, one may describe the concentration by the following equation, considering the transcription rate with which mRNA is synthesized from a gene and the decomposition rate of mRNA as well as the dilution effect: $d [mRNA] / dt = k_{transcription} [RNAP] [Promoter] - (k_{dmRNA} + μ) [mRNA]$
In equation 2, [mRNA] represents the concentration of mRNA, k_{transcription} represents a rate constant of transcription, [RNAP] represents the concentration of RNA polymerase which performs the transcription, [Promoter] represents the concentration of a promoter of the corresponding gene, k_dRNA represents a rate constant of decomposition of mRNA, and µ represents the specific growth rate.
For a protein produced by translation of mRNA, the concentration can be described by the following equation, which takes into account the translation rate with which mRNA is translated into a protein and decomposition rate of the protein as well as the dilution effect: $d [Protein] / dt = k_{translation} [Ribosome] [mRNA] - (k_{dProtein} + μ) [Protein]$
In equation 3, [Protein] represents the concentration of a protein, k_translation represents a rate constant of translation, [Ribosome] represents the concentration of ribosome which performs translation, k_dProtein represents a rate constant of protein decomposition, and µ represents the specific growth rate.
Equations 2 and 3 can also be described in various other ways.
When differential equations are prepared, values are assigned for the parameters required for simulation in the prepared differential equations. The parameters required for the simulation include rate constants, initial concentrations, and so forth. The parameters are preferably a rate constant of transcription and/or a rate constant of translation.
Parameters for respective enzymatic reactions and gene expression can be collected from literatures. However, there are many parameters for modeling of enzymatic reactions and gene expression, and therefore it is often impossible to collect all the parameters from literatures. In such a case, it is possible to estimate appropriate values from various information in literatures, or estimate them by optimization of measurement results obtained in experiments.
In order to solve such nonlinear ordinary differential equations obtained as described above, it is possible to use a mathematical calculation program such as MATLAB® (MathWorks) and MATHEMATICA® (Wolfram Research). To execute metabolic simulation, for example, ODE solvers of MATLAB® (MathWorks) can be used. As the ODE solver for solving a metabolic reaction equation or gene expression equation, ode45, ode21s and so forth are preferably used. Moreover, many kinds of metabolic simulation software have been developed as software for performing metabolic simulation, and they can be utilized (Ishii, N. et al., J. Biotechnol., 113:281-294, 2004). Examples include GEPASI (Mendes, P. Comput. Applic. Biosci., 9:563-571, 1993), SCAMP (Sauro, H.M., Comput. Appl. Biosci., 9:441-450, 1993), E-CELL (Tomita, M. et al., Bioinformatics, 15:72-84, 1999), and so forth.
A simulation of the present invention is characterized by additionally using certain conditions, as described in greater detail below. This approach enables more accurate simulation of a substance production process accompanied by growth of cells.

(a) Description of specific growth rate with time function

A differential equation that incorporates a specific growth rate of cells is included in the differential equations, and the specific growth rate is represented with a function of time.
The phrase "growth of cells" refers to a phenomenon whereby the number of cells increases in a substance production process. The number of cells usually increases with conversion of a substrate added for the substance production into cell components. When growth is represented by the number of cells, the growth rate is the rate at which the number of cells increases, and a specific growth rate is obtained by dividing the increase rate of the number of cells with the number of cells. The number of cells used in this context is one of the values serving as indexes of growth of cells, and any value may be used so long as a value having an equivalent function (for example, turbidity of culture broth) is chosen.
The function of time of the specific growth rate is preferably obtained by preparing mathematical equations that express measurement data of the specific growth rate in a production process. Growth of cells can be measured by measuring turbidity of culture broth or counting the number of cells in a diluted culture broth. A curve of cell growth experimentally measured can be approximated as a function of time. A large number of products are marketed as software for obtaining such an approximate mathematical equation, and preferred examples include TableCurve® 2D (Systat Software). The obtained time function of growth curve can be differentiated and divided with the equation of the growth curve to obtain a time equation of the specific growth rate. Examples of the growth curve obtained by measurement of turbidity (OD) and the time equation of the specific growth rate µ obtained therefrom are shown in Fig. 2.

(b) Representation of parameters with a growth rate factor

One or more of the parameters of the differential equations are represented as a growth rate factor.
The parameters change with the growth of cells. Thus, it is preferable to represent as many parameters as possible with a growth rate factor. The growth rate factor can be expressed as a function of time. A function of time of the specific growth rate is generated with measurement data relating to the specific growth rate obtained in a production process. In the alternative, the growth rate factor can be expressed as a function of the specific growth rate µ. The specific growth rate µ may be obtained by dividing the rate of increase of the number of cells by the number of cells.
For example, it is known that rate constants of transcription and translation, which are parameters of gene expression, markedly change with the growth rate of cells, and molecular numbers of RNA polymerase and ribosome, which catalyze the respective processes, markedly change with the culture rate (Bremer and Dennis, In Escherichia coli and Salmonella: Cellular and Molecular Biology/Second Edition (Neidhardt, F.C. Ed.), pp. 1553-1569, American Society for Microbiology Press, Washington, D.C., 1996). If the rate constants of transcription and translation are expressed as specific growth rate-dependent mathematical equations and incorporated into mathematical equations for gene expression, it becomes possible to prepare mathematical equations that accurately describe the expression. More specifically, in the aforementioned differential equation of mRNA (equation 2), [RNAP] can be represented with an equation of specific growth rate µ. Similarly, in the aforementioned differential equation of protein (equation 3), [Ribosome] can be represented with an equation of the specific growth rate µ. The specific growth rate µ used herein may be the same as that used in (a) mentioned above, or a specific growth rate based on another index relating to growth of cells.
As mentioned above, it is also possible to directly express values measured in a cell growth process as a function of time. For example, if a function representing a parameter is represented with a specific growth rate-dependent equation, the parameter can be converted into a time function by substituting a time function of the specific growth rate into the specific growth rate-dependent equation.

(c) Description of cell component formation rate with a growth rate factor

A formation rate with which a cell component is formed from an intracellular metabolite is incorporated into the differential equations, and the formation rate is represented as a function of a growth rate factor. As set forth above, the growth rate factor can be expressed as a function of time or as a function of the specific growth rate µ.
The phrase "cell component" refers to a major polymer compound constituting cells such as protein, RNA, DNA, lipid and lipopolysaccharide. By enzymatic conversion of a substrate into a cell component such as protein, nucleic acid and lipid, cells can obtain a required component. By enumerating biochemical reactions resulting in such a cell component, a stoichiometric matrix can be created (Savinell J.M., Palsson B.O., J. Theor. Biol., 154:421-454, 1992; Vallino, J.J. and Stephanopoulos, G., Biotechnol. Bioeng., 41:633-646, 1993). Details of creation of a stoichiometric matrix of metabolic reactions from glucose to all the amino acids in E. coli are described in W02005/001736 in detail. If a composition of the cell components is given, the formation rates of the cell components from intracellular metabolites, V_biomass, can be calculated by using the stoichiometric matrix (Pramanik, J. and Keasling J.D., Biotechnol. Bioeng., 56:398-421, 1997).
It is also known that the cell component formation rate is growth rate-dependent (Pramanik, J. and Keasling J.D., Biotechnol. Bioeng., 60:230-238, 1998). Further, the composition of cell components also depends on the growth rate, and by measuring it, it can be described with mathematical equations by using the specific growth rate (Bremer and Dennis, In Escherichia coli and Salmonella: Cellular and Molecular Biology/Second Edition (Neidhardt, F.C. Ed.), pp. 1553-1569., American Society for Microbiology Press, Washington, D.C., 1996).
With the cell component formation rate V_biomass obtained as described above, the influence on intracellular metabolites can be taken into consideration for metabolites as precursors of the cell components. Specifically, by incorporating V_biomass into V_output in the differential equation of a metabolite as a precursor of the cell component (equation 1) to perform calculation, it becomes possible to describe the material balance of the metabolite as a precursor of the cell component with a growth rate factor equation.
If the formation rate of the cell component can be measured, it is also possible to represent values measured in a production process as a time function and use it directly. Moreover, if the function representing the formation rate is represented with a specific growth rate-dependent equation, the formation rate can be converted into a time function by substituting a time function of the specific growth rate into the specific growth rate-dependent equation. The specific growth rate may be the same as that used in (a) mentioned above, or a specific growth rate based on another index relating to growth of cells.
Moreover, it is preferable to represent the composition of cell components itself with a mathematical equation using the specific growth rate of cells or its equivalent index concerning the growth.

(d) Preparation of mathematical equations representing inflow rate of metabolites from outside of cells and outflow rate of metabolites out of cells

An inflow rate of a metabolite taken up from the outside of cells and/or an outflow rate of a metabolite excreted out of cells are incorporated into the differential equations, and the inflow rate and/or the outflow rate is represented as a growth rate factor.. As set forth above, the growth rate factor can be expressed as a function of time or as a function of the specific growth rate µ.
In a production process of a substance, it is common to add an organic substance other than the substrate such as glucose as a medium component. Examples include tryptone, soybean hydrolysate, yeast extract and so forth. Substances such as amino acids derived from such medium components are also taken up into cells and affect the metabolic simulation. It is possible to describe an uptake rate of a metabolite taken up from the inside of cells V_uptake with a function of the specific growth rate or time function based on measured values of concentrations of medium components remaining in the medium. As the metabolites excreted out of the cells from the inside of the cells during a production process, substances called by-products may also be excreted other than the objective product. By also incorporating these substances into the metabolic simulation, more accurate material balance can be described. The outflow rate to the outside of the cells Vexcretion can be represented with a mathematical equation using a growth rate factor. The growth rate factor can be a function of the specific growth rate or a time function from measured values of the metabolite concentration detected in the medium. By performing calculation with incorporating V_uptake into V_input and V_excretion into V_input in the differential equation of a metabolite (equation 1), the material balance depending on the growth rate of the cells or time can be described (Fig. 3).
If the inflow rate and/or the outflow rate can be measured, it is also possible to represent values measured in a production process as a time function and use it directly. For instance, Fig. 4 shows the concentration change curve, obtained by measurement of the extracellular acetic acid concentration, and the time equation of the excretion rate obtained therefrom. Moreover, if the function representing the inflow rate and/or the outflow rate is represented with a specific growth rate-dependent equation, the inflow rate and/or the outflow rate can be converted into a time function by substituting a time function of the specific growth rate into the specific growth rate-dependent equation. The specific growth rate may be the same as that used in (a) mentioned above, or a specific growth rate based on another index relating to growth of cells.
The metabolites taken up into the cells preferably are a substrate and/or an organic substance in the medium.
The substances excreted out of the cells preferably are an objective product and/or a by-product. The substances excreted out of the cells more preferably are amino acids, organic acids and/or carbon dioxide, further preferably amino acids or organic acids.
The cells used for the production process may be those of any type, so long as those used for substance production are chosen. Examples include, for example, various cultured cells, those of mold, yeast, various bacteria, and so forth. Preferred are those of a microorganism having an ability to produce a useful compound, for example, an amino acid, nucleic acid or organic acid. As a microorganism having an ability to produce an amino acid, nucleic acid or organic acid, E. coli, Bacillus bacteria, coryneform bacteria and so forth are preferably used. More preferred are those of a microorganism having an ability to produce an amino acid and/or an ability to produce an organic acid. The microorganism is preferably Escherichia coli.
By the simulation according to the method of the present invention, it becomes possible to predict dynamic behaviors of mRNA or protein concentrations for various enzymes in addition to intracellular metabolites. Therefore, the method of the present invention can serve as a useful tool in improvement of a production process of a useful substance, typically an amino acid or nucleic acid. For example, it becomes possible to verify effect of amplification or deletion of various enzymes in a computer (in silico experiments). Moreover, easy estimation of influence of change of parameters of various enzymes such as affinity to a substrate and affinity to an inhibitor on the whole metabolism and effect of amplification or deletion of a factor controlling expression of various enzyme genes also becomes possible. These results provide an important direction for improvement of a production process, and thus also have superior industrial usefulness.
The present invention further provides a program for executing the simulation method of the present invention and a storage means storing the program.
The program of the present invention is a program for making a computer execute the simulation method of the present invention and causes a computer to function as the following means (1) to (3):

(1) a means for storing a set of differential equations concerning intracellular metabolites and gene expression and satisfying the following (a) to (d);
1. (a) the differential equations include a specific growth rate of the cells, expressed as a differential equation wherein the specific growth rate is represented as a function of time;
2. (b) all or a part of the parameters of the differential equations are represented as a function of a growth rate factor.
3. (c) the differential equations include a formation rate for formation of a cell component from an intracellular metabolite, and the formation rate is represented as a function of a growth rate factor;
4. (d) the differential equations include an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, and the inflow rate and/or the outflow rate is represented as a growth rate factor;
(2) a means for storing values of parameters in the set of differential equations required for the simulation, and
(3) a means for calculating solutions of the set of differential equations based on the stored differential equations and values of the parameters.

A flowchart of a process executed by the program of the present invention is shown in Fig. 5. Moreover, a functional block diagram of a computer on which the program of the present invention is installed is shown in Fig. 6.
The means for storing a set of differential equations is constituted by a central processing portion 1, a storing portion 2 and an input portion 3. In a routine (S1) of storing a set of differential equations, the central processing portion 1 stores data of the set of differential equations inputted from the input portion 3 in the storing portion 2. The format of the data of the set of differential equations is not particularly limited, and it may be a usual format.
A means for storing values of parameters is constituted by the central processing portion 1, the storing portion 2 and the input portion 3. In a routine (S2) of storing values of parameters, the central processing portion 1 stores values of parameters inputted from the input portion 3 in the storing portion 2.
A means for computing solutions of the set of differential equations is constituted by the central processing portion 1, the storing portion 2 and an output portion 4. In a routine (S3) of computing the solutions of the set of differential equations, the central processing portion reads out the data of the differential equations and the values of parameters from the storing portion 2, computes the solutions from them and outputs the solutions to the output portion 4.
The central processing portion 1 is, for example, a processor. The storing portion 2 is, for example, a storage device using a recording medium. The input portion 3 is, for example, an input device such as keyboard and other devices or a data receiver for data from another device. The output portion 4 is, for example, an output device such as display, or a data transmission device for transmission to other devices.
The program for causing a computer to function as the aforementioned means can be created according to a usual programming method.
Further, the program of the present invention can also be stored in a computer-readable recording medium. The "recording medium" referred to herein include any "transportable physical media" such as floppy disk (registered trademark), magneto-optical disk, ROM, EPROM, EEPROM, CD-ROM, MO and DVD, any "physical media for fixation" such as ROM, RAM and HD built in various computer systems, and "communication media" storing the program over a short period of time such as communication cables and carrier waves in the case of transmitting the program through a network, of which typical examples are LAN, WAN and Internet.
Further, the "program" is a data processing method described with any language or description mode, and the format such as source code and binary code is not limited. In addition, a "program" is not necessarily' restricted to those configured as a single program, and includes those configured as a distributed system as two or more modules and libraries and those achieving the function through cooperation with another program, of which typical example is an operation system (OS). As specific configuration for reading from the recording medium in the devices shown in the embodiment, routines for reading, routines for installation after reading and so forth, known configurations and routines can be used.
The present invention provides accurate metabolic simulation results for a substance-production process that is accompanied by growth of cells. A simulation of the invention thus enables in silico experiments, which lend practical direction to improving such a production process, typically for obtaining an amino acid or nucleic acid. Illustrative of the improvement realizable in this regard is production optimization in batch culture or fed-batch culture often used as actual industrial process.

Hereafter, the present invention will be further explained with reference to examples.

<1> System parameters

Simulation of expression of the proteins of the enzymes and transcription factors from the genes and conversion of substances by the enzymatic reactions mentioned in the central metabolic map of E. coli shown in Fig. 7 was performed. The system parameters used for the simulation and the metabolites of which quantities were included as constants are shown in Table 1. In order to use experimental values of turbidity OD (optical density) as an index of growth, the following basic data were obtained. For dry cell weight per OD, DCWperOD, the results obtained in weight measurement of dry cells obtained from 300 ml of culture broth of MG1655 strain by centrifugal separation. As for cell density per OD, celldens, cell density of the culture broth of MG1655 strain subjected to two steps of dilution was measured 5 times, this procedure was independently repeated 5 times (measurement was performed for 50 plates in total), and the average of the results was used. Cell volume, cellvol, and cell weight, cellweight, per cell were calculated by considering that cellvol of per 1 g of DCW was 0.0025 l/g (Rohwer et al., J. Biol. Chem., 275, 34909-34921, 2000). The translation rate constant was calculated from a value mentioned in literature, 11.03 (min)^-1, at µ of 0.01 (min)^-1 (Lee and Baily, Biotech. Bioeng., 26, 66-72, 1984) and a ribosome concentration. As the rate constant of proteolysis, a value mentioned in literature was used either (Miller, C.G., In Escherichia coli and Salmonella: Cellular and Molecular Biology/Second Edition (Neidhardt, F.C. Ed.), pp.680-691, American Society for Microbiology Press, Washington, D.C., 1996).

Table 1: System parameters

The system parameters and values of the quantities of metabolites included as constants are shown. The values for which titles of literature are mentioned are literature values, and "Experimental" indicates that the values were obtained by experiments.
System parameter	Name	Value	Unit	Reference
DCWperOD	dry cell weight per OD	3.91E-01	g/L	Experimental
celldens	cell density	3.99E+12	cells/L	Experimental
reacvol	culture volume	3.00E-01	L	Experimental
initOD	initilal OD	2.07E+00		Experimental
cellvol	volume of single cell	2.00E-16	L	Calculated from gDCW/cellvol (0.0025 l/g)
cellweight	weight of single cell	8.00E-14	g	Calculated from gDCW/cellvol (0.0025 l/g)
ktrans	rate constant for translation	1.21E+05	(Mmin)^-1	Biotech. Bioeng. 26, 66-72, 1984
kdeg	rate constant for protein degradation	1.67E-04	min^-1	Escherichia coli and Salmonella 44, 680-691, 1996
ATP	adenosine-5-triphosphate	4.27E-03	M	Biotech Bioeng. 79, 53-73, 2002
ADP	adenosine-5-diphosphate	5.95E-04	M	Biotech. Bioeng. 79, 53-73, 2002
AMP	adenosine-5-monophosphate	9.55E-04	M	Biotech. Bioeng. 79, 53-73, 2002
NAD	nicotinamide adenine dinucleotide	1.47E-03	M	Biotech. Bioeng. 79, 53-73, 2002
NADH	nicotinamide adenine dinucleotide reduced	1.00E-04	M	Biotech. Bioeng. 79, 53-73, 2002
NADP	dihydronicotinamide adenine dinucleotide phosphate	1.95E-04	M	Biotech. Bioeng. 79, 53-73, 2002
NADPH	dihydronicotinamide adenine dinucleotide phosphate reduced	6.20E-05	M	Biotech. Bioeng. 79, 53-73, 2002
CoA	CoA	1.23E-04	M	Metab. Eng. 4, 182-192, 2002
Pi	inorganic phosphate	1.00E-02	M	Escherichia coli and Salmonella 87, 1357-1381, 1996
ASP	aspartate	1.34E-03	M	Biochem. J. 356, 433-444, 2001
CO2	carbon dioxide	1.35E-03	M	Experimental

<2> Modeling of enzymatic reaction

The abbreviations and initial values of the metabolites used in this example are shown in Table 2. As for those obtained from literature, titles of literature are shown. In the simulation, supposing that the volumes of all cells increased in the process of growth, the total cell volume (cellvoltot) was used as a variable. The initial value of the total cell volume (cellvoltot) was calculated from the initially measured value of OD (initOD) in accordance with the following equation. $cellvoltot = cellvol \times celldens \times reacvol \times initOD$

Modeling of the enzymatic reactions was performed based on the Michaelis-Menten type equation described by Segel (Segel, Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady-State Enzyme Systems, John Wiley & Sons, New York, 1975). Names and initial values of quantities of enzymes, transcription factors and mRNA are shown in Table 3. Types, values of parameters, substrates, products and effectors of the enzymatic reactions are shown in Table 4. The enzymatic reaction formulas are shown in Table 4.

Table 2: Abbreviations and initial values of metabolites The values for which no literature name is mentioned are estimated values.

Variable	Name	Initial value	Unit	Reference
Cellvoltot	total cell volume	4.96E-04	L
Glext	external glucose	2.20E-01	M	Experimental
G6P	glucose-6-phosphate	3.48E-03	M	Biotech. Bioeng. 79, 53-73, 2002
F6P	fructose-6-phosphate	6.00E-04	M	Biotech. Bioeng. 79, 53-73, 2002
FDP	fructose-16-diphosphate	2.72E-04	M	Biotech. Bioeng. 79, 53-73, 2002
GA3P	glyceraldehyde 3-phosphate	2.18E-04	M	Biotech. Bioeng. 79, 53-73, 2002
DHAP	dihydroxyacetone phosphate	1.67E-04	M	Biotech. Bioeng. 79, 53-73, 2002
13DPG	13-bis-phosphoglycerate	8.00E-06	M	Biotech. Bioeng. 79, 53-73, 2002
3PG	3-phosphoglycerate	2.50E-03	M	Nature 305, 286-290, 1983
2PG	2-phosphoglycerate	3.99E-04	M	Biotech. Bioeng. 79, 53-73, 2002
PEP	phosphoenolpyruvate	2.67E-03	M	Biotech. Bioeng. 79, 53-73, 2002
PYR	pyruvate	2.67E-03	M	Biotech. Bioeng. 79, 53-73,2002
6PGC	6-phosphogluconate	8.08E-04	M	Biotech. Bioeng. 79, 53-73, 2002
RL5P	ribulose-5-phsophate	1.11E-04	M	Biotech. Bioeng. 79, 53-73, 2002
R5P	ribose-5-phosphate	3.98E-04	M	Biotech. Bioeng. 79, 53-73, 2002
X5P	xylulose-5-phosphate	1.38E-04	M	Biotech. Bioeng. 79, 53-73, 2002
B4P	erythrose-4-phosphate	9.80E-05	M	Biotech. Bioeng. 79, 53-73, 2002
S7P	sedoheptulose-7-phosphate	2.76E-04	M	Biotech. Bioeng. 79, 53-73, 2002
ACCoA	acetylCoA	3.00E-04	M	Anal. Biochem. 295, 129-137, 2001
OAA	oxaloacetate	6.80E-04	M	Biomol. Eng. 19, 5-15, 2002
CIT	citrate	1.50E-04	M	Biomol. Eng. 19, 5-15, 2002
ICIT	isocitrate	1.70E-04	M	Biomol. Eng. 19, 5-15, 2002
AKG	2-ketoglutarate	1.80E-04	M	Biomol. Eng. 19, 5-15, 2002
SUCCoA	succinylCoA	1.00E-04	M
SUCC	succinate	1.90E-04	M	Biomol. Eng. 19, 5-15, 2002
FUM	fumarate	1.00E-04	M
MAL	malate	6.00E-05	M	Biomol. Eng. 19, 5-15, 2002
GLX	glyoxylate	1.00E-04	M
cAMP	cyclic AMP	8.00E-06	M	Mol. Microbiol. 10, 341-350, 1993
cAMPxt	external cyclic AMP	8.00E-06	M	Mol. Microbiol. 10, 341-350, 1993
FIP	fructose 1-phosphate	1.00E-04	M
CO2xt	external CO2	0.00E+00	M

<3> Equilibirium reactions

The algebraic equations were reduced for solutions by assuming that equilibirium is established for the binding of transcription factors and effectors and activation of CYA. CRP is a transcription factor involved in catabolite suppression, and [CRP-cAMP]^tot produced by equilibration of CRP and cAMP as an effector thereof was represented as a solution of the quadratic equation shown in the row of CRP1 in Table 6. [CRP]^tot and [cAMP]^tot represent the total concentrations of intracellular CRP and cAMP, respectively. About 200 binding sites on the genome are known for CRP, and it is necessary to consider that equilibirium is established for CRP-cAMP and these binding sites. Assuming that the dissociation constant for this, K_dCRP is 4 x 10^-8 (M), the concentration of CRP-cAMP binding with the promoter on the genome can be considered a solution of the quadratic equation shown in the row of CRP2 in Table 6. It is known that Mlc suppresses a target gene in the absence of glucose, whereas it binds to non-phosphorylated IICB^Glc in the presence of glucose. [Mlc-IICB^Glc] which binds with IICB^Glc and thus is inactivated was represented as a solution of the quadratic equation shown in the row of Mlc in Table 6. Cra is a transcription factor known as an activator/repressor of many sugar metabolism-related genes. It was assumed that the Cra concentration was constant. [Cra-F1P] binding to F1P, which is an effector thereof, was represented as a solution of the quadratic equation represented in the row of Cra in Table 6. The effector of PdhR, which is a repressor of the aceEF gene coding for PDH, is PYR (Quail and Guest, Mol. Microbiol., 15, 519-529, 1995), and [PdhR-PYR] which binds with PYR and is thus inactivated was represented as a solution of the quadratic equation shown in the row of PdhR in Table 6. It has been suggested by Cortay et al. that IclR is a transcription factor that suppresses expression of the glyoxylic acid pathway, and the effector thereof is PEP (Cortay et al., EMBO J., 10, 675-679, 1991). [IclR-PEP] which binds with PEP and thus is inactivated was represented as a solution of the quadratic equation shown in the row of PclR in Table 6. Activation of CYA by phosphorylated IIA^Glc (IIA^Glc-P) is known. Although the detailed mechanism is unknown, modeling was performed by assuming that CYA and IIA^Glc-P bind to each other to form an activated CYA (CYAA). Based on the difference in CYA activity between a wild type strain and a strain deficient in crr, which codes for IIA^Glc (Reddy and Kamireddi, J. Bacteriol., 180, 732-736, 1998), the dissociation constant of CYAA, K_dCYAA, was predicted as 1.34 x 10^-4 (M). Based on this, concentration of the activated CYAA [CYAA] produced by the equilibirium of CYA and IIA^Glc-P was represented as a solution of the quadratic equation shown in the row of CYA in Table 6.

<4> Modeling of gene expression

As for the gene expression of E. coli, modeling was performed for the genes to which the transcription factors CRP, Cra, Mlc, PdhR, and IclR relate by referring to the EcoCyc database (Keseler et al., Nucleic Acids Res., 33, D334-D337, 2005). As for the transcription factors (TF) per se, modeling of expression was performed for CRP, Mlc, PdhR and IclR. A list of the equations of transcription and translation and parameters of the genes is shown in Table 7, and the equations used for the gene expression are shown in Table 8. [mRNA_gene] represents the concentration of mRNA to be transcribed, [P_gene] represents the concentration of a promoter for a gene, and [RNAP□^D] represents the concentration of RNA polymerase bound with □^D. The parameters k_gene ^base, k_gene ^TF and k_gene ^dRNA are a baseline transcription rate constant, a transcription rate constant for TF-binding gene, and a decomposition constant of mRNA, respectively. [Protein] represents the concentration of translated protein, and [Ribosome] represents the ribosome concentration. The parameters k_trans and k_deg represent a translation rate constant and a proteolysis rate constant, respectively. The NoTF equation was used for a gene of which control is not known and a gene for which control is not considered, and TF1 was used for a gene to which one transcription factor relates. For the genes to which two transcription factors relate (crp, mlc, ptsG, ptsHI, aceBAK), equation was mentioned for each gene. If the concentrations of the translated proteins of the same operon are different, it is considered that such difference is due to reduction of transcription product or difference in transcription efficiency. Therefore, T_Protein (translation efficiency coefficient of protein) was defined, and the TL1 equation was used. The sdhCDAB-sucABCD operon of E. coli contains sdhCDAB coding for SDH, and sucABCD composed of sucAB coding for the E1 and E2 subunits of the KGDH complex and sucCD coding for SCS (Cunningham and Guest, Microbiology, 144, 2113-2123, 1998). In consideration of the fact that the sucABCD operon is also transcribed from P_suc besides P_sdh, and coefficients concerning the translation efficiency, T_KGDH and T_scs, for the KGDH complex and SCS, respectively, it was described by using TL (KGDH) and TL (SCS). Because it was reported that the ratio of the products of the aceBAK operon, ICL, MSA and ICDKP, is 1:0.3:0.003 (Chung et al., J. Bacteriol., 175, 4572-4575, 1993), translation constants T_MSA and T_ICDKP were defined for the translation of MSA and ICDKP, respectively, and TL1 was used.
As concentration of a promoter, a value at µ of 0.01 (min)^-1 was estimated based on the µ-dependent intracellular gene number data described by Bremer and Dennis (Escherichia coli and Salmonella: Cellular and Molecular Biology/Second Edition (Neidhardt F.C., Ed., pp.1553-1569, American Society for Microbiology Press, Washington, D.C., 1996), and used as a constant. As for the number of binding sites on the genome of CRP, a value at µ of 0.01 (min)^-1 was calculated on the presumption that about 200 of the binding sites of CRP are uniformly distributed over the genome. As the rate constant of transcription, a value calculated so that it should give the literature value of the protein concentration or specific activity of enzyme as a constant value was used. When two or more transcription rate constants were required in control by a transcription factor, they were calculated by using data of transcription activity or protein concentration in a transcription factor-deficient strain. As the mRNA decomposition rate, if any experimental value is available from literature for a certain gene, it was used, or otherwise, measurement data based on the DNA microarray experiment of Selinger et al. (Genome Res., 13, 216-223, 2003) were used.

<5> Preparation of mathematical equation for specific growth rate µ and cell formation rate

The specific growth rate µ is an index often used for representing growth. In order to represent growth with µ as accurately as possible, the following approximate equation of OD was obtained from OD data over time from culture of a wild type strain in a S-type jar by using a curve fitting program, TableCurve 2D (Systat Software), and a time function of µ was obtained from an equation obtained by differentiating the approximate equation. The result of plotting for OD and µ based on the approximate equation is shown in Fig. 2. $OD = (2.05 + 1.53 \times 10^{- 5} t^{2} - 5.35 \times {10 -}^{10} t^{4} + 3.07 \times 10^{- 25} t^{6}) / (1 - 1.76 \times 10^{- 5} t^{2} + 1.17 \times 10^{- 10} t^{4} - 3.19 \times 10^{- 16} t^{6} + 4.19 \times 10^{- 22} t^{8})$
$μ = dOD / dt / OD$
In order to compute the cell formation rate during the growth, metabolic reactions of E. coli were defined based on the report of Chassagnole et al. (Biotechnol. Bioeng., 79, 53-73, 2002). Synthetic reactions of the cell components were defined for each of protein synthesis, RNA synthesis, DNA replication, lipid synthesis, glycogen synthesis and peptidoglycan (murein) synthesis, and stoichiometric equations were defined from ratios of components (Neidhardt and Umbarger, Escherichia coli and Salmonella: Cellular and Molecular Biology (Neidhardt, F.C. Ed., pp.13-16, American Society for Microbiology and Washington D.C., 1996; Pramanik and Keasling, Biotechnol. Bioeng., 56, 398-421, 1997) and energy required for synthesis (Stephanopoulos et al., Metabolic Engineering: Principles and Methodologies, Academic Press, San Diego, 1998). Furthermore, a stoichiometric equation was prepared for each component required for synthesis of 1 g of cells from the composition of cells (Chassagnole et al., Biotechnol. Bioeng., 79, 53-73, 2002). This equation concerning the cell formation was converted into a stoichiometric equation using intermediate metabolites used in the simulation to create the following stoichiometric equation for each intermediate metabolite required for producing 1 g of cells. $g_{-} biomass = 3.962 Pyr + 1.229 aKG + - 2.232 {CO}_{2} + 10.91 {NH}_{4} + 44.69 ATP + - 44.6 ADP + - 15.18 P + 16.09 H + 18.17 NADPH + - 18.17 NADP + 2.409 AcCoA + - 2.949 CoA + - 0.487 Fum + 2.393 OAA + 1.957 3 PG + 0.252 SO 4 + - 2.329 NADP + 2.329 NAD + 0.5402 SucCoA + - 0.4727 Suc + 0.6887 PEP + 0.3312 E 4 P + 0.4133 DHAP + 0.1023 O 2 + - 0.0432 GAP + 0.5312 R 5 P + 0.1025 F 6 P$
The amount of the intermediate metabolite required for formation of 1 g of cells was converted into a value per cell volume and minute and incorporated into the differential equations as an equation of the specific growth rate µ.

<6> Preparation of mathematical equations of RNA polymerase and ribosome

Transcription and translation in gene expression are catalyzed by RNA polymerase and ribosome, respectively. It is known that the molecular numbers of these enzymes change during the process of growth. From the data of µ-dependent intracellular molecular number described by Bremer and Dennis (Escherichia coli and Salmonella: Cellular and Molecular Biology/Second Edition (Neidhardt F.C. Ed., pp.1553-1569, American Society for Microbiology Press, Washington, D.C., 1996), mathematical equations were prepared by using approximate equations. RNA polymerase binds with the □ factor to become a holoenzyme and then function. Because □^D responsible for gene expression during the growth phase is substantially constant during the growing process, □^D-bound RNA polymerase concentration [RNAP□^D] was considered to be 1/3 of the total RNA polymerase concentration, and represented by the following equation. $[RNAP □^{D}] = 6.67 \times 10^{- 7} + 3.0 \times 10^{- 4} μ + 2.64 \times 10^{- 2} μ^{2}$
As for ribosome, by fitting the data of Bremer and Dennis (Escherichia coli and Salmonella: Cellular and Molecular Biology/Second Edition (Neidhardt F.C. Ed., pp. 1553-1569, American Society for Microbiology Press, Washington, D.C., 1996) using TableCurve 2D, the following equation was obtained. $[Ribosome] = 1.90 \times 10^{- 5} + 1.38 μ^{2} + 10.2 μ^{2.5} + 36.6 μ^{3}$

<7> Excretion to outside of cells and uptake from outside of cells

Among the substances excreted to the outside of cells, excretion of acetic acid (AcOH) and formic acid (Formate), of which amounts detected in culture of a wild type strain were large, was incorporated. Based on the measured data for extracellular concentrations of acetic acid and formic acid over time, the profiles over time were approximated to time functions, and the functions were converted into rates by differentiation and incorporated into the differential equations of ACCoA and PYR. The plot of the amount of acetic acid based on the approximated function and the rate obtained by differentiating the approximated function is shown in Fig. 4. $\begin{matrix} {AcOH}_{ex} & = (2.49 \times 10^{- 3} - 7.61 \times 10^{- 3} t - 3.38 \times 10^{- 7} t^{2} + 9.33 \times 10^{- 10} t^{3}) / (1 - 761 \times 10^{- 3} t + 2.44 \times 10^{- 5} t^{2} - 1.05 \times 10^{- 8} t^{3}) \\ {Form}_{ex} & = 4.41 \times 10^{- 4} - 1.17 \times 10^{- 9} t^{2} + 1.97 \times 10^{- 13} t^{4} + 3.93 \times 10^{- 19} t^{6} \end{matrix}$
Among the organic substances existing in medium, amino acids and so forth are taken up into the cells. Uptake of glutamic acid (Glu) and alanine (Ala), of which existing amounts in culture of a wild type strain were large, was represented with mathematical equations and incorporated. Uptake of glutamic acid and alanine contained in the initial medium was approximated with the following time functions, and the functions were converted into rates by differentiation of the functions and incorporated into the differential equations of AKG and PYR. $G 1 u_{in} = 9.2 \times 10^{- 4} - 7.26 \times 10^{- 6} t$
$Alain = 8.36 \times 10^{- 4} - 6.69 \times 10^{- 6} t$

<8> Culture of wild type strain MG1655 and metabolic flux analysis

The wild type strain MG1655 was cultured overnight in 30 ml of LB medium, and the cell were collected from the culture broth. The cells were cultured in MS medium using ¹³C glucose contained in an S-type jar under the conditions of batch culture. As for the composition of the MS medium, it had a composition of 40 g of glucose, 1 g of MgSO₄•7H₂O, 16 g of (NH₄)₂SO₄, 1 g of KH₂PO₄, 2 g of Bacto yeast extract, 0.01 g of MnSO₄•4H₂O, 0.01 g of FeSO₄•7H₂O, and 0.5 ml of GD113 (antifoaming agent) in 1 L, and as for the culture conditions, the culture was carried out in a culture volume of 0.3 L, at a temperature of 37°C and pH 7.0 with aeration by stirring. Metabolic flux analysis was performed during the growth phase (315 minutes) and the stationary phase (495 minutes). The method of metabolic flux analysis is described in International Publication No. WO2005/001736 in detail. These culture results were used for verification of the simulation. Plots of the results of the measurements of extracellular glucose concentration and extracellular CO₂ concentration are shown in Fig. 8A, B and Fig. 8B, P, respectively (broken lines). Further, the results of the metabolic flux analysis during the growth phase (315 minutes) and the stationary phase (495 minutes) are shown in Table 9. The values of metabolic fluxes converted into enzymatic reaction rates are plotted to enzymatic activity (Fig. 9A, F and I, Fig. 9B, L, O and R).

Table 9: Results of metabolic flux analysis and conversion into enzymatic reaction rates

The metabolic fluxes are represented with values standardized in mmol/10 mmol Glc. The enzymatic activity is represented with values obtained by converting a sugar consumption rate at a corresponding time into actual enzymatic reaction rate (mol/min).
Enzyme	Flux at 315 min (mmol/10mmol)	Rate at 315 min (M/min)	Flux at 495 min (mmol/10mmol)	Rate at 495 min (M/min)
IICBGIc	10.00	0.0795	10.00	0.0497
PGI	6.02	0.0479	2.94	0.0146
PFKA	6.02	0.0478	7.05	0.0351
TPIA	6.02	0.0479	7.05	0.0351
GAPA	16.47	0.1309	16.47	0.0819
PGK	16.47	0.1309	16.47	0.0819
dGPM+iGPM	14.58	0.1159	14.87	0.0739
ENO	14.58	0.1159	14.87	0.0739
PYKF	0.71	0.0057	2.52	0.0125
PDH	7.24	0.0575	10.46	0.0520
G6PD	3.98	0.0316	7.05	0.0351
6PGD	3.98	0.0316	7.05	0.0351
RPIA+RPE	3.98	0.0316	7.05	0.0351
TKTAI+TALB	0.32	0.0026	0.22	0.0011
TKTAII	0.25	0.0020	1.62	0.0081
CS	4.02	0.0320	7.77	0.0386
ACNA+ACNB	4.02	0.0320	7.77	0.0386
ICDA	4.02	0.0320	7.21	0.0358
KGDH	2.87	0.0228	6.43	0.0320
SCS	2.87	0.0228	6.43	0.0320
SDH	93.10	0.7401	21.23	0.1056
FRD	90.00	0.7154	14.09	0.0701
FUMA	3.10	0.0246	7.14	0.0355
MDH	3.10	0.0246	7.14	0.0355
FBA	6.02	0.0479	7.05	0.0351
PPSA	0.00	0.0000	0.00	0.0000
PCK	0.01	0.0001	1.49	0.0074
PPC	2.81	0.0223	3.13	0.0155
NADPME+NADME	0.00	0.0000	0.30	0.0015
ICL	0.00	0.0000	0.56	0.0028
MSA	0.00	0.0000	0.56	0.0028

<9> Simulation and verification of E. coli central metabolic model

The simulation was performed by describing differential equations using a mathematical calculation program MATLAB (MathWorks) and using ode15s as an ODE solver. The differential equations of material balance used for the simulation are shown below. Material balance of each substance is described as the sum of enzymatic reaction rate, dilution effect due to growth, synthesis rate of cell component (y_biomass(metabolite)), uptake rate of extracellular substance and rate of excretion to the outside of cells mentioned in Table 4. $d [Cellvoltot] / dt = μ * [Cellvoltot]$
$d [Glucose] / dt = - rx 1 e$
$d [G 6 P] / dt = rx 1 e - rx 2 - rx 12 - (μ * [G 6 P])$
$d [F 6 P] / dt = rx 2 + rx 29 + rx 16 + rx 17 b + rx 3 - (μ * [F 6 P]) - y_{-} biomass (F 6 P) * mu / cellvol * cellweight$
$d [FDP] / dt = rx 3 - rx 4 - rx 29 - (μ * [FDP]) - y_{-} biomass (FDP) * μ / cellvol * cellweight$
$d [DHAP] / dt = rx 4 - rx 5 - (μ * [DHAP]) - y_{-} biomass (DHAP) * μ / cellvol * cellweight$
$d [GA 3 P] / dt = rx 4 + rx 5 + rx 17 b - rx 6 - rx 16 - rx 17 - (μ * [GA 3 P]) - y_{-} biomass (GA 3 P) * μ / cellvol * cellweight$
$d [13 DPG] / dt = rx 6 - rx 7 - (μ * [13 DPG])$
$d [3 PG] / dt = rx 7 - rx 8 - rx 9 - (μ * [3 PG]) - y_{-} biomass (3 PG) * μ / cellvol * cellweight$
$d [2 PG] / dt = rx 8 + rx 9 - rx 10 - (μ * [2 PG]) PEP$
$d [PEP] / dt = rx 10 + rx 30 + rx 34 - rx 11 - rx 1 a - rx 31 - (μ * [PEP]) - y_{-} biomass (11) * μ / cellvol * cellweight$
$d [PYR] / dt = rx 11 + rx 1 a + rx 32 + rx 33 - rx 18 - rx 30 - (μ * [PYR]) - y_{-} biomass (PYR) * mu / cellvol * cellweight + Alauptake - Formin$
$d [6 PGC] / dt = rx 12 - rx 13 - (μ * [6 PGC])$
$d [RL 5 P] / dt = rx 13 - rx 14 - rx 15 - (μ * [RL 5 P])$
$d [R 5 P] / dt = rx 14 + rx 17 - (μ * [R 5 P]) - y_{-} biomass (R 5 P) * μ / cellvol * cellweight$
$d [X 5 P] / dt = rx 15 + rx 17 - rx 17 b - (μ * [X 5 P])$
$d [E 4 P] / dt = rx 16 - rx 17 b - (μ * [E 4 P]) - y_{-} biomass (E 4 P) * μ / cellvol * cellweight$
$d [S 7 P] / dt = - rx 16 - rx 17 - (μ * [S 7 P])$
$d [ACCoA] / dt = rx 18 - rx 19 - rx 36 - (μ * [ACCoA]) - AcOHin - y_{-} biomass (19) * mu / cellvol * cellweight + Formin$
$d [OAA] / dt = rx 28 + rx 31 - rx 19 - rx 34 - (μ * [OAA]) - y_{-} biomass (OAA) * μ / cellvol * cellweight$
$d [CIT] / dt = rx 19 - rx 20 - rx 21 - (μ * [CIT]) - y_{-} biomass (CIT) * μ / cellvol * cellweight$
$d [ICIT] / dt = rx 20 + rx 21 - rx 22 - rx 35 - (μ * [ICIT]) - y_{-} biomass (ICIT) * μ / cellvol * cellweight$
$d [AKG] / dt = rx 22 - rx 23 - (μ * [AKG]) - y_{-} biomass (AKG) * μ / cellvol * cellweight + Gluuptake$
$d [SUCCoA] / dt = rx 23 - rx 24 - (μ * [SUCCoA]) - y_{-} biomass (SUCCoA) * μ / cellvol * cellweight$
$d [SUCC] / dt = rx 24 + rx 35 - rx 25 - rx 26 - (μ * [SUCC]) - y_{-} biomass (SUCC) * μ / cellvol * cellweight$
$d [FUM] / dt = rx 25 + rx 26 - rx 27 - (μ * [FUM]) - y_{-} biomass (FUM) * μ / cellvol * cellweight$
$d [MAL] / dt = rx 27 + rx 36 - rx 28 - rx 32 - rx 33 - (μ * [MAL]) - y_{-} biomass (27) * μ / cellvol * cellweight$
$d [GLX] / dt = rx 35 - rx 36 - (μ * [GLX])$
$d [cAMP] / dt = rx 39 + rx 39 a - rx 40 - rx 41 - (μ * [cAMP])$
$d [cAMPex] / dt = rx 41$
$d [F 1 P] / dt = - rx 42 - (μ * [F 1 P])$
$d [CO 2] / dt = rx 13 + rx 18 + rx 22 + rx 23 + rx 32 + rx 33 + rx 34 - rx 31 - (μ * [CO 2]) - y_{-} biomass (32) * μ / cellvol * cellweight$
During the simulation, a part of the parameters were manually changed to perform the simulation. The changed parameters are shown in Tables 3, 4 and 7. Among the results of the simulation of the E. coli central metabolic model, temporal changes of major metabolites are shown in Figs. 8A and 8B. In the process of growth, the extracellular glucose concentration (Fig. 8A, B) and the extracellular CO₂ concentration (Fig. 8B, P) showed behaviors extremely close to those of the measured values (broken lines), although deviation is seen for the first half ot the culture. Among the metabolites, metabolites other than PEP (Fig. 8A, F) showed temporal changes in the culture in the presence of glucose within the numerical ranges considered physiologically reasonable. The mRNA concentrations of major genes, the protein concentrations and the enzymatic activities are shown in Figs. 9A and 9B. Since all the initial values of mRNA concentrations were set at 0, they showed common patterns that they increased with increase in the growth rate, and decreased after µ reached the maximum (Fig. 9A, A etc.). The protein concentrations showed patterns that they decreased from the initial values, then increased with the increase in the growth rate, reached the maximum slightly behind the peak of mRNA, and decreased thereafter (Fig. 9A, B etc.). As for enzymatic activities, whereas a significant deviation from the results of the metabolic flux analysis was observed for the activity of the oxidative pentose phosphate pathway enzyme, G6PD (Fig. 9B, L), profiles close to those of the values of the enzymatic reaction rates based on the flux analysis data could be obtained for the activity of PFK in the glycolysis system (Fig. 9A, F), activity of the TCA cycle enzyme, CS (Fig. 9B, 0), and activity of PEPC of the supplemental pathway (Fig. 9B, R).
The results of simulation of gene expression where RNA polymerase and ribosome concentrations were independent from µ are shown in Figs. 10A and 10B. In this simulation, the values of the RNA polymerase and ribosome concentrations were those observed at µ of 0.01 (min)^-1, and only k_ptsG ^CRP among the parameters was changed to 0.45 time of the original value to perform the simulation. The results were physiologically unreasonable, for example, the mRNA level became constant (Fig. 10A, A etc.), the protein concentration increased in the second half of the culture where µ decreased (Fig. 10A, B etc.), and so forth. These results indicate that description of µ-dependent RNA polymerase and ribosome concentrations is important for carrying out the simulation of the growing process. In order to investigate the effect of the cell formation accompanying the cell growth, simulation was performed with synthesis rates from intracellular metabolites to cell components of 0, and the results of the metabolic simulation are shown in Figs. 11A and 11B. It can be seen that many metabolites were accumulated in the cells, and the concentrations thereof are deviated from the physiological concentrations. Furthermore, in order to investigate effects of uptake and excretion of substances, simulation was performed with Glu and Ala uptake of 0 and acetic acid and formic acid excretion of 0, and the results of the metabolites are shown (Figs. 12A and 12B). It can be seen that evident differences were observed in the simulation results, such as those of AcCoA (Fig. 12B, J) and CIT (Fig. 12B, K), and they were deviated from the physiological concentrations. Thus, it was revealed that formation of cell components from intracellular metabolites, uptake and excretion of extracellular metabolites are necessary for realizing simulation similar to measurement results.

Claims

A method for effecting a simulation of a substance-production process that uses cells, wherein said simulation is based on a set of differential equations that represent intracellular metabolites and gene expression, said method comprising the steps of:
(a) including a specific growth rate of cells, expressed as a differential equation, in the set of differential equations;

(b) assigning values for parameters in the set of differential equations, wherein at least one of said parameters is represented as a growth rate factor;

(c) incorporating in the set of differential equations a formation rate for formation of a cell component from an intracellular metabolite, wherein said formation rate is represented as a growth rate factor;

(d) incorporating in the set of differential equations an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, wherein said inflow rate and said outflow rate are represented, respectively, as a growth rate factor;

(e) solving the set of differential equations; and

(f) generating data representative of the substance-production process.
The method according to claim 1, wherein the differential equations include the following equations (1) to (3): $d [Metabolite] / dt = V_{input} - V_{output} - μ [Metabolite]$
$d [mRNA] / dt = k_{transcription} [P] - (k_{dRNA} + μ) [mRNA]$
and $d [Protein] / dt = k_{translation} [mRNA] - (k_{dProtein} + μ) [Protein]$

wherein, in the Equation 1, [Metabolite] represents an intracellular concentration of a metabolite, V_input represents the sum of rates of reactions producing the metabolite, V_output represents the sum of rates of reactions consuming the metabolite, and µ represents the specific growth rate;
in the Equation 2, [mRNA] represents a concentration of mRNA, k_{transcription} represents a rate constant of transcription, [P] represents a promoter concentration, k_dRNA represents a rate constant of decomposition of mRNA, and µ represents the specific growth rate, and
in the Equation 3, [Protein] represents a concentration of a protein, k_translation represents a rate constant of translation, k_dProtein represents a rate constant of decomposition of the protein, and µ represents the specific growth rate.
The method according to claim 1 or 2, wherein the growth rate factor is a function of the specific growth rate or a function of time.
The method according to any one of claims 1 to 3, wherein the specific growth rate is represented as a function of time, and the function is obtained by generating a mathematic equation from measurement data of the specific growth rate in the production process.
The method according to any one of claims 1 to 4, wherein the growth rate factor representing the formation rate is obtained by generating a mathematic equation expressing measurement data of the formation rate in the production process.
The method according to any one of claims 1 to 5, wherein the growth rate factor representing the inflow rate and/or the outflow rate is obtained by generating a mathematic equation expressing measurement data of the inflow rate and/or the outflow rate in the production process.
The method according to any one of claims 1 to 6, wherein the metabolite taken up into the cells is a substrate and/or an organic substance in a medium.
The method according to any one of claims 1 to 7, wherein the metabolite excreted out of the cells is an objective substance and/or a by-product.
The method according to claim 8, wherein the metabolite excreted out of the cells is an amino acid, an organic acid and/or carbon dioxide.
The method according to claim 9, wherein the metabolite excreted out of the cells is an amino acid or an organic acid.
The method according to any one of claims 1 to 10, wherein at least one of said parameters is a rate constant of transcription and/or a rate constant of translation.
The method according to any one of claims 1 to 11, wherein the cells are those of a microorganism having an amino acid producing ability and/or an organic acid producing ability.
The method according to claim 12, wherein the microorganism is Escherichia coli.
The method according to any one of claims 1 to 13, wherein a composition of the cells, are represented by a mathematical equation using the specific growth rate of the cells or the cells' equivalent index concerning the growth.
A computer program product for effecting a simulation of a substance-production process that uses cells, wherein said simulation is based on a set of differential equations that represent intracellular metabolites and gene expression, comprising:
(a) computer code for including a specific growth rate of cells, expressed as a differential equation, in the set of differential equations;

(b) computer code for assigning values for parameters in the set of differential equations, wherein at least one of said parameters is represented as a growth rate factor;

(c) computer code for incorporating in the set of differential equations a formation rate for formation of a cell component from an intracellular metabolite, wherein said formation rate is represented as a growth rate factor;

(d) computer code for incorporating in the set of differential equations an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, wherein said inflow rate and said outflow rate are represented, respectively, as a growth rate factor;

(e) computer code for solving the set of differential equations; and

(f) computer code for generating data representative of the substance-production process.
A system for effecting a simulation of a substance-production process that uses cells, wherein said simulation is based on a set of differential equations that represent intracellular metabolites and gene expression, comprising:
a processor for processing information; and

a storing means, including:
(a) computer code for including a specific growth rate of cells, expressed as a differential equation, in the set of differential equations;

(b) computer code for assigning values for parameters in the set of differential equations, wherein at least one of said parameters is represented as a growth rate factor;

(c) computer code for incorporating in the set of differential equations a formation rate for formation of a cell component from an intracellular metabolite, wherein said formation rate is represented as a growth rate factor;

(d) computer code for incorporating in the set of differential equations an inflow rate of a metabolite taken up from the outside of the cells and/or an outflow rate of a metabolite excreted out of the cells from the inside of the cells, wherein said inflow rate and said outflow rate are represented, respectively, as a growth rate factor;

(e) computer code for solving the set of differential equations; and

(f) computer code for generating data representative of the substance-production process.