Sage Journals: Discover world-class research

Abstract

Measuring the concentrations of metabolites and estimating the reaction rates of each reaction step consisting of metabolic pathways are significant for an improvement in microorganisms used in maximizing the production of materials. Although the reaction pathway must be identified for such an improvement, doing so is not easy. Numerous reaction steps have been reported; however, the actual reaction steps activated vary or change according to the conditions. Furthermore, to build mathematical models for a dynamical analysis, the reaction mechanisms and parameter values must be known; however, to date, sufficient information has yet to be published for many cases. In addition, experimental observations are expensive. A new mathematical approach that is applicable to small sample data, and that requires no detailed reaction information, is strongly needed. S-system is one such model that can use smaller samples than other ordinary differential equation models. We propose a simplified S-system to apply minimal quantities of samples for a dynamic analysis of the metabolic pathways. We applied the model to the phenyl lactate production pathway of Escherichia coli. The model obtained suggests that actually activated reaction steps and feedback are inhibitions within the pathway.

Keywords

Metabolic pathway pathway dynamics small sample S-system Michaelis-Menten

Introduction

Mathematical model for metabolic pathways

The artificial and industrial uses of microorganisms for material production have a long history of more than a thousand years. Recently, genetic operations have been widely applied to improve production. Two generally considered approaches introduce enzymes that have higher activities from other organisms or species and introduce enzymes to realize metabolic pathways that do not naturally occur in the microorganisms. The former method is popular because its operation is simpler and improvements are more predictable than with the latter method. Conventional gene modifications using ultraviolet or other radiation types are easy to achieve and have been widely applied in many industries. Nevertheless, the efficiency of such improvements is quite low because gene modifications occur accidentally and uncontrollably, and progress is made serendipitously. Therefore, gene introduction is currently used along with conventional methods.

Target genes for modification are chosen based on information including the reaction rates of the respective reaction steps within the metabolic pathway and include the production materials and substrates of microorganisms, as well as changes in the reaction rates through changes in the concentrations of the metabolites that consist of the pathway. Bottleneck reaction steps and feedback loop inhibitions are suggested based on this information. The genes of enzymes used in such reactions are candidates for modification.

The rates of enzymic reactions are generally defined as the limit of changes in the compounds over time.^1,2 Several formulae are established for the types of enzymic reactions, such as inhibition schemes. The most popular is the Michaelis-Menten law^3,4 for a simple one-to-one enzymic reaction without any inhibition or catalysis using an enzyme. The reaction rate is modeled based on the following reaction scheme:

S + E ⇌ E S \to P + E

(1)

in which $S$ indicates the substrate of the reaction, $E$ signifies the enzyme, and $P$ denotes the product. A bidirectional arrow represents a reversible reaction. A 1-directional arrow signifies a 1-way reaction. The reaction rate is modeled as the following ordinary differential equation (ODE):

\frac{d [P]}{d t} = \frac{V_{\max} [S]}{K_{m} + [S]}

(2)

Therein, a pair of square brackets denotes the concentration of the compound, $t$ is the time, and $V_{\max}$ and $K_{m}$ are parameters that define the kinetic characteristic of the enzyme. All reaction steps in a metabolic pathway can be represented using the ODE above if all reactions are simple enzymic reactions and if the parameter values are defined. The pathway is then modeled as simultaneous ODEs or the ODE system. Consequently, simulations can be done of concentration changes of metabolites, stability analysis, steady-state estimation, and bottleneck finding.^5,6

Generally, finding the parameter values is difficult and expensive because it requires enzyme isolation and measurements of the reaction rates in test tubes (in vitro measurements). Although the amount of enzyme information in the literature and public databases is growing,⁷ the values of the kinetic parameters have not been sufficiently published or accumulated. An enzyme generally has different parameter values for the various conditions and species of organisms. Moreover, the parameter values generally differ between in vitro and in vivo (in living cells) conditions.^8,9 For most industrial applications, a dynamical analysis of the pathways must be conducted without reaction scheme information or kinetic parameter values.

Canonical ODE model

The ODE systems in canonical forms are applicable because they are independent of the molecular mechanism of the reaction scheme. The S-system⁶ is one such canonical ODE model. For a reaction scheme with 2 reactions, the following is used:

X_{1} + X_{2} + \dots + X_{l} \to P \to X_{l + 1} + X_{l + 2} + \dots + X_{m}

(3)

where $X_{1}$ , $X_{2} \dots$ , and $X_{l}$ and $P$ denote the substrates and product of the first reaction, respectively, and $P$ and $X_{l + 1}, \dots, X_{m}$ denote the substrate and products of the second reaction, respectively. The S-system form is represented as follows:

\frac{d [P]}{d t} = α \prod_{j = 1}^{n} {[X_{j}]}^{g_{j}} - β \prod_{j = 1}^{n} {[X_{j}]}^{h_{j}}

(4)

where $X_{j}$ indicates the concentration of the metabolite $j$ , $g_{j}$ signifies the kinetic parameter representing the influence of $X_{j}$ to the increasing processes of $P$ , $h_{j}$ denotes the influence of $X_{j}$ on the decreasing processes of $P$ , and $α$ and $β$ are rate constants. The first term of the left-hand side of the equation represents the total rate of the increasing or composing processes of $P$ . The second term is the total rate of the decreasing or decomposing processes. In addition, $l$ and $m$ in the reaction (3) are the numbers of composing and decomposing processes of $P$ , respectively. All $l$ and $m$ compounds are suffixed sequentially in Equation (4).

The S-system above is a simplified form of the general mass action law,⁶ which can be presented as follows:

\frac{d [P]}{d t} = \sum_{i = 1}^{m} α_{i} \prod_{j = 1}^{n} {[X_{i j}]}^{g_{i j}}

(5)

Therein, the $i$ suffix denotes each composing and decomposing reaction of $P$ . A simplification of the S-system summarizes the composing reactions of $P$ into a term using $α$ and the decomposing reactions into a $β$ term in Equation (4). The parameters $g_{j}$ and $h_{j}$ correspond to reaction orders in the mass action law, $g_{i j}$ , respectively, and indicate dependencies between the metabolites $P$ and $X_{j}$ . Consequently, they represent the network scheme of the reaction pathway. No direct dependence exists between $P$ and $X_{j}$ when $g_{j}$ and $h_{j}$ in Equation (4) are equal to 0. $X_{j}$ suppresses the production of $P$ when $g_{j}$ is negative.

The parameter values of $g_{j}$ , $h_{j}$ , $α$ , and $β$ can be estimated using numerical optimization methods for finding the parameter values by which the calculated time series of $P$ through a numerical integration of Equation (4) matches the observed time series of the concentration of $P$ . The determined values of $g_{j}$ and $h_{j}$ might be considered as representing orders of each reaction between $X_{j}$ and $P$ . However, parameter optimization is an inverse problem¹⁰ because several sets of different parameter values are generally found for the given observed time series data. Restrictions and limitations are effective for difficulties such as fixing some $g_{j}$ to 0 based on biological knowledge.

Method for small sample

Numerical optimizations require a sufficient number of observed samples. Smaller needs are better because observations entail a certain amount of costs. A mathematical model with fewer parameters requires fewer samples. We propose a canonical ODE model for small samples through the simplification of the second term of Equation (4), as shown below:

\frac{d [P]}{d t} = α \prod_{j = 1}^{n} {[X_{j}]}^{g_{j}} - β {[P]}^{h}

(6)

By the mass action law, the decomposition rate of a compound depends solely on its concentration in many biological processes, such as Michaelis-Menten–type reactions shown in Equation (2). We introduce this idea as an assumption in Equation (6). Although decomposition reactions are often modeled as linear ODE, such as

\frac{d [P]}{d t} = A \exp (- [P] t)

our model includes a nonlinear decomposition term because we suspect that a linear term might be too simple for the metabolite in a complex biological network system containing many unknown reactions. Our assumption is reasonable when regulation of the degradation processes by these unknown reactions is not significant.

The resulting time series of the model based on a numerical integration varies greatly through a change in the initial value. Finding the best initial value is difficult because the observed initial value frequently has errors, particularly for small sample data sets. Therefore, we compare the model and data in differential spaces. The differential of the observed data can be calculated through a numerical differentiation, and the parameter values can be evaluated by comparing the differential values with the values of Equation (6).

We evaluated the proposed method according to its application to the phenyl lactate (PL) production pathway from glucose using Escherichia coli. For a pathway that includes branches and feedback loops, we estimated the actually activated reaction steps and activities of the feedback inhibitions suggesting strategies for an improvement in production.

Method

First, we build a pathway map based on information from the literature and different databases. We then choose some metabolites in the map as observation targets. A pathway map is reconstructed using only the target metabolites. The observations are measurements of the target metabolite concentration at sampling time points with equal intervals. A simplified S-system model is defined based on the reconstructed pathway map by fixing some $g_{j}$ to 0 if the link to $P$ from $X_{j}$ does not exist in the reconstructed pathway map.

Observed time series data on the concentration of metabolites are numerically differentiated. Optimal values of $g_{j}$ , $h$ , $α$ , and $β$ are sought using a nonlinear numerical optimization method such as a genetic algorithm¹⁰ or differential evolution¹¹ to minimize the difference between differentials of the observed data through a numerical differentiation and those from Equation (6). The optimal $g_{j}$ is considered, which represents the activity of the reaction from $X_{j}$ to $P$ . A creation of the formula in Equation (6) and an optimization of the parameters are conducted for each target metabolite. We choose the PL production metabolic pathway^12–15 (Figure 1) as the application target and 6 metabolites for the observation target. Then, we reconstruct the pathway using only the observation target metabolites (Figure 2). Six sampling time points are used. Changes in the concentration of the metabolites are observed in the log phase of the cell growth (Figure 3).

Figure 1.

Main pathway of phenyl lactate metabolism of Escherichia coli. A dashed line indicates that the path consists of plural reactions. A solid 1-directional line indicates a single enzymic reaction. A bidirectional line indicates 2 reactions: forward and backward processes catalyzed by the same or different enzymes. A double-lined circle indicates a target metabolite selected for observation.

Figure 2.

Reconstructed pathway consisting of observation target metabolites.

Figure 3.

Observed time series of concentrations of target metabolites at the 6 time points. SA, PP, PA, and PL indicate shikimic acid, phenyl pyruvate, phenylalanine, and PL, respectively.

Pathway map construction and observation

The pathway from glucose to PL consists roughly of glycolysis and shikimic acid (SA) pathways. Phenyl lactate is produced from phenyl pyruvate (PP), PP is from phenylalanine (PA), and PA is from prephenic acid. In addition, prephenic acid is from chorismic acid in the SA pathway. Choosing phosphoenolpyruvate (PEP), erythrose 4-phosphate (E4P), SA, PA, PP, and PL as the observation targets, we then reconstruct the pathway for these 6 metabolites (Figure 2). The pathway includes a branching point and 2 feedback loops.

ODE models for the respective metabolite

Changes in concentration of metabolites in the reconstructed pathway are modeled mathematically using the following ODE models:

\begin{array}{l} \frac{d [S A]}{d t} = α_{S A} {[P E P]}^{g_{S A, P E P}} {[E 4 P]}^{g_{S A, E 4 P}} \\ - β_{S A} {[S A]}^{h_{S A}} \\ \frac{d [P P]}{d t} = α_{P P} {[S A]}^{g_{P P, S A}} {[P A]}^{g_{P P, P A}} {[P L]}^{g_{P P, P L}} \\ - β_{P P} {[P P]}^{h_{P P}} \\ \frac{d [P A]}{d t} = α_{P A} {[S A]}^{g_{P L, S A}} {[P P]}^{g_{P A, P P}} - β_{P A} {[P A]}^{h_{P A}} \\ \frac{d [P L]}{d t} = α_{P L} {[P P]}^{g_{P L, P P}} - β_{P L} {[P L]}^{h_{P L}} \end{array}

(7)

Actually, production of PEP and E4P is not controlled by any other metabolites, and these compounds are independent variables in the ODE system above. Their respective dynamics are not modeled. Parameter $g_{P P, S A}$ in the ODE system represents the summarized actual activity of the reaction chain to PP that is produced from SA, consisting of several reaction steps. Here, $g_{a, b}$ indicates that the rate of change in the concentration of $a$ is increased by $b$ when the sign of $g_{a, b}$ is positive. A negative value of $g_{a, b}$ means that $b$ suppresses the composition processes of $a$ . A larger $h$ signifies a higher rate of decomposition or consumption of the metabolite. A negative $h$ means that the metabolite suppresses the decomposition itself. Therein, $α$ and $β$ are fixed rate coefficients. For each metabolite, the activities of the respective reaction steps of the composition of the metabolite can be compared.

Time differential values of the metabolite concentration are calculable using the ODE system by determining all parameter values. We introduce the differential evolution algorithm¹¹ to find the parameter values that minimize the differences between the differential of the concentration values calculated using the ODE system, as shown below:

E = \sum_{t = 1}^{T} (D_{c t} - D_{o t})

(8)

Therein, $E$ indicates the summarized differences, $D_{c t}$ signifies the differentials calculated by the ODE system with the parameter values at time point $t$ , and $D_{o t}$ denotes the differential at time point t of the observed data. Here, $D_{o t}$ is calculated using the quadratic interpolation of the observed data, and $E$ is minimized by searching the parameter values of $g$ , $h$ , $α$ , and $β$ for each metabolite. The optimization algorithm introduces uniformly distributed random numbers for the initial values of the search parameters.

Result

Model interpretation

The ODE system with the determined parameter values is shown below (Figure 4):

\begin{matrix} \frac{d [S A]}{d t} = 21.1 {[P E P]}^{65.4} {[E 4 P]}^{33.9} - 99.2 {[S A]}^{- 0.0674} \\ \frac{d [P P]}{d t} = 9.38 {[S A]}^{- 9.44} {[P A]}^{8.77} {[P L]}^{- 0.390} \\ - 57.2 {[P P]}^{- 1.95} \\ \frac{d [P A]}{d t} = 0.000223 {[S A]}^{2.12} {[P P]}^{- 0.307} \\ - 90.7 {[P A]}^{0.207} \\ \frac{d [P L]}{d t} = 5.60 {[P P]}^{1.13} - 1.51 \times 10^{- 6} {[P L]}^{2.82} \end{matrix}

Figure 4.

Estimated actual reaction activities on the reconstructed pathway. Numerical values shown in the figure are $g_{j}$ in Equation (6), which represent incoming links for each metabolite. SA, PP, PA, and PL indicate shikimic acid, phenyl pyruvate, phenylalanine, and PL, respectively.

The ODE system suggests that the reaction chain to SA from PEP is more active than that from E4P. In addition, PP production is inhibited by feedback loops from both PA and PL.

There are 3 negative value parameters. Reaction steps that are represented by these parameters consist of only productive enzymic reactions in the literature.^12,13 These negative parameters may suggest unknown regulatory pathways whose effects look substrate inhibitions.

Discussion

Numerical advantages of our model

The changes in metabolite concentration over time can be modeled using the Michaelis-Menten law, as shown below:

\begin{matrix} \frac{d [S A]}{d t} = \frac{V_{S A, P E P} [P E P]}{K_{S A, P E P} + [P E P]} + \frac{V_{S A, E 4 P} [E 4 P]}{K_{S A, E 4 P} + [E 4 P]} \\ - \frac{V_{P P, S A} [S A]}{K_{P P, S A} + [S A]} - \frac{V_{P A, S A} [S A]}{K_{P A, S A} + [S A]} \\ \frac{d [P P]}{d t} = \frac{V_{P P, S A} [S A]}{K_{P P, S A} + [S A]} + \frac{V_{P P, P A} [P A]}{K_{P P, P A} + [P A]} \\ + \frac{V_{P P, P L} [P L]}{K_{P P, P L} + [P L]} \\ - \frac{V_{P A, P P} [P P]}{K_{P A, P P} + [P P]} - \frac{V_{P L, P P} [P P]}{K_{P L, P P} + [P P]} \\ \frac{d [P A]}{d t} = \frac{V_{P A, S A} [S A]}{K_{P A, S A} + [S A]} + \frac{V_{P A, P P} [P P]}{K_{P A, P P} + [P P]} \\ - \frac{V_{P P, P A} [P A]}{K_{P P, P A} + [P A]} \\ \frac{d [P L]}{d t} = \frac{V_{P L, P P} [P P]}{K_{P L, P P} + [P P]} - \frac{V_{P P, P L} [P L]}{K_{P P, P L} + [P L]} \end{matrix}

For each metabolite, there are 8 parameters for SA, 10 for PP, 6 for PA, and 4 for PL because a model of the reaction step (single enzymic reaction) has 2 parameters. Changes of concentrations in time of SA and PP cannot be modeled using the 6 sampling data in this study. There are a total of 16 parameters for the reconstructed pathway because some parameters are common (1 of the outgoing reactions of SA is an incoming reaction of PP). Some numerical optimization methods can search for the 16 parameters simultaneously. However, such a simultaneous nonlinear numerical optimization is not easy, and the difficulty increases significantly with an increase in the number of parameters (“curse of dimensionality”).

The form of the S-system model is reduction or simplification of general mass action (GMA) law model shown as follows:

\begin{matrix} \frac{d [S A]}{d t} = α_{S A} {[P E P]}^{g_{S A, P E P}} {[E 4 P]}^{g_{S A, E 4 P}} \\ - β_{1 S A} {[S A]}^{g_{P P, S A} - - β_{2 S A} g_{P A, S A}} \\ \frac{d [P P]}{d t} = α_{P P} {[S A]}^{g_{P P, S A}} {[P A]}^{g_{P P, P A}} {[P L]}^{g_{P P, P L}} \\ - β_{1 P P} {[P P]}^{g_{P A, P P} - β_{2 P P} g_{P L, P P}} \\ \frac{d [P A]}{d t} = α_{P A} {[S A]}^{g_{P L, S A}} {[P P]}^{g_{P A, P P}} - β_{P A} {[P A]}^{g_{P P, P A}} \\ \frac{d [P L]}{d t} = α_{P L} {[P P]}^{g_{P L, P P}} - β_{P L} {[P L]}^{g_{P P, P L}} \end{matrix}

For each metabolite, there are 7 S-system parameters for SA, 8 for PP, 5 for PA, and 4 for PL, or 18 in total because some of the parameters are the same. The number of parametes are reduced for SA, PP, and PA compared with Michaelis-Menten models; however, concentrations of SA and PP cannot be modeled even when using the GMA. In our proposed model, the corresponding numbers are 5, 6, 5, and 4, or 20 in all, which is not a small number; however, it is not a problem that the total number of parameters is larger than in the S-system because no parameters are common to any 2 metabolites, and the parameters of each metabolite are optimized independently of the other metabolites. Our model has fewer parameters for each metabolite, which indicates that our model is more robust against errors than the S-system and Michaelis-Menten models.

The less number of model parameters means that the model needs a less number of data. Michaelis-Menten model consists of 2 parameters for each reaction, thus $2 n$ parameter values must be determined for a system that consists of $n$ compounds, whereas the simplified S-system contains $n + 2$ . This means that the simplified S-system model needs about a half quantity of data compared with a Michaelis-Menten model.

Although the model and data do not match perfectly because the data generally include errors in the probabilistic distributions, the optimization precision of the parameters (model fitness to the data) shown in Figure 5 is apparently sufficient. In terms of biochemical engineering, the values of some parameters shown in Figure 4 are large as the reaction order. Perhaps because of fluxes of pathways other than the reconstructed pathway from which we omitted pathways other than the main reaction chain, glycolysis has many branches to other subsystems. However, the parameter values for PA (PA in Figure 4), PP, and PL might be more reliable because the reaction steps occurring naturally around these metabolites may be nearly the same as those of the reconstructed pathway.

Figure 5.

Values of differentials calculated from the observed data through a quadratic interpolation and calculated from the optimized ODE model (simplified S-system model). The dashed line shows differential values calculated from the ODE model with the optimized parameters. A solid line represents the numerical differentiation from the observed data. SA, PP, PA, and PL indicate shikimic acid, phenyl pyruvate, phenylalanine, and phenyl lactate, respectively. Michaelis-Menten model and the full S-system model cannot be obtained due to the much number of parameters of the models compared with the number of data. SA, PP, PA, and PL indicate shikimic acid, phenyl pyruvate, phenylalanine, and PL, respectively.

The parameter values of the $α$ term of a metabolite (of the incoming link of a metabolite in the pathway map) are directly comparable. Phenyl pyruvate has 3 parameters to compare, 2 of which are negative, and 1 of which is positive. Therefore, PP production, which is inhibited by SA and PL (the final product), depends mainly on PA.

Two negative feedback loops exist: PP production is inhibited by PL, and PA production is inhibited by PP. Although it can be readily imagined that the inhibition of the feedback reactions increases the production of PL, the main inhibitory effect on PP production is from the SA. Disrupting one or more genes of the reactions to PP from SA might improve the PL production.

Phenylalanine decreases gradually (Figure 3), but shows no natural decomposition, which is represented as $A \exp (- B t)$ , and might be caused by the incoming link from PP.

In conclusion, the results show that the reliability of the estimated parameter values might not be the best or even very high because the reconstructed pathway and the ODE system are simplified. These values suggest that the target genes can be modified for an industrial improvement in production using microorganisms. This case study presents several suggestions that may be useful when constrained to only a few samples or low observation costs.

Footnotes

Acknowledgements

The stock of microorganisms used in this study was provided by professor Naoki Takaya at the Graduate School of Life and Environmental Science, Tsukuba University.

Funding:

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the New Energy and Industrial Technology Development Organization, Japan.

Declaration of conflicting interests:

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Author Contributions

DT and SA conceived the mathematical analyses. DT carried out the mathematical analyses. HK and TH conceived the experimental design for observations. HK, TH and CO prepared biological and biochemical materials. HK and YH carried out the experiments and observations.

References

Heinrich

Schuster

The modeling of metabolic systems. Structure, control, and optimality. Biosystems. 1998;47:61–77.

Gavalas

GR.

Nonlinear Differential Equations of Chemically Reacting Systems. Berlin, Germany: Springer-Verlag; 1968.

Bisswanger

Enzyme Kinetics: Principles and Methods. 2nd ed. Weinheim, Germany: Wiley-VCH; 2008.

Cleland

WW.

Enzyme kinetics. Ann Rev Biochem. 1967;36:77–112.

Savageau

MA.

Biochemical Systems Analysis: A Study of Function and Design in Molecular Biology. Boston, MA: Addison-Wesley; 1977.

Voit

EO.

Computational Analysis of Biochemical Systems. Cambridge, UK: Cambridge University Press; 2000.

Placzek

Schomburg

Chang

et al . BRENDA in 2017: new perspectives and new tools in BRENDA. Nucleic Acids Res. 2017;45: D380–D388.

Davidi

Noor

Liebermeister

et al . Global characterization of in vivo enzyme catalytic rates and their correspondence to in vitro kcat measurements. Proc Natl Acad Sci U S A. 2016;113:3401–3406.

van Eunen

Kiewiet

Westerhoff

Bakker

BM.

Testing biochemistry revisited: how in vivo metabolism can be understood from in vitro enzyme kinetics. PLoS Comput Biol. 2012;8:e1002483.

10.

Tominaga

Horton

Inference of scale-free networks from gene expression time series. J Bioinfor Comp Biol. 2006;4:503–514.

11.

Storn

Price

KV.

Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces. J Global Optim. 1997;11:341–359.

12.

Kawaguchi

Uematsu

Ogino

et al . Simultaneous saccharification and fermentation of kraft pulp by recombinant Escherichia coli for phenyllactic acid production. Biochem Eng J. 2014;88:188–194.

13.

Kawaguchi

Teramura

Uematsu

et al . Phenyllactic acid production by simultaneous saccharification and fermentation of pretreated sorghum bagasse. Biores Technol. 2015;182:169–178.

14.

Koma

Yamanaka

Moriyoshi

Ohmoto

Sakai

Production of aromatic compounds by metabolically engineered Escherichia coli with an expanded shikimate pathway. Appl Envoron Microbiol. 2012;78:6203–6216.

15.

Liu

S-P

Xiao

M-R

Zhang

et al . Production of L-phenylalanine from glucose by metabolic engineering of wild type Escherichia coli W3110. Process Biochem. 2013;48:413–419.